Idealized_Hurricane.F90

Go to the documentation of this file.

!> Forcing for the idealized hurricane and SCM_idealized_hurricane examples.
module idealized_hurricane
 
! This file is part of MOM6. See LICENSE.md for the license.
 
! History
!--------
! November 2014: Origination.
! October 2018: Renamed module from SCM_idealized_hurricane to idealized_hurricane
!               This module is no longer exclusively for use in SCM mode.
!               Legacy code that can be deleted is at the bottom (currently maintained
!               only to preserve exact answers in SCM mode).
!               The T/S initializations have been removed since they are redundant
!               w/ T/S initializations in CVMix_tests (which should be moved
!               into the main state_initialization to their utility
!               for multiple example cases)..
! To do
! 1. Remove the legacy SCM_idealized_hurricane_wind_forcing code
! 2. Make the hurricane-to-background wind transition a runtime parameter
!
 
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_forcing_type, only : allocate_forcing_type, allocate_mech_forcing
use mom_grid, only : ocean_grid_type
use mom_safe_alloc, only : safe_alloc_ptr
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
use mom_verticalgrid, only : verticalgrid_type
 
implicit none ; private
 
#include <MOM_memory.h>
 
public idealized_hurricane_wind_init !Public interface to initialize the idealized
                                     ! hurricane wind profile.
public idealized_hurricane_wind_forcing !Public interface to update the idealized
                                        ! hurricane wind profile.
public scm_idealized_hurricane_wind_forcing !Public interface to the legacy idealized
                                        ! hurricane wind profile for SCM.
 
!> Container for parameters describing idealized wind structure
type, public :: idealized_hurricane_cs ; private
 
  ! Parameters used to compute Holland radial wind profile
  real    :: rho_a                !< Mean air density [kg m-3]
  real    :: pressure_ambient     !< Pressure at surface of ambient air [Pa]
  real    :: pressure_central     !< Pressure at surface at hurricane center [Pa]
  real    :: rad_max_wind         !< Radius of maximum winds [m]
  real    :: max_windspeed        !< Maximum wind speeds [m s-1]
  real    :: hurr_translation_spd !< Hurricane translation speed [m s-1]
  real    :: hurr_translation_dir !< Hurricane translation speed [m s-1]
  real    :: gustiness            !< Gustiness (optional, used in u*) [R L Z T-1 ~> Pa]
  real    :: rho0                 !< A reference ocean density [R ~> kg m-3]
  real    :: hurr_cen_y0          !< The initial y position of the hurricane
                                  !!  This experiment is conducted in a Cartesian
                                  !!  grid and this is assumed to be in meters [m]
  real    :: hurr_cen_x0          !< The initial x position of the hurricane
                                  !!  This experiment is conducted in a Cartesian
                                  !!  grid and this is assumed to be in meters [m]
  real    :: holland_a            !< Parameter 'A' from the Holland formula
  real    :: holland_b            !< Parameter 'B' from the Holland formula
  real    :: holland_axbxdp       !< 'A' x 'B' x (Pressure Ambient-Pressure central)
                                  !!  for the Holland prorfile calculation
  logical :: relative_tau         !< A logical to take difference between wind
                                  !!  and surface currents to compute the stress
 
 
  ! Parameters used if in SCM (single column model) mode
  logical :: scm_mode        !< If true this being used in Single Column Model mode
  logical :: br_bench        !< A "benchmark" configuration (which is meant to
                             !!  provide identical wind to reproduce a previous
                             !!  experiment, where that wind formula contained
                             !!  an error)
  real    :: dy_from_center  !< (Fixed) distance in y from storm center path [m]
 
  ! Par
  real :: pi      !< Mathematical constant
  real :: deg2rad !< Mathematical constant
 
end type
 
! This include declares and sets the variable "version".
#include "version_variable.h"
 
character(len=40)  :: mdl = "idealized_hurricane" !< This module's name.
 
contains
 
!> Initializes wind profile for the SCM idealized hurricane example
subroutine idealized_hurricane_wind_init(Time, G, US, param_file, CS)
  type(time_type),               intent(in) :: time   !< Model time
  type(ocean_grid_type),         intent(in) :: g      !< Grid structure
  type(unit_scale_type),         intent(in) :: us     !< A dimensional unit scaling type
  type(param_file_type),         intent(in) :: param_file !< Input parameter structure
  type(idealized_hurricane_cs),  pointer    :: cs     !< Parameter container for this module
 
  real :: dp, c
 
