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