8 use mom_cpu_clock,
only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
11 use mom_domains,
only : root_pe, to_all, scalar_pair, cgrid_ne, agrid
12 use mom_eos,
only : find_depth_of_pressure_in_cell
32 use mom_sponge,
only : set_up_sponge_field, set_up_sponge_ml_density
42 use mom_ale,
only : pressure_gradient_plm
44 use mom_eos,
only : int_specific_vol_dp, convert_temp_salt_for_teos10
55 use rgc_initialization,
only : rgc_initialize_sponges
94 use mom_ale,
only : ale_initregridding,
ale_cs, ale_initthicknesstocoord
95 use mom_ale,
only : ale_remap_scalar, ale_build_grid, ale_regrid_accelerated
102 use fms_io_mod,
only : field_size
104 implicit none ;
private
106 #include <MOM_memory.h>
108 public mom_initialize_state
115 character(len=40) :: mdl =
"MOM_state_initialization"
121 subroutine mom_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
122 restart_CS, ALE_CSp, tracer_Reg, sponge_CSp, &
123 ALE_sponge_CSp, OBC, Time_in)
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)), &
137 type(time_type),
intent(inout) :: time
144 type(
ale_cs),
pointer :: ale_csp
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")
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))
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.)
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()')
609 end subroutine mom_initialize_state
612 subroutine initialize_thickness_from_file(h, G, GV, US, param_file, file_has_thickness, &
617 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
621 logical,
intent(in) :: file_has_thickness
624 logical,
optional,
intent(in) :: just_read_params
628 real :: eta(SZI_(G),SZJ_(G),SZK_(G)+1)
629 integer :: inconsistent = 0
630 logical :: correct_thickness
632 character(len=40) :: mdl =
"initialize_thickness_from_file"
633 character(len=200) :: filename, thickness_file, inputdir, mesg
634 integer :: i, j, k, is, ie, js, je, nz
636 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
638 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
640 if (.not.just_read) &
641 call calltree_enter(trim(mdl)//
"(), MOM_state_initialization.F90")
643 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".", do_not_log=just_read)
644 inputdir = slasher(inputdir)
645 call get_param(param_file, mdl,
"THICKNESS_FILE", thickness_file, &
646 "The name of the thickness file.", &
647 fail_if_missing=.not.just_read, do_not_log=just_read)
649 filename = trim(inputdir)//trim(thickness_file)
650 if (.not.just_read)
call log_param(param_file, mdl,
"INPUTDIR/THICKNESS_FILE", filename)
652 if ((.not.just_read) .and. (.not.
file_exists(filename, g%Domain)))
call mom_error(fatal, &
653 " initialize_thickness_from_file: Unable to open "//trim(filename))
655 if (file_has_thickness)
then
657 if (just_read)
return
658 call mom_read_data(filename,
"h", h(:,:,:), g%Domain, scale=gv%m_to_H)
660 call get_param(param_file, mdl,
"ADJUST_THICKNESS", correct_thickness, &
661 "If true, all mass below the bottom removed if the "//&
662 "topography is shallower than the thickness input file "//&
663 "would indicate.", default=.false., do_not_log=just_read)
664 if (just_read)
return
666 call mom_read_data(filename,
"eta", eta(:,:,:), g%Domain, scale=us%m_to_Z)
668 if (correct_thickness)
then
669 call adjustetatofitbathymetry(g, gv, us, eta, h)
671 do k=nz,1,-1 ;
do j=js,je ;
do i=is,ie
672 if (eta(i,j,k) < (eta(i,j,k+1) + gv%Angstrom_Z))
then
673 eta(i,j,k) = eta(i,j,k+1) + gv%Angstrom_Z
674 h(i,j,k) = gv%Angstrom_H
676 h(i,j,k) = gv%Z_to_H * (eta(i,j,k) - eta(i,j,k+1))
678 enddo ;
enddo ;
enddo
680 do j=js,je ;
do i=is,ie
681 if (abs(eta(i,j,nz+1) + g%bathyT(i,j)) > 1.0*us%m_to_Z) &
682 inconsistent = inconsistent + 1
684 call sum_across_pes(inconsistent)
686 if ((inconsistent > 0) .and. (is_root_pe()))
then
687 write(mesg,
'("Thickness initial conditions are inconsistent ",'// &
688 '"with topography in ",I8," places.")') inconsistent
689 call mom_error(warning, mesg)
694 call calltree_leave(trim(mdl)//
'()')
695 end subroutine initialize_thickness_from_file
705 subroutine adjustetatofitbathymetry(G, GV, US, eta, h)
709 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1),
intent(inout) :: eta
710 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(inout) :: h
712 integer :: i, j, k, is, ie, js, je, nz, contractions, dilations
713 real :: hTolerance = 0.1
714 real :: hTmp, eTmp, dilate
715 character(len=100) :: mesg
717 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
718 htolerance = 0.1*us%m_to_Z
721 do j=js,je ;
do i=is,ie
722 if (-eta(i,j,nz+1) > g%bathyT(i,j) + htolerance)
then
723 eta(i,j,nz+1) = -g%bathyT(i,j)
724 contractions = contractions + 1
727 call sum_across_pes(contractions)
728 if ((contractions > 0) .and. (is_root_pe()))
then
729 write(mesg,
'("Thickness initial conditions were contracted ",'// &
730 '"to fit topography in ",I8," places.")') contractions
731 call mom_error(warning,
'adjustEtaToFitBathymetry: '//mesg)
736 do k=nz,1,-1 ;
do j=js,je ;
do i=is,ie
739 if (eta(i,j,k) < (eta(i,j,k+1) + gv%Angstrom_Z))
then
740 eta(i,j,k) = eta(i,j,k+1) + gv%Angstrom_Z
741 h(i,j,k) = gv%Angstrom_Z
743 h(i,j,k) = (eta(i,j,k) - eta(i,j,k+1))
745 enddo ;
enddo ;
enddo
748 do j=js,je ;
do i=is,ie
752 if (-eta(i,j,nz+1) < g%bathyT(i,j) - htolerance)
then
753 dilations = dilations + 1
754 if (eta(i,j,1) <= eta(i,j,nz+1))
then
755 do k=1,nz ; h(i,j,k) = (eta(i,j,1) + g%bathyT(i,j)) / real(nz) ;
enddo
757 dilate = (eta(i,j,1) + g%bathyT(i,j)) / (eta(i,j,1) - eta(i,j,nz+1))
758 do k=1,nz ; h(i,j,k) = h(i,j,k) * dilate ;
enddo
760 do k=nz,2,-1 ; eta(i,j,k) = eta(i,j,k+1) + h(i,j,k) ;
enddo
765 do k=1,nz ;
do j=js,je ;
do i=is,ie
766 h(i,j,k) = h(i,j,k)*gv%Z_to_H
767 enddo ;
enddo ;
enddo
769 call sum_across_pes(dilations)
770 if ((dilations > 0) .and. (is_root_pe()))
then
771 write(mesg,
'("Thickness initial conditions were dilated ",'// &
772 '"to fit topography in ",I8," places.")') dilations
773 call mom_error(warning,
'adjustEtaToFitBathymetry: '//mesg)
776 end subroutine adjustetatofitbathymetry
779 subroutine initialize_thickness_uniform(h, G, GV, param_file, just_read_params)
782 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
786 logical,
optional,
intent(in) :: just_read_params
789 character(len=40) :: mdl =
"initialize_thickness_uniform"
790 real :: e0(SZK_(G)+1)
792 real :: eta1D(SZK_(G)+1)
795 integer :: i, j, k, is, ie, js, je, nz
797 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
799 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
801 if (just_read)
return
803 call calltree_enter(trim(mdl)//
"(), MOM_state_initialization.F90")
805 if (g%max_depth<=0.)
call mom_error(fatal,
"initialize_thickness_uniform: "// &
806 "MAXIMUM_DEPTH has a non-sensical value! Was it set?")
809 e0(k) = -g%max_depth * real(k-1) / real(nz)
812 do j=js,je ;
do i=is,ie
818 eta1d(nz+1) = -g%bathyT(i,j)
821 if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z))
then
822 eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
823 h(i,j,k) = gv%Angstrom_H
825 h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
830 call calltree_leave(trim(mdl)//
'()')
831 end subroutine initialize_thickness_uniform
834 subroutine initialize_thickness_list(h, G, GV, US, param_file, just_read_params)
838 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
842 logical,
optional,
intent(in) :: just_read_params
845 character(len=40) :: mdl =
"initialize_thickness_list"
846 real :: e0(SZK_(G)+1)
848 real :: eta1D(SZK_(G)+1)
851 character(len=200) :: filename, eta_file, inputdir
852 character(len=72) :: eta_var
853 integer :: i, j, k, is, ie, js, je, nz
855 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
857 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
859 call get_param(param_file, mdl,
"INTERFACE_IC_FILE", eta_file, &
860 "The file from which horizontal mean initial conditions "//&
861 "for interface depths can be read.", fail_if_missing=.true.)
862 call get_param(param_file, mdl,
"INTERFACE_IC_VAR", eta_var, &
863 "The variable name for horizontal mean initial conditions "//&
864 "for interface depths relative to mean sea level.", &
867 if (just_read)
return
869 call calltree_enter(trim(mdl)//
"(), MOM_state_initialization.F90")
871 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
872 filename = trim(slasher(inputdir))//trim(eta_file)
873 call log_param(param_file, mdl,
"INPUTDIR/INTERFACE_IC_FILE", filename)
876 call mom_read_data(filename, eta_var, e0(:), scale=us%m_to_Z)
878 if ((abs(e0(1)) - 0.0) > 0.001)
then
880 do k=nz+1,2,-1 ; e0(k) = e0(k-1) ;
enddo
884 if (e0(2) > e0(1))
then
885 do k=1,nz ; e0(k) = -e0(k) ;
enddo
888 do j=js,je ;
do i=is,ie
894 eta1d(nz+1) = -g%bathyT(i,j)
897 if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z))
then
898 eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
899 h(i,j,k) = gv%Angstrom_H
901 h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
906 call calltree_leave(trim(mdl)//
'()')
907 end subroutine initialize_thickness_list
910 subroutine initialize_thickness_search
911 call mom_error(fatal,
" MOM_state_initialization.F90, initialize_thickness_search: NOT IMPLEMENTED")
912 end subroutine initialize_thickness_search
915 subroutine convert_thickness(h, G, GV, US, tv)
919 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
925 real,
dimension(SZI_(G),SZJ_(G)) :: &
927 real :: dz_geo(SZI_(G),SZJ_(G))
933 logical :: Boussinesq
934 integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
935 integer :: itt, max_itt
937 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
938 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
940 boussinesq = gv%Boussinesq
941 i_gearth = 1.0 / (gv%mks_g_Earth)
942 hm_rho_to_pa = gv%mks_g_Earth * gv%H_to_m
945 call mom_error(fatal,
"Not yet converting thickness with Boussinesq approx.")
947 if (
associated(tv%eqn_of_state))
then
948 do j=jsq,jeq+1 ;
do i=isq,ieq+1
949 p_bot(i,j) = 0.0 ; p_top(i,j) = 0.0
953 do i=is,ie ; p_top(i,j) = p_bot(i,j) ;
enddo
955 is, ie-is+1, tv%eqn_of_state)
957 p_bot(i,j) = p_top(i,j) + hm_rho_to_pa * (h(i,j,k) * rho(i))
962 call int_specific_vol_dp(tv%T(:,:,k), tv%S(:,:,k), p_top, p_bot, &
963 0.0, g%HI, tv%eqn_of_state, dz_geo)
964 if (itt < max_itt)
then ;
do j=js,je
966 is, ie-is+1, tv%eqn_of_state)
971 p_bot(i,j) = p_bot(i,j) + rho(i) * &
972 (hm_rho_to_pa*h(i,j,k) - dz_geo(i,j))
977 do j=js,je ;
do i=is,ie
978 h(i,j,k) = (p_bot(i,j) - p_top(i,j)) * gv%kg_m2_to_H * i_gearth
982 do k=1,nz ;
do j=js,je ;
do i=is,ie
983 h(i,j,k) = (h(i,j,k) * gv%Rlay(k)) * hm_rho_to_pa * gv%kg_m2_to_H**2
986 enddo ;
enddo ;
enddo
990 end subroutine convert_thickness
993 subroutine depress_surface(h, G, GV, US, param_file, tv, just_read_params)
997 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1001 logical,
optional,
intent(in) :: just_read_params
1004 real,
dimension(SZI_(G),SZJ_(G)) :: &
1006 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &
1009 real :: scale_factor
1011 character(len=40) :: mdl =
"depress_surface"
1012 character(len=200) :: inputdir, eta_srf_file
1013 character(len=200) :: filename, eta_srf_var
1014 logical :: just_read
1015 integer :: i, j, k, is, ie, js, je, nz
1016 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1018 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
1022 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
1023 inputdir = slasher(inputdir)
1024 call get_param(param_file, mdl,
"SURFACE_HEIGHT_IC_FILE", eta_srf_file,&
1025 "The initial condition file for the surface height.", &
1026 fail_if_missing=.not.just_read, do_not_log=just_read)
1027 call get_param(param_file, mdl,
"SURFACE_HEIGHT_IC_VAR", eta_srf_var, &
1028 "The initial condition variable for the surface height.",&
1029 default=
"SSH", do_not_log=just_read)
1030 filename = trim(inputdir)//trim(eta_srf_file)
1031 if (.not.just_read) &
1032 call log_param(param_file, mdl,
"INPUTDIR/SURFACE_HEIGHT_IC_FILE", filename)
1034 call get_param(param_file, mdl,
"SURFACE_HEIGHT_IC_SCALE", scale_factor, &
1035 "A scaling factor to convert SURFACE_HEIGHT_IC_VAR into "//&
1036 "units of m", units=
"variable", default=1.0, do_not_log=just_read)
1038 if (just_read)
return
1040 call mom_read_data(filename, eta_srf_var, eta_sfc, g%Domain, scale=scale_factor)
1043 call find_eta(h, tv, g, gv, us, eta, eta_to_m=1.0)
1045 do j=js,je ;
do i=is,ie ;
if (g%mask2dT(i,j) > 0.0)
then
1049 if (eta_sfc(i,j) > eta(i,j,1))
then
1051 if (eta_sfc(i,j) - eta(i,j,nz+1) > 10.0*(eta(i,j,1) - eta(i,j,nz+1)))
then
1053 call mom_error(warning,
"Free surface height dilation attempted "//&
1054 "to exceed 10-fold.", all_print=.true.)
1056 dilate = (eta_sfc(i,j) - eta(i,j,nz+1)) / (eta(i,j,1) - eta(i,j,nz+1))
1058 do k=1,nz ; h(i,j,k) = h(i,j,k) * dilate ;
enddo
1059 elseif (eta(i,j,1) > eta_sfc(i,j))
then
1062 if (eta(i,j,k) <= eta_sfc(i,j))
exit
1063 if (eta(i,j,k+1) >= eta_sfc(i,j))
then
1064 h(i,j,k) = gv%Angstrom_H
1066 h(i,j,k) = max(gv%Angstrom_H, h(i,j,k) * &
1067 (eta_sfc(i,j) - eta(i,j,k+1)) / (eta(i,j,k) - eta(i,j,k+1)) )
1071 endif ;
enddo ;
enddo
1073 end subroutine depress_surface
1077 subroutine trim_for_ice(PF, G, GV, US, ALE_CSp, tv, h, just_read_params)
1082 type(
ale_cs),
pointer :: ALE_CSp
1084 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1086 logical,
optional,
intent(in) :: just_read_params
1089 character(len=200) :: mdl =
"trim_for_ice"
1090 real,
dimension(SZI_(G),SZJ_(G)) :: p_surf
1091 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: S_t, S_b
1092 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: T_t, T_b
1093 character(len=200) :: inputdir, filename, p_surf_file, p_surf_var
1094 real :: scale_factor
1095 real :: min_thickness
1097 logical :: just_read
1098 logical :: use_remapping
1101 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
1103 call get_param(pf, mdl,
"SURFACE_PRESSURE_FILE", p_surf_file, &
1104 "The initial condition file for the surface height.", &
1105 fail_if_missing=.not.just_read, do_not_log=just_read)
1106 call get_param(pf, mdl,
"SURFACE_PRESSURE_VAR", p_surf_var, &
1107 "The initial condition variable for the surface height.", &
1108 units=
"kg m-2", default=
"", do_not_log=just_read)
1109 call get_param(pf, mdl,
"INPUTDIR", inputdir, default=
".", do_not_log=.true.)
1110 filename = trim(slasher(inputdir))//trim(p_surf_file)
1111 if (.not.just_read)
call log_param(pf, mdl,
"!INPUTDIR/SURFACE_HEIGHT_IC_FILE", filename)
1113 call get_param(pf, mdl,
"SURFACE_PRESSURE_SCALE", scale_factor, &
1114 "A scaling factor to convert SURFACE_PRESSURE_VAR from "//&
1115 "file SURFACE_PRESSURE_FILE into a surface pressure.", &
1116 units=
"file dependent", default=1., do_not_log=just_read)
1117 call get_param(pf, mdl,
"MIN_THICKNESS", min_thickness,
'Minimum layer thickness', &
1118 units=
'm', default=1.e-3, do_not_log=just_read, scale=us%m_to_Z)
1119 call get_param(pf, mdl,
"TRIMMING_USES_REMAPPING", use_remapping, &
1120 'When trimming the column, also remap T and S.', &
1121 default=.false., do_not_log=just_read)
1123 if (just_read)
return
1125 call mom_read_data(filename, p_surf_var, p_surf, g%Domain, scale=scale_factor)
1127 if (use_remapping)
then
1129 call initialize_remapping(remap_cs,
'PLM', boundary_extrapolation=.true.)
1133 if (
associated(ale_csp) )
then
1134 call pressure_gradient_plm(ale_csp, s_t, s_b, t_t, t_b, g, gv, tv, h, .true.)
1137 do k=1,g%ke ;
do j=g%jsc,g%jec ;
do i=g%isc,g%iec
1138 t_t(i,j,k) = tv%T(i,j,k) ; t_b(i,j,k) = tv%T(i,j,k)
1139 s_t(i,j,k) = tv%S(i,j,k) ; s_b(i,j,k) = tv%S(i,j,k)
1140 enddo ;
enddo ;
enddo
1143 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
1144 call cut_off_column_top(gv%ke, tv, gv, gv%mks_g_Earth*us%Z_to_m, g%bathyT(i,j), &
1145 min_thickness, tv%T(i,j,:), t_t(i,j,:), t_b(i,j,:), &
1146 tv%S(i,j,:), s_t(i,j,:), s_b(i,j,:), p_surf(i,j), h(i,j,:), remap_cs, &
1147 z_tol=1.0e-5*us%m_to_Z)
1150 end subroutine trim_for_ice
1155 subroutine cut_off_column_top(nk, tv, GV, G_earth, depth, min_thickness, &
1156 T, T_t, T_b, S, S_t, S_b, p_surf, h, remap_CS, z_tol)
1157 integer,
intent(in) :: nk
1160 real,
intent(in) :: G_earth
1161 real,
intent(in) :: depth
1162 real,
intent(in) :: min_thickness
1163 real,
dimension(nk),
intent(inout) :: T
1164 real,
dimension(nk),
intent(in) :: T_t
1165 real,
dimension(nk),
intent(in) :: T_b
1166 real,
dimension(nk),
intent(inout) :: S
1167 real,
dimension(nk),
intent(in) :: S_t
1168 real,
dimension(nk),
intent(in) :: S_b
1169 real,
intent(in) :: p_surf
1170 real,
dimension(nk),
intent(inout) :: h
1173 real,
optional,
intent(in) :: z_tol
1177 real,
dimension(nk+1) :: e
1178 real,
dimension(nk) :: h0, S0, T0, h1, S1, T1
1179 real :: P_t, P_b, z_out, e_top
1185 e(k) = e(k+1) + gv%H_to_Z*h(k)
1192 call find_depth_of_pressure_in_cell(t_t(k), t_b(k), s_t(k), s_b(k), e(k), e(k+1), &
1193 p_t, p_surf, gv%Rho0, g_earth, tv%eqn_of_state, &
1194 p_b, z_out, z_tol=z_tol)
1195 if (z_out>=e(k))
then
1198 elseif (z_out<=e(k+1))
then
1208 if (e_top<e(1))
then
1211 if (e(k) > e_top)
then
1214 e_top = e_top - min_thickness
1217 h(k) = gv%Z_to_H * max( min_thickness, e(k) - e(k+1) )
1218 if (e(k) < e_top)
exit
1224 if (
associated(remap_cs))
then
1230 call remapping_core_h(remap_cs, nk, h0, t0, nk, h1, t1, 1.0e-30*gv%m_to_H, 1.0e-10*gv%m_to_H)
1231 call remapping_core_h(remap_cs, nk, h0, s0, nk, h1, s1, 1.0e-30*gv%m_to_H, 1.0e-10*gv%m_to_H)
1238 end subroutine cut_off_column_top
1241 subroutine initialize_velocity_from_file(u, v, G, param_file, just_read_params)
1243 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1245 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1249 logical,
optional,
intent(in) :: just_read_params
1252 character(len=40) :: mdl =
"initialize_velocity_from_file"
1253 character(len=200) :: filename,velocity_file,inputdir
1254 logical :: just_read
1256 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
1258 if (.not.just_read)
call calltree_enter(trim(mdl)//
"(), MOM_state_initialization.F90")
1260 call get_param(param_file, mdl,
"VELOCITY_FILE", velocity_file, &
1261 "The name of the velocity initial condition file.", &
1262 fail_if_missing=.not.just_read, do_not_log=just_read)
1263 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
1264 inputdir = slasher(inputdir)
1266 if (just_read)
return
1268 filename = trim(inputdir)//trim(velocity_file)
1269 call log_param(param_file, mdl,
"INPUTDIR/VELOCITY_FILE", filename)
1271 if (.not.
file_exists(filename, g%Domain))
call mom_error(fatal, &
1272 " initialize_velocity_from_file: Unable to open "//trim(filename))
1277 call calltree_leave(trim(mdl)//
'()')
1278 end subroutine initialize_velocity_from_file
1281 subroutine initialize_velocity_zero(u, v, G, param_file, just_read_params)
1283 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1285 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1289 logical,
optional,
intent(in) :: just_read_params
1292 character(len=200) :: mdl =
"initialize_velocity_zero"
1293 logical :: just_read
1294 integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
1295 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1296 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1298 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
1300 if (.not.just_read)
call calltree_enter(trim(mdl)//
"(), MOM_state_initialization.F90")
1302 if (just_read)
return
1304 do k=1,nz ;
do j=js,je ;
do i=isq,ieq
1306 enddo ;
enddo ;
enddo
1307 do k=1,nz ;
do j=jsq,jeq ;
do i=is,ie
1309 enddo ;
enddo ;
enddo
1311 call calltree_leave(trim(mdl)//
'()')
1312 end subroutine initialize_velocity_zero
1315 subroutine initialize_velocity_uniform(u, v, G, param_file, just_read_params)
1317 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1319 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1323 logical,
optional,
intent(in) :: just_read_params
1326 integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
1327 real :: initial_u_const, initial_v_const
1328 logical :: just_read
1329 character(len=200) :: mdl =
"initialize_velocity_uniform"
1330 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1331 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1333 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
1335 call get_param(param_file, mdl,
"INITIAL_U_CONST", initial_u_const, &
1336 "A initial uniform value for the zonal flow.", &
1337 units=
"m s-1", fail_if_missing=.not.just_read, do_not_log=just_read)
1338 call get_param(param_file, mdl,
"INITIAL_V_CONST", initial_v_const, &
1339 "A initial uniform value for the meridional flow.", &
1340 units=
"m s-1", fail_if_missing=.not.just_read, do_not_log=just_read)
1342 if (just_read)
return
1344 do k=1,nz ;
do j=js,je ;
do i=isq,ieq
1345 u(i,j,k) = initial_u_const
1346 enddo ;
enddo ;
enddo
1347 do k=1,nz ;
do j=jsq,jeq ;
do i=is,ie
1348 v(i,j,k) = initial_v_const
1349 enddo ;
enddo ;
enddo
1351 end subroutine initialize_velocity_uniform
1355 subroutine initialize_velocity_circular(u, v, G, param_file, just_read_params)
1357 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1359 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1363 logical,
optional,
intent(in) :: just_read_params
1366 character(len=200) :: mdl =
"initialize_velocity_circular"
1367 real :: circular_max_u
1368 real :: dpi, psi1, psi2
1369 logical :: just_read
1370 integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
1371 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1372 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1374 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
1376 call get_param(param_file, mdl,
"CIRCULAR_MAX_U", circular_max_u, &
1377 "The amplitude of zonal flow from which to scale the "// &
1378 "circular stream function [m s-1].", &
1379 units=
"m s-1", default=0., do_not_log=just_read)
1381 if (just_read)
return
1385 do k=1,nz ;
do j=js,je ;
do i=isq,ieq
1387 psi2 = my_psi(i,j-1)
1388 u(i,j,k) = (psi1-psi2)/g%dy_Cu(i,j)
1389 enddo ;
enddo ;
enddo
1390 do k=1,nz ;
do j=jsq,jeq ;
do i=is,ie
1392 psi2 = my_psi(i-1,j)
1393 v(i,j,k) = (psi2-psi1)/g%dx_Cv(i,j)
1394 enddo ;
enddo ;
enddo
1399 real function my_psi(ig,jg)
1405 x = 2.0*(g%geoLonBu(ig,jg)-g%west_lon)/g%len_lon-1.0
1406 y = 2.0*(g%geoLatBu(ig,jg)-g%south_lat)/g%len_lat-1.0
1407 r = sqrt( x**2 + y**2 )
1409 my_psi = 0.5*(1.0 - cos(dpi*r))
1410 my_psi = my_psi * (circular_max_u*g%len_lon*1e3/dpi)
1413 end subroutine initialize_velocity_circular
1416 subroutine initialize_temp_salt_from_file(T, S, G, param_file, just_read_params)
1417 type(ocean_grid_type),
intent(in) :: G
1418 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(out) :: T
1420 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(out) :: S
1423 logical,
optional,
intent(in) :: just_read_params
1426 logical :: just_read
1427 character(len=200) :: filename, salt_filename
1428 character(len=200) :: ts_file, salt_file, inputdir
1429 character(len=40) :: mdl =
"initialize_temp_salt_from_file"
1430 character(len=64) :: temp_var, salt_var
1432 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
1434 if (.not.just_read)
call calltree_enter(trim(mdl)//
"(), MOM_state_initialization.F90")
1436 call get_param(param_file, mdl,
"TS_FILE", ts_file, &
1437 "The initial condition file for temperature.", &
1438 fail_if_missing=.not.just_read, do_not_log=just_read)
1439 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
1440 inputdir = slasher(inputdir)
1442 filename = trim(inputdir)//trim(ts_file)
1443 if (.not.just_read)
call log_param(param_file, mdl,
"INPUTDIR/TS_FILE", filename)
1444 call get_param(param_file, mdl,
"TEMP_IC_VAR", temp_var, &
1445 "The initial condition variable for potential temperature.", &
1446 default=
"PTEMP", do_not_log=just_read)
1447 call get_param(param_file, mdl,
"SALT_IC_VAR", salt_var, &
1448 "The initial condition variable for salinity.", &
1449 default=
"SALT", do_not_log=just_read)
1450 call get_param(param_file, mdl,
"SALT_FILE", salt_file, &
1451 "The initial condition file for salinity.", &
1452 default=trim(ts_file), do_not_log=just_read)
1454 if (just_read)
return
1456 if (.not.file_exists(filename, g%Domain))
call mom_error(fatal, &
1457 " initialize_temp_salt_from_file: Unable to open "//trim(filename))
1460 call mom_read_data(filename, temp_var, t(:,:,:), g%Domain)
1462 salt_filename = trim(inputdir)//trim(salt_file)
1463 if (.not.file_exists(salt_filename, g%Domain))
call mom_error(fatal, &
1464 " initialize_temp_salt_from_file: Unable to open "//trim(salt_filename))
1466 call mom_read_data(salt_filename, salt_var, s(:,:,:), g%Domain)
1468 call calltree_leave(trim(mdl)//
'()')
1469 end subroutine initialize_temp_salt_from_file
1472 subroutine initialize_temp_salt_from_profile(T, S, G, param_file, just_read_params)
1473 type(ocean_grid_type),
intent(in) :: G
1474 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(out) :: T
1476 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(out) :: S
1479 logical,
optional,
intent(in) :: just_read_params
1482 real,
dimension(SZK_(G)) :: T0, S0
1484 logical :: just_read
1485 character(len=200) :: filename, ts_file, inputdir
1486 character(len=40) :: mdl =
"initialize_temp_salt_from_profile"
1488 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
1490 if (.not.just_read)
call calltree_enter(trim(mdl)//
"(), MOM_state_initialization.F90")
1492 call get_param(param_file, mdl,
"TS_FILE", ts_file, &
1493 "The file with the reference profiles for temperature "//&
1494 "and salinity.", fail_if_missing=.not.just_read, do_not_log=just_read)
1496 if (just_read)
return
1498 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
1499 inputdir = slasher(inputdir)
1500 filename = trim(inputdir)//trim(ts_file)
1501 call log_param(param_file, mdl,
"INPUTDIR/TS_FILE", filename)
1502 if (.not.file_exists(filename))
call mom_error(fatal, &
1503 " initialize_temp_salt_from_profile: Unable to open "//trim(filename))
1506 call mom_read_data(filename,
"PTEMP", t0(:))
1507 call mom_read_data(filename,
"SALT", s0(:))
1509 do k=1,g%ke ;
do j=g%jsc,g%jec ;
do i=g%isc,g%iec
1510 t(i,j,k) = t0(k) ; s(i,j,k) = s0(k)
1511 enddo ;
enddo ;
enddo
1513 call calltree_leave(trim(mdl)//
'()')
1514 end subroutine initialize_temp_salt_from_profile
1517 subroutine initialize_temp_salt_fit(T, S, G, GV, param_file, eqn_of_state, P_Ref, just_read_params)
1518 type(ocean_grid_type),
intent(in) :: G
1519 type(verticalgrid_type),
intent(in) :: GV
1520 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(out) :: T
1522 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(out) :: S
1526 type(eos_type),
pointer :: eqn_of_state
1527 real,
intent(in) :: P_Ref
1529 logical,
optional,
intent(in) :: just_read_params
1536 real :: pres(SZK_(G))
1537 real :: drho_dT(SZK_(G))
1538 real :: drho_dS(SZK_(G))
1539 real :: rho_guess(SZK_(G))
1540 logical :: fit_salin
1541 logical :: just_read
1542 character(len=40) :: mdl =
"initialize_temp_salt_fit"
1543 integer :: i, j, k, itt, nz
1546 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
1548 if (.not.just_read)
call calltree_enter(trim(mdl)//
"(), MOM_state_initialization.F90")
1550 call get_param(param_file, mdl,
"T_REF", t_ref, &
1551 "A reference temperature used in initialization.", &
1552 units=
"degC", fail_if_missing=.not.just_read, do_not_log=just_read)
1553 call get_param(param_file, mdl,
"S_REF", s_ref, &
1554 "A reference salinity used in initialization.", units=
"PSU", &
1555 default=35.0, do_not_log=just_read)
1556 call get_param(param_file, mdl,
"FIT_SALINITY", fit_salin, &
1557 "If true, accept the prescribed temperature and fit the "//&
1558 "salinity; otherwise take salinity and fit temperature.", &
1559 default=.false., do_not_log=just_read)
1561 if (just_read)
return
1564 pres(k) = p_ref ; s0(k) = s_ref
1568 call calculate_density(t0(1),s0(1),pres(1),rho_guess(1),eqn_of_state)
1569 call calculate_density_derivs(t0,s0,pres,drho_dt,drho_ds,1,1,eqn_of_state)
1574 s0(k) = max(0.0, s0(1) + (gv%Rlay(k) - rho_guess(1)) / drho_ds(1))
1578 call calculate_density(t0,s0,pres,rho_guess,1,nz,eqn_of_state)
1579 call calculate_density_derivs(t0,s0,pres,drho_dt,drho_ds,1,nz,eqn_of_state)
1581 s0(k) = max(0.0, s0(k) + (gv%Rlay(k) - rho_guess(k)) / drho_ds(k))
1587 t0(k) = t0(1) + (gv%Rlay(k) - rho_guess(1)) / drho_dt(1)
1590 call calculate_density(t0,s0,pres,rho_guess,1,nz,eqn_of_state)
1591 call calculate_density_derivs(t0,s0,pres,drho_dt,drho_ds,1,nz,eqn_of_state)
1593 t0(k) = t0(k) + (gv%Rlay(k) - rho_guess(k)) / drho_dt(k)
1598 do k=1,nz ;
do j=g%jsd,g%jed ;
do i=g%isd,g%ied
1599 t(i,j,k) = t0(k) ; s(i,j,k) = s0(k)
1600 enddo ;
enddo ;
enddo
1602 call calltree_leave(trim(mdl)//
'()')
1603 end subroutine initialize_temp_salt_fit
1610 subroutine initialize_temp_salt_linear(T, S, G, param_file, just_read_params)
1611 type(ocean_grid_type),
intent(in) :: G
1612 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(out) :: T
1614 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(out) :: S
1618 logical,
optional,
intent(in) :: just_read_params
1624 real :: delta_S, delta_T
1625 real :: S_top, T_top
1626 real :: S_range, T_range
1628 logical :: just_read
1629 character(len=40) :: mdl =
"initialize_temp_salt_linear"
1631 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
1633 if (.not.just_read)
call calltree_enter(trim(mdl)//
"(), MOM_state_initialization.F90")
1634 call get_param(param_file, mdl,
"T_TOP", t_top, &
1635 "Initial temperature of the top surface.", &
1636 units=
"degC", fail_if_missing=.not.just_read, do_not_log=just_read)
1637 call get_param(param_file, mdl,
"T_RANGE", t_range, &
1638 "Initial temperature difference (top-bottom).", &
1639 units=
"degC", fail_if_missing=.not.just_read, do_not_log=just_read)
1640 call get_param(param_file, mdl,
"S_TOP", s_top, &
1641 "Initial salinity of the top surface.", &
1642 units=
"PSU", fail_if_missing=.not.just_read, do_not_log=just_read)
1643 call get_param(param_file, mdl,
"S_RANGE", s_range, &
1644 "Initial salinity difference (top-bottom).", &
1645 units=
"PSU", fail_if_missing=.not.just_read, do_not_log=just_read)
1647 if (just_read)
return
1656 s(:,:,k) = s_top - s_range*((real(k)-0.5)/real(g%ke))
1657 t(:,:,k) = t_top - t_range*((real(k)-0.5)/real(g%ke))
1669 call calltree_leave(trim(mdl)//
'()')
1670 end subroutine initialize_temp_salt_linear
1677 subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, param_file, CSp, ALE_CSp, Time)
1678 type(ocean_grid_type),
intent(in) :: G
1679 type(verticalgrid_type),
intent(in) :: GV
1680 type(unit_scale_type),
intent(in) :: US
1681 logical,
intent(in) :: use_temperature
1682 type(thermo_var_ptrs),
intent(in) :: tv
1685 type(sponge_cs),
pointer :: CSp
1687 type(ale_sponge_cs),
pointer :: ALE_CSp
1689 type(time_type),
intent(in) :: Time
1692 real,
allocatable,
dimension(:,:,:) :: eta
1693 real,
allocatable,
dimension(:,:,:) :: h
1695 real,
dimension (SZI_(G),SZJ_(G),SZK_(G)) :: &
1697 real,
dimension (SZI_(G),SZJ_(G)) :: &
1700 real :: Idamp(SZI_(G),SZJ_(G))
1701 real :: pres(SZI_(G))
1703 integer :: i, j, k, is, ie, js, je, nz
1704 integer :: isd, ied, jsd, jed
1705 integer,
dimension(4) :: siz
1707 character(len=40) :: potemp_var, salin_var, Idamp_var, eta_var
1708 character(len=40) :: mdl =
"initialize_sponges_file"
1709 character(len=200) :: damping_file, state_file
1710 character(len=200) :: filename, inputdir
1713 logical :: new_sponges
1716 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1717 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1719 pres(:) = 0.0 ; eta(:,:,:) = 0.0 ; tmp(:,:,:) = 0.0 ; idamp(:,:) = 0.0
1721 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
1722 inputdir = slasher(inputdir)
1723 call get_param(param_file, mdl,
"SPONGE_DAMPING_FILE", damping_file, &
1724 "The name of the file with the sponge damping rates.", &
1725 fail_if_missing=.true.)
1726 call get_param(param_file, mdl,
"SPONGE_STATE_FILE", state_file, &
1727 "The name of the file with the state to damp toward.", &
1728 default=damping_file)
1729 call get_param(param_file, mdl,
"SPONGE_PTEMP_VAR", potemp_var, &
1730 "The name of the potential temperature variable in "//&
1731 "SPONGE_STATE_FILE.", default=
"PTEMP")
1732 call get_param(param_file, mdl,
"SPONGE_SALT_VAR", salin_var, &
1733 "The name of the salinity variable in "//&
1734 "SPONGE_STATE_FILE.", default=
"SALT")
1735 call get_param(param_file, mdl,
"SPONGE_ETA_VAR", eta_var, &
1736 "The name of the interface height variable in "//&
1737 "SPONGE_STATE_FILE.", default=
"ETA")
1738 call get_param(param_file, mdl,
"SPONGE_IDAMP_VAR", idamp_var, &
1739 "The name of the inverse damping rate variable in "//&
1740 "SPONGE_DAMPING_FILE.", default=
"IDAMP")
1741 call get_param(param_file, mdl,
"USE_REGRIDDING", use_ale, do_not_log = .true.)
1743 call get_param(param_file, mdl,
"NEW_SPONGES", new_sponges, &
1744 "Set True if using the newer sponging code which "//&
1745 "performs on-the-fly regridding in lat-lon-time.",&
1746 "of sponge restoring data.", default=.false.)
1754 filename = trim(inputdir)//trim(damping_file)
1755 call log_param(param_file, mdl,
"INPUTDIR/SPONGE_DAMPING_FILE", filename)
1756 if (.not.file_exists(filename, g%Domain)) &
1757 call mom_error(fatal,
" initialize_sponges: Unable to open "//trim(filename))
1759 if (new_sponges .and. .not. use_ale) &
1760 call mom_error(fatal,
" initialize_sponges: Newer sponges are currently unavailable in layered mode ")
1762 call mom_read_data(filename,
"Idamp", idamp(:,:), g%Domain)
1768 filename = trim(inputdir)//trim(state_file)
1769 call log_param(param_file, mdl,
"INPUTDIR/SPONGE_STATE_FILE", filename)
1770 if (.not.file_exists(filename, g%Domain)) &
1771 call mom_error(fatal,
" initialize_sponges: Unable to open "//trim(filename))
1775 if (.not. use_ale)
then
1776 allocate(eta(isd:ied,jsd:jed,nz+1))
1777 call mom_read_data(filename, eta_var, eta(:,:,:), g%Domain, scale=us%m_to_Z)
1779 do j=js,je ;
do i=is,ie
1780 eta(i,j,nz+1) = -g%bathyT(i,j)
1782 do k=nz,1,-1 ;
do j=js,je ;
do i=is,ie
1783 if (eta(i,j,k) < (eta(i,j,k+1) + gv%Angstrom_Z)) &
1784 eta(i,j,k) = eta(i,j,k+1) + gv%Angstrom_Z
1785 enddo ;
enddo ;
enddo
1788 call initialize_sponge(idamp, eta, g, param_file, csp, gv)
1790 elseif (.not. new_sponges)
then
1792 call field_size(filename,eta_var,siz,no_domain=.true.)
1793 if (siz(1) /= g%ieg-g%isg+1 .or. siz(2) /= g%jeg-g%jsg+1) &
1794 call mom_error(fatal,
"initialize_sponge_file: Array size mismatch for sponge data.")
1799 allocate(eta(isd:ied,jsd:jed,nz_data+1))
1800 allocate(h(isd:ied,jsd:jed,nz_data))
1802 call mom_read_data(filename, eta_var, eta(:,:,:), g%Domain, scale=us%m_to_Z)
1804 do j=js,je ;
do i=is,ie
1805 eta(i,j,nz+1) = -g%bathyT(i,j)
1808 do k=nz,1,-1 ;
do j=js,je ;
do i=is,ie
1809 if (eta(i,j,k) < (eta(i,j,k+1) + gv%Angstrom_Z)) &
1810 eta(i,j,k) = eta(i,j,k+1) + gv%Angstrom_Z
1811 enddo ;
enddo ;
enddo
1812 do k=1,nz;
do j=js,je ;
do i=is,ie
1813 h(i,j,k) = gv%Z_to_H*(eta(i,j,k)-eta(i,j,k+1))
1814 enddo ;
enddo ;
enddo
1815 call initialize_ale_sponge(idamp, g, param_file, ale_csp, h, nz_data)
1820 call initialize_ale_sponge(idamp, g, param_file, ale_csp)
1827 if ( gv%nkml>0 .and. .not. new_sponges)
then
1831 do i=is-1,ie ; pres(i) = tv%P_Ref ;
enddo
1833 call mom_read_data(filename, potemp_var, tmp(:,:,:), g%Domain)
1834 call mom_read_data(filename, salin_var, tmp2(:,:,:), g%Domain)
1837 call calculate_density(tmp(:,j,1), tmp2(:,j,1), pres, tmp_2d(:,j), &
1838 is, ie-is+1, tv%eqn_of_state)
1841 call set_up_sponge_ml_density(tmp_2d, g, csp)
1845 if ( use_temperature .and. .not. new_sponges)
then
1846 call mom_read_data(filename, potemp_var, tmp(:,:,:), g%Domain)
1847 call set_up_sponge_field(tmp, tv%T, g, nz, csp)
1848 call mom_read_data(filename, salin_var, tmp(:,:,:), g%Domain)
1849 call set_up_sponge_field(tmp, tv%S, g, nz, csp)
1850 elseif (use_temperature)
then
1851 call set_up_ale_sponge_field(filename, potemp_var, time, g, gv, tv%T, ale_csp)
1852 call set_up_ale_sponge_field(filename, salin_var, time, g, gv, tv%S, ale_csp)
1855 end subroutine initialize_sponges_file
1859 subroutine set_velocity_depth_max(G)
1860 type(ocean_grid_type),
intent(inout) :: G
1864 do i=g%isd,g%ied-1 ;
do j=g%jsd,g%jed
1865 g%Dblock_u(i,j) = g%mask2dCu(i,j) * max(g%bathyT(i,j), g%bathyT(i+1,j))
1866 g%Dopen_u(i,j) = g%Dblock_u(i,j)
1868 do i=g%isd,g%ied ;
do j=g%jsd,g%jed-1
1869 g%Dblock_v(i,j) = g%mask2dCv(i,j) * max(g%bathyT(i,j), g%bathyT(i,j+1))
1870 g%Dopen_v(i,j) = g%Dblock_v(i,j)
1872 end subroutine set_velocity_depth_max
1876 subroutine compute_global_grid_integrals(G)
1877 type(ocean_grid_type),
intent(inout) :: G
1879 real,
dimension(G%isc:G%iec, G%jsc:G%jec) :: tmpForSumming
1882 tmpforsumming(:,:) = 0.
1883 g%areaT_global = 0.0 ; g%IareaT_global = 0.0
1884 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
1885 tmpforsumming(i,j) = g%areaT(i,j) * g%mask2dT(i,j)
1887 g%areaT_global = reproducing_sum(tmpforsumming)
1888 g%IareaT_global = 1. / g%areaT_global
1889 end subroutine compute_global_grid_integrals
1893 subroutine set_velocity_depth_min(G)
1894 type(ocean_grid_type),
intent(inout) :: G
1898 do i=g%isd,g%ied-1 ;
do j=g%jsd,g%jed
1899 g%Dblock_u(i,j) = g%mask2dCu(i,j) * min(g%bathyT(i,j), g%bathyT(i+1,j))
1900 g%Dopen_u(i,j) = g%Dblock_u(i,j)
1902 do i=g%isd,g%ied ;
do j=g%jsd,g%jed-1
1903 g%Dblock_v(i,j) = g%mask2dCv(i,j) * min(g%bathyT(i,j), g%bathyT(i,j+1))
1904 g%Dopen_v(i,j) = g%Dblock_v(i,j)
1906 end subroutine set_velocity_depth_min
1911 subroutine mom_temp_salt_initialize_from_z(h, tv, G, GV, US, PF, just_read_params)
1912 type(ocean_grid_type),
intent(inout) :: G
1913 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1915 type(thermo_var_ptrs),
intent(inout) :: tv
1917 type(verticalgrid_type),
intent(in) :: GV
1918 type(unit_scale_type),
intent(in) :: US
1921 logical,
optional,
intent(in) :: just_read_params
1925 character(len=200) :: filename
1927 character(len=200) :: tfilename
1929 character(len=200) :: sfilename
1931 character(len=200) :: shelf_file
1932 character(len=200) :: inputdir
1933 character(len=200) :: mesg, area_varname, ice_shelf_file
1935 type(eos_type),
pointer :: eos => null()
1936 type(thermo_var_ptrs) :: tv_loc
1937 type(verticalgrid_type) :: GV_loc
1939 # include "version_variable.h"
1940 character(len=40) :: mdl =
"MOM_initialize_layers_from_Z"
1942 integer :: is, ie, js, je, nz
1943 integer :: isc,iec,jsc,jec
1944 integer :: isg, ieg, jsg, jeg
1945 integer :: isd, ied, jsd, jed
1947 integer :: i, j, k, ks, np, ni, nj
1948 integer :: idbg, jdbg
1949 integer :: nkml, nkbl
1951 integer :: kd, inconsistent
1956 real,
dimension(:,:),
pointer :: shelf_area => null()
1959 real :: missing_value_temp, missing_value_salt
1960 logical :: correct_thickness
1961 character(len=40) :: potemp_var, salin_var
1962 character(len=8) :: laynum
1964 integer,
parameter :: niter=10
1965 logical :: just_read
1966 logical :: adjust_temperature = .true.
1967 real,
parameter :: missing_value = -1.e20
1968 real,
parameter :: temp_land_fill = 0.0, salt_land_fill = 35.0
1969 logical :: reentrant_x, tripolar_n,dbg
1970 logical :: debug = .false.
1973 real,
dimension(:),
allocatable :: z_edges_in, z_in, Rb
1974 real,
dimension(:,:,:),
allocatable,
target :: temp_z, salt_z, mask_z
1975 real,
dimension(:,:,:),
allocatable :: rho_z
1976 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: zi
1977 real,
dimension(SZI_(G),SZJ_(G)) :: nlevs
1978 real,
dimension(SZI_(G)) :: press
1981 real,
dimension(:),
allocatable :: hTarget
1982 real,
dimension(:,:),
allocatable :: area_shelf_h
1983 real,
dimension(:,:),
allocatable,
target :: frac_shelf_h
1984 real,
dimension(:,:,:),
allocatable,
target :: tmpT1dIn, tmpS1dIn
1985 real,
dimension(:,:,:),
allocatable :: tmp_mask_in
1986 real,
dimension(:,:,:),
allocatable :: h1
1987 real,
dimension(:,:,:),
allocatable :: dz_interface
1988 real :: zTopOfCell, zBottomOfCell
1989 type(regridding_cs) :: regridCS
1990 type(remapping_cs) :: remapCS
1992 logical :: homogenize, useALEremapping, remap_full_column, remap_general, remap_old_alg
1993 logical :: use_ice_shelf
1994 character(len=10) :: remappingScheme
1995 real :: tempAvg, saltAvg
1996 integer :: nPoints, ans
1997 integer :: id_clock_routine, id_clock_read, id_clock_interp, id_clock_fill, id_clock_ALE
1999 id_clock_routine = cpu_clock_id(
'(Initialize from Z)', grain=clock_routine)
2000 id_clock_ale = cpu_clock_id(
'(Initialize from Z) ALE', grain=clock_loop)
2002 call cpu_clock_begin(id_clock_routine)
2004 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
2005 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
2006 isg = g%isg ; ieg = g%ieg ; jsg = g%jsg ; jeg = g%jeg
2008 pi_180=atan(1.0)/45.
2010 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
2012 if (.not.just_read)
call calltree_enter(trim(mdl)//
"(), MOM_state_initialization.F90")
2013 if (.not.just_read)
call log_version(pf, mdl, version,
"")
2015 inputdir =
"." ;
call get_param(pf, mdl,
"INPUTDIR", inputdir)
2016 inputdir = slasher(inputdir)
2018 eos => tv%eqn_of_state
2022 reentrant_x = .false. ;
call get_param(pf, mdl,
"REENTRANT_X", reentrant_x,default=.true.)
2023 tripolar_n = .false. ;
call get_param(pf, mdl,
"TRIPOLAR_N", tripolar_n, default=.false.)
2024 call get_param(pf, mdl,
"MINIMUM_DEPTH", min_depth, default=0.0, scale=us%m_to_Z)
2026 call get_param(pf, mdl,
"NKML",nkml,default=0)
2027 call get_param(pf, mdl,
"NKBL",nkbl,default=0)
2029 call get_param(pf, mdl,
"TEMP_SALT_Z_INIT_FILE",filename, &
2030 "The name of the z-space input file used to initialize "//&
2031 "temperatures (T) and salinities (S). If T and S are not "//&
2032 "in the same file, TEMP_Z_INIT_FILE and SALT_Z_INIT_FILE "//&
2033 "must be set.",default=
"temp_salt_z.nc",do_not_log=just_read)
2034 call get_param(pf, mdl,
"TEMP_Z_INIT_FILE",tfilename, &
2035 "The name of the z-space input file used to initialize "//&
2036 "temperatures, only.", default=trim(filename),do_not_log=just_read)
2037 call get_param(pf, mdl,
"SALT_Z_INIT_FILE",sfilename, &
2038 "The name of the z-space input file used to initialize "//&
2039 "temperatures, only.", default=trim(filename),do_not_log=just_read)
2040 filename = trim(inputdir)//trim(filename)
2041 tfilename = trim(inputdir)//trim(tfilename)
2042 sfilename = trim(inputdir)//trim(sfilename)
2043 call get_param(pf, mdl,
"Z_INIT_FILE_PTEMP_VAR", potemp_var, &
2044 "The name of the potential temperature variable in "//&
2045 "TEMP_Z_INIT_FILE.", default=
"ptemp",do_not_log=just_read)
2046 call get_param(pf, mdl,
"Z_INIT_FILE_SALT_VAR", salin_var, &
2047 "The name of the salinity variable in "//&
2048 "SALT_Z_INIT_FILE.", default=
"salt",do_not_log=just_read)
2049 call get_param(pf, mdl,
"Z_INIT_HOMOGENIZE", homogenize, &
2050 "If True, then horizontally homogenize the interpolated "//&
2051 "initial conditions.", default=.false., do_not_log=just_read)
2052 call get_param(pf, mdl,
"Z_INIT_ALE_REMAPPING", usealeremapping, &
2053 "If True, then remap straight to model coordinate from file.", &
2054 default=.false., do_not_log=just_read)
2055 call get_param(pf, mdl,
"Z_INIT_REMAPPING_SCHEME", remappingscheme, &
2056 "The remapping scheme to use if using Z_INIT_ALE_REMAPPING "//&
2057 "is True.", default=
"PPM_IH4", do_not_log=just_read)
2058 call get_param(pf, mdl,
"Z_INIT_REMAP_GENERAL", remap_general, &
2059 "If false, only initializes to z* coordinates. "//&
2060 "If true, allows initialization directly to general coordinates.",&
2061 default=.false., do_not_log=just_read)
2062 call get_param(pf, mdl,
"Z_INIT_REMAP_FULL_COLUMN", remap_full_column, &
2063 "If false, only reconstructs profiles for valid data points. "//&
2064 "If true, inserts vanished layers below the valid data.", &
2065 default=remap_general, do_not_log=just_read)
2066 call get_param(pf, mdl,
"Z_INIT_REMAP_OLD_ALG", remap_old_alg, &
2067 "If false, uses the preferred remapping algorithm for initialization. "//&
2068 "If true, use an older, less robust algorithm for remapping.", &
2069 default=.true., do_not_log=just_read)
2070 call get_param(pf, mdl,
"ICE_SHELF", use_ice_shelf, default=.false.)
2071 if (use_ice_shelf)
then
2072 call get_param(pf, mdl,
"ICE_THICKNESS_FILE", ice_shelf_file, &
2073 "The file from which the ice bathymetry and area are read.", &
2074 fail_if_missing=.not.just_read, do_not_log=just_read)
2075 shelf_file = trim(inputdir)//trim(ice_shelf_file)
2076 if (.not.just_read)
call log_param(pf, mdl,
"INPUTDIR/THICKNESS_FILE", shelf_file)
2077 call get_param(pf, mdl,
"ICE_AREA_VARNAME", area_varname, &
2078 "The name of the area variable in ICE_THICKNESS_FILE.", &
2079 fail_if_missing=.not.just_read, do_not_log=just_read)
2081 if (.not.usealeremapping)
then
2082 call get_param(pf, mdl,
"ADJUST_THICKNESS", correct_thickness, &
2083 "If true, all mass below the bottom removed if the "//&
2084 "topography is shallower than the thickness input file "//&
2085 "would indicate.", default=.false., do_not_log=just_read)
2087 call get_param(pf, mdl,
"FIT_TO_TARGET_DENSITY_IC", adjust_temperature, &
2088 "If true, all the interior layers are adjusted to "//&
2089 "their target densities using mostly temperature "//&
2090 "This approach can be problematic, particularly in the "//&
2091 "high latitudes.", default=.true., do_not_log=just_read)
2094 call cpu_clock_end(id_clock_routine)
2099 eps_z = 1.0e-10*us%m_to_Z
2116 call horiz_interp_and_extrap_tracer(tfilename, potemp_var, 1.0, 1, &
2117 g, temp_z, mask_z, z_in, z_edges_in, missing_value_temp, reentrant_x, &
2118 tripolar_n, homogenize, m_to_z=us%m_to_Z)
2120 call horiz_interp_and_extrap_tracer(sfilename, salin_var, 1.0, 1, &
2121 g, salt_z, mask_z, z_in, z_edges_in, missing_value_salt, reentrant_x, &
2122 tripolar_n, homogenize, m_to_z=us%m_to_Z)
2127 do k=1,
size(z_edges_in,1) ; z_edges_in(k) = -z_edges_in(k) ;
enddo
2129 allocate(rho_z(isd:ied,jsd:jed,kd))
2130 allocate(area_shelf_h(isd:ied,jsd:jed))
2131 allocate(frac_shelf_h(isd:ied,jsd:jed))
2136 call convert_temp_salt_for_teos10(temp_z, salt_z, press, g, kd, mask_z, eos)
2138 do k=1,kd ;
do j=js,je
2139 call calculate_density(temp_z(:,j,k), salt_z(:,j,k), press, rho_z(:,j,k), is, ie, eos)
2148 if (use_ice_shelf)
then
2149 if (.not.file_exists(shelf_file, g%Domain))
call mom_error(fatal, &
2150 "MOM_temp_salt_initialize_from_Z: Unable to open shelf file "//trim(shelf_file))
2152 call mom_read_data(shelf_file, trim(area_varname), area_shelf_h, g%Domain)
2155 frac_shelf_h(:,:) = 0.0
2157 do j=jsd,jed ;
do i=isd,ied
2158 if (g%areaT(i,j) > 0.0) &
2159 frac_shelf_h(i,j) = area_shelf_h(i,j) / g%areaT(i,j)
2162 shelf_area => frac_shelf_h
2168 if (usealeremapping)
then
2169 call cpu_clock_begin(id_clock_ale)
2170 nkd = max(gv%ke, kd)
2177 allocate( tmp_mask_in(isd:ied,jsd:jed,nkd) ) ; tmp_mask_in(:,:,:) = 0.
2178 allocate( h1(isd:ied,jsd:jed,nkd) ) ; h1(:,:,:) = 0.
2179 allocate( tmpt1din(isd:ied,jsd:jed,nkd) ) ; tmpt1din(:,:,:) = 0.
2180 allocate( tmps1din(isd:ied,jsd:jed,nkd) ) ; tmps1din(:,:,:) = 0.
2181 do j = js, je ;
do i = is, ie
2182 if (g%mask2dT(i,j)>0.)
then
2183 ztopofcell = 0. ; zbottomofcell = 0.
2184 tmp_mask_in(i,j,1:kd) = mask_z(i,j,:)
2186 if (tmp_mask_in(i,j,k)>0. .and. k<=kd)
then
2187 zbottomofcell = max( z_edges_in(k+1), -g%bathyT(i,j) )
2188 tmpt1din(i,j,k) = temp_z(i,j,k)
2189 tmps1din(i,j,k) = salt_z(i,j,k)
2191 zbottomofcell = -g%bathyT(i,j)
2192 tmpt1din(i,j,k) = tmpt1din(i,j,k-1)
2193 tmps1din(i,j,k) = tmps1din(i,j,k-1)
2195 tmpt1din(i,j,k) = -99.9
2196 tmps1din(i,j,k) = -99.9
2198 h1(i,j,k) = gv%Z_to_H * (ztopofcell - zbottomofcell)
2199 ztopofcell = zbottomofcell
2201 h1(i,j,kd) = h1(i,j,kd) + gv%Z_to_H * max(0., ztopofcell + g%bathyT(i,j) )
2205 deallocate( tmp_mask_in )
2212 call ale_initregridding( gv, us, g%max_depth, pf, mdl, regridcs )
2214 if (.not. remap_general)
then
2216 allocate( htarget(nz) )
2217 htarget = getcoordinateresolution( regridcs )
2218 do j = js, je ;
do i = is, ie
2220 if (g%mask2dT(i,j)>0.)
then
2222 ztopofcell = 0. ; zbottomofcell = 0.
2224 zbottomofcell = max( ztopofcell - htarget(k), -g%bathyT(i,j) )
2225 h(i,j,k) = gv%Z_to_H * (ztopofcell - zbottomofcell)
2226 ztopofcell = zbottomofcell
2233 deallocate( htarget )
2237 call initialize_remapping( remapcs, remappingscheme, boundary_extrapolation=.false. )
2238 if (remap_general)
then
2239 call set_regrid_params( regridcs, min_thickness=0. )
2241 tv_loc%T => tmpt1din
2242 tv_loc%S => tmps1din
2245 allocate( dz_interface(isd:ied,jsd:jed,nkd+1) )
2246 if (use_ice_shelf)
then
2247 call regridding_main( remapcs, regridcs, g, gv_loc, h1, tv_loc, h, dz_interface, shelf_area )
2249 call regridding_main( remapcs, regridcs, g, gv_loc, h1, tv_loc, h, dz_interface )
2251 deallocate( dz_interface )
2253 call ale_remap_scalar(remapcs, g, gv, nkd, h1, tmpt1din, h, tv%T, all_cells=remap_full_column, &
2254 old_remap=remap_old_alg )
2255 call ale_remap_scalar(remapcs, g, gv, nkd, h1, tmps1din, h, tv%S, all_cells=remap_full_column, &
2256 old_remap=remap_old_alg )
2258 deallocate( tmpt1din )
2259 deallocate( tmps1din )
2261 call cpu_clock_end(id_clock_ale)
2267 nlevs = sum(mask_z,dim=3)
2271 do k=2,nz ; rb(k)=0.5*(gv%Rlay(k-1)+gv%Rlay(k)) ;
enddo
2272 rb(1) = 0.0 ; rb(nz+1) = 2.0*gv%Rlay(nz) - gv%Rlay(nz-1)
2274 zi(is:ie,js:je,:) = find_interfaces(rho_z(is:ie,js:je,:), z_in, rb, g%bathyT(is:ie,js:je), &
2275 nlevs(is:ie,js:je), nkml, nkbl, min_depth, eps_z=eps_z)
2277 if (correct_thickness)
then
2278 call adjustetatofitbathymetry(g, gv, us, zi, h)
2280 do k=nz,1,-1 ;
do j=js,je ;
do i=is,ie
2281 if (zi(i,j,k) < (zi(i,j,k+1) + gv%Angstrom_Z))
then
2282 zi(i,j,k) = zi(i,j,k+1) + gv%Angstrom_Z
2283 h(i,j,k) = gv%Angstrom_H
2285 h(i,j,k) = gv%Z_to_H * (zi(i,j,k) - zi(i,j,k+1))
2287 enddo ;
enddo ;
enddo
2289 do j=js,je ;
do i=is,ie
2290 if (abs(zi(i,j,nz+1) + g%bathyT(i,j)) > 1.0*us%m_to_Z) &
2291 inconsistent = inconsistent + 1
2293 call sum_across_pes(inconsistent)
2295 if ((inconsistent > 0) .and. (is_root_pe()))
then
2296 write(mesg,
'("Thickness initial conditions are inconsistent ",'// &
2297 '"with topography in ",I5," places.")') inconsistent
2298 call mom_error(warning, mesg)
2302 tv%T(is:ie,js:je,:) = tracer_z_init(temp_z(is:ie,js:je,:), z_edges_in, zi(is:ie,js:je,:), &
2303 nkml, nkbl, missing_value, g%mask2dT(is:ie,js:je), nz, &
2304 nlevs(is:ie,js:je),dbg,idbg,jdbg, eps_z=eps_z)
2305 tv%S(is:ie,js:je,:) = tracer_z_init(salt_z(is:ie,js:je,:), z_edges_in, zi(is:ie,js:je,:), &
2306 nkml, nkbl, missing_value, g%mask2dT(is:ie,js:je), nz, &
2307 nlevs(is:ie,js:je), eps_z=eps_z)
2310 npoints = 0 ; tempavg = 0. ; saltavg = 0.
2311 do j=js,je ;
do i=is,ie ;
if (g%mask2dT(i,j) >= 1.0)
then
2312 npoints = npoints + 1
2313 tempavg = tempavg + tv%T(i,j,k)
2314 saltavg = saltavg + tv%S(i,j,k)
2315 endif ;
enddo ;
enddo
2318 if (homogenize)
then
2319 call sum_across_pes(npoints)
2320 call sum_across_pes(tempavg)
2321 call sum_across_pes(saltavg)
2323 tempavg = tempavg / real(npoints)
2324 saltavg = saltavg / real(npoints)
2326 tv%T(:,:,k) = tempavg
2327 tv%S(:,:,k) = saltavg
2334 do k=1,nz ;
do j=js,je ;
do i=is,ie
2335 if (tv%T(i,j,k) == missing_value)
then
2336 tv%T(i,j,k) = temp_land_fill
2337 tv%S(i,j,k) = salt_land_fill
2339 enddo ;
enddo ;
enddo
2342 ks = max(0,nkml)+max(0,nkbl)+1
2344 if (adjust_temperature .and. .not. usealeremapping)
then
2345 call determine_temperature(tv%T(is:ie,js:je,:), tv%S(is:ie,js:je,:), &
2346 gv%Rlay(1:nz), tv%p_ref, niter, missing_value, h(is:ie,js:je,:), ks, eos)
2350 deallocate(z_in, z_edges_in, temp_z, salt_z, mask_z)
2351 deallocate(rho_z, area_shelf_h, frac_shelf_h)
2357 call calltree_leave(trim(mdl)//
'()')
2358 call cpu_clock_end(id_clock_routine)
2360 end subroutine mom_temp_salt_initialize_from_z
2363 subroutine mom_state_init_tests(G, GV, US, tv)
2364 type(ocean_grid_type),
intent(inout) :: G
2365 type(verticalgrid_type),
intent(in) :: GV
2366 type(unit_scale_type),
intent(in) :: US
2367 type(thermo_var_ptrs),
intent(in) :: tv
2369 integer,
parameter :: nk=5
2370 real,
dimension(nk) :: T, T_t, T_b, S, S_t, S_b, rho, h, z
2371 real,
dimension(nk+1) :: e
2373 real :: P_tot, P_t, P_b, z_out
2374 type(remapping_cs),
pointer :: remap_CS => null()
2381 e(k+1) = e(k) - h(k)
2385 z(k) = 0.5 * ( e(k) + e(k+1) )
2386 t_t(k) = 20.+(0./500.)*e(k)
2387 t(k) = 20.+(0./500.)*z(k)
2388 t_b(k) = 20.+(0./500.)*e(k+1)
2389 s_t(k) = 35.-(0./500.)*e(k)
2390 s(k) = 35.+(0./500.)*z(k)
2391 s_b(k) = 35.-(0./500.)*e(k+1)
2392 call calculate_density(0.5*(t_t(k)+t_b(k)), 0.5*(s_t(k)+s_b(k)), -gv%Rho0*gv%mks_g_Earth*z(k), &
2393 rho(k), tv%eqn_of_state)
2394 p_tot = p_tot + gv%mks_g_Earth * rho(k) * h(k)
2399 call find_depth_of_pressure_in_cell(t_t(k), t_b(k), s_t(k), s_b(k), e(k), e(k+1), &
2400 p_t, 0.5*p_tot, gv%Rho0, gv%mks_g_Earth, tv%eqn_of_state, p_b, z_out)
2401 write(0,*) k,p_t,p_b,0.5*p_tot,e(k),e(k+1),z_out
2404 write(0,*) p_b,p_tot
2407 write(0,*)
' ==================================================================== '
2410 call cut_off_column_top(nk, tv, gv, gv%mks_g_Earth, -e(nk+1), gv%Angstrom_H, &
2411 t, t_t, t_b, s, s_t, s_b, 0.5*p_tot, h, remap_cs)
2414 end subroutine mom_state_init_tests