! This include declares and sets the variable "version".
#include "version_variable.h"
 
  if (associated(cs)) then
    call mom_error(fatal, "idealized_hurricane_wind_init called "// &
                          "with an associated control structure.")
    return
  endif
 
  allocate(cs)
 
  cs%pi = 4.0*atan(1.0)
  cs%Deg2Rad = cs%pi/180.
 
  ! Read all relevant parameters and write them to the model log.
  call log_version(param_file, mdl, version, "")
 
  ! Parameters for computing a wind profile
  call get_param(param_file, mdl, "IDL_HURR_RHO_AIR", cs%rho_a, &
                 "Air density used to compute the idealized hurricane "//&
                 "wind profile.", units='kg/m3', default=1.2)
  call get_param(param_file, mdl, "IDL_HURR_AMBIENT_PRESSURE", &
                 cs%pressure_ambient, "Ambient pressure used in the "//&
                 "idealized hurricane wind profile.", units='Pa', &
                 default=101200.)
  call get_param(param_file, mdl, "IDL_HURR_CENTRAL_PRESSURE", &
                 cs%pressure_central, "Central pressure used in the "//&
                 "idealized hurricane wind profile.", units='Pa', &
                 default=96800.)
  call get_param(param_file, mdl, "IDL_HURR_RAD_MAX_WIND", &
                 cs%rad_max_wind, "Radius of maximum winds used in the "//&
                 "idealized hurricane wind profile.", units='m', &
                 default=50.e3)
  call get_param(param_file, mdl, "IDL_HURR_MAX_WIND", cs%max_windspeed, &
                 "Maximum wind speed used in the idealized hurricane"// &
                 "wind profile.", units='m/s', default=65.)
  call get_param(param_file, mdl, "IDL_HURR_TRAN_SPEED", cs%hurr_translation_spd, &
                 "Translation speed of hurricane used in the idealized "//&
                 "hurricane wind profile.", units='m/s', default=5.0)
  call get_param(param_file, mdl, "IDL_HURR_TRAN_DIR", cs%hurr_translation_dir, &
                 "Translation direction (towards) of hurricane used in the "//&
                 "idealized hurricane wind profile.", units='degrees', &
                 default=180.0)
  cs%hurr_translation_dir = cs%hurr_translation_dir * cs%Deg2Rad
  call get_param(param_file, mdl, "IDL_HURR_X0", cs%Hurr_cen_X0, &
                 "Idealized Hurricane initial X position", &
                 units='m', default=0.)
  call get_param(param_file, mdl, "IDL_HURR_Y0", cs%Hurr_cen_Y0, &
                 "Idealized Hurricane initial Y position", &
                 units='m', default=0.)
  call get_param(param_file, mdl, "IDL_HURR_TAU_CURR_REL", cs%relative_tau, &
                 "Current relative stress switch "//&
                 "used in the idealized hurricane wind profile.", &
                 units='', default=.false.)
 
  ! Parameters for SCM mode
  call get_param(param_file, mdl, "IDL_HURR_SCM_BR_BENCH", cs%BR_BENCH, &
                 "Single column mode benchmark case switch, which is "// &
                 "invoking a modification (bug) in the wind profile meant to "//&
                 "reproduce a previous implementation.", units='', default=.false.)
  call get_param(param_file, mdl, "IDL_HURR_SCM", cs%SCM_MODE, &
                 "Single Column mode switch "//&
                 "used in the SCM idealized hurricane wind profile.", &
                 units='', default=.false.)
  call get_param(param_file, mdl, "IDL_HURR_SCM_LOCY", cs%DY_from_center, &
                 "Y distance of station used in the SCM idealized hurricane "//&
                 "wind profile.", units='m', default=50.e3)
 
  ! The following parameters are model run-time parameters which are 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, "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, do_not_log=.true.)
  call get_param(param_file, mdl, "GUST_CONST", cs%gustiness, &
                 "The background gustiness in the winds.", units="Pa", &
                 default=0.0, scale=us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z, do_not_log=.true.)
 
 
  if (cs%BR_BENCH) then
    cs%rho_a = 1.2
  endif
  dp = cs%pressure_ambient - cs%pressure_central
  c = cs%max_windspeed / sqrt( dp )
  cs%Holland_B = c**2 * cs%rho_a * exp(1.0)
  cs%Holland_A = (cs%rad_max_wind)**cs%Holland_B
  cs%Holland_AxBxDP = cs%Holland_A*cs%Holland_B*dp
 
