MOM6
MOM_CVMix_KPP.F90
Go to the documentation of this file.
1 !> Provides the K-Profile Parameterization (KPP) of Large et al., 1994, via CVMix.
2 module mom_cvmix_kpp
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
6 use mom_coms, only : max_across_pes
7 use mom_debugging, only : hchksum, is_nan
8 use mom_diag_mediator, only : time_type, diag_ctrl, safe_alloc_ptr, post_data
9 use mom_diag_mediator, only : query_averaging_enabled, register_diag_field
10 use mom_error_handler, only : mom_error, mom_mesg, fatal, warning, is_root_pe
11 use mom_eos, only : eos_type, calculate_density
12 use mom_file_parser, only : get_param, log_param, log_version, param_file_type
13 use mom_file_parser, only : openparameterblock, closeparameterblock
14 use mom_grid, only : ocean_grid_type, ispointincell
15 use mom_unit_scaling, only : unit_scale_type
16 use mom_verticalgrid, only : verticalgrid_type
17 use mom_wave_interface, only : wave_parameters_cs, get_langmuir_number
18 use mom_domains, only : pass_var
19 use mpp_mod, only : mpp_pe
20 
21 use cvmix_kpp, only : cvmix_init_kpp, cvmix_put_kpp, cvmix_get_kpp_real
22 use cvmix_kpp, only : cvmix_coeffs_kpp
23 use cvmix_kpp, only : cvmix_kpp_compute_obl_depth
24 use cvmix_kpp, only : cvmix_kpp_compute_turbulent_scales
25 use cvmix_kpp, only : cvmix_kpp_compute_bulk_richardson
26 use cvmix_kpp, only : cvmix_kpp_compute_unresolved_shear
27 use cvmix_kpp, only : cvmix_kpp_params_type
28 use cvmix_kpp, only : cvmix_kpp_compute_kobl_depth
29 
30 implicit none ; private
31 
32 #include "MOM_memory.h"
33 
34 public :: kpp_init
35 public :: kpp_compute_bld
36 public :: kpp_calculate
37 public :: kpp_end
38 public :: kpp_nonlocaltransport_temp
39 public :: kpp_nonlocaltransport_saln
40 public :: kpp_get_bld
41 
42 ! Enumerated constants
43 integer, private, parameter :: nlt_shape_cvmix = 0 !< Use the CVMix profile
44 integer, private, parameter :: nlt_shape_linear = 1 !< Linear, \f$ G(\sigma) = 1-\sigma \f$
45 integer, private, parameter :: nlt_shape_parabolic = 2 !< Parabolic, \f$ G(\sigma) = (1-\sigma)^2 \f$
46 integer, private, parameter :: nlt_shape_cubic = 3 !< Cubic, \f$ G(\sigma) = 1 + (2\sigma-3) \sigma^2\f$
47 integer, private, parameter :: nlt_shape_cubic_lmd = 4 !< Original shape,
48  !! \f$ G(\sigma) = \frac{27}{4} \sigma (1-\sigma)^2 \f$
49 
50 integer, private, parameter :: sw_method_all_sw = 0 !< Use all shortwave radiation
51 integer, private, parameter :: sw_method_mxl_sw = 1 !< Use shortwave radiation absorbed in mixing layer
52 integer, private, parameter :: sw_method_lv1_sw = 2 !< Use shortwave radiation absorbed in layer 1
53 integer, private, parameter :: lt_k_constant = 1, & !< Constant enhance K through column
54  lt_k_scaled = 2, & !< Enhance K scales with G(sigma)
55  lt_k_mode_constant = 1, & !< Prescribed enhancement for K
56  lt_k_mode_vr12 = 2, & !< Enhancement for K based on
57  !! Van Roekel et al., 2012
58  lt_k_mode_rw16 = 3, & !< Enhancement for K based on
59  !! Reichl et al., 2016
60  lt_vt2_mode_constant = 1, & !< Prescribed enhancement for Vt2
61  lt_vt2_mode_vr12 = 2, & !< Enhancement for Vt2 based on
62  !! Van Roekel et al., 2012
63  lt_vt2_mode_rw16 = 3, & !< Enhancement for Vt2 based on
64  !! Reichl et al., 2016
65  lt_vt2_mode_lf17 = 4 !< Enhancement for Vt2 based on
66  !! Li and Fox-Kemper, 2017
67 
68 integer :: ncall_smooth = 0
69 
70 !> Control structure for containing KPP parameters/data
71 type, public :: kpp_cs ; private
72 
73  ! Parameters
74  real :: ri_crit !< Critical bulk Richardson number (defines OBL depth)
75  real :: vonkarman !< von Karman constant (dimensionless)
76  real :: cs !< Parameter for computing velocity scale function (dimensionless)
77  real :: cs2 !< Parameter for multiplying by non-local term
78  ! This is active for NLT_SHAPE_CUBIC_LMD only
79  logical :: enhance_diffusion !< If True, add enhanced diffusivity at base of boundary layer.
80  character(len=10) :: interptype !< Type of interpolation to compute bulk Richardson number
81  character(len=10) :: interptype2 !< Type of interpolation to compute diff and visc at OBL_depth
82  logical :: computeekman !< If True, compute Ekman depth limit for OBLdepth
83  logical :: computemoninobukhov !< If True, compute Monin-Obukhov limit for OBLdepth
84  logical :: passivemode !< If True, makes KPP passive meaning it does NOT alter the diffusivity
85  real :: deepobloffset !< If non-zero, is a distance from the bottom that the OBL can not
86  !! penetrate through [m]
87  real :: minobldepth !< If non-zero, is a minimum depth for the OBL [m]
88  real :: surf_layer_ext !< Fraction of OBL depth considered in the surface layer [nondim]
89  real :: minvtsqr !< Min for the squared unresolved velocity used in Rib CVMix calculation [m2 s-2]
90  logical :: fixedobldepth !< If True, will fix the OBL depth at fixedOBLdepth_value
91  real :: fixedobldepth_value !< value for the fixed OBL depth when fixedOBLdepth==True.
92  logical :: debug !< If True, calculate checksums and write debugging information
93  character(len=30) :: matchtechnique !< Method used in CVMix for setting diffusivity and NLT profile functions
94  integer :: nlt_shape !< MOM6 over-ride of CVMix NLT shape function
95  logical :: applynonlocaltrans !< If True, apply non-local transport to heat and scalars
96  integer :: n_smooth !< Number of times smoothing operator is applied on OBLdepth.
97  logical :: deepen_only !< If true, apply OBLdepth smoothing at a cell only if the OBLdepth gets deeper.
98  logical :: kppzerodiffusivity !< If True, will set diffusivity and viscosity from KPP to zero
99  !! for testing purposes.
100  logical :: kppisadditive !< If True, will add KPP diffusivity to initial diffusivity.
101  !! If False, will replace initial diffusivity wherever KPP diffusivity
102  !! is non-zero.
103  real :: min_thickness !< A minimum thickness used to avoid division by small numbers
104  !! in the vicinity of vanished layers.
105  ! smg: obsolete below
106  logical :: correctsurflayeravg !< If true, applies a correction to the averaging of surface layer properties
107  real :: surflayerdepth !< A guess at the depth of the surface layer (which should 0.1 of OBLdepth) [m]
108  ! smg: obsolete above
109  integer :: sw_method !< Sets method for using shortwave radiation in surface buoyancy flux
110  logical :: lt_k_enhancement !< Flags if enhancing mixing coefficients due to LT
111  integer :: lt_k_shape !< Integer for constant or shape function enhancement
112  integer :: lt_k_method !< Integer for mixing coefficients LT method
113  real :: kpp_k_enh_fac !< Factor to multiply by K if Method is CONSTANT
114  logical :: lt_vt2_enhancement !< Flags if enhancing Vt2 due to LT
115  integer :: lt_vt2_method !< Integer for Vt2 LT method
116  real :: kpp_vt2_enh_fac !< Factor to multiply by VT2 if Method is CONSTANT
117  logical :: stokes_mixing !< Flag if model is mixing down Stokes gradient
118  !! This is relavent for which current to use in RiB
119 
120  !> CVMix parameters
121  type(cvmix_kpp_params_type), pointer :: kpp_params => null()
122 
123  type(diag_ctrl), pointer :: diag => null() !< Pointer to diagnostics control structure
124  !>@{ Diagnostic handles
125  integer :: id_obldepth = -1, id_bulkri = -1
126  integer :: id_n = -1, id_n2 = -1
127  integer :: id_ws = -1, id_vt2 = -1
128  integer :: id_bulkuz2 = -1, id_bulkdrho = -1
129  integer :: id_ustar = -1, id_buoyflux = -1
130  integer :: id_qminussw = -1, id_nets = -1
131  integer :: id_sigma = -1, id_kv_kpp = -1
132  integer :: id_kt_kpp = -1, id_ks_kpp = -1
133  integer :: id_tsurf = -1, id_ssurf = -1
134  integer :: id_usurf = -1, id_vsurf = -1
135  integer :: id_kd_in = -1
136  integer :: id_nltt = -1
137  integer :: id_nlts = -1
138  integer :: id_nlt_dsdt = -1
139  integer :: id_nlt_dtdt = -1
140  integer :: id_nlt_temp_budget = -1
141  integer :: id_nlt_saln_budget = -1
142  integer :: id_enhk = -1, id_enhvt2 = -1
143  integer :: id_enhw = -1
144  integer :: id_la_sl = -1
145  integer :: id_obldepth_original = -1
146  !!@}
147 
148  ! Diagnostics arrays
149  real, allocatable, dimension(:,:) :: obldepth !< Depth (positive) of OBL [m]
150  real, allocatable, dimension(:,:) :: obldepth_original !< Depth (positive) of OBL [m] without smoothing
151  real, allocatable, dimension(:,:) :: kobl !< Level (+fraction) of OBL extent
152  real, allocatable, dimension(:,:) :: obldepthprev !< previous Depth (positive) of OBL [m]
153  real, allocatable, dimension(:,:) :: la_sl !< Langmuir number used in KPP
154  real, allocatable, dimension(:,:,:) :: drho !< Bulk difference in density [kg m-3]
155  real, allocatable, dimension(:,:,:) :: uz2 !< Square of bulk difference in resolved velocity [m2 s-2]
156  real, allocatable, dimension(:,:,:) :: bulkri !< Bulk Richardson number for each layer (dimensionless)
157  real, allocatable, dimension(:,:,:) :: sigma !< Sigma coordinate (dimensionless)
158  real, allocatable, dimension(:,:,:) :: ws !< Turbulent velocity scale for scalars [m s-1]
159  real, allocatable, dimension(:,:,:) :: n !< Brunt-Vaisala frequency [s-1]
160  real, allocatable, dimension(:,:,:) :: n2 !< Squared Brunt-Vaisala frequency [s-2]
161  real, allocatable, dimension(:,:,:) :: vt2 !< Unresolved squared turbulence velocity for bulk Ri [m2 s-2]
162  real, allocatable, dimension(:,:,:) :: kt_kpp !< Temp diffusivity from KPP [m2 s-1]
163  real, allocatable, dimension(:,:,:) :: ks_kpp !< Scalar diffusivity from KPP [m2 s-1]
164  real, allocatable, dimension(:,:,:) :: kv_kpp !< Viscosity due to KPP [m2 s-1]
165  real, allocatable, dimension(:,:) :: tsurf !< Temperature of surface layer [degC]
166  real, allocatable, dimension(:,:) :: ssurf !< Salinity of surface layer [ppt]
167  real, allocatable, dimension(:,:) :: usurf !< i-velocity of surface layer [m s-1]
168  real, allocatable, dimension(:,:) :: vsurf !< j-velocity of surface layer [m s-1]
169  real, allocatable, dimension(:,:,:) :: enhk !< Enhancement for mixing coefficient
170  real, allocatable, dimension(:,:,:) :: enhvt2 !< Enhancement for Vt2
171 
172 end type kpp_cs
173 
174 #define __DO_SAFETY_CHECKS__
175 
176 contains
177 
178 !> Initialize the CVMix KPP module and set up diagnostics
179 !! Returns True if KPP is to be used, False otherwise.
180 logical function kpp_init(paramFile, G, GV, US, diag, Time, CS, passive, Waves)
181 
182  ! Arguments
183  type(param_file_type), intent(in) :: paramfile !< File parser
184  type(ocean_grid_type), intent(in) :: g !< Ocean grid
185  type(verticalgrid_type), intent(in) :: gv !< Vertical grid structure.
186  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
187  type(diag_ctrl), target, intent(in) :: diag !< Diagnostics
188  type(time_type), intent(in) :: time !< Model time
189  type(kpp_cs), pointer :: cs !< Control structure
190  logical, optional, intent(out) :: passive !< Copy of %passiveMode
191  type(wave_parameters_cs), optional, pointer :: waves !< Wave CS
192 
193  ! Local variables
194 #include "version_variable.h"
195  character(len=40) :: mdl = 'MOM_CVMix_KPP' !< name of this module
196  character(len=20) :: string !< local temporary string
197  logical :: cs_is_one=.false. !< Logical for setting Cs based on Non-local
198  logical :: lnodgat1=.false. !< True => G'(1) = 0 (shape function)
199  !! False => compute G'(1) as in LMD94
200  if (associated(cs)) call mom_error(fatal, 'MOM_CVMix_KPP, KPP_init: '// &
201  'Control structure has already been initialized')
202 
203  ! Read parameters
204  call log_version(paramfile, mdl, version, 'This is the MOM wrapper to CVMix:KPP\n' // &
205  'See http://cvmix.github.io/')
206  call get_param(paramfile, mdl, "USE_KPP", kpp_init, &
207  "If true, turns on the [CVMix] KPP scheme of Large et al., 1994, "// &
208  "to calculate diffusivities and non-local transport in the OBL.", &
209  default=.false.)
210  ! Forego remainder of initialization if not using this scheme
211  if (.not. kpp_init) return
212  allocate(cs)
213 
214  call openparameterblock(paramfile,'KPP')
215  call get_param(paramfile, mdl, 'PASSIVE', cs%passiveMode, &
216  'If True, puts KPP into a passive-diagnostic mode.', &
217  default=.false.)
218  !BGR: Note using PASSIVE for KPP creates warning for PASSIVE from Convection
219  ! should we create a separate flag?
220  if (present(passive)) passive=cs%passiveMode ! This is passed back to the caller so
221  ! the caller knows to not use KPP output
222  call get_param(paramfile, mdl, 'APPLY_NONLOCAL_TRANSPORT', cs%applyNonLocalTrans, &
223  'If True, applies the non-local transport to heat and scalars. '// &
224  'If False, calculates the non-local transport and tendencies but '//&
225  'purely for diagnostic purposes.', &
226  default=.not. cs%passiveMode)
227  call get_param(paramfile, mdl, 'N_SMOOTH', cs%n_smooth, &
228  'The number of times the 1-1-4-1-1 Laplacian filter is applied on '// &
229  'OBL depth.', &
230  default=0)
231  if (cs%n_smooth > 0) then
232  call get_param(paramfile, mdl, 'DEEPEN_ONLY_VIA_SMOOTHING', cs%deepen_only, &
233  'If true, apply OBLdepth smoothing at a cell only if the OBLdepth '// &
234  'gets deeper via smoothing.', &
235  default=.false.)
236  endif
237  call get_param(paramfile, mdl, 'RI_CRIT', cs%Ri_crit, &
238  'Critical bulk Richardson number used to define depth of the '// &
239  'surface Ocean Boundary Layer (OBL).', &
240  units='nondim', default=0.3)
241  call get_param(paramfile, mdl, 'VON_KARMAN', cs%vonKarman, &
242  'von Karman constant.', &
243  units='nondim', default=0.40)
244  call get_param(paramfile, mdl, 'ENHANCE_DIFFUSION', cs%enhance_diffusion, &
245  'If True, adds enhanced diffusion at the based of the boundary layer.', &
246  default=.true.)
247  call get_param(paramfile, mdl, 'INTERP_TYPE', cs%interpType, &
248  'Type of interpolation to determine the OBL depth.\n'// &
249  'Allowed types are: linear, quadratic, cubic.', &
250  default='quadratic')
251  call get_param(paramfile, mdl, 'INTERP_TYPE2', cs%interpType2, &
252  'Type of interpolation to compute diff and visc at OBL_depth.\n'// &
253  'Allowed types are: linear, quadratic, cubic or LMD94.', &
254  default='LMD94')
255  call get_param(paramfile, mdl, 'COMPUTE_EKMAN', cs%computeEkman, &
256  'If True, limit OBL depth to be no deeper than Ekman depth.', &
257  default=.false.)
258  call get_param(paramfile, mdl, 'COMPUTE_MONIN_OBUKHOV', cs%computeMoninObukhov, &
259  'If True, limit the OBL depth to be no deeper than '// &
260  'Monin-Obukhov depth.', &
261  default=.false.)
262  call get_param(paramfile, mdl, 'CS', cs%cs, &
263  'Parameter for computing velocity scale function.', &
264  units='nondim', default=98.96)
265  call get_param(paramfile, mdl, 'CS2', cs%cs2, &
266  'Parameter for computing non-local term.', &
267  units='nondim', default=6.32739901508)
268  call get_param(paramfile, mdl, 'DEEP_OBL_OFFSET', cs%deepOBLoffset, &
269  'If non-zero, the distance above the bottom to which the OBL is clipped '// &
270  'if it would otherwise reach the bottom. The smaller of this and 0.1D is used.', &
271  units='m',default=0.)
272  call get_param(paramfile, mdl, 'FIXED_OBLDEPTH', cs%fixedOBLdepth, &
273  'If True, fix the OBL depth to FIXED_OBLDEPTH_VALUE '// &
274  'rather than using the OBL depth from CVMix. '// &
275  'This option is just for testing purposes.', &
276  default=.false.)
277  call get_param(paramfile, mdl, 'FIXED_OBLDEPTH_VALUE', cs%fixedOBLdepth_value, &
278  'Value for the fixed OBL depth when fixedOBLdepth==True. '// &
279  'This parameter is for just for testing purposes. '// &
280  'It will over-ride the OBLdepth computed from CVMix.', &
281  units='m',default=30.0)
282  call get_param(paramfile, mdl, 'SURF_LAYER_EXTENT', cs%surf_layer_ext, &
283  'Fraction of OBL depth considered in the surface layer.', &
284  units='nondim',default=0.10)
285  call get_param(paramfile, mdl, 'MINIMUM_OBL_DEPTH', cs%minOBLdepth, &
286  'If non-zero, a minimum depth to use for KPP OBL depth. Independent of '// &
287  'this parameter, the OBL depth is always at least as deep as the first layer.', &
288  units='m',default=0.)
289  call get_param(paramfile, mdl, 'MINIMUM_VT2', cs%minVtsqr, &
290  'Min of the unresolved velocity Vt2 used in Rib CVMix calculation.\n'// &
291  'Scaling: MINIMUM_VT2 = const1*d*N*ws, with d=1m, N=1e-5/s, ws=1e-6 m/s.', &
292  units='m2/s2',default=1e-10)
293 
294 ! smg: for removal below
295  call get_param(paramfile, mdl, 'CORRECT_SURFACE_LAYER_AVERAGE', cs%correctSurfLayerAvg, &
296  'If true, applies a correction step to the averaging of surface layer '// &
297  'properties. This option is obsolete.', default=.false.)
298  if (cs%correctSurfLayerAvg) &
299  call mom_error(fatal,'Correct surface layer average disabled in code. To recover \n'// &
300  ' feature will require code intervention.')
301  call get_param(paramfile, mdl, 'FIRST_GUESS_SURFACE_LAYER_DEPTH', cs%surfLayerDepth, &
302  'The first guess at the depth of the surface layer used for averaging '// &
303  'the surface layer properties. If =0, the top model level properties '// &
304  'will be used for the surface layer. If CORRECT_SURFACE_LAYER_AVERAGE=True, a '// &
305  'subsequent correction is applied. This parameter is obsolete', units='m', default=0.)
306 ! smg: for removal above
307 
308  call get_param(paramfile, mdl, 'NLT_SHAPE', string, &
309  'MOM6 method to set nonlocal transport profile. '// &
310  'Over-rides the result from CVMix. Allowed values are: \n'// &
311  '\t CVMix - Uses the profiles from CVMix specified by MATCH_TECHNIQUE\n'//&
312  '\t LINEAR - A linear profile, 1-sigma\n'// &
313  '\t PARABOLIC - A parablic profile, (1-sigma)^2\n'// &
314  '\t CUBIC - A cubic profile, (1-sigma)^2(1+2*sigma)\n'// &
315  '\t CUBIC_LMD - The original KPP profile', &
316  default='CVMix')
317  select case ( trim(string) )
318  case ("CVMix") ; cs%NLT_shape = nlt_shape_cvmix
319  case ("LINEAR") ; cs%NLT_shape = nlt_shape_linear
320  case ("PARABOLIC") ; cs%NLT_shape = nlt_shape_parabolic
321  case ("CUBIC") ; cs%NLT_shape = nlt_shape_cubic
322  case ("CUBIC_LMD") ; cs%NLT_shape = nlt_shape_cubic_lmd
323  case default ; call mom_error(fatal,"KPP_init: "// &
324  "Unrecognized NLT_SHAPE option"//trim(string))
325  end select
326  call get_param(paramfile, mdl, 'MATCH_TECHNIQUE', cs%MatchTechnique, &
327  'CVMix method to set profile function for diffusivity and NLT, '// &
328  'as well as matching across OBL base. Allowed values are: \n'// &
329  '\t SimpleShapes = sigma*(1-sigma)^2 for both diffusivity and NLT\n'// &
330  '\t MatchGradient = sigma*(1-sigma)^2 for NLT; diffusivity profile from matching\n'//&
331  '\t MatchBoth = match gradient for both diffusivity and NLT\n'// &
332  '\t ParabolicNonLocal = sigma*(1-sigma)^2 for diffusivity; (1-sigma)^2 for NLT', &
333  default='SimpleShapes')
334  if (cs%MatchTechnique == 'ParabolicNonLocal') then
335  ! This forces Cs2 (Cs in non-local computation) to equal 1 for parabolic non-local option.
336  ! May be used during CVMix initialization.
337  cs_is_one=.true.
338  endif
339  if (cs%MatchTechnique == 'ParabolicNonLocal' .or. cs%MatchTechnique == 'SimpleShapes') then
340  ! if gradient won't be matched, lnoDGat1=.true.
341  lnodgat1=.true.
342  endif
343 
344  ! safety check to avoid negative diff/visc
345  if (cs%MatchTechnique == 'MatchBoth' .and. (cs%interpType2 == 'cubic' .or. &
346  cs%interpType2 == 'quadratic')) then
347  call mom_error(fatal,"If MATCH_TECHNIQUE=MatchBoth, INTERP_TYPE2 must be set to \n"//&
348  "linear or LMD94 (recommended) to avoid negative viscosity and diffusivity.\n"//&
349  "Please select one of these valid options." )
350  endif
351 
352  call get_param(paramfile, mdl, 'KPP_ZERO_DIFFUSIVITY', cs%KPPzeroDiffusivity, &
353  'If True, zeroes the KPP diffusivity and viscosity; for testing purpose.',&
354  default=.false.)
355  call get_param(paramfile, mdl, 'KPP_IS_ADDITIVE', cs%KPPisAdditive, &
356  'If true, adds KPP diffusivity to diffusivity from other schemes.\n'//&
357  'If false, KPP is the only diffusivity wherever KPP is non-zero.', &
358  default=.true.)
359  call get_param(paramfile, mdl, 'KPP_SHORTWAVE_METHOD',string, &
360  'Determines contribution of shortwave radiation to KPP surface '// &
361  'buoyancy flux. Options include:\n'// &
362  ' ALL_SW: use total shortwave radiation\n'// &
363  ' MXL_SW: use shortwave radiation absorbed by mixing layer\n'// &
364  ' LV1_SW: use shortwave radiation absorbed by top model layer', &
365  default='MXL_SW')
366  select case ( trim(string) )
367  case ("ALL_SW") ; cs%SW_METHOD = sw_method_all_sw
368  case ("MXL_SW") ; cs%SW_METHOD = sw_method_mxl_sw
369  case ("LV1_SW") ; cs%SW_METHOD = sw_method_lv1_sw
370  case default ; call mom_error(fatal,"KPP_init: "// &
371  "Unrecognized KPP_SHORTWAVE_METHOD option"//trim(string))
372  end select
373  call get_param(paramfile, mdl, 'CVMix_ZERO_H_WORK_AROUND', cs%min_thickness, &
374  'A minimum thickness used to avoid division by small numbers in the vicinity '// &
375  'of vanished layers. This is independent of MIN_THICKNESS used in other parts of MOM.', &
376  units='m', default=0.)
377 
378 !/BGR: New options for including Langmuir effects
379 !/ 1. Options related to enhancing the mixing coefficient
380  call get_param(paramfile, mdl, "USE_KPP_LT_K", cs%LT_K_Enhancement, &
381  'Flag for Langmuir turbulence enhancement of turbulent'//&
382  'mixing coefficient.', units="", default=.false.)
383  call get_param(paramfile, mdl, "STOKES_MIXING", cs%STOKES_MIXING, &
384  'Flag for Langmuir turbulence enhancement of turbulent'//&
385  'mixing coefficient.', units="", default=.false.)
386  if (cs%LT_K_Enhancement) then
387  call get_param(paramfile, mdl, 'KPP_LT_K_SHAPE', string, &
388  'Vertical dependence of LT enhancement of mixing. '// &
389  'Valid options are: \n'// &
390  '\t CONSTANT = Constant value for full OBL\n'// &
391  '\t SCALED = Varies based on normalized shape function.', &
392  default='CONSTANT')
393  select case ( trim(string))
394  case ("CONSTANT") ; cs%LT_K_SHAPE = lt_k_constant
395  case ("SCALED") ; cs%LT_K_SHAPE = lt_k_scaled
396  case default ; call mom_error(fatal,"KPP_init: "//&
397  "Unrecognized KPP_LT_K_SHAPE option: "//trim(string))
398  end select
399  call get_param(paramfile, mdl, "KPP_LT_K_METHOD", string , &
400  'Method to enhance mixing coefficient in KPP. '// &
401  'Valid options are: \n'// &
402  '\t CONSTANT = Constant value (KPP_K_ENH_FAC) \n'// &
403  '\t VR12 = Function of Langmuir number based on VR12\n'// &
404  '\t RW16 = Function of Langmuir number based on RW16', &
405  default='CONSTANT')
406  select case ( trim(string))
407  case ("CONSTANT") ; cs%LT_K_METHOD = lt_k_mode_constant
408  case ("VR12") ; cs%LT_K_METHOD = lt_k_mode_vr12
409  case ("RW16") ; cs%LT_K_METHOD = lt_k_mode_rw16
410  case default ; call mom_error(fatal,"KPP_init: "//&
411  "Unrecognized KPP_LT_K_METHOD option: "//trim(string))
412  end select
413  if (cs%LT_K_METHOD==lt_k_mode_constant) then
414  call get_param(paramfile, mdl, "KPP_K_ENH_FAC",cs%KPP_K_ENH_FAC , &
415  'Constant value to enhance mixing coefficient in KPP.', &
416  default=1.0)
417  endif
418  endif
419 !/ 2. Options related to enhancing the unresolved Vt2/entrainment in Rib
420  call get_param(paramfile, mdl, "USE_KPP_LT_VT2", cs%LT_Vt2_Enhancement, &
421  'Flag for Langmuir turbulence enhancement of Vt2'//&
422  'in Bulk Richardson Number.', units="", default=.false.)
423  if (cs%LT_Vt2_Enhancement) then
424  call get_param(paramfile, mdl, "KPP_LT_VT2_METHOD",string , &
425  'Method to enhance Vt2 in KPP. '// &
426  'Valid options are: \n'// &
427  '\t CONSTANT = Constant value (KPP_VT2_ENH_FAC) \n'// &
428  '\t VR12 = Function of Langmuir number based on VR12\n'// &
429  '\t RW16 = Function of Langmuir number based on RW16\n'// &
430  '\t LF17 = Function of Langmuir number based on LF17', &
431  default='CONSTANT')
432  select case ( trim(string))
433  case ("CONSTANT") ; cs%LT_VT2_METHOD = lt_vt2_mode_constant
434  case ("VR12") ; cs%LT_VT2_METHOD = lt_vt2_mode_vr12
435  case ("RW16") ; cs%LT_VT2_METHOD = lt_vt2_mode_rw16
436  case ("LF17") ; cs%LT_VT2_METHOD = lt_vt2_mode_lf17
437  case default ; call mom_error(fatal,"KPP_init: "//&
438  "Unrecognized KPP_LT_VT2_METHOD option: "//trim(string))
439  end select
440  if (cs%LT_VT2_METHOD==lt_vt2_mode_constant) then
441  call get_param(paramfile, mdl, "KPP_VT2_ENH_FAC",cs%KPP_VT2_ENH_FAC , &
442  'Constant value to enhance VT2 in KPP.', &
443  default=1.0)
444  endif
445  endif
446 
447  call closeparameterblock(paramfile)
448  call get_param(paramfile, mdl, 'DEBUG', cs%debug, default=.false., do_not_log=.true.)
449 
450  call cvmix_init_kpp( ri_crit=cs%Ri_crit, &
451  minobldepth=cs%minOBLdepth, &
452  minvtsqr=cs%minVtsqr, &
453  vonkarman=cs%vonKarman, &
454  surf_layer_ext=cs%surf_layer_ext, &
455  interp_type=cs%interpType, &
456  interp_type2=cs%interpType2, &
457  lekman=cs%computeEkman, &
458  lmonob=cs%computeMoninObukhov, &
459  matchtechnique=cs%MatchTechnique, &
460  lenhanced_diff=cs%enhance_diffusion,&
461  lnonzero_surf_nonlocal=cs_is_one ,&
462  lnodgat1=lnodgat1 ,&
463  cvmix_kpp_params_user=cs%KPP_params )
464 
465  ! Register diagnostics
466  cs%diag => diag
467  cs%id_OBLdepth = register_diag_field('ocean_model', 'KPP_OBLdepth', diag%axesT1, time, &
468  'Thickness of the surface Ocean Boundary Layer calculated by [CVMix] KPP', 'meter', &
469  cmor_field_name='oml', cmor_long_name='ocean_mixed_layer_thickness_defined_by_mixing_scheme', &
470  cmor_units='m', cmor_standard_name='Ocean Mixed Layer Thickness Defined by Mixing Scheme')
471  ! CMOR names are placeholders; must be modified by time period
472  ! for CMOR compliance. Diag manager will be used for omlmax and
473  ! omldamax.
474  if (cs%n_smooth > 0) then
475  cs%id_OBLdepth_original = register_diag_field('ocean_model', 'KPP_OBLdepth_original', diag%axesT1, time, &
476  'Thickness of the surface Ocean Boundary Layer without smoothing calculated by [CVMix] KPP', 'meter', &
477  cmor_field_name='oml', cmor_long_name='ocean_mixed_layer_thickness_defined_by_mixing_scheme', &
478  cmor_units='m', cmor_standard_name='Ocean Mixed Layer Thickness Defined by Mixing Scheme')
479  endif
480  cs%id_BulkDrho = register_diag_field('ocean_model', 'KPP_BulkDrho', diag%axesTL, time, &
481  'Bulk difference in density used in Bulk Richardson number, as used by [CVMix] KPP', 'kg/m3')
482  cs%id_BulkUz2 = register_diag_field('ocean_model', 'KPP_BulkUz2', diag%axesTL, time, &
483  'Square of bulk difference in resolved velocity used in Bulk Richardson number via [CVMix] KPP', 'm2/s2')
484  cs%id_BulkRi = register_diag_field('ocean_model', 'KPP_BulkRi', diag%axesTL, time, &
485  'Bulk Richardson number used to find the OBL depth used by [CVMix] KPP', 'nondim')
486  cs%id_Sigma = register_diag_field('ocean_model', 'KPP_sigma', diag%axesTi, time, &
487  'Sigma coordinate used by [CVMix] KPP', 'nondim')
488  cs%id_Ws = register_diag_field('ocean_model', 'KPP_Ws', diag%axesTL, time, &
489  'Turbulent vertical velocity scale for scalars used by [CVMix] KPP', 'm/s')
490  cs%id_N = register_diag_field('ocean_model', 'KPP_N', diag%axesTi, time, &
491  '(Adjusted) Brunt-Vaisala frequency used by [CVMix] KPP', '1/s')
492  cs%id_N2 = register_diag_field('ocean_model', 'KPP_N2', diag%axesTi, time, &
493  'Square of Brunt-Vaisala frequency used by [CVMix] KPP', '1/s2')
494  cs%id_Vt2 = register_diag_field('ocean_model', 'KPP_Vt2', diag%axesTL, time, &
495  'Unresolved shear turbulence used by [CVMix] KPP', 'm2/s2')
496  cs%id_uStar = register_diag_field('ocean_model', 'KPP_uStar', diag%axesT1, time, &
497  'Friction velocity, u*, as used by [CVMix] KPP', 'm/s', conversion=us%Z_to_m*us%s_to_T)
498  cs%id_buoyFlux = register_diag_field('ocean_model', 'KPP_buoyFlux', diag%axesTi, time, &
499  'Surface (and penetrating) buoyancy flux, as used by [CVMix] KPP', 'm2/s3', conversion=us%L_to_m**2*us%s_to_T**3)
500  cs%id_QminusSW = register_diag_field('ocean_model', 'KPP_QminusSW', diag%axesT1, time, &
501  'Net temperature flux ignoring short-wave, as used by [CVMix] KPP', 'K m/s')
502  cs%id_netS = register_diag_field('ocean_model', 'KPP_netSalt', diag%axesT1, time, &
503  'Effective net surface salt flux, as used by [CVMix] KPP', 'ppt m/s')
504  cs%id_Kt_KPP = register_diag_field('ocean_model', 'KPP_Kheat', diag%axesTi, time, &
505  'Heat diffusivity due to KPP, as calculated by [CVMix] KPP', 'm2/s')
506  cs%id_Kd_in = register_diag_field('ocean_model', 'KPP_Kd_in', diag%axesTi, time, &
507  'Diffusivity passed to KPP', 'm2/s', conversion=us%Z2_T_to_m2_s)
508  cs%id_Ks_KPP = register_diag_field('ocean_model', 'KPP_Ksalt', diag%axesTi, time, &
509  'Salt diffusivity due to KPP, as calculated by [CVMix] KPP', 'm2/s')
510  cs%id_Kv_KPP = register_diag_field('ocean_model', 'KPP_Kv', diag%axesTi, time, &
511  'Vertical viscosity due to KPP, as calculated by [CVMix] KPP', 'm2/s')
512  cs%id_NLTt = register_diag_field('ocean_model', 'KPP_NLtransport_heat', diag%axesTi, time, &
513  'Non-local transport (Cs*G(sigma)) for heat, as calculated by [CVMix] KPP', 'nondim')
514  cs%id_NLTs = register_diag_field('ocean_model', 'KPP_NLtransport_salt', diag%axesTi, time, &
515  'Non-local tranpsort (Cs*G(sigma)) for scalars, as calculated by [CVMix] KPP', 'nondim')
516  cs%id_NLT_dTdt = register_diag_field('ocean_model', 'KPP_NLT_dTdt', diag%axesTL, time, &
517  'Temperature tendency due to non-local transport of heat, as calculated by [CVMix] KPP', 'K/s')
518  cs%id_NLT_dSdt = register_diag_field('ocean_model', 'KPP_NLT_dSdt', diag%axesTL, time, &
519  'Salinity tendency due to non-local transport of salt, as calculated by [CVMix] KPP', 'ppt/s')
520  cs%id_NLT_temp_budget = register_diag_field('ocean_model', 'KPP_NLT_temp_budget', diag%axesTL, time, &
521  'Heat content change due to non-local transport, as calculated by [CVMix] KPP', 'W/m^2')
522  cs%id_NLT_saln_budget = register_diag_field('ocean_model', 'KPP_NLT_saln_budget', diag%axesTL, time, &
523  'Salt content change due to non-local transport, as calculated by [CVMix] KPP', 'kg/(sec*m^2)')
524  cs%id_Tsurf = register_diag_field('ocean_model', 'KPP_Tsurf', diag%axesT1, time, &
525  'Temperature of surface layer (10% of OBL depth) as passed to [CVMix] KPP', 'C')
526  cs%id_Ssurf = register_diag_field('ocean_model', 'KPP_Ssurf', diag%axesT1, time, &
527  'Salinity of surface layer (10% of OBL depth) as passed to [CVMix] KPP', 'ppt')
528  cs%id_Usurf = register_diag_field('ocean_model', 'KPP_Usurf', diag%axesCu1, time, &
529  'i-component flow of surface layer (10% of OBL depth) as passed to [CVMix] KPP', 'm/s')
530  cs%id_Vsurf = register_diag_field('ocean_model', 'KPP_Vsurf', diag%axesCv1, time, &
531  'j-component flow of surface layer (10% of OBL depth) as passed to [CVMix] KPP', 'm/s')
532  cs%id_EnhK = register_diag_field('ocean_model', 'EnhK', diag%axesTI, time, &
533  'Langmuir number enhancement to K as used by [CVMix] KPP','nondim')
534  cs%id_EnhVt2 = register_diag_field('ocean_model', 'EnhVt2', diag%axesTL, time, &
535  'Langmuir number enhancement to Vt2 as used by [CVMix] KPP','nondim')
536  cs%id_La_SL = register_diag_field('ocean_model', 'KPP_La_SL', diag%axesT1, time, &
537  'Surface-layer Langmuir number computed in [CVMix] KPP','nondim')
538 
539  allocate( cs%N( szi_(g), szj_(g), szk_(g)+1 ) )
540  cs%N(:,:,:) = 0.
541  allocate( cs%OBLdepth( szi_(g), szj_(g) ) )
542  cs%OBLdepth(:,:) = 0.
543  allocate( cs%kOBL( szi_(g), szj_(g) ) )
544  cs%kOBL(:,:) = 0.
545  allocate( cs%La_SL( szi_(g), szj_(g) ) )
546  cs%La_SL(:,:) = 0.
547  allocate( cs%Vt2( szi_(g), szj_(g), szk_(g) ) )
548  cs%Vt2(:,:,:) = 0.
549  if (cs%id_OBLdepth_original > 0) allocate( cs%OBLdepth_original( szi_(g), szj_(g) ) )
550 
551  allocate( cs%OBLdepthprev( szi_(g), szj_(g) ) ) ; cs%OBLdepthprev(:,:) = 0.0
552  if (cs%id_BulkDrho > 0) allocate( cs%dRho( szi_(g), szj_(g), szk_(g) ) )
553  if (cs%id_BulkDrho > 0) cs%dRho(:,:,:) = 0.
554  if (cs%id_BulkUz2 > 0) allocate( cs%Uz2( szi_(g), szj_(g), szk_(g) ) )
555  if (cs%id_BulkUz2 > 0) cs%Uz2(:,:,:) = 0.
556  if (cs%id_BulkRi > 0) allocate( cs%BulkRi( szi_(g), szj_(g), szk_(g) ) )
557  if (cs%id_BulkRi > 0) cs%BulkRi(:,:,:) = 0.
558  if (cs%id_Sigma > 0) allocate( cs%sigma( szi_(g), szj_(g), szk_(g)+1 ) )
559  if (cs%id_Sigma > 0) cs%sigma(:,:,:) = 0.
560  if (cs%id_Ws > 0) allocate( cs%Ws( szi_(g), szj_(g), szk_(g) ) )
561  if (cs%id_Ws > 0) cs%Ws(:,:,:) = 0.
562  if (cs%id_N2 > 0) allocate( cs%N2( szi_(g), szj_(g), szk_(g)+1 ) )
563  if (cs%id_N2 > 0) cs%N2(:,:,:) = 0.
564  if (cs%id_Kt_KPP > 0) allocate( cs%Kt_KPP( szi_(g), szj_(g), szk_(g)+1 ) )
565  if (cs%id_Kt_KPP > 0) cs%Kt_KPP(:,:,:) = 0.
566  if (cs%id_Ks_KPP > 0) allocate( cs%Ks_KPP( szi_(g), szj_(g), szk_(g)+1 ) )
567  if (cs%id_Ks_KPP > 0) cs%Ks_KPP(:,:,:) = 0.
568  if (cs%id_Kv_KPP > 0) allocate( cs%Kv_KPP( szi_(g), szj_(g), szk_(g)+1 ) )
569  if (cs%id_Kv_KPP > 0) cs%Kv_KPP(:,:,:) = 0.
570  if (cs%id_Tsurf > 0) allocate( cs%Tsurf( szi_(g), szj_(g)) )
571  if (cs%id_Tsurf > 0) cs%Tsurf(:,:) = 0.
572  if (cs%id_Ssurf > 0) allocate( cs%Ssurf( szi_(g), szj_(g)) )
573  if (cs%id_Ssurf > 0) cs%Ssurf(:,:) = 0.
574  if (cs%id_Usurf > 0) allocate( cs%Usurf( szib_(g), szj_(g)) )
575  if (cs%id_Usurf > 0) cs%Usurf(:,:) = 0.
576  if (cs%id_Vsurf > 0) allocate( cs%Vsurf( szi_(g), szjb_(g)) )
577  if (cs%id_Vsurf > 0) cs%Vsurf(:,:) = 0.
578  if (cs%id_EnhVt2 > 0) allocate( cs%EnhVt2( szi_(g), szj_(g), szk_(g)) )
579  if (cs%id_EnhVt2 > 0) cs%EnhVt2(:,:,:) = 0.
580  if (cs%id_EnhK > 0) allocate( cs%EnhK( szi_(g), szj_(g), szk_(g)+1 ) )
581  if (cs%id_EnhK > 0) cs%EnhK(:,:,:) = 0.
582 
583 
584 end function kpp_init
585 
586 !> KPP vertical diffusivity/viscosity and non-local tracer transport
587 subroutine kpp_calculate(CS, G, GV, US, h, uStar, &
588  buoyFlux, Kt, Ks, Kv, nonLocalTransHeat,&
589  nonLocalTransScalar, waves)
590 
591  ! Arguments
592  type(kpp_cs), pointer :: cs !< Control structure
593  type(ocean_grid_type), intent(in) :: g !< Ocean grid
594  type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid
595  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
596  type(wave_parameters_cs), optional, pointer :: waves !< Wave CS
597  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer/level thicknesses [H ~> m or kg m-2]
598  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: ustar !< Surface friction velocity [Z T-1 ~> m s-1]
599  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in) :: buoyflux !< Surface buoyancy flux [L2 T-3 ~> m2 s-3]
600  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: kt !< (in) Vertical diffusivity of heat w/o KPP
601  !! (out) Vertical diffusivity including KPP
602  !! [Z2 T-1 ~> m2 s-1]
603  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: ks !< (in) Vertical diffusivity of salt w/o KPP
604  !! (out) Vertical diffusivity including KPP
605  !! [Z2 T-1 ~> m2 s-1]
606  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: kv !< (in) Vertical viscosity w/o KPP
607  !! (out) Vertical viscosity including KPP
608  !! [Z2 T-1 ~> m2 s-1]
609  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: nonlocaltransheat !< Temp non-local transport [m s-1]
610  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: nonlocaltransscalar !< scalar non-local transport [m s-1]
611 
612 ! Local variables
613  integer :: i, j, k ! Loop indices
614  real, dimension( G%ke ) :: cellheight ! Cell center heights referenced to surface [m] (negative in ocean)
615  real, dimension( G%ke+1 ) :: ifaceheight ! Interface heights referenced to surface [m] (negative in ocean)
616  real, dimension( G%ke+1, 2) :: kdiffusivity ! Vertical diffusivity at interfaces [m2 s-1]
617  real, dimension( G%ke+1 ) :: kviscosity ! Vertical viscosity at interfaces [m2 s-1]
618  real, dimension( G%ke+1, 2) :: nonlocaltrans ! Non-local transport for heat/salt at interfaces [nondim]
619 
620  real :: surffricvel, surfbuoyflux
621  real :: sigma, sigmaratio
622  real :: buoy_scale ! A unit conversion factor for buoyancy fluxes [m2 T3 L-2 s-3 ~> nondim]
623  real :: dh ! The local thickness used for calculating interface positions [m]
624  real :: hcorr ! A cumulative correction arising from inflation of vanished layers [m]
625 
626  ! For Langmuir Calculations
627  real :: langenhk ! Langmuir enhancement for mixing coefficient
628 
629 
630 #ifdef __DO_SAFETY_CHECKS__
631  if (cs%debug) then
632  call hchksum(h, "KPP in: h",g%HI,haloshift=0, scale=gv%H_to_m)
633  call hchksum(ustar, "KPP in: uStar",g%HI,haloshift=0, scale=us%Z_to_m*us%s_to_T)
634  call hchksum(buoyflux, "KPP in: buoyFlux",g%HI,haloshift=0)
635  call hchksum(kt, "KPP in: Kt",g%HI,haloshift=0, scale=us%Z2_T_to_m2_s)
636  call hchksum(ks, "KPP in: Ks",g%HI,haloshift=0, scale=us%Z2_T_to_m2_s)
637  endif
638 #endif
639 
640  nonlocaltrans(:,:) = 0.0
641 
642  if (cs%id_Kd_in > 0) call post_data(cs%id_Kd_in, kt, cs%diag)
643 
644  buoy_scale = us%L_to_m**2*us%s_to_T**3
645 
646  !$OMP parallel do default(none) firstprivate(nonLocalTrans) &
647  !$OMP private(surfFricVel, iFaceHeight, hcorr, dh, cellHeight, &
648  !$OMP surfBuoyFlux, Kdiffusivity, Kviscosity, LangEnhK, sigma, &
649  !$OMP sigmaRatio) &
650  !$OMP shared(G, GV, CS, US, uStar, h, buoy_scale, buoyFlux, Kt, &
651  !$OMP Ks, Kv, nonLocalTransHeat, nonLocalTransScalar, waves)
652  ! loop over horizontal points on processor
653  do j = g%jsc, g%jec
654  do i = g%isc, g%iec
655 
656  ! skip calling KPP for land points
657  if (g%mask2dT(i,j)==0.) cycle
658 
659  ! things independent of position within the column
660  surffricvel = us%Z_to_m*us%s_to_T * ustar(i,j)
661 
662  ifaceheight(1) = 0.0 ! BBL is all relative to the surface
663  hcorr = 0.
664  do k=1,g%ke
665 
666  ! cell center and cell bottom in meters (negative values in the ocean)
667  dh = h(i,j,k) * gv%H_to_m ! Nominal thickness to use for increment
668  dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0)
669  hcorr = min( dh - cs%min_thickness, 0. ) ! If inflating then hcorr<0
670  dh = max( dh, cs%min_thickness ) ! Limit increment dh>=min_thickness
671  cellheight(k) = ifaceheight(k) - 0.5 * dh
672  ifaceheight(k+1) = ifaceheight(k) - dh
673 
674  enddo ! k-loop finishes
675 
676  surfbuoyflux = buoy_scale*buoyflux(i,j,1) ! This is only used in kpp_compute_OBL_depth to limit
677  ! h to Monin-Obukov (default is false, ie. not used)
678 
679  ! Call CVMix/KPP to obtain OBL diffusivities, viscosities and non-local transports
680 
681  ! Unlike LMD94, we do not match to interior diffusivities. If using the original
682  ! LMD94 shape function, not matching is equivalent to matching to a zero diffusivity.
683 
684  !BGR/ Add option for use of surface buoyancy flux with total sw flux.
685  if (cs%SW_METHOD == sw_method_all_sw) then
686  surfbuoyflux = buoy_scale * buoyflux(i,j,1)
687  elseif (cs%SW_METHOD == sw_method_mxl_sw) then
688  ! We know the actual buoyancy flux into the OBL
689  surfbuoyflux = buoy_scale * (buoyflux(i,j,1) - buoyflux(i,j,int(cs%kOBL(i,j))+1))
690  elseif (cs%SW_METHOD == sw_method_lv1_sw) then
691  surfbuoyflux = buoy_scale * (buoyflux(i,j,1) - buoyflux(i,j,2))
692  endif
693 
694  ! If option "MatchBoth" is selected in CVMix, MOM should be capable of matching.
695  if (.not. (cs%MatchTechnique == 'MatchBoth')) then
696  kdiffusivity(:,:) = 0. ! Diffusivities for heat and salt [m2 s-1]
697  kviscosity(:) = 0. ! Viscosity [m2 s-1]
698  else
699  kdiffusivity(:,1) = us%Z2_T_to_m2_s * kt(i,j,:)
700  kdiffusivity(:,2) = us%Z2_T_to_m2_s * ks(i,j,:)
701  kviscosity(:) = us%Z2_T_to_m2_s * kv(i,j,:)
702  endif
703 
704  call cvmix_coeffs_kpp(kviscosity(:), & ! (inout) Total viscosity [m2 s-1]
705  kdiffusivity(:,1), & ! (inout) Total heat diffusivity [m2 s-1]
706  kdiffusivity(:,2), & ! (inout) Total salt diffusivity [m2 s-1]
707  ifaceheight, & ! (in) Height of interfaces [m]
708  cellheight, & ! (in) Height of level centers [m]
709  kviscosity(:), & ! (in) Original viscosity [m2 s-1]
710  kdiffusivity(:,1), & ! (in) Original heat diffusivity [m2 s-1]
711  kdiffusivity(:,2), & ! (in) Original salt diffusivity [m2 s-1]
712  cs%OBLdepth(i,j), & ! (in) OBL depth [m]
713  cs%kOBL(i,j), & ! (in) level (+fraction) of OBL extent
714  nonlocaltrans(:,1),& ! (out) Non-local heat transport [nondim]
715  nonlocaltrans(:,2),& ! (out) Non-local salt transport [nondim]
716  surffricvel, & ! (in) Turbulent friction velocity at surface [m s-1]
717  surfbuoyflux, & ! (in) Buoyancy flux at surface [m2 s-3]
718  g%ke, & ! (in) Number of levels to compute coeffs for
719  g%ke, & ! (in) Number of levels in array shape
720  cvmix_kpp_params_user=cs%KPP_params )
721 
722  ! safety check, Kviscosity and Kdiffusivity must be >= 0
723  do k=1, g%ke+1
724  if (kviscosity(k) < 0. .or. kdiffusivity(k,1) < 0.) then
725  call mom_error(fatal,"KPP_calculate, after CVMix_coeffs_kpp: "// &
726  "Negative vertical viscosity or diffusivity has been detected. " // &
727  "This is likely related to the choice of MATCH_TECHNIQUE and INTERP_TYPE2." //&
728  "You might consider using the default options for these parameters." )
729  endif
730  enddo
731 
732  IF (cs%LT_K_ENHANCEMENT) then
733  if (cs%LT_K_METHOD==lt_k_mode_constant) then
734  langenhk = cs%KPP_K_ENH_FAC
735  elseif (cs%LT_K_METHOD==lt_k_mode_vr12) then
736  ! Added minimum value for La_SL, so removed maximum value for LangEnhK.
737  langenhk = sqrt(1.+(1.5*cs%La_SL(i,j))**(-2) + &
738  (5.4*cs%La_SL(i,j))**(-4))
739  elseif (cs%LT_K_METHOD==lt_k_mode_rw16) then
740  !This maximum value is proposed in Reichl et al., 2016 JPO formula
741  langenhk = min(2.25, 1. + 1./cs%La_SL(i,j))
742  else
743  !This shouldn't be reached.
744  !call MOM_error(WARNING,"Unexpected behavior in MOM_CVMix_KPP, see error in LT_K_ENHANCEMENT")
745  langenhk = 1.0
746  endif
747  do k=1,g%ke
748  if (cs%LT_K_SHAPE== lt_k_constant) then
749  if (cs%id_EnhK > 0) cs%EnhK(i,j,:) = langenhk
750  kdiffusivity(k,1) = kdiffusivity(k,1) * langenhk
751  kdiffusivity(k,2) = kdiffusivity(k,2) * langenhk
752  kviscosity(k) = kviscosity(k) * langenhk
753  elseif (cs%LT_K_SHAPE == lt_k_scaled) then
754  sigma = min(1.0,-ifaceheight(k)/cs%OBLdepth(i,j))
755  sigmaratio = sigma * (1. - sigma)**2 / 0.148148037
756  if (cs%id_EnhK > 0) cs%EnhK(i,j,k) = (1.0 + (langenhk - 1.)*sigmaratio)
757  kdiffusivity(k,1) = kdiffusivity(k,1) * ( 1. + &
758  ( langenhk - 1.)*sigmaratio)
759  kdiffusivity(k,2) = kdiffusivity(k,2) * ( 1. + &
760  ( langenhk - 1.)*sigmaratio)
761  kviscosity(k) = kviscosity(k) * ( 1. + &
762  ( langenhk - 1.)*sigmaratio)
763  endif
764  enddo
765  endif
766 
767  ! Over-write CVMix NLT shape function with one of the following choices.
768  ! The CVMix code has yet to update for thse options, so we compute in MOM6.
769  ! Note that nonLocalTrans = Cs * G(sigma) (LMD94 notation), with
770  ! Cs = 6.32739901508.
771  ! Start do-loop at k=2, since k=1 is ocean surface (sigma=0)
772  ! and we do not wish to double-count the surface forcing.
773  ! Only compute nonlocal transport for 0 <= sigma <= 1.
774  ! MOM6 recommended shape is the parabolic; it gives deeper boundary layer
775  ! and no spurious extrema.
776  if (surfbuoyflux < 0.0) then
777  if (cs%NLT_shape == nlt_shape_cubic) then
778  do k = 2, g%ke
779  sigma = min(1.0,-ifaceheight(k)/cs%OBLdepth(i,j))
780  nonlocaltrans(k,1) = (1.0 - sigma)**2 * (1.0 + 2.0*sigma) !*
781  nonlocaltrans(k,2) = nonlocaltrans(k,1)
782  enddo
783  elseif (cs%NLT_shape == nlt_shape_parabolic) then
784  do k = 2, g%ke
785  sigma = min(1.0,-ifaceheight(k)/cs%OBLdepth(i,j))
786  nonlocaltrans(k,1) = (1.0 - sigma)**2 !*CS%CS2
787  nonlocaltrans(k,2) = nonlocaltrans(k,1)
788  enddo
789  elseif (cs%NLT_shape == nlt_shape_linear) then
790  do k = 2, g%ke
791  sigma = min(1.0,-ifaceheight(k)/cs%OBLdepth(i,j))
792  nonlocaltrans(k,1) = (1.0 - sigma)!*CS%CS2
793  nonlocaltrans(k,2) = nonlocaltrans(k,1)
794  enddo
795  elseif (cs%NLT_shape == nlt_shape_cubic_lmd) then
796  ! Sanity check (should agree with CVMix result using simple matching)
797  do k = 2, g%ke
798  sigma = min(1.0,-ifaceheight(k)/cs%OBLdepth(i,j))
799  nonlocaltrans(k,1) = cs%CS2 * sigma*(1.0 -sigma)**2
800  nonlocaltrans(k,2) = nonlocaltrans(k,1)
801  enddo
802  endif
803  endif
804 
805  ! we apply nonLocalTrans in subroutines
806  ! KPP_NonLocalTransport_temp and KPP_NonLocalTransport_saln
807  nonlocaltransheat(i,j,:) = nonlocaltrans(:,1) ! temp
808  nonlocaltransscalar(i,j,:) = nonlocaltrans(:,2) ! saln
809 
810  ! set the KPP diffusivity and viscosity to zero for testing purposes
811  if (cs%KPPzeroDiffusivity) then
812  kdiffusivity(:,1) = 0.0
813  kdiffusivity(:,2) = 0.0
814  kviscosity(:) = 0.0
815  endif
816 
817 
818  ! compute unresolved squared velocity for diagnostics
819  if (cs%id_Vt2 > 0) then
820 !BGR Now computing VT2 above so can modify for LT
821 ! therefore, don't repeat this operation here
822 ! CS%Vt2(i,j,:) = CVmix_kpp_compute_unresolved_shear( &
823 ! cellHeight(1:G%ke), & ! Depth of cell center [m]
824 ! ws_cntr=Ws_1d, & ! Turbulent velocity scale profile, at centers [m s-1]
825 ! N_iface=CS%N(i,j,:), & ! Buoyancy frequency at interface [s-1]
826 ! CVmix_kpp_params_user=CS%KPP_params ) ! KPP parameters
827  endif
828 
829  ! Copy 1d data into 3d diagnostic arrays
830  !/ grabbing obldepth_0d for next time step.
831  cs%OBLdepthprev(i,j)=cs%OBLdepth(i,j)
832  if (cs%id_sigma > 0) then
833  cs%sigma(i,j,:) = 0.
834  if (cs%OBLdepth(i,j)>0.) cs%sigma(i,j,:) = -ifaceheight/cs%OBLdepth(i,j)
835  endif
836  if (cs%id_Kt_KPP > 0) cs%Kt_KPP(i,j,:) = kdiffusivity(:,1)
837  if (cs%id_Ks_KPP > 0) cs%Ks_KPP(i,j,:) = kdiffusivity(:,2)
838  if (cs%id_Kv_KPP > 0) cs%Kv_KPP(i,j,:) = kviscosity(:)
839 
840  ! Update output of routine
841  if (.not. cs%passiveMode) then
842  if (cs%KPPisAdditive) then
843  do k=1, g%ke+1
844  kt(i,j,k) = kt(i,j,k) + us%m2_s_to_Z2_T * kdiffusivity(k,1)
845  ks(i,j,k) = ks(i,j,k) + us%m2_s_to_Z2_T * kdiffusivity(k,2)
846  kv(i,j,k) = kv(i,j,k) + us%m2_s_to_Z2_T * kviscosity(k)
847  if (cs%Stokes_Mixing) waves%KvS(i,j,k) = kv(i,j,k)
848  enddo
849  else ! KPP replaces prior diffusivity when former is non-zero
850  do k=1, g%ke+1
851  if (kdiffusivity(k,1) /= 0.) kt(i,j,k) = us%m2_s_to_Z2_T * kdiffusivity(k,1)
852  if (kdiffusivity(k,2) /= 0.) ks(i,j,k) = us%m2_s_to_Z2_T * kdiffusivity(k,2)
853  if (kviscosity(k) /= 0.) kv(i,j,k) = us%m2_s_to_Z2_T * kviscosity(k)
854  if (cs%Stokes_Mixing) waves%KvS(i,j,k) = kv(i,j,k)
855  enddo
856  endif
857  endif
858 
859 
860  ! end of the horizontal do-loops over the vertical columns
861  enddo ! i
862  enddo ! j
863 
864 
865 #ifdef __DO_SAFETY_CHECKS__
866  if (cs%debug) then
867  call hchksum(kt, "KPP out: Kt", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
868  call hchksum(ks, "KPP out: Ks", g%HI, haloshift=0, scale=us%Z2_T_to_m2_s)
869  endif
870 #endif
871 
872  ! send diagnostics to post_data
873  if (cs%id_OBLdepth > 0) call post_data(cs%id_OBLdepth, cs%OBLdepth, cs%diag)
874  if (cs%id_OBLdepth_original > 0) call post_data(cs%id_OBLdepth_original,cs%OBLdepth_original,cs%diag)
875  if (cs%id_sigma > 0) call post_data(cs%id_sigma, cs%sigma, cs%diag)
876  if (cs%id_Ws > 0) call post_data(cs%id_Ws, cs%Ws, cs%diag)
877  if (cs%id_Vt2 > 0) call post_data(cs%id_Vt2, cs%Vt2, cs%diag)
878  if (cs%id_uStar > 0) call post_data(cs%id_uStar, ustar, cs%diag)
879  if (cs%id_buoyFlux > 0) call post_data(cs%id_buoyFlux, buoyflux, cs%diag)
880  if (cs%id_Kt_KPP > 0) call post_data(cs%id_Kt_KPP, cs%Kt_KPP, cs%diag)
881  if (cs%id_Ks_KPP > 0) call post_data(cs%id_Ks_KPP, cs%Ks_KPP, cs%diag)
882  if (cs%id_Kv_KPP > 0) call post_data(cs%id_Kv_KPP, cs%Kv_KPP, cs%diag)
883  if (cs%id_NLTt > 0) call post_data(cs%id_NLTt, nonlocaltransheat, cs%diag)
884  if (cs%id_NLTs > 0) call post_data(cs%id_NLTs, nonlocaltransscalar,cs%diag)
885 
886 
887 end subroutine kpp_calculate
888 
889 
890 !> Compute OBL depth
891 subroutine kpp_compute_bld(CS, G, GV, US, h, Temp, Salt, u, v, EOS, uStar, buoyFlux, Waves)
892 
893  ! Arguments
894  type(kpp_cs), pointer :: cs !< Control structure
895  type(ocean_grid_type), intent(inout) :: g !< Ocean grid
896  type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid
897  type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
898  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer/level thicknesses [H ~> m or kg m-2]
899  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: temp !< potential/cons temp [degC]
900  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: salt !< Salinity [ppt]
901  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: u !< Velocity i-component [L T-1 ~> m s-1]
902  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: v !< Velocity j-component [L T-1 ~> m s-1]
903  type(eos_type), pointer :: eos !< Equation of state
904  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: ustar !< Surface friction velocity [Z T-1 ~> m s-1]
905  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in) :: buoyflux !< Surface buoyancy flux [L2 T-3 ~> m2 s-3]
906  type(wave_parameters_cs), optional, pointer :: waves !< Wave CS
907 
908  ! Local variables
909  integer :: i, j, k, km1 ! Loop indices
910  real, dimension( G%ke ) :: cellheight ! Cell center heights referenced to surface [m] (negative in ocean)
911  real, dimension( G%ke+1 ) :: ifaceheight ! Interface heights referenced to surface [m] (negative in ocean)
912  real, dimension( G%ke+1 ) :: n2_1d ! Brunt-Vaisala frequency squared, at interfaces [s-2]
913  real, dimension( G%ke ) :: ws_1d ! Profile of vertical velocity scale for scalars [m s-1]
914  real, dimension( G%ke ) :: deltarho ! delta Rho in numerator of Bulk Ri number
915  real, dimension( G%ke ) :: deltau2 ! square of delta U (shear) in denominator of Bulk Ri [m2 s-2]
916  real, dimension( G%ke ) :: surfbuoyflux2
917  real, dimension( G%ke ) :: bulkri_1d ! Bulk Richardson number for each layer
918 
919  ! for EOS calculation
920  real, dimension( 3*G%ke ) :: rho_1d
921  real, dimension( 3*G%ke ) :: pres_1d
922  real, dimension( 3*G%ke ) :: temp_1d
923  real, dimension( 3*G%ke ) :: salt_1d
924 
925  real :: surffricvel, surfbuoyflux, coriolis
926  real :: gorho ! Gravitational acceleration divided by density in MKS units [m4 s-2]
927  real :: pref, rho1, rhok, uk, vk, sigma, sigmaratio
928 
929  real :: zbottomminusoffset ! Height of bottom plus a little bit [m]
930  real :: sldepth_0d ! Surface layer depth = surf_layer_ext*OBLdepth.
931  real :: htot ! Running sum of thickness used in the surface layer average [m]
932  real :: buoy_scale ! A unit conversion factor for buoyancy fluxes [m2 T3 L-2 s-3 ~> nondim]
933  real :: delh ! Thickness of a layer [m]
934  real :: surfhtemp, surftemp ! Integral and average of temp over the surface layer
935  real :: surfhsalt, surfsalt ! Integral and average of saln over the surface layer
936  real :: surfhu, surfu ! Integral and average of u over the surface layer
937  real :: surfhv, surfv ! Integral and average of v over the surface layer
938  real :: dh ! The local thickness used for calculating interface positions [m]
939  real :: hcorr ! A cumulative correction arising from inflation of vanished layers [m]
940  integer :: kk, ksfc, ktmp
941 
942  ! For Langmuir Calculations
943  real :: langenhw ! Langmuir enhancement for turbulent velocity scale
944  real, dimension(G%ke) :: langenhvt2 ! Langmuir enhancement for unresolved shear
945  real, dimension(G%ke) :: u_h, v_h
946  real :: mld_guess, la
947  real :: surfhus, surfhvs, surfus, surfvs, wavedir, currentdir
948  real :: varup, vardn, m, varlo, varavg
949  real :: h10pct, h20pct,cmnfact, usx20pct, usy20pct, enhvt2
950  integer :: b
951  real :: wst
952 
953 
954 #ifdef __DO_SAFETY_CHECKS__
955  if (cs%debug) then
956  call hchksum(salt, "KPP in: S",g%HI,haloshift=0)
957  call hchksum(temp, "KPP in: T",g%HI,haloshift=0)
958  call hchksum(u, "KPP in: u",g%HI,haloshift=0)
959  call hchksum(v, "KPP in: v",g%HI,haloshift=0)
960  endif
961 #endif
962 
963  ! some constants
964  gorho = gv%mks_g_Earth / (us%R_to_kg_m3*gv%Rho0)
965  buoy_scale = us%L_to_m**2*us%s_to_T**3
966 
967  ! loop over horizontal points on processor
968  !$OMP parallel do default(none) private(surfFricVel, iFaceHeight, hcorr, dh, cellHeight, &
969  !$OMP surfBuoyFlux, U_H, V_H, u, v, Coriolis, pRef, SLdepth_0d, &
970  !$OMP ksfc, surfHtemp, surfHsalt, surfHu, surfHv, surfHuS, &
971  !$OMP surfHvS, hTot, delH, surftemp, surfsalt, surfu, surfv, &
972  !$OMP surfUs, surfVs, Uk, Vk, deltaU2, km1, kk, pres_1D, &
973  !$OMP Temp_1D, salt_1D, surfBuoyFlux2, MLD_GUESS, LA, rho_1D, &
974  !$OMP deltarho, N2_1d, ws_1d, LangEnhVT2, enhvt2, wst, &
975  !$OMP BulkRi_1d, zBottomMinusOffset) &
976  !$OMP shared(G, GV, CS, US, uStar, h, buoy_scale, buoyFlux, &
977  !$OMP Temp, Salt, waves, EOS, GoRho)
978  do j = g%jsc, g%jec
979  do i = g%isc, g%iec
980 
981  ! skip calling KPP for land points
982  if (g%mask2dT(i,j)==0.) cycle
983 
984  do k=1,g%ke
985  u_h(k) = 0.5 * us%L_T_to_m_s*(u(i,j,k)+u(i-1,j,k))
986  v_h(k) = 0.5 * us%L_T_to_m_s*(v(i,j,k)+v(i,j-1,k))
987  enddo
988 
989  ! things independent of position within the column
990  coriolis = 0.25*us%s_to_T*( (g%CoriolisBu(i,j) + g%CoriolisBu(i-1,j-1)) + &
991  (g%CoriolisBu(i-1,j) + g%CoriolisBu(i,j-1)) )
992  surffricvel = us%Z_to_m*us%s_to_T * ustar(i,j)
993 
994  ! Bullk Richardson number computed for each cell in a column,
995  ! assuming OBLdepth = grid cell depth. After Rib(k) is
996  ! known for the column, then CVMix interpolates to find
997  ! the actual OBLdepth. This approach avoids need to iterate
998  ! on the OBLdepth calculation. It follows that used in MOM5
999  ! and POP.
1000  ifaceheight(1) = 0.0 ! BBL is all relative to the surface
1001  pref = 0.
1002  hcorr = 0.
1003  do k=1,g%ke
1004 
1005  ! cell center and cell bottom in meters (negative values in the ocean)
1006  dh = h(i,j,k) * gv%H_to_m ! Nominal thickness to use for increment
1007  dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0)
1008  hcorr = min( dh - cs%min_thickness, 0. ) ! If inflating then hcorr<0
1009  dh = max( dh, cs%min_thickness ) ! Limit increment dh>=min_thickness
1010  cellheight(k) = ifaceheight(k) - 0.5 * dh
1011  ifaceheight(k+1) = ifaceheight(k) - dh
1012 
1013  ! find ksfc for cell where "surface layer" sits
1014  sldepth_0d = cs%surf_layer_ext*max( max(-cellheight(k),-ifaceheight(2) ), cs%minOBLdepth )
1015  ksfc = k
1016  do ktmp = 1,k
1017  if (-1.0*ifaceheight(ktmp+1) >= sldepth_0d) then
1018  ksfc = ktmp
1019  exit
1020  endif
1021  enddo
1022 
1023  ! average temp, saln, u, v over surface layer
1024  ! use C-grid average to get u,v on T-points.
1025  surfhtemp=0.0
1026  surfhsalt=0.0
1027  surfhu =0.0
1028  surfhv =0.0
1029  surfhus =0.0
1030  surfhvs =0.0
1031  htot =0.0
1032  do ktmp = 1,ksfc
1033 
1034  ! SLdepth_0d can be between cell interfaces
1035  delh = min( max(0.0, sldepth_0d - htot), h(i,j,ktmp)*gv%H_to_m )
1036 
1037  ! surface layer thickness
1038  htot = htot + delh
1039 
1040  ! surface averaged fields
1041  surfhtemp = surfhtemp + temp(i,j,ktmp) * delh
1042  surfhsalt = surfhsalt + salt(i,j,ktmp) * delh
1043  surfhu = surfhu + 0.5*us%L_T_to_m_s*(u(i,j,ktmp)+u(i-1,j,ktmp)) * delh
1044  surfhv = surfhv + 0.5*us%L_T_to_m_s*(v(i,j,ktmp)+v(i,j-1,ktmp)) * delh
1045  if (cs%Stokes_Mixing) then
1046  surfhus = surfhus + 0.5*(waves%US_x(i,j,ktmp)+waves%US_x(i-1,j,ktmp)) * delh
1047  surfhvs = surfhvs + 0.5*(waves%US_y(i,j,ktmp)+waves%US_y(i,j-1,ktmp)) * delh
1048  endif
1049 
1050  enddo
1051  surftemp = surfhtemp / htot
1052  surfsalt = surfhsalt / htot
1053  surfu = surfhu / htot
1054  surfv = surfhv / htot
1055  surfus = surfhus / htot
1056  surfvs = surfhvs / htot
1057 
1058  ! vertical shear between present layer and
1059  ! surface layer averaged surfU,surfV.
1060  ! C-grid average to get Uk and Vk on T-points.
1061  uk = 0.5*us%L_T_to_m_s*(u(i,j,k)+u(i-1,j,k)) - surfu
1062  vk = 0.5*us%L_T_to_m_s*(v(i,j,k)+v(i,j-1,k)) - surfv
1063 
1064  if (cs%Stokes_Mixing) then
1065  ! If momentum is mixed down the Stokes drift gradient, then
1066  ! the Stokes drift must be included in the bulk Richardson number
1067  ! calculation.
1068  uk = uk + (0.5*(waves%Us_x(i,j,k)+waves%US_x(i-1,j,k)) -surfus )
1069  vk = vk + (0.5*(waves%Us_y(i,j,k)+waves%Us_y(i,j-1,k)) -surfvs )
1070  endif
1071 
1072  deltau2(k) = uk**2 + vk**2
1073 
1074  ! pressure, temp, and saln for EOS
1075  ! kk+1 = surface fields
1076  ! kk+2 = k fields
1077  ! kk+3 = km1 fields
1078  km1 = max(1, k-1)
1079  kk = 3*(k-1)
1080  pres_1d(kk+1) = pref
1081  pres_1d(kk+2) = pref
1082  pres_1d(kk+3) = pref
1083  temp_1d(kk+1) = surftemp
1084  temp_1d(kk+2) = temp(i,j,k)
1085  temp_1d(kk+3) = temp(i,j,km1)
1086  salt_1d(kk+1) = surfsalt
1087  salt_1d(kk+2) = salt(i,j,k)
1088  salt_1d(kk+3) = salt(i,j,km1)
1089 
1090  ! pRef is pressure at interface between k and km1.
1091  ! iterate pRef for next pass through k-loop.
1092  pref = pref + gv%H_to_Pa * h(i,j,k)
1093 
1094  ! this difference accounts for penetrating SW
1095  surfbuoyflux2(k) = buoy_scale * (buoyflux(i,j,1) - buoyflux(i,j,k+1))
1096 
1097  enddo ! k-loop finishes
1098 
1099  if (cs%LT_K_ENHANCEMENT .or. cs%LT_VT2_ENHANCEMENT) then
1100  mld_guess = max( 1.*us%m_to_Z, abs(us%m_to_Z*cs%OBLdepthprev(i,j) ) )
1101  call get_langmuir_number(la, g, gv, us, mld_guess, ustar(i,j), i, j, &
1102  h=h(i,j,:), u_h=u_h, v_h=v_h, waves=waves)
1103  cs%La_SL(i,j)=la
1104  endif
1105 
1106 
1107  ! compute in-situ density
1108  call calculate_density(temp_1d, salt_1d, pres_1d, rho_1d, 1, 3*g%ke, eos)
1109 
1110  ! N2 (can be negative) and N (non-negative) on interfaces.
1111  ! deltaRho is non-local rho difference used for bulk Richardson number.
1112  ! CS%N is local N (with floor) used for unresolved shear calculation.
1113  do k = 1, g%ke
1114  km1 = max(1, k-1)
1115  kk = 3*(k-1)
1116  deltarho(k) = rho_1d(kk+2) - rho_1d(kk+1)
1117  n2_1d(k) = (gorho * (rho_1d(kk+2) - rho_1d(kk+3)) ) / &
1118  ((0.5*(h(i,j,km1) + h(i,j,k))+gv%H_subroundoff)*gv%H_to_m)
1119  cs%N(i,j,k) = sqrt( max( n2_1d(k), 0.) )
1120  enddo
1121  n2_1d(g%ke+1 ) = 0.0
1122  cs%N(i,j,g%ke+1 ) = 0.0
1123 
1124  ! turbulent velocity scales w_s and w_m computed at the cell centers.
1125  ! Note that if sigma > CS%surf_layer_ext, then CVMix_kpp_compute_turbulent_scales
1126  ! computes w_s and w_m velocity scale at sigma=CS%surf_layer_ext. So we only pass
1127  ! sigma=CS%surf_layer_ext for this calculation.
1128  call cvmix_kpp_compute_turbulent_scales( &
1129  cs%surf_layer_ext, & ! (in) Normalized surface layer depth; sigma = CS%surf_layer_ext
1130  -cellheight, & ! (in) Assume here that OBL depth [m] = -cellHeight(k)
1131  surfbuoyflux2, & ! (in) Buoyancy flux at surface [m2 s-3]
1132  surffricvel, & ! (in) Turbulent friction velocity at surface [m s-1]
1133  w_s=ws_1d, & ! (out) Turbulent velocity scale profile [m s-1]
1134  cvmix_kpp_params_user=cs%KPP_params )
1135 
1136  !Compute CVMix VT2
1137  cs%Vt2(i,j,:) = cvmix_kpp_compute_unresolved_shear( &
1138  zt_cntr=cellheight(1:g%ke), & ! Depth of cell center [m]
1139  ws_cntr=ws_1d, & ! Turbulent velocity scale profile, at centers [m s-1]
1140  n_iface=cs%N(i,j,:), & ! Buoyancy frequency at interface [s-1]
1141  cvmix_kpp_params_user=cs%KPP_params ) ! KPP parameters
1142 
1143  !Modify CVMix VT2
1144  IF (cs%LT_VT2_ENHANCEMENT) then
1145  IF (cs%LT_VT2_METHOD==lt_vt2_mode_constant) then
1146  do k=1,g%ke
1147  langenhvt2(k) = cs%KPP_VT2_ENH_FAC
1148  enddo
1149  elseif (cs%LT_VT2_METHOD==lt_vt2_mode_vr12) then
1150  !Introduced minimum value for La_SL, so maximum value for enhvt2 is removed.
1151  enhvt2 = sqrt(1.+(1.5*cs%La_SL(i,j))**(-2) + &
1152  (5.4*cs%La_SL(i,j))**(-4))
1153  do k=1,g%ke
1154  langenhvt2(k) = enhvt2
1155  enddo
1156  elseif (cs%LT_VT2_METHOD==lt_vt2_mode_rw16) then
1157  !Introduced minimum value for La_SL, so maximum value for enhvt2 is removed.
1158  enhvt2 = 1. + 2.3*cs%La_SL(i,j)**(-0.5)
1159  do k=1,g%ke
1160  langenhvt2(k) = enhvt2
1161  enddo
1162  elseif (cs%LT_VT2_METHOD==lt_vt2_mode_lf17) then
1163  cs%CS=cvmix_get_kpp_real('c_s',cs%KPP_params)
1164  do k=1,g%ke
1165  wst = (max(0.,-buoy_scale*buoyflux(i,j,1))*(-cellheight(k)))**(1./3.)
1166  langenhvt2(k) = sqrt((0.15*wst**3. + 0.17*surffricvel**3.* &
1167  (1.+0.49*cs%La_SL(i,j)**(-2.))) / &
1168  (0.2*ws_1d(k)**3/(cs%cs*cs%surf_layer_ext*cs%vonKarman**4.)))
1169  enddo
1170  else
1171  !This shouldn't be reached.
1172  !call MOM_error(WARNING,"Unexpected behavior in MOM_CVMix_KPP, see error in Vt2")
1173  langenhvt2(:) = 1.0
1174  endif
1175  else
1176  langenhvt2(:) = 1.0
1177  endif
1178 
1179  do k=1,g%ke
1180  cs%Vt2(i,j,k)=cs%Vt2(i,j,k)*langenhvt2(k)
1181  if (cs%id_EnhVt2 > 0) cs%EnhVt2(i,j,k)=langenhvt2(k)
1182  enddo
1183 
1184  ! Calculate Bulk Richardson number from eq (21) of LMD94
1185  bulkri_1d = cvmix_kpp_compute_bulk_richardson( &
1186  zt_cntr = cellheight(1:g%ke), & ! Depth of cell center [m]
1187  delta_buoy_cntr=gorho*deltarho, & ! Bulk buoyancy difference, Br-B(z) [s-1]
1188  delta_vsqr_cntr=deltau2, & ! Square of resolved velocity difference [m2 s-2]
1189  vt_sqr_cntr=cs%Vt2(i,j,:), &
1190  ws_cntr=ws_1d, & ! Turbulent velocity scale profile [m s-1]
1191  n_iface=cs%N(i,j,:)) ! Buoyancy frequency [s-1]
1192 
1193 
1194  surfbuoyflux = buoy_scale * buoyflux(i,j,1) ! This is only used in kpp_compute_OBL_depth to limit
1195  ! h to Monin-Obukov (default is false, ie. not used)
1196 
1197  call cvmix_kpp_compute_obl_depth( &
1198  bulkri_1d, & ! (in) Bulk Richardson number
1199  ifaceheight, & ! (in) Height of interfaces [m]
1200  cs%OBLdepth(i,j), & ! (out) OBL depth [m]
1201  cs%kOBL(i,j), & ! (out) level (+fraction) of OBL extent
1202  zt_cntr=cellheight, & ! (in) Height of cell centers [m]
1203  surf_fric=surffricvel, & ! (in) Turbulent friction velocity at surface [m s-1]
1204  surf_buoy=surfbuoyflux, & ! (in) Buoyancy flux at surface [m2 s-3]
1205  coriolis=coriolis, & ! (in) Coriolis parameter [s-1]
1206  cvmix_kpp_params_user=cs%KPP_params ) ! KPP parameters
1207 
1208  ! A hack to avoid KPP reaching the bottom. It was needed during development
1209  ! because KPP was unable to handle vanishingly small layers near the bottom.
1210  if (cs%deepOBLoffset>0.) then
1211  zbottomminusoffset = ifaceheight(g%ke+1) + min(cs%deepOBLoffset,-0.1*ifaceheight(g%ke+1))
1212  cs%OBLdepth(i,j) = min( cs%OBLdepth(i,j), -zbottomminusoffset )
1213  endif
1214 
1215  ! apply some constraints on OBLdepth
1216  if(cs%fixedOBLdepth) cs%OBLdepth(i,j) = cs%fixedOBLdepth_value
1217  cs%OBLdepth(i,j) = max( cs%OBLdepth(i,j), -ifaceheight(2) ) ! no shallower than top layer
1218  cs%OBLdepth(i,j) = min( cs%OBLdepth(i,j), -ifaceheight(g%ke+1) ) ! no deeper than bottom
1219  cs%kOBL(i,j) = cvmix_kpp_compute_kobl_depth( ifaceheight, cellheight, cs%OBLdepth(i,j) )
1220 
1221 !*************************************************************************
1222 ! smg: remove code below
1223 
1224 ! Following "correction" step has been found to be unnecessary.
1225 ! Code should be removed after further testing.
1226 ! BGR: 03/15/2018-> Restructured code (Vt2 changed to compute from call in MOM_CVMix_KPP now)
1227 ! I have not taken this restructuring into account here.
1228 ! Do we ever run with correctSurfLayerAvg?
1229 ! smg's suggested testing and removal is advised, in the meantime
1230 ! I have added warning if correctSurfLayerAvg is attempted.
1231  ! if (CS%correctSurfLayerAvg) then
1232 
1233  ! SLdepth_0d = CS%surf_layer_ext * CS%OBLdepth(i,j)
1234  ! hTot = h(i,j,1)
1235  ! surfTemp = Temp(i,j,1) ; surfHtemp = surfTemp * hTot
1236  ! surfSalt = Salt(i,j,1) ; surfHsalt = surfSalt * hTot
1237  ! surfU = 0.5*US%L_T_to_m_s*(u(i,j,1)+u(i-1,j,1)) ; surfHu = surfU * hTot
1238  ! surfV = 0.5*US%L_T_to_m_s*(v(i,j,1)+v(i,j-1,1)) ; surfHv = surfV * hTot
1239  ! pRef = 0.0
1240 
1241  ! do k = 2, G%ke
1242 
1243  ! ! Recalculate differences with surface layer
1244  ! Uk = 0.5*US%L_T_to_m_s*(u(i,j,k)+u(i-1,j,k)) - surfU
1245  ! Vk = 0.5*US%L_T_to_m_s*(v(i,j,k)+v(i,j-1,k)) - surfV
1246  ! deltaU2(k) = Uk**2 + Vk**2
1247  ! pRef = pRef + GV%H_to_Pa * h(i,j,k)
1248  ! call calculate_density(surfTemp, surfSalt, pRef, rho1, EOS)
1249  ! call calculate_density(Temp(i,j,k), Salt(i,j,k), pRef, rhoK, EOS)
1250  ! deltaRho(k) = rhoK - rho1
1251 
1252  ! ! Surface layer averaging (needed for next k+1 iteration of this loop)
1253  ! if (hTot < SLdepth_0d) then
1254  ! delH = min( max(0., SLdepth_0d - hTot), h(i,j,k)*GV%H_to_m )
1255  ! hTot = hTot + delH
1256  ! surfHtemp = surfHtemp + Temp(i,j,k) * delH ; surfTemp = surfHtemp / hTot
1257  ! surfHsalt = surfHsalt + Salt(i,j,k) * delH ; surfSalt = surfHsalt / hTot
1258  ! surfHu = surfHu + 0.5*US%L_T_to_m_s*(u(i,j,k)+u(i-1,j,k)) * delH ; surfU = surfHu / hTot
1259  ! surfHv = surfHv + 0.5*US%L_T_to_m_s*(v(i,j,k)+v(i,j-1,k)) * delH ; surfV = surfHv / hTot
1260  ! endif
1261 
1262  ! enddo
1263 
1264  ! BulkRi_1d = CVMix_kpp_compute_bulk_Richardson( &
1265  ! cellHeight(1:G%ke), & ! Depth of cell center [m]
1266  ! GoRho*deltaRho, & ! Bulk buoyancy difference, Br-B(z) [s-1]
1267  ! deltaU2, & ! Square of resolved velocity difference [m2 s-2]
1268  ! ws_cntr=Ws_1d, & ! Turbulent velocity scale profile [m s-1]
1269  ! N_iface=CS%N ) ! Buoyancy frequency [s-1]
1270 
1271  ! surfBuoyFlux = buoy_scale*buoyFlux(i,j,1) ! This is only used in kpp_compute_OBL_depth to limit
1272  ! ! h to Monin-Obukov (default is false, ie. not used)
1273 
1274  ! call CVMix_kpp_compute_OBL_depth( &
1275  ! BulkRi_1d, & ! (in) Bulk Richardson number
1276  ! iFaceHeight, & ! (in) Height of interfaces [m]
1277  ! CS%OBLdepth(i,j), & ! (out) OBL depth [m]
1278  ! CS%kOBL(i,j), & ! (out) level (+fraction) of OBL extent
1279  ! zt_cntr=cellHeight, & ! (in) Height of cell centers [m]
1280  ! surf_fric=surfFricVel, & ! (in) Turbulent friction velocity at surface [m s-1]
1281  ! surf_buoy=surfBuoyFlux, & ! (in) Buoyancy flux at surface [m2 s-3]
1282  ! Coriolis=Coriolis, & ! (in) Coriolis parameter [s-1]
1283  ! CVMix_kpp_params_user=CS%KPP_params ) ! KPP parameters
1284 
1285  ! if (CS%deepOBLoffset>0.) then
1286  ! zBottomMinusOffset = iFaceHeight(G%ke+1) + min(CS%deepOBLoffset,-0.1*iFaceHeight(G%ke+1))
1287  ! CS%OBLdepth(i,j) = min( CS%OBLdepth(i,j), -zBottomMinusOffset )
1288  ! CS%kOBL(i,j) = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, CS%OBLdepth(i,j) )
1289  ! endif
1290 
1291  ! ! apply some constraints on OBLdepth
1292  ! if(CS%fixedOBLdepth) CS%OBLdepth(i,j) = CS%fixedOBLdepth_value
1293  ! CS%OBLdepth(i,j) = max( CS%OBLdepth(i,j), -iFaceHeight(2) ) ! no shallower than top layer
1294  ! CS%OBLdepth(i,j) = min( CS%OBLdepth(i,j), -iFaceHeight(G%ke+1) ) ! no deep than bottom
1295  ! CS%kOBL(i,j) = CVMix_kpp_compute_kOBL_depth( iFaceHeight, cellHeight, CS%OBLdepth(i,j) )
1296 
1297  ! endif ! endif for "correction" step
1298 
1299 ! smg: remove code above
1300 ! **********************************************************************
1301 
1302  ! recompute wscale for diagnostics, now that we in fact know boundary layer depth
1303  !BGR consider if LTEnhancement is wanted for diagnostics
1304  if (cs%id_Ws > 0) then
1305  call cvmix_kpp_compute_turbulent_scales( &
1306  -cellheight/cs%OBLdepth(i,j), & ! (in) Normalized boundary layer coordinate
1307  cs%OBLdepth(i,j), & ! (in) OBL depth [m]
1308  surfbuoyflux, & ! (in) Buoyancy flux at surface [m2 s-3]
1309  surffricvel, & ! (in) Turbulent friction velocity at surface [m s-1]
1310  w_s=ws_1d, & ! (out) Turbulent velocity scale profile [m s-1]
1311  cvmix_kpp_params_user=cs%KPP_params) ! KPP parameters
1312  cs%Ws(i,j,:) = ws_1d(:)
1313  endif
1314 
1315  ! Diagnostics
1316  if (cs%id_N2 > 0) cs%N2(i,j,:) = n2_1d(:)
1317  if (cs%id_BulkDrho > 0) cs%dRho(i,j,:) = deltarho(:)
1318  if (cs%id_BulkRi > 0) cs%BulkRi(i,j,:) = bulkri_1d(:)
1319  if (cs%id_BulkUz2 > 0) cs%Uz2(i,j,:) = deltau2(:)
1320  if (cs%id_Tsurf > 0) cs%Tsurf(i,j) = surftemp
1321  if (cs%id_Ssurf > 0) cs%Ssurf(i,j) = surfsalt
1322  if (cs%id_Usurf > 0) cs%Usurf(i,j) = surfu
1323  if (cs%id_Vsurf > 0) cs%Vsurf(i,j) = surfv
1324 
1325  enddo
1326  enddo
1327 
1328  ! send diagnostics to post_data
1329  if (cs%id_BulkRi > 0) call post_data(cs%id_BulkRi, cs%BulkRi, cs%diag)
1330  if (cs%id_N > 0) call post_data(cs%id_N, cs%N, cs%diag)
1331  if (cs%id_N2 > 0) call post_data(cs%id_N2, cs%N2, cs%diag)
1332  if (cs%id_Tsurf > 0) call post_data(cs%id_Tsurf, cs%Tsurf, cs%diag)
1333  if (cs%id_Ssurf > 0) call post_data(cs%id_Ssurf, cs%Ssurf, cs%diag)
1334  if (cs%id_Usurf > 0) call post_data(cs%id_Usurf, cs%Usurf, cs%diag)
1335  if (cs%id_Vsurf > 0) call post_data(cs%id_Vsurf, cs%Vsurf, cs%diag)
1336  if (cs%id_BulkDrho > 0) call post_data(cs%id_BulkDrho, cs%dRho, cs%diag)
1337  if (cs%id_BulkUz2 > 0) call post_data(cs%id_BulkUz2, cs%Uz2, cs%diag)
1338  if (cs%id_EnhK > 0) call post_data(cs%id_EnhK, cs%EnhK, cs%diag)
1339  if (cs%id_EnhVt2 > 0) call post_data(cs%id_EnhVt2, cs%EnhVt2, cs%diag)
1340  if (cs%id_La_SL > 0) call post_data(cs%id_La_SL, cs%La_SL, cs%diag)
1341 
1342  ! BLD smoothing:
1343  if (cs%n_smooth > 0) call kpp_smooth_bld(cs,g,gv,h)
1344 
1345 end subroutine kpp_compute_bld
1346 
1347 
1348 !> Apply a 1-1-4-1-1 Laplacian filter one time on BLD to reduce any horizontal two-grid-point noise
1349 subroutine kpp_smooth_bld(CS,G,GV,h)
1350  ! Arguments
1351  type(kpp_cs), pointer :: CS !< Control structure
1352  type(ocean_grid_type), intent(inout) :: G !< Ocean grid
1353  type(verticalgrid_type), intent(in) :: GV !< Ocean vertical grid
1354  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer/level thicknesses [H ~> m or kg m-2]
1355 
1356  ! local
1357  real, dimension( G%ke ) :: cellHeight ! Cell center heights referenced to surface [m]
1358  ! (negative in the ocean)
1359  real, dimension( G%ke+1 ) :: iFaceHeight ! Interface heights referenced to surface [m]
1360  ! (negative in the ocean)
1361  real :: wc, ww, we, wn, ws ! averaging weights for smoothing
1362  real :: dh ! The local thickness used for calculating interface positions [m]
1363  real :: hcorr ! A cumulative correction arising from inflation of vanished layers [m]
1364  real :: pref
1365  integer :: i, j, k, s
1366 
1367  ncall_smooth = ncall_smooth+1
1368 
1369  if (cs%id_OBLdepth_original > 0) cs%OBLdepth_original = cs%OBLdepth
1370 
1371  cs%n_smooth = 2
1372  do s=1,cs%n_smooth
1373 
1374  ! Update halos
1375  !call pass_var(CS%OBLdepth, G%Domain, halo=4)
1376  if (s==1 .or. ncall_smooth<3) then
1377  call pass_var(cs%OBLdepth, g%Domain, halo=4)
1378  endif
1379 
1380  cs%OBLdepthprev = cs%OBLdepth
1381 
1382  ! apply smoothing on OBL depth
1383  do j = g%jsc-2, g%jec+2
1384  do i = g%isc-2, g%iec+2
1385 
1386  ! skip land points
1387  if (g%mask2dT(i,j)==0.) cycle
1388 
1389  ! compute weights
1390  !! fail
1391  ww = 0.125 * g%mask2dT(i-1,j)
1392  we = 0.125 * g%mask2dT(i+1,j)
1393  ws = 0.125 * g%mask2dT(i,j-1)
1394  wn = 0.125 * g%mask2dT(i,j+1)
1395 
1396  !! pass
1397  !ww = 0.25 * G%mask2dT(i-1,j) ! 0.125 * G%mask2dT(i-1,j)
1398  !we = 0.25 * G%mask2dT(i-1,j) ! 0.125 * G%mask2dT(i+1,j)
1399  !ws = 0.0 ! 0.125 * G%mask2dT(i,j-1)
1400  !wn = 0.0 ! 0.125 * G%mask2dT(i,j+1)
1401 
1402  !! fail
1403  !!ww = 0.0 * G%mask2dT(i-1,j) ! 0.125 * G%mask2dT(i-1,j)
1404  !!we = 0.0 * G%mask2dT(i-1,j) ! 0.125 * G%mask2dT(i+1,j)
1405  !!ws = 0.25 ! 0.125 * G%mask2dT(i,j-1)
1406  !!wn = 0.25 ! 0.125 * G%mask2dT(i,j+1)
1407 
1408  !! pass
1409  !!ww = 0.0 * G%mask2dT(i-1,j) ! 0.125 * G%mask2dT(i-1,j)
1410  !!we = 0.0 * G%mask2dT(i-1,j) ! 0.125 * G%mask2dT(i+1,j)
1411  !!ws = 0.25 ! 0.125 * G%mask2dT(i,j-1)
1412  !!wn = 0.0 ! 0.125 * G%mask2dT(i,j+1)
1413 
1414  !! fail
1415  !!ww = 0.0 * G%mask2dT(i-1,j) ! 0.125 * G%mask2dT(i-1,j)
1416  !!we = 0.0 * G%mask2dT(i-1,j) ! 0.125 * G%mask2dT(i+1,j)
1417  !!ws = 0.0 ! 0.125 * G%mask2dT(i,j-1)
1418  !!wn = 0.25 ! 0.125 * G%mask2dT(i,j+1)
1419 
1420  wc = 1.0 - (ww+we+wn+ws)
1421 
1422  cs%OBLdepth(i,j) = wc * cs%OBLdepthprev(i,j) &
1423  + ww * cs%OBLdepthprev(i-1,j) &
1424  + we * cs%OBLdepthprev(i+1,j) &
1425  + ws * cs%OBLdepthprev(i,j-1) &
1426  + wn * cs%OBLdepthprev(i,j+1)
1427 
1428  enddo
1429  enddo
1430 
1431  enddo ! s-loop
1432 
1433  ! Update kOBL for smoothed OBL depths
1434  do j = g%jsc, g%jec
1435  do i = g%isc, g%iec
1436 
1437  ! skip land points
1438  if (g%mask2dT(i,j)==0.) cycle
1439 
1440  ifaceheight(1) = 0.0 ! BBL is all relative to the surface
1441  hcorr = 0.
1442  do k=1,g%ke
1443 
1444  ! cell center and cell bottom in meters (negative values in the ocean)
1445  dh = h(i,j,k) * gv%H_to_m ! Nominal thickness to use for increment
1446  dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0)
1447  hcorr = min( dh - cs%min_thickness, 0. ) ! If inflating then hcorr<0
1448  dh = max( dh, cs%min_thickness ) ! Limit increment dh>=min_thickness
1449  cellheight(k) = ifaceheight(k) - 0.5 * dh
1450  ifaceheight(k+1) = ifaceheight(k) - dh
1451  enddo
1452 
1453  ! Apply OBLdepth smoothing at a cell only if the OBLdepth gets deeper via smoothing.
1454  if (cs%deepen_only) cs%OBLdepth(i,j) = max(cs%OBLdepth(i,j),cs%OBLdepth_original(i,j))
1455 
1456  ! prevent OBL depths deeper than the bathymetric depth
1457  cs%OBLdepth(i,j) = min( cs%OBLdepth(i,j), -ifaceheight(g%ke+1) ) ! no deeper than bottom
1458 
1459  cs%kOBL(i,j) = cvmix_kpp_compute_kobl_depth( ifaceheight, cellheight, cs%OBLdepth(i,j) )
1460  enddo
1461  enddo
1462 
1463 end subroutine kpp_smooth_bld
1464 
1465 
1466 
1467 !> Copies KPP surface boundary layer depth into BLD
1468 subroutine kpp_get_bld(CS, BLD, G)
1469  type(kpp_cs), pointer :: cs !< Control structure for
1470  !! this module
1471  type(ocean_grid_type), intent(in) :: g !< Grid structure
1472  real, dimension(SZI_(G),SZJ_(G)), intent(inout) :: bld!< bnd. layer depth [m]
1473  ! Local variables
1474  integer :: i,j
1475  do j = g%jsc, g%jec ; do i = g%isc, g%iec
1476  bld(i,j) = cs%OBLdepth(i,j)
1477  enddo ; enddo
1478 end subroutine kpp_get_bld
1479 
1480 !> Apply KPP non-local transport of surface fluxes for temperature.
1481 subroutine kpp_nonlocaltransport_temp(CS, G, GV, h, nonLocalTrans, surfFlux, &
1482  dt, scalar, C_p)
1483 
1484  type(kpp_cs), intent(in) :: cs !< Control structure
1485  type(ocean_grid_type), intent(in) :: g !< Ocean grid
1486  type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid
1487  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer/level thickness [H ~> m or kg m-2]
1488  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in) :: nonlocaltrans !< Non-local transport [nondim]
1489  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: surfflux !< Surface flux of scalar
1490  !! [conc H s-1 ~> conc m s-1 or conc kg m-2 s-1]
1491  real, intent(in) :: dt !< Time-step [s]
1492  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: scalar !< temperature
1493  real, intent(in) :: c_p !< Seawater specific heat capacity [J kg-1 degC-1]
1494 
1495  integer :: i, j, k
1496  real, dimension( SZI_(G), SZJ_(G), SZK_(G) ) :: dtracer
1497 
1498 
1499  dtracer(:,:,:) = 0.0
1500  !$OMP parallel do default(none) shared(dtracer, nonLocalTrans, h, G, GV, surfFlux)
1501  do k = 1, g%ke
1502  do j = g%jsc, g%jec
1503  do i = g%isc, g%iec
1504  dtracer(i,j,k) = ( nonlocaltrans(i,j,k) - nonlocaltrans(i,j,k+1) ) / &
1505  ( h(i,j,k) + gv%H_subroundoff ) * surfflux(i,j)
1506  enddo
1507  enddo
1508  enddo
1509 
1510  ! Update tracer due to non-local redistribution of surface flux
1511  if (cs%applyNonLocalTrans) then
1512  do k = 1, g%ke
1513  do j = g%jsc, g%jec
1514  do i = g%isc, g%iec
1515  scalar(i,j,k) = scalar(i,j,k) + dt * dtracer(i,j,k)
1516  enddo
1517  enddo
1518  enddo
1519  endif
1520 
1521  ! Diagnostics
1522  if (cs%id_QminusSW > 0) call post_data(cs%id_QminusSW, surfflux, cs%diag)
1523  if (cs%id_NLT_dTdt > 0) call post_data(cs%id_NLT_dTdt, dtracer, cs%diag)
1524  if (cs%id_NLT_temp_budget > 0) then
1525  dtracer(:,:,:) = 0.0
1526  do k = 1, g%ke
1527  do j = g%jsc, g%jec
1528  do i = g%isc, g%iec
1529  dtracer(i,j,k) = (nonlocaltrans(i,j,k) - nonlocaltrans(i,j,k+1)) * &
1530  surfflux(i,j) * c_p * gv%H_to_kg_m2
1531  enddo
1532  enddo
1533  enddo
1534  call post_data(cs%id_NLT_temp_budget, dtracer, cs%diag)
1535  endif
1536 
1537 end subroutine kpp_nonlocaltransport_temp
1538 
1539 
1540 !> Apply KPP non-local transport of surface fluxes for salinity.
1541 !> This routine is a useful prototype for other material tracers.
1542 subroutine kpp_nonlocaltransport_saln(CS, G, GV, h, nonLocalTrans, surfFlux, dt, scalar)
1543 
1544  type(kpp_cs), intent(in) :: cs !< Control structure
1545  type(ocean_grid_type), intent(in) :: g !< Ocean grid
1546  type(verticalgrid_type), intent(in) :: gv !< Ocean vertical grid
1547  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h !< Layer/level thickness [H ~> m or kg m-2]
1548  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in) :: nonlocaltrans !< Non-local transport [nondim]
1549  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: surfflux !< Surface flux of scalar
1550  !! [conc H s-1 ~> conc m s-1 or conc kg m-2 s-1]
1551  real, intent(in) :: dt !< Time-step [s]
1552  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(inout) :: scalar !< Scalar (scalar units [conc])
1553 
1554  integer :: i, j, k
1555  real, dimension( SZI_(G), SZJ_(G), SZK_(G) ) :: dtracer
1556 
1557 
1558  dtracer(:,:,:) = 0.0
1559  !$OMP parallel do default(none) shared(dtracer, nonLocalTrans, h, G, GV, surfFlux)
1560  do k = 1, g%ke
1561  do j = g%jsc, g%jec
1562  do i = g%isc, g%iec
1563  dtracer(i,j,k) = ( nonlocaltrans(i,j,k) - nonlocaltrans(i,j,k+1) ) / &
1564  ( h(i,j,k) + gv%H_subroundoff ) * surfflux(i,j)
1565  enddo
1566  enddo
1567  enddo
1568 
1569  ! Update tracer due to non-local redistribution of surface flux
1570  if (cs%applyNonLocalTrans) then
1571  do k = 1, g%ke
1572  do j = g%jsc, g%jec
1573  do i = g%isc, g%iec
1574  scalar(i,j,k) = scalar(i,j,k) + dt * dtracer(i,j,k)
1575  enddo
1576  enddo
1577  enddo
1578  endif
1579 
1580  ! Diagnostics
1581  if (cs%id_netS > 0) call post_data(cs%id_netS, surfflux, cs%diag)
1582  if (cs%id_NLT_dSdt > 0) call post_data(cs%id_NLT_dSdt, dtracer, cs%diag)
1583  if (cs%id_NLT_saln_budget > 0) then
1584  dtracer(:,:,:) = 0.0
1585  do k = 1, g%ke
1586  do j = g%jsc, g%jec
1587  do i = g%isc, g%iec
1588  dtracer(i,j,k) = (nonlocaltrans(i,j,k) - nonlocaltrans(i,j,k+1)) * &
1589  surfflux(i,j) * gv%H_to_kg_m2
1590  enddo
1591  enddo
1592  enddo
1593  call post_data(cs%id_NLT_saln_budget, dtracer, cs%diag)
1594  endif
1595 
1596 end subroutine kpp_nonlocaltransport_saln
1597 
1598 
1599 
1600 
1601 !> Clear pointers, deallocate memory
1602 subroutine kpp_end(CS)
1603  type(kpp_cs), pointer :: cs !< Control structure
1604 
1605  if (.not.associated(cs)) return
1606 
1607  deallocate(cs)
1608 
1609 end subroutine kpp_end
1610 
1611 !> \namespace mom_cvmix_kpp
1612 !!
1613 !! \section section_KPP The K-Profile Parameterization
1614 !!
1615 !! The K-Profile Parameterization (KPP) of Large et al., 1994, (http://dx.doi.org/10.1029/94RG01872) is
1616 !! implemented via the Community Vertical Mixing package, [CVMix](http://cvmix.github.io/),
1617 !! which is called directly by this module.
1618 !!
1619 !! The formulation and implementation of KPP is described in great detail in the
1620 !! [CVMix manual](https://github.com/CVMix/CVMix-description/raw/master/cvmix.pdf) (written by our own Steve Griffies).
1621 !!
1622 !! \subsection section_KPP_nutshell KPP in a nutshell
1623 !!
1624 !! Large et al., 1994, decompose the parameterized boundary layer turbulent flux of a scalar, \f$ s \f$, as
1625 !! \f[ \overline{w^\prime s^\prime} = -K \partial_z s + K \gamma_s(\sigma), \f]
1626 !! where \f$ \sigma = -z/h \f$ is a non-dimensional coordinate within the boundary layer of depth \f$ h \f$.
1627 !! \f$ K \f$ is the eddy diffusivity and is a function of position within the boundary layer as well as a
1628 !! function of the surface forcing:
1629 !! \f[ K = h w_s(\sigma) G(\sigma) . \f]
1630 !! Here, \f$ w_s \f$ is the vertical velocity scale of the boundary layer turbulence and \f$ G(\sigma) \f$ is
1631 !! a "shape function" which is described later.
1632 !! The last term is the "non-local transport" which involves a function \f$ \gamma_s(\sigma) \f$ that is matched
1633 !! to the forcing but is not actually needed in the final implementation.
1634 !! Instead, the entire non-local transport term can be equivalently written
1635 !! \f[ K \gamma_s(\sigma) = C_s G(\sigma) Q_s \f]
1636 !! where \f$ Q_s \f$ is the surface flux of \f$ s \f$ and \f$ C_s \f$ is a constant.
1637 !! The vertical structure of the redistribution (non-local) term is solely due to the shape function,
1638 !! \f$ G(\sigma) \f$.
1639 !! In our implementation of KPP, we allow the shape functions used for \f$ K \f$ and for the non-local transport
1640 !! to be chosen independently.
1641 !!
1642 !! [google_thread_NLT]: https://groups.google.com/forum/#!msg/CVMix-dev/i6rF-eHOtKI/Ti8BeyksrhAJ
1643 !! "Extreme values of non-local transport"
1644 !!
1645 !! The particular shape function most widely used in the atmospheric community is
1646 !! \f[ G(\sigma) = \sigma (1-\sigma)^2 \f]
1647 !! which satisfies the boundary conditions
1648 !! \f$ G(0) = 0 \f$,
1649 !! \f$ G(1) = 0 \f$,
1650 !! \f$ G^\prime(0) = 1 \f$, and
1651 !! \f$ G^\prime(1) = 0 \f$.
1652 !! Large et al, 1994, alter the function so as to match interior diffusivities but we have found that this leads
1653 !! to inconsistencies within the formulation (see google groups thread
1654 !! [Extreme values of non-local transport][google_thread_NLT]).
1655 !! Instead, we use either the above form, or even simpler forms that use alternative upper boundary conditions.
1656 !!
1657 !! The KPP boundary layer depth is a function of the bulk Richardson number, Rib.
1658 !! But to compute Rib, we need the boundary layer depth. To address this circular
1659 !! logic, we compute Rib for each vertical cell in a column, assuming the BL depth
1660 !! equals to the depth of the given grid cell. Once we have a vertical array of Rib(k),
1661 !! we then call the OBLdepth routine from CVMix to compute the actual
1662 !! OBLdepth. We optionally then "correct" the OBLdepth by cycling through once more,
1663 !! this time knowing the OBLdepth from the first pass. This "correction" step is not
1664 !! used by NCAR. It has been found in idealized MOM6 tests to not be necessary.
1665 !!
1666 !! \sa
1667 !! kpp_calculate(), kpp_applynonlocaltransport()
1668 end module mom_cvmix_kpp
mom_diag_mediator::query_averaging_enabled
logical function, public query_averaging_enabled(diag_cs, time_int, time_end)
Call this subroutine to determine whether the averaging is currently enabled. .true....
Definition: MOM_diag_mediator.F90:1850
mom_verticalgrid
Provides a transparent vertical ocean grid type and supporting routines.
Definition: MOM_verticalGrid.F90:2
mom_cvmix_kpp::kpp_get_bld
subroutine, public kpp_get_bld(CS, BLD, G)
Copies KPP surface boundary layer depth into BLD.
Definition: MOM_CVMix_KPP.F90:1469
mom_cvmix_kpp::ncall_smooth
integer ncall_smooth
Definition: MOM_CVMix_KPP.F90:68
mom_cvmix_kpp::kpp_compute_bld
subroutine, public kpp_compute_bld(CS, G, GV, US, h, Temp, Salt, u, v, EOS, uStar, buoyFlux, Waves)
Compute OBL depth.
Definition: MOM_CVMix_KPP.F90:892
mom_cvmix_kpp::sw_method_all_sw
integer, parameter, private sw_method_all_sw
Use all shortwave radiation.
Definition: MOM_CVMix_KPP.F90:50
mom_file_parser::log_version
An overloaded interface to log version information about modules.
Definition: MOM_file_parser.F90:109
mom_cvmix_kpp::kpp_end
subroutine, public kpp_end(CS)
Clear pointers, deallocate memory.
Definition: MOM_CVMix_KPP.F90:1603
mom_diag_mediator
The subroutines here provide convenient wrappers to the fms diag_manager interfaces with additional d...
Definition: MOM_diag_mediator.F90:3
mom_cvmix_kpp::lt_k_mode_rw16
integer, parameter, private lt_k_mode_rw16
Enhancement for K based on.
Definition: MOM_CVMix_KPP.F90:53
mom_cvmix_kpp::sw_method_lv1_sw
integer, parameter, private sw_method_lv1_sw
Use shortwave radiation absorbed in layer 1.
Definition: MOM_CVMix_KPP.F90:52
mom_cvmix_kpp::sw_method_mxl_sw
integer, parameter, private sw_method_mxl_sw
Use shortwave radiation absorbed in mixing layer.
Definition: MOM_CVMix_KPP.F90:51
mom_diag_mediator::register_diag_field
integer function, public register_diag_field(module_name, field_name, axes_in, init_time, long_name, units, missing_value, range, mask_variant, standard_name, verbose, do_not_log, err_msg, interp_method, tile_count, cmor_field_name, cmor_long_name, cmor_units, cmor_standard_name, cell_methods, x_cell_method, y_cell_method, v_cell_method, conversion, v_extensive)
Returns the "diag_mediator" handle for a group (native, CMOR, z-coord, ...) of diagnostics derived fr...
Definition: MOM_diag_mediator.F90:1878
mom_error_handler::mom_mesg
subroutine, public mom_mesg(message, verb, all_print)
This provides a convenient interface for writing an informative comment.
Definition: MOM_error_handler.F90:53
mom_eos
Provides subroutines for quantities specific to the equation of state.
Definition: MOM_EOS.F90:2
mom_cvmix_kpp::kpp_smooth_bld
subroutine kpp_smooth_bld(CS, G, GV, h)
Apply a 1-1-4-1-1 Laplacian filter one time on BLD to reduce any horizontal two-grid-point noise.
Definition: MOM_CVMix_KPP.F90:1350
mom_cvmix_kpp::kpp_calculate
subroutine, public kpp_calculate(CS, G, GV, US, h, uStar, buoyFlux, Kt, Ks, Kv, nonLocalTransHeat, nonLocalTransScalar, waves)
KPP vertical diffusivity/viscosity and non-local tracer transport.
Definition: MOM_CVMix_KPP.F90:590
mom_domains::pass_var
Do a halo update on an array.
Definition: MOM_domains.F90:49
mom_file_parser::param_file_type
A structure that can be parsed to read and document run-time parameters.
Definition: MOM_file_parser.F90:54
mom_file_parser::get_param
An overloaded interface to read and log the values of various types of parameters.
Definition: MOM_file_parser.F90:102
mom_cvmix_kpp::lt_k_constant
integer, parameter, private lt_k_constant
Constant enhance K through column.
Definition: MOM_CVMix_KPP.F90:53
mom_cvmix_kpp::kpp_nonlocaltransport_temp
subroutine, public kpp_nonlocaltransport_temp(CS, G, GV, h, nonLocalTrans, surfFlux, dt, scalar, C_p)
Apply KPP non-local transport of surface fluxes for temperature.
Definition: MOM_CVMix_KPP.F90:1483
mom_cvmix_kpp::kpp_nonlocaltransport_saln
subroutine, public kpp_nonlocaltransport_saln(CS, G, GV, h, nonLocalTrans, surfFlux, dt, scalar)
Apply KPP non-local transport of surface fluxes for salinity. This routine is a useful prototype for ...
Definition: MOM_CVMix_KPP.F90:1543
mom_cvmix_kpp::lt_vt2_mode_rw16
integer, parameter, private lt_vt2_mode_rw16
Enhancement for Vt2 based on.
Definition: MOM_CVMix_KPP.F90:53
mom_unit_scaling::unit_scale_type
Describes various unit conversion factors.
Definition: MOM_unit_scaling.F90:14
mom_coms
Interfaces to non-domain-oriented communication subroutines, including the MOM6 reproducing sums faci...
Definition: MOM_coms.F90:3
mom_diag_mediator::post_data
Make a diagnostic available for averaging or output.
Definition: MOM_diag_mediator.F90:70
mom_cvmix_kpp::lt_k_mode_vr12
integer, parameter, private lt_k_mode_vr12
Enhancement for K based on.
Definition: MOM_CVMix_KPP.F90:53
mom_eos::eos_type
A control structure for the equation of state.
Definition: MOM_EOS.F90:91
mom_wave_interface
Interface for surface waves.
Definition: MOM_wave_interface.F90:2
mom_wave_interface::wave_parameters_cs
Container for all surface wave related parameters.
Definition: MOM_wave_interface.F90:47
mom_cvmix_kpp::nlt_shape_cubic
integer, parameter, private nlt_shape_cubic
Cubic, .
Definition: MOM_CVMix_KPP.F90:46
mom_verticalgrid::verticalgrid_type
Describes the vertical ocean grid, including unit conversion factors.
Definition: MOM_verticalGrid.F90:24
mom_domains
Describes the decomposed MOM domain and has routines for communications across PEs.
Definition: MOM_domains.F90:2
mom_cvmix_kpp::lt_k_mode_constant
integer, parameter, private lt_k_mode_constant
Prescribed enhancement for K.
Definition: MOM_CVMix_KPP.F90:53
mom_cvmix_kpp::nlt_shape_cubic_lmd
integer, parameter, private nlt_shape_cubic_lmd
Original shape, .
Definition: MOM_CVMix_KPP.F90:47
mom_file_parser
The MOM6 facility to parse input files for runtime parameters.
Definition: MOM_file_parser.F90:2
mom_cvmix_kpp::lt_k_scaled
integer, parameter, private lt_k_scaled
Enhance K scales with G(sigma)
Definition: MOM_CVMix_KPP.F90:53
mom_grid
Provides the ocean grid type.
Definition: MOM_grid.F90:2
mom_cvmix_kpp::nlt_shape_linear
integer, parameter, private nlt_shape_linear
Linear, .
Definition: MOM_CVMix_KPP.F90:44
mom_unit_scaling
Provides a transparent unit rescaling type to facilitate dimensional consistency testing.
Definition: MOM_unit_scaling.F90:2
mom_error_handler::is_root_pe
logical function, public is_root_pe()
This returns .true. if the current PE is the root PE.
Definition: MOM_error_handler.F90:44
mom_file_parser::openparameterblock
subroutine, public openparameterblock(CS, blockName, desc)
Tags blockName onto the end of the active parameter block name.
Definition: MOM_file_parser.F90:2015
mom_debugging
Provides checksumming functions for debugging.
Definition: MOM_debugging.F90:7
mom_cvmix_kpp::kpp_init
logical function, public kpp_init(paramFile, G, GV, US, diag, Time, CS, passive, Waves)
Initialize the CVMix KPP module and set up diagnostics Returns True if KPP is to be used,...
Definition: MOM_CVMix_KPP.F90:181
mom_cvmix_kpp::lt_vt2_mode_constant
integer, parameter, private lt_vt2_mode_constant
Prescribed enhancement for Vt2.
Definition: MOM_CVMix_KPP.F90:53
mom_error_handler::mom_error
subroutine, public mom_error(level, message, all_print)
This provides a convenient interface for writing an mpp_error message with run-time filter based on a...
Definition: MOM_error_handler.F90:72
mom_cvmix_kpp::lt_vt2_mode_vr12
integer, parameter, private lt_vt2_mode_vr12
Enhancement for Vt2 based on.
Definition: MOM_CVMix_KPP.F90:53
mom_cvmix_kpp::lt_vt2_mode_lf17
integer, parameter, private lt_vt2_mode_lf17
Enhancement for Vt2 based on.
Definition: MOM_CVMix_KPP.F90:53
mom_cvmix_kpp::nlt_shape_parabolic
integer, parameter, private nlt_shape_parabolic
Parabolic, .
Definition: MOM_CVMix_KPP.F90:45
mom_file_parser::log_param
An overloaded interface to log the values of various types of parameters.
Definition: MOM_file_parser.F90:96
mom_cvmix_kpp::nlt_shape_cvmix
integer, parameter, private nlt_shape_cvmix
Use the CVMix profile.
Definition: MOM_CVMix_KPP.F90:43
mom_cvmix_kpp
Provides the K-Profile Parameterization (KPP) of Large et al., 1994, via CVMix.
Definition: MOM_CVMix_KPP.F90:2
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2
mom_cvmix_kpp::kpp_cs
Control structure for containing KPP parameters/data.
Definition: MOM_CVMix_KPP.F90:71
mom_grid::ispointincell
logical function, public ispointincell(G, i, j, x, y)
Returns true if the coordinates (x,y) are within the h-cell (i,j)
Definition: MOM_grid.F90:467
mom_grid::ocean_grid_type
Ocean grid type. See mom_grid for details.
Definition: MOM_grid.F90:26
mom_file_parser::closeparameterblock
subroutine, public closeparameterblock(CS)
Remove the lowest level of recursion from the active block name.
Definition: MOM_file_parser.F90:2033
mom_diag_mediator::diag_ctrl
The following data type a list of diagnostic fields an their variants, as well as variables that cont...
Definition: MOM_diag_mediator.F90:239
mom_wave_interface::get_langmuir_number
subroutine, public get_langmuir_number(LA, G, GV, US, HBL, ustar, i, j, H, U_H, V_H, Override_MA, Waves)
Interface to get Langmuir number based on options stored in wave structure.
Definition: MOM_wave_interface.F90:879
mom_eos::calculate_density
Calculates density of sea water from T, S and P.
Definition: MOM_EOS.F90:60