MOM_state_initialization.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Initialization functions for state variables, u, v, h, T and S.
7
8use mom_debugging, only : hchksum, qchksum, uvchksum
11use mom_coms, only : max_across_pes, min_across_pes, reproducing_sum
12use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
13use mom_cpu_clock, only : clock_routine, clock_loop
14use mom_domains, only : pass_var, pass_vector, sum_across_pes, broadcast
15use mom_domains, only : root_pe, to_all, scalar_pair, cgrid_ne, agrid
16use mom_error_handler, only : mom_mesg, mom_error, fatal, warning, is_root_pe
18use mom_file_parser, only : get_param, read_param, log_param, param_file_type
19use mom_file_parser, only : log_version
20use mom_get_input, only : directories
22use mom_interface_heights, only : find_eta, dz_to_thickness, dz_to_thickness_simple
24use mom_io, only : file_exists, field_size, mom_read_data, mom_read_vector, slasher
29use mom_restart, only : restore_state, is_new_run, copy_restart_var, copy_restart_vector
33use mom_ale_sponge, only : set_up_ale_sponge_field, set_up_ale_sponge_vel_field
34use mom_ale_sponge, only : ale_sponge_cs, initialize_ale_sponge
36use mom_time_manager, only : time_type, operator(/=)
37use mom_tracer_registry, only : tracer_registry_type
41use mom_eos, only : calculate_density, calculate_density_derivs, eos_type, eos_domain
94use mom_horizontal_regridding, only : horiz_interp_and_extrap_tracer, homogenize_field
98
99implicit none ; private
100
101#include <MOM_memory.h>
102
104
105! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
106! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
107! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
108! vary with the Boussinesq approximation, the Boussinesq variant is given first.
109
110character(len=40) :: mdl = "MOM_state_initialization" !< This module's name.
111
112contains
113
114!> Initialize temporally evolving fields, either as initial
115!! conditions or by reading them from a restart (or saves) file.
116subroutine mom_initialize_state(u, v, h, tv, Time, G, GV, US, PF, dirs, &
117 restart_CS, ALE_CSp, tracer_Reg, sponge_CSp, &
118 ALE_sponge_CSp, oda_incupd_CSp, OBC_for_remap, &
119 Time_in, frac_shelf_h, mass_shelf, OBC_for_bug)
120 type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
121 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
122 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
123 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
124 intent(out) :: u !< The zonal velocity that is being
125 !! initialized [L T-1 ~> m s-1]
126 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
127 intent(out) :: v !< The meridional velocity that is being
128 !! initialized [L T-1 ~> m s-1]
129 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
130 intent(out) :: h !< Layer thicknesses [H ~> m or kg m-2]
131 type(thermo_var_ptrs), intent(inout) :: tv !< A structure pointing to various thermodynamic
132 !! variables
133 type(time_type), intent(inout) :: time !< Time at the start of the run segment.
134 type(param_file_type), intent(in) :: pf !< A structure indicating the open file to parse
135 !! for model parameter values.
136 type(directories), intent(in) :: dirs !< A structure containing several relevant
137 !! directory paths.
138 type(mom_restart_cs), intent(inout) :: restart_cs !< MOM restart control structure
139 type(ale_cs), pointer :: ale_csp !< The ALE control structure for remapping
140 type(tracer_registry_type), pointer :: tracer_reg !< A pointer to the tracer registry
141 type(sponge_cs), pointer :: sponge_csp !< The layerwise sponge control structure.
142 type(ale_sponge_cs), pointer :: ale_sponge_csp !< The ALE sponge control structure.
143 type(ocean_obc_type), pointer :: obc_for_remap !< The open boundary condition control
144 !! structure that may be used for remapping velocities.
145 !! This must be on the unrotated grid, but only the
146 !! position and directions of the OBC faces are used.
147 type(oda_incupd_cs), pointer :: oda_incupd_csp !< The oda_incupd control structure.
148 type(time_type), optional, intent(in) :: time_in !< Time at the start of the run segment.
149 real, dimension(SZI_(G),SZJ_(G)), &
150 optional, intent(in) :: frac_shelf_h !< The fraction of the grid cell covered
151 !! by a floating ice shelf [nondim].
152 real, dimension(SZI_(G),SZJ_(G)), &
153 optional, intent(in) :: mass_shelf !< The mass per unit area of the overlying
154 !! ice shelf [R Z ~> kg m-2]
155 type(ocean_obc_type), optional, pointer :: obc_for_bug !< An open boundary condition control structure
156 !! that might be used to store OBC temperatures and
157 !! salinities if OBC_RESERVOIR_INIT_BUG is true.
158 ! Local variables
159 real :: depth_tot(szi_(g),szj_(g)) ! The nominal total depth of the ocean [Z ~> m]
160 real :: dz(szi_(g),szj_(g),szk_(gv)) ! The layer thicknesses in geopotential (z) units [Z ~> m]
161 character(len=200) :: inputdir ! The directory where NetCDF input files are.
162 character(len=200) :: config, h_config
163 real :: dt ! The baroclinic dynamics timestep for this run [T ~> s].
164
165 logical :: from_z_file, useale
166 logical :: new_sim, rotate_index
167 logical :: use_temperature, use_sponge, use_oda_incupd
168 logical :: verify_restart_time
169 logical :: obc_ts_reservoir_init_bug ! If true, set the OBC temperature and salinity reservoirs
170 ! at the startup of a new run from initial values that are set before remapping.
171 logical :: use_eos ! If true, density is calculated from T & S using an equation of state.
172 logical :: depress_sfc ! If true, remove the mass that would be displaced
173 ! by a large surface pressure by squeezing the column.
174 logical :: trim_ic_for_p_surf ! If true, remove the mass that would be displaced
175 ! by a large surface pressure, such as with an ice sheet.
176 logical :: regrid_accelerate
177 integer :: regrid_iterations
178 logical :: convert
179 logical :: just_read ! If true, only read the parameters because this
180 ! is a run from a restart file; this option
181 ! allows the use of Fatal unused parameters.
182 type(eos_type), pointer :: eos => null()
183 logical :: enable_bugs ! If true, the defaults for recently added bug-fix flags are set to
184 ! recreate the bugs, or if false bugs are only used if actively selected.
185 logical :: debug ! If true, write debugging output.
186 logical :: debug_layers = .false.
187 logical :: use_ice_shelf
188 character(len=80) :: mesg
189 ! This include declares and sets the variable "version".
190# include "version_variable.h"
191 integer :: i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz
192 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
193
194 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
195 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
196 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
197 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
198
199 call calltree_enter("MOM_initialize_state(), MOM_state_initialization.F90")
200 call log_version(pf, mdl, version, "")
201 call get_param(pf, mdl, "DEBUG", debug, default=.false.)
202
203 new_sim = is_new_run(restart_cs)
204 just_read = .not.new_sim
205
206 call get_param(pf, mdl, "INPUTDIR", inputdir, &
207 "The directory in which input files are found.", default=".")
208 inputdir = slasher(inputdir)
209
210 use_temperature = associated(tv%T)
211 useale = associated(ale_csp)
212 use_eos = associated(tv%eqn_of_state)
213 if (use_eos) eos => tv%eqn_of_state
214 use_ice_shelf = PRESENT(frac_shelf_h)
215
216 !====================================================================
217 ! Initialize temporally evolving fields, either as initial
218 ! conditions or by reading them from a restart (or saves) file.
219 !====================================================================
220
221 if (new_sim) then
222 call mom_mesg("Run initialized internally.", 3)
223
224 if (present(time_in)) time = time_in
225 ! Otherwise leave Time at its input value.
226
227 ! This initialization should not be needed. Certainly restricting it
228 ! to the computational domain helps detect possible uninitialized
229 ! data in halos which should be covered by the pass_var(h) later.
230 !do k=1,nz ; do j=js,je ; do i=is,ie
231 ! h(i,j,k) = 0.
232 !enddo
233
234 ! Initialize the layer thicknesses.
235 dz(:,:,:) = 0.0
236 endif
237
238 ! Set the nominal depth of the ocean, which might be different from the bathymetric
239 ! geopotential height, for use by the various initialization routines. G%bathyT has
240 ! already been initialized in previous calls.
241 do j=jsd,jed ; do i=isd,ied
242 depth_tot(i,j) = g%bathyT(i,j) + g%Z_ref
243 enddo ; enddo
244
245 call get_param(pf, mdl, "FATAL_INCONSISTENT_RESTART_TIME", verify_restart_time, &
246 "If true and a time_in value is provided to MOM_initialize_state, verify that "//&
247 "the time read from a restart file is the same as time_in, and issue a fatal "//&
248 "error if it is not. Otherwise, simply set the time to time_in if present.", &
249 default=.false.)
250
251 ! The remaining initialization calls are done, regardless of whether the
252 ! fields are actually initialized here (if just_read=.false.) or whether it
253 ! is just to make sure that all valid parameters are read to enable the
254 ! detection of unused parameters.
255 call get_param(pf, mdl, "INIT_LAYERS_FROM_Z_FILE", from_z_file, &
256 "If true, initialize the layer thicknesses, temperatures, and "//&
257 "salinities from a Z-space file on a latitude-longitude grid.", &
258 default=.false., do_not_log=just_read)
259
260 convert = new_sim ! Thicknesses are initialized in height units in most cases.
261 if (from_z_file) then
262 ! Initialize thickness and T/S from z-coordinate data in a file.
263 if (.NOT.use_temperature) call mom_error(fatal,"MOM_initialize_state : "//&
264 "use_temperature must be true if INIT_LAYERS_FROM_Z_FILE is true")
265
266 call mom_temp_salt_initialize_from_z(h, tv, depth_tot, g, gv, us, pf, &
267 just_read=just_read, frac_shelf_h=frac_shelf_h)
268 convert = .false.
269 else
270 ! Initialize thickness, h.
271 call get_param(pf, mdl, "THICKNESS_CONFIG", h_config, &
272 "A string that determines how the initial layer "//&
273 "thicknesses are specified for a new run: \n"//&
274 " \t file - read interface heights from the file specified \n"//&
275 " \t\t by (THICKNESS_FILE).\n"//&
276 " \t thickness_file - read thicknesses from the file specified \n"//&
277 " \t\t by (THICKNESS_FILE).\n"//&
278 " \t mass_file - read thicknesses in units of mass per unit area from the file \n"//&
279 " \t\t specified by (THICKNESS_FILE).\n"//&
280 " \t coord - determined by ALE coordinate.\n"//&
281 " \t uniform - uniform thickness layers evenly distributed \n"//&
282 " \t\t between the surface and MAXIMUM_DEPTH. \n"//&
283 " \t list - read a list of positive interface depths. \n"//&
284 " \t param - use thicknesses from parameter THICKNESS_INIT_VALUES. \n"//&
285 " \t DOME - use a slope and channel configuration for the \n"//&
286 " \t\t DOME sill-overflow test case. \n"//&
287 " \t ISOMIP - use a configuration for the \n"//&
288 " \t\t ISOMIP test case. \n"//&
289 " \t benchmark - use the benchmark test case thicknesses. \n"//&
290 " \t Neverworld - use the Neverworld test case thicknesses. \n"//&
291 " \t search - search a density profile for the interface \n"//&
292 " \t\t densities. This is not yet implemented. \n"//&
293 " \t circle_obcs - the circle_obcs test case is used. \n"//&
294 " \t DOME2D - 2D version of DOME initialization. \n"//&
295 " \t adjustment2d - 2D lock exchange thickness ICs. \n"//&
296 " \t sloshing - sloshing gravity thickness ICs. \n"//&
297 " \t seamount - no motion test with seamount ICs. \n"//&
298 " \t dumbbell - sloshing channel ICs. \n"//&
299 " \t soliton - Equatorial Rossby soliton. \n"//&
300 " \t rossby_front - a mixed layer front in thermal wind balance.\n"//&
301 " \t USER - call a user modified routine.", &
302 default="uniform", do_not_log=just_read)
303 select case (trim(h_config))
304 case ("file")
305 call initialize_thickness_from_file(dz, depth_tot, g, gv, us, pf, file_has_thickness=.false., &
306 mass_file=.false., just_read=just_read)
307 case ("thickness_file")
308 call initialize_thickness_from_file(dz, depth_tot, g, gv, us, pf, file_has_thickness=.true., &
309 mass_file=.false., just_read=just_read)
310 case ("mass_file")
311 call initialize_thickness_from_file(h, depth_tot, g, gv, us, pf, file_has_thickness=.true., &
312 mass_file=.true., just_read=just_read)
313 convert = .false.
314 case ("coord")
315 if (new_sim .and. useale) then
316 call ale_initthicknesstocoord( ale_csp, g, gv, dz, height_units=.true. )
317 elseif (new_sim) then
318 call mom_error(fatal, "MOM_initialize_state: USE_REGRIDDING must be True "//&
319 "for THICKNESS_CONFIG of 'coord'")
320 endif
321 case ("uniform"); call initialize_thickness_uniform(dz, depth_tot, g, gv, pf, &
322 just_read=just_read)
323 case ("list"); call initialize_thickness_list(dz, depth_tot, g, gv, us, pf, &
324 just_read=just_read)
325 case ("param"); call initialize_thickness_param(dz, depth_tot, g, gv, us, pf, &
326 just_read=just_read)
327 case ("DOME"); call dome_initialize_thickness(dz, depth_tot, g, gv, pf, &
328 just_read=just_read)
329 case ("ISOMIP"); call isomip_initialize_thickness(dz, depth_tot, g, gv, us, pf, tv, &
330 just_read=just_read)
331 case ("benchmark"); call benchmark_initialize_thickness(dz, depth_tot, g, gv, us, pf, &
332 tv%eqn_of_state, tv%P_Ref, just_read=just_read)
333 case ("Neverworld","Neverland"); call neverworld_initialize_thickness(dz, depth_tot, &
334 g, gv, us, pf, tv%P_Ref)
335 case ("search"); call initialize_thickness_search()
336 case ("circle_obcs"); call circle_obcs_initialize_thickness(dz, depth_tot, g, gv, us, pf, &
337 just_read=just_read)
338 case ("lock_exchange"); call lock_exchange_initialize_thickness(dz, g, gv, us, &
339 pf, just_read=just_read)
340 case ("external_gwave"); call external_gwave_initialize_thickness(dz, g, gv, us, &
341 pf, just_read=just_read)
342 case ("DOME2D"); call dome2d_initialize_thickness(dz, depth_tot, g, gv, us, pf, &
343 just_read=just_read)
344 case ("adjustment2d"); call adjustment_initialize_thickness(dz, g, gv, us, &
345 pf, just_read=just_read)
346 case ("sloshing"); call sloshing_initialize_thickness(dz, depth_tot, g, gv, us, pf, &
347 just_read=just_read)
348 case ("seamount"); call seamount_initialize_thickness(dz, depth_tot, g, gv, us, pf, &
349 just_read=just_read)
350 case ("dumbbell"); call dumbbell_initialize_thickness(dz, depth_tot, g, gv, us, pf, &
351 just_read=just_read)
352 case ("soliton"); call soliton_initialize_thickness(dz, depth_tot, g, gv, us, pf, &
353 just_read=just_read)
354 case ("phillips"); call phillips_initialize_thickness(dz, depth_tot, g, gv, us, pf, &
355 just_read=just_read)
356 case ("rossby_front")
357 call rossby_front_initialize_thickness(h, g, gv, us, pf, just_read=just_read)
358 convert = .false. ! Rossby_front initialization works directly in thickness units.
359 case ("USER"); call user_initialize_thickness(dz, g, gv, pf, &
360 just_read=just_read)
361 case default ; call mom_error(fatal, "MOM_initialize_state: "//&
362 "Unrecognized layer thickness configuration "//trim(h_config))
363 end select
364
365 ! Initialize temperature and salinity (T and S).
366 if ( use_temperature ) then
367 call get_param(pf, mdl, "TS_CONFIG", config, &
368 "A string that determines how the initial temperatures "//&
369 "and salinities are specified for a new run: \n"//&
370 " \t file - read velocities from the file specified \n"//&
371 " \t\t by (TS_FILE). \n"//&
372 " \t fit - find the temperatures that are consistent with \n"//&
373 " \t\t the layer densities and salinity S_REF. \n"//&
374 " \t TS_profile - use temperature and salinity profiles \n"//&
375 " \t\t (read from TS_FILE) to set layer densities. \n"//&
376 " \t benchmark - use the benchmark test case T & S. \n"//&
377 " \t linear - linear in logical layer space. \n"//&
378 " \t DOME2D - 2D DOME initialization. \n"//&
379 " \t ISOMIP - ISOMIP initialization. \n"//&
380 " \t adjustment2d - 2d lock exchange T/S ICs. \n"//&
381 " \t sloshing - sloshing mode T/S ICs. \n"//&
382 " \t seamount - no motion test with seamount ICs. \n"//&
383 " \t dumbbell - sloshing channel ICs. \n"//&
384 " \t rossby_front - a mixed layer front in thermal wind balance.\n"//&
385 " \t SCM_CVMix_tests - used in the SCM CVMix tests.\n"//&
386 " \t USER - call a user modified routine.", &
387 fail_if_missing=new_sim, do_not_log=just_read)
388! " \t baroclinic_zone - an analytic baroclinic zone. \n"//&
389
390 ! Check for incompatible THICKNESS_CONFIG and TS_CONFIG settings
391 if (new_sim .and. (.not.convert)) then ; select case (trim(config))
392 case ("DOME2D", "ISOMIP", "adjustment2d", "baroclinic_zone", "sloshing", &
393 "seamount", "dumbbell", "SCM_CVMix_tests", "dense")
394 call mom_error(fatal, "TS_CONFIG = "//trim(config)//" does not work with thicknesses "//&
395 "that have already been converted to thickness units, as is the case with "//&
396 "THICKNESS_CONFIG = "//trim(h_config)//".")
397 end select ; endif
398
399 select case (trim(config))
400 case ("fit"); call initialize_temp_salt_fit(tv%T, tv%S, g, gv, us, pf, &
401 eos, tv%P_Ref, just_read=just_read)
402 case ("file"); call initialize_temp_salt_from_file(tv%T, tv%S, g, gv, us, &
403 pf, just_read=just_read)
404 case ("benchmark"); call benchmark_init_temperature_salinity(tv%T, tv%S, &
405 g, gv, us, pf, eos, tv%P_Ref, just_read=just_read)
406 case ("TS_profile") ; call initialize_temp_salt_from_profile(tv%T, tv%S, &
407 g, gv, us, pf, just_read=just_read)
408 case ("linear"); call initialize_temp_salt_linear(tv%T, tv%S, g, gv, us, pf, &
409 just_read=just_read)
410 case ("DOME2D"); call dome2d_initialize_temperature_salinity (tv%T, tv%S, dz, &
411 g, gv, us, pf, just_read=just_read)
412 case ("ISOMIP"); call isomip_initialize_temperature_salinity (tv%T, tv%S, dz, &
413 depth_tot, g, gv, us, pf, eos, just_read=just_read)
414 case ("adjustment2d"); call adjustment_initialize_temperature_salinity ( tv%T, &
415 tv%S, dz, depth_tot, g, gv, us, pf, just_read=just_read)
416 case ("baroclinic_zone"); call baroclinic_zone_init_temperature_salinity( tv%T, &
417 tv%S, dz, depth_tot, g, gv, us, pf, just_read=just_read)
418 case ("sloshing"); call sloshing_initialize_temperature_salinity(tv%T, &
419 tv%S, dz, g, gv, us, pf, just_read=just_read)
420 case ("seamount"); call seamount_initialize_temperature_salinity(tv%T, &
421 tv%S, dz, g, gv, us, pf, just_read=just_read)
422 case ("dumbbell"); call dumbbell_initialize_temperature_salinity(tv%T, &
423 tv%S, dz, g, gv, us, pf, just_read=just_read)
424 case ("rossby_front")
425 if (convert .and. .not.just_read) call dz_to_thickness(dz, tv, h, g, gv, us)
426 call rossby_front_initialize_temperature_salinity ( tv%T, tv%S, h, &
427 g, gv, us, pf, just_read=just_read)
428 case ("SCM_CVMix_tests"); call scm_cvmix_tests_ts_init(tv%T, tv%S, dz, &
429 g, gv, us, pf, just_read=just_read)
430 case ("dense"); call dense_water_initialize_ts(g, gv, us, pf, tv%T, tv%S, &
431 dz, just_read=just_read)
432 case ("USER"); call user_init_temperature_salinity(tv%T, tv%S, g, gv, pf, &
433 just_read=just_read)
434 case default ; call mom_error(fatal, "MOM_initialize_state: "//&
435 "Unrecognized Temp & salt configuration "//trim(config))
436 end select
437 endif
438 endif ! not from_Z_file.
439
440 if (present(obc_for_bug)) then ; if (use_temperature .and. associated(obc_for_bug)) then
441 call get_param(pf, mdl, "ENABLE_BUGS_BY_DEFAULT", enable_bugs, &
442 default=.true., do_not_log=.true.) ! This is logged from MOM.F90.
443 ! Log this parameter later with the other OBC parameters.
444 call get_param(pf, mdl, "OBC_TS_RESERVOIR_INIT_BUG", obc_ts_reservoir_init_bug, &
445 "If true, set the OBC temperature and salinity reservoirs at the startup of a "//&
446 "new run from initial values that are set before remapping.", &
447 default=enable_bugs, do_not_log=.true.)
448 if (obc_ts_reservoir_init_bug) then
449 ! These calls should be moved down to join the OBC code, but doing so changes answers because
450 ! the temperatures and salinities can change due to the remapping and reading from the restarts.
451 call pass_var(tv%T, g%Domain, complete=.false.)
452 call pass_var(tv%S, g%Domain, complete=.true.)
453 call fill_temp_salt_segments(g, gv, us, obc_for_bug, tv)
454 endif
455 endif ; endif
456
457 ! Convert thicknesses from geometric distances in depth units to thickness units or mass-per-unit-area.
458 if (new_sim .and. convert) call dz_to_thickness(dz, tv, h, g, gv, us)
459
460 ! Handle the initial surface displacement under ice shelf
461 call get_param(pf, mdl, "DEPRESS_INITIAL_SURFACE", depress_sfc, &
462 "If true, depress the initial surface to avoid huge "//&
463 "tsunamis when a large surface pressure is applied.", &
464 default=.false., do_not_log=just_read)
465 call get_param(pf, mdl, "TRIM_IC_FOR_P_SURF", trim_ic_for_p_surf, &
466 "If true, cuts way the top of the column for initial conditions "//&
467 "at the depth where the hydrostatic pressure matches the imposed "//&
468 "surface pressure which is read from file.", default=.false., &
469 do_not_log=just_read)
470 if (depress_sfc .and. trim_ic_for_p_surf) call mom_error(fatal, "MOM_initialize_state: "//&
471 "DEPRESS_INITIAL_SURFACE and TRIM_IC_FOR_P_SURF are exclusive and cannot both be True")
472
473 if (new_sim .and. debug .and. (depress_sfc .or. trim_ic_for_p_surf)) &
474 call hchksum(h, "Pre-depress: h ", g%HI, haloshift=1, unscale=gv%H_to_MKS)
475
476 ! Remove the mass that would be displaced by an ice shelf or inverse barometer.
477 if (depress_sfc) then
478 call depress_surface(h, g, gv, us, pf, tv, just_read=just_read)
479 elseif (trim_ic_for_p_surf) then
480 call trim_for_ice(pf, g, gv, us, ale_csp, tv, h, just_read=just_read)
481 elseif (new_sim .and. use_ice_shelf .and. present(mass_shelf)) then
482 call calc_sfc_displacement(pf, g, gv, us, mass_shelf, tv, h)
483 endif
484
485 ! Perhaps we want to run the regridding coordinate generator for multiple
486 ! iterations here so the initial grid is consistent with the coordinate
487 if (useale) then
488 call get_param(pf, mdl, "REGRID_ACCELERATE_INIT", regrid_accelerate, &
489 "If true, runs REGRID_ACCELERATE_ITERATIONS iterations of the regridding "//&
490 "algorithm to push the initial grid to be consistent with the initial "//&
491 "condition. Useful only for state-based and iterative coordinates.", &
492 default=.false., do_not_log=just_read)
493 if (regrid_accelerate) then
494 call get_param(pf, mdl, "REGRID_ACCELERATE_ITERATIONS", regrid_iterations, &
495 "The number of regridding iterations to perform to generate "//&
496 "an initial grid that is consistent with the initial conditions.", &
497 default=1, do_not_log=just_read)
498
499 call get_param(pf, mdl, "DT", dt, "Timestep", &
500 units="s", scale=us%s_to_T, fail_if_missing=.true.)
501
502 if (new_sim .and. debug) &
503 call hchksum(h, "Pre-ALE_regrid: h ", g%HI, haloshift=1, unscale=gv%H_to_MKS)
504 ! In this call, OBC_for_remap is only used for the directions of OBCs when setting thicknesses at
505 ! velocity points.
506 call ale_regrid_accelerated(ale_csp, g, gv, us, h, tv, regrid_iterations, u, v, obc_for_remap, &
507 tracer_reg, dt=dt, initial=.true.)
508 endif
509 endif
510
511 ! The thicknesses in halo points might be needed to initialize the velocities.
512 if (new_sim) call pass_var(h, g%Domain)
513
514 ! Initialize velocity components, u and v
515 call get_param(pf, mdl, "VELOCITY_CONFIG", config, &
516 "A string that determines how the initial velocities "//&
517 "are specified for a new run: \n"//&
518 " \t file - read velocities from the file specified \n"//&
519 " \t\t by (VELOCITY_FILE). \n"//&
520 " \t zero - the fluid is initially at rest. \n"//&
521 " \t uniform - the flow is uniform (determined by\n"//&
522 " \t\t parameters INITIAL_U_CONST and INITIAL_V_CONST).\n"//&
523 " \t rossby_front - a mixed layer front in thermal wind balance.\n"//&
524 " \t soliton - Equatorial Rossby soliton.\n"//&
525 " \t USER - call a user modified routine.", default="zero", &
526 do_not_log=just_read)
527 select case (trim(config))
528 case ("file"); call initialize_velocity_from_file(u, v, g, gv, us, pf, just_read)
529 case ("zero"); call initialize_velocity_zero(u, v, g, gv, pf, just_read)
530 case ("uniform"); call initialize_velocity_uniform(u, v, g, gv, us, pf, just_read)
531 case ("circular"); call initialize_velocity_circular(u, v, g, gv, us, pf, just_read)
532 case ("phillips"); call phillips_initialize_velocity(u, v, g, gv, us, pf, just_read)
533 case ("rossby_front"); call rossby_front_initialize_velocity(u, v, h, &
534 g, gv, us, pf, just_read)
535 case ("soliton"); call soliton_initialize_velocity(u, v, g, gv, us, pf, just_read)
536 case ("USER"); call user_initialize_velocity(u, v, g, gv, us, pf, just_read)
537 case default ; call mom_error(fatal, "MOM_initialize_state: "//&
538 "Unrecognized velocity configuration "//trim(config))
539 end select
540
541 if (new_sim) call pass_vector(u, v, g%Domain)
542 if (debug .and. new_sim) then
543 call uvchksum("MOM_initialize_state [uv]", u, v, g%HI, haloshift=1, unscale=us%L_T_to_m_s)
544 endif
545
546 ! This is the end of the block of code that might have initialized fields
547 ! internally at the start of a new run.
548
549 ! Initialized assimilative incremental update (oda_incupd) structure and
550 ! register restart.
551 call get_param(pf, mdl, "ODA_INCUPD", use_oda_incupd, &
552 "If true, oda incremental updates will be applied "//&
553 "everywhere in the domain.", default=.false.)
554 if (use_oda_incupd) then
555 call restart_registry_lock(restart_cs, unlocked=.true.)
556 call initialize_oda_incupd_fixed(g, gv, us, oda_incupd_csp, restart_cs)
557 call restart_registry_lock(restart_cs)
558 endif
559
560 if (.not.new_sim) then ! This block restores the state from a restart file.
561 ! This line calls a subroutine that reads the initial conditions
562 ! from a previously generated file.
563 call restore_state(dirs%input_filename, dirs%restart_input_dir, time, g, restart_cs)
564 if (present(time_in)) then
565 if (verify_restart_time .and. (time /= time_in)) call mom_error(fatal, &
566 "MOM6 attempted to restart from a file from a different time than given by Time_in.")
567 time = time_in
568 endif
569 call get_param(pf, mdl, "ROTATE_INDEX", rotate_index, &
570 "Enable rotation of the horizontal indices.", &
571 default=.false., debuggingparam=.true., do_not_log=.true.)
572 if (rotate_index) then
573 ! This model is using a rotated grid, so the unrotated variables used here have not been set yet.
574 call copy_restart_var(h, "h", restart_cs, .true.)
575 call copy_restart_vector(u, v, "u", "v", restart_cs, .true.)
576 if ( use_temperature ) then
577 call copy_restart_var(tv%T, "Temp", restart_cs, .true.)
578 call copy_restart_var(tv%S, "Salt", restart_cs, .true.)
579 endif
580 endif
581 endif
582
583 if ( use_temperature ) then
584 call pass_var(tv%T, g%Domain, complete=.false.)
585 call pass_var(tv%S, g%Domain, complete=.false.)
586 endif
587 call pass_var(h, g%Domain)
588
589 if (debug) then
590 call hchksum(h, "MOM_initialize_state: h ", g%HI, haloshift=1, unscale=gv%H_to_MKS)
591 if ( use_temperature ) call hchksum(tv%T, "MOM_initialize_state: T ", g%HI, haloshift=1, unscale=us%C_to_degC)
592 if ( use_temperature ) call hchksum(tv%S, "MOM_initialize_state: S ", g%HI, haloshift=1, unscale=us%S_to_ppt)
593 if ( use_temperature .and. debug_layers) then ; do k=1,nz
594 write(mesg,'("MOM_IS: T[",I0,"]")') k
595 call hchksum(tv%T(:,:,k), mesg, g%HI, haloshift=1, unscale=us%C_to_degC)
596 write(mesg,'("MOM_IS: S[",I0,"]")') k
597 call hchksum(tv%S(:,:,k), mesg, g%HI, haloshift=1, unscale=us%S_to_ppt)
598 enddo ; endif
599 endif
600
601 call get_param(pf, mdl, "SPONGE", use_sponge, &
602 "If true, sponges may be applied anywhere in the domain. "//&
603 "The exact location and properties of those sponges are "//&
604 "specified via SPONGE_CONFIG.", default=.false.)
605 if ( use_sponge ) then
606 call get_param(pf, mdl, "SPONGE_CONFIG", config, &
607 "A string that sets how the sponges are configured: \n"//&
608 " \t file - read sponge properties from the file \n"//&
609 " \t\t specified by (SPONGE_FILE).\n"//&
610 " \t ISOMIP - apply ale sponge in the ISOMIP case \n"//&
611 " \t RGC - apply sponge in the rotating_gravity_current case \n"//&
612 " \t DOME - use a slope and channel configuration for the \n"//&
613 " \t\t DOME sill-overflow test case. \n"//&
614 " \t BFB - Sponge at the southern boundary of the domain\n"//&
615 " \t\t for buoyancy-forced basin case.\n"//&
616 " \t USER - call a user modified routine.", default="file")
617 select case (trim(config))
618 case ("DOME"); call dome_initialize_sponges(g, gv, us, tv, depth_tot, pf, sponge_csp)
619 case ("DOME2D"); call dome2d_initialize_sponges(g, gv, us, tv, depth_tot, pf, useale, &
620 sponge_csp, ale_sponge_csp)
621 case ("ISOMIP"); call isomip_initialize_sponges(g, gv, us, tv, depth_tot, pf, useale, &
622 sponge_csp, ale_sponge_csp)
623 case ("RGC"); call rgc_initialize_sponges(g, gv, us, tv, u, v, depth_tot, pf, useale, &
624 sponge_csp, ale_sponge_csp)
625 case ("USER"); call user_initialize_sponges(g, gv, use_temperature, tv, pf, sponge_csp, h)
626 case ("BFB"); call bfb_initialize_sponges_southonly(g, gv, us, use_temperature, tv, depth_tot, pf, &
627 sponge_csp, h)
628 case ("DUMBBELL"); call dumbbell_initialize_sponges(g, gv, us, tv, h, depth_tot, pf, useale, &
629 sponge_csp, ale_sponge_csp)
630 case ("phillips"); call phillips_initialize_sponges(g, gv, us, tv, pf, sponge_csp, h)
631 case ("dense"); call dense_water_initialize_sponges(g, gv, us, tv, depth_tot, pf, useale, &
632 sponge_csp, ale_sponge_csp)
633 case ("file"); call initialize_sponges_file(g, gv, us, use_temperature, tv, u, v, depth_tot, pf, &
634 sponge_csp, ale_sponge_csp, time)
635 case default ; call mom_error(fatal, "MOM_initialize_state: "//&
636 "Unrecognized sponge configuration "//trim(config))
637 end select
638 endif
639
640 ! Set-up of data Assimilation with incremental update
641 if (use_oda_incupd) then
642 call initialize_oda_incupd_file(g, gv, us, use_temperature, tv, h, u, v, &
643 pf, oda_incupd_csp, restart_cs, time)
644 endif
645
646 call calltree_leave('MOM_initialize_state()')
647
648end subroutine mom_initialize_state
649
650subroutine mom_initialize_obcs(h, tv, OBC, Time, G, GV, US, PF, restart_CS, tracer_Reg)
651 type(ocean_grid_type), intent(inout) :: g !< The ocean's grid structure.
652 type(verticalgrid_type), intent(in) :: gv !< The ocean's vertical grid structure.
653 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
654 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
655 intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2]
656 type(thermo_var_ptrs), intent(inout) :: tv !< A structure pointing to various thermodynamic
657 !! variables
658 type(ocean_obc_type), pointer :: obc !< The open boundary condition control structure.
659 type(time_type), intent(in) :: time !< Time at the start of the run segment.
660 type(param_file_type), intent(in) :: pf !< A structure indicating the open file to parse
661 !! for model parameter values.
662 type(mom_restart_cs), intent(inout) :: restart_cs !< MOM restart control structure
663 type(tracer_registry_type), pointer :: tracer_reg !< A pointer to the tracer registry
664
665 ! Local variables
666 character(len=200) :: config
667 logical :: enable_bugs ! If true, the defaults for recently added bug-fix flags are set to
668 ! recreate the bugs, or if false bugs are only used if actively selected.
669 logical :: debug ! If true, write debugging output.
670 logical :: debug_obc ! If true, do additional calls resetting values to help debug the correctness
671 ! of the open boundary condition code.
672 logical :: obc_ts_reservoir_init_bug ! If true, set the OBC temperature and salinity reservoirs
673 ! at the startup of a new run from initial values that are set before remapping.
674 logical :: obc_reservoir_init_bug ! If true, set the OBC tracer reservoirs at the startup of a new
675 ! run from the interior tracer concentrations regardless of properties that
676 ! may be explicitly specified for the reservoir concentrations.
677
678 call calltree_enter('MOM_initialize_OBCs()')
679 if (associated(obc)) then
680 call get_param(pf, mdl, "DEBUG", debug, default=.false.)
681 call get_param(pf, mdl, "OBC_DEBUGGING_TESTS", debug_obc, &
682 "If true, do additional calls resetting values to help verify the correctness "//&
683 "of the open boundary condition code.", default=.false., &
684 do_not_log=.true., old_name="DEBUG_OBC", debuggingparam=.true.)
685 call get_param(pf, mdl, "ENABLE_BUGS_BY_DEFAULT", enable_bugs, &
686 default=.true., do_not_log=.true.) ! This is logged from MOM.F90.
687 call get_param(pf, mdl, "OBC_TS_RESERVOIR_INIT_BUG", obc_ts_reservoir_init_bug, &
688 "If true, set the OBC temperature and salinity reservoirs at the startup of a "//&
689 "new run from initial values that are set before remapping.", default=enable_bugs)
690 if (associated(tv%T) .and. (.not.obc_ts_reservoir_init_bug)) then
691 ! Store the updated temperatures and salinities at the open boundaries, noting that they may
692 ! still be updated by the calls in the next 50 lines, so the code setting the tracer
693 ! reservoir values will come later in the calling routine.
694 call fill_temp_salt_segments(g, gv, us, obc, tv)
695 endif
696 call get_param(pf, mdl, "OBC_RESERVOIR_INIT_BUG", obc_reservoir_init_bug, &
697 "If true, set the OBC tracer reservoirs at the startup of a new run from the "//&
698 "interior tracer concentrations regardless of properties that may be explicitly "//&
699 "specified for the reservoir concentrations.", default=enable_bugs)
700 if (obc_reservoir_init_bug .and. associated(tv%T) .and. is_new_run(restart_cs)) then
701 ! Set up OBC%trex_x and OBC%tres_y as they have not been read from a restart file.
702 ! When OBC_RESERVOIR_INIT_BUG is false, setup_OBC_tracer_reservoirs() is called from initialize_MOM
703 ! after all tracer package initialization is finished and grid rotation has been dealt with.
704 call setup_obc_tracer_reservoirs(g, gv, obc)
705 ! Ensure that the values of the tracer reservoirs that have just been set will not be revised.
706 call set_initialized_obc_tracer_reservoirs(g, obc, restart_cs)
707 endif
708
709 ! This controls user code for setting open boundary data
710 call get_param(pf, mdl, "OBC_USER_CONFIG", config, &
711 "A string that sets how the user code is invoked to set open boundary data: \n"//&
712 " DOME - specified inflow on northern boundary\n"//&
713 " dyed_channel - supercritical with dye on the inflow boundary\n"//&
714 " dyed_obcs - circle_obcs with dyes on the open boundaries\n"//&
715 " Kelvin - barotropic Kelvin wave forcing on the western boundary\n"//&
716 " shelfwave - Flather with shelf wave forcing on western boundary\n"//&
717 " supercritical - now only needed here for the allocations\n"//&
718 " tidal_bay - Flather with tidal forcing on eastern boundary\n"//&
719 " USER - user specified", default="none")
720 if (trim(config) == "DOME") then
721 call dome_set_obc_data(obc, tv, g, gv, us, pf, tracer_reg)
722 elseif (trim(config) == "dyed_channel") then
723 call dyed_channel_set_obc_tracer_data(obc, g, gv, pf, tracer_reg)
724 obc%update_OBC = .true.
725 elseif (trim(config) == "dyed_obcs") then
726 call dyed_obcs_set_obc_data(obc, g, gv, pf, tracer_reg)
727 elseif (trim(config) == "Kelvin") then
728 obc%update_OBC = .true.
729 elseif (trim(config) == "shelfwave") then
730 obc%update_OBC = .true.
731 elseif (lowercase(trim(config)) == "supercritical") then
732 call supercritical_set_obc_data(obc, g, gv, us, pf)
733 elseif (trim(config) == "tidal_bay") then
734 obc%update_OBC = .true.
735 elseif (trim(config) == "USER") then
736 call user_set_obc_data(obc, tv, g, gv, pf, tracer_reg)
737 elseif (.not. trim(config) == "none") then
738 call mom_error(fatal, "The open boundary conditions specified by "//&
739 "OBC_USER_CONFIG = "//trim(config)//" have not been fully implemented.")
740 endif
741
742 if (debug) then
743 call hchksum(g%mask2dT, 'MOM_initialize_OBCs: mask2dT ', g%HI)
744 call uvchksum('MOM_initialize_OBCs: mask2dC[uv]', g%mask2dCu, g%mask2dCv, g%HI)
745 call qchksum(g%mask2dBu, 'MOM_initialize_OBCs: mask2dBu ', g%HI)
746 endif
747 if (debug_obc) call open_boundary_test_extern_h(g, gv, obc, h)
748
749 if (obc%use_h_res) &
750 call fill_thickness_segments(g, gv, us, obc, h)
751 endif
752
753 call calltree_leave('MOM_initialize_OBCs()')
754
755end subroutine mom_initialize_obcs
756
757!> Reads the layer thicknesses or interface heights from a file.
758subroutine initialize_thickness_from_file(h, depth_tot, G, GV, US, param_file, file_has_thickness, &
759 just_read, mass_file)
760 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
761 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
762 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
763 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
764 intent(out) :: h !< The thickness that is being initialized, in height
765 !! or thickness units, depending on the value of
766 !! mass_file [Z ~> m] or [H ~> m or kg m-2].
767 real, dimension(SZI_(G),SZJ_(G)), &
768 intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
769 type(param_file_type), intent(in) :: param_file !< A structure indicating the open file
770 !! to parse for model parameter values.
771 logical, intent(in) :: file_has_thickness !< If true, this file contains layer
772 !! thicknesses; otherwise it contains
773 !! interface heights.
774 logical, intent(in) :: just_read !< If true, this call will only read
775 !! parameters without changing h.
776 logical, intent(in) :: mass_file !< If true, this file contains layer thicknesses in
777 !! units of mass per unit area.
778
779 ! Local variables
780 real :: eta(SZI_(G),SZJ_(G),SZK_(GV)+1) ! Interface heights, in depth units [Z ~> m].
781 real :: h_rescale ! A factor by which to rescale the initial thickness variable in the input
782 ! file to convert it to units of m [various]
783 real :: eta_rescale ! A factor by which to rescale the initial interface heights to convert
784 ! them to units of m or correct sign conventions to positive upward [various]
785 real :: h_tolerance ! A parameter that controls the tolerance when adjusting the
786 ! thickness to fit the bathymetry [Z ~> m].
787 real :: tol_dz_bot ! A tolerance for detecting inconsistent bottom depths when
788 ! correct_thickness is false [Z ~> m]
789 integer :: inconsistent ! The total number of cells with in consistent topography and layer thicknesses.
790 logical :: correct_thickness
791 character(len=40) :: mdl = "initialize_thickness_from_file" ! This subroutine's name.
792 character(len=200) :: filename, thickness_file, inputdir, mesg ! Strings for file/path
793 character(len=80) :: eta_var ! The interface height variable name in the input file
794 character(len=80) :: h_var ! The thickness variable name in the input file
795 integer :: i, j, k, is, ie, js, je, nz
796
797 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
798
799 if (.not.just_read) &
800 call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
801
802 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".", do_not_log=just_read)
803 inputdir = slasher(inputdir)
804 call get_param(param_file, mdl, "THICKNESS_FILE", thickness_file, &
805 "The name of the thickness file.", &
806 fail_if_missing=.not.just_read, do_not_log=just_read)
807
808 filename = trim(thickness_file)
809 if (scan(thickness_file, "/") == 0) then ! prepend inputdir if only a filename is given
810 filename = trim(inputdir)//trim(thickness_file)
811 endif
812 if (.not.just_read) call log_param(param_file, mdl, "INPUTDIR/THICKNESS_FILE", filename)
813
814 if ((.not.just_read) .and. (.not.file_exists(filename, g%Domain))) call mom_error(fatal, &
815 " initialize_thickness_from_file: Unable to open "//trim(filename))
816
817 if (file_has_thickness) then
818 call get_param(param_file, mdl, "THICKNESS_IC_VAR", h_var, &
819 "The variable name for layer thickness initial conditions.", &
820 default="h", do_not_log=just_read)
821 call get_param(param_file, mdl, "THICKNESS_IC_RESCALE", h_rescale, &
822 'A factor by which to rescale the initial thicknesses in the input file to '//&
823 'convert them to units of kg/m2 (if THICKNESS_CONFIG="mass_file") or m.', &
824 default=1.0, units="various", do_not_log=just_read)
825 if (just_read) return ! All run-time parameters have been read, so return.
826
827 if (mass_file) then
828 h_rescale = h_rescale*gv%kg_m2_to_H
829 else
830 h_rescale = h_rescale*us%m_to_Z
831 endif
832 call mom_read_data(filename, h_var, h(:,:,:), g%Domain, scale=h_rescale)
833 else
834 call get_param(param_file, mdl, "ADJUST_THICKNESS", correct_thickness, &
835 "If true, all mass below the bottom removed if the "//&
836 "topography is shallower than the thickness input file "//&
837 "would indicate.", default=.false., do_not_log=just_read)
838 if (correct_thickness) then
839 call get_param(param_file, mdl, "THICKNESS_TOLERANCE", h_tolerance, &
840 "A parameter that controls the tolerance when adjusting the "//&
841 "thickness to fit the bathymetry. Used when ADJUST_THICKNESS=True.", &
842 units="m", default=0.1, scale=us%m_to_Z, do_not_log=just_read)
843 endif
844 call get_param(param_file, mdl, "DZ_BOTTOM_TOLERANCE", tol_dz_bot, &
845 "A tolerance for detecting inconsistent topography and input layer "//&
846 "thicknesses when ADJUST_THICKNESS is false.", &
847 units="m", default=1.0, scale=us%m_to_Z, &
848 do_not_log=(just_read.or.correct_thickness))
849 call get_param(param_file, mdl, "INTERFACE_IC_VAR", eta_var, &
850 "The variable name for initial conditions for interface heights "//&
851 "relative to mean sea level, positive upward unless otherwise rescaled.", &
852 default="eta", do_not_log=just_read)
853 call get_param(param_file, mdl, "INTERFACE_IC_RESCALE", eta_rescale, &
854 "A factor by which to rescale the initial interface heights to convert "//&
855 "them to units of m or correct sign conventions to positive upward.", &
856 default=1.0, units="various", do_not_log=just_read)
857 if (just_read) return ! All run-time parameters have been read, so return.
858
859 call mom_read_data(filename, eta_var, eta(:,:,:), g%Domain, scale=us%m_to_Z*eta_rescale)
860
861 if (correct_thickness) then
862 call adjustetatofitbathymetry(g, gv, us, eta, h, h_tolerance, dz_ref_eta=g%Z_ref)
863 else
864 do k=nz,1,-1 ; do j=js,je ; do i=is,ie
865 if (eta(i,j,k) < (eta(i,j,k+1) + gv%Angstrom_Z)) then
866 eta(i,j,k) = eta(i,j,k+1) + gv%Angstrom_Z
867 h(i,j,k) = gv%Angstrom_Z
868 else
869 h(i,j,k) = eta(i,j,k) - eta(i,j,k+1)
870 endif
871 enddo ; enddo ; enddo
872
873 inconsistent = 0
874 do j=js,je ; do i=is,ie
875 if (abs(eta(i,j,nz+1) + depth_tot(i,j)) > tol_dz_bot) &
876 inconsistent = inconsistent + 1
877 enddo ; enddo
878 call sum_across_pes(inconsistent)
879
880 if ((inconsistent > 0) .and. (is_root_pe())) then
881 write(mesg,'("Thickness initial conditions are inconsistent ",'// &
882 '"with topography in ",I0," places.")') inconsistent
883 call mom_error(warning, mesg)
884 endif
885 endif
886
887 endif
888 call calltree_leave(trim(mdl)//'()')
890
891!> Adjust interface heights to fit the bathymetry and diagnose layer thickness.
892!!
893!! If the bottom most interface is below the topography then the bottom-most
894!! layers are contracted to ANGSTROM thickness (which may be 0).
895!! If the bottom most interface is above the topography then the entire column
896!! is dilated (expanded) to fill the void.
897subroutine adjustetatofitbathymetry(G, GV, US, eta, h, ht, dZ_ref_eta)
898 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
899 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
900 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
901 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: eta !< Interface heights [Z ~> m].
902 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thicknesses [Z ~> m]
903 real, intent(in) :: ht !< Tolerance to exceed adjustment
904 !! criteria [Z ~> m]
905 real, optional, intent(in) :: dZ_ref_eta !< The difference between the
906 !! reference heights for bathyT and
907 !! eta [Z ~> m], 0 by default.
908 ! Local variables
909 integer :: i, j, k, is, ie, js, je, nz, contractions, dilations
910 real :: dilate ! A factor by which the column is dilated [nondim]
911 real :: dZ_ref ! The difference in the reference heights for G%bathyT and eta [Z ~> m]
912 character(len=100) :: mesg
913
914 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
915 dz_ref = 0.0 ; if (present(dz_ref_eta)) dz_ref = dz_ref_eta
916
917 contractions = 0
918 do j=js,je ; do i=is,ie
919 if (-eta(i,j,nz+1) > (g%bathyT(i,j) + dz_ref) + ht) then
920 eta(i,j,nz+1) = -(g%bathyT(i,j) + dz_ref)
921 contractions = contractions + 1
922 endif
923 enddo ; enddo
924 call sum_across_pes(contractions)
925 if ((contractions > 0) .and. (is_root_pe())) then
926 write(mesg,'("Thickness initial conditions were contracted ",'// &
927 '"to fit topography in ",I0," places.")') contractions
928 call mom_error(warning, 'adjustEtaToFitBathymetry: '//mesg)
929 endif
930
931 ! To preserve previous answers in non-Boussinesq cases, delay converting
932 ! thicknesses to units of H until the end of this routine.
933 do k=nz,1,-1 ; do j=js,je ; do i=is,ie
934 ! Collapse layers to thinnest possible if the thickness less than
935 ! the thinnest possible (or negative).
936 if (eta(i,j,k) < (eta(i,j,k+1) + gv%Angstrom_Z)) then
937 eta(i,j,k) = eta(i,j,k+1) + gv%Angstrom_Z
938 h(i,j,k) = gv%Angstrom_Z
939 else
940 h(i,j,k) = (eta(i,j,k) - eta(i,j,k+1))
941 endif
942 enddo ; enddo ; enddo
943
944 dilations = 0
945 do j=js,je ; do i=is,ie
946 ! The whole column is dilated to accommodate deeper topography than
947 ! the bathymetry would indicate.
948 ! This should be... if ((G%mask2dt(i,j)*(eta(i,j,1)-eta(i,j,nz+1)) > 0.0) .and. &
949 if (-eta(i,j,nz+1) < (g%bathyT(i,j) + dz_ref) - ht) then
950 dilations = dilations + 1
951 if (eta(i,j,1) <= eta(i,j,nz+1)) then
952 do k=1,nz ; h(i,j,k) = (eta(i,j,1) + (g%bathyT(i,j) + dz_ref)) / real(nz) ; enddo
953 else
954 dilate = (eta(i,j,1) + (g%bathyT(i,j) + dz_ref)) / (eta(i,j,1) - eta(i,j,nz+1))
955 do k=1,nz ; h(i,j,k) = h(i,j,k) * dilate ; enddo
956 endif
957 do k=nz,2,-1 ; eta(i,j,k) = eta(i,j,k+1) + h(i,j,k) ; enddo
958 endif
959 enddo ; enddo
960
961
962 call sum_across_pes(dilations)
963 if ((dilations > 0) .and. (is_root_pe())) then
964 write(mesg,'("Thickness initial conditions were dilated ",'// &
965 '"to fit topography in ",I0," places.")') dilations
966 call mom_error(warning, 'adjustEtaToFitBathymetry: '//mesg)
967 endif
968
969end subroutine adjustetatofitbathymetry
970
971!> Initializes thickness to be uniform
972subroutine initialize_thickness_uniform(h, depth_tot, G, GV, param_file, just_read)
973 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
974 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
975 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
976 intent(out) :: h !< The thickness that is being initialized [Z ~> m]
977 real, dimension(SZI_(G),SZJ_(G)), &
978 intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
979 type(param_file_type), intent(in) :: param_file !< A structure indicating the open file
980 !! to parse for model parameter values.
981 logical, intent(in) :: just_read !< If true, this call will only read
982 !! parameters without changing h.
983 ! Local variables
984 character(len=40) :: mdl = "initialize_thickness_uniform" ! This subroutine's name.
985 real :: e0(SZK_(GV)+1) ! The resting interface heights [Z ~> m], usually
986 ! negative because it is positive upward.
987 real :: eta1D(SZK_(GV)+1)! Interface height relative to the sea surface,
988 ! positive upward [Z ~> m].
989 integer :: i, j, k, is, ie, js, je, nz
990
991 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
992
993 if (just_read) return ! This subroutine has no run-time parameters.
994
995 call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
996
997 if (g%max_depth<=0.) call mom_error(fatal,"initialize_thickness_uniform: "// &
998 "MAXIMUM_DEPTH has a nonsensical value! Was it set?")
999
1000 do k=1,nz
1001 e0(k) = -g%max_depth * real(k-1) / real(nz)
1002 enddo
1003
1004 do j=js,je ; do i=is,ie
1005 ! This sets the initial thickness (in m) of the layers. The
1006 ! thicknesses are set to insure that: 1. each layer is at least an
1007 ! Angstrom thick, and 2. the interfaces are where they should be
1008 ! based on the resting depths and interface height perturbations,
1009 ! as long at this doesn't interfere with 1.
1010 eta1d(nz+1) = -depth_tot(i,j)
1011 do k=nz,1,-1
1012 eta1d(k) = e0(k)
1013 if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z)) then
1014 eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
1015 h(i,j,k) = gv%Angstrom_Z
1016 else
1017 h(i,j,k) = eta1d(k) - eta1d(k+1)
1018 endif
1019 enddo
1020 enddo ; enddo
1021
1022 call calltree_leave(trim(mdl)//'()')
1023end subroutine initialize_thickness_uniform
1024
1025!> Initialize thickness from a 1D list
1026subroutine initialize_thickness_list(h, depth_tot, G, GV, US, param_file, just_read)
1027 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1028 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1029 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1030 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1031 intent(out) :: h !< The thickness that is being initialized [Z ~> m]
1032 real, dimension(SZI_(G),SZJ_(G)), &
1033 intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
1034 type(param_file_type), intent(in) :: param_file !< A structure indicating the open file
1035 !! to parse for model parameter values.
1036 logical, intent(in) :: just_read !< If true, this call will only read
1037 !! parameters without changing h.
1038 ! Local variables
1039 character(len=40) :: mdl = "initialize_thickness_list" ! This subroutine's name.
1040 real :: e0(SZK_(GV)+1) ! The resting interface heights, in depth units [Z ~> m],
1041 ! usually negative because it is positive upward.
1042 real :: eta1D(SZK_(GV)+1)! Interface height relative to the sea surface
1043 ! positive upward, in depth units [Z ~> m].
1044 character(len=200) :: filename, eta_file, inputdir ! Strings for file/path
1045 character(len=72) :: eta_var ! The interface height variable name in the input file
1046 integer :: i, j, k, is, ie, js, je, nz
1047
1048 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
1049
1050 call get_param(param_file, mdl, "INTERFACE_IC_FILE", eta_file, &
1051 "The file from which horizontal mean initial conditions "//&
1052 "for interface depths can be read.", fail_if_missing=.true.)
1053 call get_param(param_file, mdl, "INTERFACE_IC_VAR", eta_var, &
1054 "The variable name for horizontal mean initial conditions "//&
1055 "for interface depths relative to mean sea level.", &
1056 default="eta")
1057
1058 if (just_read) return
1059
1060 call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
1061
1062 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
1063 filename = trim(slasher(inputdir))//trim(eta_file)
1064 call log_param(param_file, mdl, "INPUTDIR/INTERFACE_IC_FILE", filename)
1065
1066 e0(:) = 0.0
1067 call mom_read_data(filename, eta_var, e0(:), scale=us%m_to_Z)
1068
1069 if ((abs(e0(1)) - 0.0) > 0.001) then
1070 ! This list probably starts with the interior interface, so shift it up.
1071 do k=nz+1,2,-1 ; e0(k) = e0(k-1) ; enddo
1072 e0(1) = 0.0
1073 endif
1074
1075 if (e0(2) > e0(1)) then ! Switch to the convention for interface heights increasing upward.
1076 do k=1,nz ; e0(k) = -e0(k) ; enddo
1077 endif
1078
1079 do j=js,je ; do i=is,ie
1080 ! This sets the initial thickness (in m) of the layers. The
1081 ! thicknesses are set to insure that: 1. each layer is at least an
1082 ! Angstrom thick, and 2. the interfaces are where they should be
1083 ! based on the resting depths and interface height perturbations,
1084 ! as long at this doesn't interfere with 1.
1085 eta1d(nz+1) = -depth_tot(i,j)
1086 do k=nz,1,-1
1087 eta1d(k) = e0(k)
1088 if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z)) then
1089 eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
1090 h(i,j,k) = gv%Angstrom_Z
1091 else
1092 h(i,j,k) = eta1d(k) - eta1d(k+1)
1093 endif
1094 enddo
1095 enddo ; enddo
1096
1097 call calltree_leave(trim(mdl)//'()')
1098end subroutine initialize_thickness_list
1099
1100!> Initializes thickness based on a run-time parameter with nominal thickness
1101!! for each layer
1102subroutine initialize_thickness_param(h, depth_tot, G, GV, US, param_file, just_read)
1103 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1104 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1105 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1106 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1107 intent(out) :: h !< The thickness that is being initialized [Z ~> m]
1108 real, dimension(SZI_(G),SZJ_(G)), &
1109 intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
1110 type(param_file_type), intent(in) :: param_file !< A structure indicating the open file
1111 !! to parse for model parameter values.
1112 logical, intent(in) :: just_read !< If true, this call will only read
1113 !! parameters without changing h.
1114 ! Local variables
1115 character(len=40) :: mdl = "initialize_thickness_param" ! This subroutine's name.
1116 real :: e0(SZK_(GV)+1) ! The resting interface heights [Z ~> m], usually
1117 ! negative because it is positive upward.
1118 real :: eta1D(SZK_(GV)+1)! Interface height relative to the sea surface,
1119 ! positive upward [Z ~> m].
1120 real :: dz(SZK_(GV)) ! The nominal initial layer thickness [Z ~> m], usually
1121 real :: h0_def(SZK_(GV)) ! Uniform default values for dz [Z ~> m], usually
1122 integer :: i, j, k, is, ie, js, je, nz
1123
1124 call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
1125 if (g%max_depth<=0.) call mom_error(fatal, "initialize_thickness_param: "// &
1126 "MAXIMUM_DEPTH has a nonsensical value! Was it set?")
1127
1128 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
1129
1130 h0_def(:) = ( g%max_depth / real(nz) ) * us%Z_to_m
1131 call get_param(param_file, mdl, "THICKNESS_INIT_VALUES", dz, &
1132 "A list of nominal thickness for each layer to initialize with", &
1133 units="m", scale=us%m_to_Z, defaults=h0_def, do_not_log=just_read)
1134 if (just_read) return ! This subroutine has no run-time parameters.
1135
1136 e0(nz+1) = -g%max_depth
1137 do k=nz, 1, -1
1138 e0(k) = e0(k+1) + dz(k)
1139 enddo
1140
1141 do j=js,je ; do i=is,ie
1142 ! This sets the initial thickness (in m) of the layers. The
1143 ! thicknesses are set to insure that: 1. each layer is at least an
1144 ! Angstrom thick, and 2. the interfaces are where they should be
1145 ! based on the resting depths and interface height perturbations,
1146 ! as long at this doesn't interfere with 1.
1147 eta1d(nz+1) = -depth_tot(i,j)
1148 do k=nz,1,-1
1149 eta1d(k) = e0(k)
1150 if (eta1d(k) < (eta1d(k+1) + gv%Angstrom_Z)) then
1151 eta1d(k) = eta1d(k+1) + gv%Angstrom_Z
1152 h(i,j,k) = gv%Angstrom_Z
1153 else
1154 h(i,j,k) = eta1d(k) - eta1d(k+1)
1155 endif
1156 enddo
1157 enddo ; enddo
1158
1159 call calltree_leave(trim(mdl)//'()')
1160end subroutine initialize_thickness_param
1161
1162!> Search density space for location of layers (not implemented!)
1163subroutine initialize_thickness_search
1164 call mom_error(fatal," MOM_state_initialization.F90, initialize_thickness_search: NOT IMPLEMENTED")
1165end subroutine initialize_thickness_search
1166
1167!> Depress the sea-surface based on an initial condition file
1168subroutine depress_surface(h, G, GV, US, param_file, tv, just_read, z_top_shelf)
1169 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1170 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1171 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1172 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1173 intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2]
1174 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
1175 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic variables
1176 logical, intent(in) :: just_read !< If true, this call will only read
1177 !! parameters without changing h.
1178 real, dimension(SZI_(G),SZJ_(G)), &
1179 optional, intent(in) :: z_top_shelf !< Top interface position under ice shelf [Z ~> m]
1180 ! Local variables
1181 real, dimension(SZI_(G),SZJ_(G)) :: &
1182 eta_sfc ! The free surface height that the model should use [Z ~> m].
1183 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: &
1184 eta ! The free surface height that the model should use [Z ~> m].
1185 real :: dilate ! A ratio by which layers are dilated [nondim].
1186 real :: scale_factor ! A scaling factor for the eta_sfc values that are read in,
1187 ! which can be used to change units, for example, often [Z m-1 ~> 1].
1188 character(len=40) :: mdl = "depress_surface" ! This subroutine's name.
1189 character(len=200) :: inputdir, eta_srf_file ! Strings for file/path
1190 character(len=200) :: filename, eta_srf_var ! Strings for file/path
1191 integer :: i, j, k, is, ie, js, je, nz
1192 logical :: use_z_shelf
1193
1194 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
1195
1196 use_z_shelf = present(z_top_shelf)
1197
1198
1199 if (.not. use_z_shelf) then
1200 ! Read the surface height (or pressure) from a file.
1201
1202 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
1203 inputdir = slasher(inputdir)
1204 call get_param(param_file, mdl, "SURFACE_HEIGHT_IC_FILE", eta_srf_file, &
1205 "The initial condition file for the surface height.", &
1206 fail_if_missing=.not.just_read, do_not_log=just_read)
1207 call get_param(param_file, mdl, "SURFACE_HEIGHT_IC_VAR", eta_srf_var, &
1208 "The initial condition variable for the surface height.", &
1209 default="SSH", do_not_log=just_read)
1210 filename = trim(inputdir)//trim(eta_srf_file)
1211 if (.not.just_read) &
1212 call log_param(param_file, mdl, "INPUTDIR/SURFACE_HEIGHT_IC_FILE", filename)
1213
1214 call get_param(param_file, mdl, "SURFACE_HEIGHT_IC_SCALE", scale_factor, &
1215 "A scaling factor to convert SURFACE_HEIGHT_IC_VAR into units of m", &
1216 units="variable", default=1.0, scale=us%m_to_Z, do_not_log=just_read)
1217
1218 if (just_read) return ! All run-time parameters have been read, so return.
1219
1220 call mom_read_data(filename, eta_srf_var, eta_sfc, g%Domain, scale=scale_factor)
1221 else
1222 do j=js,je ; do i=is,ie
1223 eta_sfc(i,j) = z_top_shelf(i,j)
1224 enddo ; enddo
1225 endif
1226
1227 ! Convert thicknesses to interface heights.
1228 call find_eta(h, tv, g, gv, us, eta, dzref=g%Z_ref)
1229
1230 do j=js,je ; do i=is,ie ; if (g%mask2dT(i,j) > 0.0) then
1231! if (eta_sfc(i,j) < eta(i,j,nz+1)) then
1232 ! Issue a warning?
1233! endif
1234 if (eta_sfc(i,j) > eta(i,j,1)) then
1235 ! Dilate the water column to agree, but only up to 10-fold.
1236 if (eta_sfc(i,j) - eta(i,j,nz+1) > 10.0*(eta(i,j,1) - eta(i,j,nz+1))) then
1237 dilate = 10.0
1238 call mom_error(warning, "Free surface height dilation attempted "//&
1239 "to exceed 10-fold.", all_print=.true.)
1240 else
1241 dilate = (eta_sfc(i,j) - eta(i,j,nz+1)) / (eta(i,j,1) - eta(i,j,nz+1))
1242 endif
1243 do k=1,nz ; h(i,j,k) = h(i,j,k) * dilate ; enddo
1244 elseif (eta(i,j,1) > eta_sfc(i,j)) then
1245 ! Remove any mass that is above the target free surface.
1246 do k=1,nz
1247 if (eta(i,j,k) <= eta_sfc(i,j)) exit
1248 if (eta(i,j,k+1) >= eta_sfc(i,j)) then
1249 h(i,j,k) = gv%Angstrom_H
1250 else
1251 h(i,j,k) = max(gv%Angstrom_H, h(i,j,k) * &
1252 (eta_sfc(i,j) - eta(i,j,k+1)) / (eta(i,j,k) - eta(i,j,k+1)) )
1253 endif
1254 enddo
1255 endif
1256 endif ; enddo ; enddo
1257
1258end subroutine depress_surface
1259
1260!> Adjust the layer thicknesses by cutting away the top of each model column at the depth
1261!! where the hydrostatic pressure matches an imposed surface pressure read from file.
1262subroutine trim_for_ice(PF, G, GV, US, ALE_CSp, tv, h, just_read)
1263 type(param_file_type), intent(in) :: PF !< Parameter file structure
1264 type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1265 type(verticalgrid_type), intent(in) :: GV !< Vertical grid structure
1266 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1267 type(ale_cs), pointer :: ALE_CSp !< ALE control structure
1268 type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamics structure
1269 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1270 intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
1271 logical, intent(in) :: just_read !< If true, this call will only read
1272 !! parameters without changing h.
1273 ! Local variables
1274 character(len=200) :: mdl = "trim_for_ice"
1275 real, dimension(SZI_(G),SZJ_(G)) :: p_surf ! Imposed pressure on ocean at surface [R L2 T-2 ~> Pa]
1276 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: S_t, S_b ! Top and bottom edge values for reconstructions
1277 ! of salinity within each layer [S ~> ppt]
1278 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: T_t, T_b ! Top and bottom edge values for reconstructions
1279 ! of temperature within each layer [C ~> degC]
1280 character(len=200) :: inputdir, filename, p_surf_file, p_surf_var ! Strings for file/path
1281 real :: scale_factor ! A file-dependent scaling factor for the input pressure [various].
1282 real :: min_thickness ! The minimum layer thickness [H ~> m or kg m-2].
1283 real :: z_tolerance ! The tolerance with which to find the depth matching a specified pressure [Z ~> m].
1284 integer :: i, j, k
1285 integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.
1286 integer :: remap_answer_date ! The vintage of the order of arithmetic and expressions to use
1287 ! for remapping. Values below 20190101 recover the remapping
1288 ! answers from 2018, while higher values use more robust
1289 ! forms of the same remapping expressions.
1290 logical :: use_remapping ! If true, remap the initial conditions.
1291 logical :: use_frac_dp_bugfix ! If true, use bugfix. Otherwise, pressure input to EOS is negative.
1292 type(remapping_cs), pointer :: remap_CS => null()
1293
1294 call get_param(pf, mdl, "SURFACE_PRESSURE_FILE", p_surf_file, &
1295 "The initial condition file for the surface pressure exerted by ice.", &
1296 fail_if_missing=.not.just_read, do_not_log=just_read)
1297 call get_param(pf, mdl, "SURFACE_PRESSURE_VAR", p_surf_var, &
1298 "The initial condition variable for the surface pressure exerted by ice.", &
1299 default="", do_not_log=just_read)
1300 call get_param(pf, mdl, "INPUTDIR", inputdir, default=".", do_not_log=.true.)
1301 filename = trim(slasher(inputdir))//trim(p_surf_file)
1302 if (.not.just_read) call log_param(pf, mdl, "!INPUTDIR/SURFACE_HEIGHT_IC_FILE", filename)
1303
1304 call get_param(pf, mdl, "SURFACE_PRESSURE_SCALE", scale_factor, &
1305 "A scaling factor to convert SURFACE_PRESSURE_VAR from "//&
1306 "file SURFACE_PRESSURE_FILE into a surface pressure.", &
1307 units="file dependent", default=1., do_not_log=just_read)
1308 call get_param(pf, mdl, "MIN_THICKNESS", min_thickness, 'Minimum layer thickness', &
1309 units='m', default=1.e-3, scale=gv%m_to_H, do_not_log=just_read)
1310 call get_param(pf, mdl, "TRIM_IC_Z_TOLERANCE", z_tolerance, &
1311 "The tolerance with which to find the depth matching the specified "//&
1312 "surface pressure with TRIM_IC_FOR_P_SURF.", &
1313 units="m", default=1.0e-5, scale=us%m_to_Z, do_not_log=just_read)
1314 call get_param(pf, mdl, "FRAC_DP_AT_POS_NEGATIVE_P_BUGFIX", use_frac_dp_bugfix, &
1315 "If true, use bugfix in ice shelf TRIM_IC initialization. "//&
1316 "Otherwise, pressure input to density EOS is negative.", &
1317 default=.false., do_not_log=just_read)
1318 call get_param(pf, mdl, "TRIMMING_USES_REMAPPING", use_remapping, &
1319 'When trimming the column, also remap T and S.', &
1320 default=.false., do_not_log=just_read)
1321 if (use_remapping) then
1322 call get_param(pf, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
1323 "This sets the default value for the various _ANSWER_DATE parameters.", &
1324 default=99991231, do_not_log=just_read)
1325 call get_param(pf, mdl, "REMAPPING_ANSWER_DATE", remap_answer_date, &
1326 "The vintage of the expressions and order of arithmetic to use for remapping. "//&
1327 "Values below 20190101 result in the use of older, less accurate expressions "//&
1328 "that were in use at the end of 2018. Higher values result in the use of more "//&
1329 "robust and accurate forms of mathematically equivalent expressions.", &
1330 default=default_answer_date, do_not_log=just_read.or.(.not.gv%Boussinesq))
1331 if (.not.gv%Boussinesq) remap_answer_date = max(remap_answer_date, 20230701)
1332 else
1333 remap_answer_date = 20181231
1334 if (.not.gv%Boussinesq) remap_answer_date = 20230701
1335 endif
1336
1337 if (just_read) return ! All run-time parameters have been read, so return.
1338
1339 call mom_read_data(filename, p_surf_var, p_surf, g%Domain, &
1340 scale=scale_factor*us%Pa_to_RL2_T2)
1341
1342 if (use_remapping) then
1343 allocate(remap_cs)
1344 if (remap_answer_date < 20190101) then
1345 call initialize_remapping(remap_cs, 'PLM', boundary_extrapolation=.true., &
1346 h_neglect=1.0e-30*gv%m_to_H, h_neglect_edge=1.0e-10*gv%m_to_H)
1347 else
1348 call initialize_remapping(remap_cs, 'PLM', boundary_extrapolation=.true., &
1349 h_neglect=gv%H_subroundoff, h_neglect_edge=gv%H_subroundoff)
1350 endif
1351 endif
1352
1353 ! Find edge values of T and S used in reconstructions
1354 if ( associated(ale_csp) ) then ! This should only be associated if we are in ALE mode
1355 call ts_plm_edge_values(ale_csp, s_t, s_b, t_t, t_b, g, gv, tv, h, .true.)
1356 else
1357! call MOM_error(FATAL, "trim_for_ice: Does not work without ALE mode")
1358 do k=1,gv%ke ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
1359 t_t(i,j,k) = tv%T(i,j,k) ; t_b(i,j,k) = tv%T(i,j,k)
1360 s_t(i,j,k) = tv%S(i,j,k) ; s_b(i,j,k) = tv%S(i,j,k)
1361 enddo ; enddo ; enddo
1362 endif
1363
1364 do j=g%jsc,g%jec ; do i=g%isc,g%iec
1365 call cut_off_column_top(gv%ke, tv, gv, us, gv%g_Earth, g%bathyT(i,j)+g%Z_ref, min_thickness, &
1366 tv%T(i,j,:), t_t(i,j,:), t_b(i,j,:), tv%S(i,j,:), s_t(i,j,:), s_b(i,j,:), &
1367 p_surf(i,j), h(i,j,:), remap_cs, z_tol=z_tolerance, &
1368 frac_dp_bugfix=use_frac_dp_bugfix)
1369 enddo ; enddo
1370
1371end subroutine trim_for_ice
1372
1373!> Calculate the hydrostatic equilibrium position of the surface under an ice shelf
1374subroutine calc_sfc_displacement(PF, G, GV, US, mass_shelf, tv, h)
1375 type(param_file_type), intent(in) :: PF !< Parameter file structure
1376 type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1377 type(verticalgrid_type), intent(in) :: GV !< Vertical grid structure
1378 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1379 real, dimension(SZI_(G),SZJ_(G)), &
1380 intent(in) :: mass_shelf !< Ice shelf mass [R Z ~> kg m-2]
1381 type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamics structure
1382 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1383 intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
1384
1385 real :: z_top_shelf(SZI_(G),SZJ_(G)) ! The depth of the top interface under ice shelves [Z ~> m]
1386 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: &
1387 eta ! The free surface height that the model should use [Z ~> m].
1388 ! temporary arrays
1389 real, dimension(SZK_(GV)) :: rho_col ! potential density in the column for use in ice [R ~> kg m-3]
1390 real, dimension(SZK_(GV)) :: rho_h ! potential density multiplied by thickness [R Z ~> kg m-2]
1391 real, dimension(SZK_(GV)) :: h_tmp ! temporary storage for thicknesses [H ~> m]
1392 real, dimension(SZK_(GV)) :: p_ref ! pressure for density [R Z ~> kg m-2]
1393 real, dimension(SZK_(GV)+1) :: ei_tmp, ei_orig ! temporary storage for interface positions [Z ~> m]
1394 real :: z_top ! An estimate of the height of the ice-ocean interface [Z ~> m]
1395 real :: mass_disp ! The net mass of sea water that has been displaced by the shelf [R Z ~> kg m-2]
1396 real :: residual ! The difference between the displaced ocean mass and the ice shelf
1397 ! mass [R Z ~> kg m-2]
1398 real :: tol ! The initialization tolerance for ice shelf initialization [Z ~> m]
1399 integer :: is, ie, js, je, k, nz, i, j, max_iter, iter
1400
1401 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
1402
1403 call get_param(pf, mdl, "ICE_SHELF_INITIALIZATION_Z_TOLERANCE", tol, &
1404 "A initialization tolerance for the calculation of the static "// &
1405 "ice shelf displacement (m) using initial temperature and salinity profile.", &
1406 default=0.001, units="m", scale=us%m_to_Z)
1407 max_iter = 1e3
1408 call mom_mesg("Started calculating initial interface position under ice shelf ")
1409 ! Convert thicknesses to interface heights.
1410 call find_eta(h, tv, g, gv, us, eta, dzref=g%Z_ref)
1411 do j=js,je ; do i=is,ie
1412 iter = 1
1413 z_top_shelf(i,j) = 0.0
1414 p_ref(:) = tv%p_ref
1415 if ((g%mask2dT(i,j) > 0.) .and. (mass_shelf(i,j) > 0.)) then
1416 call calculate_density(tv%T(i,j,:), tv%S(i,j,:), p_ref, rho_col, tv%eqn_of_state)
1417 z_top = min(max(-1.0*mass_shelf(i,j)/rho_col(1), -g%bathyT(i,j)), 0.)
1418 h_tmp(:) = 0.0
1419 ei_tmp(1:nz+1) = eta(i,j,1:nz+1)
1420 ei_orig(1:nz+1) = eta(i,j,1:nz+1)
1421 do k=1,nz+1
1422 if (ei_tmp(k) < z_top) ei_tmp(k) = z_top
1423 enddo
1424 mass_disp = 0.0
1425 do k=1,nz
1426 h_tmp(k) = max(ei_tmp(k)-ei_tmp(k+1), gv%Angstrom_H)
1427 rho_h(k) = h_tmp(k) * rho_col(k)
1428 mass_disp = mass_disp + rho_h(k)
1429 enddo
1430 residual = mass_shelf(i,j) - mass_disp
1431 do while ((abs(residual) > tol) .and. (z_top > -g%bathyT(i,j)) .and. (iter < max_iter))
1432 z_top = min(max(z_top-(residual*0.5e-3), -g%bathyT(i,j)), 0.0)
1433 h_tmp(:) = 0.0
1434 ei_tmp(1:nz+1) = ei_orig(1:nz+1)
1435 do k=1,nz+1
1436 if (ei_tmp(k) < z_top) ei_tmp(k) = z_top
1437 enddo
1438 mass_disp = 0.0
1439 do k=1,nz
1440 h_tmp(k) = max(ei_tmp(k)-ei_tmp(k+1), gv%Angstrom_H)
1441 rho_h(k) = h_tmp(k) * rho_col(k)
1442 mass_disp = mass_disp + rho_h(k)
1443 enddo
1444 residual = mass_shelf(i,j) - mass_disp
1445 iter = iter+1
1446 enddo
1447 if (iter >= max_iter) call mom_mesg("Warning: calc_sfc_displacement too many iterations.")
1448 z_top_shelf(i,j) = z_top
1449 endif
1450 enddo ; enddo
1451 call mom_mesg("Calling depress_surface ")
1452 call depress_surface(h, g, gv, us, pf, tv, just_read=.false.,z_top_shelf=z_top_shelf)
1453 call mom_mesg("Finishing calling depress_surface ")
1454end subroutine calc_sfc_displacement
1455
1456!> Adjust the layer thicknesses by removing the top of the water column above the
1457!! depth where the hydrostatic pressure matches p_surf
1458subroutine cut_off_column_top(nk, tv, GV, US, G_earth, depth, min_thickness, T, T_t, T_b, &
1459 S, S_t, S_b, p_surf, h, remap_CS, z_tol, frac_dp_bugfix)
1460 integer, intent(in) :: nk !< Number of layers
1461 type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure
1462 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1463 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1464 real, intent(in) :: G_earth !< Gravitational acceleration [L2 Z-1 T-2 ~> m s-2]
1465 real, intent(in) :: depth !< Depth of ocean column [Z ~> m].
1466 real, intent(in) :: min_thickness !< Smallest thickness allowed [H ~> m or kg m-2].
1467 real, dimension(nk), intent(inout) :: T !< Layer mean temperature [C ~> degC]
1468 real, dimension(nk), intent(in) :: T_t !< Temperature at top of layer [C ~> degC]
1469 real, dimension(nk), intent(in) :: T_b !< Temperature at bottom of layer [C ~> degC]
1470 real, dimension(nk), intent(inout) :: S !< Layer mean salinity [S ~> ppt]
1471 real, dimension(nk), intent(in) :: S_t !< Salinity at top of layer [S ~> ppt]
1472 real, dimension(nk), intent(in) :: S_b !< Salinity at bottom of layer [S ~> ppt]
1473 real, intent(in) :: p_surf !< Imposed pressure on ocean at surface [R L2 T-2 ~> Pa]
1474 real, dimension(nk), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
1475 type(remapping_cs), pointer :: remap_CS !< Remapping structure for remapping T and S,
1476 !! if associated
1477 real, intent(in) :: z_tol !< The tolerance with which to find the depth
1478 !! matching the specified pressure [Z ~> m].
1479 logical, intent(in) :: frac_dp_bugfix !< If true, use bugfix in frac_dp_at_pos
1480
1481 ! Local variables
1482 real, dimension(nk+1) :: e ! Top and bottom edge positions for reconstructions [Z ~> m]
1483 real, dimension(nk) :: h0, h1 ! Initial and remapped layer thicknesses [H ~> m or kg m-2]
1484 real, dimension(nk) :: S0, S1 ! Initial and remapped layer salinities [S ~> ppt]
1485 real, dimension(nk) :: T0, T1 ! Initial and remapped layer temperatures [C ~> degC]
1486 real :: P_t, P_b ! Top and bottom pressures [R L2 T-2 ~> Pa]
1487 real :: z_out, e_top ! Interface height positions [Z ~> m]
1488 real :: min_dz ! The minimum thickness in depth units [Z ~> m]
1489 real :: dh_surf_rem ! The remaining thickness to remove in non-Bousinesq mode [H ~> kg m-2]
1490 integer :: k
1491
1492 ! Keep a copy of the initial thicknesses in reverse order to use in remapping
1493 do k=1,nk ; h0(k) = h(nk+1-k) ; enddo
1494
1495 if (gv%Boussinesq) then
1496 min_dz = gv%H_to_Z * min_thickness
1497 ! Calculate original interface positions
1498 e(nk+1) = -depth
1499 do k=nk,1,-1
1500 e(k) = e(k+1) + gv%H_to_Z*h(k)
1501 enddo
1502
1503 p_t = 0.
1504 e_top = e(1)
1505 do k=1,nk
1506 call find_depth_of_pressure_in_cell(t_t(k), t_b(k), s_t(k), s_b(k), e(k), e(k+1), &
1507 p_t, p_surf, gv%Rho0, g_earth, tv%eqn_of_state, &
1508 us, p_b, z_out, z_tol=z_tol, &
1509 frac_dp_bugfix=frac_dp_bugfix)
1510 if (z_out>=e(k)) then
1511 ! Imposed pressure was less that pressure at top of cell
1512 exit
1513 elseif (z_out<=e(k+1)) then
1514 ! Imposed pressure was greater than pressure at bottom of cell
1515 e_top = e(k+1)
1516 else
1517 ! Imposed pressure was fell between pressures at top and bottom of cell
1518 e_top = z_out
1519 exit
1520 endif
1521 p_t = p_b
1522 enddo
1523 if (e_top<e(1)) then
1524 ! Clip layers from the top down, if at all
1525 do k=1,nk
1526 if (e(k) > e_top) then
1527 ! Original e(K) is too high
1528 e(k) = e_top
1529 e_top = e_top - min_dz ! Next interface must be at least this deep
1530 endif
1531 ! This layer needs trimming
1532 h(k) = max( min_thickness, gv%Z_to_H * (e(k) - e(k+1)) )
1533 if (e(k) < e_top) exit ! No need to go further
1534 enddo
1535 endif
1536 else
1537 ! In non-Bousinesq mode, we are already in mass units so the calculation is much easier.
1538 if (p_surf > 0.0) then
1539 dh_surf_rem = p_surf * gv%RZ_to_H / g_earth
1540 do k=1,nk
1541 if (h(k) <= min_thickness) then ! This layer has no mass to remove.
1542 cycle
1543 elseif ((h(k) - min_thickness) < dh_surf_rem) then ! This layer should be removed entirely.
1544 dh_surf_rem = dh_surf_rem - (h(k) - min_thickness)
1545 h(k) = min_thickness
1546 else ! This is the last layer that should be removed.
1547 h(k) = h(k) - dh_surf_rem
1548 dh_surf_rem = 0.0
1549 exit
1550 endif
1551 enddo
1552 endif
1553 endif
1554
1555 ! Now we need to remap but remapping assumes the surface is at the
1556 ! same place in the two columns so we turn the column upside down.
1557 if (associated(remap_cs)) then
1558 do k=1,nk
1559 s0(k) = s(nk+1-k)
1560 t0(k) = t(nk+1-k)
1561 h1(k) = h(nk+1-k)
1562 enddo
1563 call remapping_core_h(remap_cs, nk, h0, t0, nk, h1, t1)
1564 call remapping_core_h(remap_cs, nk, h0, s0, nk, h1, s1)
1565 do k=1,nk
1566 s(k) = s1(nk+1-k)
1567 t(k) = t1(nk+1-k)
1568 enddo
1569 endif
1570
1571end subroutine cut_off_column_top
1572
1573!> Initialize horizontal velocity components from file
1574subroutine initialize_velocity_from_file(u, v, G, GV, US, param_file, just_read)
1575 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1576 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1577 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1578 intent(out) :: u !< The zonal velocity that is being initialized [L T-1 ~> m s-1]
1579 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1580 intent(out) :: v !< The meridional velocity that is being initialized [L T-1 ~> m s-1]
1581 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1582 type(param_file_type), intent(in) :: param_file !< A structure indicating the open file to
1583 !! parse for model parameter values.
1584 logical, intent(in) :: just_read !< If true, this call will only read
1585 !! parameters without changing u or v.
1586 ! Local variables
1587 character(len=40) :: mdl = "initialize_velocity_from_file" ! This subroutine's name.
1588 character(len=200) :: filename, velocity_file, inputdir ! Strings for file/path
1589 character(len=64) :: u_IC_var, v_IC_var ! Velocity component names in files
1590
1591 if (.not.just_read) call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
1592
1593 call get_param(param_file, mdl, "VELOCITY_FILE", velocity_file, &
1594 "The name of the velocity initial condition file.", &
1595 fail_if_missing=.not.just_read, do_not_log=just_read)
1596 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
1597 inputdir = slasher(inputdir)
1598
1599 filename = trim(velocity_file)
1600 if (scan(velocity_file, '/')== 0) then ! prepend inputdir if only a filename is given
1601 filename = trim(inputdir)//trim(velocity_file)
1602 endif
1603 if (.not.just_read) call log_param(param_file, mdl, "INPUTDIR/VELOCITY_FILE", filename)
1604
1605 call get_param(param_file, mdl, "U_IC_VAR", u_ic_var, &
1606 "The initial condition variable for zonal velocity in VELOCITY_FILE.", &
1607 default="u")
1608 call get_param(param_file, mdl, "V_IC_VAR", v_ic_var, &
1609 "The initial condition variable for meridional velocity in VELOCITY_FILE.", &
1610 default="v")
1611
1612 if (just_read) return ! All run-time parameters have been read, so return.
1613
1614 if (.not.file_exists(filename, g%Domain)) call mom_error(fatal, &
1615 " initialize_velocity_from_file: Unable to open "//trim(filename))
1616
1617 ! Read the velocities from a netcdf file.
1618 call mom_read_vector(filename, u_ic_var, v_ic_var, u(:,:,:), v(:,:,:), g%Domain, scale=us%m_s_to_L_T)
1619
1620 call calltree_leave(trim(mdl)//'()')
1621end subroutine initialize_velocity_from_file
1622
1623!> Initialize horizontal velocity components to zero.
1624subroutine initialize_velocity_zero(u, v, G, GV, param_file, just_read)
1625 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1626 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1627 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1628 intent(out) :: u !< The zonal velocity that is being initialized [L T-1 ~> m s-1]
1629 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1630 intent(out) :: v !< The meridional velocity that is being initialized [L T-1 ~> m s-1]
1631 type(param_file_type), intent(in) :: param_file !< A structure indicating the open file to
1632 !! parse for model parameter values.
1633 logical, intent(in) :: just_read !< If true, this call will only read
1634 !! parameters without changing h.
1635 ! Local variables
1636 character(len=200) :: mdl = "initialize_velocity_zero" ! This subroutine's name.
1637 integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
1638
1639 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
1640 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1641
1642 if (.not.just_read) call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
1643
1644 if (just_read) return ! All run-time parameters have been read, so return.
1645
1646 do k=1,nz ; do j=js,je ; do i=isq,ieq
1647 u(i,j,k) = 0.0
1648 enddo ; enddo ; enddo
1649 do k=1,nz ; do j=jsq,jeq ; do i=is,ie
1650 v(i,j,k) = 0.0
1651 enddo ; enddo ; enddo
1652
1653 call calltree_leave(trim(mdl)//'()')
1654end subroutine initialize_velocity_zero
1655
1656!> Sets the initial velocity components to uniform
1657subroutine initialize_velocity_uniform(u, v, G, GV, US, param_file, just_read)
1658 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1659 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1660 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1661 intent(out) :: u !< The zonal velocity that is being initialized [L T-1 ~> m s-1]
1662 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1663 intent(out) :: v !< The meridional velocity that is being initialized [L T-1 ~> m s-1]
1664 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1665 type(param_file_type), intent(in) :: param_file !< A structure indicating the open file to
1666 !! parse for model parameter values.
1667 logical, intent(in) :: just_read !< If true, this call will only read
1668 !! parameters without changing u or v.
1669 ! Local variables
1670 integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
1671 real :: initial_u_const, initial_v_const ! Constant initial velocities [L T-1 ~> m s-1]
1672 character(len=200) :: mdl = "initialize_velocity_uniform" ! This subroutine's name.
1673
1674 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
1675 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1676
1677 call get_param(param_file, mdl, "INITIAL_U_CONST", initial_u_const, &
1678 "A initial uniform value for the zonal flow.", &
1679 default=0.0, units="m s-1", scale=us%m_s_to_L_T, do_not_log=just_read)
1680 call get_param(param_file, mdl, "INITIAL_V_CONST", initial_v_const, &
1681 "A initial uniform value for the meridional flow.", &
1682 default=0.0, units="m s-1", scale=us%m_s_to_L_T, do_not_log=just_read)
1683
1684 if (just_read) return ! All run-time parameters have been read, so return.
1685
1686 do k=1,nz ; do j=js,je ; do i=isq,ieq
1687 u(i,j,k) = initial_u_const
1688 enddo ; enddo ; enddo
1689 do k=1,nz ; do j=jsq,jeq ; do i=is,ie
1690 v(i,j,k) = initial_v_const
1691 enddo ; enddo ; enddo
1692
1693end subroutine initialize_velocity_uniform
1694
1695!> Sets the initial velocity components to be circular with
1696!! no flow at edges of domain and center.
1697subroutine initialize_velocity_circular(u, v, G, GV, US, param_file, just_read)
1698 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1699 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1700 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1701 intent(out) :: u !< The zonal velocity that is being initialized [L T-1 ~> m s-1]
1702 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1703 intent(out) :: v !< The meridional velocity that is being initialized [L T-1 ~> m s-1]
1704 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1705 type(param_file_type), intent(in) :: param_file !< A structure indicating the open file to
1706 !! parse for model parameter values.
1707 logical, intent(in) :: just_read !< If true, this call will only read
1708 !! parameters without changing u or v.
1709 ! Local variables
1710 character(len=200) :: mdl = "initialize_velocity_circular"
1711 real :: circular_max_u ! The amplitude of the zonal flow [L T-1 ~> m s-1]
1712 real :: dpi ! A local variable storing pi = 3.14159265358979... [nondim]
1713 real :: psi1, psi2 ! Values of the streamfunction at two points [L2 T-1 ~> m2 s-1]
1714 integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
1715 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
1716 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
1717
1718 call get_param(param_file, mdl, "CIRCULAR_MAX_U", circular_max_u, &
1719 "The amplitude of zonal flow from which to scale the "// &
1720 "circular stream function [m s-1].", &
1721 units="m s-1", default=0., scale=us%m_s_to_L_T, do_not_log=just_read)
1722
1723 if (just_read) return ! All run-time parameters have been read, so return.
1724
1725 if (g%grid_unit_to_L <= 0.) call mom_error(fatal, "MOM_state_initialization.F90: "//&
1726 "initialize_velocity_circular() is only set to work with Cartesian axis units.")
1727
1728 dpi = acos(0.0)*2.0 ! pi
1729
1730 do k=1,nz ; do j=js,je ; do i=isq,ieq
1731 psi1 = my_psi(i,j)
1732 psi2 = my_psi(i,j-1)
1733 u(i,j,k) = (psi1 - psi2) / g%dy_Cu(i,j) ! *(circular_max_u*G%len_lon/(2.0*dpi))
1734 enddo ; enddo ; enddo
1735 do k=1,nz ; do j=jsq,jeq ; do i=is,ie
1736 psi1 = my_psi(i,j)
1737 psi2 = my_psi(i-1,j)
1738 v(i,j,k) = (psi2 - psi1) / g%dx_Cv(i,j) ! *(circular_max_u*G%len_lon/(2.0*dpi))
1739 enddo ; enddo ; enddo
1740
1741 contains
1742
1743 !> Returns the value of a circular stream function at (ig,jg) in [L2 T-1 ~> m2 s-1]
1744 real function my_psi(ig,jg)
1745 integer :: ig !< Global i-index
1746 integer :: jg !< Global j-index
1747 ! Local variables
1748 real :: x, y, r ! [nondim]
1749
1750 x = 2.0*(g%geoLonBu(ig,jg)-g%west_lon) / g%len_lon - 1.0 ! -1<x<1
1751 y = 2.0*(g%geoLatBu(ig,jg)-g%south_lat) / g%len_lat - 1.0 ! -1<y<1
1752 r = sqrt( (x**2) + (y**2) ) ! Circular stream function is a function of radius only
1753 r = min(1.0, r) ! Flatten stream function in corners of box
1754 my_psi = 0.5*(1.0 - cos(dpi*r))
1755 my_psi = my_psi * (circular_max_u * g%len_lon * g%grid_unit_to_L / dpi) ! len_lon is in km
1756 end function my_psi
1757
1758end subroutine initialize_velocity_circular
1759
1760!> Initializes temperature and salinity from file
1761subroutine initialize_temp_salt_from_file(T, S, G, GV, US, param_file, just_read)
1762 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1763 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1764 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T !< The potential temperature that is
1765 !! being initialized [C ~> degC]
1766 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S !< The salinity that is
1767 !! being initialized [S ~> ppt]
1768 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1769 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
1770 logical, intent(in) :: just_read !< If true, this call will only
1771 !! read parameters without changing T or S.
1772 ! Local variables
1773 character(len=200) :: filename, salt_filename ! Full paths to input files
1774 character(len=200) :: ts_file, salt_file, inputdir ! Strings for file/path
1775 character(len=40) :: mdl = "initialize_temp_salt_from_file"
1776 character(len=64) :: temp_var, salt_var ! Temperature and salinity names in files
1777
1778 if (.not.just_read) call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
1779
1780 call get_param(param_file, mdl, "TS_FILE", ts_file, &
1781 "The initial condition file for temperature.", &
1782 fail_if_missing=.not.just_read, do_not_log=just_read)
1783 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
1784 inputdir = slasher(inputdir)
1785
1786 filename = trim(ts_file)
1787 if (scan(ts_file, '/')== 0) then ! prepend inputdir if only a filename is given
1788 filename = trim(inputdir)//trim(ts_file)
1789 endif
1790 if (.not.just_read) call log_param(param_file, mdl, "INPUTDIR/TS_FILE", filename)
1791 call get_param(param_file, mdl, "TEMP_IC_VAR", temp_var, &
1792 "The initial condition variable for potential temperature.", &
1793 default="PTEMP", do_not_log=just_read)
1794 call get_param(param_file, mdl, "SALT_IC_VAR", salt_var, &
1795 "The initial condition variable for salinity.", &
1796 default="SALT", do_not_log=just_read)
1797 call get_param(param_file, mdl, "SALT_FILE", salt_file, &
1798 "The initial condition file for salinity.", &
1799 default=trim(ts_file), do_not_log=just_read)
1800
1801 if (just_read) return ! All run-time parameters have been read, so return.
1802
1803 if (.not.file_exists(filename, g%Domain)) call mom_error(fatal, &
1804 " initialize_temp_salt_from_file: Unable to open "//trim(filename))
1805
1806 ! Read the temperatures and salinities from netcdf files.
1807 call mom_read_data(filename, temp_var, t(:,:,:), g%Domain, scale=us%degC_to_C)
1808
1809 salt_filename = trim(salt_file)
1810 if (scan(salt_file, '/')== 0) then ! prepend inputdir if only a filename is given
1811 salt_filename = trim(inputdir)//trim(salt_file)
1812 endif
1813 if (.not.file_exists(salt_filename, g%Domain)) call mom_error(fatal, &
1814 " initialize_temp_salt_from_file: Unable to open "//trim(salt_filename))
1815
1816 call mom_read_data(salt_filename, salt_var, s(:,:,:), g%Domain, scale=us%ppt_to_S)
1817
1818 call calltree_leave(trim(mdl)//'()')
1819end subroutine initialize_temp_salt_from_file
1820
1821!> Initializes temperature and salinity from a 1D profile
1822subroutine initialize_temp_salt_from_profile(T, S, G, GV, US, param_file, just_read)
1823 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1824 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1825 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T !< The potential temperature that is
1826 !! being initialized [C ~> degC]
1827 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S !< The salinity that is
1828 !! being initialized [S ~> ppt]
1829 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1830 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
1831 logical, intent(in) :: just_read !< If true, this call will only read
1832 !! parameters without changing T or S.
1833 ! Local variables
1834 real, dimension(SZK_(GV)) :: T0 ! The profile of temperatures [C ~> degC]
1835 real, dimension(SZK_(GV)) :: S0 ! The profile of salinities [S ~> ppt]
1836 integer :: i, j, k
1837 character(len=200) :: filename, ts_file, inputdir ! Strings for file/path
1838 character(len=64) :: temp_var, salt_var ! Temperature and salinity names in files
1839 character(len=40) :: mdl = "initialize_temp_salt_from_profile"
1840
1841 if (.not.just_read) call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
1842
1843 call get_param(param_file, mdl, "TS_FILE", ts_file, &
1844 "The file with the reference profiles for temperature and salinity.", &
1845 fail_if_missing=.not.just_read, do_not_log=just_read)
1846 call get_param(param_file, mdl, "TEMP_IC_VAR", temp_var, &
1847 "The initial condition variable for potential temperature.", &
1848 default="PTEMP", do_not_log=just_read)
1849 call get_param(param_file, mdl, "SALT_IC_VAR", salt_var, &
1850 "The initial condition variable for salinity.", &
1851 default="SALT", do_not_log=just_read)
1852
1853 if (just_read) return ! All run-time parameters have been read, so return.
1854
1855 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
1856 inputdir = slasher(inputdir)
1857 filename = trim(inputdir)//trim(ts_file)
1858 call log_param(param_file, mdl, "INPUTDIR/TS_FILE", filename)
1859 if (.not.file_exists(filename)) call mom_error(fatal, &
1860 " initialize_temp_salt_from_profile: Unable to open "//trim(filename))
1861
1862 ! Read the temperatures and salinities from a netcdf file.
1863 call mom_read_data(filename, temp_var, t0(:), scale=us%degC_to_C)
1864 call mom_read_data(filename, salt_var, s0(:), scale=us%ppt_to_S)
1865
1866 do k=1,gv%ke ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
1867 t(i,j,k) = t0(k) ; s(i,j,k) = s0(k)
1868 enddo ; enddo ; enddo
1869
1870 call calltree_leave(trim(mdl)//'()')
1872
1873!> Initializes temperature and salinity by fitting to density
1874subroutine initialize_temp_salt_fit(T, S, G, GV, US, param_file, eqn_of_state, P_Ref, just_read)
1875 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1876 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1877 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T !< The potential temperature that is
1878 !! being initialized [C ~> degC].
1879 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S !< The salinity that is being
1880 !! initialized [S ~> ppt].
1881 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1882 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
1883 !! parameters.
1884 type(eos_type), intent(in) :: eqn_of_state !< Equation of state structure
1885 real, intent(in) :: P_Ref !< The coordinate-density reference pressure
1886 !! [R L2 T-2 ~> Pa].
1887 logical, intent(in) :: just_read !< If true, this call will only read
1888 !! parameters without changing T or S.
1889 ! Local variables
1890 real :: T0(SZK_(GV)) ! Layer potential temperatures [C ~> degC]
1891 real :: S0(SZK_(GV)) ! Layer salinities [S ~> ppt]
1892 real :: T_Ref ! Reference Temperature [C ~> degC]
1893 real :: S_Ref ! Reference Salinity [S ~> ppt]
1894 real :: pres(SZK_(GV)) ! An array of the reference pressure [R L2 T-2 ~> Pa].
1895 real :: drho_dT(SZK_(GV)) ! Derivative of density with temperature [R C-1 ~> kg m-3 degC-1].
1896 real :: drho_dS(SZK_(GV)) ! Derivative of density with salinity [R S-1 ~> kg m-3 ppt-1].
1897 real :: rho_guess(SZK_(GV)) ! Potential density at T0 & S0 [R ~> kg m-3].
1898 logical :: fit_salin ! If true, accept the prescribed temperature and fit the salinity.
1899 character(len=40) :: mdl = "initialize_temp_salt_fit" ! This subroutine's name.
1900 integer :: i, j, k, itt, nz
1901 nz = gv%ke
1902
1903 if (.not.just_read) call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
1904
1905 call get_param(param_file, mdl, "T_REF", t_ref, &
1906 "A reference temperature used in initialization.", &
1907 units="degC", scale=us%degC_to_C, fail_if_missing=.not.just_read, do_not_log=just_read)
1908 call get_param(param_file, mdl, "S_REF", s_ref, &
1909 "A reference salinity used in initialization.", &
1910 units="ppt", default=35.0, scale=us%ppt_to_S, do_not_log=just_read)
1911 call get_param(param_file, mdl, "FIT_SALINITY", fit_salin, &
1912 "If true, accept the prescribed temperature and fit the "//&
1913 "salinity; otherwise take salinity and fit temperature.", &
1914 default=.false., do_not_log=just_read)
1915
1916 if (just_read) return ! All run-time parameters have been read, so return.
1917
1918 do k=1,nz
1919 pres(k) = p_ref ; s0(k) = s_ref
1920 t0(k) = t_ref
1921 enddo
1922
1923 call calculate_density(t0(1), s0(1), pres(1), rho_guess(1), eqn_of_state)
1924 call calculate_density_derivs(t0, s0, pres, drho_dt, drho_ds, eqn_of_state, (/1,1/) )
1925
1926 if (fit_salin) then
1927 ! A first guess of the layers' temperatures.
1928 do k=nz,1,-1
1929 s0(k) = max(0.0, s0(1) + (gv%Rlay(k) - rho_guess(1)) / drho_ds(1))
1930 enddo
1931 ! Refine the guesses for each layer.
1932 do itt=1,6
1933 call calculate_density(t0, s0, pres, rho_guess, eqn_of_state)
1934 call calculate_density_derivs(t0, s0, pres, drho_dt, drho_ds, eqn_of_state)
1935 do k=1,nz
1936 s0(k) = max(0.0, s0(k) + (gv%Rlay(k) - rho_guess(k)) / drho_ds(k))
1937 enddo
1938 enddo
1939 else
1940 ! A first guess of the layers' temperatures.
1941 do k=nz,1,-1
1942 t0(k) = t0(1) + (gv%Rlay(k) - rho_guess(1)) / drho_dt(1)
1943 enddo
1944 do itt=1,6
1945 call calculate_density(t0, s0, pres, rho_guess, eqn_of_state)
1946 call calculate_density_derivs(t0, s0, pres, drho_dt, drho_ds, eqn_of_state)
1947 do k=1,nz
1948 t0(k) = t0(k) + (gv%Rlay(k) - rho_guess(k)) / drho_dt(k)
1949 enddo
1950 enddo
1951 endif
1952
1953 do k=1,nz ; do j=g%jsd,g%jed ; do i=g%isd,g%ied
1954 t(i,j,k) = t0(k) ; s(i,j,k) = s0(k)
1955 enddo ; enddo ; enddo
1956
1957 call calltree_leave(trim(mdl)//'()')
1958end subroutine initialize_temp_salt_fit
1959
1960!> Initializes T and S with linear profiles according to reference surface
1961!! layer salinity and temperature and a specified range.
1962!!
1963!! \remark Note that the linear distribution is set up with respect to the layer
1964!! number, not the physical position).
1965subroutine initialize_temp_salt_linear(T, S, G, GV, US, param_file, just_read)
1966 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1967 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1968 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T !< The potential temperature that is
1969 !! being initialized [C ~> degC]
1970 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S !< The salinity that is
1971 !! being initialized [S ~> ppt]
1972 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1973 type(param_file_type), intent(in) :: param_file !< A structure to parse for
1974 !! run-time parameters
1975 logical, intent(in) :: just_read !< If present and true,
1976 !! this call will only read parameters
1977 !! without changing T or S.
1978
1979 ! Local variables
1980 real :: S_top, S_range ! Reference salinity in the surface layer and its vertical range [S ~> ppt]
1981 real :: T_top, T_range ! Reference temperature in the surface layer and its vertical range [C ~> degC]
1982 character(len=40) :: mdl = "initialize_temp_salt_linear" ! This subroutine's name.
1983 integer :: k
1984
1985 if (.not.just_read) call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
1986 call get_param(param_file, mdl, "T_TOP", t_top, &
1987 "Initial temperature of the top surface.", &
1988 units="degC", scale=us%degC_to_C, fail_if_missing=.not.just_read, do_not_log=just_read)
1989 call get_param(param_file, mdl, "T_RANGE", t_range, &
1990 "Initial temperature difference (top-bottom).", &
1991 units="degC", scale=us%degC_to_C, fail_if_missing=.not.just_read, do_not_log=just_read)
1992 call get_param(param_file, mdl, "S_TOP", s_top, &
1993 "Initial salinity of the top surface.", &
1994 units="ppt", scale=us%ppt_to_S, fail_if_missing=.not.just_read, do_not_log=just_read)
1995 call get_param(param_file, mdl, "S_RANGE", s_range, &
1996 "Initial salinity difference (top-bottom).", &
1997 units="ppt", scale=us%ppt_to_S, fail_if_missing=.not.just_read, do_not_log=just_read)
1998
1999 if (just_read) return ! All run-time parameters have been read, so return.
2000
2001 ! Prescribe salinity and temperature, with the extrapolated top interface value prescribed.
2002 do k=1,gv%ke
2003 s(:,:,k) = s_top - s_range*((real(k)-0.5)/real(gv%ke))
2004 t(:,:,k) = t_top - t_range*((real(k)-0.5)/real(gv%ke))
2005 enddo
2006
2007 ! Prescribe salinity and temperature, but with the top layer value matching the surface value.
2008 ! S(:,:,1) = S_top ; T(:,:,1) = T_top
2009 ! do k=2,GV%ke
2010 ! S(:,:,k) = S_top - S_range * (real(k-1) / real(GV%ke-1))
2011 ! T(:,:,k) = T_top - T_range * (real(k-1) / real(GV%ke-1))
2012 ! enddo
2013
2014 call calltree_leave(trim(mdl)//'()')
2015end subroutine initialize_temp_salt_linear
2016
2017!> This subroutine sets the inverse restoration time (Idamp), and
2018!! the values towards which the interface heights and an arbitrary
2019!! number of tracers should be restored within each sponge. The
2020!! interface height is always subject to damping, and must always be
2021!! the first registered field.
2022subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, u, v, depth_tot, param_file, &
2023 Layer_CSp, ALE_CSp, Time)
2024 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
2025 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
2026 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2027 logical, intent(in) :: use_temperature !< If true, T & S are state variables.
2028 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic
2029 !! variables.
2030 real, target, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
2031 intent(in) :: u !< The zonal velocity that is being
2032 !! initialized [L T-1 ~> m s-1]
2033 real, target, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
2034 intent(in) :: v !< The meridional velocity that is being
2035 !! initialized [L T-1 ~> m s-1]
2036 real, dimension(SZI_(G),SZJ_(G)), &
2037 intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
2038 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters.
2039 type(sponge_cs), pointer :: Layer_CSp !< A pointer that is set to point to the control
2040 !! structure for this module (in layered mode).
2041 type(ale_sponge_cs), pointer :: ALE_CSp !< A pointer that is set to point to the control
2042 !! structure for this module (in ALE mode).
2043 type(time_type), intent(in) :: Time !< Time at the start of the run segment. Time_in
2044 !! overrides any value set for Time.
2045 ! Local variables
2046 real, allocatable, dimension(:,:,:) :: eta ! The target interface heights [Z ~> m].
2047 real, allocatable, dimension(:,:,:) :: dz ! The target interface thicknesses in height units [Z ~> m]
2048
2049 real, dimension (SZI_(G),SZJ_(G),SZK_(GV)) :: &
2050 tmp, & ! A temporary array for temperatures [C ~> degC] or other tracers.
2051 tmp2 ! A temporary array for salinities [S ~> ppt]
2052 real, dimension (SZI_(G),SZJ_(G)) :: &
2053 tmp_2d ! A temporary array for mixed layer densities [R ~> kg m-3]
2054 real, allocatable, dimension(:,:,:) :: tmp_T ! A temporary array for reading sponge target temperatures
2055 ! on the vertical grid of the input file [C ~> degC]
2056 real, allocatable, dimension(:,:,:) :: tmp_S ! A temporary array for reading sponge target salinities
2057 ! on the vertical grid of the input file [S ~> ppt]
2058 real, allocatable, dimension(:,:,:) :: tmp_u ! Temporary array for reading sponge target zonal
2059 ! velocities on the vertical grid of the input file [L T-1 ~> m s-1]
2060 real, allocatable, dimension(:,:,:) :: tmp_v ! Temporary array for reading sponge target meridional
2061 ! velocities on the vertical grid of the input file [L T-1 ~> m s-1]
2062
2063 real :: Idamp(SZI_(G),SZJ_(G)) ! The sponge damping rate [T-1 ~> s-1]
2064 real :: Idamp_u(SZIB_(G),SZJ_(G)) ! The sponge damping rate for velocity fields [T-1 ~> s-1]
2065 real :: Idamp_v(SZI_(G),SZJB_(G)) ! The sponge damping rate for velocity fields [T-1 ~> s-1]
2066 real :: pres(SZI_(G)) ! An array of the reference pressure [R L2 T-2 ~> Pa]
2067
2068 integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
2069 integer :: i, j, k, is, ie, js, je, nz
2070 integer :: isd, ied, jsd, jed
2071 integer, dimension(4) :: siz
2072 integer :: nz_data ! The size of the sponge source grid
2073 logical :: sponge_uv ! Apply sponges in u and v, in addition to tracers.
2074 character(len=40) :: potemp_var, salin_var, u_var, v_var, Idamp_var, Idamp_u_var, Idamp_v_var, eta_var
2075 character(len=40) :: mdl = "initialize_sponges_file"
2076 character(len=200) :: damping_file, uv_damping_file, state_file, state_uv_file ! Strings for filenames
2077 character(len=200) :: filename, inputdir ! Strings for file/path and path.
2078
2079 logical :: use_ALE ! True if ALE is being used, False if in layered mode
2080 logical :: time_space_interp_sponge ! If true use sponge data that need to be interpolated in both
2081 ! the horizontal dimension and in time prior to vertical remapping.
2082
2083 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
2084 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
2085
2086 pres(:) = 0.0 ; tmp(:,:,:) = 0.0 ; idamp(:,:) = 0.0 ; idamp_u(:,:) = 0.0 ; idamp_v(:,:) = 0.0
2087
2088 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
2089 inputdir = slasher(inputdir)
2090 call get_param(param_file, mdl, "SPONGE_DAMPING_FILE", damping_file, &
2091 "The name of the file with the sponge damping rates.", &
2092 fail_if_missing=.true.)
2093 call get_param(param_file, mdl, "SPONGE_STATE_FILE", state_file, &
2094 "The name of the file with the state to damp toward.", &
2095 default=damping_file)
2096 call get_param(param_file, mdl, "SPONGE_PTEMP_VAR", potemp_var, &
2097 "The name of the potential temperature variable in "//&
2098 "SPONGE_STATE_FILE.", default="PTEMP")
2099 call get_param(param_file, mdl, "SPONGE_SALT_VAR", salin_var, &
2100 "The name of the salinity variable in "//&
2101 "SPONGE_STATE_FILE.", default="SALT")
2102 call get_param(param_file, mdl, "SPONGE_UV", sponge_uv, &
2103 "Apply sponges in u and v, in addition to tracers.", &
2104 default=.false.)
2105 if (sponge_uv) then
2106 call get_param(param_file, mdl, "SPONGE_UV_STATE_FILE", state_uv_file, &
2107 "The name of the file with the state to damp UV toward.", &
2108 default=damping_file)
2109 call get_param(param_file, mdl, "SPONGE_U_VAR", u_var, &
2110 "The name of the zonal velocity variable in "//&
2111 "SPONGE_UV_STATE_FILE.", default="UVEL")
2112 call get_param(param_file, mdl, "SPONGE_V_VAR", v_var, &
2113 "The name of the vertical velocity variable in "//&
2114 "SPONGE_UV_STATE_FILE.", default="VVEL")
2115 endif
2116 call get_param(param_file, mdl, "SPONGE_ETA_VAR", eta_var, &
2117 "The name of the interface height variable in "//&
2118 "SPONGE_STATE_FILE.", default="ETA")
2119 call get_param(param_file, mdl, "SPONGE_IDAMP_VAR", idamp_var, &
2120 "The name of the inverse damping rate variable in "//&
2121 "SPONGE_DAMPING_FILE.", default="Idamp")
2122 if (sponge_uv) then
2123 call get_param(param_file, mdl, "SPONGE_UV_DAMPING_FILE", uv_damping_file, &
2124 "The name of the file with sponge damping rates for the velocity variables.", &
2125 default=damping_file)
2126 call get_param(param_file, mdl, "SPONGE_IDAMP_U_var", idamp_u_var, &
2127 "The name of the inverse damping rate variable in "//&
2128 "SPONGE_UV_DAMPING_FILE for the velocities.", default=idamp_var)
2129 call get_param(param_file, mdl, "SPONGE_IDAMP_V_var", idamp_v_var, &
2130 "The name of the inverse damping rate variable in "//&
2131 "SPONGE_UV_DAMPING_FILE for the velocities.", default=idamp_var)
2132 endif
2133 call get_param(param_file, mdl, "USE_REGRIDDING", use_ale, default=.false., do_not_log=.true.)
2134
2135 call get_param(param_file, mdl, "INTERPOLATE_SPONGE_TIME_SPACE", time_space_interp_sponge, &
2136 "If True, perform on-the-fly regridding in lat-lon-time of sponge restoring data.", &
2137 default=.false.)
2138
2139 ! Read in sponge damping rate for tracers
2140 filename = trim(damping_file)
2141 if (scan(damping_file, '/')== 0) then ! prepend inputdir if only a filename is given
2142 filename = trim(inputdir)//trim(damping_file)
2143 endif
2144 call log_param(param_file, mdl, "INPUTDIR/SPONGE_DAMPING_FILE", filename)
2145 if (.not.file_exists(filename, g%Domain)) &
2146 call mom_error(fatal, " initialize_sponges: Unable to open "//trim(filename))
2147
2148 if (time_space_interp_sponge .and. .not. use_ale) &
2149 call mom_error(fatal, " initialize_sponges: Time-varying sponges are currently unavailable in layered mode ")
2150
2151 call mom_read_data(filename, idamp_var, idamp(:,:), g%Domain, scale=us%T_to_s)
2152
2153 ! Read in sponge damping rate for velocities
2154 if (sponge_uv) then
2155 if (separate_idamp_for_uv()) then
2156 filename = trim(inputdir)//trim(uv_damping_file)
2157 call log_param(param_file, mdl, "INPUTDIR/SPONGE_UV_DAMPING_FILE", filename)
2158
2159 if (.not.file_exists(filename, g%Domain)) &
2160 call mom_error(fatal, " initialize_sponges: Unable to open "//trim(filename))
2161
2162 call mom_read_vector(filename, idamp_u_var,idamp_v_var,idamp_u(:,:),idamp_v(:,:), g%Domain, scale=us%T_to_s)
2163 else
2164 ! call MOM_error(FATAL, "Must provide SPONGE_IDAMP_U_var and SPONGE_IDAMP_V_var")
2165 call pass_var(idamp,g%Domain)
2166 do j=g%jsc,g%jec
2167 do i=g%iscB,g%iecB
2168 idamp_u(i,j) = 0.5*(idamp(i,j)+idamp(i+1,j))
2169 enddo
2170 enddo
2171 do j=g%jscB,g%jecB
2172 do i=g%isc,g%iec
2173 idamp_v(i,j) = 0.5*(idamp(i,j)+idamp(i,j+1))
2174 enddo
2175 enddo
2176 endif
2177 endif
2178
2179 ! Now register all of the fields which are damped in the sponge.
2180 ! By default, momentum is advected vertically within the sponge, but
2181 ! momentum is typically not damped within the sponge.
2182
2183 filename = trim(inputdir)//trim(state_file)
2184 call log_param(param_file, mdl, "INPUTDIR/SPONGE_STATE_FILE", filename)
2185 if (.not.file_exists(filename, g%Domain)) &
2186 call mom_error(fatal, " initialize_sponges: Unable to open "//trim(filename))
2187
2188
2189 if (.not. use_ale) then
2190 ! The first call to set_up_sponge_field is for the interface heights if in layered mode.
2191 allocate(eta(isd:ied,jsd:jed,nz+1), source=0.0)
2192 call mom_read_data(filename, eta_var, eta(:,:,:), g%Domain, scale=us%m_to_Z)
2193
2194 do j=js,je ; do i=is,ie
2195 eta(i,j,nz+1) = -depth_tot(i,j)
2196 enddo ; enddo
2197 do k=nz,1,-1 ; do j=js,je ; do i=is,ie
2198 if (eta(i,j,k) < (eta(i,j,k+1) + gv%Angstrom_Z)) &
2199 eta(i,j,k) = eta(i,j,k+1) + gv%Angstrom_Z
2200 enddo ; enddo ; enddo
2201 ! Set the sponge damping rates so that the model will know where to
2202 ! apply the sponges, along with the interface heights.
2203 call initialize_sponge(idamp, eta, g, param_file, layer_csp, gv)
2204 deallocate(eta)
2205
2206 if ( gv%nkml>0) then
2207 ! This call to set_up_sponge_ML_density registers the target values of the
2208 ! mixed layer density, which is used in determining which layers can be
2209 ! inflated without causing static instabilities.
2210 do i=is,ie ; pres(i) = tv%P_Ref ; enddo
2211 eosdom(:) = eos_domain(g%HI)
2212
2213 call mom_read_data(filename, potemp_var, tmp(:,:,:), g%Domain, scale=us%degC_to_C)
2214 call mom_read_data(filename, salin_var, tmp2(:,:,:), g%Domain, scale=us%ppt_to_S)
2215
2216 do j=js,je
2217 call calculate_density(tmp(:,j,1), tmp2(:,j,1), pres, tmp_2d(:,j), tv%eqn_of_state, eosdom)
2218 enddo
2219
2220 call set_up_sponge_ml_density(tmp_2d, g, layer_csp)
2221 endif
2222
2223 ! Now register all of the tracer fields which are damped in the
2224 ! sponge. By default, momentum is advected vertically within the
2225 ! sponge, but momentum is typically not damped within the sponge.
2226
2227
2228 ! The remaining calls to set_up_sponge_field can be in any order.
2229 if ( use_temperature) then
2230 call mom_read_data(filename, potemp_var, tmp(:,:,:), g%Domain, scale=us%degC_to_C)
2231 call set_up_sponge_field(tmp, tv%T, g, gv, nz, layer_csp)
2232 call mom_read_data(filename, salin_var, tmp2(:,:,:), g%Domain, scale=us%ppt_to_S)
2233 call set_up_sponge_field(tmp2, tv%S, g, gv, nz, layer_csp)
2234 endif
2235
2236! else
2237 ! Initialize sponges without supplying sponge grid
2238! if (sponge_uv) then
2239! call initialize_ALE_sponge(Idamp, G, GV, US, param_file, ALE_CSp, Idamp_u, Idamp_v)
2240! else
2241! call initialize_ALE_sponge(Idamp, G, GV, US, param_file, ALE_CSp)
2242! endif
2243 endif
2244
2245
2246 if (use_ale) then ! ALE mode
2247 if (.not. time_space_interp_sponge) then
2248 call field_size(filename,eta_var,siz,no_domain=.true.)
2249 if (siz(1) /= g%ieg-g%isg+1 .or. siz(2) /= g%jeg-g%jsg+1) &
2250 call mom_error(fatal,"initialize_sponge_file: Array size mismatch for sponge data.")
2251 nz_data = siz(3)-1
2252 allocate(eta(isd:ied,jsd:jed,nz_data+1))
2253 allocate(dz(isd:ied,jsd:jed,nz_data))
2254 call mom_read_data(filename, eta_var, eta(:,:,:), g%Domain, scale=us%m_to_Z)
2255 do j=js,je ; do i=is,ie
2256 eta(i,j,nz_data+1) = -depth_tot(i,j)
2257 enddo ; enddo
2258 do k=nz_data,1,-1 ; do j=js,je ; do i=is,ie
2259 if (eta(i,j,k) < (eta(i,j,k+1) + gv%Angstrom_Z)) &
2260 eta(i,j,k) = eta(i,j,k+1) + gv%Angstrom_Z
2261 enddo ; enddo ; enddo
2262 do k=1,nz_data ; do j=js,je ; do i=is,ie
2263 dz(i,j,k) = eta(i,j,k)-eta(i,j,k+1)
2264 enddo ; enddo ; enddo
2265 deallocate(eta)
2266
2267 if (use_temperature) then
2268 allocate(tmp_t(isd:ied,jsd:jed,nz_data))
2269 allocate(tmp_s(isd:ied,jsd:jed,nz_data))
2270 call mom_read_data(filename, potemp_var, tmp_t(:,:,:), g%Domain, scale=us%degC_to_C)
2271 call mom_read_data(filename, salin_var, tmp_s(:,:,:), g%Domain, scale=us%ppt_to_S)
2272 endif
2273
2274 if (sponge_uv) then
2275 call initialize_ale_sponge(idamp, g, gv, param_file, ale_csp, dz, nz_data, idamp_u, idamp_v, &
2276 data_h_is_z=.true.)
2277 else
2278 call initialize_ale_sponge(idamp, g, gv, param_file, ale_csp, dz, nz_data, &
2279 data_h_is_z=.true.)
2280 endif
2281 if (use_temperature) then
2282 call set_up_ale_sponge_field(tmp_t, g, gv, tv%T, ale_csp, 'temp', &
2283 sp_long_name='temperature', sp_unit='degC s-1')
2284 call set_up_ale_sponge_field(tmp_s, g, gv, tv%S, ale_csp, 'salt', &
2285 sp_long_name='salinity', sp_unit='g kg-1 s-1')
2286 deallocate(tmp_s)
2287 deallocate(tmp_t)
2288 endif
2289 deallocate(dz)
2290
2291 if (sponge_uv) then
2292 filename = trim(inputdir)//trim(state_uv_file)
2293 call log_param(param_file, mdl, "INPUTDIR/SPONGE_STATE_UV_FILE", filename)
2294 if (.not.file_exists(filename, g%Domain)) &
2295 call mom_error(fatal, " initialize_sponges: Unable to open "//trim(filename))
2296 allocate(tmp_u(g%IsdB:g%IedB,jsd:jed,nz_data))
2297 allocate(tmp_v(isd:ied,g%JsdB:g%JedB,nz_data))
2298 call mom_read_vector(filename, u_var, v_var, tmp_u(:,:,:), tmp_v(:,:,:), g%Domain, scale=us%m_s_to_L_T)
2299 call set_up_ale_sponge_vel_field(tmp_u, tmp_v, g, gv, u, v, ale_csp)
2300 deallocate(tmp_u,tmp_v)
2301 endif
2302 else
2303 ! Initialize sponges without supplying sponge grid
2304 if (sponge_uv) then
2305 call initialize_ale_sponge(idamp, g, gv, us, param_file, ale_csp, idamp_u, idamp_v)
2306 else
2307 call initialize_ale_sponge(idamp, g, gv, us, param_file, ale_csp)
2308 endif
2309 ! The remaining calls to set_up_sponge_field can be in any order.
2310 if ( use_temperature) then
2311 call set_up_ale_sponge_field(filename, potemp_var, time, g, gv, us, tv%T, ale_csp, &
2312 'temp', sp_long_name='temperature', sp_unit='degC s-1', scale=us%degC_to_C)
2313 call set_up_ale_sponge_field(filename, salin_var, time, g, gv, us, tv%S, ale_csp, &
2314 'salt', sp_long_name='salinity', sp_unit='g kg-1 s-1', scale=us%ppt_to_S)
2315 endif
2316 if (sponge_uv) then
2317 filename = trim(inputdir)//trim(state_uv_file)
2318 call log_param(param_file, mdl, "INPUTDIR/SPONGE_STATE_UV_FILE", filename)
2319 if (.not.file_exists(filename, g%Domain)) &
2320 call mom_error(fatal, " initialize_sponges: Unable to open "//trim(filename))
2321 call set_up_ale_sponge_vel_field(filename, u_var, filename, v_var, time, g, gv, us, &
2322 ale_csp, u, v, scale=us%m_s_to_L_T)
2323 endif
2324 endif
2325 endif
2326
2327 if (sponge_uv .and. .not. use_ale) call mom_error(fatal,'initialize_sponges_file: '// &
2328 'UV damping to target values only available in ALE mode')
2329
2330
2331 contains
2332
2333 ! returns true if a separate idamp is provided for u and/or v
2334 logical function separate_idamp_for_uv()
2335 separate_idamp_for_uv = (lowercase(damping_file)/=lowercase(uv_damping_file) .or. &
2336 lowercase(idamp_var)/=lowercase(idamp_u_var) .or. lowercase(idamp_var)/=lowercase(idamp_v_var))
2337 end function separate_idamp_for_uv
2338
2339end subroutine initialize_sponges_file
2340
2341subroutine initialize_oda_incupd_file(G, GV, US, use_temperature, tv, h, u, v, param_file, &
2342 oda_incupd_CSp, restart_CS, Time)
2343 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
2344 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
2345 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2346 logical, intent(in) :: use_temperature !< If true, T & S are state variables.
2347 type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic
2348 !! variables.
2349 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
2350 intent(inout) :: h !< Layer thickness [H ~> m or kg m-2] (in)
2351
2352 real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
2353 intent(in) :: u !< The zonal velocity that is being
2354 !! initialized [L T-1 ~> m s-1]
2355 real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
2356 intent(in) :: v !< The meridional velocity that is being
2357 !! initialized [L T-1 ~> m s-1]
2358 type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters.
2359 type(oda_incupd_cs), pointer :: oda_incupd_CSp !< A pointer that is set to point to the control
2360 !! structure for this module.
2361 type(mom_restart_cs), intent(in) :: restart_CS !< MOM restart control structure
2362 type(time_type), intent(in) :: Time !< Time at the start of the run segment. Time_in
2363 !! overrides any value set for Time.
2364 ! Local variables
2365 real, allocatable, dimension(:,:,:) :: hoda ! The layer thickness increment and oda layer thickness [H ~> m or kg m-2]
2366 real, allocatable, dimension(:,:,:) :: tmp_tr ! A temporary array for reading oda tracer increments
2367 ! on the vertical grid of the input file, used for both
2368 ! temperatures [C ~> degC] and salinities [S ~> ppt]
2369 real, allocatable, dimension(:,:,:) :: tmp_u ! Temporary array for reading oda zonal velocity
2370 ! increments on the vertical grid of the input file [L T-1 ~> m s-1]
2371 real, allocatable, dimension(:,:,:) :: tmp_v ! Temporary array for reading oda meridional velocity
2372 ! increments on the vertical grid of the input file [L T-1 ~> m s-1]
2373
2374 integer :: is, ie, js, je, nz
2375 integer :: isd, ied, jsd, jed
2376
2377 integer, dimension(4) :: siz
2378 integer :: nz_data ! The size of the sponge source grid
2379 logical :: oda_inc ! input files are increments (true) or full fields (false)
2380 logical :: save_inc ! save increments if using full fields
2381 logical :: uv_inc ! use u and v increments
2382 logical :: reset_ncount ! reset ncount to zero if true
2383
2384 character(len=40) :: tempinc_var, salinc_var, uinc_var, vinc_var, h_var
2385 character(len=40) :: mdl = "initialize_oda_incupd_file"
2386 character(len=200) :: inc_file, uv_inc_file ! Strings for filenames
2387 character(len=200) :: filename, inputdir ! Strings for file/path and path.
2388
2389! logical :: use_ALE ! True if ALE is being used, False if in layered mode
2390
2391 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
2392 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
2393
2394 call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
2395 inputdir = slasher(inputdir)
2396
2397 call get_param(param_file, mdl, "ODA_INCUPD_FILE", inc_file, &
2398 "The name of the file with the T,S,h increments.", &
2399 fail_if_missing=.true.)
2400 call get_param(param_file, mdl, "ODA_INCUPD_INC", oda_inc, &
2401 "INCUPD files are increments and not full fields.", &
2402 default=.true.)
2403 if (.not.oda_inc) then
2404 call get_param(param_file, mdl, "ODA_INCUPD_SAVE", save_inc, &
2405 "If true, save the increments when using full fields.", &
2406 default=.false.)
2407 endif
2408 call get_param(param_file, mdl, "ODA_INCUPD_RESET_NCOUNT", reset_ncount, &
2409 "If True, reinitialize number of updates already done, ncount.", &
2410 default=.true.)
2411 if (.not.oda_inc .and. .not.reset_ncount) &
2412 call mom_error(fatal, " initialize_oda_incupd: restarting during update "// &
2413 "necessitates increments input file")
2414
2415 call get_param(param_file, mdl, "ODA_TEMPINC_VAR", tempinc_var, &
2416 "The name of the potential temperature inc. variable in "//&
2417 "ODA_INCUPD_FILE.", default="ptemp_inc")
2418 call get_param(param_file, mdl, "ODA_SALTINC_VAR", salinc_var, &
2419 "The name of the salinity inc. variable in "//&
2420 "ODA_INCUPD_FILE.", default="sal_inc")
2421 call get_param(param_file, mdl, "ODA_THK_VAR", h_var, &
2422 "The name of the layer thickness variable in "//&
2423 "ODA_INCUPD_FILE.", default="h")
2424 call get_param(param_file, mdl, "ODA_INCUPD_UV", uv_inc, &
2425 "use U,V increments.", &
2426 default=.true.)
2427 call get_param(param_file, mdl, "ODA_INCUPD_UV_FILE", uv_inc_file, &
2428 "The name of the file with the U,V increments.", &
2429 default=inc_file)
2430 call get_param(param_file, mdl, "ODA_UINC_VAR", uinc_var, &
2431 "The name of the zonal vel. inc. variable in "//&
2432 "ODA_INCUPD_FILE.", default="u_inc")
2433 call get_param(param_file, mdl, "ODA_VINC_VAR", vinc_var, &
2434 "The name of the meridional vel. inc. variable in "//&
2435 "ODA_INCUPD_FILE.", default="v_inc")
2436
2437! call get_param(param_file, mdl, "USE_REGRIDDING", use_ALE, default=.false., do_not_log=.true.)
2438
2439 ! Read in incremental update for tracers
2440 filename = trim(inc_file)
2441 if (scan(inc_file, '/')== 0) then ! prepend inputdir if only a filename is given
2442 filename = trim(inputdir)//trim(inc_file)
2443 endif
2444 call log_param(param_file, mdl, "INPUTDIR/ODA_INCUPD_FILE", filename)
2445 if (.not.file_exists(filename, g%Domain)) &
2446 call mom_error(fatal, " initialize_oda_incupd: Unable to open "//trim(filename))
2447
2448 call field_size(filename,h_var,siz,no_domain=.true.)
2449 if (siz(1) /= g%ieg-g%isg+1 .or. siz(2) /= g%jeg-g%jsg+1) &
2450 call mom_error(fatal,"initialize_oda_incupd_file: Array size mismatch for oda data.")
2451 nz_data = siz(3)
2452 ! get h increments
2453 allocate(hoda(isd:ied,jsd:jed,nz_data))
2454 call mom_read_data(filename, h_var , hoda(:,:,:), g%Domain, scale=us%m_to_Z)
2455 call initialize_oda_incupd( g, gv, us, param_file, oda_incupd_csp, hoda, nz_data, restart_cs)
2456 deallocate(hoda)
2457
2458 ! set-up T and S increments arrays
2459 if (use_temperature) then
2460 allocate(tmp_tr(isd:ied,jsd:jed,nz_data))
2461 ! temperature inc. in array Inc(1)
2462 call mom_read_data(filename, tempinc_var, tmp_tr(:,:,:), g%Domain, scale=us%degC_to_C)
2463 call set_up_oda_incupd_field(tmp_tr, g, gv, oda_incupd_csp)
2464 ! salinity inc. in array Inc(2)
2465 call mom_read_data(filename, salinc_var, tmp_tr(:,:,:), g%Domain, scale=us%ppt_to_S)
2466 call set_up_oda_incupd_field(tmp_tr, g, gv, oda_incupd_csp)
2467 deallocate(tmp_tr)
2468 endif
2469
2470 ! set-up U and V increments arrays
2471 if (uv_inc) then
2472 filename = trim(inputdir)//trim(uv_inc_file)
2473 call log_param(param_file, mdl, "INPUTDIR/ODA_INCUPD_UV_FILE", filename)
2474 if (.not.file_exists(filename, g%Domain)) &
2475 call mom_error(fatal, " initialize_oda_incupd_uv: Unable to open "//trim(filename))
2476 allocate(tmp_u(g%IsdB:g%IedB,jsd:jed,nz_data), source=0.0)
2477 allocate(tmp_v(isd:ied,g%JsdB:g%JedB,nz_data), source=0.0)
2478 call mom_read_vector(filename, uinc_var, vinc_var, tmp_u, tmp_v, g%Domain, scale=us%m_s_to_L_T)
2479 call set_up_oda_incupd_vel_field(tmp_u, tmp_v, g, gv, oda_incupd_csp)
2480 deallocate(tmp_u, tmp_v)
2481 endif
2482
2483 ! calculate increments if input are full fields
2484 if (oda_inc) then ! input are increments
2485 if (is_root_pe()) call mom_mesg("incupd using increments fields ")
2486 else ! inputs are full fields
2487 if (is_root_pe()) call mom_mesg("incupd using full fields ")
2488 call calc_oda_increments(h, tv, u, v, g, gv, us, oda_incupd_csp)
2489 if (save_inc) then
2490 call output_oda_incupd_inc(time, g, gv, param_file, oda_incupd_csp, us)
2491 endif
2492 endif ! not oda_inc
2493
2494end subroutine initialize_oda_incupd_file
2495
2496
2497!> This subroutine sets the 4 bottom depths at velocity points to be the
2498!! maximum of the adjacent depths.
2499subroutine set_velocity_depth_max(G)
2500 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
2501 ! Local variables
2502 integer :: i, j
2503
2504 do i=g%isd,g%ied-1 ; do j=g%jsd,g%jed
2505 g%Dblock_u(i,j) = g%mask2dCu(i,j) * max(g%bathyT(i,j), g%bathyT(i+1,j))
2506 g%Dopen_u(i,j) = g%Dblock_u(i,j)
2507 enddo ; enddo
2508 do i=g%isd,g%ied ; do j=g%jsd,g%jed-1
2509 g%Dblock_v(i,j) = g%mask2dCv(i,j) * max(g%bathyT(i,j), g%bathyT(i,j+1))
2510 g%Dopen_v(i,j) = g%Dblock_v(i,j)
2511 enddo ; enddo
2512end subroutine set_velocity_depth_max
2513
2514!> This subroutine sets the 4 bottom depths at velocity points to be the
2515!! minimum of the adjacent depths.
2516subroutine set_velocity_depth_min(G)
2517 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
2518 ! Local variables
2519 integer :: i, j
2520
2521 do i=g%isd,g%ied-1 ; do j=g%jsd,g%jed
2522 g%Dblock_u(i,j) = g%mask2dCu(i,j) * min(g%bathyT(i,j), g%bathyT(i+1,j))
2523 g%Dopen_u(i,j) = g%Dblock_u(i,j)
2524 enddo ; enddo
2525 do i=g%isd,g%ied ; do j=g%jsd,g%jed-1
2526 g%Dblock_v(i,j) = g%mask2dCv(i,j) * min(g%bathyT(i,j), g%bathyT(i,j+1))
2527 g%Dopen_v(i,j) = g%Dblock_v(i,j)
2528 enddo ; enddo
2529end subroutine set_velocity_depth_min
2530
2531!> This subroutine determines the isopycnal or other coordinate interfaces and
2532!! layer potential temperatures and salinities directly from a z-space file on
2533!! a latitude-longitude grid.
2534subroutine mom_temp_salt_initialize_from_z(h, tv, depth_tot, G, GV, US, PF, just_read, frac_shelf_h)
2535 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
2536 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
2537 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
2538 intent(out) :: h !< Layer thicknesses being initialized [H ~> m or kg m-2]
2539 type(thermo_var_ptrs), intent(inout) :: tv !< A structure pointing to various thermodynamic
2540 !! variables including temperature and salinity
2541 real, dimension(SZI_(G),SZJ_(G)), &
2542 intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
2543 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2544 type(param_file_type), intent(in) :: PF !< A structure indicating the open file
2545 !! to parse for model parameter values.
2546 logical, intent(in) :: just_read !< If true, this call will only read
2547 !! parameters without changing T or S.
2548 real, dimension(SZI_(G),SZJ_(G)), &
2549 optional, intent(in) :: frac_shelf_h !< The fraction of the grid cell covered
2550 !! by a floating ice shelf [nondim].
2551
2552 ! Local variables
2553 character(len=200) :: filename !< The name of an input file containing temperature
2554 !! and salinity in z-space; by default it is also used for ice shelf area.
2555 character(len=200) :: tfilename !< The name of an input file containing temperature in z-space.
2556 character(len=200) :: sfilename !< The name of an input file containing salinity in z-space.
2557 character(len=200) :: inputdir !! The directory where NetCDF input files are.
2558 character(len=200) :: mesg
2559
2560 type(eos_type), pointer :: eos => null()
2561 type(thermo_var_ptrs) :: tv_loc ! A temporary thermo_var container
2562 type(verticalgrid_type) :: GV_loc ! A temporary vertical grid structure
2563 ! This include declares and sets the variable "version".
2564# include "version_variable.h"
2565 character(len=40) :: mdl = "MOM_initialize_layers_from_Z" ! This module's name.
2566
2567 integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
2568 integer :: is, ie, js, je, nz ! compute domain indices
2569 integer :: isg, ieg, jsg, jeg ! global extent
2570 integer :: isd, ied, jsd, jed ! data domain indices
2571
2572 integer :: i, j, k, ks
2573 integer :: nkml ! The number of layers in the mixed layer.
2574
2575 integer :: inconsistent ! The total number of cells with in consistent topography and layer thicknesses.
2576 integer :: kd ! The number of levels in the input data
2577 integer :: nkd ! number of levels to use for regridding input arrays
2578 real :: eps_Z ! A negligibly thin layer thickness [Z ~> m].
2579 real :: eps_rho ! A negligibly small density difference [R ~> kg m-3].
2580 real :: PI_180 ! for conversion from degrees to radians [radian degree-1]
2581 real :: Hmix_default ! The default initial mixed layer depth [Z ~> m].
2582 real :: Hmix_depth ! The mixed layer depth in the initial condition [Z ~> m].
2583 real :: missing_value_temp ! The missing value in the input temperature field [C ~> degC]
2584 real :: missing_value_salt ! The missing value in the input salinity field [S ~> ppt]
2585 real :: tol_temp ! The tolerance for changes in temperature during the horizontal
2586 ! interpolation from an input dataset [C ~> degC]
2587 real :: tol_sal ! The tolerance for changes in salinity during the horizontal
2588 ! interpolation from an input dataset [S ~> ppt]
2589 logical :: correct_thickness ! If true, correct the column thicknesses to match the topography
2590 real :: h_tolerance ! A parameter that controls the tolerance when adjusting the
2591 ! thickness to fit the bathymetry [Z ~> m].
2592 real :: tol_dz_bot ! A tolerance for detecting inconsistent bottom depths when
2593 ! correct_thickness is false [Z ~> m]
2594 character(len=40) :: potemp_var, salin_var
2595
2596 integer, parameter :: niter=10 ! number of iterations for t/s adjustment to layer density
2597 logical :: adjust_temperature = .true. ! fit t/s to target densities
2598 real :: temp_land_fill ! A temperature value to use for land points [C ~> degC]
2599 real :: salt_land_fill ! A salinity value to use for land points [C ~> degC]
2600
2601 ! data arrays
2602 real, dimension(:), allocatable :: z_edges_in ! Input data interface heights or depths [Z ~> m]
2603 real, dimension(:), allocatable :: z_in ! Input data cell heights or depths [Z ~> m]
2604 real, dimension(:), allocatable :: Rb ! Interface densities [R ~> kg m-3]
2605 real, dimension(:,:,:), allocatable, target :: temp_z ! Input temperatures [C ~> degC]
2606 real, dimension(:,:,:), allocatable, target :: salt_z ! Input salinities [S ~> ppt]
2607 real, dimension(:,:,:), allocatable, target :: mask_z ! 1 for valid data points [nondim]
2608 real, dimension(:,:,:), allocatable :: rho_z ! Densities in Z-space [R ~> kg m-3]
2609 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: zi ! Interface heights [Z ~> m]
2610 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dz ! Layer thicknesses in height units [Z ~> m]
2611 real, dimension(SZI_(G),SZJ_(G)) :: Z_bottom ! The (usually negative) height of the seafloor
2612 ! relative to the surface [Z ~> m].
2613 integer, dimension(SZI_(G),SZJ_(G)) :: nlevs ! The number of levels in each column with valid data
2614 real, dimension(SZI_(G)) :: press ! Pressures [R L2 T-2 ~> Pa].
2615
2616 ! Local variables for ALE remapping
2617 real, dimension(:), allocatable :: hTarget ! Target thicknesses [Z ~> m].
2618 real, dimension(:,:,:), allocatable, target :: tmpT1dIn ! Input temperatures on a model-sized grid [C ~> degC]
2619 real, dimension(:,:,:), allocatable, target :: tmpS1dIn ! Input salinities on a model-sized grid [S ~> ppt]
2620 real, dimension(:,:,:), allocatable :: tmp_mask_in ! The valid data mask on a model-sized grid [nondim]
2621 real, dimension(:,:,:), allocatable :: dz1 ! Input grid thicknesses in depth units [Z ~> m]
2622 real, dimension(:,:,:), allocatable :: h1 ! Thicknesses on the input grid [H ~> m or kg m-2].
2623 real, dimension(:,:,:), allocatable :: dz_interface ! Change in position of interface due to
2624 ! regridding [H ~> m or kg m-2]
2625 real :: dz_neglect ! A negligibly small vertical layer extent used in
2626 ! remapping cell reconstructions [Z ~> m]
2627 real :: dz_neglect_edge ! A negligibly small vertical layer extent used in
2628 ! remapping edge value calculations [Z ~> m]
2629 real :: zTopOfCell, zBottomOfCell ! Heights in Z units [Z ~> m].
2630 type(regridding_cs) :: regridCS ! Regridding parameters and work arrays
2631 type(remapping_cs) :: remapCS ! Remapping parameters and work arrays
2632
2633 logical :: homogenize, useALEremapping, remap_full_column, remap_general, remap_old_alg
2634 integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.
2635 integer :: remap_answer_date ! The vintage of the order of arithmetic and expressions to use
2636 ! for remapping. Values below 20190101 recover the remapping
2637 ! answers from 2018, while higher values use more robust
2638 ! forms of the same remapping expressions.
2639 integer :: hor_regrid_answer_date ! The vintage of the order of arithmetic and expressions to use
2640 ! for horizontal regridding. Values below 20190101 recover the
2641 ! answers from 2018, while higher values use expressions that have
2642 ! been rearranged for rotational invariance.
2643 logical :: pre_gridded
2644 logical :: separate_mixed_layer ! If true, handle the mixed layers differently.
2645 logical :: density_extrap_bug ! If true use an expression with a vertical indexing bug for
2646 ! extrapolating the densities at the bottom of unstable profiles
2647 ! from data when finding the initial interface locations in
2648 ! layered mode from a dataset of T and S.
2649 character(len=64) :: remappingScheme
2650 logical :: om4_remap_via_sub_cells ! If true, use the OM4 remapping algorithm (only used if useALEremapping)
2651 logical :: do_conv_adj, ignore
2652 logical :: use_depth_based_time_fitler, use_adjust_interface_motion
2653 integer :: id_clock_routine, id_clock_ALE
2654
2655 id_clock_routine = cpu_clock_id('(Initialize from Z)', grain=clock_routine)
2656 id_clock_ale = cpu_clock_id('(Initialize from Z) ALE', grain=clock_loop)
2657
2658 call cpu_clock_begin(id_clock_routine)
2659
2660 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
2661 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
2662 isg = g%isg ; ieg = g%ieg ; jsg = g%jsg ; jeg = g%jeg
2663
2664 pi_180=atan(1.0)/45.
2665
2666 if (.not.just_read) call calltree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
2667 if (.not.just_read) call log_version(pf, mdl, version, "")
2668
2669 call get_param(pf, mdl, "INPUTDIR", inputdir, default=".")
2670 inputdir = slasher(inputdir)
2671
2672 eos => tv%eqn_of_state
2673
2674 call get_param(pf, mdl, "TEMP_SALT_Z_INIT_FILE", filename, &
2675 "The name of the z-space input file used to initialize "//&
2676 "temperatures (T) and salinities (S). If T and S are not "//&
2677 "in the same file, TEMP_Z_INIT_FILE and SALT_Z_INIT_FILE "//&
2678 "must be set.", default="temp_salt_z.nc", do_not_log=just_read)
2679 call get_param(pf, mdl, "TEMP_Z_INIT_FILE", tfilename, &
2680 "The name of the z-space input file used to initialize "//&
2681 "temperatures, only.", default=trim(filename), do_not_log=just_read)
2682 call get_param(pf, mdl, "SALT_Z_INIT_FILE", sfilename, &
2683 "The name of the z-space input file used to initialize "//&
2684 "temperatures, only.", default=trim(filename), do_not_log=just_read)
2685 filename = trim(inputdir)//trim(filename)
2686 tfilename = trim(inputdir)//trim(tfilename)
2687 sfilename = trim(inputdir)//trim(sfilename)
2688 call get_param(pf, mdl, "Z_INIT_FILE_PTEMP_VAR", potemp_var, &
2689 "The name of the potential temperature variable in "//&
2690 "TEMP_Z_INIT_FILE.", default="ptemp", do_not_log=just_read)
2691 call get_param(pf, mdl, "Z_INIT_FILE_SALT_VAR", salin_var, &
2692 "The name of the salinity variable in "//&
2693 "SALT_Z_INIT_FILE.", default="salt", do_not_log=just_read)
2694 call get_param(pf, mdl, "Z_INIT_HOMOGENIZE", homogenize, &
2695 "If True, then horizontally homogenize the interpolated "//&
2696 "initial conditions.", default=.false., do_not_log=just_read)
2697 call get_param(pf, mdl, "Z_INIT_ALE_REMAPPING", usealeremapping, &
2698 "If True, then remap straight to model coordinate from file.", &
2699 default=.false., do_not_log=just_read)
2700 call get_param(pf, mdl, "Z_INIT_REMAPPING_SCHEME", remappingscheme, &
2701 "The remapping scheme to use if using Z_INIT_ALE_REMAPPING "//&
2702 "is True.", default="PPM_IH4", do_not_log=just_read)
2703 call get_param(pf, mdl, "Z_INIT_REMAP_GENERAL", remap_general, &
2704 "If false, only initializes to z* coordinates. "//&
2705 "If true, allows initialization directly to general coordinates.", &
2706 default=.not.(gv%Boussinesq.or.gv%semi_Boussinesq) , do_not_log=just_read)
2707 call get_param(pf, mdl, "Z_INIT_REMAP_FULL_COLUMN", remap_full_column, &
2708 "If false, only reconstructs profiles for valid data points. "//&
2709 "If true, inserts vanished layers below the valid data.", &
2710 default=remap_general, do_not_log=just_read)
2711 call get_param(pf, mdl, "Z_INIT_REMAP_OLD_ALG", remap_old_alg, &
2712 "If false, uses the preferred remapping algorithm for initialization. "//&
2713 "If true, use an older, less robust algorithm for remapping.", &
2714 default=.false., do_not_log=just_read)
2715 call get_param(pf, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
2716 "This sets the default value for the various _ANSWER_DATE parameters.", &
2717 default=99991231, do_not_log=just_read)
2718 call get_param(pf, mdl, "TEMP_SALT_INIT_VERTICAL_REMAP_ONLY", pre_gridded, &
2719 "If true, initial conditions are on the model horizontal grid. " //&
2720 "Extrapolation over missing ocean values is done using an ICE-9 "//&
2721 "procedure with vertical ALE remapping .", &
2722 default=.false., do_not_log=just_read)
2723 if (usealeremapping) then
2724 call get_param(pf, mdl, "REMAPPING_ANSWER_DATE", remap_answer_date, &
2725 "The vintage of the expressions and order of arithmetic to use for remapping. "//&
2726 "Values below 20190101 result in the use of older, less accurate expressions "//&
2727 "that were in use at the end of 2018. Higher values result in the use of more "//&
2728 "robust and accurate forms of mathematically equivalent expressions.", &
2729 default=default_answer_date, do_not_log=just_read.or.(.not.gv%Boussinesq))
2730 call get_param(pf, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
2731 do_not_log=.true., default=.true.)
2732 call get_param(pf, mdl, "Z_INIT_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
2733 "If true, use the OM4 remapping-via-subcells algorithm for initialization. "//&
2734 "See REMAPPING_USE_OM4_SUBCELLS for more details. "//&
2735 "We recommend setting this option to false.", default=om4_remap_via_sub_cells)
2736 if (.not.gv%Boussinesq) remap_answer_date = max(remap_answer_date, 20230701)
2737 endif
2738 call get_param(pf, mdl, "HOR_REGRID_ANSWER_DATE", hor_regrid_answer_date, &
2739 "The vintage of the order of arithmetic for horizontal regridding. "//&
2740 "Dates before 20190101 give the same answers as the code did in late 2018, "//&
2741 "while later versions add parentheses for rotational symmetry. "//&
2742 "Dates after 20230101 use reproducing sums for global averages.", &
2743 default=default_answer_date, do_not_log=just_read.or.(.not.gv%Boussinesq))
2744 if (.not.gv%Boussinesq) hor_regrid_answer_date = max(hor_regrid_answer_date, 20230701)
2745
2746 if (.not.usealeremapping) then
2747 call get_param(pf, mdl, "ADJUST_THICKNESS", correct_thickness, &
2748 "If true, all mass below the bottom removed if the "//&
2749 "topography is shallower than the thickness input file "//&
2750 "would indicate.", default=.false., do_not_log=just_read)
2751 call get_param(pf, mdl, "THICKNESS_TOLERANCE", h_tolerance, &
2752 "A parameter that controls the tolerance when adjusting the "//&
2753 "thickness to fit the bathymetry. Used when ADJUST_THICKNESS=True.", &
2754 units="m", default=0.1, scale=us%m_to_Z, &
2755 do_not_log=(just_read.or..not.correct_thickness))
2756 call get_param(pf, mdl, "DZ_BOTTOM_TOLERANCE", tol_dz_bot, &
2757 "A tolerance for detecting inconsistent topography and input layer "//&
2758 "thicknesses when ADJUST_THICKNESS is false.", &
2759 units="m", default=1.0, scale=us%m_to_Z, &
2760 do_not_log=(just_read.or.correct_thickness))
2761
2762 call get_param(pf, mdl, "FIT_TO_TARGET_DENSITY_IC", adjust_temperature, &
2763 "If true, all the interior layers are adjusted to "//&
2764 "their target densities using mostly temperature "//&
2765 "This approach can be problematic, particularly in the "//&
2766 "high latitudes.", default=.true., do_not_log=just_read)
2767 call get_param(pf, mdl, "Z_INIT_SEPARATE_MIXED_LAYER", separate_mixed_layer, &
2768 "If true, distribute the topmost Z_INIT_HMIX_DEPTH of water over NKML layers, "//&
2769 "and do not correct the density of the topmost NKML+NKBL layers. Otherwise "//&
2770 "all layers are initialized based on the depths of their target densities.", &
2771 default=.false., do_not_log=just_read.or.(gv%nkml==0))
2772 if (gv%nkml == 0) separate_mixed_layer = .false.
2773 call get_param(pf, mdl, "MINIMUM_DEPTH", hmix_default, &
2774 units="m", default=0.0, scale=us%m_to_Z)
2775 call get_param(pf, mdl, "Z_INIT_HMIX_DEPTH", hmix_depth, &
2776 "The mixed layer depth in the initial conditions when Z_INIT_SEPARATE_MIXED_LAYER "//&
2777 "is set to true.", units="m", default=us%Z_to_m*hmix_default, scale=us%m_to_Z, &
2778 do_not_log=(just_read .or. .not.separate_mixed_layer))
2779 ! Reusing MINIMUM_DEPTH for the default mixed layer depth may be a strange choice, but
2780 ! it reproduces previous answers.
2781 call get_param(pf, mdl, "DENSITY_INTERP_TOLERANCE", eps_rho, &
2782 "A small density tolerance used when finding depths in a density profile.", &
2783 units="kg m-3", default=1.0e-10, scale=us%kg_m3_to_R, &
2784 do_not_log=usealeremapping.or.just_read)
2785 call get_param(pf, mdl, "LAYER_Z_INIT_IC_EXTRAP_BUG", density_extrap_bug, &
2786 "If true use an expression with a vertical indexing bug for extrapolating the "//&
2787 "densities at the bottom of unstable profiles from data when finding the "//&
2788 "initial interface locations in layered mode from a dataset of T and S.", &
2789 default=.false., do_not_log=just_read)
2790 endif
2791 call get_param(pf, mdl, "LAND_FILL_TEMP", temp_land_fill, &
2792 "A value to use to fill in ocean temperatures on land points.", &
2793 units="degC", default=0.0, scale=us%degC_to_C, do_not_log=just_read)
2794 call get_param(pf, mdl, "LAND_FILL_SALIN", salt_land_fill, &
2795 "A value to use to fill in ocean salinities on land points.", &
2796 units="ppt", default=35.0, scale=us%ppt_to_S, do_not_log=just_read)
2797 call get_param(pf, mdl, "HORIZ_INTERP_TOL_TEMP", tol_temp, &
2798 "The tolerance in temperature changes between iterations when interpolating "//&
2799 "from an input dataset using horiz_interp_and_extrap_tracer. This routine "//&
2800 "converges slowly, so an overly small tolerance can get expensive.", &
2801 units="degC", default=1.0e-3, scale=us%degC_to_C, do_not_log=just_read)
2802 call get_param(pf, mdl, "HORIZ_INTERP_TOL_SALIN", tol_sal, &
2803 "The tolerance in salinity changes between iterations when interpolating "//&
2804 "from an input dataset using horiz_interp_and_extrap_tracer. This routine "//&
2805 "converges slowly, so an overly small tolerance can get expensive.", &
2806 units="ppt", default=1.0e-3, scale=us%ppt_to_S, do_not_log=just_read)
2807 call get_param(pf, mdl, "REGRID_USE_DEPTH_BASED_TIME_FILTER", use_depth_based_time_fitler, &
2808 default=.true., do_not_log=.true.)
2809 call get_param(pf, mdl, "USE_ADJUST_INTERFACE_MOTION", use_adjust_interface_motion, &
2810 default=.true., do_not_log=.true.)
2811
2812 if (just_read) then
2813 if ((.not.usealeremapping) .and. adjust_temperature) &
2814 ! This call is just here to read and log the determine_temperature parameters
2815 call determine_temperature(tv%T, tv%S, gv%Rlay(1:nz), eos, tv%P_Ref, 0, &
2816 0, g, gv, us, pf, just_read=.true.)
2817 call cpu_clock_end(id_clock_routine)
2818 return ! All run-time parameters have been read, so return.
2819 endif
2820
2821 eps_z = gv%Angstrom_Z
2822
2823 ! Read input grid coordinates for temperature and salinity field
2824 ! in z-coordinate dataset. The file is REQUIRED to contain the
2825 ! following:
2826 !
2827 ! dimension variables:
2828 ! lon (degrees_E), lat (degrees_N), depth(meters)
2829 ! variables:
2830 ! ptemp(lon,lat,depth) : degC, potential temperature
2831 ! salt (lon,lat,depth) : ppt, salinity
2832 !
2833 ! The first record will be read if there are multiple time levels.
2834 ! The observation grid MUST tile the model grid. If the model grid extends
2835 ! to the North/South Pole past the limits of the input data, they are extrapolated using the average
2836 ! value at the northernmost/southernmost latitude.
2837
2838 call horiz_interp_and_extrap_tracer(tfilename, potemp_var, 1, &
2839 g, temp_z, mask_z, z_in, z_edges_in, missing_value_temp, &
2840 scale=us%degC_to_C, homogenize=homogenize, m_to_z=us%m_to_Z, &
2841 answer_date=hor_regrid_answer_date, ongrid=pre_gridded, tr_iter_tol=tol_temp)
2842
2843 call horiz_interp_and_extrap_tracer(sfilename, salin_var, 1, &
2844 g, salt_z, mask_z, z_in, z_edges_in, missing_value_salt, &
2845 scale=us%ppt_to_S, homogenize=homogenize, m_to_z=us%m_to_Z, &
2846 answer_date=hor_regrid_answer_date, ongrid=pre_gridded, tr_iter_tol=tol_sal)
2847
2848 kd = size(z_in,1)
2849
2850 ! Convert the sign convention of Z_edges_in.
2851 do k=1,size(z_edges_in,1) ; z_edges_in(k) = -z_edges_in(k) ; enddo
2852
2853 ! Convert T&S to Absolute Salinity and Conservative Temperature if using TEOS10 or NEMO
2854 call convert_temp_salt_for_teos10(temp_z, salt_z, g%HI, kd, mask_z, eos)
2855
2856 do j=js,je ; do i=is,ie
2857 z_bottom(i,j) = -depth_tot(i,j)
2858 enddo ; enddo
2859
2860 ! Done with horizontal interpolation.
2861 ! Now remap to model coordinates
2862 if (usealeremapping) then
2863 call cpu_clock_begin(id_clock_ale)
2864 nkd = max(gv%ke, kd)
2865
2866 ! Build the source grid and copy data onto model-shaped arrays with vanished layers
2867 allocate( tmp_mask_in(isd:ied,jsd:jed,nkd), source=0.0 )
2868 allocate( dz1(isd:ied,jsd:jed,nkd), source=0.0 )
2869 allocate( h1(isd:ied,jsd:jed,nkd), source=0.0 )
2870 allocate( tmpt1din(isd:ied,jsd:jed,nkd), source=0.0 )
2871 allocate( tmps1din(isd:ied,jsd:jed,nkd), source=0.0 )
2872 do j = js, je ; do i = is, ie
2873 if (g%mask2dT(i,j) > 0.) then
2874 ztopofcell = 0. ; zbottomofcell = 0.
2875 tmp_mask_in(i,j,1:kd) = mask_z(i,j,:)
2876 do k = 1, nkd
2877 if ((tmp_mask_in(i,j,k) > 0.) .and. (k <= kd)) then
2878 zbottomofcell = max( z_edges_in(k+1), z_bottom(i,j))
2879 tmpt1din(i,j,k) = temp_z(i,j,k)
2880 tmps1din(i,j,k) = salt_z(i,j,k)
2881 elseif (k>1) then
2882 zbottomofcell = z_bottom(i,j)
2883 tmpt1din(i,j,k) = tmpt1din(i,j,k-1)
2884 tmps1din(i,j,k) = tmps1din(i,j,k-1)
2885 else ! This next block should only ever be reached over land
2886 tmpt1din(i,j,k) = temp_land_fill
2887 tmps1din(i,j,k) = salt_land_fill
2888 endif
2889 dz1(i,j,k) = (ztopofcell - zbottomofcell)
2890 ztopofcell = zbottomofcell ! Bottom becomes top for next value of k
2891 enddo
2892 dz1(i,j,kd) = dz1(i,j,kd) + max(0., ztopofcell - z_bottom(i,j) )
2893 ! The max here is in case the data data is shallower than model
2894 endif ! mask2dT
2895 enddo ; enddo
2896 deallocate( tmp_mask_in )
2897
2898 ! Convert input thicknesses to units of H. In non-Boussinesq mode this is done by inverting
2899 ! integrals of specific volume in pressure, so it can be expensive.
2900 tv_loc = tv
2901 tv_loc%T => tmpt1din
2902 tv_loc%S => tmps1din
2903 gv_loc = gv
2904 gv_loc%ke = nkd
2905 call dz_to_thickness(dz1, tv_loc, h1, g, gv_loc, us)
2906
2907 ! Build the target grid (and set the model thickness to it)
2908
2909 call ale_initregridding( g, gv, us, g%max_depth, pf, mdl, regridcs ) ! sets regridCS
2910 if (remap_general) then
2911 dz_neglect = set_h_neglect(gv, remap_answer_date, dz_neglect_edge)
2912 else
2913 dz_neglect = set_dz_neglect(gv, us, remap_answer_date, dz_neglect_edge)
2914 endif
2915 call initialize_remapping( remapcs, remappingscheme, boundary_extrapolation=.false., &
2916 om4_remap_via_sub_cells=om4_remap_via_sub_cells, answer_date=remap_answer_date, &
2917 h_neglect=dz_neglect, h_neglect_edge=dz_neglect_edge)
2918
2919 ! Now remap from source grid to target grid, first setting reconstruction parameters
2920 if (remap_general) then
2921 call set_regrid_params( regridcs, min_thickness=0., &
2922 use_adjust_interface_motion=use_adjust_interface_motion, &
2923 use_depth_based_time_filter=use_depth_based_time_fitler)
2924 allocate( dz_interface(isd:ied,jsd:jed,nkd+1), source=0.) ! Need for argument to regridding_main() but is not used
2925
2926 call regridding_preadjust_reqs(regridcs, do_conv_adj, ignore)
2927 if (do_conv_adj) call convective_adjustment(g, gv_loc, h1, tv_loc)
2928 call regridding_main( remapcs, regridcs, g, gv_loc, us, h1, tv_loc, h, dz_interface, &
2929 frac_shelf_h=frac_shelf_h )
2930
2931 deallocate( dz_interface )
2932
2933 call ale_remap_scalar(remapcs, g, gv, nkd, h1, tmpt1din, h, tv%T, all_cells=remap_full_column, &
2934 old_remap=remap_old_alg )
2935 call ale_remap_scalar(remapcs, g, gv, nkd, h1, tmps1din, h, tv%S, all_cells=remap_full_column, &
2936 old_remap=remap_old_alg )
2937 else
2938 ! This is the old way of initializing to z* coordinates only
2939 allocate( htarget(nz) )
2940 htarget = getcoordinateresolution( regridcs )
2941 do j = js, je ; do i = is, ie
2942 dz(i,j,:) = 0.
2943 if (g%mask2dT(i,j) > 0.) then
2944 ! Build the target grid combining hTarget and topography
2945 ztopofcell = 0. ; zbottomofcell = 0.
2946 do k = 1, nz
2947 zbottomofcell = max( ztopofcell - htarget(k), z_bottom(i,j))
2948 dz(i,j,k) = ztopofcell - zbottomofcell
2949 ztopofcell = zbottomofcell ! Bottom becomes top for next value of k
2950 enddo
2951 else
2952 dz(i,j,:) = 0.
2953 endif ! mask2dT
2954 enddo ; enddo
2955 deallocate( htarget )
2956
2957 dz_neglect = set_dz_neglect(gv, us, remap_answer_date, dz_neglect_edge)
2958 call ale_remap_scalar(remapcs, g, gv, nkd, dz1, tmpt1din, dz, tv%T, all_cells=remap_full_column, &
2959 old_remap=remap_old_alg)
2960 call ale_remap_scalar(remapcs, g, gv, nkd, dz1, tmps1din, dz, tv%S, all_cells=remap_full_column, &
2961 old_remap=remap_old_alg)
2962
2963 if (gv%Boussinesq .or. gv%semi_Boussinesq) then
2964 ! This is a simple conversion of the target grid to thickness units that is not
2965 ! appropriate in non-Boussinesq mode.
2966 call dz_to_thickness_simple(dz, h, g, gv, us)
2967 else
2968 ! Convert dz into thicknesses in units of H using the equation of state as appropriate.
2969 call dz_to_thickness(dz, tv, h, g, gv, us)
2970 endif
2971 endif
2972
2973 deallocate( dz1 )
2974 deallocate( h1 )
2975 deallocate( tmpt1din )
2976 deallocate( tmps1din )
2977
2978 call cpu_clock_end(id_clock_ale)
2979
2980 else ! remap to isopycnal layer space
2981
2982 ! Next find interface positions using local arrays
2983 ! nlevs contains the number of valid data points in each column
2984 nlevs = int(sum(mask_z,dim=3))
2985
2986 ! Rb contains the layer interface densities
2987 allocate(rb(nz+1))
2988 do k=2,nz ; rb(k) = 0.5*(gv%Rlay(k-1)+gv%Rlay(k)) ; enddo
2989 rb(1) = 0.0
2990 if (nz>1) then
2991 rb(nz+1) = 2.0*gv%Rlay(nz) - gv%Rlay(nz-1)
2992 else
2993 rb(nz+1) = 2.0 * gv%Rlay(1)
2994 endif
2995
2996 nkml = 0 ; if (separate_mixed_layer) nkml = gv%nkml
2997
2998 press(:) = tv%P_Ref
2999 eosdom(:) = eos_domain(g%HI)
3000 allocate(rho_z(isd:ied,jsd:jed,kd))
3001 do k=1,kd ; do j=js,je
3002 call calculate_density(temp_z(:,j,k), salt_z(:,j,k), press, rho_z(:,j,k), eos, eosdom)
3003 enddo ; enddo
3004
3005 call find_interfaces(rho_z, z_in, kd, rb, z_bottom, zi, g, gv, us, nlevs, nkml, &
3006 hmix_depth, eps_z, eps_rho, density_extrap_bug)
3007
3008 deallocate(rho_z, rb)
3009
3010 dz(:,:,:) = 0.0
3011 if (correct_thickness) then
3012 call adjustetatofitbathymetry(g, gv, us, zi, dz, h_tolerance, dz_ref_eta=g%Z_ref)
3013 else
3014 do k=nz,1,-1 ; do j=js,je ; do i=is,ie
3015 if (zi(i,j,k) < (zi(i,j,k+1) + gv%Angstrom_Z)) then
3016 zi(i,j,k) = zi(i,j,k+1) + gv%Angstrom_Z
3017 dz(i,j,k) = gv%Angstrom_Z
3018 else
3019 dz(i,j,k) = zi(i,j,k) - zi(i,j,k+1)
3020 endif
3021 enddo ; enddo ; enddo
3022 inconsistent = 0
3023 do j=js,je ; do i=is,ie
3024 if (abs(zi(i,j,nz+1) - z_bottom(i,j)) > tol_dz_bot) &
3025 inconsistent = inconsistent + 1
3026 enddo ; enddo
3027 call sum_across_pes(inconsistent)
3028
3029 if ((inconsistent > 0) .and. (is_root_pe())) then
3030 write(mesg, '("Thickness initial conditions are inconsistent ",'// &
3031 '"with topography in ",I0," places.")') inconsistent
3032 call mom_error(warning, mesg)
3033 endif
3034 endif
3035
3036 call tracer_z_init_array(temp_z, z_edges_in, kd, zi, temp_land_fill, g, nz, nlevs, eps_z, tv%T)
3037 call tracer_z_init_array(salt_z, z_edges_in, kd, zi, salt_land_fill, g, nz, nlevs, eps_z, tv%S)
3038
3039 if (homogenize) then
3040 ! Horizontally homogenize data to produce perfectly "flat" initial conditions
3041 do k=1,nz
3042 call homogenize_field(tv%T(:,:,k), g, tmp_scale=us%C_to_degC, answer_date=hor_regrid_answer_date)
3043 call homogenize_field(tv%S(:,:,k), g, tmp_scale=us%S_to_ppt, answer_date=hor_regrid_answer_date)
3044 enddo
3045 endif
3046
3047 if (adjust_temperature) then
3048 ! Finally adjust to target density
3049 ks = 1 ; if (separate_mixed_layer) ks = gv%nk_rho_varies + 1
3050 call determine_temperature(tv%T, tv%S, gv%Rlay(1:nz), eos, tv%P_Ref, niter, &
3051 ks, g, gv, us, pf, just_read)
3052 endif
3053
3054 ! Now convert dz into thicknesses in units of H.
3055 call dz_to_thickness(dz, tv, h, g, gv, us)
3056
3057 endif ! useALEremapping
3058
3059 deallocate(z_in, z_edges_in, temp_z, salt_z, mask_z)
3060
3061 call pass_var(h, g%Domain)
3062 call pass_var(tv%T, g%Domain)
3063 call pass_var(tv%S, g%Domain)
3064
3065 call calltree_leave(trim(mdl)//'()')
3066 call cpu_clock_end(id_clock_routine)
3067
3069
3070
3071!> Find interface positions corresponding to interpolated depths in a density profile
3072subroutine find_interfaces(rho, zin, nk_data, Rb, Z_bot, zi, G, GV, US, nlevs, nkml, hml, &
3073 eps_z, eps_rho, density_extrap_bug)
3074 type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
3075 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure
3076 integer, intent(in) :: nk_data !< The number of levels in the input data
3077 real, dimension(SZI_(G),SZJ_(G),nk_data), &
3078 intent(in) :: rho !< Potential density in z-space [R ~> kg m-3]
3079 real, dimension(nk_data), intent(in) :: zin !< Input data levels [Z ~> m].
3080 real, dimension(SZK_(GV)+1), intent(in) :: Rb !< target interface densities [R ~> kg m-3]
3081 real, dimension(SZI_(G),SZJ_(G)), &
3082 intent(in) :: Z_bot !< The (usually negative) height of the seafloor
3083 !! relative to the surface [Z ~> m].
3084 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
3085 intent(out) :: zi !< The returned interface heights [Z ~> m]
3086 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
3087 integer, dimension(SZI_(G),SZJ_(G)), &
3088 intent(in) :: nlevs !< number of valid points in each column
3089 integer, intent(in) :: nkml !< number of mixed layer pieces to distribute over
3090 !! a depth of hml.
3091 real, intent(in) :: hml !< mixed layer depth [Z ~> m].
3092 real, intent(in) :: eps_z !< A negligibly small layer thickness [Z ~> m].
3093 real, intent(in) :: eps_rho !< A negligibly small density difference [R ~> kg m-3].
3094 logical, intent(in) :: density_extrap_bug !< If true use an expression with an
3095 !! indexing bug for projecting the densities at
3096 !! the bottom of unstable profiles from data when
3097 !! finding the initial interface locations in
3098 !! layered mode from a dataset of T and S.
3099
3100 ! Local variables
3101 real, dimension(nk_data) :: rho_ ! A column of densities [R ~> kg m-3]
3102 real, dimension(SZK_(GV)+1) :: zi_ ! A column interface heights (negative downward) [Z ~> m].
3103 real :: slope ! The rate of change of height with density [Z R-1 ~> m4 kg-1]
3104 real :: drhodz ! A local vertical density gradient [R Z-1 ~> kg m-4]
3105 real, parameter :: zoff = 0.999 ! A small fractional adjustment to the density differences [nondim]
3106 logical :: unstable ! True if the column is statically unstable anywhere.
3107 integer :: nlevs_data ! The number of data values in a column.
3108 logical :: work_down ! This indicates whether this pass goes up or down the water column.
3109 integer :: k_int, lo_int, hi_int, mid
3110 integer :: i, j, k, is, ie, js, je, nz
3111
3112 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
3113
3114 zi(:,:,:) = 0.0
3115
3116 do j=js,je ; do i=is,ie
3117 nlevs_data = nlevs(i,j)
3118 do k=1,nlevs_data ; rho_(k) = rho(i,j,k) ; enddo
3119
3120 unstable = .true.
3121 work_down = .true.
3122 do while (unstable)
3123 ! Modify the input profile until it no longer has densities that decrease with depth.
3124 unstable = .false.
3125 if (work_down) then
3126 do k=2,nlevs_data-1 ; if (rho_(k) - rho_(k-1) < 0.0) then
3127 if (k == 2) then
3128 rho_(k-1) = rho_(k) - eps_rho
3129 else
3130 drhodz = (rho_(k+1)-rho_(k-1)) / (zin(k+1)-zin(k-1))
3131 if (drhodz < 0.0) unstable = .true.
3132 rho_(k) = rho_(k-1) + drhodz*zoff*(zin(k)-zin(k-1))
3133 endif
3134 endif ; enddo
3135 work_down = .false.
3136 else
3137 do k=nlevs_data-1,2,-1 ; if (rho_(k+1) - rho_(k) < 0.0) then
3138 if (k == nlevs_data-1) then
3139 if (density_extrap_bug) then
3140 rho_(k+1) = rho_(k-1) + eps_rho
3141 else
3142 rho_(k+1) = rho_(k) + eps_rho
3143 endif
3144 else
3145 drhodz = (rho_(k+1)-rho_(k-1)) / (zin(k+1)-zin(k-1))
3146 if (drhodz < 0.0) unstable = .true.
3147 rho_(k) = rho_(k+1) - drhodz*(zin(k+1)-zin(k))
3148 endif
3149 endif ; enddo
3150 work_down = .true.
3151 endif
3152 enddo
3153
3154 ! Find and store the interface depths.
3155 zi_(1) = 0.0
3156 if (nlevs_data < 1) then
3157 ! There is no data to use, so set the interfaces at the bottom.
3158 do k=2,nz ; zi_(k) = z_bot(i,j) ; enddo
3159 elseif (nlevs_data == 1) then
3160 ! There is data for only one input layer, so set the interfaces at the bottom or top,
3161 ! depending on how their target densities compare with the one data point.
3162 do k=2,nz
3163 if (rb(k) < rho_(1)) then ; zi_(k) = 0.0
3164 else ; zi_(k) = z_bot(i,j) ; endif
3165 enddo
3166 else
3167 do k=2,nz
3168 ! Find the value of k_int in the list of rho_ where rho_(k_int) <= Rb(K) < rho_(k_int+1).
3169 ! This might be made a little faster by exploiting the fact that Rb is
3170 ! monotonically increasing and not resetting lo_int back to 1 inside the K loop.
3171 lo_int = 1 ; hi_int = nlevs_data
3172 do while (lo_int < hi_int)
3173 mid = (lo_int+hi_int) / 2
3174 if (rb(k) < rho_(mid)) then ; hi_int = mid
3175 else ; lo_int = mid+1 ; endif
3176 enddo
3177 k_int = max(1, lo_int-1)
3178
3179 ! Linearly interpolate to find the depth, zi_, where Rb would be found.
3180 slope = (zin(k_int+1) - zin(k_int)) / max(rho_(k_int+1) - rho_(k_int), eps_rho)
3181 zi_(k) = -1.0*(zin(k_int) + slope*(rb(k)-rho_(k_int)))
3182 zi_(k) = min(max(zi_(k), z_bot(i,j)), -1.0*hml)
3183 enddo
3184 endif
3185 zi_(nz+1) = z_bot(i,j)
3186 if (nkml > 0) then ; do k=2,nkml+1
3187 zi_(k) = max(hml*((1.0-real(k))/real(nkml)), z_bot(i,j))
3188 enddo ; endif
3189 do k=nz,max(nkml+2,2),-1
3190 if (zi_(k) < zi_(k+1) + eps_z) zi_(k) = zi_(k+1) + eps_z
3191 if (zi_(k) > -1.0*hml) zi_(k) = max(-1.0*hml, z_bot(i,j))
3192 enddo
3193
3194 do k=1,nz+1
3195 zi(i,j,k) = zi_(k)
3196 enddo
3197 enddo ; enddo ! i- and j- loops
3198
3199end subroutine find_interfaces
3200
3201!> Run simple unit tests
3202subroutine mom_state_init_tests(G, GV, US, tv)
3203 type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
3204 type(verticalgrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
3205 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
3206 type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure.
3207
3208 ! Local variables
3209 integer, parameter :: nk=5
3210 real, dimension(nk) :: T, T_t, T_b ! Temperatures [C ~> degC]
3211 real, dimension(nk) :: S, S_t, S_b ! Salinities [S ~> ppt]
3212 real, dimension(nk) :: rho ! Layer density [R ~> kg m-3]
3213 real, dimension(nk) :: h ! Layer thicknesses [H ~> m or kg m-2]
3214 real, dimension(nk) :: z ! Height of layer center [Z ~> m]
3215 real, dimension(nk+1) :: e ! Interface heights [Z ~> m]
3216 real :: T_ref ! A reference temperature [C ~> degC]
3217 real :: S_ref ! A reference salinity [S ~> ppt]
3218 real :: P_tot, P_t, P_b ! Pressures [R L2 T-2 ~> Pa]
3219 real :: z_out ! Output height [Z ~> m]
3220 real :: I_z_scale ! The inverse of the height scale for prescribed gradients [Z-1 ~> m-1]
3221 real :: z_tol ! The tolerance with which to find the depth matching a specified pressure [Z ~> m].
3222 integer :: k
3223 type(remapping_cs), pointer :: remap_CS => null()
3224
3225 i_z_scale = 1.0 / (500.0*us%m_to_Z)
3226 do k = 1, nk
3227 h(k) = 100.0*gv%m_to_H
3228 enddo
3229 e(1) = 0.
3230 do k = 1, nk
3231 e(k+1) = e(k) - gv%H_to_Z * h(k)
3232 enddo
3233 p_tot = 0.
3234 t_ref = 20.0*us%degC_to_C
3235 s_ref = 35.0*us%ppt_to_S
3236 z_tol = 1.0e-5*us%m_to_Z
3237 do k = 1, nk
3238 z(k) = 0.5 * ( e(k) + e(k+1) )
3239 t_t(k) = t_ref + (0. * i_z_scale) * e(k)
3240 t(k) = t_ref + (0. * i_z_scale)*z(k)
3241 t_b(k) = t_ref + (0. * i_z_scale)*e(k+1)
3242 s_t(k) = s_ref - (0. * i_z_scale)*e(k)
3243 s(k) = s_ref + (0. * i_z_scale)*z(k)
3244 s_b(k) = s_ref - (0. * i_z_scale)*e(k+1)
3245 call calculate_density(0.5*(t_t(k)+t_b(k)), 0.5*(s_t(k)+s_b(k)), -gv%Rho0*gv%g_Earth*z(k), &
3246 rho(k), tv%eqn_of_state)
3247 p_tot = p_tot + gv%g_Earth * rho(k) * gv%H_to_Z*h(k)
3248 enddo
3249
3250 p_t = 0.
3251 do k = 1, nk
3252 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, &
3253 gv%Rho0, gv%g_Earth, tv%eqn_of_state, us, p_b, z_out, z_tol=z_tol, &
3254 frac_dp_bugfix=.false.)
3255 write(0,*) k, us%RL2_T2_to_Pa*p_t, us%RL2_T2_to_Pa*p_b, 0.5*us%RL2_T2_to_Pa*p_tot, &
3256 us%Z_to_m*e(k), us%Z_to_m*e(k+1), us%Z_to_m*z_out
3257 p_t = p_b
3258 enddo
3259 write(0,*) us%RL2_T2_to_Pa*p_b, us%RL2_T2_to_Pa*p_tot
3260
3261 write(0,*) ''
3262 write(0,*) ' ==================================================================== '
3263 write(0,*) ''
3264 write(0,*) gv%H_to_m*h(:)
3265
3266 ! For consistency with the usual call, add the following:
3267 ! if (use_remapping) then
3268 ! allocate(remap_CS)
3269 ! call initialize_remapping(remap_CS, 'PLM', boundary_extrapolation=.true., &
3270 ! h_neglect=GV%H_subroundoff, h_neglect_edge=GV%H_subroundoff)
3271 ! endif
3272 call cut_off_column_top(nk, tv, gv, us, gv%g_Earth, -e(nk+1), gv%Angstrom_H, &
3273 t, t_t, t_b, s, s_t, s_b, 0.5*p_tot, h, remap_cs, z_tol=z_tol, &
3274 frac_dp_bugfix=.false.)
3275 write(0,*) gv%H_to_m*h(:)
3276 if (associated(remap_cs)) deallocate(remap_cs)
3277
3278end subroutine mom_state_init_tests
3279
3280end module mom_state_initialization