src/ALE/MOM_ALE.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		!> This module contains the main regridding routines.
6		!!
7		!! Regridding comprises two steps:
8		!! 1. Interpolation and creation of a new grid based on target interface
9		!! densities (or any other criterion).
10		!! 2. Remapping of quantities between old grid and new grid.
11		!!
12		!! Original module written by Laurent White, 2008.06.09
13		module MOM_ALE
14
15		use MOM_debugging, only : check_column_integrals
16		use MOM_diag_mediator, only : register_diag_field, post_data, diag_ctrl
17		use MOM_diag_mediator, only : time_type, diag_update_remap_grids, query_averaging_enabled
18		use MOM_domains, only : create_group_pass, do_group_pass, group_pass_type
19		use MOM_error_handler, only : MOM_error, FATAL, WARNING
20		use MOM_error_handler, only : callTree_showQuery
21		use MOM_error_handler, only : callTree_enter, callTree_leave, callTree_waypoint
22		use MOM_hybgen_unmix, only : hybgen_unmix, init_hybgen_unmix, end_hybgen_unmix, hybgen_unmix_CS
23		use MOM_hybgen_regrid, only : hybgen_regrid_CS
24		use MOM_file_parser, only : get_param, param_file_type, log_param
25		use MOM_interface_heights,only : find_eta, calc_derived_thermo
26		use MOM_open_boundary, only : ocean_OBC_type, OBC_DIRECTION_E, OBC_DIRECTION_W
27		use MOM_open_boundary, only : OBC_DIRECTION_N, OBC_DIRECTION_S
28		use MOM_regridding, only : initialize_regridding, regridding_main, end_regridding
29		use MOM_regridding, only : uniformResolution
30		use MOM_regridding, only : inflate_vanished_layers_old
31		use MOM_regridding, only : regridding_preadjust_reqs, convective_adjustment
32		use MOM_regridding, only : set_target_densities_from_GV, set_target_densities
33		use MOM_regridding, only : regriddingCoordinateModeDoc, DEFAULT_COORDINATE_MODE
34		use MOM_regridding, only : regriddingInterpSchemeDoc, regriddingDefaultInterpScheme
35		use MOM_regridding, only : regriddingDefaultBoundaryExtrapolation
36		use MOM_regridding, only : regriddingDefaultMinThickness
37		use MOM_regridding, only : regridding_CS, set_regrid_params, write_regrid_file
38		use MOM_regridding, only : getCoordinateInterfaces
39		use MOM_regridding, only : getCoordinateUnits, getCoordinateShortName
40		use MOM_regridding, only : getStaticThickness
41		use MOM_remapping, only : initialize_remapping, end_remapping
42		use MOM_remapping, only : remapping_core_h, remapping_core_w
43		use MOM_remapping, only : remappingSchemesDoc, remappingDefaultScheme
44		use MOM_remapping, only : interpolate_column, reintegrate_column
45		use MOM_remapping, only : remapping_CS, dzFromH1H2, remapping_set_param
46		use MOM_string_functions, only : uppercase, extractWord, extract_integer
47		use MOM_tracer_registry, only : tracer_registry_type, tracer_type, MOM_tracer_chkinv
48		use MOM_unit_scaling, only : unit_scale_type
49		use MOM_variables, only : ocean_grid_type, thermo_var_ptrs
50		use MOM_verticalGrid, only : get_thickness_units, verticalGrid_type
51
52		!use regrid_consts, only : coordinateMode, DEFAULT_COORDINATE_MODE
53		use regrid_consts, only : coordinateUnits, coordinateMode, state_dependent
54		use regrid_edge_values, only : edge_values_implicit_h4
55		use PLM_functions, only : PLM_reconstruction, PLM_boundary_extrapolation
56		use PLM_functions, only : PLM_extrapolate_slope, PLM_monotonized_slope, PLM_slope_wa
57		use PPM_functions, only : PPM_reconstruction, PPM_boundary_extrapolation
58		use Recon1d_PLM_WLS, only : PLM_WLS
59
60		implicit none ; private
61		#include <MOM_memory.h>
62
63
64		!> ALE control structure
65		type, public :: ALE_CS ; private
66		logical :: remap_uv_using_old_alg !< If true, uses the old "remapping via a delta z"
67		!! method. If False, uses the new method that
68		!! remaps between grids described by h.
69		logical :: partial_cell_vel_remap !< If true, use partial cell thicknesses at velocity points
70		!! that are masked out where they extend below the shallower
71		!! of the neighboring bathymetry for remapping velocity.
72
73		real :: regrid_time_scale !< The time-scale used in blending between the current (old) grid
74		!! and the target (new) grid [T ~> s]
75
76		type(regridding_CS) :: regridCS !< Regridding parameters and work arrays
77		type(remapping_CS) :: remapCS !< Remapping parameters and work arrays
78		type(remapping_CS) :: vel_remapCS !< Remapping parameters for velocities and work arrays
79
80		type(hybgen_unmix_CS), pointer :: hybgen_unmixCS => NULL() !< Parameters for hybgen remapping
81
82		logical :: use_hybgen_unmix !< If true, use the hybgen unmixing code before regridding
83		logical :: do_conv_adj !< If true, do convective adjustment before regridding
84
85		integer :: nk !< Used only for queries, not directly by this module
86		real :: BBL_h_vel_mask !< The thickness of a bottom boundary layer within which velocities in
87		!! thin layers are zeroed out after remapping, following practice with
88		!! Hybgen remapping, or a negative value to avoid such filtering
89		!! altogether, in [H ~> m or kg m-2].
90		real :: h_vel_mask !< A thickness at velocity points below which near-bottom layers are
91		!! zeroed out after remapping, following the practice with Hybgen
92		!! remapping, or a negative value to avoid such filtering altogether,
93		!! in [H ~> m or kg m-2].
94
95		logical :: remap_after_initialization !< Indicates whether to regrid/remap after initializing the state.
96
97		integer :: answer_date !< The vintage of the expressions and order of arithmetic to use for
98		!! remapping. Values below 20190101 result in the use of older, less
99		!! accurate expressions that were in use at the end of 2018. Higher
100		!! values result in the use of more robust and accurate forms of
101		!! mathematically equivalent expressions.
102
103		logical :: conserve_ke !< Apply a correction to the baroclinic velocity after remapping to
104		!! conserve KE.
105
106		logical :: debug !< If true, write verbose checksums for debugging purposes.
107		logical :: show_call_tree !< For debugging
108
109		! for diagnostics
110		type(diag_ctrl), pointer :: diag !< structure to regulate output
111		integer, dimension(:), allocatable :: id_tracer_remap_tendency !< diagnostic id
112		integer, dimension(:), allocatable :: id_Htracer_remap_tendency !< diagnostic id
113		integer, dimension(:), allocatable :: id_Htracer_remap_tendency_2d !< diagnostic id
114		logical, dimension(:), allocatable :: do_tendency_diag !< flag for doing diagnostics
115		integer :: id_dzRegrid = -1 !< diagnostic id
116
117		! diagnostic for fields prior to applying ALE remapping
118		integer :: id_u_preale = -1 !< diagnostic id for zonal velocity before ALE.
119		integer :: id_v_preale = -1 !< diagnostic id for meridional velocity before ALE.
120		integer :: id_h_preale = -1 !< diagnostic id for layer thicknesses before ALE.
121		integer :: id_T_preale = -1 !< diagnostic id for temperatures before ALE.
122		integer :: id_S_preale = -1 !< diagnostic id for salinities before ALE.
123		integer :: id_e_preale = -1 !< diagnostic id for interface heights before ALE.
124		integer :: id_vert_remap_h = -1 !< diagnostic id for layer thicknesses used for remapping
125		integer :: id_vert_remap_h_tendency = -1 !< diagnostic id for layer thickness tendency due to ALE
126		integer :: id_remap_delta_integ_u2 = -1 !< Change in depth-integrated rho0u*2/2
127		integer :: id_remap_delta_integ_v2 = -1 !< Change in depth-integrated rho0v*2/2
128
129		end type
130
131		! Publicly available functions
132		public ALE_init
133		public ALE_end
134		public ALE_regrid
135		public ALE_offline_inputs
136		public ALE_regrid_accelerated
137		public ALE_remap_scalar
138		public ALE_remap_tracers
139		public ALE_remap_velocities
140		public ALE_remap_set_h_vel, ALE_remap_set_h_vel_via_dz
141		public ALE_remap_interface_vals
142		public ALE_remap_vertex_vals
143		public ALE_PLM_edge_values
144		public TS_PLM_edge_values
145		public TS_PPM_edge_values
146		public TS_PLM_WLS_edge_values
147		public adjustGridForIntegrity
148		public ALE_initRegridding
149		public ALE_getCoordinate
150		public ALE_getCoordinateUnits
151		public ALE_writeCoordinateFile
152		public ALE_updateVerticalGridType
153		public ALE_initThicknessToCoord
154		public ALE_update_regrid_weights
155		public pre_ALE_diagnostics
156		public pre_ALE_adjustments
157		public ALE_remap_init_conds
158		public ALE_register_diags
159		public ALE_set_extrap_boundaries
160
161		! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
162		! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
163		! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]". If the units
164		! vary with the Boussinesq approximation, the Boussinesq variant is given first.
165
166		contains
167
168		!> This routine is typically called (from initialize_MOM in file MOM.F90)
169		!! before the main time integration loop to initialize the regridding stuff.
170		!! We read the MOM_input file to register the values of different
171		!! regridding/remapping parameters.
172	1	subroutine ALE_init( param_file, G, GV, US, max_depth, CS)
173		type(param_file_type), intent(in) :: param_file !< Parameter file
174		type(ocean_grid_type), intent(in) :: G !< Grid structure
175		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
176		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
177		real, intent(in) :: max_depth !< The maximum depth of the ocean [Z ~> m].
178		type(ALE_CS), pointer :: CS !< Module control structure
179
180		! Local variables
181		character(len=40) :: mdl = "MOM_ALE" ! This module's name.
182		character(len=80) :: string, vel_string ! Temporary strings
183		real :: filter_shallow_depth, filter_deep_depth ! Depth ranges of filtering [H ~> m or kg m-2]
184		integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.
185		logical :: check_reconstruction
186		logical :: check_remapping
187		logical :: force_bounds_in_subcell
188		logical :: local_logical
189		logical :: remap_boundary_extrap
190		logical :: init_boundary_extrap
191		logical :: om4_remap_via_sub_cells
192		type(hybgen_regrid_CS), pointer :: hybgen_regridCS => NULL() ! Control structure for hybgen regridding
193		! for sharing parameters.
194		real :: h_neglect, h_neglect_edge ! small thicknesses [H ~> m or kg m-2]
195
196	1	if (associated(CS)) then
197		call MOM_error(WARNING, "ALE_init called with an associated "// &
198	0	"control structure.")
199	0	return
200		endif
201	1	allocate(CS)
202
203	1	CS%show_call_tree = callTree_showQuery()
204	1	if (CS%show_call_tree) call callTree_enter("ALE_init(), MOM_ALE.F90")
205
206		call get_param(param_file, mdl, "REMAP_UV_USING_OLD_ALG", CS%remap_uv_using_old_alg, &
207		"If true, uses the old remapping-via-a-delta-z method for "//&
208		"remapping u and v. If false, uses the new method that remaps "//&
209		"between grids described by an old and new thickness.", &
210	1	default=.false.)
211
212		! Initialize and configure regridding
213	1	call ALE_initRegridding(G, GV, US, max_depth, param_file, mdl, CS%regridCS)
214		call regridding_preadjust_reqs(CS%regridCS, CS%do_conv_adj, CS%use_hybgen_unmix, &
215	1	hybgen_CS=hybgen_regridCS)
216
217		! Initialize and configure remapping that is orchestrated by ALE.
218		call get_param(param_file, mdl, "REMAPPING_SCHEME", string, &
219		"This sets the reconstruction scheme used "//&
220		"for vertical remapping for all variables. "//&
221		"It can be one of the following schemes: \n"//&
222	1	trim(remappingSchemesDoc), default=remappingDefaultScheme)
223		call get_param(param_file, mdl, "VELOCITY_REMAPPING_SCHEME", vel_string, &
224		"This sets the reconstruction scheme used for vertical remapping "//&
225		"of velocities. By default it is the same as REMAPPING_SCHEME. "//&
226		"It can be one of the following schemes: \n"//&
227	1	trim(remappingSchemesDoc), default=trim(string))
228		call get_param(param_file, mdl, "FATAL_CHECK_RECONSTRUCTIONS", check_reconstruction, &
229		"If true, cell-by-cell reconstructions are checked for "//&
230		"consistency and if non-monotonicity or an inconsistency is "//&
231	1	"detected then a FATAL error is issued.", default=.false.)
232		call get_param(param_file, mdl, "FATAL_CHECK_REMAPPING", check_remapping, &
233		"If true, the results of remapping are checked for "//&
234		"conservation and new extrema and if an inconsistency is "//&
235	1	"detected then a FATAL error is issued.", default=.false.)
236		call get_param(param_file, mdl, "REMAP_BOUND_INTERMEDIATE_VALUES", force_bounds_in_subcell, &
237		"If true, the values on the intermediate grid used for remapping "//&
238		"are forced to be bounded, which might not be the case due to "//&
239	1	"round off.", default=.false.)
240		call get_param(param_file, mdl, "REMAP_BOUNDARY_EXTRAP", remap_boundary_extrap, &
241		"If true, values at the interfaces of boundary cells are "//&
242	1	"extrapolated instead of piecewise constant", default=.false.)
243		call get_param(param_file, mdl, "INIT_BOUNDARY_EXTRAP", init_boundary_extrap, &
244		"If true, values at the interfaces of boundary cells are "//&
245		"extrapolated instead of piecewise constant during initialization. "//&
246	1	"Defaults to REMAP_BOUNDARY_EXTRAP.", default=remap_boundary_extrap)
247		call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
248		"This sets the default value for the various _ANSWER_DATE parameters.", &
249	1	default=99991231)
250		call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
251		"This selects the remapping algorithm used in OM4 that does not use "//&
252		"the full reconstruction for the top- and lower-most sub-layers, but instead "//&
253		"assumes they are always vanished (untrue) and so just uses their edge values. "//&
254	1	"We recommend setting this option to false.", default=.true.)
255		call get_param(param_file, mdl, "REMAPPING_ANSWER_DATE", CS%answer_date, &
256		"The vintage of the expressions and order of arithmetic to use for remapping. "//&
257		"Values below 20190101 result in the use of older, less accurate expressions "//&
258		"that were in use at the end of 2018. Higher values result in the use of more "//&
259		"robust and accurate forms of mathematically equivalent expressions.", &
260	1	default=default_answer_date, do_not_log=.not.GV%Boussinesq)
261	1	if (.not.GV%Boussinesq) CS%answer_date = max(CS%answer_date, 20230701)
262
263	1	if (CS%answer_date >= 20190101) then
264	1	h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff
265	0	elseif (GV%Boussinesq) then
266	0	h_neglect = GV%m_to_H * 1.0e-30 ; h_neglect_edge = GV%m_to_H * 1.0e-10
267		else
268	0	h_neglect = GV%kg_m2_to_H * 1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H * 1.0e-10
269		endif
270
271		call initialize_remapping( CS%remapCS, string, nk=GV%ke, &
272		boundary_extrapolation=init_boundary_extrap, &
273		check_reconstruction=check_reconstruction, &
274		check_remapping=check_remapping, &
275		force_bounds_in_subcell=force_bounds_in_subcell, &
276		om4_remap_via_sub_cells=om4_remap_via_sub_cells, &
277		answer_date=CS%answer_date, &
278	1	h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)
279		call initialize_remapping( CS%vel_remapCS, vel_string, nk=GV%ke, &
280		boundary_extrapolation=init_boundary_extrap, &
281		check_reconstruction=check_reconstruction, &
282		check_remapping=check_remapping, &
283		force_bounds_in_subcell=force_bounds_in_subcell, &
284		om4_remap_via_sub_cells=om4_remap_via_sub_cells, &
285		answer_date=CS%answer_date, &
286	1	h_neglect=h_neglect, h_neglect_edge=h_neglect_edge)
287
288		call get_param(param_file, mdl, "PARTIAL_CELL_VELOCITY_REMAP", CS%partial_cell_vel_remap, &
289		"If true, use partial cell thicknesses at velocity points that are masked out "//&
290		"where they extend below the shallower of the neighboring bathymetry for "//&
291	1	"remapping velocity.", default=.false.)
292
293		call get_param(param_file, mdl, "REMAP_AFTER_INITIALIZATION", CS%remap_after_initialization, &
294		"If true, applies regridding and remapping immediately after "//&
295		"initialization so that the state is ALE consistent. This is a "//&
296		"legacy step and should not be needed if the initialization is "//&
297	1	"consistent with the coordinate mode.", default=.true.)
298
299		call get_param(param_file, mdl, "REGRID_USE_DEPTH_BASED_TIME_FILTER", local_logical, &
300		"If true, always uses depth-based time filtering code that updates the "//&
301		"generated grid using REGRID_TIME_SCALE, REGRID_FILTER_SHALLOW_DEPTH, "//&
302		"REGRID_FILTER_DEEP_DEPTH parameters. Setting to True always uses "//&
303		"filtering but setting to False bypasses calculations when filter times = 0.", &
304	1	default=.true.)
305	1	call set_regrid_params(CS%regridCS, use_depth_based_time_filter=local_logical)
306		call get_param(param_file, mdl, "REGRID_TIME_SCALE", CS%regrid_time_scale, &
307		"The time-scale used in blending between the current (old) grid "//&
308		"and the target (new) grid. A short time-scale favors the target "//&
309		"grid (0. or anything less than DT_THERM) has no memory of the old "//&
310		"grid. A very long time-scale makes the model more Lagrangian.", &
311	1	units="s", default=0., scale=US%s_to_T)
312		call get_param(param_file, mdl, "REGRID_FILTER_SHALLOW_DEPTH", filter_shallow_depth, &
313		"The depth above which no time-filtering is applied. Above this depth "//&
314		"final grid exactly matches the target (new) grid.", &
315	1	units="m", default=0., scale=GV%m_to_H)
316		call get_param(param_file, mdl, "REGRID_FILTER_DEEP_DEPTH", filter_deep_depth, &
317		"The depth below which full time-filtering is applied with time-scale "//&
318		"REGRID_TIME_SCALE. Between depths REGRID_FILTER_SHALLOW_DEPTH and "//&
319		"REGRID_FILTER_DEEP_DEPTH the filter weights adopt a cubic profile.", &
320	1	units="m", default=0., scale=GV%m_to_H)
321		call set_regrid_params(CS%regridCS, depth_of_time_filter_shallow=filter_shallow_depth, &
322	1	depth_of_time_filter_deep=filter_deep_depth)
323		call get_param(param_file, mdl, "REGRID_USE_OLD_DIRECTION", local_logical, &
324		"If true, the regridding integrates upwards from the bottom for "//&
325		"interface positions, much as the main model does. If false "//&
326		"regridding integrates downward, consistent with the remapping code.", &
327	1	default=.true., do_not_log=.true.)
328	1	call set_regrid_params(CS%regridCS, integrate_downward_for_e=.not.local_logical)
329
330		call get_param(param_file, mdl, "REMAP_VEL_MASK_BBL_THICK", CS%BBL_h_vel_mask, &
331		"A thickness of a bottom boundary layer below which velocities in thin layers "//&
332		"are zeroed out after remapping, following practice with Hybgen remapping, "//&
333		"or a negative value to avoid such filtering altogether.", &
334	1	default=-0.001, units="m", scale=GV%m_to_H)
335		call get_param(param_file, mdl, "REMAP_VEL_MASK_H_THIN", CS%h_vel_mask, &
336		"A thickness at velocity points below which near-bottom layers are zeroed out "//&
337		"after remapping, following practice with Hybgen remapping, "//&
338		"or a negative value to avoid such filtering altogether.", &
339	1	default=1.0e-6, units="m", scale=GV%m_to_H, do_not_log=(CS%BBL_h_vel_mask<=0.0))
340
341	1	if (CS%use_hybgen_unmix) &
342	0	call init_hybgen_unmix(CS%hybgen_unmixCS, GV, US, param_file, hybgen_regridCS)
343
344		call get_param(param_file, mdl, "REMAP_VEL_CONSERVE_KE", CS%conserve_ke, &
345		"If true, a correction is applied to the baroclinic component of velocity "//&
346		"after remapping so that total KE is conserved. KE may not be conserved "//&
347		"when (CS%BBL_h_vel_mask > 0.0) .and. (CS%h_vel_mask > 0.0)", &
348	1	default=.false.)
349		call get_param(param_file, "MOM", "DEBUG", CS%debug, &
350		"If true, write out verbose debugging data.", &
351	1	default=.false., debuggingParam=.true.)
352
353		! Keep a record of values for subsequent queries
354	1	CS%nk = GV%ke
355
356	1	if (CS%show_call_tree) call callTree_leave("ALE_init()")
357		end subroutine ALE_init
358
359		!> Sets the boundary extrapolation set for the remapping type.
360	1	subroutine ALE_set_extrap_boundaries( param_file, CS)
361		type(param_file_type), intent(in) :: param_file !< Parameter file
362		type(ALE_CS), pointer :: CS !< Module control structure
363
364		logical :: remap_boundary_extrap
365		call get_param(param_file, "MOM_ALE", "REMAP_BOUNDARY_EXTRAP", remap_boundary_extrap, &
366		"If true, values at the interfaces of boundary cells are "//&
367	1	"extrapolated instead of piecewise constant", default=.false.)
368	1	call remapping_set_param(CS%remapCS, boundary_extrapolation=remap_boundary_extrap)
369	1	end subroutine ALE_set_extrap_boundaries
370
371		!> Sets the remapping algorithm to that of OM4
372		!!
373		!! The remapping aglorithm used in OM4 made poor assumptions about the reconstructions
374		!! in the top/bottom layers, namely that they were always vanished and could be
375		!! represented solely by their upper/lower edge value respectively.
376		!! Passing .false. here uses the full reconstruction of those top and bottom layers
377		!! and properly sample those layers.
378	0	subroutine ALE_set_OM4_remap_algorithm( CS, om4_remap_via_sub_cells )
379		type(ALE_CS), pointer :: CS !< Module control structure
380		logical, intent(in) :: om4_remap_via_sub_cells !< If true, use OM4 remapping algorithm
381
382	0	call remapping_set_param(CS%remapCS, om4_remap_via_sub_cells=om4_remap_via_sub_cells )
383
384	0	end subroutine ALE_set_OM4_remap_algorithm
385
386		!> Initialize diagnostics for the ALE module.
387	1	subroutine ALE_register_diags(Time, G, GV, US, diag, CS)
388		type(time_type),target, intent(in) :: Time !< Time structure
389		type(ocean_grid_type), intent(in) :: G !< Grid structure
390		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
391		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
392		type(diag_ctrl), target, intent(in) :: diag !< Diagnostics control structure
393		type(ALE_CS), pointer :: CS !< Module control structure
394
395		! Local variables
396		character(len=48) :: thickness_units
397
398	1	CS%diag => diag
399	1	thickness_units = get_thickness_units(GV)
400
401		! These diagnostics of the state variables before ALE are useful for
402		! debugging the ALE code.
403		CS%id_u_preale = register_diag_field('ocean_model', 'u_preale', diag%axesCuL, Time, &
404	1	'Zonal velocity before remapping', 'm s-1', conversion=US%L_T_to_m_s)
405		CS%id_v_preale = register_diag_field('ocean_model', 'v_preale', diag%axesCvL, Time, &
406	1	'Meridional velocity before remapping', 'm s-1', conversion=US%L_T_to_m_s)
407		CS%id_h_preale = register_diag_field('ocean_model', 'h_preale', diag%axesTL, Time, &
408		'Layer Thickness before remapping', thickness_units, conversion=GV%H_to_MKS, &
409	1	v_extensive=.true.)
410		CS%id_T_preale = register_diag_field('ocean_model', 'T_preale', diag%axesTL, Time, &
411	1	'Temperature before remapping', 'degC', conversion=US%C_to_degC)
412		CS%id_S_preale = register_diag_field('ocean_model', 'S_preale', diag%axesTL, Time, &
413	1	'Salinity before remapping', 'PSU', conversion=US%S_to_ppt)
414		CS%id_e_preale = register_diag_field('ocean_model', 'e_preale', diag%axesTi, Time, &
415	1	'Interface Heights before remapping', 'm', conversion=US%Z_to_m)
416
417		CS%id_dzRegrid = register_diag_field('ocean_model', 'dzRegrid', diag%axesTi, Time, &
418	1	'Change in interface height due to ALE regridding', 'm', conversion=GV%H_to_m)
419		CS%id_vert_remap_h = register_diag_field('ocean_model', 'vert_remap_h', diag%axestl, Time, &
420		'layer thicknesses after ALE regridding and remapping', &
421	1	thickness_units, conversion=GV%H_to_MKS, v_extensive=.true.)
422		CS%id_vert_remap_h_tendency = register_diag_field('ocean_model', &
423		'vert_remap_h_tendency', diag%axestl, Time, &
424		'Layer thicknesses tendency due to ALE regridding and remapping', &
425	1	trim(thickness_units)//" s-1", conversion=GV%H_to_MKS*US%s_to_T, v_extensive=.true.)
426		CS%id_remap_delta_integ_u2 = register_diag_field('ocean_model', 'ale_u2', diag%axesCu1, Time, &
427		'Rate of change in half rho0 times depth integral of squared zonal '//&
428		'velocity by remapping. If REMAP_VEL_CONSERVE_KE is .true. then '//&
429		'this measures the change before the KE-conserving correction is applied.', &
430	1	'W m-2', conversion=US%RZ3_T3_to_W_m2GV%H_to_RZUS%L_to_Z**2)
431		CS%id_remap_delta_integ_v2 = register_diag_field('ocean_model', 'ale_v2', diag%axesCv1, Time, &
432		'Rate of change in half rho0 times depth integral of squared meridional '//&
433		'velocity by remapping. If REMAP_VEL_CONSERVE_KE is .true. then '//&
434		'this measures the change before the KE-conserving correction is applied.', &
435	1	'W m-2', conversion=US%RZ3_T3_to_W_m2GV%H_to_RZUS%L_to_Z**2)
436
437	1	end subroutine ALE_register_diags
438
439		!> Crudely adjust (initial) grid for integrity.
440		!! This routine is typically called (from initialize_MOM in file MOM.F90)
441		!! before the main time integration loop to initialize the regridding stuff.
442		!! We read the MOM_input file to register the values of different
443		!! regridding/remapping parameters.
444	1	subroutine adjustGridForIntegrity( CS, G, GV, h )
445		type(ALE_CS), intent(in) :: CS !< Regridding parameters and options
446		type(ocean_grid_type), intent(in) :: G !< Ocean grid informations
447		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
448		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Current 3D grid thickness that
449		!! are to be adjusted [H ~> m or kg m-2]
450	1	call inflate_vanished_layers_old( CS%regridCS, G, GV, h(:,:,:) )
451
452	1	end subroutine adjustGridForIntegrity
453
454
455		!> End of regridding (memory deallocation).
456		!! This routine is typically called (from MOM_end in file MOM.F90)
457		!! after the main time integration loop to deallocate the regridding stuff.
458	1	subroutine ALE_end(CS)
459		type(ALE_CS), pointer :: CS !< module control structure
460
461		! Deallocate memory used for the regridding
462	1	call end_remapping( CS%remapCS )
463
464	1	if (CS%use_hybgen_unmix) call end_hybgen_unmix( CS%hybgen_unmixCS )
465	1	call end_regridding( CS%regridCS )
466
467	1	deallocate(CS)
468
469	1	end subroutine ALE_end
470
471		!> Save any diagnostics of the state before ALE remapping. These diagnostics are
472		!! mostly used for debugging.
473	12	subroutine pre_ALE_diagnostics(G, GV, US, h, u, v, tv, CS)
474		type(ocean_grid_type), intent(in) :: G !< Ocean grid informations
475		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
476		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
477		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Current 3D grid obtained after the
478		!! last time step [H ~> m or kg m-2]
479		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: u !< Zonal velocity field [L T-1 ~> m s-1]
480		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: v !< Meridional velocity field [L T-1 ~> m s-1]
481		type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamic variable structure
482		type(ALE_CS), pointer :: CS !< Regridding parameters and options
483
484		! Local variables
485	24	real :: eta_preale(SZI_(G),SZJ_(G),SZK_(GV)+1) ! Interface heights before remapping [Z ~> m]
486
487	12	if (CS%id_u_preale > 0) call post_data(CS%id_u_preale, u, CS%diag)
488	12	if (CS%id_v_preale > 0) call post_data(CS%id_v_preale, v, CS%diag)
489	12	if (CS%id_h_preale > 0) call post_data(CS%id_h_preale, h, CS%diag)
490	12	if (CS%id_T_preale > 0) call post_data(CS%id_T_preale, tv%T, CS%diag)
491	12	if (CS%id_S_preale > 0) call post_data(CS%id_S_preale, tv%S, CS%diag)
492	12	if (CS%id_e_preale > 0) then
493		!$omp target update to(h)
494		!$omp target enter data map(alloc: eta_preale)
495	0	call find_eta(h, tv, G, GV, US, eta_preale, dZref=G%Z_ref)
496		!$omp target exit data map(from: eta_preale)
497	0	call post_data(CS%id_e_preale, eta_preale, CS%diag)
498		endif
499
500	12	end subroutine pre_ALE_diagnostics
501
502
503		!> Potentially do some preparatory work, such as convective adjustment, to clean up the model
504		!! state before regridding.
505	13	subroutine pre_ALE_adjustments(G, GV, US, h, tv, Reg, CS, u, v)
506		type(ocean_grid_type), intent(in) :: G !< Ocean grid informations
507		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
508		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
509		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Current 3D grid obtained after the
510		!! last time step [H ~> m or kg m-2]
511		type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamic variable structure
512		type(tracer_registry_type), pointer :: Reg !< Tracer registry structure
513		type(ALE_CS), pointer :: CS !< Regridding parameters and options
514		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
515		optional, intent(inout) :: u !< Zonal velocity field [L T-1 ~> m s-1]
516		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
517		optional, intent(inout) :: v !< Meridional velocity field [L T-1 ~> m s-1]
518
519		integer :: ntr
520
521		! Do column-wise convective adjustment.
522		! Tracers and velocities should probably also undergo consistent adjustments.
523	13	if (CS%do_conv_adj) call convective_adjustment(G, GV, h, tv)
524
525	13	if (CS%use_hybgen_unmix) then
526	0	ntr = 0 ; if (associated(Reg)) ntr = Reg%ntr
527	0	call hybgen_unmix(G, GV, US, CS%hybgen_unmixCS, tv, Reg, ntr, h)
528		endif
529
530	13	end subroutine pre_ALE_adjustments
531
532		!> Takes care of building a new grid. The creation of the new grid can be based on z coordinates,
533		!! target interface densities, sigma coordinates or any arbitrary coordinate system.
534	13	subroutine ALE_regrid( G, GV, US, h, h_new, dzRegrid, tv, CS, frac_shelf_h, PCM_cell)
535		type(ocean_grid_type), intent(in) :: G !< Ocean grid informations
536		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
537		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
538		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses in 3D grid before
539		!! regridding [H ~> m or kg m-2]
540		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: h_new !< Layer thicknesses in 3D grid after
541		!! regridding [H ~> m or kg m-2]
542		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(out) :: dzRegrid !< The change in grid interface positions
543		!! due to regridding, in the same units as
544		!! thicknesses [H ~> m or kg m-2]
545		type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamic variable structure
546		type(ALE_CS), pointer :: CS !< Regridding parameters and options
547		real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: frac_shelf_h !< Fractional ice shelf coverage [nondim]
548		logical, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
549		optional, intent(out) :: PCM_cell !< If true, use PCM remapping in a cell.
550
551		! Local variables
552		logical :: showCallTree
553
554	13	showCallTree = callTree_showQuery()
555
556	13	if (showCallTree) call callTree_enter("ALE_regrid(), MOM_ALE.F90")
557
558		! Build the new grid and store it in h_new. The old grid is retained as h.
559		! Both are needed for the subsequent remapping of variables.
560	8667737	dzRegrid(:,:,:) = 0.0
561		call regridding_main( CS%remapCS, CS%regridCS, G, GV, US, h, tv, h_new, dzRegrid, &
562	13	frac_shelf_h=frac_shelf_h, PCM_cell=PCM_cell)
563
564	13	if (CS%id_dzRegrid>0) then ; if (query_averaging_enabled(CS%diag)) then
565	0	call post_data(CS%id_dzRegrid, dzRegrid, CS%diag, alt_h=h_new)
566		endif ; endif
567
568	13	if (showCallTree) call callTree_leave("ALE_regrid()")
569
570	13	end subroutine ALE_regrid
571
572		!> Regrid/remap stored fields used for offline tracer integrations. These input fields are assumed to have
573		!! the same layer thicknesses at the end of the last offline interval (which should be a Zstar grid). This
574		!! routine builds a grid on the runtime specified vertical coordinate
575	0	subroutine ALE_offline_inputs(CS, G, GV, US, h, tv, Reg, uhtr, vhtr, Kd, debug, OBC)
576		type(ALE_CS), pointer :: CS !< Regridding parameters and options
577		type(ocean_grid_type), intent(in ) :: G !< Ocean grid informations
578		type(verticalGrid_type), intent(in ) :: GV !< Ocean vertical grid structure
579		type(unit_scale_type), intent(in ) :: US !< A dimensional unit scaling type
580		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: h !< Layer thicknesses [H ~> m or kg m-2]
581		type(thermo_var_ptrs), intent(inout) :: tv !< Thermodynamic variable structure
582		type(tracer_registry_type), pointer :: Reg !< Tracer registry structure
583		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(inout) :: uhtr !< Zonal mass fluxes [H L2 ~> m3 or kg]
584		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(inout) :: vhtr !< Meridional mass fluxes [H L2 ~> m3 or kg]
585		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: Kd !< Input diffusivities
586		!! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
587		logical, intent(in ) :: debug !< If true, then turn checksums
588		type(ocean_OBC_type), pointer :: OBC !< Open boundary structure
589		! Local variables
590		integer :: nk, i, j, k, isc, iec, jsc, jec
591	0	real, dimension(SZI_(G), SZJ_(G), SZK_(GV)) :: h_new ! Layer thicknesses after regridding [H ~> m or kg m-2]
592	0	real, dimension(SZI_(G), SZJ_(G), SZK_(GV)+1) :: dzRegrid ! The change in grid interface positions [H ~> m or kg m-2]
593	0	real, dimension(SZK_(GV)) :: h_src ! Source grid thicknesses at velocity points [H ~> m or kg m-2]
594	0	real, dimension(SZK_(GV)) :: h_dest ! Destination grid thicknesses at velocity points [H ~> m or kg m-2]
595	0	real, dimension(SZK_(GV)) :: temp_vec ! Transports on the destination grid [H L2 ~> m3 or kg]
596
597	0	isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; nk = GV%ke
598	0	dzRegrid(:,:,:) = 0.0
599	0	h_new(:,:,:) = 0.0
600
601	0	if (debug) call MOM_tracer_chkinv("Before ALE_offline_inputs", G, GV, h, Reg%Tr, Reg%ntr)
602
603		! Build new grid from the Zstar state onto the requested vertical coordinate. The new grid is stored
604		! in h_new. The old grid is h. Both are needed for the subsequent remapping of variables. Convective
605		! adjustment right now is not used because it is unclear what to do with vanished layers
606	0	call regridding_main( CS%remapCS, CS%regridCS, G, GV, US, h, tv, h_new, dzRegrid)
607	0	if (CS%show_call_tree) call callTree_waypoint("new grid generated (ALE_offline_inputs)")
608
609		! Remap all variables from old grid h onto new grid h_new
610	0	call ALE_remap_tracers(CS, G, GV, h, h_new, Reg, debug=CS%show_call_tree)
611	0	if (allocated(tv%SpV_avg)) tv%valid_SpV_halo = -1 ! Record that SpV_avg is no longer valid.
612	0	if (CS%show_call_tree) call callTree_waypoint("state remapped (ALE_inputs)")
613
614		! Reintegrate mass transports from Zstar to the offline vertical coordinate
615	0	do j=jsc,jec ; do i=G%iscB,G%iecB
616	0	if (G%mask2dCu(i,j)>0.) then
617	0	h_src(:) = 0.5 * (h(i,j,:) + h(i+1,j,:))
618	0	h_dest(:) = 0.5 * (h_new(i,j,:) + h_new(i+1,j,:))
619	0	call reintegrate_column(nk, h_src, uhtr(I,j,:), nk, h_dest, temp_vec)
620	0	uhtr(I,j,:) = temp_vec
621		endif
622		enddo ; enddo
623	0	do j=G%jscB,G%jecB ; do i=isc,iec
624	0	if (G%mask2dCv(i,j)>0.) then
625	0	h_src(:) = 0.5 * (h(i,j,:) + h(i,j+1,:))
626	0	h_dest(:) = 0.5 * (h_new(i,j,:) + h_new(i,j+1,:))
627	0	call reintegrate_column(nk, h_src, vhtr(I,j,:), nk, h_dest, temp_vec)
628	0	vhtr(I,j,:) = temp_vec
629		endif
630		enddo ; enddo
631
632	0	do j=jsc,jec ; do i=isc,iec
633	0	if (G%mask2dT(i,j)>0.) then
634	0	if (check_column_integrals(nk, h_src, nk, h_dest)) then
635	0	call MOM_error(FATAL, "ALE_offline_inputs: Kd interpolation columns do not match")
636		endif
637	0	call interpolate_column(nk, h(i,j,:), Kd(i,j,:), nk, h_new(i,j,:), Kd(i,j,:), .true.)
638		endif
639		enddo ; enddo
640
641	0	call ALE_remap_scalar(CS%remapCS, G, GV, nk, h, tv%T, h_new, tv%T)
642	0	call ALE_remap_scalar(CS%remapCS, G, GV, nk, h, tv%S, h_new, tv%S)
643
644	0	if (debug) call MOM_tracer_chkinv("After ALE_offline_inputs", G, GV, h_new, Reg%Tr, Reg%ntr)
645
646		! Copy over the new layer thicknesses
647	0	do k = 1,nk ; do j = jsc-1,jec+1 ; do i = isc-1,iec+1
648	0	h(i,j,k) = h_new(i,j,k)
649		enddo ; enddo ; enddo
650
651	0	if (allocated(tv%SpV_avg)) tv%valid_SpV_halo = -1 ! Record that SpV_avg is no longer valid.
652
653	0	if (CS%show_call_tree) call callTree_leave("ALE_offline_inputs()")
654	0	end subroutine ALE_offline_inputs
655
656
657		!> For a state-based coordinate, accelerate the process of regridding by
658		!! repeatedly applying the grid calculation algorithm
659	0	subroutine ALE_regrid_accelerated(CS, G, GV, US, h, tv, n_itt, u, v, OBC, Reg, dt, &
660	0	dzRegrid, initial)
661		type(ALE_CS), pointer :: CS !< ALE control structure
662		type(ocean_grid_type), intent(inout) :: G !< Ocean grid
663		type(verticalGrid_type), intent(in) :: GV !< Vertical grid
664		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
665		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
666		intent(inout) :: h !< Original thicknesses [H ~> m or kg m-2]
667		type(thermo_var_ptrs), intent(inout) :: tv !< Thermo vars (T/S/EOS)
668		integer, intent(in) :: n_itt !< Number of times to regrid
669		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
670		intent(inout) :: u !< Zonal velocity [L T-1 ~> m s-1]
671		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
672		intent(inout) :: v !< Meridional velocity [L T-1 ~> m s-1]
673		type(ocean_OBC_type), pointer :: OBC !< Open boundary structure
674		type(tracer_registry_type), &
675		optional, pointer :: Reg !< Tracer registry to remap onto new grid
676		real, optional, intent(in) :: dt !< Model timestep to provide a timescale for regridding [T ~> s]
677		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
678		optional, intent(inout) :: dzRegrid !< Final change in interface positions [H ~> m or kg m-2]
679		logical, optional, intent(in) :: initial !< Whether we're being called from an initialization
680		!! routine (and expect diagnostics to work)
681
682		! Local variables
683		integer :: i, j, itt, nz
684	0	type(thermo_var_ptrs) :: tv_local ! local/intermediate temp/salt
685		type(group_pass_type) :: pass_T_S_h ! group pass if the coordinate has a stencil
686	0	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_loc ! A working copy of layer thicknesses [H ~> m or kg m-2]
687	0	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_orig ! The original layer thicknesses [H ~> m or kg m-2]
688	0	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), target :: T ! local temporary temperatures [C ~> degC]
689	0	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), target :: S ! local temporary salinities [S ~> ppt]
690	0	real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: h_old_u ! Source grid thickness at zonal
691		! velocity points [H ~> m or kg m-2]
692	0	real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: h_old_v ! Source grid thickness at meridional
693		! velocity points [H ~> m or kg m-2]
694	0	real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)) :: h_new_u ! Destination grid thickness at zonal
695		! velocity points [H ~> m or kg m-2]
696	0	real, dimension(SZI_(G),SZJB_(G),SZK_(GV)) :: h_new_v ! Destination grid thickness at meridional
697		! velocity points [H ~> m or kg m-2]
698
699		! we have to keep track of the total dzInterface if for some reason
700		! we're using the old remapping algorithm for u/v
701	0	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: dzInterface ! Interface height changes within
702		! an iteration [H ~> m or kg m-2]
703	0	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: dzIntTotal ! Cumulative interface position changes [H ~> m or kg m-2]
704
705	0	nz = GV%ke
706
707		! initial total interface displacement due to successive regridding
708	0	if (CS%remap_uv_using_old_alg) &
709	0	dzIntTotal(:,:,:) = 0.
710
711	0	call create_group_pass(pass_T_S_h, T, G%domain)
712	0	call create_group_pass(pass_T_S_h, S, G%domain)
713	0	call create_group_pass(pass_T_S_h, h_loc, G%domain)
714
715		! copy original temp/salt and set our local tv_pointers to them
716	0	tv_local = tv
717	0	T(:,:,:) = tv%T(:,:,:)
718	0	S(:,:,:) = tv%S(:,:,:)
719	0	tv_local%T => T
720	0	tv_local%S => S
721
722		! get local copy of thickness and save original state for remapping
723	0	h_loc(:,:,:) = h(:,:,:)
724	0	h_orig(:,:,:) = h(:,:,:)
725
726		! Apply timescale to regridding (for e.g. filtered_grid_motion)
727	0	if (present(dt)) &
728	0	call ALE_update_regrid_weights(dt, CS)
729
730	0	do itt = 1, n_itt
731
732	0	call do_group_pass(pass_T_S_h, G%domain)
733
734		! generate new grid
735	0	if (CS%do_conv_adj) call convective_adjustment(G, GV, h_loc, tv_local)
736
737		! Update the layer specific volumes if necessary
738	0	if (allocated(tv_local%SpV_avg)) call calc_derived_thermo(tv_local, h, G, GV, US, halo=1)
739
740	0	call regridding_main(CS%remapCS, CS%regridCS, G, GV, US, h_loc, tv_local, h, dzInterface)
741	0	if (CS%remap_uv_using_old_alg) &
742	0	dzIntTotal(:,:,:) = dzIntTotal(:,:,:) + dzInterface(:,:,:)
743
744		! remap from original grid onto new grid
745	0	do j = G%jsc-1,G%jec+1 ; do i = G%isc-1,G%iec+1
746		call remapping_core_h(CS%remapCS, nz, h_orig(i,j,:), tv%S(i,j,:), nz, h(i,j,:), &
747	0	tv_local%S(i,j,:))
748		call remapping_core_h(CS%remapCS, nz, h_orig(i,j,:), tv%T(i,j,:), nz, h(i,j,:), &
749	0	tv_local%T(i,j,:))
750		enddo ; enddo
751
752		! starting grid for next iteration
753	0	h_loc(:,:,:) = h(:,:,:)
754		enddo
755
756		! remap all state variables (including those that weren't needed for regridding)
757	0	call ALE_remap_tracers(CS, G, GV, h_orig, h, Reg)
758
759	0	call ALE_remap_set_h_vel(CS, G, GV, h_orig, h_old_u, h_old_v, OBC)
760	0	if (CS%remap_uv_using_old_alg) then
761	0	call ALE_remap_set_h_vel_via_dz(CS, G, GV, h, h_new_u, h_new_v, OBC, h_orig, dzIntTotal)
762		else
763	0	call ALE_remap_set_h_vel(CS, G, GV, h, h_new_u, h_new_v, OBC)
764		endif
765
766	0	call ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u, v)
767
768		! save total dzregrid for diags if needed?
769	0	if (present(dzRegrid)) dzRegrid(:,:,:) = dzIntTotal(:,:,:)
770
771	0	if (allocated(tv%SpV_avg)) tv%valid_SpV_halo = -1 ! Record that SpV_avg is no longer valid.
772
773	0	end subroutine ALE_regrid_accelerated
774
775		!> This routine takes care of remapping all tracer variables between the old and the
776		!! new grids. This routine is called during initialization of the model at time=0, to
777		!! remap initial conditions to the model grid. It is also called during a
778		!! time step to update the state.
779	13	subroutine ALE_remap_tracers(CS, G, GV, h_old, h_new, Reg, debug, dt, PCM_cell)
780		type(ALE_CS), intent(in) :: CS !< ALE control structure
781		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
782		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
783		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_old !< Thickness of source grid
784		!! [H ~> m or kg m-2]
785		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_new !< Thickness of destination grid
786		!! [H ~> m or kg m-2]
787		type(tracer_registry_type), pointer :: Reg !< Tracer registry structure
788		logical, optional, intent(in) :: debug !< If true, show the call tree
789		real, optional, intent(in) :: dt !< time step for diagnostics [T ~> s]
790		logical, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
791		optional, intent(in) :: PCM_cell !< Use PCM remapping in cells where true
792
793		! Local variables
794	26	real :: tr_column(GV%ke) ! A column of updated tracer concentrations [CU ~> Conc]
795	26	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: work_conc ! The rate of change of concentrations [Conc T-1 ~> Conc s-1]
796	26	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: work_cont ! The rate of change of cell-integrated tracer
797		! content [Conc H T-1 ~> Conc m s-1 or Conc kg m-2 s-1] or
798		! cell thickness [H T-1 ~> m s-1 or kg m-2 s-1]
799	26	real, dimension(SZI_(G),SZJ_(G)) :: work_2d ! The rate of change of column-integrated tracer
800		! content [Conc H T-1 ~> Conc m s-1 or Conc kg m-2 s-1]
801	26	logical :: PCM(GV%ke) ! If true, do PCM remapping from a cell.
802		real :: Idt ! The inverse of the timestep [T-1 ~> s-1]
803	26	real :: h1(GV%ke) ! A column of source grid layer thicknesses [H ~> m or kg m-2]
804	26	real :: h2(GV%ke) ! A column of target grid layer thicknesses [H ~> m or kg m-2]
805		logical :: show_call_tree
806		type(tracer_type), pointer :: Tr => NULL()
807		integer :: i, j, k, m, nz, ntr
808
809	13	show_call_tree = .false.
810	13	if (present(debug)) show_call_tree = debug
811
812	13	if (show_call_tree) call callTree_enter("ALE_remap_tracers(), MOM_ALE.F90")
813
814	13	nz = GV%ke
815
816	13	ntr = 0 ; if (associated(Reg)) ntr = Reg%ntr
817
818	13	if (present(dt)) then
819	12	Idt = 1.0/dt
820	7895712	work_conc(:,:,:) = 0.0
821	7895712	work_cont(:,:,:) = 0.0
822		endif
823
824		! Remap all registered tracers, including temperature and salinity.
825	13	if (ntr>0) then
826	13	if (show_call_tree) call callTree_waypoint("remapping tracers (ALE_remap_tracers)")
827		!$OMP parallel do default(shared) private(h1,h2,tr_column,Tr,PCM,work_conc,work_cont,work_2d)
828	52	do m=1,ntr ! For each tracer
829	39	Tr => Reg%Tr(m)
830	283179	do j = G%jsc,G%jec ; do i = G%isc,G%iec ; if (G%mask2dT(i,j)>0.) then
831		! Build the start and final grids
832	14861496	h1(:) = h_old(i,j,:)
833	14861496	h2(:) = h_new(i,j,:)
834	195546	if (present(PCM_cell)) then
835	14861496	PCM(:) = PCM_cell(i,j,:)
836	195546	call remapping_core_h(CS%remapCS, nz, h1, Tr%t(i,j,:), nz, h2, tr_column, PCM_cell=PCM)
837		else
838	0	call remapping_core_h(CS%remapCS, nz, h1, Tr%t(i,j,:), nz, h2, tr_column)
839		endif
840
841		! Possibly underflow any very tiny tracer concentrations to 0. Note that this is not conservative!
842	195546	if (Tr%conc_underflow > 0.0) then ; do k=1,GV%ke
843	0	if (abs(tr_column(k)) < Tr%conc_underflow) tr_column(k) = 0.0
844		enddo ; endif
845
846		! Intermediate steps for tendency of tracer concentration and tracer content.
847	195546	if (present(dt)) then
848	180504	if (Tr%id_remap_conc > 0) then
849	0	do k=1,GV%ke
850	0	work_conc(i,j,k) = (tr_column(k) - Tr%t(i,j,k)) * Idt
851		enddo
852		endif
853	180504	if (Tr%id_remap_cont > 0 .or. Tr%id_remap_cont_2d > 0) then
854	0	do k=1,GV%ke
855	0	work_cont(i,j,k) = (tr_column(k)h2(k) - Tr%t(i,j,k)h1(k)) * Idt
856		enddo
857		endif
858		endif
859
860		! update tracer concentration
861	14861496	Tr%t(i,j,:) = tr_column(:)
862		endif ; enddo ; enddo
863
864		! tendency diagnostics.
865	52	if (present(dt)) then
866	36	if (Tr%id_remap_conc > 0) then
867	0	call post_data(Tr%id_remap_conc, work_conc, CS%diag)
868		endif
869	36	if (Tr%id_remap_cont > 0) then
870	0	call post_data(Tr%id_remap_cont, work_cont, CS%diag)
871		endif
872
873	36	if (Tr%id_remap_cont_2d > 0) then
874	0	do j = G%jsc,G%jec ; do i = G%isc,G%iec
875	0	work_2d(i,j) = 0.0
876	0	do k = 1,GV%ke
877	0	work_2d(i,j) = work_2d(i,j) + work_cont(i,j,k)
878		enddo
879		enddo ; enddo
880	0	call post_data(Tr%id_remap_cont_2d, work_2d, CS%diag)
881		endif
882		endif
883		enddo ! m=1,ntr
884
885		endif ! endif for ntr > 0
886
887
888	13	if (CS%id_vert_remap_h > 0) call post_data(CS%id_vert_remap_h, h_old, CS%diag)
889	13	if ((CS%id_vert_remap_h_tendency > 0) .and. present(dt)) then
890	0	do k = 1, nz ; do j = G%jsc,G%jec ; do i = G%isc,G%iec
891	0	work_cont(i,j,k) = (h_new(i,j,k) - h_old(i,j,k))*Idt
892		enddo ; enddo ; enddo
893	0	call post_data(CS%id_vert_remap_h_tendency, work_cont, CS%diag)
894		endif
895
896	13	if (show_call_tree) call callTree_leave("ALE_remap_tracers(), MOM_ALE.F90")
897
898	13	end subroutine ALE_remap_tracers
899
900		!> This routine sets the thicknesses at velocity points used for vertical remapping.
901	26	subroutine ALE_remap_set_h_vel(CS, G, GV, h_new, h_u, h_v, OBC, debug)
902		type(ALE_CS), intent(in) :: CS !< ALE control structure
903		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
904		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
905		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_new !< Thickness at tracer points of the
906		!! grid being interpolated to velocity
907		!! points [H ~> m or kg m-2]
908		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
909		intent(inout) :: h_u !< Grid thickness at zonal velocity
910		!! points [H ~> m or kg m-2]
911		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
912		intent(inout) :: h_v !< Grid thickness at meridional velocity
913		!! points [H ~> m or kg m-2]
914		type(ocean_OBC_type), pointer :: OBC !< Open boundary structure
915		logical, optional, intent(in) :: debug !< If true, show the call tree
916
917		! Local variables
918		logical :: show_call_tree
919		integer :: i, j, k
920
921	26	show_call_tree = .false.
922	26	if (present(debug)) show_call_tree = debug
923	26	if (show_call_tree) call callTree_enter("ALE_remap_set_h_vel()")
924
925		! Build the u- and v-velocity grid thicknesses for remapping.
926
927		!$OMP parallel do default(shared)
928	14275976	do k=1,GV%ke ; do j=G%jsc,G%jec ; do I=G%IscB,G%IecB ; if (G%mask2dCu(I,j)>0.) then
929	9648600	h_u(I,j,k) = 0.5*(h_new(i,j,k) + h_new(i+1,j,k))
930		endif ; enddo ; enddo ; enddo
931		!$OMP parallel do default(shared)
932	14394926	do k=1,GV%ke ; do J=G%JscB,G%JecB ; do i=G%isc,G%iec ; if (G%mask2dCv(i,J)>0.) then
933	9543300	h_v(i,J,k) = 0.5*(h_new(i,j,k) + h_new(i,j+1,k))
934		endif ; enddo ; enddo ; enddo
935
936		! Mask out blocked portions of velocity cells.
937	26	if (CS%partial_cell_vel_remap) call ALE_remap_set_h_vel_partial(CS, G, GV, h_new, h_u, h_v)
938
939		! Take open boundary conditions into account.
940	26	if (associated(OBC)) call ALE_remap_set_h_vel_OBC(G, GV, h_new, h_u, h_v, OBC)
941
942	26	if (show_call_tree) call callTree_leave("ALE_remap_set_h_vel()")
943
944	26	end subroutine ALE_remap_set_h_vel
945
946		!> This routine sets the thicknesses at velocity points used for vertical remapping using a
947		!! combination of the old grid and interface movements.
948	0	subroutine ALE_remap_set_h_vel_via_dz(CS, G, GV, h_new, h_u, h_v, OBC, h_old, dzInterface, debug)
949		type(ALE_CS), intent(in) :: CS !< ALE control structure
950		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
951		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
952		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_new !< Thickness at tracer points of the
953		!! grid being interpolated to velocity
954		!! points [H ~> m or kg m-2]
955		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
956		intent(inout) :: h_u !< Grid thickness at zonal velocity
957		!! points [H ~> m or kg m-2]
958		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
959		intent(inout) :: h_v !< Grid thickness at meridional velocity
960		!! points [H ~> m or kg m-2]
961		type(ocean_OBC_type), pointer :: OBC !< Open boundary structure
962		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
963		intent(in) :: h_old !< Thickness of source grid when generating
964		!! the destination grid via the old
965		!! algorithm [H ~> m or kg m-2]
966		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
967		intent(in) :: dzInterface !< Change in interface position
968		!! [H ~> m or kg m-2]
969		logical, optional, intent(in) :: debug !< If true, show the call tree
970
971		! Local variables
972		logical :: show_call_tree
973		integer :: i, j, k
974
975	0	show_call_tree = .false.
976	0	if (present(debug)) show_call_tree = debug
977	0	if (show_call_tree) call callTree_enter("ALE_remap_set_h_vel()")
978
979		! Build the u- and v-velocity grid thicknesses for remapping using the old grid and interface movement.
980
981		!$OMP parallel do default(shared)
982	0	do k=1,GV%ke ; do j=G%jsc,G%jec ; do I=G%IscB,G%IecB ; if (G%mask2dCu(I,j)>0.) then
983		h_u(I,j,k) = max( 0., 0.5*(h_old(i,j,k) + h_old(i+1,j,k)) + &
984		0.5 * (( dzInterface(i,j,k) + dzInterface(i+1,j,k) ) - &
985	0	( dzInterface(i,j,k+1) + dzInterface(i+1,j,k+1) )) )
986		endif ; enddo ; enddo ; enddo
987
988		!$OMP parallel do default(shared)
989	0	do k=1,GV%ke ; do J=G%JscB,G%JecB ; do i=G%isc,G%iec ; if (G%mask2dCv(i,J)>0.) then
990		h_v(i,J,k) = max( 0., 0.5*(h_old(i,j,k) + h_old(i,j+1,k)) + &
991		0.5 * (( dzInterface(i,j,k) + dzInterface(i,j+1,k) ) - &
992	0	( dzInterface(i,j,k+1) + dzInterface(i,j+1,k+1) )) )
993		endif ; enddo ; enddo ; enddo
994
995		! Mask out blocked portions of velocity cells.
996	0	if (CS%partial_cell_vel_remap) call ALE_remap_set_h_vel_partial(CS, G, GV, h_old, h_u, h_v)
997
998		! Take open boundary conditions into account.
999	0	if (associated(OBC)) call ALE_remap_set_h_vel_OBC(G, GV, h_new, h_u, h_v, OBC)
1000
1001	0	if (show_call_tree) call callTree_leave("ALE_remap_set_h_vel()")
1002
1003	0	end subroutine ALE_remap_set_h_vel_via_dz
1004
1005		!> Mask out the thicknesses at velocity points where they are below the minimum depth
1006		!! at adjacent tracer points
1007	0	subroutine ALE_remap_set_h_vel_partial(CS, G, GV, h_mask, h_u, h_v)
1008		type(ALE_CS), intent(in) :: CS !< ALE control structure
1009		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1010		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
1011		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_mask !< Thickness at tracer points
1012		!! used to apply the partial
1013		!! cell masking [H ~> m or kg m-2]
1014		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1015		intent(inout) :: h_u !< Grid thickness at zonal velocity
1016		!! points [H ~> m or kg m-2]
1017		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1018		intent(inout) :: h_v !< Grid thickness at meridional velocity
1019		!! points [H ~> m or kg m-2]
1020		! Local variables
1021	0	real, dimension(SZI_(G),SZJ_(G)) :: h_tot ! The vertically summed thicknesses [H ~> m or kg m-2]
1022		real :: h_mask_vel ! A depth below which the thicknesses at a velocity point are masked out [H ~> m or kg m-2]
1023		integer :: i, j, k
1024
1025	0	h_tot(:,:) = 0.0
1026	0	do k=1,GV%ke ; do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1
1027	0	h_tot(i,j) = h_tot(i,j) + h_mask(i,j,k)
1028		enddo ; enddo ; enddo
1029
1030		!$OMP parallel do default(shared) private(h_mask_vel)
1031	0	do j=G%jsc,G%jec ; do I=G%IscB,G%IecB ; if (G%mask2dCu(I,j)>0.) then
1032	0	h_mask_vel = min(h_tot(i,j), h_tot(i+1,j))
1033	0	call apply_partial_cell_mask(h_u(I,j,:), h_mask_vel)
1034		endif ; enddo ; enddo
1035
1036		!$OMP parallel do default(shared) private(h_mask_vel)
1037	0	do J=G%JscB,G%JecB ; do i=G%isc,G%iec ; if (G%mask2dCv(i,J)>0.) then
1038	0	h_mask_vel = min(h_tot(i,j), h_tot(i,j+1))
1039	0	call apply_partial_cell_mask(h_v(i,J,:), h_mask_vel)
1040		endif ; enddo ; enddo
1041
1042	0	end subroutine ALE_remap_set_h_vel_partial
1043
1044		! Reset thicknesses at velocity points on open boundary condition segments
1045	0	subroutine ALE_remap_set_h_vel_OBC(G, GV, h_new, h_u, h_v, OBC)
1046		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1047		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
1048		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_new !< Thickness at tracer points of the
1049		!! grid being interpolated to velocity
1050		!! points [H ~> m or kg m-2]
1051		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1052		intent(inout) :: h_u !< Grid thickness at zonal velocity
1053		!! points [H ~> m or kg m-2]
1054		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1055		intent(inout) :: h_v !< Grid thickness at meridional velocity
1056		!! points [H ~> m or kg m-2]
1057		type(ocean_OBC_type), pointer :: OBC !< Open boundary structure
1058
1059		! Local variables
1060		integer :: i, j, k, nz, is_OBC, ie_OBC, js_OBC, je_OBC
1061
1062	0	if (.not.associated(OBC)) return
1063
1064	0	nz = GV%ke
1065
1066		! Take open boundary conditions into account.
1067	0	if (OBC%u_E_OBCs_on_PE) then
1068	0	js_OBC = max(G%jsc, OBC%js_u_E_obc) ; je_OBC = min(G%jec, OBC%je_u_E_obc)
1069	0	Is_OBC = max(G%IscB, OBC%Is_u_E_obc) ; Ie_OBC = min(G%IecB, OBC%Ie_u_E_obc)
1070		!$OMP parallel do default(shared)
1071	0	do j=js_OBC,je_OBC ; do I=Is_OBC,Ie_OBC ; if (OBC%segnum_u(I,j) > 0) then ! OBC_DIRECTION_E
1072	0	do k=1,nz ; h_u(I,j,k) = h_new(i,j,k) ; enddo
1073		endif ; enddo ; enddo
1074		endif
1075	0	if (OBC%u_W_OBCs_on_PE) then
1076	0	js_OBC = max(G%jsc, OBC%js_u_W_obc) ; je_OBC = min(G%jec, OBC%je_u_W_obc)
1077	0	Is_OBC = max(G%IscB, OBC%Is_u_W_obc) ; Ie_OBC = min(G%IecB, OBC%Ie_u_W_obc)
1078		!$OMP parallel do default(shared)
1079	0	do j=js_OBC,je_OBC ; do I=Is_OBC,Ie_OBC ; if (OBC%segnum_u(I,j) < 0) then ! OBC_DIRECTION_W
1080	0	do k=1,nz ; h_u(I,j,k) = h_new(i+1,j,k) ; enddo
1081		endif ; enddo ; enddo
1082		endif
1083
1084	0	if (OBC%v_N_OBCs_on_PE) then
1085	0	Js_OBC = max(G%JscB, OBC%Js_v_N_obc) ; Je_OBC = min(G%JecB, OBC%Je_v_N_obc)
1086	0	is_OBC = max(G%isc, OBC%is_v_N_obc) ; ie_OBC = min(G%iec, OBC%ie_v_N_obc)
1087		!$OMP parallel do default(shared)
1088	0	do J=Js_OBC,Je_OBC ; do i=is_OBC,ie_OBC ; if (OBC%segnum_v(i,J) > 0) then ! OBC_DIRECTION_N
1089	0	do k=1,nz ; h_v(i,J,k) = h_new(i,j,k) ; enddo
1090		endif ; enddo ; enddo
1091		endif
1092	0	if (OBC%v_S_OBCs_on_PE) then
1093	0	Js_OBC = max(G%JscB, OBC%Js_v_S_obc) ; Je_OBC = min(G%JecB, OBC%Je_v_S_obc)
1094	0	is_OBC = max(G%isc, OBC%is_v_S_obc) ; ie_OBC = min(G%iec, OBC%ie_v_S_obc)
1095		!$OMP parallel do default(shared)
1096	0	do J=Js_OBC,Je_OBC ; do i=is_OBC,ie_OBC ; if (OBC%segnum_v(i,J) < 0) then ! OBC_DIRECTION_S
1097	0	do k=1,nz ; h_v(i,J,k) = h_new(i,j+1,k) ; enddo
1098		endif ; enddo ; enddo
1099		endif
1100
1101		end subroutine ALE_remap_set_h_vel_OBC
1102
1103		!> This routine remaps velocity components between the old and the new grids,
1104		!! with thicknesses at velocity points taken to be arithmetic averages of tracer thicknesses.
1105		!! This routine may be called during initialization of the model at time=0, to
1106		!! remap initial conditions to the model grid. It is also called during a
1107		!! time step to update the state.
1108	13	subroutine ALE_remap_velocities(CS, G, GV, h_old_u, h_old_v, h_new_u, h_new_v, u, v, debug, &
1109		dt, allow_preserve_variance)
1110		type(ALE_CS), intent(in) :: CS !< ALE control structure
1111		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1112		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
1113		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1114		intent(in) :: h_old_u !< Source grid thickness at zonal
1115		!! velocity points [H ~> m or kg m-2]
1116		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1117		intent(in) :: h_old_v !< Source grid thickness at meridional
1118		!! velocity points [H ~> m or kg m-2]
1119		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1120		intent(in) :: h_new_u !< Destination grid thickness at zonal
1121		!! velocity points [H ~> m or kg m-2]
1122		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1123		intent(in) :: h_new_v !< Destination grid thickness at meridional
1124		!! velocity points [H ~> m or kg m-2]
1125		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), &
1126		intent(inout) :: u !< Zonal velocity [L T-1 ~> m s-1]
1127		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), &
1128		intent(inout) :: v !< Meridional velocity [L T-1 ~> m s-1]
1129		logical, optional, intent(in) :: debug !< If true, show the call tree
1130		real, optional, intent(in) :: dt !< time step for diagnostics [T ~> s]
1131		logical, optional, intent(in) :: allow_preserve_variance !< If true, enables ke-conserving
1132		!! correction
1133
1134		! Local variables
1135		real :: h_mask_vel ! A depth below which the thicknesses at a velocity point are masked out [H ~> m or kg m-2]
1136	26	real :: u_src(GV%ke) ! A column of u-velocities on the source grid [L T-1 ~> m s-1]
1137	26	real :: u_tgt(GV%ke) ! A column of u-velocities on the target grid [L T-1 ~> m s-1]
1138	26	real :: v_src(GV%ke) ! A column of v-velocities on the source grid [L T-1 ~> m s-1]
1139	26	real :: v_tgt(GV%ke) ! A column of v-velocities on the target grid [L T-1 ~> m s-1]
1140	26	real :: h1(GV%ke) ! A column of source grid layer thicknesses [H ~> m or kg m-2]
1141	26	real :: h2(GV%ke) ! A column of target grid layer thicknesses [H ~> m or kg m-2]
1142		real :: rescale_coef ! Factor that scales the baroclinic velocity to conserve ke [nondim]
1143		real :: u_bt, v_bt ! Depth-averaged velocity components [L T-1 ~> m s-1]
1144		real :: ke_c_src, ke_c_tgt ! \int [u_c or v_c]^2 dz on src and tgt grids [H L2 T-2 ~> m3 s-2]
1145	26	real, dimension(SZIB_(G),SZJ_(G)) :: du2h_tot ! The rate of change of vertically integrated
1146		! 0.5 * rho0 * u**2 [R Z L2 T-3 ~> W m-2]
1147	26	real, dimension(SZI_(G),SZJB_(G)) :: dv2h_tot ! The rate of change of vertically integrated
1148		! 0.5 * rho0 * v**2 [R Z L2 T-3 ~> W m-2]
1149		real :: u2h_tot, v2h_tot ! The vertically integrated u2 and v2 [H L2 T-2 ~> m3 s-2 or kg s-2]
1150		real :: I_dt ! 1 / dt [T-1 ~> s-1]
1151		logical :: variance_option ! Contains the value of allow_preserve_variance when present, else false
1152		logical :: show_call_tree
1153		integer :: i, j, k, nz
1154
1155	13	show_call_tree = .false.
1156	13	if (present(debug)) show_call_tree = debug
1157	13	if (show_call_tree) call callTree_enter("ALE_remap_velocities()")
1158
1159		! Setup related to KE conservation
1160	13	variance_option = .false.
1161	13	if (present(allow_preserve_variance)) variance_option=allow_preserve_variance
1162	13	if (present(dt)) I_dt = 1.0 / dt
1163
1164	13	if (CS%id_remap_delta_integ_u2>0) du2h_tot(:,:) = 0.
1165	13	if (CS%id_remap_delta_integ_v2>0) dv2h_tot(:,:) = 0.
1166
1167	13	if (((CS%id_remap_delta_integ_u2>0) .or. (CS%id_remap_delta_integ_v2>0)) .and. .not.present(dt))&
1168	0	call MOM_error(FATAL, "ALE KE diagnostics requires passing dt into ALE_remap_velocities")
1169
1170	13	nz = GV%ke
1171
1172		! --- Remap u profiles from the source vertical grid onto the new target grid.
1173
1174		!$OMP parallel do default(shared) private(h1,h2,u_src,h_mask_vel,u_tgt, &
1175		!$OMP u_bt,ke_c_src,ke_c_tgt,rescale_coef, &
1176		!$OMP u2h_tot,v2h_tot)
1177	95173	do j=G%jsc,G%jec ; do I=G%IscB,G%IecB ; if (G%mask2dCu(I,j)>0.) then
1178		! Make a 1-d copy of the start and final grids and the source velocity
1179	4888624	do k=1,nz
1180	4824300	h1(k) = h_old_u(I,j,k)
1181	4824300	h2(k) = h_new_u(I,j,k)
1182	4888624	u_src(k) = u(I,j,k)
1183		enddo
1184
1185	64324	if (CS%id_remap_delta_integ_u2>0) then
1186	0	u2h_tot = 0.
1187	0	do k=1,nz
1188	0	u2h_tot = u2h_tot - h1(k) * (u_src(k)**2)
1189		enddo
1190		endif
1191
1192	64324	call remapping_core_h(CS%vel_remapCS, nz, h1, u_src, nz, h2, u_tgt)
1193
1194	64324	if (variance_option .and. CS%conserve_ke) then
1195		! Conserve ke_u by correcting baroclinic component.
1196		! Assumes total depth doesn't change during remap, and
1197		! that \int u(z) dz doesn't change during remap.
1198		! First get barotropic component
1199	0	u_bt = 0.0
1200	0	do k=1,nz
1201	0	u_bt = u_bt + h2(k) * u_tgt(k) ! Dimensions [H L T-1 ~> m2 s-1 or kg m-1 s-1]
1202		enddo
1203	0	u_bt = u_bt / (sum(h2(1:nz)) + GV%H_subroundoff) ! Dimensions return to [L T-1 ~> m s-1]
1204		! Next get baroclinic ke = \int (u-u_bt)^2 from source and target
1205	0	ke_c_src = 0.0
1206	0	ke_c_tgt = 0.0
1207	0	do k=1,nz
1208	0	ke_c_src = ke_c_src + h1(k) * (u_src(k) - u_bt)**2
1209	0	ke_c_tgt = ke_c_tgt + h2(k) * (u_tgt(k) - u_bt)**2
1210		enddo
1211		! Next rescale baroclinic component on target grid to conserve ke
1212		! The values 1.5625 = 1.25**2 and 1.25 below mean that the KE-conserving
1213		! correction cannot amplify the baroclinic part of velocity by more
1214		! than 25%. This threshold is somewhat arbitrary. It was added to
1215		! prevent unstable behavior when the amplification factor is large.
1216	0	if (ke_c_src < 1.5625 * ke_c_tgt) then
1217	0	rescale_coef = sqrt(ke_c_src / ke_c_tgt)
1218		else
1219	0	rescale_coef = 1.25
1220		endif
1221	0	do k=1,nz
1222	0	u_tgt(k) = u_bt + rescale_coef * (u_tgt(k) - u_bt)
1223		enddo
1224		endif
1225
1226	64324	if (CS%id_remap_delta_integ_u2>0) then
1227	0	do k=1,nz
1228	0	u2h_tot = u2h_tot + h2(k) * (u_tgt(k)**2)
1229		enddo
1230	0	du2h_tot(I,j) = u2h_tot * I_dt
1231		endif
1232
1233	64324	if ((CS%BBL_h_vel_mask > 0.0) .and. (CS%h_vel_mask > 0.0)) &
1234	0	call mask_near_bottom_vel(u_tgt, h2, CS%BBL_h_vel_mask, CS%h_vel_mask, nz)
1235
1236		! Copy the column of new velocities back to the 3-d array
1237	4888624	do k=1,nz
1238	4888624	u(I,j,k) = u_tgt(k)
1239		enddo !k
1240		endif ; enddo ; enddo
1241
1242	13	if (CS%id_remap_delta_integ_u2>0) call post_data(CS%id_remap_delta_integ_u2, du2h_tot, CS%diag)
1243
1244	13	if (show_call_tree) call callTree_waypoint("u remapped (ALE_remap_velocities)")
1245
1246
1247		! --- Remap v profiles from the source vertical grid onto the new target grid.
1248
1249		!$OMP parallel do default(shared) private(h1,h2,v_src,h_mask_vel,v_tgt, &
1250		!$OMP v_bt,ke_c_src,ke_c_tgt,rescale_coef, &
1251		!$OMP u2h_tot,v2h_tot)
1252	95966	do J=G%JscB,G%JecB ; do i=G%isc,G%iec ; if (G%mask2dCv(i,J)>0.) then
1253
1254	4835272	do k=1,nz
1255	4771650	h1(k) = h_old_v(i,J,k)
1256	4771650	h2(k) = h_new_v(i,J,k)
1257	4835272	v_src(k) = v(i,J,k)
1258		enddo
1259
1260	63622	if (CS%id_remap_delta_integ_v2>0) then
1261	0	v2h_tot = 0.
1262	0	do k=1,nz
1263	0	v2h_tot = v2h_tot - h1(k) * (v_src(k)**2)
1264		enddo
1265		endif
1266
1267	63622	call remapping_core_h(CS%vel_remapCS, nz, h1, v_src, nz, h2, v_tgt)
1268
1269	63622	if (variance_option .and. CS%conserve_ke) then
1270		! Conserve ke_v by correcting baroclinic component.
1271		! Assumes total depth doesn't change during remap, and
1272		! that \int v(z) dz doesn't change during remap.
1273		! First get barotropic component
1274	0	v_bt = 0.0
1275	0	do k=1,nz
1276	0	v_bt = v_bt + h2(k) * v_tgt(k) ! Dimensions [H L T-1 ~> m2 s-1 or kg m-1 s-1]
1277		enddo
1278	0	v_bt = v_bt / (sum(h2(1:nz)) + GV%H_subroundoff) ! Dimensions return to [L T-1 ~> m s-1]
1279		! Next get baroclinic ke = \int (u-u_bt)^2 from source and target
1280	0	ke_c_src = 0.0
1281	0	ke_c_tgt = 0.0
1282	0	do k=1,nz
1283	0	ke_c_src = ke_c_src + h1(k) * (v_src(k) - v_bt)**2
1284	0	ke_c_tgt = ke_c_tgt + h2(k) * (v_tgt(k) - v_bt)**2
1285		enddo
1286		! Next rescale baroclinic component on target grid to conserve ke
1287	0	if (ke_c_src < 1.5625 * ke_c_tgt) then
1288	0	rescale_coef = sqrt(ke_c_src / ke_c_tgt)
1289		else
1290	0	rescale_coef = 1.25
1291		endif
1292	0	do k=1,nz
1293	0	v_tgt(k) = v_bt + rescale_coef * (v_tgt(k) - v_bt)
1294		enddo
1295		endif
1296
1297	63622	if (CS%id_remap_delta_integ_v2>0) then
1298	0	do k=1,nz
1299	0	v2h_tot = v2h_tot + h2(k) * (v_tgt(k)**2)
1300		enddo
1301	0	dv2h_tot(I,j) = v2h_tot * I_dt
1302		endif
1303
1304	63622	if ((CS%BBL_h_vel_mask > 0.0) .and. (CS%h_vel_mask > 0.0)) then
1305	0	call mask_near_bottom_vel(v_tgt, h2, CS%BBL_h_vel_mask, CS%h_vel_mask, nz)
1306		endif
1307
1308		! Copy the column of new velocities back to the 3-d array
1309	4835272	do k=1,nz
1310	4835272	v(i,J,k) = v_tgt(k)
1311		enddo !k
1312		endif ; enddo ; enddo
1313
1314	13	if (CS%id_remap_delta_integ_v2>0) call post_data(CS%id_remap_delta_integ_v2, dv2h_tot, CS%diag)
1315
1316	13	if (show_call_tree) call callTree_waypoint("v remapped (ALE_remap_velocities)")
1317	13	if (show_call_tree) call callTree_leave("ALE_remap_velocities()")
1318
1319	13	end subroutine ALE_remap_velocities
1320
1321		!> Interpolate to find an updated array of values at interfaces after remapping.
1322	0	subroutine ALE_remap_interface_vals(CS, G, GV, h_old, h_new, int_val)
1323		type(ALE_CS), intent(in) :: CS !< ALE control structure
1324		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1325		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
1326		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_old !< Thickness of source grid
1327		!! [H ~> m or kg m-2]
1328		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_new !< Thickness of destination grid
1329		!! [H ~> m or kg m-2]
1330		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
1331		intent(inout) :: int_val !< The interface values to interpolate [A]
1332
1333	0	real :: val_src(GV%ke+1) ! A column of interface values on the source grid [A]
1334	0	real :: val_tgt(GV%ke+1) ! A column of interface values on the target grid [A]
1335	0	real :: h_src(GV%ke) ! A column of source grid layer thicknesses [H ~> m or kg m-2]
1336	0	real :: h_tgt(GV%ke) ! A column of target grid layer thicknesses [H ~> m or kg m-2]
1337		integer :: i, j, k, nz
1338
1339	0	nz = GV%ke
1340
1341	0	do j=G%jsc,G%jec ; do i=G%isc,G%iec ; if (G%mask2dT(i,j)>0.) then
1342	0	do k=1,nz
1343	0	h_src(k) = h_old(i,j,k)
1344	0	h_tgt(k) = h_new(i,j,k)
1345		enddo
1346
1347	0	do K=1,nz+1
1348	0	val_src(K) = int_val(i,j,K)
1349		enddo
1350
1351	0	call interpolate_column(nz, h_src, val_src, nz, h_tgt, val_tgt, .false.)
1352
1353	0	do K=1,nz+1
1354	0	int_val(i,j,K) = val_tgt(K)
1355		enddo
1356		endif ; enddo ; enddo
1357
1358	0	end subroutine ALE_remap_interface_vals
1359
1360		!> Interpolate to find an updated array of values at vertices of tracer cells after remapping.
1361	0	subroutine ALE_remap_vertex_vals(CS, G, GV, h_old, h_new, vert_val)
1362		type(ALE_CS), intent(in) :: CS !< ALE control structure
1363		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1364		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
1365		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_old !< Thickness of source grid
1366		!! [H ~> m or kg m-2]
1367		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h_new !< Thickness of destination grid
1368		!! [H ~> m or kg m-2]
1369		real, dimension(SZIB_(G),SZJB_(G),SZK_(GV)+1), &
1370		intent(inout) :: vert_val !< The interface values to interpolate [A]
1371
1372	0	real :: val_src(GV%ke+1) ! A column of interface values on the source grid [A]
1373	0	real :: val_tgt(GV%ke+1) ! A column of interface values on the target grid [A]
1374	0	real :: h_src(GV%ke) ! A column of source grid layer thicknesses [H ~> m or kg m-2]
1375	0	real :: h_tgt(GV%ke) ! A column of target grid layer thicknesses [H ~> m or kg m-2]
1376		real :: I_mask_sum ! The inverse of the tracer point masks surrounding a corner [nondim]
1377		integer :: i, j, k, nz
1378
1379	0	nz = GV%ke
1380
1381	0	do J=G%JscB,G%JecB ; do I=G%IscB,G%IecB
1382	0	if ((G%mask2dT(i,j) + G%mask2dT(i+1,j+1)) + (G%mask2dT(i+1,j) + G%mask2dT(i,j+1)) > 0.0 ) then
1383		I_mask_sum = 1.0 / ((G%mask2dT(i,j) + G%mask2dT(i+1,j+1)) + &
1384	0	(G%mask2dT(i+1,j) + G%mask2dT(i,j+1)))
1385
1386	0	do k=1,nz
1387		h_src(k) = ((G%mask2dT(i,j) * h_old(i,j,k) + G%mask2dT(i+1,j+1) * h_old(i+1,j+1,k)) + &
1388	0	(G%mask2dT(i+1,j) * h_old(i+1,j,k) + G%mask2dT(i,j+1) * h_old(i,j+1,k)) ) * I_mask_sum
1389		h_tgt(k) = ((G%mask2dT(i,j) * h_new(i,j,k) + G%mask2dT(i+1,j+1) * h_new(i+1,j+1,k)) + &
1390	0	(G%mask2dT(i+1,j) * h_new(i+1,j,k) + G%mask2dT(i,j+1) * h_new(i,j+1,k)) ) * I_mask_sum
1391		enddo
1392
1393	0	do K=1,nz+1
1394	0	val_src(K) = vert_val(I,J,K)
1395		enddo
1396
1397	0	call interpolate_column(nz, h_src, val_src, nz, h_tgt, val_tgt, .false.)
1398
1399	0	do K=1,nz+1
1400	0	vert_val(I,J,K) = val_tgt(K)
1401		enddo
1402		endif ; enddo ; enddo
1403
1404	0	end subroutine ALE_remap_vertex_vals
1405
1406		!> Mask out thicknesses to 0 when their running sum exceeds a specified value.
1407	0	subroutine apply_partial_cell_mask(h1, h_mask)
1408		real, dimension(:), intent(inout) :: h1 !< A column of thicknesses to be masked out after their
1409		!! running vertical sum exceeds h_mask [H ~> m or kg m-2]
1410		real, intent(in) :: h_mask !< The depth after which the thicknesses in h1 are
1411		!! masked out [H ~> m or kg m-2]
1412		! Local variables
1413		real :: h1_rsum ! The running sum of h1 [H ~> m or kg m-2]
1414		integer :: k
1415
1416	0	h1_rsum = 0.0
1417	0	do k=1,size(h1)
1418	0	if (h1(k) > h_mask - h1_rsum) then
1419		! This thickness is reduced because it extends below the shallower neighboring bathymetry.
1420	0	h1(k) = max(h_mask - h1_rsum, 0.0)
1421	0	h1_rsum = h_mask
1422		else
1423	0	h1_rsum = h1_rsum + h1(k)
1424		endif
1425		enddo
1426	0	end subroutine apply_partial_cell_mask
1427
1428
1429		!> Zero out velocities in a column in very thin layers near the seafloor
1430	0	subroutine mask_near_bottom_vel(vel, h, h_BBL, h_thin, nk)
1431		integer, intent(in) :: nk !< The number of layers in this column
1432		real, intent(inout) :: vel(nk) !< The velocity component being zeroed out [L T-1 ~> m s-1]
1433		real, intent(in) :: h(nk) !< The layer thicknesses at velocity points [H ~> m or kg m-2]
1434		real, intent(in) :: h_BBL !< The thickness of the near-bottom region over which to apply
1435		!! the filtering [H ~> m or kg m-2]
1436		real, intent(in) :: h_thin !< A layer thickness below which the filtering is applied [H ~> m or kg m-2]
1437
1438		! Local variables
1439		real :: h_from_bot ! The distance between the top of a layer and the seafloor [H ~> m or kg m-2]
1440		integer :: k
1441
1442	0	if ((h_BBL < 0.0) .or. (h_thin < 0.0)) return
1443
1444	0	h_from_bot = 0.0
1445	0	do k=nk,1,-1
1446	0	h_from_bot = h_from_bot + h(k)
1447	0	if (h_from_bot > h_BBL) return
1448		! Set the velocity to zero in thin, near-bottom layers.
1449	0	if (h(k) <= h_thin) vel(k) = 0.0
1450		enddo !k
1451
1452		end subroutine mask_near_bottom_vel
1453
1454
1455		!> Remaps a single scalar between grids described by thicknesses h_src and h_dst.
1456		!! h_dst must be dimensioned as a model array with GV%ke layers while h_src can
1457		!! have an arbitrary number of layers specified by nk_src.
1458	0	subroutine ALE_remap_scalar(CS, G, GV, nk_src, h_src, s_src, h_dst, s_dst, all_cells, old_remap)
1459		type(remapping_CS), intent(in) :: CS !< Remapping control structure
1460		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
1461		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
1462		integer, intent(in) :: nk_src !< Number of levels on source grid
1463		real, dimension(SZI_(G),SZJ_(G),nk_src), intent(in) :: h_src !< Level thickness of source grid
1464		!! [H ~> m or kg m-2] or other units
1465		!! if H_neglect is provided
1466		real, dimension(SZI_(G),SZJ_(G),nk_src), intent(in) :: s_src !< Scalar on source grid, in arbitrary units [A]
1467		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),intent(in) :: h_dst !< Level thickness of destination grid in the
1468		!! same units as h_src, often [H ~> m or kg m-2]
1469		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)),intent(inout) :: s_dst !< Scalar on destination grid, in the same
1470		!! arbitrary units as s_src [A]
1471		logical, optional, intent(in) :: all_cells !< If false, only reconstruct for
1472		!! non-vanished cells. Use all vanished
1473		!! layers otherwise (default).
1474		logical, optional, intent(in) :: old_remap !< If true, use the old "remapping_core_w"
1475		!! method, otherwise use "remapping_core_h".
1476		! Local variables
1477		integer :: i, j, k, n_points
1478	0	real :: dx(GV%ke+1) ! Change in interface position [H ~> m or kg m-2]
1479		logical :: ignore_vanished_layers, use_remapping_core_w
1480
1481	0	ignore_vanished_layers = .false.
1482	0	if (present(all_cells)) ignore_vanished_layers = .not. all_cells
1483	0	use_remapping_core_w = .false.
1484	0	if (present(old_remap)) use_remapping_core_w = old_remap
1485	0	n_points = nk_src
1486
1487		!$OMP parallel do default(shared) firstprivate(n_points,dx)
1488	0	do j = G%jsc,G%jec ; do i = G%isc,G%iec
1489	0	if (G%mask2dT(i,j) > 0.) then
1490	0	if (ignore_vanished_layers) then
1491	0	n_points = 0
1492	0	do k = 1, nk_src
1493	0	if (h_src(i,j,k)>0.) n_points = n_points + 1
1494		enddo
1495	0	s_dst(i,j,:) = 0.
1496		endif
1497	0	if (use_remapping_core_w) then
1498	0	call dzFromH1H2( n_points, h_src(i,j,1:n_points), GV%ke, h_dst(i,j,:), dx )
1499		call remapping_core_w(CS, n_points, h_src(i,j,1:n_points), s_src(i,j,1:n_points), &
1500	0	GV%ke, dx, s_dst(i,j,:))
1501		else
1502		call remapping_core_h(CS, n_points, h_src(i,j,1:n_points), s_src(i,j,1:n_points), &
1503	0	GV%ke, h_dst(i,j,:), s_dst(i,j,:))
1504		endif
1505		else
1506	0	s_dst(i,j,:) = 0.
1507		endif
1508		enddo ; enddo
1509
1510	0	end subroutine ALE_remap_scalar
1511
1512
1513		!> Calculate edge values (top and bottom of layer) for T and S consistent with a PLM reconstruction
1514		!! in the vertical direction. Boundary reconstructions are PCM unless bdry_extrap is true.
1515	24	subroutine TS_PLM_edge_values( CS, S_t, S_b, T_t, T_b, G, GV, tv, h, bdry_extrap )
1516		type(ocean_grid_type), intent(in) :: G !< ocean grid structure
1517		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
1518		type(ALE_CS), intent(inout) :: CS !< module control structure
1519		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1520		intent(inout) :: S_t !< Salinity at the top edge of each layer [S ~> ppt]
1521		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1522		intent(inout) :: S_b !< Salinity at the bottom edge of each layer [S ~> ppt]
1523		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1524		intent(inout) :: T_t !< Temperature at the top edge of each layer [C ~> degC]
1525		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1526		intent(inout) :: T_b !< Temperature at the bottom edge of each layer [C ~> degC]
1527		type(thermo_var_ptrs), intent(in) :: tv !< thermodynamics structure
1528		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1529		intent(in) :: h !< layer thickness [H ~> m or kg m-2]
1530		logical, intent(in) :: bdry_extrap !< If true, use high-order boundary
1531		!! extrapolation within boundary cells
1532
1533	24	call ALE_PLM_edge_values( CS, G, GV, h, tv%S, bdry_extrap, S_t, S_b )
1534	24	call ALE_PLM_edge_values( CS, G, GV, h, tv%T, bdry_extrap, T_t, T_b )
1535
1536	24	end subroutine TS_PLM_edge_values
1537
1538		!> Calculate edge values (top and bottom of layer) 3d scalar array.
1539		!! Boundary reconstructions are PCM unless bdry_extrap is true.
1540	48	subroutine ALE_PLM_edge_values( CS, G, GV, h, Q, bdry_extrap, Q_t, Q_b )
1541		type(ALE_CS), intent(in) :: CS !< module control structure
1542		type(ocean_grid_type), intent(in) :: G !< ocean grid structure
1543		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
1544		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1545		intent(in) :: h !< layer thickness [H ~> m or kg m-2]
1546		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1547		intent(in) :: Q !< 3d scalar array, in arbitrary units [A]
1548		logical, intent(in) :: bdry_extrap !< If true, use high-order boundary
1549		!! extrapolation within boundary cells
1550		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1551		intent(inout) :: Q_t !< Scalar at the top edge of each layer [A]
1552		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1553		intent(inout) :: Q_b !< Scalar at the bottom edge of each layer [A]
1554		! Local variables
1555		integer :: i, j, k
1556	96	real :: slp(GV%ke) ! Tracer slope times the cell width [A]
1557		real :: mslp ! Monotonized tracer slope times the cell width [A]
1558		real :: h_neglect ! Tiny thicknesses used in remapping [H ~> m or kg m-2]
1559
1560	48	if (CS%answer_date >= 20190101) then
1561	48	h_neglect = GV%H_subroundoff
1562	0	elseif (GV%Boussinesq) then
1563	0	h_neglect = GV%m_to_H*1.0e-30
1564		else
1565	0	h_neglect = GV%kg_m2_to_H*1.0e-30
1566		endif
1567
1568		!$OMP parallel do default(shared) private(slp,mslp)
1569	366096	do j = G%jsc-1,G%jec+1 ; do i = G%isc-1,G%iec+1
1570	363072	slp(1) = 0.
1571	26867328	do k = 2, GV%ke-1
1572		slp(k) = PLM_slope_wa(h(i,j,k-1), h(i,j,k), h(i,j,k+1), h_neglect, &
1573	26867328	Q(i,j,k-1), Q(i,j,k), Q(i,j,k+1))
1574		enddo
1575	363072	slp(GV%ke) = 0.
1576
1577	26867328	do k = 2, GV%ke-1
1578	26504256	mslp = PLM_monotonized_slope(Q(i,j,k-1), Q(i,j,k), Q(i,j,k+1), slp(k-1), slp(k), slp(k+1))
1579	26504256	Q_t(i,j,k) = Q(i,j,k) - 0.5 * mslp
1580	26867328	Q_b(i,j,k) = Q(i,j,k) + 0.5 * mslp
1581		enddo
1582	366048	if (bdry_extrap) then
1583	363072	mslp = - PLM_extrapolate_slope(h(i,j,2), h(i,j,1), h_neglect, Q(i,j,2), Q(i,j,1))
1584	363072	Q_t(i,j,1) = Q(i,j,1) - 0.5 * mslp
1585	363072	Q_b(i,j,1) = Q(i,j,1) + 0.5 * mslp
1586		mslp = PLM_extrapolate_slope(h(i,j,GV%ke-1), h(i,j,GV%ke), h_neglect, &
1587	363072	Q(i,j,GV%ke-1), Q(i,j,GV%ke))
1588	363072	Q_t(i,j,GV%ke) = Q(i,j,GV%ke) - 0.5 * mslp
1589	363072	Q_b(i,j,GV%ke) = Q(i,j,GV%ke) + 0.5 * mslp
1590		else
1591	0	Q_t(i,j,1) = Q(i,j,1)
1592	0	Q_b(i,j,1) = Q(i,j,1)
1593	0	Q_t(i,j,GV%ke) = Q(i,j,GV%ke)
1594	0	Q_b(i,j,GV%ke) = Q(i,j,GV%ke)
1595		endif
1596
1597		enddo ; enddo
1598
1599	48	end subroutine ALE_PLM_edge_values
1600
1601		!> Calculate edge values (top and bottom of layer) for T and S consistent with a PPM reconstruction
1602		!! in the vertical direction. Boundary reconstructions are PCM unless bdry_extrap is true.
1603	0	subroutine TS_PPM_edge_values( CS, S_t, S_b, T_t, T_b, G, GV, tv, h, bdry_extrap )
1604		type(ocean_grid_type), intent(in) :: G !< ocean grid structure
1605		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
1606		type(ALE_CS), intent(inout) :: CS !< module control structure
1607		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1608		intent(inout) :: S_t !< Salinity at the top edge of each layer [S ~> ppt]
1609		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1610		intent(inout) :: S_b !< Salinity at the bottom edge of each layer [S ~> ppt]
1611		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1612		intent(inout) :: T_t !< Temperature at the top edge of each layer [C ~> degC]
1613		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1614		intent(inout) :: T_b !< Temperature at the bottom edge of each layer [C ~> degC]
1615		type(thermo_var_ptrs), intent(in) :: tv !< thermodynamics structure
1616		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1617		intent(in) :: h !< layer thicknesses [H ~> m or kg m-2]
1618		logical, intent(in) :: bdry_extrap !< If true, use high-order boundary
1619		!! extrapolation within boundary cells
1620
1621		! Local variables
1622		integer :: i, j, k
1623	0	real :: hTmp(GV%ke) ! A 1-d copy of h [H ~> m or kg m-2]
1624	0	real :: tmp(GV%ke) ! A 1-d copy of a column of temperature [C ~> degC] or salinity [S ~> ppt]
1625		real, dimension(CS%nk,2) :: &
1626	0	ppol_E ! Edge value of polynomial in [C ~> degC] or [S ~> ppt]
1627		real, dimension(CS%nk,3) :: &
1628	0	ppol_coefs ! Coefficients of polynomial, all in [C ~> degC] or [S ~> ppt]
1629		real :: h_neglect, h_neglect_edge ! Tiny thicknesses [H ~> m or kg m-2]
1630
1631	0	if (CS%answer_date >= 20190101) then
1632	0	h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff
1633	0	elseif (GV%Boussinesq) then
1634	0	h_neglect = GV%m_to_H1.0e-30 ; h_neglect_edge = GV%m_to_H1.0e-10
1635		else
1636	0	h_neglect = GV%kg_m2_to_H1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H1.0e-10
1637		endif
1638
1639		! Determine reconstruction within each column
1640		!$OMP parallel do default(shared) private(hTmp,tmp,ppol_E,ppol_coefs)
1641	0	do j = G%jsc-1,G%jec+1 ; do i = G%isc-1,G%iec+1
1642
1643		! Build current grid
1644	0	hTmp(:) = h(i,j,:)
1645	0	tmp(:) = tv%S(i,j,:)
1646
1647		! Reconstruct salinity profile
1648	0	ppol_E(:,:) = 0.0
1649	0	ppol_coefs(:,:) = 0.0
1650		call edge_values_implicit_h4( GV%ke, hTmp, tmp, ppol_E, h_neglect=h_neglect_edge, &
1651	0	answer_date=CS%answer_date )
1652		call PPM_reconstruction( GV%ke, hTmp, tmp, ppol_E, ppol_coefs, h_neglect, &
1653	0	answer_date=CS%answer_date )
1654	0	if (bdry_extrap) &
1655	0	call PPM_boundary_extrapolation( GV%ke, hTmp, tmp, ppol_E, ppol_coefs, h_neglect )
1656
1657	0	do k = 1,GV%ke
1658	0	S_t(i,j,k) = ppol_E(k,1)
1659	0	S_b(i,j,k) = ppol_E(k,2)
1660		enddo
1661
1662		! Reconstruct temperature profile
1663	0	ppol_E(:,:) = 0.0
1664	0	ppol_coefs(:,:) = 0.0
1665	0	tmp(:) = tv%T(i,j,:)
1666	0	if (CS%answer_date < 20190101) then
1667		call edge_values_implicit_h4( GV%ke, hTmp, tmp, ppol_E, h_neglect=1.0e-10*GV%m_to_H, &
1668	0	answer_date=CS%answer_date )
1669		else
1670		call edge_values_implicit_h4( GV%ke, hTmp, tmp, ppol_E, h_neglect=GV%H_subroundoff, &
1671	0	answer_date=CS%answer_date )
1672		endif
1673		call PPM_reconstruction( GV%ke, hTmp, tmp, ppol_E, ppol_coefs, h_neglect, &
1674	0	answer_date=CS%answer_date )
1675	0	if (bdry_extrap) &
1676	0	call PPM_boundary_extrapolation(GV%ke, hTmp, tmp, ppol_E, ppol_coefs, h_neglect )
1677
1678	0	do k = 1,GV%ke
1679	0	T_t(i,j,k) = ppol_E(k,1)
1680	0	T_b(i,j,k) = ppol_E(k,2)
1681		enddo
1682
1683		enddo ; enddo
1684
1685	0	end subroutine TS_PPM_edge_values
1686
1687		!> Calculate edge values (top and bottom of layer) for T and S consistent with a PLM reconstruction
1688		!! in the vertical direction that uses weighted least squares for the slope.
1689	0	subroutine TS_PLM_WLS_edge_values(CS, S_t, S_b, T_t, T_b, G, GV, tv, h)
1690		type(ocean_grid_type), intent(in) :: G !< ocean grid structure
1691		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
1692		type(ALE_CS), intent(inout) :: CS !< module control structure
1693		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1694		intent(inout) :: S_t !< Salinity at the top edge of each layer [S ~> ppt]
1695		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1696		intent(inout) :: S_b !< Salinity at the bottom edge of each layer [S ~> ppt]
1697		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1698		intent(inout) :: T_t !< Temperature at the top edge of each layer [C ~> degC]
1699		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1700		intent(inout) :: T_b !< Temperature at the bottom edge of each layer [C ~> degC]
1701		type(thermo_var_ptrs), intent(in) :: tv !< thermodynamics structure
1702		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
1703		intent(in) :: h !< layer thickness [H ~> m or kg m-2]
1704		! Local variables
1705		integer :: i, j, k
1706	0	type(PLM_WLS) :: recon !< A PLM-WLS reconstruction
1707
1708	0	call recon%init(GV%ke, h_neglect=GV%H_subroundoff)
1709
1710		!$OMP parallel do default(shared) firstprivate(recon)
1711	0	do j = G%jsc-1,G%jec+1 ; do i = G%isc-1,G%iec+1
1712
1713	0	call recon%reconstruct(h(i,j,:), tv%T(i,j,:))
1714	0	T_t(i,j,:) = recon%ul(:)
1715	0	T_b(i,j,:) = recon%ur(:)
1716
1717	0	call recon%reconstruct(h(i,j,:), tv%S(i,j,:))
1718	0	S_t(i,j,:) = recon%ul(:)
1719	0	S_b(i,j,:) = recon%ur(:)
1720
1721		enddo ; enddo
1722
1723	0	call recon%destroy()
1724
1725	0	end subroutine TS_PLM_WLS_edge_values
1726
1727		!> Initializes regridding for the main ALE algorithm
1728	1	subroutine ALE_initRegridding(G, GV, US, max_depth, param_file, mdl, regridCS)
1729		type(ocean_grid_type), intent(in) :: G !< Grid structure
1730		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
1731		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1732		real, intent(in) :: max_depth !< The maximum depth of the ocean [Z ~> m].
1733		type(param_file_type), intent(in) :: param_file !< parameter file
1734		character(len=*), intent(in) :: mdl !< Name of calling module
1735		type(regridding_CS), intent(out) :: regridCS !< Regridding parameters and work arrays
1736		! Local variables
1737		character(len=30) :: coord_mode
1738
1739		call get_param(param_file, mdl, "REGRIDDING_COORDINATE_MODE", coord_mode, &
1740		"Coordinate mode for vertical regridding. "//&
1741		"Choose among the following possibilities: "//&
1742		trim(regriddingCoordinateModeDoc), &
1743	1	default=DEFAULT_COORDINATE_MODE, fail_if_missing=.true.)
1744
1745	1	call initialize_regridding(regridCS, G, GV, US, max_depth, param_file, mdl, coord_mode, '', '')
1746
1747	1	end subroutine ALE_initRegridding
1748
1749		!> Query the target coordinate interfaces positions
1750	0	function ALE_getCoordinate( CS )
1751		type(ALE_CS), pointer :: CS !< module control structure
1752
1753		real, dimension(CS%nk+1) :: ALE_getCoordinate !< The coordinate positions, in the appropriate units
1754		!! of the target coordinate, e.g. [Z ~> m] for z*,
1755		!! non-dimensional for sigma, etc.
1756	0	ALE_getCoordinate(:) = getCoordinateInterfaces( CS%regridCS, undo_scaling=.true. )
1757
1758		end function ALE_getCoordinate
1759
1760
1761		!> Query the target coordinate units
1762	0	function ALE_getCoordinateUnits( CS )
1763		type(ALE_CS), pointer :: CS !< module control structure
1764
1765		character(len=20) :: ALE_getCoordinateUnits
1766
1767	0	ALE_getCoordinateUnits = getCoordinateUnits( CS%regridCS )
1768
1769	0	end function ALE_getCoordinateUnits
1770
1771
1772		!> Returns true if initial conditions should be regridded and remapped
1773	1	logical function ALE_remap_init_conds( CS )
1774		type(ALE_CS), pointer :: CS !< module control structure
1775
1776	1	ALE_remap_init_conds = .false.
1777	1	if (associated(CS)) ALE_remap_init_conds = CS%remap_after_initialization
1778	1	end function ALE_remap_init_conds
1779
1780		!> Updates the weights for time filtering the new grid generated in regridding
1781	12	subroutine ALE_update_regrid_weights( dt, CS )
1782		real, intent(in) :: dt !< Time-step used between ALE calls [T ~> s]
1783		type(ALE_CS), pointer :: CS !< ALE control structure
1784		! Local variables
1785		real :: w ! An implicit weighting estimate [nondim]
1786
1787	12	if (associated(CS)) then
1788	12	w = 0.0
1789	12	if (CS%regrid_time_scale > 0.0) then
1790	0	w = CS%regrid_time_scale / (CS%regrid_time_scale + dt)
1791		endif
1792	12	call set_regrid_params(CS%regridCS, old_grid_weight=w)
1793		endif
1794
1795	12	end subroutine ALE_update_regrid_weights
1796
1797		!> Update the vertical grid type with ALE information.
1798		!! This subroutine sets information in the verticalGrid_type to be
1799		!! consistent with the use of ALE mode.
1800	1	subroutine ALE_updateVerticalGridType(CS, GV)
1801		type(ALE_CS), pointer :: CS !< ALE control structure
1802		type(verticalGrid_type), pointer :: GV !< vertical grid information
1803
1804		integer :: nk
1805
1806	1	nk = GV%ke
1807	77	GV%sInterface(1:nk+1) = getCoordinateInterfaces( CS%regridCS, undo_scaling=.true. )
1808	76	GV%sLayer(1:nk) = 0.5*( GV%sInterface(1:nk) + GV%sInterface(2:nk+1) )
1809	1	GV%zAxisUnits = getCoordinateUnits( CS%regridCS )
1810	1	GV%zAxisLongName = getCoordinateShortName( CS%regridCS )
1811	1	GV%direction = -1 ! Because of ferret in z* mode. Need method to set
1812		! as function of coordinate mode.
1813
1814	1	end subroutine ALE_updateVerticalGridType
1815
1816
1817		!> Write the vertical coordinate information into a file.
1818		!! This subroutine writes out a file containing any available data related
1819		!! to the vertical grid used by the MOM ocean model when in ALE mode.
1820	1	subroutine ALE_writeCoordinateFile( CS, GV, directory )
1821		type(ALE_CS), pointer :: CS !< module control structure
1822		type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
1823		character(len=*), intent(in) :: directory !< directory for writing grid info
1824
1825		character(len=240) :: filepath
1826
1827	1	filepath = trim(directory) // trim("Vertical_coordinate.nc")
1828
1829	1	call write_regrid_file(CS%regridCS, GV, filepath)
1830
1831	1	end subroutine ALE_writeCoordinateFile
1832
1833		!> Set h to coordinate values for fixed coordinate systems
1834	0	subroutine ALE_initThicknessToCoord( CS, G, GV, h, height_units )
1835		type(ALE_CS), intent(inout) :: CS !< module control structure
1836		type(ocean_grid_type), intent(in) :: G !< module grid structure
1837		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid structure
1838		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: h !< layer thickness in thickness units
1839		!! [H ~> m or kg m-2] or height units [Z ~> m]
1840		logical, optional, intent(in) :: height_units !< If present and true, the
1841		!! thicknesses are in height units
1842
1843		! Local variables
1844		real :: scale ! A scaling value for the thicknesses [nondim] or [H Z-1 ~> nondim or kg m-3]
1845		integer :: i, j
1846
1847	0	scale = GV%Z_to_H
1848	0	if (present(height_units)) then ; if (height_units) scale = 1.0 ; endif
1849	0	do j = G%jsd,G%jed ; do i = G%isd,G%ied
1850	0	h(i,j,:) = scale * getStaticThickness( CS%regridCS, 0., max(G%meanSL(i,j)+G%bathyT(i,j), 0.0) )
1851		enddo ; enddo
1852
1853	0	end subroutine ALE_initThicknessToCoord
1854
1855	0	end module MOM_ALE