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