MOM_neutral_diffusion.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> A column-wise toolbox for implementing neutral diffusion
7
8use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
9use mom_cpu_clock, only : clock_module, clock_routine
10use mom_domains, only : pass_var
11use mom_diag_mediator, only : diag_ctrl, time_type
12use mom_diag_mediator, only : post_data, register_diag_field
14use mom_eos, only : calculate_density, calculate_density_derivs
15use mom_eos, only : eos_linear
16use mom_error_handler, only : mom_error, fatal, warning, mom_mesg, is_root_pe
17use mom_file_parser, only : get_param, log_version, param_file_type
19use mom_grid, only : ocean_grid_type
22use mom_remapping, only : average_value_ppoly, remappingschemesdoc, remappingdefaultscheme
23use mom_tracer_registry, only : tracer_registry_type, tracer_type
33use mom_io, only : stdout, stderr
35
36implicit none ; private
37
38#include <MOM_memory.h>
39
43
44!> The control structure for the MOM_neutral_diffusion module
45type, 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
132
133! This include declares and sets the variable "version".
134#include "version_variable.h"
135character(len=40) :: mdl = "MOM_neutral_diffusion" !< module name
136
137contains
138
139!> Read parameters and allocate control structure for neutral_diffusion module.
140logical 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 if (associated(cs)) then
161 call mom_error(fatal, "neutral_diffusion_init called with associated control structure.")
162 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 default=.false., do_not_log=.true.)
168 call log_version(param_file, mdl, version, &
169 "This module implements neutral diffusion of tracers", &
170 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 default=.false.)
174
175 if (.not.neutral_diffusion_init) return
176
177 allocate(cs)
178 cs%diag => diag
179 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 "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 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 "boundary layers.", default=.false.)
197 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 "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 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 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 "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 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 default=default_answer_date)
225
226 ! Initialize and configure remapping
227 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 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 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 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 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 "We recommend setting this option to false.", default=om4_remap_via_sub_cells)
249 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 answer_date=cs%remap_answer_date )
253 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 " for neutral position", default=3)
261 if (cs%neutral_pos_method > 4 .or. cs%neutral_pos_method < 0) then
262 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 default="mid_pressure")
271 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 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 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 default=10)
284 endif
285 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 "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 default=.true.)
292 endif
293
294 if (cs%interior_only) then
295 allocate(cs%hbl(szi_(g),szj_(g)), source=0.)
296 call extract_diabatic_member(diabatic_csp, kpp_csp=cs%KPP_CSp)
297 call extract_diabatic_member(diabatic_csp, energetic_pbl_csp=cs%energetic_PBL_CSp)
298 if ( .not. ASSOCIATED(cs%energetic_PBL_CSp) .and. .not. ASSOCIATED(cs%KPP_CSp) ) then
299 call mom_error(fatal,"NDIFF_INTERIOR_ONLY is true, but no valid boundary layer scheme was found")
300 endif
301
302 if (cs%tapering) then
303 allocate(cs%coeff_l(szk_(gv)+1), source=1.)
304 allocate(cs%coeff_r(szk_(gv)+1), source=1.)
305 endif
306 endif
307
308 cs%KhTh_use_vert_struct = khth_use_ebt_struct .or. khth_use_sqg_struct
309 if (cs%KhTh_use_vert_struct) &
310 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 cs%R_to_kg_m3 = us%R_to_kg_m3
314
315! call closeParameterBlock(param_file)
316 if (cs%continuous_reconstruction) then
317 cs%nsurf = 2*gv%ke+2 ! Continuous reconstruction means that every interface has two connections
318 allocate(cs%dRdT(szi_(g),szj_(g),szk_(gv)+1), source=0.)
319 allocate(cs%dRdS(szi_(g),szj_(g),szk_(gv)+1), source=0.)
320 else
321 cs%nsurf = 4*gv%ke ! Discontinuous means that every interface has four connections
322 allocate(cs%T_i(szi_(g),szj_(g),szk_(gv),2), source=0.)
323 allocate(cs%S_i(szi_(g),szj_(g),szk_(gv),2), source=0.)
324 allocate(cs%P_i(szi_(g),szj_(g),szk_(gv),2), source=0.)
325 allocate(cs%dRdT_i(szi_(g),szj_(g),szk_(gv),2), source=0.)
326 allocate(cs%dRdS_i(szi_(g),szj_(g),szk_(gv),2), source=0.)
327 allocate(cs%ppoly_coeffs_T(szi_(g),szj_(g),szk_(gv),cs%deg+1), source=0.)
328 allocate(cs%ppoly_coeffs_S(szi_(g),szj_(g),szk_(gv),cs%deg+1), source=0.)
329 allocate(cs%ns(szi_(g),szj_(g)), source=0)
330 endif
331 ! T-points
332 allocate(cs%Tint(szi_(g),szj_(g),szk_(gv)+1), source=0.)
333 allocate(cs%Sint(szi_(g),szj_(g),szk_(gv)+1), source=0.)
334 allocate(cs%Pint(szi_(g),szj_(g),szk_(gv)+1), source=0.)
335 allocate(cs%stable_cell(szi_(g),szj_(g),szk_(gv)), source=.true.)
336 ! U-points
337 allocate(cs%uPoL(g%isd:g%ied,g%jsd:g%jed, cs%nsurf), source=0.)
338 allocate(cs%uPoR(g%isd:g%ied,g%jsd:g%jed, cs%nsurf), source=0.)
339 allocate(cs%uKoL(g%isd:g%ied,g%jsd:g%jed, cs%nsurf), source=0)
340 allocate(cs%uKoR(g%isd:g%ied,g%jsd:g%jed, cs%nsurf), source=0)
341 allocate(cs%uHeff(g%isd:g%ied,g%jsd:g%jed,cs%nsurf-1), source=0.)
342 ! V-points
343 allocate(cs%vPoL(g%isd:g%ied,g%jsd:g%jed, cs%nsurf), source=0.)
344 allocate(cs%vPoR(g%isd:g%ied,g%jsd:g%jed, cs%nsurf), source=0.)
345 allocate(cs%vKoL(g%isd:g%ied,g%jsd:g%jed, cs%nsurf), source=0)
346 allocate(cs%vKoR(g%isd:g%ied,g%jsd:g%jed, cs%nsurf), source=0)
347 allocate(cs%vHeff(g%isd:g%ied,g%jsd:g%jed,cs%nsurf-1), source=0.)
348
349end function neutral_diffusion_init
350
351!> Calculate remapping factors for u/v columns used to map adjoining columns to
352!! a shared coordinate space.
353subroutine 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 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 real, dimension(SZI_(G), SZJ_(G)) :: heff_sum ! Summed effective face thicknesses [H ~> m or kg m-2]
373 integer :: imethod
374 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 integer, dimension(SZI_(G), SZJ_(G)) :: k_top ! Index of the first layer within the boundary
377 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 integer, dimension(SZI_(G), SZJ_(G)) :: k_bot ! Index of the last layer within the boundary
380 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 pa_to_h = 1. / (gv%H_to_RZ * gv%g_Earth)
385
386 k_top(:,:) = 1 ; k_bot(:,:) = 1
387 zeta_top(:,:) = 0. ; zeta_bot(:,:) = 0.
388
389 ! Check if hbl needs to be extracted
390 if (cs%interior_only) then
391 if (associated(visc%h_ML)) then
392 cs%hbl(:,:) = visc%h_ML(:,:)
393 else
394 call mom_error(fatal, "hor_bnd_diffusion requires that visc%h_ML is associated.")
395 endif
396 call pass_var(cs%hbl, g%Domain, halo=1)
397
398 ! get k-indices and zeta
399 do j=g%jsc-1, g%jec+1 ; do i=g%isc-1,g%iec+1
400 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 zeta_bot(i,j))
403 endif
404 enddo ; enddo
405 ! TODO: add similar code for BOTTOM boundary layer
406 endif
407
408 h_neglect = gv%H_subroundoff ; h_neglect_edge = gv%H_subroundoff
409
410 if (.not. cs%continuous_reconstruction) then
411 if (cs%remap_answer_date < 20190101) then
412 if (gv%Boussinesq) then
413 h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
414 else
415 h_neglect = gv%kg_m2_to_H*1.0e-30 ; h_neglect_edge = gv%kg_m2_to_H*1.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 if (cs%ref_pres>=0.) then
422 ref_pres(:) = cs%ref_pres
423 endif
424
425 if (cs%continuous_reconstruction) then
426 cs%dRdT(:,:,:) = 0.
427 cs%dRdS(:,:,:) = 0.
428 else
429 cs%T_i(:,:,:,:) = 0.
430 cs%S_i(:,:,:,:) = 0.
431 cs%dRdT_i(:,:,:,:) = 0.
432 cs%dRdS_i(:,:,:,:) = 0.
433 cs%ns(:,:) = 0.
434 cs%stable_cell(:,:,:) = .true.
435 endif
436
437 ! Calculate pressure at interfaces and layer averaged alpha/beta
438 if (present(p_surf)) then
439 do j=g%jsc-1,g%jec+1 ; do i=g%isc-1,g%iec+1
440 cs%Pint(i,j,1) = p_surf(i,j)
441 enddo ; enddo
442 else
443 cs%Pint(:,:,1) = 0.
444 endif
445 do k=1,gv%ke ; do j=g%jsc-1,g%jec+1 ; do i=g%isc-1,g%iec+1
446 cs%Pint(i,j,k+1) = cs%Pint(i,j,k) + h(i,j,k)*(gv%g_Earth*gv%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 if (.not. cs%continuous_reconstruction) then
452 if (present(p_surf)) then
453 do j=g%jsc-1,g%jec+1 ; do i=g%isc-1,g%iec+1
454 cs%P_i(i,j,1,1) = p_surf(i,j)
455 cs%P_i(i,j,1,2) = p_surf(i,j) + h(i,j,1)*(gv%H_to_RZ*gv%g_Earth)
456 enddo ; enddo
457 else
458 do j=g%jsc-1,g%jec+1 ; do i=g%isc-1,g%iec+1
459 cs%P_i(i,j,1,1) = 0.
460 cs%P_i(i,j,1,2) = h(i,j,1)*(gv%H_to_RZ*gv%g_Earth)
461 enddo ; enddo
462 endif
463 do k=2,gv%ke ; do j=g%jsc-1,g%jec+1 ; do i=g%isc-1,g%iec+1
464 cs%P_i(i,j,k,1) = cs%P_i(i,j,k-1,2)
465 cs%P_i(i,j,k,2) = cs%P_i(i,j,k-1,2) + h(i,j,k)*(gv%H_to_RZ*gv%g_Earth)
466 enddo ; enddo ; enddo
467 endif
468
469 eosdom(:) = eos_domain(g%HI, halo=1)
470 do j = g%jsc-1, g%jec+1
471 ! Interpolate state to interface
472 do i = g%isc-1, g%iec+1
473 if (cs%continuous_reconstruction) then
474 call interface_scalar(gv%ke, h(i,j,:), t(i,j,:), cs%Tint(i,j,:), 2, h_neglect)
475 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 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 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 do k=1,gv%ke
484 cs%T_i(i,j,k,1) = evaluation_polynomial( cs%ppoly_coeffs_T(i,j,k,:), cs%deg+1, 0. )
485 cs%T_i(i,j,k,2) = evaluation_polynomial( cs%ppoly_coeffs_T(i,j,k,:), cs%deg+1, 1. )
486 cs%S_i(i,j,k,1) = evaluation_polynomial( cs%ppoly_coeffs_S(i,j,k,:), cs%deg+1, 0. )
487 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 if (cs%continuous_reconstruction) then
494 do k = 1, gv%ke+1
495 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 cs%dRdS(:,j,k), cs%EOS, eosdom)
498 enddo
499 else ! Discontinuous reconstruction
500 do k = 1, gv%ke
501 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 cs%dRdS_i(:,j,k,1), cs%EOS, eosdom)
505 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 cs%dRdS_i(:,j,k,2), cs%EOS, eosdom)
509 enddo
510 endif
511 enddo
512
513 if (.not. cs%continuous_reconstruction) then
514 do j = g%jsc-1, g%jec+1 ; do i = g%isc-1, g%iec+1
515 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 if (cs%interior_only) then
517 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 do k = 1, k_bot(i,j)
520 cs%stable_cell(i,j,k) = .false.
521 enddo
522 endif
523 enddo ; enddo
524 endif
525
526 cs%uhEff(:,:,:) = 0.
527 cs%vhEff(:,:,:) = 0.
528 cs%uPoL(:,:,:) = 0.
529 cs%vPoL(:,:,:) = 0.
530 cs%uPoR(:,:,:) = 0.
531 cs%vPoR(:,:,:) = 0.
532 cs%uKoL(:,:,:) = 1
533 cs%vKoL(:,:,:) = 1
534 cs%uKoR(:,:,:) = 1
535 cs%vKoR(:,:,:) = 1
536
537 ! Neutral surface factors at U points
538 do j = g%jsc, g%jec ; do i = g%isc-1, g%iec
539 if (g%mask2dCu(i,j) > 0.0) then
540 if (cs%continuous_reconstruction) then
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 k_bot(i,j), k_bot(i+1,j), zeta_bot(i,j), zeta_bot(i+1,j))
546 else
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 hard_fail_heff = cs%hard_fail_heff)
554 endif
555 endif
556 enddo ; enddo
557
558 ! Neutral surface factors at V points
559 do j = g%jsc-1, g%jec ; do i = g%isc, g%iec
560 if (g%mask2dCv(i,j) > 0.0) then
561 if (cs%continuous_reconstruction) then
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 k_bot(i,j), k_bot(i,j+1), zeta_bot(i,j), zeta_bot(i,j+1))
567 else
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 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 if (cs%continuous_reconstruction) then
584 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 do k = 1, cs%nsurf-1 ; do j = g%jsc, g%jec ; do i = g%isc-1, g%iec
588 cs%uhEff(i,j,k) = cs%uhEff(i,j,k) / gv%H_to_pa
589 enddo ; enddo ; enddo
590 do k = 1, cs%nsurf-1 ; do j = g%jsc-1, g%jec ; do i = g%isc, g%iec
591 cs%vhEff(i,j,k) = cs%vhEff(i,j,k) / gv%H_to_pa
592 enddo ; enddo ; enddo
593 else
594 do k = 1, cs%nsurf-1 ; do j = g%jsc, g%jec ; do i = g%isc-1, g%iec
595 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 do k = 1, cs%nsurf-1 ; do j = g%jsc-1, g%jec ; do i = g%isc, g%iec
598 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 if (cs%id_uhEff_2d>0) then
604 heff_sum(:,:) = 0.
605 do k = 1,cs%nsurf-1 ; do j=g%jsc,g%jec ; do i=g%isc-1,g%iec
606 heff_sum(i,j) = heff_sum(i,j) + cs%uhEff(i,j,k)
607 enddo ; enddo ; enddo
608 call post_data(cs%id_uhEff_2d, heff_sum, cs%diag)
609 endif
610 if (cs%id_vhEff_2d>0) then
611 heff_sum(:,:) = 0.
612 do k = 1,cs%nsurf-1 ; do j=g%jsc-1,g%jec ; do i=g%isc,g%iec
613 heff_sum(i,j) = heff_sum(i,j) + cs%vhEff(i,j,k)
614 enddo ; enddo ; enddo
615 call post_data(cs%id_vhEff_2d, heff_sum, cs%diag)
616 endif
617
619
620!> Update tracer concentration due to neutral diffusion; layer thickness unchanged by this update.
621subroutine 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 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 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 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 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 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 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 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 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 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 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 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 h_neglect = gv%H_subroundoff ; h_neglect_edge = gv%H_subroundoff
679
680 if (.not. cs%continuous_reconstruction) then
681 if (cs%remap_answer_date < 20190101) then
682 h_neglect = gv%m_to_H*1.0e-30 ; h_neglect_edge = gv%m_to_H*1.0e-10
683 endif
684 endif
685
686 if (cs%KhTh_use_vert_struct) then
687 ! Compute Coef at h points
688 cs%Coef_h(:,:,:) = 0.
689 do j = g%jsc,g%jec ; do i = g%isc,g%iec
690 if (g%mask2dT(i,j)>0.) then
691 normalize = 1.0 / ((g%mask2dCu(i-1,j)+g%mask2dCu(i,j)) + &
692 (g%mask2dCv(i,j-1)+g%mask2dCv(i,j)) + 1.0e-37)
693 do k = 1, gv%ke+1
694 cs%Coef_h(i,j,k) = normalize*g%mask2dT(i,j)*((coef_x(i-1,j,k)+coef_x(i,j,k)) + &
695 (coef_y(i,j-1,k)+coef_y(i,j,k)))
696 enddo
697 endif
698 enddo ; enddo
699 call pass_var(cs%Coef_h,g%Domain)
700 endif
701
702 nk = gv%ke
703
704 do m = 1,reg%ntr ! Loop over tracer registry
705
706 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 tracer%id_dfx_2d > 0 .or. tracer%id_dfy_2d > 0) then
711 idt = 1.0 / dt
712 tendency(:,:,:) = 0.0
713 endif
714
715 uflx(:,:,:) = 0.
716 vflx(:,:,:) = 0.
717
718 ! x-flux
719 if (cs%KhTh_use_vert_struct) then
720 if (cs%tapering) then
721 do j = g%jsc,g%jec ; do i = g%isc-1,g%iec
722 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 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 cs%coeff_r(:)*cs%Coef_h(i+1,j,:))
734 endif
735 enddo ; enddo
736 else
737 do j = g%jsc,g%jec ; do i = g%isc-1,g%iec
738 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 cs%Coef_h(i+1,j,:))
747 endif
748 enddo ; enddo
749 endif
750 else
751 if (cs%tapering) then
752 do j = g%jsc,g%jec ; do i = g%isc-1,g%iec
753 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 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 cs%coeff_r(:))
765 endif
766 enddo ; enddo
767 else
768 do j = g%jsc,g%jec ; do i = g%isc-1,g%iec
769 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 cs%remap_CS, h_neglect_edge)
777 endif
778 enddo ; enddo
779 endif
780 endif
781
782 ! y-flux
783 if (cs%KhTh_use_vert_struct) then
784 if (cs%tapering) then
785 do j = g%jsc-1,g%jec ; do i = g%isc,g%iec
786 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 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 cs%coeff_r(:)*cs%Coef_h(i,j+1,:))
799 endif
800 enddo ; enddo
801 else
802 do j = g%jsc-1,g%jec ; do i = g%isc,g%iec
803 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 cs%Coef_h(i,j+1,:))
812 endif
813 enddo ; enddo
814 endif
815 else
816 if (cs%tapering) then
817 do j = g%jsc-1,g%jec ; do i = g%isc,g%iec
818 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 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 cs%coeff_r(:))
831 endif
832 enddo ; enddo
833 else
834 do j = g%jsc-1,g%jec ; do i = g%isc,g%iec
835 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 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 if (cs%KhTh_use_vert_struct) then
851 do j = g%jsc,g%jec ; do i = g%isc,g%iec
852 if (g%mask2dT(i,j)>0.) then
853 if (cs%ndiff_answer_date <= 20240330) then
854 dtracer(:) = 0.
855 do ks = 1,cs%nsurf-1
856 k = cs%uKoL(i,j,ks)
857 dtracer(k) = dtracer(k) + uflx(i,j,ks)
858 k = cs%uKoR(i-1,j,ks)
859 dtracer(k) = dtracer(k) - uflx(i-1,j,ks)
860 k = cs%vKoL(i,j,ks)
861 dtracer(k) = dtracer(k) + vflx(i,j,ks)
862 k = cs%vKoR(i,j-1,ks)
863 dtracer(k) = dtracer(k) - vflx(i,j-1,ks)
864 enddo
865 else ! This form recovers rotational symmetry.
866 dtracer_n(:) = 0.0 ; dtracer_s(:) = 0.0 ; dtracer_e(:) = 0.0 ; dtracer_w(:) = 0.0
867 do ks = 1,cs%nsurf-1
868 k = cs%uKoL(i,j,ks)
869 dtracer_e(k) = dtracer_e(k) + uflx(i,j,ks)
870 k = cs%uKoR(i-1,j,ks)
871 dtracer_w(k) = dtracer_w(k) - uflx(i-1,j,ks)
872 k = cs%vKoL(i,j,ks)
873 dtracer_n(k) = dtracer_n(k) + vflx(i,j,ks)
874 k = cs%vKoR(i,j-1,ks)
875 dtracer_s(k) = dtracer_s(k) - vflx(i,j-1,ks)
876 enddo
877 do k = 1, gv%ke
878 dtracer(k) = (dtracer_n(k) + dtracer_s(k)) + (dtracer_e(k) + dtracer_w(k))
879 enddo
880 endif
881 do k = 1, gv%ke
882 tracer%t(i,j,k) = tracer%t(i,j,k) + dtracer(k) * &
883 ( g%IareaT(i,j) / ( h(i,j,k) + gv%H_subroundoff ) )
884 if (abs(tracer%t(i,j,k)) < tracer%conc_underflow) tracer%t(i,j,k) = 0.0
885 enddo
886
887 if (tracer%id_dfxy_conc > 0 .or. tracer%id_dfxy_cont > 0 .or. tracer%id_dfxy_cont_2d > 0 ) then
888 do k = 1, gv%ke
889 tendency(i,j,k) = dtracer(k) * g%IareaT(i,j) * idt
890 enddo
891 endif
892
893 endif
894 enddo ; enddo
895 else
896 do j = g%jsc,g%jec ; do i = g%isc,g%iec
897 if (g%mask2dT(i,j)>0.) then
898 if (cs%ndiff_answer_date <= 20240330) then
899 dtracer(:) = 0.
900 do ks = 1,cs%nsurf-1
901 k = cs%uKoL(i,j,ks)
902 dtracer(k) = dtracer(k) + coef_x(i,j,1) * uflx(i,j,ks)
903 k = cs%uKoR(i-1,j,ks)
904 dtracer(k) = dtracer(k) - coef_x(i-1,j,1) * uflx(i-1,j,ks)
905 k = cs%vKoL(i,j,ks)
906 dtracer(k) = dtracer(k) + coef_y(i,j,1) * vflx(i,j,ks)
907 k = cs%vKoR(i,j-1,ks)
908 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 dtracer_n(:) = 0.0 ; dtracer_s(:) = 0.0 ; dtracer_e(:) = 0.0 ; dtracer_w(:) = 0.0
912 do ks = 1,cs%nsurf-1
913 k = cs%uKoL(i,j,ks)
914 dtracer_e(k) = dtracer_e(k) + coef_x(i,j,1) * uflx(i,j,ks)
915 k = cs%uKoR(i-1,j,ks)
916 dtracer_w(k) = dtracer_w(k) - coef_x(i-1,j,1) * uflx(i-1,j,ks)
917 k = cs%vKoL(i,j,ks)
918 dtracer_n(k) = dtracer_n(k) + coef_y(i,j,1) * vflx(i,j,ks)
919 k = cs%vKoR(i,j-1,ks)
920 dtracer_s(k) = dtracer_s(k) - coef_y(i,j-1,1) * vflx(i,j-1,ks)
921 enddo
922 do k = 1, gv%ke
923 dtracer(k) = (dtracer_n(k) + dtracer_s(k)) + (dtracer_e(k) + dtracer_w(k))
924 enddo
925 endif
926 do k = 1, gv%ke
927 tracer%t(i,j,k) = tracer%t(i,j,k) + dtracer(k) * &
928 ( g%IareaT(i,j) / ( h(i,j,k) + gv%H_subroundoff ) )
929 if (abs(tracer%t(i,j,k)) < tracer%conc_underflow) tracer%t(i,j,k) = 0.0
930 enddo
931
932 if (tracer%id_dfxy_conc > 0 .or. tracer%id_dfxy_cont > 0 .or. tracer%id_dfxy_cont_2d > 0 ) then
933 do k = 1, gv%ke
934 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 if (tracer%conc_underflow > 0.0) then
944 do k=1,gv%ke ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
945 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 if (tracer%id_dfx_2d > 0) then
952
953 if (cs%KhTh_use_vert_struct) then
954 do j = g%jsc,g%jec ; do i = g%isc-1,g%iec
955 trans_x_2d(i,j) = 0.
956 if (g%mask2dCu(i,j)>0.) then
957 do ks = 1,cs%nsurf-1
958 trans_x_2d(i,j) = trans_x_2d(i,j) - uflx(i,j,ks)
959 enddo
960 trans_x_2d(i,j) = trans_x_2d(i,j) * idt
961 endif
962 enddo ; enddo
963 else
964 do j = g%jsc,g%jec ; do i = g%isc-1,g%iec
965 trans_x_2d(i,j) = 0.
966 if (g%mask2dCu(i,j)>0.) then
967 do ks = 1,cs%nsurf-1
968 trans_x_2d(i,j) = trans_x_2d(i,j) - coef_x(i,j,1) * uflx(i,j,ks)
969 enddo
970 trans_x_2d(i,j) = trans_x_2d(i,j) * idt
971 endif
972 enddo ; enddo
973 endif
974
975 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 if (tracer%id_dfy_2d > 0) then
981
982 if (cs%KhTh_use_vert_struct) then
983 do j = g%jsc-1,g%jec ; do i = g%isc,g%iec
984 trans_y_2d(i,j) = 0.
985 if (g%mask2dCv(i,j)>0.) then
986 do ks = 1,cs%nsurf-1
987 trans_y_2d(i,j) = trans_y_2d(i,j) - vflx(i,j,ks)
988 enddo
989 trans_y_2d(i,j) = trans_y_2d(i,j) * idt
990 endif
991 enddo ; enddo
992 else
993 do j = g%jsc-1,g%jec ; do i = g%isc,g%iec
994 trans_y_2d(i,j) = 0.
995 if (g%mask2dCv(i,j)>0.) then
996 do ks = 1,cs%nsurf-1
997 trans_y_2d(i,j) = trans_y_2d(i,j) - coef_y(i,j,1) * vflx(i,j,ks)
998 enddo
999 trans_y_2d(i,j) = trans_y_2d(i,j) * idt
1000 endif
1001 enddo ; enddo
1002 endif
1003
1004 call post_data(tracer%id_dfy_2d, trans_y_2d(:,:), cs%diag)
1005 endif
1006
1007 ! post tendency of layer-integrated tracer content
1008 if (tracer%id_dfxy_cont > 0) then
1009 call post_data(tracer%id_dfxy_cont, tendency(:,:,:), cs%diag)
1010 endif
1011
1012 ! post depth summed tendency for tracer content
1013 if (tracer%id_dfxy_cont_2d > 0) then
1014 tendency_2d(:,:) = 0.
1015 do j = g%jsc,g%jec ; do i = g%isc,g%iec
1016 do k = 1, gv%ke
1017 tendency_2d(i,j) = tendency_2d(i,j) + tendency(i,j,k)
1018 enddo
1019 enddo ; enddo
1020 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 if (tracer%id_dfxy_conc > 0) then
1027 do k = 1, gv%ke ; do j = g%jsc,g%jec ; do i = g%isc,g%iec
1028 tendency(i,j,k) = tendency(i,j,k) / ( h(i,j,k) + gv%H_subroundoff )
1029 enddo ; enddo ; enddo
1030 call post_data(tracer%id_dfxy_conc, tendency, cs%diag)
1031 endif
1032 enddo ! Loop over tracer registry
1033
1034end 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.
1038subroutine 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 coeff_l(:) = 1.0
1058 coeff_r(:) = 1.0
1059
1060 ! Calculate vertical indices containing the boundary layer depths
1061 max_bld = max(bld_l, bld_r)
1062 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 zeta_l)
1067 call boundary_k_range(surface, ne-1, h_r, min_bld, dummy1, dummy2, k_min_r, &
1068 zeta_r)
1069
1070 ! k_max
1071 call boundary_k_range(surface, ne-1, h_l, max_bld, dummy1, dummy2, k_max_l, &
1072 zeta_l)
1073 call boundary_k_range(surface, ne-1, h_r, max_bld, dummy1, dummy2, k_max_r, &
1074 zeta_r)
1075 ! left
1076 do k=1,k_min_l
1077 coeff_l(k) = 0.0
1078 enddo
1079 do k=k_min_l+1,k_max_l+1
1080 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 do k=1,k_min_r
1085 coeff_r(k) = 0.0
1086 enddo
1087 do k=k_min_r+1,k_max_r+1
1088 coeff_r(k) = (real(k - k_min_r) + 1.0)/(real(k_max_r - k_min_r) + 2.0)
1089 enddo
1090
1091end subroutine compute_tapering_coeffs
1092
1093!> Returns interface scalar, Si, for a column of layer values, S.
1094subroutine 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 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 call plm_diff(nk, h, s, 2, 1, diff)
1112 si(1) = s(1) - 0.5 * diff(1)
1113 if (i_method==1) then
1114 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 sa = s(k-1) + 0.5 * diff(k-1) ! Lower edge value of a PLM reconstruction for layer above
1118 sb = s(k) - 0.5 * diff(k) ! Upper edge value of a PLM reconstruction for layer below
1119 si(k) = 0.5 * ( sa + sb )
1120 enddo
1121 elseif (i_method==2) then
1122 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 km2 = max(1, k-2)
1126 kp1 = min(nk, k+1)
1127 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 si(nk+1) = s(nk) + 0.5 * diff(nk)
1131
1132end 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).
1137real 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 r_hk_hkp1 = hk + hkp1
1158 if (r_hk_hkp1 <= 0.) then
1159 ppm_edge = 0.5 * ( ak + akp1 )
1160 return
1161 endif
1162 r_hk_hkp1 = 1. / r_hk_hkp1
1163 if (hk<hkp1) then
1164 ppm_edge = ak + ( hk * r_hk_hkp1 ) * ( akp1 - ak )
1165 else
1166 ppm_edge = akp1 + ( hkp1 * r_hk_hkp1 ) * ( ak - akp1 )
1167 endif
1168
1169 r_2hk_hkp1 = 1. / ( ( 2. * hk + hkp1 ) + h_neglect )
1170 r_hk_2hkp1 = 1. / ( ( hk + 2. * hkp1 ) + h_neglect )
1171 f1 = 1./ ( ( hk + hkp1) + ( hkm1 + hkp2 ) )
1172 f2 = 2. * ( hkp1 * hk ) * r_hk_hkp1 * &
1173 ( ( hkm1 + hk ) * r_2hk_hkp1 - ( hkp2 + hkp1 ) * r_hk_2hkp1 )
1174 f3 = hk * ( hkm1 + hk ) * r_2hk_hkp1
1175 f4 = hkp1 * ( hkp1 + hkp2 ) * r_hk_2hkp1
1176
1177 ppm_edge = ppm_edge + f1 * ( f2 * ( akp1 - ak ) - ( f3 * pkp1 - f4 * pk ) )
1178
1179end 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)
1183real 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 dx = xr - xl
1200 xave = 0.5 * ( xr + xl )
1201 a6o3 = 2. * amean - ( al + ar ) ! a6 / 3.
1202 a6 = 3. * a6o3
1203
1204 if (dx<0.) then
1205 stop 'ppm_ave: dx<0 should not happened!'
1206 elseif (dx>1.) then
1207 stop 'ppm_ave: dx>1 should not happened!'
1208 elseif (dx==0.) then
1209 ppm_ave = al + ( ar - al ) * xr + a6 * xr * ( 1. - xr )
1210 else
1211 ppm_ave = ( al + xave * ( ( ar - al ) + a6 ) ) - a6o3 * ( xr**2 + xr * xl + xl**2 )
1212 endif
1213end 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].
1217real 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 signum = sign(a,x)
1222 if (x==0.) signum = 0.
1223
1224end 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.
1228subroutine 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 do k = 2, nk-1
1251 hkm1 = h(k-1)
1252 hk = h(k)
1253 hkp1 = h(k+1)
1254
1255 if ( ( hkp1 + hk ) * ( hkm1 + hk ) > 0.) then
1256 skm1 = s(k-1)
1257 sk = s(k)
1258 skp1 = s(k+1)
1259 if (c_method==1) then
1260 ! Simple centered diff (from White)
1261 if ( hk + 0.5 * (hkm1 + hkp1) /= 0. ) then
1262 diff_c = ( skp1 - skm1 ) * ( hk / ( hk + 0.5 * (hkm1 + hkp1) ) )
1263 else
1264 diff_c = 0.
1265 endif
1266 elseif (c_method==2) then
1267 ! Second order accurate centered finite-volume slope (from Colella and Woodward, JCP 1984)
1268 diff_c = fv_diff(hkm1, hk, hkp1, skm1, sk, skp1)
1269 elseif (c_method==3) then
1270 ! Second order accurate finite-volume least squares slope
1271 diff_c = hk * fvlsq_slope(hkm1, hk, hkp1, skm1, sk, skp1)
1272 endif
1273 ! Limit centered slope by twice the side differenced slopes
1274 diff_l = 2. * ( sk - skm1 )
1275 diff_r = 2. * ( skp1 - sk )
1276 if ( signum(1., diff_l) * signum(1., diff_r) <= 0. ) then
1277 diff(k) = 0. ! PCM for local extrema
1278 else
1279 diff(k) = sign( min( abs(diff_l), abs(diff_c), abs(diff_r) ), diff_c )
1280 endif
1281 else
1282 diff(k) = 0. ! PCM next to vanished layers
1283 endif
1284 enddo
1285 if (b_method==1) then ! PCM for top and bottom layer
1286 diff(1) = 0.
1287 diff(nk) = 0.
1288 elseif (b_method==2) then ! Linear extrapolation for top and bottom interfaces
1289 diff(1) = ( s(2) - s(1) ) * 2. * ( h(1) / ( h(1) + h(2) ) )
1290 diff(nk) = s(nk) - s(nk-1) * 2. * ( h(nk) / ( h(nk-1) + h(nk) ) )
1291 endif
1292
1293end 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.
1299real 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 h_sum = ( hkm1 + hkp1 ) + hk
1315 if (h_sum /= 0.) h_sum = 1./ h_sum
1316 hm = hkm1 + hk
1317 if (hm /= 0.) hm = 1./ hm
1318 hp = hkp1 + hk
1319 if (hp /= 0.) hp = 1./ hp
1320 fv_diff = ( hk * h_sum ) * &
1321 ( ( 2. * hkm1 + hk ) * hp * ( skp1 - sk ) &
1322 + ( 2. * hkp1 + hk ) * hm * ( sk - skm1 ) )
1323end 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].
1330real 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 xkm1 = -0.5 * ( hk + hkm1 )
1353 xkp1 = 0.5 * ( hk + hkp1 )
1354 h_sum = ( hkm1 + hkp1 ) + hk
1355 hx_sum = hkm1*xkm1 + hkp1*xkp1
1356 hxsq_sum = hkm1*(xkm1**2) + hkp1*(xkp1**2)
1357 hxy_sum = hkm1*xkm1*skm1 + hkp1*xkp1*skp1
1358 hy_sum = ( hkm1*skm1 + hkp1*skp1 ) + hk*sk
1359 det = h_sum * hxsq_sum - hx_sum**2
1360 if (det /= 0.) then
1361 !a = ( hxsq_sum * hy_sum - hx_sum*hxy_sum ) / det ! a would be mean of straight line fit
1362 fvlsq_slope = ( h_sum * hxy_sum - hx_sum*hy_sum ) / det ! Gradient of straight line fit
1363 else
1364 fvlsq_slope = 0. ! Adcroft's reciprocal rule
1365 endif
1366end function fvlsq_slope
1367
1368
1369!> Returns positions within left/right columns of combined interfaces using continuous reconstructions of T/S
1370subroutine find_neutral_surface_positions_continuous(nk, Pl, Tl, Sl, dRdTl, dRdSl, Pr, Tr, Sr, &
1371 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 ns = 2*nk+2
1412
1413 ! Initialize variables for the search
1414 kr = 1
1415 kl = 1
1416 lastp_right = 0.
1417 lastp_left = 0.
1418 lastk_right = 1
1419 lastk_left = 1
1420 reached_bottom = .false.
1421
1422 ! Check to see if we should limit the diffusion to the interior
1423 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 neutral_surfaces: do k_surface = 1, ns
1427 klm1 = max(kl-1, 1)
1428 if (klm1>nk) stop 'find_neutral_surface_positions(): klm1 went out of bounds!'
1429 krm1 = max(kr-1, 1)
1430 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 + ( drdsr(kr) + drdsl(kl) ) * ( sr(kr) - sl(kl) ) )
1435 ! Which column has the lighter surface for the current indexes, kr and kl
1436 if (.not. reached_bottom) then
1437 if (drho < 0.) then
1438 searching_left_column = .true.
1439 searching_right_column = .false.
1440 elseif (drho > 0.) then
1441 searching_right_column = .true.
1442 searching_left_column = .false.
1443 else ! dRho == 0.
1444 if (kl + kr == 2) then ! Still at surface
1445 searching_left_column = .true.
1446 searching_right_column = .false.
1447 else ! Not the surface so we simply change direction
1448 searching_left_column = .not. searching_left_column
1449 searching_right_column = .not. searching_right_column
1450 endif
1451 endif
1452 endif
1453
1454 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 + ( 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 + ( 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 if (drhotop > 0. .or. kr+kl==2) then
1466 pol(k_surface) = 0. ! The right surface is lighter than anything in layer klm1
1467 elseif (drhotop >= drhobot) then ! Left layer is unstratified
1468 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 pol(k_surface) = interpolate_for_nondim_position( drhotop, pl(klm1), drhobot, pl(klm1+1) )
1473 endif
1474 if (pol(k_surface)>=1. .and. klm1<nk) then ! >= is really ==, when PoL==1 we point to the bottom of the cell
1475 klm1 = klm1 + 1
1476 pol(k_surface) = pol(k_surface) - 1.
1477 endif
1478 if (real(klm1-lastk_left)+(pol(k_surface)-lastp_left)<0.) then
1479 pol(k_surface) = lastp_left
1480 klm1 = lastk_left
1481 endif
1482 kol(k_surface) = klm1
1483 if (kr <= nk) then
1484 por(k_surface) = 0.
1485 kor(k_surface) = kr
1486 else
1487 por(k_surface) = 1.
1488 kor(k_surface) = nk
1489 endif
1490 if (kr <= nk) then
1491 kr = kr + 1
1492 else
1493 reached_bottom = .true.
1494 searching_right_column = .true.
1495 searching_left_column = .false.
1496 endif
1497 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 ( 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 ( 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 if (drhotop >= 0. .or. kr+kl==2) then
1509 por(k_surface) = 0. ! The left surface is lighter than anything in layer krm1
1510 elseif (drhotop >= drhobot) then ! Right layer is unstratified
1511 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 por(k_surface) = interpolate_for_nondim_position( drhotop, pr(krm1), drhobot, pr(krm1+1) )
1516 endif
1517 if (por(k_surface)>=1. .and. krm1<nk) then ! >= is really ==, when PoR==1 we point to the bottom of the cell
1518 krm1 = krm1 + 1
1519 por(k_surface) = por(k_surface) - 1.
1520 endif
1521 if (real(krm1-lastk_right)+(por(k_surface)-lastp_right)<0.) then
1522 por(k_surface) = lastp_right
1523 krm1 = lastk_right
1524 endif
1525 kor(k_surface) = krm1
1526 if (kl <= nk) then
1527 pol(k_surface) = 0.
1528 kol(k_surface) = kl
1529 else
1530 pol(k_surface) = 1.
1531 kol(k_surface) = nk
1532 endif
1533 if (kl <= nk) then
1534 kl = kl + 1
1535 else
1536 reached_bottom = .true.
1537 searching_right_column = .false.
1538 searching_left_column = .true.
1539 endif
1540 else
1541 stop 'Else what?'
1542 endif
1543 if (interior_limit) then
1544 if (kol(k_surface)<=bl_kl) then
1545 kol(k_surface) = bl_kl
1546 if (pol(k_surface)<bl_zl) then
1547 pol(k_surface) = bl_zl
1548 endif
1549 endif
1550 if (kor(k_surface)<=bl_kr) then
1551 kor(k_surface) = bl_kr
1552 if (por(k_surface)<bl_zr) then
1553 por(k_surface) = bl_zr
1554 endif
1555 endif
1556 endif
1557
1558 lastk_left = kol(k_surface) ; lastp_left = pol(k_surface)
1559 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 if (k_surface>1) then
1564 hl = absolute_position(nk,ns,pl,kol,pol,k_surface) - absolute_position(nk,ns,pl,kol,pol,k_surface-1)
1565 hr = absolute_position(nk,ns,pr,kor,por,k_surface) - absolute_position(nk,ns,pr,kor,por,k_surface-1)
1566 if ( hl + hr > 0.) then
1567 heff(k_surface-1) = 2. * hl * hr / ( hl + hr ) ! Harmonic mean of layer thicknesses
1568 else
1569 heff(k_surface-1) = 0.
1570 endif
1571 endif
1572
1573 enddo neutral_surfaces
1574
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.
1580real 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 if ((ppos > pneg) .and. (drhopos - drhoneg >= 0. )) then
1589 if ( drhopos - drhoneg > 0. ) then
1590 interpolate_for_nondim_position = min( 1., max( 0., -drhoneg / ( drhopos - drhoneg ) ) )
1591 elseif (drhopos - drhoneg == 0) then
1592 if (drhoneg > 0.) then
1594 elseif (drhoneg < 0.) then
1596 else ! dRhoPos = dRhoNeg = 0
1598 endif
1599 else ! dRhoPos - dRhoNeg < 0
1601 endif
1602 elseif (ppos == pneg) then ! Handle vanished or inverted layers
1604 else ! ((Ppos < Pneg) .or. (dRhoNeg > dRhoPos) )
1605 ! Error handling for problematic cases. It is expected that this should never occur.
1606 write(mesg,*) 'dRhoNeg, Pneg, dRhoPos, Ppos', drhoneg, pneg, drhopos, ppos
1607 call mom_error(warning, 'interpolate_for_nondim_position: '//trim(mesg))
1608 ! write(stderr,*) trim(mesg)
1609 if ((ppos < pneg) .and. (drhoneg > drhopos)) then
1610 mesg = '(Ppos < Pneg) and (dRhoNeg > dRhoPos)'
1611 elseif (ppos < pneg) then
1612 mesg = 'Ppos < Pneg'
1613 elseif (drhoneg > drhopos) then
1614 mesg = trim(mesg)//'; dRhoNeg > dRhoPos'
1615 else ! This should never happen.
1616 mesg = 'Unexpected failure.'
1617 endif
1618 call mom_error(fatal, 'interpolate_for_nondim_position: '//trim(mesg))
1619 endif
1620
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
1627 Pres_l, hcol_l, Tl, Sl, ppoly_T_l, ppoly_S_l, stable_l, &
1628 Pres_r, hcol_r, Tr, Sr, ppoly_T_r, ppoly_S_r, stable_r, &
1629 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 ns = 4*nk
1686 ki_right = 1
1687 ki_left = 1
1688 kl_left = 1
1689 kl_right = 1
1690 lastp_left = 0.
1691 lastp_right = 0.
1692 reached_bottom = .false.
1693 searching_left_column = .false.
1694 searching_right_column = .false.
1695
1696 fail_heff = .true.
1697 if (PRESENT(hard_fail_heff)) fail_heff = hard_fail_heff
1698
1699 if (PRESENT(k_bot_l) .and. PRESENT(k_bot_r) .and. PRESENT(zeta_bot_l) .and. PRESENT(zeta_bot_r)) then
1700 k_init_l = k_bot_l ; k_init_r = k_bot_r
1701 p_init_l = zeta_bot_l ; p_init_r = zeta_bot_r
1702 lastp_left = zeta_bot_l ; lastp_right = zeta_bot_r
1703 kl_left = k_bot_l ; kl_right = k_bot_r
1704 else
1705 k_init_l = 1 ; k_init_r = 1
1706 p_init_l = 0. ; p_init_r = 0.
1707 endif
1708 ! Loop over each neutral surface, working from top to bottom
1709 neutral_surfaces: do k_surface = 1, ns
1710
1711 if (k_surface == ns) then
1712 pol(k_surface) = 1.
1713 por(k_surface) = 1.
1714 kol(k_surface) = nk
1715 kor(k_surface) = nk
1716 ! If the layers are unstable, then simply point the surface to the previous location
1717 elseif (.not. stable_l(kl_left)) then
1718 if (k_surface > 1) then
1719 pol(k_surface) = ki_left - 1 ! Top interface is at position = 0., Bottom is at position = 1
1720 kol(k_surface) = kl_left
1721 por(k_surface) = por(k_surface-1)
1722 kor(k_surface) = kor(k_surface-1)
1723 else
1724 por(k_surface) = p_init_r
1725 kor(k_surface) = k_init_r
1726 pol(k_surface) = p_init_l
1727 kol(k_surface) = k_init_l
1728 endif
1729 call increment_interface(nk, kl_left, ki_left, reached_bottom, searching_left_column, searching_right_column)
1730 searching_left_column = .true.
1731 searching_right_column = .false.
1732 elseif (.not. stable_r(kl_right)) then ! Check the right layer for stability
1733 if (k_surface > 1) then
1734 por(k_surface) = ki_right - 1 ! Top interface is at position = 0., Bottom is at position = 1
1735 kor(k_surface) = kl_right
1736 pol(k_surface) = pol(k_surface-1)
1737 kol(k_surface) = kol(k_surface-1)
1738 else
1739 por(k_surface) = 0.
1740 kor(k_surface) = 1
1741 pol(k_surface) = 0.
1742 kol(k_surface) = 1
1743 endif
1744 call increment_interface(nk, kl_right, ki_right, reached_bottom, searching_right_column, searching_left_column)
1745 searching_left_column = .false.
1746 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 drho)
1755 if (cs%debug) write(stdout,'(A,I0,A,E12.4,A,I0,A,I0,A,I0,A,I0)') &
1756 "k_surface=",k_surface, " dRho=",cs%R_to_kg_m3*drho, &
1757 "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 if (.not. reached_bottom) then
1760 if (drho < 0.) then
1761 searching_left_column = .true.
1762 searching_right_column = .false.
1763 elseif (drho > 0.) then
1764 searching_left_column = .false.
1765 searching_right_column = .true.
1766 else ! dRho == 0.
1767 if ( ( kl_left + kl_right == 2 ) .and. (ki_left + ki_right == 2) ) then ! Still at surface
1768 searching_left_column = .true.
1769 searching_right_column = .false.
1770 else ! Not the surface so we simply change direction
1771 searching_left_column = .not. searching_left_column
1772 searching_right_column = .not. searching_right_column
1773 endif
1774 endif
1775 endif
1776 if (searching_left_column) then
1777 ! Position of the right interface is known and all quantities are fixed
1778 por(k_surface) = ki_right - 1.
1779 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 ppoly_t_l(kl_left,:), ppoly_s_l(kl_left,:))
1785 kol(k_surface) = kl_left
1786
1787 if (cs%debug) then
1788 write(stdout,'(A,I0)') "Searching left layer ", kl_left
1789 write(stdout,'(A,I0,1X,I0)') "Searching from right: ", kl_right, ki_right
1790 write(stdout,*) "Temp/Salt Reference: ", tr(kl_right,ki_right), sr(kl_right,ki_right)
1791 write(stdout,*) "Temp/Salt Top L: ", tl(kl_left,1), sl(kl_left,1)
1792 write(stdout,*) "Temp/Salt Bot L: ", tl(kl_left,2), sl(kl_left,2)
1793 endif
1794 call increment_interface(nk, kl_right, ki_right, reached_bottom, searching_right_column, searching_left_column)
1795 lastp_left = pol(k_surface)
1796 ! If the right layer increments, then we need to reset the last position on the right
1797 if ( kl_right == (kor(k_surface) + 1) ) lastp_right = 0.
1798
1799 elseif (searching_right_column) then
1800 ! Position of the right interface is known and all quantities are fixed
1801 pol(k_surface) = ki_left - 1.
1802 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 ppoly_t_r(kl_right,:), ppoly_s_r(kl_right,:))
1808 kor(k_surface) = kl_right
1809
1810 if (cs%debug) then
1811 write(stdout,'(A,I0)') "Searching right layer ", kl_right
1812 write(stdout,'(A,I0,1X,I0)') "Searching from left: ", kl_left, ki_left
1813 write(stdout,*) "Temp/Salt Reference: ", tl(kl_left,ki_left), sl(kl_left,ki_left)
1814 write(stdout,*) "Temp/Salt Top L: ", tr(kl_right,1), sr(kl_right,1)
1815 write(stdout,*) "Temp/Salt Bot L: ", tr(kl_right,2), sr(kl_right,2)
1816 endif
1817 call increment_interface(nk, kl_left, ki_left, reached_bottom, searching_left_column, searching_right_column)
1818 lastp_right = por(k_surface)
1819 ! If the right layer increments, then we need to reset the last position on the right
1820 if ( kl_left == (kol(k_surface) + 1) ) lastp_left = 0.
1821 else
1822 stop 'Else what?'
1823 endif
1824 if (cs%debug) write(stdout,'(A,I3,A,ES16.6,A,I0,A,ES16.6)') "KoL:", kol(k_surface), " PoL:", pol(k_surface), &
1825 " KoR:", kor(k_surface), " PoR:", por(k_surface)
1826 endif
1827 ! Effective thickness
1828 if (k_surface>1) then
1829 if ( kol(k_surface) == kol(k_surface-1) .and. kor(k_surface) == kor(k_surface-1) ) then
1830 hl = (pol(k_surface) - pol(k_surface-1))*hcol_l(kol(k_surface))
1831 hr = (por(k_surface) - por(k_surface-1))*hcol_r(kor(k_surface))
1832 if (hl < 0. .or. hr < 0.) then
1833 if (fail_heff) then
1834 call mom_error(fatal,"Negative thicknesses in neutral diffusion")
1835 else
1836 if (searching_left_column) then
1837 pol(k_surface) = pol(k_surface-1)
1838 kol(k_surface) = kol(k_surface-1)
1839 elseif (searching_right_column) then
1840 por(k_surface) = por(k_surface-1)
1841 kor(k_surface) = kor(k_surface-1)
1842 endif
1843 endif
1844 elseif ( hl + hr == 0. ) then
1845 heff(k_surface-1) = 0.
1846 else
1847 heff(k_surface-1) = 2. * ( (hl * hr) / ( hl + hr ) )! Harmonic mean
1848 if ( kol(k_surface) /= kol(k_surface-1) ) then
1849 call mom_error(fatal,"Neutral sublayer spans multiple layers")
1850 endif
1851 if ( kor(k_surface) /= kor(k_surface-1) ) then
1852 call mom_error(fatal,"Neutral sublayer spans multiple layers")
1853 endif
1854 endif
1855 else
1856 heff(k_surface-1) = 0.
1857 endif
1858 endif
1859 enddo neutral_surfaces
1861
1862!> Sweep down through the column and mark as stable if the bottom interface of a cell is denser than the top
1863subroutine 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 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 t(k,1), s(k,1), max(p(k,1), cs%ref_pres), delta_rho )
1877 stable_cell(k) = (delta_rho > 0.)
1878 enddo
1879end 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]
1883real function search_other_column(CS, ksurf, pos_last, T_from, S_from, P_from, T_top, S_top, P_top, &
1884 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 if (cs%neutral_pos_method == 1 .or. cs%neutral_pos_method == 3) then
1911 call calc_delta_rho_and_derivs(cs, t_top, s_top, p_top, t_from, s_from, p_from, drhotop)
1912 call calc_delta_rho_and_derivs(cs, t_bot, s_bot, p_bot, t_from, s_from, p_from, drhobot)
1913 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 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 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 if ( (drhotop > 0.) .or. (ksurf == 1) ) then ! First interface or lighter than anything in layer
1922 pos = pos_last
1923 elseif ( drhotop > drhobot ) then ! Unstably stratified
1924 pos = 1.
1925 elseif ( drhotop < 0. .and. drhobot < 0.) then ! Denser than anything in layer
1926 pos = 1.
1927 elseif ( drhotop == 0. .and. drhobot == 0. ) then ! Perfectly unstratified
1928 pos = 1.
1929 elseif ( drhobot == 0. ) then ! Matches perfectly at the Top
1930 pos = 1.
1931 elseif ( drhotop == 0. ) then ! Matches perfectly at the Bottom
1932 pos = pos_last
1933 else ! Neutral surface within layer
1934 pos = -1
1935 endif
1936
1937 ! Can safely return if position is >= 0 otherwise will need to find the position within the layer
1938 if (pos>=0) return
1939
1940 if (cs%neutral_pos_method==1) then
1941 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 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 drdt_top, drds_top, drdt_bot, drds_bot, t_poly, s_poly )
1947 elseif (cs%neutral_pos_method == 3) then
1948 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
1951end function search_other_column
1952
1953!> Increments the interface which was just connected and also set flags if the bottom is reached
1954subroutine 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 reached_bottom = .false.
1963 if (ki == 2) then ! At the bottom interface
1964 if ((ki == 2) .and. (kl < nk) ) then ! Not at the bottom so just go to the next layer
1965 kl = kl+1
1966 ki = 1
1967 elseif ((kl == nk) .and. (ki==2)) then
1968 reached_bottom = .true.
1969 searching_this_column = .false.
1970 searching_other_column = .true.
1971 endif
1972 elseif (ki==1) then ! At the top interface
1973 ki = 2 ! Next interface is same layer, but bottom interface
1974 else
1975 call mom_error(fatal,"Unanticipated eventuality in increment_interface")
1976 endif
1977end 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
1987function find_neutral_pos_linear( CS, z0, T_ref, S_ref, dRdT_ref, dRdS_ref, &
1988 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 nterm = SIZE(ppoly_t)
2028
2029 ! Position independent quantities
2030 drdt_diff = drdt_bot - drdt_top
2031 drds_diff = drds_bot - drds_top
2032 ! Initial starting drho (used for bisection)
2033 zmin = z0 ! Lower bounding interval
2034 zmax = 1. ! Maximum bounding interval (bottom of layer)
2035 a1 = 1. - zmin
2036 a2 = zmin
2037 t_z = evaluation_polynomial( ppoly_t, nterm, zmin )
2038 s_z = evaluation_polynomial( ppoly_s, nterm, zmin )
2039 drdt_z = a1*drdt_top + a2*drdt_bot
2040 drds_z = a1*drds_top + a2*drds_bot
2041 drho_min = 0.5*((drdt_z+drdt_ref)*(t_z-t_ref) + (drds_z+drds_ref)*(s_z-s_ref))
2042
2043 t_z = evaluation_polynomial( ppoly_t, nterm, 1. )
2044 s_z = evaluation_polynomial( ppoly_s, nterm, 1. )
2045 drho_max = 0.5*((drdt_bot+drdt_ref)*(t_z-t_ref) + (drds_bot+drds_ref)*(s_z-s_ref))
2046
2047 if (drho_min >= 0.) then
2048 z = z0
2049 return
2050 elseif (drho_max == 0.) then
2051 z = 1.
2052 return
2053 endif
2054 if ( sign(1.,drho_min) == sign(1.,drho_max) ) then
2055 call mom_error(fatal, "drho_min is the same sign as dhro_max")
2056 endif
2057
2058 z = z0
2059 ztest = z0
2060 do iter = 1, cs%max_iter
2061 ! Calculate quantities at the current nondimensional position
2062 a1 = 1.-z
2063 a2 = z
2064 drdt_z = a1*drdt_top + a2*drdt_bot
2065 drds_z = a1*drds_top + a2*drds_bot
2066 t_z = evaluation_polynomial( ppoly_t, nterm, z )
2067 s_z = evaluation_polynomial( ppoly_s, nterm, z )
2068 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 if (abs(drho) <= cs%drho_tol) exit
2072 ! Update bisection bracketing intervals
2073 if (drho < 0. .and. drho > drho_min) then
2074 drho_min = drho
2075 zmin = z
2076 elseif (drho > 0. .and. drho < drho_max) then
2077 drho_max = drho
2078 zmax = z
2079 endif
2080
2081 ! Calculate a Newton step
2082 dt_dz = first_derivative_polynomial( ppoly_t, nterm, z )
2083 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 (drds_diff*(s_z - s_ref) + (drds_ref+drds_z)*ds_dz) )
2086
2087 ztest = z - drho/drho_dz
2088 ! Take a bisection if z falls out of [zmin,zmax]
2089 if (ztest < zmin .or. ztest > zmax) then
2090 if ( drho < 0. ) then
2091 ztest = 0.5*(z + zmax)
2092 else
2093 ztest = 0.5*(zmin + z)
2094 endif
2095 endif
2096
2097 ! Test to ensure we haven't stalled out
2098 if ( abs(z-ztest) <= cs%x_tol ) exit
2099 ! Reset for next iteration
2100 z = ztest
2101 enddo
2102
2103end 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
2107function 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 side = 0
2133 ! Set the first two evaluation to the endpoints of the interval
2134 b = z0 ; c = 1
2135 nterm = SIZE(ppoly_t)
2136
2137 ! Calculate drho at the minimum bound
2138 tb = evaluation_polynomial( ppoly_t, nterm, b )
2139 sb = evaluation_polynomial( ppoly_s, nterm, b )
2140 pb = p_top*(1.-b) + p_bot*b
2141 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 tc = evaluation_polynomial( ppoly_t, nterm, 1. )
2145 sc = evaluation_polynomial( ppoly_s, nterm, 1. )
2146 pc = p_bot
2147 call calc_delta_rho_and_derivs(cs, tc, sc, pc, t_ref, s_ref, p_ref, drho_c)
2148
2149 if (drho_b >= 0.) then
2150 z = z0
2151 return
2152 elseif (drho_c == 0.) then
2153 z = 1.
2154 return
2155 endif
2156 if ( sign(1.,drho_b) == sign(1.,drho_c) ) then
2157 z = z0
2158 return
2159 endif
2160
2161 do iter = 1, cs%max_iter
2162 ! Calculate new position and evaluate if we have converged
2163 a = (drho_b*c - drho_c*b)/(drho_b-drho_c)
2164 ta = evaluation_polynomial( ppoly_t, nterm, a )
2165 sa = evaluation_polynomial( ppoly_s, nterm, a )
2166 pa = p_top*(1.-a) + p_bot*a
2167 call calc_delta_rho_and_derivs(cs, ta, sa, pa, t_ref, s_ref, p_ref, drho_a)
2168 if (abs(drho_a) < cs%drho_tol) then
2169 z = a
2170 return
2171 endif
2172
2173 if (drho_a*drho_c > 0.) then
2174 if ( abs(a-c)<cs%x_tol) then
2175 z = a
2176 return
2177 endif
2178 c = a ; drho_c = drho_a
2179 if (side == -1) drho_b = 0.5*drho_b
2180 side = -1
2181 elseif ( drho_b*drho_a > 0 ) then
2182 if ( abs(a-b)<cs%x_tol) then
2183 z = a
2184 return
2185 endif
2186 b = a ; drho_b = drho_a
2187 if (side == 1) drho_c = 0.5*drho_c
2188 side = 1
2189 else
2190 z = a
2191 return
2192 endif
2193 enddo
2194
2195 z = a
2196
2197end function find_neutral_pos_full
2198
2199!> Calculate the difference in density between two points in a variety of ways
2200subroutine 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 if (cs%ref_pres > 0.) then
2222 p1 = cs%ref_pres
2223 p2 = cs%ref_pres
2224 else
2225 p1 = p1_in
2226 p2 = p2_in
2227 endif
2228
2229 ! Use the full linear equation of state to calculate the difference in density (expensive!)
2230 if (trim(cs%delta_rho_form) == 'full') then
2231 pmid = 0.5 * (p1 + p2)
2232 call calculate_density(t1, s1, pmid, rho1, cs%EOS)
2233 call calculate_density(t2, s2, pmid, rho2, cs%EOS)
2234 drho = rho1 - rho2
2235 ! Use the density derivatives at the average of pressures and the differences in temperature
2236 elseif (trim(cs%delta_rho_form) == 'mid_pressure') then
2237 pmid = 0.5 * (p1 + p2)
2238 if (cs%ref_pres>=0) pmid = cs%ref_pres
2239 call calculate_density_derivs(t1, s1, pmid, drdt1, drds1, cs%EOS)
2240 call calculate_density_derivs(t2, s2, pmid, drdt2, drds2, cs%EOS)
2241 drho = delta_rho_from_derivs( t1, s1, p1, drdt1, drds1, t2, s2, p2, drdt2, drds2)
2242 elseif (trim(cs%delta_rho_form) == 'local_pressure') then
2243 call calculate_density_derivs(t1, s1, p1, drdt1, drds1, cs%EOS)
2244 call calculate_density_derivs(t2, s2, p2, drdt2, drds2, cs%EOS)
2245 drho = delta_rho_from_derivs( t1, s1, p1, drdt1, drds1, t2, s2, p2, drdt2, drds2)
2246 else
2247 call mom_error(fatal, "delta_rho_form is not recognized")
2248 endif
2249
2250 if (PRESENT(drdt1_out)) drdt1_out = drdt1
2251 if (PRESENT(drds1_out)) drds1_out = drds1
2252 if (PRESENT(drdt2_out)) drdt2_out = drdt2
2253 if (PRESENT(drds2_out)) drds2_out = drds2
2254
2255end 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$
2261function 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 drho = 0.5 * ( (drdt1+drdt2)*(t1-t2) + (drds1+drds2)*(s1-s2))
2277
2278end function delta_rho_from_derivs
2279
2280!> Converts non-dimensional position within a layer to absolute position (for debugging)
2281function 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 k = karr(k_surface)
2294 if (k>n) stop 'absolute_position: k>nk is out of bounds!'
2295 absolute_position = pint(k) + nparr(k_surface) * ( pint(k+1) - pint(k) )
2296
2297end function absolute_position
2298
2299!> Converts non-dimensional positions within layers to absolute positions (for debugging)
2300function 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 do k_surface = 1, ns
2314 absolute_positions(k_surface) = absolute_position(n,ns,pint,karr,nparr,k_surface)
2315 enddo
2316
2317end function absolute_positions
2318
2319!> Returns a single column of neutral diffusion fluxes of a tracer.
2320subroutine neutral_surface_flux(nk, nsurf, deg, hl, hr, Tl, Tr, PiL, PiR, KoL, KoR, &
2321 hEff, Flx, continuous, h_neglect, remap_CS, h_neglect_edge, &
2322 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 real, dimension(nk+1) :: Til !< Left-column interface tracer in arbitrary concentration
2378 !! units (e.g. [C ~> degC] for temperature)
2379 real, dimension(nk+1) :: Tir !< Right-column interface tracer in arbitrary concentration
2380 !! units (e.g. [C ~> degC] for temperature)
2381 real, dimension(nk) :: aL_l !< Left-column left edge value of tracer in arbitrary concentration
2382 !! units (e.g. [C ~> degC] for temperature)
2383 real, dimension(nk) :: aR_l !< Left-column right edge value of tracer in arbitrary concentration
2384 !! units (e.g. [C ~> degC] for temperature)
2385 real, dimension(nk) :: aL_r !< Right-column left edge value of tracer in arbitrary concentration
2386 !! units (e.g. [C ~> degC] for temperature)
2387 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 real, dimension(nk,2) :: Tid_l !< Left-column interface tracer in arbitrary concentration
2392 !! units (e.g. [C ~> degC] for temperature)
2393 real, dimension(nk,2) :: Tid_r !< Right-column interface tracer in arbitrary concentration
2394 !! units (e.g. [C ~> degC] for temperature)
2395 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 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 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 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 tapering = .false.
2408 if (present(coeff_l) .and. present(coeff_r)) tapering = .true.
2409 khtr_ave = 1.0
2410
2411 ! Setup reconstruction edge values
2412 if (continuous) then
2413 call interface_scalar(nk, hl, tl, til, 2, h_neglect)
2414 call interface_scalar(nk, hr, tr, tir, 2, h_neglect)
2415 call ppm_left_right_edge_values(nk, tl, til, al_l, ar_l)
2416 call ppm_left_right_edge_values(nk, tr, tir, al_r, ar_r)
2417 else
2418 ppoly_r_coeffs_l(:,:) = 0.
2419 ppoly_r_coeffs_r(:,:) = 0.
2420 tid_l(:,:) = 0.
2421 tid_r(:,:) = 0.
2422
2423 call build_reconstructions_1d( remap_cs, nk, hl, tl, ppoly_r_coeffs_l, tid_l, &
2424 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 ppoly_r_s_r, imethod, h_neglect, h_neglect_edge )
2427 endif
2428
2429 do k_sublayer = 1, nsurf-1
2430 if (heff(k_sublayer) == 0.) then
2431 flx(k_sublayer) = 0.
2432 else
2433 if (tapering) then
2434 klb = kol(k_sublayer+1)
2435 klt = kol(k_sublayer)
2436 krb = kor(k_sublayer+1)
2437 krt = kor(k_sublayer)
2438 ! these are added in this order to preserve vertically-uniform diffusivity answers
2439 khtr_ave = 0.25 * ((coeff_l(klb) + coeff_l(klt)) + (coeff_r(krb) + coeff_r(krt)))
2440 endif
2441 if (continuous) then
2442 klb = kol(k_sublayer+1)
2443 t_left_bottom = ( 1. - pil(k_sublayer+1) ) * til(klb) + pil(k_sublayer+1) * til(klb+1)
2444 klt = kol(k_sublayer)
2445 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 al_l(klt), ar_l(klt), tl(klt))
2448
2449 krb = kor(k_sublayer+1)
2450 t_right_bottom = ( 1. - pir(k_sublayer+1) ) * tir(krb) + pir(k_sublayer+1) * tir(krb+1)
2451 krt = kor(k_sublayer)
2452 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 al_r(krt), ar_r(krt), tr(krt))
2455 dt_top = t_right_top - t_left_top
2456 dt_bottom = t_right_bottom - t_left_bottom
2457 dt_ave = 0.5 * ( dt_top + dt_bottom )
2458 dt_layer = t_right_layer - t_left_layer
2459 if (signum(1.,dt_top) * signum(1.,dt_bottom) <= 0. .or. signum(1.,dt_ave) * signum(1.,dt_layer) <= 0.) then
2460 dt_ave = 0.
2461 else
2462 dt_ave = dt_layer
2463 endif
2464 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 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 t_right_top_int, t_right_bot_int, t_right_layer)
2473
2474 dt_top = t_right_top - t_left_top
2475 dt_bottom = t_right_bottom - t_left_bottom
2476 dt_sublayer = t_right_sub - t_left_sub
2477 dt_top_int = t_right_top_int - t_left_top_int
2478 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 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 dt_bot_int >= 0.)
2489 if (down_flux) then
2490 flx(k_sublayer) = dt_sublayer * heff(k_sublayer) * khtr_ave
2491 else
2492 flx(k_sublayer) = 0.
2493 endif
2494 endif
2495 endif
2496 enddo
2497
2498end subroutine neutral_surface_flux
2499
2500!> Evaluate various parts of the reconstructions to calculate gradient-based flux limiter
2501subroutine 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 ks_top = k_sub
2533 ks_bot = k_sub + 1
2534 if ( ks(ks_top) /= ks(ks_bot) ) then
2535 call mom_error(fatal, "Neutral surfaces span more than one layer")
2536 endif
2537 kl = ks(k_sub)
2538 ! First if the neutral surfaces spans the entirety of a cell, then do not search across the discontinuity
2539 if ( (ps(ks_top) == 0.) .and. (ps(ks_bot) == 1.)) then
2540 t_top = t_int(kl,1)
2541 t_bot = t_int(kl,2)
2542 else
2543 ! Search across potential discontinuity at top
2544 if ( (kl > 1) .and. (ps(ks_top) == 0.) ) then
2545 t_top = t_int(kl-1,2)
2546 else
2547 t_top = evaluation_polynomial( t_poly(kl,:), deg+1, ps(ks_top) )
2548 endif
2549 ! Search across potential discontinuity at bottom
2550 if ( (kl < nk) .and. (ps(ks_bot) == 1.) ) then
2551 t_bot = t_int(kl+1,1)
2552 else
2553 t_bot = evaluation_polynomial( t_poly(kl,:), deg+1, ps(ks_bot) )
2554 endif
2555 endif
2556 t_sub = average_value_ppoly(nk, t_mean, t_int, t_poly, imethod, kl, ps(ks_top), ps(ks_bot))
2557 t_top_int = evaluation_polynomial( t_poly(kl,:), deg+1, ps(ks_top))
2558 t_bot_int = evaluation_polynomial( t_poly(kl,:), deg+1, ps(ks_bot))
2559 t_layer = t_mean(kl)
2560
2561end subroutine neutral_surface_t_eval
2562
2563!> Discontinuous PPM reconstructions of the left/right edge values within a cell
2564subroutine 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 do k = 1, nk
2576 al(k) = ti(k)
2577 ar(k) = ti(k+1)
2578 if ( signum(1., ar(k) - tl(k))*signum(1., tl(k) - al(k)) <= 0.0 ) then
2579 al(k) = tl(k)
2580 ar(k) = tl(k)
2581 elseif ( sign(3., ar(k) - al(k)) * ( (tl(k) - al(k)) + (tl(k) - ar(k))) > abs(ar(k) - al(k)) ) then
2582 al(k) = tl(k) + 2.0 * ( tl(k) - ar(k) )
2583 elseif ( sign(3., ar(k) - al(k)) * ( (tl(k) - al(k)) + (tl(k) - ar(k))) < -abs(ar(k) - al(k)) ) then
2584 ar(k) = tl(k) + 2.0 * ( tl(k) - al(k) )
2585 endif
2586 enddo
2587end subroutine ppm_left_right_edge_values
2588
2589!> Returns true if unit tests of neutral_diffusion functions fail. Otherwise returns false.
2590logical 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. &
2595
2597
2598!> Returns true if unit tests of neutral_diffusion functions fail. Otherwise returns false.
2599logical 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 h_neglect = 1.0e-30
2612
2613 v = verbose
2614
2615 ndiff_unit_tests_continuous = .false. ! Normally return false
2616 write(stdout,*) '==== MOM_neutral_diffusion: ndiff_unit_tests_continuous ='
2617
2619 test_fv_diff(v,1.,1.,1., 0.,1.,2., 1., 'FV: Straight line on uniform grid')
2621 test_fv_diff(v,1.,1.,0., 0.,4.,8., 7., 'FV: Vanished right cell')
2623 test_fv_diff(v,0.,1.,1., 0.,4.,8., 7., 'FV: Vanished left cell')
2625 test_fv_diff(v,1.,2.,4., 0.,3.,9., 4., 'FV: Stretched grid')
2627 test_fv_diff(v,2.,0.,2., 0.,1.,2., 0., 'FV: Vanished middle cell')
2629 test_fv_diff(v,0.,1.,0., 0.,1.,2., 2., 'FV: Vanished on both sides')
2631 test_fv_diff(v,1.,0.,0., 0.,1.,2., 0., 'FV: Two vanished cell sides')
2633 test_fv_diff(v,0.,0.,0., 0.,1.,2., 0., 'FV: All vanished cells')
2634
2636 test_fvlsq_slope(v,1.,1.,1., 0.,1.,2., 1., 'LSQ: Straight line on uniform grid')
2638 test_fvlsq_slope(v,1.,1.,0., 0.,1.,2., 1., 'LSQ: Vanished right cell')
2640 test_fvlsq_slope(v,0.,1.,1., 0.,1.,2., 1., 'LSQ: Vanished left cell')
2642 test_fvlsq_slope(v,1.,2.,4., 0.,3.,9., 2., 'LSQ: Stretched grid')
2644 test_fvlsq_slope(v,1.,0.,1., 0.,1.,2., 2., 'LSQ: Vanished middle cell')
2646 test_fvlsq_slope(v,0.,1.,0., 0.,1.,2., 0., 'LSQ: Vanished on both sides')
2648 test_fvlsq_slope(v,1.,0.,0., 0.,1.,2., 0., 'LSQ: Two vanished cell sides')
2650 test_fvlsq_slope(v,0.,0.,0., 0.,1.,2., 0., 'LSQ: All vanished cells')
2651
2652 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')
2656 test_data1d(v,5, tio, (/24.,22.5,15.,7.5,6./), 'Linear profile, linear interface temperatures')
2657 call interface_scalar(4, (/10.,10.,10.,10./), (/24.,18.,12.,6./), tio, 2, h_neglect)
2659 test_data1d(v,5, tio, (/24.,22.,15.,8.,6./), 'Linear profile, PPM interface temperatures')
2660
2662 test_ifndp(v,-1.0, 0., 1.0, 1.0, 0.5, 'Check mid-point')
2664 test_ifndp(v, 0.0, 0., 1.0, 1.0, 0.0, 'Check bottom')
2666 test_ifndp(v, 0.1, 0., 1.1, 1.0, 0.0, 'Check below')
2668 test_ifndp(v,-1.0, 0., 0.0, 1.0, 1.0, 'Check top')
2670 test_ifndp(v,-1.0, 0., -0.1, 1.0, 1.0, 'Check above')
2672 test_ifndp(v,-1.0, 0., 3.0, 1.0, 0.25, 'Check 1/4')
2674 test_ifndp(v,-3.0, 0., 1.0, 1.0, 0.75, 'Check 3/4')
2676 test_ifndp(v, 1.0, 0., 1.0, 1.0, 0.0, 'Check dRho=0 below')
2678 test_ifndp(v,-1.0, 0., -1.0, 1.0, 1.0, 'Check dRho=0 above')
2680 test_ifndp(v, 0.0, 0., 0.0, 1.0, 0.5, 'Check dRho=0 mid')
2682 test_ifndp(v,-2.0, .5, 5.0, 0.5, 0.5, 'Check dP=0')
2683
2684 ! Identical columns
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 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 'Identical columns')
2699 absolute_positions(3, 8, (/0.,10.,20.,30./), kol, pilro), &
2700 (/0.,0.,10.,10.,20.,20.,30.,30./), '... left positions')
2702 absolute_positions(3, 8, (/0.,10.,20.,30./), kor, pirlo), &
2703 (/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 pilro, pirlo, kol, kor, heff, flx, .true., h_neglect)
2708 (/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 pilro, pirlo, kol, kor, heff, flx, .true., h_neglect)
2713 (/0.,20.,0.,20.,0.,20.,0./), 'Identical columns, S flux')
2714
2715 ! Right column slightly cooler than left
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 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 'Right column slightly cooler')
2730 absolute_positions(3, 8, (/0.,10.,20.,30./), kol, pilro), &
2731 (/0.,5.,10.,15.,20.,25.,30.,30./), '... left positions')
2733 absolute_positions(3, 8, (/0.,10.,20.,30./), kor, pirlo), &
2734 (/0.,0.,5.,10.,15.,20.,25.,30./), '... right positions')
2735
2736 ! Right column slightly warmer than left
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 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 'Right column slightly warmer')
2750
2751 ! Right column somewhat cooler than left
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 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 'Right column somewhat cooler')
2765
2766 ! Right column much colder than left with no overlap
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 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 'Right column much cooler')
2780
2781 ! Right column with mixed layer
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 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 'Right column with mixed layer')
2795
2796 ! Identical columns with mixed layer
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 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 'Identical columns with mixed layer')
2810
2811 ! Right column with unstable mixed layer
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 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 'Right column with unstable mixed layer')
2825
2826 ! Left column with unstable mixed layer
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 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 'Left column with unstable mixed layer')
2840
2841 ! Two unstable mixed layers
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 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 'Two unstable mixed layers')
2855
2856 if (.not. ndiff_unit_tests_continuous) write(stdout,*) 'Pass'
2857
2858end function ndiff_unit_tests_continuous
2859
2860logical 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 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 v = verbose
2884 ndiff_unit_tests_discontinuous = .false. ! Normally return false
2885 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 allocate(cs%EOS)
2890 call eos_manual_init(cs%EOS, form_of_eos=eos_linear, drho_dt=-1., drho_ds=0.)
2891 sl(:) = 0. ; sr(:) = 0. ; ; sil(:,:) = 0. ; sir(:,:) = 0.
2892 ppoly_t_l(:,:) = 0. ; ppoly_t_r(:,:) = 0.
2893 ppoly_s_l(:,:) = 0. ; ppoly_s_r(:,:) = 0.
2894 ! Intialize any control structures needed for unit tests
2895 cs%ref_pres = -1.
2896
2897 hl = (/10.,10.,10./) ; hr = (/10.,10.,10./)
2898 pres_l(1,1) = 0. ; pres_l(1,2) = hl(1) ; pres_r(1,1) = 0. ; pres_r(1,2) = hr(1)
2899 do k = 2,nk
2900 pres_l(k,1) = pres_l(k-1,2)
2901 pres_l(k,2) = pres_l(k,1) + hl(k)
2902 pres_r(k,1) = pres_r(k-1,2)
2903 pres_r(k,2) = pres_r(k,1) + hr(k)
2904 enddo
2905 cs%delta_rho_form = 'mid_pressure'
2906 cs%neutral_pos_method = 1
2907
2908 til(1,:) = (/ 22.00, 18.00 /) ; til(2,:) = (/ 18.00, 14.00 /) ; til(3,:) = (/ 14.00, 10.00 /)
2909 tir(1,:) = (/ 22.00, 18.00 /) ; tir(2,:) = (/ 18.00, 14.00 /) ; tir(3,:) = (/ 14.00, 10.00 /)
2910 call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
2911 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 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 'Identical Columns')
2921
2922 til(1,:) = (/ 22.00, 18.00 /) ; til(2,:) = (/ 18.00, 14.00 /) ; til(3,:) = (/ 14.00, 10.00 /)
2923 tir(1,:) = (/ 20.00, 16.00 /) ; tir(2,:) = (/ 16.00, 12.00 /) ; tir(3,:) = (/ 12.00, 8.00 /)
2924 call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
2925 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 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 'Right slightly cooler')
2935
2936 til(1,:) = (/ 20.00, 16.00 /) ; til(2,:) = (/ 16.00, 12.00 /) ; til(3,:) = (/ 12.00, 8.00 /)
2937 tir(1,:) = (/ 22.00, 18.00 /) ; tir(2,:) = (/ 18.00, 14.00 /) ; tir(3,:) = (/ 14.00, 10.00 /)
2938 call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
2939 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 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 'Left slightly cooler')
2949
2950 til(1,:) = (/ 22.00, 20.00 /) ; til(2,:) = (/ 18.00, 16.00 /) ; til(3,:) = (/ 14.00, 12.00 /)
2951 tir(1,:) = (/ 32.00, 24.00 /) ; tir(2,:) = (/ 22.00, 14.00 /) ; tir(3,:) = (/ 12.00, 4.00 /)
2952 call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
2953 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 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 'Right more strongly stratified')
2963
2964 til(1,:) = (/ 22.00, 18.00 /) ; til(2,:) = (/ 18.00, 14.00 /) ; til(3,:) = (/ 14.00, 10.00 /)
2965 tir(1,:) = (/ 14.00, 14.00 /) ; tir(2,:) = (/ 14.00, 14.00 /) ; tir(3,:) = (/ 12.00, 8.00 /)
2966 call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
2967 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 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 'Deep Mixed layer on the right')
2977
2978 til(1,:) = (/ 14.00, 14.00 /) ; til(2,:) = (/ 14.00, 12.00 /) ; til(3,:) = (/ 10.00, 8.00 /)
2979 tir(1,:) = (/ 14.00, 14.00 /) ; tir(2,:) = (/ 14.00, 14.00 /) ; tir(3,:) = (/ 14.00, 14.00 /)
2980 call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
2981 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 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 'Right unstratified column')
2991
2992 til(1,:) = (/ 14.00, 14.00 /) ; til(2,:) = (/ 14.00, 12.00 /) ; til(3,:) = (/ 10.00, 8.00 /)
2993 tir(1,:) = (/ 14.00, 14.00 /) ; tir(2,:) = (/ 14.00, 14.00 /) ; tir(3,:) = (/ 12.00, 4.00 /)
2994 call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
2995 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 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 'Right unstratified column')
3005
3006 til(1,:) = (/ 14.00, 14.00 /) ; til(2,:) = (/ 14.00, 10.00 /) ; til(3,:) = (/ 10.00, 2.00 /)
3007 tir(1,:) = (/ 14.00, 14.00 /) ; tir(2,:) = (/ 14.00, 10.00 /) ; tir(3,:) = (/ 10.00, 2.00 /)
3008 call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
3009 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 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 'Identical columns with mixed layer')
3019
3020 til(1,:) = (/ 14.00, 12.00 /) ; til(2,:) = (/ 10.00, 10.00 /) ; til(3,:) = (/ 8.00, 2.00 /)
3021 tir(1,:) = (/ 14.00, 12.00 /) ; tir(2,:) = (/ 12.00, 8.00 /) ; tir(3,:) = (/ 8.00, 2.00 /)
3022 call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
3023 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 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 'Left interior unstratified')
3033
3034 til(1,:) = (/ 12.00, 12.00 /) ; til(2,:) = (/ 12.00, 10.00 /) ; til(3,:) = (/ 10.00, 6.00 /)
3035 tir(1,:) = (/ 12.00, 10.00 /) ; tir(2,:) = (/ 10.00, 12.00 /) ; tir(3,:) = (/ 8.00, 4.00 /)
3036 call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
3037 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 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 'Left mixed layer, Right unstable interior')
3047
3048 til(1,:) = (/ 14.00, 14.00 /) ; til(2,:) = (/ 10.00, 10.00 /) ; til(3,:) = (/ 8.00, 6.00 /)
3049 tir(1,:) = (/ 10.00, 14.00 /) ; tir(2,:) = (/ 16.00, 16.00 /) ; tir(3,:) = (/ 12.00, 4.00 /)
3050 call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
3051 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 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 'Left thick mixed layer, Right unstable mixed')
3061
3062 til(1,:) = (/ 8.00, 12.00 /) ; til(2,:) = (/ 12.00, 10.00 /) ; til(3,:) = (/ 8.00, 4.00 /)
3063 tir(1,:) = (/ 10.00, 14.00 /) ; tir(2,:) = (/ 14.00, 12.00 /) ; tir(3,:) = (/ 10.00, 6.00 /)
3064 call mark_unstable_cells( cs, nk, til, sil, pres_l, stable_l )
3065 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 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 'Unstable mixed layers, left cooler')
3075
3076 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 cs%max_iter = 10
3080 ! Unit tests require explicit initialization of tolerance
3081 cs%Drho_tol = 0.
3082 cs%x_tol = 0.
3084 find_neutral_pos_linear(cs, 0., 10., 35., -0.2, 0., &
3085 -0.2, 0., -0.2, 0., &
3086 (/12.,-4./), (/34.,0./)), "Temp Uniform Linearized Alpha/Beta"))
3087 ! EOS linear in S, uniform beta
3089 find_neutral_pos_linear(cs, 0., 10., 35., 0., 0.8, &
3090 0., 0.8, 0., 0.8, &
3091 (/12.,0./), (/34.,2./)), "Salt Uniform Linearized Alpha/Beta"))
3092 ! EOS linear in T/S, uniform alpha/beta
3094 find_neutral_pos_linear(cs, 0., 10., 35., -0.5, 0.5, &
3095 -0.5, 0.5, -0.5, 0.5, &
3096 (/12.,-4./), (/34.,2./)), "Temp/salt Uniform Linearized Alpha/Beta"))
3097 ! EOS linear in T, insensitive to So
3099 find_neutral_pos_linear(cs, 0., 10., 35., -0.2, 0., &
3100 -0.4, 0., -0.6, 0., &
3101 (/12.,-4./), (/34.,0./)), "Temp stratified Linearized Alpha/Beta"))
3102 ! EOS linear in S, insensitive to T
3104 find_neutral_pos_linear(cs, 0., 10., 35., 0., 0.8, &
3105 0., 1.0, 0., 0.5, &
3106 (/12.,0./), (/34.,2./)), "Salt stratified Linearized Alpha/Beta"))
3107 if (.not. ndiff_unit_tests_discontinuous) write(stdout,*) 'Pass'
3108
3110
3111!> Returns true if a test of fv_diff() fails, and conditionally writes results to stream
3112logical 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 pret = fv_diff(hkm1, hk, hkp1, skm1, sk, skp1)
3128 test_fv_diff = (pret /= ptrue)
3129
3130 if (test_fv_diff .or. verbose) then
3131 stdunit = stdout
3132 if (test_fv_diff) stdunit = stderr ! In case of wrong results, write to error stream
3133 write(stdunit,'(a)') title
3134 if (test_fv_diff) then
3135 write(stdunit,'(2(1x,a,f20.16),1x,a)') 'pRet=',pret,'pTrue=',ptrue,'WRONG!'
3136 else
3137 write(stdunit,'(2(1x,a,f20.16))') 'pRet=',pret,'pTrue=',ptrue
3138 endif
3139 endif
3140
3141end function test_fv_diff
3142
3143!> Returns true if a test of fvlsq_slope() fails, and conditionally writes results to stream
3144logical 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 pret = fvlsq_slope(hkm1, hk, hkp1, skm1, sk, skp1)
3160 test_fvlsq_slope = (pret /= ptrue)
3161
3162 if (test_fvlsq_slope .or. verbose) then
3163 stdunit = stdout
3164 if (test_fvlsq_slope) stdunit = stderr ! In case of wrong results, write to error stream
3165 write(stdunit,'(a)') title
3166 if (test_fvlsq_slope) then
3167 write(stdunit,'(2(1x,a,f20.16),1x,a)') 'pRet=',pret,'pTrue=',ptrue,'WRONG!'
3168 else
3169 write(stdunit,'(2(1x,a,f20.16))') 'pRet=',pret,'pTrue=',ptrue
3170 endif
3171 endif
3172
3173end function test_fvlsq_slope
3174
3175!> Returns true if a test of interpolate_for_nondim_position() fails, and conditionally writes results to stream
3176logical 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 pret = interpolate_for_nondim_position(rhoneg, pneg, rhopos, ppos)
3190 test_ifndp = (pret /= ptrue)
3191
3192 if (test_ifndp .or. verbose) then
3193 stdunit = stdout
3194 if (test_ifndp) stdunit = stderr ! In case of wrong results, write to error stream
3195 write(stdunit,'(a)') title
3196 if (test_ifndp) then
3197 write(stdunit,'(4(1x,a,f20.16),2(1x,a,1pe22.15),1x,a)') &
3198 '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 'r1=',rhoneg,'p1=',pneg,'r2=',rhopos,'p2=',ppos,'pRet=',pret,'pTrue=',ptrue
3202 endif
3203 endif
3204
3205end function test_ifndp
3206
3207!> Returns true if comparison of Po and Ptrue fails, and conditionally writes results to stream
3208logical 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 test_data1d = .false.
3219 do k = 1,nk
3220 if (po(k) /= ptrue(k)) test_data1d = .true.
3221 enddo
3222
3223 if (test_data1d .or. verbose) then
3224 stdunit = stdout
3225 if (test_data1d) stdunit = stderr ! In case of wrong results, write to error stream
3226 write(stdunit,'(a)') title
3227 do k = 1,nk
3228 if (po(k) /= ptrue(k)) then
3229 test_data1d = .true.
3230 write(stdunit,'(a,I0,2(1x,a,f20.16),1x,a,1pe22.15,1x,a)') &
3231 'k=',k,'Po=',po(k),'Ptrue=',ptrue(k),'err=',po(k)-ptrue(k),'WRONG!'
3232 else
3233 if (verbose) &
3234 write(stdunit,'(a,I0,2(1x,a,f20.16),1x,a,1pe22.15)') &
3235 'k=',k,'Po=',po(k),'Ptrue=',ptrue(k),'err=',po(k)-ptrue(k)
3236 endif
3237 enddo
3238 endif
3239
3240end function test_data1d
3241
3242!> Returns true if comparison of Po and Ptrue fails, and conditionally writes results to stream
3243logical 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 test_data1di = .false.
3254 do k = 1,nk
3255 if (po(k) /= ptrue(k)) test_data1di = .true.
3256 enddo
3257
3258 if (test_data1di .or. verbose) then
3259 stdunit = stdout
3260 if (test_data1di) stdunit = stderr ! In case of wrong results, write to error stream
3261 write(stdunit,'(a)') title
3262 do k = 1,nk
3263 if (po(k) /= ptrue(k)) then
3264 test_data1di = .true.
3265 write(stdunit,'(a,I0,2(1x,a,i5),1x,a)') 'k=',k,'Io=',po(k),'Itrue=',ptrue(k),'WRONG!'
3266 else
3267 if (verbose) &
3268 write(stdunit,'(a,I0,2(1x,a,i5))') 'k=',k,'Io=',po(k),'Itrue=',ptrue(k)
3269 endif
3270 enddo
3271 endif
3272
3273end 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
3277logical 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 test_nsp = .false.
3299 do k = 1,ns
3300 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 if (k < ns) then
3302 if (heff(k) /= heff0(k)) test_nsp = .true.
3303 endif
3304 enddo
3305
3306 if (test_nsp .or. verbose) then
3307 stdunit = stdout
3308 if (test_nsp) stdunit = stderr ! In case of wrong results, write to error stream
3309 write(stdunit,'(a)') title
3310 do k = 1,ns
3311 this_row_failed = compare_nsp_row(kol(k), kor(k), pl(k), pr(k), kol0(k), kor0(k), pl0(k), pr0(k))
3312 if (this_row_failed) then
3313 write(stdunit,10) k,kol(k),pl(k),kor(k),pr(k),' <-- WRONG!'
3314 write(stdunit,10) k,kol0(k),pl0(k),kor0(k),pr0(k),' <-- should be this'
3315 else
3316 write(stdunit,10) k,kol(k),pl(k),kor(k),pr(k)
3317 endif
3318 if (k < ns) then
3319 if (heff(k) /= heff0(k)) then
3320 write(stdunit,'(i3,8x,"layer hEff =",2(f20.16,a))') k,heff(k)," .neq. ",heff0(k),' <-- WRONG!'
3321 else
3322 write(stdunit,'(i3,8x,"layer hEff =",f20.16)') k,heff(k)
3323 endif
3324 endif
3325 enddo
3326 endif
3327 if (test_nsp) call mom_error(fatal,"test_nsp failed")
3328
332910 format("ks=",i3," kL=",i3," pL=",f20.16," kR=",i3," pR=",f20.16,a)
3330end function test_nsp
3331
3332!> Compares a single row, k, of output from find_neutral_surface_positions()
3333logical 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 compare_nsp_row = .false.
3344 if (kol /= kol0) compare_nsp_row = .true.
3345 if (kor /= kor0) compare_nsp_row = .true.
3346 if (pl /= pl0) compare_nsp_row = .true.
3347 if (pr /= pr0) compare_nsp_row = .true.
3348end function compare_nsp_row
3349
3350!> Compares output position from refine_nondim_position with an expected value
3351logical 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 stdunit = stdout ! Output to standard error
3359 test_rnp = abs(expected_pos - test_pos) > 2*epsilon(test_pos)
3360 if (test_rnp) then
3361 write(stdunit,'(A, f20.16, " .neq. ", f20.16, " <-- WRONG")') title, expected_pos, test_pos
3362 else
3363 write(stdunit,'(A, f20.16, " == ", f20.16)') title, expected_pos, test_pos
3364 endif
3365end function test_rnp
3366
3367!> Deallocates neutral_diffusion control structure
3368subroutine neutral_diffusion_end(CS)
3369 type(neutral_diffusion_cs), pointer :: cs !< Neutral diffusion control structure
3370
3371 if (associated(cs)) deallocate(cs)
3372
3373end subroutine neutral_diffusion_end
3374
3375end module mom_neutral_diffusion