MOM_CVMix_ddiff.F90

!> Interface to CVMix double diffusion scheme.
module mom_cvmix_ddiff

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

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_derivs
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_debugging,      only : hchksum
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_ddiff,        only : cvmix_init_ddiff, cvmix_coeffs_ddiff
use cvmix_kpp,          only : cvmix_kpp_compute_kobl_depth
implicit none ; private

#include <MOM_memory.h>

public cvmix_ddiff_init, cvmix_ddiff_end, cvmix_ddiff_is_used, compute_ddiff_coeffs

!> Control structure including parameters for CVMix double diffusion.
type, public :: cvmix_ddiff_cs ; private

  ! Parameters
  real    :: strat_param_max !< maximum value for the stratification parameter [nondim]
  real    :: kappa_ddiff_s   !< leading coefficient in formula for salt-fingering regime
                             !! for salinity diffusion [m2 s-1]
  real    :: ddiff_exp1      !< interior exponent in salt-fingering regime formula [nondim]
  real    :: ddiff_exp2      !< exterior exponent in salt-fingering regime formula [nondim]
  real    :: mol_diff        !< molecular diffusivity [m2 s-1]
  real    :: kappa_ddiff_param1 !< exterior coefficient in diffusive convection regime [nondim]
  real    :: kappa_ddiff_param2 !< middle coefficient in diffusive convection regime [nondim]
  real    :: kappa_ddiff_param3 !< interior coefficient in diffusive convection regime [nondim]
  real    :: min_thickness      !< Minimum thickness allowed [m]
  character(len=4) :: diff_conv_type !< type of diffusive convection to use. Options are Marmorino &
                                !! Caldwell 1976 ("MC76"; default) and Kelley 1988, 1990 ("K90")
  logical :: debug              !< If true, turn on debugging

end type cvmix_ddiff_cs

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

contains

!> Initialized the CVMix double diffusion module.
logical function cvmix_ddiff_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_ddiff_cs),    pointer       :: cs         !< This module's control structure.

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

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

  ! Read parameters
  call get_param(param_file, mdl, "USE_CVMIX_DDIFF", cvmix_ddiff_init, default=.false., do_not_log=.true.)
  call log_version(param_file, mdl, version, &
           "Parameterization of mixing due to double diffusion processes via CVMix", &
           all_default=.not.cvmix_ddiff_init)
  call get_param(param_file, mdl, "USE_CVMIX_DDIFF", cvmix_ddiff_init, &
                 "If true, turns on double diffusive processes via CVMix. "//&
                 "Note that double diffusive processes on viscosity are ignored "//&
                 "in CVMix, see http://cvmix.github.io/ for justification.", &
                 default=.false.)

  if (.not. cvmix_ddiff_init) return

  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_DDIFF')

  call get_param(param_file, mdl, "STRAT_PARAM_MAX", cs%strat_param_max, &
                 "The maximum value for the double dissusion stratification parameter", &
                 units="nondim", default=2.55)

  call get_param(param_file, mdl, "KAPPA_DDIFF_S", cs%kappa_ddiff_s, &
                 "Leading coefficient in formula for salt-fingering regime "//&
                 "for salinity diffusion.", units="m2 s-1", default=1.0e-4)

  call get_param(param_file, mdl, "DDIFF_EXP1", cs%ddiff_exp1, &
                 "Interior exponent in salt-fingering regime formula.", &
                 units="nondim", default=1.0)

  call get_param(param_file, mdl, "DDIFF_EXP2", cs%ddiff_exp2, &
                 "Exterior exponent in salt-fingering regime formula.", &
                 units="nondim", default=3.0)

  call get_param(param_file, mdl, "KAPPA_DDIFF_PARAM1", cs%kappa_ddiff_param1, &
                "Exterior coefficient in diffusive convection regime.", &
                 units="nondim", default=0.909)

  call get_param(param_file, mdl, "KAPPA_DDIFF_PARAM2", cs%kappa_ddiff_param2, &
                "Middle coefficient in diffusive convection regime.", &
                 units="nondim", default=4.6)

  call get_param(param_file, mdl, "KAPPA_DDIFF_PARAM3", cs%kappa_ddiff_param3, &
                "Interior coefficient in diffusive convection regime.", &
                 units="nondim", default=-0.54)

  call get_param(param_file, mdl, "MOL_DIFF", cs%mol_diff, &
                 "Molecular diffusivity used in CVMix double diffusion.", &
                 units="m2 s-1", default=1.5e-6)

  call get_param(param_file, mdl, "DIFF_CONV_TYPE", cs%diff_conv_type, &
                 "type of diffusive convection to use. Options are Marmorino \n" //&
                 "and Caldwell 1976 (MC76) and Kelley 1988, 1990 (K90).", &
                 default="MC76")

  call closeparameterblock(param_file)

  call cvmix_init_ddiff(strat_param_max=cs%strat_param_max,          &
                        kappa_ddiff_s=cs%kappa_ddiff_s,           &
                        ddiff_exp1=cs%ddiff_exp1,                 &
                        ddiff_exp2=cs%ddiff_exp2,                 &
                        mol_diff=cs%mol_diff,                     &
                        kappa_ddiff_param1=cs%kappa_ddiff_param1, &
                        kappa_ddiff_param2=cs%kappa_ddiff_param2, &
                        kappa_ddiff_param3=cs%kappa_ddiff_param3, &
                        diff_conv_type=cs%diff_conv_type)

