SCM_CVMix_tests.F90

!> Initial conditions and forcing for the single column model (SCM) CVMix test set.
module scm_cvmix_tests

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

use mom_domains,       only : pass_var, pass_vector, to_all
use mom_error_handler, only : mom_error, fatal
use mom_file_parser,   only : get_param, log_version, param_file_type
use mom_forcing_type,  only : forcing, mech_forcing
use mom_grid,          only : ocean_grid_type
use mom_verticalgrid,  only : verticalgrid_type
use mom_safe_alloc,    only : safe_alloc_ptr
use mom_unit_scaling,  only : unit_scale_type
use mom_time_manager,  only : time_type, operator(+), operator(/), time_type_to_real
use mom_unit_scaling,  only : unit_scale_type
use mom_variables,     only : thermo_var_ptrs, surface

implicit none ; private

#include <MOM_memory.h>

public scm_cvmix_tests_ts_init
public scm_cvmix_tests_surface_forcing_init
public scm_cvmix_tests_wind_forcing
public scm_cvmix_tests_buoyancy_forcing
public scm_cvmix_tests_cs

! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]".  If the units
! vary with the Boussinesq approximation, the Boussinesq variant is given first.

!> Container for surface forcing parameters
type scm_cvmix_tests_cs ; private
  logical :: usewindstress  !< True to use wind stress
  logical :: useheatflux    !< True to use heat flux
  logical :: useevaporation !< True to use evaporation
  logical :: usediurnalsw   !< True to use diurnal sw radiation
  real :: tau_x !< (Constant) Wind stress, X [Pa]
  real :: tau_y !< (Constant) Wind stress, Y [Pa]
  real :: surf_hf !< (Constant) Heat flux [degC Z T-1 ~> m degC s-1]
  real :: surf_evap !< (Constant) Evaporation rate [Z T-1 ~> m s-1]
  real :: max_sw !< maximum of diurnal sw radiation [degC Z T-1 ~> degC m s-1]
  real :: rho0 !< reference density [R ~> kg m-3]
end type

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

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

contains

!> Initializes temperature and salinity for the SCM CVMix test example
subroutine scm_cvmix_tests_ts_init(T, S, h, G, GV, US, param_file, just_read_params)
  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(out) :: t  !< Potential temperature [degC]
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(out) :: s  !< Salinity [psu]
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)  :: h  !< Layer thickness [H ~> m or kg m-2]
  type(unit_scale_type),                     intent(in)  :: us !< A dimensional unit scaling type
  type(param_file_type),                     intent(in)  :: param_file !< Input parameter structure
  logical,                         optional, intent(in)  :: just_read_params !< If present and true, this call
                                                               !! will only read parameters without changing h.
  ! Local variables
  real :: upperlayertempmld !< Upper layer Temp MLD thickness [Z ~> m].
  real :: upperlayersaltmld !< Upper layer Salt MLD thickness [Z ~> m].
  real :: upperlayertemp !< Upper layer temperature (SST if thickness 0) [degC]
  real :: upperlayersalt !< Upper layer salinity (SSS if thickness 0) [ppt]
  real :: lowerlayertemp !< Temp at top of lower layer [degC]
  real :: lowerlayersalt !< Salt at top of lower layer [ppt]
  real :: lowerlayerdtdz !< Temp gradient in lower layer [degC / Z ~> degC m-1].
  real :: lowerlayerdsdz !< Salt gradient in lower layer [ppt / Z ~> ppt m-1].
  real :: lowerlayermintemp !< Minimum temperature in lower layer [degC]
  real :: zc, dz, top, bottom ! Depths and thicknesses [Z ~> m].
  logical :: just_read    ! If true, just read parameters but set nothing.
  integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed


  just_read = .false. ; if (present(just_read_params)) just_read = just_read_params

  if (.not.just_read) call log_version(param_file, mdl, version)
  call get_param(param_file, mdl, "SCM_TEMP_MLD", upperlayertempmld, &
                 'Initial temp mixed layer depth', &
                 units='m', default=0.0, scale=us%m_to_Z, do_not_log=just_read)
  call get_param(param_file, mdl, "SCM_SALT_MLD", upperlayersaltmld, &
                 'Initial salt mixed layer depth', &
                 units='m', default=0.0, scale=us%m_to_Z, do_not_log=just_read)
  call get_param(param_file, mdl, "SCM_L1_SALT", upperlayersalt, &
                 'Layer 2 surface salinity', units='1e-3', default=35.0, do_not_log=just_read)
  call get_param(param_file, mdl, "SCM_L1_TEMP", upperlayertemp, &
                 'Layer 1 surface temperature', units='C', default=20.0, do_not_log=just_read)
  call get_param(param_file, mdl, "SCM_L2_SALT", lowerlayersalt, &
                 'Layer 2 surface salinity', units='1e-3', default=35.0, do_not_log=just_read)
  call get_param(param_file, mdl, "SCM_L2_TEMP", lowerlayertemp, &
                 'Layer 2 surface temperature', units='C', default=20.0, do_not_log=just_read)
  call get_param(param_file, mdl, "SCM_L2_DTDZ", lowerlayerdtdz,     &
                 'Initial temperature stratification in layer 2', &
                 units='C/m', default=0.0, scale=us%Z_to_m, do_not_log=just_read)
  call get_param(param_file, mdl, "SCM_L2_DSDZ", lowerlayerdsdz,  &
                 'Initial salinity stratification in layer 2', &
                 units='PPT/m', default=0.0, scale=us%Z_to_m, do_not_log=just_read)
  call get_param(param_file, mdl, "SCM_L2_MINTEMP",lowerlayermintemp, &
                 'Layer 2 minimum temperature', units='C', default=4.0, do_not_log=just_read)

  if (just_read) return ! All run-time parameters have been read, so return.

  do j=js,je ; do i=is,ie
    top = 0. ! Reference to surface
    bottom = 0.
    do k=1,nz
      bottom = bottom - h(i,j,k)*gv%H_to_Z ! Interface below layer [Z ~> m]
      zc = 0.5*( top + bottom )        ! Z of middle of layer [Z ~> m]
      dz = min(0., zc + upperlayertempmld)
      t(i,j,k) = max(lowerlayermintemp,lowerlayertemp + lowerlayerdtdz * dz)
      dz = min(0., zc + upperlayersaltmld)
      s(i,j,k) = lowerlayersalt + lowerlayerdsdz * dz
      top = bottom
    enddo ! k
  enddo ; enddo

