Line data Source code
1 : !=======================================================================
2 : !
3 : ! Reads and interpolates forcing data for atmosphere and ocean quantities.
4 : !
5 : ! authors: Elizabeth C. Hunke and William H. Lipscomb, LANL
6 : !
7 : module icedrv_forcing
8 :
9 : use icedrv_kinds
10 : use icedrv_domain_size, only: nx
11 : use icedrv_calendar, only: time, nyr, dayyr, mday, month, secday
12 : use icedrv_calendar, only: daymo, daycal, dt, yday, sec
13 : use icedrv_constants, only: nu_diag, nu_forcing, nu_open_clos
14 : use icedrv_constants, only: c0, c1, c2, c10, c100, p5, c4, c24
15 : use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted
16 : use icepack_intfc, only: icepack_query_parameters
17 : use icepack_intfc, only: icepack_sea_freezing_temperature
18 : use icepack_intfc, only: icepack_init_wave
19 : use icedrv_system, only: icedrv_system_abort
20 : use icedrv_flux, only: zlvl, Tair, potT, rhoa, uatm, vatm, wind, &
21 : strax, stray, fsw, swvdr, swvdf, swidr, swidf, Qa, flw, frain, &
22 : fsnow, sst, sss, uocn, vocn, qdp, hmix, Tf, opening, closing, sstdat
23 :
24 : implicit none
25 : private
26 : public :: init_forcing, get_forcing, interp_coeff, &
27 : interp_coeff_monthly, get_wave_spec
28 :
29 : integer (kind=int_kind), parameter :: &
30 : ntime = 8760 ! number of data points in time
31 :
32 : integer (kind=int_kind), public :: &
33 : ycycle , & ! number of years in forcing cycle
34 : fyear_init , & ! first year of data in forcing cycle
35 : fyear , & ! current year in forcing cycle
36 : fyear_final ! last year in cycle
37 :
38 : real (kind=dbl_kind), dimension(ntime) :: &
39 : fsw_data, & ! field values at temporal data points
40 : cldf_data, &
41 : fsnow_data, &
42 : Tair_data, &
43 : uatm_data, &
44 : vatm_data, &
45 : wind_data, &
46 : strax_data, &
47 : stray_data, &
48 : rhum_data, &
49 : Qa_data, &
50 : rhoa_data, &
51 : potT_data, &
52 : flw_data, &
53 : qdp_data, &
54 : sst_data, &
55 : sss_data, &
56 : uocn_data, &
57 : vocn_data, &
58 : frain_data, &
59 : swvdr_data, &
60 : swvdf_data, &
61 : swidr_data, &
62 : swidf_data, &
63 : zlvl_data, &
64 : hmix_data
65 :
66 : real (kind=dbl_kind), dimension(nx) :: &
67 : sst_temp
68 :
69 : real (kind=dbl_kind), dimension(ntime) :: &
70 : open_data, &
71 : clos_data
72 :
73 : character(char_len), public :: &
74 : atm_data_format, & ! 'bin'=binary or 'nc'=netcdf
75 : ocn_data_format, & ! 'bin'=binary or 'nc'=netcdf
76 : bgc_data_format, & ! 'bin'=binary or 'nc'=netcdf
77 : atm_data_type, & ! 'default', 'clim', 'CFS'
78 : ocn_data_type, & ! 'default', 'SHEBA'
79 : bgc_data_type, & ! 'default', 'ISPOL', 'NICE'
80 : atm_data_file, & ! atmospheric forcing data file
81 : ocn_data_file, & ! ocean forcing data file
82 : ice_data_file, & ! ice forcing data file
83 : bgc_data_file, & ! biogeochemistry forcing data file
84 : precip_units ! 'mm_per_month', 'mm_per_sec', 'mks'
85 :
86 : character(char_len_long), public :: &
87 : data_dir ! top directory for forcing data
88 :
89 : real (kind=dbl_kind), parameter, public :: &
90 : frcvdr = 0.28_dbl_kind, & ! frac of incoming sw in vis direct band
91 : frcvdf = 0.24_dbl_kind, & ! frac of incoming sw in vis diffuse band
92 : frcidr = 0.31_dbl_kind, & ! frac of incoming sw in near IR direct band
93 : frcidf = 0.17_dbl_kind ! frac of incoming sw in near IR diffuse band
94 :
95 : logical (kind=log_kind), public :: &
96 : oceanmixed_ice , & ! if true, use internal ocean mixed layer
97 : restore_ocn ! restore sst if true
98 :
99 : real (kind=dbl_kind), public :: &
100 : trest, & ! restoring time scale (sec)
101 : trestore ! restoring time scale (days)
102 :
103 : !=======================================================================
104 :
105 : contains
106 :
107 : !=======================================================================
108 :
109 42 : subroutine init_forcing
110 :
111 : ! Determine the current and final year of the forcing cycle based on
112 : ! namelist input; initialize the forcing data.
113 :
114 : integer (kind=int_kind) :: &
115 : i ! index
116 :
117 : character(len=*), parameter :: subname='(init_forcing)'
118 :
119 42 : fyear = fyear_init + mod(nyr-1,ycycle) ! current year
120 42 : fyear_final = fyear_init + ycycle - 1 ! last year in forcing cycle
121 :
122 42 : write (nu_diag,*) ' Initial forcing data year = ',fyear_init
123 42 : write (nu_diag,*) ' Final forcing data year = ',fyear_final
124 :
125 : !-------------------------------------------------------------------
126 : ! Initialize forcing data to default values
127 : !-------------------------------------------------------------------
128 :
129 : ! many default forcing values are set in init_flux_atm
130 42 : i = 1 ! use first grid box value
131 :
132 367962 : zlvl_data(:) = zlvl (i) ! atmospheric data level (m)
133 367962 : Tair_data(:) = Tair (i) ! air temperature (K)
134 367962 : potT_data(:) = potT (i) ! air potential temperature (K)
135 367962 : rhoa_data(:) = rhoa (i) ! air density (kg/m^3)
136 367962 : uatm_data(:) = uatm (i) ! wind velocity components (m/s)
137 367962 : vatm_data(:) = vatm (i)
138 367962 : wind_data(:) = wind (i) ! wind speed (m/s)
139 367962 : strax_data(:) = strax(i) ! wind stress components (N/m^2)
140 367962 : stray_data(:) = stray(i)
141 367962 : fsw_data(:) = fsw (i) ! incoming shortwave radiation (W/m^2)
142 367962 : swvdr_data(:) = swvdr(i) ! sw down, visible, direct (W/m^2)
143 367962 : swvdf_data(:) = swvdf(i) ! sw down, visible, diffuse (W/m^2)
144 367962 : swidr_data(:) = swidr(i) ! sw down, near IR, direct (W/m^2)
145 367962 : swidf_data(:) = swidf(i) ! sw down, near IR, diffuse (W/m^2)
146 367962 : Qa_data(:) = Qa (i) ! specific humidity (kg/kg)
147 367962 : flw_data(:) = flw (i) ! incoming longwave radiation (W/m^2)
148 367962 : frain_data(:) = frain(i) ! rainfall rate (kg/m^2 s)
149 367962 : fsnow_data(:) = fsnow(i) ! snowfall rate (kg/m^2 s)
150 367962 : qdp_data(:) = qdp (i) ! deep ocean heat flux (W/m^2)
151 367962 : sss_data(:) = sss (i) ! sea surface salinity
152 367962 : uocn_data(:) = uocn (i) ! ocean current components (m/s)
153 367962 : vocn_data(:) = vocn (i)
154 42 : cldf_data(:) = c0 ! cloud fraction
155 :
156 42 : if (trim(atm_data_type(1:4)) == 'CFS') call atm_CFS
157 42 : if (trim(atm_data_type(1:4)) == 'clim') call atm_climatological
158 42 : if (trim(atm_data_type(1:5)) == 'ISPOL') call atm_ISPOL
159 42 : if (trim(atm_data_type(1:4)) == 'NICE') call atm_NICE
160 42 : if (trim(ocn_data_type(1:5)) == 'SHEBA') call ice_open_clos
161 :
162 42 : if (restore_ocn) then
163 11 : if (trestore == 0) then
164 6 : trest = dt ! use data instantaneously
165 : else
166 5 : trest = real(trestore,kind=dbl_kind) * secday ! seconds
167 : end if
168 96371 : sst_data(:) = sstdat(i) ! default may be overwritten below
169 : else
170 271591 : sst_data(:) = sst (i) ! default or restart value if not restoring
171 : endif
172 :
173 42 : if (trim(ocn_data_type(1:5)) == 'ISPOL') call ocn_ISPOL
174 42 : if (trim(ocn_data_type(1:4)) == 'NICE') call ocn_NICE
175 :
176 : call prepare_forcing (Tair_data, fsw_data, &
177 : cldf_data, &
178 : frain_data, fsnow_data, &
179 : Qa_data, rhoa_data, &
180 : uatm_data, vatm_data, &
181 : strax_data, stray_data, &
182 : zlvl_data, wind_data, &
183 : swvdr_data, swvdf_data, &
184 : swidr_data, swidf_data, &
185 42 : potT_data)
186 :
187 42 : end subroutine init_forcing
188 :
189 : !=======================================================================
190 :
191 318408 : subroutine get_forcing(timestep)
192 :
193 : !ECH notes
194 : ! We will probably need to send in the time and work out what the data
195 : ! time slice is, instead of sending in the timestep. This currently assumes
196 : ! the time step and the data both start Jan 1.
197 :
198 : integer (kind=int_kind), intent(in) :: &
199 : timestep ! time step index
200 :
201 : integer (kind=int_kind) :: &
202 : i , & ! data index
203 : recslot , & ! spline slot for current record
204 : midmonth ! middle day of month
205 :
206 : integer (kind=int_kind) :: &
207 : mlast, mnext ! indices of bracketing time slices
208 :
209 : real (kind=dbl_kind) :: &
210 109812 : c1intp, c2intp ! interpolation coefficients
211 :
212 : integer (kind=int_kind) :: &
213 : recnum , & ! record number for current data value
214 : maxrec , & ! maximum number of data records
215 : dataloc , & ! = 1 for data located in middle of time interval
216 : ! = 2 for date located at end of time interval
217 : offndy ! Julian day of first data record
218 :
219 : real (kind=dbl_kind) :: &
220 109812 : sec6hr , & ! number of seconds in 6 hours
221 109812 : sec1hr , & ! number of seconds in 1 hour
222 109812 : offset ! time to first data record since 1 Jan (s)
223 :
224 : character(len=*), parameter :: subname='(get_forcing)'
225 :
226 318408 : if (trim(atm_data_type) == 'CFS') then
227 : ! calculate data index corresponding to current timestep
228 24132 : i = mod(timestep-1,ntime)+1 ! repeat forcing cycle
229 24132 : mlast = i
230 24132 : mnext = mlast
231 24132 : c1intp = c1
232 24132 : c2intp = c0
233 :
234 : ! fill all grid boxes with the same forcing data
235 120660 : Tair (:) = c1intp * Tair_data(mlast) + c2intp * Tair_data(mnext)
236 120660 : Qa (:) = c1intp * Qa_data(mlast) + c2intp * Qa_data(mnext)
237 120660 : uatm (:) = c1intp * uatm_data(mlast) + c2intp * uatm_data(mnext)
238 120660 : vatm (:) = c1intp * vatm_data(mlast) + c2intp * vatm_data(mnext)
239 120660 : fsnow(:) = c1intp * fsnow_data(mlast) + c2intp * fsnow_data(mnext)
240 120660 : flw (:) = c1intp * flw_data(mlast) + c2intp * flw_data(mnext)
241 120660 : fsw (:) = c1intp * fsw_data(mlast) + c2intp * fsw_data(mnext)
242 :
243 : ! derived (or not otherwise set)
244 120660 : potT (:) = c1intp * potT_data(mlast) + c2intp * potT_data(mnext)
245 120660 : wind (:) = c1intp * wind_data(mlast) + c2intp * wind_data(mnext)
246 120660 : strax(:) = c1intp * strax_data(mlast) + c2intp * strax_data(mnext)
247 120660 : stray(:) = c1intp * stray_data(mlast) + c2intp * stray_data(mnext)
248 120660 : rhoa (:) = c1intp * rhoa_data(mlast) + c2intp * rhoa_data(mnext)
249 120660 : frain(:) = c1intp * frain_data(mlast) + c2intp * frain_data(mnext)
250 120660 : swvdr(:) = c1intp * swvdr_data(mlast) + c2intp * swvdr_data(mnext)
251 120660 : swvdf(:) = c1intp * swvdf_data(mlast) + c2intp * swvdf_data(mnext)
252 120660 : swidr(:) = c1intp * swidr_data(mlast) + c2intp * swidr_data(mnext)
253 120660 : swidf(:) = c1intp * swidf_data(mlast) + c2intp * swidf_data(mnext)
254 :
255 294276 : elseif (trim(atm_data_type) == 'clim') then
256 199668 : midmonth = 15 ! assume data is given on 15th of every month
257 199668 : recslot = 1 ! latter half of month
258 199668 : if (mday < midmonth) recslot = 2 ! first half of month
259 199668 : if (recslot == 1) then
260 106896 : mlast = month
261 106896 : mnext = mod(month ,12) + 1
262 : else ! recslot = 2
263 92772 : mlast = mod(month+10,12) + 1
264 92772 : mnext = month
265 : endif
266 199668 : call interp_coeff_monthly(recslot, c1intp, c2intp)
267 :
268 : ! fill all grid boxes with the same forcing data
269 998340 : Tair (:) = c1intp * Tair_data(mlast) + c2intp * Tair_data(mnext)
270 998340 : Qa (:) = c1intp * Qa_data(mlast) + c2intp * Qa_data(mnext)
271 998340 : uatm (:) = c1intp * uatm_data(mlast) + c2intp * uatm_data(mnext)
272 998340 : vatm (:) = c1intp * vatm_data(mlast) + c2intp * vatm_data(mnext)
273 998340 : fsnow(:) = c1intp * fsnow_data(mlast) + c2intp * fsnow_data(mnext)
274 998340 : flw (:) = c1intp * flw_data(mlast) + c2intp * flw_data(mnext)
275 998340 : fsw (:) = c1intp * fsw_data(mlast) + c2intp * fsw_data(mnext)
276 :
277 : ! derived (or not otherwise set)
278 998340 : potT (:) = c1intp * potT_data(mlast) + c2intp * potT_data(mnext)
279 998340 : wind (:) = c1intp * wind_data(mlast) + c2intp * wind_data(mnext)
280 998340 : strax(:) = c1intp * strax_data(mlast) + c2intp * strax_data(mnext)
281 998340 : stray(:) = c1intp * stray_data(mlast) + c2intp * stray_data(mnext)
282 998340 : rhoa (:) = c1intp * rhoa_data(mlast) + c2intp * rhoa_data(mnext)
283 998340 : frain(:) = c1intp * frain_data(mlast) + c2intp * frain_data(mnext)
284 998340 : swvdr(:) = c1intp * swvdr_data(mlast) + c2intp * swvdr_data(mnext)
285 998340 : swvdf(:) = c1intp * swvdf_data(mlast) + c2intp * swvdf_data(mnext)
286 998340 : swidr(:) = c1intp * swidr_data(mlast) + c2intp * swidr_data(mnext)
287 998340 : swidf(:) = c1intp * swidf_data(mlast) + c2intp * swidf_data(mnext)
288 :
289 94608 : elseif (trim(atm_data_type) == 'ISPOL') then
290 :
291 20148 : offndy = 0 ! first data record (Julian day)
292 20148 : offset = real(offndy,dbl_kind)*secday
293 20148 : dataloc = 1 ! data located at middle of interval
294 20148 : maxrec = 365
295 20148 : recslot = 2
296 20148 : recnum = mod(int(yday)+maxrec-offndy-1,maxrec)+1
297 20148 : mlast = mod(recnum+maxrec-2,maxrec) + 1
298 20148 : mnext = mod(recnum-1, maxrec) + 1
299 : call interp_coeff (recnum, recslot, secday, dataloc, &
300 20148 : c1intp, c2intp, offset)
301 :
302 100740 : Tair (:) = c1intp * Tair_data(mlast) + c2intp * Tair_data(mnext)
303 100740 : Qa (:) = c1intp * Qa_data(mlast) + c2intp * Qa_data(mnext)
304 100740 : uatm (:) = c1intp * uatm_data(mlast) + c2intp * uatm_data(mnext)
305 100740 : vatm (:) = c1intp * vatm_data(mlast) + c2intp * vatm_data(mnext)
306 100740 : fsnow(:) = c1intp * fsnow_data(mlast) + c2intp * fsnow_data(mnext)
307 :
308 : ! derived (or not otherwise set)
309 100740 : potT (:) = c1intp * potT_data(mlast) + c2intp * potT_data(mnext)
310 100740 : wind (:) = c1intp * wind_data(mlast) + c2intp * wind_data(mnext)
311 100740 : strax(:) = c1intp * strax_data(mlast) + c2intp * strax_data(mnext)
312 100740 : stray(:) = c1intp * stray_data(mlast) + c2intp * stray_data(mnext)
313 100740 : rhoa (:) = c1intp * rhoa_data(mlast) + c2intp * rhoa_data(mnext)
314 100740 : frain(:) = c1intp * frain_data(mlast) + c2intp * frain_data(mnext)
315 :
316 20148 : sec6hr = secday/c4; ! seconds in 6 hours
317 20148 : offndy = 0
318 20148 : maxrec = 1460
319 20148 : recnum = 4*int(yday) - 3 + int(real(sec,kind=dbl_kind)/sec6hr)
320 20148 : recnum = mod(recnum+maxrec-4*offndy-1,maxrec)+1 ! data begins on 16 June 2004
321 20148 : recslot = 2
322 20148 : mlast = mod(recnum+maxrec-2,maxrec) + 1
323 20148 : mnext = mod(recnum-1, maxrec) + 1
324 : call interp_coeff (recnum, recslot, sec6hr, dataloc, &
325 20148 : c1intp, c2intp, offset)
326 :
327 100740 : fsw (:) = c1intp * fsw_data(mlast) + c2intp * fsw_data(mnext)
328 100740 : flw (:) = c1intp * flw_data(mlast) + c2intp * flw_data(mnext)
329 :
330 : ! derived
331 100740 : swvdr(:) = c1intp * swvdr_data(mlast) + c2intp * swvdr_data(mnext)
332 100740 : swvdf(:) = c1intp * swvdf_data(mlast) + c2intp * swvdf_data(mnext)
333 100740 : swidr(:) = c1intp * swidr_data(mlast) + c2intp * swidr_data(mnext)
334 100740 : swidf(:) = c1intp * swidf_data(mlast) + c2intp * swidf_data(mnext)
335 :
336 74460 : elseif (trim(atm_data_type) == 'NICE') then
337 :
338 2064 : offndy = 0 ! first data record (Julian day)
339 2064 : offset = real(offndy,dbl_kind)*secday
340 2064 : dataloc = 1 ! data located in middle of interval
341 2064 : maxrec = 365
342 2064 : recslot = 2
343 2064 : recnum = mod(int(yday)+maxrec-offndy-1,maxrec)+1
344 2064 : mlast = mod(recnum+maxrec-2,maxrec) + 1
345 2064 : mnext = mod(recnum-1, maxrec) + 1
346 : call interp_coeff (recnum, recslot, secday, dataloc, &
347 2064 : c1intp, c2intp, offset)
348 :
349 10320 : Tair (:) = c1intp * Tair_data(mlast) + c2intp * Tair_data(mnext)
350 10320 : Qa (:) = c1intp * Qa_data(mlast) + c2intp * Qa_data(mnext)
351 10320 : uatm (:) = c1intp * uatm_data(mlast) + c2intp * uatm_data(mnext)
352 10320 : vatm (:) = c1intp * vatm_data(mlast) + c2intp * vatm_data(mnext)
353 10320 : fsnow(:) = c1intp * fsnow_data(mlast) + c2intp * fsnow_data(mnext)
354 :
355 : ! derived (or not otherwise set)
356 10320 : potT (:) = c1intp * potT_data(mlast) + c2intp * potT_data(mnext)
357 10320 : wind (:) = c1intp * wind_data(mlast) + c2intp * wind_data(mnext)
358 10320 : strax(:) = c1intp * strax_data(mlast) + c2intp * strax_data(mnext)
359 10320 : stray(:) = c1intp * stray_data(mlast) + c2intp * stray_data(mnext)
360 10320 : rhoa (:) = c1intp * rhoa_data(mlast) + c2intp * rhoa_data(mnext)
361 10320 : frain(:) = c1intp * frain_data(mlast) + c2intp * frain_data(mnext)
362 :
363 2064 : sec6hr = secday/c4; ! seconds in 6 hours
364 2064 : maxrec = 1460
365 2064 : dataloc = 2 ! data located at end of interval
366 2064 : recnum = 4*int(yday) - 3 + int(real(sec,kind=dbl_kind)/sec6hr)
367 2064 : recnum = mod(recnum+maxrec-4*offndy-1,maxrec)+1
368 2064 : recslot = 2
369 2064 : mlast = mod(recnum+maxrec-2,maxrec) + 1
370 2064 : mnext = mod(recnum-1, maxrec) + 1
371 : call interp_coeff (recnum, recslot, sec6hr, dataloc, &
372 2064 : c1intp, c2intp, offset)
373 :
374 10320 : fsw (:) = c1intp * fsw_data(mlast) + c2intp * fsw_data(mnext)
375 10320 : flw (:) = c1intp * flw_data(mlast) + c2intp * flw_data(mnext)
376 :
377 : ! derived
378 10320 : swvdr(:) = c1intp * swvdr_data(mlast) + c2intp * swvdr_data(mnext)
379 10320 : swvdf(:) = c1intp * swvdf_data(mlast) + c2intp * swvdf_data(mnext)
380 10320 : swidr(:) = c1intp * swidr_data(mlast) + c2intp * swidr_data(mnext)
381 10320 : swidf(:) = c1intp * swidf_data(mlast) + c2intp * swidf_data(mnext)
382 :
383 : endif
384 :
385 : ! possible bug: is the ocean data also offset to the beginning of the field campaigns?
386 :
387 318408 : if (trim(ocn_data_type) == 'ISPOL') then
388 :
389 20148 : midmonth = 15 ! assume data is given on 15th of every month
390 20148 : recslot = 1 ! latter half of month
391 20148 : if (mday < midmonth) recslot = 2 ! first half of month
392 20148 : if (recslot == 1) then
393 11015 : mlast = month
394 11015 : mnext = mod(month ,12) + 1
395 : else ! recslot = 2
396 9133 : mlast = mod(month+10,12) + 1
397 9133 : mnext = month
398 : endif
399 20148 : call interp_coeff_monthly(recslot, c1intp, c2intp)
400 :
401 100740 : sst_temp(:) = c1intp * sst_data(mlast) + c2intp * sst_data(mnext)
402 100740 : sss (:) = c1intp * sss_data(mlast) + c2intp * sss_data(mnext)
403 100740 : uocn (:) = c1intp * uocn_data(mlast) + c2intp * uocn_data(mnext)
404 100740 : vocn (:) = c1intp * vocn_data(mlast) + c2intp * vocn_data(mnext)
405 100740 : qdp (:) = c1intp * qdp_data(mlast) + c2intp * qdp_data(mnext)
406 :
407 298260 : elseif (trim(ocn_data_type) == 'NICE') then
408 :
409 2064 : dataloc = 2 ! data located at end of interval
410 2064 : maxrec = 365
411 2064 : recslot = 2
412 2064 : recnum = int(yday)
413 2064 : mlast = mod(recnum+maxrec-2,maxrec) + 1
414 2064 : mnext = mod(recnum-1, maxrec) + 1
415 2064 : call interp_coeff ( recnum, recslot, secday, dataloc, c1intp, c2intp)
416 :
417 10320 : sst_temp(:) = c1intp * sst_data(mlast) + c2intp * sst_data(mnext)
418 10320 : sss (:) = c1intp * sss_data(mlast) + c2intp * sss_data(mnext)
419 10320 : uocn (:) = c1intp * uocn_data(mlast) + c2intp * uocn_data(mnext)
420 10320 : vocn (:) = c1intp * vocn_data(mlast) + c2intp * vocn_data(mnext)
421 10320 : qdp (:) = c1intp * qdp_data(mlast) + c2intp * qdp_data(mnext)
422 10320 : hmix (:) = c1intp * hmix_data(mlast) + c2intp * hmix_data(mnext)
423 :
424 : else
425 :
426 : ! use default values for all other data fields
427 296196 : i = mod(timestep-1,ntime)+1 ! repeat forcing cycle
428 296196 : mlast = i
429 296196 : mnext = mlast
430 296196 : c1intp = c1
431 296196 : c2intp = c0
432 :
433 1480980 : sst_temp(:) = c1intp * sst_data(mlast) + c2intp * sst_data(mnext)
434 1480980 : sss (:) = c1intp * sss_data(mlast) + c2intp * sss_data(mnext)
435 1480980 : uocn (:) = c1intp * uocn_data(mlast) + c2intp * uocn_data(mnext)
436 1480980 : vocn (:) = c1intp * vocn_data(mlast) + c2intp * vocn_data(mnext)
437 1480980 : qdp (:) = c1intp * qdp_data(mlast) + c2intp * qdp_data(mnext)
438 :
439 : endif
440 :
441 318408 : call finish_ocn_forcing(sst_temp)
442 :
443 : ! Lindsay SHEBA open/close dataset is hourly
444 318408 : if (trim(ocn_data_type) == 'SHEBA') then
445 :
446 272064 : sec1hr = secday/c24 ! seconds in 1 hour
447 272064 : maxrec = ntime
448 272064 : recnum = 24*int(yday) - 23 + int(real(sec,kind=dbl_kind)/sec1hr)
449 272064 : recslot = 2
450 272064 : dataloc = 1 ! data located at middle of interval
451 272064 : mlast = mod(recnum+maxrec-2,maxrec) + 1
452 272064 : mnext = mod(recnum-1, maxrec) + 1
453 272064 : call interp_coeff ( recnum, recslot, sec1hr, dataloc, c1intp, c2intp)
454 :
455 1360320 : opening(:) = c1intp * open_data(mlast) + c2intp * open_data(mnext)
456 1360320 : closing(:) = -(c1intp * clos_data(mlast) + c2intp * clos_data(mnext))
457 :
458 : endif
459 :
460 318408 : end subroutine get_forcing
461 :
462 : !=======================================================================
463 :
464 25 : subroutine atm_climatological
465 :
466 : real (kind=dbl_kind), dimension(12) :: &
467 : fsw_clim, & ! field values at temporal data points
468 : flw_clim, &
469 : Tair_clim, &
470 : wind_clim, &
471 : rhum_clim, &
472 : fsnow_clim
473 :
474 : real (kind=dbl_kind) :: &
475 8 : Tffresh, qqqice, TTTice, rhos
476 :
477 : character(len=*), parameter :: subname='(atm_climatological)'
478 :
479 : ! Ice station meteorology from Lindsay (1998, J. Climate), Table 1, p. 325
480 : ! zlvl = c2 ! 2-m temperatures and wind speed
481 :
482 : data fsw_clim / 0.0d0, 1.2d0, 31.5d0, 146.0d0, 263.3d0, 307.9d0, &
483 : 230.6d0, 134.7d0, 44.2d0, 2.6d0, 0.0d0, 0.0d0 /
484 : data flw_clim /164.0d0, 160.5d0, 164.1d0, 188.1d0, 245.2d0, 291.2d0, &
485 : 303.9d0, 297.0d0, 263.8d0, 210.9d0, 177.0d0, 166.0d0 /
486 : data Tair_clim /-31.4d0, -32.8d0, -31.6d0, -24.1d0, -11.0d0, -1.8d0, &
487 : -0.1d0, -1.4d0, -8.0d0, -19.5d0, -27.6d0, -31.1d0 /
488 : data rhum_clim / 78.7d0, 78.4d0, 79.6d0, 82.1d0, 86.5d0, 91.7d0, &
489 : 95.1d0, 94.3d0, 90.7d0, 83.8d0, 80.1d0, 78.7d0 /
490 : data wind_clim / 4.4d0, 4.0d0, 4.0d0, 3.9d0, 3.9d0, 4.2d0, &
491 : 4.1d0, 4.2d0, 4.5d0, 4.2d0, 3.9d0, 4.0d0 /
492 : ! data shf_clim / 9.9d0, 8.4d0, 6.6d0, 0.1d0, -5.8d0, -1.6d0, &
493 : ! 2.2d0, 1.2d0, 0.5d0, 2.0d0, 5.6d0, 7.0d0 /
494 : ! data lhf_clim / 1.3d0, 1.1d0, 1.1d0, 0.0d0, -5.9d0, -10.3d0, &
495 : ! -6.5d0, -6.7d0, -3.9d0, -0.1d0, 1.0d0, 1.1d0 /
496 :
497 : ! Semtner (1976, JPO) snowfall spec., p. 383 in m/s snow volume (.4 m/yr)
498 : data fsnow_clim/ 3.17d-9, 3.17d-9, 3.17d-9, 3.17d-9, 1.90d-8, 0.0d0, &
499 : 0.0d0, 1.63d-8, 4.89d-8, 4.89d-8, 3.17d-9, 3.17d-9 /
500 :
501 : !-----------------------------------------------------------------
502 : ! query icepack values
503 : !-----------------------------------------------------------------
504 :
505 : call icepack_query_parameters(Tffresh_out=Tffresh, qqqice_out=qqqice, &
506 25 : TTTice_out=TTTice, rhos_out=rhos)
507 25 : call icepack_warnings_flush(nu_diag)
508 25 : if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
509 0 : file=__FILE__,line= __LINE__)
510 :
511 : !-----------------------------------------------------------------
512 :
513 325 : fsw_data (1:12) = fsw_clim (1:12)
514 325 : flw_data (1:12) = flw_clim (1:12)
515 325 : rhum_data (1:12) = rhum_clim (1:12)
516 325 : wind_data (1:12) = wind_clim (1:12)
517 :
518 325 : rhoa_data (1:12) = 1.275_dbl_kind ! air density (kg/m^3)
519 325 : Tair_data (1:12) = Tair_clim (1:12) + Tffresh
520 325 : uatm_data (1:12) = wind_clim (1:12)
521 325 : vatm_data (1:12) = c0
522 :
523 : ! Qa = rhum * saturation humidity (1.275 kg/m^3 = air density)
524 : Qa_data (1:12) = (rhum_clim(1:12)/c100)*qqqice &
525 325 : * exp(-TTTice/Tair_data(1:12))/rhoa_data(1:12)
526 :
527 325 : fsnow_data(1:12) = rhos*fsnow_clim(1:12) ! convert vol -> mass flux
528 325 : frain_data(1:12) = c0
529 :
530 : ! 6 W/m2 warming of mixed layer from deep ocean
531 219025 : qdp_data(:) = -6.0 ! 2 W/m2 from deep + 4 W/m2 counteracting larger
532 : ! SH+LH with bulk transfer than in MU 71
533 :
534 25 : end subroutine atm_climatological
535 :
536 : !=======================================================================
537 :
538 3 : subroutine atm_CFS
539 :
540 : integer (kind=int_kind) :: &
541 : nt ! loop index
542 :
543 : real (kind=dbl_kind) :: &
544 1 : dlwsfc, & ! downwelling longwave (W/m2)
545 1 : dswsfc, & ! downwelling shortwave (W/m2)
546 1 : windu10, & ! wind components (m/s)
547 1 : windv10, & !
548 1 : temp2m, & ! 2m air temperature (K)
549 1 : spechum ,& ! specific humidity (kg/kg)
550 1 : precip ! precipitation (kg/m2/s)
551 :
552 : character (char_len_long) string1
553 : character (char_len_long) filename
554 : character(len=*), parameter :: subname='(atm_CFS)'
555 :
556 : ! atm_data_file = 'cfsv2_2015_220_70_01hr.txt'
557 3 : filename = trim(data_dir)//'/CFS/'//trim(atm_data_file)
558 :
559 3 : write (nu_diag,*) 'Reading ',filename
560 :
561 3 : open (nu_forcing, file=filename, form='formatted')
562 3 : read (nu_forcing, *) string1 ! headers
563 3 : read (nu_forcing, *) string1 ! units
564 :
565 26283 : do nt = 1, ntime
566 : read (nu_forcing, '(6(f10.5,1x),2(f10.8,1x))') &
567 26280 : dswsfc, dlwsfc, windu10, windv10, temp2m, spechum, precip
568 :
569 26280 : flw_data(nt) = dlwsfc
570 26280 : fsw_data(nt) = dswsfc
571 26280 : uatm_data(nt) = windu10
572 26280 : vatm_data(nt) = windv10
573 26280 : Tair_data(nt) = temp2m
574 26280 : potT_data(nt) = temp2m
575 26280 : Qa_data(nt) = spechum
576 26283 : fsnow_data(nt) = precip
577 : enddo
578 :
579 3 : close (nu_forcing)
580 :
581 3 : end subroutine atm_CFS
582 :
583 : !=======================================================================
584 :
585 42 : subroutine prepare_forcing (Tair, fsw, &
586 : cldf, &
587 : frain, fsnow, &
588 : Qa, rhoa, &
589 : uatm, vatm, &
590 : strax, stray, &
591 : zlvl, wind, &
592 : swvdr, swvdf, &
593 : swidr, swidf, &
594 : potT)
595 :
596 : ! this routine acts on the data fields prior to interpolation
597 :
598 : real (kind=dbl_kind), dimension(ntime), &
599 : intent(inout) :: &
600 : Tair , & ! air temperature (K)
601 : fsw , & ! incoming shortwave radiation (W/m^2)
602 : cldf , & ! cloud fraction
603 : frain , & ! rainfall rate (kg/m^2 s)
604 : fsnow , & ! snowfall rate (kg/m^2 s)
605 : Qa , & ! specific humidity (kg/kg)
606 : rhoa , & ! air density (kg/m^3)
607 : uatm , & ! wind velocity components (m/s)
608 : vatm , &
609 : strax , & ! wind stress components (N/m^2)
610 : stray , &
611 : zlvl , & ! atm level height (m)
612 : wind , & ! wind speed (m/s)
613 : ! flw , & ! incoming longwave radiation (W/m^2)
614 : swvdr , & ! sw down, visible, direct (W/m^2)
615 : swvdf , & ! sw down, visible, diffuse (W/m^2)
616 : swidr , & ! sw down, near IR, direct (W/m^2)
617 : swidf , & ! sw down, near IR, diffuse (W/m^2)
618 : potT ! air potential temperature (K)
619 :
620 : ! local variables
621 :
622 : integer (kind=int_kind) :: &
623 : nt
624 :
625 : logical (kind=log_kind) :: &
626 : calc_strair
627 :
628 : real (kind=dbl_kind), parameter :: &
629 : lapse_rate = 0.0065_dbl_kind ! (K/m) lapse rate over sea level
630 :
631 : real (kind=dbl_kind) :: &
632 30 : precip_factor, zlvl0, &
633 15 : Tffresh
634 :
635 : character(len=*), parameter :: subname='(prepare_forcing)'
636 :
637 42 : zlvl0 = c10 ! default
638 :
639 : !-----------------------------------------------------------------
640 : ! query icepack values
641 : !-----------------------------------------------------------------
642 :
643 42 : call icepack_query_parameters(calc_strair_out=calc_strair)
644 42 : call icepack_query_parameters(Tffresh_out=Tffresh)
645 42 : call icepack_warnings_flush(nu_diag)
646 42 : if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
647 0 : file=__FILE__,line= __LINE__)
648 :
649 : !-----------------------------------------------------------------
650 : ! convert precipitation units to kg/m^2 s
651 : !-----------------------------------------------------------------
652 42 : if (trim(precip_units) == 'mm_per_month') then
653 0 : precip_factor = 12._dbl_kind/(secday*dayyr)
654 42 : elseif (trim(precip_units) == 'mm_per_day') then
655 0 : precip_factor = c1/secday
656 42 : elseif (trim(precip_units) == 'mm_per_sec' .or. &
657 : trim(precip_units) == 'mks') then
658 42 : precip_factor = c1 ! mm/sec = kg/m^2 s
659 : endif
660 :
661 367962 : do nt = 1, ntime
662 :
663 : !-----------------------------------------------------------------
664 : ! make sure interpolated values are physically realistic
665 : !-----------------------------------------------------------------
666 367920 : cldf (nt) = max(min(cldf(nt),c1),c0)
667 367920 : fsw (nt) = max(fsw(nt),c0)
668 367920 : fsnow(nt) = max(fsnow(nt),c0)
669 367920 : rhoa (nt) = max(rhoa(nt),c0)
670 367920 : Qa (nt) = max(Qa(nt),c0)
671 :
672 : !-----------------------------------------------------------------
673 : ! calculations specific to datasets
674 : !-----------------------------------------------------------------
675 :
676 367920 : if (trim(atm_data_type) == 'CFS') then
677 : ! precip is in kg/m^2/s
678 26280 : zlvl0 = c10
679 : ! downward longwave as in Parkinson and Washington (1979)
680 : ! call longwave_parkinson_washington(Tair(nt), cldf(nt), flw(nt))
681 :
682 341640 : elseif (trim(atm_data_type) == 'clim') then
683 : ! precip is in kg/m^2/s
684 219000 : zlvl0 = c2
685 :
686 122640 : elseif (trim(atm_data_type) == 'ISPOL' .or. &
687 : trim(atm_data_type) == 'NICE') then
688 43800 : Tair(nt) = Tair(nt) - lapse_rate*8.0_dbl_kind
689 :
690 : endif ! atm_data_type
691 :
692 : ! this longwave depends on the current ice aice and sst and so can not be
693 : ! computed ahead of time
694 : ! ! longwave based on Rosati and Miyakoda, JPO 18, p. 1607 (1988)
695 : ! call longwave_rosati_miyakoda(cldf(i,j), Tsfc(i,j), &
696 : ! aice(i,j), sst(i,j), &
697 : ! Qa(i,j), Tair(i,j), &
698 : ! hm(i,j), flw(i,j))
699 :
700 : !-----------------------------------------------------------------
701 : ! Compute other fields needed by model
702 : !-----------------------------------------------------------------
703 :
704 367920 : zlvl(nt) = zlvl0
705 367920 : potT(nt) = Tair(nt)
706 :
707 : ! divide shortwave into spectral bands
708 367920 : swvdr(nt) = fsw(nt)*frcvdr ! visible direct
709 367920 : swvdf(nt) = fsw(nt)*frcvdf ! visible diffuse
710 367920 : swidr(nt) = fsw(nt)*frcidr ! near IR direct
711 367920 : swidf(nt) = fsw(nt)*frcidf ! near IR diffuse
712 :
713 : ! precipitation
714 367920 : fsnow(nt) = fsnow(nt) * precip_factor
715 :
716 : ! determine whether precip is rain or snow
717 : ! HadGEM forcing provides separate snowfall and rainfall rather
718 : ! than total precipitation
719 : ! if (trim(atm_data_type) /= 'hadgem') then
720 367920 : frain(nt) = c0
721 367920 : if (Tair(nt) >= Tffresh) then
722 8093 : frain(nt) = fsnow(nt)
723 8093 : fsnow(nt) = c0
724 : endif
725 : ! endif
726 :
727 367962 : if (calc_strair) then
728 367920 : wind (nt) = sqrt(uatm(nt)**2 + vatm(nt)**2)
729 367920 : strax(nt) = c0
730 367920 : stray(nt) = c0
731 : ! else ! strax, stray, wind are read from files
732 : endif ! calc_strair
733 :
734 : enddo ! ntime
735 :
736 42 : end subroutine prepare_forcing
737 :
738 : !=======================================================================
739 :
740 219816 : subroutine interp_coeff_monthly (recslot, c1intp, c2intp)
741 :
742 : ! Compute coefficients for interpolating monthly data to current time step.
743 :
744 : integer (kind=int_kind), intent(in) :: &
745 : recslot ! slot (1 or 2) for current record
746 :
747 : real (kind=dbl_kind), intent(inout) :: &
748 : c1intp, c2intp ! interpolation coefficients
749 :
750 : ! local variables
751 :
752 : real (kind=dbl_kind) :: &
753 72708 : tt , & ! seconds elapsed in current year
754 72708 : t1, t2 ! seconds elapsed at month midpoint
755 :
756 : real (kind=dbl_kind) :: &
757 1090620 : daymid(0:13) ! month mid-points
758 :
759 : character(len=*), parameter :: subname='(interp_coeff_monthly)'
760 :
761 3077424 : daymid(1:13) = 14._dbl_kind ! time frame ends 0 sec into day 15
762 219816 : daymid(0) = 14._dbl_kind - daymo(12) ! Dec 15, 0 sec
763 :
764 : ! make time cyclic
765 219816 : tt = mod(time/secday,dayyr)
766 :
767 : ! Find neighboring times
768 :
769 219816 : if (recslot==2) then ! first half of month
770 101905 : t2 = daycal(month) + daymid(month) ! midpoint, current month
771 101905 : if (month == 1) then
772 12446 : t1 = daymid(0) ! Dec 15 (0 sec)
773 : else
774 89459 : t1 = daycal(month-1) + daymid(month-1) ! midpoint, previous month
775 : endif
776 : else ! second half of month
777 117911 : t1 = daycal(month) + daymid(month) ! midpoint, current month
778 117911 : t2 = daycal(month+1) + daymid(month+1)! day 15 of next month (0 sec)
779 : endif
780 :
781 : ! Compute coefficients
782 219816 : c1intp = (t2 - tt) / (t2 - t1)
783 219816 : c2intp = c1 - c1intp
784 :
785 219816 : end subroutine interp_coeff_monthly
786 :
787 : !=======================================================================
788 :
789 318552 : subroutine interp_coeff (recnum, recslot, secint, dataloc, &
790 : c1intp, c2intp, offset)
791 :
792 : ! Compute coefficients for interpolating data to current time step.
793 : ! Works for any data interval that divides evenly into a
794 : ! year (daily, 6-hourly, etc.)
795 : ! Use interp_coef_monthly for monthly data.
796 :
797 : integer (kind=int_kind), intent(in) :: &
798 : recnum , & ! record number for current data value
799 : recslot , & ! spline slot for current record
800 : dataloc ! = 1 for data located in middle of time interval
801 : ! = 2 for date located at end of time interval
802 :
803 : real (kind=dbl_kind), intent(in) :: &
804 : secint ! seconds in data interval
805 :
806 : real (kind=dbl_kind), intent(inout) :: &
807 : c1intp, c2intp ! interpolation coefficients
808 :
809 : real (kind=dbl_kind), intent(in), optional :: &
810 : offset ! amount of time data is offset (s)
811 :
812 : ! local variables
813 :
814 : real (kind=dbl_kind) :: &
815 109956 : secyr ! seconds in a year
816 :
817 : real (kind=dbl_kind) :: &
818 109956 : tt , & ! seconds elapsed in current year
819 109956 : t1, t2 , & ! seconds elapsed at data points
820 109956 : rcnum ! recnum => dbl_kind
821 :
822 : character(len=*), parameter :: subname='(interp_coeff)'
823 :
824 318552 : secyr = dayyr * secday ! seconds in a year
825 318552 : if (present(offset)) then
826 44424 : tt = mod(time-offset,secyr)
827 : else
828 274128 : tt = mod(time,secyr)
829 : endif
830 :
831 : ! Find neighboring times
832 318552 : rcnum = real(recnum,kind=dbl_kind)
833 318552 : if (recslot==2) then ! current record goes in slot 2
834 318552 : if (dataloc==1) then ! data located at middle of interval
835 314424 : t2 = (rcnum-p5)*secint
836 : else ! data located at end of interval
837 4128 : t2 = rcnum*secint
838 : endif
839 318552 : t1 = t2 - secint ! - 1 interval
840 : else ! recslot = 1
841 0 : if (dataloc==1) then ! data located at middle of interval
842 0 : t1 = (rcnum-p5)*secint
843 : else
844 0 : t1 = rcnum*secint ! data located at end of interval
845 : endif
846 0 : t2 = t1 + secint ! + 1 interval
847 : endif
848 :
849 : ! Compute coefficients
850 318552 : c1intp = abs((t2 - tt) / (t2 - t1))
851 318552 : c2intp = c1 - c1intp
852 :
853 318552 : end subroutine interp_coeff
854 :
855 : !=======================================================================
856 :
857 3 : subroutine atm_ISPOL
858 : ! there is data for 366 days, but we only use 365
859 :
860 : integer (kind=int_kind) :: &
861 : i ! index
862 :
863 : real (kind=dbl_kind), dimension(366) :: &
864 367 : tair , & ! air temperature
865 367 : qa , & ! specific humidity
866 367 : uatm , & ! wind velocity components
867 367 : vatm , &
868 367 : fsnow ! snowfall rate
869 : ! aday ! not used
870 :
871 : real (kind=dbl_kind), dimension(1464) :: &
872 1465 : fsw , & ! shortwave
873 1465 : flw ! longwave
874 : ! atime ! not used
875 :
876 : character (char_len_long) filename
877 :
878 : character(len=*), parameter :: subname='(atm_ISPOL)'
879 :
880 : ! atm_data_file = 'ISPOL_atm_forcing.txt'
881 3 : filename = trim(data_dir)//'/ISPOL_2004/'//trim(atm_data_file)
882 :
883 3 : write (nu_diag,*) 'Reading ',filename
884 :
885 3 : open (nu_forcing, file=filename, form='formatted')
886 :
887 3 : read(nu_forcing,*) tair
888 3 : read(nu_forcing,*) qa
889 3 : read(nu_forcing,*) fsw
890 3 : read(nu_forcing,*) flw
891 3 : read(nu_forcing,*) uatm
892 3 : read(nu_forcing,*) vatm
893 3 : read(nu_forcing,*) fsnow
894 : ! read(nu_forcing,*) aday
895 : ! read(nu_forcing,*) atime
896 :
897 1101 : do i = 1, 366 ! daily
898 1098 : Tair_data (i) = tair (i)
899 1098 : Qa_data (i) = qa (i)
900 1098 : uatm_data (i) = uatm (i)
901 1098 : vatm_data (i) = vatm (i)
902 1101 : fsnow_data(i) = fsnow(i)
903 : end do
904 4395 : do i = 1, 1464 ! 6hr, 1464/4=366 days
905 4392 : fsw_data (i) = fsw (i)
906 4395 : flw_data (i) = flw (i)
907 : end do
908 :
909 3 : close(nu_forcing)
910 :
911 3 : end subroutine atm_ISPOL
912 :
913 : !=======================================================================
914 :
915 2 : subroutine atm_NICE
916 : ! there is data for 366 days, but we only use 365
917 :
918 : integer (kind=int_kind) :: &
919 : i ! index
920 :
921 : real (kind=dbl_kind), dimension(366) :: &
922 734 : tair , & ! air temperature
923 734 : qa , & ! specific humidity
924 734 : uatm , & ! wind velocity components
925 734 : vatm , &
926 734 : fsnow ! snowfall rate
927 : ! aday ! not used
928 :
929 : real (kind=dbl_kind), dimension(1464) :: &
930 2930 : fsw , & ! shortwave
931 2930 : flw ! longwave
932 : ! atime ! not used
933 :
934 : character (char_len_long) filename
935 :
936 : character(len=*), parameter :: subname='(atm_NICE)'
937 :
938 : ! atm_data_file = 'NICE_atm_forcing.txt'
939 2 : filename = trim(data_dir)//'/NICE_2015/'//trim(atm_data_file)
940 :
941 2 : write (nu_diag,*) 'Reading ',filename
942 :
943 2 : open (nu_forcing, file=filename, form='formatted')
944 :
945 2 : read(nu_forcing,*) tair
946 2 : read(nu_forcing,*) qa
947 2 : read(nu_forcing,*) fsw
948 2 : read(nu_forcing,*) flw
949 2 : read(nu_forcing,*) uatm
950 2 : read(nu_forcing,*) vatm
951 2 : read(nu_forcing,*) fsnow
952 : ! read(nu_forcing,*) aday
953 : ! read(nu_forcing,*) atime
954 :
955 734 : do i = 1, 366
956 732 : Tair_data (i) = tair (i)
957 732 : Qa_data (i) = qa (i)
958 732 : uatm_data (i) = uatm (i)
959 732 : vatm_data (i) = vatm (i)
960 734 : fsnow_data(i) = fsnow(i)
961 : end do
962 2930 : do i = 1, 1464
963 2928 : fsw_data (i) = fsw (i)
964 2930 : flw_data (i) = flw (i)
965 : end do
966 :
967 2 : close(nu_forcing)
968 :
969 2 : end subroutine atm_NICE
970 :
971 : !=======================================================================
972 :
973 2 : subroutine ocn_NICE
974 :
975 : integer (kind=int_kind) :: &
976 : i
977 :
978 : real (kind=dbl_kind), dimension(365) :: &
979 732 : t , & ! sea surface temperature
980 732 : s , & ! sea surface salinity
981 732 : hblt, & ! mixed layer depth
982 732 : u , & ! ocean current, x
983 732 : v , & ! ocean current, y
984 732 : dhdx, & ! sea surface slope
985 732 : dhdy, & ! sea surface slope
986 732 : qdp ! deep ocean heat flux
987 :
988 : character (char_len_long) filename
989 :
990 : character(len=*), parameter :: subname='(ocn_NICE)'
991 :
992 : ! ocn_data_file = 'oceanmixed_daily_3.txt'
993 2 : filename = trim(data_dir)//'/NICE_2015/'//trim(ocn_data_file)
994 :
995 2 : write (nu_diag,*) 'Reading ',filename
996 :
997 2 : open (nu_forcing, file=filename, form='formatted')
998 :
999 2 : read(nu_forcing,*) t
1000 2 : read(nu_forcing,*) s
1001 2 : read(nu_forcing,*) hblt
1002 2 : read(nu_forcing,*) u
1003 2 : read(nu_forcing,*) v
1004 2 : read(nu_forcing,*) dhdx ! not used for Icepack
1005 2 : read(nu_forcing,*) dhdy ! not used for Icepack
1006 2 : read(nu_forcing,*) qdp
1007 :
1008 2 : close(nu_forcing)
1009 :
1010 732 : do i = 1, 365 ! daily
1011 730 : sst_data (i) = t (i)
1012 730 : sss_data (i) = s (i)
1013 730 : hmix_data(i) = hblt(i)
1014 730 : uocn_data(i) = u (i)
1015 730 : vocn_data(i) = v (i)
1016 732 : qdp_data (i) = qdp (i)
1017 : end do
1018 :
1019 2 : end subroutine ocn_NICE
1020 :
1021 : !=======================================================================
1022 :
1023 3 : subroutine ocn_ISPOL
1024 :
1025 : integer (kind=int_kind) :: &
1026 : i
1027 :
1028 : real (kind=dbl_kind), dimension(12) :: &
1029 13 : t , & ! sea surface temperature
1030 13 : s , & ! sea surface salinity
1031 13 : hblt, & ! mixed layer depth
1032 13 : u , & ! ocean current, x
1033 13 : v , & ! ocean current, y
1034 13 : dhdx, & ! sea surface slope
1035 13 : dhdy, & ! sea surface slope
1036 13 : qdp ! deep ocean heat flux
1037 :
1038 : character (char_len_long) filename
1039 :
1040 : character(len=*), parameter :: subname='(ocn_ISPOL)'
1041 :
1042 : ! ocn_data_file = 'pop_frc.gx1v3.051202_but_hblt_from_010815_ispol.txt'
1043 3 : filename = trim(data_dir)//'/ISPOL_2004/'//trim(ocn_data_file)
1044 :
1045 3 : write (nu_diag,*) 'Reading ',filename
1046 :
1047 3 : open (nu_forcing, file=filename, form='formatted')
1048 :
1049 3 : read(nu_forcing,*) t
1050 3 : read(nu_forcing,*) s
1051 3 : read(nu_forcing,*) hblt
1052 3 : read(nu_forcing,*) u
1053 3 : read(nu_forcing,*) v
1054 3 : read(nu_forcing,*) dhdx ! not used for Icepack
1055 3 : read(nu_forcing,*) dhdy ! not used for Icepack
1056 3 : read(nu_forcing,*) qdp
1057 :
1058 3 : close(nu_forcing)
1059 :
1060 39 : do i = 1, 12 ! monthly
1061 36 : sst_data (i) = t (i)
1062 36 : sss_data (i) = s (i)
1063 36 : hmix_data(i) = hblt(i)
1064 36 : uocn_data(i) = u (i)
1065 36 : vocn_data(i) = v (i)
1066 39 : qdp_data (i) = qdp (i)
1067 : end do
1068 :
1069 3 : end subroutine ocn_ISPOL
1070 :
1071 : !=======================================================================
1072 :
1073 318408 : subroutine finish_ocn_forcing(sst_temp)
1074 :
1075 : ! Compute ocean freezing temperature Tf based on tfrz_option
1076 : ! 'minus1p8' Tf = -1.8 C (default)
1077 : ! 'linear_salt' Tf = -depressT * sss
1078 : ! 'mushy' Tf conforms with mushy layer thermo (ktherm=2)
1079 :
1080 : real (kind=dbl_kind), dimension(nx), intent(in) :: &
1081 : sst_temp
1082 :
1083 : ! local variables
1084 :
1085 : integer (kind=int_kind) :: &
1086 : i ! horizontal indices
1087 :
1088 : character(len=*), parameter :: subname='(finish_ocn_forcing)'
1089 :
1090 1592040 : do i = 1, nx
1091 1273632 : sss (i) = max (sss(i), c0)
1092 1273632 : hmix(i) = max(hmix(i), c0)
1093 1273632 : Tf (i) = icepack_sea_freezing_temperature(sss(i))
1094 1592040 : if (restore_ocn) sst(i) = sst(i) + (sst_temp(i)-sst(i))*dt/trest
1095 : enddo
1096 318408 : call icepack_warnings_flush(nu_diag)
1097 318408 : if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
1098 0 : file=__FILE__,line= __LINE__)
1099 :
1100 318408 : end subroutine finish_ocn_forcing
1101 :
1102 : !=======================================================================
1103 :
1104 34 : subroutine ice_open_clos
1105 :
1106 :
1107 : integer (kind=int_kind) :: i
1108 :
1109 11 : real (kind=dbl_kind) :: xtime
1110 :
1111 : character (char_len_long) filename
1112 :
1113 : ! ice_data_file = 'open_clos_lindsay.dat'
1114 34 : filename = trim(data_dir)//'/SHEBA/'//trim(ice_data_file)
1115 :
1116 34 : write (nu_diag,*) 'Reading ',filename
1117 :
1118 34 : open (nu_open_clos, file=filename, form='formatted')
1119 :
1120 : ! hourly data
1121 297874 : do i=1,ntime
1122 297874 : read(nu_open_clos,*) xtime, open_data(i), clos_data(i)
1123 : enddo
1124 :
1125 34 : end subroutine ice_open_clos
1126 :
1127 : !=======================================================================
1128 :
1129 0 : subroutine get_wave_spec
1130 :
1131 : use icedrv_arrays_column, only: wave_spectrum, wave_sig_ht, &
1132 : dwavefreq, wavefreq
1133 : use icedrv_domain_size, only: nfreq
1134 :
1135 : ! local variables
1136 : integer (kind=int_kind) :: &
1137 : k
1138 :
1139 : real(kind=dbl_kind), dimension(nfreq) :: &
1140 0 : wave_spectrum_profile ! wave spectrum
1141 :
1142 0 : wave_spectrum(:,:) = c0
1143 :
1144 : ! wave spectrum and frequencies
1145 : ! get hardwired frequency bin info and a dummy wave spectrum profile
1146 : call icepack_init_wave(nfreq=nfreq, &
1147 : wave_spectrum_profile=wave_spectrum_profile, &
1148 0 : wavefreq=wavefreq, dwavefreq=dwavefreq)
1149 :
1150 0 : do k = 1, nfreq
1151 0 : wave_spectrum(:,k) = wave_spectrum_profile(k)
1152 : enddo
1153 :
1154 0 : end subroutine get_wave_spec
1155 :
1156 :
1157 : !=======================================================================
1158 :
1159 : end module icedrv_forcing
1160 :
1161 : !=======================================================================
|
