+*Obsolete WIND_CONFIG = "SCM_ideal_hurr" (#770)

Changed the code to issue a fatal error message when WIND_CONFIG =
"SCM_ideal_hurr", with the message including instructions on how to recover
mathematically equivalent solutions, and eliminated the subroutine
SCM_idealized_hurricane_wind_forcing() from the Idealized_hurricane module.  The
ocean_only MOM_parameter_doc files have also been modified to reflect that
"SCM_ideal_hurr" is no longer a valid setting for WIND_CONFIG.   All answers are
bitwise identical in any cases that run, but some cases may fail during
initialization with instructions on how to fix them.

Author: Hallberg-NOAA

config_src/drivers/solo_driver/MOM_surface_forcing.F90

@@ -45,9 +45,8 @@ module MOM_surface_forcing
 use user_surface_forcing,    only : USER_surface_forcing_init, user_surface_forcing_CS
 use user_revise_forcing,     only : user_alter_forcing, user_revise_forcing_init
 use user_revise_forcing,     only : user_revise_forcing_CS
-use idealized_hurricane, only : idealized_hurricane_wind_init
-use idealized_hurricane, only : idealized_hurricane_wind_forcing, SCM_idealized_hurricane_wind_forcing
-use idealized_hurricane, only : idealized_hurricane_CS
+use idealized_hurricane,     only : idealized_hurricane_wind_forcing
+use idealized_hurricane,     only : idealized_hurricane_wind_init, idealized_hurricane_CS
 use SCM_CVmix_tests,         only : SCM_CVmix_tests_surface_forcing_init
 use SCM_CVmix_tests,         only : SCM_CVmix_tests_wind_forcing
 use SCM_CVmix_tests,         only : SCM_CVmix_tests_buoyancy_forcing
@@ -297,7 +296,8 @@ subroutine set_forcing(sfc_state, forces, fluxes, day_start, day_interval, G, US
     elseif (trim(CS%wind_config) == "ideal_hurr") then
       call idealized_hurricane_wind_forcing(sfc_state, forces, day_center, G, US, CS%idealized_hurricane_CSp)
     elseif (trim(CS%wind_config) == "SCM_ideal_hurr") then
-      call SCM_idealized_hurricane_wind_forcing(sfc_state, forces, day_center, G, US, CS%idealized_hurricane_CSp)
+      call MOM_error(FATAL, "MOM_surface_forcing (set_forcing): "//&
+                            'WIND_CONFIG = "SCM_ideal_hurr" is a depricated option.')
     elseif (trim(CS%wind_config) == "SCM_CVmix_tests") then
       call SCM_CVmix_tests_wind_forcing(sfc_state, forces, day_center, G, US, CS%SCM_CVmix_tests_CSp)
     elseif (trim(CS%wind_config) == "USER") then
@@ -1780,8 +1780,8 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, tracer_flow_C
   call get_param(param_file, mdl, "WIND_CONFIG", CS%wind_config, &
                  "The character string that indicates how wind forcing is specified.  Valid "//&
                  "options include (file), (data_override), (2gyre), (1gyre), (gyres), (zero), "//&
-                 "(const), (Neverworld), (scurves), (ideal_hurr), (SCM_ideal_hurr), "//&
-                 "(SCM_CVmix_tests) and (USER).", default="zero")
+                 "(const), (Neverworld), (scurves), (ideal_hurr), (SCM_CVmix_tests) and (USER).", &
+                 default="zero")
   if (trim(CS%wind_config) == "file") then
     call get_param(param_file, mdl, "WIND_FILE", CS%wind_file, &
                  "The file in which the wind stresses are found in "//&
@@ -1990,9 +1990,13 @@ subroutine surface_forcing_init(Time, G, US, param_file, diag, CS, tracer_flow_C
     call dumbbell_surface_forcing_init(Time, G, US, param_file, diag, CS%dumbbell_forcing_CSp)
   elseif (trim(CS%wind_config) == "MESO" .or. trim(CS%buoy_config) == "MESO" ) then
     call MESO_surface_forcing_init(Time, G, US, param_file, diag, CS%MESO_forcing_CSp)
-  elseif (trim(CS%wind_config) == "ideal_hurr" .or.&
-          trim(CS%wind_config) == "SCM_ideal_hurr") then
+  elseif (trim(CS%wind_config) == "ideal_hurr") then
     call idealized_hurricane_wind_init(Time, G, US, param_file, CS%idealized_hurricane_CSp)
+  elseif (trim(CS%wind_config) == "SCM_ideal_hurr") then
+    call MOM_error(FATAL, "MOM_surface_forcing (surface_forcing_init): "//&
+          'WIND_CONFIG = "SCM_ideal_hurr" is a depricated option.  '//&
+          'To obtain mathematically equivalent results set '//&
+          'WIND_CONFIG = "ideal_hurr", IDL_HURR_SCM = True and IDL_HURR_X0 = 6.48e+05.')
   elseif (trim(CS%wind_config) == "const") then
     call get_param(param_file, mdl, "CONST_WIND_TAUX", CS%tau_x0, &
                  "With wind_config const, this is the constant zonal wind-stress", &

src/user/Idealized_Hurricane.F90

@@ -14,9 +14,7 @@ module Idealized_hurricane
 !               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
-!
+! December 2024: Removed the legacy subroutine SCM_idealized_hurricane_wind_forcing
 
 use MOM_error_handler, only : MOM_error, FATAL
 use MOM_file_parser, only : get_param, log_version, param_file_type
@@ -37,8 +35,6 @@ module Idealized_hurricane
                                      ! 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
@@ -656,196 +652,6 @@ subroutine idealized_hurricane_wind_profile(CS, US, absf, YY, XX, UOCN, VOCN, Tx
   TY = US%L_to_Z * CS%rho_a * Cd * du10 * 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 :: 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    ! The radius of the maximum winds raised to the power given by B, used in the
-               ! wind profile expression, in [km^B]
-  real :: B    ! A power used in the wind profile expression [nondim]
-  real :: C    ! A temporary variable in units of the square root of a specific volume [sqrt(m3 kg-1)]
-  real :: rad  ! The distance from the hurricane center [L ~> m]
-  real :: radius10  ! The distance from the hurricane center to its edge [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  ! Time at which the eye crosses the origin [T ~> s]
-  real :: dP  ! The pressure difference across the hurricane [R L2 T-2 ~> Pa]
-  real :: rB  ! The distance from the center raised to the power given by B, in [m^B]
-              ! or [km^B] if BR_Bench is true.
-  real :: Cd  ! Air-sea drag coefficient [nondim]
-  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     ! The wind inflow angle (positive outward) [radians]
-  real :: Rstr     ! A function of the position normalized by the radius of maximum winds [nondim]
-  real :: A0       ! The axisymmetric inflow angle [degrees]
-  real :: A1       ! The inflow angle asymmetry [degrees]
-  real :: P1       ! The angle difference between the translation direction and the inflow direction [radians]
-  real :: Adir     ! The angle of the direction from the center to a point [radians]
-  real :: transdir   ! Translation direction [radians]
-  real :: V_TS, U_TS ! Components of the translation speed [L T-1 ~> m s-1]
-
-  ! 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., tau_mag=.true.)
-
-  ! Implementing Holland (1980) parameteric wind profile
-  !------------------------------------------------------------|
-  t0 = 129600.*US%s_to_T  ! TC 'eye' crosses (0,0) at 36 hours |
-  transdir = CS%pi        ! translation direction (-x)         |
-  !------------------------------------------------------------|
-  dP = CS%pressure_ambient - CS%pressure_central
-  if (CS%answer_date < 20190101) 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 (CS%BR_Bench) then ! rho_a reset to value used in generated wind for benchmark test
-      B = C**2 * 1.2 * exp(1.0)
-    endif
-  else
-    B = (CS%max_windspeed**2 / dP ) * CS%rho_a * exp(1.0)
-  endif
-
-  if (CS%BR_Bench) then
-    A = (US%L_to_m*CS%rad_max_wind / 1000.)**B
-  else
-    A = (US%L_to_m*CS%rad_max_wind)**B
-  endif
-  ! f_local = f(x,y), but in the SCM it is constant
-  if (CS%BR_Bench) then ! (CS%SCM_mode) then
-    f_local = CS%f_column
-  else
-    f_local = G%CoriolisBu(is,js)
-  endif
-
-  ! Calculate x position relative to hurricane center as a function of time.
-  xx = (t0 - time_type_to_real(day)*US%s_to_T) * CS%hurr_translation_spd * cos(transdir)
-  rad = sqrt((xx**2) + (CS%dy_from_center**2))
-
-  ! 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 (CS%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
-
-  ! Calculate U10 in the interior (inside of the hurricane edge radius),
-  ! while adjusting U10 to 0 outside of the ambient wind radius.
-  if (rad > 0.001*CS%rad_max_wind .AND. rad < CS%rad_edge*CS%rad_max_wind) then
-    U10 = sqrt( A*B*dP*exp(-A/rB)/(CS%rho_a*rB) + 0.25*(rkm*f_local)**2 ) - 0.5*rkm*f_local
-  elseif (rad > CS%rad_edge*CS%rad_max_wind .AND. rad < CS%rad_ambient*CS%rad_max_wind) then
-    radius10 = CS%rad_max_wind*CS%rad_edge
-    if (CS%BR_Bench) then
-      rkm = radius10/1000.
-      rB = (US%L_to_m*rkm)**B
-    else
-      rkm = radius10
-      rB = (US%L_to_m*radius10)**B
-    endif
-    if (CS%edge_taper_bug) rad = radius10
-    U10 = ( sqrt( A*B*dP*exp(-A/rB)/(CS%rho_a*rB) + 0.25*(rkm*f_local)**2 ) - 0.5*rkm*f_local) &
-          * (CS%rad_ambient - rad/CS%rad_max_wind)/(CS%rad_ambient - CS%rad_edge)
-  else
-    U10 = 0.
-  endif
-  Adir = atan2(CS%dy_from_center,xx)
-
-  ! Wind angle model following Zhang and Ulhorn (2012)
-  ! ALPH is inflow angle positive outward.
-  RSTR = min(CS%rad_edge, rad / CS%rad_max_wind)
-  A0 = CS%A0_Rnorm*RSTR + CS%A0_speed*CS%max_windspeed + CS%A0_0
-  A1 = -A0*(CS%A1_Rnorm*RSTR + CS%A1_speed*CS%hurr_translation_spd + CS%A1_0)
-  P1 = (CS%P1_Rnorm*RSTR + CS%P1_speed*CS%hurr_translation_spd + CS%P1_0) * CS%pi/180.
-  ALPH = A0 - A1*cos( (TRANSDIR - ADIR ) - P1)
-  if (rad > CS%rad_edge*CS%rad_max_wind .AND. rad < CS%rad_ambient*CS%rad_max_wind) then
-    ALPH = ALPH* (CS%rad_ambient - rad/CS%rad_max_wind) / (CS%rad_ambient - CS%rad_edge)
-  elseif (rad > CS%rad_ambient*CS%rad_max_wind) then
-    ALPH = 0.0
-  endif
-  ALPH = ALPH * CS%Deg2Rad
-
-  ! 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 [R L Z T-2 ~> 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
-    ! 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)) )
-    ! 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 - CS%pi - Alph) - Uocn + U_TS
-    dV = U10*cos(Adir - Alph) - Vocn + V_TS
-    !/----------------------------------------------------|
-    !  Add a simple drag coefficient as a function of U10 |
-    !/----------------------------------------------------|
-    du10 = sqrt((du**2) + (dv**2))
-    Cd = simple_wind_scaled_Cd(u10, du10, CS)
-
-    forces%taux(I,j) = CS%rho_a * US%L_to_Z * G%mask2dCu(I,j) * Cd*du10*dU
-  enddo ; enddo
-
-  ! 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 - CS%pi - Alph) - Uocn + U_TS
-    dV = U10*cos(Adir-Alph) - Vocn + V_TS
-    du10 = sqrt((du**2) + (dv**2))
-    Cd = simple_wind_scaled_Cd(u10, du10, CS)
-    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.
-  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) = 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 ; endif
-
-  !> Set magnitude of the wind stress [R L Z T-2 ~> Pa]
-  if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie
-    forces%tau_mag(i,j) = G%mask2dT(i,j) * (CS%gustiness + &
-            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))))
-  enddo ; enddo ; endif
-
-end subroutine SCM_idealized_hurricane_wind_forcing
-
 !> This function returns the air-sea drag coefficient using a simple function of the air-sea velocity difference.
 function simple_wind_scaled_Cd(u10, du10, CS) result(Cd)
   real,                      intent(in) :: U10  !< The 10 m wind speed [L T-1 ~> m s-1]