Neverland_surface_forcing.F90

!> Wind and buoyancy forcing for the Neverland configurations
module neverland_surface_forcing
 
! This file is part of MOM6. See LICENSE.md for the license.
 
use mom_diag_mediator, only : post_data, query_averaging_enabled
use mom_diag_mediator, only : register_diag_field, diag_ctrl, safe_alloc_ptr
use mom_domains, only : pass_var, pass_vector, agrid
use mom_error_handler, only : mom_error, fatal, warning, is_root_pe
use mom_file_parser, only : get_param, log_version, param_file_type
use mom_forcing_type, only : forcing, mech_forcing
use mom_forcing_type, only : allocate_forcing_type, allocate_mech_forcing
use mom_grid, only : ocean_grid_type
use mom_io, only : file_exists, read_data, slasher
use mom_time_manager, only : time_type, operator(+), operator(/)
use mom_unit_scaling, only : unit_scale_type
use mom_variables, only : surface
 
implicit none ; private
 
public neverland_wind_forcing
public neverland_buoyancy_forcing
public neverland_surface_forcing_init
 
!> This control structure should be used to store any run-time variables
!! associated with the Neverland forcing.
!!
!! It can be readily modified for a specific case, and because it is private there
!! will be no changes needed in other code (although they will have to be recompiled).
type, public :: neverland_surface_forcing_cs ; private
 
  logical :: use_temperature !< If true, use temperature and salinity.
  logical :: restorebuoy     !< If true, use restoring surface buoyancy forcing.
  real :: rho0               !< The density used in the Boussinesq
                             !! approximation [kg m-3].
  real :: g_earth            !< The gravitational acceleration [L2 Z-1 T-2 ~> m s-2].
  real :: flux_const         !<  The restoring rate at the surface [m s-1].
  real, dimension(:,:), pointer :: &
    buoy_restore(:,:) => null() !< The pattern to restore buoyancy to.
  character(len=200) :: inputdir !< The directory where NetCDF input files are.
  type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the
                                   !! timing of diagnostic output.
  logical :: first_call = .true. !< True until Neverland_buoyancy_forcing has been called
end type neverland_surface_forcing_cs
 
contains
 
!> Sets the surface wind stresses, forces%taux and forces%tauy for the
!! Neverland forcing configuration.
subroutine neverland_wind_forcing(sfc_state, forces, day, G, US, CS)
  type(surface),                 intent(inout) :: sfc_state !< A structure containing fields that
                                                         !! describe the surface state of the ocean.
  type(mech_forcing),            intent(inout) :: forces !< A structure with the driving mechanical forces
  type(time_type),               intent(in)    :: day    !< Time used for determining the fluxes.
  type(ocean_grid_type),         intent(inout) :: g      !< Grid structure.
  type(unit_scale_type),         intent(in)    :: us     !< A dimensional unit scaling type
  type(neverland_surface_forcing_cs), pointer  :: cs     !< Control structure for this module.
 
  ! Local variables
  integer :: i, j, is, ie, js, je, isq, ieq, jsq, jeq
  integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
  real :: x, y
  real :: pi
  real :: tau_max, off
 
  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_mech_forcing(g, forces, stress=.true.)
 
  !  Set the surface wind stresses, in units of Pa.  A positive taux
  !  accelerates the ocean to the (pseudo-)east.
 
  !  The i-loop extends to is-1 so that taux can be used later in the
  ! calculation of ustar - otherwise the lower bound would be Isq.
  pi = 4.0*atan(1.0)
  forces%taux(:,:) = 0.0
  tau_max = 0.2
  off = 0.02
  do j=js,je ; do i=is-1,ieq
!    x = (G%geoLonT(i,j)-G%west_lon)/G%len_lon
    y = (g%geoLatT(i,j)-g%south_lat)/g%len_lat
