src/parameterizations/vertical/MOM_CVMix_KPP.F90

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

1		! This file is part of MOM6, the Modular Ocean Model version 6.
2		! See the LICENSE file for licensing information.
3		! SPDX-License-Identifier: Apache-2.0
4
5		!> Provides the K-Profile Parameterization (KPP) of Large et al., 1994, via CVMix.
6		module MOM_CVMix_KPP
7
8		use MOM_coms, only : max_across_PEs
9		use MOM_debugging, only : hchksum, is_NaN
10		use MOM_diag_mediator, only : time_type, diag_ctrl, safe_alloc_ptr, post_data
11		use MOM_diag_mediator, only : query_averaging_enabled, register_diag_field
12		use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING, is_root_PE
13		use MOM_EOS, only : EOS_type, calculate_density
14		use MOM_file_parser, only : get_param, log_param, log_version, param_file_type
15		use MOM_file_parser, only : openParameterBlock, closeParameterBlock
16		use MOM_grid, only : ocean_grid_type, isPointInCell
17		use MOM_interface_heights, only : thickness_to_dz
18		use MOM_restart, only : MOM_restart_CS, register_restart_field
19		use MOM_unit_scaling, only : unit_scale_type
20		use MOM_variables, only : thermo_var_ptrs
21		use MOM_verticalGrid, only : verticalGrid_type
22		use MOM_wave_interface, only : wave_parameters_CS, Get_Langmuir_Number, get_wave_method
23		use MOM_domains, only : pass_var
24		use MOM_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
25		use MOM_cpu_clock, only : CLOCK_MODULE, CLOCK_ROUTINE
26		use MOM_tracer_types, only : tracer_type
27
28		use CVMix_kpp, only : CVMix_init_kpp, CVMix_put_kpp, CVMix_get_kpp_real
29		use CVMix_kpp, only : CVMix_coeffs_kpp
30		use CVMix_kpp, only : CVMix_kpp_compute_OBL_depth
31		use CVMix_kpp, only : CVMix_kpp_compute_turbulent_scales
32		use CVMix_kpp, only : CVMix_kpp_compute_bulk_Richardson
33		use CVMix_kpp, only : CVMix_kpp_compute_unresolved_shear
34		use CVMix_kpp, only : CVMix_kpp_params_type
35		use CVMix_kpp, only : CVMix_kpp_compute_kOBL_depth
36		use CVMix_kpp, only : CVMix_kpp_compute_StokesXi
37
38		implicit none ; private
39
40		#include "MOM_memory.h"
41
42		public :: register_KPP_restarts
43		public :: KPP_init
44		public :: KPP_compute_BLD
45		public :: KPP_calculate
46		public :: KPP_end
47		public :: KPP_NonLocalTransport_temp
48		public :: KPP_NonLocalTransport_saln
49		public :: KPP_NonLocalTransport
50		public :: KPP_get_BLD
51
52		! Enumerated constants
53		integer, private, parameter :: NLT_SHAPE_CVMix = 0 !< Use the CVMix profile
54		integer, private, parameter :: NLT_SHAPE_LINEAR = 1 !< Linear, \f$ G(\sigma) = 1-\sigma \f$
55		integer, private, parameter :: NLT_SHAPE_PARABOLIC = 2 !< Parabolic, \f$ G(\sigma) = (1-\sigma)^2 \f$
56		integer, private, parameter :: NLT_SHAPE_CUBIC = 3 !< Cubic, \f$ G(\sigma) = 1 + (2\sigma-3) \sigma^2\f$
57		integer, private, parameter :: NLT_SHAPE_CUBIC_LMD = 4 !< Original shape,
58		!! \f$ G(\sigma) = \frac{27}{4} \sigma (1-\sigma)^2 \f$
59
60		integer, private, parameter :: SW_METHOD_ALL_SW = 0 !< Use all shortwave radiation
61		integer, private, parameter :: SW_METHOD_MXL_SW = 1 !< Use shortwave radiation absorbed in mixing layer
62		integer, private, parameter :: SW_METHOD_LV1_SW = 2 !< Use shortwave radiation absorbed in layer 1
63		integer, private, parameter :: LT_K_CONSTANT = 1, & !< Constant enhance K through column
64		LT_K_SCALED = 2, & !< Enhance K scales with G(sigma)
65		LT_K_MODE_CONSTANT = 1, & !< Prescribed enhancement for K
66		LT_K_MODE_VR12 = 2, & !< Enhancement for K based on
67		!! Van Roekel et al., 2012
68		LT_K_MODE_RW16 = 3, & !< Enhancement for K based on
69		!! Reichl et al., 2016
70		LT_VT2_MODE_CONSTANT = 1, & !< Prescribed enhancement for Vt2
71		LT_VT2_MODE_VR12 = 2, & !< Enhancement for Vt2 based on
72		!! Van Roekel et al., 2012
73		LT_VT2_MODE_RW16 = 3, & !< Enhancement for Vt2 based on
74		!! Reichl et al., 2016
75		LT_VT2_MODE_LF17 = 4 !< Enhancement for Vt2 based on
76		!! Li and Fox-Kemper, 2017
77
78		!> Control structure for containing KPP parameters/data
79		type, public :: KPP_CS ; private
80
81		! Parameters
82		real :: Ri_crit !< Critical bulk Richardson number (defines OBL depth) [nondim]
83		real :: vonKarman !< von Karman constant (dimensionless) [nondim]
84		real :: cs !< Parameter for computing velocity scale function (dimensionless) [nondim]
85		real :: cs2 !< Parameter for multiplying by non-local term [nondim]
86		! This is active for NLT_SHAPE_CUBIC_LMD only
87		logical :: enhance_diffusion !< If True, add enhanced diffusivity at base of boundary layer.
88		character(len=32) :: interpType !< Type of interpolation to compute bulk Richardson number
89		character(len=32) :: interpType2 !< Type of interpolation to compute diff and visc at OBL_depth
90		logical :: StokesMOST !< If True, use Stokes similarity package
91		logical :: computeEkman !< If True, compute Ekman depth limit for OBLdepth
92		logical :: computeMoninObukhov !< If True, compute Monin-Obukhov limit for OBLdepth
93		logical :: passiveMode !< If True, makes KPP passive meaning it does NOT alter the diffusivity
94		real :: deepOBLoffset !< If non-zero, is a distance from the bottom that the OBL can not
95		!! penetrate through [Z ~> m]
96		real :: minOBLdepth !< If non-zero, is a minimum depth for the OBL [Z ~> m]
97		real :: surf_layer_ext !< Fraction of OBL depth considered in the surface layer [nondim]
98		real :: minVtsqr !< Min for the squared unresolved velocity used in Rib CVMix
99		!! calculation [L2 T-2 ~> m2 s-2]
100		logical :: fixedOBLdepth !< If True, will fix the OBL depth at fixedOBLdepth_value
101		real :: fixedOBLdepth_value !< value for the fixed OBL depth when fixedOBLdepth==True [Z ~> m]
102		logical :: debug !< If True, calculate checksums and write debugging information
103		character(len=30) :: MatchTechnique !< Method used in CVMix for setting diffusivity and NLT profile functions
104		integer :: NLT_shape !< MOM6 over-ride of CVMix NLT shape function
105		logical :: applyNonLocalTrans !< If True, apply non-local transport to all tracers
106		integer :: n_smooth !< Number of times smoothing operator is applied on OBLdepth.
107		logical :: deepen_only !< If true, apply OBLdepth smoothing at a cell only if the OBLdepth gets deeper.
108		logical :: KPPzeroDiffusivity !< If True, will set diffusivity and viscosity from KPP to zero
109		!! for testing purposes.
110		logical :: KPPisAdditive !< If True, will add KPP diffusivity to initial diffusivity.
111		!! If False, will replace initial diffusivity wherever KPP diffusivity
112		!! is non-zero.
113		real :: min_thickness !< A minimum thickness used to avoid division by small numbers
114		!! in the vicinity of vanished layers [Z ~> m]
115		integer :: SW_METHOD !< Sets method for using shortwave radiation in surface buoyancy flux
116		logical :: LT_K_Enhancement !< Flags if enhancing mixing coefficients due to LT
117		integer :: LT_K_Shape !< Integer for constant or shape function enhancement
118		integer :: LT_K_Method !< Integer for mixing coefficients LT method
119		real :: KPP_CVt2 !< Parameter for Stokes MOST convection entrainment [nondim]
120		real :: KPP_K_ENH_FAC !< Factor to multiply by K if Method is CONSTANT [nondim]
121		logical :: LT_Vt2_Enhancement !< Flags if enhancing Vt2 due to LT
122		integer :: LT_VT2_METHOD !< Integer for Vt2 LT method
123		real :: KPP_VT2_ENH_FAC !< Factor to multiply by VT2 if Method is CONSTANT [nondim]
124		real :: MLD_guess_min !< The minimum estimate of the mixed layer depth used to
125		!! calculate the Langmuir number for Langmuir turbulence
126		!! enhancement with KPP [Z ~> m]
127		logical :: STOKES_MIXING !< Flag if model is mixing down Stokes gradient
128		!! This is relevant for which current to use in RiB
129		integer :: answer_date !< The vintage of the order of arithmetic in the CVMix KPP
130		!! calculations. Values below 20240501 recover the answers
131		!! from early in 2024, while higher values use expressions
132		!! that have been refactored for rotational symmetry.
133
134		!> CVMix parameters
135		type(CVMix_kpp_params_type), pointer :: KPP_params => NULL()
136
137		type(diag_ctrl), pointer :: diag => NULL() !< Pointer to diagnostics control structure
138		!>@{ Diagnostic handles
139		integer :: id_OBLdepth = -1, id_BulkRi = -1
140		integer :: id_N = -1, id_N2 = -1
141		integer :: id_Ws = -1, id_Vt2 = -1
142		integer :: id_BulkUz2 = -1, id_BulkDrho = -1
143		integer :: id_uStar = -1, id_buoyFlux = -1
144		integer :: id_sigma = -1, id_Kv_KPP = -1
145		integer :: id_Kt_KPP = -1, id_Ks_KPP = -1
146		integer :: id_Tsurf = -1, id_Ssurf = -1
147		integer :: id_Usurf = -1, id_Vsurf = -1
148		integer :: id_Kd_in = -1
149		integer :: id_NLTt = -1
150		integer :: id_NLTs = -1
151		integer :: id_EnhK = -1, id_EnhVt2 = -1
152		integer :: id_EnhW = -1
153		integer :: id_La_SL = -1
154		integer :: id_OBLdepth_original = -1
155		integer :: id_StokesXI = -1
156		integer :: id_Lam2 = -1
157		!>@}
158
159		! Diagnostics arrays
160		real, pointer, dimension(:,:) :: OBLdepth !< Depth (positive) of ocean boundary layer (OBL) [Z ~> m]
161		real, allocatable, dimension(:,:) :: OBLdepth_original !< Depth (positive) of OBL without smoothing [Z ~> m]
162		real, allocatable, dimension(:,:) :: StokesParXI !< Stokes similarity parameter [nondim]
163		real, allocatable, dimension(:,:) :: Lam2 !< La^(-2) = Ustk0/u* [nondim]
164		real, allocatable, dimension(:,:) :: kOBL !< Level (+fraction) of OBL extent [nondim]
165		real, allocatable, dimension(:,:) :: OBLdepthprev !< previous Depth (positive) of OBL [Z ~> m]
166		real, allocatable, dimension(:,:) :: La_SL !< Langmuir number used in KPP [nondim]
167		real, allocatable, dimension(:,:,:) :: dRho !< Bulk difference in density [R ~> kg m-3]
168		real, allocatable, dimension(:,:,:) :: Uz2 !< Square of bulk difference in resolved velocity [L2 T-2 ~> m2 s-2]
169		real, allocatable, dimension(:,:,:) :: BulkRi !< Bulk Richardson number for each layer [nondim]
170		real, allocatable, dimension(:,:,:) :: sigma !< Sigma coordinate (dimensionless) [nondim]
171		real, allocatable, dimension(:,:,:) :: Ws !< Turbulent velocity scale for scalars [Z T-1 ~> m s-1]
172		real, allocatable, dimension(:,:,:) :: N !< Brunt-Vaisala frequency [T-1 ~> s-1]
173		real, allocatable, dimension(:,:,:) :: N2 !< Squared Brunt-Vaisala frequency [T-2 ~> s-2]
174		real, allocatable, dimension(:,:,:) :: Vt2 !< Unresolved squared turbulence velocity for
175		!! bulk Ri [Z2 T-2 ~> m2 s-2]
176		real, allocatable, dimension(:,:,:) :: Kt_KPP !< Temp diffusivity from KPP [Z2 T-1 ~> m2 s-1]
177		real, allocatable, dimension(:,:,:) :: Ks_KPP !< Scalar diffusivity from KPP [Z2 T-1 ~> m2 s-1]
178		real, allocatable, dimension(:,:,:) :: Kv_KPP !< Viscosity due to KPP [Z2 T-1 ~> m2 s-1]
179		real, allocatable, dimension(:,:) :: Tsurf !< Temperature of surface layer [C ~> degC]
180		real, allocatable, dimension(:,:) :: Ssurf !< Salinity of surface layer [S ~> ppt]
181		real, allocatable, dimension(:,:) :: Usurf !< i-velocity of surface layer [L T-1 ~> m s-1]
182		real, allocatable, dimension(:,:) :: Vsurf !< j-velocity of surface layer [L T-1 ~> m s-1]
183		real, allocatable, dimension(:,:,:) :: EnhK !< Enhancement for mixing coefficient [nondim]
184		real, allocatable, dimension(:,:,:) :: EnhVt2 !< Enhancement for Vt2 [nondim]
185
186		end type KPP_CS
187
188		!>@{ CPU time clocks
189		integer :: id_clock_KPP_calc, id_clock_KPP_compute_BLD, id_clock_KPP_smoothing
190		!>@}
191
192		#define __DO_SAFETY_CHECKS__
193
194		contains
195
196		!> Routine to register restarts, pass-through to children modules
197	1	subroutine register_KPP_restarts(G, param_file, restart_CSp, CS)
198		type(ocean_grid_type), intent(in) :: G !< The ocean's grid structure
199		type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
200		type(MOM_restart_CS), pointer :: restart_CSp !< MOM restart control structure
201		type(KPP_CS), pointer :: CS !< module control structure
202
203		character(len=40) :: mdl = 'MOM_CVMix_KPP' !< name of this module
204		logical :: use_kpp, fpmix
205
206	1	if (associated(CS)) call MOM_error(FATAL, 'MOM_CVMix_KPP, register_KPP_restarts: '// &
207	0	'Control structure has already been initialized')
208	1	call get_param(param_file, mdl, "USE_KPP", use_kpp, default=.false., do_not_log=.true.)
209		! Forego remainder of initialization if not using this scheme
210	1	if (.not. use_kpp) return
211	0	allocate(CS)
212
213	0	allocate(CS%OBLdepth(SZI_(G),SZJ_(G)), source=0.0)
214
215		! FPMIX is needed to decide if boundary layer depth should be added to restart file
216		call get_param(param_file, '', "FPMIX", fpmix, &
217		"If true, add non-local momentum flux increments and diffuse down the Eulerian gradient.", &
218	0	default=.false., do_not_log=.true.)
219	0	if (fpmix) call register_restart_field(CS%OBLdepth, 'KPP_OBLdepth', .false., restart_CSp)
220
221		end subroutine register_KPP_restarts
222
223		!> Initialize the CVMix KPP module and set up diagnostics
224		!! Returns True if KPP is to be used, False otherwise.
225	1	logical function KPP_init(paramFile, G, GV, US, diag, Time, CS, passive)
226
227		! Arguments
228		type(param_file_type), intent(in) :: paramFile !< File parser
229		type(ocean_grid_type), intent(in) :: G !< Ocean grid
230		type(verticalGrid_type), intent(in) :: GV !< Vertical grid structure
231		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
232		type(diag_ctrl), target, intent(in) :: diag !< Diagnostics
233		type(time_type), intent(in) :: Time !< Model time
234		type(KPP_CS), pointer :: CS !< Control structure
235		logical, optional, intent(out) :: passive !< Copy of %passiveMode
236
237		! Local variables
238		# include "version_variable.h"
239		character(len=40) :: mdl = 'MOM_CVMix_KPP' !< name of this module
240		character(len=20) :: string !< local temporary string
241		character(len=20) :: langmuir_mixing_opt = 'NONE' !< Langmuir mixing option to be passed to CVMix, e.g., LWF16
242		character(len=20) :: langmuir_entrainment_opt = 'NONE' !< Langmuir entrainment option to be
243		!! passed to CVMix, e.g., LWF16
244		integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.
245		logical :: CS_IS_ONE=.false. !< Logical for setting Cs based on Non-local
246		logical :: lnoDGat1=.false. !< True => G'(1) = 0 (shape function)
247		!! False => compute G'(1) as in LMD94
248		! Read parameters
249	1	call get_param(paramFile, mdl, "USE_KPP", KPP_init, default=.false., do_not_log=.true.)
250		call log_version(paramFile, mdl, version, 'This is the MOM wrapper to CVMix:KPP\n' // &
251	1	'See http://cvmix.github.io/', all_default=.not.KPP_init)
252		call get_param(paramFile, mdl, "USE_KPP", KPP_init, &
253		"If true, turns on the [CVMix] KPP scheme of Large et al., 1994, "// &
254		"to calculate diffusivities and non-local transport in the OBL.", &
255	1	default=.false.)
256		! Forego remainder of initialization if not using this scheme
257	1	if (.not. KPP_init) return
258
259		call get_param(paramFile, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
260		"This sets the default value for the various _ANSWER_DATE parameters.", &
261	0	default=99991231, do_not_log=.true.)
262
263	0	call openParameterBlock(paramFile,'KPP')
264		call get_param(paramFile, mdl, 'PASSIVE', CS%passiveMode, &
265		'If True, puts KPP into a passive-diagnostic mode.', &
266	0	default=.False.)
267		!BGR: Note using PASSIVE for KPP creates warning for PASSIVE from Convection
268		! should we create a separate flag?
269	0	if (present(passive)) passive=CS%passiveMode ! This is passed back to the caller so
270		! the caller knows to not use KPP output
271		call get_param(paramFile, mdl, 'APPLY_NONLOCAL_TRANSPORT', CS%applyNonLocalTrans, &
272		'If True, applies the non-local transport to all tracers. '// &
273		'If False, calculates the non-local transport and tendencies but '//&
274		'purely for diagnostic purposes.', &
275	0	default=.not. CS%passiveMode)
276		call get_param(paramFile, mdl, 'N_SMOOTH', CS%n_smooth, &
277		'The number of times the 1-1-4-1-1 Laplacian filter is applied on OBL depth.', &
278	0	default=0)
279	0	if (CS%n_smooth > G%domain%nihalo) then
280	0	call MOM_error(FATAL,'KPP smoothing number (N_SMOOTH) cannot be greater than NIHALO.')
281	0	elseif (CS%n_smooth > G%domain%njhalo) then
282	0	call MOM_error(FATAL,'KPP smoothing number (N_SMOOTH) cannot be greater than NJHALO.')
283		endif
284	0	if (CS%n_smooth > 0) then
285		call get_param(paramFile, mdl, 'DEEPEN_ONLY_VIA_SMOOTHING', CS%deepen_only, &
286		'If true, apply OBLdepth smoothing at a cell only if the OBLdepth '// &
287		'gets deeper via smoothing.', &
288	0	default=.false.)
289	0	id_clock_KPP_smoothing = cpu_clock_id('(Ocean KPP BLD smoothing)', grain=CLOCK_ROUTINE)
290		endif
291		call get_param(paramFile, mdl, 'RI_CRIT', CS%Ri_crit, &
292		'Critical bulk Richardson number used to define depth of the '// &
293		'surface Ocean Boundary Layer (OBL).', &
294	0	units='nondim', default=0.3)
295		call get_param(paramFile, mdl, 'VON_KARMAN', CS%vonKarman, &
296		'von Karman constant.', &
297	0	units='nondim', default=0.40)
298		call get_param(paramFile, mdl, 'ENHANCE_DIFFUSION', CS%enhance_diffusion, &
299		'If True, adds enhanced diffusion at the based of the boundary layer.', &
300	0	default=.true.)
301		call get_param(paramFile, mdl, 'INTERP_TYPE', CS%interpType, &
302		'Type of interpolation to determine the OBL depth.\n'// &
303		'Allowed types are: linear, quadratic, cubic.', &
304	0	default='quadratic')
305		call get_param(paramFile, mdl, 'INTERP_TYPE2', CS%interpType2, &
306		'Type of interpolation to compute diff and visc at OBL_depth.\n'// &
307		'Allowed types are: linear, quadratic, cubic or LMD94.', &
308	0	default='LMD94')
309		call get_param(paramFile, mdl, 'STOKES_MOST', CS%StokesMOST, &
310		'If True, use Stokes Similarity package.', &
311	0	default=.False.)
312		call get_param(paramFile, mdl, 'COMPUTE_EKMAN', CS%computeEkman, &
313		'If True, limit OBL depth to be no deeper than Ekman depth.', &
314	0	default=.False.)
315		call get_param(paramFile, mdl, 'COMPUTE_MONIN_OBUKHOV', CS%computeMoninObukhov, &
316		'If True, limit the OBL depth to be no deeper than '// &
317		'Monin-Obukhov depth.', &
318	0	default=.False.)
319		call get_param(paramFile, mdl, 'CS', CS%cs, &
320		'Parameter for computing velocity scale function.', &
321	0	units='nondim', default=98.96)
322		call get_param(paramFile, mdl, 'CS2', CS%cs2, &
323		'Parameter for computing non-local term.', &
324	0	units='nondim', default=6.32739901508)
325		call get_param(paramFile, mdl, 'DEEP_OBL_OFFSET', CS%deepOBLoffset, &
326		'If non-zero, the distance above the bottom to which the OBL is clipped '// &
327		'if it would otherwise reach the bottom. The smaller of this and 0.1D is used.', &
328	0	units='m', default=0., scale=US%m_to_Z)
329		call get_param(paramFile, mdl, 'FIXED_OBLDEPTH', CS%fixedOBLdepth, &
330		'If True, fix the OBL depth to FIXED_OBLDEPTH_VALUE '// &
331		'rather than using the OBL depth from CVMix. '// &
332		'This option is just for testing purposes.', &
333	0	default=.False.)
334		call get_param(paramFile, mdl, 'FIXED_OBLDEPTH_VALUE', CS%fixedOBLdepth_value, &
335		'Value for the fixed OBL depth when fixedOBLdepth==True. '// &
336		'This parameter is for just for testing purposes. '// &
337		'It will over-ride the OBLdepth computed from CVMix.', &
338	0	units='m', default=30.0, scale=US%m_to_Z)
339		call get_param(paramFile, mdl, 'SURF_LAYER_EXTENT', CS%surf_layer_ext, &
340		'Fraction of OBL depth considered in the surface layer.', &
341	0	units='nondim', default=0.10)
342		call get_param(paramFile, mdl, 'MINIMUM_OBL_DEPTH', CS%minOBLdepth, &
343		'If non-zero, a minimum depth to use for KPP OBL depth. Independent of '// &
344		'this parameter, the OBL depth is always at least as deep as the first layer.', &
345	0	units='m', default=0., scale=US%m_to_Z)
346		call get_param(paramFile, mdl, 'MINIMUM_VT2', CS%minVtsqr, &
347		'Min of the unresolved velocity Vt2 used in Rib CVMix calculation.\n'// &
348		'Scaling: MINIMUM_VT2 = const1dN*ws, with d=1m, N=1e-5/s, ws=1e-6 m/s.', &
349	0	units='m2/s2', default=1e-10, scale=US%m_s_to_L_T**2)
350
351		call get_param(paramFile, mdl, 'NLT_SHAPE', string, &
352		'MOM6 method to set nonlocal transport profile. '// &
353		'Over-rides the result from CVMix. Allowed values are: \n'// &
354		'\t CVMix - Uses the profiles from CVMix specified by MATCH_TECHNIQUE\n'//&
355		'\t LINEAR - A linear profile, 1-sigma\n'// &
356		'\t PARABOLIC - A parablic profile, (1-sigma)^2\n'// &
357		'\t CUBIC - A cubic profile, (1-sigma)^2(1+2*sigma)\n'// &
358		'\t CUBIC_LMD - The original KPP profile', &
359	0	default='CVMix')
360	0	select case ( trim(string) )
361	0	case ("CVMix") ; CS%NLT_shape = NLT_SHAPE_CVMix
362	0	case ("LINEAR") ; CS%NLT_shape = NLT_SHAPE_LINEAR
363	0	case ("PARABOLIC") ; CS%NLT_shape = NLT_SHAPE_PARABOLIC
364	0	case ("CUBIC") ; CS%NLT_shape = NLT_SHAPE_CUBIC
365	0	case ("CUBIC_LMD") ; CS%NLT_shape = NLT_SHAPE_CUBIC_LMD
366		case default ; call MOM_error(FATAL,"KPP_init: "// &
367	0	"Unrecognized NLT_SHAPE option"//trim(string))
368		end select
369		call get_param(paramFile, mdl, 'MATCH_TECHNIQUE', CS%MatchTechnique, &
370		'CVMix method to set profile function for diffusivity and NLT, '// &
371		'as well as matching across OBL base. Allowed values are: \n'// &
372		'\t SimpleShapes = sigma*(1-sigma)^2 for both diffusivity and NLT\n'// &
373		'\t MatchGradient = sigma*(1-sigma)^2 for NLT; diffusivity profile from matching\n'//&
374		'\t MatchBoth = match gradient for both diffusivity and NLT\n'// &
375		'\t ParabolicNonLocal = sigma*(1-sigma)^2 for diffusivity; (1-sigma)^2 for NLT', &
376	0	default='SimpleShapes')
377	0	if (CS%MatchTechnique == 'ParabolicNonLocal') then
378		! This forces Cs2 (Cs in non-local computation) to equal 1 for parabolic non-local option.
379		! May be used during CVMix initialization.
380	0	Cs_is_one=.true.
381		endif
382	0	if (CS%MatchTechnique == 'ParabolicNonLocal' .or. CS%MatchTechnique == 'SimpleShapes') then
383		! if gradient won't be matched, lnoDGat1=.true.
384	0	lnoDGat1=.true.
385		endif
386
387		! safety check to avoid negative diff/visc
388	0	if (CS%MatchTechnique == 'MatchBoth' .and. (CS%interpType2 == 'cubic' .or. &
389		CS%interpType2 == 'quadratic')) then
390		call MOM_error(FATAL,"If MATCH_TECHNIQUE=MatchBoth, INTERP_TYPE2 must be set to \n"//&
391		"linear or LMD94 (recommended) to avoid negative viscosity and diffusivity.\n"//&
392	0	"Please select one of these valid options." )
393		endif
394
395		call get_param(paramFile, mdl, 'KPP_ZERO_DIFFUSIVITY', CS%KPPzeroDiffusivity, &
396		'If True, zeroes the KPP diffusivity and viscosity; for testing purpose.',&
397	0	default=.False.)
398		call get_param(paramFile, mdl, 'KPP_IS_ADDITIVE', CS%KPPisAdditive, &
399		'If true, adds KPP diffusivity to diffusivity from other schemes.\n'//&
400		'If false, KPP is the only diffusivity wherever KPP is non-zero.', &
401	0	default=.True.)
402		call get_param(paramFile, mdl, 'KPP_SHORTWAVE_METHOD',string, &
403		'Determines contribution of shortwave radiation to KPP surface '// &
404		'buoyancy flux. Options include:\n'// &
405		' ALL_SW: use total shortwave radiation\n'// &
406		' MXL_SW: use shortwave radiation absorbed by mixing layer\n'// &
407		' LV1_SW: use shortwave radiation absorbed by top model layer', &
408	0	default='MXL_SW')
409	0	select case ( trim(string) )
410	0	case ("ALL_SW") ; CS%SW_METHOD = SW_METHOD_ALL_SW
411	0	case ("MXL_SW") ; CS%SW_METHOD = SW_METHOD_MXL_SW
412	0	case ("LV1_SW") ; CS%SW_METHOD = SW_METHOD_LV1_SW
413		case default ; call MOM_error(FATAL,"KPP_init: "// &
414	0	"Unrecognized KPP_SHORTWAVE_METHOD option"//trim(string))
415		end select
416		call get_param(paramFile, mdl, 'CVMix_ZERO_H_WORK_AROUND', CS%min_thickness, &
417		'A minimum thickness used to avoid division by small numbers in the vicinity '// &
418		'of vanished layers. This is independent of MIN_THICKNESS used in other parts of MOM.', &
419	0	units='m', default=0., scale=US%m_to_Z)
420
421		!/BGR: New options for including Langmuir effects
422		!/ 1. Options related to enhancing the mixing coefficient
423		call get_param(paramFile, mdl, "USE_KPP_LT_K", CS%LT_K_Enhancement, &
424		'Flag for Langmuir turbulence enhancement of turbulent '//&
425	0	'mixing coefficient.', Default=.false.)
426		call get_param(paramFile, mdl, "STOKES_MIXING", CS%Stokes_Mixing, &
427		'Flag for Langmuir turbulence enhancement of turbulent '//&
428	0	'mixing coefficient.', Default=.false.)
429	0	if (CS%LT_K_Enhancement) then
430		call get_param(paramFile, mdl, 'KPP_LT_K_SHAPE', string, &
431		'Vertical dependence of LT enhancement of mixing. '// &
432		'Valid options are: \n'// &
433		'\t CONSTANT = Constant value for full OBL\n'// &
434		'\t SCALED = Varies based on normalized shape function.', &
435	0	default='CONSTANT')
436	0	select case ( trim(string))
437	0	case ("CONSTANT") ; CS%LT_K_SHAPE = LT_K_CONSTANT
438	0	case ("SCALED") ; CS%LT_K_SHAPE = LT_K_SCALED
439		case default ; call MOM_error(FATAL,"KPP_init: "//&
440	0	"Unrecognized KPP_LT_K_SHAPE option: "//trim(string))
441		end select
442		call get_param(paramFile, mdl, "KPP_LT_K_METHOD", string , &
443		'Method to enhance mixing coefficient in KPP. '// &
444		'Valid options are: \n'// &
445		'\t CONSTANT = Constant value (KPP_K_ENH_FAC) \n'// &
446		'\t VR12 = Function of Langmuir number based on VR12\n'// &
447		'\t (Van Roekel et al. 2012)\n'// &
448		'\t (Li et al. 2016, OM) \n'// &
449		'\t RW16 = Function of Langmuir number based on RW16\n'// &
450		'\t (Reichl et al., 2016, JPO)', &
451	0	default='CONSTANT')
452	0	select case ( trim(string))
453		case ("CONSTANT")
454	0	CS%LT_K_METHOD = LT_K_MODE_CONSTANT
455	0	langmuir_mixing_opt = 'LWF16'
456		case ("VR12")
457	0	CS%LT_K_METHOD = LT_K_MODE_VR12
458	0	langmuir_mixing_opt = 'LWF16'
459		case ("RW16")
460	0	CS%LT_K_METHOD = LT_K_MODE_RW16
461	0	langmuir_mixing_opt = 'RWHGK16'
462		case default
463		call MOM_error(FATAL,"KPP_init: "//&
464	0	"Unrecognized KPP_LT_K_METHOD option: "//trim(string))
465		end select
466	0	if (CS%LT_K_METHOD==LT_K_MODE_CONSTANT) then
467		call get_param(paramFile, mdl, "KPP_K_ENH_FAC", CS%KPP_K_ENH_FAC, &
468		'Constant value to enhance mixing coefficient in KPP.', &
469	0	units="nondim", default=1.0)
470		endif
471		endif
472		!/ 2. Options related to enhancing the unresolved Vt2/entrainment in Rib
473		call get_param(paramFile, mdl, "USE_KPP_LT_VT2", CS%LT_Vt2_Enhancement, &
474		'Flag for Langmuir turbulence enhancement of Vt2 '//&
475	0	'in Bulk Richardson Number.', Default=.false.)
476	0	if (CS%LT_Vt2_Enhancement) then
477		call get_param(paramFile, mdl, "KPP_LT_VT2_METHOD",string , &
478		'Method to enhance Vt2 in KPP. '// &
479		'Valid options are: \n'// &
480		'\t CONSTANT = Constant value (KPP_VT2_ENH_FAC) \n'// &
481		'\t VR12 = Function of Langmuir number based on VR12\n'// &
482		'\t (Van Roekel et al., 2012) \n'// &
483		'\t (Li et al. 2016, OM) \n'// &
484		'\t RW16 = Function of Langmuir number based on RW16\n'// &
485		'\t (Reichl et al., 2016, JPO) \n'// &
486		'\t LF17 = Function of Langmuir number based on LF17\n'// &
487		'\t (Li and Fox-Kemper, 2017, JPO)', &
488	0	default='CONSTANT')
489	0	select case ( trim(string))
490		case ("CONSTANT")
491	0	CS%LT_VT2_METHOD = LT_VT2_MODE_CONSTANT
492	0	langmuir_entrainment_opt = 'LWF16'
493		case ("VR12")
494	0	CS%LT_VT2_METHOD = LT_VT2_MODE_VR12
495	0	langmuir_entrainment_opt = 'LWF16'
496		case ("RW16")
497	0	CS%LT_VT2_METHOD = LT_VT2_MODE_RW16
498	0	langmuir_entrainment_opt = 'RWHGK16'
499		case ("LF17")
500	0	CS%LT_VT2_METHOD = LT_VT2_MODE_LF17
501	0	langmuir_entrainment_opt = 'LF17'
502		case default
503		call MOM_error(FATAL,"KPP_init: "//&
504	0	"Unrecognized KPP_LT_VT2_METHOD option: "//trim(string))
505		end select
506	0	if (CS%LT_VT2_METHOD==LT_VT2_MODE_CONSTANT) then
507		call get_param(paramFile, mdl, "KPP_VT2_ENH_FAC", CS%KPP_VT2_ENH_FAC, &
508		'Constant value to enhance VT2 in KPP.', &
509	0	units="nondim", default=1.0)
510		endif
511		endif
512
513	0	if (CS%LT_K_ENHANCEMENT .or. CS%LT_VT2_ENHANCEMENT) then
514		call get_param(paramFile, mdl, "KPP_LT_MLD_GUESS_MIN", CS%MLD_guess_min, &
515		"The minimum estimate of the mixed layer depth used to calculate "//&
516		"the Langmuir number for Langmuir turbulence enhancement with KPP.", &
517	0	units="m", default=1.0, scale=US%m_to_Z)
518		endif
519
520		call get_param(paramFile, mdl, "KPP_CVt2", CS%KPP_CVt2, &
521		'Parameter for Stokes MOST convection entrainment', &
522	0	units="nondim", default=1.6)
523
524		call get_param(paramFile, mdl, "ANSWER_DATE", CS%answer_date, &
525		"The vintage of the order of arithmetic in the CVMix KPP calculations. Values "//&
526		"below 20240501 recover the answers from early in 2024, while higher values "//&
527		"use expressions that have been refactored for rotational symmetry.", &
528	0	default=default_answer_date)
529
530	0	call closeParameterBlock(paramFile)
531
532	0	call get_param(paramFile, mdl, 'DEBUG', CS%debug, default=.False., do_not_log=.True.)
533
534		call CVMix_init_kpp( Ri_crit=CS%Ri_crit, &
535		minOBLdepth=US%Z_to_m*CS%minOBLdepth, &
536		minVtsqr=US%L_T_to_m_s*2CS%minVtsqr, &
537		vonKarman=CS%vonKarman, &
538		surf_layer_ext=CS%surf_layer_ext, &
539		CVt2=CS%KPP_CVt2, &
540		interp_type=CS%interpType, &
541		interp_type2=CS%interpType2, &
542		lEkman=CS%computeEkman, &
543		lStokesMOST=CS%StokesMOST, &
544		lMonOb=CS%computeMoninObukhov, &
545		MatchTechnique=CS%MatchTechnique, &
546		lenhanced_diff=CS%enhance_diffusion,&
547		lnonzero_surf_nonlocal=Cs_is_one ,&
548		lnoDGat1=lnoDGat1 ,&
549		langmuir_mixing_str=langmuir_mixing_opt,&
550		langmuir_entrainment_str=langmuir_entrainment_opt,&
551	0	CVMix_kpp_params_user=CS%KPP_params )
552
553		! Register diagnostics
554	0	CS%diag => diag
555		CS%id_OBLdepth = register_diag_field('ocean_model', 'KPP_OBLdepth', diag%axesT1, Time, &
556		'Thickness of the surface Ocean Boundary Layer calculated by [CVMix] KPP', &
557		'meter', conversion=US%Z_to_m, &
558		cmor_field_name='oml', cmor_long_name='ocean_mixed_layer_thickness_defined_by_mixing_scheme', &
559	0	cmor_units='m', cmor_standard_name='Ocean Mixed Layer Thickness Defined by Mixing Scheme')
560		! CMOR names are placeholders; must be modified by time period
561		! for CMOR compliance. Diag manager will be used for omlmax and
562		! omldamax.
563	0	if (CS%n_smooth > 0) then
564		CS%id_OBLdepth_original = register_diag_field('ocean_model', 'KPP_OBLdepth_original', diag%axesT1, Time, &
565		'Thickness of the surface Ocean Boundary Layer without smoothing calculated by [CVMix] KPP', &
566		'meter', conversion=US%Z_to_m, &
567		cmor_field_name='oml', cmor_long_name='ocean_mixed_layer_thickness_defined_by_mixing_scheme', &
568	0	cmor_units='m', cmor_standard_name='Ocean Mixed Layer Thickness Defined by Mixing Scheme')
569		endif
570	0	if ( CS%StokesMOST ) then
571		CS%id_StokesXI = register_diag_field('ocean_model', 'StokesXI', diag%axesT1, Time, &
572	0	'Stokes Similarity Parameter', 'nondim')
573		CS%id_Lam2 = register_diag_field('ocean_model', 'Lam2', diag%axesT1, Time, &
574	0	'Ustk0_ustar', 'nondim')
575		endif
576		CS%id_BulkDrho = register_diag_field('ocean_model', 'KPP_BulkDrho', diag%axesTL, Time, &
577		'Bulk difference in density used in Bulk Richardson number, as used by [CVMix] KPP', &
578	0	'kg/m3', conversion=US%R_to_kg_m3)
579		CS%id_BulkUz2 = register_diag_field('ocean_model', 'KPP_BulkUz2', diag%axesTL, Time, &
580		'Square of bulk difference in resolved velocity used in Bulk Richardson number via [CVMix] KPP', &
581	0	'm2/s2', conversion=US%L_T_to_m_s**2)
582		CS%id_BulkRi = register_diag_field('ocean_model', 'KPP_BulkRi', diag%axesTL, Time, &
583	0	'Bulk Richardson number used to find the OBL depth used by [CVMix] KPP', 'nondim')
584		CS%id_Sigma = register_diag_field('ocean_model', 'KPP_sigma', diag%axesTi, Time, &
585	0	'Sigma coordinate used by [CVMix] KPP', 'nondim')
586		CS%id_Ws = register_diag_field('ocean_model', 'KPP_Ws', diag%axesTL, Time, &
587		'Turbulent vertical velocity scale for scalars used by [CVMix] KPP', &
588	0	'm/s', conversion=US%Z_to_m*US%s_to_T)
589		CS%id_N = register_diag_field('ocean_model', 'KPP_N', diag%axesTi, Time, &
590	0	'(Adjusted) Brunt-Vaisala frequency used by [CVMix] KPP', '1/s', conversion=US%s_to_T)
591		CS%id_N2 = register_diag_field('ocean_model', 'KPP_N2', diag%axesTi, Time, &
592	0	'Square of Brunt-Vaisala frequency used by [CVMix] KPP', '1/s2', conversion=US%s_to_T**2)
593		CS%id_Vt2 = register_diag_field('ocean_model', 'KPP_Vt2', diag%axesTL, Time, &
594	0	'Unresolved shear turbulence used by [CVMix] KPP', 'm2/s2', conversion=US%Z_to_m*2US%s_to_T**2)
595		CS%id_uStar = register_diag_field('ocean_model', 'KPP_uStar', diag%axesT1, Time, &
596	0	'Friction velocity, u, as used by [CVMix] KPP', 'm/s', conversion=US%Z_to_mUS%s_to_T)
597		CS%id_buoyFlux = register_diag_field('ocean_model', 'KPP_buoyFlux', diag%axesTi, Time, &
598		'Surface (and penetrating) buoyancy flux, as used by [CVMix] KPP', &
599	0	'm2/s3', conversion=US%L_to_m*2US%s_to_T**3)
600		CS%id_Kt_KPP = register_diag_field('ocean_model', 'KPP_Kheat', diag%axesTi, Time, &
601		'Heat diffusivity due to KPP, as calculated by [CVMix] KPP', &
602	0	'm2/s', conversion=US%Z2_T_to_m2_s)
603		CS%id_Kd_in = register_diag_field('ocean_model', 'KPP_Kd_in', diag%axesTi, Time, &
604	0	'Diffusivity passed to KPP', 'm2/s', conversion=GV%HZ_T_to_m2_s)
605		CS%id_Ks_KPP = register_diag_field('ocean_model', 'KPP_Ksalt', diag%axesTi, Time, &
606		'Salt diffusivity due to KPP, as calculated by [CVMix] KPP', &
607	0	'm2/s', conversion=US%Z2_T_to_m2_s)
608		CS%id_Kv_KPP = register_diag_field('ocean_model', 'KPP_Kv', diag%axesTi, Time, &
609		'Vertical viscosity due to KPP, as calculated by [CVMix] KPP', &
610	0	'm2/s', conversion=US%Z2_T_to_m2_s)
611		CS%id_NLTt = register_diag_field('ocean_model', 'KPP_NLtransport_heat', diag%axesTi, Time, &
612	0	'Non-local transport (Cs*G(sigma)) for heat, as calculated by [CVMix] KPP', 'nondim')
613		CS%id_NLTs = register_diag_field('ocean_model', 'KPP_NLtransport_salt', diag%axesTi, Time, &
614	0	'Non-local tranpsort (Cs*G(sigma)) for scalars, as calculated by [CVMix] KPP', 'nondim')
615		CS%id_Tsurf = register_diag_field('ocean_model', 'KPP_Tsurf', diag%axesT1, Time, &
616		'Temperature of surface layer (10% of OBL depth) as passed to [CVMix] KPP', &
617	0	'C', conversion=US%C_to_degC)
618		CS%id_Ssurf = register_diag_field('ocean_model', 'KPP_Ssurf', diag%axesT1, Time, &
619		'Salinity of surface layer (10% of OBL depth) as passed to [CVMix] KPP', &
620	0	'ppt', conversion=US%S_to_ppt)
621		CS%id_Usurf = register_diag_field('ocean_model', 'KPP_Usurf', diag%axesCu1, Time, &
622		'i-component flow of surface layer (10% of OBL depth) as passed to [CVMix] KPP', &
623	0	'm/s', conversion=US%L_T_to_m_s)
624		CS%id_Vsurf = register_diag_field('ocean_model', 'KPP_Vsurf', diag%axesCv1, Time, &
625		'j-component flow of surface layer (10% of OBL depth) as passed to [CVMix] KPP', &
626	0	'm/s', conversion=US%L_T_to_m_s)
627		CS%id_EnhK = register_diag_field('ocean_model', 'EnhK', diag%axesTI, Time, &
628	0	'Langmuir number enhancement to K as used by [CVMix] KPP','nondim')
629		CS%id_EnhVt2 = register_diag_field('ocean_model', 'EnhVt2', diag%axesTL, Time, &
630	0	'Langmuir number enhancement to Vt2 as used by [CVMix] KPP','nondim')
631		CS%id_La_SL = register_diag_field('ocean_model', 'KPP_La_SL', diag%axesT1, Time, &
632	0	'Surface-layer Langmuir number computed in [CVMix] KPP','nondim')
633
634	0	allocate( CS%N( SZI_(G), SZJ_(G), SZK_(GV)+1 ), source=0. )
635	0	allocate( CS%StokesParXI( SZI_(G), SZJ_(G) ), source=0. )
636	0	allocate( CS%Lam2 ( SZI_(G), SZJ_(G) ), source=0. )
637	0	allocate( CS%kOBL( SZI_(G), SZJ_(G) ), source=0. )
638	0	allocate( CS%La_SL( SZI_(G), SZJ_(G) ), source=0. )
639	0	allocate( CS%Vt2( SZI_(G), SZJ_(G), SZK_(GV) ), source=0. )
640	0	if (CS%id_OBLdepth_original > 0) allocate( CS%OBLdepth_original( SZI_(G), SZJ_(G) ) )
641
642	0	allocate( CS%OBLdepthprev( SZI_(G), SZJ_(G) ), source=0.0 )
643	0	if (CS%id_BulkDrho > 0) allocate( CS%dRho( SZI_(G), SZJ_(G), SZK_(GV) ), source=0. )
644	0	if (CS%id_BulkUz2 > 0) allocate( CS%Uz2( SZI_(G), SZJ_(G), SZK_(GV) ), source=0. )
645	0	if (CS%id_BulkRi > 0) allocate( CS%BulkRi( SZI_(G), SZJ_(G), SZK_(GV) ), source=0. )
646	0	if (CS%id_Sigma > 0) allocate( CS%sigma( SZI_(G), SZJ_(G), SZK_(GV)+1 ), source=0. )
647	0	if (CS%id_Ws > 0) allocate( CS%Ws( SZI_(G), SZJ_(G), SZK_(GV) ), source=0. )
648	0	if (CS%id_N2 > 0) allocate( CS%N2( SZI_(G), SZJ_(G), SZK_(GV)+1 ), source=0. )
649	0	if (CS%id_Kt_KPP > 0) allocate( CS%Kt_KPP( SZI_(G), SZJ_(G), SZK_(GV)+1 ), source=0. )
650	0	if (CS%id_Ks_KPP > 0) allocate( CS%Ks_KPP( SZI_(G), SZJ_(G), SZK_(GV)+1 ), source=0. )
651	0	if (CS%id_Kv_KPP > 0) allocate( CS%Kv_KPP( SZI_(G), SZJ_(G), SZK_(GV)+1 ), source=0. )
652	0	if (CS%id_Tsurf > 0) allocate( CS%Tsurf( SZI_(G), SZJ_(G) ), source=0. )
653	0	if (CS%id_Ssurf > 0) allocate( CS%Ssurf( SZI_(G), SZJ_(G) ), source=0. )
654	0	if (CS%id_Usurf > 0) allocate( CS%Usurf( SZIB_(G), SZJ_(G) ), source=0. )
655	0	if (CS%id_Vsurf > 0) allocate( CS%Vsurf( SZI_(G), SZJB_(G) ), source=0. )
656	0	if (CS%id_EnhVt2 > 0) allocate( CS%EnhVt2( SZI_(G), SZJ_(G), SZK_(GV) ), source=0. )
657	0	if (CS%id_EnhK > 0) allocate( CS%EnhK( SZI_(G), SZJ_(G), SZK_(GV)+1 ), source=0. )
658
659	0	id_clock_KPP_calc = cpu_clock_id('Ocean KPP calculate)', grain=CLOCK_MODULE)
660	0	id_clock_KPP_compute_BLD = cpu_clock_id('(Ocean KPP comp BLD)', grain=CLOCK_ROUTINE)
661
662	0	end function KPP_init
663
664		!> KPP vertical diffusivity/viscosity and non-local tracer transport
665	0	subroutine KPP_calculate(CS, G, GV, US, h, tv, uStar, buoyFlux, Kt, Ks, Kv, &
666	0	nonLocalTransHeat, nonLocalTransScalar, Waves, lamult)
667
668		! Arguments
669		type(KPP_CS), pointer :: CS !< Control structure
670		type(ocean_grid_type), intent(in) :: G !< Ocean grid
671		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid
672		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
673		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
674		type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure.
675		real, dimension(SZI_(G),SZJ_(G)), intent(in) :: uStar !< Surface friction velocity [Z T-1 ~> m s-1]
676		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: buoyFlux !< Surface buoyancy flux [L2 T-3 ~> m2 s-3]
677		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: Kt !< (in) Vertical diffusivity of heat w/o KPP
678		!! (out) Vertical diffusivity including KPP
679		!! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
680		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: Ks !< (in) Vertical diffusivity of salt w/o KPP
681		!! (out) Vertical diffusivity including KPP
682		!! [H Z T-1 ~> m2 s-1 or kg m-1 s-1]
683		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: Kv !< (in) Vertical viscosity w/o KPP
684		!! (out) Vertical viscosity including KPP
685		!! [H Z T-1 ~> m2 s-1 or Pa s]
686		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: nonLocalTransHeat !< Temp non-local transport [nondim]
687		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: nonLocalTransScalar !< scalar non-local trans. [nondim]
688		type(wave_parameters_CS), pointer :: Waves !< Wave CS for Langmuir turbulence
689		real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: lamult !< Langmuir enhancement multiplier [nondim]
690
691		! Local variables
692		integer :: i, j, k ! Loop indices
693	0	real, dimension(SZI_(G),SZK_(GV)) :: dz ! Height change across layers [Z ~> m]
694	0	real, dimension( GV%ke ) :: cellHeight ! Cell center heights referenced to surface [Z ~> m] (negative in ocean)
695	0	real, dimension( GV%ke+1 ) :: iFaceHeight ! Interface heights referenced to surface [Z ~> m] (negative in ocean)
696	0	real, dimension( GV%ke ) :: z_cell ! Cell center heights referenced to surface [m] (negative in ocean)
697	0	real, dimension( GV%ke+1 ) :: z_inter ! Cell interface heights referenced to surface [m] (negative in ocean)
698	0	real, dimension( GV%ke+1, 2) :: Kdiffusivity ! Vertical diffusivity at interfaces in MKS units [m2 s-1]
699	0	real, dimension( GV%ke+1 ) :: Kviscosity ! Vertical viscosity at interfaces in MKS units [m2 s-1]
700	0	real, dimension( GV%ke+1, 2) :: nonLocalTrans ! Non-local transport for heat/salt at interfaces [nondim]
701
702		real :: surfFricVel ! Surface friction velocity in MKS units [m s-1]
703		real :: surfBuoyFlux ! Surface buoyancy flux in MKS units [m2 s-3]
704		real :: sigma ! Fractional vertical position within the boundary layer [nondim]
705		real :: sigmaRatio ! A cubic function of sigma [nondim]
706		real :: buoy_scale ! A unit conversion factor for buoyancy fluxes [m2 T3 L-2 s-3 ~> 1]
707		real :: dh ! The local thickness used for calculating interface positions [Z ~> m]
708		real :: hcorr ! A cumulative correction arising from inflation of vanished layers [Z ~> m]
709
710		! For Langmuir Calculations
711		real :: LangEnhK ! Langmuir enhancement for mixing coefficient [nondim]
712
713	0	if (CS%Stokes_Mixing .and. .not.associated(Waves)) call MOM_error(FATAL, &
714	0	"KPP_calculate: The Waves control structure must be associated if STOKES_MIXING is True.")
715
716	0	if (CS%debug) then
717	0	call hchksum(h, "KPP in: h", G%HI, haloshift=0, unscale=GV%H_to_m)
718	0	call hchksum(uStar, "KPP in: uStar", G%HI, haloshift=0, unscale=US%Z_to_m*US%s_to_T)
719	0	call hchksum(buoyFlux, "KPP in: buoyFlux", G%HI, haloshift=0, unscale=US%L_to_m*2US%s_to_T**3)
720	0	call hchksum(Kt, "KPP in: Kt", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
721	0	call hchksum(Ks, "KPP in: Ks", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
722		endif
723
724	0	nonLocalTrans(:,:) = 0.0
725
726	0	if (CS%id_Kd_in > 0) call post_data(CS%id_Kd_in, Kt, CS%diag)
727
728	0	call cpu_clock_begin(id_clock_KPP_calc)
729	0	buoy_scale = US%L_to_m*2US%s_to_T**3
730
731		!$OMP parallel do default(none) firstprivate(nonLocalTrans) &
732		!$OMP private(surfFricVel, iFaceHeight, hcorr, dh, dz, cellHeight, &
733		!$OMP surfBuoyFlux, Kdiffusivity, Kviscosity, LangEnhK, sigma, &
734		!$OMP sigmaRatio, z_inter, z_cell) &
735		!$OMP shared(G, GV, CS, US, tv, uStar, h, buoy_scale, buoyFlux, Kt, &
736		!$OMP Ks, Kv, nonLocalTransHeat, nonLocalTransScalar, Waves, lamult)
737		! loop over horizontal points on processor
738	0	do j = G%jsc, G%jec
739
740		! Find the vertical distances across layers.
741	0	call thickness_to_dz(h, tv, dz, j, G, GV)
742
743	0	do i = G%isc, G%iec ; if (G%mask2dT(i,j) > 0.0) then
744
745		! things independent of position within the column
746	0	surfFricVel = US%Z_to_mUS%s_to_T uStar(i,j)
747
748	0	iFaceHeight(1) = 0.0 ! BBL is all relative to the surface
749	0	hcorr = 0.
750	0	do k=1,GV%ke
751
752		! cell center and cell bottom in meters (negative values in the ocean)
753	0	dh = dz(i,k) ! Nominal thickness to use for increment
754	0	dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0)
755	0	hcorr = min( dh - CS%min_thickness, 0. ) ! If inflating then hcorr<0
756	0	dh = max( dh, CS%min_thickness ) ! Limit increment dh>=min_thickness
757	0	cellHeight(k) = iFaceHeight(k) - 0.5 * dh
758	0	iFaceHeight(k+1) = iFaceHeight(k) - dh
759
760		enddo ! k-loop finishes
761
762	0	surfBuoyFlux = buoy_scale*buoyFlux(i,j,1) ! This is only used in kpp_compute_OBL_depth to limit
763		! h to Monin-Obukhov (default is false, ie. not used)
764
765		! Call CVMix/KPP to obtain OBL diffusivities, viscosities and non-local transports
766
767		! Unlike LMD94, we do not match to interior diffusivities. If using the original
768		! LMD94 shape function, not matching is equivalent to matching to a zero diffusivity.
769
770		!BGR/ Add option for use of surface buoyancy flux with total sw flux.
771	0	if (CS%SW_METHOD == SW_METHOD_ALL_SW) then
772	0	surfBuoyFlux = buoy_scale * buoyFlux(i,j,1)
773	0	elseif (CS%SW_METHOD == SW_METHOD_MXL_SW) then
774		! We know the actual buoyancy flux into the OBL
775	0	surfBuoyFlux = buoy_scale * (buoyFlux(i,j,1) - buoyFlux(i,j,int(CS%kOBL(i,j))+1))
776	0	elseif (CS%SW_METHOD == SW_METHOD_LV1_SW) then
777	0	surfBuoyFlux = buoy_scale * (buoyFlux(i,j,1) - buoyFlux(i,j,2))
778		endif
779
780		! If option "MatchBoth" is selected in CVMix, MOM should be capable of matching.
781	0	if (.not. (CS%MatchTechnique == 'MatchBoth')) then
782	0	Kdiffusivity(:,:) = 0. ! Diffusivities for heat and salt [m2 s-1]
783	0	Kviscosity(:) = 0. ! Viscosity [m2 s-1]
784		else
785	0	Kdiffusivity(:,1) = GV%HZ_T_to_m2_s * Kt(i,j,:)
786	0	Kdiffusivity(:,2) = GV%HZ_T_to_m2_s * Ks(i,j,:)
787	0	Kviscosity(:) = GV%HZ_T_to_m2_s * Kv(i,j,:)
788		endif
789
790	0	IF (CS%LT_K_ENHANCEMENT) then
791	0	if (CS%LT_K_METHOD==LT_K_MODE_CONSTANT) then
792	0	LangEnhK = CS%KPP_K_ENH_FAC
793	0	elseif (CS%LT_K_METHOD==LT_K_MODE_VR12) then
794	0	if (present(lamult)) then
795	0	LangEnhK = lamult(i,j)
796		else
797		LangEnhK = sqrt(1.+(1.5CS%La_SL(i,j))*(-2) + &
798	0	(5.4CS%La_SL(i,j))*(-4))
799		endif
800	0	elseif (CS%LT_K_METHOD==LT_K_MODE_RW16) then
801		!This maximum value is proposed in Reichl et al., 2016 JPO formula
802	0	LangEnhK = min(2.25, 1. + 1./CS%La_SL(i,j))
803		else
804		!This shouldn't be reached.
805		!call MOM_error(WARNING,"Unexpected behavior in MOM_CVMix_KPP, see error in LT_K_ENHANCEMENT")
806	0	LangEnhK = 1.0
807		endif
808
809		! diffusivities don't need to be enhanced below anymore since LangEnhK is applied within CVMix.
810		! todo: need to double check if the distinction between the two different options of LT_K_SHAPE may need to be
811		! treated specially.
812	0	do k=1,GV%ke
813	0	if (CS%LT_K_SHAPE== LT_K_CONSTANT) then
814	0	if (CS%id_EnhK > 0) CS%EnhK(i,j,:) = LangEnhK
815		!Kdiffusivity(k,1) = Kdiffusivity(k,1) * LangEnhK
816		!Kdiffusivity(k,2) = Kdiffusivity(k,2) * LangEnhK
817		!Kviscosity(k) = Kviscosity(k) * LangEnhK
818	0	elseif (CS%LT_K_SHAPE == LT_K_SCALED) then
819	0	sigma = min(1.0,-iFaceHeight(k)/CS%OBLdepth(i,j))
820	0	SigmaRatio = sigma * (1. - sigma)**2 / 0.148148037
821	0	if (CS%id_EnhK > 0) CS%EnhK(i,j,k) = (1.0 + (LangEnhK - 1.)*sigmaRatio)
822		!Kdiffusivity(k,1) = Kdiffusivity(k,1) * ( 1. + &
823		! ( LangEnhK - 1.)*sigmaRatio)
824		!Kdiffusivity(k,2) = Kdiffusivity(k,2) * ( 1. + &
825		! ( LangEnhK - 1.)*sigmaRatio)
826		!Kviscosity(k) = Kviscosity(k) * ( 1. + &
827		! ( LangEnhK - 1.)*sigmaRatio)
828		endif
829		enddo
830		endif
831
832		! Convert columns to MKS units for passing to CVMix
833	0	do k = 1, GV%ke
834	0	z_cell(k) = US%Z_to_m*cellHeight(k)
835		enddo
836	0	do K = 1, GV%ke+1
837	0	z_inter(K) = US%Z_to_m*iFaceHeight(K)
838		enddo
839
840		call CVMix_coeffs_kpp(Kviscosity(:), & ! (inout) Total viscosity [m2 s-1]
841		Kdiffusivity(:,1), & ! (inout) Total heat diffusivity [m2 s-1]
842		Kdiffusivity(:,2), & ! (inout) Total salt diffusivity [m2 s-1]
843		z_inter(:), & ! (in) Height of interfaces [m]
844		z_cell(:), & ! (in) Height of level centers [m]
845		Kviscosity(:), & ! (in) Original viscosity [m2 s-1]
846		Kdiffusivity(:,1), & ! (in) Original heat diffusivity [m2 s-1]
847		Kdiffusivity(:,2), & ! (in) Original salt diffusivity [m2 s-1]
848		US%Z_to_m*CS%OBLdepth(i,j), & ! (in) OBL depth [m]
849		CS%kOBL(i,j), & ! (in) level (+fraction) of OBL extent
850		nonLocalTrans(:,1),& ! (out) Non-local heat transport [nondim]
851		nonLocalTrans(:,2),& ! (out) Non-local salt transport [nondim]
852		surfFricVel, & ! (in) Turbulent friction velocity at surface [m s-1]
853		surfBuoyFlux, & ! (in) Buoyancy flux at surface [m2 s-3]
854		GV%ke, & ! (in) Number of levels to compute coeffs for
855		GV%ke, & ! (in) Number of levels in array shape
856		Langmuir_EFactor=LangEnhK,& ! Langmuir enhancement multiplier
857		StokesXi = CS%StokesParXI(i,j), & ! Stokes forcing parameter
858	0	CVMix_kpp_params_user=CS%KPP_params )
859
860		! safety check, Kviscosity and Kdiffusivity must be >= 0
861	0	do k=1, GV%ke+1
862	0	if (Kviscosity(k) < 0. .or. Kdiffusivity(k,1) < 0.) then
863	0	write(*,'(a,3i3)') 'interface, i, j, k = ',j, j, k
864	0	write(*,'(a,2f12.5)') 'lon,lat=', G%geoLonT(i,j), G%geoLatT(i,j)
865	0	write(*,'(a,es12.4)') 'depth, z_inter(k) =',z_inter(k)
866	0	write(*,'(a,es12.4)') 'Kviscosity(k) =',Kviscosity(k)
867	0	write(*,'(a,es12.4)') 'Kdiffusivity(k,1) =',Kdiffusivity(k,1)
868	0	write(*,'(a,es12.4)') 'Kdiffusivity(k,2) =',Kdiffusivity(k,2)
869	0	write(,'(a,es12.4)') 'OBLdepth =',US%Z_to_mCS%OBLdepth(i,j)
870	0	write(*,'(a,f8.4)') 'kOBL =',CS%kOBL(i,j)
871	0	write(,'(a,es12.4)') 'u =',surfFricVel
872	0	write(*,'(a,es12.4)') 'bottom, z_inter(GV%ke+1) =',z_inter(GV%ke+1)
873	0	write(*,'(a,es12.4)') 'CS%La_SL(i,j) =',CS%La_SL(i,j)
874	0	write(*,'(a,es12.4)') 'LangEnhK =',LangEnhK
875	0	if (present(lamult)) write(*,'(a,es12.4)') 'lamult(i,j) =',lamult(i,j)
876	0	write(,) 'Kviscosity(:) =',Kviscosity(:)
877	0	write(,) 'Kdiffusivity(:,1) =',Kdiffusivity(:,1)
878
879		call MOM_error(FATAL,"KPP_calculate, after CVMix_coeffs_kpp: "// &
880		"Negative vertical viscosity or diffusivity has been detected. " // &
881		"This is likely related to the choice of MATCH_TECHNIQUE and INTERP_TYPE2. " //&
882	0	"You might consider using the default options for these parameters." )
883		endif
884		enddo
885
886		! Over-write CVMix NLT shape function with one of the following choices.
887		! The CVMix code has yet to update for thse options, so we compute in MOM6.
888		! Note that nonLocalTrans = Cs * G(sigma) (LMD94 notation), with
889		! Cs = 6.32739901508.
890		! Start do-loop at k=2, since k=1 is ocean surface (sigma=0)
891		! and we do not wish to double-count the surface forcing.
892		! Only compute nonlocal transport for 0 <= sigma <= 1.
893		! MOM6 recommended shape is the parabolic; it gives deeper boundary layer
894		! and no spurious extrema.
895	0	if (surfBuoyFlux < 0.0) then
896	0	if (CS%NLT_shape == NLT_SHAPE_CUBIC) then
897	0	do k = 2, GV%ke
898	0	sigma = min(1.0,-iFaceHeight(k)/CS%OBLdepth(i,j))
899	0	nonLocalTrans(k,1) = (1.0 - sigma)*2 (1.0 + 2.0sigma) !
900	0	nonLocalTrans(k,2) = nonLocalTrans(k,1)
901		enddo
902	0	elseif (CS%NLT_shape == NLT_SHAPE_PARABOLIC) then
903	0	do k = 2, GV%ke
904	0	sigma = min(1.0,-iFaceHeight(k)/CS%OBLdepth(i,j))
905	0	nonLocalTrans(k,1) = (1.0 - sigma)*2 !CS%CS2
906	0	nonLocalTrans(k,2) = nonLocalTrans(k,1)
907		enddo
908	0	elseif (CS%NLT_shape == NLT_SHAPE_LINEAR) then
909	0	do k = 2, GV%ke
910	0	sigma = min(1.0,-iFaceHeight(k)/CS%OBLdepth(i,j))
911	0	nonLocalTrans(k,1) = (1.0 - sigma)!*CS%CS2
912	0	nonLocalTrans(k,2) = nonLocalTrans(k,1)
913		enddo
914	0	elseif (CS%NLT_shape == NLT_SHAPE_CUBIC_LMD) then
915		! Sanity check (should agree with CVMix result using simple matching)
916	0	do k = 2, GV%ke
917	0	sigma = min(1.0,-iFaceHeight(k)/CS%OBLdepth(i,j))
918	0	nonLocalTrans(k,1) = CS%CS2 * sigma(1.0 -sigma)*2
919	0	nonLocalTrans(k,2) = nonLocalTrans(k,1)
920		enddo
921		endif
922		endif
923
924		! we apply nonLocalTrans in subroutines
925		! KPP_NonLocalTransport_temp and KPP_NonLocalTransport_saln
926	0	nonLocalTransHeat(i,j,:) = nonLocalTrans(:,1) ! temperature
927	0	nonLocalTransScalar(i,j,:) = nonLocalTrans(:,2) ! salinity
928
929		! set the KPP diffusivity and viscosity to zero for testing purposes
930	0	if (CS%KPPzeroDiffusivity) then
931	0	Kdiffusivity(:,1) = 0.0
932	0	Kdiffusivity(:,2) = 0.0
933	0	Kviscosity(:) = 0.0
934		endif
935
936		! Copy 1d data into 3d diagnostic arrays
937		!/ grabbing obldepth_0d for next time step.
938	0	CS%OBLdepthprev(i,j) = CS%OBLdepth(i,j)
939	0	if (CS%id_sigma > 0) then
940	0	CS%sigma(i,j,:) = 0.
941	0	if (CS%OBLdepth(i,j)>0.) CS%sigma(i,j,:) = -iFaceHeight(:)/CS%OBLdepth(i,j)
942		endif
943	0	if (CS%id_Kt_KPP > 0) CS%Kt_KPP(i,j,:) = US%m2_s_to_Z2_T * Kdiffusivity(:,1)
944	0	if (CS%id_Ks_KPP > 0) CS%Ks_KPP(i,j,:) = US%m2_s_to_Z2_T * Kdiffusivity(:,2)
945	0	if (CS%id_Kv_KPP > 0) CS%Kv_KPP(i,j,:) = US%m2_s_to_Z2_T * Kviscosity(:)
946
947		! Update output of routine
948	0	if (.not. CS%passiveMode) then
949	0	if (CS%KPPisAdditive) then
950	0	do k=1, GV%ke+1
951	0	Kt(i,j,k) = Kt(i,j,k) + GV%m2_s_to_HZ_T * Kdiffusivity(k,1)
952	0	Ks(i,j,k) = Ks(i,j,k) + GV%m2_s_to_HZ_T * Kdiffusivity(k,2)
953	0	Kv(i,j,k) = Kv(i,j,k) + GV%m2_s_to_HZ_T * Kviscosity(k)
954	0	if (CS%Stokes_Mixing) Waves%KvS(i,j,k) = Kv(i,j,k)
955		enddo
956		else ! KPP replaces prior diffusivity when former is non-zero
957	0	do k=1, GV%ke+1
958	0	if (Kdiffusivity(k,1) /= 0.) Kt(i,j,k) = GV%m2_s_to_HZ_T * Kdiffusivity(k,1)
959	0	if (Kdiffusivity(k,2) /= 0.) Ks(i,j,k) = GV%m2_s_to_HZ_T * Kdiffusivity(k,2)
960	0	if (Kviscosity(k) /= 0.) Kv(i,j,k) = GV%m2_s_to_HZ_T * Kviscosity(k)
961	0	if (CS%Stokes_Mixing) Waves%KvS(i,j,k) = Kv(i,j,k)
962		enddo
963		endif
964		endif
965
966
967		! end of the horizontal do-loops over the vertical columns
968		endif ; enddo ! i
969		enddo ! j
970
971	0	call cpu_clock_end(id_clock_KPP_calc)
972
973	0	if (CS%debug) then
974	0	call hchksum(Kt, "KPP out: Kt", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
975	0	call hchksum(Ks, "KPP out: Ks", G%HI, haloshift=0, unscale=GV%HZ_T_to_m2_s)
976		endif
977
978		! send diagnostics to post_data
979	0	if (CS%id_OBLdepth > 0) call post_data(CS%id_OBLdepth, CS%OBLdepth, CS%diag)
980	0	if (CS%id_OBLdepth_original > 0) call post_data(CS%id_OBLdepth_original,CS%OBLdepth_original,CS%diag)
981	0	if (CS%id_sigma > 0) call post_data(CS%id_sigma, CS%sigma, CS%diag)
982	0	if (CS%id_Ws > 0) call post_data(CS%id_Ws, CS%Ws, CS%diag)
983	0	if (CS%id_uStar > 0) call post_data(CS%id_uStar, uStar, CS%diag)
984	0	if (CS%id_buoyFlux > 0) call post_data(CS%id_buoyFlux, buoyFlux, CS%diag)
985	0	if (CS%id_Kt_KPP > 0) call post_data(CS%id_Kt_KPP, CS%Kt_KPP, CS%diag)
986	0	if (CS%id_Ks_KPP > 0) call post_data(CS%id_Ks_KPP, CS%Ks_KPP, CS%diag)
987	0	if (CS%id_Kv_KPP > 0) call post_data(CS%id_Kv_KPP, CS%Kv_KPP, CS%diag)
988	0	if (CS%id_NLTt > 0) call post_data(CS%id_NLTt, nonLocalTransHeat, CS%diag)
989	0	if (CS%id_NLTs > 0) call post_data(CS%id_NLTs, nonLocalTransScalar,CS%diag)
990
991
992	0	end subroutine KPP_calculate
993
994
995		!> Compute OBL depth
996	0	subroutine KPP_compute_BLD(CS, G, GV, US, h, Temp, Salt, u, v, tv, uStar, buoyFlux, Waves, lamult)
997
998		! Arguments
999		type(KPP_CS), pointer :: CS !< Control structure
1000		type(ocean_grid_type), intent(inout) :: G !< Ocean grid
1001		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid
1002		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1003		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
1004		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: Temp !< potential/cons temp [C ~> degC]
1005		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: Salt !< Salinity [S ~> ppt]
1006		real, dimension(SZIB_(G),SZJ_(G),SZK_(GV)), intent(in) :: u !< Velocity i-component [L T-1 ~> m s-1]
1007		real, dimension(SZI_(G),SZJB_(G),SZK_(GV)), intent(in) :: v !< Velocity j-component [L T-1 ~> m s-1]
1008		type(thermo_var_ptrs), intent(in) :: tv !< Thermodynamics structure.
1009		real, dimension(SZI_(G),SZJ_(G)), intent(in) :: uStar !< Surface friction velocity [Z T-1 ~> m s-1]
1010		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: buoyFlux !< Surface buoyancy flux [L2 T-3 ~> m2 s-3]
1011		type(wave_parameters_CS), pointer :: Waves !< Wave CS for Langmuir turbulence
1012		real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: lamult !< Langmuir enhancement factor [nondim]
1013
1014		! Local variables
1015	0	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dz ! Height change across layers [Z ~> m]
1016		! Variables for passing to CVMix routines, often in MKS units
1017	0	real, dimension( GV%ke ) :: Ws_1d ! Profile of vertical velocity scale for scalars in MKS units [m s-1]
1018	0	real, dimension( GV%ke ) :: deltaRho ! delta Rho in numerator of Bulk Ri number [R ~> kg m-3]
1019	0	real, dimension( GV%ke ) :: deltaBuoy ! Change in Buoyancy based on deltaRho [m s-2]
1020	0	real, dimension( GV%ke ) :: deltaU2 ! square of delta U (shear) in denominator of Bulk Ri [m2 s-2]
1021	0	real, dimension( GV%ke ) :: surfBuoyFlux2 ! Surface buoyancy flux in MKS units [m2 s-3]
1022	0	real, dimension( GV%ke ) :: BulkRi_1d ! Bulk Richardson number for each layer [nondim]
1023	0	real, dimension( GV%ke ) :: Vt2_1d ! Unresolved squared turbulence velocity for bulk Ri [m2 s-2]
1024	0	real, dimension( GV%ke ) :: z_cell ! Cell center heights referenced to surface [m] (negative in ocean)
1025	0	real, dimension( GV%ke ) :: OBL_depth ! Cell center depths referenced to surface [m] (positive in ocean)
1026	0	real, dimension( GV%ke+1 ) :: z_inter ! Cell interface heights referenced to surface [m] (negative in ocean)
1027	0	real, dimension( GV%ke+1 ) :: N_col ! A column of buoyancy frequencies at interfaces in MKS units [s-1]
1028		real :: surfFricVel ! Surface friction velocity in MKS units [m s-1]
1029		real :: surfBuoyFlux ! Surface buoyancy flux in MKS units [m2 s-3]
1030		real :: Coriolis ! Coriolis parameter at tracer points in MKS units [s-1]
1031		real :: KPP_OBL_depth ! Boundary layer depth calculated by CVMix_kpp_compute_OBL_depth in MKS units [m]
1032
1033		! Variables for EOS calculations
1034	0	real, dimension( 3*GV%ke ) :: rho_1D ! A column of densities [R ~> kg m-3]
1035	0	real, dimension( 3*GV%ke ) :: pres_1D ! A column of pressures [R L2 T-2 ~> Pa]
1036	0	real, dimension( 3*GV%ke ) :: Temp_1D ! A column of temperatures [C ~> degC]
1037	0	real, dimension( 3*GV%ke ) :: Salt_1D ! A column of salinities [S ~> ppt]
1038
1039	0	real, dimension( GV%ke ) :: cellHeight ! Cell center heights referenced to surface [Z ~> m] (negative in ocean)
1040	0	real, dimension( GV%ke+1 ) :: iFaceHeight ! Interface heights referenced to surface [Z ~> m] (negative in ocean)
1041	0	real, dimension( GV%ke+1 ) :: N2_1d ! Brunt-Vaisala frequency squared, at interfaces [T-2 ~> s-2]
1042		real :: zBottomMinusOffset ! Height of bottom plus a little bit [Z ~> m]
1043		real :: GoRho ! Gravitational acceleration in MKS units divided by density [m s-2 R-1 ~> m4 kg-1 s-2]
1044		real :: GoRho_Z_L2 ! Gravitational acceleration, perhaps divided by density, times aspect ratio
1045		! rescaling [H T-2 R-1 ~> m4 kg-1 s-2 or m s-2]
1046		real :: pRef ! The interface pressure [R L2 T-2 ~> Pa]
1047		real :: Uk, Vk ! Layer velocities relative to their averages in the surface layer [L T-1 ~> m s-1]
1048		real :: SLdepth_0d ! Surface layer depth = surf_layer_ext*OBLdepth [Z ~> m]
1049		real :: hTot ! Running sum of thickness used in the surface layer average [Z ~> m]
1050		real :: I_hTot ! The inverse of hTot [Z-1 ~> m-1]
1051		real :: buoy_scale ! A unit conversion factor for buoyancy fluxes [m2 T3 L-2 s-3 ~> 1]
1052		real :: delH ! Thickness of a layer [Z ~> m]
1053		real :: surfTemp ! Average of temperature over the surface layer [C ~> degC]
1054		real :: surfHtemp ! Integral of temperature over the surface layer [Z C ~> m degC]
1055		real :: surfSalt ! Average of salinity over the surface layer [S ~> ppt]
1056		real :: surfHsalt ! Integral of salinity over the surface layer [Z S ~> m ppt]
1057		real :: surfHu, surfHv ! Integral of u and v over the surface layer [Z L T-1 ~> m2 s-1]
1058		real :: surfU, surfV ! Average of u and v over the surface layer [Z T-1 ~> m s-1]
1059		real :: dh ! The local thickness used for calculating interface positions [Z ~> m]
1060		real :: hcorr ! A cumulative correction arising from inflation of vanished layers [Z ~> m]
1061
1062		! For Langmuir Calculations
1063		real :: Vt_layer ! non-dimensional extent contribution to unresolved shear
1064		real :: LangEnhVt2 ! Langmuir enhancement for unresolved shear [nondim]
1065	0	real, dimension(GV%ke) :: U_H, V_H ! Velocities at tracer points [L T-1 ~> m s-1]
1066		real :: MLD_guess ! A guess at the mixed layer depth for calculating the Langmuir number [Z ~> m]
1067		real :: LA ! The local Langmuir number [nondim]
1068		real :: surfHuS, surfHvS ! Stokes drift velocities integrated over the boundary layer [Z L T-1 ~> m2 s-1]
1069		real :: surfUs, surfVs ! Stokes drift velocities averaged over the boundary layer [Z T-1 ~> m s-1]
1070
1071		integer :: i, j, k, km1, kk, ksfc, ktmp ! Loop indices
1072
1073	0	real, dimension(GV%ke) :: uE_H, vE_H ! Eulerian velocities h-points, centers [L T-1 ~> m s-1]
1074	0	real, dimension(GV%ke) :: uS_H, vS_H ! Stokes drift components h-points, centers [L T-1 ~> m s-1]
1075	0	real, dimension(GV%ke) :: uSbar_H, vSbar_H ! Cell Average Stokes drift h-points [L T-1 ~> m s-1]
1076	0	real, dimension(GV%ke+1) :: uS_Hi, vS_Hi ! Stokes Drift components at interfaces [L T-1 ~> m s-1]
1077		real :: uS_SLD , vS_SLD, uS_SLC , vS_SLC, uSbar_SLD, vSbar_SLD ! Stokes at/to to Surface Layer Extent
1078		! [L T-1 ~> m s-1]
1079		real :: StokesXI ! Stokes similarity parameter [nondim]
1080	0	real, dimension( GV%ke ) :: StokesXI_1d , StokesVt_1d ! Parameters of TKE production ratio [nondim]
1081		integer :: kbl ! index of cell containing boundary layer depth
1082
1083	0	if (CS%Stokes_Mixing .and. .not.associated(Waves)) call MOM_error(FATAL, &
1084	0	"KPP_compute_BLD: The Waves control structure must be associated if STOKES_MIXING is True.")
1085
1086	0	if (CS%debug) then
1087	0	call hchksum(Salt, "KPP in: S", G%HI, haloshift=0, unscale=US%S_to_ppt)
1088	0	call hchksum(Temp, "KPP in: T", G%HI, haloshift=0, unscale=US%C_to_degC)
1089	0	call hchksum(u, "KPP in: u", G%HI, haloshift=0, unscale=US%L_T_to_m_s)
1090	0	call hchksum(v, "KPP in: v", G%HI, haloshift=0, unscale=US%L_T_to_m_s)
1091		endif
1092
1093	0	call cpu_clock_begin(id_clock_KPP_compute_BLD)
1094
1095		! some constants
1096	0	GoRho = US%Z_to_mUS%s_to_T2 (GV%g_Earth_Z_T2 / GV%Rho0)
1097	0	if (GV%Boussinesq) then
1098	0	GoRho_Z_L2 = GV%Z_to_H * GV%g_Earth_Z_T2 / GV%Rho0
1099		else
1100	0	GoRho_Z_L2 = GV%g_Earth_Z_T2 * GV%RZ_to_H
1101		endif
1102	0	buoy_scale = US%L_to_m*2US%s_to_T**3
1103
1104		! Find the vertical distances across layers.
1105	0	call thickness_to_dz(h, tv, dz, G, GV, US)
1106
1107		! loop over horizontal points on processor
1108		!$OMP parallel do default(none) private(surfFricVel, iFaceHeight, hcorr, dh, cellHeight, &
1109		!$OMP surfBuoyFlux, U_H, V_H, Coriolis, pRef, SLdepth_0d, vt2_1d, &
1110		!$OMP ksfc, surfHtemp, surfHsalt, surfHu, surfHv, surfHuS, &
1111		!$OMP surfHvS, hTot, I_hTot, delH, surftemp, surfsalt, surfu, surfv, &
1112		!$OMP surfUs, surfVs, Uk, Vk, deltaU2, km1, kk, pres_1D, N_col, &
1113		!$OMP Temp_1D, salt_1D, surfBuoyFlux2, MLD_guess, LA, rho_1D, &
1114		!$OMP deltarho, deltaBuoy, N2_1d, ws_1d, LangEnhVT2,KPP_OBL_depth, z_cell, &
1115		!$OMP z_inter, OBL_depth, BulkRi_1d, zBottomMinusOffset, uE_H, vE_H, &
1116		!$OMP uS_H, vS_H, uSbar_H, vSbar_H , uS_Hi, vS_Hi, &
1117		!$OMP uS_SLD, vS_SLD, uS_SLC, vS_SLC, uSbar_SLD, vSbar_SLD, &
1118		!$OMP StokesXI, StokesXI_1d, StokesVt_1d, kbl) &
1119		!$OMP shared(G, GV, CS, US, uStar, h, dz, buoy_scale, buoyFlux, &
1120		!$OMP Temp, Salt, waves, tv, GoRho, GoRho_Z_L2, u, v, lamult, Vt_layer)
1121	0	do j = G%jsc, G%jec
1122	0	do i = G%isc, G%iec ; if (G%mask2dT(i,j) > 0.0) then
1123
1124	0	do k=1,GV%ke
1125	0	U_H(k) = 0.5 * (u(I,j,k)+u(I-1,j,k))
1126	0	V_H(k) = 0.5 * (v(i,J,k)+v(i,J-1,k))
1127		enddo
1128	0	if (CS%StokesMOST) then
1129	0	do k=1,GV%ke
1130	0	uE_H(k) = 0.5 * (u(I,j,k)+u(I-1,j,k)-Waves%US_x(I,j,k)-Waves%US_x(I-1,j,k))
1131	0	vE_H(k) = 0.5 * (v(i,J,k)+v(i,J-1,k)-Waves%US_y(i,J,k)-Waves%US_y(i,J-1,k))
1132		enddo
1133		endif
1134		! things independent of position within the column
1135		Coriolis = 0.25US%s_to_T( (G%CoriolisBu(i,j) + G%CoriolisBu(i-1,j-1)) + &
1136	0	(G%CoriolisBu(i-1,j) + G%CoriolisBu(i,j-1)) )
1137	0	surfFricVel = US%Z_to_mUS%s_to_T uStar(i,j)
1138
1139		! Bulk Richardson number computed for each cell in a column,
1140		! assuming OBLdepth = grid cell depth. After Rib(k) is
1141		! known for the column, then CVMix interpolates to find
1142		! the actual OBLdepth. This approach avoids need to iterate
1143		! on the OBLdepth calculation. It follows that used in MOM5
1144		! and POP.
1145	0	iFaceHeight(1) = 0.0 ! BBL is all relative to the surface
1146	0	pRef = 0. ; if (associated(tv%p_surf)) pRef = tv%p_surf(i,j)
1147	0	hcorr = 0.
1148
1149	0	if (CS%StokesMOST) call Compute_StokesDrift( i, j, h(i,j,1) , iFaceHeight(1), &
1150	0	uS_Hi(1), vS_Hi(1), uS_H(1), vS_H(1), uSbar_H(1), vSbar_H(1), Waves)
1151
1152	0	do k=1,GV%ke
1153		! cell center and cell bottom in meters (negative values in the ocean)
1154	0	dh = dz(i,j,k) ! Nominal thickness to use for increment
1155	0	dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0)
1156	0	hcorr = min( dh - CS%min_thickness, 0. ) ! If inflating then hcorr<0
1157	0	dh = max( dh, CS%min_thickness ) ! Limit increment dh>=min_thickness
1158	0	cellHeight(k) = iFaceHeight(k) - 0.5 * dh
1159	0	iFaceHeight(k+1) = iFaceHeight(k) - dh
1160
1161		! find ksfc for cell where "surface layer" sits
1162	0	SLdepth_0d = CS%surf_layer_ext*max( max(-cellHeight(k),-iFaceHeight(2) ), CS%minOBLdepth )
1163	0	ksfc = k
1164	0	do ktmp = 1,k
1165	0	if (-1.0*iFaceHeight(ktmp+1) >= SLdepth_0d) then
1166	0	ksfc = ktmp
1167	0	exit
1168		endif
1169		enddo
1170
1171	0	if (CS%StokesMOST) then
1172		! if k=1, want buoyFlux(i,j,1) - buoyFlux(i,j,2), otherwise
1173		! subtract average of buoyFlux(i,j,k) and buoyFlux(i,j,k+1)
1174		surfBuoyFlux = buoy_scale * &
1175	0	(buoyFlux(i,j,1) - 0.5*(buoyFlux(i,j,max(2,k))+buoyFlux(i,j,k+1)) )
1176	0	surfBuoyFlux2(k) = surfBuoyFlux
1177		call Compute_StokesDrift(i,j, iFaceHeight(k),iFaceHeight(k+1), &
1178	0	uS_Hi(k+1), vS_Hi(k+1), uS_H(k), vS_H(k), uSbar_H(k), vSbar_H(k), Waves)
1179		call Compute_StokesDrift(i,j, iFaceHeight(ksfc) , -SLdepth_0d, &
1180	0	uS_SLD , vS_SLD, uS_SLC , vS_SLC, uSbar_SLD, vSbar_SLD, Waves)
1181		call cvmix_kpp_compute_StokesXi( iFaceHeight,CellHeight,ksfc ,SLdepth_0d,surfBuoyFlux, &
1182		surfFricVel,waves%omega_w2x(i,j), uE_H, vE_H, uS_Hi, vS_Hi, uSbar_H, vSbar_H, uS_SLD,&
1183	0	vS_SLD, uSbar_SLD, vSbar_SLD, StokesXI, CVMix_kpp_params_user=CS%KPP_params )
1184	0	StokesXI_1d(k) = StokesXI
1185	0	StokesVt_1d(k) = 0.0 ! StokesXI
1186
1187		! average temperature, salinity, u and v over surface layer starting at ksfc
1188	0	delH = SLdepth_0d + iFaceHeight(ksfc)
1189	0	surfHtemp = Temp(i,j,ksfc) * delH
1190	0	surfHsalt = Salt(i,j,ksfc) * delH
1191	0	surfHu = (uE_H(ksfc) + uSbar_SLD) * delH
1192	0	surfHv = (vE_H(ksfc) + vSbar_SLD) * delH
1193	0	hTot = delH
1194	0	do ktmp = 1,ksfc-1 ! if ksfc >=2
1195	0	delH = h(i,j,ktmp)*GV%H_to_Z
1196	0	hTot = hTot + delH
1197	0	surfHtemp = surfHtemp + Temp(i,j,ktmp) * delH
1198	0	surfHsalt = surfHsalt + Salt(i,j,ktmp) * delH
1199	0	surfHu = surfHu + (uE_H(ktmp) + uSbar_H(ktmp)) * delH
1200	0	surfHv = surfHv + (vE_H(ktmp) + vSbar_H(ktmp)) * delH
1201		enddo
1202	0	I_hTot = 1./hTot
1203	0	surfTemp = surfHtemp * I_hTot
1204	0	surfSalt = surfHsalt * I_hTot
1205	0	surfU = surfHu * I_hTot
1206	0	surfV = surfHv * I_hTot
1207
1208	0	Uk = uE_H(k) + uS_H(k) - surfU
1209	0	Vk = vE_H(k) + vS_H(k) - surfV
1210
1211		else !not StokesMOST
1212	0	StokesXI_1d(k) = 0.0
1213		! average temperature, salinity, u and v over surface layer
1214		! use C-grid average to get u and v on T-points.
1215	0	surfHtemp = 0.0
1216	0	surfHsalt = 0.0
1217	0	surfHu = 0.0
1218	0	surfHv = 0.0
1219	0	surfHuS = 0.0
1220	0	surfHvS = 0.0
1221	0	hTot = 0.0
1222	0	do ktmp = 1,ksfc
1223
1224		! SLdepth_0d can be between cell interfaces
1225	0	delH = min( max(0.0, SLdepth_0d - hTot), dz(i,j,ktmp) )
1226
1227		! surface layer thickness
1228	0	hTot = hTot + delH
1229
1230		! surface averaged fields
1231	0	surfHtemp = surfHtemp + Temp(i,j,ktmp) * delH
1232	0	surfHsalt = surfHsalt + Salt(i,j,ktmp) * delH
1233	0	surfHu = surfHu + 0.5(u(i,j,ktmp)+u(i-1,j,ktmp)) delH
1234	0	surfHv = surfHv + 0.5(v(i,j,ktmp)+v(i,j-1,ktmp)) delH
1235	0	if (CS%Stokes_Mixing) then
1236	0	surfHus = surfHus + 0.5(Waves%US_x(i,j,ktmp)+Waves%US_x(i-1,j,ktmp)) delH
1237	0	surfHvs = surfHvs + 0.5(Waves%US_y(i,j,ktmp)+Waves%US_y(i,j-1,ktmp)) delH
1238		endif
1239
1240		enddo
1241		!I_hTot = 1./hTot
1242		!surfTemp = surfHtemp * I_hTot
1243		!surfSalt = surfHsalt * I_hTot
1244		!surfU = surfHu * I_hTot
1245		!surfV = surfHv * I_hTot
1246		!surfUs = surfHus * I_hTot
1247		!surfVs = surfHvs * I_hTot
1248
1249	0	surfTemp = surfHtemp / hTot
1250	0	surfSalt = surfHsalt / hTot
1251	0	surfU = surfHu / hTot
1252	0	surfV = surfHv / hTot
1253	0	surfUs = surfHus / hTot
1254	0	surfVs = surfHvs / hTot
1255		! vertical shear between present layer and surface layer averaged surfU and surfV.
1256		! C-grid average to get Uk and Vk on T-points.
1257	0	Uk = 0.5*(u(i,j,k)+u(i-1,j,k)) - surfU
1258	0	Vk = 0.5*(v(i,j,k)+v(i,j-1,k)) - surfV
1259
1260	0	if (CS%Stokes_Mixing) then
1261		! If momentum is mixed down the Stokes drift gradient, then
1262		! the Stokes drift must be included in the bulk Richardson number
1263		! calculation.
1264	0	Uk = Uk + (0.5*(Waves%Us_x(i,j,k)+Waves%US_x(i-1,j,k)) - surfUs )
1265	0	Vk = Vk + (0.5*(Waves%Us_y(i,j,k)+Waves%Us_y(i,j-1,k)) - surfVs )
1266		endif
1267
1268		! this difference accounts for penetrating SW
1269	0	surfBuoyFlux = buoy_scale * (buoyFlux(i,j,1) - buoyFlux(i,j,k+1))
1270	0	surfBuoyFlux2(k) = surfBuoyFlux
1271
1272		endif ! StokesMOST
1273
1274	0	deltaU2(k) = US%L_T_to_m_s*2 ((Uk2) + (Vk2))
1275
1276		! pressure, temperature, and salinity for calling the equation of state
1277		! kk+1 = surface fields
1278		! kk+2 = k fields
1279		! kk+3 = km1 fields
1280	0	km1 = max(1, k-1)
1281	0	kk = 3*(k-1)
1282	0	pres_1D(kk+1) = pRef
1283	0	pres_1D(kk+2) = pRef
1284	0	pres_1D(kk+3) = pRef
1285	0	Temp_1D(kk+1) = surfTemp
1286	0	Temp_1D(kk+2) = Temp(i,j,k)
1287	0	Temp_1D(kk+3) = Temp(i,j,km1)
1288	0	Salt_1D(kk+1) = surfSalt
1289	0	Salt_1D(kk+2) = Salt(i,j,k)
1290	0	Salt_1D(kk+3) = Salt(i,j,km1)
1291
1292		! pRef is pressure at interface between k and km1 [R L2 T-2 ~> Pa].
1293		! iterate pRef for next pass through k-loop.
1294	0	pRef = pRef + (GV%g_Earth * GV%H_to_RZ) * h(i,j,k)
1295
1296		enddo ! k-loop finishes
1297
1298	0	if ( (CS%LT_K_ENHANCEMENT .or. CS%LT_VT2_ENHANCEMENT) .and. .not. present(lamult)) then
1299	0	MLD_guess = max( CS%MLD_guess_min, abs(CS%OBLdepthprev(i,j) ) )
1300		call get_Langmuir_Number(LA, G, GV, US, MLD_guess, uStar(i,j), i, j, &
1301	0	dz=dz(i,j,:), U_H=U_H, V_H=V_H, WAVES=WAVES)
1302	0	CS%La_SL(i,j) = LA
1303		endif
1304
1305
1306		! compute in-situ density
1307	0	call calculate_density(Temp_1D, Salt_1D, pres_1D, rho_1D, tv%eqn_of_state)
1308
1309		! N2 (can be negative) and N (non-negative) on interfaces.
1310		! deltaRho is non-local rho difference used for bulk Richardson number.
1311		! CS%N is local N (with floor) used for unresolved shear calculation.
1312	0	do k = 1, GV%ke
1313	0	km1 = max(1, k-1)
1314	0	kk = 3*(k-1)
1315	0	deltaRho(k) = rho_1D(kk+2) - rho_1D(kk+1)
1316	0	if (GV%Boussinesq .or. GV%semi_Boussinesq) then
1317	0	deltaBuoy(k) = GoRho*(rho_1D(kk+2) - rho_1D(kk+1))
1318		else
1319		deltaBuoy(k) = (US%Z_to_mUS%s_to_T2) GV%g_Earth_Z_T2 * &
1320	0	( (rho_1D(kk+2) - rho_1D(kk+1)) / (0.5 * (rho_1D(kk+2) + rho_1D(kk+1))) )
1321		endif
1322		N2_1d(k) = (GoRho_Z_L2 * (rho_1D(kk+2) - rho_1D(kk+3)) ) / &
1323	0	((0.5*(h(i,j,km1) + h(i,j,k))+GV%H_subroundoff))
1324	0	CS%N(i,j,k) = sqrt( max( N2_1d(k), 0.) )
1325		enddo
1326	0	N2_1d(GV%ke+1 ) = 0.0
1327	0	CS%N(i,j,GV%ke+1 ) = 0.0
1328
1329		! Convert columns to MKS units for passing to CVMix
1330	0	do k = 1, GV%ke
1331	0	OBL_depth(k) = -US%Z_to_m * cellHeight(k)
1332	0	z_cell(k) = US%Z_to_m*cellHeight(k)
1333		enddo
1334	0	do K = 1, GV%ke+1
1335	0	N_col(K) = US%s_to_T*CS%N(i,j,K)
1336	0	z_inter(K) = US%Z_to_m*iFaceHeight(K)
1337		enddo
1338
1339		! CVMix_kpp_compute_turbulent_scales_1d_OBL computes w_s velocity scale at cell centers for
1340		! CVmix_kpp_compute_bulk_Richardson call to CVmix_kpp_compute_unresolved_shear
1341		! at sigma=Vt_layer (CS%surf_layer_ext or 1.0) for this calculation.
1342		! StokesVt_1d controls Stokes enhancement (= 0 for none)
1343	0	Vt_layer = 1.0 ! CS%surf_layer_ext
1344		call CVMix_kpp_compute_turbulent_scales( & ! 1d_OBL
1345		Vt_layer, & ! (in) Boundary layer extent contributing to unresolved shear
1346		OBL_depth, & ! (in) OBL depth [m]
1347		surfBuoyFlux2, & ! (in) Buoyancy flux at surface [m2 s-3]
1348		surfFricVel, & ! (in) Turbulent friction velocity at surface [m s-1]
1349		xi=StokesVt_1d, & ! (in) Stokes similarity parameter-->1/CHI(xi) enhance of Vt
1350		w_s=Ws_1d, & ! (out) Turbulent velocity scale profile [m s-1]
1351	0	CVMix_kpp_params_user=CS%KPP_params )
1352
1353		! Determine the enhancement factor for unresolved shear
1354	0	IF (CS%LT_VT2_ENHANCEMENT) then
1355	0	IF (CS%LT_VT2_METHOD==LT_VT2_MODE_CONSTANT) then
1356	0	LangEnhVT2 = CS%KPP_VT2_ENH_FAC
1357	0	elseif (CS%LT_VT2_METHOD==LT_VT2_MODE_VR12) then
1358		!Introduced minimum value for La_SL, so maximum value for enhvt2 is removed.
1359	0	if (present(lamult)) then
1360	0	LangEnhVT2 = lamult(i,j)
1361		else
1362		LangEnhVT2 = sqrt(1.+(1.5CS%La_SL(i,j))*(-2) + &
1363	0	(5.4CS%La_SL(i,j))*(-4))
1364		endif
1365		else
1366		! for other methods (e.g., LT_VT2_MODE_RW16, LT_VT2_MODE_LF17), the enhancement factor is
1367		! computed internally within CVMix using LaSL, bfsfc, and ustar to be passed to CVMix.
1368	0	LangEnhVT2 = 1.0
1369		endif
1370		else
1371	0	LangEnhVT2 = 1.0
1372		endif
1373
1374	0	surfBuoyFlux = buoy_scale * buoyFlux(i,j,1)
1375
1376		! Calculate Bulk Richardson number from eq (21) of LMD94
1377		BulkRi_1d = CVmix_kpp_compute_bulk_Richardson( &
1378		zt_cntr=z_cell, & ! Depth of cell center [m]
1379		delta_buoy_cntr=deltaBuoy, & ! Bulk buoyancy difference, Br-B(z) [m s-2]
1380		delta_Vsqr_cntr=deltaU2, & ! Square of resolved velocity difference [m2 s-2]
1381		ws_cntr=Ws_1d, & ! Turbulent velocity scale profile [m s-1]
1382		N_iface=N_col, & ! Buoyancy frequency [s-1]
1383		EFactor=LangEnhVT2, & ! Langmuir enhancement factor [nondim]
1384		LaSL=CS%La_SL(i,j), & ! surface layer averaged Langmuir number [nondim]
1385		bfsfc=surfBuoyFlux2, & ! surface buoyancy flux [m2 s-3]
1386		uStar=surfFricVel, & ! surface friction velocity [m s-1]
1387	0	CVMix_kpp_params_user=CS%KPP_params ) ! KPP parameters
1388
1389		! ! A hack to avoid KPP reaching the bottom. It was needed during development
1390		! ! because KPP was unable to handle vanishingly small layers near the bottom.
1391		! if (CS%deepOBLoffset>0.) then
1392		! zBottomMinusOffset = iFaceHeight(GV%ke+1) + min(CS%deepOBLoffset, -0.1*iFaceHeight(GV%ke+1))
1393		! CS%OBLdepth(i,j) = min( CS%OBLdepth(i,j), -zBottomMinusOffset )
1394		! endif
1395	0	zBottomMinusOffset = iFaceHeight(GV%ke+1) + min(CS%deepOBLoffset,-0.1*iFaceHeight(GV%ke+1))
1396
1397		call CVMix_kpp_compute_OBL_depth( &
1398		BulkRi_1d, & ! (in) Bulk Richardson number
1399		z_inter, & ! (in) Height of interfaces [m]
1400		KPP_OBL_depth, & ! (out) OBL depth [m]
1401		CS%kOBL(i,j), & ! (out) level (+fraction) of OBL extent
1402		zt_cntr=z_cell, & ! (in) Height of cell centers [m]
1403		surf_fric=surfFricVel, & ! (in) Turbulent friction velocity at surface [m s-1]
1404		surf_buoy=surfBuoyFlux2, & ! (in) Buoyancy flux at surface [m2 s-3]
1405		Coriolis=Coriolis, & ! (in) Coriolis parameter [s-1]
1406		Xi = StokesXI_1d, & ! (in) Stokes similarity parameter Lmob limit (1-Xi)
1407		zBottom = zBottomMinusOffset, & ! (in) Numerical limit on OBLdepth
1408	0	CVMix_kpp_params_user=CS%KPP_params ) ! KPP parameters
1409	0	CS%OBLdepth(i,j) = US%m_to_Z * KPP_OBL_depth
1410
1411	0	if (CS%StokesMOST) then
1412	0	kbl = int(CS%kOBL(i,j))
1413	0	SLdepth_0d = CS%surf_layer_ext*CS%OBLdepth(i,j)
1414	0	surfBuoyFlux = surfBuoyFlux2(kbl)
1415		! find ksfc for cell where "surface layer" sits
1416	0	ksfc = kbl
1417	0	do ktmp = 1, kbl
1418	0	if (-1.0*iFaceHeight(ktmp+1) >= SLdepth_0d) then
1419	0	ksfc = ktmp
1420	0	exit
1421		endif
1422		enddo
1423
1424		call Compute_StokesDrift(i,j, iFaceHeight(ksfc) , -SLdepth_0d, &
1425	0	uS_SLD , vS_SLD, uS_SLC , vS_SLC, uSbar_SLD, vSbar_SLD, Waves)
1426		call cvmix_kpp_compute_StokesXi( iFaceHeight,CellHeight,ksfc ,SLdepth_0d, &
1427		surfBuoyFlux, surfFricVel,waves%omega_w2x(i,j), uE_H, vE_H, uS_Hi, &
1428		vS_Hi, uSbar_H, vSbar_H, uS_SLD, vS_SLD, uSbar_SLD, vSbar_SLD, &
1429	0	StokesXI, CVMix_kpp_params_user=CS%KPP_params )
1430	0	CS%StokesParXI(i,j) = StokesXI
1431	0	CS%Lam2(i,j) = sqrt(US_Hi(1)2+VS_Hi(1)2) / MAX(surfFricVel,0.0002)
1432
1433		else !.not Stokes_MOST
1434	0	CS%StokesParXI(i,j) = 10.0
1435	0	CS%Lam2(i,j) = sqrt(US_Hi(1)2+VS_Hi(1)2) / MAX(surfFricVel,0.0002)
1436
1437		! A hack to avoid KPP reaching the bottom. It was needed during development
1438		! because KPP was unable to handle vanishingly small layers near the bottom.
1439	0	if (CS%deepOBLoffset>0.) then
1440	0	zBottomMinusOffset = iFaceHeight(GV%ke+1) + min(CS%deepOBLoffset, -0.1*iFaceHeight(GV%ke+1))
1441	0	CS%OBLdepth(i,j) = min( CS%OBLdepth(i,j), -zBottomMinusOffset )
1442		endif
1443
1444		! apply some constraints on OBLdepth
1445	0	if (CS%fixedOBLdepth) CS%OBLdepth(i,j) = CS%fixedOBLdepth_value
1446	0	CS%OBLdepth(i,j) = max( CS%OBLdepth(i,j), -iFaceHeight(2) ) ! no shallower than top layer
1447	0	CS%OBLdepth(i,j) = min( CS%OBLdepth(i,j), -iFaceHeight(GV%ke+1) ) ! no deeper than bottom
1448	0	CS%kOBL(i,j) = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, CS%OBLdepth(i,j) )
1449
1450		endif !Stokes_MOST
1451
1452		! compute unresolved squared velocity for diagnostics
1453	0	if (CS%id_Vt2 > 0) then
1454		Vt2_1d(:) = CVmix_kpp_compute_unresolved_shear( &
1455		z_cell, & ! Depth of cell center [m]
1456		ws_cntr=Ws_1d, & ! Turbulent velocity scale profile, at centers [m s-1]
1457		N_iface=N_col, & ! Buoyancy frequency at interface [s-1]
1458		EFactor=LangEnhVT2, & ! Langmuir enhancement factor [nondim]
1459		LaSL=CS%La_SL(i,j), & ! surface layer averaged Langmuir number [nondim]
1460		bfsfc=surfBuoyFlux2, & ! surface buoyancy flux [m2 s-3]
1461		uStar=surfFricVel, & ! surface friction velocity [m s-1]
1462	0	CVmix_kpp_params_user=CS%KPP_params ) ! KPP parameters
1463	0	CS%Vt2(i,j,:) = US%m_to_Z*2US%T_to_s*2 Vt2_1d(:)
1464		endif
1465
1466		! recompute wscale for diagnostics, now that we in fact know boundary layer depth
1467		!BGR consider if LTEnhancement is wanted for diagnostics
1468	0	if (CS%id_Ws > 0) then
1469		call CVMix_kpp_compute_turbulent_scales( &
1470		-cellHeight(:)/CS%OBLdepth(i,j), & ! (in) Normalized boundary layer coordinate [nondim]
1471		US%Z_to_m*CS%OBLdepth(i,j), & ! (in) OBL depth [m]
1472		surfBuoyFlux, & ! (in) Buoyancy flux at surface [m2 s-3]
1473		surfFricVel, & ! (in) Turbulent friction velocity at surface [m s-1]
1474		xi=StokesXI, & ! (in) Stokes similarity parameter-->1/CHI(xi) enhance
1475		w_s=Ws_1d, & ! (out) Turbulent velocity scale profile [m s-1]
1476	0	CVMix_kpp_params_user=CS%KPP_params) ! KPP parameters
1477	0	CS%Ws(i,j,:) = US%m_to_ZUS%T_to_sWs_1d(:)
1478		endif
1479
1480		! Diagnostics
1481	0	if (CS%id_N2 > 0) CS%N2(i,j,:) = N2_1d(:)
1482	0	if (CS%id_BulkDrho > 0) CS%dRho(i,j,:) = deltaRho(:)
1483	0	if (CS%id_BulkRi > 0) CS%BulkRi(i,j,:) = BulkRi_1d(:)
1484	0	if (CS%id_BulkUz2 > 0) CS%Uz2(i,j,:) = US%m_s_to_L_T*2 deltaU2(:)
1485	0	if (CS%id_Tsurf > 0) CS%Tsurf(i,j) = surfTemp
1486	0	if (CS%id_Ssurf > 0) CS%Ssurf(i,j) = surfSalt
1487	0	if (CS%id_Usurf > 0) CS%Usurf(i,j) = surfU
1488	0	if (CS%id_Vsurf > 0) CS%Vsurf(i,j) = surfV
1489
1490		endif ; enddo
1491		enddo
1492
1493	0	call cpu_clock_end(id_clock_KPP_compute_BLD)
1494
1495		! send diagnostics to post_data
1496	0	if (CS%id_BulkRi > 0) call post_data(CS%id_BulkRi, CS%BulkRi, CS%diag)
1497	0	if (CS%id_N > 0) call post_data(CS%id_N, CS%N, CS%diag)
1498	0	if (CS%id_N2 > 0) call post_data(CS%id_N2, CS%N2, CS%diag)
1499	0	if (CS%id_Tsurf > 0) call post_data(CS%id_Tsurf, CS%Tsurf, CS%diag)
1500	0	if (CS%id_Ssurf > 0) call post_data(CS%id_Ssurf, CS%Ssurf, CS%diag)
1501	0	if (CS%id_Usurf > 0) call post_data(CS%id_Usurf, CS%Usurf, CS%diag)
1502	0	if (CS%id_Vsurf > 0) call post_data(CS%id_Vsurf, CS%Vsurf, CS%diag)
1503	0	if (CS%id_BulkDrho > 0) call post_data(CS%id_BulkDrho, CS%dRho, CS%diag)
1504	0	if (CS%id_BulkUz2 > 0) call post_data(CS%id_BulkUz2, CS%Uz2, CS%diag)
1505	0	if (CS%id_EnhK > 0) call post_data(CS%id_EnhK, CS%EnhK, CS%diag)
1506	0	if (CS%id_EnhVt2 > 0) call post_data(CS%id_EnhVt2, CS%EnhVt2, CS%diag)
1507	0	if (CS%id_La_SL > 0) call post_data(CS%id_La_SL, CS%La_SL, CS%diag)
1508	0	if (CS%id_Vt2 > 0) call post_data(CS%id_Vt2, CS%Vt2, CS%diag)
1509
1510	0	if (CS%StokesMOST) then
1511	0	if (CS%id_StokesXI > 0) call post_data(CS%id_StokesXI, CS%StokesParXI, CS%diag)
1512	0	if (CS%id_Lam2 > 0) call post_data(CS%id_Lam2 , CS%Lam2 , CS%diag)
1513		endif
1514
1515		! BLD smoothing:
1516	0	if (CS%n_smooth > 0) call KPP_smooth_BLD(CS, G, GV, US, dz)
1517
1518	0	end subroutine KPP_compute_BLD
1519
1520
1521		!> Apply a 1-1-4-1-1 Laplacian filter one time on BLD to reduce any horizontal two-grid-point noise
1522	0	subroutine KPP_smooth_BLD(CS, G, GV, US, dz)
1523		! Arguments
1524		type(KPP_CS), pointer :: CS !< Control structure
1525		type(ocean_grid_type), intent(inout) :: G !< Ocean grid
1526		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid
1527		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1528		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: dz !< Layer thicknesses [Z ~> m]
1529
1530		! local variables
1531	0	real, dimension(SZI_(G),SZJ_(G)) :: OBLdepth_prev ! OBLdepth before s.th smoothing iteration [Z ~> m]
1532	0	real, dimension(SZI_(G),SZJ_(G)) :: total_depth ! The total depth of the water column, adjusted
1533		! for the minimum layer thickness [Z ~> m]
1534	0	real, dimension( GV%ke ) :: cellHeight ! Cell center heights referenced to surface [Z ~> m]
1535		! (negative in the ocean)
1536	0	real, dimension( GV%ke+1 ) :: iFaceHeight ! Interface heights referenced to surface [Z ~> m]
1537		! (negative in the ocean)
1538		real :: wc, ww, we, wn, ws ! averaging weights for smoothing [nondim]
1539		real :: dh ! The local thickness used for calculating interface positions [Z ~> m]
1540	0	real :: h_cor(SZI_(G)) ! A cumulative correction arising from inflation of vanished layers [Z ~> m]
1541		real :: hcorr ! A cumulative correction arising from inflation of vanished layers [Z ~> m]
1542		integer :: i, j, k, s, halo
1543
1544	0	call cpu_clock_begin(id_clock_KPP_smoothing)
1545
1546		! Find the total water column thickness first, as it is reused for each smoothing pass.
1547	0	total_depth(:,:) = 0.0
1548
1549		!$OMP parallel do default(shared) private(dh, h_cor)
1550	0	do j = G%jsc, G%jec
1551	0	h_cor(:) = 0.
1552	0	do k=1,GV%ke
1553	0	do i=G%isc,G%iec ; if (G%mask2dT(i,j) > 0.0) then
1554		! This code replicates the interface height calculations below. It could be simpler, as shown below.
1555	0	dh = dz(i,j,k) ! Nominal thickness to use for increment
1556	0	dh = dh + h_cor(i) ! Take away the accumulated error (could temporarily make dh<0)
1557	0	h_cor(i) = min( dh - CS%min_thickness, 0. ) ! If inflating then hcorr<0
1558	0	dh = max( dh, CS%min_thickness ) ! Limit increment dh>=min_thickness
1559	0	total_depth(i,j) = total_depth(i,j) + dh
1560		endif ; enddo
1561		enddo
1562		enddo
1563		! A much simpler (but answer changing) version of the total_depth calculation would be
1564		! do k=1,GV%ke ; do j=G%jsc,G%jec ; do i=G%isc,G%iec
1565		! total_depth(i,j) = total_depth(i,j) + dz(i,j,k)
1566		! enddo ; enddo ; enddo
1567
1568		! Update halos once, then march inward for each iteration
1569	0	if (CS%n_smooth > 1) call pass_var(total_depth, G%Domain, halo=CS%n_smooth, complete=.false.)
1570	0	call pass_var(CS%OBLdepth, G%Domain, halo=CS%n_smooth)
1571
1572	0	if (CS%id_OBLdepth_original > 0) CS%OBLdepth_original(:,:) = CS%OBLdepth(:,:)
1573
1574	0	do s=1,CS%n_smooth
1575
1576	0	OBLdepth_prev(:,:) = CS%OBLdepth(:,:)
1577	0	halo = CS%n_smooth - s
1578
1579		! apply smoothing on OBL depth
1580		!$OMP parallel do default(none) shared(G, GV, CS, OBLdepth_prev, total_depth, halo) &
1581		!$OMP private(wc, ww, we, wn, ws)
1582	0	do j = G%jsc-halo, G%jec+halo
1583	0	do i = G%isc-halo, G%iec+halo ; if (G%mask2dT(i,j) > 0.0) then
1584		! compute weights
1585	0	ww = 0.125 * G%mask2dT(i-1,j)
1586	0	we = 0.125 * G%mask2dT(i+1,j)
1587	0	ws = 0.125 * G%mask2dT(i,j-1)
1588	0	wn = 0.125 * G%mask2dT(i,j+1)
1589	0	wc = 1.0 - (ww+we+wn+ws)
1590
1591	0	if (CS%answer_date < 20240501) then
1592		CS%OBLdepth(i,j) = wc * OBLdepth_prev(i,j) &
1593		+ ww * OBLdepth_prev(i-1,j) &
1594		+ we * OBLdepth_prev(i+1,j) &
1595		+ ws * OBLdepth_prev(i,j-1) &
1596	0	+ wn * OBLdepth_prev(i,j+1)
1597		else
1598		CS%OBLdepth(i,j) = wc * OBLdepth_prev(i,j) &
1599		+ ((ww * OBLdepth_prev(i-1,j) + we * OBLdepth_prev(i+1,j)) &
1600	0	+ (ws * OBLdepth_prev(i,j-1) + wn * OBLdepth_prev(i,j+1)))
1601		endif
1602
1603		! Apply OBLdepth smoothing at a cell only if the OBLdepth gets deeper via smoothing.
1604	0	if (CS%deepen_only) CS%OBLdepth(i,j) = max(CS%OBLdepth(i,j), OBLdepth_prev(i,j))
1605
1606		! prevent OBL depths deeper than the bathymetric depth
1607	0	CS%OBLdepth(i,j) = min( CS%OBLdepth(i,j), total_depth(i,j) ) ! no deeper than bottom
1608		endif ; enddo
1609		enddo
1610
1611		enddo ! s-loop
1612
1613		! Determine the fractional index of the bottom of the boundary layer.
1614		!$OMP parallel do default(none) shared(G, GV, CS, dz) &
1615		!$OMP private(dh, hcorr, cellHeight, iFaceHeight)
1616	0	do j=G%jsc,G%jec ; do i=G%isc,G%iec ; if (G%mask2dT(i,j) > 0.0) then
1617
1618	0	iFaceHeight(1) = 0.0 ! BBL is all relative to the surface
1619	0	hcorr = 0.
1620	0	do k=1,GV%ke
1621		! cell center and cell bottom in meters (negative values in the ocean)
1622	0	dh = dz(i,j,k) ! Nominal thickness to use for increment
1623	0	dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0)
1624	0	hcorr = min( dh - CS%min_thickness, 0. ) ! If inflating then hcorr<0
1625	0	dh = max( dh, CS%min_thickness ) ! Limit increment dh>=min_thickness
1626	0	cellHeight(k) = iFaceHeight(k) - 0.5 * dh
1627	0	iFaceHeight(k+1) = iFaceHeight(k) - dh
1628		enddo
1629
1630	0	CS%kOBL(i,j) = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, CS%OBLdepth(i,j) )
1631		endif ; enddo ; enddo
1632
1633	0	call cpu_clock_end(id_clock_KPP_smoothing)
1634
1635	0	end subroutine KPP_smooth_BLD
1636
1637
1638		!> Copies KPP surface boundary layer depth into BLD, in units of [Z ~> m] unless other units are specified.
1639	0	subroutine KPP_get_BLD(CS, BLD, G, US, m_to_BLD_units)
1640		type(KPP_CS), pointer :: CS !< Control structure for
1641		!! this module
1642		type(ocean_grid_type), intent(in) :: G !< Grid structure
1643		type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
1644		real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: BLD !< Boundary layer depth [Z ~> m] or other units
1645		real, optional, intent(in) :: m_to_BLD_units !< A conversion factor from meters
1646		!! to the desired units for BLD [various]
1647		! Local variables
1648		real :: scale ! A dimensional rescaling factor in [nondim] or other units.
1649		integer :: i,j
1650
1651	0	scale = 1.0 ; if (present(m_to_BLD_units)) scale = US%Z_to_m*m_to_BLD_units
1652
1653		!$OMP parallel do default(none) shared(BLD, CS, G, scale)
1654	0	do j = G%jsc, G%jec ; do i = G%isc, G%iec
1655	0	BLD(i,j) = scale * CS%OBLdepth(i,j)
1656		enddo ; enddo
1657
1658	0	end subroutine KPP_get_BLD
1659
1660		!> Apply KPP non-local transport of surface fluxes for a given tracer
1661	0	subroutine KPP_NonLocalTransport(CS, G, GV, h, nonLocalTrans, surfFlux, &
1662	0	dt, diag, tr_ptr, scalar, flux_scale)
1663		type(KPP_CS), intent(in) :: CS !< Control structure
1664		type(ocean_grid_type), intent(in) :: G !< Ocean grid
1665		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid
1666		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer/level thickness [H ~> m or kg m-2]
1667		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: nonLocalTrans !< Non-local transport [nondim]
1668		real, dimension(SZI_(G),SZJ_(G)), intent(in) :: surfFlux !< Surface flux of scalar
1669		!! [conc H T-1 ~> conc m s-1 or conc kg m-2 s-1]
1670		real, intent(in) :: dt !< Time-step [T ~> s]
1671		type(diag_ctrl), target, intent(in) :: diag !< Diagnostics
1672		type(tracer_type), pointer, intent(in) :: tr_ptr !< tracer_type has diagnostic ids on it
1673		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: scalar !< Scalar (scalar units [conc])
1674		real, optional, intent(in) :: flux_scale !< Scale factor to get surfFlux
1675		!! into proper units [various]
1676
1677		integer :: i, j, k
1678	0	real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dtracer ! Rate of tracer change [conc T-1 ~> conc s-1]
1679	0	real, dimension(SZI_(G),SZJ_(G)) :: surfFlux_loc ! An optionally rescaled surface flux of the scalar
1680		! in [conc H T-1 ~> conc m s-1 or conc kg m-2 s-1] or other units
1681
1682		! term used to scale
1683	0	if (present(flux_scale)) then
1684	0	do j = G%jsc, G%jec ; do i = G%isc, G%iec
1685	0	surfFlux_loc(i,j) = surfFlux(i,j) * flux_scale
1686		enddo ; enddo
1687		else
1688	0	surfFlux_loc(:,:) = surfFlux(:,:)
1689		endif
1690
1691		! Post surface flux diagnostic
1692	0	if (tr_ptr%id_net_surfflux > 0) call post_data(tr_ptr%id_net_surfflux, surfFlux_loc(:,:), diag)
1693
1694		! Only continue if we are applying the nonlocal tendency
1695		! or the nonlocal tendency diagnostic has been requested
1696	0	if ((tr_ptr%id_NLT_tendency > 0) .or. (CS%applyNonLocalTrans)) then
1697
1698		!$OMP parallel do default(none) shared(dtracer, nonLocalTrans, h, G, GV, surfFlux_loc)
1699	0	do k = 1, GV%ke ; do j = G%jsc, G%jec ; do i = G%isc, G%iec
1700		dtracer(i,j,k) = ( nonLocalTrans(i,j,k) - nonLocalTrans(i,j,k+1) ) / &
1701	0	( h(i,j,k) + GV%H_subroundoff ) * surfFlux_loc(i,j)
1702		enddo ; enddo ; enddo
1703
1704		! Update tracer due to non-local redistribution of surface flux
1705	0	if (CS%applyNonLocalTrans) then
1706		!$OMP parallel do default(none) shared(G, GV, dt, scalar, dtracer)
1707	0	do k = 1, GV%ke ; do j = G%jsc, G%jec ; do i = G%isc, G%iec
1708	0	scalar(i,j,k) = scalar(i,j,k) + dt * dtracer(i,j,k)
1709		enddo ; enddo ; enddo
1710		endif
1711	0	if (tr_ptr%id_NLT_tendency > 0) call post_data(tr_ptr%id_NLT_tendency, dtracer, diag)
1712
1713		endif
1714
1715
1716	0	if (tr_ptr%id_NLT_budget > 0) then
1717		!$OMP parallel do default(none) shared(G, GV, dtracer, nonLocalTrans, surfFlux_loc)
1718	0	do k = 1, GV%ke ; do j = G%jsc, G%jec ; do i = G%isc, G%iec
1719		! Here dtracer has units of [Q R Z T-1 ~> W m-2].
1720	0	dtracer(i,j,k) = (nonLocalTrans(i,j,k) - nonLocalTrans(i,j,k+1)) * surfFlux_loc(i,j)
1721		enddo ; enddo ; enddo
1722	0	call post_data(tr_ptr%id_NLT_budget, dtracer(:,:,:), diag)
1723		endif
1724
1725	0	end subroutine KPP_NonLocalTransport
1726
1727
1728		!> Apply KPP non-local transport of surface fluxes for temperature.
1729	0	subroutine KPP_NonLocalTransport_temp(CS, G, GV, h, nonLocalTrans, surfFlux, dt, tr_ptr, scalar, C_p)
1730		type(KPP_CS), intent(in) :: CS !< Control structure
1731		type(ocean_grid_type), intent(in) :: G !< Ocean grid
1732		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid
1733		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer/level thickness [H ~> m or kg m-2]
1734		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: nonLocalTrans !< Non-local transport [nondim]
1735		real, dimension(SZI_(G),SZJ_(G)), intent(in) :: surfFlux !< Surface flux of temperature
1736		!! [C H T-1 ~> degC m s-1 or degC kg m-2 s-1]
1737		real, intent(in) :: dt !< Time-step [T ~> s]
1738		type(tracer_type), pointer, intent(in) :: tr_ptr !< tracer_type has diagnostic ids on it
1739		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: scalar !< temperature [C ~> degC]
1740		real, intent(in) :: C_p !< Seawater specific heat capacity
1741		!! [Q C-1 ~> J kg-1 degC-1]
1742
1743		call KPP_NonLocalTransport(CS, G, GV, h, nonLocalTrans, surfFlux, dt, CS%diag, &
1744	0	tr_ptr, scalar)
1745
1746	0	end subroutine KPP_NonLocalTransport_temp
1747
1748
1749		!> Apply KPP non-local transport of surface fluxes for salinity.
1750	0	subroutine KPP_NonLocalTransport_saln(CS, G, GV, h, nonLocalTrans, surfFlux, dt, tr_ptr, scalar)
1751		type(KPP_CS), intent(in) :: CS !< Control structure
1752		type(ocean_grid_type), intent(in) :: G !< Ocean grid
1753		type(verticalGrid_type), intent(in) :: GV !< Ocean vertical grid
1754		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer/level thickness [H ~> m or kg m-2]
1755		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(in) :: nonLocalTrans !< Non-local transport [nondim]
1756		real, dimension(SZI_(G),SZJ_(G)), intent(in) :: surfFlux !< Surface flux of salt
1757		!! [S H T-1 ~> ppt m s-1 or ppt kg m-2 s-1]
1758		real, intent(in) :: dt !< Time-step [T ~> s]
1759		real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: scalar !< Salinity [S ~> ppt]
1760		type(tracer_type), pointer, intent(in) :: tr_ptr !< tracer_type has diagnostic ids on it
1761
1762		call KPP_NonLocalTransport(CS, G, GV, h, nonLocalTrans, surfFlux, dt, CS%diag, &
1763	0	tr_ptr, scalar)
1764
1765	0	end subroutine KPP_NonLocalTransport_saln
1766
1767		!> Compute Stokes Drift components at zbot < ztop <= 0 and at k=0.5*(ztop+zbot) and
1768		!! average components from ztop to zbot <= 0
1769	0	subroutine Compute_StokesDrift(i ,j, ztop, zbot, uS_i, vS_i, uS_k, vS_k, uSbar, vSbar, waves)
1770
1771		type(wave_parameters_CS), pointer :: waves !< Wave CS for Langmuir turbulence
1772		real, intent(in) :: ztop !< cell top
1773		real, intent(in) :: zbot !< cell bottom
1774		real, intent(inout) :: uS_i !< Stokes u velocity at zbot interface
1775		real, intent(inout) :: vS_i !< Stokes v velocity at zbot interface
1776		real, intent(inout) :: uS_k !< Stokes u velocity at zk center
1777		real, intent(inout) :: vS_k !< Stokes v at zk =0.5(ztop+zbot)
1778		real, intent(inout) :: uSbar !< mean Stokes u (ztop to zbot)
1779		real, intent(inout) :: vSbar !< mean Stokes v (ztop to zbot)
1780		integer, intent(in) :: i !< Meridional index of H-point
1781		integer, intent(in) :: j !< Zonal index of H-point
1782
1783		! local variables
1784		integer :: b !< wavenumber band index
1785		real :: fexp !< an exponential function
1786		real :: WaveNum !< Wavenumber
1787
1788	0	uS_i = 0.0
1789	0	vS_i = 0.0
1790	0	uS_k = 0.0
1791	0	vS_k = 0.0
1792	0	uSbar = 0.0
1793	0	vSbar = 0.0
1794	0	do b = 1, waves%NumBands
1795	0	WaveNum = waves%WaveNum_Cen(b)
1796	0	fexp = exp(2. * WaveNum * zbot)
1797	0	uS_i = uS_i + waves%Ustk_Hb(i,j,b) * fexp
1798	0	vS_i = vS_i + waves%Vstk_Hb(i,j,b) * fexp
1799	0	fexp = exp( WaveNum * (ztop + zbot) )
1800	0	uS_k = uS_k+ waves%Ustk_Hb(i,j,b) * fexp
1801	0	vS_k = vS_k+ waves%Vstk_Hb(i,j,b) * fexp
1802	0	fexp = exp(2. * WaveNum * ztop) - exp(2. * WaveNum * zbot)
1803	0	uSbar = uSbar + 0.5 * waves%Ustk_Hb(i,j,b) * fexp / WaveNum
1804	0	vSbar = vSbar + 0.5 * waves%Vstk_Hb(i,j,b) * fexp / WaveNum
1805		enddo
1806	0	uSbar = uSbar / (ztop-zbot)
1807	0	vSbar = vSbar / (ztop-zbot)
1808
1809	0	end subroutine Compute_StokesDrift
1810
1811		!> Clear pointers, deallocate memory
1812	0	subroutine KPP_end(CS)
1813		type(KPP_CS), pointer :: CS !< Control structure
1814
1815	0	if (.not.associated(CS)) return
1816
1817	0	deallocate(CS)
1818
1819		end subroutine KPP_end
1820
1821	0	end module MOM_CVMix_KPP