MOM_CVMix_conv.F90

!> Interface to CVMix convection scheme.
module mom_cvmix_conv

! This file is part of MOM6. See LICENSE.md for the license.

use mom_debugging,      only : hchksum
use mom_diag_mediator,  only : diag_ctrl, time_type, register_diag_field
use mom_diag_mediator,  only : post_data
use mom_eos,            only : calculate_density
use mom_error_handler,  only : mom_error, is_root_pe, fatal, warning, note
use mom_file_parser,    only : openparameterblock, closeparameterblock
use mom_file_parser,    only : get_param, log_version, param_file_type
use mom_grid,           only : ocean_grid_type
use mom_unit_scaling,   only : unit_scale_type
use mom_variables,      only : thermo_var_ptrs
use mom_verticalgrid,   only : verticalgrid_type
use cvmix_convection,   only : cvmix_init_conv, cvmix_coeffs_conv
use cvmix_kpp,          only : cvmix_kpp_compute_kobl_depth

implicit none ; private

#include <MOM_memory.h>

public cvmix_conv_init, calculate_cvmix_conv, cvmix_conv_end, cvmix_conv_is_used

!> Control structure including parameters for CVMix convection.
type, public :: cvmix_conv_cs ; private

  ! Parameters
  real    :: kd_conv_const !< diffusivity constant used in convective regime [m2 s-1]
  real    :: kv_conv_const !< viscosity constant used in convective regime [m2 s-1]
  real    :: bv_sqr_conv   !< Threshold for squared buoyancy frequency
                           !! needed to trigger Brunt-Vaisala parameterization [s-2]
  real    :: min_thickness !< Minimum thickness allowed [m]
  logical :: debug         !< If true, turn on debugging

  ! Daignostic handles and pointers
  type(diag_ctrl), pointer :: diag => null() !< Pointer to diagnostics control structure
  !>@{ Diagnostics handles
  integer :: id_n2 = -1, id_kd_conv = -1, id_kv_conv = -1
  !>@}

end type cvmix_conv_cs

character(len=40)  :: mdl = "MOM_CVMix_conv"     !< This module's name.

contains

!> Initialized the CVMix convection mixing routine.
logical function cvmix_conv_init(Time, G, GV, US, param_file, diag, CS)

  type(time_type),         intent(in)    :: time       !< The current time.
  type(ocean_grid_type),   intent(in)    :: g          !< Grid structure.
  type(verticalgrid_type), intent(in)    :: gv         !< Vertical grid structure.
  type(unit_scale_type),   intent(in)    :: us         !< A dimensional unit scaling type
  type(param_file_type),   intent(in)    :: param_file !< Run-time parameter file handle
  type(diag_ctrl), target, intent(inout) :: diag       !< Diagnostics control structure.
  type(cvmix_conv_cs),     pointer       :: cs         !< This module's control structure.
  ! Local variables
  real    :: prandtl_conv !< Turbulent Prandtl number used in convective instabilities.
  logical :: useepbl      !< If True, use the ePBL boundary layer scheme.

