src/initialization/MOM_state_initialization.F90

portedportable, not yet portedexecuted, not portableexecutable, not hit by this run

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.
6		module MOM_state_initialization
7
8		use MOM_debugging, only : hchksum, qchksum, uvchksum
9		use MOM_density_integrals, only : int_specific_vol_dp
10		use MOM_density_integrals, only : find_depth_of_pressure_in_cell
11		use MOM_coms, only : max_across_PEs, min_across_PEs, reproducing_sum
12		use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
13		use MOM_cpu_clock, only : CLOCK_ROUTINE, CLOCK_LOOP
14		use MOM_domains, only : pass_var, pass_vector, sum_across_PEs, broadcast
15		use MOM_domains, only : root_PE, To_All, SCALAR_PAIR, CGRID_NE, AGRID
16		use MOM_error_handler, only : MOM_mesg, MOM_error, FATAL, WARNING, is_root_pe
17		use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint
18		use MOM_file_parser, only : get_param, read_param, log_param, param_file_type
19		use MOM_file_parser, only : log_version
20		use MOM_get_input, only : directories
21		use MOM_grid, only : ocean_grid_type, isPointInCell
22		use MOM_interface_heights, only : find_eta, dz_to_thickness, dz_to_thickness_simple
23		use MOM_interface_heights, only : calc_derived_thermo
24		use MOM_io, only : file_exists, field_size, MOM_read_data, MOM_read_vector, slasher
25		use MOM_open_boundary, only : ocean_OBC_type, open_boundary_test_extern_h
26		use MOM_open_boundary, only : fill_temp_salt_segments, setup_OBC_tracer_reservoirs
27		use MOM_open_boundary, only : fill_thickness_segments
28		use MOM_open_boundary, only : set_initialized_OBC_tracer_reservoirs
29		use MOM_restart, only : restore_state, is_new_run, copy_restart_var, copy_restart_vector
30		use MOM_restart, only : restart_registry_lock, MOM_restart_CS
31		use MOM_sponge, only : set_up_sponge_field, set_up_sponge_ML_density
32		use MOM_sponge, only : initialize_sponge, sponge_CS
33		use MOM_ALE_sponge, only : set_up_ALE_sponge_field, set_up_ALE_sponge_vel_field
34		use MOM_ALE_sponge, only : ALE_sponge_CS, initialize_ALE_sponge
35		use MOM_string_functions, only : uppercase, lowercase
36		use MOM_time_manager, only : time_type, operator(/=)
37		use MOM_tracer_registry, only : tracer_registry_type
38		use MOM_unit_scaling, only : unit_scale_type
39		use MOM_variables, only : thermo_var_ptrs
40		use MOM_verticalGrid, only : setVerticalGridAxes, verticalGrid_type
41		use MOM_EOS, only : calculate_density, calculate_density_derivs, EOS_type, EOS_domain
42		use MOM_EOS, only : convert_temp_salt_for_TEOS10
43		use user_initialization, only : user_initialize_thickness, user_initialize_velocity
44		use user_initialization, only : user_init_temperature_salinity, user_set_OBC_data
45		use user_initialization, only : user_initialize_sponges
46		use DOME_initialization, only : DOME_initialize_thickness
47		use DOME_initialization, only : DOME_set_OBC_data
48		use DOME_initialization, only : DOME_initialize_sponges
49		use ISOMIP_initialization, only : ISOMIP_initialize_thickness
50		use ISOMIP_initialization, only : ISOMIP_initialize_sponges
51		use ISOMIP_initialization, only : ISOMIP_initialize_temperature_salinity
52		use RGC_initialization, only : RGC_initialize_sponges
53		use baroclinic_zone_initialization, only : baroclinic_zone_init_temperature_salinity
54		use benchmark_initialization, only : benchmark_initialize_thickness
55		use benchmark_initialization, only : benchmark_init_temperature_salinity
56		use Neverworld_initialization, only : Neverworld_initialize_thickness
57		use circle_obcs_initialization, only : circle_obcs_initialize_thickness
58		use lock_exchange_initialization, only : lock_exchange_initialize_thickness
59		use external_gwave_initialization, only : external_gwave_initialize_thickness
60		use DOME2d_initialization, only : DOME2d_initialize_thickness
61		use DOME2d_initialization, only : DOME2d_initialize_temperature_salinity
62		use DOME2d_initialization, only : DOME2d_initialize_sponges
63		use adjustment_initialization, only : adjustment_initialize_thickness
64		use adjustment_initialization, only : adjustment_initialize_temperature_salinity
65		use sloshing_initialization, only : sloshing_initialize_thickness
66		use sloshing_initialization, only : sloshing_initialize_temperature_salinity
67		use seamount_initialization, only : seamount_initialize_thickness
68		use seamount_initialization, only : seamount_initialize_temperature_salinity
69		use dumbbell_initialization, only : dumbbell_initialize_thickness
70		use dumbbell_initialization, only : dumbbell_initialize_temperature_salinity
71		use Phillips_initialization, only : Phillips_initialize_thickness
72		use Phillips_initialization, only : Phillips_initialize_velocity
73		use Phillips_initialization, only : Phillips_initialize_sponges
74		use Rossby_front_2d_initialization, only : Rossby_front_initialize_thickness
75		use Rossby_front_2d_initialization, only : Rossby_front_initialize_temperature_salinity
76		use Rossby_front_2d_initialization, only : Rossby_front_initialize_velocity
77		use SCM_CVMix_tests, only: SCM_CVMix_tests_TS_init
78		use dyed_channel_initialization, only : dyed_channel_set_OBC_tracer_data
79		use dyed_obcs_initialization, only : dyed_obcs_set_OBC_data
80		use supercritical_initialization, only : supercritical_set_OBC_data
81		use soliton_initialization, only : soliton_initialize_velocity
82		use soliton_initialization, only : soliton_initialize_thickness
83		use BFB_initialization, only : BFB_initialize_sponges_southonly
84		use dense_water_initialization, only : dense_water_initialize_TS
85		use dense_water_initialization, only : dense_water_initialize_sponges
86		use dumbbell_initialization, only : dumbbell_initialize_sponges
87		use MOM_tracer_Z_init, only : tracer_Z_init_array, determine_temperature
88		use MOM_ALE, only : ALE_initRegridding, ALE_CS, ALE_initThicknessToCoord
89		use MOM_ALE, only : ALE_remap_scalar, ALE_regrid_accelerated, TS_PLM_edge_values
90		use MOM_regridding, only : regridding_CS, set_regrid_params, getCoordinateResolution
91		use MOM_regridding, only : regridding_main, regridding_preadjust_reqs, convective_adjustment
92		use MOM_regridding, only : set_dz_neglect, set_h_neglect
93		use MOM_remapping, only : remapping_CS, initialize_remapping, remapping_core_h
94		use MOM_horizontal_regridding, only : horiz_interp_and_extrap_tracer, homogenize_field
95		use MOM_oda_incupd, only: oda_incupd_CS, initialize_oda_incupd_fixed, initialize_oda_incupd
96		use MOM_oda_incupd, only: set_up_oda_incupd_field, set_up_oda_incupd_vel_field
97		use MOM_oda_incupd, only: calc_oda_increments, output_oda_incupd_inc
98
99		implicit none ; private
100
101		#include <MOM_memory.h>
102
103		public MOM_initialize_state, MOM_initialize_OBCs
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
110		character(len=40) :: mdl = "MOM_state_initialization" !< This module's name.
111
112		contains
113
114		!> Initialize temporally evolving fields, either as initial
115		!! conditions or by reading them from a restart (or saves) file.
116	3	subroutine 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	0	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	2	real :: depth_tot(SZI_(G),SZJ_(G)) ! The nominal total depth of the ocean [Z ~> m]
160	2	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	1	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
196	1	Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
197	1	isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
198	1	IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB
199
200	1	call callTree_enter("MOM_initialize_state(), MOM_state_initialization.F90")
201	1	call log_version(PF, mdl, version, "")
202	1	call get_param(PF, mdl, "DEBUG", debug, default=.false.)
203
204	1	new_sim = is_new_run(restart_CS)
205	1	just_read = .not.new_sim
206
207		call get_param(PF, mdl, "INPUTDIR", inputdir, &
208	1	"The directory in which input files are found.", default=".")
209	1	inputdir = slasher(inputdir)
210
211	1	use_temperature = associated(tv%T)
212	1	useALE = associated(ALE_CSp)
213	1	use_EOS = associated(tv%eqn_of_state)
214	1	if (use_EOS) eos => tv%eqn_of_state
215	1	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	1	if (new_sim) then
223	1	call MOM_mesg("Run initialized internally.", 3)
224
225	1	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	657976	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	8773	do j=jsd,jed ; do i=isd,ied
243	8772	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	1	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	1	default=.false., do_not_log=just_read)
260
261	1	convert = new_sim ! Thicknesses are initialized in height units in most cases.
262	1	if (from_Z_file) then
263		! Initialize thickness and T/S from z-coordinate data in a file.
264	0	if (.NOT.use_temperature) call MOM_error(FATAL,"MOM_initialize_state : "//&
265	0	"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	0	just_read=just_read, frac_shelf_h=frac_shelf_h)
269	0	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	1	default="uniform", do_not_log=just_read)
304	2	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	0	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	0	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	0	mass_file=.true., just_read=just_read)
314	0	convert = .false.
315		case ("coord")
316	0	if (new_sim .and. useALE) then
317	0	call ALE_initThicknessToCoord( ALE_CSp, G, GV, dz, height_units=.true. )
318	0	elseif (new_sim) then
319		call MOM_error(FATAL, "MOM_initialize_state: USE_REGRIDDING must be True "//&
320	0	"for THICKNESS_CONFIG of 'coord'")
321		endif
322		case ("uniform"); call initialize_thickness_uniform(dz, depth_tot, G, GV, PF, &
323	1	just_read=just_read)
324		case ("list"); call initialize_thickness_list(dz, depth_tot, G, GV, US, PF, &
325	0	just_read=just_read)
326		case ("param"); call initialize_thickness_param(dz, depth_tot, G, GV, US, PF, &
327	0	just_read=just_read)
328		case ("DOME"); call DOME_initialize_thickness(dz, depth_tot, G, GV, PF, &
329	0	just_read=just_read)
330		case ("ISOMIP"); call ISOMIP_initialize_thickness(dz, depth_tot, G, GV, US, PF, tv, &
331	0	just_read=just_read)
332		case ("benchmark"); call benchmark_initialize_thickness(dz, depth_tot, G, GV, US, PF, &
333	0	tv%eqn_of_state, tv%P_Ref, just_read=just_read)
334		case ("Neverworld","Neverland"); call Neverworld_initialize_thickness(dz, depth_tot, &
335	0	G, GV, US, PF, tv%P_Ref)
336	0	case ("search"); call initialize_thickness_search()
337		case ("circle_obcs"); call circle_obcs_initialize_thickness(dz, depth_tot, G, GV, US, PF, &
338	0	just_read=just_read)
339		case ("lock_exchange"); call lock_exchange_initialize_thickness(dz, G, GV, US, &
340	0	PF, just_read=just_read)
341		case ("external_gwave"); call external_gwave_initialize_thickness(dz, G, GV, US, &
342	0	PF, just_read=just_read)
343		case ("DOME2D"); call DOME2d_initialize_thickness(dz, depth_tot, G, GV, US, PF, &
344	0	just_read=just_read)
345		case ("adjustment2d"); call adjustment_initialize_thickness(dz, G, GV, US, &
346	0	PF, just_read=just_read)
347		case ("sloshing"); call sloshing_initialize_thickness(dz, depth_tot, G, GV, US, PF, &
348	0	just_read=just_read)
349		case ("seamount"); call seamount_initialize_thickness(dz, depth_tot, G, GV, US, PF, &
350	0	just_read=just_read)
351		case ("dumbbell"); call dumbbell_initialize_thickness(dz, depth_tot, G, GV, US, PF, &
352	0	just_read=just_read)
353		case ("soliton"); call soliton_initialize_thickness(dz, depth_tot, G, GV, US, PF, &
354	0	just_read=just_read)
355		case ("phillips"); call Phillips_initialize_thickness(dz, depth_tot, G, GV, US, PF, &
356	0	just_read=just_read)
357		case ("rossby_front")
358	0	call Rossby_front_initialize_thickness(h, G, GV, US, PF, just_read=just_read)
359	0	convert = .false. ! Rossby_front initialization works directly in thickness units.
360		case ("USER"); call user_initialize_thickness(dz, G, GV, PF, &
361	0	just_read=just_read)
362		case default ; call MOM_error(FATAL, "MOM_initialize_state: "//&
363	2	"Unrecognized layer thickness configuration "//trim(h_config))
364		end select
365
366		! Initialize temperature and salinity (T and S).
367	1	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	1	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	1	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	0	"THICKNESS_CONFIG = "//trim(h_config)//".")
398		end select ; endif
399
400	2	select case (trim(config))
401		case ("fit"); call initialize_temp_salt_fit(tv%T, tv%S, G, GV, US, PF, &
402	0	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	0	PF, just_read=just_read)
405		case ("benchmark"); call benchmark_init_temperature_salinity(tv%T, tv%S, &
406	0	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	0	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	0	just_read=just_read)
411		case ("DOME2D"); call DOME2d_initialize_temperature_salinity (tv%T, tv%S, dz, &
412	0	G, GV, US, PF, just_read=just_read)
413		case ("ISOMIP"); call ISOMIP_initialize_temperature_salinity (tv%T, tv%S, dz, &
414	0	depth_tot, G, GV, US, PF, eos, just_read=just_read)
415		case ("adjustment2d"); call adjustment_initialize_temperature_salinity ( tv%T, &
416	0	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	0	tv%S, dz, depth_tot, G, GV, US, PF, just_read=just_read)
419		case ("sloshing"); call sloshing_initialize_temperature_salinity(tv%T, &
420	0	tv%S, dz, G, GV, US, PF, just_read=just_read)
421		case ("seamount"); call seamount_initialize_temperature_salinity(tv%T, &
422	1	tv%S, dz, G, GV, US, PF, just_read=just_read)
423		case ("dumbbell"); call dumbbell_initialize_temperature_salinity(tv%T, &
424	0	tv%S, dz, G, GV, US, PF, just_read=just_read)
425		case ("rossby_front")
426	0	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	0	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	0	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	0	dz, just_read=just_read)
433		case ("USER"); call user_init_temperature_salinity(tv%T, tv%S, G, GV, PF, &
434	0	just_read=just_read)
435		case default ; call MOM_error(FATAL, "MOM_initialize_state: "//&
436	2	"Unrecognized Temp & salt configuration "//trim(config))
437		end select
438		endif
439		endif ! not from_Z_file.
440
441	1	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	0	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	0	"specified for the reservoir concentrations.", default=enable_bugs, do_not_log=.true.)
449	0	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	0	call pass_var(tv%T, G%Domain, complete=.false.)
453	0	call pass_var(tv%S, G%Domain, complete=.true.)
454	0	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	1	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	1	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	1	do_not_log=just_read)
471	1	if (depress_sfc .and. trim_ic_for_p_surf) call MOM_error(FATAL, "MOM_initialize_state: "//&
472	0	"DEPRESS_INITIAL_SURFACE and TRIM_IC_FOR_P_SURF are exclusive and cannot both be True")
473
474	1	if (new_sim .and. debug .and. (depress_sfc .or. trim_ic_for_p_surf)) &
475	0	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	1	if (depress_sfc) then
479	0	call depress_surface(h, G, GV, US, PF, tv, just_read=just_read)
480	1	elseif (trim_ic_for_p_surf) then
481	0	call trim_for_ice(PF, G, GV, US, ALE_CSp, tv, h, just_read=just_read)
482	1	elseif (new_sim .and. use_ice_shelf .and. present(mass_shelf)) then
483	0	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	1	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	1	default=.false., do_not_log=just_read)
494	1	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	0	default=1, do_not_log=just_read)
499
500		call get_param(PF, mdl, "DT", dt, "Timestep", &
501	0	units="s", scale=US%s_to_T, fail_if_missing=.true.)
502
503	0	if (new_sim .and. debug) &
504	0	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	0	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	1	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	1	do_not_log=just_read)
528	2	select case (trim(config))
529	0	case ("file"); call initialize_velocity_from_file(u, v, G, GV, US, PF, just_read)
530	1	case ("zero"); call initialize_velocity_zero(u, v, G, GV, PF, just_read)
531	0	case ("uniform"); call initialize_velocity_uniform(u, v, G, GV, US, PF, just_read)
532	0	case ("circular"); call initialize_velocity_circular(u, v, G, GV, US, PF, just_read)
533	0	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	0	G, GV, US, PF, just_read)
536	0	case ("soliton"); call soliton_initialize_velocity(u, v, G, GV, US, PF, just_read)
537	0	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	2	"Unrecognized velocity configuration "//trim(config))
540		end select
541
542	1	if (new_sim) call pass_vector(u, v, G%Domain)
543	1	if (debug .and. new_sim) then
544	0	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	1	"everywhere in the domain.", default=.false.)
555	1	if (use_oda_incupd) then
556	0	call restart_registry_lock(restart_CS, unlocked=.true.)
557	0	call initialize_oda_incupd_fixed(G, GV, US, oda_incupd_CSp, restart_CS)
558	0	call restart_registry_lock(restart_CS)
559		endif
560
561	1	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	0	call restore_state(dirs%input_filename, dirs%restart_input_dir, Time, G, restart_CS)
565	0	if (present(Time_in)) then
566	0	if (verify_restart_time .and. (Time /= Time_in)) call MOM_error(FATAL, &
567	0	"MOM6 attempted to restart from a file from a different time than given by Time_in.")
568	0	Time = Time_in
569		endif
570		call get_param(PF, mdl, "ROTATE_INDEX", rotate_index, &
571		"Enable rotation of the horizontal indices.", &
572	0	default=.false., debuggingParam=.true., do_not_log=.true.)
573	0	if (rotate_index) then
574		! This model is using a rotated grid, so the unrotated variables used here have not been set yet.
575	0	call copy_restart_var(h, "h", restart_CS, .true.)
576	0	call copy_restart_vector(u, v, "u", "v", restart_CS, .true.)
577	0	if ( use_temperature ) then
578	0	call copy_restart_var(tv%T, "Temp", restart_CS, .true.)
579	0	call copy_restart_var(tv%S, "Salt", restart_CS, .true.)
580		endif
581		endif
582		endif
583
584	1	if ( use_temperature ) then
585	1	call pass_var(tv%T, G%Domain, complete=.false.)
586	1	call pass_var(tv%S, G%Domain, complete=.false.)
587		endif
588	1	call pass_var(h, G%Domain)
589
590	1	if (debug) then
591	0	call hchksum(h, "MOM_initialize_state: h ", G%HI, haloshift=1, unscale=GV%H_to_MKS)
592	0	if ( use_temperature ) call hchksum(tv%T, "MOM_initialize_state: T ", G%HI, haloshift=1, unscale=US%C_to_degC)
593	0	if ( use_temperature ) call hchksum(tv%S, "MOM_initialize_state: S ", G%HI, haloshift=1, unscale=US%S_to_ppt)
594	0	if ( use_temperature .and. debug_layers) then ; do k=1,nz
595	0	write(mesg,'("MOM_IS: T[",I0,"]")') k
596	0	call hchksum(tv%T(:,:,k), mesg, G%HI, haloshift=1, unscale=US%C_to_degC)
597	0	write(mesg,'("MOM_IS: S[",I0,"]")') k
598	0	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	1	"specified via SPONGE_CONFIG.", default=.false.)
606	1	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	0	" \t USER - call a user modified routine.", default="file")
618	0	select case (trim(config))
619	0	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	0	sponge_CSp, ALE_sponge_CSp)
622		case ("ISOMIP"); call ISOMIP_initialize_sponges(G, GV, US, tv, depth_tot, PF, useALE, &
623	0	sponge_CSp, ALE_sponge_CSp)
624		case ("RGC"); call RGC_initialize_sponges(G, GV, US, tv, u, v, depth_tot, PF, useALE, &
625	0	sponge_CSp, ALE_sponge_CSp)
626	0	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	0	sponge_CSp, h)
629		case ("DUMBBELL"); call dumbbell_initialize_sponges(G, GV, US, tv, h, depth_tot, PF, useALE, &
630	0	sponge_CSp, ALE_sponge_CSp)
631	0	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	0	sponge_CSp, ALE_sponge_CSp)
634		case ("file"); call initialize_sponges_file(G, GV, US, use_temperature, tv, u, v, depth_tot, PF, &
635	0	sponge_CSp, ALE_sponge_CSp, Time)
636		case default ; call MOM_error(FATAL, "MOM_initialize_state: "//&
637	0	"Unrecognized sponge configuration "//trim(config))
638		end select
639		endif
640
641		! Set-up of data Assimilation with incremental update
642	1	if (use_oda_incupd) then
643		call initialize_oda_incupd_file(G, GV, US, use_temperature, tv, h, u, v, &
644	0	PF, oda_incupd_CSp, restart_CS, Time)
645		endif
646
647	1	call callTree_leave('MOM_initialize_state()')
648
649	1	end subroutine MOM_initialize_state
650
651	0	subroutine 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	0	call callTree_enter('MOM_initialize_OBCs()')
678	0	if (associated(OBC)) then
679	0	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	0	do_not_log=.true., old_name="DEBUG_OBC", debuggingParam=.true.)
684		call get_param(PF, mdl, "ENABLE_BUGS_BY_DEFAULT", enable_bugs, &
685	0	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	0	"specified for the reservoir concentrations.", default=enable_bugs)
690	0	if (associated(tv%T)) then
691	0	if (OBC_reservoir_init_bug) then
692	0	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	0	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	0	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	0	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	0	" USER - user specified", default="none")
717	0	if (trim(config) == "DOME") then
718	0	call DOME_set_OBC_data(OBC, tv, G, GV, US, PF, tracer_Reg)
719	0	elseif (trim(config) == "dyed_channel") then
720	0	call dyed_channel_set_OBC_tracer_data(OBC, G, GV, PF, tracer_Reg)
721	0	OBC%update_OBC = .true.
722	0	elseif (trim(config) == "dyed_obcs") then
723	0	call dyed_obcs_set_OBC_data(OBC, G, GV, PF, tracer_Reg)
724	0	elseif (trim(config) == "Kelvin") then
725	0	OBC%update_OBC = .true.
726	0	elseif (trim(config) == "shelfwave") then
727	0	OBC%update_OBC = .true.
728	0	elseif (lowercase(trim(config)) == "supercritical") then
729	0	call supercritical_set_OBC_data(OBC, G, GV, US, PF)
730	0	elseif (trim(config) == "tidal_bay") then
731	0	OBC%update_OBC = .true.
732	0	elseif (trim(config) == "USER") then
733	0	call user_set_OBC_data(OBC, tv, G, GV, PF, tracer_Reg)
734	0	elseif (.not. trim(config) == "none") then
735		call MOM_error(FATAL, "The open boundary conditions specified by "//&
736	0	"OBC_USER_CONFIG = "//trim(config)//" have not been fully implemented.")
737		endif
738
739	0	if (debug) then
740	0	call hchksum(G%mask2dT, 'MOM_initialize_OBCs: mask2dT ', G%HI)
741	0	call uvchksum('MOM_initialize_OBCs: mask2dC[uv]', G%mask2dCu, G%mask2dCv, G%HI)
742	0	call qchksum(G%mask2dBu, 'MOM_initialize_OBCs: mask2dBu ', G%HI)
743		endif
744	0	if (debug_OBC) call open_boundary_test_extern_h(G, GV, OBC, h)
745
746	0	if (OBC%use_h_res) &
747	0	call fill_thickness_segments(G, GV, US, OBC, h)
748		endif
749
750	0	call callTree_leave('MOM_initialize_OBCs()')
751
752	0	end subroutine MOM_initialize_OBCs
753
754		!> Reads the layer thicknesses or interface heights from a file.
755	0	subroutine 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	0	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	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
795
796	0	if (.not.just_read) &
797	0	call callTree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
798
799	0	call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".", do_not_log=just_read)
800	0	inputdir = slasher(inputdir)
801		call get_param(param_file, mdl, "THICKNESS_FILE", thickness_file, &
802		"The name of the thickness file.", &
803	0	fail_if_missing=.not.just_read, do_not_log=just_read)
804
805	0	filename = trim(thickness_file)
806	0	if (scan(thickness_file, "/") == 0) then ! prepend inputdir if only a filename is given
807	0	filename = trim(inputdir)//trim(thickness_file)
808		endif
809	0	if (.not.just_read) call log_param(param_file, mdl, "INPUTDIR/THICKNESS_FILE", filename)
810
811	0	if ((.not.just_read) .and. (.not.file_exists(filename, G%Domain))) call MOM_error(FATAL, &
812	0	" initialize_thickness_from_file: Unable to open "//trim(filename))
813
814	0	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	0	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	0	default=1.0, units="various", do_not_log=just_read)
822	0	if (just_read) return ! All run-time parameters have been read, so return.
823
824	0	if (mass_file) then
825	0	h_rescale = h_rescale*GV%kg_m2_to_H
826		else
827	0	h_rescale = h_rescale*US%m_to_Z
828		endif
829	0	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	0	"would indicate.", default=.false., do_not_log=just_read)
835	0	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	0	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	0	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	0	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	0	default=1.0, units="various", do_not_log=just_read)
854	0	if (just_read) return ! All run-time parameters have been read, so return.
855
856	0	call MOM_read_data(filename, eta_var, eta(:,:,:), G%Domain, scale=US%m_to_Z*eta_rescale)
857
858	0	if (correct_thickness) then
859	0	call adjustEtaToFitBathymetry(G, GV, US, eta, h, h_tolerance, dZ_ref_eta=G%Z_ref)
860		else
861	0	do k=nz,1,-1 ; do j=js,je ; do i=is,ie
862	0	if (eta(i,j,K) < (eta(i,j,K+1) + GV%Angstrom_Z)) then
863	0	eta(i,j,K) = eta(i,j,K+1) + GV%Angstrom_Z
864	0	h(i,j,k) = GV%Angstrom_Z
865		else
866	0	h(i,j,k) = eta(i,j,K) - eta(i,j,K+1)
867		endif
868		enddo ; enddo ; enddo
869
870	0	inconsistent = 0
871	0	do j=js,je ; do i=is,ie
872	0	if (abs(eta(i,j,nz+1) + depth_tot(i,j)) > tol_dz_bot) &
873	0	inconsistent = inconsistent + 1
874		enddo ; enddo
875	0	call sum_across_PEs(inconsistent)
876
877	0	if ((inconsistent > 0) .and. (is_root_pe())) then
878		write(mesg,'("Thickness initial conditions are inconsistent ",'// &
879	0	'"with topography in ",I0," places.")') inconsistent
880	0	call MOM_error(WARNING, mesg)
881		endif
882		endif
883
884		endif
885	0	call callTree_leave(trim(mdl)//'()')
886		end subroutine initialize_thickness_from_file
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.
894	0	subroutine 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	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
912	0	dZ_ref = 0.0 ; if (present(dZ_ref_eta)) dZ_ref = dZ_ref_eta
913
914	0	contractions = 0
915	0	do j=js,je ; do i=is,ie
916	0	if (-eta(i,j,nz+1) > (G%bathyT(i,j) + dZ_ref) + ht) then
917	0	eta(i,j,nz+1) = -(G%bathyT(i,j) + dZ_ref)
918	0	contractions = contractions + 1
919		endif
920		enddo ; enddo
921	0	call sum_across_PEs(contractions)
922	0	if ((contractions > 0) .and. (is_root_pe())) then
923		write(mesg,'("Thickness initial conditions were contracted ",'// &
924	0	'"to fit topography in ",I0," places.")') contractions
925	0	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	0	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	0	if (eta(i,j,K) < (eta(i,j,K+1) + GV%Angstrom_Z)) then
934	0	eta(i,j,K) = eta(i,j,K+1) + GV%Angstrom_Z
935	0	h(i,j,k) = GV%Angstrom_Z
936		else
937	0	h(i,j,k) = (eta(i,j,K) - eta(i,j,K+1))
938		endif
939		enddo ; enddo ; enddo
940
941	0	dilations = 0
942	0	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	0	if (-eta(i,j,nz+1) < (G%bathyT(i,j) + dZ_ref) - ht) then
947	0	dilations = dilations + 1
948	0	if (eta(i,j,1) <= eta(i,j,nz+1)) then
949	0	do k=1,nz ; h(i,j,k) = (eta(i,j,1) + (G%bathyT(i,j) + dZ_ref)) / real(nz) ; enddo
950		else
951	0	dilate = (eta(i,j,1) + (G%bathyT(i,j) + dZ_ref)) / (eta(i,j,1) - eta(i,j,nz+1))
952	0	do k=1,nz ; h(i,j,k) = h(i,j,k) * dilate ; enddo
953		endif
954	0	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	0	call sum_across_PEs(dilations)
960	0	if ((dilations > 0) .and. (is_root_pe())) then
961		write(mesg,'("Thickness initial conditions were dilated ",'// &
962	0	'"to fit topography in ",I0," places.")') dilations
963	0	call MOM_error(WARNING, 'adjustEtaToFitBathymetry: '//mesg)
964		endif
965
966	0	end subroutine adjustEtaToFitBathymetry
967
968		!> Initializes thickness to be uniform
969	1	subroutine 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	3	real :: e0(SZK_(GV)+1) ! The resting interface heights [Z ~> m], usually
983		! negative because it is positive upward.
984	1	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	1	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
989
990	1	if (just_read) return ! This subroutine has no run-time parameters.
991
992	1	call callTree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
993
994	1	if (G%max_depth<=0.) call MOM_error(FATAL,"initialize_thickness_uniform: "// &
995	0	"MAXIMUM_DEPTH has a nonsensical value! Was it set?")
996
997	76	do k=1,nz
998	76	e0(K) = -G%max_depth * real(k-1) / real(nz)
999		enddo
1000
1001	7261	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	7200	eta1D(nz+1) = -depth_tot(i,j)
1008	547260	do k=nz,1,-1
1009	540000	eta1D(K) = e0(K)
1010	547200	if (eta1D(K) < (eta1D(K+1) + GV%Angstrom_Z)) then
1011	353730	eta1D(K) = eta1D(K+1) + GV%Angstrom_Z
1012	353730	h(i,j,k) = GV%Angstrom_Z
1013		else
1014	186270	h(i,j,k) = eta1D(K) - eta1D(K+1)
1015		endif
1016		enddo
1017		enddo ; enddo
1018
1019	1	call callTree_leave(trim(mdl)//'()')
1020		end subroutine initialize_thickness_uniform
1021
1022		!> Initialize thickness from a 1D list
1023	0	subroutine 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	0	real :: e0(SZK_(GV)+1) ! The resting interface heights, in depth units [Z ~> m],
1038		! usually negative because it is positive upward.
1039	0	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	0	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	0	"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	0	default="eta")
1054
1055	0	if (just_read) return
1056
1057	0	call callTree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
1058
1059	0	call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
1060	0	filename = trim(slasher(inputdir))//trim(eta_file)
1061	0	call log_param(param_file, mdl, "INPUTDIR/INTERFACE_IC_FILE", filename)
1062
1063	0	e0(:) = 0.0
1064	0	call MOM_read_data(filename, eta_var, e0(:), scale=US%m_to_Z)
1065
1066	0	if ((abs(e0(1)) - 0.0) > 0.001) then
1067		! This list probably starts with the interior interface, so shift it up.
1068	0	do k=nz+1,2,-1 ; e0(K) = e0(K-1) ; enddo
1069	0	e0(1) = 0.0
1070		endif
1071
1072	0	if (e0(2) > e0(1)) then ! Switch to the convention for interface heights increasing upward.
1073	0	do k=1,nz ; e0(K) = -e0(K) ; enddo
1074		endif
1075
1076	0	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	0	eta1D(nz+1) = -depth_tot(i,j)
1083	0	do k=nz,1,-1
1084	0	eta1D(K) = e0(K)
1085	0	if (eta1D(K) < (eta1D(K+1) + GV%Angstrom_Z)) then
1086	0	eta1D(K) = eta1D(K+1) + GV%Angstrom_Z
1087	0	h(i,j,k) = GV%Angstrom_Z
1088		else
1089	0	h(i,j,k) = eta1D(K) - eta1D(K+1)
1090		endif
1091		enddo
1092		enddo ; enddo
1093
1094	0	call callTree_leave(trim(mdl)//'()')
1095		end subroutine initialize_thickness_list
1096
1097		!> Initializes thickness based on a run-time parameter with nominal thickness
1098		!! for each layer
1099	0	subroutine 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	0	real :: e0(SZK_(GV)+1) ! The resting interface heights [Z ~> m], usually
1114		! negative because it is positive upward.
1115	0	real :: eta1D(SZK_(GV)+1)! Interface height relative to the sea surface,
1116		! positive upward [Z ~> m].
1117	0	real :: dz(SZK_(GV)) ! The nominal initial layer thickness [Z ~> m], usually
1118	0	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	0	call callTree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
1122	0	if (G%max_depth<=0.) call MOM_error(FATAL, "initialize_thickness_param: "// &
1123	0	"MAXIMUM_DEPTH has a nonsensical value! Was it set?")
1124
1125	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
1126
1127	0	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	0	units="m", scale=US%m_to_Z, defaults=h0_def, do_not_log=just_read)
1131	0	if (just_read) return ! This subroutine has no run-time parameters.
1132
1133	0	e0(nz+1) = -G%max_depth
1134	0	do k=nz, 1, -1
1135	0	e0(K) = e0(K+1) + dz(k)
1136		enddo
1137
1138	0	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	0	eta1D(nz+1) = -depth_tot(i,j)
1145	0	do k=nz,1,-1
1146	0	eta1D(K) = e0(K)
1147	0	if (eta1D(K) < (eta1D(K+1) + GV%Angstrom_Z)) then
1148	0	eta1D(K) = eta1D(K+1) + GV%Angstrom_Z
1149	0	h(i,j,k) = GV%Angstrom_Z
1150		else
1151	0	h(i,j,k) = eta1D(K) - eta1D(K+1)
1152		endif
1153		enddo
1154		enddo ; enddo
1155
1156	0	call callTree_leave(trim(mdl)//'()')
1157		end subroutine initialize_thickness_param
1158
1159		!> Search density space for location of layers (not implemented!)
1160	0	subroutine initialize_thickness_search
1161	0	call MOM_error(FATAL," MOM_state_initialization.F90, initialize_thickness_search: NOT IMPLEMENTED")
1162	0	end subroutine initialize_thickness_search
1163
1164		!> Depress the sea-surface based on an initial condition file
1165	0	subroutine 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	0	eta_sfc ! The free surface height that the model should use [Z ~> m].
1180		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: &
1181	0	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	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
1192
1193	0	use_z_shelf = present(z_top_shelf)
1194
1195
1196	0	if (.not. use_z_shelf) then
1197		! Read the surface height (or pressure) from a file.
1198
1199	0	call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
1200	0	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	0	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	0	default="SSH", do_not_log=just_read)
1207	0	filename = trim(inputdir)//trim(eta_srf_file)
1208	0	if (.not.just_read) &
1209	0	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	0	units="variable", default=1.0, scale=US%m_to_Z, do_not_log=just_read)
1214
1215	0	if (just_read) return ! All run-time parameters have been read, so return.
1216
1217	0	call MOM_read_data(filename, eta_srf_var, eta_sfc, G%Domain, scale=scale_factor)
1218		else
1219	0	do j=js,je ; do i=is,ie
1220	0	eta_sfc(i,j) = z_top_shelf(i,j)
1221		enddo ; enddo
1222		endif
1223
1224		! Convert thicknesses to interface heights.
1225		!$omp target update from(h)
1226		!$omp target enter data map(alloc: eta)
1227	0	call find_eta(h, tv, G, GV, US, eta, dZref=G%Z_ref)
1228		!$omp target exit data map(from: eta)
1229
1230	0	do j=js,je ; do i=is,ie ; if (G%mask2dT(i,j) > 0.0) then
1231		! if (eta_sfc(i,j) < eta(i,j,nz+1)) then
1232		! Issue a warning?
1233		! endif
1234	0	if (eta_sfc(i,j) > eta(i,j,1)) then
1235		! Dilate the water column to agree, but only up to 10-fold.
1236	0	if (eta_sfc(i,j) - eta(i,j,nz+1) > 10.0*(eta(i,j,1) - eta(i,j,nz+1))) then
1237	0	dilate = 10.0
1238		call MOM_error(WARNING, "Free surface height dilation attempted "//&
1239	0	"to exceed 10-fold.", all_print=.true.)
1240		else
1241	0	dilate = (eta_sfc(i,j) - eta(i,j,nz+1)) / (eta(i,j,1) - eta(i,j,nz+1))
1242		endif
1243	0	do k=1,nz ; h(i,j,k) = h(i,j,k) * dilate ; enddo
1244	0	elseif (eta(i,j,1) > eta_sfc(i,j)) then
1245		! Remove any mass that is above the target free surface.
1246	0	do k=1,nz
1247	0	if (eta(i,j,K) <= eta_sfc(i,j)) exit
1248	0	if (eta(i,j,K+1) >= eta_sfc(i,j)) then
1249	0	h(i,j,k) = GV%Angstrom_H
1250		else
1251		h(i,j,k) = max(GV%Angstrom_H, h(i,j,k) * &
1252	0	(eta_sfc(i,j) - eta(i,j,K+1)) / (eta(i,j,K) - eta(i,j,K+1)) )
1253		endif
1254		enddo
1255		endif
1256		endif ; enddo ; enddo
1257
1258	0	end subroutine depress_surface
1259
1260		!> Adjust the layer thicknesses by cutting away the top of each model column at the depth
1261		!! where the hydrostatic pressure matches an imposed surface pressure read from file.
1262	0	subroutine trim_for_ice(PF, G, GV, US, ALE_CSp, tv, h, just_read)
1263		type(param_file_type), intent(in) :: PF !< Parameter file structure
1264		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1265		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
1266		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1267		type(ALE_CS), pointer :: ALE_CSp !< ALE control structure
1268		type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamics structure
1269		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1270		intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
1271		logical, intent(in) :: just_read !< If true, this call will only read
1272		!! parameters without changing h.
1273		! Local variables
1274		character(len=200) :: mdl = "trim_for_ice"
1275	0	real, dimension(SZI_(G),SZJ_(G)) :: p_surf ! Imposed pressure on ocean at surface [R L2 T-2 ~> Pa]
1276	0	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: S_t, S_b ! Top and bottom edge values for reconstructions
1277		! of salinity within each layer [S ~> ppt]
1278	0	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: T_t, T_b ! Top and bottom edge values for reconstructions
1279		! of temperature within each layer [C ~> degC]
1280		character(len=200) :: inputdir, filename, p_surf_file, p_surf_var ! Strings for file/path
1281		real :: scale_factor ! A file-dependent scaling factor for the input pressure [various].
1282		real :: min_thickness ! The minimum layer thickness [H ~> m or kg m-2].
1283		real :: z_tolerance ! The tolerance with which to find the depth matching a specified pressure [Z ~> m].
1284		integer :: i, j, k
1285		integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.
1286		integer :: remap_answer_date ! The vintage of the order of arithmetic and expressions to use
1287		! for remapping. Values below 20190101 recover the remapping
1288		! answers from 2018, while higher values use more robust
1289		! forms of the same remapping expressions.
1290		logical :: use_remapping ! If true, remap the initial conditions.
1291		logical :: use_frac_dp_bugfix ! If true, use bugfix. Otherwise, pressure input to EOS is negative.
1292		type(remapping_CS), pointer :: remap_CS => NULL()
1293
1294		call get_param(PF, mdl, "SURFACE_PRESSURE_FILE", p_surf_file, &
1295		"The initial condition file for the surface pressure exerted by ice.", &
1296	0	fail_if_missing=.not.just_read, do_not_log=just_read)
1297		call get_param(PF, mdl, "SURFACE_PRESSURE_VAR", p_surf_var, &
1298		"The initial condition variable for the surface pressure exerted by ice.", &
1299	0	default="", do_not_log=just_read)
1300	0	call get_param(PF, mdl, "INPUTDIR", inputdir, default=".", do_not_log=.true.)
1301	0	filename = trim(slasher(inputdir))//trim(p_surf_file)
1302	0	if (.not.just_read) call log_param(PF, mdl, "!INPUTDIR/SURFACE_HEIGHT_IC_FILE", filename)
1303
1304		call get_param(PF, mdl, "SURFACE_PRESSURE_SCALE", scale_factor, &
1305		"A scaling factor to convert SURFACE_PRESSURE_VAR from "//&
1306		"file SURFACE_PRESSURE_FILE into a surface pressure.", &
1307	0	units="file dependent", default=1., do_not_log=just_read)
1308		call get_param(PF, mdl, "MIN_THICKNESS", min_thickness, 'Minimum layer thickness', &
1309	0	units='m', default=1.e-3, scale=GV%m_to_H, do_not_log=just_read)
1310		call get_param(PF, mdl, "TRIM_IC_Z_TOLERANCE", z_tolerance, &
1311		"The tolerance with which to find the depth matching the specified "//&
1312		"surface pressure with TRIM_IC_FOR_P_SURF.", &
1313	0	units="m", default=1.0e-5, scale=US%m_to_Z, do_not_log=just_read)
1314		call get_param(PF, mdl, "FRAC_DP_AT_POS_NEGATIVE_P_BUGFIX", use_frac_dp_bugfix, &
1315		"If true, use bugfix in ice shelf TRIM_IC initialization. "//&
1316		"Otherwise, pressure input to density EOS is negative.", &
1317	0	default=.false., do_not_log=just_read)
1318		call get_param(PF, mdl, "TRIMMING_USES_REMAPPING", use_remapping, &
1319		'When trimming the column, also remap T and S.', &
1320	0	default=.false., do_not_log=just_read)
1321	0	if (use_remapping) then
1322		call get_param(PF, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
1323		"This sets the default value for the various _ANSWER_DATE parameters.", &
1324	0	default=99991231, do_not_log=just_read)
1325		call get_param(PF, mdl, "REMAPPING_ANSWER_DATE", remap_answer_date, &
1326		"The vintage of the expressions and order of arithmetic to use for remapping. "//&
1327		"Values below 20190101 result in the use of older, less accurate expressions "//&
1328		"that were in use at the end of 2018. Higher values result in the use of more "//&
1329		"robust and accurate forms of mathematically equivalent expressions.", &
1330	0	default=default_answer_date, do_not_log=just_read.or.(.not.GV%Boussinesq))
1331	0	if (.not.GV%Boussinesq) remap_answer_date = max(remap_answer_date, 20230701)
1332		else
1333	0	remap_answer_date = 20181231
1334	0	if (.not.GV%Boussinesq) remap_answer_date = 20230701
1335		endif
1336
1337	0	if (just_read) return ! All run-time parameters have been read, so return.
1338
1339		call MOM_read_data(filename, p_surf_var, p_surf, G%Domain, &
1340	0	scale=scale_factor*US%Pa_to_RL2_T2)
1341
1342	0	if (use_remapping) then
1343	0	allocate(remap_CS)
1344	0	if (remap_answer_date < 20190101) then
1345		call initialize_remapping(remap_CS, 'PLM', boundary_extrapolation=.true., &
1346	0	h_neglect=1.0e-30GV%m_to_H, h_neglect_edge=1.0e-10GV%m_to_H)
1347		else
1348		call initialize_remapping(remap_CS, 'PLM', boundary_extrapolation=.true., &
1349	0	h_neglect=GV%H_subroundoff, h_neglect_edge=GV%H_subroundoff)
1350		endif
1351		endif
1352
1353		! Find edge values of T and S used in reconstructions
1354	0	if ( associated(ALE_CSp) ) then ! This should only be associated if we are in ALE mode
1355	0	call TS_PLM_edge_values(ALE_CSp, S_t, S_b, T_t, T_b, G, GV, tv, h, .true.)
1356		else
1357		! call MOM_error(FATAL, "trim_for_ice: Does not work without ALE mode")
1358	0	do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec
1359	0	T_t(i,j,k) = tv%T(i,j,k) ; T_b(i,j,k) = tv%T(i,j,k)
1360	0	S_t(i,j,k) = tv%S(i,j,k) ; S_b(i,j,k) = tv%S(i,j,k)
1361		enddo ; enddo ; enddo
1362		endif
1363
1364	0	do j=G%jsc,G%jec ; do i=G%isc,G%iec
1365		call cut_off_column_top(GV%ke, tv, GV, US, GV%g_Earth, G%bathyT(i,j)+G%Z_ref, min_thickness, &
1366		tv%T(i,j,:), T_t(i,j,:), T_b(i,j,:), tv%S(i,j,:), S_t(i,j,:), S_b(i,j,:), &
1367		p_surf(i,j), h(i,j,:), remap_CS, z_tol=z_tolerance, &
1368	0	frac_dp_bugfix=use_frac_dp_bugfix)
1369		enddo ; enddo
1370
1371		end subroutine trim_for_ice
1372
1373		!> Calculate the hydrostatic equilibrium position of the surface under an ice shelf
1374	0	subroutine calc_sfc_displacement(PF, G, GV, US, mass_shelf, tv, h)
1375		type(param_file_type), intent(in) :: PF !< Parameter file structure
1376		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1377		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
1378		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1379		real, dimension(SZI_(G),SZJ_(G)), &
1380		intent(in) :: mass_shelf !< Ice shelf mass [R Z ~> kg m-2]
1381		type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamics structure
1382		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1383		intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
1384
1385	0	real :: z_top_shelf(SZI_(G),SZJ_(G)) ! The depth of the top interface under ice shelves [Z ~> m]
1386		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: &
1387	0	eta ! The free surface height that the model should use [Z ~> m].
1388		! temporary arrays
1389	0	real, dimension(SZK_(GV)) :: rho_col ! potential density in the column for use in ice [R ~> kg m-3]
1390	0	real, dimension(SZK_(GV)) :: rho_h ! potential density multiplied by thickness [R Z ~> kg m-2]
1391	0	real, dimension(SZK_(GV)) :: h_tmp ! temporary storage for thicknesses [H ~> m]
1392	0	real, dimension(SZK_(GV)) :: p_ref ! pressure for density [R Z ~> kg m-2]
1393	0	real, dimension(SZK_(GV)+1) :: ei_tmp, ei_orig ! temporary storage for interface positions [Z ~> m]
1394		real :: z_top ! An estimate of the height of the ice-ocean interface [Z ~> m]
1395		real :: mass_disp ! The net mass of sea water that has been displaced by the shelf [R Z ~> kg m-2]
1396		real :: residual ! The difference between the displaced ocean mass and the ice shelf
1397		! mass [R Z ~> kg m-2]
1398		real :: tol ! The initialization tolerance for ice shelf initialization [Z ~> m]
1399		integer :: is, ie, js, je, k, nz, i, j, max_iter, iter
1400
1401	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
1402
1403		call get_param(PF, mdl, "ICE_SHELF_INITIALIZATION_Z_TOLERANCE", tol, &
1404		"A initialization tolerance for the calculation of the static "// &
1405		"ice shelf displacement (m) using initial temperature and salinity profile.", &
1406	0	default=0.001, units="m", scale=US%m_to_Z)
1407	0	max_iter = 1e3
1408	0	call MOM_mesg("Started calculating initial interface position under ice shelf ")
1409		! Convert thicknesses to interface heights.
1410		!$omp target update to(h)
1411		!$omp target enter data map(alloc: eta)
1412	0	call find_eta(h, tv, G, GV, US, eta, dZref=G%Z_ref)
1413		!$omp target exit data map(from: eta)
1414	0	do j=js,je ; do i=is,ie
1415	0	iter = 1
1416	0	z_top_shelf(i,j) = 0.0
1417	0	p_ref(:) = tv%p_ref
1418	0	if ((G%mask2dT(i,j) > 0.) .and. (mass_shelf(i,j) > 0.)) then
1419	0	call calculate_density(tv%T(i,j,:), tv%S(i,j,:), P_Ref, rho_col, tv%eqn_of_state)
1420	0	z_top = min(max(-1.0*mass_shelf(i,j)/rho_col(1), -G%bathyT(i,j)), 0.)
1421	0	h_tmp(:) = 0.0
1422	0	ei_tmp(1:nz+1) = eta(i,j,1:nz+1)
1423	0	ei_orig(1:nz+1) = eta(i,j,1:nz+1)
1424	0	do k=1,nz+1
1425	0	if (ei_tmp(k) < z_top) ei_tmp(k) = z_top
1426		enddo
1427	0	mass_disp = 0.0
1428	0	do k=1,nz
1429	0	h_tmp(k) = max(ei_tmp(k)-ei_tmp(k+1), GV%Angstrom_H)
1430	0	rho_h(k) = h_tmp(k) * rho_col(k)
1431	0	mass_disp = mass_disp + rho_h(k)
1432		enddo
1433	0	residual = mass_shelf(i,j) - mass_disp
1434	0	do while ((abs(residual) > tol) .and. (z_top > -G%bathyT(i,j)) .and. (iter < max_iter))
1435	0	z_top = min(max(z_top-(residual*0.5e-3), -G%bathyT(i,j)), 0.0)
1436	0	h_tmp(:) = 0.0
1437	0	ei_tmp(1:nz+1) = ei_orig(1:nz+1)
1438	0	do k=1,nz+1
1439	0	if (ei_tmp(k) < z_top) ei_tmp(k) = z_top
1440		enddo
1441	0	mass_disp = 0.0
1442	0	do k=1,nz
1443	0	h_tmp(k) = max(ei_tmp(k)-ei_tmp(k+1), GV%Angstrom_H)
1444	0	rho_h(k) = h_tmp(k) * rho_col(k)
1445	0	mass_disp = mass_disp + rho_h(k)
1446		enddo
1447	0	residual = mass_shelf(i,j) - mass_disp
1448	0	iter = iter+1
1449		enddo
1450	0	if (iter >= max_iter) call MOM_mesg("Warning: calc_sfc_displacement too many iterations.")
1451	0	z_top_shelf(i,j) = z_top
1452		endif
1453		enddo ; enddo
1454	0	call MOM_mesg("Calling depress_surface ")
1455	0	call depress_surface(h, G, GV, US, PF, tv, just_read=.false.,z_top_shelf=z_top_shelf)
1456	0	call MOM_mesg("Finishing calling depress_surface ")
1457	0	end subroutine calc_sfc_displacement
1458
1459		!> Adjust the layer thicknesses by removing the top of the water column above the
1460		!! depth where the hydrostatic pressure matches p_surf
1461	0	subroutine cut_off_column_top(nk, tv, GV, US, G_earth, depth, min_thickness, T, T_t, T_b, &
1462	0	S, S_t, S_b, p_surf, h, remap_CS, z_tol, frac_dp_bugfix)
1463		integer, intent(in) :: nk !< Number of layers
1464		type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure
1465		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1466		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1467		real, intent(in) :: G_earth !< Gravitational acceleration [L2 Z-1 T-2 ~> m s-2]
1468		real, intent(in) :: depth !< Depth of ocean column [Z ~> m].
1469		real, intent(in) :: min_thickness !< Smallest thickness allowed [H ~> m or kg m-2].
1470		real, dimension(nk), intent(inout) :: T !< Layer mean temperature [C ~> degC]
1471		real, dimension(nk), intent(in) :: T_t !< Temperature at top of layer [C ~> degC]
1472		real, dimension(nk), intent(in) :: T_b !< Temperature at bottom of layer [C ~> degC]
1473		real, dimension(nk), intent(inout) :: S !< Layer mean salinity [S ~> ppt]
1474		real, dimension(nk), intent(in) :: S_t !< Salinity at top of layer [S ~> ppt]
1475		real, dimension(nk), intent(in) :: S_b !< Salinity at bottom of layer [S ~> ppt]
1476		real, intent(in) :: p_surf !< Imposed pressure on ocean at surface [R L2 T-2 ~> Pa]
1477		real, dimension(nk), intent(inout) :: h !< Layer thickness [H ~> m or kg m-2]
1478		type(remapping_CS), pointer :: remap_CS !< Remapping structure for remapping T and S,
1479		!! if associated
1480		real, intent(in) :: z_tol !< The tolerance with which to find the depth
1481		!! matching the specified pressure [Z ~> m].
1482		logical, intent(in) :: frac_dp_bugfix !< If true, use bugfix in frac_dp_at_pos
1483
1484		! Local variables
1485	0	real, dimension(nk+1) :: e ! Top and bottom edge positions for reconstructions [Z ~> m]
1486	0	real, dimension(nk) :: h0, h1 ! Initial and remapped layer thicknesses [H ~> m or kg m-2]
1487	0	real, dimension(nk) :: S0, S1 ! Initial and remapped layer salinities [S ~> ppt]
1488	0	real, dimension(nk) :: T0, T1 ! Initial and remapped layer temperatures [C ~> degC]
1489		real :: P_t, P_b ! Top and bottom pressures [R L2 T-2 ~> Pa]
1490		real :: z_out, e_top ! Interface height positions [Z ~> m]
1491		real :: min_dz ! The minimum thickness in depth units [Z ~> m]
1492		real :: dh_surf_rem ! The remaining thickness to remove in non-Bousinesq mode [H ~> kg m-2]
1493		integer :: k
1494
1495		! Keep a copy of the initial thicknesses in reverse order to use in remapping
1496	0	do k=1,nk ; h0(k) = h(nk+1-k) ; enddo
1497
1498	0	if (GV%Boussinesq) then
1499	0	min_dz = GV%H_to_Z * min_thickness
1500		! Calculate original interface positions
1501	0	e(nk+1) = -depth
1502	0	do k=nk,1,-1
1503	0	e(K) = e(K+1) + GV%H_to_Z*h(k)
1504		enddo
1505
1506	0	P_t = 0.
1507	0	e_top = e(1)
1508	0	do k=1,nk
1509		call find_depth_of_pressure_in_cell(T_t(k), T_b(k), S_t(k), S_b(k), e(K), e(K+1), &
1510		P_t, p_surf, GV%Rho0, G_earth, tv%eqn_of_state, &
1511		US, P_b, z_out, z_tol=z_tol, &
1512	0	frac_dp_bugfix=frac_dp_bugfix)
1513	0	if (z_out>=e(K)) then
1514		! Imposed pressure was less that pressure at top of cell
1515	0	exit
1516	0	elseif (z_out<=e(K+1)) then
1517		! Imposed pressure was greater than pressure at bottom of cell
1518	0	e_top = e(K+1)
1519		else
1520		! Imposed pressure was fell between pressures at top and bottom of cell
1521	0	e_top = z_out
1522	0	exit
1523		endif
1524	0	P_t = P_b
1525		enddo
1526	0	if (e_top<e(1)) then
1527		! Clip layers from the top down, if at all
1528	0	do K=1,nk
1529	0	if (e(K) > e_top) then
1530		! Original e(K) is too high
1531	0	e(K) = e_top
1532	0	e_top = e_top - min_dz ! Next interface must be at least this deep
1533		endif
1534		! This layer needs trimming
1535	0	h(k) = max( min_thickness, GV%Z_to_H * (e(K) - e(K+1)) )
1536	0	if (e(K) < e_top) exit ! No need to go further
1537		enddo
1538		endif
1539		else
1540		! In non-Bousinesq mode, we are already in mass units so the calculation is much easier.
1541	0	if (p_surf > 0.0) then
1542	0	dh_surf_rem = p_surf * GV%RZ_to_H / G_earth
1543	0	do k=1,nk
1544	0	if (h(k) <= min_thickness) then ! This layer has no mass to remove.
1545	0	cycle
1546	0	elseif ((h(k) - min_thickness) < dh_surf_rem) then ! This layer should be removed entirely.
1547	0	dh_surf_rem = dh_surf_rem - (h(k) - min_thickness)
1548	0	h(k) = min_thickness
1549		else ! This is the last layer that should be removed.
1550	0	h(k) = h(k) - dh_surf_rem
1551	0	dh_surf_rem = 0.0
1552	0	exit
1553		endif
1554		enddo
1555		endif
1556		endif
1557
1558		! Now we need to remap but remapping assumes the surface is at the
1559		! same place in the two columns so we turn the column upside down.
1560	0	if (associated(remap_CS)) then
1561	0	do k=1,nk
1562	0	S0(k) = S(nk+1-k)
1563	0	T0(k) = T(nk+1-k)
1564	0	h1(k) = h(nk+1-k)
1565		enddo
1566	0	call remapping_core_h(remap_CS, nk, h0, T0, nk, h1, T1)
1567	0	call remapping_core_h(remap_CS, nk, h0, S0, nk, h1, S1)
1568	0	do k=1,nk
1569	0	S(k) = S1(nk+1-k)
1570	0	T(k) = T1(nk+1-k)
1571		enddo
1572		endif
1573
1574	0	end subroutine cut_off_column_top
1575
1576		!> Initialize horizontal velocity components from file
1577	0	subroutine initialize_velocity_from_file(u, v, G, GV, US, param_file, just_read)
1578		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1579		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1580		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1581		intent(out) :: u !< The zonal velocity that is being initialized [L T-1 ~> m s-1]
1582		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1583		intent(out) :: v !< The meridional velocity that is being initialized [L T-1 ~> m s-1]
1584		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1585		type(param_file_type), intent(in) :: param_file !< A structure indicating the open file to
1586		!! parse for model parameter values.
1587		logical, intent(in) :: just_read !< If true, this call will only read
1588		!! parameters without changing u or v.
1589		! Local variables
1590		character(len=40) :: mdl = "initialize_velocity_from_file" ! This subroutine's name.
1591		character(len=200) :: filename, velocity_file, inputdir ! Strings for file/path
1592		character(len=64) :: u_IC_var, v_IC_var ! Velocity component names in files
1593
1594	0	if (.not.just_read) call callTree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
1595
1596		call get_param(param_file, mdl, "VELOCITY_FILE", velocity_file, &
1597		"The name of the velocity initial condition file.", &
1598	0	fail_if_missing=.not.just_read, do_not_log=just_read)
1599	0	call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
1600	0	inputdir = slasher(inputdir)
1601
1602	0	filename = trim(velocity_file)
1603	0	if (scan(velocity_file, '/')== 0) then ! prepend inputdir if only a filename is given
1604	0	filename = trim(inputdir)//trim(velocity_file)
1605		endif
1606	0	if (.not.just_read) call log_param(param_file, mdl, "INPUTDIR/VELOCITY_FILE", filename)
1607
1608		call get_param(param_file, mdl, "U_IC_VAR", u_IC_var, &
1609		"The initial condition variable for zonal velocity in VELOCITY_FILE.", &
1610	0	default="u")
1611		call get_param(param_file, mdl, "V_IC_VAR", v_IC_var, &
1612		"The initial condition variable for meridional velocity in VELOCITY_FILE.", &
1613	0	default="v")
1614
1615	0	if (just_read) return ! All run-time parameters have been read, so return.
1616
1617	0	if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, &
1618	0	" initialize_velocity_from_file: Unable to open "//trim(filename))
1619
1620		! Read the velocities from a netcdf file.
1621	0	call MOM_read_vector(filename, u_IC_var, v_IC_var, u(:,:,:), v(:,:,:), G%Domain, scale=US%m_s_to_L_T)
1622
1623	0	call callTree_leave(trim(mdl)//'()')
1624		end subroutine initialize_velocity_from_file
1625
1626		!> Initialize horizontal velocity components to zero.
1627	1	subroutine initialize_velocity_zero(u, v, G, GV, param_file, just_read)
1628		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1629		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1630		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1631		intent(out) :: u !< The zonal velocity that is being initialized [L T-1 ~> m s-1]
1632		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1633		intent(out) :: v !< The meridional velocity that is being initialized [L T-1 ~> m s-1]
1634		type(param_file_type), intent(in) :: param_file !< A structure indicating the open file to
1635		!! parse for model parameter values.
1636		logical, intent(in) :: just_read !< If true, this call will only read
1637		!! parameters without changing h.
1638		! Local variables
1639		character(len=200) :: mdl = "initialize_velocity_zero" ! This subroutine's name.
1640		integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
1641
1642	1	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
1643	1	Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
1644
1645	1	if (.not.just_read) call callTree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
1646
1647	1	if (just_read) return ! All run-time parameters have been read, so return.
1648
1649	549076	do k=1,nz ; do j=js,je ; do I=Isq,Ieq
1650	549000	u(I,j,k) = 0.0
1651		enddo ; enddo ; enddo
1652	553651	do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
1653	553575	v(i,J,k) = 0.0
1654		enddo ; enddo ; enddo
1655
1656	1	call callTree_leave(trim(mdl)//'()')
1657		end subroutine initialize_velocity_zero
1658
1659		!> Sets the initial velocity components to uniform
1660	0	subroutine initialize_velocity_uniform(u, v, G, GV, US, param_file, just_read)
1661		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1662		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1663		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1664		intent(out) :: u !< The zonal velocity that is being initialized [L T-1 ~> m s-1]
1665		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1666		intent(out) :: v !< The meridional velocity that is being initialized [L T-1 ~> m s-1]
1667		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1668		type(param_file_type), intent(in) :: param_file !< A structure indicating the open file to
1669		!! parse for model parameter values.
1670		logical, intent(in) :: just_read !< If true, this call will only read
1671		!! parameters without changing u or v.
1672		! Local variables
1673		integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
1674		real :: initial_u_const, initial_v_const ! Constant initial velocities [L T-1 ~> m s-1]
1675		character(len=200) :: mdl = "initialize_velocity_uniform" ! This subroutine's name.
1676
1677	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
1678	0	Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
1679
1680		call get_param(param_file, mdl, "INITIAL_U_CONST", initial_u_const, &
1681		"A initial uniform value for the zonal flow.", &
1682	0	default=0.0, units="m s-1", scale=US%m_s_to_L_T, do_not_log=just_read)
1683		call get_param(param_file, mdl, "INITIAL_V_CONST", initial_v_const, &
1684		"A initial uniform value for the meridional flow.", &
1685	0	default=0.0, units="m s-1", scale=US%m_s_to_L_T, do_not_log=just_read)
1686
1687	0	if (just_read) return ! All run-time parameters have been read, so return.
1688
1689	0	do k=1,nz ; do j=js,je ; do I=Isq,Ieq
1690	0	u(I,j,k) = initial_u_const
1691		enddo ; enddo ; enddo
1692	0	do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
1693	0	v(i,J,k) = initial_v_const
1694		enddo ; enddo ; enddo
1695
1696		end subroutine initialize_velocity_uniform
1697
1698		!> Sets the initial velocity components to be circular with
1699		!! no flow at edges of domain and center.
1700	0	subroutine initialize_velocity_circular(u, v, G, GV, US, param_file, just_read)
1701		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1702		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1703		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1704		intent(out) :: u !< The zonal velocity that is being initialized [L T-1 ~> m s-1]
1705		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1706		intent(out) :: v !< The meridional velocity that is being initialized [L T-1 ~> m s-1]
1707		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1708		type(param_file_type), intent(in) :: param_file !< A structure indicating the open file to
1709		!! parse for model parameter values.
1710		logical, intent(in) :: just_read !< If true, this call will only read
1711		!! parameters without changing u or v.
1712		! Local variables
1713		character(len=200) :: mdl = "initialize_velocity_circular"
1714		real :: circular_max_u ! The amplitude of the zonal flow [L T-1 ~> m s-1]
1715		real :: dpi ! A local variable storing pi = 3.14159265358979... [nondim]
1716		real :: psi1, psi2 ! Values of the streamfunction at two points [L2 T-1 ~> m2 s-1]
1717		integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
1718	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
1719	0	Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB
1720
1721		call get_param(param_file, mdl, "CIRCULAR_MAX_U", circular_max_u, &
1722		"The amplitude of zonal flow from which to scale the "// &
1723		"circular stream function [m s-1].", &
1724	0	units="m s-1", default=0., scale=US%m_s_to_L_T, do_not_log=just_read)
1725
1726	0	if (just_read) return ! All run-time parameters have been read, so return.
1727
1728	0	if (G%grid_unit_to_L <= 0.) call MOM_error(FATAL, "MOM_state_initialization.F90: "//&
1729	0	"initialize_velocity_circular() is only set to work with Cartesian axis units.")
1730
1731	0	dpi = acos(0.0)*2.0 ! pi
1732
1733	0	do k=1,nz ; do j=js,je ; do I=Isq,Ieq
1734	0	psi1 = my_psi(I,j)
1735	0	psi2 = my_psi(I,j-1)
1736	0	u(I,j,k) = (psi1 - psi2) / G%dy_Cu(I,j) ! (circular_max_uG%len_lon/(2.0*dpi))
1737		enddo ; enddo ; enddo
1738	0	do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
1739	0	psi1 = my_psi(i,J)
1740	0	psi2 = my_psi(i-1,J)
1741	0	v(i,J,k) = (psi2 - psi1) / G%dx_Cv(i,J) ! (circular_max_uG%len_lon/(2.0*dpi))
1742		enddo ; enddo ; enddo
1743
1744		contains
1745
1746		!> Returns the value of a circular stream function at (ig,jg) in [L2 T-1 ~> m2 s-1]
1747	0	real function my_psi(ig,jg)
1748		integer :: ig !< Global i-index
1749		integer :: jg !< Global j-index
1750		! Local variables
1751		real :: x, y, r ! [nondim]
1752
1753	0	x = 2.0*(G%geoLonBu(ig,jg)-G%west_lon) / G%len_lon - 1.0 ! -1<x<1
1754	0	y = 2.0*(G%geoLatBu(ig,jg)-G%south_lat) / G%len_lat - 1.0 ! -1<y<1
1755	0	r = sqrt( (x2) + (y2) ) ! Circular stream function is a function of radius only
1756	0	r = min(1.0, r) ! Flatten stream function in corners of box
1757	0	my_psi = 0.5(1.0 - cos(dpir))
1758	0	my_psi = my_psi * (circular_max_u * G%len_lon * G%grid_unit_to_L / dpi) ! len_lon is in km
1759	0	end function my_psi
1760
1761		end subroutine initialize_velocity_circular
1762
1763		!> Initializes temperature and salinity from file
1764	0	subroutine initialize_temp_salt_from_file(T, S, G, GV, US, param_file, just_read)
1765		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1766		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1767		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T !< The potential temperature that is
1768		!! being initialized [C ~> degC]
1769		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S !< The salinity that is
1770		!! being initialized [S ~> ppt]
1771		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1772		type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
1773		logical, intent(in) :: just_read !< If true, this call will only
1774		!! read parameters without changing T or S.
1775		! Local variables
1776		character(len=200) :: filename, salt_filename ! Full paths to input files
1777		character(len=200) :: ts_file, salt_file, inputdir ! Strings for file/path
1778		character(len=40) :: mdl = "initialize_temp_salt_from_file"
1779		character(len=64) :: temp_var, salt_var ! Temperature and salinity names in files
1780
1781	0	if (.not.just_read) call callTree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
1782
1783		call get_param(param_file, mdl, "TS_FILE", ts_file, &
1784		"The initial condition file for temperature.", &
1785	0	fail_if_missing=.not.just_read, do_not_log=just_read)
1786	0	call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
1787	0	inputdir = slasher(inputdir)
1788
1789	0	filename = trim(ts_file)
1790	0	if (scan(ts_file, '/')== 0) then ! prepend inputdir if only a filename is given
1791	0	filename = trim(inputdir)//trim(ts_file)
1792		endif
1793	0	if (.not.just_read) call log_param(param_file, mdl, "INPUTDIR/TS_FILE", filename)
1794		call get_param(param_file, mdl, "TEMP_IC_VAR", temp_var, &
1795		"The initial condition variable for potential temperature.", &
1796	0	default="PTEMP", do_not_log=just_read)
1797		call get_param(param_file, mdl, "SALT_IC_VAR", salt_var, &
1798		"The initial condition variable for salinity.", &
1799	0	default="SALT", do_not_log=just_read)
1800		call get_param(param_file, mdl, "SALT_FILE", salt_file, &
1801		"The initial condition file for salinity.", &
1802	0	default=trim(ts_file), do_not_log=just_read)
1803
1804	0	if (just_read) return ! All run-time parameters have been read, so return.
1805
1806	0	if (.not.file_exists(filename, G%Domain)) call MOM_error(FATAL, &
1807	0	" initialize_temp_salt_from_file: Unable to open "//trim(filename))
1808
1809		! Read the temperatures and salinities from netcdf files.
1810	0	call MOM_read_data(filename, temp_var, T(:,:,:), G%Domain, scale=US%degC_to_C)
1811
1812	0	salt_filename = trim(salt_file)
1813	0	if (scan(salt_file, '/')== 0) then ! prepend inputdir if only a filename is given
1814	0	salt_filename = trim(inputdir)//trim(salt_file)
1815		endif
1816	0	if (.not.file_exists(salt_filename, G%Domain)) call MOM_error(FATAL, &
1817	0	" initialize_temp_salt_from_file: Unable to open "//trim(salt_filename))
1818
1819	0	call MOM_read_data(salt_filename, salt_var, S(:,:,:), G%Domain, scale=US%ppt_to_S)
1820
1821	0	call callTree_leave(trim(mdl)//'()')
1822		end subroutine initialize_temp_salt_from_file
1823
1824		!> Initializes temperature and salinity from a 1D profile
1825	0	subroutine initialize_temp_salt_from_profile(T, S, G, GV, US, param_file, just_read)
1826		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1827		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1828		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T !< The potential temperature that is
1829		!! being initialized [C ~> degC]
1830		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S !< The salinity that is
1831		!! being initialized [S ~> ppt]
1832		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1833		type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
1834		logical, intent(in) :: just_read !< If true, this call will only read
1835		!! parameters without changing T or S.
1836		! Local variables
1837	0	real, dimension(SZK_(GV)) :: T0 ! The profile of temperatures [C ~> degC]
1838	0	real, dimension(SZK_(GV)) :: S0 ! The profile of salinities [S ~> ppt]
1839		integer :: i, j, k
1840		character(len=200) :: filename, ts_file, inputdir ! Strings for file/path
1841		character(len=64) :: temp_var, salt_var ! Temperature and salinity names in files
1842		character(len=40) :: mdl = "initialize_temp_salt_from_profile"
1843
1844	0	if (.not.just_read) call callTree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
1845
1846		call get_param(param_file, mdl, "TS_FILE", ts_file, &
1847		"The file with the reference profiles for temperature and salinity.", &
1848	0	fail_if_missing=.not.just_read, do_not_log=just_read)
1849		call get_param(param_file, mdl, "TEMP_IC_VAR", temp_var, &
1850		"The initial condition variable for potential temperature.", &
1851	0	default="PTEMP", do_not_log=just_read)
1852		call get_param(param_file, mdl, "SALT_IC_VAR", salt_var, &
1853		"The initial condition variable for salinity.", &
1854	0	default="SALT", do_not_log=just_read)
1855
1856	0	if (just_read) return ! All run-time parameters have been read, so return.
1857
1858	0	call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
1859	0	inputdir = slasher(inputdir)
1860	0	filename = trim(inputdir)//trim(ts_file)
1861	0	call log_param(param_file, mdl, "INPUTDIR/TS_FILE", filename)
1862	0	if (.not.file_exists(filename)) call MOM_error(FATAL, &
1863	0	" initialize_temp_salt_from_profile: Unable to open "//trim(filename))
1864
1865		! Read the temperatures and salinities from a netcdf file.
1866	0	call MOM_read_data(filename, temp_var, T0(:), scale=US%degC_to_C)
1867	0	call MOM_read_data(filename, salt_var, S0(:), scale=US%ppt_to_S)
1868
1869	0	do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec
1870	0	T(i,j,k) = T0(k) ; S(i,j,k) = S0(k)
1871		enddo ; enddo ; enddo
1872
1873	0	call callTree_leave(trim(mdl)//'()')
1874		end subroutine initialize_temp_salt_from_profile
1875
1876		!> Initializes temperature and salinity by fitting to density
1877	0	subroutine initialize_temp_salt_fit(T, S, G, GV, US, param_file, eqn_of_state, P_Ref, just_read)
1878		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
1879		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
1880		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T !< The potential temperature that is
1881		!! being initialized [C ~> degC].
1882		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S !< The salinity that is being
1883		!! initialized [S ~> ppt].
1884		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1885		type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time
1886		!! parameters.
1887		type(EOS_type), intent(in) :: eqn_of_state !< Equation of state structure
1888		real, intent(in) :: P_Ref !< The coordinate-density reference pressure
1889		!! [R L2 T-2 ~> Pa].
1890		logical, intent(in) :: just_read !< If true, this call will only read
1891		!! parameters without changing T or S.
1892		! Local variables
1893	0	real :: T0(SZK_(GV)) ! Layer potential temperatures [C ~> degC]
1894	0	real :: S0(SZK_(GV)) ! Layer salinities [S ~> ppt]
1895		real :: T_Ref ! Reference Temperature [C ~> degC]
1896		real :: S_Ref ! Reference Salinity [S ~> ppt]
1897	0	real :: pres(SZK_(GV)) ! An array of the reference pressure [R L2 T-2 ~> Pa].
1898	0	real :: drho_dT(SZK_(GV)) ! Derivative of density with temperature [R C-1 ~> kg m-3 degC-1].
1899	0	real :: drho_dS(SZK_(GV)) ! Derivative of density with salinity [R S-1 ~> kg m-3 ppt-1].
1900	0	real :: rho_guess(SZK_(GV)) ! Potential density at T0 & S0 [R ~> kg m-3].
1901		logical :: fit_salin ! If true, accept the prescribed temperature and fit the salinity.
1902		character(len=40) :: mdl = "initialize_temp_salt_fit" ! This subroutine's name.
1903		integer :: i, j, k, itt, nz
1904	0	nz = GV%ke
1905
1906	0	if (.not.just_read) call callTree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
1907
1908		call get_param(param_file, mdl, "T_REF", T_Ref, &
1909		"A reference temperature used in initialization.", &
1910	0	units="degC", scale=US%degC_to_C, fail_if_missing=.not.just_read, do_not_log=just_read)
1911		call get_param(param_file, mdl, "S_REF", S_Ref, &
1912		"A reference salinity used in initialization.", &
1913	0	units="ppt", default=35.0, scale=US%ppt_to_S, do_not_log=just_read)
1914		call get_param(param_file, mdl, "FIT_SALINITY", fit_salin, &
1915		"If true, accept the prescribed temperature and fit the "//&
1916		"salinity; otherwise take salinity and fit temperature.", &
1917	0	default=.false., do_not_log=just_read)
1918
1919	0	if (just_read) return ! All run-time parameters have been read, so return.
1920
1921	0	do k=1,nz
1922	0	pres(k) = P_Ref ; S0(k) = S_Ref
1923	0	T0(k) = T_Ref
1924		enddo
1925
1926	0	call calculate_density(T0(1), S0(1), pres(1), rho_guess(1), eqn_of_state)
1927	0	call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, eqn_of_state, (/1,1/) )
1928
1929	0	if (fit_salin) then
1930		! A first guess of the layers' temperatures.
1931	0	do k=nz,1,-1
1932	0	S0(k) = max(0.0, S0(1) + (GV%Rlay(k) - rho_guess(1)) / drho_dS(1))
1933		enddo
1934		! Refine the guesses for each layer.
1935	0	do itt=1,6
1936	0	call calculate_density(T0, S0, pres, rho_guess, eqn_of_state)
1937	0	call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, eqn_of_state)
1938	0	do k=1,nz
1939	0	S0(k) = max(0.0, S0(k) + (GV%Rlay(k) - rho_guess(k)) / drho_dS(k))
1940		enddo
1941		enddo
1942		else
1943		! A first guess of the layers' temperatures.
1944	0	do k=nz,1,-1
1945	0	T0(k) = T0(1) + (GV%Rlay(k) - rho_guess(1)) / drho_dT(1)
1946		enddo
1947	0	do itt=1,6
1948	0	call calculate_density(T0, S0, pres, rho_guess, eqn_of_state)
1949	0	call calculate_density_derivs(T0, S0, pres, drho_dT, drho_dS, eqn_of_state)
1950	0	do k=1,nz
1951	0	T0(k) = T0(k) + (GV%Rlay(k) - rho_guess(k)) / drho_dT(k)
1952		enddo
1953		enddo
1954		endif
1955
1956	0	do k=1,nz ; do j=G%jsd,G%jed ; do i=G%isd,G%ied
1957	0	T(i,j,k) = T0(k) ; S(i,j,k) = S0(k)
1958		enddo ; enddo ; enddo
1959
1960	0	call callTree_leave(trim(mdl)//'()')
1961		end subroutine initialize_temp_salt_fit
1962
1963		!> Initializes T and S with linear profiles according to reference surface
1964		!! layer salinity and temperature and a specified range.
1965		!!
1966		!! \remark Note that the linear distribution is set up with respect to the layer
1967		!! number, not the physical position).
1968	0	subroutine initialize_temp_salt_linear(T, S, G, GV, US, param_file, just_read)
1969		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
1970		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
1971		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: T !< The potential temperature that is
1972		!! being initialized [C ~> degC]
1973		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: S !< The salinity that is
1974		!! being initialized [S ~> ppt]
1975		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1976		type(param_file_type), intent(in) :: param_file !< A structure to parse for
1977		!! run-time parameters
1978		logical, intent(in) :: just_read !< If present and true,
1979		!! this call will only read parameters
1980		!! without changing T or S.
1981
1982		! Local variables
1983		real :: S_top, S_range ! Reference salinity in the surface layer and its vertical range [S ~> ppt]
1984		real :: T_top, T_range ! Reference temperature in the surface layer and its vertical range [C ~> degC]
1985		character(len=40) :: mdl = "initialize_temp_salt_linear" ! This subroutine's name.
1986		integer :: k
1987
1988	0	if (.not.just_read) call callTree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
1989		call get_param(param_file, mdl, "T_TOP", T_top, &
1990		"Initial temperature of the top surface.", &
1991	0	units="degC", scale=US%degC_to_C, fail_if_missing=.not.just_read, do_not_log=just_read)
1992		call get_param(param_file, mdl, "T_RANGE", T_range, &
1993		"Initial temperature difference (top-bottom).", &
1994	0	units="degC", scale=US%degC_to_C, fail_if_missing=.not.just_read, do_not_log=just_read)
1995		call get_param(param_file, mdl, "S_TOP", S_top, &
1996		"Initial salinity of the top surface.", &
1997	0	units="ppt", scale=US%ppt_to_S, fail_if_missing=.not.just_read, do_not_log=just_read)
1998		call get_param(param_file, mdl, "S_RANGE", S_range, &
1999		"Initial salinity difference (top-bottom).", &
2000	0	units="ppt", scale=US%ppt_to_S, fail_if_missing=.not.just_read, do_not_log=just_read)
2001
2002	0	if (just_read) return ! All run-time parameters have been read, so return.
2003
2004		! Prescribe salinity and temperature, with the extrapolated top interface value prescribed.
2005	0	do k=1,GV%ke
2006	0	S(:,:,k) = S_top - S_range*((real(k)-0.5)/real(GV%ke))
2007	0	T(:,:,k) = T_top - T_range*((real(k)-0.5)/real(GV%ke))
2008		enddo
2009
2010		! Prescribe salinity and temperature, but with the top layer value matching the surface value.
2011		! S(:,:,1) = S_top ; T(:,:,1) = T_top
2012		! do k=2,GV%ke
2013		! S(:,:,k) = S_top - S_range * (real(k-1) / real(GV%ke-1))
2014		! T(:,:,k) = T_top - T_range * (real(k-1) / real(GV%ke-1))
2015		! enddo
2016
2017	0	call callTree_leave(trim(mdl)//'()')
2018		end subroutine initialize_temp_salt_linear
2019
2020		!> This subroutine sets the inverse restoration time (Idamp), and
2021		!! the values towards which the interface heights and an arbitrary
2022		!! number of tracers should be restored within each sponge. The
2023		!! interface height is always subject to damping, and must always be
2024		!! the first registered field.
2025	0	subroutine initialize_sponges_file(G, GV, US, use_temperature, tv, u, v, depth_tot, param_file, &
2026		Layer_CSp, ALE_CSp, Time)
2027		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure.
2028		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
2029		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2030		logical, intent(in) :: use_temperature !< If true, T & S are state variables.
2031		type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic
2032		!! variables.
2033		real, target, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
2034		intent(in) :: u !< The zonal velocity that is being
2035		!! initialized [L T-1 ~> m s-1]
2036		real, target, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
2037		intent(in) :: v !< The meridional velocity that is being
2038		!! initialized [L T-1 ~> m s-1]
2039		real, dimension(SZI_(G),SZJ_(G)), &
2040		intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
2041		type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters.
2042		type(sponge_CS), pointer :: Layer_CSp !< A pointer that is set to point to the control
2043		!! structure for this module (in layered mode).
2044		type(ALE_sponge_CS), pointer :: ALE_CSp !< A pointer that is set to point to the control
2045		!! structure for this module (in ALE mode).
2046		type(time_type), intent(in) :: Time !< Time at the start of the run segment. Time_in
2047		!! overrides any value set for Time.
2048		! Local variables
2049	0	real, allocatable, dimension(:,:,:) :: eta ! The target interface heights [Z ~> m].
2050	0	real, allocatable, dimension(:,:,:) :: dz ! The target interface thicknesses in height units [Z ~> m]
2051
2052		real, dimension (SZI_(G),SZJ_(G),SZK_(GV)) :: &
2053	0	tmp, & ! A temporary array for temperatures [C ~> degC] or other tracers.
2054	0	tmp2 ! A temporary array for salinities [S ~> ppt]
2055		real, dimension (SZI_(G),SZJ_(G)) :: &
2056	0	tmp_2d ! A temporary array for mixed layer densities [R ~> kg m-3]
2057	0	real, allocatable, dimension(:,:,:) :: tmp_T ! A temporary array for reading sponge target temperatures
2058		! on the vertical grid of the input file [C ~> degC]
2059	0	real, allocatable, dimension(:,:,:) :: tmp_S ! A temporary array for reading sponge target salinities
2060		! on the vertical grid of the input file [S ~> ppt]
2061	0	real, allocatable, dimension(:,:,:) :: tmp_u ! Temporary array for reading sponge target zonal
2062		! velocities on the vertical grid of the input file [L T-1 ~> m s-1]
2063	0	real, allocatable, dimension(:,:,:) :: tmp_v ! Temporary array for reading sponge target meridional
2064		! velocities on the vertical grid of the input file [L T-1 ~> m s-1]
2065
2066	0	real :: Idamp(SZI_(G),SZJ_(G)) ! The sponge damping rate [T-1 ~> s-1]
2067	0	real :: Idamp_u(SZIB_(G),SZJ_(G)) ! The sponge damping rate for velocity fields [T-1 ~> s-1]
2068	0	real :: Idamp_v(SZI_(G),SZJB_(G)) ! The sponge damping rate for velocity fields [T-1 ~> s-1]
2069	0	real :: pres(SZI_(G)) ! An array of the reference pressure [R L2 T-2 ~> Pa]
2070
2071		integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
2072		integer :: i, j, k, is, ie, js, je, nz
2073		integer :: isd, ied, jsd, jed
2074		integer, dimension(4) :: siz
2075		integer :: nz_data ! The size of the sponge source grid
2076		logical :: sponge_uv ! Apply sponges in u and v, in addition to tracers.
2077		character(len=40) :: potemp_var, salin_var, u_var, v_var, Idamp_var, Idamp_u_var, Idamp_v_var, eta_var
2078		character(len=40) :: mdl = "initialize_sponges_file"
2079		character(len=200) :: damping_file, uv_damping_file, state_file, state_uv_file ! Strings for filenames
2080		character(len=200) :: filename, inputdir ! Strings for file/path and path.
2081
2082		logical :: use_ALE ! True if ALE is being used, False if in layered mode
2083		logical :: time_space_interp_sponge ! If true use sponge data that need to be interpolated in both
2084		! the horizontal dimension and in time prior to vertical remapping.
2085
2086	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
2087	0	isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
2088
2089	0	pres(:) = 0.0 ; tmp(:,:,:) = 0.0 ; Idamp(:,:) = 0.0 ; Idamp_u(:,:) = 0.0 ; Idamp_v(:,:) = 0.0
2090
2091	0	call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
2092	0	inputdir = slasher(inputdir)
2093		call get_param(param_file, mdl, "SPONGE_DAMPING_FILE", damping_file, &
2094		"The name of the file with the sponge damping rates.", &
2095	0	fail_if_missing=.true.)
2096		call get_param(param_file, mdl, "SPONGE_STATE_FILE", state_file, &
2097		"The name of the file with the state to damp toward.", &
2098	0	default=damping_file)
2099		call get_param(param_file, mdl, "SPONGE_PTEMP_VAR", potemp_var, &
2100		"The name of the potential temperature variable in "//&
2101	0	"SPONGE_STATE_FILE.", default="PTEMP")
2102		call get_param(param_file, mdl, "SPONGE_SALT_VAR", salin_var, &
2103		"The name of the salinity variable in "//&
2104	0	"SPONGE_STATE_FILE.", default="SALT")
2105		call get_param(param_file, mdl, "SPONGE_UV", sponge_uv, &
2106		"Apply sponges in u and v, in addition to tracers.", &
2107	0	default=.false.)
2108	0	if (sponge_uv) then
2109		call get_param(param_file, mdl, "SPONGE_UV_STATE_FILE", state_uv_file, &
2110		"The name of the file with the state to damp UV toward.", &
2111	0	default=damping_file)
2112		call get_param(param_file, mdl, "SPONGE_U_VAR", u_var, &
2113		"The name of the zonal velocity variable in "//&
2114	0	"SPONGE_UV_STATE_FILE.", default="UVEL")
2115		call get_param(param_file, mdl, "SPONGE_V_VAR", v_var, &
2116		"The name of the vertical velocity variable in "//&
2117	0	"SPONGE_UV_STATE_FILE.", default="VVEL")
2118		endif
2119		call get_param(param_file, mdl, "SPONGE_ETA_VAR", eta_var, &
2120		"The name of the interface height variable in "//&
2121	0	"SPONGE_STATE_FILE.", default="ETA")
2122		call get_param(param_file, mdl, "SPONGE_IDAMP_VAR", Idamp_var, &
2123		"The name of the inverse damping rate variable in "//&
2124	0	"SPONGE_DAMPING_FILE.", default="Idamp")
2125	0	if (sponge_uv) then
2126		call get_param(param_file, mdl, "SPONGE_UV_DAMPING_FILE", uv_damping_file, &
2127		"The name of the file with sponge damping rates for the velocity variables.", &
2128	0	default=damping_file)
2129		call get_param(param_file, mdl, "SPONGE_IDAMP_U_var", Idamp_u_var, &
2130		"The name of the inverse damping rate variable in "//&
2131	0	"SPONGE_UV_DAMPING_FILE for the velocities.", default=Idamp_var)
2132		call get_param(param_file, mdl, "SPONGE_IDAMP_V_var", Idamp_v_var, &
2133		"The name of the inverse damping rate variable in "//&
2134	0	"SPONGE_UV_DAMPING_FILE for the velocities.", default=Idamp_var)
2135		endif
2136	0	call get_param(param_file, mdl, "USE_REGRIDDING", use_ALE, default=.false., do_not_log=.true.)
2137
2138		call get_param(param_file, mdl, "INTERPOLATE_SPONGE_TIME_SPACE", time_space_interp_sponge, &
2139		"If True, perform on-the-fly regridding in lat-lon-time of sponge restoring data.", &
2140	0	default=.false.)
2141
2142		! Read in sponge damping rate for tracers
2143	0	filename = trim(damping_file)
2144	0	if (scan(damping_file, '/')== 0) then ! prepend inputdir if only a filename is given
2145	0	filename = trim(inputdir)//trim(damping_file)
2146		endif
2147	0	call log_param(param_file, mdl, "INPUTDIR/SPONGE_DAMPING_FILE", filename)
2148	0	if (.not.file_exists(filename, G%Domain)) &
2149	0	call MOM_error(FATAL, " initialize_sponges: Unable to open "//trim(filename))
2150
2151	0	if (time_space_interp_sponge .and. .not. use_ALE) &
2152	0	call MOM_error(FATAL, " initialize_sponges: Time-varying sponges are currently unavailable in layered mode ")
2153
2154	0	call MOM_read_data(filename, Idamp_var, Idamp(:,:), G%Domain, scale=US%T_to_s)
2155
2156		! Read in sponge damping rate for velocities
2157	0	if (sponge_uv) then
2158	0	if (separate_idamp_for_uv()) then
2159	0	filename = trim(inputdir)//trim(uv_damping_file)
2160	0	call log_param(param_file, mdl, "INPUTDIR/SPONGE_UV_DAMPING_FILE", filename)
2161
2162	0	if (.not.file_exists(filename, G%Domain)) &
2163	0	call MOM_error(FATAL, " initialize_sponges: Unable to open "//trim(filename))
2164
2165	0	call MOM_read_vector(filename, Idamp_u_var,Idamp_v_var,Idamp_u(:,:),Idamp_v(:,:), G%Domain, scale=US%T_to_s)
2166		else
2167		! call MOM_error(FATAL, "Must provide SPONGE_IDAMP_U_var and SPONGE_IDAMP_V_var")
2168	0	call pass_var(Idamp,G%Domain)
2169	0	do j=G%jsc,G%jec
2170	0	do i=G%iscB,G%iecB
2171	0	Idamp_u(I,j) = 0.5*(Idamp(i,j)+Idamp(i+1,j))
2172		enddo
2173		enddo
2174	0	do j=G%jscB,G%jecB
2175	0	do i=G%isc,G%iec
2176	0	Idamp_v(i,J) = 0.5*(Idamp(i,j)+Idamp(i,j+1))
2177		enddo
2178		enddo
2179		endif
2180		endif
2181
2182		! Now register all of the fields which are damped in the sponge.
2183		! By default, momentum is advected vertically within the sponge, but
2184		! momentum is typically not damped within the sponge.
2185
2186	0	filename = trim(inputdir)//trim(state_file)
2187	0	call log_param(param_file, mdl, "INPUTDIR/SPONGE_STATE_FILE", filename)
2188	0	if (.not.file_exists(filename, G%Domain)) &
2189	0	call MOM_error(FATAL, " initialize_sponges: Unable to open "//trim(filename))
2190
2191
2192	0	if (.not. use_ALE) then
2193		! The first call to set_up_sponge_field is for the interface heights if in layered mode.
2194	0	allocate(eta(isd:ied,jsd:jed,nz+1), source=0.0)
2195	0	call MOM_read_data(filename, eta_var, eta(:,:,:), G%Domain, scale=US%m_to_Z)
2196
2197	0	do j=js,je ; do i=is,ie
2198	0	eta(i,j,nz+1) = -depth_tot(i,j)
2199		enddo ; enddo
2200	0	do k=nz,1,-1 ; do j=js,je ; do i=is,ie
2201	0	if (eta(i,j,K) < (eta(i,j,K+1) + GV%Angstrom_Z)) &
2202	0	eta(i,j,K) = eta(i,j,K+1) + GV%Angstrom_Z
2203		enddo ; enddo ; enddo
2204		! Set the sponge damping rates so that the model will know where to
2205		! apply the sponges, along with the interface heights.
2206	0	call initialize_sponge(Idamp, eta, G, param_file, Layer_CSp, GV)
2207	0	deallocate(eta)
2208
2209	0	if ( GV%nkml>0) then
2210		! This call to set_up_sponge_ML_density registers the target values of the
2211		! mixed layer density, which is used in determining which layers can be
2212		! inflated without causing static instabilities.
2213	0	do i=is,ie ; pres(i) = tv%P_Ref ; enddo
2214	0	EOSdom(:) = EOS_domain(G%HI)
2215
2216	0	call MOM_read_data(filename, potemp_var, tmp(:,:,:), G%Domain, scale=US%degC_to_C)
2217	0	call MOM_read_data(filename, salin_var, tmp2(:,:,:), G%Domain, scale=US%ppt_to_S)
2218
2219	0	do j=js,je
2220	0	call calculate_density(tmp(:,j,1), tmp2(:,j,1), pres, tmp_2d(:,j), tv%eqn_of_state, EOSdom)
2221		enddo
2222
2223	0	call set_up_sponge_ML_density(tmp_2d, G, Layer_CSp)
2224		endif
2225
2226		! Now register all of the tracer fields which are damped in the
2227		! sponge. By default, momentum is advected vertically within the
2228		! sponge, but momentum is typically not damped within the sponge.
2229
2230
2231		! The remaining calls to set_up_sponge_field can be in any order.
2232	0	if ( use_temperature) then
2233	0	call MOM_read_data(filename, potemp_var, tmp(:,:,:), G%Domain, scale=US%degC_to_C)
2234	0	call set_up_sponge_field(tmp, tv%T, G, GV, nz, Layer_CSp)
2235	0	call MOM_read_data(filename, salin_var, tmp2(:,:,:), G%Domain, scale=US%ppt_to_S)
2236	0	call set_up_sponge_field(tmp2, tv%S, G, GV, nz, Layer_CSp)
2237		endif
2238
2239		! else
2240		! Initialize sponges without supplying sponge grid
2241		! if (sponge_uv) then
2242		! call initialize_ALE_sponge(Idamp, G, GV, US, param_file, ALE_CSp, Idamp_u, Idamp_v)
2243		! else
2244		! call initialize_ALE_sponge(Idamp, G, GV, US, param_file, ALE_CSp)
2245		! endif
2246		endif
2247
2248
2249	0	if (use_ALE) then ! ALE mode
2250	0	if (.not. time_space_interp_sponge) then
2251	0	call field_size(filename,eta_var,siz,no_domain=.true.)
2252	0	if (siz(1) /= G%ieg-G%isg+1 .or. siz(2) /= G%jeg-G%jsg+1) &
2253	0	call MOM_error(FATAL,"initialize_sponge_file: Array size mismatch for sponge data.")
2254	0	nz_data = siz(3)-1
2255	0	allocate(eta(isd:ied,jsd:jed,nz_data+1))
2256	0	allocate(dz(isd:ied,jsd:jed,nz_data))
2257	0	call MOM_read_data(filename, eta_var, eta(:,:,:), G%Domain, scale=US%m_to_Z)
2258	0	do j=js,je ; do i=is,ie
2259	0	eta(i,j,nz_data+1) = -depth_tot(i,j)
2260		enddo ; enddo
2261	0	do k=nz_data,1,-1 ; do j=js,je ; do i=is,ie
2262	0	if (eta(i,j,K) < (eta(i,j,K+1) + GV%Angstrom_Z)) &
2263	0	eta(i,j,K) = eta(i,j,K+1) + GV%Angstrom_Z
2264		enddo ; enddo ; enddo
2265	0	do k=1,nz_data ; do j=js,je ; do i=is,ie
2266	0	dz(i,j,k) = eta(i,j,k)-eta(i,j,k+1)
2267		enddo ; enddo ; enddo
2268	0	deallocate(eta)
2269
2270	0	if (use_temperature) then
2271	0	allocate(tmp_T(isd:ied,jsd:jed,nz_data))
2272	0	allocate(tmp_S(isd:ied,jsd:jed,nz_data))
2273	0	call MOM_read_data(filename, potemp_var, tmp_T(:,:,:), G%Domain, scale=US%degC_to_C)
2274	0	call MOM_read_data(filename, salin_var, tmp_S(:,:,:), G%Domain, scale=US%ppt_to_S)
2275		endif
2276
2277	0	if (sponge_uv) then
2278		call initialize_ALE_sponge(Idamp, G, GV, param_file, ALE_CSp, dz, nz_data, Idamp_u, Idamp_v, &
2279	0	data_h_is_Z=.true.)
2280		else
2281		call initialize_ALE_sponge(Idamp, G, GV, param_file, ALE_CSp, dz, nz_data, &
2282	0	data_h_is_Z=.true.)
2283		endif
2284	0	if (use_temperature) then
2285		call set_up_ALE_sponge_field(tmp_T, G, GV, tv%T, ALE_CSp, 'temp', &
2286	0	sp_long_name='temperature', sp_unit='degC s-1')
2287		call set_up_ALE_sponge_field(tmp_S, G, GV, tv%S, ALE_CSp, 'salt', &
2288	0	sp_long_name='salinity', sp_unit='g kg-1 s-1')
2289	0	deallocate(tmp_S)
2290	0	deallocate(tmp_T)
2291		endif
2292	0	deallocate(dz)
2293
2294	0	if (sponge_uv) then
2295	0	filename = trim(inputdir)//trim(state_uv_file)
2296	0	call log_param(param_file, mdl, "INPUTDIR/SPONGE_STATE_UV_FILE", filename)
2297	0	if (.not.file_exists(filename, G%Domain)) &
2298	0	call MOM_error(FATAL, " initialize_sponges: Unable to open "//trim(filename))
2299	0	allocate(tmp_u(G%IsdB:G%IedB,jsd:jed,nz_data))
2300	0	allocate(tmp_v(isd:ied,G%JsdB:G%JedB,nz_data))
2301	0	call MOM_read_vector(filename, u_var, v_var, tmp_u(:,:,:), tmp_v(:,:,:), G%Domain, scale=US%m_s_to_L_T)
2302	0	call set_up_ALE_sponge_vel_field(tmp_u, tmp_v, G, GV, u, v, ALE_CSp)
2303	0	deallocate(tmp_u,tmp_v)
2304		endif
2305		else
2306		! Initialize sponges without supplying sponge grid
2307	0	if (sponge_uv) then
2308	0	call initialize_ALE_sponge(Idamp, G, GV, US, param_file, ALE_CSp, Idamp_u, Idamp_v)
2309		else
2310	0	call initialize_ALE_sponge(Idamp, G, GV, US, param_file, ALE_CSp)
2311		endif
2312		! The remaining calls to set_up_sponge_field can be in any order.
2313	0	if ( use_temperature) then
2314		call set_up_ALE_sponge_field(filename, potemp_var, Time, G, GV, US, tv%T, ALE_CSp, &
2315	0	'temp', sp_long_name='temperature', sp_unit='degC s-1', scale=US%degC_to_C)
2316		call set_up_ALE_sponge_field(filename, salin_var, Time, G, GV, US, tv%S, ALE_CSp, &
2317	0	'salt', sp_long_name='salinity', sp_unit='g kg-1 s-1', scale=US%ppt_to_S)
2318		endif
2319	0	if (sponge_uv) then
2320	0	filename = trim(inputdir)//trim(state_uv_file)
2321	0	call log_param(param_file, mdl, "INPUTDIR/SPONGE_STATE_UV_FILE", filename)
2322	0	if (.not.file_exists(filename, G%Domain)) &
2323	0	call MOM_error(FATAL, " initialize_sponges: Unable to open "//trim(filename))
2324		call set_up_ALE_sponge_vel_field(filename, u_var, filename, v_var, Time, G, GV, US, &
2325	0	ALE_CSp, u, v, scale=US%m_s_to_L_T)
2326		endif
2327		endif
2328		endif
2329
2330	0	if (sponge_uv .and. .not. use_ALE) call MOM_error(FATAL,'initialize_sponges_file: '// &
2331	0	'UV damping to target values only available in ALE mode')
2332
2333
2334		contains
2335
2336		! returns true if a separate idamp is provided for u and/or v
2337	0	logical function separate_idamp_for_uv()
2338		separate_idamp_for_uv = (lowercase(damping_file)/=lowercase(uv_damping_file) .or. &
2339	0	lowercase(Idamp_var)/=lowercase(Idamp_u_var) .or. lowercase(Idamp_var)/=lowercase(Idamp_v_var))
2340	0	end function separate_idamp_for_uv
2341
2342		end subroutine initialize_sponges_file
2343
2344	0	subroutine initialize_oda_incupd_file(G, GV, US, use_temperature, tv, h, u, v, param_file, &
2345		oda_incupd_CSp, restart_CS, Time)
2346		type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
2347		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
2348		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2349		logical, intent(in) :: use_temperature !< If true, T & S are state variables.
2350		type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various thermodynamic
2351		!! variables.
2352		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
2353		intent(inout) :: h !< Layer thickness [H ~> m or kg m-2] (in)
2354
2355		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
2356		intent(in) :: u !< The zonal velocity that is being
2357		!! initialized [L T-1 ~> m s-1]
2358		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
2359		intent(in) :: v !< The meridional velocity that is being
2360		!! initialized [L T-1 ~> m s-1]
2361		type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters.
2362		type(oda_incupd_CS), pointer :: oda_incupd_CSp !< A pointer that is set to point to the control
2363		!! structure for this module.
2364		type(MOM_restart_CS), intent(in) :: restart_CS !< MOM restart control structure
2365		type(time_type), intent(in) :: Time !< Time at the start of the run segment. Time_in
2366		!! overrides any value set for Time.
2367		! Local variables
2368	0	real, allocatable, dimension(:,:,:) :: hoda ! The layer thickness increment and oda layer thickness [H ~> m or kg m-2]
2369	0	real, allocatable, dimension(:,:,:) :: tmp_tr ! A temporary array for reading oda tracer increments
2370		! on the vertical grid of the input file, used for both
2371		! temperatures [C ~> degC] and salinities [S ~> ppt]
2372	0	real, allocatable, dimension(:,:,:) :: tmp_u ! Temporary array for reading oda zonal velocity
2373		! increments on the vertical grid of the input file [L T-1 ~> m s-1]
2374	0	real, allocatable, dimension(:,:,:) :: tmp_v ! Temporary array for reading oda meridional velocity
2375		! increments on the vertical grid of the input file [L T-1 ~> m s-1]
2376
2377		integer :: is, ie, js, je, nz
2378		integer :: isd, ied, jsd, jed
2379
2380		integer, dimension(4) :: siz
2381		integer :: nz_data ! The size of the sponge source grid
2382		logical :: oda_inc ! input files are increments (true) or full fields (false)
2383		logical :: save_inc ! save increments if using full fields
2384		logical :: uv_inc ! use u and v increments
2385		logical :: reset_ncount ! reset ncount to zero if true
2386
2387		character(len=40) :: tempinc_var, salinc_var, uinc_var, vinc_var, h_var
2388		character(len=40) :: mdl = "initialize_oda_incupd_file"
2389		character(len=200) :: inc_file, uv_inc_file ! Strings for filenames
2390		character(len=200) :: filename, inputdir ! Strings for file/path and path.
2391
2392		! logical :: use_ALE ! True if ALE is being used, False if in layered mode
2393
2394	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
2395	0	isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
2396
2397	0	call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
2398	0	inputdir = slasher(inputdir)
2399
2400		call get_param(param_file, mdl, "ODA_INCUPD_FILE", inc_file, &
2401		"The name of the file with the T,S,h increments.", &
2402	0	fail_if_missing=.true.)
2403		call get_param(param_file, mdl, "ODA_INCUPD_INC", oda_inc, &
2404		"INCUPD files are increments and not full fields.", &
2405	0	default=.true.)
2406	0	if (.not.oda_inc) then
2407		call get_param(param_file, mdl, "ODA_INCUPD_SAVE", save_inc, &
2408		"If true, save the increments when using full fields.", &
2409	0	default=.false.)
2410		endif
2411		call get_param(param_file, mdl, "ODA_INCUPD_RESET_NCOUNT", reset_ncount, &
2412		"If True, reinitialize number of updates already done, ncount.", &
2413	0	default=.true.)
2414	0	if (.not.oda_inc .and. .not.reset_ncount) &
2415		call MOM_error(FATAL, " initialize_oda_incupd: restarting during update "// &
2416	0	"necessitates increments input file")
2417
2418		call get_param(param_file, mdl, "ODA_TEMPINC_VAR", tempinc_var, &
2419		"The name of the potential temperature inc. variable in "//&
2420	0	"ODA_INCUPD_FILE.", default="ptemp_inc")
2421		call get_param(param_file, mdl, "ODA_SALTINC_VAR", salinc_var, &
2422		"The name of the salinity inc. variable in "//&
2423	0	"ODA_INCUPD_FILE.", default="sal_inc")
2424		call get_param(param_file, mdl, "ODA_THK_VAR", h_var, &
2425		"The name of the layer thickness variable in "//&
2426	0	"ODA_INCUPD_FILE.", default="h")
2427		call get_param(param_file, mdl, "ODA_INCUPD_UV", uv_inc, &
2428		"use U,V increments.", &
2429	0	default=.true.)
2430		call get_param(param_file, mdl, "ODA_INCUPD_UV_FILE", uv_inc_file, &
2431		"The name of the file with the U,V increments.", &
2432	0	default=inc_file)
2433		call get_param(param_file, mdl, "ODA_UINC_VAR", uinc_var, &
2434		"The name of the zonal vel. inc. variable in "//&
2435	0	"ODA_INCUPD_FILE.", default="u_inc")
2436		call get_param(param_file, mdl, "ODA_VINC_VAR", vinc_var, &
2437		"The name of the meridional vel. inc. variable in "//&
2438	0	"ODA_INCUPD_FILE.", default="v_inc")
2439
2440		! call get_param(param_file, mdl, "USE_REGRIDDING", use_ALE, default=.false., do_not_log=.true.)
2441
2442		! Read in incremental update for tracers
2443	0	filename = trim(inc_file)
2444	0	if (scan(inc_file, '/')== 0) then ! prepend inputdir if only a filename is given
2445	0	filename = trim(inputdir)//trim(inc_file)
2446		endif
2447	0	call log_param(param_file, mdl, "INPUTDIR/ODA_INCUPD_FILE", filename)
2448	0	if (.not.file_exists(filename, G%Domain)) &
2449	0	call MOM_error(FATAL, " initialize_oda_incupd: Unable to open "//trim(filename))
2450
2451	0	call field_size(filename,h_var,siz,no_domain=.true.)
2452	0	if (siz(1) /= G%ieg-G%isg+1 .or. siz(2) /= G%jeg-G%jsg+1) &
2453	0	call MOM_error(FATAL,"initialize_oda_incupd_file: Array size mismatch for oda data.")
2454	0	nz_data = siz(3)
2455		! get h increments
2456	0	allocate(hoda(isd:ied,jsd:jed,nz_data))
2457	0	call MOM_read_data(filename, h_var , hoda(:,:,:), G%Domain, scale=US%m_to_Z)
2458	0	call initialize_oda_incupd( G, GV, US, param_file, oda_incupd_CSp, hoda, nz_data, restart_CS)
2459	0	deallocate(hoda)
2460
2461		! set-up T and S increments arrays
2462	0	if (use_temperature) then
2463	0	allocate(tmp_tr(isd:ied,jsd:jed,nz_data))
2464		! temperature inc. in array Inc(1)
2465	0	call MOM_read_data(filename, tempinc_var, tmp_tr(:,:,:), G%Domain, scale=US%degC_to_C)
2466	0	call set_up_oda_incupd_field(tmp_tr, G, GV, oda_incupd_CSp)
2467		! salinity inc. in array Inc(2)
2468	0	call MOM_read_data(filename, salinc_var, tmp_tr(:,:,:), G%Domain, scale=US%ppt_to_S)
2469	0	call set_up_oda_incupd_field(tmp_tr, G, GV, oda_incupd_CSp)
2470	0	deallocate(tmp_tr)
2471		endif
2472
2473		! set-up U and V increments arrays
2474	0	if (uv_inc) then
2475	0	filename = trim(inputdir)//trim(uv_inc_file)
2476	0	call log_param(param_file, mdl, "INPUTDIR/ODA_INCUPD_UV_FILE", filename)
2477	0	if (.not.file_exists(filename, G%Domain)) &
2478	0	call MOM_error(FATAL, " initialize_oda_incupd_uv: Unable to open "//trim(filename))
2479	0	allocate(tmp_u(G%IsdB:G%IedB,jsd:jed,nz_data), source=0.0)
2480	0	allocate(tmp_v(isd:ied,G%JsdB:G%JedB,nz_data), source=0.0)
2481	0	call MOM_read_vector(filename, uinc_var, vinc_var, tmp_u, tmp_v, G%Domain, scale=US%m_s_to_L_T)
2482	0	call set_up_oda_incupd_vel_field(tmp_u, tmp_v, G, GV, oda_incupd_CSp)
2483	0	deallocate(tmp_u, tmp_v)
2484		endif
2485
2486		! calculate increments if input are full fields
2487	0	if (oda_inc) then ! input are increments
2488	0	if (is_root_pe()) call MOM_mesg("incupd using increments fields ")
2489		else ! inputs are full fields
2490	0	if (is_root_pe()) call MOM_mesg("incupd using full fields ")
2491	0	call calc_oda_increments(h, tv, u, v, G, GV, US, oda_incupd_CSp)
2492	0	if (save_inc) then
2493	0	call output_oda_incupd_inc(Time, G, GV, param_file, oda_incupd_CSp, US)
2494		endif
2495		endif ! not oda_inc
2496
2497	0	end subroutine initialize_oda_incupd_file
2498
2499
2500		!> This subroutine sets the 4 bottom depths at velocity points to be the
2501		!! maximum of the adjacent depths.
2502	0	subroutine set_velocity_depth_max(G)
2503		type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
2504		! Local variables
2505		integer :: i, j
2506
2507	0	do I=G%isd,G%ied-1 ; do j=G%jsd,G%jed
2508	0	G%Dblock_u(I,j) = G%mask2dCu(I,j) * max(G%bathyT(i,j), G%bathyT(i+1,j))
2509	0	G%Dopen_u(I,j) = G%Dblock_u(I,j)
2510		enddo ; enddo
2511	0	do i=G%isd,G%ied ; do J=G%jsd,G%jed-1
2512	0	G%Dblock_v(I,J) = G%mask2dCv(i,J) * max(G%bathyT(i,j), G%bathyT(i,j+1))
2513	0	G%Dopen_v(I,J) = G%Dblock_v(I,J)
2514		enddo ; enddo
2515	0	end subroutine set_velocity_depth_max
2516
2517		!> This subroutine sets the 4 bottom depths at velocity points to be the
2518		!! minimum of the adjacent depths.
2519	0	subroutine set_velocity_depth_min(G)
2520		type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
2521		! Local variables
2522		integer :: i, j
2523
2524	0	do I=G%isd,G%ied-1 ; do j=G%jsd,G%jed
2525	0	G%Dblock_u(I,j) = G%mask2dCu(I,j) * min(G%bathyT(i,j), G%bathyT(i+1,j))
2526	0	G%Dopen_u(I,j) = G%Dblock_u(I,j)
2527		enddo ; enddo
2528	0	do i=G%isd,G%ied ; do J=G%jsd,G%jed-1
2529	0	G%Dblock_v(I,J) = G%mask2dCv(i,J) * min(G%bathyT(i,j), G%bathyT(i,j+1))
2530	0	G%Dopen_v(I,J) = G%Dblock_v(I,J)
2531		enddo ; enddo
2532	0	end subroutine set_velocity_depth_min
2533
2534		!> This subroutine determines the isopycnal or other coordinate interfaces and
2535		!! layer potential temperatures and salinities directly from a z-space file on
2536		!! a latitude-longitude grid.
2537	0	subroutine MOM_temp_salt_initialize_from_Z(h, tv, depth_tot, G, GV, US, PF, just_read, frac_shelf_h)
2538		type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure
2539		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
2540		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
2541		intent(out) :: h !< Layer thicknesses being initialized [H ~> m or kg m-2]
2542		type(thermo_var_ptrs), intent(inout) :: tv !< A structure pointing to various thermodynamic
2543		!! variables including temperature and salinity
2544		real, dimension(SZI_(G),SZJ_(G)), &
2545		intent(in) :: depth_tot !< The nominal total depth of the ocean [Z ~> m]
2546		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2547		type(param_file_type), intent(in) :: PF !< A structure indicating the open file
2548		!! to parse for model parameter values.
2549		logical, intent(in) :: just_read !< If true, this call will only read
2550		!! parameters without changing T or S.
2551		real, dimension(SZI_(G),SZJ_(G)), &
2552		optional, intent(in) :: frac_shelf_h !< The fraction of the grid cell covered
2553		!! by a floating ice shelf [nondim].
2554
2555		! Local variables
2556		character(len=200) :: filename !< The name of an input file containing temperature
2557		!! and salinity in z-space; by default it is also used for ice shelf area.
2558		character(len=200) :: tfilename !< The name of an input file containing temperature in z-space.
2559		character(len=200) :: sfilename !< The name of an input file containing salinity in z-space.
2560		character(len=200) :: inputdir !! The directory where NetCDF input files are.
2561		character(len=200) :: mesg
2562
2563		type(EOS_type), pointer :: eos => NULL()
2564	0	type(thermo_var_ptrs) :: tv_loc ! A temporary thermo_var container
2565	0	type(verticalGrid_type) :: GV_loc ! A temporary vertical grid structure
2566		! This include declares and sets the variable "version".
2567		# include "version_variable.h"
2568		character(len=40) :: mdl = "MOM_initialize_layers_from_Z" ! This module's name.
2569
2570		integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
2571		integer :: is, ie, js, je, nz ! compute domain indices
2572		integer :: isg, ieg, jsg, jeg ! global extent
2573		integer :: isd, ied, jsd, jed ! data domain indices
2574
2575		integer :: i, j, k, ks
2576		integer :: nkml ! The number of layers in the mixed layer.
2577
2578		integer :: inconsistent ! The total number of cells with in consistent topography and layer thicknesses.
2579		integer :: kd ! The number of levels in the input data
2580		integer :: nkd ! number of levels to use for regridding input arrays
2581		real :: eps_Z ! A negligibly thin layer thickness [Z ~> m].
2582		real :: eps_rho ! A negligibly small density difference [R ~> kg m-3].
2583		real :: PI_180 ! for conversion from degrees to radians [radian degree-1]
2584		real :: Hmix_default ! The default initial mixed layer depth [Z ~> m].
2585		real :: Hmix_depth ! The mixed layer depth in the initial condition [Z ~> m].
2586		real :: missing_value_temp ! The missing value in the input temperature field [C ~> degC]
2587		real :: missing_value_salt ! The missing value in the input salinity field [S ~> ppt]
2588		real :: tol_temp ! The tolerance for changes in temperature during the horizontal
2589		! interpolation from an input dataset [C ~> degC]
2590		real :: tol_sal ! The tolerance for changes in salinity during the horizontal
2591		! interpolation from an input dataset [S ~> ppt]
2592		logical :: correct_thickness ! If true, correct the column thicknesses to match the topography
2593		real :: h_tolerance ! A parameter that controls the tolerance when adjusting the
2594		! thickness to fit the bathymetry [Z ~> m].
2595		real :: tol_dz_bot ! A tolerance for detecting inconsistent bottom depths when
2596		! correct_thickness is false [Z ~> m]
2597		character(len=40) :: potemp_var, salin_var
2598
2599		integer, parameter :: niter=10 ! number of iterations for t/s adjustment to layer density
2600		logical :: adjust_temperature = .true. ! fit t/s to target densities
2601		real :: temp_land_fill ! A temperature value to use for land points [C ~> degC]
2602		real :: salt_land_fill ! A salinity value to use for land points [C ~> degC]
2603
2604		! data arrays
2605	0	real, dimension(:), allocatable :: z_edges_in ! Input data interface heights or depths [Z ~> m]
2606	0	real, dimension(:), allocatable :: z_in ! Input data cell heights or depths [Z ~> m]
2607	0	real, dimension(:), allocatable :: Rb ! Interface densities [R ~> kg m-3]
2608	0	real, dimension(:,:,:), allocatable, target :: temp_z ! Input temperatures [C ~> degC]
2609	0	real, dimension(:,:,:), allocatable, target :: salt_z ! Input salinities [S ~> ppt]
2610	0	real, dimension(:,:,:), allocatable, target :: mask_z ! 1 for valid data points [nondim]
2611	0	real, dimension(:,:,:), allocatable :: rho_z ! Densities in Z-space [R ~> kg m-3]
2612	0	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: zi ! Interface heights [Z ~> m]
2613	0	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dz ! Layer thicknesses in height units [Z ~> m]
2614	0	real, dimension(SZI_(G),SZJ_(G)) :: Z_bottom ! The (usually negative) height of the seafloor
2615		! relative to the surface [Z ~> m].
2616	0	integer, dimension(SZI_(G),SZJ_(G)) :: nlevs ! The number of levels in each column with valid data
2617	0	real, dimension(SZI_(G)) :: press ! Pressures [R L2 T-2 ~> Pa].
2618
2619		! Local variables for ALE remapping
2620	0	real, dimension(:), allocatable :: hTarget ! Target thicknesses [Z ~> m].
2621	0	real, dimension(:,:,:), allocatable, target :: tmpT1dIn ! Input temperatures on a model-sized grid [C ~> degC]
2622	0	real, dimension(:,:,:), allocatable, target :: tmpS1dIn ! Input salinities on a model-sized grid [S ~> ppt]
2623	0	real, dimension(:,:,:), allocatable :: tmp_mask_in ! The valid data mask on a model-sized grid [nondim]
2624	0	real, dimension(:,:,:), allocatable :: dz1 ! Input grid thicknesses in depth units [Z ~> m]
2625	0	real, dimension(:,:,:), allocatable :: h1 ! Thicknesses on the input grid [H ~> m or kg m-2].
2626	0	real, dimension(:,:,:), allocatable :: dz_interface ! Change in position of interface due to
2627		! regridding [H ~> m or kg m-2]
2628		real :: dz_neglect ! A negligibly small vertical layer extent used in
2629		! remapping cell reconstructions [Z ~> m]
2630		real :: dz_neglect_edge ! A negligibly small vertical layer extent used in
2631		! remapping edge value calculations [Z ~> m]
2632		real :: zTopOfCell, zBottomOfCell ! Heights in Z units [Z ~> m].
2633	0	type(regridding_CS) :: regridCS ! Regridding parameters and work arrays
2634		type(remapping_CS) :: remapCS ! Remapping parameters and work arrays
2635
2636		logical :: homogenize, useALEremapping, remap_full_column, remap_general, remap_old_alg
2637		integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.
2638		integer :: remap_answer_date ! The vintage of the order of arithmetic and expressions to use
2639		! for remapping. Values below 20190101 recover the remapping
2640		! answers from 2018, while higher values use more robust
2641		! forms of the same remapping expressions.
2642		integer :: hor_regrid_answer_date ! The vintage of the order of arithmetic and expressions to use
2643		! for horizontal regridding. Values below 20190101 recover the
2644		! answers from 2018, while higher values use expressions that have
2645		! been rearranged for rotational invariance.
2646		logical :: pre_gridded
2647		logical :: separate_mixed_layer ! If true, handle the mixed layers differently.
2648		logical :: density_extrap_bug ! If true use an expression with a vertical indexing bug for
2649		! extrapolating the densities at the bottom of unstable profiles
2650		! from data when finding the initial interface locations in
2651		! layered mode from a dataset of T and S.
2652		character(len=64) :: remappingScheme
2653		logical :: om4_remap_via_sub_cells ! If true, use the OM4 remapping algorithm (only used if useALEremapping)
2654		logical :: do_conv_adj, ignore
2655		logical :: use_depth_based_time_fitler, use_adjust_interface_motion
2656		integer :: id_clock_routine, id_clock_ALE
2657
2658	0	id_clock_routine = cpu_clock_id('(Initialize from Z)', grain=CLOCK_ROUTINE)
2659	0	id_clock_ALE = cpu_clock_id('(Initialize from Z) ALE', grain=CLOCK_LOOP)
2660
2661	0	call cpu_clock_begin(id_clock_routine)
2662
2663	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
2664	0	isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed
2665	0	isg = G%isg ; ieg = G%ieg ; jsg = G%jsg ; jeg = G%jeg
2666
2667	0	PI_180=atan(1.0)/45.
2668
2669	0	if (.not.just_read) call callTree_enter(trim(mdl)//"(), MOM_state_initialization.F90")
2670	0	if (.not.just_read) call log_version(PF, mdl, version, "")
2671
2672	0	call get_param(PF, mdl, "INPUTDIR", inputdir, default=".")
2673	0	inputdir = slasher(inputdir)
2674
2675	0	eos => tv%eqn_of_state
2676
2677		call get_param(PF, mdl, "TEMP_SALT_Z_INIT_FILE", filename, &
2678		"The name of the z-space input file used to initialize "//&
2679		"temperatures (T) and salinities (S). If T and S are not "//&
2680		"in the same file, TEMP_Z_INIT_FILE and SALT_Z_INIT_FILE "//&
2681	0	"must be set.", default="temp_salt_z.nc", do_not_log=just_read)
2682		call get_param(PF, mdl, "TEMP_Z_INIT_FILE", tfilename, &
2683		"The name of the z-space input file used to initialize "//&
2684	0	"temperatures, only.", default=trim(filename), do_not_log=just_read)
2685		call get_param(PF, mdl, "SALT_Z_INIT_FILE", sfilename, &
2686		"The name of the z-space input file used to initialize "//&
2687	0	"temperatures, only.", default=trim(filename), do_not_log=just_read)
2688	0	filename = trim(inputdir)//trim(filename)
2689	0	tfilename = trim(inputdir)//trim(tfilename)
2690	0	sfilename = trim(inputdir)//trim(sfilename)
2691		call get_param(PF, mdl, "Z_INIT_FILE_PTEMP_VAR", potemp_var, &
2692		"The name of the potential temperature variable in "//&
2693	0	"TEMP_Z_INIT_FILE.", default="ptemp", do_not_log=just_read)
2694		call get_param(PF, mdl, "Z_INIT_FILE_SALT_VAR", salin_var, &
2695		"The name of the salinity variable in "//&
2696	0	"SALT_Z_INIT_FILE.", default="salt", do_not_log=just_read)
2697		call get_param(PF, mdl, "Z_INIT_HOMOGENIZE", homogenize, &
2698		"If True, then horizontally homogenize the interpolated "//&
2699	0	"initial conditions.", default=.false., do_not_log=just_read)
2700		call get_param(PF, mdl, "Z_INIT_ALE_REMAPPING", useALEremapping, &
2701		"If True, then remap straight to model coordinate from file.", &
2702	0	default=.false., do_not_log=just_read)
2703		call get_param(PF, mdl, "Z_INIT_REMAPPING_SCHEME", remappingScheme, &
2704		"The remapping scheme to use if using Z_INIT_ALE_REMAPPING "//&
2705	0	"is True.", default="PPM_IH4", do_not_log=just_read)
2706		call get_param(PF, mdl, "Z_INIT_REMAP_GENERAL", remap_general, &
2707		"If false, only initializes to z* coordinates. "//&
2708		"If true, allows initialization directly to general coordinates.", &
2709	0	default=.not.(GV%Boussinesq.or.GV%semi_Boussinesq) , do_not_log=just_read)
2710		call get_param(PF, mdl, "Z_INIT_REMAP_FULL_COLUMN", remap_full_column, &
2711		"If false, only reconstructs profiles for valid data points. "//&
2712		"If true, inserts vanished layers below the valid data.", &
2713	0	default=remap_general, do_not_log=just_read)
2714		call get_param(PF, mdl, "Z_INIT_REMAP_OLD_ALG", remap_old_alg, &
2715		"If false, uses the preferred remapping algorithm for initialization. "//&
2716		"If true, use an older, less robust algorithm for remapping.", &
2717	0	default=.false., do_not_log=just_read)
2718		call get_param(PF, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
2719		"This sets the default value for the various _ANSWER_DATE parameters.", &
2720	0	default=99991231, do_not_log=just_read)
2721		call get_param(PF, mdl, "TEMP_SALT_INIT_VERTICAL_REMAP_ONLY", pre_gridded, &
2722		"If true, initial conditions are on the model horizontal grid. " //&
2723		"Extrapolation over missing ocean values is done using an ICE-9 "//&
2724		"procedure with vertical ALE remapping .", &
2725	0	default=.false., do_not_log=just_read)
2726	0	if (useALEremapping) then
2727		call get_param(PF, mdl, "REMAPPING_ANSWER_DATE", remap_answer_date, &
2728		"The vintage of the expressions and order of arithmetic to use for remapping. "//&
2729		"Values below 20190101 result in the use of older, less accurate expressions "//&
2730		"that were in use at the end of 2018. Higher values result in the use of more "//&
2731		"robust and accurate forms of mathematically equivalent expressions.", &
2732	0	default=default_answer_date, do_not_log=just_read.or.(.not.GV%Boussinesq))
2733		call get_param(PF, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
2734	0	do_not_log=.true., default=.true.)
2735		call get_param(PF, mdl, "Z_INIT_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
2736		"If true, use the OM4 remapping-via-subcells algorithm for initialization. "//&
2737		"See REMAPPING_USE_OM4_SUBCELLS for more details. "//&
2738	0	"We recommend setting this option to false.", default=om4_remap_via_sub_cells)
2739	0	if (.not.GV%Boussinesq) remap_answer_date = max(remap_answer_date, 20230701)
2740		endif
2741		call get_param(PF, mdl, "HOR_REGRID_ANSWER_DATE", hor_regrid_answer_date, &
2742		"The vintage of the order of arithmetic for horizontal regridding. "//&
2743		"Dates before 20190101 give the same answers as the code did in late 2018, "//&
2744		"while later versions add parentheses for rotational symmetry. "//&
2745		"Dates after 20230101 use reproducing sums for global averages.", &
2746	0	default=default_answer_date, do_not_log=just_read.or.(.not.GV%Boussinesq))
2747	0	if (.not.GV%Boussinesq) hor_regrid_answer_date = max(hor_regrid_answer_date, 20230701)
2748
2749	0	if (.not.useALEremapping) then
2750		call get_param(PF, mdl, "ADJUST_THICKNESS", correct_thickness, &
2751		"If true, all mass below the bottom removed if the "//&
2752		"topography is shallower than the thickness input file "//&
2753	0	"would indicate.", default=.false., do_not_log=just_read)
2754		call get_param(PF, mdl, "THICKNESS_TOLERANCE", h_tolerance, &
2755		"A parameter that controls the tolerance when adjusting the "//&
2756		"thickness to fit the bathymetry. Used when ADJUST_THICKNESS=True.", &
2757		units="m", default=0.1, scale=US%m_to_Z, &
2758	0	do_not_log=(just_read.or..not.correct_thickness))
2759		call get_param(PF, mdl, "DZ_BOTTOM_TOLERANCE", tol_dz_bot, &
2760		"A tolerance for detecting inconsistent topography and input layer "//&
2761		"thicknesses when ADJUST_THICKNESS is false.", &
2762		units="m", default=1.0, scale=US%m_to_Z, &
2763	0	do_not_log=(just_read.or.correct_thickness))
2764
2765		call get_param(PF, mdl, "FIT_TO_TARGET_DENSITY_IC", adjust_temperature, &
2766		"If true, all the interior layers are adjusted to "//&
2767		"their target densities using mostly temperature "//&
2768		"This approach can be problematic, particularly in the "//&
2769	0	"high latitudes.", default=.true., do_not_log=just_read)
2770		call get_param(PF, mdl, "Z_INIT_SEPARATE_MIXED_LAYER", separate_mixed_layer, &
2771		"If true, distribute the topmost Z_INIT_HMIX_DEPTH of water over NKML layers, "//&
2772		"and do not correct the density of the topmost NKML+NKBL layers. Otherwise "//&
2773		"all layers are initialized based on the depths of their target densities.", &
2774	0	default=.false., do_not_log=just_read.or.(GV%nkml==0))
2775	0	if (GV%nkml == 0) separate_mixed_layer = .false.
2776		call get_param(PF, mdl, "MINIMUM_DEPTH", Hmix_default, &
2777	0	units="m", default=0.0, scale=US%m_to_Z)
2778		call get_param(PF, mdl, "Z_INIT_HMIX_DEPTH", Hmix_depth, &
2779		"The mixed layer depth in the initial conditions when Z_INIT_SEPARATE_MIXED_LAYER "//&
2780		"is set to true.", units="m", default=US%Z_to_m*Hmix_default, scale=US%m_to_Z, &
2781	0	do_not_log=(just_read .or. .not.separate_mixed_layer))
2782		! Reusing MINIMUM_DEPTH for the default mixed layer depth may be a strange choice, but
2783		! it reproduces previous answers.
2784		call get_param(PF, mdl, "DENSITY_INTERP_TOLERANCE", eps_rho, &
2785		"A small density tolerance used when finding depths in a density profile.", &
2786		units="kg m-3", default=1.0e-10, scale=US%kg_m3_to_R, &
2787	0	do_not_log=useALEremapping.or.just_read)
2788		call get_param(PF, mdl, "LAYER_Z_INIT_IC_EXTRAP_BUG", density_extrap_bug, &
2789		"If true use an expression with a vertical indexing bug for extrapolating the "//&
2790		"densities at the bottom of unstable profiles from data when finding the "//&
2791		"initial interface locations in layered mode from a dataset of T and S.", &
2792	0	default=.false., do_not_log=just_read)
2793		endif
2794		call get_param(PF, mdl, "LAND_FILL_TEMP", temp_land_fill, &
2795		"A value to use to fill in ocean temperatures on land points.", &
2796	0	units="degC", default=0.0, scale=US%degC_to_C, do_not_log=just_read)
2797		call get_param(PF, mdl, "LAND_FILL_SALIN", salt_land_fill, &
2798		"A value to use to fill in ocean salinities on land points.", &
2799	0	units="ppt", default=35.0, scale=US%ppt_to_S, do_not_log=just_read)
2800		call get_param(PF, mdl, "HORIZ_INTERP_TOL_TEMP", tol_temp, &
2801		"The tolerance in temperature changes between iterations when interpolating "//&
2802		"from an input dataset using horiz_interp_and_extrap_tracer. This routine "//&
2803		"converges slowly, so an overly small tolerance can get expensive.", &
2804	0	units="degC", default=1.0e-3, scale=US%degC_to_C, do_not_log=just_read)
2805		call get_param(PF, mdl, "HORIZ_INTERP_TOL_SALIN", tol_sal, &
2806		"The tolerance in salinity changes between iterations when interpolating "//&
2807		"from an input dataset using horiz_interp_and_extrap_tracer. This routine "//&
2808		"converges slowly, so an overly small tolerance can get expensive.", &
2809	0	units="ppt", default=1.0e-3, scale=US%ppt_to_S, do_not_log=just_read)
2810		call get_param(PF, mdl, "REGRID_USE_DEPTH_BASED_TIME_FILTER", use_depth_based_time_fitler, &
2811	0	default=.true., do_not_log=.true.)
2812		call get_param(PF, mdl, "USE_ADJUST_INTERFACE_MOTION", use_adjust_interface_motion, &
2813	0	default=.true., do_not_log=.true.)
2814
2815	0	if (just_read) then
2816	0	if ((.not.useALEremapping) .and. adjust_temperature) &
2817		! This call is just here to read and log the determine_temperature parameters
2818		call determine_temperature(tv%T, tv%S, GV%Rlay(1:nz), eos, tv%P_Ref, 0, &
2819	0	0, G, GV, US, PF, just_read=.true.)
2820	0	call cpu_clock_end(id_clock_routine)
2821	0	return ! All run-time parameters have been read, so return.
2822		endif
2823
2824	0	eps_z = GV%Angstrom_Z
2825
2826		! Read input grid coordinates for temperature and salinity field
2827		! in z-coordinate dataset. The file is REQUIRED to contain the
2828		! following:
2829		!
2830		! dimension variables:
2831		! lon (degrees_E), lat (degrees_N), depth(meters)
2832		! variables:
2833		! ptemp(lon,lat,depth) : degC, potential temperature
2834		! salt (lon,lat,depth) : ppt, salinity
2835		!
2836		! The first record will be read if there are multiple time levels.
2837		! The observation grid MUST tile the model grid. If the model grid extends
2838		! to the North/South Pole past the limits of the input data, they are extrapolated using the average
2839		! value at the northernmost/southernmost latitude.
2840
2841		call horiz_interp_and_extrap_tracer(tfilename, potemp_var, 1, &
2842		G, temp_z, mask_z, z_in, z_edges_in, missing_value_temp, &
2843		scale=US%degC_to_C, homogenize=homogenize, m_to_Z=US%m_to_Z, &
2844	0	answer_date=hor_regrid_answer_date, ongrid=pre_gridded, tr_iter_tol=tol_temp)
2845
2846		call horiz_interp_and_extrap_tracer(sfilename, salin_var, 1, &
2847		G, salt_z, mask_z, z_in, z_edges_in, missing_value_salt, &
2848		scale=US%ppt_to_S, homogenize=homogenize, m_to_Z=US%m_to_Z, &
2849	0	answer_date=hor_regrid_answer_date, ongrid=pre_gridded, tr_iter_tol=tol_sal)
2850
2851	0	kd = size(z_in,1)
2852
2853		! Convert the sign convention of Z_edges_in.
2854	0	do k=1,size(Z_edges_in,1) ; Z_edges_in(k) = -Z_edges_in(k) ; enddo
2855
2856		! Convert T&S to Absolute Salinity and Conservative Temperature if using TEOS10 or NEMO
2857	0	call convert_temp_salt_for_TEOS10(temp_z, salt_z, G%HI, kd, mask_z, eos)
2858
2859	0	do j=js,je ; do i=is,ie
2860	0	Z_bottom(i,j) = -depth_tot(i,j)
2861		enddo ; enddo
2862
2863		! Done with horizontal interpolation.
2864		! Now remap to model coordinates
2865	0	if (useALEremapping) then
2866	0	call cpu_clock_begin(id_clock_ALE)
2867	0	nkd = max(GV%ke, kd)
2868
2869		! Build the source grid and copy data onto model-shaped arrays with vanished layers
2870	0	allocate( tmp_mask_in(isd:ied,jsd:jed,nkd), source=0.0 )
2871	0	allocate( dz1(isd:ied,jsd:jed,nkd), source=0.0 )
2872	0	allocate( h1(isd:ied,jsd:jed,nkd), source=0.0 )
2873	0	allocate( tmpT1dIn(isd:ied,jsd:jed,nkd), source=0.0 )
2874	0	allocate( tmpS1dIn(isd:ied,jsd:jed,nkd), source=0.0 )
2875	0	do j = js, je ; do i = is, ie
2876	0	if (G%mask2dT(i,j) > 0.) then
2877	0	zTopOfCell = 0. ; zBottomOfCell = 0.
2878	0	tmp_mask_in(i,j,1:kd) = mask_z(i,j,:)
2879	0	do k = 1, nkd
2880	0	if ((tmp_mask_in(i,j,k) > 0.) .and. (k <= kd)) then
2881	0	zBottomOfCell = max( z_edges_in(k+1), Z_bottom(i,j))
2882	0	tmpT1dIn(i,j,k) = temp_z(i,j,k)
2883	0	tmpS1dIn(i,j,k) = salt_z(i,j,k)
2884	0	elseif (k>1) then
2885	0	zBottomOfCell = Z_bottom(i,j)
2886	0	tmpT1dIn(i,j,k) = tmpT1dIn(i,j,k-1)
2887	0	tmpS1dIn(i,j,k) = tmpS1dIn(i,j,k-1)
2888		else ! This next block should only ever be reached over land
2889	0	tmpT1dIn(i,j,k) = temp_land_fill
2890	0	tmpS1dIn(i,j,k) = salt_land_fill
2891		endif
2892	0	dz1(i,j,k) = (zTopOfCell - zBottomOfCell)
2893	0	zTopOfCell = zBottomOfCell ! Bottom becomes top for next value of k
2894		enddo
2895	0	dz1(i,j,kd) = dz1(i,j,kd) + max(0., zTopOfCell - Z_bottom(i,j) )
2896		! The max here is in case the data data is shallower than model
2897		endif ! mask2dT
2898		enddo ; enddo
2899	0	deallocate( tmp_mask_in )
2900
2901		! Convert input thicknesses to units of H. In non-Boussinesq mode this is done by inverting
2902		! integrals of specific volume in pressure, so it can be expensive.
2903	0	tv_loc = tv
2904	0	tv_loc%T => tmpT1dIn
2905	0	tv_loc%S => tmpS1dIn
2906	0	GV_loc = GV
2907	0	GV_loc%ke = nkd
2908	0	call dz_to_thickness(dz1, tv_loc, h1, G, GV_loc, US)
2909
2910		! Build the target grid (and set the model thickness to it)
2911
2912	0	call ALE_initRegridding( G, GV, US, G%max_depth, PF, mdl, regridCS ) ! sets regridCS
2913	0	if (remap_general) then
2914	0	dz_neglect = set_h_neglect(GV, remap_answer_date, dz_neglect_edge)
2915		else
2916	0	dz_neglect = set_dz_neglect(GV, US, remap_answer_date, dz_neglect_edge)
2917		endif
2918		call initialize_remapping( remapCS, remappingScheme, boundary_extrapolation=.false., &
2919		om4_remap_via_sub_cells=om4_remap_via_sub_cells, answer_date=remap_answer_date, &
2920	0	h_neglect=dz_neglect, h_neglect_edge=dz_neglect_edge)
2921
2922		! Now remap from source grid to target grid, first setting reconstruction parameters
2923	0	if (remap_general) then
2924		call set_regrid_params( regridCS, min_thickness=0., &
2925		use_adjust_interface_motion=use_adjust_interface_motion, &
2926	0	use_depth_based_time_filter=use_depth_based_time_fitler)
2927	0	allocate( dz_interface(isd:ied,jsd:jed,nkd+1), source=0.) ! Need for argument to regridding_main() but is not used
2928
2929	0	call regridding_preadjust_reqs(regridCS, do_conv_adj, ignore)
2930	0	if (do_conv_adj) call convective_adjustment(G, GV_loc, h1, tv_loc)
2931		call regridding_main( remapCS, regridCS, G, GV_loc, US, h1, tv_loc, h, dz_interface, &
2932	0	frac_shelf_h=frac_shelf_h )
2933
2934	0	deallocate( dz_interface )
2935
2936		call ALE_remap_scalar(remapCS, G, GV, nkd, h1, tmpT1dIn, h, tv%T, all_cells=remap_full_column, &
2937	0	old_remap=remap_old_alg )
2938		call ALE_remap_scalar(remapCS, G, GV, nkd, h1, tmpS1dIn, h, tv%S, all_cells=remap_full_column, &
2939	0	old_remap=remap_old_alg )
2940		else
2941		! This is the old way of initializing to z* coordinates only
2942	0	allocate( hTarget(nz) )
2943	0	hTarget = getCoordinateResolution( regridCS )
2944	0	do j = js, je ; do i = is, ie
2945	0	dz(i,j,:) = 0.
2946	0	if (G%mask2dT(i,j) > 0.) then
2947		! Build the target grid combining hTarget and topography
2948	0	zTopOfCell = 0. ; zBottomOfCell = 0.
2949	0	do k = 1, nz
2950	0	zBottomOfCell = max( zTopOfCell - hTarget(k), Z_bottom(i,j))
2951	0	dz(i,j,k) = zTopOfCell - zBottomOfCell
2952	0	zTopOfCell = zBottomOfCell ! Bottom becomes top for next value of k
2953		enddo
2954		else
2955	0	dz(i,j,:) = 0.
2956		endif ! mask2dT
2957		enddo ; enddo
2958	0	deallocate( hTarget )
2959
2960	0	dz_neglect = set_dz_neglect(GV, US, remap_answer_date, dz_neglect_edge)
2961		call ALE_remap_scalar(remapCS, G, GV, nkd, dz1, tmpT1dIn, dz, tv%T, all_cells=remap_full_column, &
2962	0	old_remap=remap_old_alg)
2963		call ALE_remap_scalar(remapCS, G, GV, nkd, dz1, tmpS1dIn, dz, tv%S, all_cells=remap_full_column, &
2964	0	old_remap=remap_old_alg)
2965
2966	0	if (GV%Boussinesq .or. GV%semi_Boussinesq) then
2967		! This is a simple conversion of the target grid to thickness units that is not
2968		! appropriate in non-Boussinesq mode.
2969	0	call dz_to_thickness_simple(dz, h, G, GV, US)
2970		else
2971		! Convert dz into thicknesses in units of H using the equation of state as appropriate.
2972	0	call dz_to_thickness(dz, tv, h, G, GV, US)
2973		endif
2974		endif
2975
2976	0	deallocate( dz1 )
2977	0	deallocate( h1 )
2978	0	deallocate( tmpT1dIn )
2979	0	deallocate( tmpS1dIn )
2980
2981	0	call cpu_clock_end(id_clock_ALE)
2982
2983		else ! remap to isopycnal layer space
2984
2985		! Next find interface positions using local arrays
2986		! nlevs contains the number of valid data points in each column
2987	0	nlevs = int(sum(mask_z,dim=3))
2988
2989		! Rb contains the layer interface densities
2990	0	allocate(Rb(nz+1))
2991	0	do k=2,nz ; Rb(k) = 0.5*(GV%Rlay(k-1)+GV%Rlay(k)) ; enddo
2992	0	Rb(1) = 0.0
2993	0	if (nz>1) then
2994	0	Rb(nz+1) = 2.0*GV%Rlay(nz) - GV%Rlay(nz-1)
2995		else
2996	0	Rb(nz+1) = 2.0 * GV%Rlay(1)
2997		endif
2998
2999	0	nkml = 0 ; if (separate_mixed_layer) nkml = GV%nkml
3000
3001	0	press(:) = tv%P_Ref
3002	0	EOSdom(:) = EOS_domain(G%HI)
3003	0	allocate(rho_z(isd:ied,jsd:jed,kd))
3004	0	do k=1,kd ; do j=js,je
3005	0	call calculate_density(temp_z(:,j,k), salt_z(:,j,k), press, rho_z(:,j,k), eos, EOSdom)
3006		enddo ; enddo
3007
3008		call find_interfaces(rho_z, z_in, kd, Rb, Z_bottom, zi, G, GV, US, nlevs, nkml, &
3009	0	Hmix_depth, eps_z, eps_rho, density_extrap_bug)
3010
3011	0	deallocate(rho_z, Rb)
3012
3013	0	dz(:,:,:) = 0.0
3014	0	if (correct_thickness) then
3015	0	call adjustEtaToFitBathymetry(G, GV, US, zi, dz, h_tolerance, dZ_ref_eta=G%Z_ref)
3016		else
3017	0	do k=nz,1,-1 ; do j=js,je ; do i=is,ie
3018	0	if (zi(i,j,K) < (zi(i,j,K+1) + GV%Angstrom_Z)) then
3019	0	zi(i,j,K) = zi(i,j,K+1) + GV%Angstrom_Z
3020	0	dz(i,j,k) = GV%Angstrom_Z
3021		else
3022	0	dz(i,j,k) = zi(i,j,K) - zi(i,j,K+1)
3023		endif
3024		enddo ; enddo ; enddo
3025	0	inconsistent = 0
3026	0	do j=js,je ; do i=is,ie
3027	0	if (abs(zi(i,j,nz+1) - Z_bottom(i,j)) > tol_dz_bot) &
3028	0	inconsistent = inconsistent + 1
3029		enddo ; enddo
3030	0	call sum_across_PEs(inconsistent)
3031
3032	0	if ((inconsistent > 0) .and. (is_root_pe())) then
3033		write(mesg, '("Thickness initial conditions are inconsistent ",'// &
3034	0	'"with topography in ",I0," places.")') inconsistent
3035	0	call MOM_error(WARNING, mesg)
3036		endif
3037		endif
3038
3039	0	call tracer_z_init_array(temp_z, z_edges_in, kd, zi, temp_land_fill, G, nz, nlevs, eps_z, tv%T)
3040	0	call tracer_z_init_array(salt_z, z_edges_in, kd, zi, salt_land_fill, G, nz, nlevs, eps_z, tv%S)
3041
3042	0	if (homogenize) then
3043		! Horizontally homogenize data to produce perfectly "flat" initial conditions
3044	0	do k=1,nz
3045	0	call homogenize_field(tv%T(:,:,k), G, tmp_scale=US%C_to_degC, answer_date=hor_regrid_answer_date)
3046	0	call homogenize_field(tv%S(:,:,k), G, tmp_scale=US%S_to_ppt, answer_date=hor_regrid_answer_date)
3047		enddo
3048		endif
3049
3050	0	if (adjust_temperature) then
3051		! Finally adjust to target density
3052	0	ks = 1 ; if (separate_mixed_layer) ks = GV%nk_rho_varies + 1
3053		call determine_temperature(tv%T, tv%S, GV%Rlay(1:nz), eos, tv%P_Ref, niter, &
3054	0	ks, G, GV, US, PF, just_read)
3055		endif
3056
3057		! Now convert dz into thicknesses in units of H.
3058	0	call dz_to_thickness(dz, tv, h, G, GV, US)
3059
3060		endif ! useALEremapping
3061
3062	0	deallocate(z_in, z_edges_in, temp_z, salt_z, mask_z)
3063
3064	0	call pass_var(h, G%Domain)
3065	0	call pass_var(tv%T, G%Domain)
3066	0	call pass_var(tv%S, G%Domain)
3067
3068	0	call callTree_leave(trim(mdl)//'()')
3069	0	call cpu_clock_end(id_clock_routine)
3070
3071	0	end subroutine MOM_temp_salt_initialize_from_Z
3072
3073
3074		!> Find interface positions corresponding to interpolated depths in a density profile
3075	0	subroutine find_interfaces(rho, zin, nk_data, Rb, Z_bot, zi, G, GV, US, nlevs, nkml, hml, &
3076		eps_z, eps_rho, density_extrap_bug)
3077		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
3078		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
3079		integer, intent(in) :: nk_data !< The number of levels in the input data
3080		real, dimension(SZI_(G),SZJ_(G),nk_data), &
3081		intent(in) :: rho !< Potential density in z-space [R ~> kg m-3]
3082		real, dimension(nk_data), intent(in) :: zin !< Input data levels [Z ~> m].
3083		real, dimension(SZK_(GV)+1), intent(in) :: Rb !< target interface densities [R ~> kg m-3]
3084		real, dimension(SZI_(G),SZJ_(G)), &
3085		intent(in) :: Z_bot !< The (usually negative) height of the seafloor
3086		!! relative to the surface [Z ~> m].
3087		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
3088		intent(out) :: zi !< The returned interface heights [Z ~> m]
3089		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
3090		integer, dimension(SZI_(G),SZJ_(G)), &
3091		intent(in) :: nlevs !< number of valid points in each column
3092		integer, intent(in) :: nkml !< number of mixed layer pieces to distribute over
3093		!! a depth of hml.
3094		real, intent(in) :: hml !< mixed layer depth [Z ~> m].
3095		real, intent(in) :: eps_z !< A negligibly small layer thickness [Z ~> m].
3096		real, intent(in) :: eps_rho !< A negligibly small density difference [R ~> kg m-3].
3097		logical, intent(in) :: density_extrap_bug !< If true use an expression with an
3098		!! indexing bug for projecting the densities at
3099		!! the bottom of unstable profiles from data when
3100		!! finding the initial interface locations in
3101		!! layered mode from a dataset of T and S.
3102
3103		! Local variables
3104	0	real, dimension(nk_data) :: rho_ ! A column of densities [R ~> kg m-3]
3105	0	real, dimension(SZK_(GV)+1) :: zi_ ! A column interface heights (negative downward) [Z ~> m].
3106		real :: slope ! The rate of change of height with density [Z R-1 ~> m4 kg-1]
3107		real :: drhodz ! A local vertical density gradient [R Z-1 ~> kg m-4]
3108		real, parameter :: zoff = 0.999 ! A small fractional adjustment to the density differences [nondim]
3109		logical :: unstable ! True if the column is statically unstable anywhere.
3110		integer :: nlevs_data ! The number of data values in a column.
3111		logical :: work_down ! This indicates whether this pass goes up or down the water column.
3112		integer :: k_int, lo_int, hi_int, mid
3113		integer :: i, j, k, is, ie, js, je, nz
3114
3115	0	is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec ; nz = GV%ke
3116
3117	0	zi(:,:,:) = 0.0
3118
3119	0	do j=js,je ; do i=is,ie
3120	0	nlevs_data = nlevs(i,j)
3121	0	do k=1,nlevs_data ; rho_(k) = rho(i,j,k) ; enddo
3122
3123	0	unstable = .true.
3124	0	work_down = .true.
3125	0	do while (unstable)
3126		! Modify the input profile until it no longer has densities that decrease with depth.
3127	0	unstable = .false.
3128	0	if (work_down) then
3129	0	do k=2,nlevs_data-1 ; if (rho_(k) - rho_(k-1) < 0.0) then
3130	0	if (k == 2) then
3131	0	rho_(k-1) = rho_(k) - eps_rho
3132		else
3133	0	drhodz = (rho_(k+1)-rho_(k-1)) / (zin(k+1)-zin(k-1))
3134	0	if (drhodz < 0.0) unstable = .true.
3135	0	rho_(k) = rho_(k-1) + drhodzzoff(zin(k)-zin(k-1))
3136		endif
3137		endif ; enddo
3138	0	work_down = .false.
3139		else
3140	0	do k=nlevs_data-1,2,-1 ; if (rho_(k+1) - rho_(k) < 0.0) then
3141	0	if (k == nlevs_data-1) then
3142	0	if (density_extrap_bug) then
3143	0	rho_(k+1) = rho_(k-1) + eps_rho
3144		else
3145	0	rho_(k+1) = rho_(k) + eps_rho
3146		endif
3147		else
3148	0	drhodz = (rho_(k+1)-rho_(k-1)) / (zin(k+1)-zin(k-1))
3149	0	if (drhodz < 0.0) unstable = .true.
3150	0	rho_(k) = rho_(k+1) - drhodz*(zin(k+1)-zin(k))
3151		endif
3152		endif ; enddo
3153	0	work_down = .true.
3154		endif
3155		enddo
3156
3157		! Find and store the interface depths.
3158	0	zi_(1) = 0.0
3159	0	if (nlevs_data < 1) then
3160		! There is no data to use, so set the interfaces at the bottom.
3161	0	do K=2,nz ; zi_(K) = Z_bot(i,j) ; enddo
3162	0	elseif (nlevs_data == 1) then
3163		! There is data for only one input layer, so set the interfaces at the bottom or top,
3164		! depending on how their target densities compare with the one data point.
3165	0	do K=2,nz
3166	0	if (Rb(K) < rho_(1)) then ; zi_(K) = 0.0
3167	0	else ; zi_(K) = Z_bot(i,j) ; endif
3168		enddo
3169		else
3170	0	do K=2,nz
3171		! Find the value of k_int in the list of rho_ where rho_(k_int) <= Rb(K) < rho_(k_int+1).
3172		! This might be made a little faster by exploiting the fact that Rb is
3173		! monotonically increasing and not resetting lo_int back to 1 inside the K loop.
3174	0	lo_int = 1 ; hi_int = nlevs_data
3175	0	do while (lo_int < hi_int)
3176	0	mid = (lo_int+hi_int) / 2
3177	0	if (Rb(K) < rho_(mid)) then ; hi_int = mid
3178	0	else ; lo_int = mid+1 ; endif
3179		enddo
3180	0	k_int = max(1, lo_int-1)
3181
3182		! Linearly interpolate to find the depth, zi_, where Rb would be found.
3183	0	slope = (zin(k_int+1) - zin(k_int)) / max(rho_(k_int+1) - rho_(k_int), eps_rho)
3184	0	zi_(K) = -1.0(zin(k_int) + slope(Rb(K)-rho_(k_int)))
3185	0	zi_(K) = min(max(zi_(K), Z_bot(i,j)), -1.0*hml)
3186		enddo
3187		endif
3188	0	zi_(nz+1) = Z_bot(i,j)
3189	0	if (nkml > 0) then ; do K=2,nkml+1
3190	0	zi_(K) = max(hml*((1.0-real(K))/real(nkml)), Z_bot(i,j))
3191		enddo ; endif
3192	0	do K=nz,max(nkml+2,2),-1
3193	0	if (zi_(K) < zi_(K+1) + eps_Z) zi_(K) = zi_(K+1) + eps_Z
3194	0	if (zi_(K) > -1.0hml) zi_(K) = max(-1.0hml, Z_bot(i,j))
3195		enddo
3196
3197	0	do K=1,nz+1
3198	0	zi(i,j,K) = zi_(K)
3199		enddo
3200		enddo ; enddo ! i- and j- loops
3201
3202	0	end subroutine find_interfaces
3203
3204		!> Run simple unit tests
3205	0	subroutine MOM_state_init_tests(G, GV, US, tv)
3206		type(ocean_grid_type), intent(inout) :: G !< The ocean's grid structure.
3207		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure.
3208		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
3209		type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure.
3210
3211		! Local variables
3212		integer, parameter :: nk=5
3213		real, dimension(nk) :: T, T_t, T_b ! Temperatures [C ~> degC]
3214		real, dimension(nk) :: S, S_t, S_b ! Salinities [S ~> ppt]
3215		real, dimension(nk) :: rho ! Layer density [R ~> kg m-3]
3216		real, dimension(nk) :: h ! Layer thicknesses [H ~> m or kg m-2]
3217		real, dimension(nk) :: z ! Height of layer center [Z ~> m]
3218		real, dimension(nk+1) :: e ! Interface heights [Z ~> m]
3219		real :: T_ref ! A reference temperature [C ~> degC]
3220		real :: S_ref ! A reference salinity [S ~> ppt]
3221		real :: P_tot, P_t, P_b ! Pressures [R L2 T-2 ~> Pa]
3222		real :: z_out ! Output height [Z ~> m]
3223		real :: I_z_scale ! The inverse of the height scale for prescribed gradients [Z-1 ~> m-1]
3224		real :: z_tol ! The tolerance with which to find the depth matching a specified pressure [Z ~> m].
3225		integer :: k
3226		type(remapping_CS), pointer :: remap_CS => NULL()
3227
3228	0	I_z_scale = 1.0 / (500.0*US%m_to_Z)
3229	0	do k = 1, nk
3230	0	h(k) = 100.0*GV%m_to_H
3231		enddo
3232	0	e(1) = 0.
3233	0	do K = 1, nk
3234	0	e(K+1) = e(K) - GV%H_to_Z * h(k)
3235		enddo
3236	0	P_tot = 0.
3237	0	T_ref = 20.0*US%degC_to_C
3238	0	S_ref = 35.0*US%ppt_to_S
3239	0	z_tol = 1.0e-5*US%m_to_Z
3240	0	do k = 1, nk
3241	0	z(k) = 0.5 * ( e(K) + e(K+1) )
3242	0	T_t(k) = T_ref + (0. * I_z_scale) * e(k)
3243	0	T(k) = T_ref + (0. * I_z_scale)*z(k)
3244	0	T_b(k) = T_ref + (0. * I_z_scale)*e(k+1)
3245	0	S_t(k) = S_ref - (0. * I_z_scale)*e(k)
3246	0	S(k) = S_ref + (0. * I_z_scale)*z(k)
3247	0	S_b(k) = S_ref - (0. * I_z_scale)*e(k+1)
3248		call calculate_density(0.5(T_t(k)+T_b(k)), 0.5(S_t(k)+S_b(k)), -GV%Rho0GV%g_Earthz(k), &
3249	0	rho(k), tv%eqn_of_state)
3250	0	P_tot = P_tot + GV%g_Earth * rho(k) * GV%H_to_Z*h(k)
3251		enddo
3252
3253	0	P_t = 0.
3254	0	do k = 1, nk
3255		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, &
3256		GV%Rho0, GV%g_Earth, tv%eqn_of_state, US, P_b, z_out, z_tol=z_tol, &
3257	0	frac_dp_bugfix=.false.)
3258	0	write(0,) k, US%RL2_T2_to_PaP_t, US%RL2_T2_to_PaP_b, 0.5US%RL2_T2_to_Pa*P_tot, &
3259	0	US%Z_to_me(K), US%Z_to_me(K+1), US%Z_to_m*z_out
3260	0	P_t = P_b
3261		enddo
3262	0	write(0,) US%RL2_T2_to_PaP_b, US%RL2_T2_to_Pa*P_tot
3263
3264	0	write(0,*) ''
3265	0	write(0,*) ' ==================================================================== '
3266	0	write(0,*) ''
3267	0	write(0,) GV%H_to_mh(:)
3268
3269		! For consistency with the usual call, add the following:
3270		! if (use_remapping) then
3271		! allocate(remap_CS)
3272		! call initialize_remapping(remap_CS, 'PLM', boundary_extrapolation=.true., &
3273		! h_neglect=GV%H_subroundoff, h_neglect_edge=GV%H_subroundoff)
3274		! endif
3275		call cut_off_column_top(nk, tv, GV, US, GV%g_Earth, -e(nk+1), GV%Angstrom_H, &
3276		T, T_t, T_b, S, S_t, S_b, 0.5*P_tot, h, remap_CS, z_tol=z_tol, &
3277	0	frac_dp_bugfix=.false.)
3278	0	write(0,) GV%H_to_mh(:)
3279	0	if (associated(remap_CS)) deallocate(remap_CS)
3280
3281	0	end subroutine MOM_state_init_tests
3282
3283		end module MOM_state_initialization