!    forces%taux(I,j) = G%mask2dCu(I,j) * 0.0
 
    if (y <= 0.29) then
      forces%taux(i,j) = forces%taux(i,j) + tau_max * ( (1/0.29)*y - ( 1/(2*pi) )*sin( (2*pi*y) / 0.29 ) )
    endif
    if ((y > 0.29) .and. (y <= (0.8-off))) then
      forces%taux(i,j) = forces%taux(i,j) + tau_max *(0.35+0.65*cos(pi*(y-0.29)/(0.51-off))  )
    endif
    if ((y > (0.8-off)) .and. (y <= (1-off))) then
      forces%taux(i,j) = forces%taux(i,j) + tau_max *( 1.5*( (y-1+off) - (0.1/pi)*sin(10.0*pi*(y-0.8+off)) ) )
    endif
  enddo ; enddo
 
  do j=js-1,jeq ; do i=is,ie
    forces%tauy(i,j) = g%mask2dCv(i,j) * 0.0  ! Change this to the desired expression.
  enddo ; enddo
 
  !    Set the surface friction velocity, in units of m s-1.  ustar
  !  is always positive.
! if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie
!   !  This expression can be changed if desired, but need not be.
!   forces%ustar(i,j) = US%m_to_Z*US%T_to_s * G%mask2dT(i,j) * sqrt(CS%gust_const/CS%Rho0 + &
!      sqrt(0.5*(forces%taux(I-1,j)**2 + forces%taux(I,j)**2) + &
!           0.5*(forces%tauy(i,J-1)**2 + forces%tauy(i,J)**2))/CS%Rho0)
! enddo ; enddo ; endif
 
end subroutine neverland_wind_forcing
 
!> Returns the value of a cosine-bell function evaluated at x/L
real function cosbell(x,L)
 
  real , intent(in) :: x       !< non-dimensional position
  real , intent(in) :: l       !< non-dimensional width
  real              :: pi      !< 3.1415926... calculated as 4*atan(1)
 
  pi      = 4.0*atan(1.0)
  cosbell = 0.5 * (1 + cos(pi*min(abs(x/l),1.0)))
end function cosbell
 
!> Returns the value of a sin-spike function evaluated at x/L
real function spike(x,L)
 
  real , intent(in) :: x       !< non-dimensional position
  real , intent(in) :: l       !< non-dimensional width
  real              :: pi      !< 3.1415926... calculated as 4*atan(1)
 
  pi    = 4.0*atan(1.0)
  spike = (1 - sin(pi*min(abs(x/l),0.5)))
end function spike
 
 
!> Surface fluxes of buoyancy for the Neverland configurations.
subroutine neverland_buoyancy_forcing(sfc_state, fluxes, day, dt, G, US, CS)
  type(surface),                 intent(inout) :: sfc_state !< A structure containing fields that
                                                    !! describe the surface state of the ocean.
  type(forcing),                 intent(inout) :: fluxes !< Forcing fields.
  type(time_type),               intent(in)    :: day !< Time used for determining the fluxes.
  real,                          intent(in)    :: dt !< Forcing time step (s).
  type(ocean_grid_type),         intent(inout) :: g  !< Grid structure.
  type(unit_scale_type),         intent(in)    :: us !< A dimensional unit scaling type
  type(neverland_surface_forcing_cs), pointer  :: cs !< Control structure for this module.
  ! Local variables
  real :: buoy_rest_const  ! A constant relating density anomalies to the
                           ! restoring buoyancy flux [L2 m3 T-3 kg-1 ~> m5 s-3 kg-1].
  real :: density_restore  ! De
  integer :: i, j, is, ie, js, je
  integer :: isd, ied, jsd, jed
 
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
 
  !   When modifying the code, comment out this error message.  It is here
  ! so that the original (unmodified) version is not accidentally used.
 
  ! Allocate and zero out the forcing arrays, as necessary.  This portion is
  ! usually not changed.
  if (cs%use_temperature) then
    call mom_error(fatal, "Neverland_buoyancy_forcing: " // &
        "Temperature and salinity mode not coded!" )
  else
    ! This is the buoyancy only mode.
    call safe_alloc_ptr(fluxes%buoy, isd, ied, jsd, jed)
  endif
 
 
  ! MODIFY THE CODE IN THE FOLLOWING LOOPS TO SET THE BUOYANCY FORCING TERMS.
  if (cs%restorebuoy .and. cs%first_call) then
    call safe_alloc_ptr(cs%buoy_restore, isd, ied, jsd, jed)
    cs%first_call = .false.
    ! Set CS%buoy_restore(i,j) here
  endif
 
  if ( cs%use_temperature ) then
    call mom_error(fatal, "Neverland_buoyancy_surface_forcing: " // &
      "Temperature/salinity restoring not coded!" )
  else ! This is the buoyancy only mode.
    do j=js,je ; do i=is,ie
      !   fluxes%buoy is the buoyancy flux into the ocean [L2 T-3 ~> m2 s-3].  A positive
      ! buoyancy flux is of the same sign as heating the ocean.
      fluxes%buoy(i,j) = 0.0 * g%mask2dT(i,j)
    enddo ; enddo
  endif
 
  if (cs%restorebuoy) then
    if (cs%use_temperature) then
      call mom_error(fatal, "Neverland_buoyancy_surface_forcing: " // &
        "Temperature/salinity restoring not coded!" )
    else
      !   When modifying the code, comment out this error message.  It is here
      ! so that the original (unmodified) version is not accidentally used.
 
      ! The -1 is because density has the opposite sign to buoyancy.
      buoy_rest_const = -1.0 * (cs%G_Earth * us%m_to_Z*us%T_to_s*cs%Flux_const) / cs%Rho0
      do j=js,je ; do i=is,ie
       !   Set density_restore to an expression for the surface potential
       ! density [kg m-3] that is being restored toward.
        density_restore = 1030.0
 
        fluxes%buoy(i,j) = g%mask2dT(i,j) * buoy_rest_const * &
                          (density_restore - sfc_state%sfc_density(i,j))
      enddo ; enddo
    endif
  endif                                             ! end RESTOREBUOY
 
