MOM_KPP.F90

Go to the documentation of this file.

!> Provides the K-Profile Parameterization (KPP) of Large et al., 1994, via CVMix.
module mom_kpp

! License goes here?

use mom_coms,          only : max_across_pes
use mom_debugging,     only : hchksum, is_nan
use mom_diag_mediator, only : time_type, diag_ctrl, safe_alloc_ptr, post_data
use mom_diag_mediator, only : query_averaging_enabled, register_diag_field
use mom_error_handler, only : mom_error, mom_mesg, fatal, warning, is_root_pe
use mom_eos,           only : eos_type, calculate_density
use mom_file_parser,   only : get_param, log_param, log_version, param_file_type
use mom_file_parser,   only : openparameterblock, closeparameterblock
use mom_grid,          only : ocean_grid_type, ispointincell
use mom_verticalgrid,  only : verticalgrid_type

use cvmix_kpp, only : cvmix_init_kpp, cvmix_put_kpp, cvmix_get_kpp_real
use cvmix_kpp, only : cvmix_coeffs_kpp
use cvmix_kpp, only : cvmix_kpp_compute_obl_depth
use cvmix_kpp, only : cvmix_kpp_compute_turbulent_scales
use cvmix_kpp, only : cvmix_kpp_compute_bulk_richardson
use cvmix_kpp, only : cvmix_kpp_compute_unresolved_shear
use cvmix_kpp, only : cvmix_kpp_params_type
use cvmix_kpp, only : cvmix_kpp_compute_kobl_depth

implicit none ; private

#include "MOM_memory.h"

public :: kpp_init
public :: kpp_calculate
public :: kpp_end
public :: kpp_nonlocaltransport_temp
public :: kpp_nonlocaltransport_saln

! Enumerated constants
integer, private, parameter :: nlt_shape_cvmix     = 0 !< Use the CVmix profile
integer, private, parameter :: nlt_shape_linear    = 1 !< Linear, \f$ G(\sigma) = 1-\sigma \f$
integer, private, parameter :: nlt_shape_parabolic = 2 !< Parabolic, \f$ G(\sigma) = (1-\sigma)^2 \f$
integer, private, parameter :: nlt_shape_cubic     = 3 !< Cubic, \f$ G(\sigma) = 1 + (2\sigma-3) \sigma^2\f$
integer, private, parameter :: nlt_shape_cubic_lmd = 4 !< Original shape, \f$ G(\sigma) = \frac{27}{4} \sigma (1-\sigma)^2 \f$

integer, private, parameter :: sw_method_all_sw = 0 !< Use all shortwave radiation
integer, private, parameter :: sw_method_mxl_sw = 1 !< Use shortwave radiation absorbed in mixing layer
integer, private, parameter :: sw_method_lv1_sw = 2 !< Use shortwave radiation absorbed in layer 1

!> Control structure for containing KPP parameters/data
type, public :: kpp_cs ; private

  ! Parameters
  real    :: ri_crit                   !< Critical bulk Richardson number (defines OBL depth)
  real    :: vonkarman                 !< von Karman constant (dimensionless)
  real    :: cs                        !< Parameter for computing velocity scale function (dimensionless)
  real    :: cs2                       !< Parameter for multiplying by non-local term
                                       !   This is active for NLT_SHAPE_CUBIC_LMD only
  logical :: enhance_diffusion         !< If True, add enhanced diffusivity at base of boundary layer.
  character(len=10) :: interptype      !< Type of interpolation in determining OBL depth
  logical :: computeekman              !< If True, compute Ekman depth limit for OBLdepth
  logical :: computemoninobukhov       !< If True, compute Monin-Obukhov limit for OBLdepth
  logical :: passivemode               !< If True, makes KPP passive meaning it does NOT alter the diffusivity
  real    :: deepobloffset             !< If non-zero, is a distance from the bottom that the OBL can not penetrate through (m)
  real    :: minobldepth               !< If non-zero, is a minimum depth for the OBL (m)
  real    :: surf_layer_ext            !< Fraction of OBL depth considered in the surface layer (nondim)
  real    :: minvtsqr                  !< Min for the squared unresolved velocity used in Rib CVMix calculation (m2/s2)
  logical :: fixedobldepth             !< If True, will fix the OBL depth at fixedOBLdepth_value
  real    :: fixedobldepth_value       !< value for the fixed OBL depth when fixedOBLdepth==True.
  logical :: debug                     !< If True, calculate checksums and write debugging information
  character(len=30) :: matchtechnique  !< Method used in CVMix for setting diffusivity and NLT profile functions
  integer :: nlt_shape                 !< MOM6 over-ride of CVMix NLT shape function
  logical :: applynonlocaltrans        !< If True, apply non-local transport to heat and scalars
  logical :: kppzerodiffusivity        !< If True, will set diffusivity and viscosity from KPP to zero; for testing purposes.
  logical :: kppisadditive             !< If True, will add KPP diffusivity to initial diffusivity.
                                       !! If False, will replace initial diffusivity wherever KPP diffusivity is non-zero.
  real    :: min_thickness             !< A minimum thickness used to avoid division by small numbers in the vicinity of vanished layers.
  ! smg: obsolete below
  logical :: correctsurflayeravg       !< If true, applies a correction to the averaging of surface layer properties
  real    :: surflayerdepth            !< A guess at the depth of the surface layer (which should 0.1 of OBLdepth) (m)
  ! smg: obsolete above
  integer :: sw_method                 !<Sets method for using shortwave radiation in surface buoyancy flux

  !> CVmix parameters
  type(cvmix_kpp_params_type), pointer :: kpp_params => null()

  ! Diagnostic handles and pointers
  type(diag_ctrl), pointer :: diag => null()
  integer :: id_obldepth = -1, id_bulkri   = -1
  integer :: id_n        = -1, id_n2       = -1
  integer :: id_ws       = -1, id_vt2      = -1
  integer :: id_bulkuz2  = -1, id_bulkdrho = -1
  integer :: id_ustar    = -1, id_buoyflux = -1
  integer :: id_qminussw = -1, id_nets     = -1
  integer :: id_sigma    = -1, id_kv_kpp   = -1
  integer :: id_kt_kpp   = -1, id_ks_kpp   = -1
  integer :: id_tsurf    = -1, id_ssurf    = -1
  integer :: id_usurf    = -1, id_vsurf    = -1
  integer :: id_kd_in    = -1
  integer :: id_nltt     = -1
  integer :: id_nlts     = -1
  integer :: id_nlt_dsdt = -1
  integer :: id_nlt_dtdt = -1
  integer :: id_nlt_temp_budget = -1
  integer :: id_nlt_saln_budget = -1

  ! Diagnostics arrays
  real, allocatable, dimension(:,:)   :: obldepth  !< Depth (positive) of OBL (m)
  real, allocatable, dimension(:,:,:) :: drho      !< Bulk difference in density (kg/m3)
  real, allocatable, dimension(:,:,:) :: uz2       !< Square of bulk difference in resolved velocity (m2/s2)
  real, allocatable, dimension(:,:,:) :: bulkri    !< Bulk Richardson number for each layer (dimensionless)
  real, allocatable, dimension(:,:,:) :: sigma     !< Sigma coordinate (dimensionless)
  real, allocatable, dimension(:,:,:) :: ws        !< Turbulent velocity scale for scalars (m/s)
  real, allocatable, dimension(:,:,:) :: n         !< Brunt-Vaisala frequency (1/s)
  real, allocatable, dimension(:,:,:) :: n2        !< Squared Brunt-Vaisala frequency (1/s2)
  real, allocatable, dimension(:,:,:) :: vt2       !< Unresolved squared turbulence velocity for bulk Ri (m2/s2)
  real, allocatable, dimension(:,:,:) :: kt_kpp    !< Temp diffusivity from KPP (m2/s)
  real, allocatable, dimension(:,:,:) :: ks_kpp    !< Scalar diffusivity from KPP (m2/s)
  real, allocatable, dimension(:,:,:) :: kv_kpp    !< Viscosity due to KPP (m2/s)
  real, allocatable, dimension(:,:)   :: tsurf     !< Temperature of surface layer (C)
  real, allocatable, dimension(:,:)   :: ssurf     !< Salinity of surface layer (ppt)
  real, allocatable, dimension(:,:)   :: usurf     !< i-velocity of surface layer (m/s)
  real, allocatable, dimension(:,:)   :: vsurf     !< j-velocity of surface layer (m/s)

end type kpp_cs

! Module data used for debugging only
logical, parameter :: verbose = .false.
#define __DO_SAFETY_CHECKS__

contains

!> Initialize the CVmix KPP module and set up diagnostics
!! Returns True if KPP is to be used, False otherwise.
logical function kpp_init(paramFile, G, diag, Time, CS, passive)

  ! Arguments
  type(param_file_type),   intent(in)    :: paramFile !< File parser
  type(ocean_grid_type),   intent(in)    :: G         !< Ocean grid
  type(diag_ctrl), target, intent(in)    :: diag      !< Diagnostics
  type(time_type),         intent(in)    :: Time      !< Time
  type(kpp_cs),            pointer       :: CS        !< Control structure
  logical, optional,       intent(out)   :: passive   !< Copy of %passiveMode

  ! Local variables
