Add a legacy version of diabatic_driver + updates in KPP and tidal mixing

src/core/MOM.F90

@@ -52,6 +52,7 @@ module MOM
 use MOM_ALE,                   only : ALE_getCoordinate, ALE_getCoordinateUnits, ALE_writeCoordinateFile
 use MOM_ALE,                   only : ALE_updateVerticalGridType, ALE_remap_init_conds, ALE_register_diags
 use MOM_boundary_update,       only : call_OBC_register, OBC_register_end, update_OBC_CS
+use MOM_legacy_diabatic_driver,only : legacy_diabatic
 use MOM_diabatic_driver,       only : diabatic, diabatic_driver_init, diabatic_CS
 use MOM_diabatic_driver,       only : adiabatic, adiabatic_driver_init, diabatic_driver_end
 use MOM_diagnostics,           only : calculate_diagnostic_fields, MOM_diagnostics_init
@@ -197,6 +198,9 @@ module MOM
 
   logical :: adiabatic               !< If true, there are no diapycnal mass fluxes, and no calls
                                      !! to routines to calculate or apply diapycnal fluxes.
+  logical :: use_legacy_diabatic_driver!< If true (default), use the a legacy version of the
+                                       !! diabatic subroutine. This is temporary and is needed
+                                       !! to avoid change in answers.
   logical :: diabatic_first          !< If true, apply diabatic and thermodynamic
                                      !! processes before time stepping the dynamics.
   logical :: use_ALE_algorithm       !< If true, use the ALE algorithm rather than layered
@@ -1151,8 +1155,14 @@ subroutine step_MOM_thermo(CS, G, GV, u, v, h, tv, fluxes, dtdia, &
     endif
 
     call cpu_clock_begin(id_clock_diabatic)
-    call diabatic(u, v, h, tv, CS%Hml, fluxes, CS%visc, CS%ADp, CS%CDp, &
-                  dtdia, Time_end_thermo, G, GV, CS%diabatic_CSp, Waves=Waves)
+    if (CS%use_legacy_diabatic_driver) then
+      ! the following subroutine is legacy and will be deleted in the near future.
+      call legacy_diabatic(u, v, h, tv, CS%Hml, fluxes, CS%visc, CS%ADp, CS%CDp, &
+                          dtdia, Time_end_thermo, G, GV, CS%diabatic_CSp, Waves=Waves)
+    else
+      call diabatic(u, v, h, tv, CS%Hml, fluxes, CS%visc, CS%ADp, CS%CDp, &
+                    dtdia, Time_end_thermo, G, GV, CS%diabatic_CSp, Waves=Waves)
+    endif
     fluxes%fluxes_used = .true.
     call cpu_clock_end(id_clock_diabatic)
 
@@ -1627,6 +1637,10 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
                  "true. This assumes that KD = KDML = 0.0 and that \n"//&
                  "there is no buoyancy forcing, but makes the model \n"//&
                  "faster by eliminating subroutine calls.", default=.false.)
