src/ALE/MOM_regridding.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		!> Generates vertical grids as part of the ALE algorithm
6		module MOM_regridding
7
8		use MOM_error_handler, only : MOM_error, FATAL, WARNING, NOTE, assert
9		use MOM_file_parser, only : param_file_type, get_param, log_param
10		use MOM_io, only : file_exists, field_exists, field_size, MOM_read_data
11		use MOM_io, only : read_variable
12		use MOM_io, only : vardesc, var_desc, SINGLE_FILE
13		use MOM_io, only : MOM_netCDF_file, MOM_field
14		use MOM_io, only : create_MOM_file, MOM_write_field
15		use MOM_io, only : verify_variable_units, slasher
16		use MOM_unit_scaling, only : unit_scale_type
17		use MOM_variables, only : ocean_grid_type, thermo_var_ptrs
18		use MOM_verticalGrid, only : verticalGrid_type
19		use MOM_EOS, only : EOS_type, calculate_density
20		use MOM_domains, only : max_across_PEs, pass_var
21		use MOM_string_functions, only : uppercase, extractWord, extract_integer, extract_real
22
23		use MOM_remapping, only : remapping_CS
24		use regrid_consts, only : state_dependent, coordinateUnits
25		use regrid_consts, only : coordinateMode, DEFAULT_COORDINATE_MODE
26		use regrid_consts, only : REGRIDDING_LAYER, REGRIDDING_ZSTAR
27		use regrid_consts, only : REGRIDDING_RHO, REGRIDDING_SIGMA
28		use regrid_consts, only : REGRIDDING_ARBITRARY, REGRIDDING_SIGMA_SHELF_ZSTAR
29		use regrid_consts, only : REGRIDDING_HYCOM1, REGRIDDING_HYBGEN, REGRIDDING_ADAPTIVE
30		use regrid_interp, only : interp_CS_type
31		use regrid_interp, only : set_interp_scheme, set_interp_extrap, set_interp_answer_date
32
33		use coord_zlike, only : zlike_CS
34		use coord_zlike, only : init_coord_zlike, set_zlike_params, build_zstar_column, end_coord_zlike
35		use coord_sigma, only : sigma_CS
36		use coord_sigma, only : init_coord_sigma, set_sigma_params, build_sigma_column, end_coord_sigma
37		use coord_rho, only : init_coord_rho, rho_CS, set_rho_params, build_rho_column, end_coord_rho
38		use coord_rho, only : old_inflate_layers_1d
39		use coord_hycom, only : hycom_CS
40		use coord_hycom, only : init_coord_hycom, set_hycom_params, build_hycom1_column, end_coord_hycom
41		use coord_hycom, only : init_3d_coord_hycom
42		use coord_adapt, only : adapt_CS
43		use coord_adapt, only : init_coord_adapt, set_adapt_params, build_adapt_column, end_coord_adapt
44		use MOM_hybgen_regrid, only : hybgen_regrid, hybgen_regrid_CS, init_hybgen_regrid, end_hybgen_regrid
45		use MOM_hybgen_regrid, only : write_Hybgen_coord_file
46
47		implicit none ; private
48
49		#include <MOM_memory.h>
50
51		! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
52		! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
53		! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
54		! vary with the Boussinesq approximation, the Boussinesq variant is given first.
55
56		!> Regridding control structure
57		type, public :: regridding_CS ; private
58
59		!> This array is set by function setCoordinateResolution()
60		!! It contains the "resolution" or delta coordinate of the target
61		!! coordinate. It has the units of the target coordinate, e.g.
62		!! [Z ~> m] for z*, [nondim] for sigma, etc.
63		real, dimension(:), allocatable :: coordinateResolution
64
65		!> This is a scaling factor that restores coordinateResolution to values in
66		!! the natural units for output, perhaps [nondim]
67		real :: coord_scale = 1.0
68
69		!> This array is set by function set_target_densities()
70		!! This array is the nominal coordinate of interfaces and is the
71		!! running sum of coordinateResolution, in [R ~> kg m-3]. i.e.
72		!! target_density(k+1) = coordinateResolution(k) + coordinateResolution(k)
73		!! It is only used in "rho" or "Hycom" mode.
74		real, dimension(:), allocatable :: target_density
75
76		!> A flag to indicate that the target_density arrays has been filled with data.
77		logical :: target_density_set = .false.
78
79		!> Nominal HYCOM1 3D near-surface resolution [Z ~> m]
80		real, allocatable, dimension(:,:,:) :: coordinateResolution_3d
81
82		!> Nominal HYCOM1 3D density of interfaces [R ~> kg m-3]
83		real, allocatable, dimension(:,:,:) :: target_density_3d
84
85		!> This array is set by function set_regrid_max_depths()
86		!! It specifies the maximum depth that every interface is allowed to take [H ~> m or kg m-2].
87		real, dimension(:), allocatable :: max_interface_depths
88
89		!> This array is set by function set_regrid_max_thickness()
90		!! It specifies the maximum depth that every interface is allowed to take [H ~> m or kg m-2].
91		real, dimension(:), allocatable :: max_layer_thickness
92
93		integer :: nk !< Number of layers/levels in generated grid
94
95		!> Indicates which grid to use in the vertical (z*, sigma, target interface
96		!! densities)
97		integer :: regridding_scheme
98
99		!> Interpolation control structure
100		type(interp_CS_type) :: interp_CS
101
102		!> Minimum thickness allowed when building the new grid through regridding [H ~> m or kg m-2].
103		real :: min_thickness
104
105		!> If true, call adjust_interface_motion() after initial grid generation
106		logical :: use_adjust_interface_motion
107
108		!> Reference pressure for potential density calculations [R L2 T-2 ~> Pa]
109		real :: ref_pressure = 2.e7
110
111		!> If true, always pass through the depth-based time filtering that uses CS%old_grid_weight
112		!! If false, allows bypassing of the call if CS%old_grid_weight==0
113		logical :: use_depth_based_time_filter
114
115		!> Weight given to old coordinate when blending between new and old grids [nondim]
116		!! Used only below depth_of_time_filter_shallow, with a cubic variation
117		!! from zero to full effect between depth_of_time_filter_shallow and
118		!! depth_of_time_filter_deep.
119		real :: old_grid_weight = 0.
120
121		!> Depth above which no time-filtering of grid is applied [H ~> m or kg m-2]
122		real :: depth_of_time_filter_shallow = 0.
123
124		!> Depth below which time-filtering of grid is applied at full effect [H ~> m or kg m-2]
125		real :: depth_of_time_filter_deep = 0.
126
127		!> Fraction (between 0 and 1) of compressibility to add to potential density
128		!! profiles when interpolating for target grid positions [nondim]
129		real :: compressibility_fraction = 0.
130
131		!> If true, each interface is given a maximum depth based on a rescaling of
132		!! the indexing of coordinateResolution.
133		logical :: set_maximum_depths = .false.
134
135		!> If true, integrate for interface positions from the top downward.
136		!! If false, integrate from the bottom upward, as does the rest of the model.
137		logical :: integrate_downward_for_e = .true.
138
139		!> The vintage of the order of arithmetic and expressions to use for remapping.
140		!! Values below 20190101 recover the remapping answers from 2018.
141		!! Higher values use more robust forms of the same remapping expressions.
142		integer :: remap_answer_date = 99991231
143
144		logical :: use_hybgen_unmix = .false. !< If true, use the hybgen unmixing code before remapping
145
146		type(zlike_CS), pointer :: zlike_CS => null() !< Control structure for z-like coordinate generator
147		type(sigma_CS), pointer :: sigma_CS => null() !< Control structure for sigma coordinate generator
148		type(rho_CS), pointer :: rho_CS => null() !< Control structure for rho coordinate generator
149		type(hycom_CS), pointer :: hycom_CS => null() !< Control structure for hybrid coordinate generator
150		type(adapt_CS), pointer :: adapt_CS => null() !< Control structure for adaptive coordinate generator
151		type(hybgen_regrid_CS), pointer :: hybgen_CS => NULL() !< Control structure for hybgen regridding
152
153		end type
154
155		! The following routines are visible to the outside world
156		public initialize_regridding, end_regridding, regridding_main
157		public regridding_preadjust_reqs, convective_adjustment
158		public inflate_vanished_layers_old, check_grid_column
159		public set_regrid_params, get_regrid_size, write_regrid_file
160		public uniformResolution, setCoordinateResolution
161		public set_target_densities_from_GV, set_target_densities
162		public set_regrid_max_depths, set_regrid_max_thickness
163		public getCoordinateResolution, getCoordinateInterfaces
164		public getCoordinateUnits, getCoordinateShortName, getStaticThickness
165		public DEFAULT_COORDINATE_MODE
166		public set_h_neglect, set_dz_neglect
167		public get_zlike_CS, get_sigma_CS, get_rho_CS
168
169		!> Documentation for coordinate options
170		character(len=*), parameter, public :: regriddingCoordinateModeDoc = &
171		" LAYER - Isopycnal or stacked shallow water layers\n"//&
172		" ZSTAR, Z* - stretched geopotential z*\n"//&
173		" SIGMA_SHELF_ZSTAR - stretched geopotential z* ignoring shelf\n"//&
174		" SIGMA - terrain following coordinates\n"//&
175		" RHO - continuous isopycnal\n"//&
176		" HYCOM1 - HyCOM-like hybrid coordinate\n"//&
177		" HYBGEN - Hybrid coordinate from the Hycom hybgen code\n"//&
178		" ADAPTIVE - optimize for smooth neutral density surfaces"
179
180		!> Documentation for regridding interpolation schemes
181		character(len=*), parameter, public :: regriddingInterpSchemeDoc = &
182		" P1M_H2 (2nd-order accurate)\n"//&
183		" P1M_H4 (2nd-order accurate)\n"//&
184		" P1M_IH4 (2nd-order accurate)\n"//&
185		" PLM (2nd-order accurate)\n"//&
186		" PPM_CW (3rd-order accurate)\n"//&
187		" PPM_H4 (3rd-order accurate)\n"//&
188		" PPM_IH4 (3rd-order accurate)\n"//&
189		" P3M_IH4IH3 (4th-order accurate)\n"//&
190		" P3M_IH6IH5 (4th-order accurate)\n"//&
191		" PQM_IH4IH3 (4th-order accurate)\n"//&
192		" PQM_IH6IH5 (5th-order accurate)"
193
194		!> Default interpolation scheme
195		character(len=*), parameter, public :: regriddingDefaultInterpScheme = "P1M_H2"
196		!> Default mode for boundary extrapolation
197		logical, parameter, public :: regriddingDefaultBoundaryExtrapolation = .false.
198		!> Default minimum thickness for some coordinate generation modes [m]
199		real, parameter, public :: regriddingDefaultMinThickness = 1.e-3
200
201		!> Maximum length of parameters
202		integer, parameter :: MAX_PARAM_LENGTH = 120
203
204		#undef __DO_SAFETY_CHECKS__
205
206		contains
207
208		!> Initialization and configures a regridding control structure based on customizable run-time parameters
209	2	subroutine initialize_regridding(CS, G, GV, US, max_depth, param_file, mdl, &
210		coord_mode, param_prefix, param_suffix)
211		type(regridding_CS), intent(inout) :: CS !< Regridding control structure
212		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
213		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
214		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
215		real, intent(in) :: max_depth !< The maximum depth of the ocean [Z ~> m].
216		type(param_file_type), intent(in) :: param_file !< Parameter file
217		character(len=*), intent(in) :: mdl !< Name of calling module.
218		character(len=*), intent(in) :: coord_mode !< Coordinate mode
219		character(len=*), intent(in) :: param_prefix !< String to prefix to parameter names.
220		!! If empty, causes main model parameters to be used.
221		character(len=*), intent(in) :: param_suffix !< String to append to parameter names.
222
223		! Local variables
224		integer :: ke ! Number of levels
225		integer :: n_sigma ! Number of shallow dz's, for HYBRID_MAP or HYBRID_3D
226		integer :: np ! Number of profiles, for HYBRID_MAP
227		integer :: nceiling ! ceiling of map index, for HYBRID_MAP
228		integer :: nfloor ! floor of map index, for HYBRID_MAP
229		real :: nfrac ! fraction of map index, for HYBRID_MAP [nondim]
230		character(len=80) :: string, string2, varName ! Temporary strings
231		character(len=40) :: coord_units, coord_res_param ! Temporary strings
232		character(len=MAX_PARAM_LENGTH) :: param_name
233		character(len=200) :: inputdir, fileName, longString
234		character(len=320) :: message ! Temporary strings
235		character(len=12) :: expected_units, alt_units ! Temporary strings
236		logical :: tmpLogical, do_sum, main_parameters
237		logical :: coord_is_state_dependent, ierr
238		integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.
239		integer :: remap_answer_date ! The vintage of the remapping expressions to use.
240		integer :: regrid_answer_date ! The vintage of the regridding expressions to use.
241		real :: tmpReal ! A temporary variable used in setting other variables [various]
242		real :: P_Ref ! The coordinate variable reference pression [R L2 T-2 ~> Pa]
243		real :: maximum_depth ! The maximum depth of the ocean [m] (not in Z).
244		real :: dz_extra ! The thickness of an added layer to append to the woa09_dz profile when
245		! maximum_depth is large [m] (not in Z).
246		real :: nominalDepth ! Depth of ocean bottom in thickness units (positive downward) [H ~> m or kg m-2]
247		real :: depth_q ! A depth scale factor [nondim]
248		real :: depth_s ! The end of the shallow Z regime [m]
249		real :: depth_d ! The start of the deep Z regime [m]
250		real :: adaptTimeRatio, adaptZoomCoeff ! Temporary variables for input parameters [nondim]
251		real :: adaptBuoyCoeff, adaptAlpha ! Temporary variables for input parameters [nondim]
252		real :: adaptZoom ! The thickness of the near-surface zooming region with the adaptive coordinate [H ~> m or kg m-2]
253		real :: adaptDrho0 ! Reference density difference for stratification-dependent diffusion. [R ~> kg m-3]
254		integer :: i, j, k, nzf(4)
255	2	real, dimension(:), allocatable :: dz ! Resolution (thickness) in units of coordinate, which may be [m]
256		! or [Z ~> m] or [H ~> m or kg m-2] or [R ~> kg m-3] or other units.
257	2	real, dimension(:,:), allocatable :: dz_2d ! 2D resolution (thickness) in units of coordinate, which may be [m]
258		! or [Z ~> m] or [H ~> m or kg m-2] or [R ~> kg m-3] or other units.
259	2	real, dimension(:,:,:), allocatable :: dz_3d ! 3D resolution (thickness) in units of coordinate, which may be [m]
260		! or [Z ~> m] or [H ~> m or kg m-2] or [R ~> kg m-3] or other units.
261	2	real, dimension(:), allocatable :: dz_shallow ! Shallow resolution (thickness), for HYBRID_MAP or HYBRID_3D [m]
262	2	real, dimension(:,:), allocatable :: rho_target_2d ! 2D target density used in HYBRID mode [kg m-3]
263	2	real, dimension(:,:,:), allocatable :: rho_target_3d ! 3D target density used in HYBRID mode [kg m-3]
264	2	real, dimension(:,:), allocatable :: index_map ! Region array of indexes for HYBRID_MAP [nondim]
265	2	real, dimension(:), allocatable :: h_max ! Maximum layer thicknesses [H ~> m or kg m-2]
266	2	real, dimension(:), allocatable :: z_max ! Maximum interface depths [H ~> m or kg m-2] or other
267		! units depending on the coordinate
268	2	real, dimension(:), allocatable :: dz_max ! Thicknesses used to find maximum interface depths
269		! [H ~> m or kg m-2] or other units
270	2	real, dimension(:), allocatable :: rho_target ! Target density used in HYBRID mode [kg m-3]
271		! Thicknesses [m] that give level centers approximately corresponding to table 2 of WOA09
272		! These are approximate because the WOA09 depths are not smoothly spaced. Levels
273		! 1, 4, 5, 9, 12, 24, and 36 are 2.5, 2.5, 1.25 12.5, 37.5 and 62.5 m deeper than WOA09
274		! but all others are identical.
275		real, dimension(40) :: woa09_dz_approx = (/ 5., 10., 10., 15., 22.5, 25., 25., 25., &
276		37.5, 50., 50., 75., 100., 100., 100., 100., &
277		100., 100., 100., 100., 100., 100., 100., 175., &
278		250., 375., 500., 500., 500., 500., 500., 500., &
279		500., 500., 500., 500., 500., 500., 500., 500. /)
280		! These are the actual spacings [m] between WOA09 depths which, if used for layer thickness, places
281		! the interfaces at the WOA09 depths.
282		real, dimension(39) :: woa09_dzi = (/ 10., 10., 10., 20., 25., 25., 25., 25., &
283		50., 50., 50., 100., 100., 100., 100., 100., &
284		100., 100., 100., 100., 100., 100., 100., 250., &
285		250., 500., 500., 500., 500., 500., 500., 500., &
286		500., 500., 500., 500., 500., 500., 500. /)
287		! These are the spacings [m] between WOA23 depths from table 3 of
288		! https://www.ncei.noaa.gov/data/oceans/woa/WOA13/DOC/woa13documentation.pdf
289		real, dimension(136) :: woa23_dzi = (/ 5., 5., 5., 5., 5., 5., 5., 5., 5., 5., &
290		5., 5., 5., 5., 5., 5., 5., 5., 5., 5., &
291		25., 25., 25., 25., 25., 25., 25., 25., 25., 25., &
292		25., 25., 25., 25., 25., 25., 50., 50., 50., 50., &
293		50., 50., 50., 50., 50., 50., 50., 50., 50., 50., &
294		50., 50., 50., 50., 50., 50., 50., 50., 50., 50., &
295		50., 50., 50., 50., 50., 50., 100., 100., 100., 100., &
296		100., 100., 100., 100., 100., 100., 100., 100., 100., 100., &
297		100., 100., 100., 100., 100., 100., 100., 100., 100., 100., &
298		100., 100., 100., 100., 100., 100., 100., 100., 100., 100., &
299		100., 100., 100., 100., 100., 100., 100., 100., 100., 100., &
300		100., 100., 100., 100., 100., 100., 100., 100., 100., 100., &
301		100., 100., 100., 100., 100., 100., 100., 100., 100., 100., &
302		100., 100., 100., 100., 100., 100. /)
303
304	2	call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
305	2	inputdir = slasher(inputdir)
306
307	2	main_parameters=.false.
308	2	if (len_trim(param_prefix)==0) main_parameters=.true.
309	2	if (main_parameters .and. len_trim(param_suffix)>0) call MOM_error(FATAL,trim(mdl)//&
310	0	' initialize_regridding: Suffix provided without prefix for parameter names!')
311
312	2	CS%nk = 0
313	2	CS%regridding_scheme = coordinateMode(coord_mode)
314	2	coord_is_state_dependent = state_dependent(coord_mode)
315	2	maximum_depth = US%Z_to_m*max_depth
316
317	2	if (main_parameters) then
318		! Read coordinate units parameter (main model = REGRIDDING_COORDINATE_UNITS)
319		call get_param(param_file, mdl, "REGRIDDING_COORDINATE_UNITS", coord_units, &
320	1	"Units of the regridding coordinate.", default=coordinateUnits(coord_mode))
321		else
322	1	coord_units=coordinateUnits(coord_mode)
323		endif
324
325	2	if (coord_is_state_dependent) then
326	0	if (main_parameters) then
327	0	param_name = "INTERPOLATION_SCHEME"
328	0	string2 = regriddingDefaultInterpScheme
329		else
330	0	param_name = create_coord_param(param_prefix, "INTERP_SCHEME", param_suffix)
331	0	string2 = 'PPM_H4' ! Default for diagnostics
332		endif
333		call get_param(param_file, mdl, param_name, string, &
334		"This sets the interpolation scheme to use to "//&
335		"determine the new grid. These parameters are "//&
336		"only relevant when REGRIDDING_COORDINATE_MODE is "//&
337		"set to a function of state. Otherwise, it is not "//&
338		"used. It can be one of the following schemes: \n"//&
339	0	trim(regriddingInterpSchemeDoc), default=trim(string2))
340	0	call set_regrid_params(CS, interp_scheme=string)
341
342		call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
343		"This sets the default value for the various _ANSWER_DATE parameters.", &
344	0	default=99991231)
345		call get_param(param_file, mdl, "REMAPPING_ANSWER_DATE", remap_answer_date, &
346		"The vintage of the expressions and order of arithmetic to use for remapping. "//&
347		"Values below 20190101 result in the use of older, less accurate expressions "//&
348		"that were in use at the end of 2018. Higher values result in the use of more "//&
349		"robust and accurate forms of mathematically equivalent expressions.", &
350	0	default=default_answer_date, do_not_log=.not.GV%Boussinesq)
351	0	if (.not.GV%Boussinesq) remap_answer_date = max(remap_answer_date, 20230701)
352	0	call set_regrid_params(CS, remap_answer_date=remap_answer_date)
353		call get_param(param_file, mdl, "REGRIDDING_ANSWER_DATE", regrid_answer_date, &
354		"The vintage of the expressions and order of arithmetic to use for regridding. "//&
355		"Values below 20190101 result in the use of older, less accurate expressions "//&
356		"that were in use at the end of 2018. Higher values result in the use of more "//&
357		"robust and accurate forms of mathematically equivalent expressions.", &
358	0	default=default_answer_date, do_not_log=.not.GV%Boussinesq)
359	0	if (.not.GV%Boussinesq) regrid_answer_date = max(regrid_answer_date, 20230701)
360	0	call set_regrid_params(CS, regrid_answer_date=regrid_answer_date)
361		endif
362
363	2	if (main_parameters .and. coord_is_state_dependent) then
364		call get_param(param_file, mdl, "BOUNDARY_EXTRAPOLATION", tmpLogical, &
365		"When defined, a proper high-order reconstruction "//&
366		"scheme is used within boundary cells rather "//&
367		"than PCM. E.g., if PPM is used for remapping, a "//&
368		"PPM reconstruction will also be used within "//&
369	0	"boundary cells.", default=regriddingDefaultBoundaryExtrapolation)
370	0	call set_regrid_params(CS, boundary_extrapolation=tmpLogical)
371		else
372	2	call set_regrid_params(CS, boundary_extrapolation=.false.)
373		endif
374
375		! Read coordinate configuration parameter (main model = ALE_COORDINATE_CONFIG)
376	2	if (main_parameters) then
377	1	param_name = "ALE_COORDINATE_CONFIG"
378	1	coord_res_param = "ALE_RESOLUTION"
379	1	string2 = 'UNIFORM'
380		else
381	1	param_name = create_coord_param(param_prefix, "DEF", param_suffix)
382	1	coord_res_param = create_coord_param(param_prefix, "RES", param_suffix)
383	1	string2 = 'UNIFORM'
384	1	if ((maximum_depth>3000.) .and. (maximum_depth<9250.)) string2='WOA09' ! For convenience
385		endif
386		call get_param(param_file, mdl, param_name, string, &
387		"Determines how to specify the coordinate "//&
388		"resolution. Valid options are:\n"//&
389		" PARAM - use the vector-parameter "//trim(coord_res_param)//"\n"//&
390		" UNIFORM[:N] - uniformly distributed\n"//&
391		" FILE:string - read from a file. The string specifies\n"//&
392		" the filename and variable name, separated\n"//&
393		" by a comma or space, e.g. FILE:lev.nc,dz\n"//&
394		" or FILE:lev.nc,interfaces=zw\n"//&
395		" WOA09[:N] - the WOA09 vertical grid (approximately)\n"//&
396		" WOA09INT[:N] - layers spanned by the WOA09 depths\n"//&
397		" WOA23INT[:N] - layers spanned by the WOA23 depths\n"//&
398		" FNC1:string - FNC1:dz_min,H_total,power,precision\n"//&
399		" HYBRID:string - read from a file. The string specifies\n"//&
400		" the filename and two variable names, separated\n"//&
401		" by a comma or space, for sigma-2 and dz.\n"//&
402		" e.g. HYBRID:vgrid.nc,sigma2,dz\n"//&
403		" HYBRID_3D:string - read from a file. The string specifies\n"//&
404		" the filename and two 3D variable names, separated\n"//&
405		" by a comma or space, for sigma-2 and dz. The\n"//&
406		" latter can be FNC1:string which is used everywhere.\n"//&
407		" e.g. HYBRID_3D:vgrid.nc,sigma2,dz\n"//&
408		" HYBRID_MAP:string - read from a file. The string specifies\n"//&
409		" the filename and three variable names, separated\n"//&
410		" by a comma or space, for map, sigma-2 and dz.\n"//&
411		" Map is a spatial index array with, maxval(map)=N,\n"//&
412		" and the others are 2D arrays containing N profiles.\n"//&
413		" Map typically contains integer values, but it can\n"//&
414		" contain real values, I+w, which imply using\n"//&
415		" the weighted sum of profiles I and I+1.\n"//&
416		" Dz can be FNC1:string which is used everywhere.\n"//&
417		" e.g. HYBRID_MAP:vgrid.nc,map,sigma2,dz",&
418	2	default=trim(string2))
419		message = "The distribution of vertical resolution for the target\n"//&
420		"grid used for Eulerian-like coordinates. For example,\n"//&
421		"in z-coordinate mode, the parameter is a list of level\n"//&
422		"thicknesses (in m). In sigma-coordinate mode, the list\n"//&
423	2	"is of non-dimensional fractions of the water column."
424	2	if (index(trim(string),'UNIFORM')==1) then
425	0	if (len_trim(string)==7) then
426	0	ke = GV%ke ! Use model nk by default
427	0	tmpReal = maximum_depth
428	0	elseif (index(trim(string),'UNIFORM:')==1 .and. len_trim(string)>8) then
429		! Format is "UNIFORM:N" or "UNIFORM:N,dz"
430	0	ke = extract_integer(string(9:len_trim(string)),'',1)
431	0	tmpReal = extract_real(string(9:len_trim(string)),',',2,missing_value=maximum_depth)
432		else
433		call MOM_error(FATAL,trim(mdl)//', initialize_regridding: '// &
434	0	'Unable to interpret "'//trim(string)//'".')
435		endif
436	0	allocate(dz(ke))
437		dz(:) = uniformResolution(ke, coord_mode, tmpReal, &
438		US%R_to_kg_m3(GV%Rlay(1) + 0.5(GV%Rlay(1)-GV%Rlay(min(2,ke)))), &
439	0	US%R_to_kg_m3(GV%Rlay(ke) + 0.5(GV%Rlay(ke)-GV%Rlay(max(ke-1,1)))) )
440	0	if (main_parameters) call log_param(param_file, mdl, "!"//coord_res_param, dz, &
441	0	trim(message), units=trim(coord_units))
442	2	elseif (trim(string)=='PARAM') then
443		! Read coordinate resolution (main model = ALE_RESOLUTION)
444	0	allocate(dz(1001))
445	0	dz(:) = -1. ! Setting to <0 allows detection of unset elements
446	0	call get_param(param_file, mdl, coord_res_param, dz, "Scan", units="", do_not_log=.true.)
447	0	if (dz(1001)>=0.) call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// &
448	0	"PARAM specification is limited to 1000 values. Hack the code to use more!")
449	0	do ke=1,1000 ! Find number of defined levels
450	0	if (dz(ke+1)<0.) exit
451		enddo
452	0	deallocate(dz)
453	0	allocate(dz(ke)) ! Allocate with the correct number of levels, and re-read thicknesses
454		call get_param(param_file, mdl, coord_res_param, dz, &
455	0	trim(message), units=trim(coord_units), fail_if_missing=.true.)
456	2	elseif (index(trim(string),'FILE:')==1) then
457		! FILE:filename,var_name is assumed to be reading level thickness variables
458		! FILE:filename,interfaces=var_name reads positions
459	0	if (string(6:6)=='.' .or. string(6:6)=='/') then
460		! If we specified "FILE:./xyz" or "FILE:/xyz" then we have a relative or absolute path
461	0	fileName = trim( extractWord(trim(string(6:80)), 1) )
462		else
463		! Otherwise assume we should look for the file in INPUTDIR
464	0	fileName = trim(inputdir) // trim( extractWord(trim(string(6:80)), 1) )
465		endif
466	0	if (.not. file_exists(fileName)) call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// &
467	0	"Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(string)//")")
468
469	0	varName = trim( extractWord(trim(string(6:)), 2) )
470	0	if (len_trim(varName)==0) then
471	0	if (field_exists(fileName,'dz')) then ; varName = 'dz'
472	0	elseif (field_exists(fileName,'dsigma')) then ; varName = 'dsigma'
473	0	elseif (field_exists(fileName,'ztest')) then ; varName = 'ztest'
474		else ; call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// &
475	0	"Coordinate variable not specified and none could be guessed.")
476		endif
477		endif
478		! This check fails when the variable is a dimension variable! -AJA
479		!if (.not. field_exists(fileName,trim(varName))) call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// &
480		! "Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")")
481	0	if (CS%regridding_scheme == REGRIDDING_SIGMA) then
482	0	expected_units = 'nondim' ; alt_units = expected_units
483	0	elseif (CS%regridding_scheme == REGRIDDING_RHO) then
484	0	expected_units = 'kg m-3' ; alt_units = expected_units
485		else
486	0	expected_units = 'meters' ; alt_units = 'm'
487		endif
488	0	if (index(trim(varName),'interfaces=')==1) then
489	0	varName=trim(varName(12:))
490	0	call verify_variable_units(filename, varName, expected_units, message, ierr, alt_units)
491	0	if (ierr) call MOM_error(FATAL, trim(mdl)//", initialize_regridding: "//&
492		"Unsupported format in grid definition '"//trim(filename)//&
493	0	"'. Error message "//trim(message))
494	0	call field_size(trim(fileName), trim(varName), nzf)
495	0	ke = nzf(1)-1
496	0	if (ke < 1) call MOM_error(FATAL, trim(mdl)//" initialize_regridding via Var "//&
497		trim(varName)//"in FILE "//trim(filename)//&
498	0	" requires at least 2 target interface values.")
499	0	if (CS%regridding_scheme == REGRIDDING_RHO) then
500	0	allocate(rho_target(ke+1))
501	0	call MOM_read_data(trim(fileName), trim(varName), rho_target)
502		else
503	0	allocate(dz(ke))
504	0	allocate(z_max(ke+1))
505	0	call MOM_read_data(trim(fileName), trim(varName), z_max)
506	0	dz(:) = abs(z_max(1:ke) - z_max(2:ke+1))
507	0	deallocate(z_max)
508		endif
509		else
510		! Assume reading resolution
511	0	call field_size(trim(fileName), trim(varName), nzf)
512	0	ke = nzf(1)
513	0	allocate(dz(ke))
514	0	call MOM_read_data(trim(fileName), trim(varName), dz)
515		endif
516	0	if (main_parameters .and. (ke/=GV%ke)) then
517		call MOM_error(FATAL,trim(mdl)//', initialize_regridding: '// &
518	0	'Mismatch in number of model levels and "'//trim(string)//'".')
519		endif
520	0	if (main_parameters) call log_param(param_file, mdl, "!"//coord_res_param, dz, &
521	0	trim(message), units=coordinateUnits(coord_mode))
522	2	elseif (index(trim(string),'FNC1:')==1) then
523	1	ke = GV%ke ; allocate(dz(ke))
524	1	call dz_function1( trim(string(6:)), dz )
525	1	if (main_parameters) call log_param(param_file, mdl, "!"//coord_res_param, dz, &
526	1	trim(message), units=coordinateUnits(coord_mode))
527	1	elseif (index(trim(string),'RFNC1:')==1) then
528		! Function used for set target interface densities
529	0	ke = rho_function1( trim(string(7:)), rho_target )
530	1	elseif (index(trim(string),'HYBRID:')==1) then
531	0	ke = GV%ke
532	0	allocate(dz(ke))
533	0	allocate(rho_target(ke+1))
534		! The following assumes the FILE: syntax of above but without "FILE:" in the string
535	0	varName = trim( extractWord(trim(string(8:)), 3) )
536	0	if (varname == " ") call MOM_error(FATAL, &
537		trim(mdl)//", initialize_regridding: HYBRID "// &
538	0	"Too few arguments in ("//trim(string)//")")
539	0	fileName = trim( extractWord(trim(string(8:)), 1) )
540	0	if (fileName(1:1)/='.' .and. filename(1:1)/='/') fileName = trim(inputdir) // trim( fileName )
541	0	if (.not. file_exists(fileName)) call MOM_error(FATAL, &
542		trim(mdl)//", initialize_regridding: HYBRID "// &
543	0	"Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(string)//")")
544	0	varName = trim( extractWord(trim(string(8:)), 2) )
545	0	if (.not. field_exists(fileName,varName)) call MOM_error(FATAL, &
546		trim(mdl)//", initialize_regridding: HYBRID "// &
547	0	"Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")")
548	0	call MOM_read_data(trim(fileName), trim(varName), rho_target)
549	0	varName = trim( extractWord(trim(string(8:)), 3) )
550	0	if (varName(1:5) == 'FNC1:') then ! Use FNC1 to calculate dz
551	0	call dz_function1( trim(string((index(trim(string),'FNC1:')+5):)), dz )
552		else ! Read dz from file
553	0	if (.not. field_exists(fileName,varName)) call MOM_error(FATAL, &
554		trim(mdl)//", initialize_regridding: HYBRID "// &
555	0	"Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")")
556	0	call MOM_read_data(trim(fileName), trim(varName), dz)
557		endif
558	0	if (main_parameters) then
559		call log_param(param_file, mdl, "!"//coord_res_param, dz, &
560	0	trim(message), units=coordinateUnits(coord_mode))
561		call log_param(param_file, mdl, "!TARGET_DENSITIES", rho_target, &
562	0	'HYBRID target densities for interfaces', units="kg m-3")
563		endif
564	1	elseif (index(trim(string),'HYBRID_3D:')==1) then
565	0	ke = GV%ke
566	0	allocate(dz_3d(SZI_(G),SZJ_(G),ke), source=0.0)
567	0	allocate(rho_target_3d(SZI_(G),SZJ_(G),ke+1), source=0.0)
568		! The following assumes the FILE: syntax of above but without "FILE:" in the string
569	0	varName = trim( extractWord(trim(string(11:)), 3) )
570	0	if (varname == " ") call MOM_error(FATAL, &
571		trim(mdl)//", initialize_regridding: HYBRID_3D "// &
572	0	"Too few arguments in ("//trim(string)//")")
573	0	fileName = trim( extractWord(trim(string(11:)), 1) )
574	0	if (fileName(1:1)/='.' .and. filename(1:1)/='/') fileName = trim(inputdir) // trim( fileName )
575	0	if (.not. file_exists(fileName)) call MOM_error(FATAL, &
576		trim(mdl)//", initialize_regridding: HYBRID_3D "// &
577	0	"Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(string)//")")
578	0	varName = trim( extractWord(trim(string(11:)), 2) )
579	0	if (.not. field_exists(fileName,varName)) call MOM_error(FATAL, &
580		trim(mdl)//", initialize_regridding: HYBRID_3D "// &
581	0	"Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")")
582	0	call MOM_read_data(trim(fileName), trim(varName), rho_target_3d, G%Domain)
583	0	call pass_var(rho_target_3d, G%Domain, halo=1)
584	0	varName = trim( extractWord(trim(string(11:)), 3) )
585	0	if (varName(1:5) == 'FNC1:') then ! Use FNC1 to calculate dz_3d
586	0	allocate(dz(ke))
587	0	call dz_function1( trim(string((index(trim(string),'FNC1:')+5):)), dz )
588		! Adjust target grid to be consistent with maximum_depth
589	0	tmpReal = sum( dz(:) )
590	0	if (tmpReal < maximum_depth) then
591	0	dz(ke) = dz(ke) + ( maximum_depth - tmpReal )
592		endif
593	0	do i=G%isc-1,G%iec+1 ; do j=G%jsc-1,G%jec+1
594	0	if (G%mask2dT(i,j)>0.) then
595	0	do k=1,ke
596	0	dz_3d(i,j,k) = dz(k)
597		enddo
598		endif !mask2dT
599		enddo ; enddo
600	0	if (main_parameters) then
601		call log_param(param_file, mdl, "!"//coord_res_param, dz, &
602	0	trim(message), units=coordinateUnits(coord_mode))
603		endif
604		else ! Read dz from file
605	0	if (.not. field_exists(fileName,varName)) call MOM_error(FATAL, &
606		trim(mdl)//", initialize_regridding: HYBRID_3D "// &
607	0	"Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")")
608	0	call MOM_read_data(trim(fileName), trim(varName), dz_3d, G%Domain)
609	0	call pass_var(dz_3d, G%Domain, halo=1)
610		! set nominal 1-d dz to UNIFORM
611	0	allocate(dz(ke))
612		dz(:) = uniformResolution(ke, coord_mode, maximum_depth, &
613		US%R_to_kg_m3(GV%Rlay(1) + 0.5(GV%Rlay(1)-GV%Rlay(min(2,ke)))), &
614	0	US%R_to_kg_m3(GV%Rlay(ke) + 0.5(GV%Rlay(ke)-GV%Rlay(max(ke-1,1)))) )
615		endif !dz
616	1	elseif (index(trim(string),'HYBRID_MAP:')==1) then
617	0	ke = GV%ke
618	0	allocate(dz_3d(SZI_(G),SZJ_(G),ke), source=0.0)
619	0	allocate(rho_target_3d(SZI_(G),SZJ_(G),ke+1), source=0.0)
620	0	allocate(index_map(SZI_(G),SZJ_(G)), source=1.0)
621		! The following assumes the FILE: syntax of above but without "FILE:" in the string
622	0	varName = trim( extractWord(trim(string(12:)), 4) )
623	0	if (varname == " ") call MOM_error(FATAL, &
624		trim(mdl)//", initialize_regridding: HYBRID_3D "// &
625	0	"Too few arguments in ("//trim(string)//")")
626	0	fileName = trim( extractWord(trim(string(12:)), 1) )
627	0	if (fileName(1:1)/='.' .and. filename(1:1)/='/') fileName = trim(inputdir) // trim( fileName )
628	0	if (.not. file_exists(fileName)) call MOM_error(FATAL, &
629		trim(mdl)//", initialize_regridding: HYBRID_MAP "// &
630	0	"Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(string)//")")
631	0	varName = trim( extractWord(trim(string(12:)), 2) )
632	0	if (.not. field_exists(fileName,varName)) call MOM_error(FATAL, &
633		trim(mdl)//", initialize_regridding: HYBRID_MAP "// &
634	0	"Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")")
635	0	call MOM_read_data(trim(fileName), trim(varName), index_map, G%Domain)
636	0	call pass_var(index_map, G%Domain, halo=1)
637		!find maximum index
638	0	np = 1
639	0	do j=G%jsc, G%jec ; do i=G%isc, G%iec
640	0	np = max(np,ceiling(index_map(i,j)))
641		enddo ; enddo
642	0	call max_across_PEs(np)
643	0	write(string2,"(i3)") np
644		call MOM_error(NOTE, &
645	0	trim(mdl)//", initialize_regridding: HYBRID_MAP NP="//trim(string2))
646	0	if (np<1) call MOM_error(FATAL, &
647	0	trim(mdl)//", initialize_regridding: HYBRID_MAP to small NP from "//trim(varName))
648	0	allocate(dz_2d(ke,np))
649	0	allocate(rho_target_2d(ke+1,np))
650	0	varName = trim( extractWord(trim(string(12:)), 3) )
651	0	if (.not. field_exists(fileName,varName)) call MOM_error(FATAL, &
652		trim(mdl)//", initialize_regridding: HYBRID_MAP "// &
653	0	"Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")")
654		! MOM_read_data can't handle this array
655	0	call read_variable(trim(fileName), trim(varName), rho_target_2d)
656	0	if (main_parameters) then
657		call log_param(param_file, mdl, "!TARGET_DENSITIES", rho_target_2d(:,1), &
658	0	'HYBRID target densities for interfaces', units="kg m-3")
659		endif
660	0	do i=G%isc-1,G%iec+1 ; do j=G%jsc-1,G%jec+1
661	0	if (G%mask2dT(i,j)>0.) then
662	0	nfloor = floor(index_map(i,j))
663	0	nceiling = ceiling(index_map(i,j))
664	0	if (nfloor<1 .or. nceiling>np) then
665	0	write(0,'(a,2i5,a,g20.6)') 'HYBRID_MAP: i,j=',i,j,'index_map(i,j)=', index_map(i,j)
666		call MOM_error(FATAL, trim(mdl)//", initialize_regridding: HYBRID_MAP "// &
667	0	"index_map out of range")
668		endif
669	0	if (nfloor == nceiling) then
670	0	do k=1,ke+1
671	0	rho_target_3d(i,j,k) = rho_target_2d(k,nfloor)
672		enddo
673		else
674	0	nfrac = index_map(i,j) - nfloor !between 0.0 and 1.0
675	0	do k=1,ke+1
676		rho_target_3d(i,j,k) = (1.0-nfrac)*rho_target_2d(k,nfloor) + &
677	0	nfrac *rho_target_2d(k,nceiling)
678		enddo
679		endif !integer:else
680		endif !mask2dT
681		enddo ; enddo
682	0	varName = trim( extractWord(trim(string(12:)), 4) )
683	0	if (varName(1:5) == 'FNC1:') then ! Use FNC1 to calculate dz_3d
684	0	allocate(dz(ke))
685	0	call dz_function1( trim(string((index(trim(string),'FNC1:')+5):)), dz )
686		! Adjust target grid to be consistent with maximum_depth
687	0	tmpReal = sum( dz(:) )
688	0	if (tmpReal < maximum_depth) then
689	0	dz(ke) = dz(ke) + ( maximum_depth - tmpReal )
690		endif
691	0	do i=G%isc-1,G%iec+1 ; do j=G%jsc-1,G%jec+1
692	0	if (G%mask2dT(i,j)>0.) then
693	0	do k=1,ke
694	0	dz_3d(i,j,k) = dz(k)
695		enddo
696		endif !mask2dT
697		enddo ; enddo
698	0	if (main_parameters) then
699		call log_param(param_file, mdl, "!"//coord_res_param, dz, &
700	0	trim(message), units=coordinateUnits(coord_mode))
701		endif
702		else ! Read dz from file
703	0	if (.not. field_exists(fileName,varName)) call MOM_error(FATAL, &
704		trim(mdl)//", initialize_regridding: HYBRID_MAP "// &
705	0	"Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")")
706		! MOM_read_data can't handle this array
707	0	call read_variable(trim(fileName), trim(varName), dz_2d)
708	0	if (main_parameters) then
709		call log_param(param_file, mdl, "!"//coord_res_param, dz_2d(:,1), &
710	0	trim(message), units=coordinateUnits(coord_mode))
711		endif
712	0	do i=1,np
713	0	tmpReal = sum( dz_2d(:,i) )
714	0	if (tmpReal < maximum_depth) then
715	0	dz_2d(ke,i) = dz_2d(ke,i) + ( maximum_depth - tmpReal )
716		endif
717		enddo
718	0	allocate(dz(ke))
719	0	dz(:) = dz_2d(:,1)
720	0	if (main_parameters) then
721		call log_param(param_file, mdl, "!"//coord_res_param, dz, &
722	0	trim(message), units=coordinateUnits(coord_mode))
723		endif
724	0	do i=G%isc-1,G%iec+1 ; do j=G%jsc-1,G%jec+1
725	0	if (G%mask2dT(i,j)>0.) then
726	0	nfloor = floor(index_map(i,j))
727	0	nceiling = ceiling(index_map(i,j))
728	0	if (nfloor == nceiling) then
729	0	do k=1,ke
730	0	dz_3d(i,j,k) = dz_2d(k,nfloor)
731		enddo
732		else
733	0	nfrac = index_map(i,j) - nfloor !between 0.0 and 1.0
734	0	do k=1,ke
735		dz_3d(i,j,k) = (1.0-nfrac)*dz_2d(k,nfloor) + &
736	0	nfrac *dz_2d(k,nceiling)
737		enddo
738		endif !integer:else
739		endif !mask2dT
740		enddo ; enddo
741		endif !dz
742	0	deallocate(index_map)
743	0	deallocate(rho_target_2d)
744	0	deallocate(dz_2d)
745	1	elseif (index(trim(string),'WOA09INT')==1) then
746	0	if (len_trim(string)==8) then ! string=='WOA09INT'
747	0	tmpReal = 0. ; ke = 0 ; dz_extra = 0.
748	0	do while (tmpReal<maximum_depth)
749	0	ke = ke + 1
750	0	if (ke > size(woa09_dzi)) then
751	0	dz_extra = maximum_depth - tmpReal
752	0	exit
753		endif
754	0	tmpReal = tmpReal + woa09_dzi(ke)
755		enddo
756	0	elseif (index(trim(string),'WOA09INT:')==1) then ! string starts with 'WOA09INT:'
757	0	if (len_trim(string)==9) call MOM_error(FATAL,trim(mdl)//', initialize_regridding: '// &
758	0	'Expected string of form "WOA09INT:N" but got "'//trim(string)//'".')
759	0	ke = extract_integer(string(10:len_trim(string)),'',1)
760	0	if (ke>39 .or. ke<1) call MOM_error(FATAL,trim(mdl)//', initialize_regridding: '// &
761	0	'For "WOA05INT:N" N must 0<N<40 but got "'//trim(string)//'".')
762		endif
763	0	allocate(dz(ke))
764	0	do k=1,min(ke, size(woa09_dzi))
765	0	dz(k) = woa09_dzi(k)
766		enddo
767	0	if (ke > size(woa09_dzi)) dz(ke) = dz_extra
768	1	elseif (index(trim(string),'WOA23INT')==1) then
769	0	if (len_trim(string)==8) then ! string=='WOA23INT'
770	0	tmpReal = 0. ; ke = 0 ; dz_extra = 0.
771	0	do while (tmpReal<maximum_depth)
772	0	ke = ke + 1
773	0	if (ke > size(woa23_dzi)) then
774	0	dz_extra = maximum_depth - tmpReal
775	0	exit
776		endif
777	0	tmpReal = tmpReal + woa23_dzi(ke)
778		enddo
779	0	elseif (index(trim(string),'WOA23INT:')==1) then ! string starts with 'WOA23INT:'
780	0	if (len_trim(string)==9) call MOM_error(FATAL,trim(mdl)//', initialize_regridding: '// &
781	0	'Expected string of form "WOA23INT:N" but got "'//trim(string)//'".')
782	0	ke = extract_integer(string(10:len_trim(string)),'',1)
783	0	if (ke>39 .or. ke<1) call MOM_error(FATAL,trim(mdl)//', initialize_regridding: '// &
784	0	'For "WOA05INT:N" N must 0<N<40 but got "'//trim(string)//'".')
785		endif
786	0	allocate(dz(ke))
787	0	do k=1,min(ke, size(woa23_dzi))
788	0	dz(k) = woa23_dzi(k)
789		enddo
790	0	if (ke > size(woa23_dzi)) dz(ke) = dz_extra
791	1	elseif (index(trim(string),'WOA09')==1) then
792	1	if (len_trim(string)==5) then ! string=='WOA09'
793	1	tmpReal = 0. ; ke = 0 ; dz_extra = 0.
794	34	do while (tmpReal<maximum_depth)
795	33	ke = ke + 1
796	33	if (ke > size(woa09_dz_approx)) then
797	0	dz_extra = maximum_depth - tmpReal
798	0	exit
799		endif
800	33	tmpReal = tmpReal + woa09_dz_approx(ke)
801		enddo
802	0	elseif (index(trim(string),'WOA09:')==1) then ! string starts with 'WOA09:'
803	0	if (len_trim(string)==6) call MOM_error(FATAL,trim(mdl)//', initialize_regridding: '// &
804	0	'Expected string of form "WOA09:N" but got "'//trim(string)//'".')
805	0	ke = extract_integer(string(7:len_trim(string)),'',1)
806	0	if (ke>40 .or. ke<1) call MOM_error(FATAL,trim(mdl)//', initialize_regridding: '// &
807	0	'For "WOA05:N" N must 0<N<41 but got "'//trim(string)//'".')
808		endif
809	1	allocate(dz(ke))
810	34	do k=1,min(ke, size(woa09_dz_approx))
811	34	dz(k) = woa09_dz_approx(k)
812		enddo
813	1	if (ke > size(woa09_dz_approx)) dz(ke) = dz_extra
814		else
815		call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// &
816	0	"Unrecognized coordinate configuration"//trim(string))
817		endif
818
819	2	if (main_parameters) then
820		! This is a work around to apparently needed to work with the from_Z initialization... ???
821		if (coordinateMode(coord_mode) == REGRIDDING_ZSTAR .or. &
822		coordinateMode(coord_mode) == REGRIDDING_HYCOM1 .or. &
823	1	coordinateMode(coord_mode) == REGRIDDING_HYBGEN .or. &
824		coordinateMode(coord_mode) == REGRIDDING_ADAPTIVE) then
825	1	if (allocated(dz)) then
826		! Adjust target grid to be consistent with maximum_depth
827	76	tmpReal = sum( dz(:) )
828	1	if (tmpReal < maximum_depth) then
829	0	dz(ke) = dz(ke) + ( maximum_depth - tmpReal )
830	1	elseif (tmpReal > maximum_depth) then
831	0	if ( dz(ke) + ( maximum_depth - tmpReal ) > 0. ) then
832	0	dz(ke) = dz(ke) + ( maximum_depth - tmpReal )
833		else
834		call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// &
835	0	"MAXIMUM_DEPTH was too shallow to adjust bottom layer of DZ!"//trim(string))
836		endif
837		endif
838		endif !allocated(dz)
839		endif
840		endif
841
842	2	if (coordinateMode(coord_mode) == REGRIDDING_HYCOM1) then
843	0	allocate(dz_shallow(ke))
844		call get_param(param_file, mdl, "SHALLOW_"//trim(coord_res_param), dz_shallow, &
845		"HYBGEN-style Z-sigma-Z near surface fixed coordinate. "//&
846		"The default of all zeros turns this option off. "//&
847		"Let N_SIGMA be the number of consecutive non-zero entries, typically < NK. "//&
848		"Use SHALLOW_"//trim(coord_res_param)//" when rest depth is shallower than "//&
849		"SUM(SHALLOW_"//trim(coord_res_param)//"(1:N_SIGMA)). "//&
850		"Use "//trim(coord_res_param)//" when rest depth is deeper than "//&
851		"SUM("//trim(coord_res_param)//"(1:N_SIGMA)). "//&
852		"Otherwise use a linear sum of the two weighted by rest depth.",&
853	0	units="m", default=0.0)
854	0	n_sigma = ke
855	0	depth_s = 0.0
856	0	do k= 1,ke
857	0	depth_s = depth_s + dz_shallow(k)
858	0	if (dz_shallow(k) == 0.0) then
859	0	n_sigma = k-1
860	0	exit
861		endif
862		enddo
863	0	if (n_sigma > 0) then
864	0	if (main_parameters) call log_param(param_file, mdl, "!N_SIGMA", n_sigma, &
865		"Number of consecutive non-zero entries in SHALLOW_"//&
866	0	trim(coord_res_param)//".")
867	0	if (.not.allocated(dz_3d)) then
868	0	allocate(dz_3d(SZI_(G),SZJ_(G),ke), source=0.0)
869	0	allocate(rho_target_3d(SZI_(G),SZJ_(G),ke+1), source=0.0)
870	0	do i=G%isc-1,G%iec+1 ; do j=G%jsc-1,G%jec+1
871	0	if (G%mask2dT(i,j)>0.) then
872	0	do k=1,ke
873	0	dz_3d(i,j,k) = dz(k)
874		enddo
875	0	do k=1,ke+1
876	0	rho_target_3d(i,j,k) = rho_target(k)
877		enddo
878		endif !mask2dT
879		enddo ; enddo
880		endif
881	0	do i=G%isc-1,G%iec+1 ; do j=G%jsc-1,G%jec+1
882	0	if (G%mask2dT(i,j)>0.) then
883	0	nominalDepth = max(G%meanSL(i,j) + G%bathyT(i,j), 0.0) * US%Z_to_m
884	0	if (nominalDepth <= depth_s) then
885	0	do k= 1,n_sigma
886	0	dz_3d(i,j,k) = dz_shallow(k)
887		enddo
888	0	do k= n_sigma+1,ke
889	0	dz_3d(i,j,k) = dz_shallow(n_sigma)
890		enddo
891		else ! >depth_s
892	0	depth_d = 0.0
893	0	do k= 1,n_sigma
894	0	depth_d = depth_d + dz_3d(i,j,k)
895		enddo
896		! do nothing if nominalDepth >= depth_d
897	0	if (nominalDepth < depth_d) then
898	0	depth_q = (nominalDepth - depth_s) / (depth_d - depth_s)
899	0	do k= 1,n_sigma
900	0	dz_3d(i,j,k) = (1.0-depth_q)dz_shallow(k) + depth_qdz_3d(i,j,k)
901		enddo
902	0	do k= n_sigma+1,ke
903	0	dz_3d(i,j,k) = (1.0-depth_q)dz_shallow(n_sigma) + depth_qdz_3d(i,j,k)
904		enddo
905		endif !<depth_d and >depth_s
906		endif !nominalDepth
907		endif !mask2dT
908		enddo ; enddo
909		endif !n_sigma
910	0	deallocate(dz_shallow)
911		endif !REGRIDDING_HYCOM1
912
913	2	CS%nk=ke
914
915		! Target resolution (for fixed coordinates)
916	2	if (allocated(dz_3d)) then
917	0	allocate( CS%coordinateResolution(CS%nk), source=-1.E30 )
918	0	allocate( CS%coordinateResolution_3d(SZI_(G),SZJ_(G),CS%nk), source=-1.E30 )
919	0	allocate( CS%target_density_3d(SZI_(G),SZJ_(G),CS%nk+1), source=-1.E30*US%kg_m3_to_R )
920		else
921	110	allocate( CS%coordinateResolution(CS%nk), source=-1.E30 )
922	2	if (state_dependent(CS%regridding_scheme)) then
923		! Target values
924	0	allocate( CS%target_density(CS%nk+1), source=-1.E30*US%kg_m3_to_R )
925		endif
926		endif
927
928	2	if (allocated(dz_3d)) then
929		! set both 1d and 3d fields
930	0	call setCoordinateResolution(dz, CS, scale=US%m_to_Z)
931	0	call setCoordinateResolution_3d(dz_3d, CS, scale=US%m_to_Z)
932	0	CS%coord_scale = US%Z_to_m
933	0	deallocate(dz_3d)
934	2	elseif (allocated(dz)) then
935	2	if (coordinateMode(coord_mode) == REGRIDDING_SIGMA) then
936	0	call setCoordinateResolution(dz, CS, scale=1.0)
937	2	elseif (coordinateMode(coord_mode) == REGRIDDING_RHO) then
938	0	call setCoordinateResolution(dz, CS, scale=US%kg_m3_to_R)
939	2	elseif (coordinateMode(coord_mode) == REGRIDDING_ADAPTIVE) then
940	0	call setCoordinateResolution(dz, CS, scale=GV%m_to_H)
941	0	CS%coord_scale = GV%H_to_m
942		else
943	2	call setCoordinateResolution(dz, CS, scale=US%m_to_Z)
944	2	CS%coord_scale = US%Z_to_m
945		endif
946		endif
947
948		! set coord_scale for RHO regridding independent of allocation status of dz
949	2	if (coordinateMode(coord_mode) == REGRIDDING_RHO) then
950	0	CS%coord_scale = US%R_to_kg_m3
951		endif
952
953		! ensure CS%ref_pressure is rescaled properly
954	2	CS%ref_pressure = US%Pa_to_RL2_T2 * CS%ref_pressure
955
956	2	if (allocated(rho_target_3d)) then
957	0	call set_target_densities_3d(CS, G, US%kg_m3_to_R, rho_target_3d)
958	0	deallocate(rho_target_3d)
959	2	elseif (allocated(rho_target)) then
960	0	call set_target_densities(CS, US%kg_m3_to_R*rho_target)
961	0	deallocate(rho_target)
962	2	elseif (coordinateMode(coord_mode) == REGRIDDING_RHO) then
963	0	call set_target_densities_from_GV(GV, US, CS)
964		call log_param(param_file, mdl, "!TARGET_DENSITIES", US%R_to_kg_m3*CS%target_density(:), &
965	0	'RHO target densities for interfaces', "kg m-3")
966		endif
967
968		! initialise coordinate-specific control structure
969	2	call initCoord(CS, G, GV, US, coord_mode, param_file)
970
971	2	if (coord_is_state_dependent) then
972	0	if (main_parameters) then
973		call get_param(param_file, mdl, create_coord_param(param_prefix, "P_REF", param_suffix), &
974		P_Ref, &
975		"The pressure that is used for calculating the coordinate "//&
976		"density. (1 Pa = 1e4 dbar, so 2e7 is commonly used.) "//&
977		"This is only used if USE_EOS and ENABLE_THERMODYNAMICS are true.", &
978	0	units="Pa", default=2.0e7, scale=US%Pa_to_RL2_T2)
979		else
980		call get_param(param_file, mdl, create_coord_param(param_prefix, "P_REF", param_suffix), &
981		P_Ref, &
982		"The pressure that is used for calculating the diagnostic coordinate "//&
983		"density. (1 Pa = 1e4 dbar, so 2e7 is commonly used.) "//&
984		"This is only used for the RHO coordinate.", &
985	0	units="Pa", default=2.0e7, scale=US%Pa_to_RL2_T2)
986		endif
987		call get_param(param_file, mdl, create_coord_param(param_prefix, &
988		"REGRID_COMPRESSIBILITY_FRACTION", param_suffix), tmpReal, &
989		"When interpolating potential density profiles we can add "//&
990		"some artificial compressibility solely to make homogeneous "//&
991	0	"regions appear stratified.", units="nondim", default=0.)
992	0	call set_regrid_params(CS, compress_fraction=tmpReal, ref_pressure=P_Ref)
993		endif
994
995	2	if (main_parameters) then
996		call get_param(param_file, mdl, "MIN_THICKNESS", tmpReal, &
997		"When regridding, this is the minimum layer "//&
998		"thickness allowed.", units="m", scale=GV%m_to_H, &
999	1	default=regriddingDefaultMinThickness )
1000	1	call set_regrid_params(CS, min_thickness=tmpReal)
1001		call get_param(param_file, mdl, "USE_ADJUST_INTERFACE_MOTION", tmpLogical, &
1002		"When regridding, after the primary grid generation, call a function that ensures "//&
1003	1	"positive layer thicknesses. Historically, this was required.", default=.true.)
1004	1	call set_regrid_params(CS, use_adjust_interface_motion=tmpLogical)
1005		else
1006	1	call set_regrid_params(CS, min_thickness=0.)
1007	1	call set_regrid_params(CS, use_adjust_interface_motion=.true.)
1008	1	call set_regrid_params(CS, use_depth_based_time_filter=.true.)
1009		endif
1010
1011	2	if (main_parameters .and. coordinateMode(coord_mode) == REGRIDDING_HYCOM1) then
1012		call get_param(param_file, mdl, "HYCOM1_ONLY_IMPROVES", tmpLogical, &
1013		"When regridding, an interface is only moved if this improves "//&
1014	0	"the fit to the target density.", default=.false.)
1015	0	call set_hycom_params(CS%hycom_CS, only_improves=tmpLogical)
1016		endif
1017
1018	2	CS%use_hybgen_unmix = .false.
1019	2	if (coordinateMode(coord_mode) == REGRIDDING_HYBGEN) then
1020		call get_param(param_file, mdl, "USE_HYBGEN_UNMIX", CS%use_hybgen_unmix, &
1021		"If true, use hybgen unmixing code before regridding.", &
1022	0	default=.false.)
1023		endif
1024
1025	2	if (coordinateMode(coord_mode) == REGRIDDING_ADAPTIVE) then
1026		call get_param(param_file, mdl, "ADAPT_TIME_RATIO", adaptTimeRatio, &
1027	0	"Ratio of ALE timestep to grid timescale.", units="nondim", default=1.0e-1)
1028		call get_param(param_file, mdl, "ADAPT_ZOOM_DEPTH", adaptZoom, &
1029	0	"Depth of near-surface zooming region.", units="m", default=200.0, scale=GV%m_to_H)
1030		call get_param(param_file, mdl, "ADAPT_ZOOM_COEFF", adaptZoomCoeff, &
1031	0	"Coefficient of near-surface zooming diffusivity.", units="nondim", default=0.2)
1032		call get_param(param_file, mdl, "ADAPT_BUOY_COEFF", adaptBuoyCoeff, &
1033	0	"Coefficient of buoyancy diffusivity.", units="nondim", default=0.8)
1034		call get_param(param_file, mdl, "ADAPT_ALPHA", adaptAlpha, &
1035	0	"Scaling on optimization tendency.", units="nondim", default=1.0)
1036		call get_param(param_file, mdl, "ADAPT_DO_MIN_DEPTH", tmpLogical, &
1037		"If true, make a HyCOM-like mixed layer by preventing interfaces "//&
1038		"from being shallower than the depths specified by the regridding coordinate.", &
1039	0	default=.false.)
1040		call get_param(param_file, mdl, "ADAPT_DRHO0", adaptDrho0, &
1041		"Reference density difference for stratification-dependent diffusion.", &
1042	0	units="kg m-3", default=0.5, scale=US%kg_m3_to_R)
1043
1044		call set_regrid_params(CS, adaptTimeRatio=adaptTimeRatio, adaptZoom=adaptZoom, &
1045		adaptZoomCoeff=adaptZoomCoeff, adaptBuoyCoeff=adaptBuoyCoeff, adaptAlpha=adaptAlpha, &
1046	0	adaptDoMin=tmpLogical, adaptDrho0=adaptDrho0)
1047		endif
1048
1049	2	if (main_parameters .and. coord_is_state_dependent) then
1050		call get_param(param_file, mdl, "MAXIMUM_INT_DEPTH_CONFIG", string, &
1051		"Determines how to specify the maximum interface depths.\n"//&
1052		"Valid options are:\n"//&
1053		" NONE - there are no maximum interface depths\n"//&
1054		" PARAM - use the vector-parameter MAXIMUM_INTERFACE_DEPTHS\n"//&
1055		" FILE:string - read from a file. The string specifies\n"//&
1056		" the filename and variable name, separated\n"//&
1057		" by a comma or space, e.g. FILE:lev.nc,Z\n"//&
1058		" FNC1:string - FNC1:dz_min,H_total,power,precision",&
1059	0	default='NONE')
1060	0	message = "The list of maximum depths for each interface."
1061	0	allocate(z_max(ke+1))
1062	0	allocate(dz_max(ke))
1063	0	if ( trim(string) == "NONE") then
1064		! Do nothing.
1065	0	elseif ( trim(string) == "PARAM") then
1066		call get_param(param_file, mdl, "MAXIMUM_INTERFACE_DEPTHS", z_max, &
1067	0	trim(message), units="m", scale=GV%m_to_H, fail_if_missing=.true.)
1068	0	call set_regrid_max_depths(CS, z_max)
1069	0	elseif (index(trim(string),'FILE:')==1) then
1070	0	if (string(6:6)=='.' .or. string(6:6)=='/') then
1071		! If we specified "FILE:./xyz" or "FILE:/xyz" then we have a relative or absolute path
1072	0	fileName = trim( extractWord(trim(string(6:80)), 1) )
1073		else
1074		! Otherwise assume we should look for the file in INPUTDIR
1075	0	fileName = trim(inputdir) // trim( extractWord(trim(string(6:80)), 1) )
1076		endif
1077	0	if (.not. file_exists(fileName)) call MOM_error(FATAL,trim(mdl)// &
1078		", initialize_regridding: "// &
1079	0	"Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(string)//")")
1080
1081	0	do_sum = .false.
1082	0	varName = trim( extractWord(trim(string(6:)), 2) )
1083	0	if (.not. field_exists(fileName,varName)) call MOM_error(FATAL,trim(mdl)// &
1084		", initialize_regridding: "// &
1085	0	"Specified field not found: Looking for '"//trim(varName)//"' ("//trim(string)//")")
1086	0	if (len_trim(varName)==0) then
1087	0	if (field_exists(fileName,'z_max')) then ; varName = 'z_max'
1088	0	elseif (field_exists(fileName,'dz')) then ; varName = 'dz' ; do_sum = .true.
1089	0	elseif (field_exists(fileName,'dz_max')) then ; varName = 'dz_max' ; do_sum = .true.
1090		else ; call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// &
1091	0	"MAXIMUM_INT_DEPTHS variable not specified and none could be guessed.")
1092		endif
1093		endif
1094	0	if (do_sum) then
1095	0	call MOM_read_data(trim(fileName), trim(varName), dz_max)
1096	0	z_max(1) = 0.0 ; do K=1,ke ; z_max(K+1) = z_max(K) + dz_max(k) ; enddo
1097		else
1098	0	call MOM_read_data(trim(fileName), trim(varName), z_max)
1099		endif
1100		call log_param(param_file, mdl, "!MAXIMUM_INT_DEPTHS", z_max, &
1101	0	trim(message), units=coordinateUnits(coord_mode))
1102	0	call set_regrid_max_depths(CS, z_max, GV%m_to_H)
1103	0	elseif (index(trim(string),'FNC1:')==1) then
1104	0	call dz_function1( trim(string(6:)), dz_max )
1105	0	z_max(1) = 0.0 ; do K=1,ke ; z_max(K+1) = z_max(K) + dz_max(K) ; enddo
1106		call log_param(param_file, mdl, "!MAXIMUM_INT_DEPTHS", z_max, &
1107	0	trim(message), units=coordinateUnits(coord_mode))
1108	0	call set_regrid_max_depths(CS, z_max, GV%m_to_H)
1109		else
1110		call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// &
1111	0	"Unrecognized MAXIMUM_INT_DEPTH_CONFIG "//trim(string))
1112		endif
1113	0	deallocate(z_max)
1114	0	deallocate(dz_max)
1115
1116		! Optionally specify maximum thicknesses for each layer, enforced by moving
1117		! the interface below a layer downward.
1118		call get_param(param_file, mdl, "MAX_LAYER_THICKNESS_CONFIG", longString, &
1119		"Determines how to specify the maximum layer thicknesses.\n"//&
1120		"Valid options are:\n"//&
1121		" NONE - there are no maximum layer thicknesses\n"//&
1122		" PARAM - use the vector-parameter MAX_LAYER_THICKNESS\n"//&
1123		" FILE:string - read from a file. The string specifies\n"//&
1124		" the filename and variable name, separated\n"//&
1125		" by a comma or space, e.g. FILE:lev.nc,Z\n"//&
1126		" FNC1:string - FNC1:dz_min,H_total,power,precision",&
1127	0	default='NONE')
1128	0	message = "The list of maximum thickness for each layer."
1129	0	allocate(h_max(ke))
1130	0	if ( trim(longString) == "NONE") then
1131		! Do nothing.
1132	0	elseif ( trim(longString) == "PARAM") then
1133		call get_param(param_file, mdl, "MAX_LAYER_THICKNESS", h_max, &
1134	0	trim(message), units="m", fail_if_missing=.true., scale=GV%m_to_H)
1135	0	call set_regrid_max_thickness(CS, h_max)
1136	0	elseif (index(trim(longString),'FILE:')==1) then
1137	0	if (longString(6:6)=='.' .or. longString(6:6)=='/') then
1138		! If we specified "FILE:./xyz" or "FILE:/xyz" then we have a relative or absolute path
1139	0	fileName = trim( extractWord(trim(longString(6:200)), 1) )
1140		else
1141		! Otherwise assume we should look for the file in INPUTDIR
1142	0	fileName = trim(inputdir) // trim( extractWord(trim(longString(6:200)), 1) )
1143		endif
1144	0	if (.not. file_exists(fileName)) call MOM_error(FATAL,trim(mdl)// &
1145		", initialize_regridding: "// &
1146	0	"Specified file not found: Looking for '"//trim(fileName)//"' ("//trim(longString)//")")
1147
1148	0	varName = trim( extractWord(trim(longString(6:)), 2) )
1149	0	if (.not. field_exists(fileName,varName)) call MOM_error(FATAL,trim(mdl)// &
1150		", initialize_regridding: "// &
1151	0	"Specified field not found: Looking for '"//trim(varName)//"' ("//trim(longString)//")")
1152	0	if (len_trim(varName)==0) then
1153	0	if (field_exists(fileName,'h_max')) then ; varName = 'h_max'
1154	0	elseif (field_exists(fileName,'dz_max')) then ; varName = 'dz_max'
1155		else ; call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// &
1156	0	"MAXIMUM_INT_DEPTHS variable not specified and none could be guessed.")
1157		endif
1158		endif
1159	0	call MOM_read_data(trim(fileName), trim(varName), h_max)
1160		call log_param(param_file, mdl, "!MAX_LAYER_THICKNESS", h_max, &
1161	0	trim(message), units=coordinateUnits(coord_mode))
1162	0	call set_regrid_max_thickness(CS, h_max, GV%m_to_H)
1163	0	elseif (index(trim(longString),'FNC1:')==1) then
1164	0	call dz_function1( trim(longString(6:)), h_max )
1165		call log_param(param_file, mdl, "!MAX_LAYER_THICKNESS", h_max, &
1166	0	trim(message), units=coordinateUnits(coord_mode))
1167	0	call set_regrid_max_thickness(CS, h_max, GV%m_to_H)
1168		else
1169		call MOM_error(FATAL,trim(mdl)//", initialize_regridding: "// &
1170	0	"Unrecognized MAX_LAYER_THICKNESS_CONFIG "//trim(longString))
1171		endif
1172	0	deallocate(h_max)
1173		endif
1174
1175	2	if (allocated(dz)) deallocate(dz)
1176	4	end subroutine initialize_regridding
1177
1178
1179
1180		!> Deallocation of regridding memory
1181	1	subroutine end_regridding(CS)
1182		type(regridding_CS), intent(inout) :: CS !< Regridding control structure
1183
1184	1	if (associated(CS%zlike_CS)) call end_coord_zlike(CS%zlike_CS)
1185	1	if (associated(CS%sigma_CS)) call end_coord_sigma(CS%sigma_CS)
1186	1	if (associated(CS%rho_CS)) call end_coord_rho(CS%rho_CS)
1187	1	if (associated(CS%hycom_CS)) call end_coord_hycom(CS%hycom_CS)
1188	1	if (associated(CS%adapt_CS)) call end_coord_adapt(CS%adapt_CS)
1189	1	if (associated(CS%hybgen_CS)) call end_hybgen_regrid(CS%hybgen_CS)
1190
1191	1	deallocate( CS%coordinateResolution )
1192	1	if (allocated(CS%coordinateResolution_3d)) deallocate( CS%coordinateResolution_3d )
1193	1	if (allocated(CS%target_density_3d)) deallocate( CS%target_density_3d )
1194	1	if (allocated(CS%target_density)) deallocate( CS%target_density )
1195	1	if (allocated(CS%max_interface_depths) ) deallocate( CS%max_interface_depths )
1196	1	if (allocated(CS%max_layer_thickness) ) deallocate( CS%max_layer_thickness )
1197
1198	1	end subroutine end_regridding
1199
1200		!------------------------------------------------------------------------------
1201		!> Dispatching regridding routine for orchestrating regridding & remapping
1202	39	subroutine regridding_main( remapCS, CS, G, GV, US, h, tv, h_new, dzInterface, &
1203	13	frac_shelf_h, PCM_cell)
1204		!------------------------------------------------------------------------------
1205		! This routine takes care of (1) building a new grid and (2) remapping between
1206		! the old grid and the new grid. The creation of the new grid can be based
1207		! on z coordinates, target interface densities, sigma coordinates or any
1208		! arbitrary coordinate system.
1209		! The MOM6 interface positions are always calculated from the bottom up by
1210		! accumulating the layer thicknesses starting at z=-G%bathyT. z increases
1211		! upwards (decreasing k-index).
1212		! The new grid is defined by the change in position of those interfaces in z
1213		! dzInterface = zNew - zOld.
1214		! Thus, if the regridding inflates the top layer, hNew(1) > hOld(1), then the
1215		! second interface moves downward, zNew(2) < zOld(2), and dzInterface(2) < 0.
1216		! hNew(k) = hOld(k) - dzInterface(k+1) + dzInterface(k)
1217		! IMPORTANT NOTE:
1218		! This is the converse of the sign convention used in the remapping code!
1219		!------------------------------------------------------------------------------
1220
1221		! Arguments
1222		type(remapping_CS), intent(in) :: remapCS !< Remapping parameters and options
1223		type(regridding_CS), intent(in) :: CS !< Regridding control structure
1224		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1225		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
1226		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1227		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Current 3D grid obtained after
1228		!! the last time step [H ~> m or kg m-2]
1229		type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamical variables (T, S, ...)
1230		real, dimension(SZI_(G),SZJ_(G),CS%nk), intent(inout) :: h_new !< New 3D grid consistent with target
1231		!! coordinate [H ~> m or kg m-2]
1232		real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< The change in position of each
1233		!! interface [H ~> m or kg m-2]
1234		real, dimension(SZI_(G),SZJ_(G)), optional, intent(in ) :: frac_shelf_h !< Fractional ice shelf coverage [nondim]
1235		logical, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1236		optional, intent(out ) :: PCM_cell !< Use PCM remapping in cells where true
1237
1238		! Local variables
1239	26	real :: nom_depth_H(SZI_(G),SZJ_(G)) !< The nominal ocean depth at each point in thickness units [H ~> m or kg m-2]
1240	26	real :: tot_h(SZI_(G),SZJ_(G)) !< The total thickness of the water column [H ~> m or kg m-2]
1241	26	real :: tot_dz(SZI_(G),SZJ_(G)) !< The total distance between the top and bottom of the water column [Z ~> m]
1242		real :: Z_to_H ! A conversion factor used by some routines to convert coordinate
1243		! parameters to depth units [H Z-1 ~> nondim or kg m-3]
1244		character(len=128) :: mesg ! A string for error messages
1245		integer :: i, j, k
1246
1247	8553688	if (present(PCM_cell)) PCM_cell(:,:,:) = .false.
1248
1249	13	Z_to_H = US%Z_to_m * GV%m_to_H ! Often this is equivalent to GV%Z_to_H.
1250
1251	13	if ((allocated(tv%SpV_avg)) .and. (tv%valid_SpV_halo < 1)) then
1252	0	if (tv%valid_SpV_halo < 0) then
1253	0	mesg = "invalid values of SpV_avg."
1254		else
1255	0	mesg = "insufficiently large SpV_avg halos of width 0 but 1 is needed."
1256		endif
1257	0	call MOM_error(FATAL, "Regridding_main called in fully non-Boussinesq mode with "//trim(mesg))
1258		endif
1259
1260	13	if (allocated(tv%SpV_avg)) then ! This is the fully non-Boussinesq case
1261	0	do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1
1262	0	tot_h(i,j) = 0.0 ; tot_dz(i,j) = 0.0
1263		enddo ; enddo
1264	0	do k=1,GV%ke ; do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1
1265	0	tot_h(i,j) = tot_h(i,j) + h(i,j,k)
1266	0	tot_dz(i,j) = tot_dz(i,j) + GV%H_to_RZ * tv%SpV_avg(i,j,k) * h(i,j,k)
1267		enddo ; enddo ; enddo
1268	0	do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1
1269	0	if (tot_dz(i,j) > 0.0) then
1270	0	nom_depth_H(i,j) = max(G%meanSL(i,j) + G%bathyT(i,j), 0.0) * (tot_h(i,j) / tot_dz(i,j))
1271		else
1272	0	nom_depth_H(i,j) = 0.0
1273		endif
1274		enddo ; enddo
1275		else
1276	99151	do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1
1277	99138	nom_depth_H(i,j) = max(G%meanSL(i,j) + G%bathyT(i,j), 0.0) * Z_to_H
1278		enddo ; enddo
1279		endif
1280
1281	26	select case ( CS%regridding_scheme )
1282
1283		case ( REGRIDDING_ZSTAR )
1284	13	call build_zstar_grid( CS, G, GV, h, nom_depth_H, dzInterface, frac_shelf_h, zScale=Z_to_H )
1285	13	call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new)
1286		case ( REGRIDDING_SIGMA_SHELF_ZSTAR)
1287	0	call build_zstar_grid( CS, G, GV, h, nom_depth_H, dzInterface, zScale=Z_to_H )
1288	0	call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new)
1289		case ( REGRIDDING_SIGMA )
1290	0	call build_sigma_grid( CS, G, GV, h, nom_depth_H, dzInterface )
1291	0	call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new)
1292		case ( REGRIDDING_RHO )
1293	0	call build_rho_grid( G, GV, G%US, h, nom_depth_H, tv, dzInterface, remapCS, CS, frac_shelf_h )
1294	0	call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new)
1295		case ( REGRIDDING_HYCOM1 )
1296		call build_grid_HyCOM1( G, GV, G%US, h, nom_depth_H, tv, h_new, dzInterface, remapCS, CS, &
1297	0	frac_shelf_h, zScale=Z_to_H )
1298		case ( REGRIDDING_HYBGEN )
1299	0	call hybgen_regrid(G, GV, G%US, h, nom_depth_H, tv, CS%hybgen_CS, dzInterface, PCM_cell)
1300	0	call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new)
1301		case ( REGRIDDING_ADAPTIVE )
1302	0	call build_grid_adaptive(G, GV, G%US, h, nom_depth_H, tv, dzInterface, remapCS, CS)
1303	0	call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new)
1304
1305		case ( REGRIDDING_ARBITRARY )
1306		call MOM_error(FATAL,'MOM_regridding, regridding_main: '//&
1307	0	'Regridding mode "ARB" is not implemented.')
1308		case default
1309		call MOM_error(FATAL,'MOM_regridding, regridding_main: '//&
1310	13	'Unknown regridding scheme selected!')
1311
1312		end select ! type of grid
1313
1314		#ifdef __DO_SAFETY_CHECKS__
1315		if (CS%nk == GV%ke) then
1316		do j = G%jsc-1,G%jec+1 ; do i = G%isc-1,G%iec+1 ; if (G%mask2dT(i,j)>0.) then
1317		call check_grid_column( GV%ke, h(i,j,:), dzInterface(i,j,:), 'in regridding_main')
1318		endif ; enddo ; enddo
1319		endif
1320		#endif
1321	94393	do j=G%jsc,G%jec ; do i=G%isc,G%iec ; if (G%mask2dT(i,j) > 0.) then
1322	4953832	if (minval(h(i,j,:)) < 0.0) then
1323	0	write(0,*) 'regridding_main check_grid: i,j=', i, j, 'h_new(i,j,:)=', h_new(i,j,:)
1324	0	call MOM_error(FATAL, "regridding_main: negative thickness encountered.")
1325		endif
1326		endif ; enddo ; enddo
1327
1328	26	end subroutine regridding_main
1329
1330		!------------------------------------------------------------------------------
1331		!> This routine returns flags indicating which pre-remapping state adjustments
1332		!! are needed depending on the coordinate mode in use.
1333	1	subroutine regridding_preadjust_reqs(CS, do_conv_adj, do_hybgen_unmix, hybgen_CS)
1334
1335		! Arguments
1336		type(regridding_CS), intent(in) :: CS !< Regridding control structure
1337		logical, intent(out) :: do_conv_adj !< Convective adjustment should be done
1338		logical, intent(out) :: do_hybgen_unmix !< Hybgen unmixing should be done
1339		type(hybgen_regrid_CS), pointer, &
1340		optional, intent(out) :: hybgen_CS !< Control structure for hybgen regridding for sharing parameters.
1341
1342
1343	1	do_conv_adj = .false. ; do_hybgen_unmix = .false.
1344	2	select case ( CS%regridding_scheme )
1345
1346		case ( REGRIDDING_ZSTAR, REGRIDDING_SIGMA_SHELF_ZSTAR, REGRIDDING_SIGMA, REGRIDDING_ARBITRARY, &
1347		REGRIDDING_HYCOM1, REGRIDDING_ADAPTIVE )
1348	1	do_conv_adj = .false. ; do_hybgen_unmix = .false.
1349		case ( REGRIDDING_RHO )
1350	0	do_conv_adj = .true. ; do_hybgen_unmix = .false.
1351		case ( REGRIDDING_HYBGEN )
1352	0	do_conv_adj = .false. ; do_hybgen_unmix = CS%use_hybgen_unmix
1353		case default
1354		call MOM_error(FATAL,'MOM_regridding, regridding_preadjust_reqs: '//&
1355	1	'Unknown regridding scheme selected!')
1356		end select ! type of grid
1357
1358	1	if (present(hybgen_CS) .and. do_hybgen_unmix) hybgen_CS => CS%hybgen_CS
1359
1360	1	end subroutine regridding_preadjust_reqs
1361
1362
1363		!> Calculates h_new from h + delta_k dzInterface
1364	13	subroutine calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new)
1365		type(regridding_CS), intent(in) :: CS !< Regridding control structure
1366		type(ocean_grid_type), intent(in) :: G !< Grid structure
1367		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
1368		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Old layer thicknesses [H ~> m or kg m-2]
1369		!! or other units
1370		real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(in) :: dzInterface !< Change in interface positions
1371		!! in the same units as h [H ~> m or kg m-2]
1372		real, dimension(SZI_(G),SZJ_(G),CS%nk), intent(inout) :: h_new !< New layer thicknesses in the same
1373		!! units as h [H ~> m or kg m-2]
1374		! Local variables
1375		integer :: i, j, k, nki
1376
1377	13	nki = min(CS%nk, GV%ke)
1378
1379		!$OMP parallel do default(shared)
1380	819	do j = G%jsc-1,G%jec+1
1381	99151	do i = G%isc-1,G%iec+1
1382	99138	if (G%mask2dT(i,j)>0.) then
1383	4989400	do k=1,nki
1384	4989400	h_new(i,j,k) = max( 0., h(i,j,k) + ( dzInterface(i,j,k) - dzInterface(i,j,k+1) ) )
1385		enddo
1386	65650	if (CS%nk > GV%ke) then
1387	0	do k=nki+1, CS%nk
1388	0	h_new(i,j,k) = max( 0., dzInterface(i,j,k) - dzInterface(i,j,k+1) )
1389		enddo
1390		endif
1391		else
1392	2483832	h_new(i,j,1:nki) = h(i,j,1:nki)
1393	32682	if (CS%nk > GV%ke) h_new(i,j,nki+1:CS%nk) = 0.
1394		! On land points, why are we keeping the original h rather than setting to zero? -AJA
1395		endif
1396		enddo
1397		enddo
1398
1399	13	end subroutine calc_h_new_by_dz
1400
1401
1402		!> Check that the total thickness of new and old grids are consistent
1403	0	subroutine check_grid_column( nk, h, dzInterface, msg )
1404		integer, intent(in) :: nk !< Number of cells
1405		real, dimension(nk), intent(in) :: h !< Cell thicknesses [Z ~> m] or arbitrary units
1406		real, dimension(nk+1), intent(in) :: dzInterface !< Change in interface positions (same units as h), often [Z ~> m]
1407		character(len=*), intent(in) :: msg !< Message to append to errors
1408		! Local variables
1409		integer :: k
1410		real :: eps ! A tiny relative thickness [nondim]
1411		real :: total_h_old ! The total thickness in the old column, in [Z ~> m] or arbitrary units
1412		real :: total_h_new ! The total thickness in the updated column, in [Z ~> m] or arbitrary units
1413		real :: h_new ! A thickness in the updated column, in [Z ~> m] or arbitrary units
1414
1415	0	eps =1. ; eps = epsilon(eps)
1416
1417		! Total thickness of grid h
1418	0	total_h_old = 0.
1419	0	do k = 1,nk
1420	0	total_h_old = total_h_old + h(k)
1421		enddo
1422
1423	0	total_h_new = 0.
1424	0	do k = nk,1,-1
1425	0	h_new = h(k) + ( dzInterface(k) - dzInterface(k+1) ) ! New thickness
1426	0	if (h_new<0.) then
1427	0	write(0,*) 'k,h,hnew=',k,h(k),h_new
1428	0	write(0,*) 'dzI(k+1),dzI(k)=',dzInterface(k+1),dzInterface(k)
1429		call MOM_error( FATAL, 'MOM_regridding, check_grid_column: '//&
1430	0	'Negative layer thickness implied by re-gridding, '//trim(msg))
1431		endif
1432	0	total_h_new = total_h_new + h_new
1433
1434		enddo
1435
1436		! Conservation by implied h_new
1437	0	if (abs(total_h_new-total_h_old)>real(nk-1)0.5(total_h_old+total_h_new)*eps) then
1438	0	write(0,*) 'nk=',nk
1439	0	do k = 1,nk
1440	0	write(0,*) 'k,h,hnew=',k,h(k),h(k)+(dzInterface(k)-dzInterface(k+1))
1441		enddo
1442	0	write(0,*) 'Hold,Hnew,Hnew-Hold=',total_h_old,total_h_new,total_h_new-total_h_old
1443	0	write(0,) 'eps,(n)/2epsH=',eps,real(nk-1)0.5(total_h_old+total_h_new)eps
1444		call MOM_error( FATAL, 'MOM_regridding, check_grid_column: '//&
1445	0	'Re-gridding did NOT conserve total thickness to within roundoff '//trim(msg))
1446		endif
1447
1448		! Check that the top and bottom are intentionally moving
1449	0	if (dzInterface(1) /= 0.) call MOM_error( FATAL, &
1450	0	'MOM_regridding, check_grid_column: Non-zero dzInterface at surface! '//trim(msg))
1451	0	if (dzInterface(nk+1) /= 0.) call MOM_error( FATAL, &
1452	0	'MOM_regridding, check_grid_column: Non-zero dzInterface at bottom! '//trim(msg))
1453
1454	0	end subroutine check_grid_column
1455
1456		!> Returns the change in interface position motion after filtering and
1457		!! assuming the top and bottom interfaces do not move. The filtering is
1458		!! a function of depth, and is applied as the integrated average filtering
1459		!! over the trajectory of the interface. By design, this code can not give
1460		!! tangled interfaces provided that z_old and z_new are not already tangled.
1461	65650	subroutine filtered_grid_motion( CS, nk, z_old, z_new, dz_g )
1462		type(regridding_CS), intent(in) :: CS !< Regridding control structure
1463		integer, intent(in) :: nk !< Number of cells in source grid
1464		real, dimension(nk+1), intent(in) :: z_old !< Old grid position [H ~> m or kg m-2]
1465		real, dimension(CS%nk+1), intent(in) :: z_new !< New grid position before filtering [H ~> m or kg m-2]
1466		real, dimension(CS%nk+1), intent(inout) :: dz_g !< Change in interface positions including
1467		!! the effects of filtering [H ~> m or kg m-2]
1468		! Local variables
1469		real :: sgn ! The sign convention for downward [nondim].
1470		real :: dz_tgt ! The target grid movement of the unfiltered grid [H ~> m or kg m-2]
1471		real :: zr1 ! The old grid position of an interface relative to the surface [H ~> m or kg m-2]
1472		real :: z_old_k ! The corrected position of the old grid [H ~> m or kg m-2]
1473		real :: Aq ! A temporary variable related to the grid weights [nondim]
1474		real :: Bq ! A temporary variable used in the linear term in the quadratic expression for the
1475		! filtered grid movement [H ~> m or kg m-2]
1476		real :: z0, dz0 ! Together these give the position of an interface relative to a reference hieght
1477		! that may be adjusted for numerical accuracy in a solver [H ~> m or kg m-2]
1478		real :: F0 ! An estimated grid movement [H ~> m or kg m-2]
1479		real :: zs ! The depth at which the shallow filtering timescale applies [H ~> m or kg m-2]
1480		real :: zd ! The depth at which the deep filtering timescale applies [H ~> m or kg m-2]
1481		real :: dzwt ! The depth range over which the transition in the filtering timescale occurs [H ~> m or kg m-2]
1482		real :: Idzwt ! The Adcroft reciprocal of dzwt [H-1 ~> m-1 or m2 kg-1]
1483		real :: wtd ! The weight given to the new grid when time filtering [nondim]
1484		real :: Iwtd ! The inverse of wtd [nondim]
1485		real :: Int_zs ! A depth integral of the weights in [H ~> m or kg m-2]
1486		real :: Int_zd ! A depth integral of the weights in [H ~> m or kg m-2]
1487		real :: dInt_zs_zd ! The depth integral of the weights between the deep and shallow depths in [H ~> m or kg m-2]
1488		! For debugging:
1489	65650	real, dimension(nk+1) :: z_act ! The final grid positions after the filtered movement [H ~> m or kg m-2]
1490		! real, dimension(nk+1) :: ddz_g_s, ddz_g_d
1491		logical :: debug = .false.
1492		integer :: k
1493
1494	65650	if ((z_old(nk+1) - z_old(1)) * (z_new(CS%nk+1) - z_new(1)) < 0.0) then
1495	0	call MOM_error(FATAL, "filtered_grid_motion: z_old and z_new use different sign conventions.")
1496	65650	elseif ((z_old(nk+1) - z_old(1)) * (z_new(CS%nk+1) - z_new(1)) == 0.0) then
1497		! This is a massless column, so do nothing and return.
1498	0	do k=1,CS%nk+1 ; dz_g(k) = 0.0 ; enddo ; return
1499	65650	elseif ((z_old(nk+1) - z_old(1)) + (z_new(CS%nk+1) - z_new(1)) > 0.0) then
1500	0	sgn = 1.0
1501		else
1502	65650	sgn = -1.0
1503		endif
1504
1505	65650	if (debug) then
1506	0	do k=2,CS%nk+1
1507	0	if (sgn(z_new(k)-z_new(k-1)) < -5e-16(abs(z_new(k))+abs(z_new(k-1))) ) &
1508	0	call MOM_error(FATAL, "filtered_grid_motion: z_new is tangled.")
1509		enddo
1510	0	do k=2,nk+1
1511	0	if (sgn(z_old(k)-z_old(k-1)) < -5e-16(abs(z_old(k))+abs(z_old(k-1))) ) &
1512	0	call MOM_error(FATAL, "filtered_grid_motion: z_old is tangled.")
1513		enddo
1514		! ddz_g_s(:) = 0.0 ; ddz_g_d(:) = 0.0
1515		endif
1516
1517	65650	zs = CS%depth_of_time_filter_shallow
1518	65650	zd = CS%depth_of_time_filter_deep
1519	65650	wtd = 1.0 - CS%old_grid_weight
1520	65650	Iwtd = 1.0 / wtd
1521
1522	65650	dzwt = (zd - zs)
1523	65650	Idzwt = 0.0 ; if (abs(zd - zs) > 0.0) Idzwt = 1.0 / (zd - zs)
1524	65650	dInt_zs_zd = 0.5(1.0 + Iwtd) (zd - zs)
1525	65650	Aq = 0.5*(Iwtd - 1.0)
1526
1527	65650	dz_g(1) = 0.0
1528	65650	z_old_k = z_old(1)
1529	4989400	do k = 2,CS%nk+1
1530	4923750	if (k<=nk+1) z_old_k = z_old(k) ! This allows for virtual z_old interface at bottom of the model
1531		! zr1 is positive and increases with depth, and dz_tgt is positive downward.
1532	4923750	dz_tgt = sgn*(z_new(k) - z_old_k)
1533	4923750	zr1 = sgn*(z_old_k - z_old(1))
1534
1535		! First, handle the two simple and common cases that do not pass through
1536		! the adjustment rate transition zone.
1537	4989400	if ((zr1 > zd) .and. (zr1 + wtd * dz_tgt > zd)) then
1538	4923750	dz_g(k) = sgn * wtd * dz_tgt
1539	0	elseif ((zr1 < zs) .and. (zr1 + dz_tgt < zs)) then
1540	0	dz_g(k) = sgn * dz_tgt
1541		else
1542		! Find the new value by inverting the equation
1543		! integral(0 to dz_new) Iwt(z) dz = dz_tgt
1544		! This is trivial where Iwt is a constant, and agrees with the two limits above.
1545
1546		! Take test values at the transition points to figure out which segment
1547		! the new value will be found in.
1548	0	if (zr1 >= zd) then
1549	0	Int_zd = Iwtd*(zd - zr1)
1550	0	Int_zs = Int_zd - dInt_zs_zd
1551	0	elseif (zr1 <= zs) then
1552	0	Int_zs = (zs - zr1)
1553	0	Int_zd = dInt_zs_zd + (zs - zr1)
1554		else
1555		! Int_zd = (zd - zr1) * (Iwtd + 0.5(1.0 - Iwtd) (zd - zr1) / (zd - zs))
1556	0	Int_zd = (zd - zr1) * (Iwtd(0.5(zd+zr1) - zs) + 0.5(zd - zr1)) Idzwt
1557	0	Int_zs = (zs - zr1) * (0.5Iwtd ((zr1 - zs)) + (zd - 0.5(zr1+zs))) Idzwt
1558		! It has been verified that Int_zs = Int_zd - dInt_zs_zd to within roundoff.
1559		endif
1560
1561	0	if (dz_tgt >= Int_zd) then ! The new location is in the deep, slow region.
1562	0	dz_g(k) = sgn * ((zd-zr1) + wtd*(dz_tgt - Int_zd))
1563	0	elseif (dz_tgt <= Int_zs) then ! The new location is in the shallow region.
1564	0	dz_g(k) = sgn * ((zs-zr1) + (dz_tgt - Int_zs))
1565		else ! We need to solve a quadratic equation for z_new.
1566		! For accuracy, do the integral from the starting depth or the nearest
1567		! edge of the transition region. The results with each choice are
1568		! mathematically equivalent, but differ in roundoff, and this choice
1569		! should minimize the likelihood of inadvertently overlapping interfaces.
1570	0	if (zr1 <= zs) then ; dz0 = zs-zr1 ; z0 = zs ; F0 = dz_tgt - Int_zs
1571	0	elseif (zr1 >= zd) then ; dz0 = zd-zr1 ; z0 = zd ; F0 = dz_tgt - Int_zd
1572	0	else ; dz0 = 0.0 ; z0 = zr1 ; F0 = dz_tgt ; endif
1573
1574	0	Bq = (dzwt + 2.0Aq(z0-zs))
1575		! Solve the quadratic: Aq(zn-z0)2 + Bq(zn-z0) - F0*dzwt = 0
1576		! Note that b>=0, and the two terms in the standard form cancel for the right root.
1577	0	dz_g(k) = sgn * (dz0 + 2.0F0dzwt / (Bq + sqrt(Bq*2 + 4.0AqF0dzwt) ))
1578
1579		! if (debug) then
1580		! dz0 = zs-zr1 ; z0 = zs ; F0 = dz_tgt - Int_zs ; Bq = (dzwt + 2.0Aq(z0-zs))
1581		! ddz_g_s(k) = sgn * (dz0 + 2.0F0dzwt / (Bq + sqrt(Bq*2 + 4.0AqF0dzwt) )) - dz_g(k)
1582		! dz0 = zd-zr1 ; z0 = zd ; F0 = dz_tgt - Int_zd ; Bq = (dzwt + 2.0Aq(z0-zs))
1583		! ddz_g_d(k) = sgn * (dz0 + 2.0F0dzwt / (Bq + sqrt(Bq*2 + 4.0AqF0dzwt) )) - dz_g(k)
1584		!
1585		! if (abs(ddz_g_s(k)) > 1e-12*(abs(dz_g(k)) + abs(dz_g(k)+ddz_g_s(k)))) &
1586		! call MOM_error(WARNING, "filtered_grid_motion: Expect z_output to be tangled (sc).")
1587		! if (abs(ddz_g_d(k) - ddz_g_s(k)) > 1e-12*(abs(dz_g(k)+ddz_g_d(k)) + abs(dz_g(k)+ddz_g_s(k)))) &
1588		! call MOM_error(WARNING, "filtered_grid_motion: Expect z_output to be tangled.")
1589		! endif
1590		endif
1591
1592		endif
1593		enddo
1594		!dz_g(CS%nk+1) = 0.0
1595
1596	65650	if (debug) then
1597	0	z_old_k = z_old(1)
1598	0	do k=1,CS%nk+1
1599	0	if (k<=nk+1) z_old_k = z_old(k) ! This allows for virtual z_old interface at bottom of the model
1600	0	z_act(k) = z_old_k + dz_g(k)
1601		enddo
1602	0	do k=2,CS%nk+1
1603	0	if (sgn((z_act(k))-z_act(k-1)) < -1e-15(abs(z_act(k))+abs(z_act(k-1))) ) &
1604	0	call MOM_error(FATAL, "filtered_grid_motion: z_output is tangled.")
1605		enddo
1606		endif
1607
1608		end subroutine filtered_grid_motion
1609
1610		!> Builds a z*-coordinate grid with partial steps (Adcroft and Campin, 2004).
1611		!! z* is defined as
1612		!! z* = (z-eta)/(H+eta)H s.t. z=0 when z=eta and z*=-H when z=-H .
1613	0	subroutine build_zstar_grid( CS, G, GV, h, nom_depth_H, dzInterface, frac_shelf_h, zScale)
1614
1615		! Arguments
1616		type(regridding_CS), intent(in) :: CS !< Regridding control structure
1617		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1618		type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
1619		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1620		real, dimension(SZI_(G),SZJ_(G)), intent(in) :: nom_depth_H !< The bathymetric depth of this column
1621		!! relative to mean sea level or another locally
1622		!! valid reference height, converted to thickness
1623		!! units [H ~> m or kg m-2]
1624		real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth
1625		!! [H ~> m or kg m-2].
1626		real, dimension(SZI_(G),SZJ_(G)), optional,intent(in) :: frac_shelf_h !< Fractional
1627		!! ice shelf coverage [nondim].
1628		real, optional, intent(in) :: zScale !< Scaling factor from the target coordinate
1629		!! resolution in Z to desired units for zInterface,
1630		!! usually Z_to_H in which case it is in
1631		!! units of [H Z-1 ~> nondim or kg m-3]
1632		! Local variables
1633		real :: nominalDepth, minThickness, totalThickness ! Depths and thicknesses [H ~> m or kg m-2]
1634		#ifdef __DO_SAFETY_CHECKS__
1635		real :: dh ! The larger of the total column thickness or bathymetric depth [H ~> m or kg m-2]
1636		#endif
1637	26	real, dimension(SZK_(GV)+1) :: zOld ! Previous coordinate interface heights [H ~> m or kg m-2]
1638	26	real, dimension(CS%nk+1) :: zNew ! New coordinate interface heights [H ~> m or kg m-2]
1639		integer :: i, j, k, nz
1640		logical :: ice_shelf
1641
1642	13	nz = GV%ke
1643	13	minThickness = CS%min_thickness
1644	13	ice_shelf = present(frac_shelf_h)
1645
1646		!$OMP parallel do default(none) shared(G,GV,dzInterface,CS,nz,h,frac_shelf_h, &
1647		!$OMP ice_shelf,minThickness,zScale,nom_depth_H) &
1648		!$OMP private(nominalDepth,totalThickness, &
1649		#ifdef __DO_SAFETY_CHECKS__
1650		!$OMP dh, &
1651		#endif
1652		!$OMP zNew,zOld)
1653	819	do j = G%jsc-1,G%jec+1
1654	99151	do i = G%isc-1,G%iec+1
1655
1656	98332	if (G%mask2dT(i,j)==0.) then
1657	2516514	dzInterface(i,j,:) = 0.
1658	32682	cycle
1659		endif
1660
1661		! Local depth (positive downward)
1662	65650	nominalDepth = nom_depth_H(i,j)
1663
1664		! Determine water column thickness
1665	65650	totalThickness = 0.0
1666	4989400	do k = 1,nz
1667	4989400	totalThickness = totalThickness + h(i,j,k)
1668		enddo
1669
1670		! if (GV%Boussinesq) then
1671	65650	zOld(nz+1) = - nominalDepth
1672	4989400	do k = nz,1,-1
1673	4989400	zOld(k) = zOld(k+1) + h(i,j,k)
1674		enddo
1675		! else ! Work downward?
1676		! endif
1677
1678	65650	if (ice_shelf) then
1679	0	if (frac_shelf_h(i,j) > 0.) then ! under ice shelf
1680		call build_zstar_column(CS%zlike_CS, nominalDepth, totalThickness, zNew, &
1681		z_rigid_top=totalThickness-nominalDepth, &
1682	0	eta_orig=zOld(1), zScale=zScale)
1683		else
1684		call build_zstar_column(CS%zlike_CS, nominalDepth, totalThickness, &
1685	0	zNew, zScale=zScale)
1686		endif
1687		else
1688		call build_zstar_column(CS%zlike_CS, nominalDepth, totalThickness, &
1689	65650	zNew, zScale=zScale)
1690		endif
1691
1692		! Calculate the final change in grid position after blending new and old grids
1693	65650	if (CS%use_depth_based_time_filter .or. CS%old_grid_weight>0.) &
1694	65650	call filtered_grid_motion(CS, nz, zOld, zNew, dzInterface(i,j,:))
1695
1696		#ifdef __DO_SAFETY_CHECKS__
1697		dh = max(nominalDepth,totalThickness)
1698		if (abs(zNew(1)-zOld(1)) > (nz-1)0.5epsilon(dh)*dh) then
1699		write(0,*) 'min_thickness=',CS%min_thickness
1700		write(0,*) 'nominalDepth=',nominalDepth,'totalThickness=',totalThickness
1701		write(0,*) 'dzInterface(1) = ', dzInterface(i,j,1), epsilon(dh), nz, CS%nk
1702		do k=1,min(nz,CS%nk)+1
1703		write(0,*) k,zOld(k),zNew(k)
1704		enddo
1705		do k=min(nz,CS%nk)+2,CS%nk+1
1706		write(0,*) k,zOld(nz+1),zNew(k)
1707		enddo
1708		do k=1,min(nz,CS%nk)
1709		write(0,*) k,h(i,j,k),zNew(k)-zNew(k+1),CS%coordinateResolution(k)
1710		enddo
1711		do k=min(nz,CS%nk)+1,CS%nk
1712		write(0,*) k, 0.0, zNew(k)-zNew(k+1), CS%coordinateResolution(k)
1713		enddo
1714		call MOM_error( FATAL, &
1715		'MOM_regridding, build_zstar_grid(): top surface has moved!!!' )
1716		endif
1717		#endif
1718
1719	66456	if (CS%use_adjust_interface_motion) call adjust_interface_motion( CS, nz, h(i,j,:), dzInterface(i,j,:) )
1720
1721		enddo
1722		enddo
1723
1724	13	end subroutine build_zstar_grid
1725
1726		!------------------------------------------------------------------------------
1727		! Build sigma grid
1728		!> This routine builds a grid based on terrain-following coordinates.
1729	0	subroutine build_sigma_grid( CS, G, GV, h, nom_depth_H, dzInterface )
1730		!------------------------------------------------------------------------------
1731		! This routine builds a grid based on terrain-following coordinates.
1732		! The module parameter coordinateResolution(:) determines the resolution in
1733		! sigma coordinate, dSigma(:). sigma-coordinates are defined by
1734		! sigma = (eta-z)/(H+eta) s.t. sigma=0 at z=eta and sigma=1 at z=-H .
1735		!------------------------------------------------------------------------------
1736
1737		! Arguments
1738		type(regridding_CS), intent(in) :: CS !< Regridding control structure
1739		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1740		type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
1741		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1742		real, dimension(SZI_(G),SZJ_(G)), intent(in) :: nom_depth_H !< The bathymetric depth of this column
1743		!! relative to mean sea level or another locally
1744		!! valid reference height, converted to thickness
1745		!! units [H ~> m or kg m-2]
1746		real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth
1747		!! [H ~> m or kg m-2]
1748
1749		! Local variables
1750		real :: nominalDepth ! The nominal depth of the sea-floor in thickness units [H ~> m or kg m-2]
1751		real :: totalThickness ! The total thickness of the water column [H ~> m or kg m-2]
1752		#ifdef __DO_SAFETY_CHECKS__
1753		real :: dh ! The larger of the total column thickness or bathymetric depth [H ~> m or kg m-2]
1754		#endif
1755	0	real, dimension(SZK_(GV)+1) :: zOld ! Previous coordinate interface heights [H ~> m or kg m-2]
1756	0	real, dimension(CS%nk+1) :: zNew ! New coordinate interface heights [H ~> m or kg m-2]
1757		integer :: i, j, k, nz
1758
1759	0	nz = GV%ke
1760
1761	0	do i = G%isc-1,G%iec+1
1762	0	do j = G%jsc-1,G%jec+1
1763
1764	0	if (G%mask2dT(i,j)==0.) then
1765	0	dzInterface(i,j,:) = 0.
1766	0	cycle
1767		endif
1768
1769		! Determine water column height
1770	0	totalThickness = 0.0
1771	0	do k = 1,nz
1772	0	totalThickness = totalThickness + h(i,j,k)
1773		enddo
1774
1775		! In sigma coordinates, the bathymetric depth is only used as an arbitrary offset that
1776		! cancels out when determining coordinate motion, so referencing the column postions to
1777		! the surface is perfectly acceptable, but for preservation of previous answers the
1778		! referencing is done relative to the bottom when in Boussinesq or semi-Boussinesq mode.
1779	0	if (GV%Boussinesq .or. GV%semi_Boussinesq) then
1780	0	nominalDepth = nom_depth_H(i,j)
1781		else
1782	0	nominalDepth = totalThickness
1783		endif
1784
1785	0	call build_sigma_column(CS%sigma_CS, nominalDepth, totalThickness, zNew)
1786
1787		! Calculate the final change in grid position after blending new and old grids
1788	0	zOld(nz+1) = -nominalDepth
1789	0	do k = nz,1,-1
1790	0	zOld(k) = zOld(k+1) + h(i,j,k)
1791		enddo
1792
1793	0	if (CS%use_depth_based_time_filter .or. CS%old_grid_weight>0.) &
1794	0	call filtered_grid_motion(CS, nz, zOld, zNew, dzInterface(i,j,:))
1795
1796		#ifdef __DO_SAFETY_CHECKS__
1797		dh = max(nominalDepth,totalThickness)
1798		if (abs(zNew(1)-zOld(1)) > (CS%nk-1)0.5epsilon(dh)*dh) then
1799		write(0,*) 'min_thickness=',CS%min_thickness
1800		write(0,*) 'nominalDepth=',nominalDepth,'totalThickness=',totalThickness
1801		write(0,*) 'dzInterface(1) = ',dzInterface(i,j,1),epsilon(dh),nz,CS%nk
1802		do k=1,min(nz,CS%nk)+1
1803		write(0,*) k,zOld(k),zNew(k)
1804		enddo
1805		do k=min(nz,CS%nk)+2,CS%nk+1
1806		write(0,*) k,zOld(nz+1),zNew(k)
1807		enddo
1808		do k=1,min(nz,CS%nk)
1809		write(0,) k,h(i,j,k),zNew(k)-zNew(k+1),totalThicknessCS%coordinateResolution(k), &
1810		CS%coordinateResolution(k)
1811		enddo
1812		do k=min(nz,CS%nk)+1,CS%nk
1813		write(0,) k,0.0,zNew(k)-zNew(k+1),totalThicknessCS%coordinateResolution(k), &
1814		CS%coordinateResolution(k)
1815		enddo
1816		call MOM_error( FATAL, &
1817		'MOM_regridding, build_sigma_grid: top surface has moved!!!' )
1818		endif
1819		dzInterface(i,j,1) = 0.
1820		dzInterface(i,j,CS%nk+1) = 0.
1821		#endif
1822
1823		enddo
1824		enddo
1825
1826	0	end subroutine build_sigma_grid
1827
1828		!------------------------------------------------------------------------------
1829		! Build grid based on target interface densities
1830		!------------------------------------------------------------------------------
1831		!> This routine builds a new grid based on a given set of target interface densities.
1832	0	subroutine build_rho_grid( G, GV, US, h, nom_depth_H, tv, dzInterface, remapCS, CS, frac_shelf_h )
1833		!------------------------------------------------------------------------------
1834		! This routine builds a new grid based on a given set of target interface
1835		! densities (these target densities are computed by taking the mean value
1836		! of given layer densities). The algorithm operates as follows within each
1837		! column:
1838		! 1. Given T & S within each layer, the layer densities are computed.
1839		! 2. Based on these layer densities, a global density profile is reconstructed
1840		! (this profile is monotonically increasing and may be discontinuous)
1841		! 3. The new grid interfaces are determined based on the target interface
1842		! densities.
1843		! 4. T & S are remapped onto the new grid.
1844		! 5. Return to step 1 until convergence or until the maximum number of
1845		! iterations is reached, whichever comes first.
1846		!------------------------------------------------------------------------------
1847
1848		! Arguments
1849		type(regridding_CS), intent(in) :: CS !< Regridding control structure
1850		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1851		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
1852		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1853		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1854		real, dimension(SZI_(G),SZJ_(G)), intent(in) :: nom_depth_H !< The bathymetric depth of this column
1855		!! relative to mean sea level or another locally
1856		!! valid reference height, converted to thickness
1857		!! units [H ~> m or kg m-2]
1858		type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure
1859		real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth
1860		!! [H ~> m or kg m-2]
1861		type(remapping_CS), intent(in) :: remapCS !< The remapping control structure
1862		real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: frac_shelf_h !< Fractional ice
1863		!! shelf coverage [nondim]
1864		! Local variables
1865		integer :: nz ! The number of layers in the input grid
1866		integer :: i, j, k
1867		real :: nominalDepth ! Depth of the bottom of the ocean, positive downward [H ~> m or kg m-2]
1868	0	real, dimension(SZK_(GV)+1) :: zOld ! Previous coordinate interface heights [H ~> m or kg m-2]
1869	0	real, dimension(CS%nk+1) :: zNew ! New coordinate interface heights [H ~> m or kg m-2]
1870		real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2]
1871		real :: totalThickness ! Total thicknesses [H ~> m or kg m-2]
1872		#ifdef __DO_SAFETY_CHECKS__
1873		real :: dh ! The larger of the total column thickness or bathymetric depth [H ~> m or kg m-2]
1874		#endif
1875		logical :: ice_shelf
1876
1877	0	h_neglect = set_h_neglect(GV, CS%remap_answer_date, h_neglect_edge)
1878
1879	0	nz = GV%ke
1880	0	ice_shelf = present(frac_shelf_h)
1881
1882	0	if (.not.CS%target_density_set) call MOM_error(FATAL, "build_rho_grid: "//&
1883	0	"Target densities must be set before build_rho_grid is called.")
1884
1885		! Build grid based on target interface densities
1886	0	do j = G%jsc-1,G%jec+1
1887	0	do i = G%isc-1,G%iec+1
1888
1889	0	if (G%mask2dT(i,j)==0.) then
1890	0	dzInterface(i,j,:) = 0.
1891	0	cycle
1892		endif
1893
1894		! Determine total water column thickness
1895	0	totalThickness = 0.0
1896	0	do k=1,nz
1897	0	totalThickness = totalThickness + h(i,j,k)
1898		enddo
1899
1900		! In rho coordinates, the bathymetric depth is only used as an arbitrary offset that
1901		! cancels out when determining coordinate motion, so referencing the column postions to
1902		! the surface is perfectly acceptable, but for preservation of previous answers the
1903		! referencing is done relative to the bottom when in Boussinesq or semi-Boussinesq mode.
1904	0	if (GV%Boussinesq .or. GV%semi_Boussinesq) then
1905	0	nominalDepth = nom_depth_H(i,j)
1906		else
1907	0	nominalDepth = totalThickness
1908		endif
1909
1910		! Determine absolute interface positions
1911	0	zOld(nz+1) = - nominalDepth
1912	0	do k = nz,1,-1
1913	0	zOld(k) = zOld(k+1) + h(i,j,k)
1914		enddo
1915
1916	0	if (ice_shelf) then
1917		call build_rho_column(CS%rho_CS, nz, nominalDepth, h(i,j,:), &
1918		tv%T(i,j,:), tv%S(i,j,:), tv%eqn_of_state, zNew, &
1919		z_rigid_top=totalThickness - nominalDepth, eta_orig = zOld(1), &
1920	0	h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)
1921		else
1922		call build_rho_column(CS%rho_CS, nz, nominalDepth, h(i,j,:), &
1923		tv%T(i,j,:), tv%S(i,j,:), tv%eqn_of_state, zNew, &
1924	0	h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)
1925		endif
1926
1927	0	if (CS%integrate_downward_for_e) then
1928	0	zOld(1) = 0.
1929	0	do k = 1,nz
1930	0	zOld(k+1) = zOld(k) - h(i,j,k)
1931		enddo
1932		else
1933		! The rest of the model defines grids integrating up from the bottom
1934	0	zOld(nz+1) = - nominalDepth
1935	0	do k = nz,1,-1
1936	0	zOld(k) = zOld(k+1) + h(i,j,k)
1937		enddo
1938		endif
1939
1940		! Calculate the final change in grid position after blending new and old grids
1941	0	if (CS%use_depth_based_time_filter .or. CS%old_grid_weight>0.) &
1942	0	call filtered_grid_motion(CS, nz, zOld, zNew, dzInterface(i,j,:))
1943
1944		#ifdef __DO_SAFETY_CHECKS__
1945		do k=2,CS%nk
1946		if (zNew(k) > zOld(1)) then
1947		write(0,*) 'zOld=',zOld
1948		write(0,*) 'zNew=',zNew
1949		call MOM_error( FATAL, 'MOM_regridding, build_rho_grid: '//&
1950		'interior interface above surface!' )
1951		endif
1952		if (zNew(k) > zNew(k-1)) then
1953		write(0,*) 'zOld=',zOld
1954		write(0,*) 'zNew=',zNew
1955		call MOM_error( FATAL, 'MOM_regridding, build_rho_grid: '//&
1956		'interior interfaces cross!' )
1957		endif
1958		enddo
1959
1960		totalThickness = 0.0
1961		do k = 1,nz
1962		totalThickness = totalThickness + h(i,j,k)
1963		enddo
1964
1965		dh = max(nominalDepth, totalThickness)
1966		if (abs(zNew(1)-zOld(1)) > (nz-1)0.5epsilon(dh)*dh) then
1967		write(0,*) 'min_thickness=',CS%min_thickness
1968		write(0,*) 'nominalDepth=',nominalDepth,'totalThickness=',totalThickness
1969		write(0,*) 'zNew(1)-zOld(1) = ',zNew(1)-zOld(1),epsilon(dh),nz
1970		do k=1,min(nz,CS%nk)+1
1971		write(0,*) k,zOld(k),zNew(k)
1972		enddo
1973		do k=min(nz,CS%nk)+2,CS%nk+1
1974		write(0,*) k,zOld(nz+1),zNew(k)
1975		enddo
1976		do k=1,nz
1977		write(0,*) k,h(i,j,k),zNew(k)-zNew(k+1)
1978		enddo
1979		do k=min(nz,CS%nk)+1,CS%nk
1980		write(0,*) k, 0.0, zNew(k)-zNew(k+1), CS%coordinateResolution(k)
1981		enddo
1982		call MOM_error( FATAL, &
1983		'MOM_regridding, build_rho_grid: top surface has moved!!!' )
1984		endif
1985		#endif
1986
1987		enddo ! end loop on i
1988		enddo ! end loop on j
1989
1990	0	end subroutine build_rho_grid
1991
1992		!> Builds a simple HyCOM-like grid with the deepest location of potential
1993		!! density interpolated from the column profile and a clipping of depth for
1994		!! each interface to a fixed z* or p* grid. This should probably be (optionally?)
1995		!! changed to find the nearest location of the target density.
1996		!! \remark { Based on Bleck, 2002: An ocean-ice general circulation model framed in
1997		!! hybrid isopycnic-Cartesian coordinates, Ocean Modelling 37, 55-88.
1998		!! http://dx.doi.org/10.1016/S1463-5003(01)00012-9 }
1999	0	subroutine build_grid_HyCOM1( G, GV, US, h, nom_depth_H, tv, h_new, dzInterface, remapCS, CS, &
2000	0	frac_shelf_h, zScale )
2001		type(ocean_grid_type), intent(in) :: G !< Grid structure
2002		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
2003		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2004		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Existing model thickness [H ~> m or kg m-2]
2005		real, dimension(SZI_(G),SZJ_(G)), intent(in) :: nom_depth_H !< The bathymetric depth of this column
2006		!! relative to mean sea level or another locally
2007		!! valid reference height, converted to thickness
2008		!! units [H ~> m or kg m-2]
2009		type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure
2010		type(remapping_CS), intent(in) :: remapCS !< The remapping control structure
2011		type(regridding_CS), intent(in) :: CS !< Regridding control structure
2012		real, dimension(SZI_(G),SZJ_(G),CS%nk), intent(inout) :: h_new !< New layer thicknesses [H ~> m or kg m-2]
2013		real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< Changes in interface position
2014		!! in thickness units [H ~> m or kg m-2]
2015		real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: frac_shelf_h !< Fractional ice shelf
2016		!! coverage [nondim]
2017		real, optional, intent(in) :: zScale !< Scaling factor from the target coordinate
2018		!! resolution in Z to desired units for zInterface,
2019		!! usually Z_to_H in which case it is in
2020		!! units of [H Z-1 ~> nondim or kg m-3]
2021
2022		! Local variables
2023	0	real, dimension(SZK_(GV)+1) :: z_col ! Source interface positions relative to the surface [H ~> m or kg m-2]
2024	0	real, dimension(SZK_(GV)) :: p_col ! Layer center pressure in the input column [R L2 T-2 ~> Pa]
2025	0	real, dimension(CS%nk+1) :: z_col_new ! New interface positions relative to the surface [H ~> m or kg m-2]
2026	0	real, dimension(CS%nk+1) :: dz_col ! The realized change in z_col [H ~> m or kg m-2]
2027		real :: nominalDepth ! The nominal depth of the seafloor in thickness units [H ~> m or kg m-2]
2028		real :: h_neglect, h_neglect_edge ! Negligible thicknesses used for remapping [H ~> m or kg m-2]
2029		real :: z_top_col ! The nominal height of the sea surface or ice-ocean interface
2030		! in thickness units [H ~> m or kg m-2]
2031		real :: totalThickness ! The total thickness of the water column [H ~> m or kg m-2]
2032		logical :: ice_shelf
2033		integer :: i, j, k, nki
2034
2035	0	h_neglect = set_h_neglect(GV, CS%remap_answer_date, h_neglect_edge)
2036
2037	0	if (.not.CS%target_density_set) call MOM_error(FATAL, "build_grid_HyCOM1 : "//&
2038	0	"Target densities must be set before build_grid_HyCOM1 is called.")
2039
2040	0	nki = min(GV%ke, CS%nk)
2041	0	ice_shelf = present(frac_shelf_h)
2042
2043		! Build grid based on target interface densities
2044	0	do j = G%jsc-1,G%jec+1 ; do i = G%isc-1,G%iec+1
2045	0	if (G%mask2dT(i,j)>0.) then
2046
2047	0	nominalDepth = nom_depth_H(i,j)
2048
2049	0	if (ice_shelf) then
2050	0	totalThickness = 0.0
2051	0	do k=1,GV%ke
2052	0	totalThickness = totalThickness + h(i,j,k)
2053		enddo
2054	0	z_top_col = max(nominalDepth-totalThickness,0.0)
2055		else
2056	0	z_top_col = 0.0
2057		endif
2058
2059	0	z_col(1) = z_top_col ! Work downward rather than bottom up
2060	0	do K = 1, GV%ke
2061	0	z_col(K+1) = z_col(K) + h(i,j,k)
2062		p_col(k) = tv%P_Ref + CS%compressibility_fraction * &
2063	0	( 0.5 * ( z_col(K) + z_col(K+1) ) * (GV%H_to_RZ*GV%g_Earth) - tv%P_Ref )
2064		enddo
2065
2066		call build_hycom1_column(CS%hycom_CS, remapCS, tv%eqn_of_state, GV%ke, i, j, nominalDepth, &
2067		h(i,j,:), tv%T(i,j,:), tv%S(i,j,:), p_col, &
2068		z_col, z_col_new, zScale=zScale, &
2069	0	h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)
2070
2071		! Calculate the final change in grid position after blending new and old grids
2072	0	if (CS%use_depth_based_time_filter .or. CS%old_grid_weight>0.) &
2073	0	call filtered_grid_motion( CS, GV%ke, z_col, z_col_new, dz_col )
2074
2075		! This adjusts things robust to round-off errors
2076	0	dz_col(:) = -dz_col(:)
2077	0	if (CS%use_adjust_interface_motion) call adjust_interface_motion( CS, GV%ke, h(i,j,:), dz_col(:) )
2078
2079	0	dzInterface(i,j,1:nki+1) = dz_col(1:nki+1)
2080	0	if (nki<CS%nk) dzInterface(i,j,nki+2:CS%nk+1) = 0.
2081
2082		else ! on land
2083	0	dzInterface(i,j,:) = 0.
2084		endif ! mask2dT
2085		enddo ; enddo ! i,j
2086
2087	0	call calc_h_new_by_dz(CS, G, GV, h, dzInterface, h_new)
2088
2089	0	end subroutine build_grid_HyCOM1
2090
2091		!> This subroutine builds an adaptive grid that follows density surfaces where
2092		!! possible, subject to constraints on the smoothness of interface heights.
2093	0	subroutine build_grid_adaptive(G, GV, US, h, nom_depth_H, tv, dzInterface, remapCS, CS)
2094		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
2095		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
2096		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2097		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
2098		real, dimension(SZI_(G),SZJ_(G)), intent(in) :: nom_depth_H !< The bathymetric depth of this column
2099		!! relative to mean sea level or another locally
2100		!! valid reference height, converted to thickness
2101		!! units [H ~> m or kg m-2]
2102		type(thermo_var_ptrs), intent(in) :: tv !< A structure pointing to various
2103		!! thermodynamic variables
2104		type(regridding_CS), intent(in) :: CS !< Regridding control structure
2105		real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(inout) :: dzInterface !< The change in interface depth
2106		!! [H ~> m or kg m-2]
2107		type(remapping_CS), intent(in) :: remapCS !< The remapping control structure
2108
2109		! local variables
2110		integer :: i, j, k, nz ! indices and dimension lengths
2111	0	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: tInt ! Temperature on interfaces [C ~> degC]
2112	0	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: sInt ! Salinity on interfaces [S ~> ppt]
2113		! current interface positions and after tendency term is applied
2114		! positive downward
2115	0	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: zInt ! Interface depths [H ~> m or kg m-2]
2116	0	real, dimension(SZK_(GV)+1) :: zNext ! New interface depths [H ~> m or kg m-2]
2117
2118	0	nz = GV%ke
2119
2120		call assert((GV%ke == CS%nk), "build_grid_adaptive is only written to work "//&
2121	0	"with the same number of input and target layers.")
2122
2123		! position surface at z = 0.
2124	0	zInt(:,:,1) = 0.
2125
2126		! work on interior interfaces
2127	0	do K = 2, nz ; do j = G%jsc-2,G%jec+2 ; do i = G%isc-2,G%iec+2
2128	0	tInt(i,j,K) = 0.5 * (tv%T(i,j,k-1) + tv%T(i,j,k))
2129	0	sInt(i,j,K) = 0.5 * (tv%S(i,j,k-1) + tv%S(i,j,k))
2130	0	zInt(i,j,K) = zInt(i,j,K-1) + h(i,j,k-1) ! zInt in [H]
2131		enddo ; enddo ; enddo
2132
2133		! top and bottom temp/salt interfaces are just the layer
2134		! average values
2135	0	tInt(:,:,1) = tv%T(:,:,1) ; tInt(:,:,nz+1) = tv%T(:,:,nz)
2136	0	sInt(:,:,1) = tv%S(:,:,1) ; sInt(:,:,nz+1) = tv%S(:,:,nz)
2137
2138		! set the bottom interface depth
2139	0	zInt(:,:,nz+1) = zInt(:,:,nz) + h(:,:,nz)
2140
2141		! calculate horizontal density derivatives (alpha/beta)
2142		! between cells in a 5-point stencil, columnwise
2143	0	do j = G%jsc-1,G%jec+1 ; do i = G%isc-1,G%iec+1
2144	0	if (G%mask2dT(i,j) < 0.5) then
2145	0	dzInterface(i,j,:) = 0. ! land point, don't move interfaces, and skip
2146	0	cycle
2147		endif
2148
2149		call build_adapt_column(CS%adapt_CS, G, GV, US, tv, i, j, zInt, tInt, sInt, h, &
2150	0	nom_depth_H, zNext)
2151
2152	0	if (CS%use_depth_based_time_filter .or. CS%old_grid_weight>0.) &
2153	0	call filtered_grid_motion(CS, nz, zInt(i,j,:), zNext, dzInterface(i,j,:))
2154		! convert from depth to z
2155	0	do K = 1, nz+1 ; dzInterface(i,j,K) = -dzInterface(i,j,K) ; enddo
2156	0	if (CS%use_adjust_interface_motion) call adjust_interface_motion(CS, nz, h(i,j,:), dzInterface(i,j,:))
2157		enddo ; enddo
2158	0	end subroutine build_grid_adaptive
2159
2160		!> Adjust dz_Interface to ensure non-negative future thicknesses
2161	65650	subroutine adjust_interface_motion( CS, nk, h_old, dz_int )
2162		type(regridding_CS), intent(in) :: CS !< Regridding control structure
2163		integer, intent(in) :: nk !< Number of layers in h_old
2164		real, dimension(nk), intent(in) :: h_old !< Layer thicknesses on the old grid [H ~> m or kg m-2]
2165		real, dimension(CS%nk+1), intent(inout) :: dz_int !< Interface movements, adjusted to keep the thicknesses
2166		!! thicker than their minimum value [H ~> m or kg m-2]
2167		! Local variables
2168		real :: h_new ! A layer thickness on the new grid [H ~> m or kg m-2]
2169		real :: eps ! A tiny relative thickness [nondim]
2170		real :: h_total ! The total thickness of the old grid [H ~> m or kg m-2]
2171		real :: h_err ! An error tolerance that use used to flag unacceptably large negative layer thicknesses
2172		! that can not be explained by roundoff errors [H ~> m or kg m-2]
2173		integer :: k
2174
2175	65650	eps = 1. ; eps = epsilon(eps)
2176
2177	65650	h_total = 0. ; h_err = 0.
2178	4989400	do k = 1, min(CS%nk,nk)
2179	4923750	h_total = h_total + h_old(k)
2180	4923750	h_err = h_err + max( h_old(k), abs(dz_int(k)), abs(dz_int(k+1)) )*eps
2181	4923750	h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) )
2182	4989400	if (h_new < -3.0*h_err) then
2183	0	write(0,*) 'h<0 at k=',k,'h_old=',h_old(k), &
2184	0	'wup=',dz_int(k),'wdn=',dz_int(k+1),'dw_dz=',dz_int(k) - dz_int(k+1), &
2185	0	'h_new=',h_new,'h_err=',h_err
2186		call MOM_error( FATAL, 'MOM_regridding: adjust_interface_motion() - '//&
2187	0	'implied h<0 is larger than roundoff!')
2188		endif
2189		enddo
2190	65650	if (CS%nk>nk) then
2191	0	do k = nk+1, CS%nk
2192	0	h_err = h_err + max( abs(dz_int(k)), abs(dz_int(k+1)) )*eps
2193	0	h_new = ( dz_int(k) - dz_int(k+1) )
2194	0	if (h_new < -3.0*h_err) then
2195	0	write(0,*) 'h<0 at k=',k,'h_old was empty',&
2196	0	'wup=',dz_int(k),'wdn=',dz_int(k+1),'dw_dz=',dz_int(k) - dz_int(k+1), &
2197	0	'h_new=',h_new,'h_err=',h_err
2198		call MOM_error( FATAL, 'MOM_regridding: adjust_interface_motion() - '//&
2199	0	'implied h<0 is larger than roundoff!')
2200		endif
2201		enddo
2202		endif
2203	4923750	do k = min(CS%nk,nk),2,-1
2204	4858100	h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) )
2205	4858100	if (h_new<CS%min_thickness) &
2206	886070	dz_int(k) = ( dz_int(k+1) - h_old(k) ) + CS%min_thickness ! Implies next h_new = min_thickness
2207	4858100	h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) )
2208	4858100	if (h_new<0.) &
2209	0	dz_int(k) = ( 1. - eps ) * ( dz_int(k+1) - h_old(k) ) ! Backup in case min_thickness==0
2210	4858100	h_new = h_old(k) + ( dz_int(k) - dz_int(k+1) )
2211	4923750	if (h_new<0.) then
2212	0	write(0,*) 'h<0 at k=',k,'h_old=',h_old(k), &
2213	0	'wup=',dz_int(k),'wdn=',dz_int(k+1),'dw_dz=',dz_int(k) - dz_int(k+1), &
2214	0	'h_new=',h_new
2215	0	stop 'Still did not work!'
2216		call MOM_error( FATAL, 'MOM_regridding: adjust_interface_motion() - '//&
2217		'Repeated adjustment for roundoff h<0 failed!')
2218		endif
2219		enddo
2220		!if (dz_int(1)/=0.) stop 'MOM_regridding: adjust_interface_motion() surface moved'
2221
2222	65650	end subroutine adjust_interface_motion
2223
2224
2225		!------------------------------------------------------------------------------
2226		!> make sure all layers are at least as thick as the minimum thickness allowed
2227		!! for regridding purposes by inflating thin layers. This breaks mass conservation
2228		!! and adds mass to the model when there are excessively thin layers.
2229	1	subroutine inflate_vanished_layers_old( CS, G, GV, h )
2230		!------------------------------------------------------------------------------
2231		! This routine is called when initializing the regridding options. The
2232		! objective is to make sure all layers are at least as thick as the minimum
2233		! thickness allowed for regridding purposes (this parameter is set in the
2234		! MOM_input file or defaulted to 1.0e-3). When layers are too thin, they
2235		! are inflated up to the minimum thickness.
2236		!------------------------------------------------------------------------------
2237
2238		! Arguments
2239		type(regridding_CS), intent(in) :: CS !< Regridding control structure
2240		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
2241		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
2242		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2]
2243
2244		! Local variables
2245		integer :: i, j, k
2246	2	real :: hTmp(GV%ke) ! A copy of a 1-d column of h [H ~> m or kg m-2]
2247
2248	123	do i = G%isc-1,G%iec+1
2249	7687	do j = G%jsc-1,G%jec+1
2250
2251		! Build grid for current column
2252	574864	do k = 1,GV%ke
2253	574864	hTmp(k) = h(i,j,k)
2254		enddo
2255
2256	7564	call old_inflate_layers_1d( CS%min_thickness, GV%ke, hTmp )
2257
2258		! Save modified grid
2259	574986	do k = 1,GV%ke
2260	574864	h(i,j,k) = hTmp(k)
2261		enddo
2262
2263		enddo
2264		enddo
2265
2266	1	end subroutine inflate_vanished_layers_old
2267
2268		!------------------------------------------------------------------------------
2269		!> Achieve convective adjustment by swapping layers
2270	0	subroutine convective_adjustment(G, GV, h, tv)
2271		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
2272		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
2273		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
2274		intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2]
2275		type(thermo_var_ptrs), intent(inout) :: tv !< A structure pointing to various thermodynamic variables
2276		!------------------------------------------------------------------------------
2277		! Check each water column to see whether it is stratified. If not, sort the
2278		! layers by successive swappings of water masses (bubble sort algorithm)
2279		!------------------------------------------------------------------------------
2280
2281		! Local variables
2282		real :: T0, T1 ! temperatures of two layers [C ~> degC]
2283		real :: S0, S1 ! salinities of two layers [S ~> ppt]
2284		real :: r0, r1 ! densities of two layers [R ~> kg m-3]
2285		real :: h0, h1 ! Layer thicknesses [H ~> m or kg m-2]
2286	0	real, dimension(GV%ke) :: p_col ! A column of zero pressures [R L2 T-2 ~> Pa]
2287	0	real, dimension(GV%ke) :: densities ! Densities in the column [R ~> kg m-3]
2288		logical :: stratified
2289		integer :: i, j, k
2290
2291		!### Doing convective adjustment based on potential densities with zero pressure seems
2292		! questionable, although it does avoid ambiguous sorting. -RWH
2293	0	p_col(:) = 0.
2294
2295		! Loop on columns
2296	0	do j = G%jsc-1,G%jec+1 ; do i = G%isc-1,G%iec+1
2297
2298		! Compute densities within current water column
2299	0	call calculate_density(tv%T(i,j,:), tv%S(i,j,:), p_col, densities, tv%eqn_of_state)
2300
2301		! Repeat restratification until complete
2302	0	do
2303	0	stratified = .true.
2304	0	do k = 1,GV%ke-1
2305		! Gather information of current and next cells
2306	0	T0 = tv%T(i,j,k) ; T1 = tv%T(i,j,k+1)
2307	0	S0 = tv%S(i,j,k) ; S1 = tv%S(i,j,k+1)
2308	0	r0 = densities(k) ; r1 = densities(k+1)
2309	0	h0 = h(i,j,k) ; h1 = h(i,j,k+1)
2310		! If the density of the current cell is larger than the density
2311		! below it, we swap the cells and recalculate the densitiies
2312		! within the swapped cells
2313	0	if ( r0 > r1 ) then
2314	0	tv%T(i,j,k) = T1 ; tv%T(i,j,k+1) = T0
2315	0	tv%S(i,j,k) = S1 ; tv%S(i,j,k+1) = S0
2316	0	h(i,j,k) = h1 ; h(i,j,k+1) = h0
2317		! Recompute densities at levels k and k+1
2318		call calculate_density(tv%T(i,j,k), tv%S(i,j,k), p_col(k), &
2319	0	densities(k), tv%eqn_of_state)
2320		call calculate_density(tv%T(i,j,k+1), tv%S(i,j,k+1), p_col(k+1), &
2321	0	densities(k+1), tv%eqn_of_state )
2322		! Because p_col is has uniform values, these calculate_density calls are equivalent to
2323		! densities(k) = r1 ; densities(k+1) = r0
2324	0	stratified = .false.
2325		endif
2326		enddo ! k
2327
2328	0	if ( stratified ) exit
2329		enddo
2330
2331		enddo ; enddo ! i & j
2332
2333	0	end subroutine convective_adjustment
2334
2335
2336		!------------------------------------------------------------------------------
2337		!> Return a uniform resolution vector in the units of the coordinate
2338	0	function uniformResolution(nk,coordMode,maxDepth,rhoLight,rhoHeavy)
2339		!------------------------------------------------------------------------------
2340		! Calculate a vector of uniform resolution in the units of the coordinate
2341		!------------------------------------------------------------------------------
2342		! Arguments
2343		integer, intent(in) :: nk !< Number of cells in source grid
2344		character(len=*), intent(in) :: coordMode !< A string indicating the coordinate mode.
2345		!! See the documentation for regrid_consts
2346		!! for the recognized values.
2347		real, intent(in) :: maxDepth !< The range of the grid values in some modes, in coordinate
2348		!! dependent units that might be [m] or [kg m-3] or [nondim]
2349		!! or something else.
2350		real, intent(in) :: rhoLight !< The minimum value of the grid in RHO mode [kg m-3]
2351		real, intent(in) :: rhoHeavy !< The maximum value of the grid in RHO mode [kg m-3]
2352
2353		real :: uniformResolution(nk) !< The returned uniform resolution grid, in
2354		!! coordinate dependent units that might be [m] or
2355		!! [kg m-3] or [nondim] or something else.
2356
2357		! Local variables
2358		integer :: scheme
2359
2360	0	scheme = coordinateMode(coordMode)
2361	0	select case ( scheme )
2362
2363		case ( REGRIDDING_ZSTAR, REGRIDDING_HYCOM1, REGRIDDING_HYBGEN, &
2364		REGRIDDING_SIGMA_SHELF_ZSTAR, REGRIDDING_ADAPTIVE )
2365	0	uniformResolution(:) = maxDepth / real(nk)
2366
2367		case ( REGRIDDING_RHO )
2368	0	uniformResolution(:) = (rhoHeavy - rhoLight) / real(nk)
2369
2370		case ( REGRIDDING_SIGMA )
2371	0	uniformResolution(:) = 1. / real(nk)
2372
2373		case default
2374		call MOM_error(FATAL, "MOM_regridding, uniformResolution: "//&
2375	0	"Unrecognized choice for coordinate mode ("//trim(coordMode)//").")
2376
2377		end select ! type of grid
2378
2379		end function uniformResolution
2380
2381		!> Initialize the coordinate resolutions by calling the appropriate initialization
2382		!! routine for the specified coordinate mode.
2383	2	subroutine initCoord(CS, G, GV, US, coord_mode, param_file)
2384		type(regridding_CS), intent(inout) :: CS !< Regridding control structure
2385		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
2386		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
2387		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2388		character(len=*), intent(in) :: coord_mode !< A string indicating the coordinate mode.
2389		!! See the documentation for regrid_consts
2390		!! for the recognized values.
2391		type(param_file_type), intent(in) :: param_file !< Parameter file
2392
2393	2	select case (coordinateMode(coord_mode))
2394		case (REGRIDDING_ZSTAR)
2395	2	call init_coord_zlike(CS%zlike_CS, CS%nk, CS%coordinateResolution)
2396		case (REGRIDDING_SIGMA_SHELF_ZSTAR)
2397	0	call init_coord_zlike(CS%zlike_CS, CS%nk, CS%coordinateResolution)
2398		case (REGRIDDING_SIGMA)
2399	0	call init_coord_sigma(CS%sigma_CS, CS%nk, CS%coordinateResolution)
2400		case (REGRIDDING_RHO)
2401	0	call init_coord_rho(CS%rho_CS, CS%nk, CS%ref_pressure, CS%target_density, CS%interp_CS)
2402		case (REGRIDDING_HYCOM1)
2403	0	if (allocated(CS%coordinateResolution_3d)) then
2404		call init_3d_coord_hycom(CS%hycom_CS, G, CS%nk, &
2405		CS%coordinateResolution_3d, CS%target_density_3d, &
2406	0	CS%interp_CS)
2407		else
2408		call init_coord_hycom(CS%hycom_CS, CS%nk, CS%coordinateResolution, CS%target_density, &
2409	0	CS%interp_CS)
2410		endif
2411		case (REGRIDDING_HYBGEN)
2412	0	call init_hybgen_regrid(CS%hybgen_CS, GV, US, param_file)
2413		case (REGRIDDING_ADAPTIVE)
2414	2	call init_coord_adapt(CS%adapt_CS, CS%nk, CS%coordinateResolution, GV%m_to_H, US%kg_m3_to_R)
2415		end select
2416	2	end subroutine initCoord
2417
2418		!------------------------------------------------------------------------------
2419		!> Set the fixed resolution data
2420	2	subroutine setCoordinateResolution( dz, CS, scale )
2421		real, dimension(:), intent(in) :: dz !< A vector of vertical grid spacings, in arbitrary coordinate
2422		!! dependent units, such as [m] for a z-coordinate or [kg m-3]
2423		!! for a density coordinate.
2424		type(regridding_CS), intent(inout) :: CS !< Regridding control structure
2425		real, optional, intent(in) :: scale !< A scaling factor converting dz to the internal represetation
2426		!! of coordRes, in various units that depend on the coordinate,
2427		!! such as [Z m-1 ~> 1] for a z-coordinate or [R m3 kg-1 ~> 1] for
2428		!! a density coordinate.
2429
2430	2	if (size(dz)/=CS%nk) call MOM_error( FATAL, &
2431	0	'setCoordinateResolution: inconsistent number of levels' )
2432
2433	2	if (present(scale)) then
2434	110	CS%coordinateResolution(:) = scale*dz(:)
2435		else
2436	0	CS%coordinateResolution(:) = dz(:)
2437		endif
2438
2439	2	end subroutine setCoordinateResolution
2440
2441		!> Set the 3d fixed resolution data
2442	0	subroutine setCoordinateResolution_3d( dz_3d, CS, scale )
2443		real, dimension(:,:,:), intent(in) :: dz_3d !< A vector of vertical grid spacings, in arbitrary coordinate
2444		!! dependent units, such as [m] for a z-coordinate or [kg m-3]
2445		!! for a density coordinate.
2446		type(regridding_CS), intent(inout) :: CS !< Regridding control structure
2447		real, optional, intent(in) :: scale !< A scaling factor converting dz to coordRes [Z m-1 ~> 1]
2448
2449	0	if (.not.allocated(CS%coordinateResolution_3d)) &
2450		call MOM_error(FATAL,'setCoordinateResolution_3d: '//&
2451	0	'CS%coordinateResolution_3d not allocated.')
2452
2453	0	if (present(scale)) then
2454	0	CS%coordinateResolution_3d(:,:,:) = scale*dz_3d(:,:,:)
2455		else
2456	0	CS%coordinateResolution_3d(:,:,:) = dz_3d(:,:,:)
2457		endif
2458
2459	0	end subroutine setCoordinateResolution_3d
2460
2461		!> Set target densities based on the old Rlay variable
2462	0	subroutine set_target_densities_from_GV( GV, US, CS )
2463		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
2464		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2465		type(regridding_CS), intent(inout) :: CS !< Regridding control structure
2466		! Local variables
2467		integer :: k, nz
2468
2469	0	nz = CS%nk
2470	0	if (nz == 1) then ! Set a broad range of bounds. Regridding may not be meaningful in this case.
2471	0	CS%target_density(1) = 0.0
2472	0	CS%target_density(2) = 2.0*GV%Rlay(1)
2473		else
2474	0	CS%target_density(1) = (GV%Rlay(1) + 0.5*(GV%Rlay(1)-GV%Rlay(2)))
2475	0	CS%target_density(nz+1) = (GV%Rlay(nz) + 0.5*(GV%Rlay(nz)-GV%Rlay(nz-1)))
2476	0	do k=2,nz
2477	0	CS%target_density(k) = CS%target_density(k-1) + CS%coordinateResolution(k)
2478		enddo
2479		endif
2480	0	CS%target_density_set = .true.
2481
2482	0	end subroutine set_target_densities_from_GV
2483
2484		!> Set target densities based on vector of interface values
2485	0	subroutine set_target_densities_3d( CS, G, scale, rho_int_3d )
2486		type(regridding_CS), intent(inout) :: CS !< Regridding control structure
2487		type(ocean_grid_type),intent(in) :: G !< Ocean grid structure
2488		real, intent(in) :: scale !< A scaling factor converting densities [R m3 kg-1 ~> 1]
2489		real, dimension(SZI_(G),SZJ_(G),CS%nk+1), intent(in) :: rho_int_3d !< Interface densities [kg m-3]
2490
2491	0	if (.not.allocated(CS%target_density_3d)) &
2492		call MOM_error(FATAL,'set_target_densities_3d: '//&
2493	0	'CS%target_density_3d not allocated.')
2494
2495	0	CS%target_density_3d(:,:,:) = scale * rho_int_3d(:,:,:)
2496	0	CS%target_density_set = .true.
2497
2498	0	end subroutine set_target_densities_3d
2499
2500		!> Set target densities based on vector of interface values
2501	0	subroutine set_target_densities( CS, rho_int )
2502		type(regridding_CS), intent(inout) :: CS !< Regridding control structure
2503		real, dimension(CS%nk+1), intent(in) :: rho_int !< Interface densities [R ~> kg m-3]
2504
2505	0	if (size(CS%target_density)/=size(rho_int)) then
2506	0	call MOM_error(FATAL, "set_target_densities inconsistent args!")
2507		endif
2508
2509	0	CS%target_density(:) = rho_int(:)
2510	0	CS%target_density_set = .true.
2511
2512	0	end subroutine set_target_densities
2513
2514		!> Set maximum interface depths based on a vector of input values.
2515	0	subroutine set_regrid_max_depths( CS, max_depths, units_to_H )
2516		type(regridding_CS), intent(inout) :: CS !< Regridding control structure
2517		real, dimension(CS%nk+1), intent(in) :: max_depths !< Maximum interface depths, in arbitrary units, often [m]
2518		real, optional, intent(in) :: units_to_H !< A conversion factor for max_depths into H units,
2519		!! often in [H m-1 ~> 1 or kg m-3]
2520		! Local variables
2521		real :: val_to_H ! A conversion factor from the units for max_depths into H units, often [H m-1 ~> 1 or kg m-3]
2522		! if units_to_H is present, or [nondim] if it is absent.
2523		integer :: K
2524
2525	0	if (.not.allocated(CS%max_interface_depths)) allocate(CS%max_interface_depths(1:CS%nk+1))
2526
2527	0	val_to_H = 1.0 ; if (present(units_to_H)) val_to_H = units_to_H
2528	0	if (max_depths(CS%nk+1) < max_depths(1)) val_to_H = -1.0*val_to_H
2529
2530		! Check for sign reversals in the depths.
2531	0	if (max_depths(CS%nk+1) < max_depths(1)) then
2532	0	do K=1,CS%nk
2533	0	if (max_depths(K+1) > max_depths(K)) &
2534	0	call MOM_error(FATAL, "Unordered list of maximum depths sent to set_regrid_max_depths!")
2535		enddo
2536		else
2537	0	do K=1,CS%nk
2538	0	if (max_depths(K+1) < max_depths(K)) &
2539	0	call MOM_error(FATAL, "Unordered list of maximum depths sent to set_regrid_max_depths.")
2540		enddo
2541		endif
2542
2543	0	do K=1,CS%nk+1
2544	0	CS%max_interface_depths(K) = val_to_H * max_depths(K)
2545		enddo
2546
2547		! set max depths for coordinate
2548	0	select case (CS%regridding_scheme)
2549		case (REGRIDDING_HYCOM1)
2550	0	call set_hycom_params(CS%hycom_CS, max_interface_depths=CS%max_interface_depths)
2551		end select
2552	0	end subroutine set_regrid_max_depths
2553
2554		!> Set maximum layer thicknesses based on a vector of input values.
2555	0	subroutine set_regrid_max_thickness( CS, max_h, units_to_H )
2556		type(regridding_CS), intent(inout) :: CS !< Regridding control structure
2557		real, dimension(CS%nk+1), intent(in) :: max_h !< Maximum layer thicknesses, in arbitrary units, often [m]
2558		real, optional, intent(in) :: units_to_H !< A conversion factor for max_h into H units,
2559		!! often [H m-1 ~> 1 or kg m-3]
2560		! Local variables
2561		real :: val_to_H ! A conversion factor from the units for max_h into H units, often [H m-1 ~> 1 or kg m-3]
2562		! if units_to_H is present, or [nondim] if it is absent.
2563		integer :: k
2564
2565	0	if (.not.allocated(CS%max_layer_thickness)) allocate(CS%max_layer_thickness(1:CS%nk))
2566
2567	0	val_to_H = 1.0 ; if (present( units_to_H)) val_to_H = units_to_H
2568
2569	0	do k=1,CS%nk
2570	0	CS%max_layer_thickness(k) = val_to_H * max_h(k)
2571		enddo
2572
2573		! set max thickness for coordinate
2574	0	select case (CS%regridding_scheme)
2575		case (REGRIDDING_HYCOM1)
2576	0	call set_hycom_params(CS%hycom_CS, max_layer_thickness=CS%max_layer_thickness)
2577		end select
2578	0	end subroutine set_regrid_max_thickness
2579
2580
2581		!> Write the vertical coordinate information into a file.
2582		!! This subroutine writes out a file containing any available data related
2583		!! to the vertical grid used by the MOM ocean model when in ALE mode.
2584	1	subroutine write_regrid_file( CS, GV, filepath )
2585		type(regridding_CS), intent(in) :: CS !< Regridding control structure
2586		type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
2587		character(len=*), intent(in) :: filepath !< The full path to the file to write
2588
2589	31	type(vardesc) :: vars(2)
2590	6	type(MOM_field) :: fields(2)
2591	1	type(MOM_netCDF_file) :: IO_handle ! The I/O handle of the fileset
2592	2	real :: ds(GV%ke), dsi(GV%ke+1) ! The labeling layer and interface coordinates for output
2593		! in axes in files, in coordinate-dependent units that can
2594		! be obtained from getCoordinateUnits [various]
2595
2596	1	if (CS%regridding_scheme == REGRIDDING_HYBGEN) then
2597	0	call write_Hybgen_coord_file(GV, CS%hybgen_CS, filepath)
2598	0	return
2599		endif
2600
2601	76	ds(:) = CS%coord_scale * CS%coordinateResolution(:)
2602	1	dsi(1) = 0.5*ds(1)
2603	75	dsi(2:GV%ke) = 0.5*( ds(1:GV%ke-1) + ds(2:GV%ke) )
2604	1	dsi(GV%ke+1) = 0.5*ds(GV%ke)
2605
2606		vars(1) = var_desc('ds', getCoordinateUnits( CS ), &
2607	6	'Layer Coordinate Thickness', '1', 'L', '1')
2608		vars(2) = var_desc('ds_interface', getCoordinateUnits( CS ), &
2609	6	'Layer Center Coordinate Separation', '1', 'i', '1')
2610
2611		call create_MOM_file(IO_handle, trim(filepath), vars, 2, fields, &
2612	1	SINGLE_FILE, GV=GV)
2613	1	call MOM_write_field(IO_handle, fields(1), ds)
2614	1	call MOM_write_field(IO_handle, fields(2), dsi)
2615	1	call IO_handle%close()
2616
2617	2	end subroutine write_regrid_file
2618
2619		!> Set appropriate values for the negligible thicknesses used for remapping based on an answer date.
2620	193	function set_h_neglect(GV, remap_answer_date, h_neglect_edge) result(h_neglect)
2621		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
2622		integer, intent(in) :: remap_answer_date !< The vintage of the expressions to use
2623		!! for remapping. Values below 20190101 recover the
2624		!! remapping answers from 2018. Higher values use more
2625		!! robust forms of the same remapping algorithms.
2626		real, intent(out) :: h_neglect_edge !< A negligibly small thickness used in
2627		!! remapping edge value calculations [H ~> m or kg m-2]
2628		real :: h_neglect !< A negligibly small thickness used in
2629		!! remapping cell reconstructions [H ~> m or kg m-2]
2630
2631	193	if (remap_answer_date >= 20190101) then
2632	193	h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff
2633	0	elseif (GV%Boussinesq) then
2634	0	h_neglect = GV%m_to_H1.0e-30 ; h_neglect_edge = GV%m_to_H1.0e-10
2635		else
2636	0	h_neglect = GV%kg_m2_to_H1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H1.0e-10
2637		endif
2638	193	end function set_h_neglect
2639
2640		!> Set appropriate values for the negligible vertical layer extents used for remapping based on an answer date.
2641	0	function set_dz_neglect(GV, US, remap_answer_date, dz_neglect_edge) result(dz_neglect)
2642		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
2643		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
2644		integer, intent(in) :: remap_answer_date !< The vintage of the expressions to use
2645		!! for remapping. Values below 20190101 recover the
2646		!! remapping answers from 2018. Higher values use more
2647		!! robust forms of the same remapping algorithms.
2648		real, intent(out) :: dz_neglect_edge !< A negligibly small vertical layer extent
2649		!! used in remapping edge value calculations [Z ~> m]
2650		real :: dz_neglect !< A negligibly small vertical layer extent
2651		!! used in remapping cell reconstructions [Z ~> m]
2652
2653	0	if (remap_answer_date >= 20190101) then
2654	0	dz_neglect = GV%dZ_subroundoff ; dz_neglect_edge = GV%dZ_subroundoff
2655	0	elseif (GV%Boussinesq) then
2656	0	dz_neglect = US%m_to_Z1.0e-30 ; dz_neglect_edge = US%m_to_Z1.0e-10
2657		else
2658	0	dz_neglect = GV%kg_m2_to_H * (GV%H_to_mUS%m_to_Z) 1.0e-30
2659	0	dz_neglect_edge = GV%kg_m2_to_H * (GV%H_to_mUS%m_to_Z) 1.0e-10
2660		endif
2661	0	end function set_dz_neglect
2662
2663		!------------------------------------------------------------------------------
2664		!> Query the fixed resolution data
2665	0	function getCoordinateResolution( CS, undo_scaling )
2666		type(regridding_CS), intent(in) :: CS !< Regridding control structure
2667		logical, optional, intent(in) :: undo_scaling !< If present and true, undo any internal
2668		!! rescaling of the resolution data.
2669		real, dimension(CS%nk) :: getCoordinateResolution !< The resolution or delta of the target coordinate,
2670		!! in units that depend on the coordinate [various]
2671
2672		logical :: unscale
2673	0	unscale = .false. ; if (present(undo_scaling)) unscale = undo_scaling
2674
2675	0	if (unscale) then
2676	0	getCoordinateResolution(:) = CS%coord_scale * CS%coordinateResolution(:)
2677		else
2678	0	getCoordinateResolution(:) = CS%coordinateResolution(:)
2679		endif
2680
2681		end function getCoordinateResolution
2682
2683		!> Query the target coordinate interface positions
2684	2	function getCoordinateInterfaces( CS, undo_scaling )
2685		type(regridding_CS), intent(in) :: CS !< Regridding control structure
2686		logical, optional, intent(in) :: undo_scaling !< If present and true, undo any internal
2687		!! rescaling of the resolution data.
2688		real, dimension(CS%nk+1) :: getCoordinateInterfaces !< Interface positions in target coordinate,
2689		!! in units that depend on the coordinate [various]
2690
2691		integer :: k
2692		logical :: unscale
2693	2	unscale = .false. ; if (present(undo_scaling)) unscale = undo_scaling
2694
2695		! When using a coordinate with target densities, we need to get the actual
2696		! densities, rather than computing the interfaces based on resolution
2697	2	if (CS%regridding_scheme == REGRIDDING_RHO) then
2698	0	if (.not. CS%target_density_set) &
2699		call MOM_error(FATAL, 'MOM_regridding, getCoordinateInterfaces: '//&
2700	0	'target densities not set!')
2701
2702	0	if (unscale) then
2703	0	getCoordinateInterfaces(:) = CS%coord_scale * CS%target_density(:)
2704		else
2705	0	getCoordinateInterfaces(:) = CS%target_density(:)
2706		endif
2707		else
2708	2	if (unscale) then
2709	2	getCoordinateInterfaces(1) = 0.
2710	110	do k = 1, CS%nk
2711		getCoordinateInterfaces(K+1) = getCoordinateInterfaces(K) - &
2712	110	CS%coord_scale * CS%coordinateResolution(k)
2713		enddo
2714		else
2715	0	getCoordinateInterfaces(1) = 0.
2716	0	do k = 1, CS%nk
2717		getCoordinateInterfaces(K+1) = getCoordinateInterfaces(K) - &
2718	0	CS%coordinateResolution(k)
2719		enddo
2720		endif
2721		! The following line has an "abs()" to allow ferret users to reference
2722		! data by index. It is a temporary work around... :( -AJA
2723	112	getCoordinateInterfaces(:) = abs( getCoordinateInterfaces(:) )
2724		endif
2725
2726		end function getCoordinateInterfaces
2727
2728		!------------------------------------------------------------------------------
2729		!> Query the target coordinate units
2730	3	function getCoordinateUnits( CS )
2731		type(regridding_CS), intent(in) :: CS !< Regridding control structure
2732		character(len=20) :: getCoordinateUnits
2733
2734	6	select case ( CS%regridding_scheme )
2735		case ( REGRIDDING_ZSTAR, REGRIDDING_HYCOM1, REGRIDDING_HYBGEN, &
2736		REGRIDDING_ADAPTIVE )
2737	3	getCoordinateUnits = 'meter'
2738		case ( REGRIDDING_SIGMA_SHELF_ZSTAR )
2739	0	getCoordinateUnits = 'meter/fraction'
2740		case ( REGRIDDING_SIGMA )
2741	0	getCoordinateUnits = 'fraction'
2742		case ( REGRIDDING_RHO )
2743	0	getCoordinateUnits = 'kg/m3'
2744		case ( REGRIDDING_ARBITRARY )
2745	0	getCoordinateUnits = 'unknown'
2746		case default
2747		call MOM_error(FATAL,'MOM_regridding, getCoordinateUnits: '//&
2748	3	'Unknown regridding scheme selected!')
2749		end select ! type of grid
2750
2751	3	end function getCoordinateUnits
2752
2753		!------------------------------------------------------------------------------
2754		!> Query the short name of the coordinate
2755	1	function getCoordinateShortName( CS )
2756		type(regridding_CS), intent(in) :: CS !< Regridding control structure
2757		character(len=20) :: getCoordinateShortName
2758
2759	2	select case ( CS%regridding_scheme )
2760		case ( REGRIDDING_ZSTAR )
2761		!getCoordinateShortName = 'z*'
2762		! The following line is a temporary work around... :( -AJA
2763	1	getCoordinateShortName = 'pseudo-depth, -z*'
2764		case ( REGRIDDING_SIGMA_SHELF_ZSTAR )
2765	0	getCoordinateShortName = 'pseudo-depth, -z*/sigma'
2766		case ( REGRIDDING_SIGMA )
2767	0	getCoordinateShortName = 'sigma'
2768		case ( REGRIDDING_RHO )
2769	0	getCoordinateShortName = 'rho'
2770		case ( REGRIDDING_ARBITRARY )
2771	0	getCoordinateShortName = 'coordinate'
2772		case ( REGRIDDING_HYCOM1 )
2773	0	getCoordinateShortName = 'z-rho'
2774		case ( REGRIDDING_HYBGEN )
2775	0	getCoordinateShortName = 'hybrid'
2776		case ( REGRIDDING_ADAPTIVE )
2777	0	getCoordinateShortName = 'adaptive'
2778		case default
2779		call MOM_error(FATAL,'MOM_regridding, getCoordinateShortName: '//&
2780	1	'Unknown regridding scheme selected!')
2781		end select ! type of grid
2782
2783	1	end function getCoordinateShortName
2784
2785		!> Can be used to set any of the parameters for MOM_regridding.
2786	23	subroutine set_regrid_params( CS, boundary_extrapolation, min_thickness, old_grid_weight, &
2787		use_depth_based_time_filter, depth_of_time_filter_shallow, depth_of_time_filter_deep, &
2788		interp_scheme, use_adjust_interface_motion, compress_fraction, ref_pressure, &
2789		integrate_downward_for_e, remap_answers_2018, remap_answer_date, regrid_answer_date, &
2790		adaptTimeRatio, adaptZoom, adaptZoomCoeff, adaptBuoyCoeff, &
2791		adaptAlpha, adaptDoMin, adaptDrho0)
2792		type(regridding_CS), intent(inout) :: CS !< Regridding control structure
2793		logical, optional, intent(in) :: boundary_extrapolation !< Extrapolate in boundary cells
2794		real, optional, intent(in) :: min_thickness !< Minimum thickness allowed when building the
2795		!! new grid [H ~> m or kg m-2]
2796		real, optional, intent(in) :: old_grid_weight !< Weight given to old coordinate when time-filtering grid [nondim]
2797		logical, optional, intent(in) :: use_depth_based_time_filter !< Allow depth-based time filtering
2798		real, optional, intent(in) :: depth_of_time_filter_shallow !< Depth to start cubic [H ~> m or kg m-2]
2799		real, optional, intent(in) :: depth_of_time_filter_deep !< Depth to end cubic [H ~> m or kg m-2]
2800		character(len=*), optional, intent(in) :: interp_scheme !< Interpolation method for state-dependent coordinates
2801		logical, optional, intent(in) :: use_adjust_interface_motion !< Call adjust_interface_motion()
2802		real, optional, intent(in) :: compress_fraction !< Fraction of compressibility to add to potential density [nondim]
2803		real, optional, intent(in) :: ref_pressure !< The reference pressure for density-dependent
2804		!! coordinates [R L2 T-2 ~> Pa]
2805		logical, optional, intent(in) :: integrate_downward_for_e !< If true, integrate for interface positions downward
2806		!! from the top.
2807		logical, optional, intent(in) :: remap_answers_2018 !< If true, use the order of arithmetic and expressions
2808		!! that recover the remapping answers from 2018. Otherwise
2809		!! use more robust but mathematically equivalent expressions.
2810		integer, optional, intent(in) :: remap_answer_date !< The vintage of the expressions to use for remapping
2811		integer, optional, intent(in) :: regrid_answer_date !< The vintage of the expressions to use for regridding
2812		real, optional, intent(in) :: adaptTimeRatio !< Ratio of the ALE timestep to the grid timescale [nondim].
2813		real, optional, intent(in) :: adaptZoom !< Depth of near-surface zooming region [H ~> m or kg m-2].
2814		real, optional, intent(in) :: adaptZoomCoeff !< Coefficient of near-surface zooming diffusivity [nondim].
2815		real, optional, intent(in) :: adaptBuoyCoeff !< Coefficient of buoyancy diffusivity [nondim].
2816		real, optional, intent(in) :: adaptAlpha !< Scaling factor on optimization tendency [nondim].
2817		logical, optional, intent(in) :: adaptDoMin !< If true, make a HyCOM-like mixed layer by
2818		!! preventing interfaces from being shallower than
2819		!! the depths specified by the regridding coordinate.
2820		real, optional, intent(in) :: adaptDrho0 !< Reference density difference for stratification-dependent
2821		!! diffusion. [R ~> kg m-3]
2822
2823	0	if (present(interp_scheme)) call set_interp_scheme(CS%interp_CS, interp_scheme)
2824	23	if (present(boundary_extrapolation)) call set_interp_extrap(CS%interp_CS, boundary_extrapolation)
2825	23	if (present(regrid_answer_date)) call set_interp_answer_date(CS%interp_CS, regrid_answer_date)
2826
2827	23	if (present(old_grid_weight)) then
2828	12	if (old_grid_weight<0. .or. old_grid_weight>1.) &
2829	0	call MOM_error(FATAL,'MOM_regridding, set_regrid_params: Weight is out side the range 0..1!')
2830	12	CS%old_grid_weight = old_grid_weight
2831		endif
2832	23	if (present(use_depth_based_time_filter)) CS%use_depth_based_time_filter = &
2833	2	use_depth_based_time_filter
2834	23	if (present(depth_of_time_filter_shallow)) CS%depth_of_time_filter_shallow = &
2835	1	depth_of_time_filter_shallow
2836	23	if (present(depth_of_time_filter_deep)) CS%depth_of_time_filter_deep = &
2837	1	depth_of_time_filter_deep
2838	23	if (present(depth_of_time_filter_shallow) .or. present(depth_of_time_filter_deep)) then
2839	1	if (CS%depth_of_time_filter_deep<CS%depth_of_time_filter_shallow) call MOM_error(FATAL, &
2840		'MOM_regridding, '//&
2841	0	'set_regrid_params: depth_of_time_filter_deep<depth_of_time_filter_shallow!')
2842		endif
2843
2844	23	if (present(min_thickness)) CS%min_thickness = min_thickness
2845	23	if (present(use_adjust_interface_motion)) CS%use_adjust_interface_motion = use_adjust_interface_motion
2846	23	if (present(compress_fraction)) CS%compressibility_fraction = compress_fraction
2847	23	if (present(ref_pressure)) CS%ref_pressure = ref_pressure
2848	23	if (present(integrate_downward_for_e)) CS%integrate_downward_for_e = integrate_downward_for_e
2849	23	if (present(remap_answers_2018)) then
2850	0	if (remap_answers_2018) then
2851	0	CS%remap_answer_date = 20181231
2852		else
2853	0	CS%remap_answer_date = 20190101
2854		endif
2855		endif
2856	23	if (present(remap_answer_date)) CS%remap_answer_date = remap_answer_date
2857
2858	46	select case (CS%regridding_scheme)
2859		case (REGRIDDING_ZSTAR)
2860	23	if (present(min_thickness)) call set_zlike_params(CS%zlike_CS, min_thickness=min_thickness)
2861		case (REGRIDDING_SIGMA_SHELF_ZSTAR)
2862	0	if (present(min_thickness)) call set_zlike_params(CS%zlike_CS, min_thickness=min_thickness)
2863		case (REGRIDDING_SIGMA)
2864	0	if (present(min_thickness)) call set_sigma_params(CS%sigma_CS, min_thickness=min_thickness)
2865		case (REGRIDDING_RHO)
2866	0	if (present(min_thickness)) call set_rho_params(CS%rho_CS, min_thickness=min_thickness)
2867	0	if (present(ref_pressure)) call set_rho_params(CS%rho_CS, ref_pressure=ref_pressure)
2868	0	if (present(integrate_downward_for_e)) &
2869	0	call set_rho_params(CS%rho_CS, integrate_downward_for_e=integrate_downward_for_e)
2870	0	if (associated(CS%rho_CS) .and. (present(interp_scheme) .or. &
2871		present(boundary_extrapolation))) &
2872	0	call set_rho_params(CS%rho_CS, interp_CS=CS%interp_CS)
2873		case (REGRIDDING_HYCOM1)
2874	0	if (associated(CS%hycom_CS) .and. (present(interp_scheme) .or. &
2875		present(boundary_extrapolation))) &
2876	0	call set_hycom_params(CS%hycom_CS, interp_CS=CS%interp_CS)
2877		case (REGRIDDING_HYBGEN)
2878		! Do nothing for now.
2879		case (REGRIDDING_ADAPTIVE)
2880	0	if (present(adaptTimeRatio)) call set_adapt_params(CS%adapt_CS, adaptTimeRatio=adaptTimeRatio)
2881	0	if (present(adaptZoom)) call set_adapt_params(CS%adapt_CS, adaptZoom=adaptZoom)
2882	0	if (present(adaptZoomCoeff)) call set_adapt_params(CS%adapt_CS, adaptZoomCoeff=adaptZoomCoeff)
2883	0	if (present(adaptBuoyCoeff)) call set_adapt_params(CS%adapt_CS, adaptBuoyCoeff=adaptBuoyCoeff)
2884	0	if (present(adaptAlpha)) call set_adapt_params(CS%adapt_CS, adaptAlpha=adaptAlpha)
2885	0	if (present(adaptDoMin)) call set_adapt_params(CS%adapt_CS, adaptDoMin=adaptDoMin)
2886	23	if (present(adaptDrho0)) call set_adapt_params(CS%adapt_CS, adaptDrho0=adaptDrho0)
2887		end select
2888
2889	23	end subroutine set_regrid_params
2890
2891		!> Returns the number of levels/layers in the regridding control structure
2892	1	integer function get_regrid_size(CS)
2893		type(regridding_CS), intent(inout) :: CS !< Regridding control structure
2894
2895	1	get_regrid_size = CS%nk
2896
2897	1	end function get_regrid_size
2898
2899		!> This returns a copy of the zlike_CS stored in the regridding control structure.
2900	974650	function get_zlike_CS(CS)
2901		type(regridding_CS), intent(in) :: CS !< Regridding control structure
2902		type(zlike_CS) :: get_zlike_CS
2903
2904	974650	get_zlike_CS = CS%zlike_CS
2905	1949300	end function get_zlike_CS
2906
2907		!> This returns a copy of the sigma_CS stored in the regridding control structure.
2908	0	function get_sigma_CS(CS)
2909		type(regridding_CS), intent(in) :: CS !< Regridding control structure
2910		type(sigma_CS) :: get_sigma_CS
2911
2912	0	get_sigma_CS = CS%sigma_CS
2913	0	end function get_sigma_CS
2914
2915		!> This returns a copy of the rho_CS stored in the regridding control structure.
2916	0	function get_rho_CS(CS)
2917		type(regridding_CS), intent(in) :: CS !< Regridding control structure
2918		type(rho_CS) :: get_rho_CS
2919
2920	0	get_rho_CS = CS%rho_CS
2921	0	end function get_rho_CS
2922
2923		!------------------------------------------------------------------------------
2924		!> Return coordinate-derived thicknesses for fixed coordinate systems
2925	0	function getStaticThickness( CS, SSH, depth )
2926		type(regridding_CS), intent(in) :: CS !< Regridding control structure
2927		real, intent(in) :: SSH !< The sea surface height, in the same units as depth, often [Z ~> m]
2928		real, intent(in) :: depth !< The maximum depth of the grid, often [Z ~> m]
2929		real, dimension(CS%nk) :: getStaticThickness !< The returned thicknesses in the units of
2930		!! depth, often [Z ~> m]
2931		! Local
2932		integer :: k
2933		real :: z, dz ! Vertical positions and grid spacing [Z ~> m]
2934
2935	0	select case ( CS%regridding_scheme )
2936		case ( REGRIDDING_ZSTAR, REGRIDDING_SIGMA_SHELF_ZSTAR, REGRIDDING_HYCOM1, &
2937		REGRIDDING_HYBGEN, REGRIDDING_ADAPTIVE )
2938	0	if (depth>0.) then
2939	0	z = ssh
2940	0	do k = 1, CS%nk
2941	0	dz = CS%coordinateResolution(k) * ( 1. + ssh/depth ) ! Nominal dz*
2942	0	dz = max(dz, 0.) ! Avoid negative incase ssh=-depth
2943	0	dz = min(dz, depth - z) ! Clip if below topography
2944	0	z = z + dz ! Bottom of layer
2945	0	getStaticThickness(k) = dz
2946		enddo
2947		else
2948	0	getStaticThickness(:) = 0. ! On land ...
2949		endif
2950		case ( REGRIDDING_SIGMA )
2951	0	getStaticThickness(:) = CS%coordinateResolution(:) * ( depth + ssh )
2952		case ( REGRIDDING_RHO )
2953	0	getStaticThickness(:) = 0. ! Not applicable
2954		case ( REGRIDDING_ARBITRARY )
2955	0	getStaticThickness(:) = 0. ! Not applicable
2956		case default
2957		call MOM_error(FATAL,'MOM_regridding, getStaticThickness: '//&
2958	0	'Unknown regridding scheme selected!')
2959		end select ! type of grid
2960
2961		end function getStaticThickness
2962
2963		!> Parses a string and generates a dz(:) profile that goes like k**power.
2964	1	subroutine dz_function1( string, dz )
2965		character(len=*), intent(in) :: string !< String with list of parameters in form
2966		!! dz_min, H_total, power, precision
2967		real, dimension(:), intent(inout) :: dz !< Profile of nominal thicknesses [m] or other units
2968		! Local variables
2969		integer :: nk, k
2970		real :: dz_min ! minimum grid spacing [m] or other units
2971		real :: power ! A power to raise the relative position in index space [nondim]
2972		real :: prec ! The precision with which positions are returned [m] or other units
2973		real :: H_total ! The sum of the nominal thicknesses [m] or other units
2974
2975	1	nk = size(dz) ! Number of cells
2976	1	prec = -1024.
2977	1	read( string, *) dz_min, H_total, power, prec
2978	1	if (prec == -1024.) call MOM_error(FATAL,"dz_function1: "// &
2979	0	"Problem reading FNC1: string ="//trim(string))
2980		! Create profile of ( dz - dz_min )
2981	76	do k = 1, nk
2982	76	dz(k) = (real(k-1)/real(nk-1))**power
2983		enddo
2984	151	dz(:) = ( H_total - real(nk) * dz_min ) * ( dz(:) / sum(dz) ) ! Rescale to so total is H_total
2985	76	dz(:) = anint( dz(:) / prec ) * prec ! Rounds to precision prec
2986	151	dz(:) = ( H_total - real(nk) * dz_min ) * ( dz(:) / sum(dz) ) ! Rescale to so total is H_total
2987	76	dz(:) = anint( dz(:) / prec ) * prec ! Rounds to precision prec
2988	76	dz(nk) = dz(nk) + ( H_total - sum( dz(:) + dz_min ) ) ! Adjust bottommost layer
2989	76	dz(:) = anint( dz(:) / prec ) * prec ! Rounds to precision prec
2990	76	dz(:) = dz(:) + dz_min ! Finally add in the constant dz_min
2991
2992	1	end subroutine dz_function1
2993
2994		!> Construct the name of a parameter for a specific coordinate based on param_prefix and param_suffix. For the main,
2995		!! prognostic coordinate this will simply return the parameter name (e.g. P_REF)
2996	2	function create_coord_param(param_prefix, param_name, param_suffix) result(coord_param)
2997		character(len=*) :: param_name !< The base name of the parameter (e.g. the one used for the main coordinate)
2998		character(len=*) :: param_prefix !< String to prefix to parameter names.
2999		character(len=*) :: param_suffix !< String to append to parameter names.
3000		character(len=MAX_PARAM_LENGTH) :: coord_param !< Parameter name prepended by param_prefix
3001		!! and appended with param_suffix
3002		integer :: out_length
3003
3004	2	if (len_trim(param_prefix) + len_trim(param_suffix) == 0) then
3005	0	coord_param = param_name
3006		else
3007		! Note the +2 is because of two underscores
3008	2	out_length = len_trim(param_name)+len_trim(param_prefix)+len_trim(param_suffix)+2
3009	2	if (out_length > MAX_PARAM_LENGTH) then
3010	0	call MOM_error(FATAL,"Coordinate parameter is too long; increase MAX_PARAM_LENGTH")
3011		endif
3012	2	coord_param = TRIM(param_prefix)//"_"//TRIM(param_name)//"_"//TRIM(param_suffix)
3013		endif
3014
3015	2	end function create_coord_param
3016
3017		!> Parses a string and generates a rho_target(:) profile with refined resolution downward
3018		!! and returns the number of levels
3019	0	integer function rho_function1( string, rho_target )
3020		character(len=*), intent(in) :: string !< String with list of parameters in form
3021		!! dz_min, H_total, power, precision
3022		real, dimension(:), allocatable, intent(inout) :: rho_target !< Profile of interface densities [kg m-3]
3023		! Local variables
3024		integer :: nki, k, nk
3025		real :: dx ! Fractional distance from interface nki [nondim]
3026		real :: ddx ! Change in dx between interfaces [nondim]
3027		real :: rho_1, rho_2 ! Density of the top two layers in a profile [kg m-3]
3028		real :: rho_3 ! Density in the third layer, below which the density increase linearly
3029		! in subsequent layers [kg m-3]
3030		real :: drho ! Change in density over the linear region [kg m-3]
3031		real :: rho_4 ! The densest density in this profile [kg m-3], which might be very large.
3032		real :: drho_min ! A minimal fractional density difference [nondim]?
3033
3034	0	read( string, *) nk, rho_1, rho_2, rho_3, drho, rho_4, drho_min
3035	0	allocate(rho_target(nk+1))
3036	0	nki = nk + 1 - 4 ! Number of interfaces minus 4 specified values
3037	0	rho_target(1) = rho_1
3038	0	rho_target(2) = rho_2
3039	0	dx = 0.
3040	0	do k = 0, nki
3041	0	ddx = max( drho_min, real(nki-k)/real(nki*nki) )
3042	0	dx = dx + ddx
3043	0	rho_target(3+k) = rho_3 + (2. * drho) * dx
3044		enddo
3045	0	rho_target(nki+4) = rho_4
3046
3047	0	rho_function1 = nk
3048
3049	0	end function rho_function1
3050
3051		!> \namespace mom_regridding
3052		!!
3053		!! A vertical grid is defined solely by the cell thicknesses, \f$h\f$.
3054		!! Most calculations in this module start with the coordinate at the bottom
3055		!! of the column set to -depth, and use a increasing value of coordinate with
3056		!! decreasing k. This is consistent with the rest of MOM6 that uses position,
3057		!! \f$z\f$ which is a negative quantity for most of the ocean.
3058		!!
3059		!! A change in grid is define through a change in position of the interfaces:
3060		!! \f[
3061		!! z^n_{k+1/2} = z^{n-1}_{k+1/2} + \Delta z_{k+1/2}
3062		!! \f]
3063		!! with the positive upward coordinate convention
3064		!! \f[
3065		!! z_{k-1/2} = z_{k+1/2} + h_k
3066		!! \f]
3067		!! so that
3068		!! \f[
3069		!! h^n_k = h^{n-1}_k + ( \Delta z_{k-1/2} - \Delta z_{k+1/2} )
3070		!! \f]
3071		!!
3072		!! Original date of creation: 2008.06.09 by L. White
3073
3074	0	end module MOM_regridding