8 use mom_cpu_clock,
only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
103 use fms_io_mod,
only : field_size
105 implicit none ;
private
107 #include <MOM_memory.h>
116 character(len=40) ::
mdl =
"MOM_state_initialization"
122 subroutine mom_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
123 restart_CS, ALE_CSp, tracer_Reg, sponge_CSp, &
124 ALE_sponge_CSp, OBC, Time_in)
128 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
131 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
134 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
138 type(time_type),
intent(inout) :: time
145 type(
ale_cs),
pointer :: ale_csp
150 type(time_type),
optional,
intent(in) :: time_in
153 character(len=200) :: filename
154 character(len=200) :: filename2
155 character(len=200) :: inputdir
156 character(len=200) :: config
162 logical :: from_z_file, useale
164 integer :: write_geom
165 logical :: use_temperature, use_sponge, use_obc
167 logical :: depress_sfc
169 logical :: trim_ic_for_p_surf
171 logical :: regrid_accelerate
172 integer :: regrid_iterations
178 type(
eos_type),
pointer :: eos => null()
181 logical :: debug_layers = .false.
182 character(len=80) :: mesg
184 #include "version_variable.h"
185 integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz
186 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
188 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
189 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
190 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
191 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
193 call calltree_enter(
"MOM_initialize_state(), MOM_state_initialization.F90")
195 call get_param(pf,
mdl,
"DEBUG", debug, default=.false.)
196 call get_param(pf,
mdl,
"DEBUG_OBC", debug_obc, default=.false.)
200 just_read = .not.new_sim
203 "The directory in which input files are found.", default=
".")
204 inputdir = slasher(inputdir)
206 use_temperature =
associated(tv%T)
207 useale =
associated(ale_csp)
208 use_eos =
associated(tv%eqn_of_state)
209 use_obc =
associated(obc)
210 if (use_eos) eos => tv%eqn_of_state
218 call mom_mesg(
"Run initialized internally.", 3)
220 if (
present(time_in)) time = time_in
235 call get_param(pf,
mdl,
"INIT_LAYERS_FROM_Z_FILE", from_z_file, &
236 "If true, initialize the layer thicknesses, temperatures, "//&
237 "and salinities from a Z-space file on a latitude-longitude "//&
238 "grid.", default=.false., do_not_log=just_read)
240 if (from_z_file)
then
242 if (.NOT.use_temperature)
call mom_error(fatal,
"MOM_initialize_state : "//&
243 "use_temperature must be true if INIT_LAYERS_FROM_Z_FILE is true")
250 "A string that determines how the initial layer "//&
251 "thicknesses are specified for a new run: \n"//&
252 " \t file - read interface heights from the file specified \n"//&
253 " \t thickness_file - read thicknesses from the file specified \n"//&
254 " \t\t by (THICKNESS_FILE).\n"//&
255 " \t coord - determined by ALE coordinate.\n"//&
256 " \t uniform - uniform thickness layers evenly distributed \n"//&
257 " \t\t between the surface and MAXIMUM_DEPTH. \n"//&
258 " \t list - read a list of positive interface depths. \n"//&
259 " \t DOME - use a slope and channel configuration for the \n"//&
260 " \t\t DOME sill-overflow test case. \n"//&
261 " \t ISOMIP - use a configuration for the \n"//&
262 " \t\t ISOMIP test case. \n"//&
263 " \t benchmark - use the benchmark test case thicknesses. \n"//&
264 " \t Neverland - use the Neverland test case thicknesses. \n"//&
265 " \t search - search a density profile for the interface \n"//&
266 " \t\t densities. This is not yet implemented. \n"//&
267 " \t circle_obcs - the circle_obcs test case is used. \n"//&
268 " \t DOME2D - 2D version of DOME initialization. \n"//&
269 " \t adjustment2d - 2D lock exchange thickness ICs. \n"//&
270 " \t sloshing - sloshing gravity thickness ICs. \n"//&
271 " \t seamount - no motion test with seamount ICs. \n"//&
272 " \t dumbbell - sloshing channel ICs. \n"//&
273 " \t soliton - Equatorial Rossby soliton. \n"//&
274 " \t rossby_front - a mixed layer front in thermal wind balance.\n"//&
275 " \t USER - call a user modified routine.", &
276 fail_if_missing=new_sim, do_not_log=just_read)
277 select case (trim(config))
280 case (
"thickness_file")
283 if (new_sim .and. useale)
then
284 call ale_initthicknesstocoord( ale_csp, g, gv, h )
285 elseif (new_sim)
then
286 call mom_error(fatal,
"MOM_initialize_state: USE_REGRIDDING must be True "//&
287 "for THICKNESS_CONFIG of 'coord'")
290 just_read_params=just_read)
292 just_read_params=just_read)
293 case (
"DOME");
call dome_initialize_thickness(h, g, gv, pf, &
294 just_read_params=just_read)
295 case (
"ISOMIP");
call isomip_initialize_thickness(h, g, gv, us, pf, tv, &
296 just_read_params=just_read)
297 case (
"benchmark");
call benchmark_initialize_thickness(h, g, gv, us, pf, &
298 tv%eqn_of_state, tv%P_Ref, just_read_params=just_read)
300 tv%eqn_of_state, tv%P_Ref)
303 just_read_params=just_read)
305 pf, just_read_params=just_read)
307 pf, just_read_params=just_read)
308 case (
"DOME2D");
call dome2d_initialize_thickness(h, g, gv, us, pf, &
309 just_read_params=just_read)
310 case (
"adjustment2d");
call adjustment_initialize_thickness(h, g, gv, us, &
311 pf, just_read_params=just_read)
312 case (
"sloshing");
call sloshing_initialize_thickness(h, g, gv, us, pf, &
313 just_read_params=just_read)
314 case (
"seamount");
call seamount_initialize_thickness(h, g, gv, us, pf, &
315 just_read_params=just_read)
316 case (
"dumbbell");
call dumbbell_initialize_thickness(h, g, gv, us, pf, &
317 just_read_params=just_read)
319 case (
"phillips");
call phillips_initialize_thickness(h, g, gv, us, pf, &
320 just_read_params=just_read)
321 case (
"rossby_front");
call rossby_front_initialize_thickness(h, g, gv, us, &
322 pf, just_read_params=just_read)
323 case (
"USER");
call user_initialize_thickness(h, g, gv, pf, &
324 just_read_params=just_read)
325 case default ;
call mom_error(fatal,
"MOM_initialize_state: "//&
326 "Unrecognized layer thickness configuration "//trim(config))
330 if ( use_temperature )
then
332 "A string that determines how the initial tempertures "//&
333 "and salinities are specified for a new run: \n"//&
334 " \t file - read velocities from the file specified \n"//&
335 " \t\t by (TS_FILE). \n"//&
336 " \t fit - find the temperatures that are consistent with \n"//&
337 " \t\t the layer densities and salinity S_REF. \n"//&
338 " \t TS_profile - use temperature and salinity profiles \n"//&
339 " \t\t (read from TS_FILE) to set layer densities. \n"//&
340 " \t benchmark - use the benchmark test case T & S. \n"//&
341 " \t linear - linear in logical layer space. \n"//&
342 " \t DOME2D - 2D DOME initialization. \n"//&
343 " \t ISOMIP - ISOMIP initialization. \n"//&
344 " \t adjustment2d - 2d lock exchange T/S ICs. \n"//&
345 " \t sloshing - sloshing mode T/S ICs. \n"//&
346 " \t seamount - no motion test with seamount ICs. \n"//&
347 " \t dumbbell - sloshing channel ICs. \n"//&
348 " \t rossby_front - a mixed layer front in thermal wind balance.\n"//&
349 " \t SCM_CVMix_tests - used in the SCM CVMix tests.\n"//&
350 " \t USER - call a user modified routine.", &
351 fail_if_missing=new_sim, do_not_log=just_read)
353 select case (trim(config))
355 eos, tv%P_Ref, just_read_params=just_read)
357 pf, just_read_params=just_read)
359 g, gv, us, pf, eos, tv%P_Ref, just_read_params=just_read)
361 g, pf, just_read_params=just_read)
363 just_read_params=just_read)
364 case (
"DOME2D");
call dome2d_initialize_temperature_salinity ( tv%T, &
365 tv%S, h, g, gv, pf, eos, just_read_params=just_read)
367 tv%S, h, g, gv, us, pf, eos, just_read_params=just_read)
369 tv%S, h, g, gv, pf, eos, just_read_params=just_read)
371 tv%S, h, g, gv, us, pf, just_read_params=just_read)
373 tv%S, h, g, gv, pf, eos, just_read_params=just_read)
375 tv%S, h, g, gv, pf, eos, just_read_params=just_read)
376 case (
"dumbbell");
call dumbbell_initialize_temperature_salinity(tv%T, &
377 tv%S, h, g, gv, pf, eos, just_read_params=just_read)
378 case (
"rossby_front");
call rossby_front_initialize_temperature_salinity ( tv%T, &
379 tv%S, h, g, gv, pf, eos, just_read_params=just_read)
381 g, gv, us, pf, just_read_params=just_read)
382 case (
"dense");
call dense_water_initialize_ts(g, gv, pf, eos, tv%T, tv%S, &
383 h, just_read_params=just_read)
384 case (
"USER");
call user_init_temperature_salinity(tv%T, tv%S, g, pf, eos, &
385 just_read_params=just_read)
386 case default ;
call mom_error(fatal,
"MOM_initialize_state: "//&
387 "Unrecognized Temp & salt configuration "//trim(config))
391 if (use_temperature .and. use_obc) &
392 call fill_temp_salt_segments(g, obc, tv)
395 if (new_sim)
call pass_var(h, g%Domain)
399 "A string that determines how the initial velocities "//&
400 "are specified for a new run: \n"//&
401 " \t file - read velocities from the file specified \n"//&
402 " \t\t by (VELOCITY_FILE). \n"//&
403 " \t zero - the fluid is initially at rest. \n"//&
404 " \t uniform - the flow is uniform (determined by\n"//&
405 " \t\t parameters INITIAL_U_CONST and INITIAL_V_CONST).\n"//&
406 " \t rossby_front - a mixed layer front in thermal wind balance.\n"//&
407 " \t soliton - Equatorial Rossby soliton.\n"//&
408 " \t USER - call a user modified routine.", default=
"zero", &
409 do_not_log=just_read)
410 select case (trim(config))
412 just_read_params=just_read)
414 just_read_params=just_read)
416 just_read_params=just_read)
418 just_read_params=just_read)
419 case (
"phillips");
call phillips_initialize_velocity(u, v, g, gv, us, pf, &
420 just_read_params=just_read)
422 g, gv, us, pf, just_read_params=just_read)
423 case (
"soliton");
call soliton_initialize_velocity(u, v, h, g, us)
424 case (
"USER");
call user_initialize_velocity(u, v, g, us, pf, &
425 just_read_params=just_read)
426 case default ;
call mom_error(fatal,
"MOM_initialize_state: "//&
427 "Unrecognized velocity configuration "//trim(config))
431 if (debug .and. new_sim)
then
432 call uvchksum(
"MOM_initialize_state [uv]", u, v, g%HI, haloshift=1, scale=us%m_s_to_L_T)
437 call get_param(pf,
mdl,
"CONVERT_THICKNESS_UNITS", convert, &
438 "If true, convert the thickness initial conditions from "//&
439 "units of m to kg m-2 or vice versa, depending on whether "//&
440 "BOUSSINESQ is defined. This does not apply if a restart "//&
441 "file is read.", default=.not.gv%Boussinesq, do_not_log=just_read)
443 if (new_sim .and. convert .and. .not.gv%Boussinesq) &
448 call get_param(pf,
mdl,
"DEPRESS_INITIAL_SURFACE", depress_sfc, &
449 "If true, depress the initial surface to avoid huge "//&
450 "tsunamis when a large surface pressure is applied.", &
451 default=.false., do_not_log=just_read)
452 call get_param(pf,
mdl,
"TRIM_IC_FOR_P_SURF", trim_ic_for_p_surf, &
453 "If true, cuts way the top of the column for initial conditions "//&
454 "at the depth where the hydrostatic pressure matches the imposed "//&
455 "surface pressure which is read from file.", default=.false., &
456 do_not_log=just_read)
457 if (depress_sfc .and. trim_ic_for_p_surf)
call mom_error(fatal,
"MOM_initialize_state: "//&
458 "DEPRESS_INITIAL_SURFACE and TRIM_IC_FOR_P_SURF are exclusive and cannot both be True")
459 if (new_sim .and. debug .and. (depress_sfc .or. trim_ic_for_p_surf)) &
460 call hchksum(h,
"Pre-depress: h ", g%HI, haloshift=1, scale=gv%H_to_m)
461 if (depress_sfc)
call depress_surface(h, g, gv, us, pf, tv, just_read_params=just_read)
462 if (trim_ic_for_p_surf)
call trim_for_ice(pf, g, gv, us, ale_csp, tv, h, just_read_params=just_read)
467 call get_param(pf,
mdl,
"REGRID_ACCELERATE_INIT", regrid_accelerate, &
468 "If true, runs REGRID_ACCELERATE_ITERATIONS iterations of the regridding "//&
469 "algorithm to push the initial grid to be consistent with the initial "//&
470 "condition. Useful only for state-based and iterative coordinates.", &
471 default=.false., do_not_log=just_read)
472 if (regrid_accelerate)
then
473 call get_param(pf,
mdl,
"REGRID_ACCELERATE_ITERATIONS", regrid_iterations, &
474 "The number of regridding iterations to perform to generate "//&
475 "an initial grid that is consistent with the initial conditions.", &
476 default=1, do_not_log=just_read)
478 call get_param(pf,
mdl,
"DT", dt,
"Timestep", fail_if_missing=.true., scale=us%s_to_T)
480 if (new_sim .and. debug) &
481 call hchksum(h,
"Pre-ALE_regrid: h ", g%HI, haloshift=1, scale=gv%H_to_m)
482 call ale_regrid_accelerated(ale_csp, g, gv, h, tv, regrid_iterations, u, v, obc, tracer_reg, &
483 dt=dt, initial=.true.)
489 if (.not.new_sim)
then
492 call restore_state(dirs%input_filename, dirs%restart_input_dir, time, &
494 if (
present(time_in)) time = time_in
495 if ((gv%m_to_H_restart /= 0.0) .and. (gv%m_to_H_restart /= gv%m_to_H))
then
496 h_rescale = gv%m_to_H / gv%m_to_H_restart
497 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
499 if ( (us%s_to_T_restart * us%m_to_L_restart /= 0.0) .and. &
500 ((us%m_to_L * us%s_to_T_restart) /= (us%m_to_L_restart * us%s_to_T)) )
then
501 vel_rescale = (us%m_to_L * us%s_to_T_restart) / (us%m_to_L_restart * us%s_to_T)
502 do k=1,nz ;
do j=jsd,jed ;
do i=isdb,iedb ; u(i,j,k) = vel_rescale * u(i,j,k) ;
enddo ;
enddo ;
enddo
503 do k=1,nz ;
do j=jsdb,jedb ;
do i=isd,ied ; v(i,j,k) = vel_rescale * v(i,j,k) ;
enddo ;
enddo ;
enddo
507 if ( use_temperature )
then
508 call pass_var(tv%T, g%Domain, complete=.false.)
509 call pass_var(tv%S, g%Domain, complete=.false.)
514 call hchksum(h,
"MOM_initialize_state: h ", g%HI, haloshift=1, scale=gv%H_to_m)
515 if ( use_temperature )
call hchksum(tv%T,
"MOM_initialize_state: T ", g%HI, haloshift=1)
516 if ( use_temperature )
call hchksum(tv%S,
"MOM_initialize_state: S ", g%HI, haloshift=1)
517 if ( use_temperature .and. debug_layers)
then ;
do k=1,nz
518 write(mesg,
'("MOM_IS: T[",I2,"]")') k
519 call hchksum(tv%T(:,:,k), mesg, g%HI, haloshift=1)
520 write(mesg,
'("MOM_IS: S[",I2,"]")') k
521 call hchksum(tv%S(:,:,k), mesg, g%HI, haloshift=1)
526 "If true, sponges may be applied anywhere in the domain. "//&
527 "The exact location and properties of those sponges are "//&
528 "specified via SPONGE_CONFIG.", default=.false.)
529 if ( use_sponge )
then
531 "A string that sets how the sponges are configured: \n"//&
532 " \t file - read sponge properties from the file \n"//&
533 " \t\t specified by (SPONGE_FILE).\n"//&
534 " \t ISOMIP - apply ale sponge in the ISOMIP case \n"//&
535 " \t RGC - apply sponge in the rotating_gravity_current case \n"//&
536 " \t DOME - use a slope and channel configuration for the \n"//&
537 " \t\t DOME sill-overflow test case. \n"//&
538 " \t BFB - Sponge at the southern boundary of the domain\n"//&
539 " \t\t for buoyancy-forced basin case.\n"//&
540 " \t USER - call a user modified routine.", default=
"file")
541 select case (trim(config))
544 sponge_csp, ale_sponge_csp)
545 case (
"ISOMIP");
call isomip_initialize_sponges(g, gv, us, tv, pf, useale, &
546 sponge_csp, ale_sponge_csp)
548 sponge_csp, ale_sponge_csp)
553 sponge_csp, ale_sponge_csp)
556 sponge_csp, ale_sponge_csp)
558 sponge_csp, ale_sponge_csp, time)
559 case default ;
call mom_error(fatal,
"MOM_initialize_state: "//&
560 "Unrecognized sponge configuration "//trim(config))
565 call open_boundary_init(g, gv, us, pf, obc, restart_cs)
568 if (
associated(obc))
then
570 if (.not. obc%needs_IO_for_data) &
574 "A string that sets how the user code is invoked to set open boundary data: \n"//&
575 " DOME - specified inflow on northern boundary\n"//&
576 " dyed_channel - supercritical with dye on the inflow boundary\n"//&
577 " dyed_obcs - circle_obcs with dyes on the open boundaries\n"//&
578 " Kelvin - barotropic Kelvin wave forcing on the western boundary\n"//&
579 " shelfwave - Flather with shelf wave forcing on western boundary\n"//&
580 " supercritical - now only needed here for the allocations\n"//&
581 " tidal_bay - Flather with tidal forcing on eastern boundary\n"//&
582 " USER - user specified", default=
"none")
583 if (trim(config) ==
"DOME")
then
584 call dome_set_obc_data(obc, tv, g, gv, us, pf, tracer_reg)
585 elseif (trim(config) ==
"dyed_channel")
then
587 obc%update_OBC = .true.
588 elseif (trim(config) ==
"dyed_obcs")
then
590 elseif (trim(config) ==
"Kelvin")
then
591 obc%update_OBC = .true.
592 elseif (trim(config) ==
"shelfwave")
then
593 obc%update_OBC = .true.
594 elseif (
lowercase(trim(config)) ==
"supercritical")
then
596 elseif (trim(config) ==
"tidal_bay")
then
597 obc%update_OBC = .true.
598 elseif (trim(config) ==
"USER")
then
599 call user_set_obc_data(obc, tv, g, pf, tracer_reg)
600 elseif (.not. trim(config) ==
"none")
then
601 call mom_error(fatal,
"The open boundary conditions specified by "//&
602 "OBC_USER_CONFIG = "//trim(config)//
" have not been fully implemented.")
604 if (open_boundary_query(obc, apply_open_obc=.true.))
then
605 call set_tracer_data(obc, tv, h, g, pf, tracer_reg)
612 if (debug.and.
associated(obc))
then
613 call hchksum(g%mask2dT,
'MOM_initialize_state: mask2dT ', g%HI)
614 call uvchksum(
'MOM_initialize_state: mask2dC[uv]', g%mask2dCu, &
616 call qchksum(g%mask2dBu,
'MOM_initialize_state: mask2dBu ', g%HI)
619 if (debug_obc)
call open_boundary_test_extern_h(g, gv, obc, h)
630 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
634 logical,
intent(in) :: file_has_thickness
637 logical,
optional,
intent(in) :: just_read_params
641 real :: eta(SZI_(G),SZJ_(G),SZK_(G)+1)
642 integer :: inconsistent = 0
643 logical :: correct_thickness
645 character(len=40) :: mdl =
"initialize_thickness_from_file"
646 character(len=200) :: filename, thickness_file, inputdir, mesg
647 integer :: i, j, k, is, ie, js, je, nz
649 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
651 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
653 if (.not.just_read) &
654 call calltree_enter(trim(mdl)//
"(), MOM_state_initialization.F90")
656 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".", do_not_log=just_read)
657 inputdir = slasher(inputdir)
658 call get_param(param_file, mdl,
"THICKNESS_FILE", thickness_file, &
659 "The name of the thickness file.", &
660 fail_if_missing=.not.just_read, do_not_log=just_read)
662 filename = trim(inputdir)//trim(thickness_file)
663 if (.not.just_read)
call log_param(param_file, mdl,
"INPUTDIR/THICKNESS_FILE", filename)
665 if ((.not.just_read) .and. (.not.
file_exists(filename, g%Domain)))
call mom_error(fatal, &
666 " initialize_thickness_from_file: Unable to open "//trim(filename))
668 if (file_has_thickness)
then
670 if (just_read)
return
671 call mom_read_data(filename,
"h", h(:,:,:), g%Domain, scale=gv%m_to_H)
673 call get_param(param_file, mdl,
"ADJUST_THICKNESS", correct_thickness, &
674 "If true, all mass below the bottom removed if the "//&
675 "topography is shallower than the thickness input file "//&
676 "would indicate.", default=.false., do_not_log=just_read)
677 if (just_read)
return
679 call mom_read_data(filename,
"eta", eta(:,:,:), g%Domain, scale=us%m_to_Z)
681 if (correct_thickness)
then
684 do k=nz,1,-1 ;
do j=js,je ;
do i=is,ie
685 if (eta(i,j,k) < (eta(i,j,k+1) + gv%Angstrom_Z))
then
686 eta(i,j,k) = eta(i,j,k+1) + gv%Angstrom_Z
687 h(i,j,k) = gv%Angstrom_H
689 h(i,j,k) = gv%Z_to_H * (eta(i,j,k) - eta(i,j,k+1))
691 enddo ;
enddo ;
enddo
693 do j=js,je ;
do i=is,ie
694 if (abs(eta(i,j,nz+1) + g%bathyT(i,j)) > 1.0*us%m_to_Z) &
695 inconsistent = inconsistent + 1
697 call sum_across_pes(inconsistent)
699 if ((inconsistent > 0) .and. (is_root_pe()))
then
700 write(mesg,
'("Thickness initial conditions are inconsistent ",'// &
701 '"with topography in ",I8," places.")') inconsistent
702 call mom_error(warning, mesg)
722 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1),
intent(inout) :: eta
723 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(inout) :: h
725 integer :: i, j, k, is, ie, js, je, nz, contractions, dilations
726 real :: hTolerance = 0.1
727 real :: hTmp, eTmp, dilate
728 character(len=100) :: mesg
730 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
731 htolerance = 0.1*us%m_to_Z
734 do j=js,je ;
do i=is,ie
735 if (-eta(i,j,nz+1) > g%bathyT(i,j) + htolerance)
then
736 eta(i,j,nz+1) = -g%bathyT(i,j)
737 contractions = contractions + 1
740 call sum_across_pes(contractions)
741 if ((contractions > 0) .and. (is_root_pe()))
then
742 write(mesg,
'("Thickness initial conditions were contracted ",'// &
743 '"to fit topography in ",I8," places.")') contractions
744 call mom_error(warning,
'adjustEtaToFitBathymetry: '//mesg)
749 do k=nz,1,-1 ;
do j=js,je ;
do i=is,ie
752 if (eta(i,j,k) < (eta(i,j,k+1) + gv%Angstrom_Z))
then
753 eta(i,j,k) = eta(i,j,k+1) + gv%Angstrom_Z
754 h(i,j,k) = gv%Angstrom_Z
756 h(i,j,k) = (eta(i,j,k) - eta(i,j,k+1))
758 enddo ;
enddo ;
enddo
761 do j=js,je ;
do i=is,ie
765 if (-eta(i,j,nz+1) < g%bathyT(i,j) - htolerance)
then
766 dilations = dilations + 1
767 if (eta(i,j,1) <= eta(i,j,nz+1))
then
768 do k=1,nz ; h(i,j,k) = (eta(i,j,1) + g%bathyT(i,j)) / real(nz) ;
enddo
770 dilate = (eta(i,j,1) + g%bathyT(i,j)) / (eta(i,j,1) - eta(i,j,nz+1))
771 do k=1,nz ; h(i,j,k) = h(i,j,k) * dilate ;
enddo
773 do k=nz,2,-1 ; eta(i,j,k) = eta(i,j,k+1) + h(i,j,k) ;
enddo
778 do k=1,nz ;
do j=js,je ;
do i=is,ie
779 h(i,j,k) = h(i,j,k)*gv%Z_to_H
780 enddo ;
enddo ;
enddo
782 call sum_across_pes(dilations)
783 if ((dilations > 0) .and. (is_root_pe()))
then
784 write(mesg,
'("Thickness initial conditions were dilated ",'// &
785 '"to fit topography in ",I8," places.")') dilations
786 call mom_error(warning,
'adjustEtaToFitBathymetry: '//mesg)
795 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
799 logical,
optional,
intent(in) :: just_read_params
802 character(len=40) :: mdl =
"initialize_thickness_uniform"
803 real :: e0(SZK_(G)+1)
805 real :: eta1D(SZK_(G)+1)
808 integer :: i, j, k, is, ie, js, je, nz
810 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
812 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
814 if (just_read)
return
816 call calltree_enter(trim(mdl)//
"(), MOM_state_initialization.F90")
818 if (g%max_depth<=0.)
call mom_error(fatal,
"initialize_thickness_uniform: "// &
819 "MAXIMUM_DEPTH has a non-sensical value! Was it set?")
822 e0(k) = -g%max_depth * real(k-1) / real(nz)
825 do j=js,je ;
do i=is,ie
831 eta1d(nz+1) = -g%bathyT(i,j)
834 if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z))
then
835 eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
836 h(i,j,k) = gv%Angstrom_H
838 h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
851 real,
dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
855 logical,
optional,
intent(in) :: just_read_params
858 character(len=40) :: mdl =
"initialize_thickness_list"
859 real :: e0(SZK_(G)+1)
861 real :: eta1D(SZK_(G)+1)
864 character(len=200) :: filename, eta_file, inputdir
865 character(len=72) :: eta_var
866 integer :: i, j, k, is, ie, js, je, nz
868 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
870 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
872 call get_param(param_file, mdl,
"INTERFACE_IC_FILE", eta_file, &
873 "The file from which horizontal mean initial conditions "//&
874 "for interface depths can be read.", fail_if_missing=.true.)
875 call get_param(param_file, mdl,
"INTERFACE_IC_VAR", eta_var, &
876 "The variable name for horizontal mean initial conditions "//&
877 "for interface depths relative to mean sea level.", &
880 if (just_read)
return
882 call calltree_enter(trim(mdl)//
"(), MOM_state_initialization.F90")
884 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
885 filename = trim(slasher(inputdir))//trim(eta_file)
886 call log_param(param_file, mdl,
"INPUTDIR/INTERFACE_IC_FILE", filename)
889 call mom_read_data(filename, eta_var, e0(:), scale=us%m_to_Z)
891 if ((abs(e0(1)) - 0.0) > 0.001)
then
893 do k=nz+1,2,-1 ; e0(k) = e0(k-1) ;
enddo
897 if (e0(2) > e0(1))
then
898 do k=1,nz ; e0(k) = -e0(k) ;
enddo
901 do j=js,je ;
do i=is,ie
907 eta1d(nz+1) = -g%bathyT(i,j)
910 if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z))
then
911 eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
912 h(i,j,k) = gv%Angstrom_H
914 h(i,j,k) = gv%Z_to_H * (eta1d(k) - eta1d(k+1))
924 call mom_error(fatal,
" MOM_state_initialization.F90, initialize_thickness_search: NOT IMPLEMENTED")
932 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
938 real,
dimension(SZI_(G),SZJ_(G)) :: &
940 real :: dz_geo(SZI_(G),SZJ_(G))
946 logical :: Boussinesq
947 integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
948 integer :: itt, max_itt
950 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
951 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
953 boussinesq = gv%Boussinesq
954 i_gearth = 1.0 / (gv%mks_g_Earth)
955 hm_rho_to_pa = gv%mks_g_Earth * gv%H_to_m
958 call mom_error(fatal,
"Not yet converting thickness with Boussinesq approx.")
960 if (
associated(tv%eqn_of_state))
then
961 do j=jsq,jeq+1 ;
do i=isq,ieq+1
962 p_bot(i,j) = 0.0 ; p_top(i,j) = 0.0
966 do i=is,ie ; p_top(i,j) = p_bot(i,j) ;
enddo
968 is, ie-is+1, tv%eqn_of_state)
970 p_bot(i,j) = p_top(i,j) + hm_rho_to_pa * (h(i,j,k) * rho(i))
976 0.0, g%HI, tv%eqn_of_state, dz_geo)
977 if (itt < max_itt)
then ;
do j=js,je
979 is, ie-is+1, tv%eqn_of_state)
984 p_bot(i,j) = p_bot(i,j) + rho(i) * &
985 (hm_rho_to_pa*h(i,j,k) - dz_geo(i,j))
990 do j=js,je ;
do i=is,ie
991 h(i,j,k) = (p_bot(i,j) - p_top(i,j)) * gv%kg_m2_to_H * i_gearth
995 do k=1,nz ;
do j=js,je ;
do i=is,ie
996 h(i,j,k) = h(i,j,k) * (gv%Rlay(k) / gv%Rho0)
997 enddo ;
enddo ;
enddo
1004 subroutine depress_surface(h, G, GV, US, param_file, tv, just_read_params)
1008 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1012 logical,
optional,
intent(in) :: just_read_params
1015 real,
dimension(SZI_(G),SZJ_(G)) :: &
1017 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &
1020 real :: scale_factor
1022 character(len=40) :: mdl =
"depress_surface"
1023 character(len=200) :: inputdir, eta_srf_file
1024 character(len=200) :: filename, eta_srf_var
1025 logical :: just_read
1026 integer :: i, j, k, is, ie, js, je, nz
1027 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1029 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
1033 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
1034 inputdir = slasher(inputdir)
1035 call get_param(param_file, mdl,
"SURFACE_HEIGHT_IC_FILE", eta_srf_file,&
1036 "The initial condition file for the surface height.", &
1037 fail_if_missing=.not.just_read, do_not_log=just_read)
1038 call get_param(param_file, mdl,
"SURFACE_HEIGHT_IC_VAR", eta_srf_var, &
1039 "The initial condition variable for the surface height.",&
1040 default=
"SSH", do_not_log=just_read)
1041 filename = trim(inputdir)//trim(eta_srf_file)
1042 if (.not.just_read) &
1043 call log_param(param_file, mdl,
"INPUTDIR/SURFACE_HEIGHT_IC_FILE", filename)
1045 call get_param(param_file, mdl,
"SURFACE_HEIGHT_IC_SCALE", scale_factor, &
1046 "A scaling factor to convert SURFACE_HEIGHT_IC_VAR into "//&
1047 "units of m", units=
"variable", default=1.0, do_not_log=just_read)
1049 if (just_read)
return
1051 call mom_read_data(filename, eta_srf_var, eta_sfc, g%Domain, scale=scale_factor)
1054 call find_eta(h, tv, g, gv, us, eta, eta_to_m=1.0)
1056 do j=js,je ;
do i=is,ie ;
if (g%mask2dT(i,j) > 0.0)
then
1060 if (eta_sfc(i,j) > eta(i,j,1))
then
1062 if (eta_sfc(i,j) - eta(i,j,nz+1) > 10.0*(eta(i,j,1) - eta(i,j,nz+1)))
then
1064 call mom_error(warning,
"Free surface height dilation attempted "//&
1065 "to exceed 10-fold.", all_print=.true.)
1067 dilate = (eta_sfc(i,j) - eta(i,j,nz+1)) / (eta(i,j,1) - eta(i,j,nz+1))
1069 do k=1,nz ; h(i,j,k) = h(i,j,k) * dilate ;
enddo
1070 elseif (eta(i,j,1) > eta_sfc(i,j))
then
1073 if (eta(i,j,k) <= eta_sfc(i,j))
exit
1074 if (eta(i,j,k+1) >= eta_sfc(i,j))
then
1075 h(i,j,k) = gv%Angstrom_H
1077 h(i,j,k) = max(gv%Angstrom_H, h(i,j,k) * &
1078 (eta_sfc(i,j) - eta(i,j,k+1)) / (eta(i,j,k) - eta(i,j,k+1)) )
1082 endif ;
enddo ;
enddo
1088 subroutine trim_for_ice(PF, G, GV, US, ALE_CSp, tv, h, just_read_params)
1093 type(
ale_cs),
pointer :: ALE_CSp
1095 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1097 logical,
optional,
intent(in) :: just_read_params
1100 character(len=200) :: mdl =
"trim_for_ice"
1101 real,
dimension(SZI_(G),SZJ_(G)) :: p_surf
1102 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: S_t, S_b
1103 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)) :: T_t, T_b
1104 character(len=200) :: inputdir, filename, p_surf_file, p_surf_var
1105 real :: scale_factor
1106 real :: min_thickness
1108 logical :: just_read
1109 logical :: use_remapping
1112 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
1114 call get_param(pf, mdl,
"SURFACE_PRESSURE_FILE", p_surf_file, &
1115 "The initial condition file for the surface height.", &
1116 fail_if_missing=.not.just_read, do_not_log=just_read)
1117 call get_param(pf, mdl,
"SURFACE_PRESSURE_VAR", p_surf_var, &
1118 "The initial condition variable for the surface height.", &
1119 units=
"kg m-2", default=
"", do_not_log=just_read)
1120 call get_param(pf, mdl,
"INPUTDIR", inputdir, default=
".", do_not_log=.true.)
1121 filename = trim(slasher(inputdir))//trim(p_surf_file)
1122 if (.not.just_read)
call log_param(pf, mdl,
"!INPUTDIR/SURFACE_HEIGHT_IC_FILE", filename)
1124 call get_param(pf, mdl,
"SURFACE_PRESSURE_SCALE", scale_factor, &
1125 "A scaling factor to convert SURFACE_PRESSURE_VAR from "//&
1126 "file SURFACE_PRESSURE_FILE into a surface pressure.", &
1127 units=
"file dependent", default=1., do_not_log=just_read)
1128 call get_param(pf, mdl,
"MIN_THICKNESS", min_thickness,
'Minimum layer thickness', &
1129 units=
'm', default=1.e-3, do_not_log=just_read, scale=us%m_to_Z)
1130 call get_param(pf, mdl,
"TRIMMING_USES_REMAPPING", use_remapping, &
1131 'When trimming the column, also remap T and S.', &
1132 default=.false., do_not_log=just_read)
1134 if (just_read)
return
1136 call mom_read_data(filename, p_surf_var, p_surf, g%Domain, scale=scale_factor)
1138 if (use_remapping)
then
1140 call initialize_remapping(remap_cs,
'PLM', boundary_extrapolation=.true.)
1144 if (
associated(ale_csp) )
then
1145 call pressure_gradient_plm(ale_csp, s_t, s_b, t_t, t_b, g, gv, tv, h, .true.)
1148 do k=1,g%ke ;
do j=g%jsc,g%jec ;
do i=g%isc,g%iec
1149 t_t(i,j,k) = tv%T(i,j,k) ; t_b(i,j,k) = tv%T(i,j,k)
1150 s_t(i,j,k) = tv%S(i,j,k) ; s_b(i,j,k) = tv%S(i,j,k)
1151 enddo ;
enddo ;
enddo
1154 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
1155 call cut_off_column_top(gv%ke, tv, gv, us, gv%mks_g_Earth*us%Z_to_m, g%bathyT(i,j), &
1156 min_thickness, tv%T(i,j,:), t_t(i,j,:), t_b(i,j,:), &
1157 tv%S(i,j,:), s_t(i,j,:), s_b(i,j,:), p_surf(i,j), h(i,j,:), remap_cs, &
1158 z_tol=1.0e-5*us%m_to_Z)
1167 T, T_t, T_b, S, S_t, S_b, p_surf, h, remap_CS, z_tol)
1168 integer,
intent(in) :: nk
1172 real,
intent(in) :: G_earth
1173 real,
intent(in) :: depth
1174 real,
intent(in) :: min_thickness
1175 real,
dimension(nk),
intent(inout) :: T
1176 real,
dimension(nk),
intent(in) :: T_t
1177 real,
dimension(nk),
intent(in) :: T_b
1178 real,
dimension(nk),
intent(inout) :: S
1179 real,
dimension(nk),
intent(in) :: S_t
1180 real,
dimension(nk),
intent(in) :: S_b
1181 real,
intent(in) :: p_surf
1182 real,
dimension(nk),
intent(inout) :: h
1185 real,
optional,
intent(in) :: z_tol
1189 real,
dimension(nk+1) :: e
1190 real,
dimension(nk) :: h0, S0, T0, h1, S1, T1
1191 real :: P_t, P_b, z_out, e_top
1197 e(k) = e(k+1) + gv%H_to_Z*h(k)
1204 call find_depth_of_pressure_in_cell(t_t(k), t_b(k), s_t(k), s_b(k), e(k), e(k+1), &
1205 p_t, p_surf, us%R_to_kg_m3*gv%Rho0, g_earth, tv%eqn_of_state, &
1206 p_b, z_out, z_tol=z_tol)
1207 if (z_out>=e(k))
then
1210 elseif (z_out<=e(k+1))
then
1220 if (e_top<e(1))
then
1223 if (e(k) > e_top)
then
1226 e_top = e_top - min_thickness
1229 h(k) = gv%Z_to_H * max( min_thickness, e(k) - e(k+1) )
1230 if (e(k) < e_top)
exit
1236 if (
associated(remap_cs))
then
1242 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)
1243 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)
1255 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1257 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1262 logical,
optional,
intent(in) :: just_read_params
1265 character(len=40) :: mdl =
"initialize_velocity_from_file"
1266 character(len=200) :: filename,velocity_file,inputdir
1267 logical :: just_read
1269 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
1271 if (.not.just_read)
call calltree_enter(trim(mdl)//
"(), MOM_state_initialization.F90")
1273 call get_param(param_file, mdl,
"VELOCITY_FILE", velocity_file, &
1274 "The name of the velocity initial condition file.", &
1275 fail_if_missing=.not.just_read, do_not_log=just_read)
1276 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
1277 inputdir = slasher(inputdir)
1279 if (just_read)
return
1281 filename = trim(inputdir)//trim(velocity_file)
1282 call log_param(param_file, mdl,
"INPUTDIR/VELOCITY_FILE", filename)
1284 if (.not.
file_exists(filename, g%Domain))
call mom_error(fatal, &
1285 " initialize_velocity_from_file: Unable to open "//trim(filename))
1288 call mom_read_vector(filename,
"u",
"v", u(:,:,:), v(:,:,:), g%Domain, scale=us%m_s_to_L_T)
1296 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1298 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1302 logical,
optional,
intent(in) :: just_read_params
1305 character(len=200) :: mdl =
"initialize_velocity_zero"
1306 logical :: just_read
1307 integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
1308 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1309 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1311 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
1313 if (.not.just_read)
call calltree_enter(trim(mdl)//
"(), MOM_state_initialization.F90")
1315 if (just_read)
return
1317 do k=1,nz ;
do j=js,je ;
do i=isq,ieq
1319 enddo ;
enddo ;
enddo
1320 do k=1,nz ;
do j=jsq,jeq ;
do i=is,ie
1322 enddo ;
enddo ;
enddo
1330 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1332 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1337 logical,
optional,
intent(in) :: just_read_params
1340 integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
1341 real :: initial_u_const, initial_v_const
1342 logical :: just_read
1343 character(len=200) :: mdl =
"initialize_velocity_uniform"
1344 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1345 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1347 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
1349 call get_param(param_file, mdl,
"INITIAL_U_CONST", initial_u_const, &
1350 "A initial uniform value for the zonal flow.", &
1351 units=
"m s-1", scale=us%m_s_to_L_T, fail_if_missing=.not.just_read, do_not_log=just_read)
1352 call get_param(param_file, mdl,
"INITIAL_V_CONST", initial_v_const, &
1353 "A initial uniform value for the meridional flow.", &
1354 units=
"m s-1", scale=us%m_s_to_L_T, fail_if_missing=.not.just_read, do_not_log=just_read)
1356 if (just_read)
return
1358 do k=1,nz ;
do j=js,je ;
do i=isq,ieq
1359 u(i,j,k) = initial_u_const
1360 enddo ;
enddo ;
enddo
1361 do k=1,nz ;
do j=jsq,jeq ;
do i=is,ie
1362 v(i,j,k) = initial_v_const
1363 enddo ;
enddo ;
enddo
1371 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1373 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1378 logical,
optional,
intent(in) :: just_read_params
1381 character(len=200) :: mdl =
"initialize_velocity_circular"
1382 real :: circular_max_u
1383 real :: dpi, psi1, psi2
1384 logical :: just_read
1385 integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
1386 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1387 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1389 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
1391 call get_param(param_file, mdl,
"CIRCULAR_MAX_U", circular_max_u, &
1392 "The amplitude of zonal flow from which to scale the "// &
1393 "circular stream function [m s-1].", &
1394 units=
"m s-1", default=0., scale=us%L_T_to_m_s, do_not_log=just_read)
1396 if (just_read)
return
1400 do k=1,nz ;
do j=js,je ;
do i=isq,ieq
1403 u(i,j,k) = (psi1-psi2) / (g%US%L_to_m*g%dy_Cu(i,j))
1404 enddo ;
enddo ;
enddo
1405 do k=1,nz ;
do j=jsq,jeq ;
do i=is,ie
1408 v(i,j,k) = (psi2-psi1) / (g%US%L_to_m*g%dx_Cv(i,j))
1409 enddo ;
enddo ;
enddo
1414 real function my_psi(ig,jg)
1420 x = 2.0*(g%geoLonBu(ig,jg)-g%west_lon) / g%len_lon - 1.0
1421 y = 2.0*(g%geoLatBu(ig,jg)-g%south_lat) / g%len_lat - 1.0
1422 r = sqrt( x**2 + y**2 )
1424 my_psi = 0.5*(1.0 - cos(dpi*r))
1425 my_psi = my_psi * (circular_max_u*g%len_lon*1e3 / dpi)
1432 type(ocean_grid_type),
intent(in) :: G
1433 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(out) :: T
1435 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(out) :: S
1438 logical,
optional,
intent(in) :: just_read_params
1441 logical :: just_read
1442 character(len=200) :: filename, salt_filename
1443 character(len=200) :: ts_file, salt_file, inputdir
1444 character(len=40) :: mdl =
"initialize_temp_salt_from_file"
1445 character(len=64) :: temp_var, salt_var
1447 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
1449 if (.not.just_read)
call calltree_enter(trim(mdl)//
"(), MOM_state_initialization.F90")
1451 call get_param(param_file, mdl,
"TS_FILE", ts_file, &
1452 "The initial condition file for temperature.", &
1453 fail_if_missing=.not.just_read, do_not_log=just_read)
1454 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
1455 inputdir = slasher(inputdir)
1457 filename = trim(inputdir)//trim(ts_file)
1458 if (.not.just_read)
call log_param(param_file, mdl,
"INPUTDIR/TS_FILE", filename)
1459 call get_param(param_file, mdl,
"TEMP_IC_VAR", temp_var, &
1460 "The initial condition variable for potential temperature.", &
1461 default=
"PTEMP", do_not_log=just_read)
1462 call get_param(param_file, mdl,
"SALT_IC_VAR", salt_var, &
1463 "The initial condition variable for salinity.", &
1464 default=
"SALT", do_not_log=just_read)
1465 call get_param(param_file, mdl,
"SALT_FILE", salt_file, &
1466 "The initial condition file for salinity.", &
1467 default=trim(ts_file), do_not_log=just_read)
1469 if (just_read)
return
1471 if (.not.file_exists(filename, g%Domain))
call mom_error(fatal, &
1472 " initialize_temp_salt_from_file: Unable to open "//trim(filename))
1475 call mom_read_data(filename, temp_var, t(:,:,:), g%Domain)
1477 salt_filename = trim(inputdir)//trim(salt_file)
1478 if (.not.file_exists(salt_filename, g%Domain))
call mom_error(fatal, &
1479 " initialize_temp_salt_from_file: Unable to open "//trim(salt_filename))
1481 call mom_read_data(salt_filename, salt_var, s(:,:,:), g%Domain)
1483 call calltree_leave(trim(mdl)//
'()')
1488 type(ocean_grid_type),
intent(in) :: G
1489 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(out) :: T
1491 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(out) :: S
1494 logical,
optional,
intent(in) :: just_read_params
1497 real,
dimension(SZK_(G)) :: T0, S0
1499 logical :: just_read
1500 character(len=200) :: filename, ts_file, inputdir
1501 character(len=40) :: mdl =
"initialize_temp_salt_from_profile"
1503 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
1505 if (.not.just_read)
call calltree_enter(trim(mdl)//
"(), MOM_state_initialization.F90")
1507 call get_param(param_file, mdl,
"TS_FILE", ts_file, &
1508 "The file with the reference profiles for temperature "//&
1509 "and salinity.", fail_if_missing=.not.just_read, do_not_log=just_read)
1511 if (just_read)
return
1513 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
1514 inputdir = slasher(inputdir)
1515 filename = trim(inputdir)//trim(ts_file)
1516 call log_param(param_file, mdl,
"INPUTDIR/TS_FILE", filename)
1517 if (.not.file_exists(filename))
call mom_error(fatal, &
1518 " initialize_temp_salt_from_profile: Unable to open "//trim(filename))
1521 call mom_read_data(filename,
"PTEMP", t0(:))
1522 call mom_read_data(filename,
"SALT", s0(:))
1524 do k=1,g%ke ;
do j=g%jsc,g%jec ;
do i=g%isc,g%iec
1525 t(i,j,k) = t0(k) ; s(i,j,k) = s0(k)
1526 enddo ;
enddo ;
enddo
1528 call calltree_leave(trim(mdl)//
'()')
1533 type(ocean_grid_type),
intent(in) :: G
1534 type(verticalgrid_type),
intent(in) :: GV
1535 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(out) :: T
1537 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(out) :: S
1539 type(unit_scale_type),
intent(in) :: US
1542 type(eos_type),
pointer :: eqn_of_state
1543 real,
intent(in) :: P_Ref
1545 logical,
optional,
intent(in) :: just_read_params
1552 real :: pres(SZK_(G))
1553 real :: drho_dT(SZK_(G))
1554 real :: drho_dS(SZK_(G))
1555 real :: rho_guess(SZK_(G))
1556 logical :: fit_salin
1557 logical :: just_read
1558 character(len=40) :: mdl =
"initialize_temp_salt_fit"
1559 integer :: i, j, k, itt, nz
1562 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
1564 if (.not.just_read)
call calltree_enter(trim(mdl)//
"(), MOM_state_initialization.F90")
1566 call get_param(param_file, mdl,
"T_REF", t_ref, &
1567 "A reference temperature used in initialization.", &
1568 units=
"degC", fail_if_missing=.not.just_read, do_not_log=just_read)
1569 call get_param(param_file, mdl,
"S_REF", s_ref, &
1570 "A reference salinity used in initialization.", units=
"PSU", &
1571 default=35.0, do_not_log=just_read)
1572 call get_param(param_file, mdl,
"FIT_SALINITY", fit_salin, &
1573 "If true, accept the prescribed temperature and fit the "//&
1574 "salinity; otherwise take salinity and fit temperature.", &
1575 default=.false., do_not_log=just_read)
1577 if (just_read)
return
1580 pres(k) = p_ref ; s0(k) = s_ref
1584 call calculate_density(t0(1),s0(1),pres(1),rho_guess(1),eqn_of_state, scale=us%kg_m3_to_R)
1585 call calculate_density_derivs(t0,s0,pres,drho_dt,drho_ds,1,1,eqn_of_state, scale=us%kg_m3_to_R)
1590 s0(k) = max(0.0, s0(1) + (gv%Rlay(k) - rho_guess(1)) / drho_ds(1))
1594 call calculate_density(t0,s0,pres,rho_guess,1,nz,eqn_of_state, scale=us%kg_m3_to_R)
1595 call calculate_density_derivs(t0,s0,pres,drho_dt,drho_ds,1,nz,eqn_of_state, scale=us%kg_m3_to_R)
1597 s0(k) = max(0.0, s0(k) + (gv%Rlay(k) - rho_guess(k)) / drho_ds(k))
1603 t0(k) = t0(1) + (gv%Rlay(k) - rho_guess(1)) / drho_dt(1)
1606 call calculate_density(t0,s0,pres,rho_guess,1,nz,eqn_of_state, scale=us%kg_m3_to_R)
1607 call calculate_density_derivs(t0,s0,pres,drho_dt,drho_ds,1,nz,eqn_of_state, scale=us%kg_m3_to_R)
1609 t0(k) = t0(k) + (gv%Rlay(k) - rho_guess(k)) / drho_dt(k)
1614 do k=1,nz ;
do j=g%jsd,g%jed ;
do i=g%isd,g%ied
1615 t(i,j,k) = t0(k) ; s(i,j,k) = s0(k)
1616 enddo ;
enddo ;
enddo
1618 call calltree_leave(trim(mdl)//
'()')
1627 type(ocean_grid_type),
intent(in) :: G
1628 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(out) :: T
1630 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(out) :: S
1634 logical,
optional,
intent(in) :: just_read_params
1640 real :: delta_S, delta_T
1641 real :: S_top, T_top
1642 real :: S_range, T_range
1644 logical :: just_read
1645 character(len=40) :: mdl =
"initialize_temp_salt_linear"
1647 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
1649 if (.not.just_read)
call calltree_enter(trim(mdl)//
"(), MOM_state_initialization.F90")
1650 call get_param(param_file, mdl,
"T_TOP", t_top, &
1651 "Initial temperature of the top surface.", &
1652 units=
"degC", fail_if_missing=.not.just_read, do_not_log=just_read)
1653 call get_param(param_file, mdl,
"T_RANGE", t_range, &
1654 "Initial temperature difference (top-bottom).", &
1655 units=
"degC", fail_if_missing=.not.just_read, do_not_log=just_read)
1656 call get_param(param_file, mdl,
"S_TOP", s_top, &
1657 "Initial salinity of the top surface.", &
1658 units=
"PSU", fail_if_missing=.not.just_read, do_not_log=just_read)
1659 call get_param(param_file, mdl,
"S_RANGE", s_range, &
1660 "Initial salinity difference (top-bottom).", &
1661 units=
"PSU", fail_if_missing=.not.just_read, do_not_log=just_read)
1663 if (just_read)
return
1672 s(:,:,k) = s_top - s_range*((real(k)-0.5)/real(g%ke))
1673 t(:,:,k) = t_top - t_range*((real(k)-0.5)/real(g%ke))
1685 call calltree_leave(trim(mdl)//
'()')
1694 type(ocean_grid_type),
intent(in) :: G
1695 type(verticalgrid_type),
intent(in) :: GV
1696 type(unit_scale_type),
intent(in) :: US
1697 logical,
intent(in) :: use_temperature
1698 type(thermo_var_ptrs),
intent(in) :: tv
1701 type(sponge_cs),
pointer :: CSp
1703 type(ale_sponge_cs),
pointer :: ALE_CSp
1705 type(time_type),
intent(in) :: Time
1708 real,
allocatable,
dimension(:,:,:) :: eta
1709 real,
allocatable,
dimension(:,:,:) :: h
1711 real,
dimension (SZI_(G),SZJ_(G),SZK_(G)) :: &
1713 real,
dimension (SZI_(G),SZJ_(G)) :: &
1716 real :: Idamp(SZI_(G),SZJ_(G))
1717 real :: pres(SZI_(G))
1719 integer :: i, j, k, is, ie, js, je, nz
1720 integer :: isd, ied, jsd, jed
1721 integer,
dimension(4) :: siz
1723 character(len=40) :: potemp_var, salin_var, Idamp_var, eta_var
1724 character(len=40) :: mdl =
"initialize_sponges_file"
1725 character(len=200) :: damping_file, state_file
1726 character(len=200) :: filename, inputdir
1729 logical :: new_sponges
1732 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1733 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
1735 pres(:) = 0.0 ; tmp(:,:,:) = 0.0 ; idamp(:,:) = 0.0
1737 call get_param(param_file, mdl,
"INPUTDIR", inputdir, default=
".")
1738 inputdir = slasher(inputdir)
1739 call get_param(param_file, mdl,
"SPONGE_DAMPING_FILE", damping_file, &
1740 "The name of the file with the sponge damping rates.", &
1741 fail_if_missing=.true.)
1742 call get_param(param_file, mdl,
"SPONGE_STATE_FILE", state_file, &
1743 "The name of the file with the state to damp toward.", &
1744 default=damping_file)
1745 call get_param(param_file, mdl,
"SPONGE_PTEMP_VAR", potemp_var, &
1746 "The name of the potential temperature variable in "//&
1747 "SPONGE_STATE_FILE.", default=
"PTEMP")
1748 call get_param(param_file, mdl,
"SPONGE_SALT_VAR", salin_var, &
1749 "The name of the salinity variable in "//&
1750 "SPONGE_STATE_FILE.", default=
"SALT")
1751 call get_param(param_file, mdl,
"SPONGE_ETA_VAR", eta_var, &
1752 "The name of the interface height variable in "//&
1753 "SPONGE_STATE_FILE.", default=
"ETA")
1754 call get_param(param_file, mdl,
"SPONGE_IDAMP_VAR", idamp_var, &
1755 "The name of the inverse damping rate variable in "//&
1756 "SPONGE_DAMPING_FILE.", default=
"IDAMP")
1757 call get_param(param_file, mdl,
"USE_REGRIDDING", use_ale, do_not_log = .true.)
1759 call get_param(param_file, mdl,
"NEW_SPONGES", new_sponges, &
1760 "Set True if using the newer sponging code which "//&
1761 "performs on-the-fly regridding in lat-lon-time.",&
1762 "of sponge restoring data.", default=.false.)
1770 filename = trim(inputdir)//trim(damping_file)
1771 call log_param(param_file, mdl,
"INPUTDIR/SPONGE_DAMPING_FILE", filename)
1772 if (.not.file_exists(filename, g%Domain)) &
1773 call mom_error(fatal,
" initialize_sponges: Unable to open "//trim(filename))
1775 if (new_sponges .and. .not. use_ale) &
1776 call mom_error(fatal,
" initialize_sponges: Newer sponges are currently unavailable in layered mode ")
1778 call mom_read_data(filename,
"Idamp", idamp(:,:), g%Domain)
1784 filename = trim(inputdir)//trim(state_file)
1785 call log_param(param_file, mdl,
"INPUTDIR/SPONGE_STATE_FILE", filename)
1786 if (.not.file_exists(filename, g%Domain)) &
1787 call mom_error(fatal,
" initialize_sponges: Unable to open "//trim(filename))
1791 if (.not. use_ale)
then
1792 allocate(eta(isd:ied,jsd:jed,nz+1)); eta(:,:,:) = 0.0
1793 call mom_read_data(filename, eta_var, eta(:,:,:), g%Domain, scale=us%m_to_Z)
1795 do j=js,je ;
do i=is,ie
1796 eta(i,j,nz+1) = -g%bathyT(i,j)
1798 do k=nz,1,-1 ;
do j=js,je ;
do i=is,ie
1799 if (eta(i,j,k) < (eta(i,j,k+1) + gv%Angstrom_Z)) &
1800 eta(i,j,k) = eta(i,j,k+1) + gv%Angstrom_Z
1801 enddo ;
enddo ;
enddo
1804 call initialize_sponge(idamp, eta, g, param_file, csp, gv)
1806 elseif (.not. new_sponges)
then
1808 call field_size(filename,eta_var,siz,no_domain=.true.)
1809 if (siz(1) /= g%ieg-g%isg+1 .or. siz(2) /= g%jeg-g%jsg+1) &
1810 call mom_error(fatal,
"initialize_sponge_file: Array size mismatch for sponge data.")
1815 allocate(eta(isd:ied,jsd:jed,nz_data+1))
1816 allocate(h(isd:ied,jsd:jed,nz_data))
1818 call mom_read_data(filename, eta_var, eta(:,:,:), g%Domain, scale=us%m_to_Z)
1820 do j=js,je ;
do i=is,ie
1821 eta(i,j,nz+1) = -g%bathyT(i,j)
1824 do k=nz,1,-1 ;
do j=js,je ;
do i=is,ie
1825 if (eta(i,j,k) < (eta(i,j,k+1) + gv%Angstrom_Z)) &
1826 eta(i,j,k) = eta(i,j,k+1) + gv%Angstrom_Z
1827 enddo ;
enddo ;
enddo
1828 do k=1,nz;
do j=js,je ;
do i=is,ie
1829 h(i,j,k) = gv%Z_to_H*(eta(i,j,k)-eta(i,j,k+1))
1830 enddo ;
enddo ;
enddo
1831 call initialize_ale_sponge(idamp, g, param_file, ale_csp, h, nz_data)
1836 call initialize_ale_sponge(idamp, g, param_file, ale_csp)
1843 if ( gv%nkml>0 .and. .not. new_sponges)
then
1847 do i=is-1,ie ; pres(i) = tv%P_Ref ;
enddo
1849 call mom_read_data(filename, potemp_var, tmp(:,:,:), g%Domain)
1850 call mom_read_data(filename, salin_var, tmp2(:,:,:), g%Domain)
1853 call calculate_density(tmp(:,j,1), tmp2(:,j,1), pres, tmp_2d(:,j), &
1854 is, ie-is+1, tv%eqn_of_state, scale=us%kg_m3_to_R)
1857 call set_up_sponge_ml_density(tmp_2d, g, csp)
1861 if ( use_temperature .and. .not. new_sponges)
then
1862 call mom_read_data(filename, potemp_var, tmp(:,:,:), g%Domain)
1863 call set_up_sponge_field(tmp, tv%T, g, nz, csp)
1864 call mom_read_data(filename, salin_var, tmp(:,:,:), g%Domain)
1865 call set_up_sponge_field(tmp, tv%S, g, nz, csp)
1866 elseif (use_temperature)
then
1867 call set_up_ale_sponge_field(filename, potemp_var, time, g, gv, us, tv%T, ale_csp)
1868 call set_up_ale_sponge_field(filename, salin_var, time, g, gv, us, tv%S, ale_csp)
1876 type(ocean_grid_type),
intent(inout) :: G
1880 do i=g%isd,g%ied-1 ;
do j=g%jsd,g%jed
1881 g%Dblock_u(i,j) = g%mask2dCu(i,j) * max(g%bathyT(i,j), g%bathyT(i+1,j))
1882 g%Dopen_u(i,j) = g%Dblock_u(i,j)
1884 do i=g%isd,g%ied ;
do j=g%jsd,g%jed-1
1885 g%Dblock_v(i,j) = g%mask2dCv(i,j) * max(g%bathyT(i,j), g%bathyT(i,j+1))
1886 g%Dopen_v(i,j) = g%Dblock_v(i,j)
1893 type(ocean_grid_type),
intent(inout) :: G
1894 type(unit_scale_type),
intent(in) :: US
1896 real,
dimension(G%isc:G%iec, G%jsc:G%jec) :: tmpForSumming
1900 area_scale = us%L_to_m**2
1901 tmpforsumming(:,:) = 0.
1902 g%areaT_global = 0.0 ; g%IareaT_global = 0.0
1903 do j=g%jsc,g%jec ;
do i=g%isc,g%iec
1904 tmpforsumming(i,j) = area_scale*g%areaT(i,j) * g%mask2dT(i,j)
1906 g%areaT_global = reproducing_sum(tmpforsumming)
1907 g%IareaT_global = 1. / (g%areaT_global)
1913 type(ocean_grid_type),
intent(inout) :: G
1917 do i=g%isd,g%ied-1 ;
do j=g%jsd,g%jed
1918 g%Dblock_u(i,j) = g%mask2dCu(i,j) * min(g%bathyT(i,j), g%bathyT(i+1,j))
1919 g%Dopen_u(i,j) = g%Dblock_u(i,j)
1921 do i=g%isd,g%ied ;
do j=g%jsd,g%jed-1
1922 g%Dblock_v(i,j) = g%mask2dCv(i,j) * min(g%bathyT(i,j), g%bathyT(i,j+1))
1923 g%Dopen_v(i,j) = g%Dblock_v(i,j)
1931 type(ocean_grid_type),
intent(inout) :: G
1932 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1934 type(thermo_var_ptrs),
intent(inout) :: tv
1936 type(verticalgrid_type),
intent(in) :: GV
1937 type(unit_scale_type),
intent(in) :: US
1940 logical,
optional,
intent(in) :: just_read_params
1944 character(len=200) :: filename
1946 character(len=200) :: tfilename
1948 character(len=200) :: sfilename
1950 character(len=200) :: shelf_file
1951 character(len=200) :: inputdir
1952 character(len=200) :: mesg, area_varname, ice_shelf_file
1954 type(eos_type),
pointer :: eos => null()
1955 type(thermo_var_ptrs) :: tv_loc
1956 type(verticalgrid_type) :: GV_loc
1958 # include "version_variable.h"
1959 character(len=40) :: mdl =
"MOM_initialize_layers_from_Z"
1961 integer :: is, ie, js, je, nz
1962 integer :: isc,iec,jsc,jec
1963 integer :: isg, ieg, jsg, jeg
1964 integer :: isd, ied, jsd, jed
1966 integer :: i, j, k, ks, np, ni, nj
1967 integer :: idbg, jdbg
1968 integer :: nkml, nkbl
1970 integer :: kd, inconsistent
1976 real,
dimension(:,:),
pointer :: shelf_area => null()
1979 real :: missing_value_temp, missing_value_salt
1980 logical :: correct_thickness
1981 character(len=40) :: potemp_var, salin_var
1982 character(len=8) :: laynum
1984 integer,
parameter :: niter=10
1985 logical :: just_read
1986 logical :: adjust_temperature = .true.
1987 real,
parameter :: missing_value = -1.e20
1988 real,
parameter :: temp_land_fill = 0.0, salt_land_fill = 35.0
1989 logical :: reentrant_x, tripolar_n,dbg
1990 logical :: debug = .false.
1993 real,
dimension(:),
allocatable :: z_edges_in, z_in
1994 real,
dimension(:),
allocatable :: Rb
1995 real,
dimension(:,:,:),
allocatable,
target :: temp_z, salt_z, mask_z
1996 real,
dimension(:,:,:),
allocatable :: rho_z
1997 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: zi
1998 real,
dimension(SZI_(G),SZJ_(G)) :: nlevs
1999 real,
dimension(SZI_(G)) :: press
2002 real,
dimension(:),
allocatable :: hTarget
2003 real,
dimension(:,:),
allocatable :: area_shelf_h
2004 real,
dimension(:,:),
allocatable,
target :: frac_shelf_h
2005 real,
dimension(:,:,:),
allocatable,
target :: tmpT1dIn, tmpS1dIn
2006 real,
dimension(:,:,:),
allocatable :: tmp_mask_in
2007 real,
dimension(:,:,:),
allocatable :: h1
2008 real,
dimension(:,:,:),
allocatable :: dz_interface
2009 real :: zTopOfCell, zBottomOfCell
2010 type(regridding_cs) :: regridCS
2011 type(remapping_cs) :: remapCS
2013 logical :: homogenize, useALEremapping, remap_full_column, remap_general, remap_old_alg
2014 logical :: use_ice_shelf
2015 character(len=10) :: remappingScheme
2016 real :: tempAvg, saltAvg
2017 integer :: nPoints, ans
2018 integer :: id_clock_routine, id_clock_read, id_clock_interp, id_clock_fill, id_clock_ALE
2020 id_clock_routine = cpu_clock_id(
'(Initialize from Z)', grain=clock_routine)
2021 id_clock_ale = cpu_clock_id(
'(Initialize from Z) ALE', grain=clock_loop)
2023 call cpu_clock_begin(id_clock_routine)
2025 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
2026 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
2027 isg = g%isg ; ieg = g%ieg ; jsg = g%jsg ; jeg = g%jeg
2029 pi_180=atan(1.0)/45.
2031 just_read = .false. ;
if (
present(just_read_params)) just_read = just_read_params
2033 if (.not.just_read)
call calltree_enter(trim(mdl)//
"(), MOM_state_initialization.F90")
2034 if (.not.just_read)
call log_version(pf, mdl, version,
"")
2036 inputdir =
"." ;
call get_param(pf, mdl,
"INPUTDIR", inputdir)
2037 inputdir = slasher(inputdir)
2039 eos => tv%eqn_of_state
2043 reentrant_x = .false. ;
call get_param(pf, mdl,
"REENTRANT_X", reentrant_x,default=.true.)
2044 tripolar_n = .false. ;
call get_param(pf, mdl,
"TRIPOLAR_N", tripolar_n, default=.false.)
2045 call get_param(pf, mdl,
"MINIMUM_DEPTH", min_depth, default=0.0, scale=us%m_to_Z)
2047 call get_param(pf, mdl,
"NKML",nkml,default=0)
2048 call get_param(pf, mdl,
"NKBL",nkbl,default=0)
2050 call get_param(pf, mdl,
"TEMP_SALT_Z_INIT_FILE",filename, &
2051 "The name of the z-space input file used to initialize "//&
2052 "temperatures (T) and salinities (S). If T and S are not "//&
2053 "in the same file, TEMP_Z_INIT_FILE and SALT_Z_INIT_FILE "//&
2054 "must be set.",default=
"temp_salt_z.nc",do_not_log=just_read)
2055 call get_param(pf, mdl,
"TEMP_Z_INIT_FILE",tfilename, &
2056 "The name of the z-space input file used to initialize "//&
2057 "temperatures, only.", default=trim(filename),do_not_log=just_read)
2058 call get_param(pf, mdl,
"SALT_Z_INIT_FILE",sfilename, &
2059 "The name of the z-space input file used to initialize "//&
2060 "temperatures, only.", default=trim(filename),do_not_log=just_read)
2061 filename = trim(inputdir)//trim(filename)
2062 tfilename = trim(inputdir)//trim(tfilename)
2063 sfilename = trim(inputdir)//trim(sfilename)
2064 call get_param(pf, mdl,
"Z_INIT_FILE_PTEMP_VAR", potemp_var, &
2065 "The name of the potential temperature variable in "//&
2066 "TEMP_Z_INIT_FILE.", default=
"ptemp",do_not_log=just_read)
2067 call get_param(pf, mdl,
"Z_INIT_FILE_SALT_VAR", salin_var, &
2068 "The name of the salinity variable in "//&
2069 "SALT_Z_INIT_FILE.", default=
"salt",do_not_log=just_read)
2070 call get_param(pf, mdl,
"Z_INIT_HOMOGENIZE", homogenize, &
2071 "If True, then horizontally homogenize the interpolated "//&
2072 "initial conditions.", default=.false., do_not_log=just_read)
2073 call get_param(pf, mdl,
"Z_INIT_ALE_REMAPPING", usealeremapping, &
2074 "If True, then remap straight to model coordinate from file.", &
2075 default=.false., do_not_log=just_read)
2076 call get_param(pf, mdl,
"Z_INIT_REMAPPING_SCHEME", remappingscheme, &
2077 "The remapping scheme to use if using Z_INIT_ALE_REMAPPING "//&
2078 "is True.", default=
"PPM_IH4", do_not_log=just_read)
2079 call get_param(pf, mdl,
"Z_INIT_REMAP_GENERAL", remap_general, &
2080 "If false, only initializes to z* coordinates. "//&
2081 "If true, allows initialization directly to general coordinates.",&
2082 default=.false., do_not_log=just_read)
2083 call get_param(pf, mdl,
"Z_INIT_REMAP_FULL_COLUMN", remap_full_column, &
2084 "If false, only reconstructs profiles for valid data points. "//&
2085 "If true, inserts vanished layers below the valid data.", &
2086 default=remap_general, do_not_log=just_read)
2087 call get_param(pf, mdl,
"Z_INIT_REMAP_OLD_ALG", remap_old_alg, &
2088 "If false, uses the preferred remapping algorithm for initialization. "//&
2089 "If true, use an older, less robust algorithm for remapping.", &
2090 default=.true., do_not_log=just_read)
2091 call get_param(pf, mdl,
"ICE_SHELF", use_ice_shelf, default=.false.)
2092 if (use_ice_shelf)
then
2093 call get_param(pf, mdl,
"ICE_THICKNESS_FILE", ice_shelf_file, &
2094 "The file from which the ice bathymetry and area are read.", &
2095 fail_if_missing=.not.just_read, do_not_log=just_read)
2096 shelf_file = trim(inputdir)//trim(ice_shelf_file)
2097 if (.not.just_read)
call log_param(pf, mdl,
"INPUTDIR/THICKNESS_FILE", shelf_file)
2098 call get_param(pf, mdl,
"ICE_AREA_VARNAME", area_varname, &
2099 "The name of the area variable in ICE_THICKNESS_FILE.", &
2100 fail_if_missing=.not.just_read, do_not_log=just_read)
2102 if (.not.usealeremapping)
then
2103 call get_param(pf, mdl,
"ADJUST_THICKNESS", correct_thickness, &
2104 "If true, all mass below the bottom removed if the "//&
2105 "topography is shallower than the thickness input file "//&
2106 "would indicate.", default=.false., do_not_log=just_read)
2108 call get_param(pf, mdl,
"FIT_TO_TARGET_DENSITY_IC", adjust_temperature, &
2109 "If true, all the interior layers are adjusted to "//&
2110 "their target densities using mostly temperature "//&
2111 "This approach can be problematic, particularly in the "//&
2112 "high latitudes.", default=.true., do_not_log=just_read)
2115 call cpu_clock_end(id_clock_routine)
2120 eps_z = 1.0e-10*us%m_to_Z
2121 eps_rho = 1.0e-10*us%kg_m3_to_R
2138 call horiz_interp_and_extrap_tracer(tfilename, potemp_var, 1.0, 1, &
2139 g, temp_z, mask_z, z_in, z_edges_in, missing_value_temp, reentrant_x, &
2140 tripolar_n, homogenize, m_to_z=us%m_to_Z)
2142 call horiz_interp_and_extrap_tracer(sfilename, salin_var, 1.0, 1, &
2143 g, salt_z, mask_z, z_in, z_edges_in, missing_value_salt, reentrant_x, &
2144 tripolar_n, homogenize, m_to_z=us%m_to_Z)
2149 do k=1,
size(z_edges_in,1) ; z_edges_in(k) = -z_edges_in(k) ;
enddo
2151 allocate(rho_z(isd:ied,jsd:jed,kd))
2152 allocate(area_shelf_h(isd:ied,jsd:jed))
2153 allocate(frac_shelf_h(isd:ied,jsd:jed))
2158 call convert_temp_salt_for_teos10(temp_z, salt_z, press, g, kd, mask_z, eos)
2160 do k=1,kd ;
do j=js,je
2161 call calculate_density(temp_z(:,j,k), salt_z(:,j,k), press, rho_z(:,j,k), is, ie, &
2162 eos, scale=us%kg_m3_to_R)
2171 if (use_ice_shelf)
then
2172 if (.not.file_exists(shelf_file, g%Domain))
call mom_error(fatal, &
2173 "MOM_temp_salt_initialize_from_Z: Unable to open shelf file "//trim(shelf_file))
2175 call mom_read_data(shelf_file, trim(area_varname), area_shelf_h, g%Domain)
2178 frac_shelf_h(:,:) = 0.0
2180 do j=jsd,jed ;
do i=isd,ied
2181 if (g%areaT(i,j) > 0.0) &
2182 frac_shelf_h(i,j) = area_shelf_h(i,j) / (us%L_to_m**2*g%areaT(i,j))
2185 shelf_area => frac_shelf_h
2191 if (usealeremapping)
then
2192 call cpu_clock_begin(id_clock_ale)
2193 nkd = max(gv%ke, kd)
2200 allocate( tmp_mask_in(isd:ied,jsd:jed,nkd) ) ; tmp_mask_in(:,:,:) = 0.
2201 allocate( h1(isd:ied,jsd:jed,nkd) ) ; h1(:,:,:) = 0.
2202 allocate( tmpt1din(isd:ied,jsd:jed,nkd) ) ; tmpt1din(:,:,:) = 0.
2203 allocate( tmps1din(isd:ied,jsd:jed,nkd) ) ; tmps1din(:,:,:) = 0.
2204 do j = js, je ;
do i = is, ie
2205 if (g%mask2dT(i,j)>0.)
then
2206 ztopofcell = 0. ; zbottomofcell = 0.
2207 tmp_mask_in(i,j,1:kd) = mask_z(i,j,:)
2209 if (tmp_mask_in(i,j,k)>0. .and. k<=kd)
then
2210 zbottomofcell = max( z_edges_in(k+1), -g%bathyT(i,j) )
2211 tmpt1din(i,j,k) = temp_z(i,j,k)
2212 tmps1din(i,j,k) = salt_z(i,j,k)
2214 zbottomofcell = -g%bathyT(i,j)
2215 tmpt1din(i,j,k) = tmpt1din(i,j,k-1)
2216 tmps1din(i,j,k) = tmps1din(i,j,k-1)
2218 tmpt1din(i,j,k) = -99.9
2219 tmps1din(i,j,k) = -99.9
2221 h1(i,j,k) = gv%Z_to_H * (ztopofcell - zbottomofcell)
2222 ztopofcell = zbottomofcell
2224 h1(i,j,kd) = h1(i,j,kd) + gv%Z_to_H * max(0., ztopofcell + g%bathyT(i,j) )
2228 deallocate( tmp_mask_in )
2235 call ale_initregridding( gv, us, g%max_depth, pf, mdl, regridcs )
2237 if (.not. remap_general)
then
2239 allocate( htarget(nz) )
2240 htarget = getcoordinateresolution( regridcs )
2241 do j = js, je ;
do i = is, ie
2243 if (g%mask2dT(i,j)>0.)
then
2245 ztopofcell = 0. ; zbottomofcell = 0.
2247 zbottomofcell = max( ztopofcell - htarget(k), -g%bathyT(i,j) )
2248 h(i,j,k) = gv%Z_to_H * (ztopofcell - zbottomofcell)
2249 ztopofcell = zbottomofcell
2256 deallocate( htarget )
2260 call initialize_remapping( remapcs, remappingscheme, boundary_extrapolation=.false. )
2261 if (remap_general)
then
2262 call set_regrid_params( regridcs, min_thickness=0. )
2264 tv_loc%T => tmpt1din
2265 tv_loc%S => tmps1din
2268 allocate( dz_interface(isd:ied,jsd:jed,nkd+1) )
2269 if (use_ice_shelf)
then
2270 call regridding_main( remapcs, regridcs, g, gv_loc, h1, tv_loc, h, dz_interface, shelf_area )
2272 call regridding_main( remapcs, regridcs, g, gv_loc, h1, tv_loc, h, dz_interface )
2274 deallocate( dz_interface )
2276 call ale_remap_scalar(remapcs, g, gv, nkd, h1, tmpt1din, h, tv%T, all_cells=remap_full_column, &
2277 old_remap=remap_old_alg )
2278 call ale_remap_scalar(remapcs, g, gv, nkd, h1, tmps1din, h, tv%S, all_cells=remap_full_column, &
2279 old_remap=remap_old_alg )
2281 deallocate( tmpt1din )
2282 deallocate( tmps1din )
2284 call cpu_clock_end(id_clock_ale)
2290 nlevs = sum(mask_z,dim=3)
2294 do k=2,nz ; rb(k) = 0.5*(gv%Rlay(k-1)+gv%Rlay(k)) ;
enddo
2295 rb(1) = 0.0 ; rb(nz+1) = 2.0*gv%Rlay(nz) - gv%Rlay(nz-1)
2297 zi(is:ie,js:je,:) = find_interfaces(rho_z(is:ie,js:je,:), z_in, rb, g%bathyT(is:ie,js:je), &
2298 nlevs(is:ie,js:je), nkml, nkbl, min_depth, eps_z=eps_z, &
2301 if (correct_thickness)
then
2304 do k=nz,1,-1 ;
do j=js,je ;
do i=is,ie
2305 if (zi(i,j,k) < (zi(i,j,k+1) + gv%Angstrom_Z))
then
2306 zi(i,j,k) = zi(i,j,k+1) + gv%Angstrom_Z
2307 h(i,j,k) = gv%Angstrom_H
2309 h(i,j,k) = gv%Z_to_H * (zi(i,j,k) - zi(i,j,k+1))
2311 enddo ;
enddo ;
enddo
2313 do j=js,je ;
do i=is,ie
2314 if (abs(zi(i,j,nz+1) + g%bathyT(i,j)) > 1.0*us%m_to_Z) &
2315 inconsistent = inconsistent + 1
2317 call sum_across_pes(inconsistent)
2319 if ((inconsistent > 0) .and. (is_root_pe()))
then
2320 write(mesg,
'("Thickness initial conditions are inconsistent ",'// &
2321 '"with topography in ",I5," places.")') inconsistent
2326 tv%T(is:ie,js:je,:) = tracer_z_init(temp_z(is:ie,js:je,:), z_edges_in, zi(is:ie,js:je,:), &
2327 nkml, nkbl, missing_value, g%mask2dT(is:ie,js:je), nz, &
2328 nlevs(is:ie,js:je),dbg,idbg,jdbg, eps_z=eps_z)
2329 tv%S(is:ie,js:je,:) = tracer_z_init(salt_z(is:ie,js:je,:), z_edges_in, zi(is:ie,js:je,:), &
2330 nkml, nkbl, missing_value, g%mask2dT(is:ie,js:je), nz, &
2331 nlevs(is:ie,js:je), eps_z=eps_z)
2334 npoints = 0 ; tempavg = 0. ; saltavg = 0.
2335 do j=js,je ;
do i=is,ie ;
if (g%mask2dT(i,j) >= 1.0)
then
2336 npoints = npoints + 1
2337 tempavg = tempavg + tv%T(i,j,k)
2338 saltavg = saltavg + tv%S(i,j,k)
2339 endif ;
enddo ;
enddo
2342 if (homogenize)
then
2343 call sum_across_pes(npoints)
2344 call sum_across_pes(tempavg)
2345 call sum_across_pes(saltavg)
2347 tempavg = tempavg / real(npoints)
2348 saltavg = saltavg / real(npoints)
2350 tv%T(:,:,k) = tempavg
2351 tv%S(:,:,k) = saltavg
2358 do k=1,nz ;
do j=js,je ;
do i=is,ie
2359 if (tv%T(i,j,k) == missing_value)
then
2360 tv%T(i,j,k) = temp_land_fill
2361 tv%S(i,j,k) = salt_land_fill
2363 enddo ;
enddo ;
enddo
2366 ks = max(0,nkml)+max(0,nkbl)+1
2368 if (adjust_temperature .and. .not. usealeremapping)
then
2369 call determine_temperature(tv%T(is:ie,js:je,:), tv%S(is:ie,js:je,:), &
2370 us%R_to_kg_m3*gv%Rlay(1:nz), tv%p_ref, niter, missing_value, h(is:ie,js:je,:), ks, eos)
2374 deallocate(z_in, z_edges_in, temp_z, salt_z, mask_z)
2375 deallocate(rho_z, area_shelf_h, frac_shelf_h)
2381 call calltree_leave(trim(mdl)//
'()')
2382 call cpu_clock_end(id_clock_routine)
2388 type(ocean_grid_type),
intent(inout) :: G
2389 type(verticalgrid_type),
intent(in) :: GV
2390 type(unit_scale_type),
intent(in) :: US
2391 type(thermo_var_ptrs),
intent(in) :: tv
2393 integer,
parameter :: nk=5
2394 real,
dimension(nk) :: T, T_t, T_b, S, S_t, S_b, rho, h, z
2395 real,
dimension(nk+1) :: e
2397 real :: P_tot, P_t, P_b, z_out
2398 type(remapping_cs),
pointer :: remap_CS => null()
2405 e(k+1) = e(k) - h(k)
2409 z(k) = 0.5 * ( e(k) + e(k+1) )
2410 t_t(k) = 20.+(0./500.)*e(k)
2411 t(k) = 20.+(0./500.)*z(k)
2412 t_b(k) = 20.+(0./500.)*e(k+1)
2413 s_t(k) = 35.-(0./500.)*e(k)
2414 s(k) = 35.+(0./500.)*z(k)
2415 s_b(k) = 35.-(0./500.)*e(k+1)
2416 call calculate_density(0.5*(t_t(k)+t_b(k)), 0.5*(s_t(k)+s_b(k)), -us%R_to_kg_m3*gv%Rho0*gv%mks_g_Earth*z(k), &
2417 rho(k), tv%eqn_of_state)
2418 p_tot = p_tot + gv%mks_g_Earth * rho(k) * h(k)
2423 call find_depth_of_pressure_in_cell(t_t(k), t_b(k), s_t(k), s_b(k), e(k), e(k+1), p_t, 0.5*p_tot, &
2424 us%R_to_kg_m3*gv%Rho0, gv%mks_g_Earth, tv%eqn_of_state, p_b, z_out)
2425 write(0,*) k,p_t,p_b,0.5*p_tot,e(k),e(k+1),z_out
2428 write(0,*) p_b,p_tot
2431 write(0,*)
' ==================================================================== '
2435 t, t_t, t_b, s, s_t, s_b, 0.5*p_tot, h, remap_cs)