+  call get_param(param_file, "MOM", "USE_LEGACY_DIABATIC_DRIVER", CS%use_legacy_diabatic_driver, &
+                 "If true, use the a legacy version of the diabatic subroutine. \n"//&
+                 "This is temporary and is needed avoid change in answers.", &
+                 default=.true.)
   call get_param(param_file, "MOM", "DO_DYNAMICS", CS%do_dynamics, &
                  "If False, skips the dynamics calls that update u & v, as well as \n"//&
                  "the gravity wave adjustment to h. This is a fragile feature and \n"//&
@@ -2364,6 +2378,9 @@ subroutine initialize_MOM(Time, Time_init, param_file, dirs, CS, restart_CSp, &
   if (associated(CS%visc%Kv_shear)) &
     call pass_var(CS%visc%Kv_shear, G%Domain, To_All+Omit_Corners, halo=1)
 
+  if (associated(CS%visc%Kv_slow)) &
+    call pass_var(CS%visc%Kv_slow, G%Domain, To_All+Omit_Corners, halo=1)
+
   call cpu_clock_end(id_clock_pass_init)
 
   call register_obsolete_diagnostics(param_file, CS%diag)

src/core/MOM_variables.F90

@@ -233,6 +233,9 @@ module MOM_variables
                          !! convection etc).
     TKE_turb => NULL()  !< The turbulent kinetic energy per unit mass defined
                         !! at the interfaces between each layer, in m2 s-2.
+    logical :: add_Kv_slow !< If True, adds Kv_slow when calculating the
+                        !! 'coupling coefficient' (a[k]) at the interfaces.
+                        !! This is done in find_coupling_coef.
 end type vertvisc_type
 
 !> The BT_cont_type structure contains information about the summed layer

src/parameterizations/vertical/MOM_CVMix_conv.F90

@@ -212,6 +212,7 @@ subroutine calculate_CVMix_conv(h, tv, G, GV, CS, hbl)
         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))
 
       call CVMix_coeffs_conv(Mdiff_out=CS%kv_conv(i,j,:), &

src/parameterizations/vertical/MOM_CVMix_ddiff.F90

@@ -0,0 +1,301 @@
+!> 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_variables,      only : thermo_var_ptrs
+use MOM_error_handler,  only : MOM_error, is_root_pe, FATAL, WARNING, NOTE
+use MOM_file_parser,    only : openParameterBlock, closeParameterBlock
+use MOM_debugging,      only : hchksum
+use MOM_grid,           only : ocean_grid_type
+use MOM_verticalGrid,   only : verticalGrid_type
+use MOM_file_parser,    only : get_param, log_version, param_file_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
+
+  ! 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 (m^2/s)
+  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 (m^2/s)
+  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
+
+  ! Daignostic handles and pointers
+  type(diag_ctrl), pointer :: diag => NULL()
+  integer :: id_KT_extra = -1, id_KS_extra = -1, id_R_rho = -1
+
+  ! Diagnostics arrays
+  real, allocatable, dimension(:,:,:) :: KT_extra  !< double diffusion diffusivity for temp (m2/s)
+  real, allocatable, dimension(:,:,:) :: KS_extra  !< double diffusion diffusivity for salt (m2/s)
+  real, allocatable, dimension(:,:,:) :: R_rho     !< double-diffusion density ratio (nondim)
+
+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, 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(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 log_version(param_file, mdl, version, &
+    "Parameterization of mixing due to double diffusion processes via CVMix")
+  call get_param(param_file, mdl, "USE_CVMIX_DDIFF", CVMix_ddiff_init, &
+                 "If true, turns on double diffusive processes via CVMix.  \n"// &
+                 "Note that double diffusive processes on viscosity are ignored \n"// &
+                 "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 \n"// &
+                 "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)
+
+  ! Register diagnostics
+  CS%diag => diag
+
+  CS%id_KT_extra = register_diag_field('ocean_model','KT_extra',diag%axesTi,Time, &
+         'Double-diffusive diffusivity for temperature', 'm2 s-1')
+
+  CS%id_KS_extra = register_diag_field('ocean_model','KS_extra',diag%axesTi,Time, &
+         'Double-diffusive diffusivity for salinity', 'm2 s-1')
+
+  CS%id_R_rho = register_diag_field('ocean_model','R_rho',diag%axesTi,Time, &
+         'Double-diffusion density ratio', 'nondim')
+  if (CS%id_R_rho > 0) &
+     allocate(CS%R_rho( SZI_(G), SZJ_(G), SZK_(G)+1)); CS%R_rho(:,:,:) = 0.0
+
+  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, j, Kd_T, Kd_S, CS)
+
+  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_(G)),   intent(in)  :: h    !< Layer thickness, in m or kg m-2.
+  type(thermo_var_ptrs),                      intent(in)  :: tv   !< Thermodynamics structure.
+  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: Kd_T !< Interface double diffusion diapycnal
+                                                                  !! diffusivity for temp (m2/sec).
+  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(out) :: Kd_S !< Interface double diffusion diapycnal
+                                                                  !! diffusivity for salt (m2/sec).
+  type(CVMix_ddiff_cs),                           pointer :: CS   !< The control structure returned
+                                                                  !! by a previous call to CVMix_ddiff_init.
+  integer,                                   intent(in)   :: j    !< Meridional grid indice.
+!  real, dimension(:,:),                 optional, pointer :: hbl  !< Depth of ocean boundary layer (m)
+
+  ! local variables
+  real, dimension(SZK_(G)) :: &
+    cellHeight, &  !< Height of cell centers (m)
+    dRho_dT,    &  !< partial derivatives of density wrt temp (kg m-3 degC-1)
+    dRho_dS,    &  !< partial derivatives of density wrt saln (kg m-3 PPT-1)
+    pres_int,   &  !< pressure at each interface (Pa)
+    temp_int,   &  !< temp and at interfaces (degC)
+    salt_int,   &  !< salt at at interfaces
+    alpha_dT,   &  !< alpha*dT across interfaces
+    beta_dS,    &  !< beta*dS across interfaces
+    dT,         &  !< temp. difference between adjacent layers (degC)
+    dS             !< salt difference between adjacent layers
+
+  real, dimension(SZK_(G)+1) :: iFaceHeight !< Height of interfaces (m)
+  integer :: kOBL                        !< level of OBL extent
+  real :: pref, g_o_rho0, rhok, rhokm1, dz, 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
+
+  ! set Kd_T and Kd_S to zero to avoid passing values from previous call
+  Kd_T(:,j,:) = 0.0; Kd_S(:,j,:) = 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
+
+    pRef = 0.
+    pres_int(1) = pRef
+    ! 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) = pRef + GV%H_to_Pa * 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))
+      pRef = pRef + GV%H_to_Pa * h(i,j,k-1)
+    enddo ! k-loop finishes
+
+    call calculate_density_derivs(temp_int(:), salt_int(:), pres_int(:), drho_dT(:), drho_dS(:), 1, G%ke, 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*drho_dT(k) * dT(k)
+      beta_dS(k)  = drho_dS(k) * dS(k)
+    enddo
+
+    if (CS%id_R_rho > 0.0)  then
+      do k=1,G%ke
+        CS%R_rho(i,j,k) = alpha_dT(k)/beta_dS(k)
+        ! avoid NaN's
+        if(CS%R_rho(i,j,k) /= CS%R_rho(i,j,k)) CS%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))
+
+    call CVMix_coeffs_ddiff(Tdiff_out=Kd_T(i,j,:), &
+                            Sdiff_out=Kd_S(i,j,:), &
+                            strat_param_num=alpha_dT(:), &
+                            strat_param_denom=beta_dS(:), &
+                            nlev=G%ke,    &
+                            max_nlev=G%ke)
+
+    ! 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
+
+  deallocate(CS)
+
+end subroutine CVMix_ddiff_end
+
+
+end module MOM_CVMix_ddiff