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