#include "version_variable.h"
  character(len=40) :: mdl = 'MOM_KPP' ! name of this module
  character(len=20) :: string          ! local temporary string
  logical :: CS_IS_ONE=.false.

  if (associated(cs)) call mom_error(fatal, 'MOM_KPP, KPP_init: '// &
           'Control structure has already been initialized')
  allocate(cs)

  ! Read parameters
  call log_version(paramfile, mdl, version, 'This is the MOM wrapper to CVmix:KPP\n' // &
            'See http://code.google.com/p/cvmix/')
  call get_param(paramfile, mdl, "USE_KPP", kpp_init, &
                 "If true, turns on the [CVmix] KPP scheme of Large et al., 1994,\n"// &
                 "to calculate diffusivities and non-local transport in the OBL.",     &
                 default=.false.)
  ! Forego remainder of initialization if not using this scheme
  if (.not. kpp_init) return

  call openparameterblock(paramfile,'KPP')
  call get_param(paramfile, mdl, 'PASSIVE', cs%passiveMode,           &
                 'If True, puts KPP into a passive-diagnostic mode.', &
                  default=.false.)
  if (present(passive)) passive=cs%passiveMode ! This is passed back to the caller so
                                               ! the caller knows to not use KPP output
  call get_param(paramfile, mdl, 'APPLY_NONLOCAL_TRANSPORT', cs%applyNonLocalTrans,  &
                 'If True, applies the non-local transport to heat and scalars.\n'//  &
                 'If False, calculates the non-local transport and tendencies but\n'//&
                 'purely for diagnostic purposes.',                                   &
                 default=.not. cs%passiveMode)
  call get_param(paramfile, mdl, 'RI_CRIT', cs%Ri_crit,                            &
                 'Critical bulk Richardson number used to define depth of the\n'// &
                 'surface Ocean Boundary Layer (OBL).',                            &
                 units='nondim', default=0.3)
  call get_param(paramfile, mdl, 'VON_KARMAN', cs%vonKarman, &
                 'von Karman constant.',                     &
                 units='nondim', default=0.40)
  call get_param(paramfile, mdl, 'ENHANCE_DIFFUSION', cs%enhance_diffusion,              &
                 'If True, adds enhanced diffusion at the based of the boundary layer.', &
                 default=.true.)
  call get_param(paramfile, mdl, 'INTERP_TYPE', cs%interpType,           &
                 'Type of interpolation to determine the OBL depth.\n'// &
                 'Allowed types are: linear, quadratic, cubic.',         &
                 default='cubic')
  call get_param(paramfile, mdl, 'COMPUTE_EKMAN', cs%computeEkman,             &
                 'If True, limit OBL depth to be no deeper than Ekman depth.', &
                 default=.false.)
  call get_param(paramfile, mdl, 'COMPUTE_MONIN_OBUKHOV', cs%computeMoninObukhov, &
                 'If True, limit the OBL depth to be no deeper than\n'//          &
                 'Monin-Obukhov depth.',                                          &
                 default=.false.)
  call get_param(paramfile, mdl, 'CS', cs%cs,                        &
                 'Parameter for computing velocity scale function.', &
                 units='nondim', default=98.96)
  call get_param(paramfile, mdl, 'CS2', cs%cs2,                        &
                 'Parameter for computing non-local term.', &
                 units='nondim', default=6.32739901508)
  call get_param(paramfile, mdl, 'DEEP_OBL_OFFSET', cs%deepOBLoffset,                             &
                 'If non-zero, the distance above the bottom to which the OBL is clipped\n'//     &
                 'if it would otherwise reach the bottom. The smaller of this and 0.1D is used.', &
                 units='m',default=0.)
  call get_param(paramfile, mdl, 'FIXED_OBLDEPTH', cs%fixedOBLdepth,       &
                 'If True, fix the OBL depth to FIXED_OBLDEPTH_VALUE\n'//  &
                 'rather than using the OBL depth from CVMix.\n'//         &
                 'This option is just for testing purposes.',              &
                 default=.false.)
  call get_param(paramfile, mdl, 'FIXED_OBLDEPTH_VALUE', cs%fixedOBLdepth_value,  &
                 'Value for the fixed OBL depth when fixedOBLdepth==True. \n'//   &
                 'This parameter is for just for testing purposes. \n'//          &
                 'It will over-ride the OBLdepth computed from CVMix.',           &
                 units='m',default=30.0)
  call get_param(paramfile, mdl, 'SURF_LAYER_EXTENT', cs%surf_layer_ext,   &
                 'Fraction of OBL depth considered in the surface layer.', &
                 units='nondim',default=0.10)
  call get_param(paramfile, mdl, 'MINIMUM_OBL_DEPTH', cs%minOBLdepth,                            &
                 'If non-zero, a minimum depth to use for KPP OBL depth. Independent of\n'//     &
                 'this parameter, the OBL depth is always at least as deep as the first layer.', &
                 units='m',default=0.)
  call get_param(paramfile, mdl, 'MINIMUM_VT2', cs%minVtsqr,                                   &
                 'Min of the unresolved velocity Vt2 used in Rib CVMix calculation.     \n'//  &
                 'Scaling: MINIMUM_VT2 = const1*d*N*ws, with d=1m, N=1e-5/s, ws=1e-6 m/s.',    &
                 units='m2/s2',default=1e-10)

! smg: for removal below
  call get_param(paramfile, mdl, 'CORRECT_SURFACE_LAYER_AVERAGE', cs%correctSurfLayerAvg,   &
                 'If true, applies a correction step to the averaging of surface layer\n'// &
                 'properties. This option is obsolete.', default=.false.)
  call get_param(paramfile, mdl, 'FIRST_GUESS_SURFACE_LAYER_DEPTH', cs%surfLayerDepth,              &
                 'The first guess at the depth of the surface layer used for averaging\n'//         &
                 'the surface layer properties. If =0, the top model level properties\n'//          &
                 'will be used for the surface layer. If CORRECT_SURFACE_LAYER_AVERAGE=True, a\n'// &
                 'subsequent correction is applied. This parameter is obsolete', units='m', default=0.)