end function cvmix_ddiff_init

!> Subroutine for computing vertical diffusion coefficients for the
!! double diffusion mixing parameterization.
subroutine compute_ddiff_coeffs(h, tv, G, GV, US, j, Kd_T, Kd_S, CS, R_rho)

  type(ocean_grid_type),                      intent(in)    :: g    !< Grid structure.
  type(verticalgrid_type),                    intent(in)    :: gv   !< Vertical grid structure.
  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(unit_scale_type),                      intent(in)    :: us   !< A dimensional unit scaling type
  integer,                                    intent(in)    :: j    !< Meridional grid index to work on.
  ! Kd_T and Kd_S are intent inout because only one j-row is set here, but they are essentially outputs.
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: kd_t !< Interface double diffusion diapycnal
                                                                    !! diffusivity for temp [Z2 T-1 ~> m2 s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), intent(inout) :: kd_s !< Interface double diffusion diapycnal
                                                                    !! diffusivity for salt [Z2 T-1 ~> m2 s-1].
  type(cvmix_ddiff_cs),                       pointer       :: cs   !< The control structure returned
                                                                    !! by a previous call to CVMix_ddiff_init.
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)+1), &
                                    optional, intent(inout) :: r_rho !< The density ratios at interfaces [nondim].

  ! Local variables
  real, dimension(SZK_(GV)) :: &
    cellheight, &  !< Height of cell centers [m]
    drho_dt,    &  !< partial derivatives of density wrt temp [R degC-1 ~> kg m-3 degC-1]
    drho_ds,    &  !< partial derivatives of density wrt saln [R ppt-1 ~> kg m-3 ppt-1]
    pres_int,   &  !< pressure at each interface [R L2 T-2 ~> Pa]
    temp_int,   &  !< temp and at interfaces [degC]
    salt_int,   &  !< salt at at interfaces [ppt]
    alpha_dt,   &  !< alpha*dT across interfaces [kg m-3]
    beta_ds,    &  !< beta*dS across interfaces [kg m-3]
    dt,         &  !< temp. difference between adjacent layers [degC]
    ds             !< salt difference between adjacent layers [ppt]
  real, dimension(SZK_(GV)+1) :: &
    kd1_t,      &  !< Diapycanal diffusivity of temperature [m2 s-1].
    kd1_s          !< Diapycanal diffusivity of salinity [m2 s-1].

  real, dimension(SZK_(GV)+1) :: ifaceheight !< Height of interfaces [m]
  integer :: kobl                        !< level of OBL extent
  real :: dh, hcorr
  integer :: i, k

  ! initialize dummy variables
  pres_int(:) = 0.0; temp_int(:) = 0.0; salt_int(:) = 0.0
  alpha_dt(:) = 0.0; beta_ds(:) = 0.0; drho_dt(:) = 0.0
  drho_ds(:) = 0.0; dt(:) = 0.0; ds(:) = 0.0


  ! GMM, I am leaving some code commented below. We need to pass BLD to
  ! this soubroutine to avoid adding diffusivity above that. This needs
  ! to be done once we re-structure the order of the calls.
  !if (.not. associated(hbl)) then
  !  allocate(hbl(SZI_(G), SZJ_(G)));
  !  hbl(:,:) = 0.0
  !endif

  do i = g%isc, g%iec

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

    pres_int(1) = 0. ;  if (associated(tv%p_surf)) pres_int(1) = tv%p_surf(i,j)
    ! we don't have SST and SSS, so let's use values at top-most layer
    temp_int(1) = tv%T(i,j,1); salt_int(1) = tv%S(i,j,1)
    do k=2,g%ke
      ! pressure at interface
      pres_int(k) = pres_int(k-1) + (gv%g_Earth * gv%H_to_RZ) * h(i,j,k-1)
      ! temp and salt at interface
      ! for temp: (t1*h1 + t2*h2)/(h1+h2)
      temp_int(k) = (tv%T(i,j,k-1)*h(i,j,k-1) + tv%T(i,j,k)*h(i,j,k)) / (h(i,j,k-1)+h(i,j,k))
      salt_int(k) = (tv%S(i,j,k-1)*h(i,j,k-1) + tv%S(i,j,k)*h(i,j,k)) / (h(i,j,k-1)+h(i,j,k))
      ! dT and dS
      dt(k) = (tv%T(i,j,k-1)-tv%T(i,j,k))
      ds(k) = (tv%S(i,j,k-1)-tv%S(i,j,k))
    enddo ! k-loop finishes

    call calculate_density_derivs(temp_int, salt_int, pres_int, drho_dt, drho_ds, tv%eqn_of_state)

    ! The "-1.0" below is needed so that the following criteria is satisfied:
    ! if ((alpha_dT > beta_dS) .and. (beta_dS > 0.0)) then "salt finger"
    ! if ((alpha_dT < 0.) .and. (beta_dS < 0.) .and. (alpha_dT > beta_dS)) then "diffusive convection"
    do k=1,g%ke
      alpha_dt(k) = -1.0*us%R_to_kg_m3*drho_dt(k) * dt(k)
      beta_ds(k)  = us%R_to_kg_m3*drho_ds(k) * ds(k)
    enddo

    if (present(r_rho))  then
      do k=1,g%ke
        ! Set R_rho using Adcroft's rule of reciprocals.
        r_rho(i,j,k) = 0.0 ; if (abs(beta_ds(k)) > 0.0) r_rho(i,j,k) = alpha_dt(k) / beta_ds(k)
        ! avoid NaN's again for safety, perhaps unnecessarily.
        if (r_rho(i,j,k) /= r_rho(i,j,k)) r_rho(i,j,k) = 0.0
      enddo
    endif

    ifaceheight(1) = 0.0 ! BBL is all relative to the surface
    hcorr = 0.0
    ! compute heights at cell center and interfaces
    do k=1,g%ke
      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
    enddo

    ! gets index of the level and interface above hbl
    !kOBL = CVmix_kpp_compute_kOBL_depth(iFaceHeight, cellHeight,hbl(i,j))

    kd1_t(:) = 0.0 ; kd1_s(:) = 0.0
    call cvmix_coeffs_ddiff(tdiff_out=kd1_t(:), &
                            sdiff_out=kd1_s(:), &
                            strat_param_num=alpha_dt(:), &
                            strat_param_denom=beta_ds(:), &
                            nlev=g%ke,    &
                            max_nlev=g%ke)
    do k=1,g%ke+1
      kd_t(i,j,k) = us%m2_s_to_Z2_T * kd1_t(k)
      kd_s(i,j,k) = us%m2_s_to_Z2_T * kd1_s(k)
    enddo

    ! Do not apply mixing due to convection within the boundary layer
    !do k=1,kOBL
    !  Kd_T(i,j,k) = 0.0
    !  Kd_S(i,j,k) = 0.0
    !enddo

  enddo ! i-loop

end subroutine compute_ddiff_coeffs

!> Reads the parameter "USE_CVMIX_DDIFF" 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_ddiff_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_DDIFF", cvmix_ddiff_is_used, &
                 default=.false., do_not_log = .true.)

end function cvmix_ddiff_is_used

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

  deallocate(cs)

end subroutine cvmix_ddiff_end

end module mom_cvmix_ddiff