end subroutine scm_cvmix_tests_ts_init

!> Initializes surface forcing for the CVMix test case suite
subroutine scm_cvmix_tests_surface_forcing_init(Time, G, param_file, CS)
  type(time_type),          intent(in) :: time       !< Model time
  type(ocean_grid_type),    intent(in) :: g          !< Grid structure
  type(param_file_type),    intent(in) :: param_file !< Input parameter structure
  type(scm_cvmix_tests_cs), pointer    :: cs         !< Parameter container


  ! This include declares and sets the variable "version".
# include "version_variable.h"
  type(unit_scale_type), pointer :: us => null() !< A dimensional unit scaling type

  us => g%US

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

  ! Read all relevant parameters and write them to the model log.
  call log_version(param_file, mdl, version, "")
  call get_param(param_file, mdl, "SCM_USE_WIND_STRESS",              &
                 cs%UseWindStress, "Wind Stress switch "//            &
                 "used in the SCM CVMix surface forcing.",            &
                 units='', default=.false.)
  call get_param(param_file, mdl, "SCM_USE_HEAT_FLUX",                &
                 cs%UseHeatFlux, "Heat flux switch "//                &
                 "used in the SCM CVMix test surface forcing.",       &
                 units='', default=.false.)
  call get_param(param_file, mdl, "SCM_USE_EVAPORATION",              &
                 cs%UseEvaporation, "Evaporation switch "//           &
                 "used in the SCM CVMix test surface forcing.",       &
                 units='', default=.false.)
  call get_param(param_file, mdl, "SCM_USE_DIURNAL_SW",               &
                 cs%UseDiurnalSW, "Diurnal sw radation switch "//     &
                 "used in the SCM CVMix test surface forcing.",       &
                 units='', default=.false.)
  if (cs%UseWindStress) then
    call get_param(param_file, mdl, "SCM_TAU_X",                      &
                 cs%tau_x, "Constant X-dir wind stress "//            &
                 "used in the SCM CVMix test surface forcing.",       &
                 units='N/m2', scale=us%kg_m2s_to_RZ_T*us%m_s_to_L_T, fail_if_missing=.true.)
    call get_param(param_file, mdl, "SCM_TAU_Y",                      &
                 cs%tau_y, "Constant y-dir wind stress "//            &
                 "used in the SCM CVMix test surface forcing.",       &
                 units='N/m2', scale=us%kg_m2s_to_RZ_T*us%m_s_to_L_T, fail_if_missing=.true.)
  endif
  if (cs%UseHeatFlux) then
    call get_param(param_file, mdl, "SCM_HEAT_FLUX",                  &
                 cs%surf_HF, "Constant surface heat flux "//          &
                 "used in the SCM CVMix test surface forcing.",       &
                 units='m K/s', scale=us%m_to_Z*us%T_to_s, fail_if_missing=.true.)
  endif
  if (cs%UseEvaporation) then
    call get_param(param_file, mdl, "SCM_EVAPORATION",                &
                 cs%surf_evap, "Constant surface evaporation "//      &
                 "used in the SCM CVMix test surface forcing.",       &
                 units='m/s', scale=us%m_to_Z*us%T_to_s, fail_if_missing=.true.)
  endif
  if (cs%UseDiurnalSW) then
    call get_param(param_file, mdl, "SCM_DIURNAL_SW_MAX",             &
                 cs%Max_sw, "Maximum diurnal sw radiation "//         &
                 "used in the SCM CVMix test surface forcing.",       &
                 units='m K/s', scale=us%m_to_Z*us%T_to_s, fail_if_missing=.true.)
  endif
  call get_param(param_file, mdl, "RHO_0", cs%Rho0, &
                 "The mean ocean density used with BOUSSINESQ true to "//&
                 "calculate accelerations and the mass for conservation "//&
                 "properties, or with BOUSSINSEQ false to convert some "//&
                 "parameters from vertical units of m to kg m-2.", &
                 units="kg m-3", default=1035.0, scale=us%kg_m3_to_R)