! smg: for removal above

  call get_param(paramfile, mdl, 'NLT_SHAPE', string, &
                 'MOM6 method to set nonlocal transport profile.\n'//                          &
                 'Over-rides the result from CVMix.  Allowed values are: \n'//                 &
                 '\t CVMIX     - Uses the profiles from CVmix specified by MATCH_TECHNIQUE\n'//&
                 '\t LINEAR    - A linear profile, 1-sigma\n'//                                &
                 '\t PARABOLIC - A parablic profile, (1-sigma)^2\n'//                          &
                 '\t CUBIC     - A cubic profile, (1-sigma)^2(1+2*sigma)\n'//                  &
                 '\t CUBIC_LMD - The original KPP profile',                                    &
                 default='CVMIX')
  select case ( trim(string) )
    case ("CVMIX")     ; cs%NLT_shape = nlt_shape_cvmix
    case ("LINEAR")    ; cs%NLT_shape = nlt_shape_linear
    case ("PARABOLIC") ; cs%NLT_shape = nlt_shape_parabolic
    case ("CUBIC")     ; cs%NLT_shape = nlt_shape_cubic
    case ("CUBIC_LMD") ; cs%NLT_shape = nlt_shape_cubic_lmd
    case default ; call mom_error(fatal,"KPP_init: "// &
                   "Unrecognized NLT_SHAPE option"//trim(string))
  end select
  call get_param(paramfile, mdl, 'MATCH_TECHNIQUE', cs%MatchTechnique,                                    &
                 'CVMix method to set profile function for diffusivity and NLT,\n'//                      &
                 'as well as matching across OBL base. Allowed values are: \n'//                          &
                 '\t SimpleShapes      = sigma*(1-sigma)^2 for both diffusivity and NLT\n'//              &
                 '\t MatchGradient     = sigma*(1-sigma)^2 for NLT; diffusivity profile from matching\n'//&
                 '\t MatchBoth         = match gradient for both diffusivity and NLT\n'//                 &
                 '\t ParabolicNonLocal = sigma*(1-sigma)^2 for diffusivity; (1-sigma)^2 for NLT',         &
                 default='SimpleShapes')
  if (cs%MatchTechnique.eq.'ParabolicNonLocal') then
     ! This forces Cs2 (Cs in non-local computation) to equal 1 for parabolic non-local option.
     !  May be used during CVmix initialization.
     cs_is_one=.true.
  endif
  call get_param(paramfile, mdl, 'KPP_ZERO_DIFFUSIVITY', cs%KPPzeroDiffusivity,            &
                 'If True, zeroes the KPP diffusivity and viscosity; for testing purpose.',&
                 default=.false.)
  call get_param(paramfile, mdl, 'KPP_IS_ADDITIVE', cs%KPPisAdditive,                &
                 'If true, adds KPP diffusivity to diffusivity from other schemes.'//&
                 'If false, KPP is the only diffusivity wherever KPP is non-zero.',  &
                 default=.true.)
  call get_param(paramfile, mdl, 'KPP_SHORTWAVE_METHOD',string,                      &
                 'Determines contribution of shortwave radiation to KPP surface '// &
                 'buoyancy flux.  Options include:\n'//                             &
                 '  ALL_SW: use total shortwave radiation\n'//                      &
                 '  MXL_SW:  use shortwave radiation absorbed by mixing layer\n'//  &
                 '  LV1_SW:  use shortwave radiation absorbed by top model layer',  &
                 default='MXL_SW')
  select case ( trim(string) )
    case ("ALL_SW") ; cs%SW_METHOD = sw_method_all_sw
    case ("MXL_SW") ; cs%SW_METHOD = sw_method_mxl_sw
    case ("LV1_SW") ; cs%SW_METHOD = sw_method_lv1_sw
    case default ; call mom_error(fatal,"KPP_init: "// &
                   "Unrecognized KPP_SHORTWAVE_METHOD option"//trim(string))
  end select
  call get_param(paramfile, mdl, 'CVMIX_ZERO_H_WORK_AROUND', cs%min_thickness,                           &
                 'A minimum thickness used to avoid division by small numbers in the vicinity\n'//       &
                 'of vanished layers. This is independent of MIN_THICKNESS used in other parts of MOM.', &
                 units='m', default=0.)

  call closeparameterblock(paramfile)
  call get_param(paramfile, mdl, 'DEBUG', cs%debug, default=.false., do_not_log=.true.)

  call cvmix_init_kpp( ri_crit=cs%Ri_crit,                 &
                       minobldepth=cs%minOBLdepth,         &
                       minvtsqr=cs%minVtsqr,               &
                       vonkarman=cs%vonKarman,             &
                       surf_layer_ext=cs%surf_layer_ext,   &
                       interp_type=cs%interpType,          &
                       lekman=cs%computeEkman,             &
                       lmonob=cs%computeMoninObukhov,      &
                       matchtechnique=cs%MatchTechnique,   &
                       lenhanced_diff=cs%enhance_diffusion,&
                       lnonzero_surf_nonlocal=cs_is_one   ,&
                       cvmix_kpp_params_user=cs%KPP_params )

  ! Register diagnostics
  cs%diag => diag
  cs%id_OBLdepth = register_diag_field('ocean_model', 'KPP_OBLdepth', diag%axesT1, time, &
      'Thickness of the surface Ocean Boundary Layer calculated by [CVmix] KPP', 'meter', &
      cmor_field_name='oml', cmor_long_name='ocean_mixed_layer_thickness_defined_by_mixing_scheme', &
      cmor_units='m', cmor_standard_name='Ocean Mixed Layer Thickness Defined by Mixing Scheme')
      ! CMOR names are placeholders; must be modified by time period
      ! for CMOR compliance. Diag manager will be used for omlmax and
      ! omldamax.
  cs%id_BulkDrho = register_diag_field('ocean_model', 'KPP_BulkDrho', diag%axesTL, time, &
      'Bulk difference in density used in Bulk Richardson number, as used by [CVmix] KPP', 'kg/m3')
  cs%id_BulkUz2 = register_diag_field('ocean_model', 'KPP_BulkUz2', diag%axesTL, time, &
      'Square of bulk difference in resolved velocity used in Bulk Richardson number via [CVmix] KPP', 'm2/s2')
  cs%id_BulkRi = register_diag_field('ocean_model', 'KPP_BulkRi', diag%axesTL, time, &
      'Bulk Richardson number used to find the OBL depth used by [CVmix] KPP', 'nondim')
  cs%id_Sigma = register_diag_field('ocean_model', 'KPP_sigma', diag%axesTi, time, &
      'Sigma coordinate used by [CVmix] KPP', 'nondim')
  cs%id_Ws = register_diag_field('ocean_model', 'KPP_Ws', diag%axesTL, time, &
      'Turbulent vertical velocity scale for scalars used by [CVmix] KPP', 'm/s')
  cs%id_N = register_diag_field('ocean_model', 'KPP_N', diag%axesTi, time, &
      '(Adjusted) Brunt-Vaisala frequency used by [CVmix] KPP', '1/s')
  cs%id_N2 = register_diag_field('ocean_model', 'KPP_N2', diag%axesTi, time, &
      'Square of Brunt-Vaisala frequency used by [CVmix] KPP', '1/s2')
  cs%id_Vt2 = register_diag_field('ocean_model', 'KPP_Vt2', diag%axesTL, time, &
      'Unresolved shear turbulence used by [CVmix] KPP', 'm2/s2')
  cs%id_uStar = register_diag_field('ocean_model', 'KPP_uStar', diag%axesT1, time, &
      'Friction velocity, u*, as used by [CVmix] KPP', 'm/s')
  cs%id_buoyFlux = register_diag_field('ocean_model', 'KPP_buoyFlux', diag%axesTi, time, &
      'Surface (and penetrating) buoyancy flux, as used by [CVmix] KPP', 'm2/s3')
  cs%id_QminusSW = register_diag_field('ocean_model', 'KPP_QminusSW', diag%axesT1, time, &
      'Net temperature flux ignoring short-wave, as used by [CVmix] KPP', 'K m/s')
  cs%id_netS = register_diag_field('ocean_model', 'KPP_netSalt', diag%axesT1, time, &
      'Effective net surface salt flux, as used by [CVmix] KPP', 'ppt m/s')
  cs%id_Kt_KPP = register_diag_field('ocean_model', 'KPP_Kheat', diag%axesTi, time, &
      'Heat diffusivity due to KPP, as calculated by [CVmix] KPP', 'm2/s')
  cs%id_Kd_in = register_diag_field('ocean_model', 'KPP_Kd_in', diag%axesTi, time, &
      'Diffusivity passed to KPP', 'm2/s')
  cs%id_Ks_KPP = register_diag_field('ocean_model', 'KPP_Ksalt', diag%axesTi, time, &
      'Salt diffusivity due to KPP, as calculated by [CVmix] KPP', 'm2/s')
  cs%id_Kv_KPP = register_diag_field('ocean_model', 'KPP_Kv', diag%axesTi, time, &
      'Vertical viscosity due to KPP, as calculated by [CVmix] KPP', 'm2/s')
  cs%id_NLTt = register_diag_field('ocean_model', 'KPP_NLtransport_heat', diag%axesTi, time, &
      'Non-local transport (Cs*G(sigma)) for heat, as calculated by [CVmix] KPP', 'nondim')
  cs%id_NLTs = register_diag_field('ocean_model', 'KPP_NLtransport_salt', diag%axesTi, time, &
      'Non-local tranpsort (Cs*G(sigma)) for scalars, as calculated by [CVmix] KPP', 'nondim')
  cs%id_NLT_dTdt = register_diag_field('ocean_model', 'KPP_NLT_dTdt', diag%axesTL, time, &
      'Temperature tendency due to non-local transport of heat, as calculated by [CVmix] KPP', 'K/s')
  cs%id_NLT_dSdt = register_diag_field('ocean_model', 'KPP_NLT_dSdt', diag%axesTL, time, &
      'Salinity tendency due to non-local transport of salt, as calculated by [CVmix] KPP', 'ppt/s')
  cs%id_NLT_temp_budget = register_diag_field('ocean_model', 'KPP_NLT_temp_budget', diag%axesTL, time, &
      'Heat content change due to non-local transport, as calculated by [CVmix] KPP', 'W/m^2')
  cs%id_NLT_saln_budget = register_diag_field('ocean_model', 'KPP_NLT_saln_budget', diag%axesTL, time, &
      'Salt content change due to non-local transport, as calculated by [CVmix] KPP', 'kg/(sec*m^2)')
  cs%id_Tsurf = register_diag_field('ocean_model', 'KPP_Tsurf', diag%axesT1, time, &
      'Temperature of surface layer (10% of OBL depth) as passed to [CVmix] KPP', 'C')
  cs%id_Ssurf = register_diag_field('ocean_model', 'KPP_Ssurf', diag%axesT1, time, &
      'Salinity of surface layer (10% of OBL depth) as passed to [CVmix] KPP', 'ppt')
  cs%id_Usurf = register_diag_field('ocean_model', 'KPP_Usurf', diag%axesCu1, time, &
      'i-component flow of surface layer (10% of OBL depth) as passed to [CVmix] KPP', 'm/s')
  cs%id_Vsurf = register_diag_field('ocean_model', 'KPP_Vsurf', diag%axesCv1, time, &
      'j-component flow of surface layer (10% of OBL depth) as passed to [CVmix] KPP', 'm/s')

  if (cs%id_OBLdepth > 0) allocate( cs%OBLdepth( szi_(g), szj_(g) ) )
  if (cs%id_OBLdepth > 0) cs%OBLdepth(:,:) = 0.
  if (cs%id_BulkDrho > 0) allocate( cs%dRho( szi_(g), szj_(g), szk_(g) ) )
  if (cs%id_BulkDrho > 0) cs%dRho(:,:,:) = 0.
  if (cs%id_BulkUz2 > 0)  allocate( cs%Uz2( szi_(g), szj_(g), szk_(g) ) )
  if (cs%id_BulkUz2 > 0)  cs%Uz2(:,:,:) = 0.
  if (cs%id_BulkRi > 0)   allocate( cs%BulkRi( szi_(g), szj_(g), szk_(g) ) )
  if (cs%id_BulkRi > 0)   cs%BulkRi(:,:,:) = 0.
  if (cs%id_Sigma > 0)    allocate( cs%sigma( szi_(g), szj_(g), szk_(g)+1 ) )
  if (cs%id_Sigma > 0)    cs%sigma(:,:,:) = 0.
  if (cs%id_Ws > 0)       allocate( cs%Ws( szi_(g), szj_(g), szk_(g) ) )
  if (cs%id_Ws > 0)       cs%Ws(:,:,:) = 0.
  if (cs%id_N > 0)        allocate( cs%N( szi_(g), szj_(g), szk_(g)+1 ) )
  if (cs%id_N > 0)        cs%N(:,:,:) = 0.
  if (cs%id_N2 > 0)       allocate( cs%N2( szi_(g), szj_(g), szk_(g)+1 ) )
  if (cs%id_N2 > 0)       cs%N2(:,:,:) = 0.
  if (cs%id_Vt2 > 0)      allocate( cs%Vt2( szi_(g), szj_(g), szk_(g) ) )
  if (cs%id_Vt2 > 0)      cs%Vt2(:,:,:) = 0.
  if (cs%id_Kt_KPP > 0)   allocate( cs%Kt_KPP( szi_(g), szj_(g), szk_(g)+1 ) )
  if (cs%id_Kt_KPP > 0)   cs%Kt_KPP(:,:,:) = 0.
  if (cs%id_Ks_KPP > 0)   allocate( cs%Ks_KPP( szi_(g), szj_(g), szk_(g)+1 ) )
  if (cs%id_Ks_KPP > 0)   cs%Ks_KPP(:,:,:) = 0.
  if (cs%id_Kv_KPP > 0)   allocate( cs%Kv_KPP( szi_(g), szj_(g), szk_(g)+1 ) )
  if (cs%id_Kv_KPP > 0)   cs%Kv_KPP(:,:,:) = 0.
  if (cs%id_Tsurf > 0)    allocate( cs%Tsurf( szi_(g), szj_(g)) )
  if (cs%id_Tsurf > 0)    cs%Tsurf(:,:) = 0.
  if (cs%id_Ssurf > 0)    allocate( cs%Ssurf( szi_(g), szj_(g)) )
  if (cs%id_Ssurf > 0)    cs%Ssurf(:,:) = 0.
  if (cs%id_Usurf > 0)    allocate( cs%Usurf( szib_(g), szj_(g)) )
  if (cs%id_Usurf > 0)    cs%Usurf(:,:) = 0.
  if (cs%id_Vsurf > 0)    allocate( cs%Vsurf( szi_(g), szjb_(g)) )
  if (cs%id_Vsurf > 0)    cs%Vsurf(:,:) = 0.

end function kpp_init



!> KPP vertical diffusivity/viscosity and non-local tracer transport
subroutine kpp_calculate(CS, G, GV, h, Temp, Salt, u, v, EOS, uStar, &
                         buoyFlux, Kt, Ks, Kv, nonLocalTransHeat,&
                         nonLocalTransScalar)

  ! Arguments
  type(kpp_cs),                           pointer       :: CS             !< Control structure
  type(ocean_grid_type),                  intent(in)    :: G              !< Ocean grid
  type(verticalgrid_type),                intent(in)    :: GV             !< Ocean vertical grid
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h              !< Layer/level thicknesses (units of H)
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: Temp           !< potential/cons temp (deg C)
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: Salt           !< Salinity (ppt)
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in)    :: u              !< Velocity i-component (m/s)
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in)    :: v              !< Velocity j-component (m/s)
  type(eos_type),                         pointer       :: EOS            !< Equation of state
  real, dimension(SZI_(G),SZJ_(G)),         intent(in)    :: uStar          !< Surface friction velocity (m/s)
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in)    :: buoyFlux !< Surface buoyancy flux (m2/s3)
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: Kt       !< (in)  Vertical diffusivity of heat w/o KPP (m2/s)
                                                                          !< (out) Vertical diffusivity including KPP (m2/s)
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: Ks       !< (in)  Vertical diffusivity of salt w/o KPP (m2/s)
                                                                          !< (out) Vertical diffusivity including KPP (m2/s)
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: Kv       !< (in)  Vertical viscosity w/o KPP (m2/s)
                                                                          !< (out) Vertical viscosity including KPP (m2/s)
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: nonLocalTransHeat   !< Temp non-local transport (m/s)
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: nonLocalTransScalar !< scalar non-local transport (m/s)

  ! Local variables
  integer :: i, j, k, km1,kp1                    ! Loop indices
  real, dimension( G%ke )     :: cellHeight      ! Cell center heights referenced to surface (m) (negative in ocean)
  real, dimension( G%ke+1 )   :: iFaceHeight     ! Interface heights referenced to surface (m) (negative in ocean)
  real, dimension( G%ke+1 )   :: N2_1d           ! Brunt-Vaisala frequency squared, at interfaces (1/s2)
  real, dimension( G%ke+1 )   :: N_1d            ! Brunt-Vaisala frequency at interfaces (1/s) (floored at 0)
  real, dimension( G%ke )     :: Ws_1d           ! Profile of vertical velocity scale for scalars (m/s)
  real, dimension( G%ke )     :: Wm_1d           ! Profile of vertical velocity scale for momentum (m/s)
  real, dimension( G%ke )     :: Vt2_1d          ! Unresolved velocity for bulk Ri calculation/diagnostic (m2/s2)
  real, dimension( G%ke )     :: BulkRi_1d       ! Bulk Richardson number for each layer
  real, dimension( G%ke )     :: deltaRho        ! delta Rho in numerator of Bulk Ri number
  real, dimension( G%ke )     :: deltaU2         ! square of delta U (shear) in denominator of Bulk Ri (m2/s2)
  real, dimension( G%ke+1, 2) :: Kdiffusivity    ! Vertical diffusivity at interfaces (m2/s)
  real, dimension( G%ke+1 )   :: Kviscosity      ! Vertical viscosity at interfaces (m2/s)
  real, dimension( G%ke+1, 2) :: nonLocalTrans   ! Non-local transport for heat/salt at interfaces (non-dimensional)
  real, dimension( G%ke )     :: surfBuoyFlux2

  ! for EOS calculation
  real, dimension( 3*G%ke )   :: rho_1D
  real, dimension( 3*G%ke )   :: pres_1D
  real, dimension( 3*G%ke )   :: Temp_1D
  real, dimension( 3*G%ke )   :: Salt_1D

  real :: kOBL, OBLdepth_0d, surfFricVel, surfBuoyFlux, Coriolis
  real :: GoRho, pRef, rho1, rhoK, rhoKm1, Uk, Vk, sigma

  real :: zBottomMinusOffset   ! Height of bottom plus a little bit (m)
  real :: SLdepth_0d           ! Surface layer depth = surf_layer_ext*OBLdepth.
  real :: hTot                 ! Running sum of thickness used in the surface layer average (m)
  real :: delH                 ! Thickness of a layer (m)
  real :: surfHtemp, surfTemp  ! Integral and average of temp over the surface layer
  real :: surfHsalt, surfSalt  ! Integral and average of saln over the surface layer
  real :: surfHu, surfU        ! Integral and average of u over the surface layer
  real :: surfHv, surfV        ! Integral and average of v over the surface layer
  real :: dh    ! The local thickness used for calculating interface positions (m)
  real :: hcorr ! A cumulative correction arising from inflation of vanished layers (m)
  integer :: kk, ksfc, ktmp

