MOM_CVMix_KPP.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Provides the K-Profile Parameterization (KPP) of Large et al., 1994, via CVMix.
6module mom_cvmix_kpp
7
8use mom_coms, only : max_across_pes
9use mom_debugging, only : hchksum, is_nan
10use mom_diag_mediator, only : time_type, diag_ctrl, safe_alloc_ptr, post_data
12use mom_error_handler, only : mom_error, mom_mesg, fatal, warning, is_root_pe
13use mom_eos, only : eos_type, calculate_density
14use mom_file_parser, only : get_param, log_param, log_version, param_file_type
17use mom_interface_heights, only : thickness_to_dz
18use mom_restart, only : mom_restart_cs, register_restart_field
23use mom_domains, only : pass_var
24use mom_cpu_clock, only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
25use mom_cpu_clock, only : clock_module, clock_routine
27
28use cvmix_kpp, only : cvmix_init_kpp, cvmix_put_kpp, cvmix_get_kpp_real
29use cvmix_kpp, only : cvmix_coeffs_kpp
30use cvmix_kpp, only : cvmix_kpp_compute_obl_depth
31use cvmix_kpp, only : cvmix_kpp_compute_turbulent_scales
32use cvmix_kpp, only : cvmix_kpp_compute_bulk_richardson
33use cvmix_kpp, only : cvmix_kpp_compute_unresolved_shear
34use cvmix_kpp, only : cvmix_kpp_params_type
35use cvmix_kpp, only : cvmix_kpp_compute_kobl_depth
36use cvmix_kpp, only : cvmix_kpp_compute_stokesxi
37
38implicit none ; private
39
40#include "MOM_memory.h"
41
43public :: kpp_init
44public :: kpp_compute_bld
45public :: kpp_calculate
46public :: kpp_end
50public :: kpp_get_bld
51
52! Enumerated constants
53integer, private, parameter :: nlt_shape_cvmix = 0 !< Use the CVMix profile
54integer, private, parameter :: nlt_shape_linear = 1 !< Linear, \f$ G(\sigma) = 1-\sigma \f$
55integer, private, parameter :: nlt_shape_parabolic = 2 !< Parabolic, \f$ G(\sigma) = (1-\sigma)^2 \f$
56integer, private, parameter :: nlt_shape_cubic = 3 !< Cubic, \f$ G(\sigma) = 1 + (2\sigma-3) \sigma^2\f$
57integer, private, parameter :: nlt_shape_cubic_lmd = 4 !< Original shape,
58 !! \f$ G(\sigma) = \frac{27}{4} \sigma (1-\sigma)^2 \f$
59
60integer, private, parameter :: sw_method_all_sw = 0 !< Use all shortwave radiation
61integer, private, parameter :: sw_method_mxl_sw = 1 !< Use shortwave radiation absorbed in mixing layer
62integer, private, parameter :: sw_method_lv1_sw = 2 !< Use shortwave radiation absorbed in layer 1
63integer, 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
79type, 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
186end type kpp_cs
187
188!>@{ CPU time clocks
189integer :: id_clock_kpp_calc, id_clock_kpp_compute_bld, id_clock_kpp_smoothing
190!>@}
191
192#define __DO_SAFETY_CHECKS__
193
194contains
195
196!> Routine to register restarts, pass-through to children modules
197subroutine 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 if (associated(cs)) call mom_error(fatal, 'MOM_CVMix_KPP, register_KPP_restarts: '// &
207 'Control structure has already been initialized')
208 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 if (.not. use_kpp) return
211 allocate(cs)
212
213 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 default=.false., do_not_log=.true.)
219 if (fpmix) call register_restart_field(cs%OBLdepth, 'KPP_OBLdepth', .false., restart_csp)
220
221end 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.
225logical 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 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 '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 default=.false.)
256 ! Forego remainder of initialization if not using this scheme
257 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 default=99991231, do_not_log=.true.)
262
263 call openparameterblock(paramfile,'KPP')
264 call get_param(paramfile, mdl, 'PASSIVE', cs%passiveMode, &
265 'If True, puts KPP into a passive-diagnostic mode.', &
266 default=.false.)
267 !BGR: Note using PASSIVE for KPP creates warning for PASSIVE from Convection
268 ! should we create a separate flag?
269 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 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 default=0)
279 if (cs%n_smooth > g%domain%nihalo) then
280 call mom_error(fatal,'KPP smoothing number (N_SMOOTH) cannot be greater than NIHALO.')
281 elseif (cs%n_smooth > g%domain%njhalo) then
282 call mom_error(fatal,'KPP smoothing number (N_SMOOTH) cannot be greater than NJHALO.')
283 endif
284 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 default=.false.)
289 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 units='nondim', default=0.3)
295 call get_param(paramfile, mdl, 'VON_KARMAN', cs%vonKarman, &
296 'von Karman constant.', &
297 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 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 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 default='LMD94')
309 call get_param(paramfile, mdl, 'STOKES_MOST', cs%StokesMOST, &
310 'If True, use Stokes Similarity package.', &
311 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 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 default=.false.)
319 call get_param(paramfile, mdl, 'CS', cs%cs, &
320 'Parameter for computing velocity scale function.', &
321 units='nondim', default=98.96)
322 call get_param(paramfile, mdl, 'CS2', cs%cs2, &
323 'Parameter for computing non-local term.', &
324 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 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 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 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 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 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 = const1*d*N*ws, with d=1m, N=1e-5/s, ws=1e-6 m/s.', &
349 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 default='CVMix')
360 select case ( trim(string) )
361 case ("CVMix") ; cs%NLT_shape = nlt_shape_cvmix
362 case ("LINEAR") ; cs%NLT_shape = nlt_shape_linear
363 case ("PARABOLIC") ; cs%NLT_shape = nlt_shape_parabolic
364 case ("CUBIC") ; cs%NLT_shape = nlt_shape_cubic
365 case ("CUBIC_LMD") ; cs%NLT_shape = nlt_shape_cubic_lmd
366 case default ; call mom_error(fatal,"KPP_init: "// &
367 "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 default='SimpleShapes')
377 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 cs_is_one=.true.
381 endif
382 if (cs%MatchTechnique == 'ParabolicNonLocal' .or. cs%MatchTechnique == 'SimpleShapes') then
383 ! if gradient won't be matched, lnoDGat1=.true.
384 lnodgat1=.true.
385 endif
386
387 ! safety check to avoid negative diff/visc
388 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 "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 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 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 default='MXL_SW')
409 select case ( trim(string) )
410 case ("ALL_SW") ; cs%SW_METHOD = sw_method_all_sw
411 case ("MXL_SW") ; cs%SW_METHOD = sw_method_mxl_sw
412 case ("LV1_SW") ; cs%SW_METHOD = sw_method_lv1_sw
413 case default ; call mom_error(fatal,"KPP_init: "// &
414 "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 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 'mixing coefficient.', default=.false.)
426 call get_param(paramfile, mdl, "STOKES_MIXING", cs%Stokes_Mixing, &
427 'Flag for Langmuir turbulence enhancement of turbulent '//&
428 'mixing coefficient.', default=.false.)
429 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 default='CONSTANT')
436 select case ( trim(string))
437 case ("CONSTANT") ; cs%LT_K_SHAPE = lt_k_constant
438 case ("SCALED") ; cs%LT_K_SHAPE = lt_k_scaled
439 case default ; call mom_error(fatal,"KPP_init: "//&
440 "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 default='CONSTANT')
452 select case ( trim(string))
453 case ("CONSTANT")
454 cs%LT_K_METHOD = lt_k_mode_constant
455 langmuir_mixing_opt = 'LWF16'
456 case ("VR12")
457 cs%LT_K_METHOD = lt_k_mode_vr12
458 langmuir_mixing_opt = 'LWF16'
459 case ("RW16")
460 cs%LT_K_METHOD = lt_k_mode_rw16
461 langmuir_mixing_opt = 'RWHGK16'
462 case default
463 call mom_error(fatal,"KPP_init: "//&
464 "Unrecognized KPP_LT_K_METHOD option: "//trim(string))
465 end select
466 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 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 'in Bulk Richardson Number.', default=.false.)
476 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 default='CONSTANT')
489 select case ( trim(string))
490 case ("CONSTANT")
491 cs%LT_VT2_METHOD = lt_vt2_mode_constant
492 langmuir_entrainment_opt = 'LWF16'
493 case ("VR12")
494 cs%LT_VT2_METHOD = lt_vt2_mode_vr12
495 langmuir_entrainment_opt = 'LWF16'
496 case ("RW16")
497 cs%LT_VT2_METHOD = lt_vt2_mode_rw16
498 langmuir_entrainment_opt = 'RWHGK16'
499 case ("LF17")
500 cs%LT_VT2_METHOD = lt_vt2_mode_lf17
501 langmuir_entrainment_opt = 'LF17'
502 case default
503 call mom_error(fatal,"KPP_init: "//&
504 "Unrecognized KPP_LT_VT2_METHOD option: "//trim(string))
505 end select
506 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 units="nondim", default=1.0)
510 endif
511 endif
512
513 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 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 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 default=default_answer_date)
529
530 call closeparameterblock(paramfile)
531
532 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**2*cs%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 cvmix_kpp_params_user=cs%KPP_params )
552
553 ! Register diagnostics
554 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 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 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 cmor_units='m', cmor_standard_name='Ocean Mixed Layer Thickness Defined by Mixing Scheme')
569 endif
570 if ( cs%StokesMOST ) then
571 cs%id_StokesXI = register_diag_field('ocean_model', 'StokesXI', diag%axesT1, time, &
572 'Stokes Similarity Parameter', 'nondim')
573 cs%id_Lam2 = register_diag_field('ocean_model', 'Lam2', diag%axesT1, time, &
574 '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 '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 '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 '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 '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 '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 '(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 '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 'Unresolved shear turbulence used by [CVMix] KPP', 'm2/s2', conversion=us%Z_to_m**2*us%s_to_T**2)
595 cs%id_uStar = register_diag_field('ocean_model', 'KPP_uStar', diag%axesT1, time, &
596 'Friction velocity, u*, as used by [CVMix] KPP', 'm/s', conversion=us%Z_to_m*us%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 'm2/s3', conversion=us%L_to_m**2*us%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 '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 '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 '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 '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 '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 '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 '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 '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 '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 'm/s', conversion=us%L_T_to_m_s)
627 cs%id_EnhK = register_diag_field('ocean_model', 'EnhK', diag%axesTI, time, &
628 '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 '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 'Surface-layer Langmuir number computed in [CVMix] KPP','nondim')
633
634 allocate( cs%N( szi_(g), szj_(g), szk_(gv)+1 ), source=0. )
635 allocate( cs%StokesParXI( szi_(g), szj_(g) ), source=0. )
636 allocate( cs%Lam2 ( szi_(g), szj_(g) ), source=0. )
637 allocate( cs%kOBL( szi_(g), szj_(g) ), source=0. )
638 allocate( cs%La_SL( szi_(g), szj_(g) ), source=0. )
639 allocate( cs%Vt2( szi_(g), szj_(g), szk_(gv) ), source=0. )
640 if (cs%id_OBLdepth_original > 0) allocate( cs%OBLdepth_original( szi_(g), szj_(g) ) )
641
642 allocate( cs%OBLdepthprev( szi_(g), szj_(g) ), source=0.0 )
643 if (cs%id_BulkDrho > 0) allocate( cs%dRho( szi_(g), szj_(g), szk_(gv) ), source=0. )
644 if (cs%id_BulkUz2 > 0) allocate( cs%Uz2( szi_(g), szj_(g), szk_(gv) ), source=0. )
645 if (cs%id_BulkRi > 0) allocate( cs%BulkRi( szi_(g), szj_(g), szk_(gv) ), source=0. )
646 if (cs%id_Sigma > 0) allocate( cs%sigma( szi_(g), szj_(g), szk_(gv)+1 ), source=0. )
647 if (cs%id_Ws > 0) allocate( cs%Ws( szi_(g), szj_(g), szk_(gv) ), source=0. )
648 if (cs%id_N2 > 0) allocate( cs%N2( szi_(g), szj_(g), szk_(gv)+1 ), source=0. )
649 if (cs%id_Kt_KPP > 0) allocate( cs%Kt_KPP( szi_(g), szj_(g), szk_(gv)+1 ), source=0. )
650 if (cs%id_Ks_KPP > 0) allocate( cs%Ks_KPP( szi_(g), szj_(g), szk_(gv)+1 ), source=0. )
651 if (cs%id_Kv_KPP > 0) allocate( cs%Kv_KPP( szi_(g), szj_(g), szk_(gv)+1 ), source=0. )
652 if (cs%id_Tsurf > 0) allocate( cs%Tsurf( szi_(g), szj_(g) ), source=0. )
653 if (cs%id_Ssurf > 0) allocate( cs%Ssurf( szi_(g), szj_(g) ), source=0. )
654 if (cs%id_Usurf > 0) allocate( cs%Usurf( szib_(g), szj_(g) ), source=0. )
655 if (cs%id_Vsurf > 0) allocate( cs%Vsurf( szi_(g), szjb_(g) ), source=0. )
656 if (cs%id_EnhVt2 > 0) allocate( cs%EnhVt2( szi_(g), szj_(g), szk_(gv) ), source=0. )
657 if (cs%id_EnhK > 0) allocate( cs%EnhK( szi_(g), szj_(g), szk_(gv)+1 ), source=0. )
658
659 id_clock_kpp_calc = cpu_clock_id('Ocean KPP calculate)', grain=clock_module)
660 id_clock_kpp_compute_bld = cpu_clock_id('(Ocean KPP comp BLD)', grain=clock_routine)
661
662end function kpp_init
663
664!> KPP vertical diffusivity/viscosity and non-local tracer transport
665subroutine kpp_calculate(CS, G, GV, US, h, tv, uStar, buoyFlux, Kt, Ks, Kv, &
666 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 real, dimension(SZI_(G),SZK_(GV)) :: dz ! Height change across layers [Z ~> m]
694 real, dimension( GV%ke ) :: cellheight ! Cell center heights referenced to surface [Z ~> m] (negative in ocean)
695 real, dimension( GV%ke+1 ) :: ifaceheight ! Interface heights referenced to surface [Z ~> m] (negative in ocean)
696 real, dimension( GV%ke ) :: z_cell ! Cell center heights referenced to surface [m] (negative in ocean)
697 real, dimension( GV%ke+1 ) :: z_inter ! Cell interface heights referenced to surface [m] (negative in ocean)
698 real, dimension( GV%ke+1, 2) :: kdiffusivity ! Vertical diffusivity at interfaces in MKS units [m2 s-1]
699 real, dimension( GV%ke+1 ) :: kviscosity ! Vertical viscosity at interfaces in MKS units [m2 s-1]
700 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 if (cs%Stokes_Mixing .and. .not.associated(waves)) call mom_error(fatal, &
714 "KPP_calculate: The Waves control structure must be associated if STOKES_MIXING is True.")
715
716 if (cs%debug) then
717 call hchksum(h, "KPP in: h", g%HI, haloshift=0, unscale=gv%H_to_m)
718 call hchksum(ustar, "KPP in: uStar", g%HI, haloshift=0, unscale=us%Z_to_m*us%s_to_T)
719 call hchksum(buoyflux, "KPP in: buoyFlux", g%HI, haloshift=0, unscale=us%L_to_m**2*us%s_to_T**3)
720 call hchksum(kt, "KPP in: Kt", g%HI, haloshift=0, unscale=gv%HZ_T_to_m2_s)
721 call hchksum(ks, "KPP in: Ks", g%HI, haloshift=0, unscale=gv%HZ_T_to_m2_s)
722 endif
723
724 nonlocaltrans(:,:) = 0.0
725
726 if (cs%id_Kd_in > 0) call post_data(cs%id_Kd_in, kt, cs%diag)
727
728 call cpu_clock_begin(id_clock_kpp_calc)
729 buoy_scale = us%L_to_m**2*us%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 do j = g%jsc, g%jec
739
740 ! Find the vertical distances across layers.
741 call thickness_to_dz(h, tv, dz, j, g, gv)
742
743 do i = g%isc, g%iec ; if (g%mask2dT(i,j) > 0.0) then
744
745 ! things independent of position within the column
746 surffricvel = us%Z_to_m*us%s_to_T * ustar(i,j)
747
748 ifaceheight(1) = 0.0 ! BBL is all relative to the surface
749 hcorr = 0.
750 do k=1,gv%ke
751
752 ! cell center and cell bottom in meters (negative values in the ocean)
753 dh = dz(i,k) ! Nominal thickness to use for increment
754 dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0)
755 hcorr = min( dh - cs%min_thickness, 0. ) ! If inflating then hcorr<0
756 dh = max( dh, cs%min_thickness ) ! Limit increment dh>=min_thickness
757 cellheight(k) = ifaceheight(k) - 0.5 * dh
758 ifaceheight(k+1) = ifaceheight(k) - dh
759
760 enddo ! k-loop finishes
761
762 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 if (cs%SW_METHOD == sw_method_all_sw) then
772 surfbuoyflux = buoy_scale * buoyflux(i,j,1)
773 elseif (cs%SW_METHOD == sw_method_mxl_sw) then
774 ! We know the actual buoyancy flux into the OBL
775 surfbuoyflux = buoy_scale * (buoyflux(i,j,1) - buoyflux(i,j,int(cs%kOBL(i,j))+1))
776 elseif (cs%SW_METHOD == sw_method_lv1_sw) then
777 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 if (.not. (cs%MatchTechnique == 'MatchBoth')) then
782 kdiffusivity(:,:) = 0. ! Diffusivities for heat and salt [m2 s-1]
783 kviscosity(:) = 0. ! Viscosity [m2 s-1]
784 else
785 kdiffusivity(:,1) = gv%HZ_T_to_m2_s * kt(i,j,:)
786 kdiffusivity(:,2) = gv%HZ_T_to_m2_s * ks(i,j,:)
787 kviscosity(:) = gv%HZ_T_to_m2_s * kv(i,j,:)
788 endif
789
790 IF (cs%LT_K_ENHANCEMENT) then
791 if (cs%LT_K_METHOD==lt_k_mode_constant) then
792 langenhk = cs%KPP_K_ENH_FAC
793 elseif (cs%LT_K_METHOD==lt_k_mode_vr12) then
794 if (present(lamult)) then
795 langenhk = lamult(i,j)
796 else
797 langenhk = sqrt(1.+(1.5*cs%La_SL(i,j))**(-2) + &
798 (5.4*cs%La_SL(i,j))**(-4))
799 endif
800 elseif (cs%LT_K_METHOD==lt_k_mode_rw16) then
801 !This maximum value is proposed in Reichl et al., 2016 JPO formula
802 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 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 do k=1,gv%ke
813 if (cs%LT_K_SHAPE== lt_k_constant) then
814 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 elseif (cs%LT_K_SHAPE == lt_k_scaled) then
819 sigma = min(1.0,-ifaceheight(k)/cs%OBLdepth(i,j))
820 sigmaratio = sigma * (1. - sigma)**2 / 0.148148037
821 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 do k = 1, gv%ke
834 z_cell(k) = us%Z_to_m*cellheight(k)
835 enddo
836 do k = 1, gv%ke+1
837 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 cvmix_kpp_params_user=cs%KPP_params )
859
860 ! safety check, Kviscosity and Kdiffusivity must be >= 0
861 do k=1, gv%ke+1
862 if (kviscosity(k) < 0. .or. kdiffusivity(k,1) < 0.) then
863 write(*,'(a,3i3)') 'interface, i, j, k = ',j, j, k
864 write(*,'(a,2f12.5)') 'lon,lat=', g%geoLonT(i,j), g%geoLatT(i,j)
865 write(*,'(a,es12.4)') 'depth, z_inter(k) =',z_inter(k)
866 write(*,'(a,es12.4)') 'Kviscosity(k) =',kviscosity(k)
867 write(*,'(a,es12.4)') 'Kdiffusivity(k,1) =',kdiffusivity(k,1)
868 write(*,'(a,es12.4)') 'Kdiffusivity(k,2) =',kdiffusivity(k,2)
869 write(*,'(a,es12.4)') 'OBLdepth =',us%Z_to_m*cs%OBLdepth(i,j)
870 write(*,'(a,f8.4)') 'kOBL =',cs%kOBL(i,j)
871 write(*,'(a,es12.4)') 'u* =',surffricvel
872 write(*,'(a,es12.4)') 'bottom, z_inter(GV%ke+1) =',z_inter(gv%ke+1)
873 write(*,'(a,es12.4)') 'CS%La_SL(i,j) =',cs%La_SL(i,j)
874 write(*,'(a,es12.4)') 'LangEnhK =',langenhk
875 if (present(lamult)) write(*,'(a,es12.4)') 'lamult(i,j) =',lamult(i,j)
876 write(*,*) 'Kviscosity(:) =',kviscosity(:)
877 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 "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 if (surfbuoyflux < 0.0) then
896 if (cs%NLT_shape == nlt_shape_cubic) then
897 do k = 2, gv%ke
898 sigma = min(1.0,-ifaceheight(k)/cs%OBLdepth(i,j))
899 nonlocaltrans(k,1) = (1.0 - sigma)**2 * (1.0 + 2.0*sigma) !*
900 nonlocaltrans(k,2) = nonlocaltrans(k,1)
901 enddo
902 elseif (cs%NLT_shape == nlt_shape_parabolic) then
903 do k = 2, gv%ke
904 sigma = min(1.0,-ifaceheight(k)/cs%OBLdepth(i,j))
905 nonlocaltrans(k,1) = (1.0 - sigma)**2 !*CS%CS2
906 nonlocaltrans(k,2) = nonlocaltrans(k,1)
907 enddo
908 elseif (cs%NLT_shape == nlt_shape_linear) then
909 do k = 2, gv%ke
910 sigma = min(1.0,-ifaceheight(k)/cs%OBLdepth(i,j))
911 nonlocaltrans(k,1) = (1.0 - sigma)!*CS%CS2
912 nonlocaltrans(k,2) = nonlocaltrans(k,1)
913 enddo
914 elseif (cs%NLT_shape == nlt_shape_cubic_lmd) then
915 ! Sanity check (should agree with CVMix result using simple matching)
916 do k = 2, gv%ke
917 sigma = min(1.0,-ifaceheight(k)/cs%OBLdepth(i,j))
918 nonlocaltrans(k,1) = cs%CS2 * sigma*(1.0 -sigma)**2
919 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 nonlocaltransheat(i,j,:) = nonlocaltrans(:,1) ! temperature
927 nonlocaltransscalar(i,j,:) = nonlocaltrans(:,2) ! salinity
928
929 ! set the KPP diffusivity and viscosity to zero for testing purposes
930 if (cs%KPPzeroDiffusivity) then
931 kdiffusivity(:,1) = 0.0
932 kdiffusivity(:,2) = 0.0
933 kviscosity(:) = 0.0
934 endif
935
936 ! Copy 1d data into 3d diagnostic arrays
937 !/ grabbing obldepth_0d for next time step.
938 cs%OBLdepthprev(i,j) = cs%OBLdepth(i,j)
939 if (cs%id_sigma > 0) then
940 cs%sigma(i,j,:) = 0.
941 if (cs%OBLdepth(i,j)>0.) cs%sigma(i,j,:) = -ifaceheight(:)/cs%OBLdepth(i,j)
942 endif
943 if (cs%id_Kt_KPP > 0) cs%Kt_KPP(i,j,:) = us%m2_s_to_Z2_T * kdiffusivity(:,1)
944 if (cs%id_Ks_KPP > 0) cs%Ks_KPP(i,j,:) = us%m2_s_to_Z2_T * kdiffusivity(:,2)
945 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 if (.not. cs%passiveMode) then
949 if (cs%KPPisAdditive) then
950 do k=1, gv%ke+1
951 kt(i,j,k) = kt(i,j,k) + gv%m2_s_to_HZ_T * kdiffusivity(k,1)
952 ks(i,j,k) = ks(i,j,k) + gv%m2_s_to_HZ_T * kdiffusivity(k,2)
953 kv(i,j,k) = kv(i,j,k) + gv%m2_s_to_HZ_T * kviscosity(k)
954 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 do k=1, gv%ke+1
958 if (kdiffusivity(k,1) /= 0.) kt(i,j,k) = gv%m2_s_to_HZ_T * kdiffusivity(k,1)
959 if (kdiffusivity(k,2) /= 0.) ks(i,j,k) = gv%m2_s_to_HZ_T * kdiffusivity(k,2)
960 if (kviscosity(k) /= 0.) kv(i,j,k) = gv%m2_s_to_HZ_T * kviscosity(k)
961 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 call cpu_clock_end(id_clock_kpp_calc)
972
973 if (cs%debug) then
974 call hchksum(kt, "KPP out: Kt", g%HI, haloshift=0, unscale=gv%HZ_T_to_m2_s)
975 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 if (cs%id_OBLdepth > 0) call post_data(cs%id_OBLdepth, cs%OBLdepth, cs%diag)
980 if (cs%id_OBLdepth_original > 0) call post_data(cs%id_OBLdepth_original,cs%OBLdepth_original,cs%diag)
981 if (cs%id_sigma > 0) call post_data(cs%id_sigma, cs%sigma, cs%diag)
982 if (cs%id_Ws > 0) call post_data(cs%id_Ws, cs%Ws, cs%diag)
983 if (cs%id_uStar > 0) call post_data(cs%id_uStar, ustar, cs%diag)
984 if (cs%id_buoyFlux > 0) call post_data(cs%id_buoyFlux, buoyflux, cs%diag)
985 if (cs%id_Kt_KPP > 0) call post_data(cs%id_Kt_KPP, cs%Kt_KPP, cs%diag)
986 if (cs%id_Ks_KPP > 0) call post_data(cs%id_Ks_KPP, cs%Ks_KPP, cs%diag)
987 if (cs%id_Kv_KPP > 0) call post_data(cs%id_Kv_KPP, cs%Kv_KPP, cs%diag)
988 if (cs%id_NLTt > 0) call post_data(cs%id_NLTt, nonlocaltransheat, cs%diag)
989 if (cs%id_NLTs > 0) call post_data(cs%id_NLTs, nonlocaltransscalar,cs%diag)
990
991
992end subroutine kpp_calculate
993
994
995!> Compute OBL depth
996subroutine 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 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 real, dimension( GV%ke ) :: ws_1d ! Profile of vertical velocity scale for scalars in MKS units [m s-1]
1018 real, dimension( GV%ke ) :: deltarho ! delta Rho in numerator of Bulk Ri number [R ~> kg m-3]
1019 real, dimension( GV%ke ) :: deltabuoy ! Change in Buoyancy based on deltaRho [m s-2]
1020 real, dimension( GV%ke ) :: deltau2 ! square of delta U (shear) in denominator of Bulk Ri [m2 s-2]
1021 real, dimension( GV%ke ) :: surfbuoyflux2 ! Surface buoyancy flux in MKS units [m2 s-3]
1022 real, dimension( GV%ke ) :: bulkri_1d ! Bulk Richardson number for each layer [nondim]
1023 real, dimension( GV%ke ) :: vt2_1d ! Unresolved squared turbulence velocity for bulk Ri [m2 s-2]
1024 real, dimension( GV%ke ) :: z_cell ! Cell center heights referenced to surface [m] (negative in ocean)
1025 real, dimension( GV%ke ) :: obl_depth ! Cell center depths referenced to surface [m] (positive in ocean)
1026 real, dimension( GV%ke+1 ) :: z_inter ! Cell interface heights referenced to surface [m] (negative in ocean)
1027 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 real, dimension( 3*GV%ke ) :: rho_1d ! A column of densities [R ~> kg m-3]
1035 real, dimension( 3*GV%ke ) :: pres_1d ! A column of pressures [R L2 T-2 ~> Pa]
1036 real, dimension( 3*GV%ke ) :: temp_1d ! A column of temperatures [C ~> degC]
1037 real, dimension( 3*GV%ke ) :: salt_1d ! A column of salinities [S ~> ppt]
1038
1039 real, dimension( GV%ke ) :: cellheight ! Cell center heights referenced to surface [Z ~> m] (negative in ocean)
1040 real, dimension( GV%ke+1 ) :: ifaceheight ! Interface heights referenced to surface [Z ~> m] (negative in ocean)
1041 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 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 real, dimension(GV%ke) :: ue_h, ve_h ! Eulerian velocities h-points, centers [L T-1 ~> m s-1]
1074 real, dimension(GV%ke) :: us_h, vs_h ! Stokes drift components h-points, centers [L T-1 ~> m s-1]
1075 real, dimension(GV%ke) :: usbar_h, vsbar_h ! Cell Average Stokes drift h-points [L T-1 ~> m s-1]
1076 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 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 if (cs%Stokes_Mixing .and. .not.associated(waves)) call mom_error(fatal, &
1084 "KPP_compute_BLD: The Waves control structure must be associated if STOKES_MIXING is True.")
1085
1086 if (cs%debug) then
1087 call hchksum(salt, "KPP in: S", g%HI, haloshift=0, unscale=us%S_to_ppt)
1088 call hchksum(temp, "KPP in: T", g%HI, haloshift=0, unscale=us%C_to_degC)
1089 call hchksum(u, "KPP in: u", g%HI, haloshift=0, unscale=us%L_T_to_m_s)
1090 call hchksum(v, "KPP in: v", g%HI, haloshift=0, unscale=us%L_T_to_m_s)
1091 endif
1092
1093 call cpu_clock_begin(id_clock_kpp_compute_bld)
1094
1095 ! some constants
1096 gorho = us%Z_to_m*us%s_to_T**2 * (gv%g_Earth_Z_T2 / gv%Rho0)
1097 if (gv%Boussinesq) then
1098 gorho_z_l2 = gv%Z_to_H * gv%g_Earth_Z_T2 / gv%Rho0
1099 else
1100 gorho_z_l2 = gv%g_Earth_Z_T2 * gv%RZ_to_H
1101 endif
1102 buoy_scale = us%L_to_m**2*us%s_to_T**3
1103
1104 ! Find the vertical distances across layers.
1105 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 do j = g%jsc, g%jec
1122 do i = g%isc, g%iec ; if (g%mask2dT(i,j) > 0.0) then
1123
1124 do k=1,gv%ke
1125 u_h(k) = 0.5 * (u(i,j,k)+u(i-1,j,k))
1126 v_h(k) = 0.5 * (v(i,j,k)+v(i,j-1,k))
1127 enddo
1128 if (cs%StokesMOST) then
1129 do k=1,gv%ke
1130 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 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.25*us%s_to_T*( (g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j-1)) + &
1136 (g%CoriolisBu(i-1,j) + g%CoriolisBu(i,j-1)) )
1137 surffricvel = us%Z_to_m*us%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 ifaceheight(1) = 0.0 ! BBL is all relative to the surface
1146 pref = 0. ; if (associated(tv%p_surf)) pref = tv%p_surf(i,j)
1147 hcorr = 0.
1148
1149 if (cs%StokesMOST) call compute_stokesdrift( i, j, h(i,j,1) , ifaceheight(1), &
1150 us_hi(1), vs_hi(1), us_h(1), vs_h(1), usbar_h(1), vsbar_h(1), waves)
1151
1152 do k=1,gv%ke
1153 ! cell center and cell bottom in meters (negative values in the ocean)
1154 dh = dz(i,j,k) ! Nominal thickness to use for increment
1155 dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0)
1156 hcorr = min( dh - cs%min_thickness, 0. ) ! If inflating then hcorr<0
1157 dh = max( dh, cs%min_thickness ) ! Limit increment dh>=min_thickness
1158 cellheight(k) = ifaceheight(k) - 0.5 * dh
1159 ifaceheight(k+1) = ifaceheight(k) - dh
1160
1161 ! find ksfc for cell where "surface layer" sits
1162 sldepth_0d = cs%surf_layer_ext*max( max(-cellheight(k),-ifaceheight(2) ), cs%minOBLdepth )
1163 ksfc = k
1164 do ktmp = 1,k
1165 if (-1.0*ifaceheight(ktmp+1) >= sldepth_0d) then
1166 ksfc = ktmp
1167 exit
1168 endif
1169 enddo
1170
1171 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 (buoyflux(i,j,1) - 0.5*(buoyflux(i,j,max(2,k))+buoyflux(i,j,k+1)) )
1176 surfbuoyflux2(k) = surfbuoyflux
1177 call compute_stokesdrift(i,j, ifaceheight(k),ifaceheight(k+1), &
1178 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 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 vs_sld, usbar_sld, vsbar_sld, stokesxi, cvmix_kpp_params_user=cs%KPP_params )
1184 stokesxi_1d(k) = stokesxi
1185 stokesvt_1d(k) = 0.0 ! StokesXI
1186
1187 ! average temperature, salinity, u and v over surface layer starting at ksfc
1188 delh = sldepth_0d + ifaceheight(ksfc)
1189 surfhtemp = temp(i,j,ksfc) * delh
1190 surfhsalt = salt(i,j,ksfc) * delh
1191 surfhu = (ue_h(ksfc) + usbar_sld) * delh
1192 surfhv = (ve_h(ksfc) + vsbar_sld) * delh
1193 htot = delh
1194 do ktmp = 1,ksfc-1 ! if ksfc >=2
1195 delh = h(i,j,ktmp)*gv%H_to_Z
1196 htot = htot + delh
1197 surfhtemp = surfhtemp + temp(i,j,ktmp) * delh
1198 surfhsalt = surfhsalt + salt(i,j,ktmp) * delh
1199 surfhu = surfhu + (ue_h(ktmp) + usbar_h(ktmp)) * delh
1200 surfhv = surfhv + (ve_h(ktmp) + vsbar_h(ktmp)) * delh
1201 enddo
1202 i_htot = 1./htot
1203 surftemp = surfhtemp * i_htot
1204 surfsalt = surfhsalt * i_htot
1205 surfu = surfhu * i_htot
1206 surfv = surfhv * i_htot
1207
1208 uk = ue_h(k) + us_h(k) - surfu
1209 vk = ve_h(k) + vs_h(k) - surfv
1210
1211 else !not StokesMOST
1212 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 surfhtemp = 0.0
1216 surfhsalt = 0.0
1217 surfhu = 0.0
1218 surfhv = 0.0
1219 surfhus = 0.0
1220 surfhvs = 0.0
1221 htot = 0.0
1222 do ktmp = 1,ksfc
1223
1224 ! SLdepth_0d can be between cell interfaces
1225 delh = min( max(0.0, sldepth_0d - htot), dz(i,j,ktmp) )
1226
1227 ! surface layer thickness
1228 htot = htot + delh
1229
1230 ! surface averaged fields
1231 surfhtemp = surfhtemp + temp(i,j,ktmp) * delh
1232 surfhsalt = surfhsalt + salt(i,j,ktmp) * delh
1233 surfhu = surfhu + 0.5*(u(i,j,ktmp)+u(i-1,j,ktmp)) * delh
1234 surfhv = surfhv + 0.5*(v(i,j,ktmp)+v(i,j-1,ktmp)) * delh
1235 if (cs%Stokes_Mixing) then
1236 surfhus = surfhus + 0.5*(waves%US_x(i,j,ktmp)+waves%US_x(i-1,j,ktmp)) * delh
1237 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 surftemp = surfhtemp / htot
1250 surfsalt = surfhsalt / htot
1251 surfu = surfhu / htot
1252 surfv = surfhv / htot
1253 surfus = surfhus / htot
1254 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 uk = 0.5*(u(i,j,k)+u(i-1,j,k)) - surfu
1258 vk = 0.5*(v(i,j,k)+v(i,j-1,k)) - surfv
1259
1260 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 uk = uk + (0.5*(waves%Us_x(i,j,k)+waves%US_x(i-1,j,k)) - surfus )
1265 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 surfbuoyflux = buoy_scale * (buoyflux(i,j,1) - buoyflux(i,j,k+1))
1270 surfbuoyflux2(k) = surfbuoyflux
1271
1272 endif ! StokesMOST
1273
1274 deltau2(k) = us%L_T_to_m_s**2 * ((uk**2) + (vk**2))
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 km1 = max(1, k-1)
1281 kk = 3*(k-1)
1282 pres_1d(kk+1) = pref
1283 pres_1d(kk+2) = pref
1284 pres_1d(kk+3) = pref
1285 temp_1d(kk+1) = surftemp
1286 temp_1d(kk+2) = temp(i,j,k)
1287 temp_1d(kk+3) = temp(i,j,km1)
1288 salt_1d(kk+1) = surfsalt
1289 salt_1d(kk+2) = salt(i,j,k)
1290 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 pref = pref + (gv%g_Earth * gv%H_to_RZ) * h(i,j,k)
1295
1296 enddo ! k-loop finishes
1297
1298 if ( (cs%LT_K_ENHANCEMENT .or. cs%LT_VT2_ENHANCEMENT) .and. .not. present(lamult)) then
1299 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 dz=dz(i,j,:), u_h=u_h, v_h=v_h, waves=waves)
1302 cs%La_SL(i,j) = la
1303 endif
1304
1305
1306 ! compute in-situ density
1307 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 do k = 1, gv%ke
1313 km1 = max(1, k-1)
1314 kk = 3*(k-1)
1315 deltarho(k) = rho_1d(kk+2) - rho_1d(kk+1)
1316 if (gv%Boussinesq .or. gv%semi_Boussinesq) then
1317 deltabuoy(k) = gorho*(rho_1d(kk+2) - rho_1d(kk+1))
1318 else
1319 deltabuoy(k) = (us%Z_to_m*us%s_to_T**2) * gv%g_Earth_Z_T2 * &
1320 ( (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.5*(h(i,j,km1) + h(i,j,k))+gv%H_subroundoff))
1324 cs%N(i,j,k) = sqrt( max( n2_1d(k), 0.) )
1325 enddo
1326 n2_1d(gv%ke+1 ) = 0.0
1327 cs%N(i,j,gv%ke+1 ) = 0.0
1328
1329 ! Convert columns to MKS units for passing to CVMix
1330 do k = 1, gv%ke
1331 obl_depth(k) = -us%Z_to_m * cellheight(k)
1332 z_cell(k) = us%Z_to_m*cellheight(k)
1333 enddo
1334 do k = 1, gv%ke+1
1335 n_col(k) = us%s_to_T*cs%N(i,j,k)
1336 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 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 cvmix_kpp_params_user=cs%KPP_params )
1352
1353 ! Determine the enhancement factor for unresolved shear
1354 IF (cs%LT_VT2_ENHANCEMENT) then
1355 IF (cs%LT_VT2_METHOD==lt_vt2_mode_constant) then
1356 langenhvt2 = cs%KPP_VT2_ENH_FAC
1357 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 if (present(lamult)) then
1360 langenhvt2 = lamult(i,j)
1361 else
1362 langenhvt2 = sqrt(1.+(1.5*cs%La_SL(i,j))**(-2) + &
1363 (5.4*cs%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 langenhvt2 = 1.0
1369 endif
1370 else
1371 langenhvt2 = 1.0
1372 endif
1373
1374 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 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 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 cvmix_kpp_params_user=cs%KPP_params ) ! KPP parameters
1409 cs%OBLdepth(i,j) = us%m_to_Z * kpp_obl_depth
1410
1411 if (cs%StokesMOST) then
1412 kbl = int(cs%kOBL(i,j))
1413 sldepth_0d = cs%surf_layer_ext*cs%OBLdepth(i,j)
1414 surfbuoyflux = surfbuoyflux2(kbl)
1415 ! find ksfc for cell where "surface layer" sits
1416 ksfc = kbl
1417 do ktmp = 1, kbl
1418 if (-1.0*ifaceheight(ktmp+1) >= sldepth_0d) then
1419 ksfc = ktmp
1420 exit
1421 endif
1422 enddo
1423
1424 call compute_stokesdrift(i,j, ifaceheight(ksfc) , -sldepth_0d, &
1425 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 stokesxi, cvmix_kpp_params_user=cs%KPP_params )
1430 cs%StokesParXI(i,j) = stokesxi
1431 cs%Lam2(i,j) = sqrt(us_hi(1)**2+vs_hi(1)**2) / max(surffricvel,0.0002)
1432
1433 else !.not Stokes_MOST
1434 cs%StokesParXI(i,j) = 10.0
1435 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 if (cs%deepOBLoffset>0.) then
1440 zbottomminusoffset = ifaceheight(gv%ke+1) + min(cs%deepOBLoffset, -0.1*ifaceheight(gv%ke+1))
1441 cs%OBLdepth(i,j) = min( cs%OBLdepth(i,j), -zbottomminusoffset )
1442 endif
1443
1444 ! apply some constraints on OBLdepth
1445 if (cs%fixedOBLdepth) cs%OBLdepth(i,j) = cs%fixedOBLdepth_value
1446 cs%OBLdepth(i,j) = max( cs%OBLdepth(i,j), -ifaceheight(2) ) ! no shallower than top layer
1447 cs%OBLdepth(i,j) = min( cs%OBLdepth(i,j), -ifaceheight(gv%ke+1) ) ! no deeper than bottom
1448 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 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 cvmix_kpp_params_user=cs%KPP_params ) ! KPP parameters
1463 cs%Vt2(i,j,:) = us%m_to_Z**2*us%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 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 cvmix_kpp_params_user=cs%KPP_params) ! KPP parameters
1477 cs%Ws(i,j,:) = us%m_to_Z*us%T_to_s*ws_1d(:)
1478 endif
1479
1480 ! Diagnostics
1481 if (cs%id_N2 > 0) cs%N2(i,j,:) = n2_1d(:)
1482 if (cs%id_BulkDrho > 0) cs%dRho(i,j,:) = deltarho(:)
1483 if (cs%id_BulkRi > 0) cs%BulkRi(i,j,:) = bulkri_1d(:)
1484 if (cs%id_BulkUz2 > 0) cs%Uz2(i,j,:) = us%m_s_to_L_T**2 * deltau2(:)
1485 if (cs%id_Tsurf > 0) cs%Tsurf(i,j) = surftemp
1486 if (cs%id_Ssurf > 0) cs%Ssurf(i,j) = surfsalt
1487 if (cs%id_Usurf > 0) cs%Usurf(i,j) = surfu
1488 if (cs%id_Vsurf > 0) cs%Vsurf(i,j) = surfv
1489
1490 endif ; enddo
1491 enddo
1492
1493 call cpu_clock_end(id_clock_kpp_compute_bld)
1494
1495 ! send diagnostics to post_data
1496 if (cs%id_BulkRi > 0) call post_data(cs%id_BulkRi, cs%BulkRi, cs%diag)
1497 if (cs%id_N > 0) call post_data(cs%id_N, cs%N, cs%diag)
1498 if (cs%id_N2 > 0) call post_data(cs%id_N2, cs%N2, cs%diag)
1499 if (cs%id_Tsurf > 0) call post_data(cs%id_Tsurf, cs%Tsurf, cs%diag)
1500 if (cs%id_Ssurf > 0) call post_data(cs%id_Ssurf, cs%Ssurf, cs%diag)
1501 if (cs%id_Usurf > 0) call post_data(cs%id_Usurf, cs%Usurf, cs%diag)
1502 if (cs%id_Vsurf > 0) call post_data(cs%id_Vsurf, cs%Vsurf, cs%diag)
1503 if (cs%id_BulkDrho > 0) call post_data(cs%id_BulkDrho, cs%dRho, cs%diag)
1504 if (cs%id_BulkUz2 > 0) call post_data(cs%id_BulkUz2, cs%Uz2, cs%diag)
1505 if (cs%id_EnhK > 0) call post_data(cs%id_EnhK, cs%EnhK, cs%diag)
1506 if (cs%id_EnhVt2 > 0) call post_data(cs%id_EnhVt2, cs%EnhVt2, cs%diag)
1507 if (cs%id_La_SL > 0) call post_data(cs%id_La_SL, cs%La_SL, cs%diag)
1508 if (cs%id_Vt2 > 0) call post_data(cs%id_Vt2, cs%Vt2, cs%diag)
1509
1510 if (cs%StokesMOST) then
1511 if (cs%id_StokesXI > 0) call post_data(cs%id_StokesXI, cs%StokesParXI, cs%diag)
1512 if (cs%id_Lam2 > 0) call post_data(cs%id_Lam2 , cs%Lam2 , cs%diag)
1513 endif
1514
1515 ! BLD smoothing:
1516 if (cs%n_smooth > 0) call kpp_smooth_bld(cs, g, gv, us, dz)
1517
1518end 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
1522subroutine 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 real, dimension(SZI_(G),SZJ_(G)) :: OBLdepth_prev ! OBLdepth before s.th smoothing iteration [Z ~> m]
1532 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 real, dimension( GV%ke ) :: cellHeight ! Cell center heights referenced to surface [Z ~> m]
1535 ! (negative in the ocean)
1536 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 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 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 total_depth(:,:) = 0.0
1548
1549 !$OMP parallel do default(shared) private(dh, h_cor)
1550 do j = g%jsc, g%jec
1551 h_cor(:) = 0.
1552 do k=1,gv%ke
1553 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 dh = dz(i,j,k) ! Nominal thickness to use for increment
1556 dh = dh + h_cor(i) ! Take away the accumulated error (could temporarily make dh<0)
1557 h_cor(i) = min( dh - cs%min_thickness, 0. ) ! If inflating then hcorr<0
1558 dh = max( dh, cs%min_thickness ) ! Limit increment dh>=min_thickness
1559 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 if (cs%n_smooth > 1) call pass_var(total_depth, g%Domain, halo=cs%n_smooth, complete=.false.)
1570 call pass_var(cs%OBLdepth, g%Domain, halo=cs%n_smooth)
1571
1572 if (cs%id_OBLdepth_original > 0) cs%OBLdepth_original(:,:) = cs%OBLdepth(:,:)
1573
1574 do s=1,cs%n_smooth
1575
1576 obldepth_prev(:,:) = cs%OBLdepth(:,:)
1577 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 do j = g%jsc-halo, g%jec+halo
1583 do i = g%isc-halo, g%iec+halo ; if (g%mask2dT(i,j) > 0.0) then
1584 ! compute weights
1585 ww = 0.125 * g%mask2dT(i-1,j)
1586 we = 0.125 * g%mask2dT(i+1,j)
1587 ws = 0.125 * g%mask2dT(i,j-1)
1588 wn = 0.125 * g%mask2dT(i,j+1)
1589 wc = 1.0 - (ww+we+wn+ws)
1590
1591 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 + 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 + (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 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 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 do j=g%jsc,g%jec ; do i=g%isc,g%iec ; if (g%mask2dT(i,j) > 0.0) then
1617
1618 ifaceheight(1) = 0.0 ! BBL is all relative to the surface
1619 hcorr = 0.
1620 do k=1,gv%ke
1621 ! cell center and cell bottom in meters (negative values in the ocean)
1622 dh = dz(i,j,k) ! Nominal thickness to use for increment
1623 dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0)
1624 hcorr = min( dh - cs%min_thickness, 0. ) ! If inflating then hcorr<0
1625 dh = max( dh, cs%min_thickness ) ! Limit increment dh>=min_thickness
1626 cellheight(k) = ifaceheight(k) - 0.5 * dh
1627 ifaceheight(k+1) = ifaceheight(k) - dh
1628 enddo
1629
1630 cs%kOBL(i,j) = cvmix_kpp_compute_kobl_depth( ifaceheight, cellheight, cs%OBLdepth(i,j) )
1631 endif ; enddo ; enddo
1632
1633 call cpu_clock_end(id_clock_kpp_smoothing)
1634
1635end 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.
1639subroutine 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 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 do j = g%jsc, g%jec ; do i = g%isc, g%iec
1655 bld(i,j) = scale * cs%OBLdepth(i,j)
1656 enddo ; enddo
1657
1658end subroutine kpp_get_bld
1659
1660!> Apply KPP non-local transport of surface fluxes for a given tracer
1661subroutine kpp_nonlocaltransport(CS, G, GV, h, nonLocalTrans, surfFlux, &
1662 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 real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dtracer ! Rate of tracer change [conc T-1 ~> conc s-1]
1679 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 if (present(flux_scale)) then
1684 do j = g%jsc, g%jec ; do i = g%isc, g%iec
1685 surfflux_loc(i,j) = surfflux(i,j) * flux_scale
1686 enddo ; enddo
1687 else
1688 surfflux_loc(:,:) = surfflux(:,:)
1689 endif
1690
1691 ! Post surface flux diagnostic
1692 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 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 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 ( 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 if (cs%applyNonLocalTrans) then
1706 !$OMP parallel do default(none) shared(G, GV, dt, scalar, dtracer)
1707 do k = 1, gv%ke ; do j = g%jsc, g%jec ; do i = g%isc, g%iec
1708 scalar(i,j,k) = scalar(i,j,k) + dt * dtracer(i,j,k)
1709 enddo ; enddo ; enddo
1710 endif
1711 if (tr_ptr%id_NLT_tendency > 0) call post_data(tr_ptr%id_NLT_tendency, dtracer, diag)
1712
1713 endif
1714
1715
1716 if (tr_ptr%id_NLT_budget > 0) then
1717 !$OMP parallel do default(none) shared(G, GV, dtracer, nonLocalTrans, surfFlux_loc)
1718 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 dtracer(i,j,k) = (nonlocaltrans(i,j,k) - nonlocaltrans(i,j,k+1)) * surfflux_loc(i,j)
1721 enddo ; enddo ; enddo
1722 call post_data(tr_ptr%id_NLT_budget, dtracer(:,:,:), diag)
1723 endif
1724
1725end subroutine kpp_nonlocaltransport
1726
1727
1728!> Apply KPP non-local transport of surface fluxes for temperature.
1729subroutine 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 tr_ptr, scalar)
1745
1746end subroutine kpp_nonlocaltransport_temp
1747
1748
1749!> Apply KPP non-local transport of surface fluxes for salinity.
1750subroutine 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 tr_ptr, scalar)
1764
1765end 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
1769subroutine 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 us_i = 0.0
1789 vs_i = 0.0
1790 us_k = 0.0
1791 vs_k = 0.0
1792 usbar = 0.0
1793 vsbar = 0.0
1794 do b = 1, waves%NumBands
1795 wavenum = waves%WaveNum_Cen(b)
1796 fexp = exp(2. * wavenum * zbot)
1797 us_i = us_i + waves%Ustk_Hb(i,j,b) * fexp
1798 vs_i = vs_i + waves%Vstk_Hb(i,j,b) * fexp
1799 fexp = exp( wavenum * (ztop + zbot) )
1800 us_k = us_k+ waves%Ustk_Hb(i,j,b) * fexp
1801 vs_k = vs_k+ waves%Vstk_Hb(i,j,b) * fexp
1802 fexp = exp(2. * wavenum * ztop) - exp(2. * wavenum * zbot)
1803 usbar = usbar + 0.5 * waves%Ustk_Hb(i,j,b) * fexp / wavenum
1804 vsbar = vsbar + 0.5 * waves%Vstk_Hb(i,j,b) * fexp / wavenum
1805 enddo
1806 usbar = usbar / (ztop-zbot)
1807 vsbar = vsbar / (ztop-zbot)
1808
1809end subroutine compute_stokesdrift
1810
1811!> Clear pointers, deallocate memory
1812subroutine kpp_end(CS)
1813 type(kpp_cs), pointer :: cs !< Control structure
1814
1815 if (.not.associated(cs)) return
1816
1817 deallocate(cs)
1818
1819end subroutine kpp_end
1820
1821end module mom_cvmix_kpp