src/tracer/MOM_neutral_diffusion.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		!> A column-wise toolbox for implementing neutral diffusion
6		module MOM_neutral_diffusion
7
8		use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
9		use MOM_cpu_clock, only : CLOCK_MODULE, CLOCK_ROUTINE
10		use MOM_domains, only : pass_var
11		use MOM_diag_mediator, only : diag_ctrl, time_type
12		use MOM_diag_mediator, only : post_data, register_diag_field
13		use MOM_EOS, only : EOS_type, EOS_manual_init, EOS_domain
14		use MOM_EOS, only : calculate_density, calculate_density_derivs
15		use MOM_EOS, only : EOS_LINEAR
16		use MOM_error_handler, only : MOM_error, FATAL, WARNING, MOM_mesg, is_root_pe
17		use MOM_file_parser, only : get_param, log_version, param_file_type
18		use MOM_file_parser, only : openParameterBlock, closeParameterBlock
19		use MOM_grid, only : ocean_grid_type
20		use MOM_remapping, only : remapping_CS, initialize_remapping
21		use MOM_remapping, only : extract_member_remapping_CS, build_reconstructions_1d
22		use MOM_remapping, only : average_value_ppoly, remappingSchemesDoc, remappingDefaultScheme
23		use MOM_tracer_registry, only : tracer_registry_type, tracer_type
24		use MOM_unit_scaling, only : unit_scale_type
25		use MOM_variables, only : vertvisc_type
26		use MOM_verticalGrid, only : verticalGrid_type
27		use polynomial_functions, only : evaluation_polynomial, first_derivative_polynomial
28		use PPM_functions, only : PPM_reconstruction, PPM_boundary_extrapolation
29		use regrid_edge_values, only : edge_values_implicit_h4
30		use MOM_CVMix_KPP, only : KPP_get_BLD, KPP_CS
31		use MOM_energetic_PBL, only : energetic_PBL_get_MLD, energetic_PBL_CS
32		use MOM_diabatic_driver, only : diabatic_CS, extract_diabatic_member
33		use MOM_io, only : stdout, stderr
34		use MOM_hor_bnd_diffusion, only : boundary_k_range, SURFACE, BOTTOM
35
36		implicit none ; private
37
38		#include <MOM_memory.h>
39
40		public neutral_diffusion, neutral_diffusion_init, neutral_diffusion_end
41		public neutral_diffusion_calc_coeffs
42		public neutral_diffusion_unit_tests
43
44		!> The control structure for the MOM_neutral_diffusion module
45		type, public :: neutral_diffusion_CS ; private
46		integer :: nkp1 !< Number of interfaces for a column = nk + 1
47		integer :: nsurf !< Number of neutral surfaces
48		integer :: deg = 2 !< Degree of polynomial used for reconstructions
49		logical :: continuous_reconstruction = .true. !< True if using continuous PPM reconstruction at interfaces
50		logical :: debug = .false. !< If true, write verbose debugging messages
51		logical :: hard_fail_heff !< Bring down the model if a problem with heff is detected
52		integer :: max_iter !< Maximum number of iterations if refine_position is defined
53		real :: drho_tol !< Convergence criterion representing density difference from true neutrality [R ~> kg m-3]
54		real :: x_tol !< Convergence criterion for how small an update of the position can be [nondim]
55		real :: ref_pres !< Reference pressure, negative if using locally referenced neutral
56		!! density [R L2 T-2 ~> Pa]
57		logical :: interior_only !< If true, only applies neutral diffusion in the ocean interior.
58		!! That is, the algorithm will exclude the surface and bottom boundary layers.
59		logical :: tapering = .false. !< If true, neutral diffusion linearly decays towards zero within a
60		!! transition zone defined using boundary layer depths. Only available when
61		!! interior_only=true.
62		logical :: KhTh_use_vert_struct !< If true, uses vertical structure
63		!! for tracer diffusivity.
64		logical :: use_unmasked_transport_bug !< If true, use an older form for the accumulation of
65		!! neutral-diffusion transports that were unmasked, as used prior to Jan 2018.
66		real, allocatable, dimension(:,:) :: hbl !< Boundary layer depth [H ~> m or kg m-2]
67		! Coefficients used to apply tapering from neutral to horizontal direction
68		real, allocatable, dimension(:) :: coeff_l !< Non-dimensional coefficient in the left column,
69		!! at cell interfaces [nondim]
70		real, allocatable, dimension(:) :: coeff_r !< Non-dimensional coefficient in the right column,
71		!! at cell interfaces [nondim]
72		! Array used when KhTh_use_vert_struct is true
73		real, allocatable, dimension(:,:,:) :: Coef_h !< Coef_x and Coef_y averaged at t-points [L2 ~> m2]
74		! Positions of neutral surfaces in both the u, v directions
75		real, allocatable, dimension(:,:,:) :: uPoL !< Non-dimensional position with left layer uKoL-1, u-point [nondim]
76		real, allocatable, dimension(:,:,:) :: uPoR !< Non-dimensional position with right layer uKoR-1, u-point [nondim]
77		integer, allocatable, dimension(:,:,:) :: uKoL !< Index of left interface corresponding to neutral surface,
78		!! at a u-point
79		integer, allocatable, dimension(:,:,:) :: uKoR !< Index of right interface corresponding to neutral surface,
80		!! at a u-point
81		real, allocatable, dimension(:,:,:) :: uHeff !< Effective thickness at u-point [H ~> m or kg m-2]
82		real, allocatable, dimension(:,:,:) :: vPoL !< Non-dimensional position with left layer uKoL-1, v-point [nondim]
83		real, allocatable, dimension(:,:,:) :: vPoR !< Non-dimensional position with right layer uKoR-1, v-point [nondim]
84		integer, allocatable, dimension(:,:,:) :: vKoL !< Index of left interface corresponding to neutral surface,
85		!! at a v-point
86		integer, allocatable, dimension(:,:,:) :: vKoR !< Index of right interface corresponding to neutral surface,
87		!! at a v-point
88		real, allocatable, dimension(:,:,:) :: vHeff !< Effective thickness at v-point [H ~> m or kg m-2]
89		! Coefficients of polynomial reconstructions for temperature and salinity
90		real, allocatable, dimension(:,:,:,:) :: ppoly_coeffs_T !< Polynomial coefficients of the
91		!! sub-gridscale temperatures [C ~> degC]
92		real, allocatable, dimension(:,:,:,:) :: ppoly_coeffs_S !< Polynomial coefficients of the
93		!! sub-gridscale salinity [S ~> ppt]
94		! Variables needed for continuous reconstructions
95		real, allocatable, dimension(:,:,:) :: dRdT !< dRho/dT [R C-1 ~> kg m-3 degC-1] at interfaces
96		real, allocatable, dimension(:,:,:) :: dRdS !< dRho/dS [R S-1 ~> kg m-3 ppt-1] at interfaces
97		real, allocatable, dimension(:,:,:) :: Tint !< Interface T [C ~> degC]
98		real, allocatable, dimension(:,:,:) :: Sint !< Interface S [S ~> ppt]
99		real, allocatable, dimension(:,:,:) :: Pint !< Interface pressure [R L2 T-2 ~> Pa]
100		! Variables needed for discontinuous reconstructions
101		real, allocatable, dimension(:,:,:,:) :: T_i !< Top edge reconstruction of temperature [C ~> degC]
102		real, allocatable, dimension(:,:,:,:) :: S_i !< Top edge reconstruction of salinity [S ~> ppt]
103		real, allocatable, dimension(:,:,:,:) :: P_i !< Interface pressures [R L2 T-2 ~> Pa]
104		real, allocatable, dimension(:,:,:,:) :: dRdT_i !< dRho/dT [R C-1 ~> kg m-3 degC-1] at top edge
105		real, allocatable, dimension(:,:,:,:) :: dRdS_i !< dRho/dS [R S-1 ~> kg m-3 ppt-1] at top edge
106		integer, allocatable, dimension(:,:) :: ns !< Number of interfaces in a column
107		logical, allocatable, dimension(:,:,:) :: stable_cell !< True if the cell is stably stratified wrt to the next cell
108		real :: R_to_kg_m3 = 1.0 !< A rescaling factor translating density to kg m-3 for
109		!! use in diagnostic messages [kg m-3 R-1 ~> 1].
110		type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to
111		!! regulate the timing of diagnostic output.
112		integer :: neutral_pos_method !< Method to find the position of a neutral surface within the layer
113		character(len=40) :: delta_rho_form !< Determine which (if any) approximation is made to the
114		!! equation describing the difference in density
115
116		integer :: id_uhEff_2d = -1 !< Diagnostic IDs
117		integer :: id_vhEff_2d = -1 !< Diagnostic IDs
118
119		type(EOS_type), pointer :: EOS => NULL() !< Equation of state parameters
120		type(remapping_CS) :: remap_CS !< Remapping control structure used to create sublayers
121		integer :: remap_answer_date !< The vintage of the order of arithmetic and expressions to use
122		!! for remapping. Values below 20190101 recover the remapping
123		!! answers from 2018, while higher values use more robust
124		!! forms of the same remapping expressions.
125		integer :: ndiff_answer_date !< The vintage of the order of arithmetic to use for the neutral
126		!! diffusion. Values of 20240330 or below recover the answers
127		!! from the original form of this code, while higher values use
128		!! mathematically equivalent expressions that recover rotational symmetry.
129		type(KPP_CS), pointer :: KPP_CSp => NULL() !< KPP control structure needed to get BLD
130		type(energetic_PBL_CS), pointer :: energetic_PBL_CSp => NULL()!< ePBL control structure needed to get MLD
131		end type neutral_diffusion_CS
132
133		! This include declares and sets the variable "version".
134		#include "version_variable.h"
135		character(len=40) :: mdl = "MOM_neutral_diffusion" !< module name
136
137		contains
138
139		!> Read parameters and allocate control structure for neutral_diffusion module.
140	1	logical function neutral_diffusion_init(Time, G, GV, US, param_file, diag, EOS, diabatic_CSp, CS)
141		type(time_type), target, intent(in) :: Time !< Time structure
142		type(ocean_grid_type), intent(in) :: G !< Grid structure
143		type(verticalGrid_type), intent(in) :: GV !< The ocean's vertical grid structure
144		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
145		type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure
146		type(param_file_type), intent(in) :: param_file !< Parameter file structure
147		type(EOS_type), target, intent(in) :: EOS !< Equation of state
148		type(diabatic_CS), pointer :: diabatic_CSp!< diabatic control structure needed to get BLD
149		type(neutral_diffusion_CS), pointer :: CS !< Neutral diffusion control structure
150
151		! Local variables
152		character(len=80) :: string ! Temporary strings
153		integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.
154		logical :: debug ! If true, write verbose checksums for debugging purposes.
155		logical :: boundary_extrap ! Indicate whether high-order boundary
156		!! extrapolation should be used within boundary cells.
157		logical :: om4_remap_via_sub_cells ! If true, use the OM4 remapping algorithm
158		logical :: KhTh_use_ebt_struct, KhTh_use_sqg_struct
159
160	1	if (associated(CS)) then
161	0	call MOM_error(FATAL, "neutral_diffusion_init called with associated control structure.")
162	0	return
163		endif
164
165		! Log this module and master switch for turning it on/off
166		call get_param(param_file, mdl, "USE_NEUTRAL_DIFFUSION", neutral_diffusion_init, &
167	1	default=.false., do_not_log=.true.)
168		call log_version(param_file, mdl, version, &
169		"This module implements neutral diffusion of tracers", &
170	1	all_default=.not.neutral_diffusion_init)
171		call get_param(param_file, mdl, "USE_NEUTRAL_DIFFUSION", neutral_diffusion_init, &
172		"If true, enables the neutral diffusion module.", &
173	1	default=.false.)
174
175	1	if (.not.neutral_diffusion_init) return
176
177	0	allocate(CS)
178	0	CS%diag => diag
179	0	CS%EOS => EOS
180		! call openParameterBlock(param_file,'NEUTRAL_DIFF')
181
182		! Read all relevant parameters and write them to the model log.
183		call get_param(param_file, mdl, "NDIFF_CONTINUOUS", CS%continuous_reconstruction, &
184		"If true, uses a continuous reconstruction of T and S when "//&
185		"finding neutral surfaces along which diffusion will happen. "//&
186		"If false, a PPM discontinuous reconstruction of T and S "//&
187		"is done which results in a higher order routine but exacts "//&
188	0	"a higher computational cost.", default=.true.)
189		call get_param(param_file, mdl, "NDIFF_REF_PRES", CS%ref_pres, &
190		"The reference pressure (Pa) used for the derivatives of "//&
191		"the equation of state. If negative (default), local pressure is used.", &
192	0	units="Pa", default=-1., scale=US%Pa_to_RL2_T2)
193		call get_param(param_file, mdl, "NDIFF_INTERIOR_ONLY", CS%interior_only, &
194		"If true, only applies neutral diffusion in the ocean interior. "//&
195		"That is, the algorithm will exclude the surface and bottom "//&
196	0	"boundary layers.", default=.false.)
197	0	if (CS%interior_only) then
198		call get_param(param_file, mdl, "NDIFF_TAPERING", CS%tapering, &
199		"If true, neutral diffusion linearly decays to zero within "//&
200		"a transition zone defined using boundary layer depths. "//&
201	0	"Only applicable when NDIFF_INTERIOR_ONLY=True", default=.false.)
202		endif
203		call get_param(param_file, mdl, "KHTR_USE_EBT_STRUCT", KhTh_use_ebt_struct, &
204		"If true, uses the equivalent barotropic structure "//&
205		"as the vertical structure of the tracer diffusivity.",&
206	0	default=.false.,do_not_log=.true.)
207		call get_param(param_file, mdl, "KHTR_USE_SQG_STRUCT", KhTh_use_sqg_struct, &
208		"If true, uses the surface geostrophic structure "//&
209		"as the vertical structure of the tracer diffusivity.",&
210	0	default=.false.,do_not_log=.true.)
211		call get_param(param_file, mdl, "NDIFF_USE_UNMASKED_TRANSPORT_BUG", CS%use_unmasked_transport_bug, &
212		"If true, use an older form for the accumulation of neutral-diffusion "//&
213		"transports that were unmasked, as used prior to Jan 2018. This is not "//&
214	0	"recommended.", default=.false.)
215
216		call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
217		"This sets the default value for the various _ANSWER_DATE parameters.", &
218	0	default=99991231)
219		call get_param(param_file, mdl, "NDIFF_ANSWER_DATE", CS%ndiff_answer_date, &
220		"The vintage of the order of arithmetic to use for the neutral diffusion. "//&
221		"Values of 20240330 or below recover the answers from the original form of the "//&
222		"neutral diffusion code, while higher values use mathematically equivalent "//&
223		"expressions that recover rotational symmetry.", &
224	0	default=default_answer_date)
225
226		! Initialize and configure remapping
227	0	if ( .not.CS%continuous_reconstruction ) then
228		call get_param(param_file, mdl, "NDIFF_BOUNDARY_EXTRAP", boundary_extrap, &
229		"Extrapolate at the top and bottommost cells, otherwise \n"// &
230		"assume boundaries are piecewise constant", &
231	0	default=.false.)
232		call get_param(param_file, mdl, "NDIFF_REMAPPING_SCHEME", string, &
233		"This sets the reconstruction scheme used "//&
234		"for vertical remapping for all variables. "//&
235		"It can be one of the following schemes: "//&
236	0	trim(remappingSchemesDoc), default=remappingDefaultScheme)
237		call get_param(param_file, mdl, "REMAPPING_ANSWER_DATE", CS%remap_answer_date, &
238		"The vintage of the expressions and order of arithmetic to use for remapping. "//&
239		"Values below 20190101 result in the use of older, less accurate expressions "//&
240		"that were in use at the end of 2018. Higher values result in the use of more "//&
241		"robust and accurate forms of mathematically equivalent expressions.", &
242	0	default=default_answer_date, do_not_log=.not.GV%Boussinesq)
243		call get_param(param_file, mdl, "REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
244	0	do_not_log=.true., default=.true.)
245		call get_param(param_file, mdl, "NDIFF_REMAPPING_USE_OM4_SUBCELLS", om4_remap_via_sub_cells, &
246		"If true, use the OM4 remapping-via-subcells algorithm for neutral diffusion. "//&
247		"See REMAPPING_USE_OM4_SUBCELLS for more details. "//&
248	0	"We recommend setting this option to false.", default=om4_remap_via_sub_cells)
249	0	if (.not.GV%Boussinesq) CS%remap_answer_date = max(CS%remap_answer_date, 20230701)
250		call initialize_remapping( CS%remap_CS, string, boundary_extrapolation=boundary_extrap, &
251		om4_remap_via_sub_cells=om4_remap_via_sub_cells, &
252	0	answer_date=CS%remap_answer_date )
253	0	call extract_member_remapping_CS(CS%remap_CS, degree=CS%deg)
254		call get_param(param_file, mdl, "NEUTRAL_POS_METHOD", CS%neutral_pos_method, &
255		"Method used to find the neutral position \n"// &
256		"1. Delta_rho varies linearly, find 0 crossing \n"// &
257		"2. Alpha and beta vary linearly from top to bottom, \n"// &
258		" Newton's method for neutral position \n"// &
259		"3. Full nonlinear equation of state, use regula falsi \n"// &
260	0	" for neutral position", default=3)
261	0	if (CS%neutral_pos_method > 4 .or. CS%neutral_pos_method < 0) then
262	0	call MOM_error(FATAL,"Invalid option for NEUTRAL_POS_METHOD")
263		endif
264
265		call get_param(param_file, mdl, "DELTA_RHO_FORM", CS%delta_rho_form, &
266		"Determine how the difference in density is calculated \n"// &
267		" full : Difference of in-situ densities \n"// &
268		" no_pressure: Calculated from dRdT, dRdS, but no \n"// &
269		" pressure dependence", &
270	0	default="mid_pressure")
271	0	if (CS%neutral_pos_method > 1) then
272		call get_param(param_file, mdl, "NDIFF_DRHO_TOL", CS%drho_tol, &
273		"Sets the convergence criterion for finding the neutral "// &
274		"position within a layer in kg m-3.", &
275	0	units="kg m-3", default=1.e-10, scale=US%kg_m3_to_R)
276		call get_param(param_file, mdl, "NDIFF_X_TOL", CS%x_tol, &
277		"Sets the convergence criterion for a change in nondimensional "// &
278		"position within a layer.", &
279	0	units="nondim", default=0.)
280		call get_param(param_file, mdl, "NDIFF_MAX_ITER", CS%max_iter, &
281		"The maximum number of iterations to be done before "// &
282		"exiting the iterative loop to find the neutral surface", &
283	0	default=10)
284		endif
285	0	call get_param(param_file, mdl, "DEBUG", debug, default=.false., do_not_log=.true.)
286		call get_param(param_file, mdl, "NDIFF_DEBUG", CS%debug, &
287		"Turns on verbose output for discontinuous neutral "//&
288	0	"diffusion routines.", default=debug)
289		call get_param(param_file, mdl, "HARD_FAIL_HEFF", CS%hard_fail_heff, &
290		"Bring down the model if a problem with heff is detected",&
291	0	default=.true.)
292		endif
293
294	0	if (CS%interior_only) then
295	0	allocate(CS%hbl(SZI_(G),SZJ_(G)), source=0.)
296	0	call extract_diabatic_member(diabatic_CSp, KPP_CSp=CS%KPP_CSp)
297	0	call extract_diabatic_member(diabatic_CSp, energetic_PBL_CSp=CS%energetic_PBL_CSp)
298	0	if ( .not. ASSOCIATED(CS%energetic_PBL_CSp) .and. .not. ASSOCIATED(CS%KPP_CSp) ) then
299	0	call MOM_error(FATAL,"NDIFF_INTERIOR_ONLY is true, but no valid boundary layer scheme was found")
300		endif
301
302	0	if (CS%tapering) then
303	0	allocate(CS%coeff_l(SZK_(GV)+1), source=1.)
304	0	allocate(CS%coeff_r(SZK_(GV)+1), source=1.)
305		endif
306		endif
307
308	0	CS%KhTh_use_vert_struct = KhTh_use_ebt_struct .or. KhTh_use_sqg_struct
309	0	if (CS%KhTh_use_vert_struct) &
310	0	allocate(CS%Coef_h(G%isd:G%ied,G%jsd:G%jed,SZK_(GV)+1), source=0.)
311
312		! Store a rescaling factor for use in diagnostic messages.
313	0	CS%R_to_kg_m3 = US%R_to_kg_m3
314
315		! call closeParameterBlock(param_file)
316	0	if (CS%continuous_reconstruction) then
317	0	CS%nsurf = 2*GV%ke+2 ! Continuous reconstruction means that every interface has two connections
318	0	allocate(CS%dRdT(SZI_(G),SZJ_(G),SZK_(GV)+1), source=0.)
319	0	allocate(CS%dRdS(SZI_(G),SZJ_(G),SZK_(GV)+1), source=0.)
320		else
321	0	CS%nsurf = 4*GV%ke ! Discontinuous means that every interface has four connections
322	0	allocate(CS%T_i(SZI_(G),SZJ_(G),SZK_(GV),2), source=0.)
323	0	allocate(CS%S_i(SZI_(G),SZJ_(G),SZK_(GV),2), source=0.)
324	0	allocate(CS%P_i(SZI_(G),SZJ_(G),SZK_(GV),2), source=0.)
325	0	allocate(CS%dRdT_i(SZI_(G),SZJ_(G),SZK_(GV),2), source=0.)
326	0	allocate(CS%dRdS_i(SZI_(G),SZJ_(G),SZK_(GV),2), source=0.)
327	0	allocate(CS%ppoly_coeffs_T(SZI_(G),SZJ_(G),SZK_(GV),CS%deg+1), source=0.)
328	0	allocate(CS%ppoly_coeffs_S(SZI_(G),SZJ_(G),SZK_(GV),CS%deg+1), source=0.)
329	0	allocate(CS%ns(SZI_(G),SZJ_(G)), source=0)
330		endif
331		! T-points
332	0	allocate(CS%Tint(SZI_(G),SZJ_(G),SZK_(GV)+1), source=0.)
333	0	allocate(CS%Sint(SZI_(G),SZJ_(G),SZK_(GV)+1), source=0.)
334	0	allocate(CS%Pint(SZI_(G),SZJ_(G),SZK_(GV)+1), source=0.)
335	0	allocate(CS%stable_cell(SZI_(G),SZJ_(G),SZK_(GV)), source=.true.)
336		! U-points
337	0	allocate(CS%uPoL(G%isd:G%ied,G%jsd:G%jed, CS%nsurf), source=0.)
338	0	allocate(CS%uPoR(G%isd:G%ied,G%jsd:G%jed, CS%nsurf), source=0.)
339	0	allocate(CS%uKoL(G%isd:G%ied,G%jsd:G%jed, CS%nsurf), source=0)
340	0	allocate(CS%uKoR(G%isd:G%ied,G%jsd:G%jed, CS%nsurf), source=0)
341	0	allocate(CS%uHeff(G%isd:G%ied,G%jsd:G%jed,CS%nsurf-1), source=0.)
342		! V-points
343	0	allocate(CS%vPoL(G%isd:G%ied,G%jsd:G%jed, CS%nsurf), source=0.)
344	0	allocate(CS%vPoR(G%isd:G%ied,G%jsd:G%jed, CS%nsurf), source=0.)
345	0	allocate(CS%vKoL(G%isd:G%ied,G%jsd:G%jed, CS%nsurf), source=0)
346	0	allocate(CS%vKoR(G%isd:G%ied,G%jsd:G%jed, CS%nsurf), source=0)
347	0	allocate(CS%vHeff(G%isd:G%ied,G%jsd:G%jed,CS%nsurf-1), source=0.)
348
349	0	end function neutral_diffusion_init
350
351		!> Calculate remapping factors for u/v columns used to map adjoining columns to
352		!! a shared coordinate space.
353	0	subroutine neutral_diffusion_calc_coeffs(G, GV, US, h, T, S, visc, CS, p_surf)
354		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
355		type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
356		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
357		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
358		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: T !< Potential temperature [C ~> degC]
359		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: S !< Salinity [S ~> ppt]
360		type(vertvisc_type), intent(in) :: visc !< Structure with vertical viscosities,
361		!! boundary layer properties and related fields
362		type(neutral_diffusion_CS), pointer :: CS !< Neutral diffusion control structure
363		real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: p_surf !< Surface pressure to include in pressures used
364		!! for equation of state calculations [R L2 T-2 ~> Pa]
365
366		! Local variables
367		integer, dimension(2) :: EOSdom ! The i-computational domain for the equation of state
368		integer :: i, j, k
369		! Variables used for reconstructions
370	0	real, dimension(SZK_(GV),2) :: ppoly_r_S ! Reconstruction slopes that are unused here, in units of a vertical
371		! gradient, which for temperature would be [C H-1 ~> degC m-1 or degC m2 kg-1].
372	0	real, dimension(SZI_(G), SZJ_(G)) :: hEff_sum ! Summed effective face thicknesses [H ~> m or kg m-2]
373		integer :: iMethod
374	0	real, dimension(SZI_(G)) :: ref_pres ! Reference pressure used to calculate alpha/beta [R L2 T-2 ~> Pa]
375		real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2]
376	0	integer, dimension(SZI_(G), SZJ_(G)) :: k_top ! Index of the first layer within the boundary
377	0	real, dimension(SZI_(G), SZJ_(G)) :: zeta_top ! Distance from the top of a layer to the intersection of the
378		! top extent of the boundary layer (0 at top, 1 at bottom) [nondim]
379	0	integer, dimension(SZI_(G), SZJ_(G)) :: k_bot ! Index of the last layer within the boundary
380	0	real, dimension(SZI_(G), SZJ_(G)) :: zeta_bot ! Distance of the lower layer to the boundary layer depth [nondim]
381		real :: pa_to_H ! A conversion factor from rescaled pressure to thickness
382		! (H) units [H T2 R-1 Z-2 ~> m Pa-1 or s2 m-1]
383
384	0	pa_to_H = 1. / (GV%H_to_RZ * GV%g_Earth)
385
386	0	k_top(:,:) = 1 ; k_bot(:,:) = 1
387	0	zeta_top(:,:) = 0. ; zeta_bot(:,:) = 0.
388
389		! Check if hbl needs to be extracted
390	0	if (CS%interior_only) then
391	0	if (associated(visc%h_ML)) then
392	0	CS%hbl(:,:) = visc%h_ML(:,:)
393		else
394	0	call MOM_error(FATAL, "hor_bnd_diffusion requires that visc%h_ML is associated.")
395		endif
396	0	call pass_var(CS%hbl, G%Domain, halo=1)
397
398		! get k-indices and zeta
399	0	do j=G%jsc-1, G%jec+1 ; do i=G%isc-1,G%iec+1
400	0	if (G%mask2dT(i,j) > 0.0) then
401		call boundary_k_range(SURFACE, G%ke, h(i,j,:), CS%hbl(i,j), k_top(i,j), zeta_top(i,j), k_bot(i,j), &
402	0	zeta_bot(i,j))
403		endif
404		enddo ; enddo
405		! TODO: add similar code for BOTTOM boundary layer
406		endif
407
408	0	h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff
409
410	0	if (.not. CS%continuous_reconstruction) then
411	0	if (CS%remap_answer_date < 20190101) then
412	0	if (GV%Boussinesq) then
413	0	h_neglect = GV%m_to_H1.0e-30 ; h_neglect_edge = GV%m_to_H1.0e-10
414		else
415	0	h_neglect = GV%kg_m2_to_H1.0e-30 ; h_neglect_edge = GV%kg_m2_to_H1.0e-10
416		endif
417		endif
418		endif
419
420		! If doing along isopycnal diffusion (as opposed to neutral diffusion, set the reference pressure)
421	0	if (CS%ref_pres>=0.) then
422	0	ref_pres(:) = CS%ref_pres
423		endif
424
425	0	if (CS%continuous_reconstruction) then
426	0	CS%dRdT(:,:,:) = 0.
427	0	CS%dRdS(:,:,:) = 0.
428		else
429	0	CS%T_i(:,:,:,:) = 0.
430	0	CS%S_i(:,:,:,:) = 0.
431	0	CS%dRdT_i(:,:,:,:) = 0.
432	0	CS%dRdS_i(:,:,:,:) = 0.
433	0	CS%ns(:,:) = 0.
434	0	CS%stable_cell(:,:,:) = .true.
435		endif
436
437		! Calculate pressure at interfaces and layer averaged alpha/beta
438	0	if (present(p_surf)) then
439	0	do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1
440	0	CS%Pint(i,j,1) = p_surf(i,j)
441		enddo ; enddo
442		else
443	0	CS%Pint(:,:,1) = 0.
444		endif
445	0	do k=1,GV%ke ; do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1
446	0	CS%Pint(i,j,k+1) = CS%Pint(i,j,k) + h(i,j,k)(GV%g_EarthGV%H_to_RZ)
447		enddo ; enddo ; enddo
448
449		! Pressures at the interfaces, this is redundant as P_i(k,1) = P_i(k-1,2) however retain this
450		! for now to ensure consistency of indexing for discontinuous reconstructions
451	0	if (.not. CS%continuous_reconstruction) then
452	0	if (present(p_surf)) then
453	0	do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1
454	0	CS%P_i(i,j,1,1) = p_surf(i,j)
455	0	CS%P_i(i,j,1,2) = p_surf(i,j) + h(i,j,1)(GV%H_to_RZGV%g_Earth)
456		enddo ; enddo
457		else
458	0	do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1
459	0	CS%P_i(i,j,1,1) = 0.
460	0	CS%P_i(i,j,1,2) = h(i,j,1)(GV%H_to_RZGV%g_Earth)
461		enddo ; enddo
462		endif
463	0	do k=2,GV%ke ; do j=G%jsc-1,G%jec+1 ; do i=G%isc-1,G%iec+1
464	0	CS%P_i(i,j,k,1) = CS%P_i(i,j,k-1,2)
465	0	CS%P_i(i,j,k,2) = CS%P_i(i,j,k-1,2) + h(i,j,k)(GV%H_to_RZGV%g_Earth)
466		enddo ; enddo ; enddo
467		endif
468
469	0	EOSdom(:) = EOS_domain(G%HI, halo=1)
470	0	do j = G%jsc-1, G%jec+1
471		! Interpolate state to interface
472	0	do i = G%isc-1, G%iec+1
473	0	if (CS%continuous_reconstruction) then
474	0	call interface_scalar(GV%ke, h(i,j,:), T(i,j,:), CS%Tint(i,j,:), 2, h_neglect)
475	0	call interface_scalar(GV%ke, h(i,j,:), S(i,j,:), CS%Sint(i,j,:), 2, h_neglect)
476		else
477		call build_reconstructions_1d( CS%remap_CS, GV%ke, h(i,j,:), T(i,j,:), CS%ppoly_coeffs_T(i,j,:,:), &
478	0	CS%T_i(i,j,:,:), ppoly_r_S, iMethod, h_neglect, h_neglect_edge )
479		call build_reconstructions_1d( CS%remap_CS, GV%ke, h(i,j,:), S(i,j,:), CS%ppoly_coeffs_S(i,j,:,:), &
480	0	CS%S_i(i,j,:,:), ppoly_r_S, iMethod, h_neglect, h_neglect_edge )
481		! In the current ALE formulation, interface values are not exactly at the 0. or 1. of the
482		! polynomial reconstructions
483	0	do k=1,GV%ke
484	0	CS%T_i(i,j,k,1) = evaluation_polynomial( CS%ppoly_coeffs_T(i,j,k,:), CS%deg+1, 0. )
485	0	CS%T_i(i,j,k,2) = evaluation_polynomial( CS%ppoly_coeffs_T(i,j,k,:), CS%deg+1, 1. )
486	0	CS%S_i(i,j,k,1) = evaluation_polynomial( CS%ppoly_coeffs_S(i,j,k,:), CS%deg+1, 0. )
487	0	CS%S_i(i,j,k,2) = evaluation_polynomial( CS%ppoly_coeffs_S(i,j,k,:), CS%deg+1, 1. )
488		enddo
489		endif
490		enddo
491
492		! Continuous reconstruction
493	0	if (CS%continuous_reconstruction) then
494	0	do k = 1, GV%ke+1
495	0	if (CS%ref_pres<0) ref_pres(:) = CS%Pint(:,j,k)
496		call calculate_density_derivs(CS%Tint(:,j,k), CS%Sint(:,j,k), ref_pres, CS%dRdT(:,j,k), &
497	0	CS%dRdS(:,j,k), CS%EOS, EOSdom)
498		enddo
499		else ! Discontinuous reconstruction
500	0	do k = 1, GV%ke
501	0	if (CS%ref_pres<0) ref_pres(:) = CS%Pint(:,j,k)
502		! Calculate derivatives for the top interface
503		call calculate_density_derivs(CS%T_i(:,j,k,1), CS%S_i(:,j,k,1), ref_pres, CS%dRdT_i(:,j,k,1), &
504	0	CS%dRdS_i(:,j,k,1), CS%EOS, EOSdom)
505	0	if (CS%ref_pres<0) ref_pres(:) = CS%Pint(:,j,k+1)
506		! Calculate derivatives at the bottom interface
507		call calculate_density_derivs(CS%T_i(:,j,k,2), CS%S_i(:,j,k,2), ref_pres, CS%dRdT_i(:,j,k,2), &
508	0	CS%dRdS_i(:,j,k,2), CS%EOS, EOSdom)
509		enddo
510		endif
511		enddo
512
513	0	if (.not. CS%continuous_reconstruction) then
514	0	do j = G%jsc-1, G%jec+1 ; do i = G%isc-1, G%iec+1
515	0	call mark_unstable_cells( CS, GV%ke, CS%T_i(i,j,:,:), CS%S_i(i,j,:,:), CS%P_i(i,j,:,:), CS%stable_cell(i,j,:) )
516	0	if (CS%interior_only) then
517	0	if (.not. CS%stable_cell(i,j,k_bot(i,j))) zeta_bot(i,j) = -1.
518		! set values in the surface and bottom boundary layer to false.
519	0	do k = 1, k_bot(i,j)
520	0	CS%stable_cell(i,j,k) = .false.
521		enddo
522		endif
523		enddo ; enddo
524		endif
525
526	0	CS%uhEff(:,:,:) = 0.
527	0	CS%vhEff(:,:,:) = 0.
528	0	CS%uPoL(:,:,:) = 0.
529	0	CS%vPoL(:,:,:) = 0.
530	0	CS%uPoR(:,:,:) = 0.
531	0	CS%vPoR(:,:,:) = 0.
532	0	CS%uKoL(:,:,:) = 1
533	0	CS%vKoL(:,:,:) = 1
534	0	CS%uKoR(:,:,:) = 1
535	0	CS%vKoR(:,:,:) = 1
536
537		! Neutral surface factors at U points
538	0	do j = G%jsc, G%jec ; do I = G%isc-1, G%iec
539	0	if (G%mask2dCu(I,j) > 0.0) then
540	0	if (CS%continuous_reconstruction) then
541		call find_neutral_surface_positions_continuous(GV%ke, &
542		CS%Pint(i,j,:), CS%Tint(i,j,:), CS%Sint(i,j,:), CS%dRdT(i,j,:), CS%dRdS(i,j,:), &
543		CS%Pint(i+1,j,:), CS%Tint(i+1,j,:), CS%Sint(i+1,j,:), CS%dRdT(i+1,j,:), CS%dRdS(i+1,j,:), &
544		CS%uPoL(I,j,:), CS%uPoR(I,j,:), CS%uKoL(I,j,:), CS%uKoR(I,j,:), CS%uhEff(I,j,:), &
545	0	k_bot(I,j), k_bot(I+1,j), zeta_bot(I,j), zeta_bot(I+1,j))
546		else
547		call find_neutral_surface_positions_discontinuous(CS, GV%ke, &
548		CS%P_i(i,j,:,:), h(i,j,:), CS%T_i(i,j,:,:), CS%S_i(i,j,:,:), CS%ppoly_coeffs_T(i,j,:,:), &
549		CS%ppoly_coeffs_S(i,j,:,:),CS%stable_cell(i,j,:), &
550		CS%P_i(i+1,j,:,:), h(i+1,j,:), CS%T_i(i+1,j,:,:), CS%S_i(i+1,j,:,:), CS%ppoly_coeffs_T(i+1,j,:,:), &
551		CS%ppoly_coeffs_S(i+1,j,:,:), CS%stable_cell(i+1,j,:), &
552		CS%uPoL(I,j,:), CS%uPoR(I,j,:), CS%uKoL(I,j,:), CS%uKoR(I,j,:), CS%uhEff(I,j,:), &
553	0	hard_fail_heff = CS%hard_fail_heff)
554		endif
555		endif
556		enddo ; enddo
557
558		! Neutral surface factors at V points
559	0	do J = G%jsc-1, G%jec ; do i = G%isc, G%iec
560	0	if (G%mask2dCv(i,J) > 0.0) then
561	0	if (CS%continuous_reconstruction) then
562		call find_neutral_surface_positions_continuous(GV%ke, &
563		CS%Pint(i,j,:), CS%Tint(i,j,:), CS%Sint(i,j,:), CS%dRdT(i,j,:), CS%dRdS(i,j,:), &
564		CS%Pint(i,j+1,:), CS%Tint(i,j+1,:), CS%Sint(i,j+1,:), CS%dRdT(i,j+1,:), CS%dRdS(i,j+1,:), &
565		CS%vPoL(i,J,:), CS%vPoR(i,J,:), CS%vKoL(i,J,:), CS%vKoR(i,J,:), CS%vhEff(i,J,:), &
566	0	k_bot(i,J), k_bot(i,J+1), zeta_bot(i,J), zeta_bot(i,J+1))
567		else
568		call find_neutral_surface_positions_discontinuous(CS, GV%ke, &
569		CS%P_i(i,j,:,:), h(i,j,:), CS%T_i(i,j,:,:), CS%S_i(i,j,:,:), CS%ppoly_coeffs_T(i,j,:,:), &
570		CS%ppoly_coeffs_S(i,j,:,:),CS%stable_cell(i,j,:), &
571		CS%P_i(i,j+1,:,:), h(i,j+1,:), CS%T_i(i,j+1,:,:), CS%S_i(i,j+1,:,:), CS%ppoly_coeffs_T(i,j+1,:,:), &
572		CS%ppoly_coeffs_S(i,j+1,:,:), CS%stable_cell(i,j+1,:), &
573		CS%vPoL(I,j,:), CS%vPoR(I,j,:), CS%vKoL(I,j,:), CS%vKoR(I,j,:), CS%vhEff(I,j,:), &
574	0	hard_fail_heff = CS%hard_fail_heff)
575		endif
576		endif
577		enddo ; enddo
578
579		! Continuous reconstructions calculate hEff as the difference between the pressures of the
580		! neutral surfaces which need to be reconverted to thickness units. The discontinuous version
581		! calculates hEff from the nondimensional fraction of the layer spanned by adjacent neutral
582		! surfaces, so hEff is already in thickness units.
583	0	if (CS%continuous_reconstruction) then
584	0	if (CS%use_unmasked_transport_bug) then
585		! This option is not recommended but needed to recover answers prior to Jan 2018.
586		! It is independent of the other 2018 answers flags.
587	0	do k = 1, CS%nsurf-1 ; do j = G%jsc, G%jec ; do I = G%isc-1, G%iec
588	0	CS%uhEff(I,j,k) = CS%uhEff(I,j,k) / GV%H_to_pa
589		enddo ; enddo ; enddo
590	0	do k = 1, CS%nsurf-1 ; do J = G%jsc-1, G%jec ; do i = G%isc, G%iec
591	0	CS%vhEff(I,j,k) = CS%vhEff(I,j,k) / GV%H_to_pa
592		enddo ; enddo ; enddo
593		else
594	0	do k = 1, CS%nsurf-1 ; do j = G%jsc, G%jec ; do I = G%isc-1, G%iec
595	0	if (G%mask2dCu(I,j) > 0.0) CS%uhEff(I,j,k) = CS%uhEff(I,j,k) * pa_to_H
596		enddo ; enddo ; enddo
597	0	do k = 1, CS%nsurf-1 ; do J = G%jsc-1, G%jec ; do i = G%isc, G%iec
598	0	if (G%mask2dCv(i,J) > 0.0) CS%vhEff(i,J,k) = CS%vhEff(i,J,k) * pa_to_H
599		enddo ; enddo ; enddo
600		endif
601		endif
602
603	0	if (CS%id_uhEff_2d>0) then
604	0	hEff_sum(:,:) = 0.
605	0	do k = 1,CS%nsurf-1 ; do j=G%jsc,G%jec ; do i=G%isc-1,G%iec
606	0	hEff_sum(i,j) = hEff_sum(i,j) + CS%uhEff(i,j,k)
607		enddo ; enddo ; enddo
608	0	call post_data(CS%id_uhEff_2d, hEff_sum, CS%diag)
609		endif
610	0	if (CS%id_vhEff_2d>0) then
611	0	hEff_sum(:,:) = 0.
612	0	do k = 1,CS%nsurf-1 ; do j=G%jsc-1,G%jec ; do i=G%isc,G%iec
613	0	hEff_sum(i,j) = hEff_sum(i,j) + CS%vhEff(i,j,k)
614		enddo ; enddo ; enddo
615	0	call post_data(CS%id_vhEff_2d, hEff_sum, CS%diag)
616		endif
617
618	0	end subroutine neutral_diffusion_calc_coeffs
619
620		!> Update tracer concentration due to neutral diffusion; layer thickness unchanged by this update.
621	0	subroutine neutral_diffusion(G, GV, h, Coef_x, Coef_y, dt, Reg, US, CS)
622		type(ocean_grid_type), intent(in) :: G !< Ocean grid structure
623		type(verticalGrid_type), intent(in) :: GV !< ocean vertical grid structure
624		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
625		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: Coef_x !< dt * Kh * dy / dx at u-points [L2 ~> m2]
626		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)+1), intent(in) :: Coef_y !< dt * Kh * dx / dy at v-points [L2 ~> m2]
627		real, intent(in) :: dt !< Tracer time step * I_numitts [T ~> s]
628		!! (I_numitts is in tracer_hordiff)
629		type(tracer_registry_type), pointer :: Reg !< Tracer registry
630		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
631		type(neutral_diffusion_CS), pointer :: CS !< Neutral diffusion control structure
632
633		! Local variables
634	0	real, dimension(SZIB_(G),SZJ_(G),CS%nsurf-1) :: uFlx ! Zonal flux of tracer in units that vary between a
635		! thickness times a concentration ([C H ~> degC m or degC kg m-2] for temperature) or a
636		! volume or mass times a concentration ([C H L2 ~> degC m3 or degC kg] for temperature),
637		! depending on the setting of CS%KhTh_use_vert_struct.
638	0	real, dimension(SZI_(G),SZJB_(G),CS%nsurf-1) :: vFlx ! Meridional flux of tracer in units that vary between a
639		! thickness times a concentration ([C H ~> degC m or degC kg m-2] for temperature) or a
640		! volume or mass times a concentration ([C H L2 ~> degC m3 or degC kg] for temperature),
641		! depending on the setting of CS%KhTh_use_vert_struct.
642	0	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: tendency ! tendency array for diagnostics
643		! [H conc T-1 ~> m conc s-1 or kg m-2 conc s-1]
644		! For temperature these units are
645		! [C H T-1 ~> degC m s-1 or degC kg m-2 s-1].
646	0	real, dimension(SZI_(G),SZJ_(G)) :: tendency_2d ! Depth integrated content tendency for diagnostics
647		! [H conc T-1 ~> m conc s-1 or kg m-2 conc s-1].
648		! For temperature these units are
649		! [C H T-1 ~> degC m s-1 or degC kg m-2 s-1].
650	0	real, dimension(SZIB_(G),SZJ_(G)) :: trans_x_2d ! Depth integrated diffusive tracer x-transport
651		! diagnostic. For temperature this has units of
652		! [C H L2 ~> degC m3 or degC kg].
653	0	real, dimension(SZI_(G),SZJB_(G)) :: trans_y_2d ! depth integrated diffusive tracer y-transport
654		! diagnostic. For temperature this has units of
655		! [C H L2 ~> degC m3 or degC kg].
656	0	real, dimension(SZK_(GV)) :: dTracer ! Change in tracer concentration due to neutral diffusion
657		! [H L2 conc ~> m3 conc or kg conc]. For temperature
658		! these units are [C H L2 ~> degC m3 or degC kg].
659	0	real, dimension(SZK_(GV)) :: dTracer_N ! Change in tracer concentration due to neutral diffusion
660		! into a cell via its logically northern face, in
661		! [H L2 conc ~> m3 conc or kg conc].
662	0	real, dimension(SZK_(GV)) :: dTracer_S ! Change in tracer concentration due to neutral diffusion
663		! into a cell via its logically southern face, in
664		! [H L2 conc ~> m3 conc or kg conc].
665	0	real, dimension(SZK_(GV)) :: dTracer_E ! Change in tracer concentration due to neutral diffusion
666		! into a cell via its logically eastern face, in
667		! [H L2 conc ~> m3 conc or kg conc].
668	0	real, dimension(SZK_(GV)) :: dTracer_W ! Change in tracer concentration due to neutral diffusion
669		! into a cell via its logically western face, in
670		! [H L2 conc ~> m3 conc or kg conc].
671		real :: normalize ! normalization used for averaging Coef_x and Coef_y to t-points [nondim].
672
673		type(tracer_type), pointer :: Tracer => NULL() ! Pointer to the current tracer
674
675		integer :: i, j, k, m, ks, nk
676		real :: Idt ! The inverse of the time step [T-1 ~> s-1]
677		real :: h_neglect, h_neglect_edge ! Negligible thicknesses [H ~> m or kg m-2]
678	0	h_neglect = GV%H_subroundoff ; h_neglect_edge = GV%H_subroundoff
679
680	0	if (.not. CS%continuous_reconstruction) then
681	0	if (CS%remap_answer_date < 20190101) then
682	0	h_neglect = GV%m_to_H1.0e-30 ; h_neglect_edge = GV%m_to_H1.0e-10
683		endif
684		endif
685
686	0	if (CS%KhTh_use_vert_struct) then
687		! Compute Coef at h points
688	0	CS%Coef_h(:,:,:) = 0.
689	0	do j = G%jsc,G%jec ; do i = G%isc,G%iec
690	0	if (G%mask2dT(i,j)>0.) then
691		normalize = 1.0 / ((G%mask2dCu(I-1,j)+G%mask2dCu(I,j)) + &
692	0	(G%mask2dCv(i,J-1)+G%mask2dCv(i,J)) + 1.0e-37)
693	0	do k = 1, GV%ke+1
694		CS%Coef_h(i,j,k) = normalizeG%mask2dT(i,j)((Coef_x(I-1,j,k)+Coef_x(I,j,k)) + &
695	0	(Coef_y(i,J-1,k)+Coef_y(i,J,k)))
696		enddo
697		endif
698		enddo ; enddo
699	0	call pass_var(CS%Coef_h,G%Domain)
700		endif
701
702	0	nk = GV%ke
703
704	0	do m = 1,Reg%ntr ! Loop over tracer registry
705
706	0	tracer => Reg%Tr(m)
707
708		! for diagnostics
709		if (tracer%id_dfxy_conc > 0 .or. tracer%id_dfxy_cont > 0 .or. tracer%id_dfxy_cont_2d > 0 .or. &
710	0	tracer%id_dfx_2d > 0 .or. tracer%id_dfy_2d > 0) then
711	0	Idt = 1.0 / dt
712	0	tendency(:,:,:) = 0.0
713		endif
714
715	0	uFlx(:,:,:) = 0.
716	0	vFlx(:,:,:) = 0.
717
718		! x-flux
719	0	if (CS%KhTh_use_vert_struct) then
720	0	if (CS%tapering) then
721	0	do j = G%jsc,G%jec ; do I = G%isc-1,G%iec
722	0	if (G%mask2dCu(I,j)>0.) then
723		! compute coeff_l and coeff_r and pass them to neutral_surface_flux
724		call compute_tapering_coeffs(G%ke+1, CS%hbl(I,j), CS%hbl(I+1,j), CS%coeff_l(:), CS%coeff_r(:), &
725	0	h(I,j,:), h(I+1,j,:))
726		call neutral_surface_flux(nk, CS%nsurf, CS%deg, h(i,j,:), h(i+1,j,:), &
727		tracer%t(i,j,:), tracer%t(i+1,j,:), &
728		CS%uPoL(I,j,:), CS%uPoR(I,j,:), &
729		CS%uKoL(I,j,:), CS%uKoR(I,j,:), &
730		CS%uhEff(I,j,:), uFlx(I,j,:), &
731		CS%continuous_reconstruction, h_neglect, &
732		CS%remap_CS, h_neglect_edge, CS%coeff_l(:)*CS%Coef_h(i,j,:), &
733	0	CS%coeff_r(:)*CS%Coef_h(i+1,j,:))
734		endif
735		enddo ; enddo
736		else
737	0	do j = G%jsc,G%jec ; do I = G%isc-1,G%iec
738	0	if (G%mask2dCu(I,j)>0.) then
739		call neutral_surface_flux(nk, CS%nsurf, CS%deg, h(i,j,:), h(i+1,j,:), &
740		tracer%t(i,j,:), tracer%t(i+1,j,:), &
741		CS%uPoL(I,j,:), CS%uPoR(I,j,:), &
742		CS%uKoL(I,j,:), CS%uKoR(I,j,:), &
743		CS%uhEff(I,j,:), uFlx(I,j,:), &
744		CS%continuous_reconstruction, h_neglect, &
745		CS%remap_CS, h_neglect_edge, CS%Coef_h(i,j,:), &
746	0	CS%Coef_h(i+1,j,:))
747		endif
748		enddo ; enddo
749		endif
750		else
751	0	if (CS%tapering) then
752	0	do j = G%jsc,G%jec ; do I = G%isc-1,G%iec
753	0	if (G%mask2dCu(I,j)>0.) then
754		! compute coeff_l and coeff_r and pass them to neutral_surface_flux
755		call compute_tapering_coeffs(G%ke+1, CS%hbl(I,j), CS%hbl(I+1,j), CS%coeff_l(:), CS%coeff_r(:), &
756	0	h(I,j,:), h(I+1,j,:))
757		call neutral_surface_flux(nk, CS%nsurf, CS%deg, h(i,j,:), h(i+1,j,:), &
758		tracer%t(i,j,:), tracer%t(i+1,j,:), &
759		CS%uPoL(I,j,:), CS%uPoR(I,j,:), &
760		CS%uKoL(I,j,:), CS%uKoR(I,j,:), &
761		CS%uhEff(I,j,:), uFlx(I,j,:), &
762		CS%continuous_reconstruction, h_neglect, &
763		CS%remap_CS, h_neglect_edge, CS%coeff_l(:), &
764	0	CS%coeff_r(:))
765		endif
766		enddo ; enddo
767		else
768	0	do j = G%jsc,G%jec ; do I = G%isc-1,G%iec
769	0	if (G%mask2dCu(I,j)>0.) then
770		call neutral_surface_flux(nk, CS%nsurf, CS%deg, h(i,j,:), h(i+1,j,:), &
771		tracer%t(i,j,:), tracer%t(i+1,j,:), &
772		CS%uPoL(I,j,:), CS%uPoR(I,j,:), &
773		CS%uKoL(I,j,:), CS%uKoR(I,j,:), &
774		CS%uhEff(I,j,:), uFlx(I,j,:), &
775		CS%continuous_reconstruction, h_neglect, &
776	0	CS%remap_CS, h_neglect_edge)
777		endif
778		enddo ; enddo
779		endif
780		endif
781
782		! y-flux
783	0	if (CS%KhTh_use_vert_struct) then
784	0	if (CS%tapering) then
785	0	do J = G%jsc-1,G%jec ; do i = G%isc,G%iec
786	0	if (G%mask2dCv(i,J)>0.) then
787		! compute coeff_l and coeff_r and pass them to neutral_surface_flux
788		call compute_tapering_coeffs(G%ke+1, CS%hbl(i,J), CS%hbl(i,J+1), CS%coeff_l(:), CS%coeff_r(:), &
789	0	h(i,J,:), h(i,J+1,:))
790
791		call neutral_surface_flux(nk, CS%nsurf, CS%deg, h(i,j,:), h(i,j+1,:), &
792		tracer%t(i,j,:), tracer%t(i,j+1,:), &
793		CS%vPoL(i,J,:), CS%vPoR(i,J,:), &
794		CS%vKoL(i,J,:), CS%vKoR(i,J,:), &
795		CS%vhEff(i,J,:), vFlx(i,J,:), &
796		CS%continuous_reconstruction, h_neglect, &
797		CS%remap_CS, h_neglect_edge, CS%coeff_l(:)*CS%Coef_h(i,j,:), &
798	0	CS%coeff_r(:)*CS%Coef_h(i,j+1,:))
799		endif
800		enddo ; enddo
801		else
802	0	do J = G%jsc-1,G%jec ; do i = G%isc,G%iec
803	0	if (G%mask2dCv(i,J)>0.) then
804		call neutral_surface_flux(nk, CS%nsurf, CS%deg, h(i,j,:), h(i,j+1,:), &
805		tracer%t(i,j,:), tracer%t(i,j+1,:), &
806		CS%vPoL(i,J,:), CS%vPoR(i,J,:), &
807		CS%vKoL(i,J,:), CS%vKoR(i,J,:), &
808		CS%vhEff(i,J,:), vFlx(i,J,:), &
809		CS%continuous_reconstruction, h_neglect, &
810		CS%remap_CS, h_neglect_edge, CS%Coef_h(i,j,:), &
811	0	CS%Coef_h(i,j+1,:))
812		endif
813		enddo ; enddo
814		endif
815		else
816	0	if (CS%tapering) then
817	0	do J = G%jsc-1,G%jec ; do i = G%isc,G%iec
818	0	if (G%mask2dCv(i,J)>0.) then
819		! compute coeff_l and coeff_r and pass them to neutral_surface_flux
820		call compute_tapering_coeffs(G%ke+1, CS%hbl(i,J), CS%hbl(i,J+1), CS%coeff_l(:), CS%coeff_r(:), &
821	0	h(i,J,:), h(i,J+1,:))
822
823		call neutral_surface_flux(nk, CS%nsurf, CS%deg, h(i,j,:), h(i,j+1,:), &
824		tracer%t(i,j,:), tracer%t(i,j+1,:), &
825		CS%vPoL(i,J,:), CS%vPoR(i,J,:), &
826		CS%vKoL(i,J,:), CS%vKoR(i,J,:), &
827		CS%vhEff(i,J,:), vFlx(i,J,:), &
828		CS%continuous_reconstruction, h_neglect, &
829		CS%remap_CS, h_neglect_edge, CS%coeff_l(:), &
830	0	CS%coeff_r(:))
831		endif
832		enddo ; enddo
833		else
834	0	do J = G%jsc-1,G%jec ; do i = G%isc,G%iec
835	0	if (G%mask2dCv(i,J)>0.) then
836		call neutral_surface_flux(nk, CS%nsurf, CS%deg, h(i,j,:), h(i,j+1,:), &
837		tracer%t(i,j,:), tracer%t(i,j+1,:), &
838		CS%vPoL(i,J,:), CS%vPoR(i,J,:), &
839		CS%vKoL(i,J,:), CS%vKoR(i,J,:), &
840		CS%vhEff(i,J,:), vFlx(i,J,:), &
841		CS%continuous_reconstruction, h_neglect, &
842	0	CS%remap_CS, h_neglect_edge)
843		endif
844		enddo ; enddo
845		endif
846		endif
847
848		! Update the tracer concentration from divergence of neutral diffusive flux components, noting
849		! that uFlx and vFlx use an unexpected sign convention.
850	0	if (CS%KhTh_use_vert_struct) then
851	0	do j = G%jsc,G%jec ; do i = G%isc,G%iec
852	0	if (G%mask2dT(i,j)>0.) then
853	0	if (CS%ndiff_answer_date <= 20240330) then
854	0	dTracer(:) = 0.
855	0	do ks = 1,CS%nsurf-1
856	0	k = CS%uKoL(I,j,ks)
857	0	dTracer(k) = dTracer(k) + uFlx(I,j,ks)
858	0	k = CS%uKoR(I-1,j,ks)
859	0	dTracer(k) = dTracer(k) - uFlx(I-1,j,ks)
860	0	k = CS%vKoL(i,J,ks)
861	0	dTracer(k) = dTracer(k) + vFlx(i,J,ks)
862	0	k = CS%vKoR(i,J-1,ks)
863	0	dTracer(k) = dTracer(k) - vFlx(i,J-1,ks)
864		enddo
865		else ! This form recovers rotational symmetry.
866	0	dTracer_N(:) = 0.0 ; dTracer_S(:) = 0.0 ; dTracer_E(:) = 0.0 ; dTracer_W(:) = 0.0
867	0	do ks = 1,CS%nsurf-1
868	0	k = CS%uKoL(I,j,ks)
869	0	dTracer_E(k) = dTracer_E(k) + uFlx(I,j,ks)
870	0	k = CS%uKoR(I-1,j,ks)
871	0	dTracer_W(k) = dTracer_W(k) - uFlx(I-1,j,ks)
872	0	k = CS%vKoL(i,J,ks)
873	0	dTracer_N(k) = dTracer_N(k) + vFlx(i,J,ks)
874	0	k = CS%vKoR(i,J-1,ks)
875	0	dTracer_S(k) = dTracer_S(k) - vFlx(i,J-1,ks)
876		enddo
877	0	do k = 1, GV%ke
878	0	dTracer(k) = (dTracer_N(k) + dTracer_S(k)) + (dTracer_E(k) + dTracer_W(k))
879		enddo
880		endif
881	0	do k = 1, GV%ke
882		tracer%t(i,j,k) = tracer%t(i,j,k) + dTracer(k) * &
883	0	( G%IareaT(i,j) / ( h(i,j,k) + GV%H_subroundoff ) )
884	0	if (abs(tracer%t(i,j,k)) < tracer%conc_underflow) tracer%t(i,j,k) = 0.0
885		enddo
886
887	0	if (tracer%id_dfxy_conc > 0 .or. tracer%id_dfxy_cont > 0 .or. tracer%id_dfxy_cont_2d > 0 ) then
888	0	do k = 1, GV%ke
889	0	tendency(i,j,k) = dTracer(k) * G%IareaT(i,j) * Idt
890		enddo
891		endif
892
893		endif
894		enddo ; enddo
895		else
896	0	do j = G%jsc,G%jec ; do i = G%isc,G%iec
897	0	if (G%mask2dT(i,j)>0.) then
898	0	if (CS%ndiff_answer_date <= 20240330) then
899	0	dTracer(:) = 0.
900	0	do ks = 1,CS%nsurf-1
901	0	k = CS%uKoL(I,j,ks)
902	0	dTracer(k) = dTracer(k) + Coef_x(I,j,1) * uFlx(I,j,ks)
903	0	k = CS%uKoR(I-1,j,ks)
904	0	dTracer(k) = dTracer(k) - Coef_x(I-1,j,1) * uFlx(I-1,j,ks)
905	0	k = CS%vKoL(i,J,ks)
906	0	dTracer(k) = dTracer(k) + Coef_y(i,J,1) * vFlx(i,J,ks)
907	0	k = CS%vKoR(i,J-1,ks)
908	0	dTracer(k) = dTracer(k) - Coef_y(i,J-1,1) * vFlx(i,J-1,ks)
909		enddo
910		else ! This form recovers rotational symmetry.
911	0	dTracer_N(:) = 0.0 ; dTracer_S(:) = 0.0 ; dTracer_E(:) = 0.0 ; dTracer_W(:) = 0.0
912	0	do ks = 1,CS%nsurf-1
913	0	k = CS%uKoL(I,j,ks)
914	0	dTracer_E(k) = dTracer_E(k) + Coef_x(I,j,1) * uFlx(I,j,ks)
915	0	k = CS%uKoR(I-1,j,ks)
916	0	dTracer_W(k) = dTracer_W(k) - Coef_x(I-1,j,1) * uFlx(I-1,j,ks)
917	0	k = CS%vKoL(i,J,ks)
918	0	dTracer_N(k) = dTracer_N(k) + Coef_y(i,J,1) * vFlx(i,J,ks)
919	0	k = CS%vKoR(i,J-1,ks)
920	0	dTracer_S(k) = dTracer_S(k) - Coef_y(i,J-1,1) * vFlx(i,J-1,ks)
921		enddo
922	0	do k = 1, GV%ke
923	0	dTracer(k) = (dTracer_N(k) + dTracer_S(k)) + (dTracer_E(k) + dTracer_W(k))
924		enddo
925		endif
926	0	do k = 1, GV%ke
927		tracer%t(i,j,k) = tracer%t(i,j,k) + dTracer(k) * &
928	0	( G%IareaT(i,j) / ( h(i,j,k) + GV%H_subroundoff ) )
929	0	if (abs(tracer%t(i,j,k)) < tracer%conc_underflow) tracer%t(i,j,k) = 0.0
930		enddo
931
932	0	if (tracer%id_dfxy_conc > 0 .or. tracer%id_dfxy_cont > 0 .or. tracer%id_dfxy_cont_2d > 0 ) then
933	0	do k = 1, GV%ke
934	0	tendency(i,j,k) = dTracer(k) * G%IareaT(i,j) * Idt
935		enddo
936		endif
937
938		endif
939		enddo ; enddo
940		endif
941
942		! Do user controlled underflow of the tracer concentrations.
943	0	if (tracer%conc_underflow > 0.0) then
944	0	do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec
945	0	if (abs(tracer%t(i,j,k)) < tracer%conc_underflow) tracer%t(i,j,k) = 0.0
946		enddo ; enddo ; enddo
947		endif
948
949		! Diagnose vertically summed zonal flux, giving zonal tracer transport from ndiff.
950		! Note sign corresponds to downgradient flux convention.
951	0	if (tracer%id_dfx_2d > 0) then
952
953	0	if (CS%KhTh_use_vert_struct) then
954	0	do j = G%jsc,G%jec ; do I = G%isc-1,G%iec
955	0	trans_x_2d(I,j) = 0.
956	0	if (G%mask2dCu(I,j)>0.) then
957	0	do ks = 1,CS%nsurf-1
958	0	trans_x_2d(I,j) = trans_x_2d(I,j) - uFlx(I,j,ks)
959		enddo
960	0	trans_x_2d(I,j) = trans_x_2d(I,j) * Idt
961		endif
962		enddo ; enddo
963		else
964	0	do j = G%jsc,G%jec ; do I = G%isc-1,G%iec
965	0	trans_x_2d(I,j) = 0.
966	0	if (G%mask2dCu(I,j)>0.) then
967	0	do ks = 1,CS%nsurf-1
968	0	trans_x_2d(I,j) = trans_x_2d(I,j) - Coef_x(I,j,1) * uFlx(I,j,ks)
969		enddo
970	0	trans_x_2d(I,j) = trans_x_2d(I,j) * Idt
971		endif
972		enddo ; enddo
973		endif
974
975	0	call post_data(tracer%id_dfx_2d, trans_x_2d(:,:), CS%diag)
976		endif
977
978		! Diagnose vertically summed merid flux, giving meridional tracer transport from ndiff.
979		! Note sign corresponds to downgradient flux convention.
980	0	if (tracer%id_dfy_2d > 0) then
981
982	0	if (CS%KhTh_use_vert_struct) then
983	0	do J = G%jsc-1,G%jec ; do i = G%isc,G%iec
984	0	trans_y_2d(i,J) = 0.
985	0	if (G%mask2dCv(i,J)>0.) then
986	0	do ks = 1,CS%nsurf-1
987	0	trans_y_2d(i,J) = trans_y_2d(i,J) - vFlx(i,J,ks)
988		enddo
989	0	trans_y_2d(i,J) = trans_y_2d(i,J) * Idt
990		endif
991		enddo ; enddo
992		else
993	0	do J = G%jsc-1,G%jec ; do i = G%isc,G%iec
994	0	trans_y_2d(i,J) = 0.
995	0	if (G%mask2dCv(i,J)>0.) then
996	0	do ks = 1,CS%nsurf-1
997	0	trans_y_2d(i,J) = trans_y_2d(i,J) - Coef_y(i,J,1) * vFlx(i,J,ks)
998		enddo
999	0	trans_y_2d(i,J) = trans_y_2d(i,J) * Idt
1000		endif
1001		enddo ; enddo
1002		endif
1003
1004	0	call post_data(tracer%id_dfy_2d, trans_y_2d(:,:), CS%diag)
1005		endif
1006
1007		! post tendency of layer-integrated tracer content
1008	0	if (tracer%id_dfxy_cont > 0) then
1009	0	call post_data(tracer%id_dfxy_cont, tendency(:,:,:), CS%diag)
1010		endif
1011
1012		! post depth summed tendency for tracer content
1013	0	if (tracer%id_dfxy_cont_2d > 0) then
1014	0	tendency_2d(:,:) = 0.
1015	0	do j = G%jsc,G%jec ; do i = G%isc,G%iec
1016	0	do k = 1, GV%ke
1017	0	tendency_2d(i,j) = tendency_2d(i,j) + tendency(i,j,k)
1018		enddo
1019		enddo ; enddo
1020	0	call post_data(tracer%id_dfxy_cont_2d, tendency_2d(:,:), CS%diag)
1021		endif
1022
1023		! post tendency of tracer concentration; this step must be
1024		! done after posting tracer content tendency, since we alter
1025		! the tendency array.
1026	0	if (tracer%id_dfxy_conc > 0) then
1027	0	do k = 1, GV%ke ; do j = G%jsc,G%jec ; do i = G%isc,G%iec
1028	0	tendency(i,j,k) = tendency(i,j,k) / ( h(i,j,k) + GV%H_subroundoff )
1029		enddo ; enddo ; enddo
1030	0	call post_data(tracer%id_dfxy_conc, tendency, CS%diag)
1031		endif
1032		enddo ! Loop over tracer registry
1033
1034	0	end subroutine neutral_diffusion
1035
1036		!> Computes linear tapering coefficients at interfaces of the left and right columns
1037		!! within a region defined by the boundary layer depths in the two columns.
1038	0	subroutine compute_tapering_coeffs(ne, bld_l, bld_r, coeff_l, coeff_r, h_l, h_r)
1039		integer, intent(in) :: ne !< Number of interfaces
1040		real, intent(in) :: bld_l !< Boundary layer depth, left column [H ~> m or kg m-2]
1041		real, intent(in) :: bld_r !< Boundary layer depth, right column [H ~> m or kg m-2]
1042		real, dimension(ne-1), intent(in) :: h_l !< Layer thickness, left column [H ~> m or kg m-2]
1043		real, dimension(ne-1), intent(in) :: h_r !< Layer thickness, right column [H ~> m or kg m-2]
1044		real, dimension(ne), intent(inout) :: coeff_l !< Tapering coefficient, left column [nondim]
1045		real, dimension(ne), intent(inout) :: coeff_r !< Tapering coefficient, right column [nondim]
1046
1047		! Local variables
1048		real :: min_bld ! Minimum of the boundary layer depth in two adjacent columns [H ~> m or kg m-2]
1049		real :: max_bld ! Maximum of the boundary layer depth in two adjacent columns [H ~> m or kg m-2]
1050		integer :: dummy1 ! dummy integer
1051		real :: dummy2 ! dummy real [nondim]
1052		integer :: k_min_l, k_min_r, k_max_l, k_max_r ! Min/max vertical indices in two adjacent columns
1053		real :: zeta_l, zeta_r ! dummy variables [nondim]
1054		integer :: k ! vertical index
1055
1056		! Initialize coefficients
1057	0	coeff_l(:) = 1.0
1058	0	coeff_r(:) = 1.0
1059
1060		! Calculate vertical indices containing the boundary layer depths
1061	0	max_bld = MAX(bld_l, bld_r)
1062	0	min_bld = MIN(bld_l, bld_r)
1063
1064		! k_min
1065		call boundary_k_range(SURFACE, ne-1, h_l, min_bld, dummy1, dummy2, k_min_l, &
1066	0	zeta_l)
1067		call boundary_k_range(SURFACE, ne-1, h_r, min_bld, dummy1, dummy2, k_min_r, &
1068	0	zeta_r)
1069
1070		! k_max
1071		call boundary_k_range(SURFACE, ne-1, h_l, max_bld, dummy1, dummy2, k_max_l, &
1072	0	zeta_l)
1073		call boundary_k_range(SURFACE, ne-1, h_r, max_bld, dummy1, dummy2, k_max_r, &
1074	0	zeta_r)
1075		! left
1076	0	do k=1,k_min_l
1077	0	coeff_l(k) = 0.0
1078		enddo
1079	0	do k=k_min_l+1,k_max_l+1
1080	0	coeff_l(k) = (real(k - k_min_l) + 1.0)/(real(k_max_l - k_min_l) + 2.0)
1081		enddo
1082
1083		! right
1084	0	do k=1,k_min_r
1085	0	coeff_r(k) = 0.0
1086		enddo
1087	0	do k=k_min_r+1,k_max_r+1
1088	0	coeff_r(k) = (real(k - k_min_r) + 1.0)/(real(k_max_r - k_min_r) + 2.0)
1089		enddo
1090
1091	0	end subroutine compute_tapering_coeffs
1092
1093		!> Returns interface scalar, Si, for a column of layer values, S.
1094	0	subroutine interface_scalar(nk, h, S, Si, i_method, h_neglect)
1095		integer, intent(in) :: nk !< Number of levels
1096		real, dimension(nk), intent(in) :: h !< Layer thickness [H ~> m or kg m-2]
1097		real, dimension(nk), intent(in) :: S !< Layer scalar (or concentrations) in arbitrary
1098		!! concentration units (e.g. [C ~> degC] for temperature)
1099		real, dimension(nk+1), intent(inout) :: Si !< Interface scalar (or concentrations) in arbitrary
1100		!! concentration units (e.g. [C ~> degC] for temperature)
1101		integer, intent(in) :: i_method !< =1 use average of PLM edges
1102		!! =2 use continuous PPM edge interpolation
1103		real, intent(in) :: h_neglect !< A negligibly small thickness [H ~> m or kg m-2]
1104		! Local variables
1105		integer :: k, km2, kp1
1106	0	real, dimension(nk) :: diff ! Difference in scalar concentrations between layer centers in arbitrary
1107		! concentration units (e.g. [C ~> degC] for temperature)
1108		real :: Sb, Sa ! Values of scalar concentrations at the upper and lower edges of a layer in arbitrary
1109		! concentration units (e.g. [C ~> degC] for temperature)
1110
1111	0	call PLM_diff(nk, h, S, 2, 1, diff)
1112	0	Si(1) = S(1) - 0.5 * diff(1)
1113	0	if (i_method==1) then
1114	0	do k = 2, nk
1115		! Average of the two edge values (will be bounded and,
1116		! when slopes are unlimited, notionally second-order accurate)
1117	0	Sa = S(k-1) + 0.5 * diff(k-1) ! Lower edge value of a PLM reconstruction for layer above
1118	0	Sb = S(k) - 0.5 * diff(k) ! Upper edge value of a PLM reconstruction for layer below
1119	0	Si(k) = 0.5 * ( Sa + Sb )
1120		enddo
1121	0	elseif (i_method==2) then
1122	0	do k = 2, nk
1123		! PPM quasi-fourth order interpolation for edge values following
1124		! equation 1.6 in Colella & Woodward, 1984: JCP 54, 174-201.
1125	0	km2 = max(1, k-2)
1126	0	kp1 = min(nk, k+1)
1127	0	Si(k) = ppm_edge(h(km2), h(k-1), h(k), h(kp1), S(k-1), S(k), diff(k-1), diff(k), h_neglect)
1128		enddo
1129		endif
1130	0	Si(nk+1) = S(nk) + 0.5 * diff(nk)
1131
1132	0	end subroutine interface_scalar
1133
1134		!> Returns the PPM quasi-fourth order edge value at k+1/2 following
1135		!! equation 1.6 in Colella & Woodward, 1984: JCP 54, 174-201.
1136		!! The returned units are the same as those of Ak (e.g. [C ~> degC] for temperature).
1137	0	real function ppm_edge(hkm1, hk, hkp1, hkp2, Ak, Akp1, Pk, Pkp1, h_neglect)
1138		real, intent(in) :: hkm1 !< Width of cell k-1 in [H ~> m or kg m-2] or other units
1139		real, intent(in) :: hk !< Width of cell k in [H ~> m or kg m-2] or other units
1140		real, intent(in) :: hkp1 !< Width of cell k+1 in [H ~> m or kg m-2] or other units
1141		real, intent(in) :: hkp2 !< Width of cell k+2 in [H ~> m or kg m-2] or other units
1142		real, intent(in) :: Ak !< Average scalar value of cell k in arbitrary concentration
1143		!! units (e.g. [C ~> degC] for temperature)
1144		real, intent(in) :: Akp1 !< Average scalar value of cell k+1 in arbitrary concentration
1145		!! units (e.g. [C ~> degC] for temperature)
1146		real, intent(in) :: Pk !< PLM slope for cell k in arbitrary concentration
1147		!! units (e.g. [C ~> degC] for temperature)
1148		real, intent(in) :: Pkp1 !< PLM slope for cell k+1 in arbitrary concentration
1149		!! units (e.g. [C ~> degC] for temperature)
1150		real, intent(in) :: h_neglect !< A negligibly small thickness [H ~> m or kg m-2]
1151
1152		! Local variables
1153		real :: R_hk_hkp1, R_2hk_hkp1, R_hk_2hkp1 ! Reciprocals of combinations of thicknesses [H-1 ~> m-1 or m2 kg-1]
1154		real :: f1 ! A work variable with units of an inverse cell width [H-1 ~> m-1 or m2 kg-1]
1155		real :: f2, f3, f4 ! Work variables with units of the cell width [H ~> m or kg m-2]
1156
1157	0	R_hk_hkp1 = hk + hkp1
1158	0	if (R_hk_hkp1 <= 0.) then
1159	0	ppm_edge = 0.5 * ( Ak + Akp1 )
1160	0	return
1161		endif
1162	0	R_hk_hkp1 = 1. / R_hk_hkp1
1163	0	if (hk<hkp1) then
1164	0	ppm_edge = Ak + ( hk * R_hk_hkp1 ) * ( Akp1 - Ak )
1165		else
1166	0	ppm_edge = Akp1 + ( hkp1 * R_hk_hkp1 ) * ( Ak - Akp1 )
1167		endif
1168
1169	0	R_2hk_hkp1 = 1. / ( ( 2. * hk + hkp1 ) + h_neglect )
1170	0	R_hk_2hkp1 = 1. / ( ( hk + 2. * hkp1 ) + h_neglect )
1171	0	f1 = 1./ ( ( hk + hkp1) + ( hkm1 + hkp2 ) )
1172		f2 = 2. * ( hkp1 * hk ) * R_hk_hkp1 * &
1173	0	( ( hkm1 + hk ) * R_2hk_hkp1 - ( hkp2 + hkp1 ) * R_hk_2hkp1 )
1174	0	f3 = hk * ( hkm1 + hk ) * R_2hk_hkp1
1175	0	f4 = hkp1 * ( hkp1 + hkp2 ) * R_hk_2hkp1
1176
1177	0	ppm_edge = ppm_edge + f1 * ( f2 * ( Akp1 - Ak ) - ( f3 * Pkp1 - f4 * Pk ) )
1178
1179	0	end function ppm_edge
1180
1181		!> Returns the average of a PPM reconstruction between two fractional positions in the same
1182		!! arbitrary concentration units as aMean (e.g. usually [C ~> degC] for temperature)
1183	0	real function ppm_ave(xL, xR, aL, aR, aMean)
1184		real, intent(in) :: xL !< Fraction position of left bound (0,1) [nondim]
1185		real, intent(in) :: xR !< Fraction position of right bound (0,1) [nondim]
1186		real, intent(in) :: aL !< Left edge scalar value, at x=0, in arbitrary concentration
1187		!! units (e.g. usually [C ~> degC] for temperature)
1188		real, intent(in) :: aR !< Right edge scalar value, at x=1 in arbitrary concentration
1189		!! units (e.g. usually [C ~> degC] for temperature)
1190		real, intent(in) :: aMean !< Average scalar value of cell in arbitrary concentration
1191		!! units (e.g. usually [C ~> degC] for temperature)
1192
1193		! Local variables
1194		real :: dx ! Distance between the bounds [nondim]
1195		real :: xave ! Average fractional position [nondim]
1196		real :: a6, a6o3 ! Terms proportional to the normalized scalar curvature in the same arbitrary
1197		! concentration units as aMean (e.g. usually [C ~> degC] for temperature)
1198
1199	0	dx = xR - xL
1200	0	xave = 0.5 * ( xR + xL )
1201	0	a6o3 = 2. * aMean - ( aL + aR ) ! a6 / 3.
1202	0	a6 = 3. * a6o3
1203
1204	0	if (dx<0.) then
1205	0	stop 'ppm_ave: dx<0 should not happened!'
1206	0	elseif (dx>1.) then
1207	0	stop 'ppm_ave: dx>1 should not happened!'
1208	0	elseif (dx==0.) then
1209	0	ppm_ave = aL + ( aR - aL ) * xR + a6 * xR * ( 1. - xR )
1210		else
1211	0	ppm_ave = ( aL + xave * ( ( aR - aL ) + a6 ) ) - a6o3 * ( xR*2 + xR xL + xL**2 )
1212		endif
1213	0	end function ppm_ave
1214
1215		!> A true signum function that returns either -abs(a), when x<0; or abs(a) when x>0; or 0 when x=0.
1216		!! The returned units are the same as those of a [arbitrary].
1217	0	real function signum(a,x)
1218		real, intent(in) :: a !< The magnitude argument in arbitrary units [arbitrary]
1219		real, intent(in) :: x !< The sign (or zero) argument [arbitrary]
1220
1221	0	signum = sign(a,x)
1222	0	if (x==0.) signum = 0.
1223
1224	0	end function signum
1225
1226		!> Returns PLM slopes for a column where the slopes are the difference in value across each cell.
1227		!! The limiting follows equation 1.8 in Colella & Woodward, 1984: JCP 54, 174-201.
1228	0	subroutine PLM_diff(nk, h, S, c_method, b_method, diff)
1229		integer, intent(in) :: nk !< Number of levels
1230		real, dimension(nk), intent(in) :: h !< Layer thickness [H ~> m or kg m-2] or other units
1231		real, dimension(nk), intent(in) :: S !< Layer salinity (conc, e.g. ppt) or other tracer
1232		!! concentration in arbitrary units [A ~> a]
1233		integer, intent(in) :: c_method !< Method to use for the centered difference
1234		integer, intent(in) :: b_method !< =1, use PCM in first/last cell, =2 uses linear extrapolation
1235		real, dimension(nk), intent(inout) :: diff !< Scalar difference across layer (conc, e.g. ppt)
1236		!! in the same arbitrary units as S [A ~> a],
1237		!! determined by the following values for c_method:
1238		!! 1. Second order finite difference (not recommended)
1239		!! 2. Second order finite volume (used in original PPM)
1240		!! 3. Finite-volume weighted least squares linear fit
1241		!! \todo The use of c_method to choose a scheme is inefficient
1242		!! and should eventually be moved up the call tree.
1243
1244		! Local variables
1245		integer :: k
1246		real :: hkm1, hk, hkp1 ! Successive layer thicknesses [H ~> m or kg m-2] or other units
1247		real :: Skm1, Sk, Skp1 ! Successive layer tracer concentrations in the same arbitrary units as S [A ~> a]
1248		real :: diff_l, diff_r, diff_c ! Differences in tracer concentrations in arbitrary units [A ~> a]
1249
1250	0	do k = 2, nk-1
1251	0	hkm1 = h(k-1)
1252	0	hk = h(k)
1253	0	hkp1 = h(k+1)
1254
1255	0	if ( ( hkp1 + hk ) * ( hkm1 + hk ) > 0.) then
1256	0	Skm1 = S(k-1)
1257	0	Sk = S(k)
1258	0	Skp1 = S(k+1)
1259	0	if (c_method==1) then
1260		! Simple centered diff (from White)
1261	0	if ( hk + 0.5 * (hkm1 + hkp1) /= 0. ) then
1262	0	diff_c = ( Skp1 - Skm1 ) * ( hk / ( hk + 0.5 * (hkm1 + hkp1) ) )
1263		else
1264	0	diff_c = 0.
1265		endif
1266	0	elseif (c_method==2) then
1267		! Second order accurate centered finite-volume slope (from Colella and Woodward, JCP 1984)
1268	0	diff_c = fv_diff(hkm1, hk, hkp1, Skm1, Sk, Skp1)
1269	0	elseif (c_method==3) then
1270		! Second order accurate finite-volume least squares slope
1271	0	diff_c = hk * fvlsq_slope(hkm1, hk, hkp1, Skm1, Sk, Skp1)
1272		endif
1273		! Limit centered slope by twice the side differenced slopes
1274	0	diff_l = 2. * ( Sk - Skm1 )
1275	0	diff_r = 2. * ( Skp1 - Sk )
1276	0	if ( signum(1., diff_l) * signum(1., diff_r) <= 0. ) then
1277	0	diff(k) = 0. ! PCM for local extrema
1278		else
1279	0	diff(k) = sign( min( abs(diff_l), abs(diff_c), abs(diff_r) ), diff_c )
1280		endif
1281		else
1282	0	diff(k) = 0. ! PCM next to vanished layers
1283		endif
1284		enddo
1285	0	if (b_method==1) then ! PCM for top and bottom layer
1286	0	diff(1) = 0.
1287	0	diff(nk) = 0.
1288	0	elseif (b_method==2) then ! Linear extrapolation for top and bottom interfaces
1289	0	diff(1) = ( S(2) - S(1) ) * 2. * ( h(1) / ( h(1) + h(2) ) )
1290	0	diff(nk) = S(nk) - S(nk-1) * 2. * ( h(nk) / ( h(nk-1) + h(nk) ) )
1291		endif
1292
1293	0	end subroutine PLM_diff
1294
1295		!> Returns the cell-centered second-order finite volume (unlimited PLM) slope using three
1296		!! consecutive cell widths and average values. Slope is returned as a difference across
1297		!! the central cell (i.e. units of scalar S, e.g. [C ~> degC] for temperature).
1298		!! Discretization follows equation 1.7 in Colella & Woodward, 1984: JCP 54, 174-201.
1299	0	real function fv_diff(hkm1, hk, hkp1, Skm1, Sk, Skp1)
1300		real, intent(in) :: hkm1 !< Left cell width [H ~> m or kg m-2] or other arbitrary units
1301		real, intent(in) :: hk !< Center cell width [H ~> m or kg m-2] or other arbitrary units
1302		real, intent(in) :: hkp1 !< Right cell width [H ~> m or kg m-2] or other arbitrary units
1303		real, intent(in) :: Skm1 !< Left cell average value in arbitrary concentration
1304		!! units (e.g. [C ~> degC] for temperature)
1305		real, intent(in) :: Sk !< Center cell average value in arbitrary concentration
1306		!! units (e.g. [C ~> degC] for temperature)
1307		real, intent(in) :: Skp1 !< Right cell average value in arbitrary concentration
1308		!! units (e.g. [C ~> degC] for temperature)
1309
1310		! Local variables
1311		real :: h_sum, hp, hm ! At first sums of thicknesses [H ~> m or kg m-2], then changed into
1312		! their reciprocals [H-1 ~> m-1 or m2 kg-1]
1313
1314	0	h_sum = ( hkm1 + hkp1 ) + hk
1315	0	if (h_sum /= 0.) h_sum = 1./ h_sum
1316	0	hm = hkm1 + hk
1317	0	if (hm /= 0.) hm = 1./ hm
1318	0	hp = hkp1 + hk
1319	0	if (hp /= 0.) hp = 1./ hp
1320		fv_diff = ( hk * h_sum ) * &
1321		( ( 2. * hkm1 + hk ) * hp * ( Skp1 - Sk ) &
1322	0	+ ( 2. * hkp1 + hk ) * hm * ( Sk - Skm1 ) )
1323	0	end function fv_diff
1324
1325
1326		!> Returns the cell-centered second-order weighted least squares slope
1327		!! using three consecutive cell widths and average values. Slope is returned
1328		!! as a gradient (i.e. units of scalar S over width units). For example, for temperature
1329		!! fvlsq_slope would usually be returned in units of [C H-1 ~> degC m-1 or degC m2 kg-1].
1330	0	real function fvlsq_slope(hkm1, hk, hkp1, Skm1, Sk, Skp1)
1331		real, intent(in) :: hkm1 !< Left cell width [H ~> m or kg m-2] or other arbitrary units
1332		real, intent(in) :: hk !< Center cell width [H ~> m or kg m-2] or other arbitrary units
1333		real, intent(in) :: hkp1 !< Right cell width [H ~> m or kg m-2] or other arbitrary units
1334		real, intent(in) :: Skm1 !< Left cell average value in arbitrary concentration
1335		!! units (e.g. [C ~> degC] for temperature)
1336		real, intent(in) :: Sk !< Center cell average value often in arbitrary concentration
1337		!! units (e.g. [C ~> degC] for temperature)
1338		real, intent(in) :: Skp1 !< Right cell average value often in arbitrary concentration
1339		!! units (e.g. [C ~> degC] for temperature)
1340
1341		! Local variables
1342		real :: xkm1, xkp1 ! Distances between layer centers [H ~> m or kg m-2] or other arbitrary units
1343		real :: h_sum ! Sum of the successive cell widths [H ~> m or kg m-2] or other arbitrary units
1344		real :: hx_sum ! Thicknesses times distances [H2 ~> m2 or kg2 m-4]
1345		real :: hxsq_sum ! Thicknesses times squared distances [H3 ~> m3 or kg3 m-6]
1346		real :: det ! The denominator in the weighted slope calculation [H4 ~> m4 or kg4 m-8]
1347		real :: hxy_sum ! Sum of layer concentrations times thicknesses and distances in units that
1348		! depend on those of Sk (e.g. [C H2 ~> degC m2 or degC kg2 m-4] for temperature)
1349		real :: hy_sum ! Sum of layer concentrations times thicknesses in units that depend on
1350		! those of Sk (e.g. [C H ~> degC m or degC kg m-2] for temperature)
1351
1352	0	xkm1 = -0.5 * ( hk + hkm1 )
1353	0	xkp1 = 0.5 * ( hk + hkp1 )
1354	0	h_sum = ( hkm1 + hkp1 ) + hk
1355	0	hx_sum = hkm1xkm1 + hkp1xkp1
1356	0	hxsq_sum = hkm1(xkm12) + hkp1(xkp1**2)
1357	0	hxy_sum = hkm1xkm1Skm1 + hkp1xkp1Skp1
1358	0	hy_sum = ( hkm1Skm1 + hkp1Skp1 ) + hk*Sk
1359	0	det = h_sum * hxsq_sum - hx_sum**2
1360	0	if (det /= 0.) then
1361		!a = ( hxsq_sum * hy_sum - hx_sum*hxy_sum ) / det ! a would be mean of straight line fit
1362	0	fvlsq_slope = ( h_sum * hxy_sum - hx_sum*hy_sum ) / det ! Gradient of straight line fit
1363		else
1364	0	fvlsq_slope = 0. ! Adcroft's reciprocal rule
1365		endif
1366	0	end function fvlsq_slope
1367
1368
1369		!> Returns positions within left/right columns of combined interfaces using continuous reconstructions of T/S
1370	0	subroutine find_neutral_surface_positions_continuous(nk, Pl, Tl, Sl, dRdTl, dRdSl, Pr, Tr, Sr, &
1371	0	dRdTr, dRdSr, PoL, PoR, KoL, KoR, hEff, bl_kl, bl_kr, bl_zl, bl_zr)
1372		integer, intent(in) :: nk !< Number of levels
1373		real, dimension(nk+1), intent(in) :: Pl !< Left-column interface pressure [R L2 T-2 ~> Pa] or other units
1374		real, dimension(nk+1), intent(in) :: Tl !< Left-column interface potential temperature [C ~> degC]
1375		real, dimension(nk+1), intent(in) :: Sl !< Left-column interface salinity [S ~> ppt]
1376		real, dimension(nk+1), intent(in) :: dRdTl !< Left-column dRho/dT [R C-1 ~> kg m-3 degC-1]
1377		real, dimension(nk+1), intent(in) :: dRdSl !< Left-column dRho/dS [R S-1 ~> kg m-3 ppt-1]
1378		real, dimension(nk+1), intent(in) :: Pr !< Right-column interface pressure [R L2 T-2 ~> Pa] or other units
1379		real, dimension(nk+1), intent(in) :: Tr !< Right-column interface potential temperature [C ~> degC]
1380		real, dimension(nk+1), intent(in) :: Sr !< Right-column interface salinity [S ~> ppt]
1381		real, dimension(nk+1), intent(in) :: dRdTr !< Left-column dRho/dT [R C-1 ~> kg m-3 degC-1]
1382		real, dimension(nk+1), intent(in) :: dRdSr !< Left-column dRho/dS [R S-1 ~> kg m-3 ppt-1]
1383		real, dimension(2*nk+2), intent(inout) :: PoL !< Fractional position of neutral surface within
1384		!! layer KoL of left column [nondim]
1385		real, dimension(2*nk+2), intent(inout) :: PoR !< Fractional position of neutral surface within
1386		!! layer KoR of right column [nondim]
1387		integer, dimension(2*nk+2), intent(inout) :: KoL !< Index of first left interface above neutral surface
1388		integer, dimension(2*nk+2), intent(inout) :: KoR !< Index of first right interface above neutral surface
1389		real, dimension(2*nk+1), intent(inout) :: hEff !< Effective thickness between two neutral surfaces
1390		!! [R L2 T-2 ~> Pa] or other units following Pl and Pr.
1391		integer, optional, intent(in) :: bl_kl !< Layer index of the boundary layer (left)
1392		integer, optional, intent(in) :: bl_kr !< Layer index of the boundary layer (right)
1393		real, optional, intent(in) :: bl_zl !< Fractional position of the boundary layer (left) [nondim]
1394		real, optional, intent(in) :: bl_zr !< Fractional position of the boundary layer (right) [nondim]
1395
1396		! Local variables
1397		integer :: ns ! Number of neutral surfaces
1398		integer :: k_surface ! Index of neutral surface
1399		integer :: kl ! Index of left interface
1400		integer :: kr ! Index of right interface
1401		logical :: searching_left_column ! True if searching for the position of a right interface in the left column
1402		logical :: searching_right_column ! True if searching for the position of a left interface in the right column
1403		logical :: reached_bottom ! True if one of the bottom-most interfaces has been used as the target
1404		integer :: krm1, klm1
1405		real :: dRho, dRhoTop, dRhoBot ! Potential density differences at various points [R ~> kg m-3]
1406		real :: hL, hR ! Pressure thicknesses [R L2 T-2 ~> Pa]
1407		integer :: lastK_left, lastK_right ! Layers used during the last iteration
1408		real :: lastP_left, lastP_right ! Fractional positions during the last iteration [nondim]
1409		logical :: interior_limit
1410
1411	0	ns = 2*nk+2
1412
1413		! Initialize variables for the search
1414	0	kr = 1
1415	0	kl = 1
1416	0	lastP_right = 0.
1417	0	lastP_left = 0.
1418	0	lastK_right = 1
1419	0	lastK_left = 1
1420	0	reached_bottom = .false.
1421
1422		! Check to see if we should limit the diffusion to the interior
1423	0	interior_limit = PRESENT(bl_kl) .and. PRESENT(bl_kr) .and. PRESENT(bl_zr) .and. PRESENT(bl_zl)
1424
1425		! Loop over each neutral surface, working from top to bottom
1426	0	neutral_surfaces: do k_surface = 1, ns
1427	0	klm1 = max(kl-1, 1)
1428	0	if (klm1>nk) stop 'find_neutral_surface_positions(): klm1 went out of bounds!'
1429	0	krm1 = max(kr-1, 1)
1430	0	if (krm1>nk) stop 'find_neutral_surface_positions(): krm1 went out of bounds!'
1431
1432		! Potential density difference, rho(kr) - rho(kl)
1433		dRho = 0.5 * ( ( dRdTr(kr) + dRdTl(kl) ) * ( Tr(kr) - Tl(kl) ) &
1434	0	+ ( dRdSr(kr) + dRdSl(kl) ) * ( Sr(kr) - Sl(kl) ) )
1435		! Which column has the lighter surface for the current indexes, kr and kl
1436	0	if (.not. reached_bottom) then
1437	0	if (dRho < 0.) then
1438	0	searching_left_column = .true.
1439	0	searching_right_column = .false.
1440	0	elseif (dRho > 0.) then
1441	0	searching_right_column = .true.
1442	0	searching_left_column = .false.
1443		else ! dRho == 0.
1444	0	if (kl + kr == 2) then ! Still at surface
1445	0	searching_left_column = .true.
1446	0	searching_right_column = .false.
1447		else ! Not the surface so we simply change direction
1448	0	searching_left_column = .not. searching_left_column
1449	0	searching_right_column = .not. searching_right_column
1450		endif
1451		endif
1452		endif
1453
1454	0	if (searching_left_column) then
1455		! Interpolate for the neutral surface position within the left column, layer klm1
1456		! Potential density difference, rho(kl-1) - rho(kr) (should be negative)
1457		dRhoTop = 0.5 * ( ( dRdTl(klm1) + dRdTr(kr) ) * ( Tl(klm1) - Tr(kr) ) &
1458	0	+ ( dRdSl(klm1) + dRdSr(kr) ) * ( Sl(klm1) - Sr(kr) ) )
1459		! Potential density difference, rho(kl) - rho(kr) (will be positive)
1460		dRhoBot = 0.5 * ( ( dRdTl(klm1+1) + dRdTr(kr) ) * ( Tl(klm1+1) - Tr(kr) ) &
1461	0	+ ( dRdSl(klm1+1) + dRdSr(kr) ) * ( Sl(klm1+1) - Sr(kr) ) )
1462
1463		! Because we are looking left, the right surface, kr, is lighter than klm1+1 and should be denser than klm1
1464		! unless we are still at the top of the left column (kl=1)
1465	0	if (dRhoTop > 0. .or. kr+kl==2) then
1466	0	PoL(k_surface) = 0. ! The right surface is lighter than anything in layer klm1
1467	0	elseif (dRhoTop >= dRhoBot) then ! Left layer is unstratified
1468	0	PoL(k_surface) = 1.
1469		else
1470		! Linearly interpolate for the position between Pl(kl-1) and Pl(kl) where the density difference
1471		! between right and left is zero. The Pl here are only used to handle massless layers.
1472	0	PoL(k_surface) = interpolate_for_nondim_position( dRhoTop, Pl(klm1), dRhoBot, Pl(klm1+1) )
1473		endif
1474	0	if (PoL(k_surface)>=1. .and. klm1<nk) then ! >= is really ==, when PoL==1 we point to the bottom of the cell
1475	0	klm1 = klm1 + 1
1476	0	PoL(k_surface) = PoL(k_surface) - 1.
1477		endif
1478	0	if (real(klm1-lastK_left)+(PoL(k_surface)-lastP_left)<0.) then
1479	0	PoL(k_surface) = lastP_left
1480	0	klm1 = lastK_left
1481		endif
1482	0	KoL(k_surface) = klm1
1483	0	if (kr <= nk) then
1484	0	PoR(k_surface) = 0.
1485	0	KoR(k_surface) = kr
1486		else
1487	0	PoR(k_surface) = 1.
1488	0	KoR(k_surface) = nk
1489		endif
1490	0	if (kr <= nk) then
1491	0	kr = kr + 1
1492		else
1493	0	reached_bottom = .true.
1494	0	searching_right_column = .true.
1495	0	searching_left_column = .false.
1496		endif
1497	0	elseif (searching_right_column) then
1498		! Interpolate for the neutral surface position within the right column, layer krm1
1499		! Potential density difference, rho(kr-1) - rho(kl) (should be negative)
1500		dRhoTop = 0.5 * ( ( dRdTr(krm1) + dRdTl(kl) ) * ( Tr(krm1) - Tl(kl) ) + &
1501	0	( dRdSr(krm1) + dRdSl(kl) ) * ( Sr(krm1) - Sl(kl) ) )
1502		! Potential density difference, rho(kr) - rho(kl) (will be positive)
1503		dRhoBot = 0.5 * ( ( dRdTr(krm1+1) + dRdTl(kl) ) * ( Tr(krm1+1) - Tl(kl) ) + &
1504	0	( dRdSr(krm1+1) + dRdSl(kl) ) * ( Sr(krm1+1) - Sl(kl) ) )
1505
1506		! Because we are looking right, the left surface, kl, is lighter than krm1+1 and should be denser than krm1
1507		! unless we are still at the top of the right column (kr=1)
1508	0	if (dRhoTop >= 0. .or. kr+kl==2) then
1509	0	PoR(k_surface) = 0. ! The left surface is lighter than anything in layer krm1
1510	0	elseif (dRhoTop >= dRhoBot) then ! Right layer is unstratified
1511	0	PoR(k_surface) = 1.
1512		else
1513		! Linearly interpolate for the position between Pr(kr-1) and Pr(kr) where the density difference
1514		! between right and left is zero. The Pr here are only used to handle massless layers.
1515	0	PoR(k_surface) = interpolate_for_nondim_position( dRhoTop, Pr(krm1), dRhoBot, Pr(krm1+1) )
1516		endif
1517	0	if (PoR(k_surface)>=1. .and. krm1<nk) then ! >= is really ==, when PoR==1 we point to the bottom of the cell
1518	0	krm1 = krm1 + 1
1519	0	PoR(k_surface) = PoR(k_surface) - 1.
1520		endif
1521	0	if (real(krm1-lastK_right)+(PoR(k_surface)-lastP_right)<0.) then
1522	0	PoR(k_surface) = lastP_right
1523	0	krm1 = lastK_right
1524		endif
1525	0	KoR(k_surface) = krm1
1526	0	if (kl <= nk) then
1527	0	PoL(k_surface) = 0.
1528	0	KoL(k_surface) = kl
1529		else
1530	0	PoL(k_surface) = 1.
1531	0	KoL(k_surface) = nk
1532		endif
1533	0	if (kl <= nk) then
1534	0	kl = kl + 1
1535		else
1536	0	reached_bottom = .true.
1537	0	searching_right_column = .false.
1538	0	searching_left_column = .true.
1539		endif
1540		else
1541	0	stop 'Else what?'
1542		endif
1543	0	if (interior_limit) then
1544	0	if (KoL(k_surface)<=bl_kl) then
1545	0	KoL(k_surface) = bl_kl
1546	0	if (PoL(k_surface)<bl_zl) then
1547	0	PoL(k_surface) = bl_zl
1548		endif
1549		endif
1550	0	if (KoR(k_surface)<=bl_kr) then
1551	0	KoR(k_surface) = bl_kr
1552	0	if (PoR(k_surface)<bl_zr) then
1553	0	PoR(k_surface) = bl_zr
1554		endif
1555		endif
1556		endif
1557
1558	0	lastK_left = KoL(k_surface) ; lastP_left = PoL(k_surface)
1559	0	lastK_right = KoR(k_surface) ; lastP_right = PoR(k_surface)
1560		! Effective thickness
1561		! NOTE: This would be better expressed in terms of the layers thicknesses rather
1562		! than as differences of position - AJA
1563	0	if (k_surface>1) then
1564	0	hL = absolute_position(nk,ns,Pl,KoL,PoL,k_surface) - absolute_position(nk,ns,Pl,KoL,PoL,k_surface-1)
1565	0	hR = absolute_position(nk,ns,Pr,KoR,PoR,k_surface) - absolute_position(nk,ns,Pr,KoR,PoR,k_surface-1)
1566	0	if ( hL + hR > 0.) then
1567	0	hEff(k_surface-1) = 2. * hL * hR / ( hL + hR ) ! Harmonic mean of layer thicknesses
1568		else
1569	0	hEff(k_surface-1) = 0.
1570		endif
1571		endif
1572
1573		enddo neutral_surfaces
1574
1575	0	end subroutine find_neutral_surface_positions_continuous
1576
1577		!> Returns the non-dimensional position between Pneg and Ppos where the
1578		!! interpolated density difference equals zero [nondim].
1579		!! The result is always bounded to be between 0 and 1.
1580	0	real function interpolate_for_nondim_position(dRhoNeg, Pneg, dRhoPos, Ppos)
1581		real, intent(in) :: dRhoNeg !< Negative density difference [R ~> kg m-3]
1582		real, intent(in) :: Pneg !< Position of negative density difference [R L2 T-2 ~> Pa] or [nondim]
1583		real, intent(in) :: dRhoPos !< Positive density difference [R ~> kg m-3]
1584		real, intent(in) :: Ppos !< Position of positive density difference [R L2 T-2 ~> Pa] or [nondim]
1585
1586		character(len=120) :: mesg
1587
1588	0	if ((Ppos > Pneg) .and. (dRhoPos - dRhoNeg >= 0. )) then
1589	0	if ( dRhoPos - dRhoNeg > 0. ) then
1590	0	interpolate_for_nondim_position = min( 1., max( 0., -dRhoNeg / ( dRhoPos - dRhoNeg ) ) )
1591	0	elseif (dRhoPos - dRhoNeg == 0) then
1592	0	if (dRhoNeg > 0.) then
1593	0	interpolate_for_nondim_position = 0.
1594	0	elseif (dRhoNeg < 0.) then
1595	0	interpolate_for_nondim_position = 1.
1596		else ! dRhoPos = dRhoNeg = 0
1597	0	interpolate_for_nondim_position = 0.5
1598		endif
1599		else ! dRhoPos - dRhoNeg < 0
1600	0	interpolate_for_nondim_position = 0.5
1601		endif
1602	0	elseif (Ppos == Pneg) then ! Handle vanished or inverted layers
1603	0	interpolate_for_nondim_position = 0.5
1604		else ! ((Ppos < Pneg) .or. (dRhoNeg > dRhoPos) )
1605		! Error handling for problematic cases. It is expected that this should never occur.
1606	0	write(mesg,*) 'dRhoNeg, Pneg, dRhoPos, Ppos', dRhoNeg, Pneg, dRhoPos, Ppos
1607	0	call MOM_error(WARNING, 'interpolate_for_nondim_position: '//trim(mesg))
1608		! write(stderr,*) trim(mesg)
1609	0	if ((Ppos < Pneg) .and. (dRhoNeg > dRhoPos)) then
1610	0	mesg = '(Ppos < Pneg) and (dRhoNeg > dRhoPos)'
1611	0	elseif (Ppos < Pneg) then
1612	0	mesg = 'Ppos < Pneg'
1613	0	elseif (dRhoNeg > dRhoPos) then
1614	0	mesg = trim(mesg)//'; dRhoNeg > dRhoPos'
1615		else ! This should never happen.
1616	0	mesg = 'Unexpected failure.'
1617		endif
1618	0	call MOM_error(FATAL, 'interpolate_for_nondim_position: '//trim(mesg))
1619		endif
1620
1621	0	end function interpolate_for_nondim_position
1622
1623		!> Higher order version of find_neutral_surface_positions. Returns positions within left/right columns
1624		!! of combined interfaces using intracell reconstructions of T/S. Note that the polynomial reconstructions
1625		!! of T and S are optional to aid with unit testing, but will always be passed otherwise
1626	0	subroutine find_neutral_surface_positions_discontinuous(CS, nk, &
1627	0	Pres_l, hcol_l, Tl, Sl, ppoly_T_l, ppoly_S_l, stable_l, &
1628	0	Pres_r, hcol_r, Tr, Sr, ppoly_T_r, ppoly_S_r, stable_r, &
1629	0	PoL, PoR, KoL, KoR, hEff, zeta_bot_L, zeta_bot_R, k_bot_L, k_bot_R, hard_fail_heff)
1630
1631		type(neutral_diffusion_CS), intent(inout) :: CS !< Neutral diffusion control structure
1632		integer, intent(in) :: nk !< Number of levels
1633		real, dimension(nk,2), intent(in) :: Pres_l !< Left-column interface pressure [R L2 T-2 ~> Pa]
1634		real, dimension(nk), intent(in) :: hcol_l !< Left-column layer thicknesses [H ~> m or kg m-2]
1635		!! or other units
1636		real, dimension(nk,2), intent(in) :: Tl !< Left-column top interface potential
1637		!! temperature [C ~> degC]
1638		real, dimension(nk,2), intent(in) :: Sl !< Left-column top interface salinity [S ~> ppt]
1639		real, dimension(:,:), intent(in) :: ppoly_T_l !< Left-column coefficients of T reconstruction [C ~> degC]
1640		real, dimension(:,:), intent(in) :: ppoly_S_l !< Left-column coefficients of S reconstruction [S ~> ppt]
1641		logical, dimension(nk), intent(in) :: stable_l !< True where the left-column is stable
1642		real, dimension(nk,2), intent(in) :: Pres_r !< Right-column interface pressure [R L2 T-2 ~> Pa]
1643		real, dimension(nk), intent(in) :: hcol_r !< Left-column layer thicknesses [H ~> m or kg m-2]
1644		!! or other units
1645		real, dimension(nk,2), intent(in) :: Tr !< Right-column top interface potential
1646		!! temperature [C ~> degC]
1647		real, dimension(nk,2), intent(in) :: Sr !< Right-column top interface salinity [S ~> ppt]
1648		real, dimension(:,:), intent(in) :: ppoly_T_r !< Right-column coefficients of T
1649		!! reconstruction [C ~> degC]
1650		real, dimension(:,:), intent(in) :: ppoly_S_r !< Right-column coefficients of S reconstruction [S ~> ppt]
1651		logical, dimension(nk), intent(in) :: stable_r !< True where the right-column is stable
1652		real, dimension(4*nk), intent(inout) :: PoL !< Fractional position of neutral surface within
1653		!! layer KoL of left column [nondim]
1654		real, dimension(4*nk), intent(inout) :: PoR !< Fractional position of neutral surface within
1655		!! layer KoR of right column [nondim]
1656		integer, dimension(4*nk), intent(inout) :: KoL !< Index of first left interface above neutral surface
1657		integer, dimension(4*nk), intent(inout) :: KoR !< Index of first right interface above neutral surface
1658		real, dimension(4*nk-1), intent(inout) :: hEff !< Effective thickness between two neutral surfaces
1659		!! [H ~> m or kg m-2] or other units taken from hcol_l
1660		real, optional, intent(in) :: zeta_bot_L!< Non-dimensional distance to where the boundary layer
1661		!! intersects the cell (left) [nondim]
1662		real, optional, intent(in) :: zeta_bot_R!< Non-dimensional distance to where the boundary layer
1663		!! intersects the cell (right) [nondim]
1664
1665		integer, optional, intent(in) :: k_bot_L !< k-index for the boundary layer (left) [nondim]
1666		integer, optional, intent(in) :: k_bot_R !< k-index for the boundary layer (right) [nondim]
1667		logical, optional, intent(in) :: hard_fail_heff !< If true (default) bring down the model if the
1668		!! neutral surfaces ever cross
1669		! Local variables
1670		integer :: ns ! Number of neutral surfaces
1671		integer :: k_surface ! Index of neutral surface
1672		integer :: kl_left, kl_right ! Index of layers on the left/right
1673		integer :: ki_left, ki_right ! Index of interfaces on the left/right
1674		logical :: searching_left_column ! True if searching for the position of a right interface in the left column
1675		logical :: searching_right_column ! True if searching for the position of a left interface in the right column
1676		logical :: reached_bottom ! True if one of the bottom-most interfaces has been used as the target
1677		logical :: fail_heff ! Fail if negative thickness are encountered. By default this
1678		! is true, but it can take its value from hard_fail_heff.
1679		real :: dRho ! A density difference between columns [R ~> kg m-3]
1680		real :: hL, hR ! Left and right layer thicknesses [H ~> m or kg m-2] or units from hcol_l
1681		real :: lastP_left, lastP_right ! Previous positions for left and right [nondim]
1682		integer :: k_init_L, k_init_R ! Starting indices layers for left and right
1683		real :: p_init_L, p_init_R ! Starting positions for left and right [nondim]
1684		! Initialize variables for the search
1685	0	ns = 4*nk
1686	0	ki_right = 1
1687	0	ki_left = 1
1688	0	kl_left = 1
1689	0	kl_right = 1
1690	0	lastP_left = 0.
1691	0	lastP_right = 0.
1692	0	reached_bottom = .false.
1693	0	searching_left_column = .false.
1694	0	searching_right_column = .false.
1695
1696	0	fail_heff = .true.
1697	0	if (PRESENT(hard_fail_heff)) fail_heff = hard_fail_heff
1698
1699	0	if (PRESENT(k_bot_L) .and. PRESENT(k_bot_R) .and. PRESENT(zeta_bot_L) .and. PRESENT(zeta_bot_R)) then
1700	0	k_init_L = k_bot_L ; k_init_R = k_bot_R
1701	0	p_init_L = zeta_bot_L ; p_init_R = zeta_bot_R
1702	0	lastP_left = zeta_bot_L ; lastP_right = zeta_bot_R
1703	0	kl_left = k_bot_L ; kl_right = k_bot_R
1704		else
1705	0	k_init_L = 1 ; k_init_R = 1
1706	0	p_init_L = 0. ; p_init_R = 0.
1707		endif
1708		! Loop over each neutral surface, working from top to bottom
1709	0	neutral_surfaces: do k_surface = 1, ns
1710
1711	0	if (k_surface == ns) then
1712	0	PoL(k_surface) = 1.
1713	0	PoR(k_surface) = 1.
1714	0	KoL(k_surface) = nk
1715	0	KoR(k_surface) = nk
1716		! If the layers are unstable, then simply point the surface to the previous location
1717	0	elseif (.not. stable_l(kl_left)) then
1718	0	if (k_surface > 1) then
1719	0	PoL(k_surface) = ki_left - 1 ! Top interface is at position = 0., Bottom is at position = 1
1720	0	KoL(k_surface) = kl_left
1721	0	PoR(k_surface) = PoR(k_surface-1)
1722	0	KoR(k_surface) = KoR(k_surface-1)
1723		else
1724	0	PoR(k_surface) = p_init_R
1725	0	KoR(k_surface) = k_init_R
1726	0	PoL(k_surface) = p_init_L
1727	0	KoL(k_Surface) = k_init_L
1728		endif
1729	0	call increment_interface(nk, kl_left, ki_left, reached_bottom, searching_left_column, searching_right_column)
1730	0	searching_left_column = .true.
1731	0	searching_right_column = .false.
1732	0	elseif (.not. stable_r(kl_right)) then ! Check the right layer for stability
1733	0	if (k_surface > 1) then
1734	0	PoR(k_surface) = ki_right - 1 ! Top interface is at position = 0., Bottom is at position = 1
1735	0	KoR(k_surface) = kl_right
1736	0	PoL(k_surface) = PoL(k_surface-1)
1737	0	KoL(k_surface) = KoL(k_surface-1)
1738		else
1739	0	PoR(k_surface) = 0.
1740	0	KoR(k_surface) = 1
1741	0	PoL(k_surface) = 0.
1742	0	KoL(k_surface) = 1
1743		endif
1744	0	call increment_interface(nk, kl_right, ki_right, reached_bottom, searching_right_column, searching_left_column)
1745	0	searching_left_column = .false.
1746	0	searching_right_column = .true.
1747		else ! Layers are stable so need to figure out whether we need to search right or left
1748		! For convenience, the left column uses the searched "from" interface variables, and the right column
1749		! uses the searched 'to'. These will get reset in subsequent calc_delta_rho calls
1750
1751		call calc_delta_rho_and_derivs(CS, &
1752		Tr(kl_right, ki_right), Sr(kl_right, ki_right), Pres_r(kl_right,ki_right), &
1753		Tl(kl_left, ki_left), Sl(kl_left, ki_left) , Pres_l(kl_left,ki_left), &
1754	0	dRho)
1755	0	if (CS%debug) write(stdout,'(A,I0,A,E12.4,A,I0,A,I0,A,I0,A,I0)') &
1756	0	"k_surface=",k_surface, " dRho=",CS%R_to_kg_m3*dRho, &
1757	0	"kl_left=",kl_left, " ki_left=",ki_left, " kl_right=",kl_right, " ki_right=",ki_right
1758		! Which column has the lighter surface for the current indexes, kr and kl
1759	0	if (.not. reached_bottom) then
1760	0	if (dRho < 0.) then
1761	0	searching_left_column = .true.
1762	0	searching_right_column = .false.
1763	0	elseif (dRho > 0.) then
1764	0	searching_left_column = .false.
1765	0	searching_right_column = .true.
1766		else ! dRho == 0.
1767	0	if ( ( kl_left + kl_right == 2 ) .and. (ki_left + ki_right == 2) ) then ! Still at surface
1768	0	searching_left_column = .true.
1769	0	searching_right_column = .false.
1770		else ! Not the surface so we simply change direction
1771	0	searching_left_column = .not. searching_left_column
1772	0	searching_right_column = .not. searching_right_column
1773		endif
1774		endif
1775		endif
1776	0	if (searching_left_column) then
1777		! Position of the right interface is known and all quantities are fixed
1778	0	PoR(k_surface) = ki_right - 1.
1779	0	KoR(k_surface) = kl_right
1780		PoL(k_surface) = search_other_column(CS, k_surface, lastP_left, &
1781		Tr(kl_right, ki_right), Sr(kl_right, ki_right), Pres_r(kl_right, ki_right), &
1782		Tl(kl_left,1), Sl(kl_left,1), Pres_l(kl_left,1), &
1783		Tl(kl_left,2), Sl(kl_left,2), Pres_l(kl_left,2), &
1784	0	ppoly_T_l(kl_left,:), ppoly_S_l(kl_left,:))
1785	0	KoL(k_surface) = kl_left
1786
1787	0	if (CS%debug) then
1788	0	write(stdout,'(A,I0)') "Searching left layer ", kl_left
1789	0	write(stdout,'(A,I0,1X,I0)') "Searching from right: ", kl_right, ki_right
1790	0	write(stdout,*) "Temp/Salt Reference: ", Tr(kl_right,ki_right), Sr(kl_right,ki_right)
1791	0	write(stdout,*) "Temp/Salt Top L: ", Tl(kl_left,1), Sl(kl_left,1)
1792	0	write(stdout,*) "Temp/Salt Bot L: ", Tl(kl_left,2), Sl(kl_left,2)
1793		endif
1794	0	call increment_interface(nk, kl_right, ki_right, reached_bottom, searching_right_column, searching_left_column)
1795	0	lastP_left = PoL(k_surface)
1796		! If the right layer increments, then we need to reset the last position on the right
1797	0	if ( kl_right == (KoR(k_surface) + 1) ) lastP_right = 0.
1798
1799	0	elseif (searching_right_column) then
1800		! Position of the right interface is known and all quantities are fixed
1801	0	PoL(k_surface) = ki_left - 1.
1802	0	KoL(k_surface) = kl_left
1803		PoR(k_surface) = search_other_column(CS, k_surface, lastP_right, &
1804		Tl(kl_left, ki_left), Sl(kl_left, ki_left), Pres_l(kl_left, ki_left), &
1805		Tr(kl_right,1), Sr(kl_right,1), Pres_r(kl_right,1), &
1806		Tr(kl_right,2), Sr(kl_right,2), Pres_r(kl_right,2), &
1807	0	ppoly_T_r(kl_right,:), ppoly_S_r(kl_right,:))
1808	0	KoR(k_surface) = kl_right
1809
1810	0	if (CS%debug) then
1811	0	write(stdout,'(A,I0)') "Searching right layer ", kl_right
1812	0	write(stdout,'(A,I0,1X,I0)') "Searching from left: ", kl_left, ki_left
1813	0	write(stdout,*) "Temp/Salt Reference: ", Tl(kl_left,ki_left), Sl(kl_left,ki_left)
1814	0	write(stdout,*) "Temp/Salt Top L: ", Tr(kl_right,1), Sr(kl_right,1)
1815	0	write(stdout,*) "Temp/Salt Bot L: ", Tr(kl_right,2), Sr(kl_right,2)
1816		endif
1817	0	call increment_interface(nk, kl_left, ki_left, reached_bottom, searching_left_column, searching_right_column)
1818	0	lastP_right = PoR(k_surface)
1819		! If the right layer increments, then we need to reset the last position on the right
1820	0	if ( kl_left == (KoL(k_surface) + 1) ) lastP_left = 0.
1821		else
1822	0	stop 'Else what?'
1823		endif
1824	0	if (CS%debug) write(stdout,'(A,I3,A,ES16.6,A,I0,A,ES16.6)') "KoL:", KoL(k_surface), " PoL:", PoL(k_surface), &
1825	0	" KoR:", KoR(k_surface), " PoR:", PoR(k_surface)
1826		endif
1827		! Effective thickness
1828	0	if (k_surface>1) then
1829	0	if ( KoL(k_surface) == KoL(k_surface-1) .and. KoR(k_surface) == KoR(k_surface-1) ) then
1830	0	hL = (PoL(k_surface) - PoL(k_surface-1))*hcol_l(KoL(k_surface))
1831	0	hR = (PoR(k_surface) - PoR(k_surface-1))*hcol_r(KoR(k_surface))
1832	0	if (hL < 0. .or. hR < 0.) then
1833	0	if (fail_heff) then
1834	0	call MOM_error(FATAL,"Negative thicknesses in neutral diffusion")
1835		else
1836	0	if (searching_left_column) then
1837	0	PoL(k_surface) = PoL(k_surface-1)
1838	0	KoL(k_surface) = KoL(k_surface-1)
1839	0	elseif (searching_right_column) then
1840	0	PoR(k_surface) = PoR(k_surface-1)
1841	0	KoR(k_surface) = KoR(k_surface-1)
1842		endif
1843		endif
1844	0	elseif ( hL + hR == 0. ) then
1845	0	hEff(k_surface-1) = 0.
1846		else
1847	0	hEff(k_surface-1) = 2. * ( (hL * hR) / ( hL + hR ) )! Harmonic mean
1848	0	if ( KoL(k_surface) /= KoL(k_surface-1) ) then
1849	0	call MOM_error(FATAL,"Neutral sublayer spans multiple layers")
1850		endif
1851	0	if ( KoR(k_surface) /= KoR(k_surface-1) ) then
1852	0	call MOM_error(FATAL,"Neutral sublayer spans multiple layers")
1853		endif
1854		endif
1855		else
1856	0	hEff(k_surface-1) = 0.
1857		endif
1858		endif
1859		enddo neutral_surfaces
1860	0	end subroutine find_neutral_surface_positions_discontinuous
1861
1862		!> Sweep down through the column and mark as stable if the bottom interface of a cell is denser than the top
1863	0	subroutine mark_unstable_cells(CS, nk, T, S, P, stable_cell)
1864		type(neutral_diffusion_CS), intent(inout) :: CS !< Neutral diffusion control structure
1865		integer, intent(in) :: nk !< Number of levels in a column
1866		real, dimension(nk,2), intent(in) :: T !< Temperature at interfaces [C ~> degC]
1867		real, dimension(nk,2), intent(in) :: S !< Salinity at interfaces [S ~> ppt]
1868		real, dimension(nk,2), intent(in) :: P !< Pressure at interfaces [R L2 T-2 ~> Pa]
1869		logical, dimension(nk), intent( out) :: stable_cell !< True if this cell is unstably stratified
1870
1871		integer :: k
1872		real :: delta_rho ! A density difference [R ~> kg m-3]
1873
1874	0	do k = 1,nk
1875		call calc_delta_rho_and_derivs( CS, T(k,2), S(k,2), max(P(k,2), CS%ref_pres), &
1876	0	T(k,1), S(k,1), max(P(k,1), CS%ref_pres), delta_rho )
1877	0	stable_cell(k) = (delta_rho > 0.)
1878		enddo
1879	0	end subroutine mark_unstable_cells
1880
1881		!> Searches the "other" (searched) column for the position of the neutral surface, returning
1882		!! the fractional postion within the layer [nondim]
1883	0	real function search_other_column(CS, ksurf, pos_last, T_from, S_from, P_from, T_top, S_top, P_top, &
1884	0	T_bot, S_bot, P_bot, T_poly, S_poly ) result(pos)
1885		type(neutral_diffusion_CS), intent(in ) :: CS !< Neutral diffusion control structure
1886		integer, intent(in ) :: ksurf !< Current index of neutral surface
1887		real, intent(in ) :: pos_last !< Last position within the current layer, used as the lower
1888		!! bound in the root finding algorithm [nondim]
1889		real, intent(in ) :: T_from !< Temperature at the searched from interface [C ~> degC]
1890		real, intent(in ) :: S_from !< Salinity at the searched from interface [S ~> ppt]
1891		real, intent(in ) :: P_from !< Pressure at the searched from interface [R L2 T-2 ~> Pa]
1892		real, intent(in ) :: T_top !< Temperature at the searched to top interface [C ~> degC]
1893		real, intent(in ) :: S_top !< Salinity at the searched to top interface [S ~> ppt]
1894		real, intent(in ) :: P_top !< Pressure at the searched to top interface [R L2 T-2 ~> Pa]
1895		!! interface [R L2 T-2 ~> Pa]
1896		real, intent(in ) :: T_bot !< Temperature at the searched to bottom interface [C ~> degC]
1897		real, intent(in ) :: S_bot !< Salinity at the searched to bottom interface [S ~> ppt]
1898		real, intent(in ) :: P_bot !< Pressure at the searched to bottom
1899		!! interface [R L2 T-2 ~> Pa]
1900		real, dimension(:), intent(in ) :: T_poly !< Temperature polynomial reconstruction
1901		!! coefficients [C ~> degC]
1902		real, dimension(:), intent(in ) :: S_poly !< Salinity polynomial reconstruction
1903		!! coefficients [S ~> ppt]
1904		! Local variables
1905		real :: dRhotop, dRhobot ! Density differences [R ~> kg m-3]
1906		real :: dRdT_top, dRdT_bot, dRdT_from ! Partial derivatives of density with temperature [R C-1 ~> kg m-3 degC-1]
1907		real :: dRdS_top, dRdS_bot, dRdS_from ! Partial derivatives of density with salinity [R S-1 ~> kg m-3 ppt-1]
1908
1909		! Calculate the difference in density at the tops or the bottom
1910	0	if (CS%neutral_pos_method == 1 .or. CS%neutral_pos_method == 3) then
1911	0	call calc_delta_rho_and_derivs(CS, T_top, S_top, P_top, T_from, S_from, P_from, dRhoTop)
1912	0	call calc_delta_rho_and_derivs(CS, T_bot, S_bot, P_bot, T_from, S_from, P_from, dRhoBot)
1913	0	elseif (CS%neutral_pos_method == 2) then
1914		call calc_delta_rho_and_derivs(CS, T_top, S_top, P_top, T_from, S_from, P_from, dRhoTop, &
1915	0	dRdT_top, dRdS_top, dRdT_from, dRdS_from)
1916		call calc_delta_rho_and_derivs(CS, T_bot, S_bot, P_bot, T_from, S_from, P_from, dRhoBot, &
1917	0	dRdT_bot, dRdS_bot, dRdT_from, dRdS_from)
1918		endif
1919
1920		! Handle all the special cases EXCEPT if it connects within the layer
1921	0	if ( (dRhoTop > 0.) .or. (ksurf == 1) ) then ! First interface or lighter than anything in layer
1922	0	pos = pos_last
1923	0	elseif ( dRhoTop > dRhoBot ) then ! Unstably stratified
1924	0	pos = 1.
1925	0	elseif ( dRhoTop < 0. .and. dRhoBot < 0.) then ! Denser than anything in layer
1926	0	pos = 1.
1927	0	elseif ( dRhoTop == 0. .and. dRhoBot == 0. ) then ! Perfectly unstratified
1928	0	pos = 1.
1929	0	elseif ( dRhoBot == 0. ) then ! Matches perfectly at the Top
1930	0	pos = 1.
1931	0	elseif ( dRhoTop == 0. ) then ! Matches perfectly at the Bottom
1932	0	pos = pos_last
1933		else ! Neutral surface within layer
1934	0	pos = -1
1935		endif
1936
1937		! Can safely return if position is >= 0 otherwise will need to find the position within the layer
1938	0	if (pos>=0) return
1939
1940	0	if (CS%neutral_pos_method==1) then
1941	0	pos = interpolate_for_nondim_position( dRhoTop, P_top, dRhoBot, P_bot )
1942		! For the 'Linear' case of finding the neutral position, the reference pressure to use is the average
1943		! of the midpoint of the layer being searched and the interface being searched from
1944	0	elseif (CS%neutral_pos_method == 2) then
1945		pos = find_neutral_pos_linear( CS, pos_last, T_from, S_from, dRdT_from, dRdS_from, &
1946	0	dRdT_top, dRdS_top, dRdT_bot, dRdS_bot, T_poly, S_poly )
1947	0	elseif (CS%neutral_pos_method == 3) then
1948	0	pos = find_neutral_pos_full( CS, pos_last, T_from, S_from, P_from, P_top, P_bot, T_poly, S_poly)
1949		endif
1950
1951	0	end function search_other_column
1952
1953		!> Increments the interface which was just connected and also set flags if the bottom is reached
1954	0	subroutine increment_interface(nk, kl, ki, reached_bottom, searching_this_column, searching_other_column)
1955		integer, intent(in ) :: nk !< Number of vertical levels
1956		integer, intent(inout) :: kl !< Current layer (potentially updated)
1957		integer, intent(inout) :: ki !< Current interface
1958		logical, intent(inout) :: reached_bottom !< Updated when kl == nk and ki == 2
1959		logical, intent(inout) :: searching_this_column !< Updated when kl == nk and ki == 2
1960		logical, intent(inout) :: searching_other_column !< Updated when kl == nk and ki == 2
1961
1962	0	reached_bottom = .false.
1963	0	if (ki == 2) then ! At the bottom interface
1964	0	if ((ki == 2) .and. (kl < nk) ) then ! Not at the bottom so just go to the next layer
1965	0	kl = kl+1
1966	0	ki = 1
1967	0	elseif ((kl == nk) .and. (ki==2)) then
1968	0	reached_bottom = .true.
1969	0	searching_this_column = .false.
1970	0	searching_other_column = .true.
1971		endif
1972	0	elseif (ki==1) then ! At the top interface
1973	0	ki = 2 ! Next interface is same layer, but bottom interface
1974		else
1975	0	call MOM_error(FATAL,"Unanticipated eventuality in increment_interface")
1976		endif
1977	0	end subroutine increment_interface
1978
1979		!> Search a layer to find where delta_rho = 0 based on a linear interpolation of alpha and beta of the top and bottom
1980		!! being searched and polynomial reconstructions of T and S. Compressibility is not needed because either, we are
1981		!! assuming incompressibility in the equation of state for this module or alpha and beta are calculated having been
1982		!! displaced to the average pressures of the two pressures We need Newton's method because the T and S reconstructions
1983		!! make delta_rho a polynomial function of z if using PPM or higher. If Newton's method would search fall out of the
1984		!! interval [0,1], a bisection step would be taken instead. Also this linearization of alpha, beta means that second
1985		!! derivatives of the EOS are not needed. Note that delta in variable names below refers to horizontal differences and
1986		!! 'd' refers to vertical differences
1987	0	function find_neutral_pos_linear( CS, z0, T_ref, S_ref, dRdT_ref, dRdS_ref, &
1988	0	dRdT_top, dRdS_top, dRdT_bot, dRdS_bot, ppoly_T, ppoly_S ) result( z )
1989		type(neutral_diffusion_CS),intent(in) :: CS !< Control structure with parameters for this module
1990		real, intent(in) :: z0 !< Lower bound of position, also serves as the
1991		!! initial guess [nondim]
1992		real, intent(in) :: T_ref !< Temperature at the searched from interface [C ~> degC]
1993		real, intent(in) :: S_ref !< Salinity at the searched from interface [S ~> ppt]
1994		real, intent(in) :: dRdT_ref !< dRho/dT at the searched from interface
1995		!! [R C-1 ~> kg m-3 degC-1]
1996		real, intent(in) :: dRdS_ref !< dRho/dS at the searched from interface
1997		!! [R S-1 ~> kg m-3 ppt-1]
1998		real, intent(in) :: dRdT_top !< dRho/dT at top of layer being searched
1999		!! [R C-1 ~> kg m-3 degC-1]
2000		real, intent(in) :: dRdS_top !< dRho/dS at top of layer being searched
2001		!! [R S-1 ~> kg m-3 ppt-1]
2002		real, intent(in) :: dRdT_bot !< dRho/dT at bottom of layer being searched
2003		!! [R C-1 ~> kg m-3 degC-1]
2004		real, intent(in) :: dRdS_bot !< dRho/dS at bottom of layer being searched
2005		!! [R S-1 ~> kg m-3 ppt-1]
2006		real, dimension(:), intent(in) :: ppoly_T !< Coefficients of the polynomial reconstruction of T within
2007		!! the layer to be searched [C ~> degC].
2008		real, dimension(:), intent(in) :: ppoly_S !< Coefficients of the polynomial reconstruction of S within
2009		!! the layer to be searched [S ~> ppt].
2010		real :: z !< Position where drho = 0 [nondim]
2011		! Local variables
2012		real :: dRdT_diff ! Difference in the partial derivative of density with temperature across the
2013		! layer [R C-1 ~> kg m-3 degC-1]
2014		real :: dRdS_diff ! Difference in the partial derivative of density with salinity across the
2015		! layer [R S-1 ~> kg m-3 ppt-1]
2016		real :: drho, drho_dz ! Density anomaly and its derivative with fractional position [R ~> kg m-3]
2017		real :: dRdT_z ! Partial derivative of density with temperature at a point [R C-1 ~> kg m-3 degC-1]
2018		real :: dRdS_z ! Partial derivative of density with salinity at a point [R S-1 ~> kg m-3 ppt-1]
2019		real :: T_z, dT_dz ! Temperature at a point and its derivative with fractional position [C ~> degC]
2020		real :: S_z, dS_dz ! Salinity at a point and its derivative with fractional position [S ~> ppt]
2021		real :: drho_min, drho_max ! Bounds on density differences [R ~> kg m-3]
2022		real :: ztest, zmin, zmax ! Fractional positions in the cell [nondim]
2023		real :: a1, a2 ! Fractional weights of the top and bottom values [nondim]
2024		integer :: iter
2025		integer :: nterm
2026
2027	0	nterm = SIZE(ppoly_T)
2028
2029		! Position independent quantities
2030	0	dRdT_diff = dRdT_bot - dRdT_top
2031	0	dRdS_diff = dRdS_bot - dRdS_top
2032		! Initial starting drho (used for bisection)
2033	0	zmin = z0 ! Lower bounding interval
2034	0	zmax = 1. ! Maximum bounding interval (bottom of layer)
2035	0	a1 = 1. - zmin
2036	0	a2 = zmin
2037	0	T_z = evaluation_polynomial( ppoly_T, nterm, zmin )
2038	0	S_z = evaluation_polynomial( ppoly_S, nterm, zmin )
2039	0	dRdT_z = a1dRdT_top + a2dRdT_bot
2040	0	dRdS_z = a1dRdS_top + a2dRdS_bot
2041	0	drho_min = 0.5((dRdT_z+dRdT_ref)(T_z-T_ref) + (dRdS_z+dRdS_ref)*(S_z-S_ref))
2042
2043	0	T_z = evaluation_polynomial( ppoly_T, nterm, 1. )
2044	0	S_z = evaluation_polynomial( ppoly_S, nterm, 1. )
2045	0	drho_max = 0.5((dRdT_bot+dRdT_ref)(T_z-T_ref) + (dRdS_bot+dRdS_ref)*(S_z-S_ref))
2046
2047	0	if (drho_min >= 0.) then
2048	0	z = z0
2049	0	return
2050	0	elseif (drho_max == 0.) then
2051	0	z = 1.
2052	0	return
2053		endif
2054	0	if ( SIGN(1.,drho_min) == SIGN(1.,drho_max) ) then
2055	0	call MOM_error(FATAL, "drho_min is the same sign as dhro_max")
2056		endif
2057
2058	0	z = z0
2059	0	ztest = z0
2060	0	do iter = 1, CS%max_iter
2061		! Calculate quantities at the current nondimensional position
2062	0	a1 = 1.-z
2063	0	a2 = z
2064	0	dRdT_z = a1dRdT_top + a2dRdT_bot
2065	0	dRdS_z = a1dRdS_top + a2dRdS_bot
2066	0	T_z = evaluation_polynomial( ppoly_T, nterm, z )
2067	0	S_z = evaluation_polynomial( ppoly_S, nterm, z )
2068	0	drho = 0.5((dRdT_z+dRdT_ref)(T_z-T_ref) + (dRdS_z+dRdS_ref)*(S_z-S_ref))
2069
2070		! Check for convergence
2071	0	if (ABS(drho) <= CS%drho_tol) exit
2072		! Update bisection bracketing intervals
2073	0	if (drho < 0. .and. drho > drho_min) then
2074	0	drho_min = drho
2075	0	zmin = z
2076	0	elseif (drho > 0. .and. drho < drho_max) then
2077	0	drho_max = drho
2078	0	zmax = z
2079		endif
2080
2081		! Calculate a Newton step
2082	0	dT_dz = first_derivative_polynomial( ppoly_T, nterm, z )
2083	0	dS_dz = first_derivative_polynomial( ppoly_S, nterm, z )
2084		drho_dz = 0.5( (dRdT_diff(T_z - T_ref) + (dRdT_ref+dRdT_z)*dT_dz) + &
2085	0	(dRdS_diff(S_z - S_ref) + (dRdS_ref+dRdS_z)dS_dz) )
2086
2087	0	ztest = z - drho/drho_dz
2088		! Take a bisection if z falls out of [zmin,zmax]
2089	0	if (ztest < zmin .or. ztest > zmax) then
2090	0	if ( drho < 0. ) then
2091	0	ztest = 0.5*(z + zmax)
2092		else
2093	0	ztest = 0.5*(zmin + z)
2094		endif
2095		endif
2096
2097		! Test to ensure we haven't stalled out
2098	0	if ( abs(z-ztest) <= CS%x_tol ) exit
2099		! Reset for next iteration
2100	0	z = ztest
2101		enddo
2102
2103	0	end function find_neutral_pos_linear
2104
2105		!> Use the full equation of state to calculate the difference in locally referenced potential density. The derivatives
2106		!! in this case are not trivial to calculate, so instead we use a regula falsi method
2107	0	function find_neutral_pos_full( CS, z0, T_ref, S_ref, P_ref, P_top, P_bot, ppoly_T, ppoly_S ) result( z )
2108		type(neutral_diffusion_CS),intent(in) :: CS !< Control structure with parameters for this module
2109		real, intent(in) :: z0 !< Lower bound of position, also serves as the
2110		!! initial guess [nondim]
2111		real, intent(in) :: T_ref !< Temperature at the searched from interface [C ~> degC]
2112		real, intent(in) :: S_ref !< Salinity at the searched from interface [S ~> ppt]
2113		real, intent(in) :: P_ref !< Pressure at the searched from interface [R L2 T-2 ~> Pa]
2114		real, intent(in) :: P_top !< Pressure at top of layer being searched [R L2 T-2 ~> Pa]
2115		real, intent(in) :: P_bot !< Pressure at bottom of layer being searched [R L2 T-2 ~> Pa]
2116		real, dimension(:), intent(in) :: ppoly_T !< Coefficients of the polynomial reconstruction of T within
2117		!! the layer to be searched [C ~> degC]
2118		real, dimension(:), intent(in) :: ppoly_S !< Coefficients of the polynomial reconstruction of T within
2119		!! the layer to be searched [S ~> ppt]
2120		real :: z !< Position where drho = 0 [nondim]
2121		! Local variables
2122		integer :: iter
2123		integer :: nterm
2124
2125		real :: drho_a, drho_b, drho_c ! Density differences [R ~> kg m-3]
2126		real :: a, b, c ! Fractional positions [nondim]
2127		real :: Ta, Tb, Tc ! Temperatures [C ~> degC]
2128		real :: Sa, Sb, Sc ! Salinities [S ~> ppt]
2129		real :: Pa, Pb, Pc ! Pressures [R L2 T-2 ~> Pa]
2130		integer :: side
2131
2132	0	side = 0
2133		! Set the first two evaluation to the endpoints of the interval
2134	0	b = z0 ; c = 1
2135	0	nterm = SIZE(ppoly_T)
2136
2137		! Calculate drho at the minimum bound
2138	0	Tb = evaluation_polynomial( ppoly_T, nterm, b )
2139	0	Sb = evaluation_polynomial( ppoly_S, nterm, b )
2140	0	Pb = P_top(1.-b) + P_botb
2141	0	call calc_delta_rho_and_derivs(CS, Tb, Sb, Pb, T_ref, S_ref, P_ref, drho_b)
2142
2143		! Calculate drho at the maximum bound
2144	0	Tc = evaluation_polynomial( ppoly_T, nterm, 1. )
2145	0	Sc = evaluation_polynomial( ppoly_S, nterm, 1. )
2146	0	Pc = P_Bot
2147	0	call calc_delta_rho_and_derivs(CS, Tc, Sc, Pc, T_ref, S_ref, P_ref, drho_c)
2148
2149	0	if (drho_b >= 0.) then
2150	0	z = z0
2151	0	return
2152	0	elseif (drho_c == 0.) then
2153	0	z = 1.
2154	0	return
2155		endif
2156	0	if ( SIGN(1.,drho_b) == SIGN(1.,drho_c) ) then
2157	0	z = z0
2158	0	return
2159		endif
2160
2161	0	do iter = 1, CS%max_iter
2162		! Calculate new position and evaluate if we have converged
2163	0	a = (drho_bc - drho_cb)/(drho_b-drho_c)
2164	0	Ta = evaluation_polynomial( ppoly_T, nterm, a )
2165	0	Sa = evaluation_polynomial( ppoly_S, nterm, a )
2166	0	Pa = P_top(1.-a) + P_bota
2167	0	call calc_delta_rho_and_derivs(CS, Ta, Sa, Pa, T_ref, S_ref, P_ref, drho_a)
2168	0	if (ABS(drho_a) < CS%drho_tol) then
2169	0	z = a
2170	0	return
2171		endif
2172
2173	0	if (drho_a*drho_c > 0.) then
2174	0	if ( ABS(a-c)<CS%x_tol) then
2175	0	z = a
2176	0	return
2177		endif
2178	0	c = a ; drho_c = drho_a
2179	0	if (side == -1) drho_b = 0.5*drho_b
2180	0	side = -1
2181	0	elseif ( drho_b*drho_a > 0 ) then
2182	0	if ( ABS(a-b)<CS%x_tol) then
2183	0	z = a
2184	0	return
2185		endif
2186	0	b = a ; drho_b = drho_a
2187	0	if (side == 1) drho_c = 0.5*drho_c
2188	0	side = 1
2189		else
2190	0	z = a
2191	0	return
2192		endif
2193		enddo
2194
2195	0	z = a
2196
2197	0	end function find_neutral_pos_full
2198
2199		!> Calculate the difference in density between two points in a variety of ways
2200	0	subroutine calc_delta_rho_and_derivs(CS, T1, S1, p1_in, T2, S2, p2_in, drho, &
2201		drdt1_out, drds1_out, drdt2_out, drds2_out )
2202		type(neutral_diffusion_CS) :: CS !< Neutral diffusion control structure
2203		real, intent(in ) :: T1 !< Temperature at point 1 [C ~> degC]
2204		real, intent(in ) :: S1 !< Salinity at point 1 [S ~> ppt]
2205		real, intent(in ) :: p1_in !< Pressure at point 1 [R L2 T-2 ~> Pa]
2206		real, intent(in ) :: T2 !< Temperature at point 2 [C ~> degC]
2207		real, intent(in ) :: S2 !< Salinity at point 2 [S ~> ppt]
2208		real, intent(in ) :: p2_in !< Pressure at point 2 [R L2 T-2 ~> Pa]
2209		real, intent( out) :: drho !< Difference in density between the two points [R ~> kg m-3]
2210		real, optional, intent( out) :: dRdT1_out !< drho_dt at point 1 [R C-1 ~> kg m-3 degC-1]
2211		real, optional, intent( out) :: dRdS1_out !< drho_ds at point 1 [R S-1 ~> kg m-3 ppt-1]
2212		real, optional, intent( out) :: dRdT2_out !< drho_dt at point 2 [R C-1 ~> kg m-3 degC-1]
2213		real, optional, intent( out) :: dRdS2_out !< drho_ds at point 2 [R S-1 ~> kg m-3 ppt-1]
2214		! Local variables
2215		real :: rho1, rho2 ! Densities [R ~> kg m-3]
2216		real :: p1, p2, pmid ! Pressures [R L2 T-2 ~> Pa]
2217		real :: drdt1, drdt2 ! Partial derivatives of density with temperature [R C-1 ~> kg m-3 degC-1]
2218		real :: drds1, drds2 ! Partial derivatives of density with salinity [R S-1 ~> kg m-3 ppt-1]
2219
2220		! Use the same reference pressure or the in-situ pressure
2221	0	if (CS%ref_pres > 0.) then
2222	0	p1 = CS%ref_pres
2223	0	p2 = CS%ref_pres
2224		else
2225	0	p1 = p1_in
2226	0	p2 = p2_in
2227		endif
2228
2229		! Use the full linear equation of state to calculate the difference in density (expensive!)
2230	0	if (TRIM(CS%delta_rho_form) == 'full') then
2231	0	pmid = 0.5 * (p1 + p2)
2232	0	call calculate_density(T1, S1, pmid, rho1, CS%EOS)
2233	0	call calculate_density(T2, S2, pmid, rho2, CS%EOS)
2234	0	drho = rho1 - rho2
2235		! Use the density derivatives at the average of pressures and the differences in temperature
2236	0	elseif (TRIM(CS%delta_rho_form) == 'mid_pressure') then
2237	0	pmid = 0.5 * (p1 + p2)
2238	0	if (CS%ref_pres>=0) pmid = CS%ref_pres
2239	0	call calculate_density_derivs(T1, S1, pmid, drdt1, drds1, CS%EOS)
2240	0	call calculate_density_derivs(T2, S2, pmid, drdt2, drds2, CS%EOS)
2241	0	drho = delta_rho_from_derivs( T1, S1, p1, drdt1, drds1, T2, S2, p2, drdt2, drds2)
2242	0	elseif (TRIM(CS%delta_rho_form) == 'local_pressure') then
2243	0	call calculate_density_derivs(T1, S1, p1, drdt1, drds1, CS%EOS)
2244	0	call calculate_density_derivs(T2, S2, p2, drdt2, drds2, CS%EOS)
2245	0	drho = delta_rho_from_derivs( T1, S1, p1, drdt1, drds1, T2, S2, p2, drdt2, drds2)
2246		else
2247	0	call MOM_error(FATAL, "delta_rho_form is not recognized")
2248		endif
2249
2250	0	if (PRESENT(drdt1_out)) drdt1_out = drdt1
2251	0	if (PRESENT(drds1_out)) drds1_out = drds1
2252	0	if (PRESENT(drdt2_out)) drdt2_out = drdt2
2253	0	if (PRESENT(drds2_out)) drds2_out = drds2
2254
2255	0	end subroutine calc_delta_rho_and_derivs
2256
2257		!> Calculate delta rho from derivatives and gradients of properties
2258		!! \f$ \Delta \rho = \frac{1}{2}\left[ (\alpha_1 + \alpha_2)*(T_1-T_2) +
2259		!! (\beta_1 + \beta_2)*(S_1-S_2) +
2260		!! (\gamma^{-1}_1 + \gamma^{-1}_2)*(P_1-P_2) \right] \f$
2261	0	function delta_rho_from_derivs( T1, S1, P1, dRdT1, dRdS1, &
2262		T2, S2, P2, dRdT2, dRdS2 ) result (drho)
2263		real :: T1 !< Temperature at point 1 [C ~> degC]
2264		real :: S1 !< Salinity at point 1 [S ~> ppt]
2265		real :: P1 !< Pressure at point 1 [R L2 T-2 ~> Pa]
2266		real :: dRdT1 !< The partial derivative of density with temperature at point 1 [R C-1 ~> kg m-3 degC-1]
2267		real :: dRdS1 !< The partial derivative of density with salinity at point 1 [R S-1 ~> kg m-3 ppt-1]
2268		real :: T2 !< Temperature at point 2 [C ~> degC]
2269		real :: S2 !< Salinity at point 2 [S ~> ppt]
2270		real :: P2 !< Pressure at point 2 [R L2 T-2 ~> Pa]
2271		real :: dRdT2 !< The partial derivative of density with temperature at point 2 [R C-1 ~> kg m-3 degC-1]
2272		real :: dRdS2 !< The partial derivative of density with salinity at point 2 [R S-1 ~> kg m-3 ppt-1]
2273		! Local variables
2274		real :: drho ! The density difference [R ~> kg m-3]
2275
2276	0	drho = 0.5 * ( (dRdT1+dRdT2)(T1-T2) + (dRdS1+dRdS2)(S1-S2))
2277
2278	0	end function delta_rho_from_derivs
2279
2280		!> Converts non-dimensional position within a layer to absolute position (for debugging)
2281	0	function absolute_position(n,ns,Pint,Karr,NParr,k_surface)
2282		integer, intent(in) :: n !< Number of levels
2283		integer, intent(in) :: ns !< Number of neutral surfaces
2284		real, intent(in) :: Pint(n+1) !< Position of interfaces [R L2 T-2 ~> Pa] or other units
2285		integer, intent(in) :: Karr(ns) !< Index of interface above position
2286		real, intent(in) :: NParr(ns) !< Non-dimensional position within layer Karr(:) [nondim]
2287		integer, intent(in) :: k_surface !< k-interface to query
2288		real :: absolute_position !< The absolute position of a location [R L2 T-2 ~> Pa]
2289		!! or other units following Pint
2290		! Local variables
2291		integer :: k
2292
2293	0	k = Karr(k_surface)
2294	0	if (k>n) stop 'absolute_position: k>nk is out of bounds!'
2295	0	absolute_position = Pint(k) + NParr(k_surface) * ( Pint(k+1) - Pint(k) )
2296
2297	0	end function absolute_position
2298
2299		!> Converts non-dimensional positions within layers to absolute positions (for debugging)
2300	0	function absolute_positions(n,ns,Pint,Karr,NParr)
2301		integer, intent(in) :: n !< Number of levels
2302		integer, intent(in) :: ns !< Number of neutral surfaces
2303		real, intent(in) :: Pint(n+1) !< Position of interface [R L2 T-2 ~> Pa] or other units
2304		integer, intent(in) :: Karr(ns) !< Indexes of interfaces about positions
2305		real, intent(in) :: NParr(ns) !< Non-dimensional positions within layers Karr(:) [nondim]
2306
2307		real, dimension(ns) :: absolute_positions !< Absolute positions [R L2 T-2 ~> Pa]
2308		!! or other units following Pint
2309
2310		! Local variables
2311		integer :: k_surface
2312
2313	0	do k_surface = 1, ns
2314	0	absolute_positions(k_surface) = absolute_position(n,ns,Pint,Karr,NParr,k_surface)
2315		enddo
2316
2317		end function absolute_positions
2318
2319		!> Returns a single column of neutral diffusion fluxes of a tracer.
2320	0	subroutine neutral_surface_flux(nk, nsurf, deg, hl, hr, Tl, Tr, PiL, PiR, KoL, KoR, &
2321	0	hEff, Flx, continuous, h_neglect, remap_CS, h_neglect_edge, &
2322	0	coeff_l, coeff_r)
2323		integer, intent(in) :: nk !< Number of levels
2324		integer, intent(in) :: nsurf !< Number of neutral surfaces
2325		integer, intent(in) :: deg !< Degree of polynomial reconstructions
2326		real, dimension(nk), intent(in) :: hl !< Left-column layer thickness [H ~> m or kg m-2]
2327		real, dimension(nk), intent(in) :: hr !< Right-column layer thickness [H ~> m or kg m-2]
2328		real, dimension(nk), intent(in) :: Tl !< Left-column layer tracer in arbitrary concentration
2329		!! units (e.g. [C ~> degC] for temperature)
2330		real, dimension(nk), intent(in) :: Tr !< Right-column layer tracer in arbitrary concentration
2331		!! units (e.g. [C ~> degC] for temperature)
2332		real, dimension(nsurf), intent(in) :: PiL !< Fractional position of neutral surface
2333		!! within layer KoL of left column [nondim]
2334		real, dimension(nsurf), intent(in) :: PiR !< Fractional position of neutral surface
2335		!! within layer KoR of right column [nondim]
2336		integer, dimension(nsurf), intent(in) :: KoL !< Index of first left interface above neutral surface
2337		integer, dimension(nsurf), intent(in) :: KoR !< Index of first right interface above neutral surface
2338		real, dimension(nsurf-1), intent(in) :: hEff !< Effective thickness between two neutral
2339		!! surfaces [H ~> m or kg m-2]
2340		real, dimension(nsurf-1), intent(inout) :: Flx !< Flux of tracer between pairs of neutral layers
2341		!! in units (conc H or conc H L2) that depend on
2342		!! the presence and units of coeff_l and coeff_r.
2343		!! If the tracer is temperature, this could have
2344		!! units of [C H ~> degC m or degC kg m-2] or
2345		!! [C H L2 ~> degC m3 or degC kg] if coeff_l has
2346		!! units of [L2 ~> m2]
2347		logical, intent(in) :: continuous !< True if using continuous reconstruction
2348		real, intent(in) :: h_neglect !< A negligibly small width for the purpose
2349		!! of cell reconstructions [H ~> m or kg m-2]
2350		type(remapping_CS), optional, intent(in) :: remap_CS !< Remapping control structure used
2351		!! to create sublayers
2352		real, optional, intent(in) :: h_neglect_edge !< A negligibly small width used for edge value
2353		!! calculations if continuous is false [H ~> m or kg m-2]
2354		real, dimension(nk+1), optional, intent(in) :: coeff_l !< Left-column diffusivity [L2 ~> m2] or [nondim]
2355		real, dimension(nk+1), optional, intent(in) :: coeff_r !< Right-column diffusivity [L2 ~> m2] or [nondim]
2356
2357		! Local variables
2358		integer :: k_sublayer, klb, klt, krb, krt
2359		real :: T_right_sub, T_left_sub ! Tracer concentrations averaged over sub-intervals in the right and left
2360		! columns in arbitrary concentration units (e.g. [C ~> degC] for temperature).
2361		real :: T_right_layer, T_left_layer ! Tracer concentrations averaged over layers in the right and left
2362		! columns in arbitrary concentration units (e.g. [C ~> degC] for temperature).
2363		real :: T_right_top, T_right_bottom, T_right_top_int, T_right_bot_int ! Tracer concentrations
2364		! at various positions in the right column in arbitrary
2365		! concentration units (e.g. [C ~> degC] for temperature).
2366		real :: T_left_top, T_left_bottom, T_left_top_int, T_left_bot_int ! Tracer concentrations
2367		! at various positions in the left column in arbitrary
2368		! concentration units (e.g. [C ~> degC] for temperature).
2369		real :: dT_layer, dT_ave, dT_sublayer ! Differences in vertically averaged tracer concentrations
2370		! over various portions of the right and left columns in arbitrary
2371		! concentration units (e.g. [C ~> degC] for temperature).
2372		real :: dT_top, dT_bottom, dT_top_int, dT_bot_int ! Differences in tracer concentrations
2373		! at various positions between the right and left columns in arbitrary
2374		! concentration units (e.g. [C ~> degC] for temperature).
2375		real :: khtr_ave ! An averaged diffusivity in normalized units [nondim] if coeff_l and coeff_r are
2376		! absent or in units copied from coeff_l and coeff_r [L2 ~> m2] or [nondim]
2377	0	real, dimension(nk+1) :: Til !< Left-column interface tracer in arbitrary concentration
2378		!! units (e.g. [C ~> degC] for temperature)
2379	0	real, dimension(nk+1) :: Tir !< Right-column interface tracer in arbitrary concentration
2380		!! units (e.g. [C ~> degC] for temperature)
2381	0	real, dimension(nk) :: aL_l !< Left-column left edge value of tracer in arbitrary concentration
2382		!! units (e.g. [C ~> degC] for temperature)
2383	0	real, dimension(nk) :: aR_l !< Left-column right edge value of tracer in arbitrary concentration
2384		!! units (e.g. [C ~> degC] for temperature)
2385	0	real, dimension(nk) :: aL_r !< Right-column left edge value of tracer in arbitrary concentration
2386		!! units (e.g. [C ~> degC] for temperature)
2387	0	real, dimension(nk) :: aR_r !< Right-column right edge value of tracer in arbitrary concentration
2388		!! units (e.g. [C ~> degC] for temperature)
2389		! Discontinuous reconstruction
2390		integer :: iMethod
2391	0	real, dimension(nk,2) :: Tid_l !< Left-column interface tracer in arbitrary concentration
2392		!! units (e.g. [C ~> degC] for temperature)
2393	0	real, dimension(nk,2) :: Tid_r !< Right-column interface tracer in arbitrary concentration
2394		!! units (e.g. [C ~> degC] for temperature)
2395	0	real, dimension(nk,deg+1) :: ppoly_r_coeffs_l ! Coefficients of the polynomial descriptions of
2396		! sub-gridscale tracer concentrations in the left column, in arbitrary
2397		! concentration units (e.g. [C ~> degC] for temperature)
2398	0	real, dimension(nk,deg+1) :: ppoly_r_coeffs_r ! Coefficients of the polynomial descriptions of
2399		! sub-gridscale tracer concentrations in the right column, in arbitrary
2400		! concentration units (e.g. [C ~> degC] for temperature)
2401	0	real, dimension(nk,deg+1) :: ppoly_r_S_l ! Reconstruction slopes that are unused here, in units of a vertical
2402		! gradient, which for temperature would be [C H-1 ~> degC m-1 or degC m2 kg-1].
2403	0	real, dimension(nk,deg+1) :: ppoly_r_S_r ! Reconstruction slopes that are unused here, in units of a vertical
2404		! gradient, which for temperature would be [C H-1 ~> degC m-1 or degC m2 kg-1].
2405		logical :: down_flux, tapering
2406
2407	0	tapering = .false.
2408	0	if (present(coeff_l) .and. present(coeff_r)) tapering = .true.
2409	0	khtr_ave = 1.0
2410
2411		! Setup reconstruction edge values
2412	0	if (continuous) then
2413	0	call interface_scalar(nk, hl, Tl, Til, 2, h_neglect)
2414	0	call interface_scalar(nk, hr, Tr, Tir, 2, h_neglect)
2415	0	call ppm_left_right_edge_values(nk, Tl, Til, aL_l, aR_l)
2416	0	call ppm_left_right_edge_values(nk, Tr, Tir, aL_r, aR_r)
2417		else
2418	0	ppoly_r_coeffs_l(:,:) = 0.
2419	0	ppoly_r_coeffs_r(:,:) = 0.
2420	0	Tid_l(:,:) = 0.
2421	0	Tid_r(:,:) = 0.
2422
2423		call build_reconstructions_1d( remap_CS, nk, hl, Tl, ppoly_r_coeffs_l, Tid_l, &
2424	0	ppoly_r_S_l, iMethod, h_neglect, h_neglect_edge )
2425		call build_reconstructions_1d( remap_CS, nk, hr, Tr, ppoly_r_coeffs_r, Tid_r, &
2426	0	ppoly_r_S_r, iMethod, h_neglect, h_neglect_edge )
2427		endif
2428
2429	0	do k_sublayer = 1, nsurf-1
2430	0	if (hEff(k_sublayer) == 0.) then
2431	0	Flx(k_sublayer) = 0.
2432		else
2433	0	if (tapering) then
2434	0	klb = KoL(k_sublayer+1)
2435	0	klt = KoL(k_sublayer)
2436	0	krb = KoR(k_sublayer+1)
2437	0	krt = KoR(k_sublayer)
2438		! these are added in this order to preserve vertically-uniform diffusivity answers
2439	0	khtr_ave = 0.25 * ((coeff_l(klb) + coeff_l(klt)) + (coeff_r(krb) + coeff_r(krt)))
2440		endif
2441	0	if (continuous) then
2442	0	klb = KoL(k_sublayer+1)
2443	0	T_left_bottom = ( 1. - PiL(k_sublayer+1) ) * Til(klb) + PiL(k_sublayer+1) * Til(klb+1)
2444	0	klt = KoL(k_sublayer)
2445	0	T_left_top = ( 1. - PiL(k_sublayer) ) * Til(klt) + PiL(k_sublayer) * Til(klt+1)
2446		T_left_layer = ppm_ave(PiL(k_sublayer), PiL(k_sublayer+1) + real(klb-klt), &
2447	0	aL_l(klt), aR_l(klt), Tl(klt))
2448
2449	0	krb = KoR(k_sublayer+1)
2450	0	T_right_bottom = ( 1. - PiR(k_sublayer+1) ) * Tir(krb) + PiR(k_sublayer+1) * Tir(krb+1)
2451	0	krt = KoR(k_sublayer)
2452	0	T_right_top = ( 1. - PiR(k_sublayer) ) * Tir(krt) + PiR(k_sublayer) * Tir(krt+1)
2453		T_right_layer = ppm_ave(PiR(k_sublayer), PiR(k_sublayer+1) + real(krb-krt), &
2454	0	aL_r(krt), aR_r(krt), Tr(krt))
2455	0	dT_top = T_right_top - T_left_top
2456	0	dT_bottom = T_right_bottom - T_left_bottom
2457	0	dT_ave = 0.5 * ( dT_top + dT_bottom )
2458	0	dT_layer = T_right_layer - T_left_layer
2459	0	if (signum(1.,dT_top) * signum(1.,dT_bottom) <= 0. .or. signum(1.,dT_ave) * signum(1.,dT_layer) <= 0.) then
2460	0	dT_ave = 0.
2461		else
2462	0	dT_ave = dT_layer
2463		endif
2464	0	Flx(k_sublayer) = dT_ave * hEff(k_sublayer) * khtr_ave
2465		else ! Discontinuous reconstruction
2466		! Calculate tracer values on left and right side of the neutral surface
2467		call neutral_surface_T_eval(nk, nsurf, k_sublayer, KoL, PiL, Tl, Tid_l, deg, iMethod, &
2468		ppoly_r_coeffs_l, T_left_top, T_left_bottom, T_left_sub, &
2469	0	T_left_top_int, T_left_bot_int, T_left_layer)
2470		call neutral_surface_T_eval(nk, nsurf, k_sublayer, KoR, PiR, Tr, Tid_r, deg, iMethod, &
2471		ppoly_r_coeffs_r, T_right_top, T_right_bottom, T_right_sub, &
2472	0	T_right_top_int, T_right_bot_int, T_right_layer)
2473
2474	0	dT_top = T_right_top - T_left_top
2475	0	dT_bottom = T_right_bottom - T_left_bottom
2476	0	dT_sublayer = T_right_sub - T_left_sub
2477	0	dT_top_int = T_right_top_int - T_left_top_int
2478	0	dT_bot_int = T_right_bot_int - T_left_bot_int
2479		! Enforcing the below criterion incorrectly zero out fluxes
2480		!dT_layer = T_right_layer - T_left_layer
2481
2482		down_flux = dT_top <= 0. .and. dT_bottom <= 0. .and. &
2483		dT_sublayer <= 0. .and. dT_top_int <= 0. .and. &
2484	0	dT_bot_int <= 0.
2485		down_flux = down_flux .or. &
2486		(dT_top >= 0. .and. dT_bottom >= 0. .and. &
2487		dT_sublayer >= 0. .and. dT_top_int >= 0. .and. &
2488	0	dT_bot_int >= 0.)
2489	0	if (down_flux) then
2490	0	Flx(k_sublayer) = dT_sublayer * hEff(k_sublayer) * khtr_ave
2491		else
2492	0	Flx(k_sublayer) = 0.
2493		endif
2494		endif
2495		endif
2496		enddo
2497
2498	0	end subroutine neutral_surface_flux
2499
2500		!> Evaluate various parts of the reconstructions to calculate gradient-based flux limiter
2501	0	subroutine neutral_surface_T_eval(nk, ns, k_sub, Ks, Ps, T_mean, T_int, deg, iMethod, T_poly, &
2502		T_top, T_bot, T_sub, T_top_int, T_bot_int, T_layer)
2503		integer, intent(in ) :: nk !< Number of cell averages
2504		integer, intent(in ) :: ns !< Number of neutral surfaces
2505		integer, intent(in ) :: k_sub !< Index of current neutral layer
2506		integer, dimension(ns), intent(in ) :: Ks !< List of the layers associated with each neutral surface
2507		real, dimension(ns), intent(in ) :: Ps !< List of the positions within a layer of each surface [nondim]
2508		real, dimension(nk), intent(in ) :: T_mean !< Layer average of tracer in arbitrary concentration
2509		!! units (e.g. [C ~> degC] for temperature)
2510		real, dimension(nk,2), intent(in ) :: T_int !< Layer interface values of tracer from reconstruction
2511		!! in concentration units (e.g. [C ~> degC] for temperature)
2512		integer, intent(in ) :: deg !< Degree of reconstruction polynomial (e.g. 1 is linear)
2513		integer, intent(in ) :: iMethod !< Method of integration to use
2514		real, dimension(nk,deg+1), intent(in ) :: T_poly !< Coefficients of polynomial reconstructions in arbitrary
2515		!! concentration units (e.g. [C ~> degC] for temperature)
2516		real, intent( out) :: T_top !< Tracer value at top (across discontinuity if necessary) in
2517		!! concentration units (e.g. [C ~> degC] for temperature)
2518		real, intent( out) :: T_bot !< Tracer value at bottom (across discontinuity if necessary)
2519		!! in concentration units (e.g. [C ~> degC] for temperature)
2520		real, intent( out) :: T_sub !< Average of the tracer value over the sublayer in arbitrary
2521		!! concentration units (e.g. [C ~> degC] for temperature)
2522		real, intent( out) :: T_top_int !< Tracer value at the top interface of a neutral layer in
2523		!! concentration units (e.g. [C ~> degC] for temperature)
2524		real, intent( out) :: T_bot_int !< Tracer value at the bottom interface of a neutral layer in
2525		!! concentration units (e.g. [C ~> degC] for temperature)
2526		real, intent( out) :: T_layer !< Cell-average tracer concentration in a layer that
2527		!! the reconstruction belongs to in concentration
2528		!! units (e.g. [C ~> degC] for temperature)
2529
2530		integer :: kl, ks_top, ks_bot
2531
2532	0	ks_top = k_sub
2533	0	ks_bot = k_sub + 1
2534	0	if ( Ks(ks_top) /= Ks(ks_bot) ) then
2535	0	call MOM_error(FATAL, "Neutral surfaces span more than one layer")
2536		endif
2537	0	kl = Ks(k_sub)
2538		! First if the neutral surfaces spans the entirety of a cell, then do not search across the discontinuity
2539	0	if ( (Ps(ks_top) == 0.) .and. (Ps(ks_bot) == 1.)) then
2540	0	T_top = T_int(kl,1)
2541	0	T_bot = T_int(kl,2)
2542		else
2543		! Search across potential discontinuity at top
2544	0	if ( (kl > 1) .and. (Ps(ks_top) == 0.) ) then
2545	0	T_top = T_int(kl-1,2)
2546		else
2547	0	T_top = evaluation_polynomial( T_poly(kl,:), deg+1, Ps(ks_top) )
2548		endif
2549		! Search across potential discontinuity at bottom
2550	0	if ( (kl < nk) .and. (Ps(ks_bot) == 1.) ) then
2551	0	T_bot = T_int(kl+1,1)
2552		else
2553	0	T_bot = evaluation_polynomial( T_poly(kl,:), deg+1, Ps(ks_bot) )
2554		endif
2555		endif
2556	0	T_sub = average_value_ppoly(nk, T_mean, T_int, T_poly, iMethod, kl, Ps(ks_top), Ps(ks_bot))
2557	0	T_top_int = evaluation_polynomial( T_poly(kl,:), deg+1, Ps(ks_top))
2558	0	T_bot_int = evaluation_polynomial( T_poly(kl,:), deg+1, Ps(ks_bot))
2559	0	T_layer = T_mean(kl)
2560
2561	0	end subroutine neutral_surface_T_eval
2562
2563		!> Discontinuous PPM reconstructions of the left/right edge values within a cell
2564	0	subroutine ppm_left_right_edge_values(nk, Tl, Ti, aL, aR)
2565		integer, intent(in) :: nk !< Number of levels
2566		real, dimension(nk), intent(in) :: Tl !< Layer tracer (conc, e.g. degC) in arbitrary units [A ~> a]
2567		real, dimension(nk+1), intent(in) :: Ti !< Interface tracer (conc, e.g. degC) in arbitrary units [A ~> a]
2568		real, dimension(nk), intent(inout) :: aL !< Left edge value of tracer (conc, e.g. degC)
2569		!! in the same arbitrary units as Tl and Ti [A ~> a]
2570		real, dimension(nk), intent(inout) :: aR !< Right edge value of tracer (conc, e.g. degC)
2571		!! in the same arbitrary units as Tl and Ti [A ~> a]
2572
2573		integer :: k
2574		! Setup reconstruction edge values
2575	0	do k = 1, nk
2576	0	aL(k) = Ti(k)
2577	0	aR(k) = Ti(k+1)
2578	0	if ( signum(1., aR(k) - Tl(k))*signum(1., Tl(k) - aL(k)) <= 0.0 ) then
2579	0	aL(k) = Tl(k)
2580	0	aR(k) = Tl(k)
2581	0	elseif ( sign(3., aR(k) - aL(k)) * ( (Tl(k) - aL(k)) + (Tl(k) - aR(k))) > abs(aR(k) - aL(k)) ) then
2582	0	aL(k) = Tl(k) + 2.0 * ( Tl(k) - aR(k) )
2583	0	elseif ( sign(3., aR(k) - aL(k)) * ( (Tl(k) - aL(k)) + (Tl(k) - aR(k))) < -abs(aR(k) - aL(k)) ) then
2584	0	aR(k) = Tl(k) + 2.0 * ( Tl(k) - aL(k) )
2585		endif
2586		enddo
2587	0	end subroutine ppm_left_right_edge_values
2588
2589		!> Returns true if unit tests of neutral_diffusion functions fail. Otherwise returns false.
2590	0	logical function neutral_diffusion_unit_tests(verbose)
2591		logical, intent(in) :: verbose !< If true, write results to stdout
2592
2593		neutral_diffusion_unit_tests = .false. .or. &
2594	0	ndiff_unit_tests_continuous(verbose) .or. ndiff_unit_tests_discontinuous(verbose)
2595
2596	0	end function neutral_diffusion_unit_tests
2597
2598		!> Returns true if unit tests of neutral_diffusion functions fail. Otherwise returns false.
2599	0	logical function ndiff_unit_tests_continuous(verbose)
2600		logical, intent(in) :: verbose !< If true, write results to stdout
2601		! Local variables
2602		integer, parameter :: nk = 4
2603		real, dimension(nk+1) :: Tio ! Test interface temperatures [degC]
2604		real, dimension(2*nk+2) :: PiLRo, PiRLo ! Fractional test positions [nondim]
2605		integer, dimension(2*nk+2) :: KoL, KoR ! Test indexes
2606		real, dimension(2*nk+1) :: hEff ! Test positions in arbitrary units [arbitrary]
2607		real, dimension(2*nk+1) :: Flx ! Test flux in the arbitrary units of hEff times [degC]
2608		logical :: v
2609		real :: h_neglect ! A negligible thickness in arbitrary units [arbitrary]
2610
2611	0	h_neglect = 1.0e-30
2612
2613	0	v = verbose
2614
2615	0	ndiff_unit_tests_continuous = .false. ! Normally return false
2616	0	write(stdout,*) '==== MOM_neutral_diffusion: ndiff_unit_tests_continuous ='
2617
2618		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2619	0	test_fv_diff(v,1.,1.,1., 0.,1.,2., 1., 'FV: Straight line on uniform grid')
2620		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2621	0	test_fv_diff(v,1.,1.,0., 0.,4.,8., 7., 'FV: Vanished right cell')
2622		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2623	0	test_fv_diff(v,0.,1.,1., 0.,4.,8., 7., 'FV: Vanished left cell')
2624		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2625	0	test_fv_diff(v,1.,2.,4., 0.,3.,9., 4., 'FV: Stretched grid')
2626		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2627	0	test_fv_diff(v,2.,0.,2., 0.,1.,2., 0., 'FV: Vanished middle cell')
2628		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2629	0	test_fv_diff(v,0.,1.,0., 0.,1.,2., 2., 'FV: Vanished on both sides')
2630		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2631	0	test_fv_diff(v,1.,0.,0., 0.,1.,2., 0., 'FV: Two vanished cell sides')
2632		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2633	0	test_fv_diff(v,0.,0.,0., 0.,1.,2., 0., 'FV: All vanished cells')
2634
2635		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2636	0	test_fvlsq_slope(v,1.,1.,1., 0.,1.,2., 1., 'LSQ: Straight line on uniform grid')
2637		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2638	0	test_fvlsq_slope(v,1.,1.,0., 0.,1.,2., 1., 'LSQ: Vanished right cell')
2639		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2640	0	test_fvlsq_slope(v,0.,1.,1., 0.,1.,2., 1., 'LSQ: Vanished left cell')
2641		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2642	0	test_fvlsq_slope(v,1.,2.,4., 0.,3.,9., 2., 'LSQ: Stretched grid')
2643		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2644	0	test_fvlsq_slope(v,1.,0.,1., 0.,1.,2., 2., 'LSQ: Vanished middle cell')
2645		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2646	0	test_fvlsq_slope(v,0.,1.,0., 0.,1.,2., 0., 'LSQ: Vanished on both sides')
2647		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2648	0	test_fvlsq_slope(v,1.,0.,0., 0.,1.,2., 0., 'LSQ: Two vanished cell sides')
2649		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2650	0	test_fvlsq_slope(v,0.,0.,0., 0.,1.,2., 0., 'LSQ: All vanished cells')
2651
2652	0	call interface_scalar(4, (/10.,10.,10.,10./), (/24.,18.,12.,6./), Tio, 1, h_neglect)
2653		!ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2654		! test_data1d(5, Tio, (/27.,21.,15.,9.,3./), 'Linear profile, interface temperatures')
2655		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2656	0	test_data1d(v,5, Tio, (/24.,22.5,15.,7.5,6./), 'Linear profile, linear interface temperatures')
2657	0	call interface_scalar(4, (/10.,10.,10.,10./), (/24.,18.,12.,6./), Tio, 2, h_neglect)
2658		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2659	0	test_data1d(v,5, Tio, (/24.,22.,15.,8.,6./), 'Linear profile, PPM interface temperatures')
2660
2661		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2662	0	test_ifndp(v,-1.0, 0., 1.0, 1.0, 0.5, 'Check mid-point')
2663		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2664	0	test_ifndp(v, 0.0, 0., 1.0, 1.0, 0.0, 'Check bottom')
2665		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2666	0	test_ifndp(v, 0.1, 0., 1.1, 1.0, 0.0, 'Check below')
2667		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2668	0	test_ifndp(v,-1.0, 0., 0.0, 1.0, 1.0, 'Check top')
2669		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2670	0	test_ifndp(v,-1.0, 0., -0.1, 1.0, 1.0, 'Check above')
2671		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2672	0	test_ifndp(v,-1.0, 0., 3.0, 1.0, 0.25, 'Check 1/4')
2673		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2674	0	test_ifndp(v,-3.0, 0., 1.0, 1.0, 0.75, 'Check 3/4')
2675		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2676	0	test_ifndp(v, 1.0, 0., 1.0, 1.0, 0.0, 'Check dRho=0 below')
2677		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2678	0	test_ifndp(v,-1.0, 0., -1.0, 1.0, 1.0, 'Check dRho=0 above')
2679		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2680	0	test_ifndp(v, 0.0, 0., 0.0, 1.0, 0.5, 'Check dRho=0 mid')
2681		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. &
2682	0	test_ifndp(v,-2.0, .5, 5.0, 0.5, 0.5, 'Check dP=0')
2683
2684		! Identical columns
2685		call find_neutral_surface_positions_continuous(3, &
2686		(/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S
2687		(/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
2688		(/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Right positions, T and S
2689		(/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
2690	0	PiLRo, PiRLo, KoL, KoR, hEff)
2691		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, KoL, KoR, PiLRo, PiRLo, hEff, &
2692		(/1,1,2,2,3,3,3,3/), & ! KoL
2693		(/1,1,2,2,3,3,3,3/), & ! KoR
2694		(/0.,0.,0.,0.,0.,0.,1.,1./), & ! pL
2695		(/0.,0.,0.,0.,0.,0.,1.,1./), & ! pR
2696		(/0.,10.,0.,10.,0.,10.,0./), & ! hEff
2697	0	'Identical columns')
2698		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_data1d(v, 8, &
2699		absolute_positions(3, 8, (/0.,10.,20.,30./), KoL, PiLRo), &
2700	0	(/0.,0.,10.,10.,20.,20.,30.,30./), '... left positions')
2701		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_data1d(v, 8, &
2702		absolute_positions(3, 8, (/0.,10.,20.,30./), KoR, PiRLo), &
2703	0	(/0.,0.,10.,10.,20.,20.,30.,30./), '... right positions')
2704		call neutral_surface_flux(3, 2*3+2, 2, (/10.,10.,10./), (/10.,10.,10./), & ! nk, hL, hR
2705		(/20.,16.,12./), (/20.,16.,12./), & ! Tl, Tr
2706	0	PiLRo, PiRLo, KoL, KoR, hEff, Flx, .true., h_neglect)
2707		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_data1d(v, 7, Flx, &
2708	0	(/0.,0.,0.,0.,0.,0.,0./), 'Identical columns, rho flux (=0)')
2709		call neutral_surface_flux(3, 2*3+2, 2, (/10.,10.,10./), (/10.,10.,10./), & ! nk, hL, hR
2710		(/-1.,-1.,-1./), (/1.,1.,1./), & ! Sl, Sr
2711	0	PiLRo, PiRLo, KoL, KoR, hEff, Flx, .true., h_neglect)
2712		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_data1d(v, 7, Flx, &
2713	0	(/0.,20.,0.,20.,0.,20.,0./), 'Identical columns, S flux')
2714
2715		! Right column slightly cooler than left
2716		call find_neutral_surface_positions_continuous(3, &
2717		(/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S
2718		(/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
2719		(/0.,10.,20.,30./), (/20.,16.,12.,8./), (/0.,0.,0.,0./), & ! Right positions, T and S
2720		(/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
2721	0	PiLRo, PiRLo, KoL, KoR, hEff)
2722		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, KoL, KoR, PiLRo, PiRLo, hEff, &
2723		(/1,1,2,2,3,3,3,3/), & ! kL
2724		(/1,1,1,2,2,3,3,3/), & ! kR
2725		(/0.,0.5,0.,0.5,0.,0.5,1.,1./), & ! pL
2726		(/0.,0.,0.5,0.,0.5,0.,0.5,1./), & ! pR
2727		(/0.,5.,5.,5.,5.,5.,0./), & ! hEff
2728	0	'Right column slightly cooler')
2729		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_data1d(v, 8, &
2730		absolute_positions(3, 8, (/0.,10.,20.,30./), KoL, PiLRo), &
2731	0	(/0.,5.,10.,15.,20.,25.,30.,30./), '... left positions')
2732		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_data1d(v, 8, &
2733		absolute_positions(3, 8, (/0.,10.,20.,30./), KoR, PiRLo), &
2734	0	(/0.,0.,5.,10.,15.,20.,25.,30./), '... right positions')
2735
2736		! Right column slightly warmer than left
2737		call find_neutral_surface_positions_continuous(3, &
2738		(/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S
2739		(/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
2740		(/0.,10.,20.,30./), (/24.,20.,16.,12./), (/0.,0.,0.,0./), & ! Right positions, T and S
2741		(/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
2742	0	PiLRo, PiRLo, KoL, KoR, hEff)
2743		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, KoL, KoR, PiLRo, PiRLo, hEff, &
2744		(/1,1,1,2,2,3,3,3/), & ! kL
2745		(/1,1,2,2,3,3,3,3/), & ! kR
2746		(/0.,0.,0.5,0.,0.5,0.,0.5,1./), & ! pL
2747		(/0.,0.5,0.,0.5,0.,0.5,1.,1./), & ! pR
2748		(/0.,5.,5.,5.,5.,5.,0./), & ! hEff
2749	0	'Right column slightly warmer')
2750
2751		! Right column somewhat cooler than left
2752		call find_neutral_surface_positions_continuous(3, &
2753		(/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S
2754		(/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
2755		(/0.,10.,20.,30./), (/16.,12.,8.,4./), (/0.,0.,0.,0./), & ! Right positions, T and S
2756		(/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
2757	0	PiLRo, PiRLo, KoL, KoR, hEff)
2758		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, KoL, KoR, PiLRo, PiRLo, hEff, &
2759		(/1,2,2,3,3,3,3,3/), & ! kL
2760		(/1,1,1,1,2,2,3,3/), & ! kR
2761		(/0.,0.,0.5,0.,0.5,1.,1.,1./), & ! pL
2762		(/0.,0.,0.,0.5,0.,0.5,0.,1./), & ! pR
2763		(/0.,0.,5.,5.,5.,0.,0./), & ! hEff
2764	0	'Right column somewhat cooler')
2765
2766		! Right column much colder than left with no overlap
2767		call find_neutral_surface_positions_continuous(3, &
2768		(/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S
2769		(/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
2770		(/0.,10.,20.,30./), (/9.,7.,5.,3./), (/0.,0.,0.,0./), & ! Right positions, T and S
2771		(/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
2772	0	PiLRo, PiRLo, KoL, KoR, hEff)
2773		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, KoL, KoR, PiLRo, PiRLo, hEff, &
2774		(/1,2,3,3,3,3,3,3/), & ! kL
2775		(/1,1,1,1,1,2,3,3/), & ! kR
2776		(/0.,0.,0.,1.,1.,1.,1.,1./), & ! pL
2777		(/0.,0.,0.,0.,0.,0.,0.,1./), & ! pR
2778		(/0.,0.,0.,0.,0.,0.,0./), & ! hEff
2779	0	'Right column much cooler')
2780
2781		! Right column with mixed layer
2782		call find_neutral_surface_positions_continuous(3, &
2783		(/0.,10.,20.,30./), (/22.,18.,14.,10./), (/0.,0.,0.,0./), & ! Left positions, T and S
2784		(/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
2785		(/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), & ! Right positions, T and S
2786		(/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
2787	0	PiLRo, PiRLo, KoL, KoR, hEff)
2788		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, KoL, KoR, PiLRo, PiRLo, hEff, &
2789		(/1,2,3,3,3,3,3,3/), & ! kL
2790		(/1,1,1,1,2,3,3,3/), & ! kR
2791		(/0.,0.,0.,0.,0.,1.,1.,1./), & ! pL
2792		(/0.,0.,0.,0.,0.,0.,0.,1./), & ! pR
2793		(/0.,0.,0.,0.,10.,0.,0./), & ! hEff
2794	0	'Right column with mixed layer')
2795
2796		! Identical columns with mixed layer
2797		call find_neutral_surface_positions_continuous(3, &
2798		(/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), & ! Left positions, T and S
2799		(/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
2800		(/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), & ! Right positions, T and S
2801		(/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
2802	0	PiLRo, PiRLo, KoL, KoR, hEff)
2803		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, KoL, KoR, PiLRo, PiRLo, hEff, &
2804		(/1,1,2,2,3,3,3,3/), & ! kL
2805		(/1,1,2,2,3,3,3,3/), & ! kR
2806		(/0.,0.,0.,0.,0.,0.,1.,1./), & ! pL
2807		(/0.,0.,0.,0.,0.,0.,1.,1./), & ! pR
2808		(/0.,10.,0.,10.,0.,10.,0./), & ! hEff
2809	0	'Identical columns with mixed layer')
2810
2811		! Right column with unstable mixed layer
2812		call find_neutral_surface_positions_continuous(3, &
2813		(/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), & ! Left positions, T and S
2814		(/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
2815		(/0.,10.,20.,30./), (/10.,14.,12.,4./), (/0.,0.,0.,0./), & ! Right positions, T and S
2816		(/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
2817	0	PiLRo, PiRLo, KoL, KoR, hEff)
2818		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, KoL, KoR, PiLRo, PiRLo, hEff, &
2819		(/1,2,3,3,3,3,3,3/), & ! kL
2820		(/1,1,1,2,3,3,3,3/), & ! kR
2821		(/0.,0.,0.,0.,0.,0.,.75,1./), & ! pL
2822		(/0.,0.,0.,0.,0.,0.25,1.,1./), & ! pR
2823		(/0.,0.,0.,0.,0.,7.5,0./), & ! hEff
2824	0	'Right column with unstable mixed layer')
2825
2826		! Left column with unstable mixed layer
2827		call find_neutral_surface_positions_continuous(3, &
2828		(/0.,10.,20.,30./), (/10.,14.,12.,4./), (/0.,0.,0.,0./), & ! Left positions, T and S
2829		(/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
2830		(/0.,10.,20.,30./), (/14.,14.,10.,2./), (/0.,0.,0.,0./), & ! Right positions, T and S
2831		(/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
2832	0	PiLRo, PiRLo, KoL, KoR, hEff)
2833		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, KoL, KoR, PiLRo, PiRLo, hEff, &
2834		(/1,1,1,2,3,3,3,3/), & ! kL
2835		(/1,2,3,3,3,3,3,3/), & ! kR
2836		(/0.,0.,0.,0.,0.,0.25,1.,1./), & ! pL
2837		(/0.,0.,0.,0.,0.,0.,.75,1./), & ! pR
2838		(/0.,0.,0.,0.,0.,7.5,0./), & ! hEff
2839	0	'Left column with unstable mixed layer')
2840
2841		! Two unstable mixed layers
2842		call find_neutral_surface_positions_continuous(3, &
2843		(/0.,10.,20.,30./), (/8.,12.,10.,2./), (/0.,0.,0.,0./), & ! Left positions, T and S
2844		(/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Left dRdT and dRdS
2845		(/0.,10.,20.,30./), (/10.,14.,12.,4./), (/0.,0.,0.,0./), & ! Right positions, T and S
2846		(/-1.,-1.,-1.,-1./), (/1.,1.,1.,1./), &! Right dRdT and dRdS
2847	0	PiLRo, PiRLo, KoL, KoR, hEff)
2848		ndiff_unit_tests_continuous = ndiff_unit_tests_continuous .or. test_nsp(v, 8, KoL, KoR, PiLRo, PiRLo, hEff, &
2849		(/1,1,1,1,2,3,3,3/), & ! kL
2850		(/1,2,3,3,3,3,3,3/), & ! kR
2851		(/0.,0.,0.,0.,0.,0.,0.75,1./), & ! pL
2852		(/0.,0.,0.,0.5,0.5,0.5,1.,1./), & ! pR
2853		(/0.,0.,0.,0.,0.,6.,0./), & ! hEff
2854	0	'Two unstable mixed layers')
2855
2856	0	if (.not. ndiff_unit_tests_continuous) write(stdout,*) 'Pass'
2857
2858	0	end function ndiff_unit_tests_continuous
2859
2860	0	logical function ndiff_unit_tests_discontinuous(verbose)
2861		logical, intent(in) :: verbose !< It true, write results to stdout
2862		! Local variables
2863		integer, parameter :: nk = 3
2864		integer, parameter :: ns = nk*4
2865		real, dimension(nk) :: Sl, Sr ! Salinities [ppt] and temperatures [degC]
2866		real, dimension(nk) :: hl, hr ! Thicknesses in pressure units [R L2 T-2 ~> Pa] or other
2867		! arbitrary units [arbitrary]
2868		real, dimension(nk,2) :: TiL, SiL, TiR, SiR ! Cell edge salinities [ppt] and temperatures [degC]
2869		real, dimension(nk,2) :: Pres_l, Pres_r ! Interface pressures [R L2 T-2 ~> Pa]
2870		integer, dimension(ns) :: KoL, KoR ! Index of the layer where the interface is found in the
2871		! left and right columns
2872		real, dimension(ns) :: PoL, PoR ! Fractional position of neutral surface within layer KoL
2873		! of the left column or KoR of the right column [nondim]
2874		real, dimension(ns-1) :: hEff ! Effective thickness between two neutral surfaces
2875		! in the same units as hl and hr [arbitrary]
2876	0	type(neutral_diffusion_CS) :: CS !< Neutral diffusion control structure
2877		real, dimension(nk,2) :: ppoly_T_l, ppoly_T_r ! Linear reconstruction for T [degC]
2878		real, dimension(nk,2) :: ppoly_S_l, ppoly_S_r ! Linear reconstruction for S [ppt]
2879		logical, dimension(nk) :: stable_l, stable_r
2880		integer :: k
2881		logical :: v
2882
2883	0	v = verbose
2884	0	ndiff_unit_tests_discontinuous = .false. ! Normally return false
2885	0	write(stdout,*) '==== MOM_neutral_diffusion: ndiff_unit_tests_discontinuous ='
2886
2887		! Unit tests for find_neutral_surface_positions_discontinuous
2888		! Salinity is 0 for all these tests
2889	0	allocate(CS%EOS)
2890	0	call EOS_manual_init(CS%EOS, form_of_EOS=EOS_LINEAR, dRho_dT=-1., dRho_dS=0.)
2891	0	Sl(:) = 0. ; Sr(:) = 0. ; ; SiL(:,:) = 0. ; SiR(:,:) = 0.
2892	0	ppoly_T_l(:,:) = 0. ; ppoly_T_r(:,:) = 0.
2893	0	ppoly_S_l(:,:) = 0. ; ppoly_S_r(:,:) = 0.
2894		! Intialize any control structures needed for unit tests
2895	0	CS%ref_pres = -1.
2896
2897	0	hL = (/10.,10.,10./) ; hR = (/10.,10.,10./)
2898	0	Pres_l(1,1) = 0. ; Pres_l(1,2) = hL(1) ; Pres_r(1,1) = 0. ; Pres_r(1,2) = hR(1)
2899	0	do k = 2,nk
2900	0	Pres_l(k,1) = Pres_l(k-1,2)
2901	0	Pres_l(k,2) = Pres_l(k,1) + hL(k)
2902	0	Pres_r(k,1) = Pres_r(k-1,2)
2903	0	Pres_r(k,2) = Pres_r(k,1) + hR(k)
2904		enddo
2905	0	CS%delta_rho_form = 'mid_pressure'
2906	0	CS%neutral_pos_method = 1
2907
2908	0	TiL(1,:) = (/ 22.00, 18.00 /) ; TiL(2,:) = (/ 18.00, 14.00 /) ; TiL(3,:) = (/ 14.00, 10.00 /)
2909	0	TiR(1,:) = (/ 22.00, 18.00 /) ; TiR(2,:) = (/ 18.00, 14.00 /) ; TiR(3,:) = (/ 14.00, 10.00 /)
2910	0	call mark_unstable_cells( CS, nk, Til, Sil, Pres_l, stable_l )
2911	0	call mark_unstable_cells( CS, nk, Tir, Sir, Pres_r, stable_r )
2912		call find_neutral_surface_positions_discontinuous(CS, nk, Pres_l, hL, TiL, SiL, ppoly_T_l, ppoly_S_l, stable_l, &
2913	0	Pres_r, hR, TiR, SiR, ppoly_T_r, ppoly_S_r, stable_r, PoL, PoR, KoL, KoR, hEff)
2914		ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, KoL, KoR, PoL, PoR, hEff, &
2915		(/ 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3 /), & ! KoL
2916		(/ 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3 /), & ! KoR
2917		(/ 0.00, 0.00, 1.00, 1.00, 0.00, 0.00, 1.00, 1.00, 0.00, 0.00, 1.00, 1.00 /), & ! PoL
2918		(/ 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 1.00, 1.00 /), & ! PoR
2919		(/ 0.00, 10.00, 0.00, 0.00, 0.00, 10.00, 0.00, 0.00, 0.00, 10.00, 0.00 /), & ! hEff
2920	0	'Identical Columns')
2921
2922	0	TiL(1,:) = (/ 22.00, 18.00 /) ; TiL(2,:) = (/ 18.00, 14.00 /) ; TiL(3,:) = (/ 14.00, 10.00 /)
2923	0	TiR(1,:) = (/ 20.00, 16.00 /) ; TiR(2,:) = (/ 16.00, 12.00 /) ; TiR(3,:) = (/ 12.00, 8.00 /)
2924	0	call mark_unstable_cells( CS, nk, Til, Sil, Pres_l, stable_l )
2925	0	call mark_unstable_cells( CS, nk, Tir, Sir, Pres_r, stable_r )
2926		call find_neutral_surface_positions_discontinuous(CS, nk, Pres_l, hL, TiL, SiL, ppoly_T_l, ppoly_S_l, stable_l, &
2927	0	Pres_r, hR, TiR, SiR, ppoly_T_r, ppoly_S_r, stable_r, PoL, PoR, KoL, KoR, hEff)
2928		ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, KoL, KoR, PoL, PoR, hEff, &
2929		(/ 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3 /), & ! KoL
2930		(/ 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3 /), & ! KoR
2931		(/ 0.00, 0.50, 1.00, 0.00, 0.50, 0.50, 1.00, 0.00, 0.50, 0.50, 1.00, 1.00 /), & ! PoL
2932		(/ 0.00, 0.00, 0.50, 0.50, 1.00, 0.00, 0.50, 0.50, 1.00, 0.00, 0.50, 1.00 /), & ! PoR
2933		(/ 0.00, 5.00, 0.00, 5.00, 0.00, 5.00, 0.00, 5.00, 0.00, 5.00, 0.00 /), & ! hEff
2934	0	'Right slightly cooler')
2935
2936	0	TiL(1,:) = (/ 20.00, 16.00 /) ; TiL(2,:) = (/ 16.00, 12.00 /) ; TiL(3,:) = (/ 12.00, 8.00 /)
2937	0	TiR(1,:) = (/ 22.00, 18.00 /) ; TiR(2,:) = (/ 18.00, 14.00 /) ; TiR(3,:) = (/ 14.00, 10.00 /)
2938	0	call mark_unstable_cells( CS, nk, Til, Sil, Pres_l, stable_l )
2939	0	call mark_unstable_cells( CS, nk, Tir, Sir, Pres_r, stable_r )
2940		call find_neutral_surface_positions_discontinuous(CS, nk, Pres_l, hL, TiL, SiL, ppoly_T_l, ppoly_S_l, stable_l, &
2941	0	Pres_r, hR, TiR, SiR, ppoly_T_r, ppoly_S_r, stable_r, PoL, PoR, KoL, KoR, hEff)
2942		ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, KoL, KoR, PoL, PoR, hEff, &
2943		(/ 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3 /), & ! KoL
2944		(/ 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3 /), & ! KoR
2945		(/ 0.00, 0.00, 0.50, 0.50, 1.00, 0.00, 0.50, 0.50, 1.00, 0.00, 0.50, 1.00 /), & ! PoL
2946		(/ 0.00, 0.50, 1.00, 0.00, 0.50, 0.50, 1.00, 0.00, 0.50, 0.50, 1.00, 1.00 /), & ! PoR
2947		(/ 0.00, 5.00, 0.00, 5.00, 0.00, 5.00, 0.00, 5.00, 0.00, 5.00, 0.00 /), & ! hEff
2948	0	'Left slightly cooler')
2949
2950	0	TiL(1,:) = (/ 22.00, 20.00 /) ; TiL(2,:) = (/ 18.00, 16.00 /) ; TiL(3,:) = (/ 14.00, 12.00 /)
2951	0	TiR(1,:) = (/ 32.00, 24.00 /) ; TiR(2,:) = (/ 22.00, 14.00 /) ; TiR(3,:) = (/ 12.00, 4.00 /)
2952	0	call mark_unstable_cells( CS, nk, Til, Sil, Pres_l, stable_l )
2953	0	call mark_unstable_cells( CS, nk, Tir, Sir, Pres_r, stable_r )
2954		call find_neutral_surface_positions_discontinuous(CS, nk, Pres_l, hL, TiL, SiL, ppoly_T_l, ppoly_S_l, stable_l, &
2955	0	Pres_r, hR, TiR, SiR, ppoly_T_r, ppoly_S_r, stable_r, PoL, PoR, KoL, KoR, hEff)
2956		ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, KoL, KoR, PoL, PoR, hEff, &
2957		(/ 1, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 3 /), & ! KoL
2958		(/ 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3 /), & ! KoR
2959		(/ 0.00, 0.00, 0.00, 0.00, 1.00, 0.00, 1.00, 0.00, 0.00, 1.00, 1.00, 1.00 /), & ! PoL
2960		(/ 0.00, 1.00, 0.00, 0.00, 0.25, 0.50, 0.75, 1.00, 0.00, 0.00, 0.00, 1.00 /), & ! PoR
2961		(/ 0.00, 0.00, 0.00, 4.00, 0.00, 4.00, 0.00, 0.00, 0.00, 0.00, 0.00 /), & ! hEff
2962	0	'Right more strongly stratified')
2963
2964	0	TiL(1,:) = (/ 22.00, 18.00 /) ; TiL(2,:) = (/ 18.00, 14.00 /) ; TiL(3,:) = (/ 14.00, 10.00 /)
2965	0	TiR(1,:) = (/ 14.00, 14.00 /) ; TiR(2,:) = (/ 14.00, 14.00 /) ; TiR(3,:) = (/ 12.00, 8.00 /)
2966	0	call mark_unstable_cells( CS, nk, Til, Sil, Pres_l, stable_l )
2967	0	call mark_unstable_cells( CS, nk, Tir, Sir, Pres_r, stable_r )
2968		call find_neutral_surface_positions_discontinuous(CS, nk, Pres_l, hL, TiL, SiL, ppoly_T_l, ppoly_S_l, stable_l, &
2969	0	Pres_r, hR, TiR, SiR, ppoly_T_r, ppoly_S_r, stable_r, PoL, PoR, KoL, KoR, hEff)
2970		ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, KoL, KoR, PoL, PoR, hEff, &
2971		(/ 1, 1, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3 /), & ! KoL
2972		(/ 1, 1, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3 /), & ! KoR
2973		(/ 0.00, 0.00, 0.00, 0.00, 0.00, 1.00, 0.00, 1.00, 0.00, 0.50, 1.00, 1.00 /), & ! PoL
2974		(/ 0.00, 1.00, 0.00, 1.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.50, 1.00 /), & ! PoR
2975		(/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 5.00, 0.00 /), & ! hEff
2976	0	'Deep Mixed layer on the right')
2977
2978	0	TiL(1,:) = (/ 14.00, 14.00 /) ; TiL(2,:) = (/ 14.00, 12.00 /) ; TiL(3,:) = (/ 10.00, 8.00 /)
2979	0	TiR(1,:) = (/ 14.00, 14.00 /) ; TiR(2,:) = (/ 14.00, 14.00 /) ; TiR(3,:) = (/ 14.00, 14.00 /)
2980	0	call mark_unstable_cells( CS, nk, Til, Sil, Pres_l, stable_l )
2981	0	call mark_unstable_cells( CS, nk, Tir, Sir, Pres_r, stable_r )
2982		call find_neutral_surface_positions_discontinuous(CS, nk, Pres_l, hL, TiL, SiL, ppoly_T_l, ppoly_S_l, stable_l, &
2983	0	Pres_r, hR, TiR, SiR, ppoly_T_r, ppoly_S_r, stable_r, PoL, PoR, KoL, KoR, hEff)
2984		ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, KoL, KoR, PoL, PoR, hEff, &
2985		(/ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3 /), & ! KoL
2986		(/ 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 3, 3 /), & ! KoR
2987		(/ 0.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00 /), & ! PoL
2988		(/ 0.00, 0.00, 0.00, 1.00, 0.00, 1.00, 0.00, 1.00, 1.00, 1.00, 1.00, 1.00 /), & ! PoR
2989		(/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 /), & ! hEff
2990	0	'Right unstratified column')
2991
2992	0	TiL(1,:) = (/ 14.00, 14.00 /) ; TiL(2,:) = (/ 14.00, 12.00 /) ; TiL(3,:) = (/ 10.00, 8.00 /)
2993	0	TiR(1,:) = (/ 14.00, 14.00 /) ; TiR(2,:) = (/ 14.00, 14.00 /) ; TiR(3,:) = (/ 12.00, 4.00 /)
2994	0	call mark_unstable_cells( CS, nk, Til, Sil, Pres_l, stable_l )
2995	0	call mark_unstable_cells( CS, nk, Tir, Sir, Pres_r, stable_r )
2996		call find_neutral_surface_positions_discontinuous(CS, nk, Pres_l, hL, TiL, SiL, ppoly_T_l, ppoly_S_l, stable_l, &
2997	0	Pres_r, hR, TiR, SiR, ppoly_T_r, ppoly_S_r, stable_r, PoL, PoR, KoL, KoR, hEff)
2998		ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, KoL, KoR, PoL, PoR, hEff, &
2999		(/ 1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3 /), & ! KoL
3000		(/ 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 3, 3 /), & ! KoR
3001		(/ 0.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.00, 1.00, 1.00, 0.00, 1.00, 1.00 /), & ! PoL
3002		(/ 0.00, 0.00, 0.00, 1.00, 0.00, 1.00, 0.00, 0.00, 0.00, 0.25, 0.50, 1.00 /), & ! PoR
3003		(/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 4.00, 0.00 /), & ! hEff
3004	0	'Right unstratified column')
3005
3006	0	TiL(1,:) = (/ 14.00, 14.00 /) ; TiL(2,:) = (/ 14.00, 10.00 /) ; TiL(3,:) = (/ 10.00, 2.00 /)
3007	0	TiR(1,:) = (/ 14.00, 14.00 /) ; TiR(2,:) = (/ 14.00, 10.00 /) ; TiR(3,:) = (/ 10.00, 2.00 /)
3008	0	call mark_unstable_cells( CS, nk, Til, Sil, Pres_l, stable_l )
3009	0	call mark_unstable_cells( CS, nk, Tir, Sir, Pres_r, stable_r )
3010		call find_neutral_surface_positions_discontinuous(CS, nk, Pres_l, hL, TiL, SiL, ppoly_T_l, ppoly_S_l, stable_l, &
3011	0	Pres_r, hR, TiR, SiR, ppoly_T_r, ppoly_S_r, stable_r, PoL, PoR, KoL, KoR, hEff)
3012		ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, KoL, KoR, PoL, PoR, hEff, &
3013		(/ 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3 /), & ! KoL
3014		(/ 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3 /), & ! KoR
3015		(/ 0.00, 1.00, 1.00, 1.00, 0.00, 0.00, 1.00, 1.00, 0.00, 0.00, 1.00, 1.00 /), & ! PoL
3016		(/ 0.00, 0.00, 0.00, 1.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 1.00, 1.00 /), & ! PoR
3017		(/ 0.00, 0.00, 0.00, 0.00, 0.00, 10.00, 0.00, 0.00, 0.00, 10.00, 0.00 /), & ! hEff
3018	0	'Identical columns with mixed layer')
3019
3020	0	TiL(1,:) = (/ 14.00, 12.00 /) ; TiL(2,:) = (/ 10.00, 10.00 /) ; TiL(3,:) = (/ 8.00, 2.00 /)
3021	0	TiR(1,:) = (/ 14.00, 12.00 /) ; TiR(2,:) = (/ 12.00, 8.00 /) ; TiR(3,:) = (/ 8.00, 2.00 /)
3022	0	call mark_unstable_cells( CS, nk, Til, Sil, Pres_l, stable_l )
3023	0	call mark_unstable_cells( CS, nk, Tir, Sir, Pres_r, stable_r )
3024		call find_neutral_surface_positions_discontinuous(CS, nk, Pres_l, hL, TiL, SiL, ppoly_T_l, ppoly_S_l, stable_l, &
3025	0	Pres_r, hR, TiR, SiR, ppoly_T_r, ppoly_S_r, stable_r, PoL, PoR, KoL, KoR, hEff)
3026		ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, KoL, KoR, PoL, PoR, hEff, &
3027		(/ 1, 1, 1, 1, 2, 2, 3, 3, 3, 3, 3, 3 /), & ! KoL
3028		(/ 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3 /), & ! KoR
3029		(/ 0.00, 0.00, 1.00, 1.00, 0.00, 1.00, 0.00, 0.00, 0.00, 0.00, 1.00, 1.00 /), & ! PoL
3030		(/ 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 0.00, 1.00, 1.00, 0.00, 1.00, 1.00 /), & ! PoR
3031		(/ 0.00, 10.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 10.00, 0.00 /), & ! hEff
3032	0	'Left interior unstratified')
3033
3034	0	TiL(1,:) = (/ 12.00, 12.00 /) ; TiL(2,:) = (/ 12.00, 10.00 /) ; TiL(3,:) = (/ 10.00, 6.00 /)
3035	0	TiR(1,:) = (/ 12.00, 10.00 /) ; TiR(2,:) = (/ 10.00, 12.00 /) ; TiR(3,:) = (/ 8.00, 4.00 /)
3036	0	call mark_unstable_cells( CS, nk, Til, Sil, Pres_l, stable_l )
3037	0	call mark_unstable_cells( CS, nk, Tir, Sir, Pres_r, stable_r )
3038		call find_neutral_surface_positions_discontinuous(CS, nk, Pres_l, hL, TiL, SiL, ppoly_T_l, ppoly_S_l, stable_l, &
3039	0	Pres_r, hR, TiR, SiR, ppoly_T_r, ppoly_S_r, stable_r, PoL, PoR, KoL, KoR, hEff)
3040		ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, KoL, KoR, PoL, PoR, hEff, &
3041		(/ 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3 /), & ! KoL
3042		(/ 1, 1, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3 /), & ! KoR
3043		(/ 0.00, 1.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 0.00, 0.50, 1.00, 1.00 /), & ! PoL
3044		(/ 0.00, 0.00, 0.00, 0.00, 1.00, 1.00, 0.00, 1.00, 0.00, 0.00, 0.50, 1.00 /), & ! PoR
3045		(/ 0.00, 0.00, 0.00, 10.00, 0.00, 0.00, 0.00, 0.00, 0.00, 5.00, 0.00 /), & ! hEff
3046	0	'Left mixed layer, Right unstable interior')
3047
3048	0	TiL(1,:) = (/ 14.00, 14.00 /) ; TiL(2,:) = (/ 10.00, 10.00 /) ; TiL(3,:) = (/ 8.00, 6.00 /)
3049	0	TiR(1,:) = (/ 10.00, 14.00 /) ; TiR(2,:) = (/ 16.00, 16.00 /) ; TiR(3,:) = (/ 12.00, 4.00 /)
3050	0	call mark_unstable_cells( CS, nk, Til, Sil, Pres_l, stable_l )
3051	0	call mark_unstable_cells( CS, nk, Tir, Sir, Pres_r, stable_r )
3052		call find_neutral_surface_positions_discontinuous(CS, nk, Pres_l, hL, TiL, SiL, ppoly_T_l, ppoly_S_l, stable_l, &
3053	0	Pres_r, hR, TiR, SiR, ppoly_T_r, ppoly_S_r, stable_r, PoL, PoR, KoL, KoR, hEff)
3054		ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, KoL, KoR, PoL, PoR, hEff, &
3055		(/ 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3 /), & ! KoL
3056		(/ 1, 1, 1, 1, 1, 1, 2, 2, 3, 3, 3, 3 /), & ! KoR
3057		(/ 0.00, 1.00, 0.00, 1.00, 1.00, 1.00, 1.00, 1.00, 0.00, 0.00, 1.00, 1.00 /), & ! PoL
3058		(/ 0.00, 0.00, 0.00, 0.00, 0.00, 1.00, 0.00, 1.00, 0.00, 0.50, 0.75, 1.00 /), & ! PoR
3059		(/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 4.00, 0.00 /), & ! hEff
3060	0	'Left thick mixed layer, Right unstable mixed')
3061
3062	0	TiL(1,:) = (/ 8.00, 12.00 /) ; TiL(2,:) = (/ 12.00, 10.00 /) ; TiL(3,:) = (/ 8.00, 4.00 /)
3063	0	TiR(1,:) = (/ 10.00, 14.00 /) ; TiR(2,:) = (/ 14.00, 12.00 /) ; TiR(3,:) = (/ 10.00, 6.00 /)
3064	0	call mark_unstable_cells( CS, nk, Til, Sil, Pres_l, stable_l )
3065	0	call mark_unstable_cells( CS, nk, Tir, Sir, Pres_r, stable_r )
3066		call find_neutral_surface_positions_discontinuous(CS, nk, Pres_l, hL, TiL, SiL, ppoly_T_l, ppoly_S_l, stable_l, &
3067	0	Pres_r, hR, TiR, SiR, ppoly_T_r, ppoly_S_r, stable_r, PoL, PoR, KoL, KoR, hEff)
3068		ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. test_nsp(v, 12, KoL, KoR, PoL, PoR, hEff, &
3069		(/ 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3 /), & ! KoL
3070		(/ 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3 /), & ! KoR
3071		(/ 0.00, 1.00, 1.00, 1.00, 0.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.50, 1.00 /), & ! PoL
3072		(/ 0.00, 0.00, 0.00, 1.00, 0.00, 1.00, 1.00, 0.00, 0.00, 0.50, 1.00, 1.00 /), & ! PoR
3073		(/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 5.00, 0.00 /), & ! hEff
3074	0	'Unstable mixed layers, left cooler')
3075
3076	0	call EOS_manual_init(CS%EOS, form_of_EOS = EOS_LINEAR, dRho_dT = -1., dRho_dS = 2.)
3077		! Tests for linearized version of searching the layer for neutral surface position
3078		! EOS linear in T, uniform alpha
3079	0	CS%max_iter = 10
3080		! Unit tests require explicit initialization of tolerance
3081	0	CS%Drho_tol = 0.
3082	0	CS%x_tol = 0.
3083		ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5, &
3084		find_neutral_pos_linear(CS, 0., 10., 35., -0.2, 0., &
3085		-0.2, 0., -0.2, 0., &
3086	0	(/12.,-4./), (/34.,0./)), "Temp Uniform Linearized Alpha/Beta"))
3087		! EOS linear in S, uniform beta
3088		ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5, &
3089		find_neutral_pos_linear(CS, 0., 10., 35., 0., 0.8, &
3090		0., 0.8, 0., 0.8, &
3091	0	(/12.,0./), (/34.,2./)), "Salt Uniform Linearized Alpha/Beta"))
3092		! EOS linear in T/S, uniform alpha/beta
3093		ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5, &
3094		find_neutral_pos_linear(CS, 0., 10., 35., -0.5, 0.5, &
3095		-0.5, 0.5, -0.5, 0.5, &
3096	0	(/12.,-4./), (/34.,2./)), "Temp/salt Uniform Linearized Alpha/Beta"))
3097		! EOS linear in T, insensitive to So
3098		ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5, &
3099		find_neutral_pos_linear(CS, 0., 10., 35., -0.2, 0., &
3100		-0.4, 0., -0.6, 0., &
3101	0	(/12.,-4./), (/34.,0./)), "Temp stratified Linearized Alpha/Beta"))
3102		! EOS linear in S, insensitive to T
3103		ndiff_unit_tests_discontinuous = ndiff_unit_tests_discontinuous .or. (test_rnp(0.5, &
3104		find_neutral_pos_linear(CS, 0., 10., 35., 0., 0.8, &
3105		0., 1.0, 0., 0.5, &
3106	0	(/12.,0./), (/34.,2./)), "Salt stratified Linearized Alpha/Beta"))
3107	0	if (.not. ndiff_unit_tests_discontinuous) write(stdout,*) 'Pass'
3108
3109	0	end function ndiff_unit_tests_discontinuous
3110
3111		!> Returns true if a test of fv_diff() fails, and conditionally writes results to stream
3112	0	logical function test_fv_diff(verbose, hkm1, hk, hkp1, Skm1, Sk, Skp1, Ptrue, title)
3113		logical, intent(in) :: verbose !< If true, write results to stdout
3114		real, intent(in) :: hkm1 !< Left cell width [nondim]
3115		real, intent(in) :: hk !< Center cell width [nondim]
3116		real, intent(in) :: hkp1 !< Right cell width [nondim]
3117		real, intent(in) :: Skm1 !< Left cell average value in arbitrary units [arbitrary]
3118		real, intent(in) :: Sk !< Center cell average value in arbitrary units [arbitrary]
3119		real, intent(in) :: Skp1 !< Right cell average value in arbitrary units [arbitrary]
3120		real, intent(in) :: Ptrue !< True answer in arbitrary units [arbitrary]
3121		character(len=*), intent(in) :: title !< Title for messages
3122
3123		! Local variables
3124		integer :: stdunit
3125		real :: Pret ! Returned normalized gradient in arbitrary units [arbitrary]
3126
3127	0	Pret = fv_diff(hkm1, hk, hkp1, Skm1, Sk, Skp1)
3128	0	test_fv_diff = (Pret /= Ptrue)
3129
3130	0	if (test_fv_diff .or. verbose) then
3131	0	stdunit = stdout
3132	0	if (test_fv_diff) stdunit = stderr ! In case of wrong results, write to error stream
3133	0	write(stdunit,'(a)') title
3134	0	if (test_fv_diff) then
3135	0	write(stdunit,'(2(1x,a,f20.16),1x,a)') 'pRet=',Pret,'pTrue=',Ptrue,'WRONG!'
3136		else
3137	0	write(stdunit,'(2(1x,a,f20.16))') 'pRet=',Pret,'pTrue=',Ptrue
3138		endif
3139		endif
3140
3141	0	end function test_fv_diff
3142
3143		!> Returns true if a test of fvlsq_slope() fails, and conditionally writes results to stream
3144	0	logical function test_fvlsq_slope(verbose, hkm1, hk, hkp1, Skm1, Sk, Skp1, Ptrue, title)
3145		logical, intent(in) :: verbose !< If true, write results to stdout
3146		real, intent(in) :: hkm1 !< Left cell width in arbitrary units [B ~> b]
3147		real, intent(in) :: hk !< Center cell width in arbitrary units [B ~> b]
3148		real, intent(in) :: hkp1 !< Right cell width in arbitrary units [B ~> b]
3149		real, intent(in) :: Skm1 !< Left cell average value in arbitrary units [A ~> a]
3150		real, intent(in) :: Sk !< Center cell average value in arbitrary units [A ~> a]
3151		real, intent(in) :: Skp1 !< Right cell average value in arbitrary units [A ~> a]
3152		real, intent(in) :: Ptrue !< True answer in arbitrary units [A B-1 ~> a b-1]
3153		character(len=*), intent(in) :: title !< Title for messages
3154
3155		! Local variables
3156		integer :: stdunit
3157		real :: Pret ! Returned slope value [A B-1 ~> a b-1]
3158
3159	0	Pret = fvlsq_slope(hkm1, hk, hkp1, Skm1, Sk, Skp1)
3160	0	test_fvlsq_slope = (Pret /= Ptrue)
3161
3162	0	if (test_fvlsq_slope .or. verbose) then
3163	0	stdunit = stdout
3164	0	if (test_fvlsq_slope) stdunit = stderr ! In case of wrong results, write to error stream
3165	0	write(stdunit,'(a)') title
3166	0	if (test_fvlsq_slope) then
3167	0	write(stdunit,'(2(1x,a,f20.16),1x,a)') 'pRet=',Pret,'pTrue=',Ptrue,'WRONG!'
3168		else
3169	0	write(stdunit,'(2(1x,a,f20.16))') 'pRet=',Pret,'pTrue=',Ptrue
3170		endif
3171		endif
3172
3173	0	end function test_fvlsq_slope
3174
3175		!> Returns true if a test of interpolate_for_nondim_position() fails, and conditionally writes results to stream
3176	0	logical function test_ifndp(verbose, rhoNeg, Pneg, rhoPos, Ppos, Ptrue, title)
3177		logical, intent(in) :: verbose !< If true, write results to stdout
3178		real, intent(in) :: rhoNeg !< Lighter density [R ~> kg m-3]
3179		real, intent(in) :: Pneg !< Interface position of lighter density [nondim]
3180		real, intent(in) :: rhoPos !< Heavier density [R ~> kg m-3]
3181		real, intent(in) :: Ppos !< Interface position of heavier density [nondim]
3182		real, intent(in) :: Ptrue !< True answer [nondim]
3183		character(len=*), intent(in) :: title !< Title for messages
3184
3185		! Local variables
3186		integer :: stdunit
3187		real :: Pret ! Interpolated fractional position [nondim]
3188
3189	0	Pret = interpolate_for_nondim_position(rhoNeg, Pneg, rhoPos, Ppos)
3190	0	test_ifndp = (Pret /= Ptrue)
3191
3192	0	if (test_ifndp .or. verbose) then
3193	0	stdunit = stdout
3194	0	if (test_ifndp) stdunit = stderr ! In case of wrong results, write to error stream
3195	0	write(stdunit,'(a)') title
3196	0	if (test_ifndp) then
3197		write(stdunit,'(4(1x,a,f20.16),2(1x,a,1pe22.15),1x,a)') &
3198	0	'r1=',rhoNeg,'p1=',Pneg,'r2=',rhoPos,'p2=',Ppos,'pRet=',Pret,'pTrue=',Ptrue,'WRONG!'
3199		else
3200		write(stdunit,'(4(1x,a,f20.16),2(1x,a,1pe22.15))') &
3201	0	'r1=',rhoNeg,'p1=',Pneg,'r2=',rhoPos,'p2=',Ppos,'pRet=',Pret,'pTrue=',Ptrue
3202		endif
3203		endif
3204
3205	0	end function test_ifndp
3206
3207		!> Returns true if comparison of Po and Ptrue fails, and conditionally writes results to stream
3208	0	logical function test_data1d(verbose, nk, Po, Ptrue, title)
3209		logical, intent(in) :: verbose !< If true, write results to stdout
3210		integer, intent(in) :: nk !< Number of layers
3211		real, dimension(nk), intent(in) :: Po !< Calculated answer [arbitrary]
3212		real, dimension(nk), intent(in) :: Ptrue !< True answer [arbitrary]
3213		character(len=*), intent(in) :: title !< Title for messages
3214
3215		! Local variables
3216		integer :: k, stdunit
3217
3218	0	test_data1d = .false.
3219	0	do k = 1,nk
3220	0	if (Po(k) /= Ptrue(k)) test_data1d = .true.
3221		enddo
3222
3223	0	if (test_data1d .or. verbose) then
3224	0	stdunit = stdout
3225	0	if (test_data1d) stdunit = stderr ! In case of wrong results, write to error stream
3226	0	write(stdunit,'(a)') title
3227	0	do k = 1,nk
3228	0	if (Po(k) /= Ptrue(k)) then
3229	0	test_data1d = .true.
3230		write(stdunit,'(a,I0,2(1x,a,f20.16),1x,a,1pe22.15,1x,a)') &
3231	0	'k=',k,'Po=',Po(k),'Ptrue=',Ptrue(k),'err=',Po(k)-Ptrue(k),'WRONG!'
3232		else
3233	0	if (verbose) &
3234		write(stdunit,'(a,I0,2(1x,a,f20.16),1x,a,1pe22.15)') &
3235	0	'k=',k,'Po=',Po(k),'Ptrue=',Ptrue(k),'err=',Po(k)-Ptrue(k)
3236		endif
3237		enddo
3238		endif
3239
3240	0	end function test_data1d
3241
3242		!> Returns true if comparison of Po and Ptrue fails, and conditionally writes results to stream
3243	0	logical function test_data1di(verbose, nk, Po, Ptrue, title)
3244		logical, intent(in) :: verbose !< If true, write results to stdout
3245		integer, intent(in) :: nk !< Number of layers
3246		integer, dimension(nk), intent(in) :: Po !< Calculated answer [arbitrary]
3247		integer, dimension(nk), intent(in) :: Ptrue !< True answer [arbitrary]
3248		character(len=*), intent(in) :: title !< Title for messages
3249
3250		! Local variables
3251		integer :: k, stdunit
3252
3253	0	test_data1di = .false.
3254	0	do k = 1,nk
3255	0	if (Po(k) /= Ptrue(k)) test_data1di = .true.
3256		enddo
3257
3258	0	if (test_data1di .or. verbose) then
3259	0	stdunit = stdout
3260	0	if (test_data1di) stdunit = stderr ! In case of wrong results, write to error stream
3261	0	write(stdunit,'(a)') title
3262	0	do k = 1,nk
3263	0	if (Po(k) /= Ptrue(k)) then
3264	0	test_data1di = .true.
3265	0	write(stdunit,'(a,I0,2(1x,a,i5),1x,a)') 'k=',k,'Io=',Po(k),'Itrue=',Ptrue(k),'WRONG!'
3266		else
3267	0	if (verbose) &
3268	0	write(stdunit,'(a,I0,2(1x,a,i5))') 'k=',k,'Io=',Po(k),'Itrue=',Ptrue(k)
3269		endif
3270		enddo
3271		endif
3272
3273	0	end function test_data1di
3274
3275		!> Returns true if output of find_neutral_surface_positions() does not match correct values,
3276		!! and conditionally writes results to stream
3277	0	logical function test_nsp(verbose, ns, KoL, KoR, pL, pR, hEff, KoL0, KoR0, pL0, pR0, hEff0, title)
3278		logical, intent(in) :: verbose !< If true, write results to stdout
3279		integer, intent(in) :: ns !< Number of surfaces
3280		integer, dimension(ns), intent(in) :: KoL !< Index of first left interface above neutral surface
3281		integer, dimension(ns), intent(in) :: KoR !< Index of first right interface above neutral surface
3282		real, dimension(ns), intent(in) :: pL !< Fractional position of neutral surface within layer
3283		!! KoL of left column [nondim]
3284		real, dimension(ns), intent(in) :: pR !< Fractional position of neutral surface within layer
3285		!! KoR of right column [nondim]
3286		real, dimension(ns-1), intent(in) :: hEff !< Effective thickness between two neutral surfaces [R L2 T-2 ~> Pa]
3287		integer, dimension(ns), intent(in) :: KoL0 !< Correct value for KoL
3288		integer, dimension(ns), intent(in) :: KoR0 !< Correct value for KoR
3289		real, dimension(ns), intent(in) :: pL0 !< Correct value for pL [nondim]
3290		real, dimension(ns), intent(in) :: pR0 !< Correct value for pR [nondim]
3291		real, dimension(ns-1), intent(in) :: hEff0 !< Correct value for hEff [R L2 T-2 ~> Pa]
3292		character(len=*), intent(in) :: title !< Title for messages
3293
3294		! Local variables
3295		integer :: k, stdunit
3296		logical :: this_row_failed
3297
3298	0	test_nsp = .false.
3299	0	do k = 1,ns
3300	0	test_nsp = test_nsp .or. compare_nsp_row(KoL(k), KoR(k), pL(k), pR(k), KoL0(k), KoR0(k), pL0(k), pR0(k))
3301	0	if (k < ns) then
3302	0	if (hEff(k) /= hEff0(k)) test_nsp = .true.
3303		endif
3304		enddo
3305
3306	0	if (test_nsp .or. verbose) then
3307	0	stdunit = stdout
3308	0	if (test_nsp) stdunit = stderr ! In case of wrong results, write to error stream
3309	0	write(stdunit,'(a)') title
3310	0	do k = 1,ns
3311	0	this_row_failed = compare_nsp_row(KoL(k), KoR(k), pL(k), pR(k), KoL0(k), KoR0(k), pL0(k), pR0(k))
3312	0	if (this_row_failed) then
3313	0	write(stdunit,10) k,KoL(k),pL(k),KoR(k),pR(k),' <-- WRONG!'
3314	0	write(stdunit,10) k,KoL0(k),pL0(k),KoR0(k),pR0(k),' <-- should be this'
3315		else
3316	0	write(stdunit,10) k,KoL(k),pL(k),KoR(k),pR(k)
3317		endif
3318	0	if (k < ns) then
3319	0	if (hEff(k) /= hEff0(k)) then
3320	0	write(stdunit,'(i3,8x,"layer hEff =",2(f20.16,a))') k,hEff(k)," .neq. ",hEff0(k),' <-- WRONG!'
3321		else
3322	0	write(stdunit,'(i3,8x,"layer hEff =",f20.16)') k,hEff(k)
3323		endif
3324		endif
3325		enddo
3326		endif
3327	0	if (test_nsp) call MOM_error(FATAL,"test_nsp failed")
3328
3329		10 format("ks=",i3," kL=",i3," pL=",f20.16," kR=",i3," pR=",f20.16,a)
3330	0	end function test_nsp
3331
3332		!> Compares a single row, k, of output from find_neutral_surface_positions()
3333	0	logical function compare_nsp_row(KoL, KoR, pL, pR, KoL0, KoR0, pL0, pR0)
3334		integer, intent(in) :: KoL !< Index of first left interface above neutral surface
3335		integer, intent(in) :: KoR !< Index of first right interface above neutral surface
3336		real, intent(in) :: pL !< Fractional position of neutral surface within layer KoL of left column [nondim]
3337		real, intent(in) :: pR !< Fractional position of neutral surface within layer KoR of right column [nondim]
3338		integer, intent(in) :: KoL0 !< Correct value for KoL
3339		integer, intent(in) :: KoR0 !< Correct value for KoR
3340		real, intent(in) :: pL0 !< Correct value for pL [nondim]
3341		real, intent(in) :: pR0 !< Correct value for pR [nondim]
3342
3343	0	compare_nsp_row = .false.
3344	0	if (KoL /= KoL0) compare_nsp_row = .true.
3345	0	if (KoR /= KoR0) compare_nsp_row = .true.
3346	0	if (pL /= pL0) compare_nsp_row = .true.
3347	0	if (pR /= pR0) compare_nsp_row = .true.
3348	0	end function compare_nsp_row
3349
3350		!> Compares output position from refine_nondim_position with an expected value
3351	0	logical function test_rnp(expected_pos, test_pos, title)
3352		real, intent(in) :: expected_pos !< The expected position [arbitrary]
3353		real, intent(in) :: test_pos !< The position returned by the code [arbitrary]
3354		character(len=*), intent(in) :: title !< A label for this test
3355		! Local variables
3356		integer :: stdunit
3357
3358	0	stdunit = stdout ! Output to standard error
3359	0	test_rnp = ABS(expected_pos - test_pos) > 2*EPSILON(test_pos)
3360	0	if (test_rnp) then
3361	0	write(stdunit,'(A, f20.16, " .neq. ", f20.16, " <-- WRONG")') title, expected_pos, test_pos
3362		else
3363	0	write(stdunit,'(A, f20.16, " == ", f20.16)') title, expected_pos, test_pos
3364		endif
3365	0	end function test_rnp
3366
3367		!> Deallocates neutral_diffusion control structure
3368	1	subroutine neutral_diffusion_end(CS)
3369		type(neutral_diffusion_CS), pointer :: CS !< Neutral diffusion control structure
3370
3371	1	if (associated(CS)) deallocate(CS)
3372
3373	1	end subroutine neutral_diffusion_end
3374
3375	0	end module MOM_neutral_diffusion