end subroutine neverland_buoyancy_forcing
 
!> Initializes the Neverland control structure.
subroutine neverland_surface_forcing_init(Time, G, US, param_file, diag, CS)
  type(time_type),         intent(in) :: time       !< The current model time.
  type(ocean_grid_type),   intent(in) :: g          !< The ocean's grid structure.
  type(unit_scale_type),   intent(in) :: us         !< A dimensional unit scaling type
  type(param_file_type),   intent(in) :: param_file !< A structure indicating the open file to parse for
                                                    !! model parameter values.
  type(diag_ctrl), target, intent(in) :: diag       !< A structure that is used to regulate diagnostic output.
  type(neverland_surface_forcing_cs), pointer :: cs !< A pointer that is set to point to the control structure
                                                    !! for this module
  ! This include declares and sets the variable "version".
#include "version_variable.h"
  ! Local variables
  character(len=40) :: mdl = "Neverland_surface_forcing" ! This module's name.
 
  if (associated(cs)) then
    call mom_error(warning, "Neverland_surface_forcing_init called with an associated "// &
                             "control structure.")
    return
  endif
  allocate(cs)
  cs%diag => diag
 
  ! Read all relevant parameters and write them to the model log.
  call log_version(param_file, mdl, version, "")
  call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", cs%use_temperature, &
                 "If true, Temperature and salinity are used as state "//&
                 "variables.", default=.true.)
 
  call get_param(param_file, mdl, "G_EARTH", cs%G_Earth, &
                 "The gravitational acceleration of the Earth.", &
                 units="m s-2", default = 9.80, scale=us%m_to_L**2*us%Z_to_m*us%T_to_s**2)
  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)
! call get_param(param_file, mdl, "GUST_CONST", CS%gust_const, &
!                "The background gustiness in the winds.", units="Pa", &
!                default=0.02)
 
  call get_param(param_file, mdl, "RESTOREBUOY", cs%restorebuoy, &
                 "If true, the buoyancy fluxes drive the model back "//&
                 "toward some specified surface state with a rate "//&
                 "given by FLUXCONST.", default= .false.)
 
  if (cs%restorebuoy) then
    call get_param(param_file, mdl, "FLUXCONST", cs%flux_const, &
                 "The constant that relates the restoring surface fluxes "//&
                 "to the relative surface anomalies (akin to a piston "//&
                 "velocity).  Note the non-MKS units.", units="m day-1", &
                 fail_if_missing=.true.)
    ! Convert CS%flux_const from m day-1 to m s-1.
    cs%flux_const = cs%flux_const / 86400.0
  endif
 
end subroutine neverland_surface_forcing_init
 
end module neverland_surface_forcing