#ifdef __DO_SAFETY_CHECKS__
  if (cs%debug) then
    call hchksum(h, "KPP in: h",g%HI,haloshift=0, scale=gv%H_to_m)
    call hchksum(temp, "KPP in: T",g%HI,haloshift=0)
    call hchksum(salt, "KPP in: S",g%HI,haloshift=0)
    call hchksum(u, "KPP in: u",g%HI,haloshift=0)
    call hchksum(v, "KPP in: v",g%HI,haloshift=0)
    call hchksum(ustar, "KPP in: uStar",g%HI,haloshift=0)
    call hchksum(buoyflux, "KPP in: buoyFlux",g%HI,haloshift=0)
    call hchksum(kt, "KPP in: Kt",g%HI,haloshift=0)
    call hchksum(ks, "KPP in: Ks",g%HI,haloshift=0)
  endif
#endif

  ! some constants
  gorho = gv%g_Earth / gv%Rho0
  nonlocaltrans(:,:) = 0.0

  if (cs%id_Kd_in > 0) call post_data(cs%id_Kd_in, kt, cs%diag)

!$OMP parallel do default(none) shared(G,GV,CS,EOS,uStar,Temp,Salt,u,v,h,GoRho,       &
!$OMP                                  buoyFlux, nonLocalTransHeat,                   &
!$OMP                                  nonLocalTransScalar,Kt,Ks,Kv)                  &
!$OMP                     firstprivate(nonLocalTrans)                                 &
!$OMP                          private(Coriolis,surfFricVel,SLdepth_0d,hTot,surfTemp, &
!$OMP                                  surfHtemp,surfSalt,surfHsalt,surfU,            &
!$OMP                                  surfHu,surfV,surfHv,iFaceHeight,               &
!$OMP                                  pRef,km1,cellHeight,Uk,Vk,deltaU2,             &
!$OMP                                  rho1,rhoK,rhoKm1,deltaRho,N2_1d,N_1d,delH,     &
!$OMP                                  surfBuoyFlux,Ws_1d,Vt2_1d,BulkRi_1d,           &
!$OMP                                  OBLdepth_0d,zBottomMinusOffset,Kdiffusivity,   &
!$OMP                                  Kviscosity,sigma,kOBL,kk,pres_1D,Temp_1D,      &
!$OMP                                  Salt_1D,rho_1D,surfBuoyFlux2,ksfc,dh,hcorr)

  ! loop over horizontal points on processor
  do j = g%jsc, g%jec
    do i = g%isc, g%iec

      ! skip calling KPP for land points
      if (g%mask2dT(i,j)==0.) cycle

      ! things independent of position within the column
      coriolis = 0.25*( (g%CoriolisBu(i,j)   + g%CoriolisBu(i-1,j-1)) &
                       +(g%CoriolisBu(i-1,j) + g%CoriolisBu(i,j-1)) )
      surffricvel = ustar(i,j)

      ! Bullk Richardson number computed for each cell in a column,
      ! assuming OBLdepth = grid cell depth. After Rib(k) is
      ! known for the column, then CVMix interpolates to find
      ! the actual OBLdepth. This approach avoids need to iterate
      ! on the OBLdepth calculation. It follows that used in MOM5
      ! and POP.
      ifaceheight(1) = 0.0 ! BBL is all relative to the surface
      pref = 0.
      hcorr = 0.
      do k=1,g%ke

        ! cell center and cell bottom in meters (negative values in the ocean)
        dh = h(i,j,k) * gv%H_to_m ! Nominal thickness to use for increment
        dh = dh + hcorr ! Take away the accumulated error (could temporarily make dh<0)
        hcorr = min( dh - cs%min_thickness, 0. ) ! If inflating then hcorr<0
        dh = max( dh, cs%min_thickness ) ! Limit increment dh>=min_thickness
        cellheight(k)    = ifaceheight(k) - 0.5 * dh
        ifaceheight(k+1) = ifaceheight(k) - dh

        ! find ksfc for cell where "surface layer" sits
        sldepth_0d = cs%surf_layer_ext*max( max(-cellheight(k),-ifaceheight(2) ), cs%minOBLdepth )
        ksfc = k
        do ktmp = 1,k
          if (-1.0*ifaceheight(ktmp+1) >= sldepth_0d) then
            ksfc = ktmp
            exit
          endif
        enddo

        ! average temp, saln, u, v over surface layer
        ! use C-grid average to get u,v on T-points.
        surfhtemp=0.0
        surfhsalt=0.0
        surfhu   =0.0
        surfhv   =0.0
        htot     =0.0
        do ktmp = 1,ksfc

          ! SLdepth_0d can be between cell interfaces
          delh = min( max(0.0, sldepth_0d - htot), h(i,j,ktmp)*gv%H_to_m )

          ! surface layer thickness
          htot = htot + delh

          ! surface averaged fields
          surfhtemp = surfhtemp + temp(i,j,ktmp) * delh
          surfhsalt = surfhsalt + salt(i,j,ktmp) * delh
          surfhu    = surfhu + 0.5*(u(i,j,ktmp)+u(i-1,j,ktmp)) * delh
          surfhv    = surfhv + 0.5*(v(i,j,ktmp)+v(i,j-1,ktmp)) * delh

        enddo
        surftemp = surfhtemp / htot
        surfsalt = surfhsalt / htot
        surfu    = surfhu    / htot
        surfv    = surfhv    / htot

        ! vertical shear between present layer and
        ! surface layer averaged surfU,surfV.
        ! C-grid average to get Uk and Vk on T-points.
        uk         = 0.5*(u(i,j,k)+u(i-1,j,k)) - surfu
        vk         = 0.5*(v(i,j,k)+v(i,j-1,k)) - surfv
        deltau2(k) = uk**2 + vk**2

        ! pressure, temp, and saln for EOS
        ! kk+1 = surface fields
        ! kk+2 = k fields
        ! kk+3 = km1 fields
        km1  = max(1, k-1)
        kk   = 3*(k-1)
        pres_1d(kk+1) = pref
        pres_1d(kk+2) = pref
        pres_1d(kk+3) = pref
        temp_1d(kk+1) = surftemp
        temp_1d(kk+2) = temp(i,j,k)
        temp_1d(kk+3) = temp(i,j,km1)
        salt_1d(kk+1) = surfsalt
        salt_1d(kk+2) = salt(i,j,k)
        salt_1d(kk+3) = salt(i,j,km1)

        ! pRef is pressure at interface between k and km1.
        ! iterate pRef for next pass through k-loop.
        pref = pref + gv%H_to_Pa * h(i,j,k)

        ! this difference accounts for penetrating SW
        surfbuoyflux2(k) = buoyflux(i,j,1) - buoyflux(i,j,k+1)

      enddo ! k-loop finishes

      ! compute in-situ density
      call calculate_density(temp_1d, salt_1d, pres_1d, rho_1d, 1, 3*g%ke, eos)

      ! N2 (can be negative) and N (non-negative) on interfaces.
      ! deltaRho is non-local rho difference used for bulk Richardson number.
      ! N_1d is local N (with floor) used for unresolved shear calculation.
      do k = 1, g%ke
        km1 = max(1, k-1)
        kk = 3*(k-1)
        deltarho(k) = rho_1d(kk+2) - rho_1d(kk+1)
        n2_1d(k)    = (gorho * (rho_1d(kk+2) - rho_1d(kk+3)) ) / &
                      ((0.5*(h(i,j,km1) + h(i,j,k))+gv%H_subroundoff)*gv%H_to_m)
        n_1d(k)     = sqrt( max( n2_1d(k), 0.) )
      enddo
      n2_1d(g%ke+1 ) = 0.0
      n_1d(g%ke+1 )  = 0.0

      ! turbulent velocity scales w_s and w_m computed at the cell centers.
      ! Note that if sigma > CS%surf_layer_ext, then CVmix_kpp_compute_turbulent_scales
      ! computes w_s and w_m velocity scale at sigma=CS%surf_layer_ext. So we only pass
      ! sigma=CS%surf_layer_ext for this calculation.
      call cvmix_kpp_compute_turbulent_scales( &
        cs%surf_layer_ext, & ! (in)  Normalized surface layer depth; sigma = CS%surf_layer_ext
        -cellheight,       & ! (in)  Assume here that OBL depth (m) = -cellHeight(k)
        surfbuoyflux2,     & ! (in)  Buoyancy flux at surface (m2/s3)
        surffricvel,       & ! (in)  Turbulent friction velocity at surface (m/s)
        w_s=ws_1d,         & ! (out) Turbulent velocity scale profile (m/s)
        cvmix_kpp_params_user=cs%KPP_params )

      ! Calculate Bulk Richardson number from eq (21) of LMD94
      bulkri_1d = cvmix_kpp_compute_bulk_richardson( &
                  cellheight(1:g%ke),                & ! Depth of cell center (m)
                  gorho*deltarho,                    & ! Bulk buoyancy difference, Br-B(z) (1/s)
                  deltau2,                           & ! Square of resolved velocity difference (m2/s2)
                  ws_cntr=ws_1d,                     & ! Turbulent velocity scale profile (m/s)
                  n_iface=n_1d)                        ! Buoyancy frequency (1/s)


      surfbuoyflux = buoyflux(i,j,1) ! This is only used in kpp_compute_OBL_depth to limit
                                     ! h to Monin-Obukov (default is false, ie. not used)

      call cvmix_kpp_compute_obl_depth( &
        bulkri_1d,              & ! (in) Bulk Richardson number
        ifaceheight,            & ! (in) Height of interfaces (m)
        obldepth_0d,            & ! (out) OBL depth (m)
        kobl,                   & ! (out) level (+fraction) of OBL extent
        zt_cntr=cellheight,     & ! (in) Height of cell centers (m)
        surf_fric=surffricvel,  & ! (in) Turbulent friction velocity at surface (m/s)
        surf_buoy=surfbuoyflux, & ! (in) Buoyancy flux at surface (m2/s3)
        coriolis=coriolis,      & ! (in) Coriolis parameter (1/s)
        cvmix_kpp_params_user=cs%KPP_params ) ! KPP parameters

      ! A hack to avoid KPP reaching the bottom. It was needed during development
      ! because KPP was unable to handle vanishingly small layers near the bottom.
      if (cs%deepOBLoffset>0.) then
        zbottomminusoffset = ifaceheight(g%ke+1) + min(cs%deepOBLoffset,-0.1*ifaceheight(g%ke+1))
        obldepth_0d = min( obldepth_0d, -zbottomminusoffset )
      endif

      ! apply some constraints on OBLdepth
      if(cs%fixedOBLdepth)  obldepth_0d = cs%fixedOBLdepth_value
      obldepth_0d = max( obldepth_0d, -ifaceheight(2) )      ! no shallower than top layer
      obldepth_0d = min( obldepth_0d, -ifaceheight(g%ke+1) ) ! no deeper than bottom
      kobl        = cvmix_kpp_compute_kobl_depth( ifaceheight, cellheight, obldepth_0d )

