From efc47bcadb83130c65cafcda0cf8a8f96404ab7c Mon Sep 17 00:00:00 2001 From: Robert Hallberg Date: Fri, 5 Apr 2024 09:09:08 -0400 Subject: [PATCH] +(*)Refactor Idealize_Hurricane Refactored the Idealized_Hurricane module to clean up strange or unscaled variable units and eliminate dimensional scaling factors with the latest answer dates. This includes the introduction of 18 new runtime variables to replace hard-coded dimensional constants and two runtime logical flags to reproduce existing bugs. An inconsistency in the sign convention for the distance from the hurricane center with idealized_hurricane_wind_forcing in (the probably not yet used) single column mode was also corrected. Also added descriptions with units for all the variables in this module. By default or with appropriate parameter settings all answers are bitwise identical. --- src/user/Idealized_Hurricane.F90 | 564 +++++++++++++++++++++---------- 1 file changed, 380 insertions(+), 184 deletions(-) diff --git a/src/user/Idealized_Hurricane.F90 b/src/user/Idealized_Hurricane.F90 index 0c9d5cd330..c47366b23c 100644 --- a/src/user/Idealized_Hurricane.F90 +++ b/src/user/Idealized_Hurricane.F90 @@ -13,10 +13,9 @@ module Idealized_hurricane ! 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).. +! 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 @@ -49,6 +48,10 @@ module Idealized_hurricane 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 :: rad_edge !< Radius of the edge of the hurricane, normalized by + !! the radius of maximum winds [nondim] + real :: rad_ambient !< Radius at which the winds are at their ambient background values, + !! normalized by the radius of maximum winds [nondim] 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] @@ -60,34 +63,60 @@ module Idealized_hurricane 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 integer :: answer_date !< The vintage of the expressions in the idealized hurricane !! test case. Values below 20190101 recover the answers !! from the end of 2018, while higher values use expressions !! that are rescalable and respect rotational symmetry. + ! Parameters used in a simple wind-speed dependent expression for C_drag + real :: Cd_calm !< The drag coefficient with weak relative winds [nondim] + real :: calm_speed !< The relative wind speed below which the drag coefficient takes its + !! calm value [L T-1 ~> m s-1] + real :: Cd_windy !< The drag coefficient with strong relative winds [nondim] + real :: windy_speed !< The relative wind speed below which the drag coefficient takes its + !! windy value [L T-1 ~> m s-1] + real :: dCd_dU10 !< The partial derivative of the drag coefficient times 1000 with the 10 m + !! wind speed for intermediate wind speeds [T L-1 ~> s m-1] + real :: Cd_intercept !< The zero-wind intercept times 1000 of the linear fit for the drag + !! coefficient for the intermediate speeds where there is a linear + !! dependence on the 10 m wind speed [nondim] + + ! Parameters used to set the inflow angle as a function of radius and maximum wind speed + real :: A0_0 !< The zero-radius, zero-speed intercept of the axisymmetric inflow angle [degrees] + real :: A0_Rnorm !< The normalized radius dependence of the axisymmetric inflow angle [degrees] + real :: A0_speed !< The maximum wind speed dependence of the axisymmetric inflow angle + !! [degrees T L-1 ~> degrees s m-1] + real :: A1_0 !< The zero-radius, zero-speed intercept of the normalized inflow angle + !! asymmetry [degrees] + real :: A1_Rnorm !< The normalized radius dependence of the normalized inflow angle asymmetry [degrees] + real :: A1_speed !< The translation speed dependence of the normalized inflow angle asymmetry + !! [degrees T L-1 ~> degrees s m-1] + real :: P1_0 !< The zero-radius, zero-speed intercept of the angle difference between the + !! translation direction and the inflow direction [degrees] + real :: P1_Rnorm !< The normalized radius dependence of the angle difference between the + !! translation direction and the inflow direction [degrees] + real :: P1_speed !< The translation speed dependence of the angle difference between the + !! translation direction and the inflow direction [degrees T L-1 ~> degrees s m-1] ! 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) + logical :: SCM_mode !< If true this being used in Single Column Model mode + logical :: edge_taper_bug !< If true and SCM_mode is true, use a bug that does all of the tapering + !! and inflow angle calculations for radii between RAD_EDGE and RAD_AMBIENT + !! as though they were at RAD_EDGE. + real :: f_column !< Coriolis parameter used in the single column mode idealized + !! hurricane wind profile [T-1 ~> s-1] + 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 + real :: pi !< The circumference of a circle divided by its diameter [nondim] + real :: Deg2Rad !< The conversion factor from degrees to radians [radian degree-1] 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 @@ -102,8 +131,11 @@ subroutine idealized_hurricane_wind_init(Time, G, US, param_file, CS) ! Local variables real :: dP ! The pressure difference across the hurricane [R L2 T-2 ~> Pa] - real :: C ! A temporary variable [nondim] + real :: C ! A temporary variable in units of the square root of a specific volume [sqrt(m3 kg-1)] integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags. + logical :: continuous_Cd ! If true, use a continuous form for the simple drag coefficient as a + ! function of wind speed with the idealized hurricane. When this is false, the + ! linear shape for the mid-range wind speeds is specified separately. ! This include declares and sets the variable "version". # include "version_variable.h" @@ -132,16 +164,22 @@ subroutine idealized_hurricane_wind_init(Time, G, US, param_file, CS) 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%Pa_to_RL2_T2) - 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.", & + 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_RAD_EDGE", CS%rad_edge, & + "Radius of the edge of the hurricane, normalized by the radius of maximum winds.", & + units='nondim', default=10.0) + call get_param(param_file, mdl, "IDL_HURR_RAD_AMBIENT", CS%rad_ambient, & + "Radius at which the winds are at their ambient background values, "//& + "normalized by the radius of maximum winds.", & + units='nondim', default=CS%rad_edge+2.0) 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) + "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) + "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.", & @@ -156,17 +194,67 @@ subroutine idealized_hurricane_wind_init(Time, G, US, param_file, CS) "Current relative stress switch used in the idealized hurricane wind profile.", & default=.false.) + call get_param(param_file, mdl, "IDL_HURR_AXI_INFLOW_0", CS%A0_0, & + "The zero-radius asymmetry, zero-speed intercept of the axisymmetric inflow "//& + "angle for the parametric idealized hurricane.", & + default=-14.33, units="degrees") + call get_param(param_file, mdl, "IDL_HURR_AXI_INFLOW_RNORM", CS%A0_Rnorm, & + "The normalized radius dependence of the axisymmetric inflow angle "//& + "for the parametric idealized hurricane.", & + default=-0.9, units="degrees") + call get_param(param_file, mdl, "IDL_HURR_AXI_INFLOW_MAX_SPEED", CS%A0_speed, & + "The maximum wind speed dependence of the axisymmetric inflow angle "//& + "for the parametric idealized hurricane.", & + default=-0.09, units="degrees s m-1", scale=US%L_T_to_m_s) + call get_param(param_file, mdl, "IDL_HURR_ASYM_INFLOW_0", CS%A1_0, & + "The zero-radius, zero-speed intercept of the normalized inflow angle asymmetry "//& + "for the parametric idealized hurricane.", & + default=0.14, units="degrees") + call get_param(param_file, mdl, "IDL_HURR_ASYM_INFLOW_RNORM", CS%A1_Rnorm, & + "The normalized radius dependence of the normalized inflow angle asymmetry "//& + "for the parametric idealized hurricane.", & + default=0.04, units="degrees") + call get_param(param_file, mdl, "IDL_HURR_ASYM_INFLOW_TR_SPEED", CS%A1_speed, & + "The translation speed dependence of the normalized inflow angle asymmetry "//& + "for the parametric idealized hurricane.", & + default=0.05, units="degrees s m-1", scale=US%L_T_to_m_s) + call get_param(param_file, mdl, "IDL_HURR_INFLOW_DANGLE_0", CS%P1_0, & + "The zero-radius, zero-speed intercept of the angle difference between the "//& + "translation direction and the inflow direction "//& + "for the parametric idealized hurricane.", & + default=85.31, units="degrees") + call get_param(param_file, mdl, "IDL_HURR_INFLOW_DANGLE_RNORM", CS%P1_Rnorm, & + "The normalized radius dependence of the angle difference between the "//& + "translation direction and the inflow direction "//& + "for the parametric idealized hurricane.", & + default=6.88, units="degrees") + call get_param(param_file, mdl, "IDL_HURR_INFLOW_DANGLE_TR_SPEED", CS%P1_speed, & + "The translation speed dependence of the angle difference between the "//& + "translation direction and the inflow direction"//& + "for the parametric idealized hurricane.", & + default=-9.60, units="degrees s m-1", scale=US%L_T_to_m_s) + ! Parameters for SCM mode - call get_param(param_file, mdl, "IDL_HURR_SCM_BR_BENCH", CS%BR_BENCH, & + 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.", default=.false.) - call get_param(param_file, mdl, "IDL_HURR_SCM", CS%SCM_MODE, & + call get_param(param_file, mdl, "IDL_HURR_SCM", CS%SCM_mode, & "Single Column mode switch used in the SCM idealized hurricane wind profile.", & default=.false.) + call get_param(param_file, mdl, "IDL_HURR_SCM_EDGE_TAPER_BUG", CS%edge_taper_bug, & + "If true and IDL_HURR_SCM is true, use a bug that does all of the tapering and "//& + "inflow angle calculations for radii between RAD_EDGE and RAD_AMBIENT as though "//& + "they were at RAD_EDGE.", & + default=CS%SCM_mode, do_not_log=.not.CS%SCM_mode) !### Change the default to false. + if (.not.CS%SCM_mode) CS%edge_taper_bug = .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) + "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, "IDL_HURR_SCM_CORIOLIS", CS%f_column, & + "Coriolis parameter used in the single column mode idealized hurricane wind profile.", & + units='s-1', default=5.5659e-05, scale=US%T_to_s, do_not_log=.not.CS%BR_Bench) ! (CS%SCM_mode) + call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, & "This sets the default value for the various _ANSWER_DATE parameters.", & default=99991231) @@ -176,6 +264,48 @@ subroutine idealized_hurricane_wind_init(Time, G, US, param_file, CS) "values use expressions that are rescalable and respect rotational symmetry.", & default=default_answer_date) + ! Parameters for the simple Cdrag expression + call get_param(param_file, mdl, "IDL_HURR_CD_CALM", CS%Cd_calm, & + "The drag coefficient with weak relative winds "//& + "for the simple drag coefficient expression used with the idealized hurricane.", & + units='nondim', default=1.2e-3) + call get_param(param_file, mdl, "IDL_HURR_CD_CALM_SPEED", CS%calm_speed, & + "The relative wind speed below which the drag coefficient takes its calm value "//& + "for the simple drag coefficient expression used with the idealized hurricane.", & + units='m s-1', default=11.0, scale=US%m_s_to_L_T) + call get_param(param_file, mdl, "IDL_HURR_CD_WINDY", CS%Cd_windy, & + "The drag coefficient with strong relative winds "//& + "for the simple drag coefficient expression used with the idealized hurricane.", & + units='nondim', default=1.8e-3) + call get_param(param_file, mdl, "IDL_HURR_CD_WINDY_SPEED", CS%windy_speed, & + "The relative wind speed below which the drag coefficient takes its windy value "//& + "for the simple drag coefficient expression used with the idealized hurricane.", & + units='m s-1', default=20.0, scale=US%m_s_to_L_T) + call get_param(param_file, mdl, "IDL_HURR_CD_CONTINUOUS", continuous_Cd, & + "If true, use a continuous form for the simple drag coefficient as a function of "//& + "wind speed with the idealized hurricane. When this is false, the linear shape "//& + "for the mid-range wind speeds is specified separately.", & + default=.false.) + call get_param(param_file, mdl, "IDL_HURR_CD_DCD_DU10", CS%dCd_dU10, & + "The partial derivative of the drag coefficient times 1000 with the 10 m wind speed "//& + "for the simple drag coefficient expression used with the idealized hurricane.", & + units="s m-1", default=0.065, scale=US%L_T_to_m_s, do_not_log=continuous_Cd) + call get_param(param_file, mdl, "IDL_HURR_CD_INTERCEPT", CS%Cd_intercept, & + "The zero-wind intercept times 1000 of the linear fit for the drag coefficient "//& + "for the intermediate speeds where there is a linear dependence on the 10 m wind speed "//& + "for the simple drag coefficient expression used with the idealized hurricane.", & + units="nondim", default=0.49, do_not_log=continuous_Cd) + if (continuous_Cd) then + if (CS%windy_speed > CS%calm_speed) then + CS%dCd_dU10 = (CS%Cd_windy - CS%Cd_calm) / (CS%windy_speed - CS%calm_speed) + CS%Cd_intercept = CS%Cd_calm - CS%dCd_dU10 * CS%calm_speed + else + CS%dCd_dU10 = 0.0 + CS%Cd_intercept = CS%Cd_windy + endif + endif + + ! 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. @@ -189,9 +319,9 @@ subroutine idealized_hurricane_wind_init(Time, G, US, param_file, CS) "The background gustiness in the winds.", & units="Pa", default=0.0, scale=US%kg_m2s_to_RZ_T*US%m_s_to_L_T, do_not_log=.true.) - if (CS%BR_BENCH) then - CS%rho_a = 1.2*US%kg_m3_to_R - endif + if (CS%rad_edge >= CS%rad_ambient) call MOM_error(FATAL, & + "idealized_hurricane_wind_init: IDL_HURR_RAD_AMBIENT must be larger than IDL_HURR_RAD_EDGE.") + dP = CS%pressure_ambient - CS%pressure_central if (CS%answer_date < 20190101) then C = CS%max_windspeed / sqrt( US%R_to_kg_m3 * dP ) @@ -199,8 +329,6 @@ subroutine idealized_hurricane_wind_init(Time, G, US, param_file, CS) 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 @@ -225,6 +353,7 @@ subroutine idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, CS) 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 :: km_to_L !< The conversion factor from the units of latitude to L [L km-1 ~> 1e3] real :: rel_tau_fac !< A factor that is set to 0 to disable !! current relative stress calculation [nondim] @@ -234,6 +363,8 @@ subroutine idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, CS) isd = G%isd ; ied = G%ied ; jsd = G%jsd ; jed = G%jed IsdB = G%IsdB ; IedB = G%IedB ; JsdB = G%JsdB ; JedB = G%JedB + km_to_L = 1.0e3*US%m_to_L + ! Allocate the forcing arrays, if necessary. call allocate_mech_forcing(G, forces, stress=.true., ustar=.true., tau_mag=.true.) @@ -252,7 +383,7 @@ subroutine idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, CS) 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 = CS%f_column fbench_fac = 0.0 else fbench = 0.0 @@ -267,17 +398,17 @@ subroutine idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, CS) 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 + 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 + YY = CS%dy_from_center - YC + 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 + YY = G%geoLatCu(I,j)*km_to_L - YC + XX = G%geoLonCu(I,j)*km_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 @@ -297,11 +428,11 @@ subroutine idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, CS) 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 + YY = CS%dy_from_center - YC + 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 + YY = G%geoLatCv(i,J)*km_to_L - YC + XX = G%geoLonCv(i,J)*km_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 @@ -347,30 +478,41 @@ subroutine idealized_hurricane_wind_profile(CS, US, absf, YY, XX, UOCN, VOCN, Tx ! 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 :: radius10 ! The distance from the hurricane center to its edge [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 + real :: Cd ! The drag coefficient [nondim] + ! These variables with weird units are only used with pre-20240501 expressions + real :: radiusB ! A rescaled radius in m raised to the variable power CS%Holland_B [m^B] + real :: Holland_A ! Parameter 'A' from the Holland formula, in units of m raised to Holland_B [m^B] + real :: Holland_AxBxDP ! 'A' x 'B' x (Pressure Ambient-Pressure central) + ! for the Holland profile calculation [m^B R L2 T-2 ~> m^B Pa] + real :: tmp ! A temporary variable [m^B R L T-1 ~> m^B kg m-2 s-1] + ! These variables are used with expressions from 20240501 or later + real :: dP ! The pressure difference across the hurricane [R L2 T-2 ~> Pa] + real :: tmpA ! A temporary variable [R L2 T-2 ~> Pa] + real :: tmpB ! A temporary variable [R L T-1 ~> kg m-2 s-1] + real :: rad_max_rad_B ! The radius of maximum wind divided by the distance from the center raised + ! to the power of Holland_B [nondim] + real :: rad_rad_max ! The radius normalized by the radius of maximum winds [nondim] !Wind angle variables - real :: Alph !< The resulting inflow angle (positive outward) - real :: Rstr - real :: A0 - real :: A1 - real :: P1 - real :: Adir + 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 :: 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 + ! Implementing Holland (1980) parametric wind profile radius = SQRT(XX**2 + YY**2) + rad_rad_max = radius / CS%rad_max_wind - !/ 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 @@ -382,17 +524,24 @@ subroutine idealized_hurricane_wind_profile(CS, US, absf, YY, XX, UOCN, VOCN, Tx ! 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. + ! Calculate U10 in the interior (inside of the hurricane edge radius), + ! while adjusting U10 to 0 outside of the ambient wind radius. if (CS%answer_date < 20190101) 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) + & + radiusB = (US%L_to_m*radius)**CS%Holland_B + Holland_A = (US%L_to_m*CS%rad_max_wind)**CS%Holland_B + if ( (radius > 0.001*CS%rad_max_wind) .and. (radius < CS%rad_edge*CS%rad_max_wind) ) then + Holland_AxBxDP = Holland_A*CS%Holland_B*(CS%pressure_ambient - CS%pressure_central) + U10 = sqrt(Holland_AxBxDP*exp(-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. + elseif ( (radius > CS%rad_edge*CS%rad_max_wind) .and. (radius < CS%rad_ambient*CS%rad_max_wind) ) then + if (CS%edge_taper_bug) then ! This recreates a bug that was in SCM_idealized_hurricane_wind_forcing. + radius = CS%rad_edge * CS%rad_max_wind + rad_rad_max = CS%rad_edge + endif + + radius10 = CS%rad_max_wind*CS%rad_edge if (CS%BR_Bench) then radius_km = radius10/1000. else @@ -400,24 +549,64 @@ subroutine idealized_hurricane_wind_profile(CS, US, absf, YY, XX, UOCN, VOCN, Tx endif radiusB = (US%L_to_m*radius10)**CS%Holland_B - U10 = (sqrt(CS%Holland_AxBxDp*exp(-CS%Holland_A/radiusB) / (CS%rho_a*radiusB) + & + Holland_AxBxDP = Holland_A*CS%Holland_B*(CS%pressure_ambient - CS%pressure_central) + U10 = (sqrt(Holland_AxBxDp*exp(-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. + * (CS%rad_ambient - radius/CS%rad_max_wind) / (CS%rad_ambient - CS%rad_edge) 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 + elseif (CS%answer_date < 20240501) then + ! This is mathematically equivalent to that is above but more accurate. + radiusB = (US%L_to_m*radius)**CS%Holland_B + Holland_A = (US%L_to_m*CS%rad_max_wind)**CS%Holland_B + if ( (radius > 0.001*CS%rad_max_wind) .and. (radius < CS%rad_edge*CS%rad_max_wind) ) then + Holland_AxBxDP = Holland_A*CS%Holland_B*(CS%pressure_ambient - CS%pressure_central) 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 + U10 = (Holland_AxBxDP * exp(-Holland_A/radiusB)) / & + ( tmp + sqrt(Holland_AxBxDP*exp(-Holland_A/radiusB) * (CS%rho_a*radiusB) + tmp**2) ) + elseif ( (radius > CS%rad_edge*CS%rad_max_wind) .and. (radius < CS%rad_ambient*CS%rad_max_wind) ) then + if (CS%edge_taper_bug) then ! This recreates a bug that was in SCM_idealized_hurricane_wind_forcing. + radius = CS%rad_edge * CS%rad_max_wind + rad_rad_max = CS%rad_edge + endif + + radius_km = CS%rad_edge * 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 + radiusB = (CS%rad_edge*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) ) + Holland_AxBxDP = Holland_A*CS%Holland_B*(CS%pressure_ambient - CS%pressure_central) + U10 = ((CS%rad_ambient/(CS%rad_ambient - CS%rad_edge)) - & + radius/((CS%rad_ambient - CS%rad_edge)*CS%rad_max_wind)) * & + (Holland_AxBxDp*exp(-Holland_A/radiusB) ) / & + ( tmp + sqrt(Holland_AxBxDp*exp(-Holland_A/radiusB) * (CS%rho_a*radiusB) + tmp**2) ) + else + U10 = 0.0 + endif + else + ! This is mathematically equivalent to the expressions above, but allows for full + ! dimensional consistency testing. + dP = CS%pressure_ambient - CS%pressure_central + if ( (rad_rad_max > 0.001) .and. (rad_rad_max <= CS%rad_edge) ) then + rad_max_rad_B = (rad_rad_max)**(-CS%Holland_B) + tmpA = (rad_max_rad_B*CS%Holland_B) * dp + tmpB = (0.5*radius_km*absf) * CS%rho_a + U10 = ( tmpA * exp(-rad_max_rad_B) ) / & + ( tmpB + sqrt( (tmpA * CS%rho_a) * exp(-rad_max_rad_B) + tmpB**2) ) + elseif ( (rad_rad_max > CS%rad_edge) .and. (rad_rad_max < CS%rad_ambient) ) then + if (CS%edge_taper_bug) then ! This recreates a bug that was in SCM_idealized_hurricane_wind_forcing. + radius = CS%rad_edge * CS%rad_max_wind + rad_rad_max = CS%rad_edge + endif + + radius_km = CS%rad_edge * CS%rad_max_wind + if (CS%BR_Bench) radius_km = radius_km * 0.001 + rad_max_rad_B = CS%rad_edge**(-CS%Holland_B) + tmpA = (rad_max_rad_B*CS%Holland_B) * dp + tmpB = (0.5*radius_km*absf) * CS%rho_a + U10 = ((CS%rad_ambient - rad_rad_max) * ( tmpA * exp(-rad_max_rad_B) )) / & + ((CS%rad_ambient - CS%rad_edge) * & + ( tmpB + sqrt((tmpA * CS%rho_a) * exp(-rad_max_rad_B) + tmpB**2) ) ) else U10 = 0.0 endif @@ -429,45 +618,42 @@ subroutine idealized_hurricane_wind_profile(CS, US, absf, YY, XX, UOCN, VOCN, Tx ! 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 + RSTR = min(CS%rad_edge, rad_rad_max) + if (CS%answer_date < 20240501) then + 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%Deg2Rad + ALPH = A0 - A1*cos(CS%hurr_translation_dir-Adir-P1) + if ( (radius > CS%rad_edge*CS%rad_max_wind) .and. (radius < CS%rad_ambient*CS%rad_max_wind) ) then + ALPH = ALPH*(CS%rad_ambient - rad_rad_max) / (CS%rad_ambient - CS%rad_edge) + elseif (radius > CS%rad_ambient*CS%rad_max_wind) then ! This should be >= to avoid a jump at CS%rad_ambient + ALPH = 0.0 + endif + ALPH = ALPH * CS%Deg2Rad + else + 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%Deg2Rad + ALPH = (A0 - A1*cos((CS%hurr_translation_dir- Adir) - P1) ) * CS%Deg2Rad + if (rad_rad_max > CS%rad_edge) & + ALPH = ALPH * (max(CS%rad_ambient - rad_rad_max, 0.0) / (CS%rad_ambient - CS%rad_edge)) 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 + 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%answer_date < 20190101) 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 + Cd = simple_wind_scaled_Cd(u10, du10, CS) ! 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 - + TX = US%L_to_Z * CS%rho_a * Cd * du10 * dU + 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. @@ -484,24 +670,34 @@ subroutine SCM_idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, C ! 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 :: 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 !for location + 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 - real :: Cd ! Air-sea drag coefficient + 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,Rstr, A0, A1, P1, Adir, transdir + ! 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] - 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 @@ -509,46 +705,46 @@ subroutine SCM_idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, C 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.) - 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) | - !------------------------------------------------------| + !------------------------------------------------------------| + 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 (BR_Bench) then ! rho_a reset to value used in generated wind for benchmark test + 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 - 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) + 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 + 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 - !/ 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) + + ! 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) - !/ 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 + if (CS%BR_Bench) then rkm = rad/1000. rB = (US%L_to_m*rkm)**B else @@ -556,43 +752,42 @@ subroutine SCM_idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, C 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. + + ! 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 = rad - rB = (US%L_to_m*rad)**B + rkm = radius10 + rB = (US%L_to_m*radius10)**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. + 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) - !/ 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. + 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 > 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 + 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 * Deg2Rad - !/BR + ALPH = ALPH * CS%Deg2Rad + ! Prepare for wind calculation ! X_TS is component of translation speed added to wind vector ! due to background steering wind. @@ -604,55 +799,33 @@ subroutine SCM_idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, C ! 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 + dU = U10*sin(Adir - CS%pi - 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%answer_date < 20190101) 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 + 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 - !/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 + dU = U10*sin(Adir - CS%pi - 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%answer_date < 20190101) 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 + 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 @@ -673,4 +846,27 @@ subroutine SCM_idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, C 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] + real, intent(in) :: du10 !< The magnitude of the difference between the 10 m wind + !! and the ocean flow [L T-1 ~> m s-1] + type(idealized_hurricane_CS), pointer :: CS !< Container for SCM parameters + real :: Cd ! Air-sea drag coefficient [nondim] + + ! Note that these expressions are discontinuous at dU10 = 11 and 20 m s-1. + if (dU10 < CS%calm_speed) then + Cd = CS%Cd_calm + elseif (dU10 < CS%windy_speed) then + if (CS%answer_date < 20190101) then + Cd = (CS%Cd_intercept + CS%dCd_dU10 * U10 )*0.001 + else + Cd = (CS%Cd_intercept + CS%dCd_dU10 * dU10 )*0.001 + endif + else + Cd = CS%Cd_windy + endif + +end function simple_wind_scaled_Cd + end module idealized_hurricane