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