!*************************************************************************
! smg: remove code below

! Following "correction" step has been found to be unnecessary.
! Code should be removed after further testing.
      if (cs%correctSurfLayerAvg) then

        sldepth_0d = cs%surf_layer_ext * obldepth_0d
        htot      = h(i,j,1)
        surftemp  = temp(i,j,1) ; surfhtemp = surftemp * htot
        surfsalt  = salt(i,j,1) ; surfhsalt = surfsalt * htot
        surfu     = 0.5*(u(i,j,1)+u(i-1,j,1)) ; surfhu = surfu * htot
        surfv     = 0.5*(v(i,j,1)+v(i,j-1,1)) ; surfhv = surfv * htot
        pref      = 0.0

        do k = 2, g%ke

          ! Recalculate differences with surface layer
          uk = 0.5*(u(i,j,k)+u(i-1,j,k)) - surfu
          vk = 0.5*(v(i,j,k)+v(i,j-1,k)) - surfv
          deltau2(k) = uk**2 + vk**2
          pref = pref + gv%H_to_Pa * h(i,j,k)
          call calculate_density(surftemp, surfsalt, pref, rho1, eos)
          call calculate_density(temp(i,j,k), salt(i,j,k), pref, rhok, eos)
          deltarho(k) = rhok - rho1

          ! Surface layer averaging (needed for next k+1 iteration of this loop)
          if (htot < sldepth_0d) then
            delh = min( max(0., sldepth_0d - htot), h(i,j,k)*gv%H_to_m )
            htot = htot + delh
            surfhtemp = surfhtemp + temp(i,j,k) * delh ; surftemp = surfhtemp / htot
            surfhsalt = surfhsalt + salt(i,j,k) * delh ; surfsalt = surfhsalt / htot
            surfhu = surfhu + 0.5*(u(i,j,k)+u(i-1,j,k)) * delh ; surfu = surfhu / htot
            surfhv = surfhv + 0.5*(v(i,j,k)+v(i,j-1,k)) * delh ; surfv = surfhv / htot
          endif

        enddo

        bulkri_1d = cvmix_kpp_compute_bulk_richardson( &
                    cellheight(1:g%ke),                & ! Depth of cell center (m)
                    gorho*deltarho,                    & ! Bulk buoyancy difference, Br-B(z) (1/s)
                    deltau2,                           & ! Square of resolved velocity difference (m2/s2)
                    ws_cntr=ws_1d,                     & ! Turbulent velocity scale profile (m/s)
                    n_iface=n_1d )                       ! Buoyancy frequency (1/s)

        surfbuoyflux = buoyflux(i,j,1) ! This is only used in kpp_compute_OBL_depth to limit
                                       ! h to Monin-Obukov (default is false, ie. not used)

        call cvmix_kpp_compute_obl_depth( &
          bulkri_1d,              & ! (in) Bulk Richardson number
          ifaceheight,            & ! (in) Height of interfaces (m)
          obldepth_0d,            & ! (out) OBL depth (m)
          kobl,                   & ! (out) level (+fraction) of OBL extent
          zt_cntr=cellheight,     & ! (in) Height of cell centers (m)
          surf_fric=surffricvel,  & ! (in) Turbulent friction velocity at surface (m/s)
          surf_buoy=surfbuoyflux, & ! (in) Buoyancy flux at surface (m2/s3)
          coriolis=coriolis,      & ! (in) Coriolis parameter (1/s)
          cvmix_kpp_params_user=cs%KPP_params ) ! KPP parameters

        if (cs%deepOBLoffset>0.) then
          zbottomminusoffset = ifaceheight(g%ke+1) + min(cs%deepOBLoffset,-0.1*ifaceheight(g%ke+1))
          obldepth_0d = min( obldepth_0d, -zbottomminusoffset )
          kobl = cvmix_kpp_compute_kobl_depth( ifaceheight, cellheight, obldepth_0d )
        endif

        ! apply some constraints on OBLdepth
        if(cs%fixedOBLdepth)  obldepth_0d = cs%fixedOBLdepth_value
        obldepth_0d = max( obldepth_0d, -ifaceheight(2) )      ! no shallower than top layer
        obldepth_0d = min( obldepth_0d, -ifaceheight(g%ke+1) ) ! no deep than bottom
        kobl        = cvmix_kpp_compute_kobl_depth( ifaceheight, cellheight, obldepth_0d )

      endif   ! endif for "correction" step

