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