SCM_CVmix_tests.F90

Go to the documentation of this file.

!> 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_error_handler, only : mom_error, fatal
use mom_file_parser, only : get_param, log_version, param_file_type
use mom_forcing_type, only : forcing, allocate_forcing_type
use mom_grid, only : ocean_grid_type
use mom_verticalgrid, only: verticalgrid_type
use mom_safe_alloc, only : safe_alloc_ptr
use mom_time_manager, only : time_type, operator(+), operator(/), get_time,&
                             time_type_to_real
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

!> 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 (m K s^{-1})
  real :: surf_evap !< (Constant) Evaporation rate (m/s)
  real :: max_sw !< maximum of diurnal sw radiation (m K s^{-1})
  real,public :: rho0 !< reference density copied for easy passing
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, param_file, just_read_params)
  real, dimension(NIMEM_,NJMEM_, NKMEM_), intent(out) :: T !< Potential temperature (degC)
  real, dimension(NIMEM_,NJMEM_, NKMEM_), intent(out) :: S !< Salinity (psu)
  real, dimension(NIMEM_,NJMEM_, NKMEM_), intent(in)  :: h !< Layer thickness (m or Pa)
  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 !< 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 :: eta(szk_(g)+1) ! The 1-d nominal positions of the interfaces.
  real :: UpperLayerTempMLD !< Upper layer Temp MLD thickness (m)
  real :: UpperLayerSaltMLD !< Upper layer Salt MLD thickness (m)
  real :: UpperLayerTemp !< Upper layer temperature (SST if thickness 0) (deg C)
  real :: UpperLayerSalt !< Upper layer salinity (SSS if thickness 0) (PPT)
  real :: LowerLayerTemp !< Temp at top of lower layer (deg C)
  real :: LowerLayerSalt !< Salt at top of lower layer (PPT)
  real :: LowerLayerdTdz !< Temp gradient in lower layer (deg C m^{-1})
  real :: LowerLayerdSdz !< Salt gradient in lower layer (PPT m^{-1})
  real :: zC, DZ
  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, 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, 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.00, 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.00, 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
    eta(1) = 0. ! Reference to surface
    do k=1,nz
      eta(k+1) = eta(k) - h(i,j,k)*gv%H_to_m ! Interface below layer (in m)
      zc = 0.5*( eta(k) + eta(k+1) )        ! Z of middle of layer (in m)
      dz = min(0., zc+upperlayertempmld)
      if (dz.ge.0.0) then ! in Layer 1
        t(i,j,k) = upperlayertemp
      else ! in Layer 2
        t(i,j,k) = lowerlayertemp + lowerlayerdtdz * dz
      endif
      dz = min(0., zc+upperlayersaltmld)
      if (dz.ge.0.0) then ! in Layer 1
        s(i,j,k) = upperlayersalt
      else ! in Layer 2
        s(i,j,k) = lowerlayersalt + lowerlayerdsdz * dz
      endif
    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       !< 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"

  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', 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', 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', 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', 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', fail_if_missing=.true.)
  endif

end subroutine scm_cvmix_tests_surface_forcing_init

subroutine scm_cvmix_tests_wind_forcing(state, fluxes, day, G, CS)
  type(surface),                    intent(in)    :: state  !< Surface state structure
  type(forcing),                    intent(inout) :: fluxes !< Surface fluxes structure
  type(time_type),                  intent(in)    :: day    !< Time in days
  type(ocean_grid_type),            intent(inout) :: G      !< Grid structure
  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

  ! Allocate the forcing arrays, if necessary.
  call allocate_forcing_type(g, fluxes, stress=.true., ustar=.true. )

  do j=js,je ; do i=isq,ieq
    fluxes%taux(i,j) = cs%tau_x
  enddo ; enddo
  do j=jsq,jeq ; do i=is,ie
    fluxes%tauy(i,j) = cs%tau_y
  enddo ; enddo
  mag_tau = sqrt(cs%tau_x*cs%tau_x + cs%tau_y*cs%tau_y)
  if (associated(fluxes%ustar)) then ; do j=js,je ; do i=is,ie
    fluxes%ustar(i,j) = sqrt(  mag_tau / cs%Rho0 )
  enddo ; enddo ; endif

end subroutine scm_cvmix_tests_wind_forcing


subroutine scm_cvmix_tests_buoyancy_forcing(state, fluxes, day, G, CS)
  type(surface),                    intent(in)    :: state  !< Surface state structure
  type(forcing),                    intent(inout) :: fluxes !< Surface fluxes structure
  type(time_type),                  intent(in)    :: day    !< Time in days (seconds?)
  type(ocean_grid_type),            intent(inout) :: G      !< Grid structure
  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

  ! Allocate the forcing arrays, if necessary.
  call allocate_forcing_type(g, fluxes, water=.true., heat=.true. )

  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
    ! This therefore must be converted to mass flux
    ! by multiplying by density
      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 K/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.5))) * cs%RHO0 * fluxes%C_p
    enddo; enddo
  endif

end subroutine scm_cvmix_tests_buoyancy_forcing

end module scm_cvmix_tests