! smg: remove code above
! **********************************************************************


      ! Call CVMix/KPP to obtain OBL diffusivities, viscosities and non-local transports

      ! Unlike LMD94, we do not match to interior diffusivities. If using the original
      ! LMD94 shape function, not matching is equivalent to matching to a zero diffusivity.

      !BGR/ Add option for use of surface buoyancy flux with total sw flux.
      if (cs%SW_METHOD .eq. sw_method_all_sw) then
         surfbuoyflux = buoyflux(i,j,1)
      elseif (cs%SW_METHOD .eq. sw_method_mxl_sw) then
         surfbuoyflux  = buoyflux(i,j,1) - buoyflux(i,j,int(kobl)+1) ! We know the actual buoyancy flux into the OBL
      elseif (cs%SW_METHOD .eq. sw_method_lv1_sw) then
         surfbuoyflux  = buoyflux(i,j,1) - buoyflux(i,j,2)
      endif

      ! If option "MatchBoth" is selected in CVMix, MOM should be capable of matching.
      if (.not. (cs%MatchTechnique.eq.'MatchBoth')) then
         kdiffusivity(:,:) = 0. ! Diffusivities for heat and salt (m2/s)
         kviscosity(:)     = 0. ! Viscosity (m2/s)
      else
         kdiffusivity(:,1) = kt(i,j,:)
         kdiffusivity(:,2) = ks(i,j,:)
         kviscosity(:)=kv(i,j,:)
      endif

      call cvmix_coeffs_kpp(kviscosity,        & ! (inout) Total viscosity (m2/s)
                            kdiffusivity(:,1), & ! (inout) Total heat diffusivity (m2/s)
                            kdiffusivity(:,2), & ! (inout) Total salt diffusivity (m2/s)
                            ifaceheight,       & ! (in) Height of interfaces (m)
                            cellheight,        & ! (in) Height of level centers (m)
                            kviscosity,        & ! (in) Original viscosity (m2/s)
                            kdiffusivity(:,1), & ! (in) Original heat diffusivity (m2/s)
                            kdiffusivity(:,2), & ! (in) Original salt diffusivity (m2/s)
                            obldepth_0d,       & ! (in) OBL depth (m)
                            kobl,              & ! (in) level (+fraction) of OBL extent
                            nonlocaltrans(:,1),& ! (out) Non-local heat transport (non-dimensional)
                            nonlocaltrans(:,2),& ! (out) Non-local salt transport (non-dimensional)
                            surffricvel,       & ! (in) Turbulent friction velocity at surface (m/s)
                            surfbuoyflux,      & ! (in) Buoyancy flux at surface (m2/s3)
                            g%ke,              & ! (in) Number of levels to compute coeffs for
                            g%ke,              & ! (in) Number of levels in array shape
                            cvmix_kpp_params_user=cs%KPP_params )


      ! Over-write CVMix NLT shape function with one of the following choices.
      ! The CVMix code has yet to update for thse options, so we compute in MOM6.
      ! Note that nonLocalTrans = Cs * G(sigma) (LMD94 notation), with
      ! Cs = 6.32739901508.
      ! Start do-loop at k=2, since k=1 is ocean surface (sigma=0)
      ! and we do not wish to double-count the surface forcing.
      ! Only compute nonlocal transport for 0 <= sigma <= 1.
      ! MOM6 recommended shape is the parabolic; it gives deeper boundary layer
      ! and no spurious extrema.
      if (surfbuoyflux < 0.0) then
        if (cs%NLT_shape == nlt_shape_cubic) then
          do k = 2, g%ke
            sigma = min(1.0,-ifaceheight(k)/obldepth_0d)
            nonlocaltrans(k,1) = (1.0 - sigma)**2 * (1.0 + 2.0*sigma) !*
            nonlocaltrans(k,2) = nonlocaltrans(k,1)
          enddo
        elseif (cs%NLT_shape == nlt_shape_parabolic) then
          do k = 2, g%ke
            sigma = min(1.0,-ifaceheight(k)/obldepth_0d)
            nonlocaltrans(k,1) = (1.0 - sigma)**2 !*CS%CS2
            nonlocaltrans(k,2) = nonlocaltrans(k,1)
          enddo
        elseif (cs%NLT_shape == nlt_shape_linear) then
          do k = 2, g%ke
            sigma = min(1.0,-ifaceheight(k)/obldepth_0d)
            nonlocaltrans(k,1) = (1.0 - sigma)!*CS%CS2
            nonlocaltrans(k,2) = nonlocaltrans(k,1)
          enddo
        elseif (cs%NLT_shape == nlt_shape_cubic_lmd) then
          ! Sanity check (should agree with CVMix result using simple matching)
          do k = 2, g%ke
            sigma = min(1.0,-ifaceheight(k)/obldepth_0d)
            nonlocaltrans(k,1) = cs%CS2 * sigma*(1.0 -sigma)**2
            nonlocaltrans(k,2) = nonlocaltrans(k,1)
          enddo
        endif
      endif

      ! we apply nonLocalTrans in subroutines
      ! KPP_NonLocalTransport_temp and KPP_NonLocalTransport_saln
      nonlocaltransheat(i,j,:)   = nonlocaltrans(:,1) ! temp
      nonlocaltransscalar(i,j,:) = nonlocaltrans(:,2) ! saln

      ! set the KPP diffusivity and viscosity to zero for testing purposes
      if(cs%KPPzeroDiffusivity) then
         kdiffusivity(:,1) = 0.0
         kdiffusivity(:,2) = 0.0
         kviscosity(:)     = 0.0
      endif

      ! recompute wscale for diagnostics, now that we in fact know boundary layer depth
      if (cs%id_Ws > 0) then
          call cvmix_kpp_compute_turbulent_scales( &
            -cellheight/obldepth_0d,               & ! (in)  Normalized boundary layer coordinate
            obldepth_0d,                           & ! (in)  OBL depth (m)
            surfbuoyflux,                          & ! (in)  Buoyancy flux at surface (m2/s3)
            surffricvel,                           & ! (in)  Turbulent friction velocity at surface (m/s)
            w_s=ws_1d,                             & ! (out) Turbulent velocity scale profile (m/s)
            cvmix_kpp_params_user=cs%KPP_params    & !       KPP parameters
            )
          cs%Ws(i,j,:) = ws_1d(:)
      endif

      ! compute unresolved squared velocity for diagnostics
      if (cs%id_Vt2 > 0) then
        vt2_1d(:) = cvmix_kpp_compute_unresolved_shear( &
                    cellheight(1:g%ke),                 & ! Depth of cell center (m)
                    ws_cntr=ws_1d,                      & ! Turbulent velocity scale profile, at centers (m/s)
                    n_iface=n_1d,                       & ! Buoyancy frequency at interface (1/s)
                    cvmix_kpp_params_user=cs%KPP_params ) ! KPP parameters
        cs%Vt2(i,j,:) = vt2_1d(:)
      endif

      ! Copy 1d data into 3d diagnostic arrays
      if (cs%id_OBLdepth > 0) cs%OBLdepth(i,j) = obldepth_0d
      if (cs%id_BulkDrho > 0) cs%dRho(i,j,:)   = deltarho(:)
      if (cs%id_BulkUz2 > 0)  cs%Uz2(i,j,:)    = deltau2(:)
      if (cs%id_BulkRi > 0)   cs%BulkRi(i,j,:) = bulkri_1d(:)
      if (cs%id_sigma > 0) then
        cs%sigma(i,j,:)  = 0.
        if (obldepth_0d>0.)   cs%sigma(i,j,:)  = -ifaceheight/obldepth_0d
      endif
      if (cs%id_N      > 0)   cs%N(i,j,:)      = n_1d(:)
      if (cs%id_N2     > 0)   cs%N2(i,j,:)     = n2_1d(:)
      if (cs%id_Kt_KPP > 0)   cs%Kt_KPP(i,j,:) = kdiffusivity(:,1)
      if (cs%id_Ks_KPP > 0)   cs%Ks_KPP(i,j,:) = kdiffusivity(:,2)
      if (cs%id_Kv_KPP > 0)   cs%Kv_KPP(i,j,:) = kviscosity(:)
      if (cs%id_Tsurf  > 0)   cs%Tsurf(i,j)    = surftemp
      if (cs%id_Ssurf  > 0)   cs%Ssurf(i,j)    = surfsalt
      if (cs%id_Usurf  > 0)   cs%Usurf(i,j)    = surfu
      if (cs%id_Vsurf  > 0)   cs%Vsurf(i,j)    = surfv


       ! Update output of routine
      if (.not. cs%passiveMode) then
        if (cs%KPPisAdditive) then
          do k=1, g%ke+1
            kt(i,j,k) = kt(i,j,k) + kdiffusivity(k,1)
            ks(i,j,k) = ks(i,j,k) + kdiffusivity(k,2)
            kv(i,j,k) = kv(i,j,k) + kviscosity(k)
          enddo
        else ! KPP replaces prior diffusivity when former is non-zero
          do k=1, g%ke+1
            if (kdiffusivity(k,1) /= 0.) kt(i,j,k) = kdiffusivity(k,1)
            if (kdiffusivity(k,2) /= 0.) ks(i,j,k) = kdiffusivity(k,2)
            if (kviscosity(k) /= 0.) kv(i,j,k) = kviscosity(k)
          enddo
        endif
      endif


    ! end of the horizontal do-loops over the vertical columns
    enddo ! i
  enddo ! j


