SCM_idealized_hurricane.F90

Go to the documentation of this file.

!> Initial conditions and forcing for the single column model (SCM) idealized
!! hurricane example.
module scm_idealized_hurricane

! 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_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
use mom_verticalgrid, only : verticalgrid_type

implicit none ; private

#include <MOM_memory.h>

public scm_idealized_hurricane_ts_init
public scm_idealized_hurricane_wind_init
public scm_idealized_hurricane_wind_forcing
public scm_idealized_hurricane_cs

!> Container for parameters describing idealized wind structure
type scm_idealized_hurricane_cs ; private
  real :: rho_a  !< Air density
  real :: p_n    !< Ambient pressure
  real :: p_c    !< Central pressure
  real :: r_max  !< Radius of maximum winds
  real :: u_max  !< Maximum wind speeds
  real :: yy     !< Distance (positive north) of storm center
  real :: tran_speed !< Hurricane translation speed
  real :: gust_const !< Gustiness (used in u*)
  real :: rho0   !< A reference ocean density in kg/m3
end type

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

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

contains

!> Initializes temperature and salinity for the SCM idealized hurricane example
subroutine scm_idealized_hurricane_ts_init(T, S, h, G, GV, 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_(G)), intent(out) :: T !< Potential temperature (degC)
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(out) :: S !< Salinity (psu)
  real, dimension(SZI_(G),SZJ_(G), SZK_(G)), intent(in)  :: h !< Layer thickness (m or Pa)
  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 :: S_ref, SST_ref, dTdZ, MLD
  real :: zC
  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_S_REF",s_ref, &
                 'Reference salinity', units='1e-3',default=35.0, do_not_log=just_read)
  call get_param(param_file, mdl,"SCM_SST_REF",sst_ref, &
                 'Reference surface temperature', units='C', &
                 fail_if_missing=.not.just_read, do_not_log=just_read)
  call get_param(param_file, mdl,"SCM_DTDZ",dtdz,                         &
                 'Initial temperature stratification below mixed layer', &
                 units='C/m', fail_if_missing=.not.just_read, do_not_log=just_read)
  call get_param(param_file, mdl,"SCM_MLD",mld, &
                 'Initial mixed layer depth', units='m', &
                 fail_if_missing=.not.just_read, 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)
      t(i,j,k) = sst_ref + dtdz*min(0., zc+mld)
      s(i,j,k) = s_ref
    enddo ! k
  enddo ; enddo

end subroutine scm_idealized_hurricane_ts_init

!> Initializes wind profile for the SCM idealized hurricane example
subroutine scm_idealized_hurricane_wind_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_idealized_hurricane_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_idealized_hurricane_wind_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_RHO_AIR", cs%rho_a,            &
                 "Air density "//                                     &
                 "used in the SCM idealized hurricane wind profile.", &
                 units='kg/m3', default=1.2)
  call get_param(param_file, mdl, "SCM_AMBIENT_PRESSURE", cs%p_n,     &
                 "Ambient pressure "//                                &
                 "used in the SCM idealized hurricane wind profile.", &
                 units='Pa', default=101200.)
  call get_param(param_file, mdl, "SCM_CENTRAL_PRESSURE", cs%p_c,     &
                 "Central pressure "//                                &
                 "used in the SCM idealized hurricane wind profile.", &
                 units='Pa', default=96800.)
  call get_param(param_file, mdl, "SCM_RADIUS_MAX_WINDS", cs%r_max,   &
                 "Radius of maximum winds "//                         &
                 "used in the SCM idealized hurricane wind profile.", &
                 units='m', default=50.e3)
  call get_param(param_file, mdl, "SCM_MAX_WIND_SPEED", cs%U_max,     &
                 "Maximum wind speed "//                              &
                 "used in the SCM idealized hurricane wind profile.", &
                 units='m/s', default=65.)
  call get_param(param_file, mdl, "SCM_YY", cs%YY,     &
                 "Y distance of station "//                           &
                 "used in the SCM idealized hurricane wind profile.", &
                 units='m', default=50.e3)
  call get_param(param_file, mdl, "SCM_TRAN_SPEED", cs%TRAN_SPEED,     &
                 "Translation speed of hurricane"//                   &
                 "used in the SCM idealized hurricane wind profile.", &
                 units='m/s', default=5.0)
  call get_param(param_file, mdl, "RHO_0", cs%Rho0, &
                 "The mean ocean density used with BOUSSINESQ true to \n"//&
                 "calculate accelerations and the mass for conservation \n"//&
                 "properties, or with BOUSSINSEQ false to convert some \n"//&
                 "parameters from vertical units of m to kg m-2.", &
                 units="kg m-3", default=1035.0)
  ! The following parameter is a model run-time parameter which is used
  ! and logged elsewhere and so should not be logged here. The default
  ! value should be consistent with the rest of the model.
  call get_param(param_file, mdl, "GUST_CONST", cs%gust_const, &
                 "The background gustiness in the winds.", units="Pa", &
                 default=0.00, do_not_log=.true.)