! This include declares and sets the variable "version".
#include "version_variable.h"

  if (associated(cs)) then
    call mom_error(warning, "CVMix_conv_init called with an associated "// &
                            "control structure.")
    return
  endif
  allocate(cs)

  ! Read parameters
  call get_param(param_file, mdl, "USE_CVMix_CONVECTION", cvmix_conv_init, default=.false., do_not_log=.true.)
  call log_version(param_file, mdl, version, &
           "Parameterization of enhanced mixing due to convection via CVMix", &
           all_default=.not.cvmix_conv_init)
  call get_param(param_file, mdl, "USE_CVMix_CONVECTION", cvmix_conv_init, &
                 "If true, turns on the enhanced mixing due to convection "//&
                 "via CVMix. This scheme increases diapycnal diffs./viscs. "//&
                 "at statically unstable interfaces. Relevant parameters are "//&
                 "contained in the CVMix_CONVECTION% parameter block.", &
                 default=.false.)

  if (.not. cvmix_conv_init) return

  call get_param(param_file, mdl, "ENERGETICS_SFC_PBL", useepbl, default=.false., &
                do_not_log=.true.)

  ! Warn user if EPBL is being used, since in this case mixing due to convection will
  ! be aplied in the boundary layer
  if (useepbl) then
     call mom_error(warning, 'MOM_CVMix_conv_init: '// &
           'CVMix convection may not be properly applied when ENERGETICS_SFC_PBL = True'//&
           'as convective mixing might occur in the boundary layer.')
  endif

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

  call get_param(param_file, mdl, 'MIN_THICKNESS', cs%min_thickness, default=0.001, do_not_log=.true.)

  call openparameterblock(param_file,'CVMix_CONVECTION')

  call get_param(param_file, mdl, "PRANDTL_CONV", prandtl_conv, &
                 "The turbulent Prandtl number applied to convective "//&
                 "instabilities (i.e., used to convert KD_CONV into KV_CONV)", &
                 units="nondim", default=1.0)

  call get_param(param_file, mdl, 'KD_CONV', cs%kd_conv_const, &
                 "Diffusivity used in convective regime. Corresponding viscosity "//&
                 "(KV_CONV) will be set to KD_CONV * PRANDTL_TURB.", &
                 units='m2/s', default=1.00)

  call get_param(param_file, mdl, 'BV_SQR_CONV', cs%bv_sqr_conv, &
                 "Threshold for squared buoyancy frequency needed to trigger "//&
                 "Brunt-Vaisala parameterization.", &
                 units='1/s^2', default=0.0)

  call closeparameterblock(param_file)

  ! set kv_conv_const based on kd_conv_const and prandtl_conv
  cs%kv_conv_const = cs%kd_conv_const * prandtl_conv

  ! Register diagnostics
  cs%diag => diag
  cs%id_N2 = register_diag_field('ocean_model', 'N2_conv', diag%axesTi, time, &
      'Square of Brunt-Vaisala frequency used by MOM_CVMix_conv module', '1/s2', conversion=us%s_to_T**2)
  cs%id_kd_conv = register_diag_field('ocean_model', 'kd_conv', diag%axesTi, time, &
      'Additional diffusivity added by MOM_CVMix_conv module', 'm2/s', conversion=us%Z2_T_to_m2_s)
  cs%id_kv_conv = register_diag_field('ocean_model', 'kv_conv', diag%axesTi, time, &
      'Additional viscosity added by MOM_CVMix_conv module', 'm2/s', conversion=us%Z2_T_to_m2_s)

  call cvmix_init_conv(convect_diff=cs%kd_conv_const, &
                       convect_visc=cs%kv_conv_const, &
                       lbruntvaisala=.true.,    &
                       bvsqr_convect=cs%bv_sqr_conv)

end function cvmix_conv_init