#ifdef __DO_SAFETY_CHECKS__
  if (cs%debug) then
    call hchksum(kt, "KPP out: Kt",g%HI,haloshift=0)
    call hchksum(ks, "KPP out: Ks",g%HI,haloshift=0)
  endif
#endif

  ! send diagnostics to post_data
  if (cs%id_OBLdepth > 0) call post_data(cs%id_OBLdepth, cs%OBLdepth,        cs%diag)
  if (cs%id_BulkDrho > 0) call post_data(cs%id_BulkDrho, cs%dRho,            cs%diag)
  if (cs%id_BulkUz2  > 0) call post_data(cs%id_BulkUz2,  cs%Uz2,             cs%diag)
  if (cs%id_BulkRi   > 0) call post_data(cs%id_BulkRi,   cs%BulkRi,          cs%diag)
  if (cs%id_sigma    > 0) call post_data(cs%id_sigma,    cs%sigma,           cs%diag)
  if (cs%id_Ws       > 0) call post_data(cs%id_Ws,       cs%Ws,              cs%diag)
  if (cs%id_N        > 0) call post_data(cs%id_N,        cs%N,               cs%diag)
  if (cs%id_N2       > 0) call post_data(cs%id_N2,       cs%N2,              cs%diag)
  if (cs%id_Vt2      > 0) call post_data(cs%id_Vt2,      cs%Vt2,             cs%diag)
  if (cs%id_uStar    > 0) call post_data(cs%id_uStar,    ustar,              cs%diag)
  if (cs%id_buoyFlux > 0) call post_data(cs%id_buoyFlux, buoyflux,           cs%diag)
  if (cs%id_Kt_KPP   > 0) call post_data(cs%id_Kt_KPP,   cs%Kt_KPP,          cs%diag)
  if (cs%id_Ks_KPP   > 0) call post_data(cs%id_Ks_KPP,   cs%Ks_KPP,          cs%diag)
  if (cs%id_Kv_KPP   > 0) call post_data(cs%id_Kv_KPP,   cs%Kv_KPP,          cs%diag)
  if (cs%id_NLTt     > 0) call post_data(cs%id_NLTt,     nonlocaltransheat,  cs%diag)
  if (cs%id_NLTs     > 0) call post_data(cs%id_NLTs,     nonlocaltransscalar,cs%diag)
  if (cs%id_Tsurf    > 0) call post_data(cs%id_Tsurf,    cs%Tsurf,           cs%diag)
  if (cs%id_Ssurf    > 0) call post_data(cs%id_Ssurf,    cs%Ssurf,           cs%diag)
  if (cs%id_Usurf    > 0) call post_data(cs%id_Usurf,    cs%Usurf,           cs%diag)
  if (cs%id_Vsurf    > 0) call post_data(cs%id_Vsurf,    cs%Vsurf,           cs%diag)

end subroutine kpp_calculate



