MOM_diffConvection.F90

Go to the documentation of this file.

module mom_diffconvection

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

implicit none ; private

#include "MOM_memory.h"

public :: diffconvection_init, diffconvection_calculate, diffconvection_end

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

  ! Parameters
  real    :: kd_convection ! The value of diffusivity to add at statically unstable interfaces (m2/s)
  logical :: debug         ! If true, turn on debugging
  logical :: passivemode   ! If true, make the motions but go nowhere

  ! Daignostic handles and pointers
  type(diag_ctrl), pointer :: diag => null()
  integer :: id_n2 = -1, id_kd_conv = -1

  ! Diagnostics arrays
  real, allocatable, dimension(:,:,:) :: n2      ! Brunt-Vaisala frequency (1/s2)
  real, allocatable, dimension(:,:,:) :: kd_conv ! Diffusivity added by convection (m2/s)

end type diffconvection_cs

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

contains

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

! 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(diffconvection_cs), pointer       :: CS        !< Control structure
! Local variables
#include "version_variable.h"
  character(len=40) :: mdl = 'MOM_diffConvection' ! This module's name.

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

! Read parameters
  call log_version(paramfile, mdl, version, &
            'This module implements enhanced diffusivity as a\n' // &
            'function of static stability, N^2.')
  call get_param(paramfile, mdl, "USE_CONVECTION", diffconvection_init, &
                 "If true, turns on the diffusive convection scheme that\n"// &
                 "increases diapycnal diffusivities at statically unstable\n"// &
                 "interfaces. Relevant parameters are contained in the\n"// &
                 "CONVECTION% parameter block.", &
                 default=.false.)

  call openparameterblock(paramfile,'CONVECTION')
  call get_param(paramfile, mdl, 'PASSIVE', cs%passiveMode,           &
                 'If True, puts KPP into a passive-diagnostic mode.', &
                 default=.false.)
  call get_param(paramfile, mdl, 'KD_CONV', cs%Kd_convection,                  &
                 'DIffusivity used in statically unstable regions of column.', &
                 units='m2/s', default=1.00)
  call closeparameterblock(paramfile)
  call get_param(paramfile, mdl, 'DEBUG', cs%debug, default=.false., do_not_log=.true.)

! Forego remainder of initialization if not using this scheme
  if (.not. diffconvection_init) return

! Register diagnostics
  cs%diag => diag
  cs%id_N2 = register_diag_field('ocean_model', 'Conv_N2', diag%axesTi, time, &
      'Square of Brunt-Vaisala frequency used by diffConvection module', '1/s2')
  if (cs%id_N2 > 0) allocate( cs%N2( szi_(g), szj_(g), szk_(g)+1 ) )
  cs%id_Kd_conv = register_diag_field('ocean_model', 'Conv_Kd', diag%axesTi, time, &
      'Additional diffusivity added by diffConvection module', 'm2/s')
  if (cs%id_Kd_conv > 0) allocate( cs%Kd_conv( szi_(g), szj_(g), szk_(g)+1 ) )

  if (cs%id_N2 > 0) cs%N2(:,:,:) = 0.
  if (cs%id_Kd_conv > 0) cs%Kd_conv(:,:,:) = 0.

end function diffconvection_init


subroutine diffconvection_calculate(CS, G, GV, h, Temp, Salt, EOS, Kd_int)
!< Calculates diffusivity and non-local transport for KPP parameterization

! Arguments
  type(diffconvection_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  !< Pot. temperature (degrees C)
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: Salt  !< Salinity (ppt)
  type(eos_type),                            pointer       :: EOS   !< Equation of state
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1), intent(inout) :: Kd_int !< (in) Vertical diffusivity on interfaces (m2/s)
                                                                 !! (out) Modified vertical diffusivity (m2/s)
! Local variables
  integer :: i, j, k
  real, dimension( G%ke+1 ) :: N2_1d ! Brunt-Vaisala frequency squared, at interfaces (1/s2)
  real, dimension( G%ke+1 ) :: Kd_1d ! Vertical diffusivity at interfaces (m2/s)
  real :: GoRho, pRef, rhoK, rhoKm1

  gorho = gv%g_Earth / gv%Rho0

  n2_1d( 1 ) = 0.
  n2_1d( g%ke+1 ) = 0.
  kd_1d( 1 ) = 0.
  kd_1d( g%ke+1 ) = 0.
  do j=g%jsc,g%jec ; do i=g%isc,g%iec
    ! This k-loop calculates external quantities independent of any iterations
    ! Start at bottom of top level
    pref = 0. ! Ignore atmospheric pressure
    do k = 2, g%ke
      ! Pressure at interface K is incremented by mass of level above
      pref = pref + gv%g_Earth * gv%Rho0 * h(i,j,k-1) * gv%H_to_m ! Boussinesq approximation!!!! ?????
      ! Compute Brunt-Vaisala frequency (static stability) on interfaces
      call calculate_density(temp(i,j,k),   salt(i,j,k),   pref, rhok,   eos)
      call calculate_density(temp(i,j,k-1), salt(i,j,k-1), pref, rhokm1, eos)
      n2_1d(k) = gorho * (rhok - rhokm1) / &
              (0.5*(h(i,j,k-1) + h(i,j,k)) + gv%H_subroundoff) ! Can be negative
      kd_1d(k) = 0.
      if (n2_1d(k) < 0.) kd_1d(k) = cs%Kd_convection
    enddo ! k

    if (.not. cs%passiveMode) kd_int(i,j,:) = kd_int(i,j,:) + kd_1d(:)

    if (cs%id_N2 > 0) cs%N2(i,j,:) = n2_1d(:)
    if (cs%id_Kd_conv > 0) cs%Kd_conv(i,j,:) = kd_1d(:)

  enddo ; enddo ! j

  if (cs%id_N2 > 0) call post_data(cs%id_N2, cs%N2, cs%diag)
  if (cs%id_Kd_conv > 0) call post_data(cs%id_Kd_conv, cs%Kd_conv, cs%diag)

end subroutine diffconvection_calculate


subroutine diffconvection_end(CS)
! Clear pointers, dealocate memory
  type(diffconvection_cs), pointer :: CS ! Control structure

  if (cs%id_N2 > 0) deallocate(cs%N2, cs%diag)
  if (cs%id_Kd_conv > 0) deallocate(cs%Kd_conv, cs%diag)
  deallocate(cs)
end subroutine diffconvection_end

end module mom_diffconvection