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