!> Subroutine for calculating enhanced diffusivity/viscosity
!! due to convection via CVMix
subroutine calculate_cvmix_conv(h, tv, G, GV, US, CS, hbl, Kd, Kv, Kd_aux)

  type(ocean_grid_type),                     intent(in)  :: g  !< Grid structure.
  type(verticalgrid_type),                   intent(in)  :: gv !< Vertical grid structure.
  type(unit_scale_type),                     intent(in)  :: us !< A dimensional unit scaling type
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)  :: h  !< Layer thickness [H ~> m or kg m-2].
  type(thermo_var_ptrs),                     intent(in)  :: tv !< Thermodynamics structure.
  type(cvmix_conv_cs),                       pointer     :: cs !< The control structure returned
                                                                !! by a previous call to CVMix_conv_init.
  real, dimension(SZI_(G),SZJ_(G)),          intent(in)  :: hbl !< Depth of ocean boundary layer [Z ~> m]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), &
                                   optional, intent(inout) :: kd !< Diapycnal diffusivity at each interface that
                                                                 !! will be incremented here [Z2 T-1 ~> m2 s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), &
                                   optional, intent(inout) :: kv !< Viscosity at each interface that will be
                                                                 !! incremented here [Z2 T-1 ~> m2 s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), &
                                   optional, intent(inout) :: kd_aux !< A second diapycnal diffusivity at each
                                                                 !! interface that will also be incremented
                                                                 !! here [Z2 T-1 ~> m2 s-1].

  ! local variables
  real, dimension(SZK_(GV)) :: rho_lwr !< Adiabatic Water Density, this is a dummy
                                       !! variable since here convection is always
                                       !! computed based on Brunt Vaisala.
  real, dimension(SZK_(GV)) :: rho_1d  !< water density in a column, this is also
                                       !! a dummy variable, same reason as above.
  real, dimension(SZK_(GV)+1) :: n2    !< Squared buoyancy frequency [s-2]
  real, dimension(SZK_(GV)+1) :: kv_col !< Viscosities at interfaces in the column [m2 s-1]
  real, dimension(SZK_(GV)+1) :: kd_col !< Diffusivities at interfaces in the column [m2 s-1]
  real, dimension(SZK_(GV)+1) :: ifaceheight !< Height of interfaces [m]
  real, dimension(SZK_(GV))   :: cellheight  !< Height of cell centers [m]
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1) :: &
    kd_conv, &                         !< Diffusivity added by convection for diagnostics [Z2 T-1 ~> m2 s-1]
    kv_conv, &                         !< Viscosity added by convection for diagnostics [Z2 T-1 ~> m2 s-1]
    n2_3d                              !< Squared buoyancy frequency for diagnostics [N-2 ~> s-2]
  integer :: kobl                      !< level of OBL extent
  real :: g_o_rho0  ! Gravitational acceleration divided by density times unit convserion factors
                    ! [Z s-2 R-1 ~> m4 s-2 kg-1]
  real :: pref      ! Interface pressures [R L2 T-2 ~> Pa]
  real :: rhok, rhokm1 ! In situ densities of the layers above and below at the interface pressure [R ~> kg m-3]
  real :: hbl_kpp   ! The depth of the ocean boundary as used by KPP [m]
  real :: dz        ! A thickness [Z ~> m]
  real :: dh, hcorr ! Two thicknesses [m]
  integer :: i, j, k

  g_o_rho0 = us%L_to_Z**2*us%s_to_T**2 * gv%g_Earth / gv%Rho0

  ! initialize dummy variables
  rho_lwr(:) = 0.0 ; rho_1d(:) = 0.0

  ! set N2 to zero at the top- and bottom-most interfaces
  n2(1) = 0.0 ; n2(gv%ke+1) = 0.0

  if (cs%id_N2 > 0) n2_3d(:,:,:) = 0.0
  if (cs%id_kv_conv > 0) kv_conv(:,:,:) = 0.0
  if (cs%id_kd_conv > 0) kd_conv(:,:,:) = 0.0

  do j = g%jsc, g%jec
    do i = g%isc, g%iec

      ! skip calling at land points
      !if (G%mask2dT(i,j) == 0.) cycle

      pref = 0. ; if (associated(tv%p_surf)) pref = tv%p_surf(i,j)
      ! Compute Brunt-Vaisala frequency (static stability) on interfaces
      do k=2,gv%ke

        ! pRef is pressure at interface between k and km1 [R L2 T-2 ~> Pa].
        pref = pref + (gv%H_to_RZ*gv%g_Earth) * h(i,j,k)
        call calculate_density(tv%t(i,j,k), tv%s(i,j,k), pref, rhok, tv%eqn_of_state)
        call calculate_density(tv%t(i,j,k-1), tv%s(i,j,k-1), pref, rhokm1, tv%eqn_of_state)

        dz = ((0.5*(h(i,j,k-1) + h(i,j,k))+gv%H_subroundoff)*gv%H_to_Z)
        n2(k) = g_o_rho0 * (rhok - rhokm1) / dz ! Can be negative

      enddo

      ifaceheight(1) = 0.0 ! BBL is all relative to the surface
      hcorr = 0.0
      ! compute heights at cell center and interfaces
      do k=1,gv%ke
        dh = h(i,j,k) * gv%H_to_m ! Nominal thickness to use for increment, in the units used by CVMix.
        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
      enddo

      ! gets index of the level and interface above hbl
      hbl_kpp = us%Z_to_m*hbl(i,j)  ! Convert to the units used by CVMix.
      kobl = cvmix_kpp_compute_kobl_depth(ifaceheight, cellheight, hbl_kpp)

      kv_col(:) = 0.0 ; kd_col(:) = 0.0
      call cvmix_coeffs_conv(mdiff_out=kv_col(:), &
                             tdiff_out=kd_col(:), &
                             nsqr=n2(:), &
                             dens=rho_1d(:), &
                             dens_lwr=rho_lwr(:), &
                             nlev=gv%ke,    &
                             max_nlev=gv%ke, &
                             obl_ind=kobl)

      if (present(kd)) then
        ! Increment the diffusivity outside of the boundary layer.
        do k=max(1,kobl+1),gv%ke+1
          kd(i,j,k) = kd(i,j,k) + us%m2_s_to_Z2_T * kd_col(k)
        enddo
      endif
      if (present(kd_aux)) then
        ! Increment the other diffusivity outside of the boundary layer.
        do k=max(1,kobl+1),gv%ke+1
          kd_aux(i,j,k) = kd_aux(i,j,k) + us%m2_s_to_Z2_T * kd_col(k)
        enddo
      endif

      if (present(kv)) then
        ! Increment the viscosity outside of the boundary layer.
        do k=max(1,kobl+1),gv%ke+1
          kv(i,j,k) = kv(i,j,k) + us%m2_s_to_Z2_T * kv_col(k)
        enddo
      endif

      ! Store 3-d arrays for diagnostics.
      if (cs%id_kv_conv > 0) then
        ! Do not apply mixing due to convection within the boundary layer
        do k=max(1,kobl+1),gv%ke+1
          kv_conv(i,j,k) = us%m2_s_to_Z2_T * kv_col(k)
        enddo
      endif
      if (cs%id_kd_conv > 0) then
        ! Do not apply mixing due to convection within the boundary layer
        do k=max(1,kobl+1),gv%ke+1
          kd_conv(i,j,k) = us%m2_s_to_Z2_T * kd_col(k)
        enddo
      endif

      if (cs%id_N2 > 0) then ; do k=2,gv%ke ; n2_3d(i,j,k) = us%T_to_s**2*n2(k) ; enddo ; endif

    enddo
  enddo

  if (cs%debug) then
    ! if (CS%id_N2 > 0) call hchksum(N2_3d, "MOM_CVMix_conv: N2",G%HI,haloshift=0)
    ! if (CS%id_kd_conv > 0) &
    !   call hchksum(Kd_conv, "MOM_CVMix_conv: Kd_conv", G%HI, haloshift=0, scale=US%Z2_T_to_m2_s)
    ! if (CS%id_kv_conv > 0) &
    !   call hchksum(Kv_conv, "MOM_CVMix_conv: Kv_conv", G%HI, haloshift=0, scale=US%m2_s_to_Z2_T)
    if (present(kd)) call hchksum(kv, "MOM_CVMix_conv: Kd", g%HI, haloshift=0, scale=us%m2_s_to_Z2_T)
    if (present(kv)) call hchksum(kv, "MOM_CVMix_conv: Kv", g%HI, haloshift=0, scale=us%m2_s_to_Z2_T)
  endif

  ! send diagnostics to post_data
  if (cs%id_N2 > 0) call post_data(cs%id_N2, n2_3d, cs%diag)
  if (cs%id_kd_conv > 0) call post_data(cs%id_kd_conv, kd_conv, cs%diag)
  if (cs%id_kv_conv > 0) call post_data(cs%id_kv_conv, kv_conv, cs%diag)

end subroutine calculate_cvmix_conv

!> Reads the parameter "USE_CVMix_CONVECTION" and returns state.
!! This function allows other modules to know whether this parameterization will
!! be used without needing to duplicate the log entry.
logical function cvmix_conv_is_used(param_file)
  type(param_file_type), intent(in) :: param_file !< A structure to parse for run-time parameters
  call get_param(param_file, mdl, "USE_CVMix_CONVECTION", cvmix_conv_is_used, &
                 default=.false., do_not_log = .true.)

end function cvmix_conv_is_used

!> Clear pointers and dealocate memory
subroutine cvmix_conv_end(CS)
  type(cvmix_conv_cs), pointer :: cs !< Control structure for this module that
                                     !! will be deallocated in this subroutine

  if (.not. associated(cs)) return

  deallocate(cs)

end subroutine cvmix_conv_end

end module mom_cvmix_conv