Line data Source code
1 : #ifdef ncdf
2 : #define USE_NETCDF
3 : #endif
4 : !=======================================================================
5 : !
6 : ! Reads and interpolates forcing data for atmosphere and ocean quantities.
7 : !
8 : ! authors: Elizabeth C. Hunke and William H. Lipscomb, LANL
9 : !
10 : ! 2004 WHL: Block structure added
11 : ! 2005 WHL: ECMWF option added
12 : ! 2006 ECH: LY option added
13 : ! 2006 WHL: Module name changed from ice_flux_in
14 : ! 2006 ECH: Fixed bugs, rearranged routines, edited comments, etc.
15 : ! Added NCAR ocean forcing file
16 : ! Converted to free source form (F90)
17 : ! 2007: netcdf version of read_data added by Alison McLaren, Met Office
18 : !
19 : module ice_forcing
20 :
21 : use ice_kinds_mod
22 : use ice_boundary, only: ice_HaloUpdate
23 : use ice_blocks, only: nx_block, ny_block
24 : use ice_domain, only: halo_info
25 : use ice_domain_size, only: ncat, max_blocks, nx_global, ny_global, nfreq
26 : use ice_communicate, only: my_task, master_task
27 : use ice_calendar, only: istep, istep1, &
28 : msec, mday, mmonth, myear, yday, daycal, & ! LCOV_EXCL_LINE
29 : daymo, days_per_year, compute_days_between
30 : use ice_fileunits, only: nu_diag, nu_forcing
31 : use ice_exit, only: abort_ice
32 : use ice_read_write, only: ice_open, ice_read, &
33 : ice_get_ncvarsize, ice_read_vec_nc, & ! LCOV_EXCL_LINE
34 : ice_open_nc, ice_read_nc, ice_close_nc
35 : use ice_timers, only: ice_timer_start, ice_timer_stop, timer_readwrite, &
36 : timer_bound, timer_forcing
37 : use ice_arrays_column, only: oceanmixed_ice, restore_bgc
38 : use ice_constants, only: c0, c1, c2, c3, c4, c5, c8, c10, c12, c15, c20, &
39 : c180, c360, c365, c1000, c3600
40 : use ice_constants, only: p001, p01, p1, p2, p25, p5, p6
41 : use ice_constants, only: cm_to_m
42 : use ice_constants, only: field_loc_center, field_type_scalar, &
43 : field_type_vector, field_loc_NEcorner
44 : use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted
45 : use icepack_intfc, only: icepack_sea_freezing_temperature
46 : use icepack_intfc, only: icepack_init_wave, icepack_init_parameters
47 : use icepack_intfc, only: icepack_query_tracer_indices, icepack_query_parameters
48 :
49 : implicit none
50 : private
51 : public :: init_forcing_atmo, init_forcing_ocn, alloc_forcing, &
52 : get_forcing_atmo, get_forcing_ocn, get_wave_spec, & ! LCOV_EXCL_LINE
53 : read_clim_data, read_clim_data_nc, & ! LCOV_EXCL_LINE
54 : interpolate_data, interp_coeff_monthly, & ! LCOV_EXCL_LINE
55 : read_data_nc_point, interp_coeff, & ! LCOV_EXCL_LINE
56 : init_snowtable
57 :
58 : integer (kind=int_kind), public :: &
59 : ycycle , & ! number of years in forcing cycle, set by namelist ! LCOV_EXCL_LINE
60 : fyear_init , & ! first year of data in forcing cycle, set by namelist ! LCOV_EXCL_LINE
61 : fyear , & ! current year in forcing cycle, varying during the run ! LCOV_EXCL_LINE
62 : fyear_final ! last year in cycle, computed at init
63 :
64 : character (char_len_long) :: & ! input data file names
65 : uwind_file, & ! this is also used a generic file containing all fields for JRA55 ! LCOV_EXCL_LINE
66 : vwind_file, & ! LCOV_EXCL_LINE
67 : wind_file, & ! LCOV_EXCL_LINE
68 : strax_file, & ! LCOV_EXCL_LINE
69 : stray_file, & ! LCOV_EXCL_LINE
70 : tair_file, & ! LCOV_EXCL_LINE
71 : humid_file, & ! LCOV_EXCL_LINE
72 : rhoa_file, & ! LCOV_EXCL_LINE
73 : fsw_file, & ! LCOV_EXCL_LINE
74 : flw_file, & ! LCOV_EXCL_LINE
75 : rain_file, & ! LCOV_EXCL_LINE
76 : sst_file, & ! LCOV_EXCL_LINE
77 : sss_file, & ! LCOV_EXCL_LINE
78 : sublim_file, & ! LCOV_EXCL_LINE
79 : snow_file
80 :
81 : character (char_len_long), dimension(:), allocatable, public :: & ! input data file names
82 : topmelt_file, & ! LCOV_EXCL_LINE
83 : botmelt_file
84 :
85 : real (kind=dbl_kind), public :: &
86 : c1intp, c2intp ! interpolation coefficients
87 :
88 : integer (kind=int_kind) :: &
89 : oldrecnum = 0 , & ! old record number (save between steps) ! LCOV_EXCL_LINE
90 : oldrecnum4X = 0 !
91 :
92 : real (kind=dbl_kind), dimension(:,:,:), allocatable, public :: &
93 : cldf ! cloud fraction
94 :
95 : real (kind=dbl_kind), dimension(:,:,:,:), allocatable, public :: &
96 : fsw_data, & ! field values at 2 temporal data points ! LCOV_EXCL_LINE
97 : cldf_data, & ! LCOV_EXCL_LINE
98 : fsnow_data, & ! LCOV_EXCL_LINE
99 : Tair_data, & ! LCOV_EXCL_LINE
100 : uatm_data, & ! LCOV_EXCL_LINE
101 : vatm_data, & ! LCOV_EXCL_LINE
102 : wind_data, & ! LCOV_EXCL_LINE
103 : strax_data, & ! LCOV_EXCL_LINE
104 : stray_data, & ! LCOV_EXCL_LINE
105 : Qa_data, & ! LCOV_EXCL_LINE
106 : rhoa_data, & ! LCOV_EXCL_LINE
107 : flw_data, & ! LCOV_EXCL_LINE
108 : sst_data, & ! LCOV_EXCL_LINE
109 : sss_data, & ! LCOV_EXCL_LINE
110 : uocn_data, & ! LCOV_EXCL_LINE
111 : vocn_data, & ! LCOV_EXCL_LINE
112 : sublim_data, & ! LCOV_EXCL_LINE
113 : frain_data
114 :
115 : real (kind=dbl_kind), dimension(:,:,:,:,:), allocatable, public :: &
116 : topmelt_data, & ! LCOV_EXCL_LINE
117 : botmelt_data
118 :
119 : real (kind=dbl_kind), dimension(:,:,:,:,:), allocatable :: &
120 : wave_spectrum_data ! field values at 2 temporal data points
121 :
122 : character(char_len), public :: &
123 : atm_data_format, & ! 'bin'=binary or 'nc'=netcdf ! LCOV_EXCL_LINE
124 : ocn_data_format, & ! 'bin'=binary or 'nc'=netcdf ! LCOV_EXCL_LINE
125 : atm_data_type, & ! 'default', 'monthly', 'ncar', 'box2001' ! LCOV_EXCL_LINE
126 : ! 'hadgem', 'oned', 'calm', 'uniform'
127 : ! 'JRA55' or 'JRA55do'
128 : bgc_data_type, & ! 'default', 'clim'
129 : ocn_data_type, & ! 'default', 'clim', 'ncar', 'oned', 'calm', 'box2001' ! LCOV_EXCL_LINE
130 : ! 'hadgem_sst' or 'hadgem_sst_uvocn', 'uniform'
131 : ice_data_type, & ! 'latsst', 'box2001', 'boxslotcyl', etc
132 : ice_data_conc, & ! 'p5','p8','p9','c1','parabolic', 'box2001', etc ! LCOV_EXCL_LINE
133 : ice_data_dist, & ! 'box2001','gauss', 'uniform', etc ! LCOV_EXCL_LINE
134 : precip_units ! 'mm_per_month', 'mm_per_sec', 'mks','m_per_sec'
135 :
136 : logical (kind=log_kind), public :: &
137 : rotate_wind ! rotate wind/stress to computational grid from true north directed
138 :
139 : character(char_len_long), public :: &
140 : atm_data_dir , & ! top directory for atmospheric data ! LCOV_EXCL_LINE
141 : ocn_data_dir , & ! top directory for ocean data ! LCOV_EXCL_LINE
142 : wave_spec_dir, & ! dir name for wave spectrum ! LCOV_EXCL_LINE
143 : wave_spec_file,& ! file name for wave spectrum ! LCOV_EXCL_LINE
144 : oceanmixed_file ! file name for ocean forcing data
145 :
146 : integer (kind=int_kind), parameter :: &
147 : nfld = 8 ! number of fields to search for in forcing file
148 :
149 : ! as in the dummy atm (latm)
150 : real (kind=dbl_kind), parameter, public :: &
151 : frcvdr = 0.28_dbl_kind, & ! frac of incoming sw in vis direct band ! LCOV_EXCL_LINE
152 : frcvdf = 0.24_dbl_kind, & ! frac of incoming sw in vis diffuse band ! LCOV_EXCL_LINE
153 : frcidr = 0.31_dbl_kind, & ! frac of incoming sw in near IR direct band ! LCOV_EXCL_LINE
154 : frcidf = 0.17_dbl_kind ! frac of incoming sw in near IR diffuse band
155 :
156 : real (kind=dbl_kind), dimension (:,:,:,:,:), allocatable, public :: &
157 : ocn_frc_m ! ocn data for 12 months
158 :
159 : logical (kind=log_kind), public :: &
160 : restore_ocn ! restore sst if true
161 :
162 : integer (kind=int_kind), public :: &
163 : trestore ! restoring time scale (days)
164 :
165 : real (kind=dbl_kind), public :: &
166 : trest ! restoring time scale (sec)
167 :
168 : logical (kind=log_kind), public :: &
169 : debug_forcing ! prints forcing debugging output if true
170 :
171 : real (dbl_kind), dimension(:), allocatable, public :: &
172 : jday_atm ! jday time vector from atm forcing files
173 :
174 : integer (kind=int_kind), public :: &
175 : Njday_atm ! Number of atm forcing timesteps
176 :
177 : character (len=char_len_long), public :: &
178 : snw_filename ! filename for snow lookup table
179 :
180 : character (char_len), public :: &
181 : snw_rhos_fname , & ! snow table 1d rhos field name ! LCOV_EXCL_LINE
182 : snw_Tgrd_fname , & ! snow table 1d Tgrd field name ! LCOV_EXCL_LINE
183 : snw_T_fname , & ! snow table 1d T field name ! LCOV_EXCL_LINE
184 : snw_tau_fname , & ! snow table 3d tau field name ! LCOV_EXCL_LINE
185 : snw_kappa_fname, & ! snow table 3d kappa field name ! LCOV_EXCL_LINE
186 : snw_drdt0_fname ! snow table 3d drdt0 field name
187 :
188 : ! PRIVATE:
189 :
190 : real (dbl_kind), parameter :: &
191 : mixed_layer_depth_default = c20 ! default mixed layer depth in m
192 :
193 : logical (kind=log_kind), parameter :: &
194 : local_debug = .false. ! local debug flag
195 :
196 : !=======================================================================
197 :
198 : contains
199 :
200 : !=======================================================================
201 : !
202 : ! Allocate space for all variables
203 : !
204 37 : subroutine alloc_forcing
205 : integer (int_kind) :: ierr
206 : character(len=*), parameter :: subname = '(alloc_forcing)'
207 :
208 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
209 :
210 : allocate ( &
211 : cldf(nx_block,ny_block, max_blocks), & ! cloud fraction ! LCOV_EXCL_LINE
212 : fsw_data(nx_block,ny_block,2,max_blocks), & ! field values at 2 temporal data points ! LCOV_EXCL_LINE
213 : cldf_data(nx_block,ny_block,2,max_blocks), & ! LCOV_EXCL_LINE
214 : fsnow_data(nx_block,ny_block,2,max_blocks), & ! LCOV_EXCL_LINE
215 : Tair_data(nx_block,ny_block,2,max_blocks), & ! LCOV_EXCL_LINE
216 : uatm_data(nx_block,ny_block,2,max_blocks), & ! LCOV_EXCL_LINE
217 : vatm_data(nx_block,ny_block,2,max_blocks), & ! LCOV_EXCL_LINE
218 : wind_data(nx_block,ny_block,2,max_blocks), & ! LCOV_EXCL_LINE
219 : strax_data(nx_block,ny_block,2,max_blocks), & ! LCOV_EXCL_LINE
220 : stray_data(nx_block,ny_block,2,max_blocks), & ! LCOV_EXCL_LINE
221 : Qa_data(nx_block,ny_block,2,max_blocks), & ! LCOV_EXCL_LINE
222 : rhoa_data(nx_block,ny_block,2,max_blocks), & ! LCOV_EXCL_LINE
223 : flw_data(nx_block,ny_block,2,max_blocks), & ! LCOV_EXCL_LINE
224 : sst_data(nx_block,ny_block,2,max_blocks), & ! LCOV_EXCL_LINE
225 : sss_data(nx_block,ny_block,2,max_blocks), & ! LCOV_EXCL_LINE
226 : uocn_data(nx_block,ny_block,2,max_blocks), & ! LCOV_EXCL_LINE
227 : vocn_data(nx_block,ny_block,2,max_blocks), & ! LCOV_EXCL_LINE
228 : sublim_data(nx_block,ny_block,2,max_blocks), & ! LCOV_EXCL_LINE
229 : frain_data(nx_block,ny_block,2,max_blocks), & ! LCOV_EXCL_LINE
230 : topmelt_data(nx_block,ny_block,2,max_blocks,ncat), & ! LCOV_EXCL_LINE
231 : botmelt_data(nx_block,ny_block,2,max_blocks,ncat), & ! LCOV_EXCL_LINE
232 : ocn_frc_m(nx_block,ny_block, max_blocks,nfld,12), & ! ocn data for 12 months ! LCOV_EXCL_LINE
233 : topmelt_file(ncat), & ! LCOV_EXCL_LINE
234 : botmelt_file(ncat), & ! LCOV_EXCL_LINE
235 : wave_spectrum_data(nx_block,ny_block,nfreq,2,max_blocks), & ! LCOV_EXCL_LINE
236 37 : stat=ierr)
237 37 : if (ierr/=0) call abort_ice('(alloc_forcing): Out of Memory')
238 :
239 : ! initialize this, not set in box2001 (and some other forcings?)
240 :
241 111936 : cldf = c0
242 :
243 37 : end subroutine alloc_forcing
244 :
245 : !=======================================================================
246 :
247 37 : subroutine init_forcing_atmo
248 :
249 : ! Determine the current and final year of the forcing cycle based on
250 : ! namelist input; initialize the atmospheric forcing data filenames.
251 :
252 : use ice_calendar, only: use_leap_years
253 :
254 : integer (kind=int_kind) :: modadj ! adjustment for mod function
255 : character(len=*), parameter :: subname = '(init_forcing_atmo)'
256 :
257 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
258 :
259 37 : modadj = abs((min(0,myear-fyear_init)/ycycle+1)*ycycle)
260 37 : fyear = fyear_init + mod(myear-fyear_init+modadj,ycycle)
261 37 : fyear_final = fyear_init + ycycle - 1 ! last year in forcing cycle
262 :
263 : if (local_debug .and. my_task == master_task) then
264 : write(nu_diag,*) subname,'fdbg fyear = ',fyear,fyear_init,fyear_final
265 : write(nu_diag,*) subname,'fdbg atm_data_type = ',trim(atm_data_type)
266 : endif
267 :
268 37 : if (trim(atm_data_type) /= 'default' .and. &
269 : my_task == master_task) then
270 7 : write (nu_diag,*) ' Initial forcing data year = ',fyear_init
271 7 : write (nu_diag,*) ' Final forcing data year = ',fyear_final
272 : endif
273 :
274 37 : if (trim(atm_data_type) == 'hadgem' .and. &
275 : trim(precip_units) /= 'mks') then
276 0 : if (my_task == master_task) then
277 0 : write (nu_diag,*) 'WARNING: HadGEM atmospheric data chosen with wrong precip_units'
278 0 : write (nu_diag,*) 'WARNING: Changing precip_units to mks (i.e. kg/m2 s).'
279 : endif
280 : call abort_ice(error_message=subname//' HadGEM precip_units error', &
281 0 : file=__FILE__, line=__LINE__)
282 : endif
283 :
284 37 : if (use_leap_years .and. (index(trim(atm_data_type),'JRA55') == 0 .and. &
285 : trim(atm_data_type) /= 'hycom' .and. & ! LCOV_EXCL_LINE
286 : trim(atm_data_type) /= 'box2001')) then
287 0 : write(nu_diag,*) 'use_leap_years option is currently only supported for'
288 0 : write(nu_diag,*) 'JRA55, JRA55do, default , and box2001 atmospheric data'
289 0 : call abort_ice(error_message=subname, file=__FILE__, line=__LINE__)
290 : endif
291 :
292 : !-------------------------------------------------------------------
293 : ! Get filenames for input forcing data
294 : !-------------------------------------------------------------------
295 :
296 : ! default forcing values from init_flux_atm
297 37 : if (trim(atm_data_type) == 'ncar') then
298 0 : call NCAR_files(fyear)
299 37 : elseif (index(trim(atm_data_type),'JRA55') > 0) then
300 21 : call JRA55_files(fyear)
301 16 : elseif (trim(atm_data_type) == 'hadgem') then
302 0 : call hadgem_files(fyear)
303 16 : elseif (trim(atm_data_type) == 'monthly') then
304 0 : call monthly_files(fyear)
305 16 : elseif (trim(atm_data_type) == 'oned') then
306 0 : call oned_files
307 16 : elseif (trim(atm_data_type) == 'ISPOL') then
308 0 : call ISPOL_files
309 16 : elseif (trim(atm_data_type) == 'box2001') then
310 16 : call box2001_data_atm
311 0 : elseif (trim(atm_data_type) == 'uniform_northeast') then
312 0 : call uniform_data_atm('NE')
313 0 : elseif (trim(atm_data_type) == 'uniform_north') then
314 0 : call uniform_data_atm('N')
315 0 : elseif (trim(atm_data_type) == 'uniform_east') then
316 0 : call uniform_data_atm('E')
317 0 : elseif (trim(atm_data_type) == 'uniform_south') then
318 0 : call uniform_data_atm('S')
319 0 : elseif (trim(atm_data_type) == 'uniform_west') then
320 0 : call uniform_data_atm('W')
321 0 : elseif (trim(atm_data_type) == 'calm') then
322 0 : call uniform_data_atm('N',c0) ! direction does not matter when c0
323 0 : elseif (trim(atm_data_type) == 'hycom') then
324 0 : call hycom_atm_files
325 0 : elseif (trim(atm_data_type) == 'default') then
326 : ! don't need to do anything more
327 : else
328 : call abort_ice (error_message=subname//' ERROR atm_data_type unknown = '// &
329 0 : trim(atm_data_type), file=__FILE__, line=__LINE__)
330 : endif
331 :
332 37 : end subroutine init_forcing_atmo
333 :
334 : !=======================================================================
335 :
336 37 : subroutine init_forcing_ocn(dt)
337 :
338 : ! Set sea surface salinity and freezing temperature to annual mean value
339 : ! using a 12-month climatology.
340 : ! Read sst data for current month, and adjust sst based on freezing
341 : ! temperature. No interpolation in time.
342 :
343 : ! Note: SST is subsequently prognosed if CICE is run
344 : ! with a mixed layer ocean (oceanmixed_ice = T), and can be
345 : ! restored to data (restore_ocn = T).
346 :
347 : use ice_blocks, only: nx_block, ny_block
348 : use ice_domain, only: nblocks
349 : use ice_domain_size, only: max_blocks
350 : use ice_flux, only: sss, sst, Tf
351 :
352 : real (kind=dbl_kind), intent(in) :: &
353 : dt ! time step
354 :
355 : ! local variables
356 :
357 : integer (kind=int_kind) :: &
358 : i, j, iblk , & ! horizontal indices ! LCOV_EXCL_LINE
359 : k , & ! month index ! LCOV_EXCL_LINE
360 : fid , & ! file id for netCDF file ! LCOV_EXCL_LINE
361 : nbits
362 :
363 : logical (kind=log_kind) :: diag
364 :
365 8 : real (kind=dbl_kind) :: secday
366 :
367 : character (char_len) :: &
368 : fieldname ! field name in netcdf file
369 :
370 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
371 27874 : work1
372 :
373 : character(len=*), parameter :: subname = '(init_forcing_ocn)'
374 :
375 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
376 :
377 37 : call icepack_query_parameters(secday_out=secday)
378 37 : call icepack_warnings_flush(nu_diag)
379 37 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
380 0 : file=__FILE__, line=__LINE__)
381 :
382 37 : call alloc_forcing()
383 :
384 223996 : sst_data(:,:,:,:) = c0
385 223996 : sss_data(:,:,:,:) = c0
386 223996 : uocn_data(:,:,:,:) = c0
387 223996 : vocn_data(:,:,:,:) = c0
388 :
389 37 : nbits = 64 ! double precision data
390 :
391 37 : if (restore_ocn .or. restore_bgc) then
392 0 : if (trestore == 0) then
393 0 : trest = dt ! use data instantaneously
394 : else
395 0 : trest = real(trestore,kind=dbl_kind) * secday ! seconds
396 : endif
397 : endif
398 :
399 : !-------------------------------------------------------------------
400 : ! Sea surface salinity (SSS)
401 : ! initialize to annual climatology created from monthly data
402 : !-------------------------------------------------------------------
403 :
404 37 : if (trim(ocn_data_type) == 'clim') then
405 :
406 0 : sss_file = trim(ocn_data_dir)//'/sss.mm.100x116.da' ! gx3 only
407 :
408 0 : if (my_task == master_task) then
409 0 : write (nu_diag,*) ' '
410 0 : write (nu_diag,*) 'SSS climatology computed from:'
411 0 : write (nu_diag,*) trim(sss_file)
412 : endif
413 :
414 0 : if (my_task == master_task) &
415 0 : call ice_open (nu_forcing, sss_file, nbits)
416 :
417 0 : sss(:,:,:) = c0
418 :
419 0 : do k = 1,12 ! loop over 12 months
420 : call ice_read (nu_forcing, k, work1, 'rda8', debug_forcing, &
421 0 : field_loc_center, field_type_scalar)
422 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j)
423 0 : do iblk = 1, nblocks
424 0 : do j = 1, ny_block
425 0 : do i = 1, nx_block
426 0 : sss(i,j,iblk) = sss(i,j,iblk) + work1(i,j,iblk)
427 : enddo
428 : enddo
429 : enddo
430 : !$OMP END PARALLEL DO
431 : enddo ! k
432 :
433 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j)
434 0 : do iblk = 1, nblocks
435 0 : do j = 1, ny_block
436 0 : do i = 1, nx_block
437 0 : sss(i,j,iblk) = sss(i,j,iblk) / c12 ! annual average
438 0 : sss(i,j,iblk) = max(sss(i,j,iblk),c0)
439 : enddo
440 : enddo
441 : enddo
442 : !$OMP END PARALLEL DO
443 :
444 0 : call ocn_freezing_temperature
445 :
446 0 : if (my_task == master_task) close(nu_forcing)
447 :
448 : !-------------------------------------------------------------------
449 : ! Sea surface temperature (SST)
450 : ! initialize to data for current month
451 : !-------------------------------------------------------------------
452 :
453 0 : if (nx_global == 320) then ! gx1
454 0 : sst_file = trim(ocn_data_dir)//'/sst_clim_hurrell.dat'
455 : else ! gx3
456 0 : sst_file = trim(ocn_data_dir)//'/sst.mm.100x116.da'
457 : endif
458 :
459 0 : if (my_task == master_task) then
460 0 : write (nu_diag,*) ' '
461 0 : write (nu_diag,*) 'Initial SST file:', trim(sst_file)
462 : endif
463 :
464 0 : if (my_task == master_task) &
465 0 : call ice_open (nu_forcing, sst_file, nbits)
466 :
467 : call ice_read (nu_forcing, mmonth, sst, 'rda8', debug_forcing, &
468 0 : field_loc_center, field_type_scalar)
469 :
470 0 : if (my_task == master_task) close(nu_forcing)
471 :
472 : ! Make sure sst is not less than freezing temperature Tf
473 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j)
474 0 : do iblk = 1, nblocks
475 0 : do j = 1, ny_block
476 0 : do i = 1, nx_block
477 0 : sst(i,j,iblk) = max(sst(i,j,iblk),Tf(i,j,iblk))
478 : enddo
479 : enddo
480 : enddo
481 : !$OMP END PARALLEL DO
482 :
483 37 : elseif (trim(ocn_data_type) == 'hadgem_sst' .or. &
484 : trim(ocn_data_type) == 'hadgem_sst_uvocn') then
485 :
486 0 : diag = .true. ! write diagnostic information
487 :
488 0 : sst_file = trim (ocn_data_dir)//'/MONTHLY/sst.1997.nc'
489 :
490 0 : if (my_task == master_task) then
491 :
492 0 : write (nu_diag,*) ' '
493 0 : write (nu_diag,*) 'Initial SST file:', trim(sst_file)
494 :
495 0 : call ice_open_nc(sst_file,fid)
496 :
497 : endif
498 :
499 0 : fieldname='sst'
500 0 : call ice_read_nc(fid,mmonth,fieldname,sst,diag)
501 :
502 0 : if (my_task == master_task) call ice_close_nc(fid)
503 :
504 : ! Make sure sst is not less than freezing temperature Tf
505 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j)
506 0 : do iblk = 1, nblocks
507 0 : do j = 1, ny_block
508 0 : do i = 1, nx_block
509 0 : sst(i,j,iblk) = max(sst(i,j,iblk),Tf(i,j,iblk))
510 : enddo
511 : enddo
512 : enddo
513 : !$OMP END PARALLEL DO
514 :
515 37 : elseif (trim(ocn_data_type) == 'ncar') then
516 0 : call ocn_data_ncar_init
517 : ! call ocn_data_ncar_init_3D
518 :
519 37 : elseif (trim(ocn_data_type) == 'hycom') then
520 0 : call ocn_data_hycom_init
521 :
522 37 : elseif (trim(ocn_data_type) == 'box2001') then
523 16 : call box2001_data_ocn
524 :
525 : ! uniform forcing options
526 21 : elseif (trim(ocn_data_type) == 'uniform_northeast') then
527 0 : call uniform_data_ocn('NE',p1)
528 21 : elseif (trim(ocn_data_type) == 'uniform_east') then
529 0 : call uniform_data_ocn('E',p1)
530 21 : elseif (trim(ocn_data_type) == 'uniform_north') then
531 0 : call uniform_data_ocn('N',p1)
532 21 : elseif (trim(ocn_data_type) == 'calm') then
533 0 : call uniform_data_ocn('N',c0) ! directon does not matter for c0
534 21 : elseif (trim(ocn_data_type) == 'default') then
535 : ! don't need to do anything more
536 : else
537 : call abort_ice (error_message=subname//' ERROR ocn_data_type unknown = '// &
538 0 : trim(ocn_data_type), file=__FILE__, line=__LINE__)
539 : endif
540 :
541 37 : end subroutine init_forcing_ocn
542 :
543 : !=======================================================================
544 :
545 0 : subroutine ocn_freezing_temperature
546 :
547 : ! Compute ocean freezing temperature Tf based on tfrz_option
548 : ! 'minus1p8' Tf = -1.8 C (default)
549 : ! 'linear_salt' Tf = -depressT * sss
550 : ! 'mushy' Tf conforms with mushy layer thermo (ktherm=2)
551 :
552 : use ice_blocks, only: nx_block, ny_block
553 : use ice_domain, only: nblocks
554 : use ice_flux, only: sss, Tf
555 :
556 : ! local variables
557 :
558 : integer (kind=int_kind) :: &
559 : i, j, iblk ! horizontal indices
560 :
561 : character(len=*), parameter :: subname = '(ocn_freezing_temperature)'
562 :
563 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
564 :
565 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j)
566 0 : do iblk = 1, nblocks
567 0 : do j = 1, ny_block
568 0 : do i = 1, nx_block
569 0 : Tf(i,j,iblk) = icepack_sea_freezing_temperature(sss(i,j,iblk))
570 : enddo
571 : enddo
572 : enddo
573 : !$OMP END PARALLEL DO
574 :
575 0 : call icepack_warnings_flush(nu_diag)
576 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
577 0 : file=__FILE__, line=__LINE__)
578 :
579 0 : end subroutine ocn_freezing_temperature
580 :
581 : !=======================================================================
582 :
583 5784 : subroutine get_forcing_atmo
584 :
585 : ! Get atmospheric forcing data and interpolate as necessary
586 :
587 : use ice_blocks, only: block, get_block
588 : use ice_domain, only: nblocks, blocks_ice
589 : use ice_flux, only: Tair, fsw, flw, frain, fsnow, Qa, rhoa, &
590 : uatm, vatm, strax, stray, zlvl, wind, swvdr, swvdf, swidr, swidf, & ! LCOV_EXCL_LINE
591 : potT, sst
592 : use ice_state, only: aice, trcr
593 : use ice_grid, only: ANGLET, hm
594 :
595 : integer (kind=int_kind) :: &
596 : iblk, & ! block index ! LCOV_EXCL_LINE
597 : ilo,ihi,jlo,jhi, & ! beginning and end of physical domain ! LCOV_EXCL_LINE
598 : modadj, & ! adjustment to make mod a postive number ! LCOV_EXCL_LINE
599 : fyear_old, & ! fyear setting on last timestep ! LCOV_EXCL_LINE
600 : nt_Tsfc
601 :
602 : type (block) :: &
603 : this_block ! block information for current block
604 :
605 : character(len=*), parameter :: subname = '(get_forcing_atmo)'
606 :
607 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
608 :
609 5784 : call ice_timer_start(timer_forcing)
610 :
611 5784 : fyear_old = fyear
612 5784 : modadj = abs((min(0,myear-fyear_init)/ycycle+1)*ycycle)
613 5784 : fyear = fyear_init + mod(myear-fyear_init+modadj,ycycle)
614 5784 : if (trim(atm_data_type) /= 'default' .and. &
615 : (istep <= 1 .or. fyear /= fyear_old)) then
616 37 : if (my_task == master_task) then
617 7 : write (nu_diag,*) ' Set current forcing data year = ',fyear
618 : endif
619 : endif
620 :
621 5784 : call icepack_query_tracer_indices(nt_Tsfc_out=nt_Tsfc)
622 5784 : call icepack_warnings_flush(nu_diag)
623 5784 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
624 0 : file=__FILE__, line=__LINE__)
625 :
626 : !-------------------------------------------------------------------
627 : ! Read and interpolate atmospheric data
628 : !-------------------------------------------------------------------
629 :
630 : if (local_debug .and. my_task == master_task) then
631 : write(nu_diag,*) subname,'fdbg fyear = ',fyear
632 : write(nu_diag,*) subname,'fdbg atm_data_type = ',trim(atm_data_type)
633 : endif
634 :
635 5784 : if (trim(atm_data_type) == 'ncar') then
636 0 : call ncar_data
637 5784 : elseif (index(trim(atm_data_type),'JRA55') > 0) then
638 2904 : call JRA55_data
639 2880 : elseif (trim(atm_data_type) == 'hadgem') then
640 0 : call hadgem_data
641 2880 : elseif (trim(atm_data_type) == 'monthly') then
642 0 : call monthly_data
643 2880 : elseif (trim(atm_data_type) == 'oned') then
644 0 : call oned_data
645 2880 : elseif (trim(atm_data_type) == 'box2001') then
646 2880 : call box2001_data_atm
647 0 : elseif (trim(atm_data_type) == 'uniform_northeast') then
648 0 : call uniform_data_atm('NE')
649 0 : elseif (trim(atm_data_type) == 'uniform_north') then
650 0 : call uniform_data_atm('N')
651 0 : elseif (trim(atm_data_type) == 'uniform_east') then
652 0 : call uniform_data_atm('E')
653 0 : elseif (trim(atm_data_type) == 'uniform_south') then
654 0 : call uniform_data_atm('S')
655 0 : elseif (trim(atm_data_type) == 'uniform_west') then
656 0 : call uniform_data_atm('W')
657 0 : elseif (trim(atm_data_type) == 'calm') then
658 0 : call uniform_data_atm('N',c0) ! direction does not matter when c0
659 0 : elseif (trim(atm_data_type) == 'hycom') then
660 0 : call hycom_atm_data
661 : !elseif (trim(atm_data_type) == 'uniform_northeast') then
662 : !elseif (trim(atm_data_type) == 'uniform_east') then
663 : !elseif (trim(atm_data_type) == 'uniform_north') then
664 : else ! default values set in init_flux
665 0 : return
666 : endif
667 :
668 : !-------------------------------------------------------------------
669 : ! Convert forcing data to fields needed by ice model
670 : !-------------------------------------------------------------------
671 :
672 2880 : !$OMP PARALLEL DO PRIVATE(iblk,ilo,ihi,jlo,jhi,this_block)
673 8688 : do iblk = 1, nblocks
674 :
675 5784 : this_block = get_block(blocks_ice(iblk),iblk)
676 5784 : ilo = this_block%ilo
677 5784 : ihi = this_block%ihi
678 5784 : jlo = this_block%jlo
679 5784 : jhi = this_block%jhi
680 :
681 : call prepare_forcing (nx_block, ny_block, &
682 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
683 : hm (:,:,iblk), & ! LCOV_EXCL_LINE
684 : Tair (:,:,iblk), & ! LCOV_EXCL_LINE
685 : fsw (:,:,iblk), & ! LCOV_EXCL_LINE
686 : cldf (:,:,iblk), & ! LCOV_EXCL_LINE
687 : flw (:,:,iblk), & ! LCOV_EXCL_LINE
688 : frain (:,:,iblk), & ! LCOV_EXCL_LINE
689 : fsnow (:,:,iblk), & ! LCOV_EXCL_LINE
690 : Qa (:,:,iblk), & ! LCOV_EXCL_LINE
691 : rhoa (:,:,iblk), & ! LCOV_EXCL_LINE
692 : uatm (:,:,iblk), & ! LCOV_EXCL_LINE
693 : vatm (:,:,iblk), & ! LCOV_EXCL_LINE
694 : strax (:,:,iblk), & ! LCOV_EXCL_LINE
695 : stray (:,:,iblk), & ! LCOV_EXCL_LINE
696 : zlvl (:,:,iblk), & ! LCOV_EXCL_LINE
697 : wind (:,:,iblk), & ! LCOV_EXCL_LINE
698 : swvdr (:,:,iblk), & ! LCOV_EXCL_LINE
699 : swvdf (:,:,iblk), & ! LCOV_EXCL_LINE
700 : swidr (:,:,iblk), & ! LCOV_EXCL_LINE
701 : swidf (:,:,iblk), & ! LCOV_EXCL_LINE
702 : potT (:,:,iblk), & ! LCOV_EXCL_LINE
703 : ANGLET(:,:,iblk), & ! LCOV_EXCL_LINE
704 : trcr (:,:,nt_Tsfc,iblk), & ! LCOV_EXCL_LINE
705 : sst (:,:,iblk), & ! LCOV_EXCL_LINE
706 11568 : aice (:,:,iblk) )
707 :
708 : enddo ! iblk
709 : !$OMP END PARALLEL DO
710 :
711 5784 : call ice_timer_start(timer_bound)
712 : call ice_HaloUpdate (swvdr, halo_info, &
713 5784 : field_loc_center, field_type_scalar)
714 : call ice_HaloUpdate (swvdf, halo_info, &
715 5784 : field_loc_center, field_type_scalar)
716 : call ice_HaloUpdate (swidr, halo_info, &
717 5784 : field_loc_center, field_type_scalar)
718 : call ice_HaloUpdate (swidf, halo_info, &
719 5784 : field_loc_center, field_type_scalar)
720 5784 : call ice_timer_stop(timer_bound)
721 :
722 5784 : call ice_timer_stop(timer_forcing)
723 :
724 5784 : end subroutine get_forcing_atmo
725 :
726 : !=======================================================================
727 :
728 5784 : subroutine get_forcing_ocn (dt)
729 :
730 : ! Read and interpolate annual climatologies of SSS and SST.
731 : ! Restore model SST to data if desired.
732 : ! Interpolate ocean fields to U grid if necessary.
733 :
734 : real (kind=dbl_kind), intent(in) :: &
735 : dt ! time step
736 :
737 : character(len=*), parameter :: subname = '(get_forcing_ocn)'
738 :
739 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
740 :
741 5784 : call ice_timer_start(timer_forcing)
742 :
743 : if (local_debug .and. my_task == master_task) then
744 : write(nu_diag,*) subname,'fdbg fyear = ',fyear
745 : write(nu_diag,*) subname,'fdbg ocn_data_type = ',trim(ocn_data_type)
746 : endif
747 :
748 5784 : if (trim(ocn_data_type) == 'clim') then
749 0 : call ocn_data_clim(dt)
750 5784 : elseif (trim(ocn_data_type) == 'ncar' .or. &
751 : trim(ocn_data_type) == 'ISPOL') then
752 0 : call ocn_data_ncar(dt)
753 5784 : elseif (trim(ocn_data_type) == 'hadgem_sst' .or. &
754 : trim(ocn_data_type) == 'hadgem_sst_uvocn') then
755 0 : call ocn_data_hadgem(dt)
756 5784 : elseif (trim(ocn_data_type) == 'oned') then
757 0 : call ocn_data_oned
758 5784 : elseif (trim(ocn_data_type) == 'hycom') then
759 : ! call ocn_data_hycom(dt)
760 : !MHRI: NOT IMPLEMENTED YET
761 5784 : elseif (trim(ocn_data_type) == 'box2001') then
762 2880 : call box2001_data_ocn
763 : ! uniform forcing options
764 2904 : elseif (trim(ocn_data_type) == 'uniform_northeast') then
765 : ! tcraig, not time varying
766 0 : call uniform_data_ocn('NE',p1)
767 2904 : elseif (trim(ocn_data_type) == 'uniform_east') then
768 0 : call uniform_data_ocn('E',p1)
769 2904 : elseif (trim(ocn_data_type) == 'uniform_north') then
770 0 : call uniform_data_ocn('N',p1)
771 2904 : elseif (trim(ocn_data_type) == 'calm') then
772 0 : call uniform_data_ocn('N',c0) ! directon does not matter for c0
773 : endif
774 :
775 5784 : call ice_timer_stop(timer_forcing)
776 :
777 5784 : end subroutine get_forcing_ocn
778 :
779 : !=======================================================================
780 :
781 0 : subroutine read_data (flag, recd, yr, ixm, ixx, ixp, &
782 : maxrec, data_file, field_data, & ! LCOV_EXCL_LINE
783 : field_loc, field_type)
784 :
785 : ! If data is at the beginning of a one-year record, get data from
786 : ! the previous year.
787 : ! If data is at the end of a one-year record, get data from the
788 : ! following year.
789 : ! If no earlier data exists (beginning of fyear_init), then
790 : ! (1) For monthly data, get data from the end of fyear_final.
791 : ! (2) For more frequent data, let the ixm value equal the
792 : ! first value of the year.
793 : ! If no later data exists (end of fyear_final), then
794 : ! (1) For monthly data, get data from the beginning of fyear_init.
795 : ! (2) For more frequent data, let the ixp value
796 : ! equal the last value of the year.
797 : ! In other words, we assume persistence when daily or 6-hourly
798 : ! data is missing, and we assume periodicity when monthly data
799 : ! is missing.
800 :
801 : use ice_diagnostics, only: debug_model_step
802 :
803 : logical (kind=log_kind), intent(in) :: flag
804 :
805 : integer (kind=int_kind), intent(in) :: &
806 : recd , & ! baseline record number ! LCOV_EXCL_LINE
807 : yr , & ! year of forcing data ! LCOV_EXCL_LINE
808 : ixm, ixx, ixp , & ! record numbers of 3 data values ! LCOV_EXCL_LINE
809 : ! relative to recd
810 : maxrec ! maximum record value
811 :
812 : real (kind=dbl_kind), dimension(nx_block,ny_block,2,max_blocks), intent(inout) :: &
813 : field_data ! 2 values needed for interpolation
814 :
815 : integer (kind=int_kind), intent(in) :: &
816 : field_loc, & ! location of field on staggered grid ! LCOV_EXCL_LINE
817 : field_type ! type of field (scalar, vector, angle)
818 :
819 : ! local variables
820 :
821 : character (char_len_long) :: &
822 : data_file ! data file to be read
823 :
824 : integer (kind=int_kind) :: &
825 : nbits , & ! = 32 for single precision, 64 for double ! LCOV_EXCL_LINE
826 : nrec , & ! record number to read ! LCOV_EXCL_LINE
827 : n2, n4 , & ! like ixm and ixp, but ! LCOV_EXCL_LINE
828 : ! adjusted at beginning and end of data
829 : arg ! value of time argument in field_data
830 :
831 : character(len=*), parameter :: subname = '(read_data)'
832 :
833 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
834 :
835 0 : call ice_timer_start(timer_readwrite) ! reading/writing
836 :
837 0 : nbits = 64 ! double precision data
838 :
839 0 : if (istep1 > debug_model_step) debug_forcing = .true. !! debugging
840 :
841 0 : if (my_task==master_task .and. (debug_forcing)) then
842 0 : write(nu_diag,*) ' ', trim(data_file)
843 : endif
844 :
845 0 : if (flag) then
846 :
847 : !-----------------------------------------------------------------
848 : ! Initialize record counters
849 : ! (n2, n4 will change only at the very beginning or end of
850 : ! a forcing cycle.)
851 : !-----------------------------------------------------------------
852 0 : n2 = ixm
853 0 : n4 = ixp
854 0 : arg = 0
855 :
856 : !-----------------------------------------------------------------
857 : ! read data
858 : !-----------------------------------------------------------------
859 :
860 0 : if (ixm /= -99) then
861 : ! currently in first half of data interval
862 0 : if (ixx <= 1) then
863 0 : if (yr > fyear_init) then ! get data from previous year
864 0 : call file_year (data_file, yr-1)
865 : else ! yr = fyear_init, no prior data exists
866 0 : if (maxrec > 12) then ! extrapolate from first record
867 0 : if (ixx == 1) n2 = ixx
868 : else ! go to end of fyear_final
869 0 : call file_year (data_file, fyear_final)
870 : endif
871 : endif ! yr > fyear_init
872 : endif ! ixx <= 1
873 :
874 0 : call ice_open (nu_forcing, data_file, nbits)
875 :
876 0 : arg = 1
877 0 : nrec = recd + n2
878 0 : call ice_read (nu_forcing, nrec, field_data(:,:,arg,:), &
879 0 : 'rda8', debug_forcing, field_loc, field_type)
880 :
881 0 : if (ixx==1 .and. my_task == master_task) close(nu_forcing)
882 : endif ! ixm ne -99
883 :
884 : ! always read ixx data from data file for current year
885 0 : call file_year (data_file, yr)
886 0 : call ice_open (nu_forcing, data_file, nbits)
887 :
888 0 : arg = arg + 1
889 0 : nrec = recd + ixx
890 0 : call ice_read (nu_forcing, nrec, field_data(:,:,arg,:), &
891 0 : 'rda8', debug_forcing, field_loc, field_type)
892 :
893 0 : if (ixp /= -99) then
894 : ! currently in latter half of data interval
895 0 : if (ixx==maxrec) then
896 0 : if (yr < fyear_final) then ! get data from following year
897 0 : if (my_task == master_task) close(nu_forcing)
898 0 : call file_year (data_file, yr+1)
899 0 : call ice_open (nu_forcing, data_file, nbits)
900 : else ! yr = fyear_final, no more data exists
901 0 : if (maxrec > 12) then ! extrapolate from ixx
902 0 : n4 = ixx
903 : else ! go to beginning of fyear_init
904 0 : if (my_task == master_task) close(nu_forcing)
905 0 : call file_year (data_file, fyear_init)
906 :
907 0 : call ice_open (nu_forcing, data_file, nbits)
908 :
909 : endif
910 : endif ! yr < fyear_final
911 : endif ! ixx = maxrec
912 :
913 0 : arg = arg + 1
914 0 : nrec = recd + n4
915 0 : call ice_read (nu_forcing, nrec, field_data(:,:,arg,:), &
916 0 : 'rda8', debug_forcing, field_loc, field_type)
917 : endif ! ixp /= -99
918 :
919 0 : if (my_task == master_task) close(nu_forcing)
920 :
921 : endif ! flag
922 :
923 0 : call ice_timer_stop(timer_readwrite) ! reading/writing
924 :
925 0 : end subroutine read_data
926 :
927 : !=======================================================================
928 :
929 0 : subroutine read_data_nc (flag, recd, yr, ixm, ixx, ixp, &
930 : maxrec, data_file, fieldname, field_data, & ! LCOV_EXCL_LINE
931 : field_loc, field_type)
932 :
933 : ! If data is at the beginning of a one-year record, get data from
934 : ! the previous year.
935 : ! If data is at the end of a one-year record, get data from the
936 : ! following year.
937 : ! If no earlier data exists (beginning of fyear_init), then
938 : ! (1) For monthly data, get data from the end of fyear_final.
939 : ! (2) For more frequent data, let the ixm value equal the
940 : ! first value of the year.
941 : ! If no later data exists (end of fyear_final), then
942 : ! (1) For monthly data, get data from the beginning of fyear_init.
943 : ! (2) For more frequent data, let the ixp value
944 : ! equal the last value of the year.
945 : ! In other words, we assume persistence when daily or 6-hourly
946 : ! data is missing, and we assume periodicity when monthly data
947 : ! is missing.
948 : !
949 : ! Adapted by Alison McLaren, Met Office from read_data
950 :
951 : use ice_diagnostics, only: debug_model_step
952 :
953 : logical (kind=log_kind), intent(in) :: flag
954 :
955 : integer (kind=int_kind), intent(in) :: &
956 : recd , & ! baseline record number ! LCOV_EXCL_LINE
957 : yr , & ! year of forcing data ! LCOV_EXCL_LINE
958 : ixm, ixx, ixp , & ! record numbers of 3 data values ! LCOV_EXCL_LINE
959 : ! relative to recd
960 : maxrec ! maximum record value
961 :
962 : character (char_len_long) :: &
963 : data_file ! data file to be read
964 :
965 : character (char_len), intent(in) :: &
966 : fieldname ! field name in netCDF file
967 :
968 : integer (kind=int_kind), intent(in) :: &
969 : field_loc, & ! location of field on staggered grid ! LCOV_EXCL_LINE
970 : field_type ! type of field (scalar, vector, angle)
971 :
972 : real (kind=dbl_kind), dimension(nx_block,ny_block,2,max_blocks), intent(out) :: &
973 : field_data ! 2 values needed for interpolation
974 :
975 : ! local variables
976 :
977 : integer (kind=int_kind) :: &
978 : nrec , & ! record number to read ! LCOV_EXCL_LINE
979 : n2, n4 , & ! like ixm and ixp, but ! LCOV_EXCL_LINE
980 : ! adjusted at beginning and end of data
981 : arg , & ! value of time argument in field_data
982 : fid ! file id for netCDF routines
983 :
984 : character(len=*), parameter :: subname = '(read_data_nc)'
985 :
986 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
987 :
988 0 : call ice_timer_start(timer_readwrite) ! reading/writing
989 :
990 0 : if (istep1 > debug_model_step) debug_forcing = .true. !! debugging
991 :
992 0 : if (my_task==master_task .and. (debug_forcing)) then
993 0 : write(nu_diag,*) ' ', trim(data_file)
994 : endif
995 :
996 0 : if (flag) then
997 :
998 : !-----------------------------------------------------------------
999 : ! Initialize record counters
1000 : ! (n2, n4 will change only at the very beginning or end of
1001 : ! a forcing cycle.)
1002 : !-----------------------------------------------------------------
1003 0 : n2 = ixm
1004 0 : n4 = ixp
1005 0 : arg = 0
1006 :
1007 : !-----------------------------------------------------------------
1008 : ! read data
1009 : !-----------------------------------------------------------------
1010 :
1011 0 : if (ixm /= -99) then
1012 : ! currently in first half of data interval
1013 0 : if (ixx <= 1) then
1014 0 : if (yr > fyear_init) then ! get data from previous year
1015 0 : call file_year (data_file, yr-1)
1016 : else ! yr = fyear_init, no prior data exists
1017 0 : if (maxrec > 12) then ! extrapolate from first record
1018 0 : if (ixx == 1) n2 = ixx
1019 : else ! go to end of fyear_final
1020 0 : call file_year (data_file, fyear_final)
1021 : endif
1022 : endif ! yr > fyear_init
1023 : endif ! ixx <= 1
1024 :
1025 0 : call ice_open_nc (data_file, fid)
1026 :
1027 0 : arg = 1
1028 0 : nrec = recd + n2
1029 :
1030 : call ice_read_nc &
1031 : (fid, nrec, fieldname, field_data(:,:,arg,:), debug_forcing, & ! LCOV_EXCL_LINE
1032 0 : field_loc, field_type)
1033 :
1034 0 : if (ixx==1) call ice_close_nc(fid)
1035 : endif ! ixm ne -99
1036 :
1037 : ! always read ixx data from data file for current year
1038 0 : call file_year (data_file, yr)
1039 0 : call ice_open_nc (data_file, fid)
1040 :
1041 0 : arg = arg + 1
1042 0 : nrec = recd + ixx
1043 :
1044 : call ice_read_nc &
1045 : (fid, nrec, fieldname, field_data(:,:,arg,:), debug_forcing, & ! LCOV_EXCL_LINE
1046 0 : field_loc, field_type)
1047 :
1048 0 : if (ixp /= -99) then
1049 : ! currently in latter half of data interval
1050 0 : if (ixx==maxrec) then
1051 0 : if (yr < fyear_final) then ! get data from following year
1052 0 : call ice_close_nc(fid)
1053 0 : call file_year (data_file, yr+1)
1054 0 : call ice_open_nc (data_file, fid)
1055 : else ! yr = fyear_final, no more data exists
1056 0 : if (maxrec > 12) then ! extrapolate from ixx
1057 0 : n4 = ixx
1058 : else ! go to beginning of fyear_init
1059 0 : call ice_close_nc(fid)
1060 0 : call file_year (data_file, fyear_init)
1061 0 : call ice_open_nc (data_file, fid)
1062 :
1063 : endif
1064 : endif ! yr < fyear_final
1065 : endif ! ixx = maxrec
1066 :
1067 0 : arg = arg + 1
1068 0 : nrec = recd + n4
1069 :
1070 : call ice_read_nc &
1071 : (fid, nrec, fieldname, field_data(:,:,arg,:), debug_forcing, & ! LCOV_EXCL_LINE
1072 0 : field_loc, field_type)
1073 : endif ! ixp /= -99
1074 :
1075 0 : call ice_close_nc(fid)
1076 :
1077 : endif ! flag
1078 :
1079 0 : call ice_timer_stop(timer_readwrite) ! reading/writing
1080 :
1081 0 : end subroutine read_data_nc
1082 :
1083 : !=======================================================================
1084 :
1085 0 : subroutine read_data_nc_hycom (flag, recd, &
1086 : data_file, fieldname, field_data, & ! LCOV_EXCL_LINE
1087 : field_loc, field_type)
1088 :
1089 : ! Data is assumed to cover the entire time period of simulation.
1090 : ! It is not bounded by start of year nor end of year
1091 : ! Data must be accesible both before and after (or on) the point in time
1092 : ! Assume increasing timeaxis within the forcing files, but they do not
1093 : ! have to be equal spaced. Read time vector from "MT" in "init_hycom"
1094 : !
1095 : ! Adapted by Mads Hvid Ribergaard, DMI from read_data_nc
1096 :
1097 : use ice_diagnostics, only: debug_model_step
1098 : use ice_timers, only: ice_timer_start, ice_timer_stop, timer_readwrite
1099 :
1100 : logical (kind=log_kind), intent(in) :: flag
1101 :
1102 : integer (kind=int_kind), intent(in) :: &
1103 : recd ! baseline record number
1104 :
1105 : character (char_len_long) :: &
1106 : data_file ! data file to be read
1107 : character (char_len), intent(in) :: &
1108 : fieldname ! field name in netCDF file
1109 :
1110 : integer (kind=int_kind), intent(in) :: &
1111 : field_loc, & ! location of field on staggered grid ! LCOV_EXCL_LINE
1112 : field_type ! type of field (scalar, vector, angle)
1113 :
1114 : real (kind=dbl_kind), dimension(nx_block,ny_block,2,max_blocks), &
1115 : intent(out) :: & ! LCOV_EXCL_LINE
1116 : field_data ! 2 values needed for interpolation
1117 :
1118 : ! local variables
1119 : integer (kind=int_kind) :: &
1120 : fid ! file id for netCDF routines
1121 :
1122 : character(len=*), parameter :: subname = '(read_data_nc_hycom)'
1123 :
1124 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
1125 :
1126 0 : call ice_timer_start(timer_readwrite) ! reading/writing
1127 :
1128 0 : if (istep1 > debug_model_step) debug_forcing = .true. !! debugging
1129 :
1130 0 : if (my_task==master_task .and. (debug_forcing)) then
1131 0 : write(nu_diag,*) ' ', trim(data_file)
1132 : endif
1133 :
1134 0 : if (flag) then
1135 :
1136 0 : call ice_open_nc (data_file, fid)
1137 : !-----------------------------------------------------------------
1138 : ! read data
1139 : !-----------------------------------------------------------------
1140 : call ice_read_nc &
1141 : (fid, recd , fieldname, field_data(:,:,1,:), debug_forcing, & ! LCOV_EXCL_LINE
1142 0 : field_loc, field_type)
1143 :
1144 : call ice_read_nc &
1145 : (fid, recd+1, fieldname, field_data(:,:,2,:), debug_forcing, & ! LCOV_EXCL_LINE
1146 0 : field_loc, field_type)
1147 :
1148 0 : call ice_close_nc(fid)
1149 :
1150 : endif ! flag
1151 :
1152 0 : call ice_timer_stop(timer_readwrite) ! reading/writing
1153 :
1154 0 : end subroutine read_data_nc_hycom
1155 :
1156 : !=======================================================================
1157 :
1158 0 : subroutine read_clim_data (readflag, recd, ixm, ixx, ixp, &
1159 : data_file, field_data, & ! LCOV_EXCL_LINE
1160 : field_loc, field_type)
1161 :
1162 : ! Read data needed for interpolation, as in read_data.
1163 : ! Assume a one-year cycle of climatological data, so that there is
1164 : ! no need to get data from other years or to extrapolate data beyond
1165 : ! the forcing time period.
1166 :
1167 : use ice_diagnostics, only: debug_model_step
1168 :
1169 : logical (kind=log_kind),intent(in) :: readflag
1170 :
1171 : integer (kind=int_kind), intent(in) :: &
1172 : recd , & ! baseline record number ! LCOV_EXCL_LINE
1173 : ixm,ixx,ixp ! record numbers of 3 data values
1174 : ! relative to recd
1175 :
1176 : character (char_len_long), intent(in) :: data_file
1177 :
1178 : integer (kind=int_kind), intent(in) :: &
1179 : field_loc, & ! location of field on staggered grid ! LCOV_EXCL_LINE
1180 : field_type ! type of field (scalar, vector, angle)
1181 :
1182 : real (kind=dbl_kind), dimension(nx_block,ny_block,2,max_blocks), intent(inout) :: &
1183 : field_data ! 2 values needed for interpolation
1184 :
1185 : ! local variables
1186 :
1187 : integer (kind=int_kind) :: &
1188 : nbits , & ! = 32 for single precision, 64 for double ! LCOV_EXCL_LINE
1189 : nrec , & ! record number to read ! LCOV_EXCL_LINE
1190 : arg ! value of time argument in field_data
1191 :
1192 : character(len=*), parameter :: subname = '(read_clim_data)'
1193 :
1194 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
1195 :
1196 0 : call ice_timer_start(timer_readwrite) ! reading/writing
1197 :
1198 0 : nbits = 64 ! double precision data
1199 :
1200 0 : if (istep1 > debug_model_step) debug_forcing = .true. !! debugging
1201 :
1202 0 : if (my_task==master_task .and. (debug_forcing)) &
1203 0 : write(nu_diag,*) ' ', trim(data_file)
1204 :
1205 0 : if (readflag) then
1206 :
1207 : !-----------------------------------------------------------------
1208 : ! read data
1209 : !-----------------------------------------------------------------
1210 :
1211 0 : call ice_open (nu_forcing, data_file, nbits)
1212 :
1213 0 : arg = 0
1214 0 : if (ixm /= -99) then
1215 0 : arg = 1
1216 0 : nrec = recd + ixm
1217 0 : call ice_read (nu_forcing, nrec, field_data(:,:,arg,:), &
1218 0 : 'rda8', debug_forcing, field_loc, field_type)
1219 : endif
1220 :
1221 0 : arg = arg + 1
1222 0 : nrec = recd + ixx
1223 0 : call ice_read (nu_forcing, nrec, field_data(:,:,arg,:), &
1224 0 : 'rda8', debug_forcing, field_loc, field_type)
1225 :
1226 0 : if (ixp /= -99) then
1227 0 : arg = arg + 1
1228 0 : nrec = recd + ixp
1229 0 : call ice_read (nu_forcing, nrec, field_data(:,:,arg,:), &
1230 0 : 'rda8', debug_forcing, field_loc, field_type)
1231 : endif
1232 :
1233 0 : if (my_task == master_task) close (nu_forcing)
1234 : endif ! readflag
1235 :
1236 0 : call ice_timer_stop(timer_readwrite) ! reading/writing
1237 :
1238 0 : end subroutine read_clim_data
1239 :
1240 : !=======================================================================
1241 :
1242 0 : subroutine read_clim_data_nc (readflag, recd, ixm, ixx, ixp, &
1243 : data_file, fieldname, field_data, & ! LCOV_EXCL_LINE
1244 : field_loc, field_type)
1245 :
1246 : ! Read data needed for interpolation, as in read_data.
1247 : ! Assume a one-year cycle of climatological data, so that there is
1248 : ! no need to get data from other years or to extrapolate data beyond
1249 : ! the forcing time period.
1250 :
1251 : use ice_diagnostics, only: debug_model_step
1252 :
1253 : logical (kind=log_kind),intent(in) :: readflag
1254 :
1255 : integer (kind=int_kind), intent(in) :: &
1256 : recd , & ! baseline record number ! LCOV_EXCL_LINE
1257 : ixm,ixx,ixp ! record numbers of 3 data values
1258 : ! relative to recd
1259 :
1260 : character (char_len_long), intent(in) :: data_file
1261 :
1262 : character (char_len), intent(in) :: &
1263 : fieldname ! field name in netCDF file
1264 :
1265 : integer (kind=int_kind), intent(in) :: &
1266 : field_loc, & ! location of field on staggered grid ! LCOV_EXCL_LINE
1267 : field_type ! type of field (scalar, vector, angle)
1268 :
1269 : real (kind=dbl_kind), dimension(nx_block,ny_block,2,max_blocks), intent(out) :: &
1270 : field_data ! 2 values needed for interpolation
1271 :
1272 : ! local variables
1273 :
1274 : integer (kind=int_kind) :: &
1275 : nrec , & ! record number to read ! LCOV_EXCL_LINE
1276 : arg , & ! value of time argument in field_data ! LCOV_EXCL_LINE
1277 : fid ! file id for netCDF routines
1278 :
1279 : character(len=*), parameter :: subname = '(read_clim_data_nc)'
1280 :
1281 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
1282 :
1283 0 : call ice_timer_start(timer_readwrite) ! reading/writing
1284 :
1285 0 : if (istep1 > debug_model_step) debug_forcing = .true. !! debugging
1286 :
1287 0 : if (my_task==master_task .and. (debug_forcing)) &
1288 0 : write(nu_diag,*) ' ', trim(data_file)
1289 :
1290 0 : if (readflag) then
1291 :
1292 : !-----------------------------------------------------------------
1293 : ! read data
1294 : !-----------------------------------------------------------------
1295 :
1296 0 : call ice_open_nc (data_file, fid)
1297 :
1298 0 : arg = 0
1299 0 : if (ixm /= -99) then
1300 0 : arg = 1
1301 0 : nrec = recd + ixm
1302 : call ice_read_nc &
1303 : (fid, nrec, fieldname, field_data(:,:,arg,:), & ! LCOV_EXCL_LINE
1304 0 : debug_forcing, field_loc, field_type)
1305 : endif
1306 :
1307 0 : arg = arg + 1
1308 0 : nrec = recd + ixx
1309 : call ice_read_nc &
1310 : (fid, nrec, fieldname, field_data(:,:,arg,:), & ! LCOV_EXCL_LINE
1311 0 : debug_forcing, field_loc, field_type)
1312 :
1313 0 : if (ixp /= -99) then
1314 0 : arg = arg + 1
1315 0 : nrec = recd + ixp
1316 : call ice_read_nc &
1317 : (fid, nrec, fieldname, field_data(:,:,arg,:), & ! LCOV_EXCL_LINE
1318 0 : debug_forcing, field_loc, field_type)
1319 : endif
1320 :
1321 0 : if (my_task == master_task) call ice_close_nc (fid)
1322 : endif ! readflag
1323 :
1324 0 : call ice_timer_stop(timer_readwrite) ! reading/writing
1325 :
1326 0 : end subroutine read_clim_data_nc
1327 :
1328 : !=======================================================================
1329 :
1330 0 : subroutine interp_coeff_monthly (recslot)
1331 :
1332 : ! Compute coefficients for interpolating monthly data to current time step.
1333 :
1334 : integer (kind=int_kind), intent(in) :: &
1335 : recslot ! slot (1 or 2) for current record
1336 :
1337 : ! local variables
1338 :
1339 : real (kind=dbl_kind) :: &
1340 : secday , & ! seconds in day ! LCOV_EXCL_LINE
1341 : tt , & ! days elapsed in current year ! LCOV_EXCL_LINE
1342 0 : t1, t2 ! days elapsed at month midpoint
1343 :
1344 : real (kind=dbl_kind) :: &
1345 0 : daymid(0:13) ! month mid-points
1346 :
1347 : character(len=*), parameter :: subname = '(interp_coeff_monthly)'
1348 :
1349 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
1350 :
1351 0 : call icepack_query_parameters(secday_out=secday)
1352 0 : call icepack_warnings_flush(nu_diag)
1353 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
1354 0 : file=__FILE__, line=__LINE__)
1355 :
1356 0 : daymid(1:13) = 14._dbl_kind ! time frame ends 0 sec into day 15
1357 0 : daymid(0) = 14._dbl_kind - daymo(12) ! Dec 15, 0 sec
1358 :
1359 : ! compute days since Jan 1, 00h, yday is the day counter for the year
1360 0 : tt = real(yday-1,kind=dbl_kind) + real(msec,kind=dbl_kind)/secday
1361 :
1362 : ! Find neighboring times
1363 :
1364 0 : if (recslot==2) then ! first half of month
1365 0 : t2 = daycal(mmonth) + daymid(mmonth) ! midpoint, current month
1366 0 : if (mmonth == 1) then
1367 0 : t1 = daymid(0) ! Dec 15 (0 sec)
1368 : else
1369 0 : t1 = daycal(mmonth-1) + daymid(mmonth-1) ! midpoint, previous month
1370 : endif
1371 : else ! second half of month
1372 0 : t1 = daycal(mmonth) + daymid(mmonth) ! midpoint, current month
1373 0 : t2 = daycal(mmonth+1) + daymid(mmonth+1)! day 15 of next month (0 sec)
1374 : endif
1375 :
1376 0 : if (tt < t1 .or. tt > t2) then
1377 0 : write(nu_diag,*) subname,' ERROR in tt',tt,t1,t2
1378 : call abort_ice (error_message=subname//' ERROR in tt', &
1379 0 : file=__FILE__, line=__LINE__)
1380 : endif
1381 :
1382 : ! Compute coefficients
1383 0 : c1intp = (t2 - tt) / (t2 - t1)
1384 0 : c2intp = c1 - c1intp
1385 :
1386 0 : end subroutine interp_coeff_monthly
1387 :
1388 : !=======================================================================
1389 :
1390 0 : subroutine interp_coeff (recnum, recslot, secint, dataloc)
1391 :
1392 : ! Compute coefficients for interpolating data to current time step.
1393 : ! Works for any data interval that divides evenly into a
1394 : ! year (daily, 6-hourly, etc.)
1395 : ! Use interp_coef_monthly for monthly data.
1396 :
1397 : integer (kind=int_kind), intent(in) :: &
1398 : recnum , & ! record number for current data value ! LCOV_EXCL_LINE
1399 : recslot , & ! spline slot for current record ! LCOV_EXCL_LINE
1400 : dataloc ! = 1 for data located in middle of time interval
1401 : ! = 2 for date located at end of time interval
1402 :
1403 : real (kind=dbl_kind), intent(in) :: &
1404 : secint ! seconds in data interval
1405 :
1406 : ! local variables
1407 :
1408 : real (kind=dbl_kind) :: &
1409 0 : secday ! seconds in a day
1410 :
1411 : real (kind=dbl_kind) :: &
1412 : tt , & ! seconds elapsed in current year ! LCOV_EXCL_LINE
1413 : t1, t2 , & ! seconds elapsed at data points ! LCOV_EXCL_LINE
1414 0 : rcnum ! recnum => dbl_kind
1415 :
1416 : character(len=*), parameter :: subname = '(interp_coeff)'
1417 :
1418 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
1419 :
1420 0 : call icepack_query_parameters(secday_out=secday)
1421 0 : call icepack_warnings_flush(nu_diag)
1422 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
1423 0 : file=__FILE__, line=__LINE__)
1424 :
1425 : ! compute seconds since Jan 1, 00h, yday is the day counter for the year
1426 0 : tt = real(yday-1,kind=dbl_kind)*secday + real(msec,kind=dbl_kind)
1427 :
1428 : ! Find neighboring times
1429 0 : rcnum = real(recnum,kind=dbl_kind)
1430 0 : if (recslot==2) then ! current record goes in slot 2
1431 0 : if (dataloc==1) then ! data located at middle of interval
1432 0 : t2 = (rcnum-p5)*secint
1433 : else ! data located at end of interval
1434 0 : t2 = rcnum*secint
1435 : endif
1436 0 : t1 = t2 - secint ! - 1 interval
1437 : else ! recslot = 1
1438 0 : if (dataloc==1) then ! data located at middle of interval
1439 0 : t1 = (rcnum-p5)*secint
1440 : else
1441 0 : t1 = rcnum*secint ! data located at end of interval
1442 : endif
1443 0 : t2 = t1 + secint ! + 1 interval
1444 : endif
1445 :
1446 : ! Compute coefficients
1447 0 : c1intp = abs((t2 - tt) / (t2 - t1))
1448 0 : c2intp = c1 - c1intp
1449 :
1450 : if (local_debug .and. my_task == master_task) then
1451 : write(nu_diag,*) subname,'fdbg yday,sec = ',yday,msec
1452 : write(nu_diag,*) subname,'fdbg tt = ',tt
1453 : write(nu_diag,*) subname,'fdbg c12intp = ',c1intp,c2intp
1454 : endif
1455 :
1456 0 : end subroutine interp_coeff
1457 :
1458 : !=======================================================================
1459 :
1460 0 : subroutine interp_coeff2 (tt, t1, t2)
1461 :
1462 : ! Compute coefficients for interpolating data to current time step.
1463 : ! Works for any data interval using decimal daynumbers
1464 :
1465 : ! local variables
1466 : real (kind=dbl_kind), intent(in) :: &
1467 : tt , & ! current decimal daynumber ! LCOV_EXCL_LINE
1468 : t1, t2 ! first+last decimal daynumber
1469 : character(len=*), parameter :: subname = '(interp_coeff2)'
1470 :
1471 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
1472 :
1473 : ! Compute coefficients
1474 0 : c1intp = abs((t2 - tt) / (t2 - t1))
1475 0 : c2intp = c1 - c1intp
1476 :
1477 0 : end subroutine interp_coeff2
1478 :
1479 : !=======================================================================
1480 :
1481 11616 : subroutine interpolate_data (field_data, field)
1482 :
1483 : ! Linear interpolation
1484 :
1485 : ! author: Elizabeth C. Hunke, LANL
1486 :
1487 : use ice_domain, only: nblocks
1488 :
1489 : real (kind=dbl_kind), dimension(nx_block,ny_block,2,max_blocks), intent(in) :: &
1490 : field_data ! 2 values used for interpolation
1491 :
1492 : real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks), intent(out) :: &
1493 : field ! interpolated field
1494 :
1495 : ! local variables
1496 :
1497 : integer (kind=int_kind) :: i,j, iblk
1498 :
1499 : character(len=*), parameter :: subname = '(interpolate data)'
1500 :
1501 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
1502 :
1503 11520 : !$OMP PARALLEL DO PRIVATE(iblk,i,j)
1504 192 : do iblk = 1, nblocks
1505 23040 : do j = 1, ny_block
1506 1166880 : do i = 1, nx_block
1507 : field(i,j,iblk) = c1intp * field_data(i,j,1,iblk) &
1508 1166784 : + c2intp * field_data(i,j,2,iblk)
1509 : enddo
1510 : enddo
1511 : enddo
1512 : !$OMP END PARALLEL DO
1513 :
1514 11616 : end subroutine interpolate_data
1515 :
1516 : !=======================================================================
1517 :
1518 0 : subroutine interpolate_wavespec_data (field_data, field)
1519 :
1520 : ! Linear interpolation
1521 :
1522 : ! author: Elizabeth C. Hunke, LANL
1523 :
1524 : use ice_domain, only: nblocks
1525 :
1526 : real (kind=dbl_kind), dimension(nx_block,ny_block,nfreq,2,max_blocks), intent(in) :: &
1527 : field_data ! 2 values used for interpolation
1528 :
1529 : real (kind=dbl_kind), dimension(nx_block,ny_block,nfreq,max_blocks), intent(out) :: &
1530 : field ! interpolated field
1531 :
1532 : ! local variables
1533 :
1534 : integer (kind=int_kind) :: i,j, iblk, freq
1535 :
1536 : character(len=*), parameter :: subname = '(interpolate data)'
1537 :
1538 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j)
1539 0 : do iblk = 1, nblocks
1540 0 : do j = 1, ny_block
1541 0 : do i = 1, nx_block
1542 0 : do freq = 1, nfreq
1543 : field(i,j,freq,iblk) = c1intp * field_data(i,j,freq,1,iblk) &
1544 0 : + c2intp * field_data(i,j,freq,2,iblk)
1545 : enddo
1546 : enddo
1547 : enddo
1548 : enddo
1549 : !$OMP END PARALLEL DO
1550 :
1551 0 : end subroutine interpolate_wavespec_data
1552 :
1553 :
1554 : !=======================================================================
1555 :
1556 8733 : subroutine file_year (data_file, yr)
1557 :
1558 : ! Construct the correct name of the atmospheric data file
1559 : ! to be read, given the year and assuming the naming convention
1560 : ! that filenames end with 'yyyy.dat' or 'yyyy.r' or 'yyyy.nc'.
1561 :
1562 : character (char_len_long), intent(inout) :: data_file
1563 :
1564 : integer (kind=int_kind), intent(in) :: yr
1565 :
1566 : character (char_len_long) :: tmpname
1567 :
1568 : integer (kind=int_kind) :: i
1569 :
1570 : character(len=*), parameter :: subname = '(file_year)'
1571 :
1572 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
1573 :
1574 8733 : if (trim(atm_data_type) == 'hadgem') then ! netcdf
1575 0 : i = index(data_file,'.nc') - 5
1576 0 : tmpname = data_file
1577 0 : write(data_file,'(a,i4.4,a)') tmpname(1:i), yr, '.nc'
1578 8733 : elseif (index(trim(atm_data_type),'JRA55') > 0) then ! netcdf
1579 8733 : i = index(data_file,'.nc') - 5
1580 8733 : tmpname = data_file
1581 8733 : write(data_file,'(a,i4.4,a)') tmpname(1:i), yr, '.nc'
1582 : else ! LANL/NCAR naming convention
1583 0 : i = index(data_file,'.dat') - 5
1584 0 : tmpname = data_file
1585 0 : write(data_file,'(a,i4.4,a)') tmpname(1:i), yr, '.dat'
1586 : endif
1587 :
1588 8733 : end subroutine file_year
1589 :
1590 : !=======================================================================
1591 :
1592 22944 : subroutine prepare_forcing (nx_block, ny_block, &
1593 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
1594 : hm, & ! LCOV_EXCL_LINE
1595 : Tair, fsw, & ! LCOV_EXCL_LINE
1596 : cldf, flw, & ! LCOV_EXCL_LINE
1597 : frain, fsnow, & ! LCOV_EXCL_LINE
1598 : Qa, rhoa, & ! LCOV_EXCL_LINE
1599 : uatm, vatm, & ! LCOV_EXCL_LINE
1600 : strax, stray, & ! LCOV_EXCL_LINE
1601 : zlvl, wind, & ! LCOV_EXCL_LINE
1602 : swvdr, swvdf, & ! LCOV_EXCL_LINE
1603 : swidr, swidf, & ! LCOV_EXCL_LINE
1604 : potT, ANGLET, & ! LCOV_EXCL_LINE
1605 : Tsfc, sst, & ! LCOV_EXCL_LINE
1606 : aice)
1607 :
1608 : integer (kind=int_kind), intent(in) :: &
1609 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
1610 : ilo,ihi,jlo,jhi ! beginning and end of physical domain
1611 :
1612 : real (kind=dbl_kind), dimension(nx_block,ny_block), intent(in) :: &
1613 : ANGLET , & ! ANGLE converted to T-cells ! LCOV_EXCL_LINE
1614 : Tsfc , & ! ice skin temperature ! LCOV_EXCL_LINE
1615 : sst , & ! sea surface temperature ! LCOV_EXCL_LINE
1616 : aice , & ! ice area fraction ! LCOV_EXCL_LINE
1617 : hm ! land mask
1618 :
1619 : real (kind=dbl_kind), dimension(nx_block,ny_block), intent(inout) :: &
1620 : fsw , & ! incoming shortwave radiation (W/m^2) ! LCOV_EXCL_LINE
1621 : cldf , & ! cloud fraction ! LCOV_EXCL_LINE
1622 : frain , & ! rainfall rate (kg/m^2 s) ! LCOV_EXCL_LINE
1623 : fsnow , & ! snowfall rate (kg/m^2 s) ! LCOV_EXCL_LINE
1624 : Tair , & ! air temperature (K) ! LCOV_EXCL_LINE
1625 : Qa , & ! specific humidity (kg/kg) ! LCOV_EXCL_LINE
1626 : rhoa , & ! air density (kg/m^3) ! LCOV_EXCL_LINE
1627 : uatm , & ! wind velocity components (m/s) ! LCOV_EXCL_LINE
1628 : vatm , & ! LCOV_EXCL_LINE
1629 : strax , & ! wind stress components (N/m^2) ! LCOV_EXCL_LINE
1630 : stray , & ! LCOV_EXCL_LINE
1631 : zlvl , & ! atm level height (m) ! LCOV_EXCL_LINE
1632 : wind , & ! wind speed (m/s) ! LCOV_EXCL_LINE
1633 : flw , & ! incoming longwave radiation (W/m^2) ! LCOV_EXCL_LINE
1634 : swvdr , & ! sw down, visible, direct (W/m^2) ! LCOV_EXCL_LINE
1635 : swvdf , & ! sw down, visible, diffuse (W/m^2) ! LCOV_EXCL_LINE
1636 : swidr , & ! sw down, near IR, direct (W/m^2) ! LCOV_EXCL_LINE
1637 : swidf , & ! sw down, near IR, diffuse (W/m^2) ! LCOV_EXCL_LINE
1638 : potT ! air potential temperature (K)
1639 :
1640 : ! local variables
1641 :
1642 : integer (kind=int_kind) :: &
1643 : i, j
1644 :
1645 5760 : real (kind=dbl_kind) :: workx, worky, &
1646 5760 : precip_factor, zlvl0, secday, Tffresh, puny
1647 :
1648 : logical (kind=log_kind) :: calc_strair
1649 :
1650 : character(len=*), parameter :: subname = '(prepare_forcing)'
1651 :
1652 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
1653 :
1654 22944 : call icepack_query_parameters(Tffresh_out=Tffresh, puny_out=puny)
1655 22944 : call icepack_query_parameters(secday_out=secday)
1656 22944 : call icepack_query_parameters(calc_strair_out=calc_strair)
1657 22944 : call icepack_warnings_flush(nu_diag)
1658 22944 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
1659 0 : file=__FILE__, line=__LINE__)
1660 :
1661 707688 : do j = jlo, jhi
1662 13826928 : do i = ilo, ihi
1663 :
1664 13119240 : zlvl0 = c10 ! default
1665 :
1666 : !-----------------------------------------------------------------
1667 : ! make sure interpolated values are physically realistic
1668 : !-----------------------------------------------------------------
1669 13119240 : cldf (i,j) = max(min(cldf(i,j),c1),c0)
1670 13119240 : fsw (i,j) = max(fsw(i,j),c0)
1671 13119240 : fsnow(i,j) = max(fsnow(i,j),c0)
1672 13119240 : rhoa (i,j) = max(rhoa(i,j),c0)
1673 13803984 : Qa (i,j) = max(Qa(i,j),c0)
1674 :
1675 : ! if (rhoa(i,j) .lt. puny) rhoa(i,j) = 1.3_dbl_kind
1676 : ! if (Tair(i,j) .lt. puny) Tair(i,j) = Tffresh
1677 : ! if (Qa(i,j) .lt. puny) Qa(i,j) = 0.0035_dbl_kind
1678 : enddo ! i
1679 : enddo ! j
1680 :
1681 : !-----------------------------------------------------------------
1682 : ! calculations specific to datasets
1683 : !-----------------------------------------------------------------
1684 :
1685 22944 : if (trim(atm_data_type) == 'ncar') then
1686 :
1687 : ! precip is in mm/month
1688 :
1689 0 : zlvl0 = c10
1690 :
1691 0 : do j = jlo, jhi
1692 0 : do i = ilo, ihi
1693 : ! correct known biases in NCAR data (as in CESM latm)
1694 0 : Qa (i,j) = Qa (i,j) * 0.94_dbl_kind
1695 0 : fsw(i,j) = fsw(i,j) * 0.92_dbl_kind
1696 :
1697 : ! downward longwave as in Parkinson and Washington (1979)
1698 0 : call longwave_parkinson_washington(Tair(i,j), cldf(i,j), &
1699 0 : flw(i,j))
1700 : enddo
1701 : enddo
1702 :
1703 22944 : elseif (trim(atm_data_type) == 'oned') then ! rectangular grid
1704 :
1705 : ! precip is in kg/m^2/s
1706 :
1707 0 : zlvl0 = c10
1708 :
1709 0 : do j = jlo, jhi
1710 0 : do i = ilo, ihi
1711 :
1712 : !-----------------------------------------------------------------
1713 : ! compute downward longwave as in Parkinson and Washington (1979)
1714 : !-----------------------------------------------------------------
1715 :
1716 : ! downward longwave as in Parkinson and Washington (1979)
1717 0 : call longwave_parkinson_washington(Tair(i,j), cldf(i,j), &
1718 0 : flw(i,j))
1719 :
1720 : ! longwave based on Rosati and Miyakoda, JPO 18, p. 1607 (1988)
1721 : ! call longwave_rosati_miyakoda(cldf(i,j), Tsfc(i,j), &
1722 : ! aice(i,j), sst(i,j), & ! LCOV_EXCL_LINE
1723 : ! Qa(i,j), Tair(i,j), & ! LCOV_EXCL_LINE
1724 : ! hm(i,j), flw(i,j))
1725 : enddo
1726 : enddo
1727 :
1728 : endif ! atm_data_type
1729 :
1730 : !-----------------------------------------------------------------
1731 : ! Compute other fields needed by model
1732 : !-----------------------------------------------------------------
1733 :
1734 : ! convert precipitation units to kg/m^2 s
1735 22944 : if (trim(precip_units) == 'mm_per_month') then
1736 5760 : precip_factor = c12/(secday*real(days_per_year,kind=dbl_kind))
1737 17184 : elseif (trim(precip_units) == 'mm_per_day') then
1738 0 : precip_factor = c1/secday
1739 17184 : elseif (trim(precip_units) == 'mm_per_sec' .or. &
1740 : trim(precip_units) == 'mks') then
1741 17184 : precip_factor = c1 ! mm/sec = kg/m^2 s
1742 0 : elseif (trim(precip_units) == 'm_per_sec') then
1743 0 : precip_factor = c1000
1744 : endif
1745 :
1746 707688 : do j = jlo, jhi
1747 13826928 : do i = ilo, ihi
1748 :
1749 13119240 : zlvl(i,j) = zlvl0
1750 13119240 : potT(i,j) = Tair(i,j)
1751 :
1752 : ! divide shortwave into spectral bands
1753 13119240 : swvdr(i,j) = fsw(i,j)*frcvdr ! visible direct
1754 13119240 : swvdf(i,j) = fsw(i,j)*frcvdf ! visible diffuse
1755 13119240 : swidr(i,j) = fsw(i,j)*frcidr ! near IR direct
1756 13119240 : swidf(i,j) = fsw(i,j)*frcidf ! near IR diffuse
1757 :
1758 : ! convert precipitation units to kg/m^2 s
1759 13803984 : fsnow(i,j) = fsnow(i,j) * precip_factor
1760 : enddo ! i
1761 : enddo ! j
1762 :
1763 : ! determine whether precip is rain or snow
1764 : ! HadGEM forcing provides separate snowfall and rainfall rather
1765 : ! than total precipitation
1766 22944 : if (trim(atm_data_type) /= 'hadgem') then
1767 :
1768 707688 : do j = jlo, jhi
1769 13826928 : do i = ilo, ihi
1770 13119240 : frain(i,j) = c0
1771 13803984 : if (Tair(i,j) >= Tffresh) then
1772 4301404 : frain(i,j) = fsnow(i,j)
1773 4301404 : fsnow(i,j) = c0
1774 : endif
1775 : enddo ! i
1776 : enddo ! j
1777 :
1778 : endif
1779 :
1780 22944 : if (calc_strair) then
1781 :
1782 22944 : if (rotate_wind) then
1783 707688 : do j = jlo, jhi
1784 13826928 : do i = ilo, ihi
1785 13119240 : wind(i,j) = sqrt(uatm(i,j)**2 + vatm(i,j)**2)
1786 : !-----------------------------------------------------------------
1787 : ! Rotate zonal/meridional vectors to local coordinates.
1788 : ! Velocity comes in on T grid, but is oriented geographically ---
1789 : ! need to rotate to pop-grid FIRST using ANGLET
1790 : ! then interpolate to the U-cell centers (otherwise we
1791 : ! interpolate across the pole).
1792 : ! Use ANGLET which is on the T grid !
1793 : ! Atmo variables are needed in T cell centers in subroutine
1794 : ! atmo_boundary_layer, and are interpolated to the U grid later as
1795 : ! necessary.
1796 : !-----------------------------------------------------------------
1797 13119240 : workx = uatm(i,j) ! wind velocity, m/s
1798 13119240 : worky = vatm(i,j)
1799 4176000 : uatm (i,j) = workx*cos(ANGLET(i,j)) & ! convert to POP grid
1800 13119240 : + worky*sin(ANGLET(i,j)) ! note uatm, vatm, wind
1801 4176000 : vatm (i,j) = worky*cos(ANGLET(i,j)) & ! are on the T-grid here
1802 13803984 : - workx*sin(ANGLET(i,j))
1803 : enddo ! i
1804 : enddo ! j
1805 : else ! not rotated
1806 0 : do j = jlo, jhi
1807 0 : do i = ilo, ihi
1808 0 : wind(i,j) = sqrt(uatm(i,j)**2 + vatm(i,j)**2)
1809 : enddo ! i
1810 : enddo ! j
1811 : endif ! rotated
1812 :
1813 : else ! strax, stray, wind are read from files
1814 :
1815 0 : if (rotate_wind) then
1816 0 : do j = jlo, jhi
1817 0 : do i = ilo, ihi
1818 0 : workx = strax(i,j) ! wind stress
1819 0 : worky = stray(i,j)
1820 0 : strax(i,j) = workx*cos(ANGLET(i,j)) & ! convert to POP grid
1821 0 : + worky*sin(ANGLET(i,j)) ! note strax, stray, wind
1822 0 : stray(i,j) = worky*cos(ANGLET(i,j)) & ! are on the T-grid here
1823 0 : - workx*sin(ANGLET(i,j))
1824 : enddo ! i
1825 : enddo ! j
1826 : else ! not rotated
1827 : ! wind (speed) is already read from file, so all is in place
1828 : endif ! rotated
1829 :
1830 : endif ! calc_strair
1831 :
1832 22944 : end subroutine prepare_forcing
1833 :
1834 : !=======================================================================
1835 :
1836 0 : subroutine longwave_parkinson_washington(Tair, cldf, flw)
1837 :
1838 : ! compute downward longwave as in Parkinson and Washington (1979)
1839 : ! (for now)
1840 : ! Parkinson, C. L. and W. M. Washington (1979),
1841 : ! Large-scale numerical-model of sea ice,
1842 : ! JGR, 84, 311-337, doi:10.1029/JC084iC01p00311
1843 :
1844 : real(kind=dbl_kind), intent(in) :: &
1845 : Tair , & ! air temperature (K) ! LCOV_EXCL_LINE
1846 : cldf ! cloud fraction
1847 :
1848 : real(kind=dbl_kind), intent(out) :: &
1849 : flw ! incoming longwave radiation (W/m^2)
1850 :
1851 : real(kind=dbl_kind) :: &
1852 0 : Tffresh, stefan_boltzmann
1853 :
1854 : character(len=*), parameter :: subname = '(longwave_parkinson_washington)'
1855 :
1856 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
1857 :
1858 : call icepack_query_parameters(Tffresh_out=Tffresh, &
1859 0 : stefan_boltzmann_out=stefan_boltzmann)
1860 0 : call icepack_warnings_flush(nu_diag)
1861 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
1862 0 : file=__FILE__, line=__LINE__)
1863 :
1864 : flw = stefan_boltzmann*Tair**4 &
1865 : * (c1 - 0.261_dbl_kind & ! LCOV_EXCL_LINE
1866 : * exp(-7.77e-4_dbl_kind*(Tffresh - Tair)**2)) & ! LCOV_EXCL_LINE
1867 0 : * (c1 + 0.275_dbl_kind*cldf)
1868 :
1869 0 : end subroutine longwave_parkinson_washington
1870 :
1871 : !=======================================================================
1872 :
1873 0 : subroutine longwave_rosati_miyakoda(cldf, Tsfc, &
1874 : aice, sst, & ! LCOV_EXCL_LINE
1875 : Qa, Tair, & ! LCOV_EXCL_LINE
1876 : hm, flw)
1877 :
1878 : ! based on
1879 : ! Rosati, A. and K. Miyakoda (1988),
1880 : ! A general-circulation model for upper ocean simulation,
1881 : ! J. Physical Oceanography, 18, 1601-1626,
1882 : ! doi:10.1175/1520-0485(1988)018<1601:AGCMFU>2.0.CO;2
1883 :
1884 : real(kind=dbl_kind), intent(in) :: &
1885 : cldf , & ! cloud fraction ! LCOV_EXCL_LINE
1886 : Tsfc , & ! ice skin temperature ! LCOV_EXCL_LINE
1887 : aice , & ! ice area fraction ! LCOV_EXCL_LINE
1888 : sst , & ! sea surface temperature ! LCOV_EXCL_LINE
1889 : Qa , & ! specific humidity (kg/kg) ! LCOV_EXCL_LINE
1890 : Tair , & ! air temperature (K) ! LCOV_EXCL_LINE
1891 : hm ! land mask
1892 :
1893 : real(kind=dbl_kind), intent(out) :: &
1894 : flw ! incoming longwave radiation (W/m^2)
1895 :
1896 : real(kind=dbl_kind) :: &
1897 : fcc , & ! cloudiness modification ! LCOV_EXCL_LINE
1898 : sstk , & ! ice/ocean surface temperature (K) ! LCOV_EXCL_LINE
1899 : rtea , & ! square root of the vapour pressure ! LCOV_EXCL_LINE
1900 : ptem , & ! potential air temperature (K) ! LCOV_EXCL_LINE
1901 0 : qlwm
1902 :
1903 : real(kind=dbl_kind) :: &
1904 0 : Tffresh, stefan_boltzmann, emissivity
1905 :
1906 : character(len=*), parameter :: subname = '(longwave_rosati_miyakoda)'
1907 :
1908 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
1909 :
1910 : call icepack_query_parameters(Tffresh_out=Tffresh, &
1911 : stefan_boltzmann_out=stefan_boltzmann, & ! LCOV_EXCL_LINE
1912 0 : emissivity_out=emissivity)
1913 0 : call icepack_warnings_flush(nu_diag)
1914 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
1915 0 : file=__FILE__, line=__LINE__)
1916 :
1917 0 : fcc = c1 - 0.8_dbl_kind * cldf
1918 : sstk = (Tsfc * aice &
1919 0 : + sst * (c1 - aice)) + Tffresh
1920 : rtea = sqrt(c1000*Qa / &
1921 0 : (0.622_dbl_kind+0.378_dbl_kind*Qa))
1922 0 : ptem = Tair ! get this from stability?
1923 : qlwm = ptem * ptem * ptem &
1924 : * ( ptem*(0.39_dbl_kind-0.05_dbl_kind*rtea)*fcc & ! LCOV_EXCL_LINE
1925 0 : + c4*(sstk-ptem) )
1926 0 : flw = emissivity*stefan_boltzmann * ( sstk**4 - qlwm )
1927 0 : flw = flw * hm ! land mask
1928 :
1929 0 : end subroutine longwave_rosati_miyakoda
1930 :
1931 : !=======================================================================
1932 : ! NCAR atmospheric forcing
1933 : !=======================================================================
1934 :
1935 0 : subroutine ncar_files (yr)
1936 :
1937 : ! Construct filenames based on the LANL naming conventions for NCAR data.
1938 : ! Edit for other directory structures or filenames.
1939 : ! Note: The year number in these filenames does not matter, because
1940 : ! subroutine file_year will insert the correct year.
1941 : ! Note: atm_data_dir may have NCAR_bulk or not
1942 : !
1943 : ! atm_data_type should be 'ncar'
1944 : ! atm_dat_dir should be ${CICE_DATA_root}/forcing/$grid/[NCAR_bulk,'']
1945 : ! atm_data_dir should be set to ${CICE_DATA_root}/forcing/$grid/[JRA55,JRA55do,'']
1946 : ! NCAR_bulk at the end of the atm_data_dir is optional to provide backwards
1947 : ! compatibility and if not included, will be appended automaticaly.
1948 : ! The grid is typically gx1, gx3, tx1, or similar.
1949 :
1950 : integer (kind=int_kind), intent(in) :: &
1951 : yr ! current forcing year
1952 :
1953 : character (char_len_long) :: &
1954 : atm_data_dir_extra ! atm_dat_dir extra if needed
1955 :
1956 : integer (kind=int_kind) :: &
1957 : strind ! string index
1958 :
1959 : character(len=*), parameter :: subname = '(ncar_files)'
1960 :
1961 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
1962 :
1963 : ! decide whether NCAR_bulk is part of atm_data_dir and set atm_data_dir_extra
1964 0 : atm_data_dir_extra = '/NCAR_bulk'
1965 0 : strind = index(trim(atm_data_dir),'NCAR_bulk')
1966 0 : if (strind > 0) then
1967 0 : atm_data_dir_extra = ''
1968 : endif
1969 :
1970 0 : fsw_file = trim(atm_data_dir)//trim(atm_data_dir_extra)//'/MONTHLY/swdn.1996.dat'
1971 0 : call file_year(fsw_file,yr)
1972 :
1973 0 : flw_file = trim(atm_data_dir)//trim(atm_data_dir_extra)//'/MONTHLY/cldf.1996.dat'
1974 0 : call file_year(flw_file,yr)
1975 :
1976 0 : rain_file = trim(atm_data_dir)//trim(atm_data_dir_extra)//'/MONTHLY/prec.1996.dat'
1977 0 : call file_year(rain_file,yr)
1978 :
1979 0 : uwind_file = trim(atm_data_dir)//trim(atm_data_dir_extra)//'/4XDAILY/u_10.1996.dat'
1980 0 : call file_year(uwind_file,yr)
1981 :
1982 0 : vwind_file = trim(atm_data_dir)//trim(atm_data_dir_extra)//'/4XDAILY/v_10.1996.dat'
1983 0 : call file_year(vwind_file,yr)
1984 :
1985 0 : tair_file = trim(atm_data_dir)//trim(atm_data_dir_extra)//'/4XDAILY/t_10.1996.dat'
1986 0 : call file_year(tair_file,yr)
1987 :
1988 0 : humid_file = trim(atm_data_dir)//trim(atm_data_dir_extra)//'/4XDAILY/q_10.1996.dat'
1989 0 : call file_year(humid_file,yr)
1990 :
1991 0 : rhoa_file = trim(atm_data_dir)//trim(atm_data_dir_extra)//'/4XDAILY/dn10.1996.dat'
1992 0 : call file_year(rhoa_file,yr)
1993 :
1994 0 : if (my_task == master_task) then
1995 0 : write (nu_diag,*) ' '
1996 0 : write (nu_diag,*) 'Forcing data year =', fyear
1997 0 : write (nu_diag,*) 'Atmospheric data files:'
1998 0 : write (nu_diag,'(3a)') trim(fsw_file)
1999 0 : write (nu_diag,'(3a)') trim(flw_file)
2000 0 : write (nu_diag,'(3a)') trim(rain_file)
2001 0 : write (nu_diag,'(3a)') trim(uwind_file)
2002 0 : write (nu_diag,'(3a)') trim(vwind_file)
2003 0 : write (nu_diag,'(3a)') trim(tair_file)
2004 0 : write (nu_diag,'(3a)') trim(humid_file)
2005 0 : write (nu_diag,'(3a)') trim(rhoa_file)
2006 : endif ! master_task
2007 :
2008 0 : end subroutine ncar_files
2009 :
2010 : !=======================================================================
2011 :
2012 0 : subroutine ncar_data
2013 :
2014 : use ice_flux, only: fsw, fsnow, Tair, uatm, vatm, rhoa, Qa
2015 :
2016 : integer (kind=int_kind) :: &
2017 : ixm,ixx,ixp , & ! record numbers for neighboring months ! LCOV_EXCL_LINE
2018 : recnum , & ! record number ! LCOV_EXCL_LINE
2019 : maxrec , & ! maximum record number ! LCOV_EXCL_LINE
2020 : recslot , & ! spline slot for current record ! LCOV_EXCL_LINE
2021 : dataloc , & ! = 1 for data located in middle of time interval ! LCOV_EXCL_LINE
2022 : ! = 2 for date located at end of time interval
2023 : midmonth ! middle day of month
2024 :
2025 : real (kind=dbl_kind) :: &
2026 : secday, & ! number of seconds in day ! LCOV_EXCL_LINE
2027 0 : sec6hr ! number of seconds in 6 hours
2028 :
2029 : logical (kind=log_kind) :: readm, read6
2030 :
2031 : character(len=*), parameter :: subname = '(ncar_data)'
2032 :
2033 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
2034 :
2035 0 : call icepack_query_parameters(secday_out=secday)
2036 0 : call icepack_warnings_flush(nu_diag)
2037 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
2038 0 : file=__FILE__, line=__LINE__)
2039 :
2040 : !-------------------------------------------------------------------
2041 : ! monthly data
2042 : !
2043 : ! Assume that monthly data values are located in the middle of the
2044 : ! month.
2045 : !-------------------------------------------------------------------
2046 :
2047 0 : midmonth = 15 ! data is given on 15th of every month
2048 : ! midmonth = fix(p5 * real(daymo(mmonth))) ! exact middle
2049 :
2050 : ! Compute record numbers for surrounding months
2051 0 : maxrec = 12
2052 0 : ixm = mod(mmonth+maxrec-2,maxrec) + 1
2053 0 : ixp = mod(mmonth, maxrec) + 1
2054 0 : if (mday >= midmonth) ixm = -99 ! other two points will be used
2055 0 : if (mday < midmonth) ixp = -99
2056 :
2057 : ! Determine whether interpolation will use values 1:2 or 2:3
2058 : ! recslot = 2 means we use values 1:2, with the current value (2)
2059 : ! in the second slot
2060 : ! recslot = 1 means we use values 2:3, with the current value (2)
2061 : ! in the first slot
2062 0 : recslot = 1 ! latter half of month
2063 0 : if (mday < midmonth) recslot = 2 ! first half of month
2064 :
2065 : ! Find interpolation coefficients
2066 0 : call interp_coeff_monthly (recslot)
2067 :
2068 : ! Read 2 monthly values
2069 0 : readm = .false.
2070 0 : if (istep==1 .or. (mday==midmonth .and. msec==0)) readm = .true.
2071 :
2072 0 : if (trim(atm_data_format) == 'bin') then
2073 : call read_data (readm, 0, fyear, ixm, mmonth, ixp, &
2074 : maxrec, fsw_file, fsw_data, & ! LCOV_EXCL_LINE
2075 0 : field_loc_center, field_type_scalar)
2076 : call read_data (readm, 0, fyear, ixm, mmonth, ixp, &
2077 : maxrec, flw_file, cldf_data, & ! LCOV_EXCL_LINE
2078 0 : field_loc_center, field_type_scalar)
2079 : call read_data (readm, 0, fyear, ixm, mmonth, ixp, &
2080 : maxrec, rain_file, fsnow_data, & ! LCOV_EXCL_LINE
2081 0 : field_loc_center, field_type_scalar)
2082 : else
2083 : call abort_ice (error_message=subname//'nonbinary atm_data_format unavailable', &
2084 0 : file=__FILE__, line=__LINE__)
2085 : ! The routine exists, for example:
2086 : ! call read_data_nc (readm, 0, fyear, ixm, mmonth, ixp, &
2087 : ! maxrec, fsw_file, 'fsw', fsw_data, & ! LCOV_EXCL_LINE
2088 : ! field_loc_center, field_type_scalar)
2089 : ! call read_data_nc (readm, 0, fyear, ixm, mmonth, ixp, &
2090 : ! maxrec, flw_file, 'cldf',cldf_data, & ! LCOV_EXCL_LINE
2091 : ! field_loc_center, field_type_scalar)
2092 : ! call read_data_nc (readm, 0, fyear, ixm, mmonth, ixp, &
2093 : ! maxrec, rain_file,'prec',fsnow_data, & ! LCOV_EXCL_LINE
2094 : ! field_loc_center, field_type_scalar)
2095 : endif
2096 :
2097 : ! Interpolate to current time step
2098 0 : call interpolate_data (fsw_data, fsw)
2099 0 : call interpolate_data (cldf_data, cldf)
2100 0 : call interpolate_data (fsnow_data, fsnow)
2101 :
2102 : !-------------------------------------------------------------------
2103 : ! 6-hourly data
2104 : !
2105 : ! Assume that the 6-hourly value is located at the end of the
2106 : ! 6-hour period. This is the convention for NCEP reanalysis data.
2107 : ! E.g. record 1 gives conditions at 6 am GMT on 1 January.
2108 : !-------------------------------------------------------------------
2109 :
2110 0 : dataloc = 2 ! data located at end of interval
2111 0 : sec6hr = secday/c4 ! seconds in 6 hours
2112 0 : maxrec = 1460 ! 365*4
2113 :
2114 : ! current record number
2115 0 : recnum = 4*int(yday) - 3 + int(real(msec,kind=dbl_kind)/sec6hr)
2116 :
2117 : ! Compute record numbers for surrounding data
2118 :
2119 0 : ixm = mod(recnum+maxrec-2,maxrec) + 1
2120 0 : ixx = mod(recnum-1, maxrec) + 1
2121 : ! ixp = mod(recnum, maxrec) + 1
2122 :
2123 : ! Compute interpolation coefficients
2124 : ! If data is located at the end of the time interval, then the
2125 : ! data value for the current record always goes in slot 2.
2126 :
2127 0 : recslot = 2
2128 0 : ixp = -99
2129 0 : call interp_coeff (recnum, recslot, sec6hr, dataloc)
2130 :
2131 : ! Read
2132 0 : read6 = .false.
2133 0 : if (istep==1 .or. oldrecnum /= recnum) read6 = .true.
2134 :
2135 0 : if (trim(atm_data_format) == 'bin') then
2136 : call read_data (read6, 0, fyear, ixm, ixx, ixp, &
2137 : maxrec, tair_file, Tair_data, & ! LCOV_EXCL_LINE
2138 0 : field_loc_center, field_type_scalar)
2139 : call read_data (read6, 0, fyear, ixm, ixx, ixp, &
2140 : maxrec, uwind_file, uatm_data, & ! LCOV_EXCL_LINE
2141 0 : field_loc_center, field_type_vector)
2142 : call read_data (read6, 0, fyear, ixm, ixx, ixp, &
2143 : maxrec, vwind_file, vatm_data, & ! LCOV_EXCL_LINE
2144 0 : field_loc_center, field_type_vector)
2145 : call read_data (read6, 0, fyear, ixm, ixx, ixp, &
2146 : maxrec, rhoa_file, rhoa_data, & ! LCOV_EXCL_LINE
2147 0 : field_loc_center, field_type_scalar)
2148 : call read_data (read6, 0, fyear, ixm, ixx, ixp, &
2149 : maxrec, humid_file, Qa_data, & ! LCOV_EXCL_LINE
2150 0 : field_loc_center, field_type_scalar)
2151 : else
2152 : call abort_ice (error_message=subname//'nonbinary atm_data_format unavailable', &
2153 0 : file=__FILE__, line=__LINE__)
2154 : endif
2155 :
2156 : ! Interpolate
2157 0 : call interpolate_data (Tair_data, Tair)
2158 0 : call interpolate_data (uatm_data, uatm)
2159 0 : call interpolate_data (vatm_data, vatm)
2160 0 : call interpolate_data (rhoa_data, rhoa)
2161 0 : call interpolate_data (Qa_data, Qa)
2162 :
2163 : ! Save record number for next time step
2164 0 : oldrecnum = recnum
2165 :
2166 0 : end subroutine ncar_data
2167 :
2168 : !=======================================================================
2169 :
2170 21 : subroutine JRA55_files(yr)
2171 :
2172 : ! find the JRA55 files:
2173 : ! This subroutine finds the JRA55 atm forcing files based on settings
2174 : ! in atm_data_type and atm_data_dir. Because the filenames are not
2175 : ! entirely consistent, we need a flexible method.
2176 : !
2177 : ! atm_data_type could be JRA55 or JRA55do with/without _grid appended
2178 : ! atm_data_dir could contain JRA55 or JRA55do or not
2179 : ! actual files could have grid in name in two location or not at all
2180 : !
2181 : ! The files will generally be of the format
2182 : ! $atm_data_type/[JRA55,JRA55do,'']/8XDAILY/[JRA55,JRA55do][_$grid,'']_03hr_forcing[_$grid,'']_$year.nc
2183 : ! The options defined by cnt try several versions of paths/filenames
2184 : ! As a user,
2185 : ! atm_data_type should be set to JRA55, JRA55do, JRA55_xxx, or JRA55do_xxx
2186 : ! where xxx can be any set of characters. The _xxx if included will be ignored.
2187 : ! Historically, these were set to JRA55_gx1 and so forth but the _gx1 is no longer needed
2188 : ! but this is still allowed for backwards compatibility. atm_data_type_prefix
2189 : ! is atm_data_type with _ and everything after _ removed.
2190 : ! atm_data_dir should be set to ${CICE_DATA_root}/forcing/$grid/[JRA55,JRA55do,'']
2191 : ! The [JRA55,JRA55do] at the end of the atm_data_dir is optional to provide backwards
2192 : ! compatibility and if not included, will be appended automaticaly using
2193 : ! the atm_data_type_prefix value. The grid is typically gx1, gx3, tx1, or similar.
2194 : ! In general, we recommend using the following format
2195 : ! atm_data_type = [JRA55,JRA55do]
2196 : ! atm_data_dir = ${CICE_DATA_root}/forcing/$grid
2197 :
2198 : integer (kind=int_kind), intent(in) :: &
2199 : yr ! current forcing year
2200 :
2201 : ! local variables
2202 : character(len=16) :: &
2203 : grd ! gx3, gx1, tx1
2204 :
2205 : character(len=64) :: &
2206 : atm_data_type_prefix ! atm_data_type prefix
2207 :
2208 : integer (kind=int_kind) :: &
2209 : cnt , & ! search for files ! LCOV_EXCL_LINE
2210 : strind ! string index
2211 :
2212 : logical :: &
2213 : exists ! file existance
2214 :
2215 : character(len=*), parameter :: subname = '(JRA55_files)'
2216 :
2217 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
2218 :
2219 : ! this could be JRA55[do] or JRA55[do]_grid, drop the _grid if set
2220 21 : atm_data_type_prefix = trim(atm_data_type)
2221 21 : strind = index(trim(atm_data_type),'_')
2222 21 : if (strind > 0) then
2223 0 : atm_data_type_prefix = atm_data_type(1:strind-1)
2224 : endif
2225 :
2226 : ! check for grid version using fortran INDEX intrinsic
2227 21 : if (index(trim(atm_data_dir),'gx1') > 0) then
2228 0 : grd = 'gx1'
2229 21 : else if (index(trim(atm_data_dir),'gx3') > 0) then
2230 21 : grd = 'gx3'
2231 0 : else if (index(trim(atm_data_dir),'tx1') > 0) then
2232 0 : grd = 'tx1'
2233 : else
2234 0 : call abort_ice(error_message=subname//' unknown grid type')
2235 : endif
2236 :
2237 : ! cnt represents the possible file format options and steps thru them until one is found
2238 21 : exists = .false.
2239 21 : cnt = 1
2240 42 : do while (.not.exists .and. cnt <= 6)
2241 21 : if (cnt == 1) uwind_file = trim(atm_data_dir)//'/'//trim(atm_data_type_prefix)// &
2242 21 : '/8XDAILY/'//trim(atm_data_type_prefix)//'_'//trim(grd)//'_03hr_forcing_2005.nc'
2243 :
2244 21 : if (cnt == 2) uwind_file = trim(atm_data_dir)//'/'//trim(atm_data_type_prefix)// &
2245 0 : '/8XDAILY/'//trim(atm_data_type_prefix)//'_03hr_forcing_'//trim(grd)//'_2005.nc'
2246 :
2247 21 : if (cnt == 3) uwind_file = trim(atm_data_dir)//'/'//trim(atm_data_type_prefix)// &
2248 0 : '/8XDAILY/'//trim(atm_data_type_prefix)// '_03hr_forcing_2005.nc'
2249 :
2250 21 : if (cnt == 4) uwind_file = trim(atm_data_dir)// &
2251 0 : '/8XDAILY/'//trim(atm_data_type_prefix)//'_'//trim(grd)//'_03hr_forcing_2005.nc'
2252 :
2253 21 : if (cnt == 5) uwind_file = trim(atm_data_dir)// &
2254 0 : '/8XDAILY/'//trim(atm_data_type_prefix)//'_03hr_forcing_'//trim(grd)//'_2005.nc'
2255 :
2256 21 : if (cnt == 6) uwind_file = trim(atm_data_dir)// &
2257 0 : '/8XDAILY/'//trim(atm_data_type_prefix)// '_03hr_forcing_2005.nc'
2258 :
2259 21 : call file_year(uwind_file,yr)
2260 21 : INQUIRE(FILE=uwind_file,EXIST=exists)
2261 : ! if (my_task == master_task) then
2262 : ! write(nu_diag,*) subname,cnt,exists,trim(uwind_file)
2263 : ! endif
2264 21 : cnt = cnt + 1
2265 : enddo
2266 :
2267 21 : if (.not.exists) then
2268 0 : call abort_ice(error_message=subname//' could not find forcing file')
2269 : endif
2270 :
2271 21 : if (my_task == master_task) then
2272 5 : write (nu_diag,'(2a)') ' '
2273 5 : write (nu_diag,'(2a)') subname,'Atmospheric data files:'
2274 5 : write (nu_diag,'(2a)') subname,trim(uwind_file)
2275 : endif
2276 :
2277 21 : end subroutine JRA55_files
2278 :
2279 : !=======================================================================
2280 :
2281 2904 : subroutine JRA55_data
2282 :
2283 : use ice_blocks, only: block, get_block
2284 : use ice_global_reductions, only: global_minval, global_maxval
2285 : use ice_domain, only: nblocks, distrb_info
2286 : use ice_flux, only: fsnow, Tair, uatm, vatm, Qa, fsw, flw
2287 : use ice_grid, only: hm, tmask, umask
2288 : use ice_state, only: aice
2289 : use ice_calendar, only: days_per_year
2290 :
2291 : integer (kind=int_kind) :: &
2292 : ncid , & ! netcdf file id ! LCOV_EXCL_LINE
2293 : i, j, n1 , & ! LCOV_EXCL_LINE
2294 : lfyear , & ! local year value ! LCOV_EXCL_LINE
2295 : recnum , & ! record number ! LCOV_EXCL_LINE
2296 : maxrec , & ! maximum record number ! LCOV_EXCL_LINE
2297 : iblk ! block index
2298 :
2299 : integer (kind=int_kind), save :: &
2300 : frec_info(2,2) = -99 ! remember prior values to reduce reading
2301 : ! first dim is yr, recnum
2302 : ! second dim is data1 data2
2303 :
2304 : real (kind=dbl_kind) :: &
2305 : sec3hr , & ! number of seconds in 3 hours ! LCOV_EXCL_LINE
2306 : secday , & ! number of seconds in day ! LCOV_EXCL_LINE
2307 : eps, tt , & ! for interpolation coefficients ! LCOV_EXCL_LINE
2308 : Tffresh , & ! LCOV_EXCL_LINE
2309 1440 : vmin, vmax
2310 :
2311 : character(len=64) :: fieldname !netcdf field name
2312 : character (char_len_long) :: uwind_file_old
2313 : character(len=*), parameter :: subname = '(JRA55_data)'
2314 :
2315 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
2316 :
2317 2904 : call icepack_query_parameters(Tffresh_out=Tffresh)
2318 2904 : call icepack_query_parameters(secday_out=secday)
2319 2904 : call icepack_warnings_flush(nu_diag)
2320 2904 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
2321 0 : file=__FILE__, line=__LINE__)
2322 :
2323 2904 : sec3hr = secday/c8 ! seconds in 3 hours
2324 2904 : maxrec = days_per_year * 8
2325 :
2326 : if (local_debug .and. my_task == master_task) then
2327 : write(nu_diag,*) subname,'fdbg dpy, maxrec = ',days_per_year,maxrec
2328 : endif
2329 :
2330 : !-------------------------------------------------------------------
2331 : ! 3-hourly data
2332 : ! states are instantaneous, 1st record is 00z Jan 1
2333 : ! fluxes are 3 hour averages, 1st record is 00z-03z Jan 1
2334 : ! interpolate states, do not interpolate fluxes
2335 : !-------------------------------------------------------------------
2336 : ! File is NETCDF with winds in NORTH and EAST direction
2337 : ! file variable names are:
2338 : ! glbrad (shortwave W/m^2), 3 hr average
2339 : ! dlwsfc (longwave W/m^2), 3 hr average
2340 : ! wndewd (eastward wind m/s), instantaneous
2341 : ! wndnwd (northward wind m/s), instantaneous
2342 : ! airtmp (air temperature K), instantaneous
2343 : ! spchmd (specific humidity kg/kg), instantaneous
2344 : ! ttlpcp (precipitation kg/m s-1), 3 hr average
2345 : !-------------------------------------------------------------------
2346 :
2347 2904 : uwind_file_old = uwind_file
2348 2904 : if (uwind_file /= uwind_file_old .and. my_task == master_task) then
2349 0 : write(nu_diag,'(2a)') subname,' reading forcing file = ',trim(uwind_file)
2350 : endif
2351 :
2352 2904 : call ice_open_nc(uwind_file,ncid)
2353 :
2354 8712 : do n1 = 1, 2
2355 :
2356 5808 : lfyear = fyear
2357 5808 : call file_year(uwind_file,lfyear)
2358 5808 : if (n1 == 1) then
2359 2904 : recnum = 8*int(yday) - 7 + int(real(msec,kind=dbl_kind)/sec3hr)
2360 2904 : if (my_task == master_task .and. (recnum <= 2 .or. recnum >= maxrec-1)) then
2361 20 : write(nu_diag,'(3a)') subname,' reading forcing file 1st ts = ',trim(uwind_file)
2362 : endif
2363 2904 : elseif (n1 == 2) then
2364 2904 : recnum = 8*int(yday) - 7 + int(real(msec,kind=dbl_kind)/sec3hr) + 1
2365 2904 : if (recnum > maxrec) then
2366 0 : lfyear = fyear + 1 ! next year
2367 0 : if (lfyear > fyear_final) lfyear = fyear_init
2368 0 : recnum = 1
2369 0 : call file_year(uwind_file,lfyear)
2370 0 : if (my_task == master_task) then
2371 0 : write(nu_diag,'(3a)') subname,' reading forcing file 2nd ts = ',trim(uwind_file)
2372 : endif
2373 0 : call ice_close_nc(ncid)
2374 0 : call ice_open_nc(uwind_file,ncid)
2375 : endif
2376 : endif
2377 :
2378 : if (local_debug .and. my_task == master_task) then
2379 : write(nu_diag,*) subname,'fdbg read recnum = ',recnum,n1
2380 : endif
2381 :
2382 : ! to reduce reading, check whether it's the same data as last read
2383 :
2384 5808 : if (lfyear /= frec_info(1,n1) .or. recnum /= frec_info(2,n1)) then
2385 :
2386 : ! check whether we can copy values from 2 to 1, should be faster than reading
2387 : ! can only do this from 2 to 1 or 1 to 2 without setting up a temporary
2388 : ! it's more likely that the values from data2 when time advances are needed in data1
2389 : ! compare n1=1 year/record with data from last timestep at n1=2
2390 :
2391 1978 : if (n1 == 1 .and. lfyear == frec_info(1,2) .and. recnum == frec_info(2,2)) then
2392 3119648 : Tair_data(:,:,1,:) = Tair_data(:,:,2,:)
2393 3119648 : uatm_data(:,:,1,:) = uatm_data(:,:,2,:)
2394 3119648 : vatm_data(:,:,1,:) = vatm_data(:,:,2,:)
2395 3119648 : Qa_data(:,:,1,:) = Qa_data(:,:,2,:)
2396 3119648 : fsw_data(:,:,1,:) = fsw_data(:,:,2,:)
2397 3119648 : flw_data(:,:,1,:) = flw_data(:,:,2,:)
2398 3119648 : fsnow_data(:,:,1,:) = fsnow_data(:,:,2,:)
2399 : else
2400 :
2401 1010 : fieldname = 'airtmp'
2402 0 : call ice_read_nc(ncid,recnum,fieldname,Tair_data(:,:,n1,:),local_debug, &
2403 : field_loc=field_loc_center, & ! LCOV_EXCL_LINE
2404 1010 : field_type=field_type_scalar)
2405 :
2406 1010 : fieldname = 'wndewd'
2407 0 : call ice_read_nc(ncid,recnum,fieldname,uatm_data(:,:,n1,:),local_debug, &
2408 : field_loc=field_loc_center, & ! LCOV_EXCL_LINE
2409 1010 : field_type=field_type_scalar)
2410 :
2411 1010 : fieldname = 'wndnwd'
2412 0 : call ice_read_nc(ncid,recnum,fieldname,vatm_data(:,:,n1,:),local_debug, &
2413 : field_loc=field_loc_center, & ! LCOV_EXCL_LINE
2414 1010 : field_type=field_type_scalar)
2415 :
2416 1010 : fieldname = 'spchmd'
2417 0 : call ice_read_nc(ncid,recnum,fieldname,Qa_data(:,:,n1,:),local_debug, &
2418 : field_loc=field_loc_center, & ! LCOV_EXCL_LINE
2419 1010 : field_type=field_type_scalar)
2420 :
2421 1010 : fieldname = 'glbrad'
2422 0 : call ice_read_nc(ncid,recnum,fieldname,fsw_data(:,:,n1,:),local_debug, &
2423 : field_loc=field_loc_center, & ! LCOV_EXCL_LINE
2424 1010 : field_type=field_type_scalar)
2425 :
2426 1010 : fieldname = 'dlwsfc'
2427 0 : call ice_read_nc(ncid,recnum,fieldname,flw_data(:,:,n1,:),local_debug, &
2428 : field_loc=field_loc_center, & ! LCOV_EXCL_LINE
2429 1010 : field_type=field_type_scalar)
2430 :
2431 1010 : fieldname = 'ttlpcp'
2432 0 : call ice_read_nc(ncid,recnum,fieldname,fsnow_data(:,:,n1,:),local_debug, &
2433 : field_loc=field_loc_center, & ! LCOV_EXCL_LINE
2434 1010 : field_type=field_type_scalar)
2435 : endif ! copy data from n1=2 from last timestep to n1=1
2436 : endif ! input data is same as last timestep
2437 :
2438 5808 : frec_info(1,n1) = lfyear
2439 8712 : frec_info(2,n1) = recnum
2440 :
2441 : enddo ! n1
2442 :
2443 2904 : call ice_close_nc(ncid)
2444 :
2445 : ! reset uwind_file to original year
2446 2904 : call file_year(uwind_file,fyear)
2447 :
2448 : ! Compute interpolation coefficients
2449 2904 : eps = 1.0e-6
2450 2904 : tt = real(mod(msec,nint(sec3hr)),kind=dbl_kind)
2451 2904 : c2intp = tt / sec3hr
2452 2904 : if (c2intp < c0 .and. c2intp > c0-eps) c2intp = c0
2453 2904 : if (c2intp > c1 .and. c2intp < c1+eps) c2intp = c1
2454 2904 : c1intp = 1.0_dbl_kind - c2intp
2455 2904 : if (c2intp < c0 .or. c2intp > c1) then
2456 0 : write(nu_diag,*) subname,' ERROR: c2intp = ',c2intp
2457 : call abort_ice (error_message=subname//' ERROR: c2intp out of range', &
2458 0 : file=__FILE__, line=__LINE__)
2459 : endif
2460 : if (local_debug .and. my_task == master_task) then
2461 : write(nu_diag,*) subname,'fdbg c12intp = ',c1intp,c2intp
2462 : endif
2463 :
2464 : ! Interpolate
2465 2904 : call interpolate_data (Tair_data, Tair)
2466 2904 : call interpolate_data (uatm_data, uatm)
2467 2904 : call interpolate_data (vatm_data, vatm)
2468 2904 : call interpolate_data (Qa_data, Qa)
2469 : ! use 3 hr average for heat flux and precip fields, no interpolation
2470 : ! call interpolate_data (fsw_data, fsw)
2471 : ! call interpolate_data (flw_data, flw)
2472 : ! call interpolate_data (fsnow_data, fsnow)
2473 9358944 : fsw(:,:,:) = fsw_data(:,:,1,:)
2474 9358944 : flw(:,:,:) = flw_data(:,:,1,:)
2475 9358944 : fsnow(:,:,:) = fsnow_data(:,:,1,:)
2476 :
2477 2880 : !$OMP PARALLEL DO PRIVATE(iblk,i,j)
2478 48 : do iblk = 1, nblocks
2479 : ! limit summer Tair values where ice is present
2480 2856 : do j = 1, ny_block
2481 291720 : do i = 1, nx_block
2482 291696 : if (aice(i,j,iblk) > p1) Tair(i,j,iblk) = min(Tair(i,j,iblk), Tffresh+p1)
2483 : enddo
2484 : enddo
2485 :
2486 5760 : do j = 1, ny_block
2487 291720 : do i = 1, nx_block
2488 288864 : Qa (i,j,iblk) = Qa (i,j,iblk) * hm(i,j,iblk)
2489 288864 : Tair(i,j,iblk) = Tair(i,j,iblk) * hm(i,j,iblk)
2490 288864 : uatm(i,j,iblk) = uatm(i,j,iblk) * hm(i,j,iblk)
2491 288864 : vatm(i,j,iblk) = vatm(i,j,iblk) * hm(i,j,iblk)
2492 288864 : fsw (i,j,iblk) = fsw (i,j,iblk) * hm(i,j,iblk)
2493 288864 : flw (i,j,iblk) = flw (i,j,iblk) * hm(i,j,iblk)
2494 291696 : fsnow(i,j,iblk) = fsnow (i,j,iblk) * hm(i,j,iblk)
2495 : enddo
2496 : enddo
2497 :
2498 : enddo ! iblk
2499 : !$OMP END PARALLEL DO
2500 :
2501 2904 : if (debug_forcing .or. local_debug) then
2502 0 : if (my_task.eq.master_task) write (nu_diag,*) subname,'fdbg JRA55_bulk_data'
2503 0 : vmin = global_minval(fsw,distrb_info,tmask)
2504 0 : vmax = global_maxval(fsw,distrb_info,tmask)
2505 0 : if (my_task.eq.master_task) write (nu_diag,*) subname,'fdbg fsw',vmin,vmax
2506 0 : vmin = global_minval(flw,distrb_info,tmask)
2507 0 : vmax = global_maxval(flw,distrb_info,tmask)
2508 0 : if (my_task.eq.master_task) write (nu_diag,*) subname,'fdbg flw',vmin,vmax
2509 0 : vmin =global_minval(fsnow,distrb_info,tmask)
2510 0 : vmax =global_maxval(fsnow,distrb_info,tmask)
2511 0 : if (my_task.eq.master_task) write (nu_diag,*) subname,'fdbg fsnow',vmin,vmax
2512 0 : vmin = global_minval(Tair,distrb_info,tmask)
2513 0 : vmax = global_maxval(Tair,distrb_info,tmask)
2514 0 : if (my_task.eq.master_task) write (nu_diag,*) subname,'fdbg Tair',vmin,vmax
2515 0 : vmin = global_minval(uatm,distrb_info,umask)
2516 0 : vmax = global_maxval(uatm,distrb_info,umask)
2517 0 : if (my_task.eq.master_task) write (nu_diag,*) subname,'fdbg uatm',vmin,vmax
2518 0 : vmin = global_minval(vatm,distrb_info,umask)
2519 0 : vmax = global_maxval(vatm,distrb_info,umask)
2520 0 : if (my_task.eq.master_task) write (nu_diag,*) subname,'fdbg vatm',vmin,vmax
2521 0 : vmin = global_minval(Qa,distrb_info,tmask)
2522 0 : vmax = global_maxval(Qa,distrb_info,tmask)
2523 0 : if (my_task.eq.master_task) write (nu_diag,*) subname,'fdbg Qa',vmin,vmax
2524 : endif ! debug_forcing
2525 :
2526 5808 : end subroutine JRA55_data
2527 :
2528 : !=======================================================================
2529 : !
2530 : ! AOMIP shortwave forcing
2531 : ! standard calculation using solar declination angle
2532 : ! then shortwave is reduced using a function of cloud fraction
2533 :
2534 0 : subroutine compute_shortwave(nx_block, ny_block, &
2535 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
2536 0 : TLON, TLAT, hm, Qa, cldf, fsw)
2537 :
2538 : !---!-------------------------------------------------------------------
2539 : !---!-------------------------------------------------------------------
2540 :
2541 : integer (kind=int_kind), intent(in) :: &
2542 : nx_block, ny_block, & ! block dimensions ! LCOV_EXCL_LINE
2543 : ilo,ihi,jlo,jhi ! beginning and end of physical domain
2544 :
2545 : real (kind=dbl_kind), dimension(nx_block,ny_block), intent(in) :: &
2546 : TLON, TLAT , & ! longitude, latitude ! LCOV_EXCL_LINE
2547 : Qa , & ! specific humidity ! LCOV_EXCL_LINE
2548 : cldf , & ! cloud fraction ! LCOV_EXCL_LINE
2549 : hm ! land mask
2550 :
2551 : real (kind=dbl_kind), dimension(nx_block,ny_block), intent(inout) :: &
2552 : fsw ! shortwave
2553 :
2554 : real (kind=dbl_kind) :: &
2555 : hour_angle, & ! LCOV_EXCL_LINE
2556 : solar_time, & ! LCOV_EXCL_LINE
2557 : declin , & ! LCOV_EXCL_LINE
2558 : cosZ , & ! LCOV_EXCL_LINE
2559 : e, d , & ! LCOV_EXCL_LINE
2560 : sw0 , & ! LCOV_EXCL_LINE
2561 : secday , & ! LCOV_EXCL_LINE
2562 : pi , & ! LCOV_EXCL_LINE
2563 : lontmp , & ! LCOV_EXCL_LINE
2564 0 : deg2rad
2565 :
2566 : integer (kind=int_kind) :: &
2567 : i, j
2568 :
2569 : character(len=*), parameter :: subname = '(compute_shortwave)'
2570 :
2571 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
2572 :
2573 0 : call icepack_query_parameters(secday_out=secday, pi_out=pi)
2574 0 : call icepack_warnings_flush(nu_diag)
2575 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
2576 0 : file=__FILE__, line=__LINE__)
2577 :
2578 0 : do j=jlo,jhi
2579 0 : do i=ilo,ihi
2580 0 : deg2rad = pi/c180
2581 : ! solar_time = mod(real(msec,kind=dbl_kind),secday)/c3600 &
2582 : ! + c12*sin(p5*TLON(i,j))
2583 :
2584 : ! Convert longitude to range of -180 to 180 for LST calculation
2585 :
2586 0 : lontmp = mod(TLON(i,j)/deg2rad,c360)
2587 0 : if (lontmp .gt. c180) lontmp = lontmp - c360
2588 0 : if (lontmp .lt. -c180) lontmp = lontmp + c360
2589 :
2590 : solar_time = mod(real(msec,kind=dbl_kind),secday)/c3600 &
2591 0 : + lontmp/c15
2592 0 : if (solar_time .ge. 24._dbl_kind) solar_time = solar_time - 24._dbl_kind
2593 0 : hour_angle = (c12 - solar_time)*pi/c12
2594 : declin = 23.44_dbl_kind*cos((172._dbl_kind-yday) &
2595 0 : * c2*pi/c365)*deg2rad ! use dayyr instead of c365???
2596 0 : cosZ = sin(TLAT(i,j))*sin(declin) &
2597 0 : + cos(TLAT(i,j))*cos(declin)*cos(hour_angle)
2598 0 : cosZ = max(cosZ,c0)
2599 0 : e = 1.e5*Qa(i,j)/(0.622_dbl_kind + 0.378_dbl_kind*Qa(i,j))
2600 0 : d = (cosZ+2.7_dbl_kind)*e*1.e-5_dbl_kind+1.085_dbl_kind*cosZ+p1
2601 0 : sw0 = 1353._dbl_kind*cosZ**2/d
2602 0 : sw0 = max(sw0,c0)
2603 :
2604 : ! total downward shortwave for cice
2605 0 : Fsw(i,j) = sw0*(c1-p6*cldf(i,j)**3)
2606 0 : Fsw(i,j) = Fsw(i,j)*hm(i,j)
2607 : enddo
2608 : enddo
2609 :
2610 0 : end subroutine compute_shortwave
2611 :
2612 : !=======================================================================
2613 : !
2614 : ! prevents humidity from being super-saturated
2615 :
2616 0 : subroutine Qa_fixLY(nx_block, ny_block, Tair, Qa)
2617 :
2618 : integer (kind=int_kind), intent(in) :: &
2619 : nx_block, ny_block ! block dimensions
2620 :
2621 : real (kind=dbl_kind), dimension(nx_block,ny_block), intent(in) :: &
2622 : Tair ! air temperature
2623 :
2624 : real (kind=dbl_kind), dimension(nx_block,ny_block), intent(inout) :: &
2625 : Qa ! specific humidity
2626 :
2627 : real (kind=dbl_kind), dimension (nx_block,ny_block) :: &
2628 0 : worka
2629 :
2630 : real (kind=dbl_kind) :: &
2631 0 : Tffresh, puny
2632 :
2633 : character(len=*), parameter :: subname = '(Qa_fixLY)'
2634 :
2635 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
2636 :
2637 0 : call icepack_query_parameters(Tffresh_out=Tffresh, puny_out=puny)
2638 0 : call icepack_warnings_flush(nu_diag)
2639 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
2640 0 : file=__FILE__, line=__LINE__)
2641 :
2642 0 : worka = Tair - Tffresh
2643 0 : worka = c2 + (0.7859_dbl_kind + 0.03477_dbl_kind*worka) &
2644 : /(c1 + 0.00412_dbl_kind*worka) & ! 2+ converts ea mb -> Pa ! LCOV_EXCL_LINE
2645 0 : + 0.00422_dbl_kind*worka ! for ice
2646 : ! vapor pressure
2647 0 : worka = (c10**worka) ! saturated
2648 0 : worka = max(worka,puny) ! puny over land to prevent division by zero
2649 : ! specific humidity
2650 0 : worka = 0.622_dbl_kind*worka/(1.e5_dbl_kind-0.378_dbl_kind*worka)
2651 :
2652 0 : Qa = min(Qa, worka)
2653 :
2654 0 : end subroutine Qa_fixLY
2655 :
2656 : !=======================================================================
2657 : ! HadGEM or HadGAM atmospheric forcing
2658 : !=======================================================================
2659 :
2660 0 : subroutine hadgem_files (yr)
2661 :
2662 : ! Construct filenames based on selected model options
2663 : !
2664 : ! Note: The year number in these filenames does not matter, because
2665 : ! subroutine file_year will insert the correct year.
2666 : !
2667 : ! author: Alison McLaren, Met Office
2668 :
2669 : integer (kind=int_kind), intent(in) :: &
2670 : yr ! current forcing year
2671 :
2672 : integer (kind=int_kind) :: &
2673 : n ! thickness category index
2674 :
2675 : logical (kind=log_kind) :: calc_strair, calc_Tsfc
2676 :
2677 : character(len=*), parameter :: subname = '(hadgem_files)'
2678 :
2679 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
2680 :
2681 : call icepack_query_parameters(calc_strair_out=calc_strair, &
2682 0 : calc_Tsfc_out=calc_Tsfc)
2683 0 : call icepack_warnings_flush(nu_diag)
2684 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
2685 0 : file=__FILE__, line=__LINE__)
2686 :
2687 : ! -----------------------------------------------------------
2688 : ! Rainfall and snowfall
2689 : ! -----------------------------------------------------------
2690 :
2691 : snow_file = &
2692 0 : trim(atm_data_dir)//'/MONTHLY/snowfall.1996.nc'
2693 0 : call file_year(snow_file,yr)
2694 :
2695 : rain_file = &
2696 0 : trim(atm_data_dir)//'/MONTHLY/rainfall.1996.nc'
2697 0 : call file_year(rain_file,yr)
2698 :
2699 0 : if (my_task == master_task) then
2700 0 : write (nu_diag,*) ' '
2701 0 : write (nu_diag,*) 'Atmospheric data files:'
2702 0 : write (nu_diag,*) trim(rain_file)
2703 0 : write (nu_diag,*) trim(snow_file)
2704 : endif
2705 :
2706 0 : if (calc_strair) then
2707 :
2708 : ! --------------------------------------------------------
2709 : ! Wind velocity
2710 : ! --------------------------------------------------------
2711 :
2712 : uwind_file = &
2713 0 : trim(atm_data_dir)//'/MONTHLY/u_10.1996.nc'
2714 0 : call file_year(uwind_file,yr)
2715 :
2716 : vwind_file = &
2717 0 : trim(atm_data_dir)//'/MONTHLY/v_10.1996.nc'
2718 0 : call file_year(vwind_file,yr)
2719 :
2720 0 : if (my_task == master_task) then
2721 0 : write (nu_diag,*) trim(uwind_file)
2722 0 : write (nu_diag,*) trim(vwind_file)
2723 : endif
2724 :
2725 : else
2726 :
2727 : ! --------------------------------------------------------
2728 : ! Wind stress
2729 : ! --------------------------------------------------------
2730 :
2731 : strax_file = &
2732 0 : trim(atm_data_dir)//'/MONTHLY/taux.1996.nc'
2733 0 : call file_year(strax_file,yr)
2734 :
2735 : stray_file = &
2736 0 : trim(atm_data_dir)//'/MONTHLY/tauy.1996.nc'
2737 0 : call file_year(stray_file,yr)
2738 :
2739 0 : if (my_task == master_task) then
2740 0 : write (nu_diag,*) trim(strax_file)
2741 0 : write (nu_diag,*) trim(stray_file)
2742 : endif
2743 :
2744 0 : if (calc_Tsfc .or. oceanmixed_ice) then
2745 :
2746 : ! --------------------------------------------------
2747 : ! Wind speed
2748 : ! --------------------------------------------------
2749 :
2750 : wind_file = &
2751 0 : trim(atm_data_dir)//'/MONTHLY/wind_10.1996.nc'
2752 0 : call file_year(wind_file,yr)
2753 :
2754 0 : if (my_task == master_task) then
2755 0 : write (nu_diag,*) trim(wind_file)
2756 : endif
2757 :
2758 : endif ! calc_Tsfc or oceanmixed_ice
2759 :
2760 : endif ! calc_strair
2761 :
2762 : ! --------------------------------------------------------------
2763 : ! Atmosphere properties. Even if these fields are not
2764 : ! being used to force the ice (i.e. calc_Tsfc=.false.), they
2765 : ! are still needed to generate forcing for mixed layer model or
2766 : ! to calculate wind stress
2767 : ! --------------------------------------------------------------
2768 :
2769 0 : if (calc_Tsfc .or. oceanmixed_ice .or. calc_strair) then
2770 :
2771 : fsw_file = &
2772 0 : trim(atm_data_dir)//'/MONTHLY/SW_incoming.1996.nc'
2773 0 : call file_year(fsw_file,yr)
2774 :
2775 : flw_file = &
2776 0 : trim(atm_data_dir)//'/MONTHLY/LW_incoming.1996.nc'
2777 0 : call file_year(flw_file,yr)
2778 :
2779 : tair_file = &
2780 0 : trim(atm_data_dir)//'/MONTHLY/t_10.1996.nc'
2781 0 : call file_year(tair_file,yr)
2782 :
2783 : humid_file = &
2784 0 : trim(atm_data_dir)//'/MONTHLY/q_10.1996.nc'
2785 0 : call file_year(humid_file,yr)
2786 :
2787 : rhoa_file = &
2788 0 : trim(atm_data_dir)//'/MONTHLY/rho_10.1996.nc'
2789 0 : call file_year(rhoa_file,yr)
2790 :
2791 0 : if (my_task == master_task) then
2792 0 : write (nu_diag,*) trim(fsw_file)
2793 0 : write (nu_diag,*) trim(flw_file)
2794 0 : write (nu_diag,*) trim(tair_file)
2795 0 : write (nu_diag,*) trim(humid_file)
2796 0 : write (nu_diag,*) trim(rhoa_file)
2797 : endif ! master_task
2798 :
2799 : endif ! calc_Tsfc or oceanmixed_ice or calc_strair
2800 :
2801 0 : if (.not. calc_Tsfc) then
2802 :
2803 : ! ------------------------------------------------------
2804 : ! Sublimation, topmelt and botmelt
2805 : ! ------------------------------------------------------
2806 :
2807 0 : do n = 1, ncat
2808 :
2809 : ! 'topmelt' = fsurf - fcondtop.
2810 0 : write(topmelt_file(n), '(a,i1,a)') &
2811 0 : trim(atm_data_dir)//'/MONTHLY/topmeltn',n,'.1996.nc'
2812 0 : call file_year(topmelt_file(n),yr)
2813 :
2814 : ! 'botmelt' = fcondtop.
2815 0 : write(botmelt_file(n), '(a,i1,a)') &
2816 0 : trim(atm_data_dir)//'/MONTHLY/botmeltn',n,'.1996.nc'
2817 0 : call file_year(botmelt_file(n),yr)
2818 :
2819 : enddo
2820 :
2821 : ! 'sublim' = - flat / Lsub.
2822 : sublim_file = &
2823 0 : trim(atm_data_dir)//'/MONTHLY/sublim.1996.nc'
2824 0 : call file_year(sublim_file,yr)
2825 :
2826 0 : if (my_task == master_task) then
2827 0 : do n = 1, ncat
2828 0 : write (nu_diag,*) trim(topmelt_file(n))
2829 0 : write (nu_diag,*) trim(botmelt_file(n))
2830 : enddo
2831 0 : write (nu_diag,*) trim(sublim_file)
2832 :
2833 : endif
2834 :
2835 : endif ! .not. calc_Tsfc
2836 :
2837 0 : end subroutine hadgem_files
2838 :
2839 : !=======================================================================
2840 :
2841 : ! read HadGEM or HadGAM atmospheric data
2842 :
2843 0 : subroutine hadgem_data
2844 :
2845 : ! authors: Alison McLaren, Met Office
2846 :
2847 : use ice_domain, only: nblocks
2848 : use ice_flux, only: fsnow, frain, uatm, vatm, strax, stray, wind, &
2849 : fsw, flw, Tair, rhoa, Qa, fcondtopn_f, fsurfn_f, flatn_f
2850 :
2851 : integer (kind=int_kind) :: &
2852 : i, j , & ! horizontal indices ! LCOV_EXCL_LINE
2853 : n , & ! thickness category index ! LCOV_EXCL_LINE
2854 : iblk , & ! block index ! LCOV_EXCL_LINE
2855 : ixm,ixp , & ! record numbers for neighboring months ! LCOV_EXCL_LINE
2856 : maxrec , & ! maximum record number ! LCOV_EXCL_LINE
2857 : recslot , & ! spline slot for current record ! LCOV_EXCL_LINE
2858 : midmonth ! middle day of month
2859 :
2860 : logical (kind=log_kind) :: readm
2861 :
2862 : real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: &
2863 : topmelt, & ! temporary fields ! LCOV_EXCL_LINE
2864 : botmelt, & ! LCOV_EXCL_LINE
2865 0 : sublim
2866 :
2867 : character (char_len) :: &
2868 : fieldname ! field name in netcdf file
2869 :
2870 : real (kind=dbl_kind) :: &
2871 0 : Lsub
2872 :
2873 : logical (kind=log_kind) :: &
2874 : calc_strair, & ! LCOV_EXCL_LINE
2875 : calc_Tsfc
2876 :
2877 : character(len=*), parameter :: subname = '(hadgem_data)'
2878 :
2879 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
2880 :
2881 0 : call icepack_query_parameters(Lsub_out=Lsub)
2882 : call icepack_query_parameters(calc_strair_out=calc_strair, &
2883 0 : calc_Tsfc_out=calc_Tsfc)
2884 0 : call icepack_warnings_flush(nu_diag)
2885 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
2886 0 : file=__FILE__, line=__LINE__)
2887 :
2888 : !-------------------------------------------------------------------
2889 : ! monthly data
2890 : !
2891 : ! Assume that monthly data values are located in the middle of the
2892 : ! month.
2893 : !-------------------------------------------------------------------
2894 :
2895 0 : midmonth = 15 ! data is given on 15th of every month
2896 : ! midmonth = fix(p5 * real(daymo(mmonth))) ! exact middle
2897 :
2898 : ! Compute record numbers for surrounding months
2899 0 : maxrec = 12
2900 0 : ixm = mod(mmonth+maxrec-2,maxrec) + 1
2901 0 : ixp = mod(mmonth, maxrec) + 1
2902 0 : if (mday >= midmonth) ixm = -99 ! other two points will be used
2903 0 : if (mday < midmonth) ixp = -99
2904 :
2905 : ! Determine whether interpolation will use values 1:2 or 2:3
2906 : ! recslot = 2 means we use values 1:2, with the current value (2)
2907 : ! in the second slot
2908 : ! recslot = 1 means we use values 2:3, with the current value (2)
2909 : ! in the first slot
2910 0 : recslot = 1 ! latter half of month
2911 0 : if (mday < midmonth) recslot = 2 ! first half of month
2912 :
2913 : ! Find interpolation coefficients
2914 0 : call interp_coeff_monthly (recslot)
2915 :
2916 : ! Read 2 monthly values
2917 0 : readm = .false.
2918 0 : if (istep==1 .or. (mday==midmonth .and. msec==0)) readm = .true.
2919 :
2920 : ! -----------------------------------------------------------
2921 : ! Rainfall and snowfall
2922 : ! -----------------------------------------------------------
2923 :
2924 0 : fieldname='rainfall'
2925 : call read_data_nc (readm, 0, fyear, ixm, mmonth, ixp, &
2926 : maxrec, rain_file, fieldname, frain_data, & ! LCOV_EXCL_LINE
2927 0 : field_loc_center, field_type_scalar)
2928 0 : fieldname='snowfall'
2929 : call read_data_nc (readm, 0, fyear, ixm, mmonth, ixp, &
2930 : maxrec, snow_file, fieldname, fsnow_data, & ! LCOV_EXCL_LINE
2931 0 : field_loc_center, field_type_scalar)
2932 :
2933 : ! Interpolate to current time step
2934 0 : call interpolate_data (fsnow_data, fsnow)
2935 0 : call interpolate_data (frain_data, frain)
2936 :
2937 0 : if (calc_strair) then
2938 :
2939 : ! --------------------------------------------------------
2940 : ! Wind velocity
2941 : ! --------------------------------------------------------
2942 :
2943 0 : fieldname='u_10'
2944 : call read_data_nc (readm, 0, fyear, ixm, mmonth, ixp, &
2945 : maxrec, uwind_file, fieldname, uatm_data, & ! LCOV_EXCL_LINE
2946 0 : field_loc_center, field_type_vector)
2947 0 : fieldname='v_10'
2948 : call read_data_nc (readm, 0, fyear, ixm, mmonth, ixp, &
2949 : maxrec, vwind_file, fieldname, vatm_data, & ! LCOV_EXCL_LINE
2950 0 : field_loc_center, field_type_vector)
2951 :
2952 : ! Interpolate to current time step
2953 0 : call interpolate_data (uatm_data, uatm)
2954 0 : call interpolate_data (vatm_data, vatm)
2955 :
2956 : else
2957 :
2958 : ! --------------------------------------------------------
2959 : ! Wind stress
2960 : ! --------------------------------------------------------
2961 :
2962 0 : fieldname='taux'
2963 : call read_data_nc (readm, 0, fyear, ixm, mmonth, ixp, &
2964 : maxrec, strax_file, fieldname, strax_data, & ! LCOV_EXCL_LINE
2965 0 : field_loc_center, field_type_vector)
2966 0 : fieldname='tauy'
2967 : call read_data_nc (readm, 0, fyear, ixm, mmonth, ixp, &
2968 : maxrec, stray_file, fieldname, stray_data, & ! LCOV_EXCL_LINE
2969 0 : field_loc_center, field_type_vector)
2970 :
2971 : ! Interpolate to current time step
2972 0 : call interpolate_data (strax_data, strax)
2973 0 : call interpolate_data (stray_data, stray)
2974 :
2975 0 : if (calc_Tsfc .or. oceanmixed_ice) then
2976 :
2977 : ! --------------------------------------------------
2978 : ! Wind speed
2979 : ! --------------------------------------------------
2980 :
2981 0 : fieldname='wind_10'
2982 : call read_data_nc (readm, 0, fyear, ixm, mmonth, ixp, &
2983 : maxrec, wind_file, fieldname, wind_data, & ! LCOV_EXCL_LINE
2984 0 : field_loc_center, field_type_scalar)
2985 :
2986 : ! Interpolate to current time step
2987 0 : call interpolate_data (wind_data, wind)
2988 :
2989 : endif ! calc_Tsfc or oceanmixed_ice
2990 :
2991 : endif ! calc_strair
2992 :
2993 : ! -----------------------------------------------------------
2994 : ! SW incoming, LW incoming, air temperature, density and
2995 : ! humidity at 10m.
2996 : !
2997 : ! Even if these fields are not being used to force the ice
2998 : ! (i.e. calc_Tsfc=.false.), they are still needed to generate
2999 : ! forcing for mixed layer model or to calculate wind stress
3000 : ! -----------------------------------------------------------
3001 :
3002 0 : if (calc_Tsfc .or. oceanmixed_ice .or. calc_strair) then
3003 :
3004 0 : fieldname='SW_incoming'
3005 : call read_data_nc (readm, 0, fyear, ixm, mmonth, ixp, &
3006 : maxrec, fsw_file, fieldname, fsw_data, & ! LCOV_EXCL_LINE
3007 0 : field_loc_center, field_type_scalar)
3008 0 : fieldname='LW_incoming'
3009 : call read_data_nc (readm, 0, fyear, ixm, mmonth, ixp, &
3010 : maxrec, flw_file, fieldname, flw_data, & ! LCOV_EXCL_LINE
3011 0 : field_loc_center, field_type_scalar)
3012 0 : fieldname='t_10'
3013 : call read_data_nc (readm, 0, fyear, ixm, mmonth, ixp, &
3014 : maxrec, tair_file, fieldname, Tair_data, & ! LCOV_EXCL_LINE
3015 0 : field_loc_center, field_type_scalar)
3016 0 : fieldname='rho_10'
3017 : call read_data_nc (readm, 0, fyear, ixm, mmonth, ixp, &
3018 : maxrec, rhoa_file, fieldname, rhoa_data, & ! LCOV_EXCL_LINE
3019 0 : field_loc_center, field_type_scalar)
3020 0 : fieldname='q_10'
3021 : call read_data_nc (readm, 0, fyear, ixm, mmonth, ixp, &
3022 : maxrec, humid_file, fieldname, Qa_data, & ! LCOV_EXCL_LINE
3023 0 : field_loc_center, field_type_scalar)
3024 :
3025 : ! Interpolate onto current timestep
3026 :
3027 0 : call interpolate_data (fsw_data, fsw)
3028 0 : call interpolate_data (flw_data, flw)
3029 0 : call interpolate_data (Tair_data, Tair)
3030 0 : call interpolate_data (rhoa_data, rhoa)
3031 0 : call interpolate_data (Qa_data, Qa)
3032 :
3033 : endif ! calc_Tsfc or oceanmixed_ice or calc_strair
3034 :
3035 0 : if (.not. calc_Tsfc) then
3036 :
3037 : ! ------------------------------------------------------
3038 : ! Sublimation, topmelt and botmelt
3039 : ! ------------------------------------------------------
3040 :
3041 0 : fieldname='sublim'
3042 : call read_data_nc (readm, 0, fyear, ixm, mmonth, ixp, &
3043 : maxrec, sublim_file, fieldname, sublim_data, & ! LCOV_EXCL_LINE
3044 0 : field_loc_center, field_type_scalar)
3045 :
3046 : ! Interpolate to current time step
3047 0 : call interpolate_data (sublim_data, sublim)
3048 :
3049 0 : do n = 1, ncat
3050 0 : write(fieldname, '(a,i1)') 'topmeltn',n
3051 : call read_data_nc (readm, 0, fyear, ixm, mmonth, ixp, &
3052 : maxrec, topmelt_file(n), fieldname, topmelt_data(:,:,:,:,n), & ! LCOV_EXCL_LINE
3053 0 : field_loc_center, field_type_scalar)
3054 :
3055 0 : write(fieldname, '(a,i1)') 'botmeltn',n
3056 : call read_data_nc (readm, 0, fyear, ixm, mmonth, ixp, &
3057 : maxrec, botmelt_file(n), fieldname, botmelt_data(:,:,:,:,n), & ! LCOV_EXCL_LINE
3058 0 : field_loc_center, field_type_scalar)
3059 :
3060 0 : call interpolate_data (topmelt_data(:,:,:,:,n), topmelt)
3061 0 : call interpolate_data (botmelt_data(:,:,:,:,n), botmelt)
3062 :
3063 : !--------------------------------------------------------
3064 : ! Convert from UM variables to CICE variables
3065 : ! topmelt = fsurf - fcondtop
3066 : ! botmelt = fcondtop (as zero layer)
3067 : !
3068 : ! Convert UM sublimation data into CICE LH flux
3069 : ! (sublim = - flatn / Lsub) and have same value for all
3070 : ! categories
3071 : !--------------------------------------------------------
3072 :
3073 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j)
3074 0 : do iblk = 1, nblocks
3075 0 : do j = 1, ny_block
3076 0 : do i = 1, nx_block
3077 0 : fcondtopn_f(i,j,n,iblk) = botmelt(i,j,iblk)
3078 : fsurfn_f(i,j,n,iblk) = topmelt(i,j,iblk) &
3079 0 : + botmelt(i,j,iblk)
3080 0 : flatn_f(i,j,n,iblk) = - sublim(i,j,iblk)*Lsub
3081 : enddo
3082 : enddo
3083 : enddo
3084 : !$OMP END PARALLEL DO
3085 :
3086 : enddo ! ncat
3087 :
3088 : endif ! .not. calc_Tsfc
3089 :
3090 0 : end subroutine hadgem_data
3091 :
3092 : !=======================================================================
3093 : ! monthly forcing
3094 : !=======================================================================
3095 :
3096 0 : subroutine monthly_files (yr)
3097 :
3098 : ! Construct filenames based on the LANL naming conventions for NCAR data.
3099 : ! Edit for other directory structures or filenames.
3100 : ! Note: The year number in these filenames does not matter, because
3101 : ! subroutine file_year will insert the correct year.
3102 :
3103 : ! author: Elizabeth C. Hunke, LANL
3104 :
3105 : integer (kind=int_kind), intent(in) :: &
3106 : yr ! current forcing year
3107 :
3108 : character(len=*), parameter :: subname = '(monthly_files)'
3109 :
3110 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
3111 :
3112 : flw_file = &
3113 0 : trim(atm_data_dir)//'/MONTHLY/cldf.omip.dat'
3114 :
3115 : rain_file = &
3116 0 : trim(atm_data_dir)//'/MONTHLY/prec.nmyr.dat'
3117 :
3118 : tair_file = &
3119 0 : trim(atm_data_dir)//'/MONTHLY/t_10.1996.dat'
3120 0 : call file_year(tair_file,yr)
3121 :
3122 : humid_file = &
3123 0 : trim(atm_data_dir)//'/MONTHLY/q_10.1996.dat'
3124 0 : call file_year(humid_file,yr)
3125 :
3126 : ! stress/speed is used instead of wind components
3127 : strax_file = &
3128 0 : trim(atm_data_dir)//'/MONTHLY/strx.1996.dat'
3129 0 : call file_year(strax_file,yr)
3130 :
3131 : stray_file = &
3132 0 : trim(atm_data_dir)//'/MONTHLY/stry.1996.dat'
3133 0 : call file_year(stray_file,yr)
3134 :
3135 : wind_file = &
3136 0 : trim(atm_data_dir)//'/MONTHLY/wind.1996.dat'
3137 0 : call file_year(wind_file,yr)
3138 :
3139 0 : if (my_task == master_task) then
3140 0 : write (nu_diag,*) ' '
3141 0 : write (nu_diag,*) 'Forcing data year = ', fyear
3142 0 : write (nu_diag,*) 'Atmospheric data files:'
3143 0 : write (nu_diag,*) trim(flw_file)
3144 0 : write (nu_diag,*) trim(rain_file)
3145 0 : write (nu_diag,*) trim(tair_file)
3146 0 : write (nu_diag,*) trim(humid_file)
3147 0 : write (nu_diag,*) trim(uwind_file)
3148 0 : write (nu_diag,*) trim(vwind_file)
3149 : endif ! master_task
3150 :
3151 0 : end subroutine monthly_files
3152 :
3153 : !=======================================================================
3154 : ! read monthly atmospheric data
3155 :
3156 0 : subroutine monthly_data
3157 :
3158 : use ice_blocks, only: block, get_block
3159 : use ice_global_reductions, only: global_minval, global_maxval
3160 : use ice_domain, only: nblocks, distrb_info, blocks_ice
3161 : use ice_flux, only: fsnow, Tair, Qa, wind, strax, stray, fsw
3162 : use ice_grid, only: hm, tlon, tlat, tmask, umask
3163 :
3164 : integer (kind=int_kind) :: &
3165 : i, j , & ! LCOV_EXCL_LINE
3166 : ixm,ixp , & ! record numbers for neighboring months ! LCOV_EXCL_LINE
3167 : maxrec , & ! maximum record number ! LCOV_EXCL_LINE
3168 : recslot , & ! spline slot for current record ! LCOV_EXCL_LINE
3169 : midmonth , & ! middle day of month ! LCOV_EXCL_LINE
3170 : iblk , & ! block index ! LCOV_EXCL_LINE
3171 : ilo,ihi,jlo,jhi ! beginning and end of physical domain
3172 :
3173 : real (kind=dbl_kind) :: &
3174 0 : vmin, vmax
3175 :
3176 : logical (kind=log_kind) :: readm
3177 :
3178 : type (block) :: &
3179 : this_block ! block information for current block
3180 :
3181 : character(len=*), parameter :: subname = '(monthly_data)'
3182 :
3183 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
3184 :
3185 : !-------------------------------------------------------------------
3186 : ! monthly data
3187 : !
3188 : ! Assume that monthly data values are located in the middle of the
3189 : ! month.
3190 : !-------------------------------------------------------------------
3191 :
3192 0 : midmonth = 15 ! data is given on 15th of every month
3193 : ! midmonth = fix(p5 * real(daymo(mmonth))) ! exact middle
3194 :
3195 : ! Compute record numbers for surrounding months
3196 0 : maxrec = 12
3197 0 : ixm = mod(mmonth+maxrec-2,maxrec) + 1
3198 0 : ixp = mod(mmonth, maxrec) + 1
3199 0 : if (mday >= midmonth) ixm = -99 ! other two points will be used
3200 0 : if (mday < midmonth) ixp = -99
3201 :
3202 : ! Determine whether interpolation will use values 1:2 or 2:3
3203 : ! recslot = 2 means we use values 1:2, with the current value (2)
3204 : ! in the second slot
3205 : ! recslot = 1 means we use values 2:3, with the current value (2)
3206 : ! in the first slot
3207 0 : recslot = 1 ! latter half of month
3208 0 : if (mday < midmonth) recslot = 2 ! first half of month
3209 :
3210 : ! Find interpolation coefficients
3211 0 : call interp_coeff_monthly (recslot)
3212 :
3213 : ! Read 2 monthly values
3214 0 : readm = .false.
3215 0 : if (istep==1 .or. (mday==midmonth .and. msec==0)) readm = .true.
3216 :
3217 : call read_clim_data (readm, 0, ixm, mmonth, ixp, &
3218 : flw_file, cldf_data, & ! LCOV_EXCL_LINE
3219 0 : field_loc_center, field_type_scalar)
3220 : call read_clim_data (readm, 0, ixm, mmonth, ixp, &
3221 : rain_file, fsnow_data, & ! LCOV_EXCL_LINE
3222 0 : field_loc_center, field_type_scalar)
3223 : call read_clim_data (readm, 0, ixm, mmonth, ixp, &
3224 : tair_file, Tair_data, & ! LCOV_EXCL_LINE
3225 0 : field_loc_center, field_type_scalar)
3226 : call read_clim_data (readm, 0, ixm, mmonth, ixp, &
3227 : humid_file, Qa_data, & ! LCOV_EXCL_LINE
3228 0 : field_loc_center, field_type_scalar)
3229 : call read_clim_data (readm, 0, ixm, mmonth, ixp, &
3230 : wind_file, wind_data, & ! LCOV_EXCL_LINE
3231 0 : field_loc_center, field_type_scalar)
3232 : call read_clim_data (readm, 0, ixm, mmonth, ixp, &
3233 : strax_file, strax_data, & ! LCOV_EXCL_LINE
3234 0 : field_loc_center, field_type_vector)
3235 : call read_clim_data (readm, 0, ixm, mmonth, ixp, &
3236 : stray_file, stray_data, & ! LCOV_EXCL_LINE
3237 0 : field_loc_center, field_type_vector)
3238 :
3239 0 : call interpolate_data (cldf_data, cldf)
3240 0 : call interpolate_data (fsnow_data, fsnow) ! units mm/s = kg/m^2/s
3241 0 : call interpolate_data (Tair_data, Tair)
3242 0 : call interpolate_data (Qa_data, Qa)
3243 0 : call interpolate_data (wind_data, wind)
3244 0 : call interpolate_data (strax_data, strax)
3245 0 : call interpolate_data (stray_data, stray)
3246 :
3247 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,ilo,ihi,jlo,jhi,this_block)
3248 0 : do iblk = 1, nblocks
3249 : call Qa_fixLY(nx_block, ny_block, &
3250 : Tair (:,:,iblk), & ! LCOV_EXCL_LINE
3251 0 : Qa (:,:,iblk))
3252 :
3253 0 : do j = 1, ny_block
3254 0 : do i = 1, nx_block
3255 0 : Qa (i,j,iblk) = Qa (i,j,iblk) * hm(i,j,iblk)
3256 0 : Tair (i,j,iblk) = Tair (i,j,iblk) * hm(i,j,iblk)
3257 0 : wind (i,j,iblk) = wind (i,j,iblk) * hm(i,j,iblk)
3258 0 : strax(i,j,iblk) = strax(i,j,iblk) * hm(i,j,iblk)
3259 0 : stray(i,j,iblk) = stray(i,j,iblk) * hm(i,j,iblk)
3260 : enddo
3261 : enddo
3262 :
3263 : ! AOMIP
3264 0 : this_block = get_block(blocks_ice(iblk),iblk)
3265 0 : ilo = this_block%ilo
3266 0 : ihi = this_block%ihi
3267 0 : jlo = this_block%jlo
3268 0 : jhi = this_block%jhi
3269 :
3270 : call compute_shortwave(nx_block, ny_block, &
3271 : ilo, ihi, jlo, jhi, & ! LCOV_EXCL_LINE
3272 : TLON (:,:,iblk), & ! LCOV_EXCL_LINE
3273 : TLAT (:,:,iblk), & ! LCOV_EXCL_LINE
3274 : hm (:,:,iblk), & ! LCOV_EXCL_LINE
3275 : Qa (:,:,iblk), & ! LCOV_EXCL_LINE
3276 : cldf (:,:,iblk), & ! LCOV_EXCL_LINE
3277 0 : fsw (:,:,iblk))
3278 :
3279 : enddo ! iblk
3280 : !$OMP END PARALLEL DO
3281 :
3282 0 : if (debug_forcing) then
3283 0 : if (my_task == master_task) write (nu_diag,*) 'LY_bulk_data'
3284 0 : vmin = global_minval(fsw,distrb_info,tmask)
3285 0 : vmax = global_maxval(fsw,distrb_info,tmask)
3286 0 : if (my_task.eq.master_task) &
3287 0 : write (nu_diag,*) 'fsw',vmin,vmax
3288 0 : vmin = global_minval(cldf,distrb_info,tmask)
3289 0 : vmax = global_maxval(cldf,distrb_info,tmask)
3290 0 : if (my_task.eq.master_task) &
3291 0 : write (nu_diag,*) 'cldf',vmin,vmax
3292 0 : vmin =global_minval(fsnow,distrb_info,tmask)
3293 0 : vmax =global_maxval(fsnow,distrb_info,tmask)
3294 0 : if (my_task.eq.master_task) &
3295 0 : write (nu_diag,*) 'fsnow',vmin,vmax
3296 0 : vmin = global_minval(Tair,distrb_info,tmask)
3297 0 : vmax = global_maxval(Tair,distrb_info,tmask)
3298 0 : if (my_task.eq.master_task) &
3299 0 : write (nu_diag,*) 'Tair',vmin,vmax
3300 0 : vmin = global_minval(wind,distrb_info,umask)
3301 0 : vmax = global_maxval(wind,distrb_info,umask)
3302 0 : if (my_task.eq.master_task) &
3303 0 : write (nu_diag,*) 'wind',vmin,vmax
3304 0 : vmin = global_minval(strax,distrb_info,umask)
3305 0 : vmax = global_maxval(strax,distrb_info,umask)
3306 0 : if (my_task.eq.master_task) &
3307 0 : write (nu_diag,*) 'strax',vmin,vmax
3308 0 : vmin = global_minval(stray,distrb_info,umask)
3309 0 : vmax = global_maxval(stray,distrb_info,umask)
3310 0 : if (my_task.eq.master_task) &
3311 0 : write (nu_diag,*) 'stray',vmin,vmax
3312 0 : vmin = global_minval(Qa,distrb_info,tmask)
3313 0 : vmax = global_maxval(Qa,distrb_info,tmask)
3314 0 : if (my_task.eq.master_task) &
3315 0 : write (nu_diag,*) 'Qa',vmin,vmax
3316 :
3317 : endif ! debug_forcing
3318 :
3319 0 : end subroutine monthly_data
3320 :
3321 : !=======================================================================
3322 : ! Oned atmospheric data
3323 : !=======================================================================
3324 :
3325 0 : subroutine oned_data
3326 :
3327 : use ice_flux, only: uatm, vatm, Tair, fsw, fsnow, Qa, rhoa, frain
3328 :
3329 : ! local parameters
3330 :
3331 : character (char_len_long) :: &
3332 : met_file, & ! netcdf filename ! LCOV_EXCL_LINE
3333 : fieldname ! field name in netcdf file
3334 :
3335 : integer (kind=int_kind) :: &
3336 : fid ! file id for netCDF file
3337 :
3338 : real (kind=dbl_kind):: &
3339 0 : work ! temporary variable
3340 :
3341 : logical (kind=log_kind) :: diag
3342 :
3343 : integer (kind=int_kind) :: &
3344 : status ! status flag
3345 :
3346 : real (kind=dbl_kind) :: & ! used to determine specific humidity
3347 : Temp , & ! air temperature (K) ! LCOV_EXCL_LINE
3348 : rh , & ! relative humidity (%) ! LCOV_EXCL_LINE
3349 : Psat , & ! saturation vapour pressure (hPa) ! LCOV_EXCL_LINE
3350 0 : ws ! saturation mixing ratio
3351 :
3352 : real (kind=dbl_kind), parameter :: & ! coefficients for Hyland-Wexler Qa
3353 : ps1 = 0.58002206e4_dbl_kind, & ! (K) ! LCOV_EXCL_LINE
3354 : ps2 = 1.3914993_dbl_kind, & ! ! LCOV_EXCL_LINE
3355 : ps3 = 0.48640239e-1_dbl_kind, & ! (K^-1) ! LCOV_EXCL_LINE
3356 : ps4 = 0.41764768e-4_dbl_kind, & ! (K^-2) ! LCOV_EXCL_LINE
3357 : ps5 = 0.14452093e-7_dbl_kind, & ! (K^-3) ! LCOV_EXCL_LINE
3358 : ps6 = 6.5459673_dbl_kind, & ! ! LCOV_EXCL_LINE
3359 : ws1 = 621.97_dbl_kind, & ! for saturation mixing ratio ! LCOV_EXCL_LINE
3360 : Pair = 1020._dbl_kind ! Sea level pressure (hPa)
3361 :
3362 : character(len=*), parameter :: subname = '(oned_data)'
3363 :
3364 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
3365 :
3366 0 : diag = .false. ! write diagnostic information
3367 :
3368 0 : if (trim(atm_data_format) == 'nc') then ! read nc file
3369 :
3370 : ! hourly data beginning Jan 1, 1989, 01:00
3371 : ! HARDWIRED for dt = 1 hour!
3372 0 : met_file = uwind_file
3373 0 : call ice_open_nc(met_file,fid)
3374 :
3375 0 : fieldname='Uatm'
3376 0 : call ice_read_nc(fid,istep1,fieldname,work,diag)
3377 0 : uatm(:,:,:) = work
3378 :
3379 0 : fieldname='Vatm'
3380 0 : call ice_read_nc(fid,istep1,fieldname,work,diag)
3381 0 : vatm(:,:,:) = work
3382 :
3383 0 : fieldname='Tair'
3384 0 : call ice_read_nc(fid,istep1,fieldname,work,diag)
3385 0 : Temp = work
3386 0 : Tair(:,:,:) = Temp
3387 :
3388 0 : call ice_close_nc(fid)
3389 :
3390 : ! hourly solar data beginning Jan 1, 1989, 01:00
3391 0 : met_file = fsw_file
3392 0 : call ice_open_nc(met_file,fid)
3393 :
3394 0 : fieldname='fsw'
3395 0 : call ice_read_nc(fid,istep1,fieldname,work,diag)
3396 0 : fsw(:,:,:) = work
3397 :
3398 0 : call ice_close_nc(fid)
3399 :
3400 : ! hourly interpolated monthly data beginning Jan 1, 1989, 01:00
3401 0 : met_file = humid_file
3402 0 : call ice_open_nc(met_file,fid)
3403 :
3404 0 : fieldname='rh'
3405 0 : call ice_read_nc(fid,istep1,fieldname,work,diag)
3406 0 : rh = work
3407 :
3408 0 : fieldname='fsnow'
3409 0 : call ice_read_nc(fid,istep1,fieldname,work,diag)
3410 0 : fsnow(:,:,:) = work
3411 :
3412 0 : call ice_close_nc(fid)
3413 :
3414 : !-------------------------------------------------------------------
3415 : ! Find specific humidity using Hyland-Wexler formulation
3416 : ! Hyland, R.W. and A. Wexler, Formulations for the Thermodynamic
3417 : ! Properties of the saturated phases of H20 from 173.15K to 473.15K,
3418 : ! ASHRAE Trans, 89(2A), 500-519, 1983
3419 : !-------------------------------------------------------------------
3420 :
3421 : Psat = exp(-ps1/Temp + ps2 - ps3*Temp + ps4*Temp**2 - ps5 * Temp**3 &
3422 0 : + ps6 * log(Temp))*p01 ! saturation vapour pressure
3423 0 : ws = ws1 * Psat/(Pair - Psat) ! saturation mixing ratio
3424 0 : Qa(:,:,:) = rh * ws * p01/(c1 + rh * ws * p01) * p001
3425 : ! specific humidity (kg/kg)
3426 : endif ! atm_data_format
3427 :
3428 : ! flw calculated in prepare_forcing
3429 0 : rhoa (:,:,:) = 1.3_dbl_kind ! air density (kg/m^3)
3430 0 : cldf (:,:,:) = p25 ! cloud fraction
3431 0 : frain(:,:,:) = c0 ! this is available in hourlymet_rh file
3432 :
3433 0 : end subroutine oned_data
3434 :
3435 : !=======================================================================
3436 :
3437 0 : subroutine oned_files
3438 :
3439 : character(len=*), parameter :: subname = '(oned_files)'
3440 :
3441 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
3442 :
3443 : fsw_file = &
3444 0 : trim(atm_data_dir)//'/hourlysolar_brw1989_5yr.nc'
3445 :
3446 : rain_file = &
3447 0 : trim(atm_data_dir)//'/hourlymet_rh_5yr.nc'
3448 :
3449 : uwind_file = &
3450 0 : trim(atm_data_dir)//'/hourlymet_brw1989_5yr.nc'
3451 :
3452 : vwind_file = &
3453 0 : trim(atm_data_dir)//'/hourlymet_brw1989_5yr.nc'
3454 :
3455 : tair_file = &
3456 0 : trim(atm_data_dir)//'/hourlymet_brw1989_5yr.nc'
3457 :
3458 : humid_file = &
3459 0 : trim(atm_data_dir)//'/hourlymet_rh_5yr.nc'
3460 :
3461 0 : if (my_task == master_task) then
3462 0 : write (nu_diag,*) ' '
3463 0 : write (nu_diag,*) 'Atmospheric data files:'
3464 0 : write (nu_diag,*) trim(fsw_file)
3465 0 : write (nu_diag,*) trim(rain_file)
3466 0 : write (nu_diag,*) trim(uwind_file)
3467 0 : write (nu_diag,*) trim(vwind_file)
3468 0 : write (nu_diag,*) trim(tair_file)
3469 0 : write (nu_diag,*) trim(humid_file)
3470 : endif ! master_task
3471 :
3472 0 : end subroutine oned_files
3473 :
3474 : !=======================================================================
3475 : ! Climatological ocean forcing
3476 : !=======================================================================
3477 :
3478 0 : subroutine ocn_data_clim (dt)
3479 :
3480 : ! Interpolate monthly sss, sst data to timestep.
3481 : ! Restore prognostic sst to data.
3482 : ! Interpolate fields from U grid to T grid if necessary.
3483 :
3484 : ! author: Elizabeth C. Hunke and William H. Lipscomb, LANL
3485 :
3486 : use ice_domain, only: nblocks
3487 : use ice_flux, only: sss, sst
3488 :
3489 : real (kind=dbl_kind), intent(in) :: &
3490 : dt ! time step
3491 :
3492 : ! local variables
3493 :
3494 : integer (kind=int_kind) :: &
3495 : i, j, iblk , & ! horizontal indices ! LCOV_EXCL_LINE
3496 : ixm,ixp , & ! record numbers for neighboring months ! LCOV_EXCL_LINE
3497 : maxrec , & ! maximum record number ! LCOV_EXCL_LINE
3498 : recslot , & ! spline slot for current record ! LCOV_EXCL_LINE
3499 : midmonth ! middle day of month
3500 :
3501 : real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: &
3502 0 : sstdat ! data value toward which SST is restored
3503 :
3504 : logical (kind=log_kind) :: readm
3505 :
3506 : character(len=*), parameter :: subname = '(ocn_data_clim)'
3507 :
3508 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
3509 :
3510 0 : if (my_task == master_task .and. istep == 1) then
3511 0 : if (trim(ocn_data_type)=='clim') then
3512 0 : write (nu_diag,*) ' '
3513 0 : write (nu_diag,*) 'SSS data interpolated to timestep:'
3514 0 : write (nu_diag,*) trim(sss_file)
3515 0 : write (nu_diag,*) ' '
3516 0 : write (nu_diag,*) 'SST data interpolated to timestep:'
3517 0 : write (nu_diag,*) trim(sst_file)
3518 0 : if (restore_ocn) write (nu_diag,*) &
3519 0 : 'SST restoring timescale (days) =', trestore
3520 : endif
3521 : endif ! my_task, istep
3522 :
3523 : !-------------------------------------------------------------------
3524 : ! monthly data
3525 : !
3526 : ! Assume that monthly data values are located in the middle of the
3527 : ! month.
3528 : !-------------------------------------------------------------------
3529 :
3530 0 : if (trim(ocn_data_type)=='clim') then
3531 :
3532 0 : midmonth = 15 ! data is given on 15th of every month
3533 : !!! midmonth = fix(p5 * real(daymo(mmonth))) ! exact middle
3534 :
3535 : ! Compute record numbers for surrounding months
3536 0 : maxrec = 12
3537 0 : ixm = mod(mmonth+maxrec-2,maxrec) + 1
3538 0 : ixp = mod(mmonth, maxrec) + 1
3539 0 : if (mday >= midmonth) ixm = -99 ! other two points will be used
3540 0 : if (mday < midmonth) ixp = -99
3541 :
3542 : ! Determine whether interpolation will use values 1:2 or 2:3
3543 : ! recslot = 2 means we use values 1:2, with the current value (2)
3544 : ! in the second slot
3545 : ! recslot = 1 means we use values 2:3, with the current value (2)
3546 : ! in the first slot
3547 0 : recslot = 1 ! latter half of month
3548 0 : if (mday < midmonth) recslot = 2 ! first half of month
3549 :
3550 : ! Find interpolation coefficients
3551 0 : call interp_coeff_monthly (recslot)
3552 :
3553 0 : readm = .false.
3554 0 : if (istep==1 .or. (mday==midmonth .and. msec==0)) readm = .true.
3555 :
3556 : !-------------------------------------------------------------------
3557 : ! Read two monthly SSS values and interpolate.
3558 : ! Note: SSS is restored instantaneously to data.
3559 : !-------------------------------------------------------------------
3560 :
3561 : call read_clim_data (readm, 0, ixm, mmonth, ixp, &
3562 : sss_file, sss_data, & ! LCOV_EXCL_LINE
3563 0 : field_loc_center, field_type_scalar)
3564 0 : call interpolate_data (sss_data, sss)
3565 :
3566 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j)
3567 0 : do iblk = 1, nblocks
3568 0 : do j = 1, ny_block
3569 0 : do i = 1, nx_block
3570 0 : sss(i,j,iblk) = max(sss(i,j,iblk), c0)
3571 : enddo
3572 : enddo
3573 : enddo
3574 : !$OMP END PARALLEL DO
3575 :
3576 0 : call ocn_freezing_temperature
3577 : endif
3578 :
3579 : !-------------------------------------------------------------------
3580 : ! Read two monthly SST values and interpolate.
3581 : ! Restore toward interpolated value.
3582 : !-------------------------------------------------------------------
3583 :
3584 0 : if (trim(ocn_data_type)=='clim') then
3585 : call read_clim_data (readm, 0, ixm, mmonth, ixp, &
3586 : sst_file, sst_data, & ! LCOV_EXCL_LINE
3587 0 : field_loc_center, field_type_scalar)
3588 0 : call interpolate_data (sst_data, sstdat)
3589 :
3590 0 : if (restore_ocn) then
3591 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j)
3592 0 : do iblk = 1, nblocks
3593 0 : do j = 1, ny_block
3594 0 : do i = 1, nx_block
3595 : sst(i,j,iblk) = sst(i,j,iblk) &
3596 0 : + (sstdat(i,j,iblk)-sst(i,j,iblk))*dt/trest
3597 : enddo
3598 : enddo
3599 : enddo
3600 : !$OMP END PARALLEL DO
3601 : endif
3602 : endif
3603 :
3604 0 : end subroutine ocn_data_clim
3605 :
3606 : !=======================================================================
3607 : ! NCAR CESM M-configuration (AIO) ocean forcing
3608 : !=======================================================================
3609 :
3610 0 : subroutine ocn_data_ncar_init
3611 :
3612 : ! Reads NCAR pop ocean forcing data set 'pop_frc_gx1v3_010815.nc'
3613 : !
3614 : ! List of ocean forcing fields: Note that order is important!
3615 : ! (order is determined by field list in vname).
3616 : !
3617 : ! For ocean mixed layer-----------------------------units
3618 : !
3619 : ! 1 sst------temperature---------------------------(C)
3620 : ! 2 sss------salinity------------------------------(ppt)
3621 : ! 3 hbl------depth---------------------------------(m)
3622 : ! 4 u--------surface u current---------------------(m/s)
3623 : ! 5 v--------surface v current---------------------(m/s)
3624 : ! 6 dhdx-----surface tilt x direction--------------(m/m)
3625 : ! 7 dhdy-----surface tilt y direction--------------(m/m)
3626 : ! 8 qdp------ocean sub-mixed layer heat flux-------(W/m2)
3627 : !
3628 : ! Fields 4, 5, 6, 7 are on the U-grid; 1, 2, 3, and 8 are
3629 : ! on the T-grid.
3630 :
3631 : ! authors: Bruce Briegleb, NCAR
3632 : ! Elizabeth Hunke, LANL
3633 :
3634 : use ice_blocks, only: nx_block, ny_block
3635 : use ice_domain_size, only: max_blocks
3636 : #ifdef USE_NETCDF
3637 : use netcdf
3638 : #endif
3639 :
3640 : integer (kind=int_kind) :: &
3641 : n , & ! field index ! LCOV_EXCL_LINE
3642 : m , & ! month index ! LCOV_EXCL_LINE
3643 : nrec, & ! record number for direct access ! LCOV_EXCL_LINE
3644 : nbits
3645 :
3646 : character(char_len) :: &
3647 : vname(nfld) ! variable names to search for in file
3648 : data vname / &
3649 : 'T', 'S', 'hblt', 'U', 'V', & ! LCOV_EXCL_LINE
3650 : 'dhdx', 'dhdy', 'qdp' /
3651 :
3652 : integer (kind=int_kind) :: &
3653 : status , & ! status flag ! LCOV_EXCL_LINE
3654 : fid , & ! file id ! LCOV_EXCL_LINE
3655 : dimid , & ! dimension id ! LCOV_EXCL_LINE
3656 : nlat , & ! number of longitudes of data ! LCOV_EXCL_LINE
3657 : nlon ! number of latitudes of data
3658 :
3659 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
3660 0 : work1
3661 :
3662 : character(len=*), parameter :: subname = '(ocn_data_ncar_init)'
3663 :
3664 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
3665 :
3666 0 : if (my_task == master_task) then
3667 :
3668 0 : write (nu_diag,*) 'WARNING: evp_prep calculates surface tilt'
3669 0 : write (nu_diag,*) 'WARNING: stress from geostrophic currents,'
3670 0 : write (nu_diag,*) 'WARNING: not data from ocean forcing file.'
3671 0 : write (nu_diag,*) 'WARNING: Alter ice_dyn_evp.F90 if desired.'
3672 :
3673 0 : if (restore_ocn) write (nu_diag,*) &
3674 0 : 'SST restoring timescale = ',trestore,' days'
3675 :
3676 0 : sst_file = trim(ocn_data_dir)//'/'//trim(oceanmixed_file) ! not just sst
3677 :
3678 : !---------------------------------------------------------------
3679 : ! Read in ocean forcing data from an existing file
3680 : !---------------------------------------------------------------
3681 0 : write (nu_diag,*) 'ocean mixed layer forcing data file = ', &
3682 0 : trim(sst_file)
3683 :
3684 : endif ! master_task
3685 :
3686 0 : if (trim(ocn_data_format) == 'nc') then
3687 : #ifdef USE_NETCDF
3688 0 : if (my_task == master_task) then
3689 0 : call ice_open_nc(sst_file, fid)
3690 :
3691 : ! status = nf90_inq_dimid(fid,'nlon',dimid)
3692 0 : status = nf90_inq_dimid(fid,'ni',dimid)
3693 0 : status = nf90_inquire_dimension(fid,dimid,len=nlon)
3694 :
3695 : ! status = nf90_inq_dimid(fid,'nlat',dimid)
3696 0 : status = nf90_inq_dimid(fid,'nj',dimid)
3697 0 : status = nf90_inquire_dimension(fid,dimid,len=nlat)
3698 :
3699 0 : if( nlon .ne. nx_global ) then
3700 : call abort_ice (error_message=subname//'ice: ocn frc file nlon ne nx_global', &
3701 0 : file=__FILE__, line=__LINE__)
3702 : endif
3703 0 : if( nlat .ne. ny_global ) then
3704 : call abort_ice (error_message=subname//'ice: ocn frc file nlat ne ny_global', &
3705 0 : file=__FILE__, line=__LINE__)
3706 : endif
3707 :
3708 : endif ! master_task
3709 :
3710 : ! Read in ocean forcing data for all 12 months
3711 0 : do n=1,nfld
3712 0 : do m=1,12
3713 :
3714 : ! Note: netCDF does single to double conversion if necessary
3715 : ! if (n >= 4 .and. n <= 7) then
3716 : ! call ice_read_nc(fid, m, vname(n), work1, debug_forcing, &
3717 : ! field_loc_NEcorner, field_type_vector)
3718 : ! else
3719 0 : call ice_read_nc(fid, m, vname(n), work1, debug_forcing, &
3720 0 : field_loc_center, field_type_scalar)
3721 : ! endif
3722 :
3723 0 : ocn_frc_m(:,:,:,n,m) = work1(:,:,:)
3724 :
3725 : enddo ! month loop
3726 : enddo ! field loop
3727 :
3728 0 : if (my_task == master_task) call ice_close_nc(fid)
3729 : #else
3730 : call abort_ice(subname//'ERROR: USE_NETCDF cpp not defined for '//trim(sst_file), &
3731 : file=__FILE__, line=__LINE__)
3732 : #endif
3733 :
3734 : else ! binary format
3735 :
3736 0 : nbits = 64
3737 0 : call ice_open (nu_forcing, sst_file, nbits)
3738 :
3739 0 : nrec = 0
3740 0 : do n=1,nfld
3741 0 : do m=1,12
3742 0 : nrec = nrec + 1
3743 0 : if (n >= 4 .and. n <= 7) then
3744 : call ice_read (nu_forcing, nrec, work1, 'rda8', debug_forcing, &
3745 0 : field_loc_NEcorner, field_type_vector)
3746 : else
3747 : call ice_read (nu_forcing, nrec, work1, 'rda8', debug_forcing, &
3748 0 : field_loc_center, field_type_scalar)
3749 : endif
3750 0 : ocn_frc_m(:,:,:,n,m) = work1(:,:,:)
3751 : enddo ! month loop
3752 : enddo ! field loop
3753 0 : close (nu_forcing)
3754 :
3755 : endif
3756 :
3757 : !echmod - currents cause Fram outflow to be too large
3758 : ! ocn_frc_m(:,:,:,4,:) = c0
3759 : ! ocn_frc_m(:,:,:,5,:) = c0
3760 : !echmod
3761 :
3762 0 : end subroutine ocn_data_ncar_init
3763 :
3764 : !=======================================================================
3765 :
3766 0 : subroutine ocn_data_ncar_init_3D
3767 :
3768 : ! Reads NCAR pop ocean forcing data set 'oceanmixed_ice_depth.nc'
3769 : !
3770 : ! List of ocean forcing fields: Note that order is important!
3771 : ! (order is determined by field list in vname).
3772 : !
3773 : ! For ocean mixed layer-----------------------------units
3774 : !
3775 : ! 1 sst------temperature---------------------------(C)
3776 : ! 2 sss------salinity------------------------------(ppt)
3777 : ! 3 hbl------depth---------------------------------(m)
3778 : ! 4 u--------surface u current---------------------(m/s)
3779 : ! 5 v--------surface v current---------------------(m/s)
3780 : ! 6 dhdx-----surface tilt x direction--------------(m/m)
3781 : ! 7 dhdy-----surface tilt y direction--------------(m/m)
3782 : ! 8 qdp------ocean sub-mixed layer heat flux-------(W/m2)
3783 : !
3784 : ! All fields are on the T-grid.
3785 : !
3786 : ! authors: Bruce Briegleb, NCAR
3787 : ! Elizabeth Hunke, LANL
3788 :
3789 : use ice_blocks, only: nx_block, ny_block
3790 : use ice_domain_size, only: max_blocks
3791 : use ice_grid, only: grid_average_X2Y, ANGLET
3792 : use ice_read_write, only: ice_read_nc_uv
3793 : #ifdef USE_NETCDF
3794 : use netcdf
3795 : #endif
3796 :
3797 : #ifdef USE_NETCDF
3798 : integer (kind=int_kind) :: &
3799 : n , & ! field index ! LCOV_EXCL_LINE
3800 : m , & ! month index ! LCOV_EXCL_LINE
3801 : nzlev ! z level of currents
3802 :
3803 : character(char_len) :: &
3804 : vname(nfld) ! variable names to search for in file
3805 : data vname / &
3806 : 'T', 'S', 'hblt', 'U', 'V', & ! LCOV_EXCL_LINE
3807 : 'dhdx', 'dhdy', 'qdp' /
3808 :
3809 : integer (kind=int_kind) :: &
3810 : fid , & ! file id ! LCOV_EXCL_LINE
3811 : dimid ! dimension id
3812 :
3813 : integer (kind=int_kind) :: &
3814 : status , & ! status flag ! LCOV_EXCL_LINE
3815 : nlat , & ! number of longitudes of data ! LCOV_EXCL_LINE
3816 : nlon ! number of latitudes of data
3817 :
3818 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
3819 0 : work1, work2
3820 : #endif
3821 :
3822 : character(len=*), parameter :: subname = '(ocn_data_ncar_init_3D)'
3823 :
3824 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
3825 :
3826 0 : if (my_task == master_task) then
3827 :
3828 0 : write (nu_diag,*) 'WARNING: evp_prep calculates surface tilt'
3829 0 : write (nu_diag,*) 'WARNING: stress from geostrophic currents,'
3830 0 : write (nu_diag,*) 'WARNING: not data from ocean forcing file.'
3831 0 : write (nu_diag,*) 'WARNING: Alter ice_dyn_evp.F if desired.'
3832 :
3833 0 : if (restore_ocn) write (nu_diag,*) &
3834 0 : 'SST restoring timescale = ',trestore,' days'
3835 :
3836 0 : sst_file = trim(ocn_data_dir)//'/'//trim(oceanmixed_file) ! not just sst
3837 :
3838 : !---------------------------------------------------------------
3839 : ! Read in ocean forcing data from an existing file
3840 : !---------------------------------------------------------------
3841 0 : write (nu_diag,*) 'ocean mixed layer forcing data file = ', &
3842 0 : trim(sst_file)
3843 0 : write (nu_diag,*)
3844 :
3845 : endif ! master_task
3846 :
3847 0 : if (trim(ocn_data_format) == 'nc') then
3848 : #ifdef USE_NETCDF
3849 0 : if (my_task == master_task) then
3850 0 : call ice_open_nc(sst_file, fid)
3851 :
3852 : ! status = nf90_inq_dimid(fid,'nlon',dimid)
3853 0 : status = nf90_inq_dimid(fid,'ni',dimid)
3854 0 : status = nf90_inquire_dimension(fid,dimid,len=nlon)
3855 :
3856 : ! status = nf90_inq_dimid(fid,'nlat',dimid)
3857 0 : status = nf90_inq_dimid(fid,'nj',dimid)
3858 0 : status = nf90_inquire_dimension(fid,dimid,len=nlat)
3859 :
3860 0 : if( nlon .ne. nx_global ) then
3861 : call abort_ice (error_message=subname//'ice: ocn frc file nlon ne nx_global', &
3862 0 : file=__FILE__, line=__LINE__)
3863 : endif
3864 0 : if( nlat .ne. ny_global ) then
3865 : call abort_ice (error_message=subname//'ice: ocn frc file nlat ne ny_global', &
3866 0 : file=__FILE__, line=__LINE__)
3867 : endif
3868 :
3869 : endif ! master_task
3870 :
3871 : ! Read in ocean forcing data for all 12 months
3872 0 : do n=1,nfld
3873 0 : do m=1,12
3874 :
3875 : ! Note: netCDF does single to double conversion if necessary
3876 0 : if (n == 4 .or. n == 5) then ! 3D currents
3877 0 : nzlev = 1 ! surface currents
3878 0 : call ice_read_nc_uv(fid, m, nzlev, vname(n), work1, debug_forcing, &
3879 0 : field_loc_center, field_type_scalar)
3880 : else
3881 0 : call ice_read_nc(fid, m, vname(n), work1, debug_forcing, &
3882 0 : field_loc_center, field_type_scalar)
3883 : endif
3884 :
3885 : ! the land mask used in ocean_mixed_depth.nc does not
3886 : ! match our gx1v3 mask (hm)
3887 0 : where (work1(:,:,:) < -900.) work1(:,:,:) = c0
3888 :
3889 0 : ocn_frc_m(:,:,:,n,m) = work1(:,:,:)
3890 :
3891 : enddo ! month loop
3892 : enddo ! field loop
3893 :
3894 0 : if (my_task == master_task) call ice_close_nc(fid)
3895 :
3896 : ! Rotate vector quantities and shift to U-grid
3897 0 : do n=4,6,2
3898 0 : do m=1,12
3899 :
3900 0 : work1(:,:,:) = ocn_frc_m(:,:,:,n ,m)
3901 0 : work2(:,:,:) = ocn_frc_m(:,:,:,n+1,m)
3902 0 : ocn_frc_m(:,:,:,n ,m) = work1(:,:,:)*cos(ANGLET(:,:,:)) &
3903 0 : + work2(:,:,:)*sin(ANGLET(:,:,:))
3904 0 : ocn_frc_m(:,:,:,n+1,m) = work2(:,:,:)*cos(ANGLET(:,:,:)) &
3905 0 : - work1(:,:,:)*sin(ANGLET(:,:,:))
3906 :
3907 0 : work1(:,:,:) = ocn_frc_m(:,:,:,n ,m)
3908 0 : work2(:,:,:) = ocn_frc_m(:,:,:,n+1,m)
3909 0 : call grid_average_X2Y('A',work1,'T',ocn_frc_m(:,:,:,n ,m),'U')
3910 0 : call grid_average_X2Y('A',work2,'T',ocn_frc_m(:,:,:,n+1,m),'U')
3911 :
3912 : enddo ! month loop
3913 : enddo ! field loop
3914 :
3915 : #else
3916 : call abort_ice(subname//'ERROR: USE_NETCDF cpp not defined', &
3917 : file=__FILE__, line=__LINE__)
3918 : #endif
3919 :
3920 : else ! binary format
3921 :
3922 : call abort_ice (error_message=subname//'new ocean forcing is netcdf only', &
3923 0 : file=__FILE__, line=__LINE__)
3924 :
3925 : endif
3926 :
3927 0 : end subroutine ocn_data_ncar_init_3D
3928 :
3929 : !=======================================================================
3930 :
3931 0 : subroutine ocn_data_ncar(dt)
3932 :
3933 : ! Interpolate monthly ocean data to timestep.
3934 : ! Restore sst if desired. sst is updated with surface fluxes in ice_ocean.F.
3935 :
3936 : use ice_blocks, only: nx_block, ny_block
3937 : use ice_global_reductions, only: global_minval, global_maxval
3938 : use ice_domain, only: nblocks, distrb_info
3939 : use ice_domain_size, only: max_blocks
3940 : use ice_flux, only: sss, sst, Tf, uocn, vocn, ss_tltx, ss_tlty, &
3941 : qdp, hmix
3942 : use ice_restart_shared, only: restart
3943 : use ice_grid, only: hm, tmask, umask
3944 :
3945 : real (kind=dbl_kind), intent(in) :: &
3946 : dt ! time step
3947 :
3948 : integer (kind=int_kind) :: &
3949 : i, j, n, iblk , & ! LCOV_EXCL_LINE
3950 : ixm,ixp , & ! record numbers for neighboring months ! LCOV_EXCL_LINE
3951 : maxrec , & ! maximum record number ! LCOV_EXCL_LINE
3952 : recslot , & ! spline slot for current record ! LCOV_EXCL_LINE
3953 : midmonth ! middle day of month
3954 :
3955 : real (kind=dbl_kind) :: &
3956 0 : vmin, vmax
3957 :
3958 : real (kind=dbl_kind), dimension (nx_block,ny_block,max_blocks) :: &
3959 0 : work1
3960 :
3961 : character(len=*), parameter :: subname = '(ocn_data_ncar)'
3962 :
3963 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
3964 :
3965 : !-------------------------------------------------------------------
3966 : ! monthly data
3967 : !
3968 : ! Assume that monthly data values are located in the middle of the
3969 : ! month.
3970 : !-------------------------------------------------------------------
3971 :
3972 0 : midmonth = 15 ! data is given on 15th of every month
3973 : ! midmonth = fix(p5 * real(daymo(mmonth),kind=dbl_kind)) ! exact middle
3974 :
3975 : ! Compute record numbers for surrounding months
3976 0 : maxrec = 12
3977 0 : ixm = mod(mmonth+maxrec-2,maxrec) + 1
3978 0 : ixp = mod(mmonth, maxrec) + 1
3979 0 : if (mday >= midmonth) ixm = -99 ! other two points will be used
3980 0 : if (mday < midmonth) ixp = -99
3981 :
3982 : ! Determine whether interpolation will use values 1:2 or 2:3
3983 : ! recslot = 2 means we use values 1:2, with the current value (2)
3984 : ! in the second slot
3985 : ! recslot = 1 means we use values 2:3, with the current value (2)
3986 : ! in the first slot
3987 0 : recslot = 1 ! latter half of month
3988 0 : if (mday < midmonth) recslot = 2 ! first half of month
3989 :
3990 : ! Find interpolation coefficients
3991 0 : call interp_coeff_monthly (recslot)
3992 :
3993 0 : sst_data(:,:,:,:) = c0
3994 0 : do n = nfld, 1, -1
3995 0 : do iblk = 1, nblocks
3996 : ! use sst_data arrays as temporary work space until n=1
3997 0 : if (ixm /= -99) then ! first half of month
3998 0 : sst_data(:,:,1,iblk) = ocn_frc_m(:,:,iblk,n,ixm)
3999 0 : sst_data(:,:,2,iblk) = ocn_frc_m(:,:,iblk,n,mmonth)
4000 : else ! second half of month
4001 0 : sst_data(:,:,1,iblk) = ocn_frc_m(:,:,iblk,n,mmonth)
4002 0 : sst_data(:,:,2,iblk) = ocn_frc_m(:,:,iblk,n,ixp)
4003 : endif
4004 : enddo
4005 :
4006 0 : call interpolate_data (sst_data,work1)
4007 : ! masking by hm is necessary due to NaNs in the data file
4008 0 : do j = 1, ny_block
4009 0 : do i = 1, nx_block
4010 0 : if (n == 2) sss (i,j,:) = c0
4011 0 : if (n == 3) hmix (i,j,:) = c0
4012 0 : if (n == 4) uocn (i,j,:) = c0
4013 0 : if (n == 5) vocn (i,j,:) = c0
4014 0 : if (n == 6) ss_tltx(i,j,:) = c0
4015 0 : if (n == 7) ss_tlty(i,j,:) = c0
4016 0 : if (n == 8) qdp (i,j,:) = c0
4017 0 : do iblk = 1, nblocks
4018 0 : if (hm(i,j,iblk) == c1) then
4019 0 : if (n == 2) sss (i,j,iblk) = work1(i,j,iblk)
4020 0 : if (n == 3) hmix (i,j,iblk) = max(mixed_layer_depth_default,work1(i,j,iblk))
4021 0 : if (n == 4) uocn (i,j,iblk) = work1(i,j,iblk)
4022 0 : if (n == 5) vocn (i,j,iblk) = work1(i,j,iblk)
4023 0 : if (n == 6) ss_tltx(i,j,iblk) = work1(i,j,iblk)
4024 0 : if (n == 7) ss_tlty(i,j,iblk) = work1(i,j,iblk)
4025 0 : if (n == 8) qdp (i,j,iblk) = work1(i,j,iblk)
4026 : endif
4027 : enddo
4028 : enddo
4029 : enddo
4030 : enddo
4031 :
4032 0 : do j = 1, ny_block
4033 0 : do i = 1, nx_block
4034 0 : sss (i,j,:) = max (sss(i,j,:), c0)
4035 0 : hmix(i,j,:) = max(hmix(i,j,:), c0)
4036 : enddo
4037 : enddo
4038 :
4039 0 : call ocn_freezing_temperature
4040 :
4041 0 : if (restore_ocn) then
4042 0 : do j = 1, ny_block
4043 0 : do i = 1, nx_block
4044 0 : sst(i,j,:) = sst(i,j,:) + (work1(i,j,:)-sst(i,j,:))*dt/trest
4045 : enddo
4046 : enddo
4047 : ! else sst is only updated in ice_ocean.F
4048 : endif
4049 :
4050 : ! initialize sst properly on first step
4051 0 : if (istep1 <= 1 .and. .not. (restart)) then
4052 0 : call interpolate_data (sst_data,sst)
4053 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j)
4054 0 : do iblk = 1, nblocks
4055 0 : do j = 1, ny_block
4056 0 : do i = 1, nx_block
4057 0 : if (hm(i,j,iblk) == c1) then
4058 0 : sst(i,j,iblk) = max (sst(i,j,iblk), Tf(i,j,iblk))
4059 : else
4060 0 : sst(i,j,iblk) = c0
4061 : endif
4062 : enddo
4063 : enddo
4064 : enddo
4065 : !$OMP END PARALLEL DO
4066 : endif
4067 :
4068 0 : if (debug_forcing) then
4069 0 : if (my_task == master_task) &
4070 0 : write (nu_diag,*) 'ocn_data_ncar'
4071 0 : vmin = global_minval(Tf,distrb_info,tmask)
4072 0 : vmax = global_maxval(Tf,distrb_info,tmask)
4073 0 : if (my_task.eq.master_task) &
4074 0 : write (nu_diag,*) 'Tf',vmin,vmax
4075 0 : vmin = global_minval(sst,distrb_info,tmask)
4076 0 : vmax = global_maxval(sst,distrb_info,tmask)
4077 0 : if (my_task.eq.master_task) &
4078 0 : write (nu_diag,*) 'sst',vmin,vmax
4079 0 : vmin = global_minval(sss,distrb_info,tmask)
4080 0 : vmax = global_maxval(sss,distrb_info,tmask)
4081 0 : if (my_task.eq.master_task) &
4082 0 : write (nu_diag,*) 'sss',vmin,vmax
4083 0 : vmin = global_minval(hmix,distrb_info,tmask)
4084 0 : vmax = global_maxval(hmix,distrb_info,tmask)
4085 0 : if (my_task.eq.master_task) &
4086 0 : write (nu_diag,*) 'hmix',vmin,vmax
4087 0 : vmin = global_minval(uocn,distrb_info,umask)
4088 0 : vmax = global_maxval(uocn,distrb_info,umask)
4089 0 : if (my_task.eq.master_task) &
4090 0 : write (nu_diag,*) 'uocn',vmin,vmax
4091 0 : vmin = global_minval(vocn,distrb_info,umask)
4092 0 : vmax = global_maxval(vocn,distrb_info,umask)
4093 0 : if (my_task.eq.master_task) &
4094 0 : write (nu_diag,*) 'vocn',vmin,vmax
4095 0 : vmin = global_minval(ss_tltx,distrb_info,umask)
4096 0 : vmax = global_maxval(ss_tltx,distrb_info,umask)
4097 0 : if (my_task.eq.master_task) &
4098 0 : write (nu_diag,*) 'ss_tltx',vmin,vmax
4099 0 : vmin = global_minval(ss_tlty,distrb_info,umask)
4100 0 : vmax = global_maxval(ss_tlty,distrb_info,umask)
4101 0 : if (my_task.eq.master_task) &
4102 0 : write (nu_diag,*) 'ss_tlty',vmin,vmax
4103 0 : vmin = global_minval(qdp,distrb_info,tmask)
4104 0 : vmax = global_maxval(qdp,distrb_info,tmask)
4105 0 : if (my_task.eq.master_task) &
4106 0 : write (nu_diag,*) 'qdp',vmin,vmax
4107 : endif
4108 :
4109 0 : end subroutine ocn_data_ncar
4110 :
4111 : !=======================================================================
4112 : ! ocean data for oned configuration
4113 : ! Current (released) values are the same as the defaults (ice_flux.F90)
4114 :
4115 0 : subroutine ocn_data_oned
4116 :
4117 : use ice_flux, only: sss, sst, Tf, uocn, vocn, ss_tltx, ss_tlty, &
4118 : qdp, hmix, frzmlt
4119 :
4120 : character(len=*), parameter :: subname = '(ocn_data_oned)'
4121 :
4122 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
4123 :
4124 0 : sss (:,:,:) = 34.0_dbl_kind ! sea surface salinity (ppt)
4125 :
4126 0 : call ocn_freezing_temperature
4127 :
4128 0 : sst (:,:,:) = Tf(:,:,:) ! sea surface temp (C)
4129 0 : uocn (:,:,:) = c0 ! surface ocean currents (m/s)
4130 0 : vocn (:,:,:) = c0
4131 0 : ss_tltx(:,:,:) = c0 ! sea surface tilt (m/m)
4132 0 : ss_tlty(:,:,:) = c0
4133 0 : frzmlt (:,:,:) = c0 ! freezing/melting potential (W/m^2)
4134 0 : qdp (:,:,:) = c0 ! deep ocean heat flux (W/m^2)
4135 0 : hmix (:,:,:) = mixed_layer_depth_default ! ocean mixed layer depth
4136 :
4137 0 : end subroutine ocn_data_oned
4138 :
4139 : !=======================================================================
4140 :
4141 0 : subroutine ocn_data_hadgem(dt)
4142 :
4143 : ! Reads in HadGEM ocean forcing data as required from netCDF files
4144 : ! Current options (selected by ocn_data_type)
4145 : ! hadgem_sst: read in sst only
4146 : ! hadgem_sst_uvocn: read in sst plus uocn and vocn
4147 :
4148 : ! authors: Ann Keen, Met Office
4149 :
4150 : use ice_domain, only: nblocks
4151 : use ice_domain_size, only: max_blocks
4152 : use ice_flux, only: sst, uocn, vocn
4153 : use ice_grid, only: ANGLET
4154 :
4155 : real (kind=dbl_kind), intent(in) :: &
4156 : dt ! time step
4157 :
4158 : integer (kind=int_kind) :: &
4159 : i, j, iblk , & ! LCOV_EXCL_LINE
4160 : ixm,ixp , & ! record numbers for neighboring months ! LCOV_EXCL_LINE
4161 : maxrec , & ! maximum record number ! LCOV_EXCL_LINE
4162 : recslot , & ! spline slot for current record ! LCOV_EXCL_LINE
4163 : midmonth ! middle day of month
4164 :
4165 : real (kind=dbl_kind), dimension(nx_block,ny_block,max_blocks) :: &
4166 0 : sstdat ! data value toward which SST is restored
4167 :
4168 0 : real (kind=dbl_kind) :: workx, worky
4169 :
4170 : logical (kind=log_kind) :: readm
4171 :
4172 : character (char_len) :: &
4173 : fieldname ! field name in netcdf file
4174 :
4175 : character (char_len_long) :: &
4176 : filename ! name of netCDF file
4177 :
4178 : character(len=*), parameter :: subname = '(ocn_data_hadgem)'
4179 :
4180 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
4181 :
4182 : !-------------------------------------------------------------------
4183 : ! monthly data
4184 : !
4185 : ! Assume that monthly data values are located in the middle of the
4186 : ! month.
4187 : !-------------------------------------------------------------------
4188 :
4189 0 : midmonth = 15 ! data is given on 15th of every month
4190 : ! midmonth = fix(p5 * real(daymo(mmonth))) ! exact middle
4191 :
4192 : ! Compute record numbers for surrounding months
4193 0 : maxrec = 12
4194 0 : ixm = mod(mmonth+maxrec-2,maxrec) + 1
4195 0 : ixp = mod(mmonth, maxrec) + 1
4196 0 : if (mday >= midmonth) ixm = -99 ! other two points will be used
4197 0 : if (mday < midmonth) ixp = -99
4198 :
4199 : ! Determine whether interpolation will use values 1:2 or 2:3
4200 : ! recslot = 2 means we use values 1:2, with the current value (2)
4201 : ! in the second slot
4202 : ! recslot = 1 means we use values 2:3, with the current value (2)
4203 : ! in the first slot
4204 0 : recslot = 1 ! latter half of month
4205 0 : if (mday < midmonth) recslot = 2 ! first half of month
4206 :
4207 : ! Find interpolation coefficients
4208 0 : call interp_coeff_monthly (recslot)
4209 :
4210 : ! Read 2 monthly values
4211 0 : readm = .false.
4212 0 : if (istep==1 .or. (mday==midmonth .and. msec==0)) readm = .true.
4213 :
4214 0 : if (my_task == master_task .and. istep == 1) then
4215 0 : write (nu_diag,*) ' '
4216 0 : write (nu_diag,*) 'SST data interpolated to timestep:'
4217 0 : write (nu_diag,*) trim(ocn_data_dir)//'/MONTHLY/sst.1997.nc'
4218 0 : if (restore_ocn) write (nu_diag,*) &
4219 0 : 'SST restoring timescale (days) =', trestore
4220 0 : if (trim(ocn_data_type)=='hadgem_sst_uvocn') then
4221 0 : write (nu_diag,*) ' '
4222 0 : write (nu_diag,*) 'uocn and vocn interpolated to timestep:'
4223 0 : write (nu_diag,*) trim(ocn_data_dir)//'/MONTHLY/uocn.1997.nc'
4224 0 : write (nu_diag,*) trim(ocn_data_dir)//'/MONTHLY/vocn.1997.nc'
4225 : endif
4226 : endif ! my_task, istep
4227 :
4228 : ! -----------------------------------------------------------
4229 : ! SST
4230 : ! -----------------------------------------------------------
4231 0 : sst_file = trim(ocn_data_dir)//'/MONTHLY/sst.1997.nc'
4232 0 : fieldname='sst'
4233 : call read_data_nc (readm, 0, fyear, ixm, mmonth, ixp, &
4234 : maxrec, sst_file, fieldname, sst_data, & ! LCOV_EXCL_LINE
4235 0 : field_loc_center, field_type_scalar)
4236 :
4237 : ! Interpolate to current time step
4238 0 : call interpolate_data (sst_data, sstdat)
4239 :
4240 : ! Restore SSTs if required
4241 0 : if (restore_ocn) then
4242 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j)
4243 0 : do iblk = 1, nblocks
4244 0 : do j = 1, ny_block
4245 0 : do i = 1, nx_block
4246 : sst(i,j,iblk) = sst(i,j,iblk) &
4247 0 : + (sstdat(i,j,iblk)-sst(i,j,iblk))*dt/trest
4248 : enddo
4249 : enddo
4250 : enddo
4251 : !$OMP END PARALLEL DO
4252 : endif
4253 :
4254 : ! -----------------------------------------------------------
4255 : ! Ocean currents
4256 : ! --------------
4257 : ! Values read in are on T grid and oriented geographically, hence
4258 : ! vectors need to be rotated to model grid and then interpolated
4259 : ! to U grid.
4260 : ! Also need to be converted from cm s-1 (UM) to m s-1 (CICE)
4261 : ! -----------------------------------------------------------
4262 :
4263 0 : if (trim(ocn_data_type)=='hadgem_sst_uvocn') then
4264 :
4265 0 : filename = trim(ocn_data_dir)//'/MONTHLY/uocn.1997.nc'
4266 0 : fieldname='uocn'
4267 : call read_data_nc (readm, 0, fyear, ixm, mmonth, ixp, &
4268 : maxrec, filename, fieldname, uocn_data, & ! LCOV_EXCL_LINE
4269 0 : field_loc_center, field_type_vector)
4270 :
4271 : ! Interpolate to current time step
4272 0 : call interpolate_data (uocn_data, uocn)
4273 :
4274 0 : filename = trim(ocn_data_dir)//'/MONTHLY/vocn.1997.nc'
4275 0 : fieldname='vocn'
4276 : call read_data_nc (readm, 0, fyear, ixm, mmonth, ixp, &
4277 : maxrec, filename, fieldname, vocn_data, & ! LCOV_EXCL_LINE
4278 0 : field_loc_center, field_type_vector)
4279 :
4280 : ! Interpolate to current time step
4281 0 : call interpolate_data (vocn_data, vocn)
4282 :
4283 : !-----------------------------------------------------------------
4284 : ! Rotate zonal/meridional vectors to local coordinates,
4285 : ! and change units
4286 : !-----------------------------------------------------------------
4287 :
4288 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j,workx,worky)
4289 0 : do iblk = 1, nblocks
4290 0 : do j = 1, ny_block
4291 0 : do i = 1, nx_block
4292 :
4293 0 : workx = uocn(i,j,iblk)
4294 0 : worky = vocn(i,j,iblk)
4295 : uocn(i,j,iblk) = workx*cos(ANGLET(i,j,iblk)) &
4296 0 : + worky*sin(ANGLET(i,j,iblk))
4297 : vocn(i,j,iblk) = worky*cos(ANGLET(i,j,iblk)) &
4298 0 : - workx*sin(ANGLET(i,j,iblk))
4299 :
4300 0 : uocn(i,j,iblk) = uocn(i,j,iblk) * cm_to_m
4301 0 : vocn(i,j,iblk) = vocn(i,j,iblk) * cm_to_m
4302 :
4303 : enddo ! i
4304 : enddo ! j
4305 : enddo ! nblocks
4306 : !$OMP END PARALLEL DO
4307 :
4308 : !-----------------------------------------------------------------
4309 : ! Interpolate to U grid
4310 : !-----------------------------------------------------------------
4311 :
4312 : ! tcraig, this is now computed in dynamics for consistency
4313 :
4314 : endif ! ocn_data_type = hadgem_sst_uvocn
4315 :
4316 0 : end subroutine ocn_data_hadgem
4317 :
4318 : !=======================================================================
4319 :
4320 0 : subroutine ocn_data_hycom_init
4321 : ! Read SSS+SST from a HYCOM file converted to NetCDF format.
4322 : ! HYCOM binary2NetCDF: hcdata2ncdf2d (or hcdata2ncdf3z)
4323 : ! + rename/link file
4324 : use ice_blocks, only: nx_block, ny_block
4325 : use ice_domain, only: nblocks
4326 : use ice_flux, only: sss, sst, Tf
4327 :
4328 : integer (kind=int_kind) :: &
4329 : i, j, iblk , & ! horizontal indices ! LCOV_EXCL_LINE
4330 : fid ! file id for netCDF file
4331 :
4332 : character (char_len) :: &
4333 : fieldname ! field name in netcdf file
4334 :
4335 : character(len=*), parameter :: subname = '(ocn_data_hycom_init)'
4336 :
4337 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
4338 :
4339 0 : if (trim(ocn_data_type) == 'hycom') then
4340 0 : sss_file = trim(ocn_data_dir)//'ice.restart.surf.nc'
4341 :
4342 0 : if (my_task == master_task) then
4343 0 : write (nu_diag,*)' '
4344 0 : write (nu_diag,*)'Initial SSS file: ',trim(sss_file)
4345 : endif
4346 :
4347 0 : fieldname = 'sss'
4348 0 : call ice_open_nc (sss_file, fid)
4349 : call ice_read_nc (fid, 1 , fieldname, sss, debug_forcing, &
4350 0 : field_loc_center, field_type_scalar)
4351 0 : call ice_close_nc(fid)
4352 :
4353 0 : call ocn_freezing_temperature
4354 :
4355 0 : sst_file = trim(ocn_data_dir)//'ice.restart.surf.nc'
4356 :
4357 0 : if (my_task == master_task) then
4358 0 : write (nu_diag,*)' '
4359 0 : write (nu_diag,*)'Initial SST file: ',trim(sst_file)
4360 : endif
4361 :
4362 0 : fieldname = 'sst'
4363 0 : call ice_open_nc (sst_file, fid)
4364 : call ice_read_nc (fid, 1 , fieldname, sst, debug_forcing, &
4365 0 : field_loc_center, field_type_scalar)
4366 0 : call ice_close_nc(fid)
4367 :
4368 : ! Make sure sst is not less than freezing temperature Tf
4369 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j)
4370 0 : do iblk = 1, nblocks
4371 0 : do j = 1, ny_block
4372 0 : do i = 1, nx_block
4373 0 : sst(i,j,iblk) = max(sst(i,j,iblk),Tf(i,j,iblk))
4374 : enddo
4375 : enddo
4376 : enddo
4377 : !$OMP END PARALLEL DO
4378 : endif
4379 :
4380 0 : end subroutine ocn_data_hycom_init
4381 :
4382 : !=======================================================================
4383 :
4384 0 : subroutine hycom_atm_files
4385 :
4386 : use ice_broadcast, only: broadcast_array, broadcast_scalar
4387 :
4388 : integer (kind = int_kind) :: &
4389 : fid ! File id
4390 : character (char_len) :: &
4391 : varname ! variable name in netcdf file
4392 : character(len=*), parameter :: subname = '(hycom_atm_files)'
4393 :
4394 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
4395 :
4396 0 : fsw_file = trim(atm_data_dir)//'/forcing.shwflx.nc'
4397 0 : flw_file = trim(atm_data_dir)//'/forcing.radflx.nc'
4398 0 : rain_file = trim(atm_data_dir)//'/forcing.precip.nc'
4399 0 : uwind_file = trim(atm_data_dir)//'/forcing.wndewd.nc'
4400 0 : vwind_file = trim(atm_data_dir)//'/forcing.wndnwd.nc'
4401 0 : tair_file = trim(atm_data_dir)//'/forcing.airtmp.nc'
4402 0 : humid_file = trim(atm_data_dir)//'/forcing.vapmix.nc'
4403 :
4404 : ! Read time vector from "tair_file"
4405 0 : call ice_open_nc(tair_file, fid)
4406 0 : varname='MT'
4407 0 : call ice_get_ncvarsize(fid,varname,Njday_atm)
4408 0 : call broadcast_scalar(Njday_atm,master_task)
4409 0 : allocate(jday_atm(Njday_atm))
4410 0 : call ice_read_vec_nc(fid,Njday_atm, varname,jday_atm, .true.)
4411 0 : call ice_close_nc(fid)
4412 0 : call broadcast_array(jday_atm ,master_task)
4413 :
4414 : ! Write diag info
4415 0 : if (my_task == master_task) then
4416 0 : write (nu_diag,*) ' '
4417 0 : write (nu_diag,*) 'CICE: Atm. (hycomdate) Start = ',jday_atm(1)
4418 0 : write (nu_diag,*) 'CICE: Atm. (hycomdate) End = ',jday_atm(Njday_atm)
4419 0 : write (nu_diag,*) 'CICE: Total Atm timesteps = ',Njday_atm
4420 0 : write (nu_diag,*) 'CICE: Atmospheric forcing files:'
4421 0 : write (nu_diag,*) trim(fsw_file)
4422 0 : write (nu_diag,*) trim(flw_file)
4423 0 : write (nu_diag,*) trim(rain_file)
4424 0 : write (nu_diag,*) trim(uwind_file)
4425 0 : write (nu_diag,*) trim(vwind_file)
4426 0 : write (nu_diag,*) trim(tair_file)
4427 0 : write (nu_diag,*) trim(humid_file)
4428 : endif ! master_task
4429 :
4430 0 : end subroutine hycom_atm_files
4431 :
4432 : !=======================================================================
4433 :
4434 0 : subroutine hycom_atm_data
4435 :
4436 : use ice_flux, only: fsw, fsnow, Tair, uatm, vatm, Qa, flw
4437 : use ice_domain, only: nblocks
4438 :
4439 : integer (kind=int_kind) :: &
4440 : recnum ! record number
4441 :
4442 : real (kind=dbl_kind) :: &
4443 0 : hcdate ! current time in HYCOM jday units
4444 :
4445 : logical (kind=log_kind) :: read6
4446 :
4447 : character (char_len) :: &
4448 : fieldname ! field name in netcdf file
4449 :
4450 : integer (kind=int_kind) :: &
4451 : i, j, iblk ! horizontal indices
4452 :
4453 0 : real (kind=dbl_kind) :: Tffresh, secday
4454 :
4455 : character(len=*), parameter :: subname = '(hycom_atm_data)'
4456 :
4457 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
4458 :
4459 0 : call icepack_query_parameters(Tffresh_out=Tffresh)
4460 0 : call icepack_query_parameters(secday_out=secday)
4461 :
4462 : ! current time in HYCOM jday units (HYCOM ref year: 1900,12,31,000000)
4463 0 : hcdate = real(compute_days_between(1900,12,31,myear,mmonth,mday)) + msec/secday
4464 :
4465 : ! Init recnum try
4466 0 : recnum=min(max(oldrecnum,1),Njday_atm-1)
4467 :
4468 : ! Find correct time in ATM data ... assume cont. incr. time-axis
4469 0 : do while ( recnum<Njday_atm )
4470 0 : if ( hcdate>=jday_atm(recnum) .and. &
4471 0 : hcdate<=jday_atm(recnum+1) ) exit
4472 0 : if ( abs(hcdate-jday_atm(recnum))<p001 ) exit ! Accept within tolerance = 0.001 days
4473 0 : if ( abs(hcdate-jday_atm(recnum+1))<p001 ) exit ! Accept within tolerance = 0.001 days
4474 0 : recnum=recnum+1
4475 : enddo
4476 :
4477 : ! Last Atm date might be the same as last CICE date.
4478 0 : recnum=min(recnum,Njday_atm-1)
4479 :
4480 : ! Check if current time do not exceed last forcing time
4481 : ! + check forcing is available before (or at) current forcing time
4482 0 : if ( hcdate>jday_atm(recnum+1)+p001 .or. hcdate<jday_atm(recnum)-p001) then
4483 : write (nu_diag,*) &
4484 0 : 'ERROR: CICE: Atm forcing not available at hcdate =',hcdate
4485 : write (nu_diag,*) &
4486 0 : 'ERROR: CICE: myear, yday ,msec = ',myear, yday, msec
4487 0 : call abort_ice ('ERROR: CICE stopped')
4488 : endif
4489 :
4490 : ! Compute interpolation coefficients
4491 0 : call interp_coeff2 (hcdate, jday_atm(recnum), jday_atm(recnum+1) )
4492 :
4493 :
4494 : ! Read
4495 0 : read6 = .false.
4496 0 : if (istep==1 .or. oldrecnum /= recnum) read6 = .true.
4497 :
4498 0 : if (trim(atm_data_format) == 'nc') then
4499 :
4500 0 : if (read6 .and. my_task == master_task) write(nu_diag,*) &
4501 0 : 'CICE: Atm. read: = ',jday_atm(recnum), jday_atm(recnum+1)
4502 :
4503 0 : fieldname = 'airtmp'
4504 : call read_data_nc_hycom (read6, recnum, &
4505 : tair_file, fieldname, Tair_data, & ! LCOV_EXCL_LINE
4506 0 : field_loc_center, field_type_scalar)
4507 0 : fieldname = 'wndewd'
4508 : call read_data_nc_hycom (read6, recnum, &
4509 : uwind_file, fieldname, uatm_data, & ! LCOV_EXCL_LINE
4510 0 : field_loc_center, field_type_vector)
4511 0 : fieldname = 'wndnwd'
4512 : call read_data_nc_hycom (read6, recnum, &
4513 : vwind_file, fieldname, vatm_data, & ! LCOV_EXCL_LINE
4514 0 : field_loc_center, field_type_vector)
4515 0 : fieldname = 'vapmix'
4516 : call read_data_nc_hycom (read6, recnum, &
4517 : humid_file, fieldname, Qa_data, & ! LCOV_EXCL_LINE
4518 0 : field_loc_center, field_type_scalar)
4519 0 : fieldname = 'shwflx'
4520 : call read_data_nc_hycom (read6, recnum, &
4521 : fsw_file, fieldname, fsw_data, & ! LCOV_EXCL_LINE
4522 0 : field_loc_center, field_type_scalar)
4523 0 : fieldname = 'radflx'
4524 : call read_data_nc_hycom (read6, recnum, &
4525 : flw_file, fieldname, flw_data, & ! LCOV_EXCL_LINE
4526 0 : field_loc_center, field_type_scalar)
4527 0 : fieldname = 'precip'
4528 : call read_data_nc_hycom (read6, recnum, &
4529 : rain_file, fieldname, fsnow_data,& ! LCOV_EXCL_LINE
4530 0 : field_loc_center, field_type_scalar)
4531 :
4532 : else
4533 0 : call abort_ice(subname//'ERROR: atm_data_format unavailable for hycom')
4534 : endif
4535 :
4536 : ! Interpolate
4537 0 : if (debug_forcing) then
4538 0 : if (my_task == master_task) then
4539 0 : write(nu_diag,*)'CICE: Atm. interpolate: = ',&
4540 0 : hcdate,c1intp,c2intp
4541 : endif
4542 : endif
4543 0 : call interpolate_data (Tair_data, Tair)
4544 0 : call interpolate_data (uatm_data, uatm)
4545 0 : call interpolate_data (vatm_data, vatm)
4546 0 : call interpolate_data ( fsw_data, fsw)
4547 0 : call interpolate_data ( flw_data, flw)
4548 0 : call interpolate_data ( Qa_data, Qa)
4549 0 : call interpolate_data (fsnow_data, fsnow)
4550 :
4551 : ! Adjust data forcing to CICE units
4552 0 : !$OMP PARALLEL DO PRIVATE(iblk,i,j)
4553 0 : do iblk = 1, nblocks
4554 0 : do j = 1, ny_block
4555 0 : do i = 1, nx_block
4556 : ! Air temperature: Degrees --> Kelvin
4557 0 : Tair(i,j,iblk) = Tair(i,j,iblk) + Tffresh
4558 : enddo ! i
4559 : enddo ! j
4560 : enddo ! nblocks
4561 : !$OMP END PARALLEL DO
4562 :
4563 : ! Save record number for next time step
4564 0 : oldrecnum = recnum
4565 :
4566 0 : end subroutine hycom_atm_data
4567 :
4568 : !=======================================================================
4569 : !
4570 0 : subroutine read_data_nc_point (flag, recd, yr, ixm, ixx, ixp, &
4571 : maxrec, data_file, fieldname, field_data, & ! LCOV_EXCL_LINE
4572 : field_loc, field_type)
4573 : !
4574 : ! If data is at the beginning of a one-year record, get data from
4575 : ! the previous year.
4576 : ! If data is at the end of a one-year record, get data from the
4577 : ! following year.
4578 : ! If no earlier data exists (beginning of fyear_init), then
4579 : ! (1) For monthly data, get data from the end of fyear_final.
4580 : ! (2) For more frequent data, let the ixm value equal the
4581 : ! first value of the year.
4582 : ! If no later data exists (end of fyear_final), then
4583 : ! (1) For monthly data, get data from the beginning of fyear_init.
4584 : ! (2) For more frequent data, let the ixp value
4585 : ! equal the last value of the year.
4586 : ! In other words, we assume persistence when daily or 6-hourly
4587 : ! data is missing, and we assume periodicity when monthly data
4588 : ! is missing.
4589 : !
4590 : use ice_diagnostics, only: debug_model_step
4591 :
4592 : logical (kind=log_kind), intent(in) :: flag
4593 :
4594 : integer (kind=int_kind), intent(in) :: &
4595 : recd , & ! baseline record number ! LCOV_EXCL_LINE
4596 : yr , & ! year of forcing data ! LCOV_EXCL_LINE
4597 : ixm, ixx, ixp , & ! record numbers of 3 data values ! LCOV_EXCL_LINE
4598 : ! relative to recd
4599 : maxrec ! maximum record value
4600 :
4601 : character (char_len_long), intent(in) :: &
4602 : data_file ! data file to be read
4603 :
4604 : character (char_len), intent(in) :: &
4605 : fieldname ! field name in netCDF file
4606 :
4607 : integer (kind=int_kind), intent(in) :: &
4608 : field_loc, & ! location of field on staggered grid ! LCOV_EXCL_LINE
4609 : field_type ! type of field (scalar, vector, angle)
4610 :
4611 : real (kind=dbl_kind), dimension(2), intent(inout) :: &
4612 : field_data ! 2 values needed for interpolation
4613 :
4614 : integer (kind=int_kind) :: &
4615 : nrec , & ! record number to read ! LCOV_EXCL_LINE
4616 : n2, n4 , & ! like ixm and ixp, but ! LCOV_EXCL_LINE
4617 : ! adjusted at beginning and end of data
4618 : arg , & ! value of time argument in field_data
4619 : fid ! file id for netCDF routines
4620 :
4621 : character(len=*), parameter :: subname = '(read_data_nc_point)'
4622 :
4623 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
4624 :
4625 0 : call ice_timer_start(timer_readwrite) ! reading/writing
4626 :
4627 0 : field_data = c0 ! to satisfy intent(out) attribute
4628 :
4629 0 : if (istep1 > debug_model_step) debug_forcing = .true. !! debugging
4630 :
4631 0 : if (my_task==master_task .and. (debug_forcing)) then
4632 0 : write(nu_diag,*) ' ', trim(data_file)
4633 : endif
4634 :
4635 0 : if (flag) then
4636 :
4637 : !-----------------------------------------------------------------
4638 : ! Initialize record counters
4639 : ! (n2, n4 will change only at the very beginning or end of
4640 : ! a forcing cycle.)
4641 : !-----------------------------------------------------------------
4642 0 : n2 = ixm
4643 0 : n4 = ixp
4644 0 : arg = 0
4645 :
4646 : !-----------------------------------------------------------------
4647 : ! read data
4648 : !-----------------------------------------------------------------
4649 :
4650 0 : if (ixm /= -99) then
4651 : ! currently in first half of data interval
4652 0 : if (ixx <= 1) then
4653 0 : if (yr > fyear_init) then ! get data from previous year
4654 : !call file_year (data_file, yr-1)
4655 : else ! yr = fyear_init, no prior data exists
4656 0 : if (maxrec > 12) then ! extrapolate from first record
4657 0 : if (ixx == 1) n2 = ixx
4658 : else ! go to end of fyear_final
4659 : ! call file_year (data_file, fyear_final)
4660 : endif
4661 : endif ! yr > fyear_init
4662 : endif ! ixx <= 1
4663 :
4664 : ! write(nu_diag,*) '!! read_data_nc !!!', trim(data_file)
4665 : ! write(nu_diag,*) 'istep ', istep
4666 : ! write(nu_diag,*) 'fyear_final ', fyear_final
4667 : ! write(nu_diag,*) 'fyear_init ', fyear_init
4668 : ! write(nu_diag,*) 'ixm, ixx, ixp ', ixm, ixx, ixp
4669 : ! write(nu_diag,*) 'maxrec ', maxrec
4670 : ! write(nu_diag,*) 'fieldname ', fieldname
4671 :
4672 0 : call ice_open_nc (data_file, fid)
4673 :
4674 0 : arg = 1
4675 0 : nrec = recd + n2
4676 :
4677 : call ice_read_nc &
4678 : (fid, nrec, fieldname, field_data(arg), debug_forcing, & ! LCOV_EXCL_LINE
4679 0 : field_loc, field_type)
4680 :
4681 : !if (ixx==1) call ice_close_nc(fid)
4682 0 : call ice_close_nc(fid)
4683 : endif ! ixm ne -99
4684 :
4685 : ! always read ixx data from data file for current year
4686 : ! call file_year (data_file, yr)
4687 0 : call ice_open_nc (data_file, fid)
4688 :
4689 0 : arg = arg + 1
4690 0 : nrec = recd + ixx
4691 :
4692 : call ice_read_nc &
4693 : (fid, nrec, fieldname, field_data(arg), debug_forcing, & ! LCOV_EXCL_LINE
4694 0 : field_loc, field_type)
4695 :
4696 0 : if (ixp /= -99) then
4697 : ! currently in latter half of data interval
4698 0 : if (ixx==maxrec) then
4699 0 : if (yr < fyear_final) then ! get data from following year
4700 0 : call ice_close_nc(fid)
4701 : !call file_year (data_file, yr+1)
4702 0 : call ice_open_nc (data_file, fid)
4703 : else ! yr = fyear_final, no more data exists
4704 0 : if (maxrec > 12) then ! extrapolate from ixx
4705 0 : n4 = ixx
4706 : else ! go to beginning of fyear_init
4707 0 : call ice_close_nc(fid)
4708 : ! call file_year (data_file, fyear_init)
4709 0 : call ice_open_nc (data_file, fid)
4710 :
4711 : endif
4712 : endif ! yr < fyear_final
4713 : endif ! ixx = maxrec
4714 :
4715 0 : arg = arg + 1
4716 0 : nrec = recd + n4
4717 :
4718 : call ice_read_nc &
4719 : (fid, nrec, fieldname, field_data(arg), debug_forcing, & ! LCOV_EXCL_LINE
4720 0 : field_loc, field_type)
4721 : endif ! ixp /= -99
4722 :
4723 0 : call ice_close_nc(fid)
4724 :
4725 : endif ! flag
4726 :
4727 0 : call ice_timer_stop(timer_readwrite) ! reading/writing
4728 :
4729 0 : end subroutine read_data_nc_point
4730 :
4731 : !=======================================================================
4732 :
4733 0 : subroutine ISPOL_files
4734 :
4735 : character(len=*), parameter :: subname = '(ISPOL_files)'
4736 :
4737 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
4738 :
4739 : fsw_file = &
4740 0 : trim(atm_data_dir)//'/fsw_sfc_4Xdaily.nc'
4741 :
4742 : flw_file = &
4743 0 : trim(atm_data_dir)//'/flw_sfc_4Xdaily.nc'
4744 :
4745 : rain_file = &
4746 0 : trim(atm_data_dir)//'/fsnow_sfc_daily_mod3.nc'
4747 :
4748 : uwind_file = &
4749 0 : trim(atm_data_dir)//'/uatm_10m_daily.nc'
4750 :
4751 : vwind_file = &
4752 0 : trim(atm_data_dir)//'/vatm_10m_daily.nc'
4753 :
4754 : tair_file = &
4755 0 : trim(atm_data_dir)//'/Tair_2m_daily.nc'
4756 :
4757 : humid_file = &
4758 0 : trim(atm_data_dir)//'/Qa_2m_daily.nc'
4759 :
4760 0 : if (my_task == master_task) then
4761 0 : write (nu_diag,*) ' '
4762 0 : write (nu_diag,*) 'Atmospheric data files:'
4763 0 : write (nu_diag,*) trim(fsw_file)
4764 0 : write (nu_diag,*) trim(flw_file)
4765 0 : write (nu_diag,*) trim(rain_file)
4766 0 : write (nu_diag,*) trim(uwind_file)
4767 0 : write (nu_diag,*) trim(vwind_file)
4768 0 : write (nu_diag,*) trim(tair_file)
4769 0 : write (nu_diag,*) trim(humid_file)
4770 : endif ! master_task
4771 :
4772 0 : end subroutine ISPOL_files
4773 :
4774 : !=======================================================================
4775 :
4776 0 : subroutine ISPOL_data
4777 :
4778 : ! Defines atmospheric data fields for Antarctic Weddell sea location
4779 :
4780 : ! authors: Nicole Jeffery, LANL
4781 : !
4782 : use ice_flux, only: uatm, vatm, Tair, fsw, Qa, rhoa, &
4783 : frain, fsnow, flw
4784 :
4785 : !local parameters
4786 :
4787 : character (char_len_long) :: &
4788 : met_file, & ! netcdf filename ! LCOV_EXCL_LINE
4789 : fieldname ! field name in netcdf file
4790 :
4791 : real (kind=dbl_kind), dimension(2), save :: &
4792 : Tair_data_p , & ! air temperature (K) for interpolation ! LCOV_EXCL_LINE
4793 : Qa_data_p, fsnow_data_p, & ! LCOV_EXCL_LINE
4794 : fsw_data_p, flw_data_p, & ! LCOV_EXCL_LINE
4795 : uatm_data_p, vatm_data_p
4796 :
4797 : real (kind=dbl_kind), parameter :: & ! coefficients for Hyland-Wexler Qa
4798 : ps1 = 0.58002206e4_dbl_kind, & ! (K) ! LCOV_EXCL_LINE
4799 : ps2 = 1.3914993_dbl_kind, & ! ! LCOV_EXCL_LINE
4800 : ps3 = 0.48640239e-1_dbl_kind, & ! (K^-1) ! LCOV_EXCL_LINE
4801 : ps4 = 0.41764768e-4_dbl_kind, & ! (K^-2) ! LCOV_EXCL_LINE
4802 : ps5 = 0.14452093e-7_dbl_kind, & ! (K^-3) ! LCOV_EXCL_LINE
4803 : ps6 = 6.5459673_dbl_kind, & ! ! LCOV_EXCL_LINE
4804 : ws1 = 621.97_dbl_kind, & ! for saturation mixing ratio ! LCOV_EXCL_LINE
4805 : Pair = 1020._dbl_kind, & ! Sea level pressure (hPa) ! LCOV_EXCL_LINE
4806 : lapse_rate = 0.0065_dbl_kind ! (K/m) lapse rate over sea level
4807 :
4808 : ! for interpolation of hourly data
4809 : integer (kind=int_kind) :: &
4810 : ixm,ixx,ixp , & ! record numbers for neighboring months ! LCOV_EXCL_LINE
4811 : maxrec , & ! maximum record number ! LCOV_EXCL_LINE
4812 : recslot , & ! spline slot for current record ! LCOV_EXCL_LINE
4813 : dataloc ! = 1 for data located in middle of time interval
4814 : ! = 2 for date located at end of time interval
4815 : real (kind=dbl_kind) :: &
4816 : secday , & ! LCOV_EXCL_LINE
4817 0 : Qa_pnt
4818 :
4819 : real (kind=dbl_kind) :: &
4820 0 : sec1hr ! number of seconds in 1 hour
4821 :
4822 : logical (kind=log_kind) :: read1
4823 :
4824 : integer (kind=int_kind) :: &
4825 : recnum , & ! record number ! LCOV_EXCL_LINE
4826 : recnum4X ! record number
4827 :
4828 : character(len=*), parameter :: subname = '(ISPOL_data)'
4829 :
4830 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
4831 :
4832 0 : call icepack_query_parameters(secday_out=secday)
4833 0 : call icepack_warnings_flush(nu_diag)
4834 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
4835 0 : file=__FILE__, line=__LINE__)
4836 :
4837 0 : if (trim(atm_data_format) == 'nc') then ! read nc file
4838 :
4839 : !-------------------------------------------------------------------
4840 : ! data from NCEP_DOE Reanalysis 2 and Bareiss et al 2008
4841 : ! daily data located at the end of the 24-hour period.
4842 : !-------------------------------------------------------------------
4843 :
4844 0 : dataloc = 2 ! data located at end of interval
4845 0 : sec1hr = secday ! seconds in day
4846 0 : maxrec = 366 !
4847 :
4848 : ! current record number
4849 0 : recnum = int(yday)
4850 :
4851 : ! Compute record numbers for surrounding data (2 on each side)
4852 0 : ixm = mod(recnum+maxrec-2,maxrec) + 1
4853 0 : ixx = mod(recnum-1, maxrec) + 1
4854 : ! ixp = mod(recnum, maxrec) + 1
4855 :
4856 : ! Compute interpolation coefficients
4857 : ! If data is located at the end of the time interval, then the
4858 : ! data value for the current record goes in slot 2
4859 :
4860 0 : recslot = 2
4861 0 : ixp = -99
4862 0 : call interp_coeff (recnum, recslot, sec1hr, dataloc)
4863 :
4864 0 : read1 = .false.
4865 0 : if (istep==1 .or. oldrecnum .ne. recnum) read1 = .true.
4866 :
4867 : ! Daily 2m Air temperature 1991
4868 :
4869 0 : met_file = tair_file
4870 0 : fieldname='Tair'
4871 :
4872 : call read_data_nc_point(read1, 0, fyear, ixm, ixx, ixp, &
4873 : maxrec, met_file, fieldname, Tair_data_p, & ! LCOV_EXCL_LINE
4874 0 : field_loc_center, field_type_scalar)
4875 :
4876 0 : Tair(:,:,:) = c1intp * Tair_data_p(1) &
4877 : + c2intp * Tair_data_p(2) & ! LCOV_EXCL_LINE
4878 0 : - lapse_rate*8.0_dbl_kind
4879 :
4880 0 : met_file = humid_file
4881 0 : fieldname='Qa'
4882 :
4883 : call read_data_nc_point(read1, 0, fyear, ixm, ixx, ixp, &
4884 : maxrec, met_file, fieldname, Qa_data_p, & ! LCOV_EXCL_LINE
4885 0 : field_loc_center, field_type_scalar)
4886 :
4887 : Qa_pnt= c1intp * Qa_data_p(1) &
4888 0 : + c2intp * Qa_data_p(2)
4889 0 : Qa(:,:,:) = Qa_pnt
4890 :
4891 0 : met_file = uwind_file
4892 0 : fieldname='uatm'
4893 :
4894 : call read_data_nc_point(read1, 0, fyear, ixm, ixx, ixp, &
4895 : maxrec, met_file, fieldname, uatm_data_p, & ! LCOV_EXCL_LINE
4896 0 : field_loc_center, field_type_scalar)
4897 :
4898 0 : uatm(:,:,:) = c1intp * uatm_data_p(1) &
4899 0 : + c2intp * uatm_data_p(2)
4900 :
4901 0 : met_file = vwind_file
4902 0 : fieldname='vatm'
4903 :
4904 : call read_data_nc_point(read1, 0, fyear, ixm, ixx, ixp, &
4905 : maxrec, met_file, fieldname, vatm_data_p, & ! LCOV_EXCL_LINE
4906 0 : field_loc_center, field_type_scalar)
4907 :
4908 0 : vatm(:,:,:) = c1intp * vatm_data_p(1) &
4909 0 : + c2intp * vatm_data_p(2)
4910 :
4911 0 : met_file = rain_file
4912 0 : fieldname='fsnow'
4913 :
4914 : call read_data_nc_point(read1, 0, fyear, ixm, ixx, ixp, &
4915 : maxrec, met_file, fieldname, fsnow_data_p, & ! LCOV_EXCL_LINE
4916 0 : field_loc_center, field_type_scalar)
4917 :
4918 0 : fsnow(:,:,:) = (c1intp * fsnow_data_p(1) + &
4919 0 : c2intp * fsnow_data_p(2))
4920 :
4921 : !-----------------------------
4922 : !fsw and flw are every 6 hours
4923 : !------------------------------
4924 0 : dataloc = 2 ! data located at end of interval
4925 0 : sec1hr = secday/c4 ! seconds in 6 hours
4926 0 : maxrec = 1460 ! 366*4
4927 :
4928 : ! current record number
4929 0 : recnum4X = 4*int(yday) - 3 + int(real(msec,kind=dbl_kind)/sec1hr)
4930 :
4931 : ! Compute record numbers for surrounding data (2 on each side)
4932 0 : ixm = mod(recnum4X+maxrec-2,maxrec) + 1
4933 0 : ixx = mod(recnum4X-1, maxrec) + 1
4934 :
4935 : ! Compute interpolation coefficients
4936 : ! If data is located at the end of the time interval, then the
4937 : ! data value for the current record goes in slot 2
4938 :
4939 0 : recslot = 2
4940 0 : ixp = -99
4941 0 : call interp_coeff (recnum4X, recslot, sec1hr, dataloc)
4942 :
4943 0 : read1 = .false.
4944 0 : if (istep==1 .or. oldrecnum4X .ne. recnum4X) read1 = .true.
4945 :
4946 0 : met_file = fsw_file
4947 0 : fieldname='fsw'
4948 :
4949 : call read_data_nc_point(read1, 0, fyear, ixm, ixx, ixp, &
4950 : maxrec, met_file, fieldname, fsw_data_p, & ! LCOV_EXCL_LINE
4951 0 : field_loc_center, field_type_scalar)
4952 :
4953 0 : fsw(:,:,:) = c1intp * fsw_data_p(1) &
4954 0 : + c2intp * fsw_data_p(2)
4955 :
4956 0 : met_file = flw_file
4957 0 : fieldname='flw'
4958 :
4959 : call read_data_nc_point(read1, 0, fyear, ixm, ixx, ixp, &
4960 : maxrec, met_file, fieldname, flw_data_p, & ! LCOV_EXCL_LINE
4961 0 : field_loc_center, field_type_scalar)
4962 :
4963 0 : flw(:,:,:) = c1intp * flw_data_p(1) &
4964 0 : + c2intp * flw_data_p(2)
4965 : endif !nc
4966 :
4967 : !flw given cldf and Tair calculated in prepare_forcing
4968 :
4969 : !-----------------------------
4970 : ! fixed data
4971 : ! May not be needed
4972 : !-----------------------------
4973 0 : rhoa (:,:,:) = 1.3_dbl_kind ! air density (kg/m^3)
4974 0 : cldf(:,:,:) = c1 !0.25_dbl_kind ! cloud fraction
4975 0 : frain(:,:,:) = c0 ! this is available in hourlymet_rh file
4976 :
4977 : ! Save record number for next time step
4978 0 : oldrecnum = recnum
4979 0 : oldrecnum4X = recnum4X
4980 :
4981 0 : end subroutine ISPOL_data
4982 :
4983 : !=======================================================================
4984 :
4985 0 : subroutine ocn_data_ispol_init
4986 :
4987 : ! Reads NCAR pop ocean forcing data set 'pop_frc_gx1v3_010815.nc'
4988 : ! at the ISPOL location -67.4677N, 310.4375E
4989 : !
4990 : ! For ocean mixed layer-----------------------------units
4991 : !
4992 : ! 1 sst------temperature---------------------------(C)
4993 : ! 2 sss------salinity------------------------------(ppt)
4994 : ! 3 hbl------depth---------------------------------(m)
4995 : ! 4 u--------surface u current---------------------(m/s)
4996 : ! 5 v--------surface v current---------------------(m/s)
4997 : ! 6 dhdx-----surface tilt x direction--------------(m/m)
4998 : ! 7 dhdy-----surface tilt y direction--------------(m/m)
4999 : ! 8 qdp------ocean sub-mixed layer heat flux-------(W/m2)
5000 : !
5001 : ! Fields 4, 5, 6, 7 are on the U-grid; 1, 2, 3, and 8 are
5002 : ! on the T-grid.
5003 : !
5004 : ! authors: Nicole Jeffery, LANL
5005 : !
5006 : use ice_gather_scatter
5007 : use ice_read_write
5008 :
5009 : integer (kind=int_kind) :: &
5010 : n , & ! field index ! LCOV_EXCL_LINE
5011 : m ! month index
5012 :
5013 : character(char_len) :: &
5014 : vname(nfld) ! variable names to search for in file
5015 : data vname / &
5016 : 'T', 'S', 'hblt', 'U', 'V', & ! LCOV_EXCL_LINE
5017 : 'dhdx', 'dhdy', 'qdp' /
5018 :
5019 : real (kind=dbl_kind) :: &
5020 0 : work
5021 :
5022 : integer (kind=int_kind) :: &
5023 : fid ! file id
5024 :
5025 : character(len=*), parameter :: subname = '(ocn_data_ispol_init)'
5026 :
5027 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
5028 :
5029 0 : if (my_task == master_task) then
5030 :
5031 0 : if (restore_ocn) write (nu_diag,*) &
5032 0 : 'SST restoring timescale = ',trestore,' days'
5033 :
5034 0 : sst_file = trim(ocn_data_dir)//'/'//trim(oceanmixed_file) ! not just sst
5035 :
5036 : !---------------------------------------------------------------
5037 : ! Read in ocean forcing data from an existing file
5038 : !---------------------------------------------------------------
5039 0 : write (nu_diag,*) 'ocean mixed layer forcing data file = ', &
5040 0 : sst_file
5041 :
5042 : endif ! master_task
5043 :
5044 0 : if (trim(ocn_data_format) == 'nc') then
5045 0 : if (my_task == master_task) then
5046 0 : call ice_open_nc(sst_file, fid)
5047 : endif ! master_task
5048 :
5049 : ! Read in ocean forcing data for all 12 months
5050 0 : do n=1,nfld
5051 0 : do m=1,12
5052 : ! Note: netCDF does single to double conversion if necessary
5053 0 : if (n >= 4 .and. n <= 7) then
5054 0 : call ice_read_nc(fid, m, vname(n), work, debug_forcing, &
5055 0 : field_loc_NEcorner, field_type_vector)
5056 : else
5057 0 : call ice_read_nc(fid, m, vname(n), work, debug_forcing, &
5058 0 : field_loc_center, field_type_scalar)
5059 : endif
5060 0 : ocn_frc_m(:,:,:,n,m) = work
5061 : enddo ! month loop
5062 : enddo ! field loop
5063 :
5064 0 : if (my_task == master_task) call ice_close_nc(fid)
5065 :
5066 : else ! binary format
5067 : call abort_ice (error_message=subname//'new ocean forcing is netcdf only', &
5068 0 : file=__FILE__, line=__LINE__)
5069 : endif
5070 :
5071 : !echmod - currents cause Fram outflow to be too large
5072 0 : ocn_frc_m(:,:,:,4,:) = c0
5073 0 : ocn_frc_m(:,:,:,5,:) = c0
5074 : !echmod
5075 :
5076 0 : end subroutine ocn_data_ispol_init
5077 :
5078 : !=======================================================================
5079 : !
5080 2896 : subroutine box2001_data_atm
5081 :
5082 : ! wind fields as in Hunke, JCP 2001
5083 : ! these are defined at the u point
5084 : ! authors: Elizabeth Hunke, LANL
5085 :
5086 : use ice_domain, only: nblocks, blocks_ice
5087 : use ice_calendar, only: timesecs
5088 : use ice_blocks, only: block, get_block, nx_block, ny_block, nghost
5089 : use ice_flux, only: uatm, vatm, wind, rhoa, strax, stray
5090 : use ice_state, only: aice
5091 :
5092 : ! local parameters
5093 :
5094 : integer (kind=int_kind) :: &
5095 : iblk, i,j ! loop indices
5096 :
5097 : integer (kind=int_kind) :: &
5098 : iglob(nx_block), & ! global indices ! LCOV_EXCL_LINE
5099 5792 : jglob(ny_block) ! global indices
5100 :
5101 : type (block) :: &
5102 : this_block ! block information for current block
5103 :
5104 : real (kind=dbl_kind) :: &
5105 0 : secday, pi , puny, period, pi2, tau
5106 :
5107 : character(len=*), parameter :: subname = '(box2001_data_atm)'
5108 :
5109 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
5110 :
5111 2896 : call icepack_query_parameters(pi_out=pi, pi2_out=pi2, puny_out=puny)
5112 2896 : call icepack_query_parameters(secday_out=secday)
5113 :
5114 2896 : period = c4*secday
5115 :
5116 8688 : do iblk = 1, nblocks
5117 205616 : do j = 1, ny_block
5118 6898272 : do i = 1, nx_block
5119 :
5120 6695552 : this_block = get_block(blocks_ice(iblk),iblk)
5121 234344320 : iglob = this_block%i_glob
5122 234344320 : jglob = this_block%j_glob
5123 :
5124 : ! wind components
5125 0 : uatm(i,j,iblk) = c5 + (sin(pi2*timesecs/period)-c3) &
5126 : * sin(pi2*real(iglob(i), kind=dbl_kind) & ! LCOV_EXCL_LINE
5127 : /real(nx_global,kind=dbl_kind)) & ! LCOV_EXCL_LINE
5128 : * sin(pi *real(jglob(j), kind=dbl_kind) & ! LCOV_EXCL_LINE
5129 6695552 : /real(ny_global,kind=dbl_kind))
5130 0 : vatm(i,j,iblk) = c5 + (sin(pi2*timesecs/period)-c3) &
5131 : * sin(pi *real(iglob(i), kind=dbl_kind) & ! LCOV_EXCL_LINE
5132 : /real(nx_global,kind=dbl_kind)) & ! LCOV_EXCL_LINE
5133 : * sin(pi2*real(jglob(j), kind=dbl_kind) & ! LCOV_EXCL_LINE
5134 6695552 : /real(ny_global,kind=dbl_kind))
5135 : ! wind stress
5136 6695552 : wind(i,j,iblk) = sqrt(uatm(i,j,iblk)**2 + vatm(i,j,iblk)**2)
5137 6695552 : tau = rhoa(i,j,iblk) * 0.0012_dbl_kind * wind(i,j,iblk)
5138 :
5139 6695552 : strax(i,j,iblk) = aice(i,j,iblk) * tau * uatm(i,j,iblk)
5140 6892480 : stray(i,j,iblk) = aice(i,j,iblk) * tau * vatm(i,j,iblk)
5141 :
5142 : ! initialization test
5143 : ! Diagonal wind vectors 1
5144 : !uatm(i,j,iblk) = c1 *real(j-nghost, kind=dbl_kind) &
5145 : ! / real(ny_global,kind=dbl_kind)
5146 : !vatm(i,j,iblk) = c1 *real(j-nghost, kind=dbl_kind) &
5147 : ! / real(ny_global,kind=dbl_kind)
5148 :
5149 : ! Diagonal wind vectors 2
5150 : !uatm(i,j,iblk) = c1 *real(i-nghost, kind=dbl_kind) &
5151 : ! / real(nx_global,kind=dbl_kind)
5152 : !vatm(i,j,iblk) = -c1 *real(i-nghost, kind=dbl_kind) &
5153 : ! / real(nx_global,kind=dbl_kind)
5154 :
5155 : ! Wind in x direction
5156 : ! uatm(i,j,iblk) = c1 *real(i-nghost, kind=dbl_kind) &
5157 : ! / real(nx_global,kind=dbl_kind)
5158 : ! vatm(i,j,iblk) = c0
5159 :
5160 : ! Wind in y direction
5161 : ! uatm(i,j,iblk) = c0
5162 : ! vatm(i,j,iblk) = c1 *real(j-nghost, kind=dbl_kind) &
5163 : ! / real(ny_global,kind=dbl_kind)
5164 : ! initialization test
5165 :
5166 : enddo
5167 : enddo
5168 : enddo ! nblocks
5169 :
5170 2896 : end subroutine box2001_data_atm
5171 :
5172 : !=======================================================================
5173 : !
5174 2896 : subroutine box2001_data_ocn
5175 :
5176 : ! current fields as in Hunke, JCP 2001
5177 : ! these are defined at the u point
5178 : ! authors: Elizabeth Hunke, LANL
5179 :
5180 : use ice_domain, only: nblocks, blocks_ice
5181 : use ice_blocks, only: block, get_block, nx_block, ny_block, nghost
5182 : use ice_flux, only: uocn, vocn
5183 : use ice_grid, only: uvm
5184 :
5185 : ! local parameters
5186 :
5187 : integer (kind=int_kind) :: &
5188 : iblk, i,j ! loop indices
5189 :
5190 : integer (kind=int_kind) :: &
5191 : iglob(nx_block), & ! global indices ! LCOV_EXCL_LINE
5192 5792 : jglob(ny_block) ! global indices
5193 :
5194 : type (block) :: &
5195 : this_block ! block information for current block
5196 :
5197 : character(len=*), parameter :: subname = '(box2001_data_ocn)'
5198 :
5199 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
5200 :
5201 8688 : do iblk = 1, nblocks
5202 205616 : do j = 1, ny_block
5203 6898272 : do i = 1, nx_block
5204 :
5205 6695552 : this_block = get_block(blocks_ice(iblk),iblk)
5206 234344320 : iglob = this_block%i_glob
5207 234344320 : jglob = this_block%j_glob
5208 :
5209 : ! ocean current
5210 : ! constant in time, could be initialized in ice_flux.F90
5211 0 : uocn(i,j,iblk) = p2*real(jglob(j), kind=dbl_kind) &
5212 6695552 : / real(ny_global,kind=dbl_kind) - p1
5213 0 : vocn(i,j,iblk) = -p2*real(iglob(i), kind=dbl_kind) &
5214 6695552 : / real(nx_global,kind=dbl_kind) + p1
5215 :
5216 6695552 : uocn(i,j,iblk) = uocn(i,j,iblk) * uvm(i,j,iblk)
5217 6892480 : vocn(i,j,iblk) = vocn(i,j,iblk) * uvm(i,j,iblk)
5218 :
5219 : enddo
5220 : enddo
5221 : enddo ! nblocks
5222 :
5223 2896 : end subroutine box2001_data_ocn
5224 :
5225 : !=======================================================================
5226 : !
5227 0 : subroutine uniform_data_atm(dir,spd)
5228 : ! uniform wind fields in some direction
5229 :
5230 : use ice_domain, only: nblocks
5231 : use ice_blocks, only: nx_block, ny_block, nghost
5232 : use ice_flux, only: uatm, vatm, wind, rhoa, strax, stray
5233 : use ice_state, only: aice
5234 :
5235 : character(len=*), intent(in) :: dir
5236 : real(kind=dbl_kind), intent(in), optional :: spd ! velocity
5237 :
5238 : ! local parameters
5239 :
5240 : integer (kind=int_kind) :: &
5241 : iblk, i,j ! loop indices
5242 :
5243 : real (kind=dbl_kind) :: &
5244 : tau, & ! LCOV_EXCL_LINE
5245 0 : atm_val ! value to use for atm speed
5246 :
5247 : character(len=*), parameter :: subname = '(uniform_data_atm)'
5248 :
5249 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
5250 :
5251 : ! check for optional spd
5252 0 : if (present(spd)) then
5253 0 : atm_val = spd
5254 : else
5255 0 : atm_val = c5 ! default
5256 : endif
5257 :
5258 : ! wind components
5259 0 : if (dir == 'NE') then
5260 0 : uatm = atm_val
5261 0 : vatm = atm_val
5262 0 : elseif (dir == 'N') then
5263 0 : uatm = c0
5264 0 : vatm = atm_val
5265 0 : elseif (dir == 'E') then
5266 0 : uatm = atm_val
5267 0 : vatm = c0
5268 0 : elseif (dir == 'S') then
5269 0 : uatm = c0
5270 0 : vatm = -atm_val
5271 0 : elseif (dir == 'W') then
5272 0 : uatm = -atm_val
5273 0 : vatm = c0
5274 : else
5275 : call abort_ice (subname//'ERROR: dir unknown, dir = '//trim(dir), &
5276 0 : file=__FILE__, line=__LINE__)
5277 : endif
5278 :
5279 0 : do iblk = 1, nblocks
5280 0 : do j = 1, ny_block
5281 0 : do i = 1, nx_block
5282 :
5283 : ! wind stress
5284 0 : wind(i,j,iblk) = sqrt(uatm(i,j,iblk)**2 + vatm(i,j,iblk)**2)
5285 0 : tau = rhoa(i,j,iblk) * 0.0012_dbl_kind * wind(i,j,iblk)
5286 0 : strax(i,j,iblk) = aice(i,j,iblk) * tau * uatm(i,j,iblk)
5287 0 : stray(i,j,iblk) = aice(i,j,iblk) * tau * vatm(i,j,iblk)
5288 :
5289 : enddo
5290 : enddo
5291 : enddo ! nblocks
5292 :
5293 0 : end subroutine uniform_data_atm
5294 : !=======================================================================
5295 :
5296 : !
5297 0 : subroutine uniform_data_ocn(dir,spd)
5298 :
5299 : ! uniform current fields in some direction
5300 :
5301 : use ice_flux, only: uocn, vocn
5302 :
5303 : character(len=*), intent(in) :: dir
5304 :
5305 : real(kind=dbl_kind), intent(in), optional :: spd ! velocity
5306 :
5307 : ! local parameters
5308 :
5309 : real(kind=dbl_kind) :: &
5310 0 : ocn_val ! value to use for ocean currents
5311 :
5312 : character(len=*), parameter :: subname = '(uniform_data_ocn)'
5313 :
5314 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
5315 :
5316 0 : if (present(spd)) then
5317 0 : ocn_val = spd
5318 : else
5319 0 : ocn_val = p1 ! default
5320 : endif
5321 :
5322 : ! ocn components
5323 0 : if (dir == 'NE') then
5324 0 : uocn = ocn_val
5325 0 : vocn = ocn_val
5326 0 : elseif (dir == 'N') then
5327 0 : uocn = c0
5328 0 : vocn = ocn_val
5329 0 : elseif (dir == 'E') then
5330 0 : uocn = ocn_val
5331 0 : vocn = c0
5332 : else
5333 : call abort_ice (subname//'ERROR: dir unknown, dir = '//trim(dir), &
5334 0 : file=__FILE__, line=__LINE__)
5335 : endif
5336 :
5337 0 : end subroutine uniform_data_ocn
5338 : !=======================================================================
5339 :
5340 0 : subroutine get_wave_spec
5341 :
5342 : use ice_read_write, only: ice_read_nc_xyf
5343 : use ice_arrays_column, only: wave_spectrum, &
5344 : dwavefreq, wavefreq
5345 : use ice_constants, only: c0
5346 : use ice_timers, only: ice_timer_start, ice_timer_stop, timer_fsd
5347 :
5348 : ! local variables
5349 : integer (kind=int_kind) :: &
5350 : fid ! file id for netCDF routines
5351 :
5352 : real(kind=dbl_kind), dimension(nfreq) :: &
5353 0 : wave_spectrum_profile ! wave spectrum
5354 :
5355 : character(char_len) :: wave_spec_type
5356 : logical (kind=log_kind) :: wave_spec
5357 : character(len=*), parameter :: subname = '(get_wave_spec)'
5358 :
5359 : if (local_debug .and. my_task == master_task) write(nu_diag,*) subname,'fdbg start'
5360 :
5361 0 : call ice_timer_start(timer_fsd)
5362 :
5363 : call icepack_query_parameters(wave_spec_out=wave_spec, &
5364 0 : wave_spec_type_out=wave_spec_type)
5365 0 : call icepack_warnings_flush(nu_diag)
5366 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
5367 0 : file=__FILE__, line=__LINE__)
5368 :
5369 : ! if no wave data is provided, wave_spectrum is zero everywhere
5370 0 : wave_spectrum(:,:,:,:) = c0
5371 0 : debug_forcing = .false.
5372 :
5373 : ! wave spectrum and frequencies
5374 0 : if (wave_spec) then
5375 : ! get hardwired frequency bin info and a dummy wave spectrum profile
5376 : ! the latter is used if wave_spec_type == profile
5377 : call icepack_init_wave(nfreq, &
5378 : wave_spectrum_profile, & ! LCOV_EXCL_LINE
5379 0 : wavefreq, dwavefreq)
5380 :
5381 : ! read more realistic data from a file
5382 0 : if ((trim(wave_spec_type) == 'constant').OR.(trim(wave_spec_type) == 'random')) then
5383 0 : if (trim(wave_spec_file(1:4)) == 'unkn') then
5384 : call abort_ice (subname//'ERROR: wave_spec_file '//trim(wave_spec_file), &
5385 0 : file=__FILE__, line=__LINE__)
5386 : else
5387 : #ifdef USE_NETCDF
5388 0 : call wave_spec_data
5389 : #else
5390 : write (nu_diag,*) "wave spectrum file not available, requires cpp USE_NETCDF"
5391 : write (nu_diag,*) "wave spectrum file not available, using default profile"
5392 : call abort_ice (subname//'ERROR: wave_spec_file '//trim(wave_spec_file), &
5393 : file=__FILE__, line=__LINE__)
5394 : #endif
5395 : endif
5396 : endif
5397 : endif
5398 :
5399 0 : call ice_timer_stop(timer_fsd)
5400 :
5401 0 : end subroutine get_wave_spec
5402 :
5403 : !=======================================================================
5404 : !
5405 : ! Read in wave spectrum forcing as a function of time. 6 hourly
5406 : ! LR started working from JRA55_data routine
5407 : ! Changed fields, and changed 3 hourly to 6 hourly
5408 : !
5409 0 : subroutine wave_spec_data
5410 :
5411 : use ice_blocks, only: block, get_block
5412 : use ice_global_reductions, only: global_minval, global_maxval
5413 : use ice_domain, only: nblocks, distrb_info, blocks_ice
5414 : use ice_arrays_column, only: wave_spectrum, &
5415 : dwavefreq, wavefreq
5416 : use ice_read_write, only: ice_read_nc_xyf
5417 : use ice_grid, only: hm, tlon, tlat, tmask, umask
5418 : use ice_calendar, only: days_per_year, use_leap_years
5419 :
5420 : integer (kind=int_kind) :: &
5421 : ncid , & ! netcdf file id ! LCOV_EXCL_LINE
5422 : i, j, freq , & ! LCOV_EXCL_LINE
5423 : ixm,ixx,ixp , & ! record numbers for neighboring months ! LCOV_EXCL_LINE
5424 : recnum , & ! record number ! LCOV_EXCL_LINE
5425 : maxrec , & ! maximum record number ! LCOV_EXCL_LINE
5426 : recslot , & ! spline slot for current record ! LCOV_EXCL_LINE
5427 : midmonth , & ! middle day of month ! LCOV_EXCL_LINE
5428 : dataloc , & ! = 1 for data located in middle of time interval ! LCOV_EXCL_LINE
5429 : ! = 2 for date located at end of time interval
5430 : iblk , & ! block index
5431 : ilo,ihi,jlo,jhi, & ! beginning and end of physical domain ! LCOV_EXCL_LINE
5432 : yr ! current forcing year
5433 :
5434 : real (kind=dbl_kind) :: &
5435 : sec6hr , & ! number of seconds in 3 hours ! LCOV_EXCL_LINE
5436 : secday , & ! number of seconds in day ! LCOV_EXCL_LINE
5437 : vmin, vmax
5438 :
5439 : logical (kind=log_kind) :: readm, read6,debug_n_d
5440 :
5441 : type (block) :: &
5442 : this_block ! block information for current block
5443 :
5444 : real(kind=dbl_kind), dimension(nfreq) :: &
5445 0 : wave_spectrum_profile ! wave spectrum
5446 :
5447 : character(len=64) :: fieldname !netcdf field name
5448 : character(char_len_long) :: spec_file
5449 : character(char_len) :: wave_spec_type
5450 : logical (kind=log_kind) :: wave_spec
5451 : character(len=*), parameter :: subname = '(wave_spec_data)'
5452 :
5453 :
5454 :
5455 0 : debug_n_d = .false. !usually false
5456 :
5457 0 : call icepack_query_parameters(secday_out=secday)
5458 0 : call icepack_warnings_flush(nu_diag)
5459 0 : if (icepack_warnings_aborted()) call abort_ice(error_message=subname, &
5460 0 : file=__FILE__, line=__LINE__)
5461 :
5462 : call icepack_init_wave(nfreq, &
5463 : wave_spectrum_profile, & ! LCOV_EXCL_LINE
5464 0 : wavefreq, dwavefreq)
5465 :
5466 :
5467 : !spec_file = trim(ocn_data_dir)//'/'//trim(wave_spec_file)
5468 0 : spec_file = trim(wave_spec_file)
5469 0 : wave_spectrum_data = c0
5470 0 : wave_spectrum = c0
5471 0 : yr = fyear ! current year
5472 : !-------------------------------------------------------------------
5473 : ! 6-hourly data
5474 : !
5475 : ! Assume that the 6-hourly value is located at the end of the
5476 : ! 6-hour period. This is the convention for NCEP reanalysis data.
5477 : ! E.g. record 1 gives conditions at 6 am GMT on 1 January.
5478 : !-------------------------------------------------------------------
5479 :
5480 0 : dataloc = 2 ! data located at end of interval
5481 0 : sec6hr = secday/c4 ! seconds in 6 hours
5482 : !maxrec = 2920 ! 365*8; for leap years = 366*8
5483 :
5484 0 : if (use_leap_years) days_per_year = 366 !overrides setting of 365 in ice_calendar
5485 0 : maxrec = days_per_year*4
5486 :
5487 0 : if(days_per_year == 365 .and. (mod(yr, 4) == 0)) then
5488 0 : call abort_ice('days_per_year should be set to 366 for leap years')
5489 : end if
5490 :
5491 : ! current record number
5492 0 : recnum = 4*int(yday) - 3 + int(real(msec,kind=dbl_kind)/sec6hr)
5493 :
5494 : ! Compute record numbers for surrounding data (2 on each side)
5495 :
5496 0 : ixm = mod(recnum+maxrec-2,maxrec) + 1
5497 0 : ixx = mod(recnum-1, maxrec) + 1
5498 :
5499 : ! Compute interpolation coefficients
5500 : ! If data is located at the end of the time interval, then the
5501 : ! data value for the current record goes in slot 2
5502 :
5503 0 : recslot = 2
5504 0 : ixp = -99
5505 0 : call interp_coeff (recnum, recslot, sec6hr, dataloc)
5506 :
5507 : ! Read
5508 0 : read6 = .false.
5509 0 : if (istep==1 .or. oldrecnum .ne. recnum) read6 = .true.
5510 : !-------------------------------------------------------------------
5511 : ! File is NETCDF
5512 : ! file variable names are:
5513 : ! efreq (wave spectrum, energy as a function of wave frequency UNITS)
5514 : !-------------------------------------------------------------------
5515 0 : call ice_open_nc(spec_file,ncid)
5516 :
5517 0 : call ice_read_nc_xyf(ncid,recnum,'efreq',wave_spectrum_data(:,:,:,1,:),debug_n_d, &
5518 : field_loc=field_loc_center, & ! LCOV_EXCL_LINE
5519 0 : field_type=field_type_scalar)
5520 0 : call ice_read_nc_xyf(ncid,recnum,'efreq',wave_spectrum_data(:,:,:,2,:),debug_n_d, &
5521 : field_loc=field_loc_center, & ! LCOV_EXCL_LINE
5522 0 : field_type=field_type_scalar)
5523 0 : call ice_close_nc(ncid)
5524 :
5525 :
5526 : ! Interpolate
5527 0 : call interpolate_wavespec_data (wave_spectrum_data, wave_spectrum)
5528 :
5529 : ! Save record number
5530 0 : oldrecnum = recnum
5531 :
5532 : if (local_debug) then
5533 : if (my_task == master_task) write (nu_diag,*) &
5534 : 'wave_spec_data ',spec_file
5535 : if (my_task.eq.master_task) &
5536 : write (nu_diag,*) 'maxrec',maxrec
5537 : write (nu_diag,*) 'days_per_year', days_per_year
5538 :
5539 : endif ! local debug
5540 :
5541 0 : end subroutine wave_spec_data
5542 :
5543 : !=======================================================================
5544 :
5545 : ! initial snow aging lookup table
5546 : !
5547 : ! Dry snow metamorphism table
5548 : ! snicar_drdt_bst_fit_60_c070416.nc
5549 : ! Flanner (file metadata units mislabelled)
5550 : ! drdsdt0 (10^-6 m/hr) tau (10^-6 m)
5551 : !
5552 0 : subroutine init_snowtable
5553 :
5554 : use ice_broadcast, only: broadcast_array, broadcast_scalar
5555 : integer (kind=int_kind) :: &
5556 : idx_T_max , & ! Table dimensions ! LCOV_EXCL_LINE
5557 : idx_rhos_max, & ! LCOV_EXCL_LINE
5558 : idx_Tgrd_max
5559 : real (kind=dbl_kind), allocatable :: &
5560 : snowage_rhos (:), & ! LCOV_EXCL_LINE
5561 : snowage_Tgrd (:), & ! LCOV_EXCL_LINE
5562 : snowage_T (:), & ! LCOV_EXCL_LINE
5563 : snowage_tau (:,:,:), & ! LCOV_EXCL_LINE
5564 : snowage_kappa(:,:,:), & ! LCOV_EXCL_LINE
5565 0 : snowage_drdt0(:,:,:)
5566 :
5567 : ! local variables
5568 :
5569 : logical (kind=log_kind) :: diag = .false.
5570 :
5571 : integer (kind=int_kind) :: &
5572 : fid ! file id for netCDF file
5573 :
5574 : character (char_len) :: &
5575 : snw_aging_table, & ! aging table setting ! LCOV_EXCL_LINE
5576 : fieldname ! field name in netcdf file
5577 :
5578 : character(len=*), parameter :: subname = '(init_snowtable)'
5579 :
5580 : !-----------------------------------------------------------------
5581 : ! read table of snow aging parameters
5582 : !-----------------------------------------------------------------
5583 :
5584 : call icepack_query_parameters(snw_aging_table_out=snw_aging_table, &
5585 0 : isnw_rhos_out=idx_rhos_max, isnw_Tgrd_out=idx_Tgrd_max, isnw_T_out=idx_T_max)
5586 :
5587 0 : if (my_task == master_task) then
5588 0 : write (nu_diag,*) ' '
5589 0 : write (nu_diag,*) 'Snow aging file:', trim(snw_filename)
5590 : endif
5591 :
5592 0 : if (snw_aging_table == 'snicar') then
5593 : ! just read the 3d data and pass it in
5594 :
5595 0 : call ice_open_nc(snw_filename,fid)
5596 :
5597 0 : allocate(snowage_tau (idx_rhos_max, idx_Tgrd_max, idx_T_max))
5598 0 : allocate(snowage_kappa(idx_rhos_max, idx_Tgrd_max, idx_T_max))
5599 0 : allocate(snowage_drdt0(idx_rhos_max, idx_Tgrd_max, idx_T_max))
5600 :
5601 0 : fieldname = trim(snw_tau_fname)
5602 : call ice_read_nc(fid,fieldname,snowage_tau, diag, &
5603 0 : idx_rhos_max,idx_Tgrd_max,idx_T_max)
5604 0 : fieldname = trim(snw_kappa_fname)
5605 : call ice_read_nc(fid,fieldname,snowage_kappa,diag, &
5606 0 : idx_rhos_max,idx_Tgrd_max,idx_T_max)
5607 0 : fieldname = trim(snw_drdt0_fname)
5608 : call ice_read_nc(fid,fieldname,snowage_drdt0,diag, &
5609 0 : idx_rhos_max,idx_Tgrd_max,idx_T_max)
5610 :
5611 0 : call ice_close_nc(fid)
5612 :
5613 0 : call broadcast_array(snowage_tau , master_task)
5614 0 : call broadcast_array(snowage_kappa, master_task)
5615 0 : call broadcast_array(snowage_drdt0, master_task)
5616 :
5617 0 : if (my_task == master_task) then
5618 0 : write(nu_diag,*) subname,' '
5619 0 : write(nu_diag,*) subname,' Successfully read snow aging properties:'
5620 0 : write(nu_diag,*) subname,' snw_aging_table = ',trim(snw_aging_table)
5621 0 : write(nu_diag,*) subname,' idx_rhos_max = ',idx_rhos_max
5622 0 : write(nu_diag,*) subname,' idx_Tgrd_max = ',idx_Tgrd_max
5623 0 : write(nu_diag,*) subname,' idx_T_max = ',idx_T_max
5624 0 : write(nu_diag,*) subname,' Data at rhos, Tgrd, T at first index '
5625 0 : write(nu_diag,*) subname,' snoage_tau (1,1,1) = ',snowage_tau (1,1,1)
5626 0 : write(nu_diag,*) subname,' snoage_kappa (1,1,1) = ',snowage_kappa(1,1,1)
5627 0 : write(nu_diag,*) subname,' snoage_drdt0 (1,1,1) = ',snowage_drdt0(1,1,1)
5628 0 : write(nu_diag,*) subname,' Data at rhos, Tgrd, T at max index'
5629 0 : write(nu_diag,*) subname,' snoage_tau (max,max,max) = ',snowage_tau (idx_rhos_max, idx_Tgrd_max, idx_T_max)
5630 0 : write(nu_diag,*) subname,' snoage_kappa (max,max,max) = ',snowage_kappa(idx_rhos_max, idx_Tgrd_max, idx_T_max)
5631 0 : write(nu_diag,*) subname,' snoage_drdt0 (max,max,max) = ',snowage_drdt0(idx_rhos_max, idx_Tgrd_max, idx_T_max)
5632 : endif
5633 :
5634 : call icepack_init_parameters( &
5635 : snowage_tau_in = snowage_tau, & ! LCOV_EXCL_LINE
5636 : snowage_kappa_in = snowage_kappa, & ! LCOV_EXCL_LINE
5637 0 : snowage_drdt0_in = snowage_drdt0 )
5638 :
5639 0 : deallocate(snowage_tau)
5640 0 : deallocate(snowage_kappa)
5641 0 : deallocate(snowage_drdt0)
5642 :
5643 : else
5644 : ! read everything and pass it in
5645 :
5646 0 : call ice_open_nc(snw_filename,fid)
5647 :
5648 0 : fieldname = trim(snw_rhos_fname)
5649 0 : call ice_get_ncvarsize(fid,fieldname,idx_rhos_max)
5650 0 : fieldname = trim(snw_Tgrd_fname)
5651 0 : call ice_get_ncvarsize(fid,fieldname,idx_Tgrd_max)
5652 0 : fieldname = trim(snw_T_fname)
5653 0 : call ice_get_ncvarsize(fid,fieldname,idx_T_max)
5654 :
5655 0 : call broadcast_scalar(idx_rhos_max, master_task)
5656 0 : call broadcast_scalar(idx_Tgrd_max, master_task)
5657 0 : call broadcast_scalar(idx_T_max , master_task)
5658 :
5659 0 : allocate(snowage_rhos (idx_rhos_max))
5660 0 : allocate(snowage_Tgrd (idx_Tgrd_max))
5661 0 : allocate(snowage_T (idx_T_max))
5662 0 : allocate(snowage_tau (idx_rhos_max, idx_Tgrd_max, idx_T_max))
5663 0 : allocate(snowage_kappa(idx_rhos_max, idx_Tgrd_max, idx_T_max))
5664 0 : allocate(snowage_drdt0(idx_rhos_max, idx_Tgrd_max, idx_T_max))
5665 :
5666 0 : fieldname = trim(snw_rhos_fname)
5667 : call ice_read_nc(fid,fieldname,snowage_rhos, diag, &
5668 0 : idx_rhos_max)
5669 0 : fieldname = trim(snw_Tgrd_fname)
5670 : call ice_read_nc(fid,fieldname,snowage_Tgrd, diag, &
5671 0 : idx_Tgrd_max)
5672 0 : fieldname = trim(snw_T_fname)
5673 : call ice_read_nc(fid,fieldname,snowage_T, diag, &
5674 0 : idx_T_max)
5675 :
5676 0 : fieldname = trim(snw_tau_fname)
5677 : call ice_read_nc(fid,fieldname,snowage_tau, diag, &
5678 0 : idx_rhos_max,idx_Tgrd_max,idx_T_max)
5679 0 : fieldname = trim(snw_kappa_fname)
5680 : call ice_read_nc(fid,fieldname,snowage_kappa,diag, &
5681 0 : idx_rhos_max,idx_Tgrd_max,idx_T_max)
5682 0 : fieldname = trim(snw_drdt0_fname)
5683 : call ice_read_nc(fid,fieldname,snowage_drdt0,diag, &
5684 0 : idx_rhos_max,idx_Tgrd_max,idx_T_max)
5685 :
5686 0 : call ice_close_nc(fid)
5687 :
5688 0 : call broadcast_array(snowage_rhos , master_task)
5689 0 : call broadcast_array(snowage_Tgrd , master_task)
5690 0 : call broadcast_array(snowage_T , master_task)
5691 0 : call broadcast_array(snowage_tau , master_task)
5692 0 : call broadcast_array(snowage_kappa, master_task)
5693 0 : call broadcast_array(snowage_drdt0, master_task)
5694 :
5695 0 : if (my_task == master_task) then
5696 0 : write(nu_diag,*) subname,' '
5697 0 : write(nu_diag,*) subname,' Successfully read snow aging properties:'
5698 0 : write(nu_diag,*) subname,' idx_rhos_max = ',idx_rhos_max
5699 0 : write(nu_diag,*) subname,' idx_Tgrd_max = ',idx_Tgrd_max
5700 0 : write(nu_diag,*) subname,' idx_T_max = ',idx_T_max
5701 0 : write(nu_diag,*) subname,' Data at rhos, Tgrd, T = ',snowage_rhos(1),snowage_Tgrd(1),snowage_T(1)
5702 0 : write(nu_diag,*) subname,' snoage_tau (1,1,1) = ',snowage_tau (1,1,1)
5703 0 : write(nu_diag,*) subname,' snoage_kappa (1,1,1) = ',snowage_kappa(1,1,1)
5704 0 : write(nu_diag,*) subname,' snoage_drdt0 (1,1,1) = ',snowage_drdt0(1,1,1)
5705 0 : write(nu_diag,*) subname,' Data at rhos, Tgrd, T = ', &
5706 0 : snowage_rhos(idx_rhos_max),snowage_Tgrd(idx_Tgrd_max),snowage_T(idx_T_max)
5707 0 : write(nu_diag,*) subname,' snoage_tau (max,max,max) = ',snowage_tau (idx_rhos_max, idx_Tgrd_max, idx_T_max)
5708 0 : write(nu_diag,*) subname,' snoage_kappa (max,max,max) = ',snowage_kappa(idx_rhos_max, idx_Tgrd_max, idx_T_max)
5709 0 : write(nu_diag,*) subname,' snoage_drdt0 (max,max,max) = ',snowage_drdt0(idx_rhos_max, idx_Tgrd_max, idx_T_max)
5710 : endif
5711 :
5712 : call icepack_init_parameters( &
5713 : isnw_t_in = idx_T_max, & ! LCOV_EXCL_LINE
5714 : isnw_Tgrd_in = idx_Tgrd_max, & ! LCOV_EXCL_LINE
5715 : isnw_rhos_in = idx_rhos_max, & ! LCOV_EXCL_LINE
5716 : snowage_rhos_in = snowage_rhos, & ! LCOV_EXCL_LINE
5717 : snowage_Tgrd_in = snowage_Tgrd, & ! LCOV_EXCL_LINE
5718 : snowage_T_in = snowage_T, & ! LCOV_EXCL_LINE
5719 : snowage_tau_in = snowage_tau, & ! LCOV_EXCL_LINE
5720 : snowage_kappa_in = snowage_kappa, & ! LCOV_EXCL_LINE
5721 0 : snowage_drdt0_in = snowage_drdt0 )
5722 :
5723 0 : deallocate(snowage_rhos)
5724 0 : deallocate(snowage_Tgrd)
5725 0 : deallocate(snowage_T)
5726 0 : deallocate(snowage_tau)
5727 0 : deallocate(snowage_kappa)
5728 0 : deallocate(snowage_drdt0)
5729 :
5730 : endif
5731 :
5732 0 : end subroutine init_snowtable
5733 :
5734 : !=======================================================================
5735 :
5736 : end module ice_forcing
5737 :
5738 : !=======================================================================
|