!> Apply KPP non-local transport of surface fluxes for temperature.
subroutine kpp_nonlocaltransport_temp(CS, G, GV, h, nonLocalTrans, surfFlux, &
                                      dt, scalar, C_p)

  type(kpp_cs),                               intent(in)    :: CS            !< Control structure
  type(ocean_grid_type),                      intent(in)    :: G             !< Ocean grid
  type(verticalgrid_type),                    intent(in)    :: GV            !< Ocean vertical grid
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   intent(in)    :: h             !< Layer/level thickness (units of H)
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in)    :: nonLocalTrans !< Non-local transport (non-dimensional)
  real, dimension(SZI_(G),SZJ_(G)),           intent(in)    :: surfFlux      !< Surface flux of scalar (H/s * scalar)
  real,                                       intent(in)    :: dt            !< Time-step (s)
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   intent(inout) :: scalar        !< temperature
  real,                                       intent(in)    :: C_p           !< Seawater specific heat capacity (J/(kg*K))

  integer :: i, j, k
  real, dimension( SZI_(G), SZJ_(G), SZK_(G) ) :: dtracer


  dtracer(:,:,:) = 0.0
!$OMP parallel do default(none) shared(G,GV,dtracer,nonLocalTrans,h,surfFlux,CS,scalar,dt)
  do k = 1, g%ke
    do j = g%jsc, g%jec
      do i = g%isc, g%iec
        dtracer(i,j,k) = ( nonlocaltrans(i,j,k) - nonlocaltrans(i,j,k+1) ) / &
                         ( h(i,j,k) + gv%H_subroundoff ) * surfflux(i,j)
      enddo
    enddo
  enddo

  !  Update tracer due to non-local redistribution of surface flux
  if (cs%applyNonLocalTrans) then
    do k = 1, g%ke
      do j = g%jsc, g%jec
        do i = g%isc, g%iec
          scalar(i,j,k) = scalar(i,j,k) + dt * dtracer(i,j,k)
        enddo
      enddo
    enddo
  endif

  ! Diagnostics
  if (cs%id_QminusSW        > 0) call post_data(cs%id_QminusSW, surfflux, cs%diag)
  if (cs%id_NLT_dTdt        > 0) call post_data(cs%id_NLT_dTdt, dtracer,  cs%diag)
  if (cs%id_NLT_temp_budget > 0) then
    dtracer(:,:,:) = 0.0
    do k = 1, g%ke
      do j = g%jsc, g%jec
        do i = g%isc, g%iec
          dtracer(i,j,k) = (nonlocaltrans(i,j,k) - nonlocaltrans(i,j,k+1)) * &
                           surfflux(i,j) * c_p * gv%H_to_kg_m2
        enddo
      enddo
    enddo
    call post_data(cs%id_NLT_temp_budget, dtracer, cs%diag)
  endif

end subroutine kpp_nonlocaltransport_temp


!> Apply KPP non-local transport of surface fluxes for salinity.
!> This routine is a useful prototype for other material tracers.
subroutine kpp_nonlocaltransport_saln(CS, G, GV, h, nonLocalTrans, surfFlux, dt, scalar)

  type(kpp_cs),                               intent(in)    :: CS            !< Control structure
  type(ocean_grid_type),                      intent(in)    :: G             !< Ocean grid
  type(verticalgrid_type),                    intent(in)    :: GV            !< Ocean vertical grid
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   intent(in)    :: h             !< Layer/level thickness (units of H)
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(in)    :: nonLocalTrans !< Non-local transport (non-dimensional)
  real, dimension(SZI_(G),SZJ_(G)),           intent(in)    :: surfFlux      !< Surface flux of scalar (H/s * scalar)
  real,                                       intent(in)    :: dt            !< Time-step (s)
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   intent(inout) :: scalar        !< Scalar (scalar units)

  integer :: i, j, k
  real, dimension( SZI_(G), SZJ_(G), SZK_(G) ) :: dtracer


  dtracer(:,:,:) = 0.0
!$OMP parallel do default(none) shared(G,GV,dtracer,nonLocalTrans,h,surfFlux,CS,scalar,dt)
  do k = 1, g%ke
    do j = g%jsc, g%jec
      do i = g%isc, g%iec
        dtracer(i,j,k) = ( nonlocaltrans(i,j,k) - nonlocaltrans(i,j,k+1) ) / &
                         ( h(i,j,k) + gv%H_subroundoff ) * surfflux(i,j)
      enddo
    enddo
  enddo

  !  Update tracer due to non-local redistribution of surface flux
  if (cs%applyNonLocalTrans) then
    do k = 1, g%ke
      do j = g%jsc, g%jec
        do i = g%isc, g%iec
          scalar(i,j,k) = scalar(i,j,k) + dt * dtracer(i,j,k)
        enddo
      enddo
    enddo
  endif

  ! Diagnostics
  if (cs%id_netS            > 0) call post_data(cs%id_netS,     surfflux, cs%diag)
  if (cs%id_NLT_dSdt        > 0) call post_data(cs%id_NLT_dSdt, dtracer,  cs%diag)
  if (cs%id_NLT_saln_budget > 0) then
    dtracer(:,:,:) = 0.0
    do k = 1, g%ke
      do j = g%jsc, g%jec
        do i = g%isc, g%iec
          dtracer(i,j,k) = (nonlocaltrans(i,j,k) - nonlocaltrans(i,j,k+1)) * &
                           surfflux(i,j) * gv%H_to_kg_m2
        enddo
      enddo
    enddo
    call post_data(cs%id_NLT_saln_budget, dtracer, cs%diag)
  endif

end subroutine kpp_nonlocaltransport_saln




!> Clear pointers, deallocate memory
subroutine kpp_end(CS)
  type(kpp_cs), pointer :: CS !< Control structure

  deallocate(cs)
end subroutine kpp_end

!> \namespace mom_kpp
!!
!! \section section_KPP The K-Profile Parameterization
!!
!! The K-Profile Parameterization (KPP) of Large et al., 1994, (http://dx.doi.org/10.1029/94RG01872) is
!! implemented via the Community Vertical Mixing package, [CVmix](https://code.google.com/p/cvmix),
!! which is called directly by this module.
!!
!! The formulation and implementation of KPP is described in great detail in the
!! [CVMix manual](https://cvmix.googlecode.com/svn/trunk/manual/cvmix.pdf) (written by our own Stephen Griffies).
!!
!! \subsection section_KPP_nutshell KPP in a nutshell
!!
!! Large et al., 1994, decompose the parameterized boundary layer turbulent flux of a scalar, \f$ s \f$, as
!! \f[ \overline{w^\prime s^\prime} = -K \partial_z s + K \gamma_s(\sigma), \f]
!! where \f$ \sigma = -z/h \f$ is a non-dimensional coordinate within the boundary layer of depth \f$ h \f$.
!! \f$ K \f$ is the eddy diffusivity and is a function of position within the boundary layer as well as a
!! function of the surface forcing:
!! \f[ K = h w_s(\sigma) G(\sigma) . \f]
!! Here, \f$ w_s \f$ is the vertical velocity scale of the boundary layer turbulence and \f$ G(\sigma) \f$ is
!! a "shape function" which is described later.
!! The last term is the "non-local transport" which involves a function \f$ \gamma_s(\sigma) \f$ that is matched
!! to the forcing but is not actually needed in the final implementation.
!! Instead, the entire non-local transport term can be equivalently written
!! \f[ K \gamma_s(\sigma) = C_s G(\sigma) Q_s \f]
!! where \f$ Q_s \f$ is the surface flux of \f$ s \f$ and \f$ C_s \f$ is a constant.
!! The vertical structure of the redistribution (non-local) term is solely due  to the shape function, \f$ G(\sigma) \f$.
!! In our implementation of KPP, we allow the shape functions used for \f$ K \f$ and for the non-local transport
!! to be chosen independently.
!!
!! [google_thread_NLT]: https://groups.google.com/forum/#!msg/cvmix-dev/i6rF-eHOtKI/Ti8BeyksrhAJ "Extreme values of non-local transport"
!!
!! The particular shape function most widely used in the atmospheric community is
!! \f[ G(\sigma) = \sigma (1-\sigma)^2 \f]
!! which satisfies the boundary conditions
!!  \f$ G(0) = 0 \f$,
!!  \f$ G(1) = 0 \f$,
!!  \f$ G^\prime(0) = 1 \f$, and
!!  \f$ G^\prime(1) = 0 \f$.
!! Large et al, 1994, alter the function so as to match interior diffusivities but we have found that this leads
!! to inconsistencies within the formulation (see google groups thread [Extreme values of non-local transport][google_thread_NLT]).
!! Instead, we use either the above form, or even simpler forms that use alternative upper boundary conditions.
!!
!! The KPP boundary layer depth is a function of the bulk Richardson number, Rib.
!! But to compute Rib, we need the boundary layer depth.  To address this circular
!! logic, we compute Rib for each vertical cell in a column, assuming the BL depth
!! equals to the depth of the given grid cell.  Once we have a vertical array of Rib(k),
!! we then call the OBLdepth routine from CVMix to compute the actual
!! OBLdepth. We optionally then "correct" the OBLdepth by cycling through once more,
!! this time knowing the OBLdepth from the first pass. This "correction" step is not
!! used by NCAR. It has been found in idealized MOM6 tests to not be necessary.
!!
!! \sa
!! kpp_calculate(), kpp_applynonlocaltransport()
end module mom_kpp