Line data Source code
1 : !=======================================================================
2 :
3 : ! parameter and variable initializations
4 : !
5 : ! authors Elizabeth C. Hunke, LANL
6 :
7 : module icedrv_init
8 :
9 : use icedrv_kinds
10 : use icedrv_constants, only: nu_diag, ice_stdout, nu_diag_out, nu_nml
11 : use icedrv_constants, only: c0, c1, c2, c3, p2, p5
12 : use icedrv_domain_size, only: nx
13 : use icepack_intfc, only: icepack_init_parameters
14 : use icepack_intfc, only: icepack_init_fsd
15 : use icepack_intfc, only: icepack_init_tracer_flags
16 : use icepack_intfc, only: icepack_init_tracer_sizes
17 : use icepack_intfc, only: icepack_init_tracer_indices
18 : use icepack_intfc, only: icepack_init_trcr
19 : use icepack_intfc, only: icepack_query_parameters
20 : use icepack_intfc, only: icepack_query_tracer_flags
21 : use icepack_intfc, only: icepack_query_tracer_sizes
22 : use icepack_intfc, only: icepack_query_tracer_indices
23 : use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted
24 : use icedrv_system, only: icedrv_system_abort
25 :
26 : implicit none
27 : private
28 : public :: input_data, init_grid2, init_state, init_fsd
29 :
30 : character(len=char_len_long), public :: &
31 : ice_ic ! method of ice cover initialization
32 : ! 'default' or 'none' => conditions specified in code
33 : ! restart = .true. overwrites default initial
34 : ! condition using filename given by ice_ic
35 :
36 : real (kind=dbl_kind), dimension (nx), public :: &
37 : TLON , & ! longitude of temp pts (radians)
38 : TLAT ! latitude of temp pts (radians)
39 :
40 : logical (kind=log_kind), &
41 : dimension (nx), public :: &
42 : tmask , & ! land/boundary mask, thickness (T-cell)
43 : lmask_n, & ! northern hemisphere mask
44 : lmask_s ! southern hemisphere mask
45 :
46 : !=======================================================================
47 :
48 : contains
49 :
50 : !=======================================================================
51 :
52 : ! Namelist variables, set to default values; may be altered
53 : ! at run time
54 : !
55 : ! author Elizabeth C. Hunke, LANL
56 :
57 46 : subroutine input_data
58 :
59 : use icedrv_diagnostics, only: diag_file, nx_names
60 : use icedrv_domain_size, only: nilyr, nslyr, nblyr, max_ntrcr, ncat
61 : use icedrv_domain_size, only: n_iso, n_aero, nfsd
62 : use icedrv_calendar, only: year_init, istep0
63 : use icedrv_calendar, only: dumpfreq, diagfreq, dump_last
64 : use icedrv_calendar, only: npt, dt, ndtd, days_per_year, use_leap_years
65 : use icedrv_restart_shared, only: restart, restart_dir, restart_file
66 : use icedrv_flux, only: update_ocn_f, l_mpond_fresh, cpl_bgc
67 : use icedrv_flux, only: default_season
68 : use icedrv_forcing, only: precip_units, fyear_init, ycycle
69 : use icedrv_forcing, only: atm_data_type, ocn_data_type, bgc_data_type
70 : use icedrv_forcing, only: atm_data_file, ocn_data_file, bgc_data_file
71 : use icedrv_forcing, only: ice_data_file
72 : use icedrv_forcing, only: atm_data_format, ocn_data_format, bgc_data_format
73 : use icedrv_forcing, only: data_dir
74 : use icedrv_forcing, only: oceanmixed_ice, restore_ocn, trestore
75 :
76 : ! local variables
77 :
78 : character (32) :: &
79 : nml_filename = 'icepack_in' ! namelist input file name
80 :
81 : integer (kind=int_kind) :: &
82 : nml_error, & ! namelist i/o error flag
83 : n ! loop index
84 :
85 : character (len=char_len) :: diag_file_names
86 :
87 34 : real (kind=dbl_kind) :: ustar_min, albicev, albicei, albsnowv, albsnowi, &
88 85 : ahmax, R_ice, R_pnd, R_snw, dT_mlt, rsnw_mlt, ksno, &
89 85 : mu_rdg, hs0, dpscale, rfracmin, rfracmax, pndaspect, hs1, hp1, &
90 68 : a_rapid_mode, Rac_rapid_mode, aspect_rapid_mode, dSdt_slow_mode, &
91 51 : phi_c_slow_mode, phi_i_mushy, kalg, emissivity
92 :
93 : integer (kind=int_kind) :: ktherm, kstrength, krdg_partic, krdg_redist, &
94 : natmiter, kitd, kcatbound
95 :
96 : character (len=char_len) :: shortwave, albedo_type, conduct, fbot_xfer_type, &
97 : tfrz_option, frzpnd, atmbndy, wave_spec_type
98 :
99 : logical (kind=log_kind) :: sw_redist
100 17 : real (kind=dbl_kind) :: sw_frac, sw_dtemp
101 :
102 : ! Flux convergence tolerance
103 17 : real (kind=dbl_kind) :: atmiter_conv
104 :
105 : logical (kind=log_kind) :: calc_Tsfc, formdrag, highfreq, calc_strair
106 : logical (kind=log_kind) :: conserv_check
107 :
108 : integer (kind=int_kind) :: ntrcr
109 : logical (kind=log_kind) :: tr_iage, tr_FY, tr_lvl, tr_pond
110 : logical (kind=log_kind) :: tr_iso, tr_aero, tr_fsd
111 : logical (kind=log_kind) :: tr_pond_cesm, tr_pond_lvl, tr_pond_topo, wave_spec
112 : integer (kind=int_kind) :: nt_Tsfc, nt_sice, nt_qice, nt_qsno, nt_iage, nt_FY
113 : integer (kind=int_kind) :: nt_alvl, nt_vlvl, nt_apnd, nt_hpnd, nt_ipnd, &
114 : nt_aero, nt_fsd, nt_isosno, nt_isoice
115 :
116 17 : real (kind=real_kind) :: rpcesm, rplvl, rptopo
117 34 : real (kind=dbl_kind) :: Cf, puny
118 : character(len=*), parameter :: subname='(input_data)'
119 :
120 : !-----------------------------------------------------------------
121 : ! Namelist variables
122 : !-----------------------------------------------------------------
123 :
124 : namelist /setup_nml/ &
125 : days_per_year, use_leap_years, year_init, istep0, &
126 : dt, npt, ndtd, dump_last, &
127 : ice_ic, restart, restart_dir, restart_file, &
128 : dumpfreq, diagfreq, diag_file, cpl_bgc, &
129 : conserv_check
130 :
131 : namelist /grid_nml/ &
132 : kcatbound
133 :
134 : namelist /thermo_nml/ &
135 : kitd, ktherm, ksno, conduct, &
136 : a_rapid_mode, Rac_rapid_mode, aspect_rapid_mode, &
137 : dSdt_slow_mode, phi_c_slow_mode, phi_i_mushy, &
138 : sw_redist, sw_frac, sw_dtemp
139 :
140 : namelist /dynamics_nml/ &
141 : kstrength, krdg_partic, krdg_redist, mu_rdg, &
142 : Cf
143 :
144 : namelist /shortwave_nml/ &
145 : shortwave, albedo_type, &
146 : albicev, albicei, albsnowv, albsnowi, &
147 : ahmax, R_ice, R_pnd, R_snw, &
148 : dT_mlt, rsnw_mlt, kalg
149 :
150 : namelist /ponds_nml/ &
151 : hs0, dpscale, frzpnd, &
152 : rfracmin, rfracmax, pndaspect, hs1, &
153 : hp1
154 :
155 : namelist /forcing_nml/ &
156 : atmbndy, calc_strair, calc_Tsfc, &
157 : update_ocn_f, l_mpond_fresh, ustar_min, &
158 : fbot_xfer_type, oceanmixed_ice, emissivity, &
159 : formdrag, highfreq, natmiter, &
160 : atmiter_conv, &
161 : tfrz_option, default_season, wave_spec_type, &
162 : precip_units, fyear_init, ycycle, &
163 : atm_data_type, ocn_data_type, bgc_data_type, &
164 : atm_data_file, ocn_data_file, bgc_data_file, &
165 : ice_data_file, &
166 : atm_data_format, ocn_data_format, bgc_data_format, &
167 : data_dir, trestore, restore_ocn
168 :
169 : namelist /tracer_nml/ &
170 : tr_iage, &
171 : tr_FY, &
172 : tr_lvl, &
173 : tr_pond_cesm, &
174 : tr_pond_lvl, &
175 : tr_pond_topo, &
176 : tr_aero, &
177 : tr_fsd, &
178 : tr_iso
179 :
180 : !-----------------------------------------------------------------
181 : ! query Icepack values
182 : !-----------------------------------------------------------------
183 :
184 : call icepack_query_parameters(ustar_min_out=ustar_min, Cf_out=Cf, &
185 : albicev_out=albicev, albicei_out=albicei, ksno_out = ksno, &
186 : albsnowv_out=albsnowv, albsnowi_out=albsnowi, &
187 : natmiter_out=natmiter, ahmax_out=ahmax, shortwave_out=shortwave, &
188 : atmiter_conv_out = atmiter_conv, &
189 : albedo_type_out=albedo_type, R_ice_out=R_ice, R_pnd_out=R_pnd, &
190 : R_snw_out=R_snw, dT_mlt_out=dT_mlt, rsnw_mlt_out=rsnw_mlt, &
191 : kstrength_out=kstrength, krdg_partic_out=krdg_partic, &
192 : krdg_redist_out=krdg_redist, mu_rdg_out=mu_rdg, &
193 : atmbndy_out=atmbndy, calc_strair_out=calc_strair, &
194 : formdrag_out=formdrag, highfreq_out=highfreq, &
195 : emissivity_out=emissivity, &
196 : kitd_out=kitd, kcatbound_out=kcatbound, hs0_out=hs0, &
197 : dpscale_out=dpscale, frzpnd_out=frzpnd, &
198 : rfracmin_out=rfracmin, rfracmax_out=rfracmax, &
199 : pndaspect_out=pndaspect, hs1_out=hs1, hp1_out=hp1, &
200 : ktherm_out=ktherm, calc_Tsfc_out=calc_Tsfc, &
201 : update_ocn_f_out = update_ocn_f, &
202 : conduct_out=conduct, a_rapid_mode_out=a_rapid_mode, &
203 : Rac_rapid_mode_out=Rac_rapid_mode, &
204 : aspect_rapid_mode_out=aspect_rapid_mode, &
205 : dSdt_slow_mode_out=dSdt_slow_mode, &
206 : phi_c_slow_mode_out=phi_c_slow_mode, &
207 : phi_i_mushy_out=phi_i_mushy, conserv_check_out=conserv_check, &
208 : tfrz_option_out=tfrz_option, kalg_out=kalg, &
209 : fbot_xfer_type_out=fbot_xfer_type, puny_out=puny, &
210 : wave_spec_type_out=wave_spec_type, &
211 46 : sw_redist_out=sw_redist, sw_frac_out=sw_frac, sw_dtemp_out=sw_dtemp)
212 46 : call icepack_warnings_flush(nu_diag)
213 46 : if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
214 0 : file=__FILE__, line=__LINE__)
215 :
216 : !-----------------------------------------------------------------
217 : ! default values
218 : !-----------------------------------------------------------------
219 :
220 46 : days_per_year = 365 ! number of days in a year
221 46 : use_leap_years= .false.! if true, use leap years (Feb 29)
222 46 : year_init = 0 ! initial year
223 46 : istep0 = 0 ! no. of steps taken in previous integrations,
224 : ! real (dumped) or imagined (to set calendar)
225 46 : dt = 3600.0_dbl_kind ! time step, s
226 46 : npt = 99999 ! total number of time steps (dt)
227 46 : diagfreq = 24 ! how often diag output is written
228 46 : diag_file = 'ice_diag' ! history file name prefix
229 46 : cpl_bgc = .false. !
230 46 : dumpfreq='y' ! restart frequency option
231 46 : dump_last=.false. ! restart at end of run
232 46 : restart = .false. ! if true, read restart files for initialization
233 46 : restart_dir = './' ! write to executable dir for default
234 46 : restart_file = 'iced' ! restart file name prefix
235 46 : ice_ic = 'default' ! initial conditions are specified in the code
236 : ! otherwise, the filename for reading restarts
237 46 : ndtd = 1 ! dynamic time steps per thermodynamic time step
238 46 : l_mpond_fresh = .false. ! logical switch for including meltpond freshwater
239 : ! flux feedback to ocean model
240 46 : default_season = 'winter' ! default forcing data, if data is not read in
241 46 : fyear_init = 1998 ! initial forcing year
242 46 : ycycle = 1 ! number of years in forcing cycle
243 46 : atm_data_format = 'bin' ! file format ('bin'=binary or 'nc'=netcdf)
244 46 : atm_data_type = 'default' ! source of atmospheric forcing data
245 46 : atm_data_file = ' ' ! atmospheric forcing data file
246 46 : precip_units = 'mks' ! 'mm_per_month' or
247 : ! 'mm_per_sec' = 'mks' = kg/m^2 s
248 46 : oceanmixed_ice = .false. ! if true, use internal ocean mixed layer
249 46 : wave_spec_type = 'none' ! type of wave spectrum forcing
250 46 : ocn_data_format = 'bin' ! file format ('bin'=binary or 'nc'=netcdf)
251 46 : ocn_data_type = 'default' ! source of ocean forcing data
252 46 : ocn_data_file = ' ' ! ocean forcing data file
253 46 : ice_data_file = ' ' ! ice forcing data file (opening, closing)
254 46 : bgc_data_format = 'bin' ! file format ('bin'=binary or 'nc'=netcdf)
255 46 : bgc_data_type = 'default' ! source of BGC forcing data
256 46 : bgc_data_file = ' ' ! biogeochemistry forcing data file
257 46 : data_dir = ' ' ! root location of data files
258 46 : restore_ocn = .false. ! restore sst if true
259 46 : trestore = 90 ! restoring timescale, days (0 instantaneous)
260 :
261 : ! extra tracers
262 46 : tr_iage = .false. ! ice age
263 46 : tr_FY = .false. ! ice age
264 46 : tr_lvl = .false. ! level ice
265 46 : tr_pond_cesm = .false. ! CESM melt ponds
266 46 : tr_pond_lvl = .false. ! level-ice melt ponds
267 46 : tr_pond_topo = .false. ! explicit melt ponds (topographic)
268 46 : tr_aero = .false. ! aerosols
269 46 : tr_fsd = .false. ! floe size distribution
270 46 : tr_iso = .false. ! isotopes
271 :
272 : !-----------------------------------------------------------------
273 : ! read from input file
274 : !-----------------------------------------------------------------
275 :
276 46 : open (nu_nml, file=nml_filename, status='old',iostat=nml_error)
277 46 : if (nml_error /= 0) then
278 0 : nml_error = -1
279 : else
280 46 : nml_error = 1
281 : endif
282 :
283 92 : do while (nml_error > 0)
284 46 : print*,'Reading namelist file ',nml_filename
285 :
286 46 : print*,'Reading setup_nml'
287 46 : read(nu_nml, nml=setup_nml,iostat=nml_error)
288 46 : if (nml_error /= 0) exit
289 :
290 46 : print*,'Reading grid_nml'
291 46 : read(nu_nml, nml=grid_nml,iostat=nml_error)
292 46 : if (nml_error /= 0) exit
293 :
294 46 : print*,'Reading tracer_nml'
295 46 : read(nu_nml, nml=tracer_nml,iostat=nml_error)
296 46 : if (nml_error /= 0) exit
297 :
298 46 : print*,'Reading thermo_nml'
299 46 : read(nu_nml, nml=thermo_nml,iostat=nml_error)
300 46 : if (nml_error /= 0) exit
301 :
302 46 : print*,'Reading shortwave_nml'
303 46 : read(nu_nml, nml=shortwave_nml,iostat=nml_error)
304 46 : if (nml_error /= 0) exit
305 :
306 46 : print*,'Reading ponds_nml'
307 46 : read(nu_nml, nml=ponds_nml,iostat=nml_error)
308 46 : if (nml_error /= 0) exit
309 :
310 46 : print*,'Reading forcing_nml'
311 46 : read(nu_nml, nml=forcing_nml,iostat=nml_error)
312 92 : if (nml_error /= 0) exit
313 : end do
314 46 : if (nml_error == 0) close(nu_nml)
315 46 : if (nml_error /= 0) then
316 0 : write(ice_stdout,*) 'error reading namelist'
317 0 : call icedrv_system_abort(file=__FILE__,line=__LINE__)
318 : endif
319 46 : close(nu_nml)
320 :
321 : !-----------------------------------------------------------------
322 : ! set up diagnostics output and resolve conflicts
323 : !-----------------------------------------------------------------
324 :
325 46 : write(ice_stdout,*) 'Diagnostic output will be in files '
326 46 : write(ice_stdout,*)' ','icepack.runlog.timestamp'
327 :
328 230 : do n = 1,nx
329 230 : write(nx_names(n),'(a,i2.2)') 'point_',n
330 : enddo
331 46 : nx_names(1) = 'icefree'
332 46 : nx_names(2) = 'slab'
333 46 : nx_names(3) = 'full_ITD'
334 46 : nx_names(4) = 'land'
335 :
336 230 : do n = 1,nx
337 184 : diag_file_names=' '
338 184 : write(diag_file_names,'(a,a,a)') trim(diag_file),'.',trim(nx_names(n))
339 184 : write(ice_stdout,*)' ',trim(diag_file_names)
340 230 : open(nu_diag_out+n-1, file=diag_file_names, status='unknown')
341 : end do
342 :
343 46 : write(nu_diag,*) '-----------------------------------'
344 46 : write(nu_diag,*) ' ICEPACK model diagnostic output '
345 46 : write(nu_diag,*) '-----------------------------------'
346 46 : write(nu_diag,*) ' '
347 :
348 3 : if (ncat == 1 .and. kitd == 1) then
349 0 : write (nu_diag,*) 'Remapping the ITD is not allowed for ncat=1.'
350 0 : write (nu_diag,*) 'Use kitd = 0 (delta function ITD) with kcatbound = 0'
351 0 : write (nu_diag,*) 'or for column configurations use kcatbound = -1'
352 0 : call icedrv_system_abort(file=__FILE__,line=__LINE__)
353 : endif
354 :
355 43 : if (ncat /= 1 .and. kcatbound == -1) then
356 0 : write (nu_diag,*) 'WARNING: ITD required for ncat > 1'
357 0 : write (nu_diag,*) 'WARNING: Setting kitd and kcatbound to default values'
358 0 : kitd = 1
359 0 : kcatbound = 0
360 : endif
361 :
362 46 : rpcesm = c0
363 46 : rplvl = c0
364 46 : rptopo = c0
365 46 : if (tr_pond_cesm) rpcesm = c1
366 46 : if (tr_pond_lvl ) rplvl = c1
367 46 : if (tr_pond_topo) rptopo = c1
368 :
369 46 : tr_pond = .false. ! explicit melt ponds
370 46 : if (rpcesm + rplvl + rptopo > puny) tr_pond = .true.
371 :
372 46 : if (rpcesm + rplvl + rptopo > c1 + puny) then
373 0 : write (nu_diag,*) 'WARNING: Must use only one melt pond scheme'
374 0 : call icedrv_system_abort(file=__FILE__,line=__LINE__)
375 : endif
376 :
377 46 : if (tr_pond_lvl .and. .not. tr_lvl) then
378 0 : write (nu_diag,*) 'WARNING: tr_pond_lvl=T but tr_lvl=F'
379 0 : write (nu_diag,*) 'WARNING: Setting tr_lvl=T'
380 0 : tr_lvl = .true.
381 : endif
382 :
383 46 : if (tr_pond_lvl .and. abs(hs0) > puny) then
384 0 : write (nu_diag,*) 'WARNING: tr_pond_lvl=T and hs0/=0'
385 0 : write (nu_diag,*) 'WARNING: Setting hs0=0'
386 0 : hs0 = c0
387 : endif
388 :
389 46 : if (tr_pond_cesm .and. trim(frzpnd) /= 'cesm') then
390 2 : write (nu_diag,*) 'WARNING: tr_pond_cesm=T'
391 2 : write (nu_diag,*) 'WARNING: frzpnd, dpscale not used'
392 2 : frzpnd = 'cesm'
393 : endif
394 :
395 46 : if (trim(shortwave) /= 'dEdd' .and. tr_pond .and. calc_tsfc) then
396 0 : write (nu_diag,*) 'WARNING: Must use dEdd shortwave'
397 0 : write (nu_diag,*) 'WARNING: with tr_pond and calc_tsfc=T.'
398 0 : write (nu_diag,*) 'WARNING: Setting shortwave = dEdd'
399 0 : shortwave = 'dEdd'
400 : endif
401 :
402 44 : if (tr_iso .and. n_iso==0) then
403 0 : write (nu_diag,*) 'WARNING: isotopes activated but'
404 0 : write (nu_diag,*) 'WARNING: not allocated in tracer array.'
405 0 : write (nu_diag,*) 'WARNING: Activate in compilation script.'
406 0 : call icedrv_system_abort(file=__FILE__,line=__LINE__)
407 : endif
408 :
409 5 : if (tr_aero .and. n_aero==0) then
410 0 : write (nu_diag,*) 'WARNING: aerosols activated but'
411 0 : write (nu_diag,*) 'WARNING: not allocated in tracer array.'
412 0 : write (nu_diag,*) 'WARNING: Activate in compilation script.'
413 0 : call icedrv_system_abort(file=__FILE__,line=__LINE__)
414 : endif
415 :
416 46 : if (tr_aero .and. trim(shortwave) /= 'dEdd') then
417 0 : write (nu_diag,*) 'WARNING: aerosols activated but dEdd'
418 0 : write (nu_diag,*) 'WARNING: shortwave is not.'
419 0 : write (nu_diag,*) 'WARNING: Setting shortwave = dEdd'
420 0 : shortwave = 'dEdd'
421 : endif
422 :
423 46 : rfracmin = min(max(rfracmin,c0),c1)
424 46 : rfracmax = min(max(rfracmax,c0),c1)
425 :
426 46 : if (ktherm == 2 .and. .not. calc_Tsfc) then
427 0 : write (nu_diag,*) 'WARNING: ktherm = 2 and calc_Tsfc = F'
428 0 : write (nu_diag,*) 'WARNING: Setting calc_Tsfc = T'
429 0 : calc_Tsfc = .true.
430 : endif
431 :
432 46 : if (ktherm == 1 .and. trim(tfrz_option) /= 'linear_salt') then
433 : write (nu_diag,*) &
434 0 : 'WARNING: ktherm = 1 and tfrz_option = ',trim(tfrz_option)
435 : write (nu_diag,*) &
436 0 : 'WARNING: For consistency, set tfrz_option = linear_salt'
437 : endif
438 46 : if (ktherm == 2 .and. trim(tfrz_option) /= 'mushy') then
439 : write (nu_diag,*) &
440 0 : 'WARNING: ktherm = 2 and tfrz_option = ',trim(tfrz_option)
441 : write (nu_diag,*) &
442 0 : 'WARNING: For consistency, set tfrz_option = mushy'
443 : endif
444 :
445 46 : if (formdrag) then
446 3 : if (trim(atmbndy) == 'constant') then
447 0 : write (nu_diag,*) 'WARNING: atmbndy = constant not allowed with formdrag'
448 0 : write (nu_diag,*) 'WARNING: Setting atmbndy = default'
449 0 : atmbndy = 'default'
450 : endif
451 :
452 3 : if (.not. calc_strair) then
453 0 : write (nu_diag,*) 'WARNING: formdrag=T but calc_strair=F'
454 0 : write (nu_diag,*) 'WARNING: Setting calc_strair=T'
455 0 : calc_strair = .true.
456 : endif
457 :
458 3 : if (tr_pond_cesm) then
459 0 : write (nu_diag,*) 'ERROR: formdrag=T but frzpnd=cesm'
460 0 : call icedrv_system_abort(file=__FILE__,line=__LINE__)
461 : endif
462 :
463 3 : if (.not. tr_lvl) then
464 3 : write (nu_diag,*) 'WARNING: formdrag=T but tr_lvl=F'
465 3 : write (nu_diag,*) 'WARNING: Setting tr_lvl=T'
466 3 : tr_lvl = .true.
467 : endif
468 : endif
469 :
470 46 : if (trim(fbot_xfer_type) == 'Cdn_ocn' .and. .not. formdrag) then
471 0 : write (nu_diag,*) 'WARNING: formdrag=F but fbot_xfer_type=Cdn_ocn'
472 0 : write (nu_diag,*) 'WARNING: Setting fbot_xfer_type = constant'
473 0 : fbot_xfer_type = 'constant'
474 : endif
475 :
476 46 : wave_spec = .false.
477 46 : if (tr_fsd .and. (trim(wave_spec_type) /= 'none')) wave_spec = .true.
478 :
479 : !-----------------------------------------------------------------
480 : ! spew
481 : !-----------------------------------------------------------------
482 :
483 46 : write(nu_diag,*) ' Document ice_in namelist parameters:'
484 46 : write(nu_diag,*) ' ==================================== '
485 46 : write(nu_diag,*) ' '
486 46 : write(nu_diag,1020) ' days_per_year = ', days_per_year
487 46 : write(nu_diag,1010) ' use_leap_years = ', use_leap_years
488 46 : write(nu_diag,1020) ' year_init = ', year_init
489 46 : write(nu_diag,1020) ' istep0 = ', istep0
490 46 : write(nu_diag,1000) ' dt = ', dt
491 46 : write(nu_diag,1020) ' npt = ', npt
492 46 : write(nu_diag,1020) ' diagfreq = ', diagfreq
493 46 : write(nu_diag,1030) ' dumpfreq = ', &
494 92 : trim(dumpfreq)
495 46 : write(nu_diag,1010) ' dump_last = ', dump_last
496 46 : write(nu_diag,1010) ' restart = ', restart
497 46 : write(nu_diag,*) ' restart_dir = ', &
498 92 : trim(restart_dir)
499 46 : write(nu_diag,*) ' restart_file = ', &
500 92 : trim(restart_file)
501 46 : write(nu_diag,*) ' ice_ic = ', &
502 92 : trim(ice_ic)
503 46 : write(nu_diag,1010) ' conserv_check = ', conserv_check
504 46 : write(nu_diag,1020) ' kitd = ', kitd
505 46 : write(nu_diag,1020) ' kcatbound = ', &
506 92 : kcatbound
507 46 : write(nu_diag,1020) ' ndtd = ', ndtd
508 46 : write(nu_diag,1020) ' kstrength = ', kstrength
509 46 : write(nu_diag,1020) ' krdg_partic = ', &
510 92 : krdg_partic
511 46 : write(nu_diag,1020) ' krdg_redist = ', &
512 92 : krdg_redist
513 46 : if (krdg_redist == 1) &
514 46 : write(nu_diag,1000) ' mu_rdg = ', mu_rdg
515 46 : if (kstrength == 1) &
516 46 : write(nu_diag,1000) ' Cf = ', Cf
517 46 : write(nu_diag,1000) ' ksno = ', ksno
518 46 : write(nu_diag,1030) ' shortwave = ', &
519 92 : trim(shortwave)
520 46 : if (cpl_bgc) then
521 0 : write(nu_diag,1000) ' BGC coupling is switched ON'
522 : else
523 46 : write(nu_diag,1000) ' BGC coupling is switched OFF'
524 : endif
525 :
526 46 : if (trim(shortwave) == 'dEdd') then
527 40 : write(nu_diag,1000) ' R_ice = ', R_ice
528 40 : write(nu_diag,1000) ' R_pnd = ', R_pnd
529 40 : write(nu_diag,1000) ' R_snw = ', R_snw
530 40 : write(nu_diag,1000) ' dT_mlt = ', dT_mlt
531 40 : write(nu_diag,1000) ' rsnw_mlt = ', rsnw_mlt
532 40 : write(nu_diag,1000) ' kalg = ', kalg
533 40 : write(nu_diag,1000) ' hp1 = ', hp1
534 40 : write(nu_diag,1000) ' hs0 = ', hs0
535 : else
536 6 : write(nu_diag,1030) ' albedo_type = ', &
537 12 : trim(albedo_type)
538 6 : write(nu_diag,1000) ' albicev = ', albicev
539 6 : write(nu_diag,1000) ' albicei = ', albicei
540 6 : write(nu_diag,1000) ' albsnowv = ', albsnowv
541 6 : write(nu_diag,1000) ' albsnowi = ', albsnowi
542 6 : write(nu_diag,1000) ' ahmax = ', ahmax
543 : endif
544 :
545 46 : write(nu_diag,1000) ' rfracmin = ', rfracmin
546 46 : write(nu_diag,1000) ' rfracmax = ', rfracmax
547 46 : if (tr_pond_lvl) then
548 30 : write(nu_diag,1000) ' hs1 = ', hs1
549 30 : write(nu_diag,1000) ' dpscale = ', dpscale
550 30 : write(nu_diag,1030) ' frzpnd = ', trim(frzpnd)
551 : endif
552 46 : if (tr_pond .and. .not. tr_pond_lvl) &
553 7 : write(nu_diag,1000) ' pndaspect = ', pndaspect
554 :
555 46 : write(nu_diag,1020) ' ktherm = ', ktherm
556 46 : if (ktherm == 1) &
557 9 : write(nu_diag,1030) ' conduct = ', conduct
558 46 : write(nu_diag,1005) ' emissivity = ', emissivity
559 46 : if (ktherm == 2) then
560 31 : write(nu_diag,1005) ' a_rapid_mode = ', a_rapid_mode
561 31 : write(nu_diag,1005) ' Rac_rapid_mode = ', Rac_rapid_mode
562 31 : write(nu_diag,1005) ' aspect_rapid_mode = ', aspect_rapid_mode
563 31 : write(nu_diag,1005) ' dSdt_slow_mode = ', dSdt_slow_mode
564 31 : write(nu_diag,1005) ' phi_c_slow_mode = ', phi_c_slow_mode
565 31 : write(nu_diag,1005) ' phi_i_mushy = ', phi_i_mushy
566 : endif
567 46 : write(nu_diag,1010) ' sw_redist = ', sw_redist
568 46 : write(nu_diag,1005) ' sw_frac = ', sw_frac
569 46 : write(nu_diag,1005) ' sw_dtemp = ', sw_dtemp
570 :
571 46 : write(nu_diag,1030) ' atmbndy = ', &
572 92 : trim(atmbndy)
573 46 : write(nu_diag,1010) ' formdrag = ', formdrag
574 46 : write(nu_diag,1010) ' highfreq = ', highfreq
575 46 : write(nu_diag,1020) ' natmiter = ', natmiter
576 46 : write(nu_diag,1005) ' atmiter_conv = ', atmiter_conv
577 46 : write(nu_diag,1010) ' calc_strair = ', calc_strair
578 46 : write(nu_diag,1010) ' calc_Tsfc = ', calc_Tsfc
579 :
580 46 : write(nu_diag,*) ' atm_data_type = ', &
581 92 : trim(atm_data_type)
582 46 : write(nu_diag,*) ' ocn_data_type = ', &
583 92 : trim(ocn_data_type)
584 46 : write(nu_diag,*) ' bgc_data_type = ', &
585 92 : trim(bgc_data_type)
586 :
587 46 : write(nu_diag,*) ' atm_data_file = ', &
588 92 : trim(atm_data_file)
589 46 : write(nu_diag,*) ' ocn_data_file = ', &
590 92 : trim(ocn_data_file)
591 46 : write(nu_diag,*) ' bgc_data_file = ', &
592 92 : trim(bgc_data_file)
593 46 : write(nu_diag,*) ' ice_data_file = ', &
594 92 : trim(ice_data_file)
595 :
596 46 : if (trim(atm_data_type)=='default') &
597 9 : write(nu_diag,*) ' default_season = ', trim(default_season)
598 :
599 46 : write(nu_diag,1010) ' update_ocn_f = ', update_ocn_f
600 46 : write(nu_diag,1010) ' wave_spec = ', wave_spec
601 46 : if (wave_spec) &
602 3 : write(nu_diag,*) ' wave_spec_type = ', wave_spec_type
603 46 : write(nu_diag,1010) ' l_mpond_fresh = ', l_mpond_fresh
604 46 : write(nu_diag,1005) ' ustar_min = ', ustar_min
605 46 : write(nu_diag,*) ' fbot_xfer_type = ', &
606 92 : trim(fbot_xfer_type)
607 46 : write(nu_diag,1010) ' oceanmixed_ice = ', oceanmixed_ice
608 46 : write(nu_diag,*) ' tfrz_option = ', &
609 92 : trim(tfrz_option)
610 46 : write(nu_diag,1010) ' restore_ocn = ', restore_ocn
611 46 : if (restore_ocn) &
612 11 : write(nu_diag,1005) ' trestore = ', trestore
613 :
614 : ! tracers
615 46 : write(nu_diag,1010) ' tr_iage = ', tr_iage
616 46 : write(nu_diag,1010) ' tr_FY = ', tr_FY
617 46 : write(nu_diag,1010) ' tr_lvl = ', tr_lvl
618 46 : write(nu_diag,1010) ' tr_pond_cesm = ', tr_pond_cesm
619 46 : write(nu_diag,1010) ' tr_pond_lvl = ', tr_pond_lvl
620 46 : write(nu_diag,1010) ' tr_pond_topo = ', tr_pond_topo
621 46 : write(nu_diag,1010) ' tr_aero = ', tr_aero
622 46 : write(nu_diag,1010) ' tr_fsd = ', tr_fsd
623 :
624 46 : nt_Tsfc = 1 ! index tracers, starting with Tsfc = 1
625 46 : ntrcr = 1 ! count tracers, starting with Tsfc = 1
626 :
627 46 : nt_qice = ntrcr + 1
628 46 : ntrcr = ntrcr + nilyr ! qice in nilyr layers
629 :
630 46 : nt_qsno = ntrcr + 1
631 46 : ntrcr = ntrcr + nslyr ! qsno in nslyr layers
632 :
633 46 : nt_sice = ntrcr + 1
634 46 : ntrcr = ntrcr + nilyr ! sice in nilyr layers
635 :
636 46 : nt_iage = max_ntrcr
637 46 : if (tr_iage) then
638 2 : ntrcr = ntrcr + 1
639 2 : nt_iage = ntrcr ! chronological ice age
640 : endif
641 :
642 46 : nt_FY = max_ntrcr
643 46 : if (tr_FY) then
644 2 : ntrcr = ntrcr + 1
645 2 : nt_FY = ntrcr ! area of first year ice
646 : endif
647 :
648 46 : nt_alvl = max_ntrcr
649 46 : nt_vlvl = max_ntrcr
650 46 : if (tr_lvl) then
651 42 : ntrcr = ntrcr + 1
652 42 : nt_alvl = ntrcr ! area of level ice
653 42 : ntrcr = ntrcr + 1
654 42 : nt_vlvl = ntrcr ! volume of level ice
655 : endif
656 :
657 46 : nt_apnd = max_ntrcr
658 46 : nt_hpnd = max_ntrcr
659 46 : nt_ipnd = max_ntrcr
660 46 : if (tr_pond) then ! all explicit melt pond schemes
661 37 : ntrcr = ntrcr + 1
662 37 : nt_apnd = ntrcr
663 37 : ntrcr = ntrcr + 1
664 37 : nt_hpnd = ntrcr
665 37 : if (tr_pond_lvl) then
666 30 : ntrcr = ntrcr + 1 ! refrozen pond ice lid thickness
667 30 : nt_ipnd = ntrcr ! on level-ice ponds (if frzpnd='hlid')
668 : endif
669 37 : if (tr_pond_topo) then
670 5 : ntrcr = ntrcr + 1 !
671 5 : nt_ipnd = ntrcr ! refrozen pond ice lid thickness
672 : endif
673 : endif
674 :
675 46 : nt_fsd = max_ntrcr
676 46 : if (tr_fsd) then
677 4 : nt_fsd = ntrcr + 1 ! floe size distribution
678 4 : ntrcr = ntrcr + nfsd
679 : end if
680 :
681 46 : nt_isosno = max_ntrcr
682 46 : nt_isoice = max_ntrcr
683 46 : if (tr_iso) then ! isotopes
684 2 : nt_isosno = ntrcr + 1
685 2 : ntrcr = ntrcr + n_iso ! n_iso species in snow
686 2 : nt_isoice = ntrcr + 1
687 2 : ntrcr = ntrcr + n_iso ! n_iso species in ice
688 : end if
689 :
690 46 : nt_aero = max_ntrcr - 4*n_aero
691 46 : if (tr_aero) then
692 2 : nt_aero = ntrcr + 1
693 2 : ntrcr = ntrcr + 4*n_aero ! 4 dEdd layers, n_aero species
694 : endif
695 :
696 46 : if (ntrcr > max_ntrcr-1) then
697 0 : write(nu_diag,*) 'max_ntrcr-1 < number of namelist tracers'
698 0 : write(nu_diag,*) 'max_ntrcr-1 = ',max_ntrcr-1,' ntrcr = ',ntrcr
699 0 : call icedrv_system_abort(file=__FILE__,line=__LINE__)
700 : endif
701 :
702 46 : write(nu_diag,*) ' '
703 46 : write(nu_diag,1020) 'max_ntrcr = ', max_ntrcr
704 46 : write(nu_diag,1020) 'ntrcr = ', ntrcr
705 46 : write(nu_diag,*) ' '
706 46 : write(nu_diag,1020) 'nt_sice = ', nt_sice
707 46 : write(nu_diag,1020) 'nt_qice = ', nt_qice
708 46 : write(nu_diag,1020) 'nt_qsno = ', nt_qsno
709 46 : write(nu_diag,*)' '
710 46 : write(nu_diag,1020) 'ncat = ', ncat
711 46 : write(nu_diag,1020) 'nilyr = ', nilyr
712 46 : write(nu_diag,1020) 'nslyr = ', nslyr
713 46 : write(nu_diag,1020) 'nblyr = ', nblyr
714 46 : write(nu_diag,1020) 'nfsd = ', nfsd
715 46 : write(nu_diag,1020) 'n_iso = ', n_iso
716 46 : write(nu_diag,1020) 'n_aero = ', n_aero
717 46 : write(nu_diag,*)' '
718 46 : write(nu_diag,1020) 'nx = ', nx
719 46 : write(nu_diag,*)' '
720 :
721 : 1000 format (a30,2x,f9.2) ! a30 to align formatted, unformatted statements
722 : 1005 format (a30,2x,f10.6) ! float
723 : 1010 format (a30,2x,l6) ! logical
724 : 1020 format (a30,2x,i6) ! integer
725 : 1030 format (a30, a8) ! character
726 : 1040 format (a30,2x,6i6) ! integer
727 : 1050 format (a30,2x,6a6) ! character
728 :
729 46 : if (formdrag) then
730 3 : if (nt_apnd==0) then
731 0 : write(nu_diag,*)'ERROR: nt_apnd:',nt_apnd
732 0 : call icedrv_system_abort(file=__FILE__,line=__LINE__)
733 3 : elseif (nt_hpnd==0) then
734 0 : write(nu_diag,*)'ERROR: nt_hpnd:',nt_hpnd
735 0 : call icedrv_system_abort(file=__FILE__,line=__LINE__)
736 3 : elseif (nt_ipnd==0) then
737 0 : write(nu_diag,*)'ERROR: nt_ipnd:',nt_ipnd
738 0 : call icedrv_system_abort(file=__FILE__,line=__LINE__)
739 3 : elseif (nt_alvl==0) then
740 0 : write(nu_diag,*)'ERROR: nt_alvl:',nt_alvl
741 0 : call icedrv_system_abort(file=__FILE__,line=__LINE__)
742 3 : elseif (nt_vlvl==0) then
743 0 : write(nu_diag,*)'ERROR: nt_vlvl:',nt_vlvl
744 0 : call icedrv_system_abort(file=__FILE__,line=__LINE__)
745 : endif
746 : endif
747 :
748 : !-----------------------------------------------------------------
749 : ! set Icepack values
750 : !-----------------------------------------------------------------
751 :
752 : call icepack_init_parameters(ustar_min_in=ustar_min, Cf_in=Cf, &
753 : albicev_in=albicev, albicei_in=albicei, ksno_in=ksno, &
754 : albsnowv_in=albsnowv, albsnowi_in=albsnowi, &
755 : natmiter_in=natmiter, ahmax_in=ahmax, shortwave_in=shortwave, &
756 : atmiter_conv_in = atmiter_conv, &
757 : albedo_type_in=albedo_type, R_ice_in=R_ice, R_pnd_in=R_pnd, &
758 : R_snw_in=R_snw, dT_mlt_in=dT_mlt, rsnw_mlt_in=rsnw_mlt, &
759 : kstrength_in=kstrength, krdg_partic_in=krdg_partic, &
760 : krdg_redist_in=krdg_redist, mu_rdg_in=mu_rdg, &
761 : atmbndy_in=atmbndy, calc_strair_in=calc_strair, &
762 : formdrag_in=formdrag, highfreq_in=highfreq, &
763 : emissivity_in=emissivity, &
764 : kitd_in=kitd, kcatbound_in=kcatbound, hs0_in=hs0, &
765 : dpscale_in=dpscale, frzpnd_in=frzpnd, &
766 : rfracmin_in=rfracmin, rfracmax_in=rfracmax, &
767 : pndaspect_in=pndaspect, hs1_in=hs1, hp1_in=hp1, &
768 : ktherm_in=ktherm, calc_Tsfc_in=calc_Tsfc, &
769 : conduct_in=conduct, a_rapid_mode_in=a_rapid_mode, &
770 : Rac_rapid_mode_in=Rac_rapid_mode, &
771 : aspect_rapid_mode_in=aspect_rapid_mode, &
772 : dSdt_slow_mode_in=dSdt_slow_mode, &
773 : phi_c_slow_mode_in=phi_c_slow_mode, &
774 : phi_i_mushy_in=phi_i_mushy, conserv_check_in=conserv_check, &
775 : tfrz_option_in=tfrz_option, kalg_in=kalg, &
776 : fbot_xfer_type_in=fbot_xfer_type, &
777 : wave_spec_type_in=wave_spec_type, wave_spec_in=wave_spec, &
778 46 : sw_redist_in=sw_redist, sw_frac_in=sw_frac, sw_dtemp_in=sw_dtemp)
779 : call icepack_init_tracer_sizes(ntrcr_in=ntrcr, &
780 : ncat_in=ncat, nilyr_in=nilyr, nslyr_in=nslyr, nblyr_in=nblyr, &
781 46 : nfsd_in=nfsd, n_iso_in=n_iso, n_aero_in=n_aero)
782 : call icepack_init_tracer_flags(tr_iage_in=tr_iage, &
783 : tr_FY_in=tr_FY, tr_lvl_in=tr_lvl, tr_aero_in=tr_aero, &
784 : tr_iso_in=tr_iso, &
785 : tr_pond_in=tr_pond, tr_pond_cesm_in=tr_pond_cesm, &
786 : tr_pond_lvl_in=tr_pond_lvl, &
787 46 : tr_pond_topo_in=tr_pond_topo, tr_fsd_in=tr_fsd)
788 : call icepack_init_tracer_indices(nt_Tsfc_in=nt_Tsfc, &
789 : nt_sice_in=nt_sice, nt_qice_in=nt_qice, &
790 : nt_qsno_in=nt_qsno, nt_iage_in=nt_iage, &
791 : nt_fy_in=nt_fy, nt_alvl_in=nt_alvl, nt_vlvl_in=nt_vlvl, &
792 : nt_apnd_in=nt_apnd, nt_hpnd_in=nt_hpnd, nt_ipnd_in=nt_ipnd, &
793 : nt_aero_in=nt_aero, nt_fsd_in=nt_fsd, &
794 46 : nt_isosno_in=nt_isosno, nt_isoice_in=nt_isoice)
795 :
796 46 : call icepack_warnings_flush(nu_diag)
797 46 : if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
798 0 : file=__FILE__,line= __LINE__)
799 :
800 46 : end subroutine input_data
801 :
802 : !=======================================================================
803 :
804 : ! Horizontal grid initialization:
805 : !
806 : ! author: Elizabeth C. Hunke, LANL
807 :
808 46 : subroutine init_grid2
809 :
810 : integer :: i
811 17 : real (kind=dbl_kind) :: pi, puny
812 : character(len=*), parameter :: subname='(init_grid2)'
813 :
814 : !-----------------------------------------------------------------
815 : ! query Icepack values
816 : !-----------------------------------------------------------------
817 :
818 46 : call icepack_query_parameters(pi_out=pi,puny_out=puny)
819 46 : call icepack_warnings_flush(nu_diag)
820 46 : if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
821 0 : file=__FILE__, line=__LINE__)
822 :
823 : !-----------------------------------------------------------------
824 : ! lat, lon, cell widths, angle, land mask
825 : !-----------------------------------------------------------------
826 :
827 230 : TLAT(:) = p5*pi ! pi/2, North pole
828 46 : TLON(:) = c0
829 :
830 184 : do i = 2, nx
831 184 : TLAT(i) = TLAT(i-1) - p5*pi/180._dbl_kind ! half-deg increments
832 : enddo
833 :
834 230 : tmask(:) = .true.
835 46 : tmask(nx) = .false. ! land in last grid cell
836 :
837 : !-----------------------------------------------------------------
838 : ! create hemisphere masks
839 : !-----------------------------------------------------------------
840 :
841 46 : lmask_n(:) = .false.
842 46 : lmask_s(:) = .false.
843 :
844 230 : do i = 1, nx
845 184 : if (TLAT(i) >= -puny) lmask_n(i) = .true. ! N. Hem.
846 230 : if (TLAT(i) < -puny) lmask_s(i) = .true. ! S. Hem.
847 : enddo
848 :
849 46 : end subroutine init_grid2
850 :
851 : !=======================================================================
852 :
853 : ! Initialize state for the itd model
854 : !
855 : ! authors: C. M. Bitz, UW
856 : ! William H. Lipscomb, LANL
857 :
858 46 : subroutine init_state
859 :
860 : use icepack_intfc, only: icepack_aggregate
861 : use icedrv_domain_size, only: ncat, nilyr, nslyr, nblyr, max_ntrcr
862 : use icedrv_domain_size, only: n_iso, n_aero, nfsd
863 : use icedrv_flux, only: sst, Tf, Tair, salinz, Tmltz
864 : use icedrv_state, only: trcr_depend, aicen, trcrn, vicen, vsnon
865 : use icedrv_state, only: aice0, aice, vice, vsno, trcr, aice_init
866 : use icedrv_state, only: n_trcr_strata, nt_strata, trcr_base
867 :
868 : integer (kind=int_kind) :: &
869 : i , & ! horizontal indes
870 : k , & ! vertical index
871 : it ! tracer index
872 :
873 : logical (kind=log_kind) :: &
874 : heat_capacity ! from icepack
875 :
876 : integer (kind=int_kind) :: ntrcr
877 : logical (kind=log_kind) :: tr_iage, tr_FY, tr_lvl, tr_aero, tr_fsd, tr_iso
878 : logical (kind=log_kind) :: tr_pond_cesm, tr_pond_lvl, tr_pond_topo
879 : integer (kind=int_kind) :: nt_Tsfc, nt_sice, nt_qice, nt_qsno, nt_iage, nt_fy
880 : integer (kind=int_kind) :: nt_alvl, nt_vlvl, nt_apnd, nt_hpnd, &
881 : nt_ipnd, nt_aero, nt_fsd, nt_isosno, nt_isoice
882 :
883 : character(len=*), parameter :: subname='(init_state)'
884 :
885 : !-----------------------------------------------------------------
886 : ! query Icepack values
887 : !-----------------------------------------------------------------
888 :
889 46 : call icepack_query_parameters(heat_capacity_out=heat_capacity)
890 46 : call icepack_query_tracer_sizes(ntrcr_out=ntrcr)
891 : call icepack_query_tracer_flags(tr_iage_out=tr_iage, &
892 : tr_FY_out=tr_FY, tr_lvl_out=tr_lvl, tr_aero_out=tr_aero, &
893 : tr_iso_out=tr_iso, &
894 : tr_pond_cesm_out=tr_pond_cesm, tr_pond_lvl_out=tr_pond_lvl, &
895 46 : tr_pond_topo_out=tr_pond_topo, tr_fsd_out=tr_fsd)
896 : call icepack_query_tracer_indices(nt_Tsfc_out=nt_Tsfc, &
897 : nt_sice_out=nt_sice, nt_qice_out=nt_qice, &
898 : nt_qsno_out=nt_qsno, nt_iage_out=nt_iage, nt_fy_out=nt_fy, &
899 : nt_alvl_out=nt_alvl, nt_vlvl_out=nt_vlvl, &
900 : nt_apnd_out=nt_apnd, nt_hpnd_out=nt_hpnd, &
901 : nt_ipnd_out=nt_ipnd, &
902 : nt_isosno_out=nt_isosno, nt_isoice_out=nt_isoice, &
903 46 : nt_aero_out=nt_aero, nt_fsd_out=nt_fsd)
904 46 : call icepack_warnings_flush(nu_diag)
905 46 : if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
906 0 : file=__FILE__,line= __LINE__)
907 :
908 : !-----------------------------------------------------------------
909 : ! Check number of layers in ice and snow.
910 : !-----------------------------------------------------------------
911 : if (nilyr < 1) then
912 : write (nu_diag,*) 'nilyr =', nilyr
913 : write (nu_diag,*) 'Must have at least one ice layer'
914 : call icedrv_system_abort(file=__FILE__,line=__LINE__)
915 : endif
916 :
917 : if (nslyr < 1) then
918 : write (nu_diag,*) 'nslyr =', nslyr
919 : write (nu_diag,*) 'Must have at least one snow layer'
920 : call icedrv_system_abort(file=__FILE__,line=__LINE__)
921 : endif
922 :
923 46 : if (.not.heat_capacity) then
924 :
925 6 : write (nu_diag,*) 'WARNING - Zero-layer thermodynamics'
926 :
927 : if (nilyr > 1) then
928 0 : write (nu_diag,*) 'nilyr =', nilyr
929 : write (nu_diag,*) &
930 0 : 'Must have nilyr = 1 if ktherm = 0'
931 0 : call icedrv_system_abort(file=__FILE__,line=__LINE__)
932 : endif
933 :
934 : if (nslyr > 1) then
935 0 : write (nu_diag,*) 'nslyr =', nslyr
936 : write (nu_diag,*) &
937 0 : 'Must have nslyr = 1 if heat_capacity = F'
938 0 : call icedrv_system_abort(file=__FILE__,line=__LINE__)
939 : endif
940 :
941 : endif ! heat_capacity = F
942 :
943 : !-----------------------------------------------------------------
944 : ! Set tracer types
945 : !-----------------------------------------------------------------
946 :
947 46 : trcr_depend(nt_Tsfc) = 0 ! ice/snow surface temperature
948 332 : do k = 1, nilyr
949 286 : trcr_depend(nt_sice + k - 1) = 1 ! volume-weighted ice salinity
950 332 : trcr_depend(nt_qice + k - 1) = 1 ! volume-weighted ice enthalpy
951 : enddo
952 104 : do k = 1, nslyr
953 104 : trcr_depend(nt_qsno + k - 1) = 2 ! volume-weighted snow enthalpy
954 : enddo
955 46 : if (tr_iage) trcr_depend(nt_iage) = 1 ! volume-weighted ice age
956 46 : if (tr_FY) trcr_depend(nt_FY) = 0 ! area-weighted first-year ice area
957 46 : if (tr_lvl) trcr_depend(nt_alvl) = 0 ! level ice area
958 46 : if (tr_lvl) trcr_depend(nt_vlvl) = 1 ! level ice volume
959 46 : if (tr_pond_cesm) then
960 2 : trcr_depend(nt_apnd) = 0 ! melt pond area
961 2 : trcr_depend(nt_hpnd) = 2+nt_apnd ! melt pond depth
962 : endif
963 46 : if (tr_pond_lvl) then
964 30 : trcr_depend(nt_apnd) = 2+nt_alvl ! melt pond area
965 30 : trcr_depend(nt_hpnd) = 2+nt_apnd ! melt pond depth
966 30 : trcr_depend(nt_ipnd) = 2+nt_apnd ! refrozen pond lid
967 : endif
968 46 : if (tr_pond_topo) then
969 5 : trcr_depend(nt_apnd) = 0 ! melt pond area
970 5 : trcr_depend(nt_hpnd) = 2+nt_apnd ! melt pond depth
971 5 : trcr_depend(nt_ipnd) = 2+nt_apnd ! refrozen pond lid
972 : endif
973 46 : if (tr_fsd) then
974 41 : do it = 1, nfsd
975 41 : trcr_depend(nt_fsd + it - 1) = 0 ! area-weighted floe size distribution
976 : enddo
977 : endif
978 46 : if (tr_iso) then ! isotopes
979 8 : do it = 1, n_iso
980 6 : trcr_depend(nt_isosno) = 2 ! snow
981 8 : trcr_depend(nt_isoice) = 1 ! ice
982 : enddo
983 : endif
984 46 : if (tr_aero) then ! volume-weighted aerosols
985 4 : do it = 1, n_aero
986 2 : trcr_depend(nt_aero+(it-1)*4 ) = 2 ! snow
987 2 : trcr_depend(nt_aero+(it-1)*4+1) = 2 ! snow
988 2 : trcr_depend(nt_aero+(it-1)*4+2) = 1 ! ice
989 4 : trcr_depend(nt_aero+(it-1)*4+3) = 1 ! ice
990 : enddo
991 : endif
992 :
993 1833 : do it = 1, ntrcr
994 : ! mask for base quantity on which tracers are carried
995 1787 : if (trcr_depend(it) == 0) then ! area
996 158 : trcr_base(it,1) = c1
997 1629 : elseif (trcr_depend(it) == 1) then ! ice volume
998 627 : trcr_base(it,2) = c1
999 1002 : elseif (trcr_depend(it) == 2) then ! snow volume
1000 216 : trcr_base(it,3) = c1
1001 : else
1002 786 : trcr_base(it,1) = c1 ! default: ice area
1003 786 : trcr_base(it,2) = c0
1004 786 : trcr_base(it,3) = c0
1005 : endif
1006 :
1007 : ! initialize number of underlying tracer layers
1008 1787 : n_trcr_strata(it) = 0
1009 : ! default indices of underlying tracer layers
1010 1787 : nt_strata (it,1) = 0
1011 1833 : nt_strata (it,2) = 0
1012 : enddo
1013 :
1014 46 : if (tr_pond_cesm) then
1015 2 : n_trcr_strata(nt_hpnd) = 1 ! melt pond depth
1016 2 : nt_strata (nt_hpnd,1) = nt_apnd ! on melt pond area
1017 : endif
1018 46 : if (tr_pond_lvl) then
1019 30 : n_trcr_strata(nt_apnd) = 1 ! melt pond area
1020 30 : nt_strata (nt_apnd,1) = nt_alvl ! on level ice area
1021 30 : n_trcr_strata(nt_hpnd) = 2 ! melt pond depth
1022 30 : nt_strata (nt_hpnd,2) = nt_apnd ! on melt pond area
1023 30 : nt_strata (nt_hpnd,1) = nt_alvl ! on level ice area
1024 30 : n_trcr_strata(nt_ipnd) = 2 ! refrozen pond lid
1025 30 : nt_strata (nt_ipnd,2) = nt_apnd ! on melt pond area
1026 30 : nt_strata (nt_ipnd,1) = nt_alvl ! on level ice area
1027 : endif
1028 46 : if (tr_pond_topo) then
1029 5 : n_trcr_strata(nt_hpnd) = 1 ! melt pond depth
1030 5 : nt_strata (nt_hpnd,1) = nt_apnd ! on melt pond area
1031 5 : n_trcr_strata(nt_ipnd) = 1 ! refrozen pond lid
1032 5 : nt_strata (nt_ipnd,1) = nt_apnd ! on melt pond area
1033 : endif
1034 :
1035 : !-----------------------------------------------------------------
1036 : ! Set state variables
1037 : !-----------------------------------------------------------------
1038 :
1039 : call set_state_var (nx, &
1040 : Tair (:), sst (:), &
1041 : Tf (:), &
1042 : salinz(:,:), Tmltz(:,:), &
1043 : aicen (:,:), trcrn(:,:,:), &
1044 46 : vicen (:,:), vsnon(:,:))
1045 :
1046 : !-----------------------------------------------------------------
1047 : ! compute aggregate ice state and open water area
1048 : !-----------------------------------------------------------------
1049 :
1050 : !$OMP PARALLEL DO PRIVATE(i)
1051 230 : do i = 1, nx
1052 184 : aice(i) = c0
1053 184 : vice(i) = c0
1054 184 : vsno(i) = c0
1055 9612 : do it = 1, max_ntrcr
1056 9612 : trcr(i,it) = c0
1057 : enddo
1058 :
1059 184 : if (tmask(i)) &
1060 : call icepack_aggregate(ncat=ncat, &
1061 0 : trcrn=trcrn(i,1:ntrcr,:), &
1062 : aicen=aicen(i,:), &
1063 : vicen=vicen(i,:), &
1064 : vsnon=vsnon(i,:), &
1065 0 : trcr=trcr (i,1:ntrcr), &
1066 51 : aice=aice (i), &
1067 51 : vice=vice (i), &
1068 51 : vsno=vsno (i), &
1069 51 : aice0=aice0(i), &
1070 : ntrcr=ntrcr, &
1071 0 : trcr_depend=trcr_depend(1:ntrcr), &
1072 0 : trcr_base=trcr_base (1:ntrcr,:), &
1073 0 : n_trcr_strata=n_trcr_strata(1:ntrcr), &
1074 240 : nt_strata=nt_strata (1:ntrcr,:))
1075 :
1076 230 : aice_init(i) = aice(i)
1077 :
1078 : enddo
1079 : !$OMP END PARALLEL DO
1080 :
1081 46 : call icepack_warnings_flush(nu_diag)
1082 46 : if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
1083 0 : file=__FILE__, line=__LINE__)
1084 :
1085 46 : end subroutine init_state
1086 :
1087 : !=======================================================================
1088 :
1089 : ! Initialize state in each ice thickness category
1090 : !
1091 : ! authors: Elizabeth Hunke, LANL
1092 :
1093 46 : subroutine set_state_var (nx, &
1094 92 : Tair, sst, &
1095 46 : Tf, &
1096 92 : salinz, Tmltz, &
1097 92 : aicen, trcrn, &
1098 46 : vicen, vsnon)
1099 :
1100 : use icedrv_arrays_column, only: hin_max
1101 : use icedrv_domain_size, only: nilyr, nslyr, max_ntrcr, ncat, nfsd
1102 : use icedrv_arrays_column, only: floe_rad_c, floe_binwidth
1103 :
1104 :
1105 : integer (kind=int_kind), intent(in) :: &
1106 : nx ! number of grid cells
1107 :
1108 : real (kind=dbl_kind), dimension (nx), intent(in) :: &
1109 : Tair ! air temperature (K)
1110 :
1111 : ! ocean values may be redefined here, unlike in CICE
1112 : real (kind=dbl_kind), dimension (nx), intent(inout) :: &
1113 : Tf , & ! freezing temperature (C)
1114 : sst ! sea surface temperature (C)
1115 :
1116 : real (kind=dbl_kind), dimension (nx,nilyr), &
1117 : intent(in) :: &
1118 : salinz , & ! initial salinity profile
1119 : Tmltz ! initial melting temperature profile
1120 :
1121 : real (kind=dbl_kind), dimension (nx,ncat), &
1122 : intent(out) :: &
1123 : aicen , & ! concentration of ice
1124 : vicen , & ! volume per unit area of ice (m)
1125 : vsnon ! volume per unit area of snow (m)
1126 :
1127 : real (kind=dbl_kind), dimension (nx,max_ntrcr,ncat), &
1128 : intent(out) :: &
1129 : trcrn ! ice tracers
1130 : ! 1: surface temperature of ice/snow (C)
1131 :
1132 : ! local variables
1133 :
1134 : integer (kind=int_kind) :: &
1135 : i , & ! horizontal indices
1136 : k , & ! ice layer index
1137 : n , & ! thickness category index
1138 : it ! tracer index
1139 :
1140 : real (kind=dbl_kind) :: &
1141 34 : Tsfc, sum, hbar, &
1142 34 : rhos, Lfresh, puny
1143 :
1144 : real (kind=dbl_kind), dimension(ncat) :: &
1145 179 : ainit, hinit ! initial area, thickness
1146 :
1147 : real (kind=dbl_kind), dimension(nilyr) :: &
1148 124 : qin ! ice enthalpy (J/m3)
1149 :
1150 : real (kind=dbl_kind), dimension(nslyr) :: &
1151 38 : qsn ! snow enthalpy (J/m3)
1152 :
1153 : real (kind=dbl_kind), parameter :: &
1154 : hsno_init = 0.25_dbl_kind ! initial snow thickness (m)
1155 :
1156 : logical (kind=log_kind) :: tr_brine, tr_lvl, tr_fsd
1157 : integer (kind=int_kind) :: nt_Tsfc, nt_qice, nt_qsno, nt_sice, nt_fsd
1158 : integer (kind=int_kind) :: nt_fbri, nt_alvl, nt_vlvl
1159 :
1160 : character(len=*), parameter :: subname='(set_state_var)'
1161 :
1162 : !-----------------------------------------------------------------
1163 : ! query Icepack values
1164 : !-----------------------------------------------------------------
1165 :
1166 : call icepack_query_tracer_flags(tr_brine_out=tr_brine, tr_lvl_out=tr_lvl, &
1167 46 : tr_fsd_out=tr_fsd)
1168 : call icepack_query_tracer_indices( nt_Tsfc_out=nt_Tsfc, nt_qice_out=nt_qice, &
1169 : nt_qsno_out=nt_qsno, nt_sice_out=nt_sice, nt_fsd_out=nt_fsd, &
1170 46 : nt_fbri_out=nt_fbri, nt_alvl_out=nt_alvl, nt_vlvl_out=nt_vlvl)
1171 46 : call icepack_query_parameters(rhos_out=rhos, Lfresh_out=Lfresh, puny_out=puny)
1172 46 : call icepack_warnings_flush(nu_diag)
1173 46 : if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
1174 0 : file=__FILE__,line= __LINE__)
1175 :
1176 : !-----------------------------------------------------------------
1177 : ! Initialize state variables.
1178 : ! If restarting, these values are overwritten.
1179 : !-----------------------------------------------------------------
1180 :
1181 264 : do n = 1, ncat
1182 1090 : do i = 1, nx
1183 872 : aicen(i,n) = c0
1184 872 : vicen(i,n) = c0
1185 872 : vsnon(i,n) = c0
1186 872 : trcrn(i,nt_Tsfc,n) = Tf(i) ! surface temperature
1187 : if (max_ntrcr >= 2) then
1188 46324 : do it = 2, max_ntrcr
1189 46324 : trcrn(i,it,n) = c0
1190 : enddo
1191 : endif
1192 872 : if (tr_lvl) trcrn(i,nt_alvl,n) = c1
1193 872 : if (tr_lvl) trcrn(i,nt_vlvl,n) = c1
1194 872 : if (tr_brine) trcrn(i,nt_fbri,n) = c1
1195 6544 : do k = 1, nilyr
1196 6544 : trcrn(i,nt_sice+k-1,n) = salinz(i,k)
1197 : enddo
1198 2202 : do k = 1, nslyr
1199 1984 : trcrn(i,nt_qsno+k-1,n) = -rhos * Lfresh
1200 : enddo
1201 : enddo
1202 218 : ainit(n) = c0
1203 264 : hinit(n) = c0
1204 : enddo
1205 :
1206 : !-----------------------------------------------------------------
1207 : ! For Icepack testing, the grid vector is populated with several
1208 : ! different ice distributions, including ice-free, a single-
1209 : ! thickness slab, a full thickness distribution (as in CICE),
1210 : ! and land
1211 : !-----------------------------------------------------------------
1212 :
1213 46 : i = 1 ! ice-free
1214 : ! already initialized above
1215 :
1216 : !-----------------------------------------------------------------
1217 :
1218 46 : i = 2 ! 2-m slab, no snow
1219 : if (3 <= ncat) then
1220 43 : n = 3
1221 43 : ainit(n) = c1 ! assumes we are using the default ITD boundaries
1222 43 : hinit(n) = c2
1223 : else
1224 3 : ainit(ncat) = c1
1225 3 : hinit(ncat) = c2
1226 : endif
1227 264 : do n = 1, ncat
1228 : ! ice volume, snow volume
1229 218 : aicen(i,n) = ainit(n)
1230 218 : vicen(i,n) = hinit(n) * ainit(n) ! m
1231 218 : vsnon(i,n) = c0
1232 : ! tracers
1233 81 : call icepack_init_trcr(Tair = Tair(i), &
1234 81 : Tf = Tf(i), &
1235 : Sprofile = salinz(i,:), &
1236 : Tprofile = Tmltz(i,:), &
1237 : Tsfc = Tsfc, &
1238 : nilyr=nilyr, nslyr=nslyr, &
1239 218 : qin=qin(:), qsn=qsn(:))
1240 :
1241 : ! floe size distribution
1242 218 : if (tr_fsd) call icepack_init_fsd(nfsd=nfsd, ice_ic=ice_ic, &
1243 : floe_rad_c=floe_rad_c, &
1244 : floe_binwidth=floe_binwidth, &
1245 20 : afsd=trcrn(i,nt_fsd:nt_fsd+nfsd-1,n))
1246 : ! surface temperature
1247 218 : trcrn(i,nt_Tsfc,n) = Tsfc ! deg C
1248 : ! ice enthalpy, salinity
1249 1636 : do k = 1, nilyr
1250 1418 : trcrn(i,nt_qice+k-1,n) = qin(k)
1251 1636 : trcrn(i,nt_sice+k-1,n) = salinz(i,k)
1252 : enddo
1253 : ! snow enthalpy
1254 496 : do k = 1, nslyr
1255 496 : trcrn(i,nt_qsno+k-1,n) = qsn(k)
1256 : enddo ! nslyr
1257 : ! brine fraction
1258 264 : if (tr_brine) trcrn(i,nt_fbri,n) = c1
1259 : enddo ! ncat
1260 46 : call icepack_warnings_flush(nu_diag)
1261 46 : if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
1262 0 : file=__FILE__, line=__LINE__)
1263 :
1264 : !-----------------------------------------------------------------
1265 :
1266 46 : i = 3 ! full thickness distribution
1267 : ! initial category areas in cells with ice
1268 46 : hbar = c3 ! initial ice thickness with greatest area
1269 : ! Note: the resulting average ice thickness
1270 : ! tends to be less than hbar due to the
1271 : ! nonlinear distribution of ice thicknesses
1272 :
1273 46 : sum = c0
1274 264 : do n = 1, ncat
1275 218 : if (n < ncat) then
1276 172 : hinit(n) = p5*(hin_max(n-1) + hin_max(n)) ! m
1277 : else ! n=ncat
1278 46 : hinit(n) = (hin_max(n-1) + c1) ! m
1279 : endif
1280 : ! parabola, max at h=hbar, zero at h=0, 2*hbar
1281 218 : ainit(n) = max(c0, (c2*hbar*hinit(n) - hinit(n)**2))
1282 264 : sum = sum + ainit(n)
1283 : enddo
1284 264 : do n = 1, ncat
1285 264 : ainit(n) = ainit(n) / (sum + puny/ncat) ! normalize
1286 : enddo
1287 :
1288 264 : do n = 1, ncat
1289 : ! ice volume, snow volume
1290 218 : aicen(i,n) = ainit(n)
1291 218 : vicen(i,n) = hinit(n) * ainit(n) ! m
1292 218 : vsnon(i,n) = min(aicen(i,n)*hsno_init,p2*vicen(i,n))
1293 : ! tracers
1294 81 : call icepack_init_trcr(Tair = Tair(i), &
1295 81 : Tf = Tf(i), &
1296 : Sprofile = salinz(i,:), &
1297 : Tprofile = Tmltz(i,:), &
1298 : Tsfc = Tsfc, &
1299 : nilyr=nilyr, nslyr=nslyr, &
1300 218 : qin=qin(:), qsn=qsn(:))
1301 : ! floe size distribution
1302 218 : if (tr_fsd) call icepack_init_fsd(nfsd=nfsd, ice_ic=ice_ic, &
1303 : floe_rad_c=floe_rad_c, &
1304 : floe_binwidth=floe_binwidth, &
1305 20 : afsd=trcrn(i,nt_fsd:nt_fsd+nfsd-1,n))
1306 :
1307 : ! surface temperature
1308 218 : trcrn(i,nt_Tsfc,n) = Tsfc ! deg C
1309 : ! ice enthalpy, salinity
1310 1636 : do k = 1, nilyr
1311 1418 : trcrn(i,nt_qice+k-1,n) = qin(k)
1312 1636 : trcrn(i,nt_sice+k-1,n) = salinz(i,k)
1313 : enddo
1314 : ! snow enthalpy
1315 496 : do k = 1, nslyr
1316 496 : trcrn(i,nt_qsno+k-1,n) = qsn(k)
1317 : enddo ! nslyr
1318 : ! brine fraction
1319 264 : if (tr_brine) trcrn(i,nt_fbri,n) = c1
1320 : enddo ! ncat
1321 46 : call icepack_warnings_flush(nu_diag)
1322 46 : if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
1323 0 : file=__FILE__, line=__LINE__)
1324 :
1325 : !-----------------------------------------------------------------
1326 :
1327 : ! land
1328 : ! already initialized above (tmask = 0)
1329 46 : i = 4
1330 46 : sst(i) = c0
1331 46 : Tf(i) = c0
1332 :
1333 46 : end subroutine set_state_var
1334 :
1335 : !=======================================================================
1336 :
1337 : ! Initialize floe size distribution tracer (call prior to reading restart data)
1338 :
1339 46 : subroutine init_fsd
1340 :
1341 : use icedrv_arrays_column, only: wavefreq, dwavefreq, wave_sig_ht, &
1342 : wave_spectrum, d_afsd_newi, d_afsd_latg, d_afsd_latm, &
1343 : d_afsd_wave, d_afsd_weld
1344 :
1345 46 : wavefreq (:) = c0
1346 46 : dwavefreq (:) = c0
1347 46 : wave_sig_ht (:) = c0
1348 46 : wave_spectrum (:,:) = c0
1349 46 : d_afsd_newi (:,:) = c0
1350 46 : d_afsd_latg (:,:) = c0
1351 46 : d_afsd_latm (:,:) = c0
1352 46 : d_afsd_wave (:,:) = c0
1353 46 : d_afsd_weld (:,:) = c0
1354 :
1355 46 : end subroutine init_fsd
1356 :
1357 : !=======================================================================
1358 :
1359 : end module icedrv_init
1360 :
1361 : !=======================================================================
|