end subroutine scm_cvmix_tests_surface_forcing_init

subroutine scm_cvmix_tests_wind_forcing(sfc_state, forces, day, G, US, CS)
  type(surface),            intent(in)    :: sfc_state  !< Surface state structure
  type(mech_forcing),       intent(inout) :: forces !< A structure with the driving mechanical forces
  type(time_type),          intent(in)    :: day    !< Time in days
  type(ocean_grid_type),    intent(inout) :: g      !< Grid structure
  type(unit_scale_type),    intent(in)    :: us     !< A dimensional unit scaling type
  type(scm_cvmix_tests_cs), pointer       :: cs     !< Container for SCM parameters
  ! Local variables
  integer :: i, j, is, ie, js, je, isq, ieq, jsq, jeq
  integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
  real    :: mag_tau
  ! Bounds for loops and memory allocation
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

  do j=js,je ; do i=isq,ieq
    forces%taux(i,j) = cs%tau_x
  enddo ; enddo
  do j=jsq,jeq ; do i=is,ie
    forces%tauy(i,j) = cs%tau_y
  enddo ; enddo
  call pass_vector(forces%taux, forces%tauy, g%Domain, to_all)

  mag_tau = sqrt(cs%tau_x*cs%tau_x + cs%tau_y*cs%tau_y)
  if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie
    forces%ustar(i,j) = sqrt( us%L_to_Z * mag_tau / (cs%Rho0) )
  enddo ; enddo ; endif

end subroutine scm_cvmix_tests_wind_forcing


subroutine scm_cvmix_tests_buoyancy_forcing(sfc_state, fluxes, day, G, US, CS)
  type(surface),            intent(in)    :: sfc_state  !< Surface state structure
  type(forcing),            intent(inout) :: fluxes !< Surface fluxes structure
  type(time_type),          intent(in)    :: day    !< Current model time
  type(ocean_grid_type),    intent(inout) :: g      !< Grid structure
  type(unit_scale_type),    intent(in)    :: us     !< A dimensional unit scaling type
  type(scm_cvmix_tests_cs), pointer       :: cs     !< Container for SCM parameters

  ! Local variables
  integer :: i, j, is, ie, js, je, isq, ieq, jsq, jeq
  integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
  real :: pi

  pi = 4.0*atan(1.0)

  ! Bounds for loops and memory allocation
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

  if (cs%UseHeatFlux) then
    ! Note CVMix test inputs give Heat flux in [m K/s]
    ! therefore must convert to W/m2 by multiplying
    ! by Rho0*Cp
    do j=jsq,jeq ; do i=is,ie
      fluxes%sens(i,j) = cs%surf_HF * cs%Rho0 * fluxes%C_p
    enddo ; enddo
  endif

  if (cs%UseEvaporation) then
    do j=jsq,jeq ; do i=is,ie
    ! Note CVMix test inputs give evaporation in [m s-1]
    ! This therefore must be converted to mass flux in [R Z T-1 ~> kg m-2 s-1]
    ! by multiplying by density and some unit conversion factors.
      fluxes%evap(i,j) = cs%surf_evap * cs%Rho0
    enddo ; enddo
  endif

  if (cs%UseDiurnalSW) then
    do j=jsq,jeq ; do i=is,ie
    ! Note CVMix test inputs give max sw rad in [m degC/s]
    ! therefore must convert to W/m2 by multiplying by Rho0*Cp
    ! Note diurnal cycle peaks at Noon.
      fluxes%sw(i,j) = cs%Max_sw *  max(0.0, cos(2*pi*(time_type_to_real(day)/86400.0 - 0.5))) * cs%RHO0 * fluxes%C_p
    enddo ; enddo
  endif

end subroutine scm_cvmix_tests_buoyancy_forcing

end module scm_cvmix_tests