end subroutine scm_idealized_hurricane_wind_init

subroutine scm_idealized_hurricane_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_idealized_hurricane_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 :: pie, Deg2Rad
  real :: U10, A, B, C, r, f,du10,rkm ! For wind profile expression
  real :: xx, t0 !for location
  real :: dp, rB
  real :: Cd ! Air-sea drag coefficient
  real :: Uocn, Vocn ! Surface ocean velocity components
  real :: dU, dV ! Air-sea differential motion
  !Wind angle variables
  real :: Alph,Rstr, A0, A1, P1, Adir, transdir, V_TS, U_TS
  logical :: BR_Bench
  ! 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.)

  pie = 4.0*atan(1.0) ; deg2rad = pie/180.

  !/ BR
  ! Implementing Holland (1980) parameteric wind profile
  !------------------------------------------------------|
  br_bench = .true.   !true if comparing to LES runs     |
  t0 = 129600.        !TC 'eye' crosses (0,0) at 36 hours|
  transdir = pie      !translation direction (-x)        |
  !------------------------------------------------------|
  dp = cs%p_n - cs%p_c
  c = cs%U_max / sqrt( dp )
  b = c**2 * cs%rho_a * exp(1.0)
  if (br_bench) then
     ! rho_a reset to value used in generated wind for benchmark test
     b = c**2 * 1.2 * exp(1.0)
  endif
  a = (cs%r_max/1000.)**b
  f =g%CoriolisBu(is,js) ! f=f(x,y) but in the SCM is constant
  if (br_bench) then
     ! f reset to value used in generated wind for benchmark test
     f = 5.5659e-05
  endif
  !/ BR
  ! Calculate x position as a function of time.
  xx = ( t0 - time_type_to_real(day)) * cs%tran_speed * cos(transdir)
  r = sqrt(xx**2.+cs%YY**2.)
  !/ BR
  ! rkm - r converted to km for Holland prof.
  !       used in km due to error, correct implementation should
  !       not need rkm, but to match winds w/ experiment this must
  !       be maintained.  Causes winds far from storm center to be a
  !       couple of m/s higher than the correct Holland prof.
  if (br_bench) then
     rkm = r/1000.
     rb = (rkm)**b
  else
     ! if not comparing to benchmark, then use correct Holland prof.
     rkm = r
     rb = r**b
  endif
  !/ BR
  ! Calculate U10 in the interior (inside of 10x radius of maximum wind),
  ! while adjusting U10 to 0 outside of 12x radius of maximum wind.
  ! Note that rho_a is set to 1.2 following generated wind for experiment
  if (r/cs%r_max.gt.0.001 .AND. r/cs%r_max.lt.10.) then
     u10 = sqrt( a*b*dp*exp(-a/rb)/(1.2*rb) + 0.25*(rkm*f)**2 ) - 0.5*rkm*f
  elseif (r/cs%r_max.gt.10. .AND. r/cs%r_max.lt.12.) then
     r=cs%r_max*10.
     if (br_bench) then
        rkm = r/1000.
        rb=rkm**b
     else
        rkm = r
        rb = r**b
     endif
     u10 = ( sqrt( a*b*dp*exp(-a/rb)/(1.2*rb) + 0.25*(rkm*f)**2 ) - 0.5*rkm*f) &
           * (12. - r/cs%r_max)/2.
  else
     u10 = 0.
  end if
  adir = atan2(cs%YY,xx)

  !/ BR
  ! Wind angle model following Zhang and Ulhorn (2012)
  ! ALPH is inflow angle positive outward.
  rstr = min(10.,r / cs%r_max)
  a0 = -0.9*rstr -0.09*cs%U_max -14.33
  a1 = -a0 *(0.04*rstr +0.05*cs%tran_speed+0.14)
  p1 = (6.88*rstr -9.60*cs%tran_speed+85.31)*pie/180.
  alph = a0 - a1*cos( (transdir - adir ) - p1)
  if (r/cs%r_max.gt.10. .AND. r/cs%r_max.lt.12.) then
     alph = alph* (12. - r/cs%r_max)/2.
  elseif (r/cs%r_max.gt.12.) then
     alph = 0.0
  endif
  alph = alph * deg2rad

  !/BR
  ! Prepare for wind calculation
  ! X_TS is component of translation speed added to wind vector
  ! due to background steering wind.
  u_ts = cs%tran_speed/2.*cos(transdir)
  v_ts = cs%tran_speed/2.*sin(transdir)

  ! 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.
  do j=js,je ; do i=is-1,ieq
    !/BR
    ! Turn off surface current for stress calculation to be
    ! consistent with test case.
    uocn = 0.!state%u(I,j)
    vocn = 0.!0.25*( (state%v(i,J) + state%v(i+1,J-1)) &
             !    +(state%v(i+1,J) + state%v(i,J-1)) )
    !/BR
    ! Wind vector calculated from location/direction (sin/cos flipped b/c
    ! cyclonic wind is 90 deg. phase shifted from position angle).
    du = u10*sin(adir-pie-alph) - uocn + u_ts
    dv = u10*cos(adir-alph) - vocn + v_ts
    !/----------------------------------------------------|
    !BR
    !  Add a simple drag coefficient as a function of U10 |
    !/----------------------------------------------------|
    du10=sqrt(du**2+dv**2)
    if (du10.LT.11.) then
       cd = 1.2e-3
    elseif (du10.LT.20.) then
       cd = (0.49 + 0.065 * u10 )*0.001
    else
       cd = 0.0018
    endif
    fluxes%taux(i,j) = cs%rho_a * g%mask2dCu(i,j) * cd*sqrt(du**2+dv**2)*du
  enddo ; enddo
  !/BR
  ! See notes above
  do j=js-1,jeq ; do i=is,ie
    uocn = 0.!0.25*( (state%u(I,j) + state%u(I-1,j+1)) &
             !    +(state%u(I-1,j) + state%u(I,j+1)) )
    vocn = 0.!state%v(i,J)
    du = u10*sin(adir-pie-alph) - uocn + u_ts
    dv = u10*cos(adir-alph) - vocn + v_ts
    du10=sqrt(du**2+dv**2)
    if (du10.LT.11.) then
       cd = 1.2e-3
    elseif (du10.LT.20.) then
       cd = (0.49 + 0.065 * u10 )*0.001
    else
       cd = 0.0018
    endif
    fluxes%tauy(i,j) = cs%rho_a * g%mask2dCv(i,j) * cd*du10*dv
  enddo ; enddo
  ! Set the surface friction velocity, in units of m s-1. ustar is always positive.
  do j=js,je ; do i=is,ie
    !  This expression can be changed if desired, but need not be.
    fluxes%ustar(i,j) = g%mask2dT(i,j) * sqrt(cs%gust_const/cs%Rho0 + &
       sqrt(0.5*(fluxes%taux(i-1,j)**2 + fluxes%taux(i,j)**2) + &
            0.5*(fluxes%tauy(i,j-1)**2 + fluxes%tauy(i,j)**2))/cs%Rho0)
  enddo ; enddo

end subroutine scm_idealized_hurricane_wind_forcing

end module scm_idealized_hurricane