Idealized_Hurricane.F90

!> 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 [R ~> kg m-3]
  real    :: pressure_ambient     !< Pressure at surface of ambient air [R L2 T-2 ~> Pa]
  real    :: pressure_central     !< Pressure at surface at hurricane center [R L2 T-2 ~> Pa]
  real    :: rad_max_wind         !< Radius of maximum winds [L ~> m]
  real    :: max_windspeed        !< Maximum wind speeds [L T-1 ~> m s-1]
  real    :: hurr_translation_spd !< Hurricane translation speed [L T-1 ~> m s-1]
  real    :: hurr_translation_dir !< Hurricane translation direction [radians]
  real    :: gustiness            !< Gustiness (optional, used in u*) [R L Z T-2 ~> 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 [L ~> 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 [L ~> m]
  real    :: holland_a            !< Parameter 'A' from the Holland formula [nondim]
  real    :: holland_b            !< Parameter 'B' from the Holland formula [nondim]
  real    :: holland_axbxdp       !< 'A' x 'B' x (Pressure Ambient-Pressure central)
                                  !! for the Holland prorfile calculation [R L2 T-2 ~> Pa]
  logical :: relative_tau         !< A logical to take difference between wind
                                  !! and surface currents to compute the stress
  logical :: answers_2018         !< If true, use expressions driving the idealized hurricane test
                                  !! case that recover the answers from the end of 2018.  Otherwise use
                                  !! expressions that are rescalable and respect rotational symmetry.

  ! 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 [L ~> 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

  ! Local variables
  real :: dp  ! The pressure difference across the hurricane [R L2 T-2 ~> Pa]
  real :: c
  logical :: default_2018_answers ! The default setting for the various 2018_ANSWERS flags.

  ! 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, scale=us%kg_m3_to_R)
  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., scale=us%m_s_to_L_T**2*us%kg_m3_to_R)
  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., scale=us%m_s_to_L_T**2*us%kg_m3_to_R)
  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, scale=us%m_to_L)
  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., scale=us%m_s_to_L_T)
  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, scale=us%m_s_to_L_T)
  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, scale=cs%Deg2Rad)
  call get_param(param_file, mdl, "IDL_HURR_X0", cs%Hurr_cen_X0, &
                 "Idealized Hurricane initial X position", &
                 units='m', default=0., scale=us%m_to_L)
  call get_param(param_file, mdl, "IDL_HURR_Y0", cs%Hurr_cen_Y0, &
                 "Idealized Hurricane initial Y position", &
                 units='m', default=0., scale=us%m_to_L)
  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, scale=us%m_to_L)
  call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, &
                 "This sets the default value for the various _2018_ANSWERS parameters.", &
                 default=.false.)
  call get_param(param_file, mdl, "IDL_HURR_2018_ANSWERS", cs%answers_2018, &
                 "If true, use expressions driving the idealized hurricane test case that recover "//&
                 "the answers from the end of 2018.  Otherwise use expressions that are rescalable "//&
                 "and respect rotational symmetry.", default=default_2018_answers)

  ! 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*us%kg_m3_to_R
  endif
  dp = cs%pressure_ambient - cs%pressure_central
  if (cs%answers_2018) then
    c = cs%max_windspeed / sqrt( us%R_to_kg_m3 * dp )
    cs%Holland_B = c**2 * us%R_to_kg_m3*cs%rho_a * exp(1.0)
  else
    cs%Holland_B = cs%max_windspeed**2 * cs%rho_a * exp(1.0) / dp
  endif
  cs%Holland_A = (us%L_to_m*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(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(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 components [R L Z T-2 ~> Pa]
  real :: uocn, vocn  !< Surface ocean velocity components [L T-1 ~> m s-1]
  real :: yy, xx      !< storm relative position [L ~> m]
  real :: xc, yc      !< Storm center location [L ~> m]
  real :: f_local     !< Local Coriolis parameter [T-1 ~> s-1]
  real :: fbench      !< The benchmark 'f' value [T-1 ~> s-1]
  real :: fbench_fac  !< A factor that is set to 0 to use the
                      !!  benchmark 'f' value [nondim]
  real :: rel_tau_fac !< A factor that is set to 0 to disable
                      !!  current relative stress calculation [nondim]

  ! 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)*us%s_to_T * cs%hurr_translation_spd * &
       cos(cs%hurr_translation_dir))
  yc = cs%Hurr_cen_Y0 + (time_type_to_real(day)*us%s_to_T * 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 * us%T_to_s
    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 = sfc_state%u(i,j) * rel_tau_fac
      if (cs%answers_2018) then
        vocn = 0.25*(sfc_state%v(i,j)+sfc_state%v(i+1,j-1)&
                    +sfc_state%v(i+1,j)+sfc_state%v(i,j-1))*rel_tau_fac
      else
        vocn =0.25*((sfc_state%v(i,j)+sfc_state%v(i+1,j-1)) +&
                    (sfc_state%v(i+1,j)+sfc_state%v(i,j-1))) * rel_tau_fac
      endif
      f_local = abs(0.5*(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
        yy = g%geoLatCu(i,j)*1000.*us%m_to_L - yc
        xx = g%geoLonCu(i,j)*1000.*us%m_to_L - xc
      endif
      call idealized_hurricane_wind_profile(cs, us, f_local, yy, xx, uocn, vocn, tx, ty)
      forces%taux(i,j) = g%mask2dCu(i,j) * tx
    enddo
  enddo
  !> Computes tauy
  do j=js-1,jeq
    do i=is,ie
      if (cs%answers_2018) then
        uocn = 0.25*(sfc_state%u(i,j)+sfc_state%u(i-1,j+1) + &
                     sfc_state%u(i-1,j)+sfc_state%u(i,j+1))*rel_tau_fac
      else
        uocn = 0.25*((sfc_state%u(i,j)+sfc_state%u(i-1,j+1)) + &
                     (sfc_state%u(i-1,j)+sfc_state%u(i,j+1))) * rel_tau_fac
      endif
      vocn = sfc_state%v(i,j) * rel_tau_fac
      f_local = abs(0.5*(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
        yy = g%geoLatCv(i,j)*1000.*us%m_to_L - yc
        xx = g%geoLonCv(i,j)*1000.*us%m_to_L - xc
      endif
      call idealized_hurricane_wind_profile(cs, us, f_local, yy, xx, uocn, vocn, tx, ty)
      forces%tauy(i,j) = g%mask2dCv(i,j) * 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, US, 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 idealized hurricane parameters
  type(unit_scale_type),        intent(in)    :: US     !< A dimensional unit scaling type
  real, intent(in)  :: absf !< Input Coriolis magnitude [T-1 ~> s-1]
  real, intent(in)  :: YY   !< Location in m relative to center y [L ~> m]
  real, intent(in)  :: XX   !< Location in m relative to center x [L ~> m]
  real, intent(in)  :: UOCN !< X surface current [L T-1 ~> m s-1]
  real, intent(in)  :: VOCN !< Y surface current [L T-1 ~> m s-1]
  real, intent(out) :: Tx   !< X stress [R L Z T-2 ~> Pa]
  real, intent(out) :: Ty   !< Y stress [R L Z T-2 ~> Pa]

  ! Local variables

  ! Wind profile terms
  real :: U10  ! The 10 m wind speed [L T-1 ~> m s-1]
  real :: radius    ! The distance from the hurricane center [L ~> m]
  real :: radius10  ! 10 times the distance from the hurricane center [L ~> m]
  real :: radius_km ! The distance from the hurricane center, perhaps in km [L ~> m] or [1000 L ~> km]
  real :: radiusB
  real :: tmp  ! A temporary variable [R L T-1 ~> kg m-2 s-1]
  real :: du10 ! The magnitude of the difference between the 10 m wind and the ocean flow [L T-1 ~> m s-1]
  real :: du   ! The difference between the zonal 10 m wind and the zonal ocean flow [L T-1 ~> m s-1]
  real :: dv   ! The difference between the meridional 10 m wind and the zonal ocean flow [L T-1 ~> m s-1]
  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 ! Meridional hurricane translation speed [L T-1 ~> m s-1]
  real :: U_TS ! Zonal hurricane translation speed [L T-1 ~> m s-1]

  ! 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 = (us%L_to_m*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 (cs%answers_2018) then
    if ( (radius > 0.001*cs%rad_max_wind) .and. (radius < 10.*cs%rad_max_wind) ) 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 > 10.*cs%rad_max_wind) .and. (radius < 15.*cs%rad_max_wind) ) then
      radius10 = cs%rad_max_wind*10.
      if (cs%BR_Bench) then
        radius_km = radius10/1000.
      else
        radius_km = radius10
      endif
      radiusb = (us%L_to_m*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
  else  ! This is mathematically equivalent to that is above but more accurate.
    if ( (radius > 0.001*cs%rad_max_wind) .and. (radius < 10.*cs%rad_max_wind) ) then
      tmp = ( 0.5*radius_km*absf) * (cs%rho_a*radiusb)
      u10 = (cs%Holland_AxBxDP * exp(-cs%Holland_A/radiusb)) / &
            ( tmp + sqrt(cs%Holland_AxBxDP*exp(-cs%Holland_A/radiusb) * (cs%rho_a*radiusb) + tmp**2) )
    elseif ( (radius > 10.*cs%rad_max_wind) .and. (radius < 15.*cs%rad_max_wind) ) then
      radius_km = 10.0 * cs%rad_max_wind
      if (cs%BR_Bench) radius_km = radius_km/1000.
      radiusb = (10.0*us%L_to_m*cs%rad_max_wind)**cs%Holland_B
      tmp = ( 0.5*radius_km*absf) * (cs%rho_a*radiusb)
      u10 = (3.0 - radius/(5.0*cs%rad_max_wind)) * (cs%Holland_AxBxDp*exp(-cs%Holland_A/radiusb) ) / &
            ( tmp + sqrt(cs%Holland_AxBxDp*exp(-cs%Holland_A/radiusb) * (cs%rho_a*radiusb) + tmp**2) )
    else
      u10 = 0.0
    endif
  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*us%L_T_to_m_s*cs%max_windspeed - 14.33
  a1 = -a0*(0.04*rstr + 0.05*us%L_T_to_m_s*cs%hurr_translation_spd + 0.14)
  p1 = (6.88*rstr - 9.60*us%L_T_to_m_s*cs%hurr_translation_spd + 85.31) * cs%Deg2Rad
  alph = a0 - a1*cos(cs%hurr_translation_dir-adir-p1)
  if ( (radius > 10.*cs%rad_max_wind) .and.&
       (radius < 15.*cs%rad_max_wind) ) then
     alph = alph*(15.0 - radius/cs%rad_max_wind)/5.
  elseif (radius > 15.*cs%rad_max_wind) then
     alph = 0.0
  endif
  alph = alph * cs%Deg2Rad

  ! Calculate translation speed components
  u_ts = cs%hurr_translation_spd * 0.5*cos(cs%hurr_translation_dir)
  v_ts = cs%hurr_translation_spd * 0.5*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 < 11.0*us%m_s_to_L_T) then
    cd = 1.2e-3
  elseif (du10 < 20.0*us%m_s_to_L_T) then
    if (cs%answers_2018) then
      cd = (0.49 + 0.065*us%L_T_to_m_s*u10)*1.e-3
    else
      cd = (0.49 + 0.065*us%L_T_to_m_s*du10)*1.e-3
    endif
  else
    cd = 1.8e-3
  endif

  ! Compute stress vector
  tx = us%L_to_Z * cs%rho_a * cd * sqrt(du**2 + dv**2) * du
  ty = us%L_to_Z * 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(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(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 :: du10 ! The magnitude of the difference between the 10 m wind and the ocean flow [L T-1 ~> m s-1]
  real :: u10  ! The 10 m wind speed [L T-1 ~> m s-1]
  real :: a, b, c ! For wind profile expression
  real :: rad  ! The distance from the hurricane center [L ~> m]
  real :: rkm  ! The distance from the hurricane center, sometimes scaled to km [L ~> m] or [1000 L ~> km]
  real :: f_local  ! The local Coriolis parameter [T-1 ~> s-1]
  real :: xx  ! x-position [L ~> m]
  real :: t0 !for location
  real :: dp  ! The pressure difference across the hurricane [R L2 T-2 ~> Pa]
  real :: rb
  real :: cd ! Air-sea drag coefficient
  real :: uocn, vocn ! Surface ocean velocity components [L T-1 ~> m s-1]
  real :: du, dv ! Air-sea differential motion [L T-1 ~> m s-1]
  !Wind angle variables
  real :: alph,rstr, a0, a1, p1, adir, transdir
  real :: v_ts, u_ts ! Components of the translation speed [L T-1 ~> m s-1]
  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
  if (cs%answers_2018) then
    c = cs%max_windspeed / sqrt( us%R_to_kg_m3*dp )
    b = c**2 * us%R_to_kg_m3*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
  elseif (br_bench) then ! rho_a reset to value used in generated wind for benchmark test
    b = (cs%max_windspeed**2 / dp ) * 1.2*us%kg_m3_to_R * exp(1.0)
  else
    b = (cs%max_windspeed**2 /dp ) * cs%rho_a * exp(1.0)
  endif

  a = (us%L_to_m*cs%rad_max_wind / 1000.)**b
  f_local = 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_local = 5.5659e-05*us%T_to_s
  endif
  !/ BR
  ! Calculate x position as a function of time.
  xx = us%s_to_T*( t0 - time_type_to_real(day)) * cs%hurr_translation_spd * cos(transdir)
  rad = sqrt(xx**2 + cs%dy_from_center**2)
  !/ BR
  ! rkm - rad 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 = rad/1000.
     rb = (us%L_to_m*rkm)**b
  else
     ! if not comparing to benchmark, then use correct Holland prof.
     rkm = rad
     rb = (us%L_to_m*rad)**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 (rad > 0.001*cs%rad_max_wind .AND. rad < 10.*cs%rad_max_wind) then
    u10 = sqrt( a*b*dp*exp(-a/rb)/(1.2*us%kg_m3_to_R*rb) + 0.25*(rkm*f_local)**2 ) - 0.5*rkm*f_local
  elseif (rad > 10.*cs%rad_max_wind .AND. rad < 12.*cs%rad_max_wind) then
    rad=(cs%rad_max_wind)*10.
    if (br_bench) then
       rkm = rad/1000.
       rb = (us%L_to_m*rkm)**b
    else
       rkm = rad
       rb = (us%L_to_m*rad)**b
    endif
    u10 = ( sqrt( a*b*dp*exp(-a/rb)/(1.2*us%kg_m3_to_R*rb) + 0.25*(rkm*f_local)**2 ) - 0.5*rkm*f_local) &
          * (12. - rad/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., rad / cs%rad_max_wind)
  a0 = -0.9*rstr - 0.09*us%L_T_to_m_s*cs%max_windspeed - 14.33
  a1 = -a0 *(0.04*rstr + 0.05*us%L_T_to_m_s*cs%hurr_translation_spd + 0.14)
  p1 = (6.88*rstr - 9.60*us%L_T_to_m_s*cs%hurr_translation_spd + 85.31)*pie/180.
  alph = a0 - a1*cos( (transdir - adir ) - p1)
  if (rad > 10.*cs%rad_max_wind .AND. rad < 12.*cs%rad_max_wind) then
    alph = alph* (12. - rad/cs%rad_max_wind)/2.
  elseif (rad > 12.*cs%rad_max_wind) 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*0.5*cos(transdir)
  v_ts = cs%hurr_translation_spd*0.5*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. ! sfc_state%u(I,j)
    vocn = 0. ! 0.25*( (sfc_state%v(i,J) + sfc_state%v(i+1,J-1)) + &
              !        (sfc_state%v(i+1,J) + sfc_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.0*us%m_s_to_L_T) then
      cd = 1.2e-3
    elseif (du10 < 20.0*us%m_s_to_L_T) then
      if (cs%answers_2018) then
        cd = (0.49 + 0.065 * us%L_T_to_m_s*u10 )*0.001
      else
        cd = (0.49 + 0.065 * us%L_T_to_m_s*du10 )*0.001
      endif
    else
      cd = 0.0018
    endif
    forces%taux(i,j) = cs%rho_a * us%L_to_Z * g%mask2dCu(i,j) * cd*du10*du
  enddo ; enddo
  !/BR
  ! See notes above
  do j=js-1,jeq ; do i=is,ie
    uocn = 0. ! 0.25*( (sfc_state%u(I,j) + sfc_state%u(I-1,j+1)) + &
              !        (sfc_state%u(I-1,j) + sfc_state%u(I,j+1)) )
    vocn = 0. ! sfc_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.0*us%m_s_to_L_T) then
      cd = 1.2e-3
    elseif (du10 < 20.0*us%m_s_to_L_T) then
      if (cs%answers_2018) then
        cd = (0.49 + 0.065 * us%L_T_to_m_s*u10 )*0.001
      else
        cd = (0.49 + 0.065 * us%L_T_to_m_s*du10 )*0.001
      endif
    else
      cd = 0.0018
    endif
    forces%tauy(i,j) = cs%rho_a * us%L_to_Z * g%mask2dCv(i,j) * cd*du10*dv
  enddo ; enddo
  ! Set the surface friction velocity [Z T-1 ~> 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