Initialize temporally evolving fields, either as initial conditions or by reading them from a restart (or saves) file.
124 type(ocean_grid_type),
intent(inout) :: G
125 type(verticalGrid_type),
intent(in) :: GV
126 type(unit_scale_type),
intent(in) :: US
127 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
130 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
133 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
135 type(thermo_var_ptrs),
intent(inout) :: tv
137 type(time_type),
intent(inout) :: Time
138 type(param_file_type),
intent(in) :: PF
140 type(directories),
intent(in) :: dirs
142 type(MOM_restart_CS),
pointer :: restart_CS
144 type(ALE_CS),
pointer :: ALE_CSp
145 type(tracer_registry_type),
pointer :: tracer_Reg
146 type(sponge_CS),
pointer :: sponge_CSp
147 type(ALE_sponge_CS),
pointer :: ALE_sponge_CSp
148 type(ocean_OBC_type),
pointer :: OBC
149 type(time_type),
optional,
intent(in) :: Time_in
152 character(len=200) :: filename
153 character(len=200) :: filename2
154 character(len=200) :: inputdir
155 character(len=200) :: config
159 logical :: from_Z_file, useALE
161 integer :: write_geom
162 logical :: use_temperature, use_sponge, use_OBC
164 logical :: depress_sfc
166 logical :: trim_ic_for_p_surf
168 logical :: regrid_accelerate
169 integer :: regrid_iterations
175 type(EOS_type),
pointer :: eos => null()
178 logical :: debug_layers = .false.
179 character(len=80) :: mesg
181 #include "version_variable.h"
182 integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
183 integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
185 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
186 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
187 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
188 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
190 call calltree_enter(
"MOM_initialize_state(), MOM_state_initialization.F90")
191 call log_version(pf, mdl, version,
"")
192 call get_param(pf, mdl,
"DEBUG", debug, default=.false.)
193 call get_param(pf, mdl,
"DEBUG_OBC", debug_obc, default=.false.)
195 new_sim = determine_is_new_run(dirs%input_filename, dirs%restart_input_dir, &
197 just_read = .not.new_sim
199 call get_param(pf, mdl,
"INPUTDIR", inputdir, &
200 "The directory in which input files are found.", default=
".")
201 inputdir = slasher(inputdir)
203 use_temperature =
associated(tv%T)
204 useale =
associated(ale_csp)
205 use_eos =
associated(tv%eqn_of_state)
206 use_obc =
associated(obc)
207 if (use_eos) eos => tv%eqn_of_state
215 call mom_mesg(
"Run initialized internally.", 3)
217 if (
present(time_in)) time = time_in
232 call get_param(pf, mdl,
"INIT_LAYERS_FROM_Z_FILE", from_z_file, &
233 "If true, initialize the layer thicknesses, temperatures, "//&
234 "and salinities from a Z-space file on a latitude-longitude "//&
235 "grid.", default=.false., do_not_log=just_read)
237 if (from_z_file)
then
239 if (.NOT.use_temperature)
call mom_error(fatal,
"MOM_initialize_state : "//&
240 "use_temperature must be true if INIT_LAYERS_FROM_Z_FILE is true")
242 call mom_temp_salt_initialize_from_z(h, tv, g, gv, us, pf, just_read_params=just_read)
246 call get_param(pf, mdl,
"THICKNESS_CONFIG", config, &
247 "A string that determines how the initial layer "//&
248 "thicknesses are specified for a new run: \n"//&
249 " \t file - read interface heights from the file specified \n"//&
250 " \t thickness_file - read thicknesses from the file specified \n"//&
251 " \t\t by (THICKNESS_FILE).\n"//&
252 " \t coord - determined by ALE coordinate.\n"//&
253 " \t uniform - uniform thickness layers evenly distributed \n"//&
254 " \t\t between the surface and MAXIMUM_DEPTH. \n"//&
255 " \t list - read a list of positive interface depths. \n"//&
256 " \t DOME - use a slope and channel configuration for the \n"//&
257 " \t\t DOME sill-overflow test case. \n"//&
258 " \t ISOMIP - use a configuration for the \n"//&
259 " \t\t ISOMIP test case. \n"//&
260 " \t benchmark - use the benchmark test case thicknesses. \n"//&
261 " \t Neverland - use the Neverland test case thicknesses. \n"//&
262 " \t search - search a density profile for the interface \n"//&
263 " \t\t densities. This is not yet implemented. \n"//&
264 " \t circle_obcs - the circle_obcs test case is used. \n"//&
265 " \t DOME2D - 2D version of DOME initialization. \n"//&
266 " \t adjustment2d - 2D lock exchange thickness ICs. \n"//&
267 " \t sloshing - sloshing gravity thickness ICs. \n"//&
268 " \t seamount - no motion test with seamount ICs. \n"//&
269 " \t dumbbell - sloshing channel ICs. \n"//&
270 " \t soliton - Equatorial Rossby soliton. \n"//&
271 " \t rossby_front - a mixed layer front in thermal wind balance.\n"//&
272 " \t USER - call a user modified routine.", &
273 fail_if_missing=new_sim, do_not_log=just_read)
274 select case (trim(config))
276 call initialize_thickness_from_file(h, g, gv, us, pf, .false., just_read_params=just_read)
277 case (
"thickness_file")
278 call initialize_thickness_from_file(h, g, gv, us, pf, .true., just_read_params=just_read)
280 if (new_sim .and. useale)
then
281 call ale_initthicknesstocoord( ale_csp, g, gv, h )
282 elseif (new_sim)
then
283 call mom_error(fatal,
"MOM_initialize_state: USE_REGRIDDING must be True "//&
284 "for THICKNESS_CONFIG of 'coord'")
286 case (
"uniform");
call initialize_thickness_uniform(h, g, gv, pf, &
287 just_read_params=just_read)
288 case (
"list");
call initialize_thickness_list(h, g, gv, us, pf, &
289 just_read_params=just_read)
290 case (
"DOME");
call dome_initialize_thickness(h, g, gv, pf, &
291 just_read_params=just_read)
292 case (
"ISOMIP");
call isomip_initialize_thickness(h, g, gv, us, pf, tv, &
293 just_read_params=just_read)
294 case (
"benchmark");
call benchmark_initialize_thickness(h, g, gv, us, pf, &
295 tv%eqn_of_state, tv%P_Ref, just_read_params=just_read)
296 case (
"Neverland");
call neverland_initialize_thickness(h, g, gv, us, pf, &
297 tv%eqn_of_state, tv%P_Ref)
298 case (
"search");
call initialize_thickness_search
299 case (
"circle_obcs");
call circle_obcs_initialize_thickness(h, g, gv, pf, &
300 just_read_params=just_read)
301 case (
"lock_exchange");
call lock_exchange_initialize_thickness(h, g, gv, us, &
302 pf, just_read_params=just_read)
303 case (
"external_gwave");
call external_gwave_initialize_thickness(h, g, gv, us, &
304 pf, just_read_params=just_read)
305 case (
"DOME2D");
call dome2d_initialize_thickness(h, g, gv, us, pf, &
306 just_read_params=just_read)
307 case (
"adjustment2d");
call adjustment_initialize_thickness(h, g, gv, us, &
308 pf, just_read_params=just_read)
309 case (
"sloshing");
call sloshing_initialize_thickness(h, g, gv, us, pf, &
310 just_read_params=just_read)
311 case (
"seamount");
call seamount_initialize_thickness(h, g, gv, us, pf, &
312 just_read_params=just_read)
313 case (
"dumbbell");
call dumbbell_initialize_thickness(h, g, gv, us, pf, &
314 just_read_params=just_read)
315 case (
"soliton");
call soliton_initialize_thickness(h, g, gv, us)
316 case (
"phillips");
call phillips_initialize_thickness(h, g, gv, us, pf, &
317 just_read_params=just_read)
318 case (
"rossby_front");
call rossby_front_initialize_thickness(h, g, gv, &
319 pf, just_read_params=just_read)
320 case (
"USER");
call user_initialize_thickness(h, g, gv, pf, &
321 just_read_params=just_read)
322 case default ;
call mom_error(fatal,
"MOM_initialize_state: "//&
323 "Unrecognized layer thickness configuration "//trim(config))
327 if ( use_temperature )
then
328 call get_param(pf, mdl,
"TS_CONFIG", config, &
329 "A string that determines how the initial tempertures "//&
330 "and salinities are specified for a new run: \n"//&
331 " \t file - read velocities from the file specified \n"//&
332 " \t\t by (TS_FILE). \n"//&
333 " \t fit - find the temperatures that are consistent with \n"//&
334 " \t\t the layer densities and salinity S_REF. \n"//&
335 " \t TS_profile - use temperature and salinity profiles \n"//&
336 " \t\t (read from TS_FILE) to set layer densities. \n"//&
337 " \t benchmark - use the benchmark test case T & S. \n"//&
338 " \t linear - linear in logical layer space. \n"//&
339 " \t DOME2D - 2D DOME initialization. \n"//&
340 " \t ISOMIP - ISOMIP initialization. \n"//&
341 " \t adjustment2d - 2d lock exchange T/S ICs. \n"//&
342 " \t sloshing - sloshing mode T/S ICs. \n"//&
343 " \t seamount - no motion test with seamount ICs. \n"//&
344 " \t dumbbell - sloshing channel ICs. \n"//&
345 " \t rossby_front - a mixed layer front in thermal wind balance.\n"//&
346 " \t SCM_CVMix_tests - used in the SCM CVMix tests.\n"//&
347 " \t USER - call a user modified routine.", &
348 fail_if_missing=new_sim, do_not_log=just_read)
350 select case (trim(config))
351 case (
"fit");
call initialize_temp_salt_fit(tv%T, tv%S, g, gv, pf, &
352 eos, tv%P_Ref, just_read_params=just_read)
353 case (
"file");
call initialize_temp_salt_from_file(tv%T, tv%S, g, &
354 pf, just_read_params=just_read)
355 case (
"benchmark");
call benchmark_init_temperature_salinity(tv%T, tv%S, &
356 g, gv, pf, eos, tv%P_Ref, just_read_params=just_read)
357 case (
"TS_profile") ;
call initialize_temp_salt_from_profile(tv%T, tv%S, &
358 g, pf, just_read_params=just_read)
359 case (
"linear");
call initialize_temp_salt_linear(tv%T, tv%S, g, pf, &
360 just_read_params=just_read)
361 case (
"DOME2D");
call dome2d_initialize_temperature_salinity ( tv%T, &
362 tv%S, h, g, gv, pf, eos, just_read_params=just_read)
363 case (
"ISOMIP");
call isomip_initialize_temperature_salinity ( tv%T, &
364 tv%S, h, g, gv, pf, eos, just_read_params=just_read)
365 case (
"adjustment2d");
call adjustment_initialize_temperature_salinity ( tv%T, &
366 tv%S, h, g, gv, pf, eos, just_read_params=just_read)
367 case (
"baroclinic_zone");
call baroclinic_zone_init_temperature_salinity( tv%T, &
368 tv%S, h, g, gv, us, pf, just_read_params=just_read)
369 case (
"sloshing");
call sloshing_initialize_temperature_salinity(tv%T, &
370 tv%S, h, g, gv, pf, eos, just_read_params=just_read)
371 case (
"seamount");
call seamount_initialize_temperature_salinity(tv%T, &
372 tv%S, h, g, gv, pf, eos, just_read_params=just_read)
373 case (
"dumbbell");
call dumbbell_initialize_temperature_salinity(tv%T, &
374 tv%S, h, g, gv, pf, eos, just_read_params=just_read)
375 case (
"rossby_front");
call rossby_front_initialize_temperature_salinity ( tv%T, &
376 tv%S, h, g, gv, pf, eos, just_read_params=just_read)
377 case (
"SCM_CVMix_tests");
call scm_cvmix_tests_ts_init(tv%T, tv%S, h, &
378 g, gv, us, pf, just_read_params=just_read)
379 case (
"dense");
call dense_water_initialize_ts(g, gv, pf, eos, tv%T, tv%S, &
380 h, just_read_params=just_read)
381 case (
"USER");
call user_init_temperature_salinity(tv%T, tv%S, g, pf, eos, &
382 just_read_params=just_read)
383 case default ;
call mom_error(fatal,
"MOM_initialize_state: "//&
384 "Unrecognized Temp & salt configuration "//trim(config))
388 if (use_temperature .and. use_obc) &
389 call fill_temp_salt_segments(g, obc, tv)
392 if (new_sim)
call pass_var(h, g%Domain)
395 call get_param(pf, mdl,
"VELOCITY_CONFIG", config, &
396 "A string that determines how the initial velocities "//&
397 "are specified for a new run: \n"//&
398 " \t file - read velocities from the file specified \n"//&
399 " \t\t by (VELOCITY_FILE). \n"//&
400 " \t zero - the fluid is initially at rest. \n"//&
401 " \t uniform - the flow is uniform (determined by\n"//&
402 " \t\t parameters INITIAL_U_CONST and INITIAL_V_CONST).\n"//&
403 " \t rossby_front - a mixed layer front in thermal wind balance.\n"//&
404 " \t soliton - Equatorial Rossby soliton.\n"//&
405 " \t USER - call a user modified routine.", default=
"zero", &
406 do_not_log=just_read)
407 select case (trim(config))
408 case (
"file");
call initialize_velocity_from_file(u, v, g, pf, &
409 just_read_params=just_read)
410 case (
"zero");
call initialize_velocity_zero(u, v, g, pf, &
411 just_read_params=just_read)
412 case (
"uniform");
call initialize_velocity_uniform(u, v, g, pf, &
413 just_read_params=just_read)
414 case (
"circular");
call initialize_velocity_circular(u, v, g, pf, &
415 just_read_params=just_read)
416 case (
"phillips");
call phillips_initialize_velocity(u, v, g, gv, us, pf, &
417 just_read_params=just_read)
418 case (
"rossby_front");
call rossby_front_initialize_velocity(u, v, h, &
419 g, gv, us, pf, just_read_params=just_read)
420 case (
"soliton");
call soliton_initialize_velocity(u, v, h, g)
421 case (
"USER");
call user_initialize_velocity(u, v, g, pf, &
422 just_read_params=just_read)
423 case default ;
call mom_error(fatal,
"MOM_initialize_state: "//&
424 "Unrecognized velocity configuration "//trim(config))
427 if (new_sim)
call pass_vector(u, v, g%Domain)
428 if (debug .and. new_sim)
then
429 call uvchksum(
"MOM_initialize_state [uv]", u, v, g%HI, haloshift=1)
434 call get_param(pf, mdl,
"CONVERT_THICKNESS_UNITS", convert, &
435 "If true, convert the thickness initial conditions from "//&
436 "units of m to kg m-2 or vice versa, depending on whether "//&
437 "BOUSSINESQ is defined. This does not apply if a restart "//&
438 "file is read.", default=.not.gv%Boussinesq, do_not_log=just_read)
440 if (new_sim .and. convert .and. .not.gv%Boussinesq) &
442 call convert_thickness(h, g, gv, us, tv)
445 call get_param(pf, mdl,
"DEPRESS_INITIAL_SURFACE", depress_sfc, &
446 "If true, depress the initial surface to avoid huge "//&
447 "tsunamis when a large surface pressure is applied.", &
448 default=.false., do_not_log=just_read)
449 call get_param(pf, mdl,
"TRIM_IC_FOR_P_SURF", trim_ic_for_p_surf, &
450 "If true, cuts way the top of the column for initial conditions "//&
451 "at the depth where the hydrostatic pressure matches the imposed "//&
452 "surface pressure which is read from file.", default=.false., &
453 do_not_log=just_read)
454 if (depress_sfc .and. trim_ic_for_p_surf)
call mom_error(fatal,
"MOM_initialize_state: "//&
455 "DEPRESS_INITIAL_SURFACE and TRIM_IC_FOR_P_SURF are exclusive and cannot both be True")
456 if (new_sim .and. debug .and. (depress_sfc .or. trim_ic_for_p_surf)) &
457 call hchksum(h,
"Pre-depress: h ", g%HI, haloshift=1, scale=gv%H_to_m)
458 if (depress_sfc)
call depress_surface(h, g, gv, us, pf, tv, just_read_params=just_read)
459 if (trim_ic_for_p_surf)
call trim_for_ice(pf, g, gv, us, ale_csp, tv, h, just_read_params=just_read)
464 call get_param(pf, mdl,
"REGRID_ACCELERATE_INIT", regrid_accelerate, &
465 "If true, runs REGRID_ACCELERATE_ITERATIONS iterations of the regridding "//&
466 "algorithm to push the initial grid to be consistent with the initial "//&
467 "condition. Useful only for state-based and iterative coordinates.", &
468 default=.false., do_not_log=just_read)
469 if (regrid_accelerate)
then
470 call get_param(pf, mdl,
"REGRID_ACCELERATE_ITERATIONS", regrid_iterations, &
471 "The number of regridding iterations to perform to generate "//&
472 "an initial grid that is consistent with the initial conditions.", &
473 default=1, do_not_log=just_read)
475 call get_param(pf, mdl,
"DT", dt,
"Timestep", fail_if_missing=.true.)
477 if (new_sim .and. debug) &
478 call hchksum(h,
"Pre-ALE_regrid: h ", g%HI, haloshift=1, scale=gv%H_to_m)
479 call ale_regrid_accelerated(ale_csp, g, gv, h, tv, regrid_iterations, u, v, tracer_reg, &
480 dt=dt, initial=.true.)
486 if (.not.new_sim)
then
489 call restore_state(dirs%input_filename, dirs%restart_input_dir, time, &
491 if (
present(time_in)) time = time_in
492 if ((gv%m_to_H_restart /= 0.0) .and. (gv%m_to_H_restart /= gv%m_to_H))
then
493 h_rescale = gv%m_to_H / gv%m_to_H_restart
494 do k=1,nz ;
do j=js,je ;
do i=is,ie ; h(i,j,k) = h_rescale * h(i,j,k) ;
enddo ;
enddo ;
enddo
498 if ( use_temperature )
then
499 call pass_var(tv%T, g%Domain, complete=.false.)
500 call pass_var(tv%S, g%Domain, complete=.false.)
502 call pass_var(h, g%Domain)
505 call hchksum(h,
"MOM_initialize_state: h ", g%HI, haloshift=1, scale=gv%H_to_m)
506 if ( use_temperature )
call hchksum(tv%T,
"MOM_initialize_state: T ", g%HI, haloshift=1)
507 if ( use_temperature )
call hchksum(tv%S,
"MOM_initialize_state: S ", g%HI, haloshift=1)
508 if ( use_temperature .and. debug_layers)
then ;
do k=1,nz
509 write(mesg,
'("MOM_IS: T[",I2,"]")') k
510 call hchksum(tv%T(:,:,k), mesg, g%HI, haloshift=1)
511 write(mesg,
'("MOM_IS: S[",I2,"]")') k
512 call hchksum(tv%S(:,:,k), mesg, g%HI, haloshift=1)
516 call get_param(pf, mdl,
"SPONGE", use_sponge, &
517 "If true, sponges may be applied anywhere in the domain. "//&
518 "The exact location and properties of those sponges are "//&
519 "specified via SPONGE_CONFIG.", default=.false.)
520 if ( use_sponge )
then
521 call get_param(pf, mdl,
"SPONGE_CONFIG", config, &
522 "A string that sets how the sponges are configured: \n"//&
523 " \t file - read sponge properties from the file \n"//&
524 " \t\t specified by (SPONGE_FILE).\n"//&
525 " \t ISOMIP - apply ale sponge in the ISOMIP case \n"//&
526 " \t RGC - apply sponge in the rotating_gravity_current case \n"//&
527 " \t DOME - use a slope and channel configuration for the \n"//&
528 " \t\t DOME sill-overflow test case. \n"//&
529 " \t BFB - Sponge at the southern boundary of the domain\n"//&
530 " \t\t for buoyancy-forced basin case.\n"//&
531 " \t USER - call a user modified routine.", default=
"file")
532 select case (trim(config))
533 case (
"DOME");
call dome_initialize_sponges(g, gv, us, tv, pf, sponge_csp)
534 case (
"DOME2D");
call dome2d_initialize_sponges(g, gv, tv, pf, useale, &
535 sponge_csp, ale_sponge_csp)
536 case (
"ISOMIP");
call isomip_initialize_sponges(g, gv, us, tv, pf, useale, &
537 sponge_csp, ale_sponge_csp)
538 case(
"RGC");
call rgc_initialize_sponges(g, gv, tv, u, v, pf, useale, &
539 sponge_csp, ale_sponge_csp)
540 case (
"USER");
call user_initialize_sponges(g, gv, use_temperature, tv, pf, sponge_csp, h)
541 case (
"BFB");
call bfb_initialize_sponges_southonly(g, gv, us, use_temperature, tv, pf, &
543 case (
"DUMBBELL");
call dumbbell_initialize_sponges(g, gv, us, tv, pf, useale, &
544 sponge_csp, ale_sponge_csp)
545 case (
"phillips");
call phillips_initialize_sponges(g, gv, us, tv, pf, sponge_csp, h)
546 case (
"dense");
call dense_water_initialize_sponges(g, gv, tv, pf, useale, &
547 sponge_csp, ale_sponge_csp)
548 case (
"file");
call initialize_sponges_file(g, gv, us, use_temperature, tv, pf, &
549 sponge_csp, ale_sponge_csp, time)
550 case default ;
call mom_error(fatal,
"MOM_initialize_state: "//&
551 "Unrecognized sponge configuration "//trim(config))
556 call open_boundary_init(g, pf, obc)
559 if (
associated(obc))
then
560 call get_param(pf, mdl,
"OBC_USER_CONFIG", config, &
561 "A string that sets how the user code is invoked to set open boundary data: \n"//&
562 " DOME - specified inflow on northern boundary\n"//&
563 " dyed_channel - supercritical with dye on the inflow boundary\n"//&
564 " dyed_obcs - circle_obcs with dyes on the open boundaries\n"//&
565 " Kelvin - barotropic Kelvin wave forcing on the western boundary\n"//&
566 " shelfwave - Flather with shelf wave forcing on western boundary\n"//&
567 " supercritical - now only needed here for the allocations\n"//&
568 " tidal_bay - Flather with tidal forcing on eastern boundary\n"//&
569 " USER - user specified", default=
"none")
570 if (trim(config) ==
"DOME")
then
571 call dome_set_obc_data(obc, tv, g, gv, us, pf, tracer_reg)
572 elseif (trim(config) ==
"dyed_channel")
then
573 call dyed_channel_set_obc_tracer_data(obc, g, gv, pf, tracer_reg)
574 obc%update_OBC = .true.
575 elseif (trim(config) ==
"dyed_obcs")
then
576 call dyed_obcs_set_obc_data(obc, g, gv, pf, tracer_reg)
577 elseif (trim(config) ==
"Kelvin")
then
578 obc%update_OBC = .true.
579 elseif (trim(config) ==
"shelfwave")
then
580 obc%update_OBC = .true.
581 elseif (lowercase(trim(config)) ==
"supercritical")
then
582 call supercritical_set_obc_data(obc, g, pf)
583 elseif (trim(config) ==
"tidal_bay")
then
584 obc%update_OBC = .true.
585 elseif (trim(config) ==
"USER")
then
586 call user_set_obc_data(obc, tv, g, pf, tracer_reg)
587 elseif (.not. trim(config) ==
"none")
then
588 call mom_error(fatal,
"The open boundary conditions specified by "//&
589 "OBC_USER_CONFIG = "//trim(config)//
" have not been fully implemented.")
591 if (open_boundary_query(obc, apply_open_obc=.true.))
then
592 call set_tracer_data(obc, tv, h, g, pf, tracer_reg)
599 if (debug.and.
associated(obc))
then
600 call hchksum(g%mask2dT,
'MOM_initialize_state: mask2dT ', g%HI)
601 call uvchksum(
'MOM_initialize_state: mask2dC[uv]', g%mask2dCu, &
603 call qchksum(g%mask2dBu,
'MOM_initialize_state: mask2dBu ', g%HI)
606 if (debug_obc)
call open_boundary_test_extern_h(g, obc, h)
607 call calltree_leave(
'MOM_initialize_state()')