end subroutine idealized_hurricane_wind_init
 
!> Computes the surface wind for the idealized hurricane test cases
subroutine idealized_hurricane_wind_forcing(state, forces, day, G, US, CS)
  type(surface),                intent(in)    :: 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(idealized_hurricane_cs), pointer       :: cs     !< Container for idealized hurricane parameters
 
  ! Local variables
  integer :: i, j, is, ie, js, je, isq, ieq, jsq, jeq
  integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
 
  real :: tx,ty       !< wind stress
  real :: uocn, vocn  !< Surface ocean velocity components
  real :: lat, lon    !< Grid location
  real :: yy, xx      !< storm relative position
  real :: xc, yc      !< Storm center location
  real :: f           !< Coriolis
  real :: fbench      !< The benchmark 'f' value
  real :: fbench_fac  !< A factor that is set to 0 to use the
                      !!  benchmark 'f' value
  real :: rel_tau_fac !< A factor that is set to 0 to disable
                      !!  current relative stress calculation
 
  ! 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_mech_forcing(g, forces, stress=.true., ustar=.true.)
 
  if (cs%relative_tau) then
     rel_tau_fac = 1.
  else
     rel_tau_fac = 0. !Multiplied to 0 surface current
  endif
 
  !> Compute storm center location
  xc = cs%Hurr_cen_X0 + (time_type_to_real(day)*cs%hurr_translation_spd*&
       cos(cs%hurr_translation_dir))
  yc = cs%Hurr_cen_Y0 + (time_type_to_real(day)*cs%hurr_translation_spd*&
       sin(cs%hurr_translation_dir))
 
 
  if (cs%BR_Bench) then
     ! f reset to value used in generated wind for benchmark test
     fbench = 5.5659e-05
     fbench_fac = 0.0
  else
     fbench = 0.0
     fbench_fac = 1.0
  endif
 
  !> Computes taux
  do j=js,je
    do i=is-1,ieq
      uocn = state%u(i,j)*rel_tau_fac
      vocn = 0.25*(state%v(i,j)+state%v(i+1,j-1)&
             +state%v(i+1,j)+state%v(i,j-1))*rel_tau_fac
      f = abs(0.5*us%s_to_T*(g%CoriolisBu(i,j)+g%CoriolisBu(i,j-1)))*fbench_fac + fbench
      ! Calculate position as a function of time.
      if (cs%SCM_mode) then
        yy = yc + cs%dy_from_center
        xx = xc
      else
        lat = g%geoLatCu(i,j)*1000. ! Convert Lat from km to m.
        lon = g%geoLonCu(i,j)*1000. ! Convert Lon from km to m.
        yy = lat - yc
        xx = lon - xc
      endif
      call idealized_hurricane_wind_profile(cs,f,yy,xx,uocn,vocn,tx,ty)
      forces%taux(i,j) = g%mask2dCu(i,j) * us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z * tx
    enddo
  enddo
  !> Computes tauy
  do j=js-1,jeq
    do i=is,ie
      uocn = 0.25*(state%u(i,j)+state%u(i-1,j+1)&
            +state%u(i-1,j)+state%u(i,j+1))*rel_tau_fac
      vocn = state%v(i,j)*rel_tau_fac
      f = abs(0.5*us%s_to_T*(g%CoriolisBu(i-1,j)+g%CoriolisBu(i,j)))*fbench_fac + fbench
      ! Calculate position as a function of time.
      if (cs%SCM_mode) then
        yy = yc + cs%dy_from_center
        xx = xc
      else
        lat = g%geoLatCv(i,j)*1000. ! Convert Lat from km to m.
        lon = g%geoLonCv(i,j)*1000. ! Convert Lon from km to m.
        yy = lat - yc
        xx = lon - xc
      endif
      call idealized_hurricane_wind_profile(cs, f, yy, xx, uocn, vocn, tx, ty)
      forces%tauy(i,j) = g%mask2dCv(i,j) * us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z * ty
    enddo
  enddo
 
  !> Get Ustar
  do j=js,je
    do i=is,ie
      !  This expression can be changed if desired, but need not be.
      forces%ustar(i,j) = g%mask2dT(i,j) * sqrt(us%L_to_Z * (cs%gustiness/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
 
  return
end subroutine idealized_hurricane_wind_forcing
 
!> Calculate the wind speed at a location as a function of time.
subroutine idealized_hurricane_wind_profile(CS, absf, YY, XX, UOCN, VOCN, Tx, Ty)
  ! Author: Brandon Reichl
  ! Date: Nov-20-2014
  !       Aug-14-2018 Generalized for non-SCM configuration
 
  ! Input parameters
  type(idealized_hurricane_cs), &
        pointer     :: CS   !< Container for SCM parameters
  real, intent(in)  :: absf !<Input coriolis magnitude
  real, intent(in)  :: YY   !< Location in m relative to center y
  real, intent(in)  :: XX   !< Location in m relative to center x
  real, intent(in)  :: UOCN !< X surface current
  real, intent(in)  :: VOCN !< Y surface current
  real, intent(out) :: Tx   !< X stress
  real, intent(out) :: Ty   !< Y stress
 
  ! Local variables
 
  ! Wind profile terms
  real :: U10
  real :: radius
  real :: radius10
  real :: radius_km
  real :: radiusB
  real :: fcor
  real :: du10
  real :: du
  real :: dv
  real :: CD
 
  !Wind angle variables
  real :: Alph !< The resulting inflow angle (positive outward)
  real :: Rstr
  real :: A0
  real :: A1
  real :: P1
  real :: Adir
  real :: V_TS
  real :: U_TS
 
  ! Implementing Holland (1980) parameteric wind profile
 
  radius = sqrt(xx**2 + yy**2)
 
  !/ BGR
  ! 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 (cs%BR_Bench) then
    radius_km = radius/1000.
  else
    ! if not comparing to benchmark, then use correct Holland prof.
    radius_km = radius
  endif
  radiusb = (radius)**cs%Holland_B
 
  !/
  ! Calculate U10 in the interior (inside of 10x radius of maximum wind),
  ! while adjusting U10 to 0 outside of 12x radius of maximum wind.
  if ( (radius/cs%rad_max_wind .gt. 0.001) .and. &
       (radius/cs%rad_max_wind .lt. 10.) ) then
    u10 = sqrt(cs%Holland_AxBxDP*exp(-cs%Holland_A/radiusb)/(cs%rho_A*radiusb)&
                +0.25*(radius_km*absf)**2) - 0.5*radius_km*absf
  elseif ( (radius/cs%rad_max_wind .gt. 10.) .and. &
           (radius/cs%rad_max_wind .lt. 15.) ) then
 
    radius10 = cs%rad_max_wind*10.
 
    if (cs%BR_Bench) then
      radius_km = radius10/1000.
    else
      radius_km = radius10
    endif
    radiusb=radius10**cs%Holland_B
 
    u10 = (sqrt(cs%Holland_AxBxDp*exp(-cs%Holland_A/radiusb)/(cs%rho_A*radiusb)&
                  +0.25*(radius_km*absf)**2)-0.5*radius_km*absf) &
           * (15.-radius/cs%rad_max_wind)/5.
  else
    u10 = 0.
  endif
  adir = atan2(yy,xx)
  !\
 
  ! Wind angle model following Zhang and Ulhorn (2012)
  ! ALPH is inflow angle positive outward.
  rstr = min(10.,radius / cs%rad_max_wind)
  a0 = -0.9*rstr - 0.09*cs%max_windspeed - 14.33
  a1 = -a0*(0.04*rstr + 0.05*cs%Hurr_translation_spd + 0.14)
  p1 = (6.88*rstr - 9.60*cs%Hurr_translation_spd + 85.31) * cs%Deg2Rad
  alph = a0 - a1*cos(cs%hurr_translation_dir-adir-p1)
  if ( (radius/cs%rad_max_wind.gt.10.) .and.&
       (radius/cs%rad_max_wind.lt.15.) ) then
     alph = alph*(15.0-radius/cs%rad_max_wind)/5.
  elseif (radius/cs%rad_max_wind.gt.15.) then
     alph = 0.0
  endif
  alph = alph * cs%Deg2Rad
 
  ! Calculate translation speed components
  u_ts = cs%hurr_translation_spd/2.*cos(cs%hurr_translation_dir)
  v_ts = cs%hurr_translation_spd/2.*sin(cs%hurr_translation_dir)
 
  ! Set output (relative) winds
  du = u10*sin(adir-cs%Pi-alph) - uocn + u_ts
  dv = u10*cos(adir-alph) - vocn + v_ts
 
  !  Use a simple drag coefficient as a function of U10 (from Sullivan et al., 2010)
  du10 = sqrt(du**2+dv**2)
  if (du10.lt.11.) then
     cd = 1.2e-3
  elseif (du10.lt.20.0) then
     cd = (0.49 + 0.065*u10)*1.e-3
  else
     cd = 1.8e-3
  endif
 
  ! Compute stress vector
  tx = cs%rho_A * cd * sqrt(du**2 + dv**2) * du
  ty = cs%rho_A * cd * sqrt(du**2 + dv**2) * dv
 
end subroutine idealized_hurricane_wind_profile
 
!> This subroutine is primarily needed as a legacy for reproducing answers.
!! It is included as an additional subroutine rather than padded into the previous
!! routine with flags to ease its eventual removal.  Its functionality is replaced
!! with the new routines and it can be deleted when answer changes are acceptable.
subroutine scm_idealized_hurricane_wind_forcing(state, forces, day, G, US, CS)
  type(surface),                intent(in)    :: 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(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_mech_forcing(g, forces, 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%pressure_ambient - cs%pressure_central
  c = cs%max_windspeed / 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%rad_max_wind/1000.)**b
  f = us%s_to_T*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  !### A constant value in s-1.
  endif
  !/ BR
  ! Calculate x position as a function of time.
  xx = ( t0 - time_type_to_real(day)) * cs%hurr_translation_spd * cos(transdir)
  r = sqrt(xx**2 + cs%DY_from_center**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%rad_max_wind > 0.001 .AND. r/cs%rad_max_wind < 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%rad_max_wind > 10. .AND. r/cs%rad_max_wind < 12.) then
     r=cs%rad_max_wind*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%rad_max_wind)/2.
  else
     u10 = 0.
  endif
  adir = atan2(cs%DY_from_center,xx)
 
  !/ BR
  ! Wind angle model following Zhang and Ulhorn (2012)
  ! ALPH is inflow angle positive outward.
  rstr = min(10.,r / cs%rad_max_wind)
  a0 = -0.9*rstr -0.09*cs%max_windspeed - 14.33
  a1 = -a0 *(0.04*rstr +0.05*cs%hurr_translation_spd+0.14)
  p1 = (6.88*rstr -9.60*cs%hurr_translation_spd+85.31)*pie/180.
  alph = a0 - a1*cos( (transdir - adir ) - p1)
  if (r/cs%rad_max_wind > 10. .AND. r/cs%rad_max_wind < 12.) then
     alph = alph* (12. - r/cs%rad_max_wind)/2.
  elseif (r/cs%rad_max_wind > 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%hurr_translation_spd/2.*cos(transdir)
  v_ts = cs%hurr_translation_spd/2.*sin(transdir)
 
  ! Set the surface wind stresses, in [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 < 11.) then
       cd = 1.2e-3
    elseif (du10 < 20.) then
       cd = (0.49 + 0.065 * u10 )*0.001
    else
       cd = 0.0018
    endif
    forces%taux(i,j) = cs%rho_a * us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z * &
                       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 < 11.) then
       cd = 1.2e-3
    elseif (du10 < 20.) then
       cd = (0.49 + 0.065 * u10 )*0.001
    else
       cd = 0.0018
    endif
    forces%tauy(i,j) = cs%rho_a * us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z * &
                       g%mask2dCv(i,j) * cd*du10*dv
  enddo ; enddo
  ! Set the surface friction velocity [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.
    forces%ustar(i,j) = g%mask2dT(i,j) * sqrt(us%L_to_Z * (cs%gustiness/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
 
end subroutine scm_idealized_hurricane_wind_forcing
 
end module idealized_hurricane