Idealized_Hurricane.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Forcing for the idealized hurricane and SCM_idealized_hurricane examples.
7
8! History
9!--------
10! November 2014: Origination.
11! October 2018: Renamed module from SCM_idealized_hurricane to idealized_hurricane
12! This module is no longer exclusively for use in SCM mode.
13! Legacy code that can be deleted is at the bottom (currently maintained
14! only to preserve exact answers in SCM mode).
15! The T/S initializations have been removed since they are redundant
16! w/ T/S initializations in CVMix_tests (which should be moved
17! into the main state_initialization to their utility
18! for multiple example cases).
19! December 2024: Removed the legacy subroutine SCM_idealized_hurricane_wind_forcing
20
21use mom_error_handler, only : mom_error, fatal
22use mom_file_parser, only : get_param, log_version, param_file_type
24use mom_forcing_type, only : allocate_mech_forcing
25use mom_grid, only : ocean_grid_type
26use mom_safe_alloc, only : safe_alloc_ptr
27use mom_time_manager, only : time_type, operator(+), operator(/), time_to_real
31
32implicit none ; private
33
34#include <MOM_memory.h>
35
36public idealized_hurricane_wind_init !Public interface to initialize the idealized
37 ! hurricane wind profile.
38public idealized_hurricane_wind_forcing !Public interface to update the idealized
39 ! hurricane wind profile.
40
41!> Container for parameters describing idealized wind structure
42type, public :: idealized_hurricane_cs ; private
43
44 ! Parameters used to compute Holland radial wind profile
45 real :: rho_a !< Mean air density [R ~> kg m-3]
46 real :: pressure_ambient !< Pressure at surface of ambient air [R L2 T-2 ~> Pa]
47 real :: pressure_central !< Pressure at surface at hurricane center [R L2 T-2 ~> Pa]
48 real :: rad_max_wind !< Radius of maximum winds [L ~> m]
49 real :: rad_edge !< Radius of the edge of the hurricane, normalized by
50 !! the radius of maximum winds [nondim]
51 real :: rad_ambient !< Radius at which the winds are at their ambient background values,
52 !! normalized by the radius of maximum winds [nondim]
53 real :: max_windspeed !< Maximum wind speeds [L T-1 ~> m s-1]
54 real :: hurr_translation_spd !< Hurricane translation speed [L T-1 ~> m s-1]
55 real :: hurr_translation_dir !< Hurricane translation direction [radians]
56 real :: gustiness !< Gustiness (used in u*) [R Z2 T-2 ~> Pa]
57 real :: rho0 !< A reference ocean density [R ~> kg m-3]
58 real :: hurr_cen_y0 !< The initial y position of the hurricane
59 !! This experiment is conducted in a Cartesian
60 !! grid and this is assumed to be in meters [L ~> m]
61 real :: hurr_cen_x0 !< The initial x position of the hurricane
62 !! This experiment is conducted in a Cartesian
63 !! grid and this is assumed to be in meters [L ~> m]
64 real :: holland_b !< Parameter 'B' from the Holland formula [nondim]
65 logical :: relative_tau !< A logical to take difference between wind
66 !! and surface currents to compute the stress
67 integer :: answer_date !< The vintage of the expressions in the idealized hurricane
68 !! test case. Values below 20190101 recover the answers
69 !! from the end of 2018, while higher values use expressions
70 !! that are rescalable and respect rotational symmetry.
71 ! Parameters used in a simple wind-speed dependent expression for C_drag
72 real :: cd_calm !< The drag coefficient with weak relative winds [nondim]
73 real :: calm_speed !< The relative wind speed below which the drag coefficient takes its
74 !! calm value [L T-1 ~> m s-1]
75 real :: cd_windy !< The drag coefficient with strong relative winds [nondim]
76 real :: windy_speed !< The relative wind speed below which the drag coefficient takes its
77 !! windy value [L T-1 ~> m s-1]
78 real :: dcd_du10 !< The partial derivative of the drag coefficient times 1000 with the 10 m
79 !! wind speed for intermediate wind speeds [T L-1 ~> s m-1]
80 real :: cd_intercept !< The zero-wind intercept times 1000 of the linear fit for the drag
81 !! coefficient for the intermediate speeds where there is a linear
82 !! dependence on the 10 m wind speed [nondim]
83
84 ! Parameters used to set the inflow angle as a function of radius and maximum wind speed
85 real :: a0_0 !< The zero-radius, zero-speed intercept of the axisymmetric inflow angle [degrees]
86 real :: a0_rnorm !< The normalized radius dependence of the axisymmetric inflow angle [degrees]
87 real :: a0_speed !< The maximum wind speed dependence of the axisymmetric inflow angle
88 !! [degrees T L-1 ~> degrees s m-1]
89 real :: a1_0 !< The zero-radius, zero-speed intercept of the normalized inflow angle
90 !! asymmetry [degrees]
91 real :: a1_rnorm !< The normalized radius dependence of the normalized inflow angle asymmetry [degrees]
92 real :: a1_speed !< The translation speed dependence of the normalized inflow angle asymmetry
93 !! [degrees T L-1 ~> degrees s m-1]
94 real :: p1_0 !< The zero-radius, zero-speed intercept of the angle difference between the
95 !! translation direction and the inflow direction [degrees]
96 real :: p1_rnorm !< The normalized radius dependence of the angle difference between the
97 !! translation direction and the inflow direction [degrees]
98 real :: p1_speed !< The translation speed dependence of the angle difference between the
99 !! translation direction and the inflow direction [degrees T L-1 ~> degrees s m-1]
100
101 ! Parameters used if in SCM (single column model) mode
102 logical :: scm_mode !< If true this being used in Single Column Model mode
103 logical :: edge_taper_bug !< If true and SCM_mode is true, use a bug that does all of the tapering
104 !! and inflow angle calculations for radii between RAD_EDGE and RAD_AMBIENT
105 !! as though they were at RAD_EDGE.
106 real :: f_column !< Coriolis parameter used in the single column mode idealized
107 !! hurricane wind profile [T-1 ~> s-1]
108 logical :: br_bench !< A "benchmark" configuration (which is meant to
109 !! provide identical wind to reproduce a previous
110 !! experiment, where that wind formula contained an error)
111 real :: dy_from_center !< (Fixed) distance in y from storm center path [L ~> m]
112
113 real :: pi !< The circumference of a circle divided by its diameter [nondim]
114 real :: deg2rad !< The conversion factor from degrees to radians [radian degree-1]
115
116end type
117
118character(len=40) :: mdl = "idealized_hurricane" !< This module's name.
119
120contains
121
122!> Initializes wind profile for the SCM idealized hurricane example
123subroutine idealized_hurricane_wind_init(Time, G, US, param_file, CS)
124 type(time_type), intent(in) :: time !< Model time
125 type(ocean_grid_type), intent(in) :: g !< Grid structure
126 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
127 type(param_file_type), intent(in) :: param_file !< Input parameter structure
128 type(idealized_hurricane_cs), pointer :: cs !< Parameter container for this module
129
130 ! Local variables
131 real :: dp ! The pressure difference across the hurricane [R L2 T-2 ~> Pa]
132 real :: c ! A temporary variable in units of the square root of a specific volume [sqrt(m3 kg-1)]
133 integer :: default_answer_date ! The default setting for the various ANSWER_DATE flags.
134 logical :: continuous_cd ! If true, use a continuous form for the simple drag coefficient as a
135 ! function of wind speed with the idealized hurricane. When this is false, the
136 ! linear shape for the mid-range wind speeds is specified separately.
137
138 ! This include declares and sets the variable "version".
139# include "version_variable.h"
140
141 if (associated(cs)) then
142 call mom_error(fatal, "idealized_hurricane_wind_init called "// &
143 "with an associated control structure.")
144 return
145 endif
146
147 allocate(cs)
148
149 cs%pi = 4.0*atan(1.0)
150 cs%Deg2Rad = cs%pi/180.
151
152 ! Read all relevant parameters and write them to the model log.
153 call log_version(param_file, mdl, version, "")
154
155 ! Parameters for computing a wind profile
156 call get_param(param_file, mdl, "IDL_HURR_RHO_AIR", cs%rho_a, &
157 "Air density used to compute the idealized hurricane wind profile.", &
158 units='kg/m3', default=1.2, scale=us%kg_m3_to_R)
159 call get_param(param_file, mdl, "IDL_HURR_AMBIENT_PRESSURE", cs%pressure_ambient, &
160 "Ambient pressure used in the idealized hurricane wind profile.", &
161 units='Pa', default=101200., scale=us%Pa_to_RL2_T2)
162 call get_param(param_file, mdl, "IDL_HURR_CENTRAL_PRESSURE", cs%pressure_central, &
163 "Central pressure used in the idealized hurricane wind profile.", &
164 units='Pa', default=96800., scale=us%Pa_to_RL2_T2)
165 call get_param(param_file, mdl, "IDL_HURR_RAD_MAX_WIND", cs%rad_max_wind, &
166 "Radius of maximum winds used in the idealized hurricane wind profile.", &
167 units='m', default=50.e3, scale=us%m_to_L)
168 call get_param(param_file, mdl, "IDL_HURR_RAD_EDGE", cs%rad_edge, &
169 "Radius of the edge of the hurricane, normalized by the radius of maximum winds.", &
170 units='nondim', default=10.0)
171 call get_param(param_file, mdl, "IDL_HURR_RAD_AMBIENT", cs%rad_ambient, &
172 "Radius at which the winds are at their ambient background values, "//&
173 "normalized by the radius of maximum winds.", &
174 units='nondim', default=cs%rad_edge+2.0)
175 call get_param(param_file, mdl, "IDL_HURR_MAX_WIND", cs%max_windspeed, &
176 "Maximum wind speed used in the idealized hurricane wind profile.", &
177 units='m/s', default=65., scale=us%m_s_to_L_T)
178 call get_param(param_file, mdl, "IDL_HURR_TRAN_SPEED", cs%hurr_translation_spd, &
179 "Translation speed of hurricane used in the idealized hurricane wind profile.", &
180 units='m/s', default=5.0, scale=us%m_s_to_L_T)
181 call get_param(param_file, mdl, "IDL_HURR_TRAN_DIR", cs%hurr_translation_dir, &
182 "Translation direction (towards) of hurricane used in the "//&
183 "idealized hurricane wind profile.", &
184 units='degrees', default=180.0, scale=cs%Deg2Rad)
185 call get_param(param_file, mdl, "IDL_HURR_X0", cs%Hurr_cen_X0, &
186 "Idealized Hurricane initial X position", &
187 units='m', default=0., scale=us%m_to_L)
188 call get_param(param_file, mdl, "IDL_HURR_Y0", cs%Hurr_cen_Y0, &
189 "Idealized Hurricane initial Y position", &
190 units='m', default=0., scale=us%m_to_L)
191 call get_param(param_file, mdl, "IDL_HURR_TAU_CURR_REL", cs%relative_tau, &
192 "Current relative stress switch used in the idealized hurricane wind profile.", &
193 default=.false.)
194
195 call get_param(param_file, mdl, "IDL_HURR_AXI_INFLOW_0", cs%A0_0, &
196 "The zero-radius asymmetry, zero-speed intercept of the axisymmetric inflow "//&
197 "angle for the parametric idealized hurricane.", &
198 default=-14.33, units="degrees")
199 call get_param(param_file, mdl, "IDL_HURR_AXI_INFLOW_RNORM", cs%A0_Rnorm, &
200 "The normalized radius dependence of the axisymmetric inflow angle "//&
201 "for the parametric idealized hurricane.", &
202 default=-0.9, units="degrees")
203 call get_param(param_file, mdl, "IDL_HURR_AXI_INFLOW_MAX_SPEED", cs%A0_speed, &
204 "The maximum wind speed dependence of the axisymmetric inflow angle "//&
205 "for the parametric idealized hurricane.", &
206 default=-0.09, units="degrees s m-1", scale=us%L_T_to_m_s)
207 call get_param(param_file, mdl, "IDL_HURR_ASYM_INFLOW_0", cs%A1_0, &
208 "The zero-radius, zero-speed intercept of the normalized inflow angle asymmetry "//&
209 "for the parametric idealized hurricane.", &
210 default=0.14, units="degrees")
211 call get_param(param_file, mdl, "IDL_HURR_ASYM_INFLOW_RNORM", cs%A1_Rnorm, &
212 "The normalized radius dependence of the normalized inflow angle asymmetry "//&
213 "for the parametric idealized hurricane.", &
214 default=0.04, units="degrees")
215 call get_param(param_file, mdl, "IDL_HURR_ASYM_INFLOW_TR_SPEED", cs%A1_speed, &
216 "The translation speed dependence of the normalized inflow angle asymmetry "//&
217 "for the parametric idealized hurricane.", &
218 default=0.05, units="degrees s m-1", scale=us%L_T_to_m_s)
219 call get_param(param_file, mdl, "IDL_HURR_INFLOW_DANGLE_0", cs%P1_0, &
220 "The zero-radius, zero-speed intercept of the angle difference between the "//&
221 "translation direction and the inflow direction "//&
222 "for the parametric idealized hurricane.", &
223 default=85.31, units="degrees")
224 call get_param(param_file, mdl, "IDL_HURR_INFLOW_DANGLE_RNORM", cs%P1_Rnorm, &
225 "The normalized radius dependence of the angle difference between the "//&
226 "translation direction and the inflow direction "//&
227 "for the parametric idealized hurricane.", &
228 default=6.88, units="degrees")
229 call get_param(param_file, mdl, "IDL_HURR_INFLOW_DANGLE_TR_SPEED", cs%P1_speed, &
230 "The translation speed dependence of the angle difference between the "//&
231 "translation direction and the inflow direction "//&
232 "for the parametric idealized hurricane.", &
233 default=-9.60, units="degrees s m-1", scale=us%L_T_to_m_s)
234
235 ! Parameters for SCM mode
236 call get_param(param_file, mdl, "IDL_HURR_SCM_BR_BENCH", cs%BR_Bench, &
237 "Single column mode benchmark case switch, which is "// &
238 "invoking a modification (bug) in the wind profile meant to "//&
239 "reproduce a previous implementation.", default=.false.)
240 call get_param(param_file, mdl, "IDL_HURR_SCM", cs%SCM_mode, &
241 "Single Column mode switch used in the SCM idealized hurricane wind profile.", &
242 default=.false.)
243 call get_param(param_file, mdl, "IDL_HURR_SCM_EDGE_TAPER_BUG", cs%edge_taper_bug, &
244 "If true and IDL_HURR_SCM is true, use a bug that does all of the tapering and "//&
245 "inflow angle calculations for radii between RAD_EDGE and RAD_AMBIENT as though "//&
246 "they were at RAD_EDGE.", &
247 default=.false., do_not_log=.not.cs%SCM_mode)
248 if (.not.cs%SCM_mode) cs%edge_taper_bug = .false.
249 call get_param(param_file, mdl, "IDL_HURR_SCM_LOCY", cs%dy_from_center, &
250 "Y distance of station used in the SCM idealized hurricane wind profile.", &
251 units='m', default=50.e3, scale=us%m_to_L)
252 call get_param(param_file, mdl, "IDL_HURR_SCM_CORIOLIS", cs%f_column, &
253 "Coriolis parameter used in the single column mode idealized hurricane wind profile.", &
254 units='s-1', default=5.5659e-05, scale=us%T_to_s, do_not_log=.not.cs%BR_Bench) ! (CS%SCM_mode)
255
256 call get_param(param_file, mdl, "DEFAULT_ANSWER_DATE", default_answer_date, &
257 "This sets the default value for the various _ANSWER_DATE parameters.", &
258 default=99991231)
259 call get_param(param_file, mdl, "IDL_HURR_ANSWER_DATE", cs%answer_date, &
260 "The vintage of the expressions in the idealized hurricane test case. "//&
261 "Values below 20190101 recover the answers from the end of 2018, while higher "//&
262 "values use expressions that are rescalable and respect rotational symmetry.", &
263 default=default_answer_date)
264
265 ! Parameters for the simple Cdrag expression
266 call get_param(param_file, mdl, "IDL_HURR_CD_CALM", cs%Cd_calm, &
267 "The drag coefficient with weak relative winds "//&
268 "for the simple drag coefficient expression used with the idealized hurricane.", &
269 units='nondim', default=1.2e-3)
270 call get_param(param_file, mdl, "IDL_HURR_CD_CALM_SPEED", cs%calm_speed, &
271 "The relative wind speed below which the drag coefficient takes its calm value "//&
272 "for the simple drag coefficient expression used with the idealized hurricane.", &
273 units='m s-1', default=11.0, scale=us%m_s_to_L_T)
274 call get_param(param_file, mdl, "IDL_HURR_CD_WINDY", cs%Cd_windy, &
275 "The drag coefficient with strong relative winds "//&
276 "for the simple drag coefficient expression used with the idealized hurricane.", &
277 units='nondim', default=1.8e-3)
278 call get_param(param_file, mdl, "IDL_HURR_CD_WINDY_SPEED", cs%windy_speed, &
279 "The relative wind speed below which the drag coefficient takes its windy value "//&
280 "for the simple drag coefficient expression used with the idealized hurricane.", &
281 units='m s-1', default=20.0, scale=us%m_s_to_L_T)
282 call get_param(param_file, mdl, "IDL_HURR_CD_CONTINUOUS", continuous_cd, &
283 "If true, use a continuous form for the simple drag coefficient as a function of "//&
284 "wind speed with the idealized hurricane. When this is false, the linear shape "//&
285 "for the mid-range wind speeds is specified separately.", &
286 default=.false.)
287 call get_param(param_file, mdl, "IDL_HURR_CD_DCD_DU10", cs%dCd_dU10, &
288 "The partial derivative of the drag coefficient times 1000 with the 10 m wind speed "//&
289 "for the simple drag coefficient expression used with the idealized hurricane.", &
290 units="s m-1", default=0.065, scale=us%L_T_to_m_s, do_not_log=continuous_cd)
291 call get_param(param_file, mdl, "IDL_HURR_CD_INTERCEPT", cs%Cd_intercept, &
292 "The zero-wind intercept times 1000 of the linear fit for the drag coefficient "//&
293 "for the intermediate speeds where there is a linear dependence on the 10 m wind speed "//&
294 "for the simple drag coefficient expression used with the idealized hurricane.", &
295 units="nondim", default=0.49, do_not_log=continuous_cd)
296 if (continuous_cd) then
297 if (cs%windy_speed > cs%calm_speed) then
298 cs%dCd_dU10 = (cs%Cd_windy - cs%Cd_calm) / (cs%windy_speed - cs%calm_speed)
299 cs%Cd_intercept = cs%Cd_calm - cs%dCd_dU10 * cs%calm_speed
300 else
301 cs%dCd_dU10 = 0.0
302 cs%Cd_intercept = cs%Cd_windy
303 endif
304 endif
305
306
307 ! The following parameters are model run-time parameters which are used
308 ! and logged elsewhere and so should not be logged here. The default
309 ! value should be consistent with the rest of the model.
310 call get_param(param_file, mdl, "RHO_0", cs%Rho0, &
311 "The mean ocean density used with BOUSSINESQ true to "//&
312 "calculate accelerations and the mass for conservation "//&
313 "properties, or with BOUSSINESQ false to convert some "//&
314 "parameters from vertical units of m to kg m-2.", &
315 units="kg m-3", default=1035.0, scale=us%kg_m3_to_R, do_not_log=.true.)
316 call get_param(param_file, mdl, "GUST_CONST", cs%gustiness, &
317 "The background gustiness in the winds.", &
318 units="Pa", default=0.0, scale=us%Pa_to_RLZ_T2*us%L_to_Z, do_not_log=.true.)
319
320 if (cs%rad_edge >= cs%rad_ambient) call mom_error(fatal, &
321 "idealized_hurricane_wind_init: IDL_HURR_RAD_AMBIENT must be larger than IDL_HURR_RAD_EDGE.")
322
323 dp = cs%pressure_ambient - cs%pressure_central
324 if (cs%answer_date < 20190101) then
325 c = cs%max_windspeed / sqrt( us%R_to_kg_m3 * dp )
326 cs%Holland_B = c**2 * us%R_to_kg_m3*cs%rho_a * exp(1.0)
327 else
328 cs%Holland_B = cs%max_windspeed**2 * cs%rho_a * exp(1.0) / dp
329 endif
330
332
333!> Computes the surface wind for the idealized hurricane test cases
334subroutine idealized_hurricane_wind_forcing(sfc_state, forces, day, G, US, CS)
335 type(surface), intent(in) :: sfc_state !< Surface state structure
336 type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces
337 type(time_type), intent(in) :: day !< Time in days
338 type(ocean_grid_type), intent(inout) :: g !< Grid structure
339 type(unit_scale_type), intent(in) :: us !< A dimensional unit scaling type
340 type(idealized_hurricane_cs), pointer :: cs !< Container for idealized hurricane parameters
341
342 ! Local variables
343 integer :: i, j, is, ie, js, je, isq, ieq, jsq, jeq
344 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
345
346 real :: tx, ty !< wind stress components [R L Z T-2 ~> Pa]
347 real :: uocn, vocn !< Surface ocean velocity components [L T-1 ~> m s-1]
348 real :: yy, xx !< storm relative position [L ~> m]
349 real :: xc, yc !< Storm center location [L ~> m]
350 real :: f_local !< Local Coriolis parameter [T-1 ~> s-1]
351 real :: fbench !< The benchmark 'f' value [T-1 ~> s-1]
352 real :: fbench_fac !< A factor that is set to 0 to use the
353 !! benchmark 'f' value [nondim]
354 real :: rel_tau_fac !< A factor that is set to 0 to disable
355 !! current relative stress calculation [nondim]
356
357 ! Bounds for loops and memory allocation
358 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
359 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
360 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
361 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
362
363 if ((g%grid_unit_to_L <= 0.) .and. (.not.cs%SCM_mode)) call mom_error(fatal, "Idealized_Hurricane.F90: " //&
364 "idealized_hurricane_wind_forcing is only set to work with Cartesian axis units.")
365
366 ! Allocate the forcing arrays, if necessary.
367 call allocate_mech_forcing(g, forces, stress=.true., ustar=.true., tau_mag=.true.)
368
369 if (cs%relative_tau) then
370 rel_tau_fac = 1.
371 else
372 rel_tau_fac = 0. !Multiplied to 0 surface current
373 endif
374
375 !> Compute storm center location
376 xc = cs%Hurr_cen_X0 + (time_to_real(day, scale=us%s_to_T) * cs%hurr_translation_spd * &
377 cos(cs%hurr_translation_dir))
378 yc = cs%Hurr_cen_Y0 + (time_to_real(day, scale=us%s_to_T) * cs%hurr_translation_spd * &
379 sin(cs%hurr_translation_dir))
380
381 if (cs%BR_Bench) then
382 ! f reset to value used in generated wind for benchmark test
383 fbench = cs%f_column
384 fbench_fac = 0.0
385 else
386 fbench = 0.0
387 fbench_fac = 1.0
388 endif
389
390 !> Computes taux
391 do j=js,je
392 do i=is-1,ieq
393 uocn = sfc_state%u(i,j) * rel_tau_fac
394 if (cs%answer_date < 20190101) then
395 vocn = 0.25*(sfc_state%v(i,j)+sfc_state%v(i+1,j-1)&
396 +sfc_state%v(i+1,j)+sfc_state%v(i,j-1))*rel_tau_fac
397 else
398 vocn = 0.25*((sfc_state%v(i,j)+sfc_state%v(i+1,j-1)) +&
399 (sfc_state%v(i+1,j)+sfc_state%v(i,j-1))) * rel_tau_fac
400 endif
401 f_local = abs(0.5*(g%CoriolisBu(i,j)+g%CoriolisBu(i,j-1)))*fbench_fac + fbench
402 ! Calculate position as a function of time.
403 if (cs%SCM_mode) then
404 yy = cs%dy_from_center - yc
405 xx = -xc
406 else
407 yy = g%geoLatCu(i,j) * g%grid_unit_to_L - yc
408 xx = g%geoLonCu(i,j) * g%grid_unit_to_L - xc
409 endif
410 call idealized_hurricane_wind_profile(cs, us, f_local, yy, xx, uocn, vocn, tx, ty)
411 forces%taux(i,j) = g%mask2dCu(i,j) * tx
412 enddo
413 enddo
414 !> Computes tauy
415 do j=js-1,jeq
416 do i=is,ie
417 if (cs%answer_date < 20190101) then
418 uocn = 0.25*(sfc_state%u(i,j)+sfc_state%u(i-1,j+1) + &
419 sfc_state%u(i-1,j)+sfc_state%u(i,j+1))*rel_tau_fac
420 else
421 uocn = 0.25*((sfc_state%u(i,j)+sfc_state%u(i-1,j+1)) + &
422 (sfc_state%u(i-1,j)+sfc_state%u(i,j+1))) * rel_tau_fac
423 endif
424 vocn = sfc_state%v(i,j) * rel_tau_fac
425 f_local = abs(0.5*(g%CoriolisBu(i-1,j)+g%CoriolisBu(i,j)))*fbench_fac + fbench
426 ! Calculate position as a function of time.
427 if (cs%SCM_mode) then
428 yy = cs%dy_from_center - yc
429 xx = -xc
430 else
431 yy = g%geoLatCv(i,j) * g%grid_unit_to_L - yc
432 xx = g%geoLonCv(i,j) * g%grid_unit_to_L - xc
433 endif
434 call idealized_hurricane_wind_profile(cs, us, f_local, yy, xx, uocn, vocn, tx, ty)
435 forces%tauy(i,j) = g%mask2dCv(i,j) * ty
436 enddo
437 enddo
438
439 !> Get Ustar
440 if (associated(forces%ustar)) then ; do j=js,je ; do i=is,ie
441 ! This expression can be changed if desired, but need not be.
442 forces%ustar(i,j) = g%mask2dT(i,j) * sqrt(cs%gustiness/cs%Rho0 + &
443 us%L_to_Z * sqrt(0.5*((forces%taux(i-1,j)**2) + (forces%taux(i,j)**2)) + &
444 0.5*((forces%tauy(i,j-1)**2) + (forces%tauy(i,j)**2)))/cs%Rho0)
445 enddo ; enddo ; endif
446
447 !> Get tau_mag [R Z2 T-2 ~> Pa]
448 if (associated(forces%tau_mag)) then ; do j=js,je ; do i=is,ie
449 forces%tau_mag(i,j) = g%mask2dT(i,j) * (cs%gustiness + &
450 us%L_to_Z * sqrt(0.5*((forces%taux(i-1,j)**2) + (forces%taux(i,j)**2)) + &
451 0.5*((forces%tauy(i,j-1)**2) + (forces%tauy(i,j)**2))))
452 enddo ; enddo ; endif
453
455
456!> Calculate the wind speed at a location as a function of time.
457subroutine idealized_hurricane_wind_profile(CS, US, absf, YY, XX, UOCN, VOCN, Tx, Ty)
458 ! Author: Brandon Reichl
459 ! Date: Nov-20-2014
460 ! Aug-14-2018 Generalized for non-SCM configuration
461
462 ! Input parameters
463 type(idealized_hurricane_cs), pointer :: CS !< Container for idealized hurricane parameters
464 type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
465 real, intent(in) :: absf !< Input Coriolis magnitude [T-1 ~> s-1]
466 real, intent(in) :: YY !< Location in m relative to center y [L ~> m]
467 real, intent(in) :: XX !< Location in m relative to center x [L ~> m]
468 real, intent(in) :: UOCN !< X surface current [L T-1 ~> m s-1]
469 real, intent(in) :: VOCN !< Y surface current [L T-1 ~> m s-1]
470 real, intent(out) :: Tx !< X stress [R L Z T-2 ~> Pa]
471 real, intent(out) :: Ty !< Y stress [R L Z T-2 ~> Pa]
472
473 ! Local variables
474
475 ! Wind profile terms
476 real :: U10 ! The 10 m wind speed [L T-1 ~> m s-1]
477 real :: radius ! The distance from the hurricane center [L ~> m]
478 real :: radius10 ! The distance from the hurricane center to its edge [L ~> m]
479 real :: radius_km ! The distance from the hurricane center, perhaps in km [L ~> m] or [1000 L ~> km]
480 real :: du10 ! The magnitude of the difference between the 10 m wind and the ocean flow [L T-1 ~> m s-1]
481 real :: du ! The difference between the zonal 10 m wind and the zonal ocean flow [L T-1 ~> m s-1]
482 real :: dv ! The difference between the meridional 10 m wind and the zonal ocean flow [L T-1 ~> m s-1]
483 real :: Cd ! The drag coefficient [nondim]
484 ! These variables with weird units are only used with pre-20240501 expressions
485 real :: radiusB ! A rescaled radius in m raised to the variable power CS%Holland_B [m^B]
486 real :: Holland_A ! Parameter 'A' from the Holland formula, in units of m raised to Holland_B [m^B]
487 real :: Holland_AxBxDP ! 'A' x 'B' x (Pressure Ambient-Pressure central)
488 ! for the Holland profile calculation [m^B R L2 T-2 ~> m^B Pa]
489 real :: tmp ! A temporary variable [m^B R L T-1 ~> m^B kg m-2 s-1]
490 ! These variables are used with expressions from 20240501 or later
491 real :: dP ! The pressure difference across the hurricane [R L2 T-2 ~> Pa]
492 real :: tmpA ! A temporary variable [R L2 T-2 ~> Pa]
493 real :: tmpB ! A temporary variable [R L T-1 ~> kg m-2 s-1]
494 real :: rad_max_rad_B ! The radius of maximum wind divided by the distance from the center raised
495 ! to the power of Holland_B [nondim]
496 real :: rad_rad_max ! The radius normalized by the radius of maximum winds [nondim]
497
498 !Wind angle variables
499 real :: Alph ! The wind inflow angle (positive outward) [radians]
500 real :: Rstr ! A function of the position normalized by the radius of maximum winds [nondim]
501 real :: A0 ! The axisymmetric inflow angle [degrees]
502 real :: A1 ! The inflow angle asymmetry [degrees]
503 real :: P1 ! The angle difference between the translation direction and the inflow direction [radians]
504 real :: Adir ! The angle of the direction from the center to a point [radians]
505 real :: V_TS ! Meridional hurricane translation speed [L T-1 ~> m s-1]
506 real :: U_TS ! Zonal hurricane translation speed [L T-1 ~> m s-1]
507
508 ! Implementing Holland (1980) parametric wind profile
509
510 radius = sqrt((xx**2) + (yy**2))
511 rad_rad_max = radius / cs%rad_max_wind
512
513 ! rkm - r converted to km for Holland prof.
514 ! used in km due to error, correct implementation should
515 ! not need rkm, but to match winds w/ experiment this must
516 ! be maintained. Causes winds far from storm center to be a
517 ! couple of m/s higher than the correct Holland prof.
518 if (cs%BR_Bench) then
519 radius_km = radius/1000.
520 else
521 ! if not comparing to benchmark, then use correct Holland prof.
522 radius_km = radius
523 endif
524
525 !/
526 ! Calculate U10 in the interior (inside of the hurricane edge radius),
527 ! while adjusting U10 to 0 outside of the ambient wind radius.
528 if (cs%answer_date < 20190101) then
529 radiusb = (us%L_to_m*radius)**cs%Holland_B
530 holland_a = (us%L_to_m*cs%rad_max_wind)**cs%Holland_B
531 if ( (radius > 0.001*cs%rad_max_wind) .and. (radius < cs%rad_edge*cs%rad_max_wind) ) then
532 holland_axbxdp = holland_a*cs%Holland_B*(cs%pressure_ambient - cs%pressure_central)
533 u10 = sqrt(holland_axbxdp*exp(-holland_a/radiusb) / (cs%rho_a*radiusb) + &
534 0.25*(radius_km*absf)**2) - 0.5*radius_km*absf
535 elseif ( (radius > cs%rad_edge*cs%rad_max_wind) .and. (radius < cs%rad_ambient*cs%rad_max_wind) ) then
536 if (cs%edge_taper_bug) then ! This recreates a bug that was in SCM_idealized_hurricane_wind_forcing.
537 radius = cs%rad_edge * cs%rad_max_wind
538 rad_rad_max = cs%rad_edge
539 endif
540
541 radius10 = cs%rad_max_wind*cs%rad_edge
542 if (cs%BR_Bench) then
543 radius_km = radius10/1000.
544 else
545 radius_km = radius10
546 endif
547 radiusb = (us%L_to_m*radius10)**cs%Holland_B
548
549 holland_axbxdp = holland_a*cs%Holland_B*(cs%pressure_ambient - cs%pressure_central)
550 u10 = (sqrt(holland_axbxdp*exp(-holland_a/radiusb) / (cs%rho_a*radiusb) + &
551 0.25*(radius_km*absf)**2) - 0.5*radius_km*absf) &
552 * (cs%rad_ambient - radius/cs%rad_max_wind) / (cs%rad_ambient - cs%rad_edge)
553 else
554 u10 = 0.
555 endif
556 elseif (cs%answer_date < 20240501) then
557 ! This is mathematically equivalent to that is above but more accurate.
558 radiusb = (us%L_to_m*radius)**cs%Holland_B
559 holland_a = (us%L_to_m*cs%rad_max_wind)**cs%Holland_B
560 if ( (radius > 0.001*cs%rad_max_wind) .and. (radius < cs%rad_edge*cs%rad_max_wind) ) then
561 holland_axbxdp = holland_a*cs%Holland_B*(cs%pressure_ambient - cs%pressure_central)
562 tmp = ( 0.5*radius_km*absf) * (cs%rho_a*radiusb)
563 u10 = (holland_axbxdp * exp(-holland_a/radiusb)) / &
564 ( tmp + sqrt(holland_axbxdp*exp(-holland_a/radiusb) * (cs%rho_a*radiusb) + tmp**2) )
565 elseif ( (radius > cs%rad_edge*cs%rad_max_wind) .and. (radius < cs%rad_ambient*cs%rad_max_wind) ) then
566 if (cs%edge_taper_bug) then ! This recreates a bug that was in SCM_idealized_hurricane_wind_forcing.
567 radius = cs%rad_edge * cs%rad_max_wind
568 rad_rad_max = cs%rad_edge
569 endif
570
571 radius_km = cs%rad_edge * cs%rad_max_wind
572 if (cs%BR_Bench) radius_km = radius_km/1000.
573 radiusb = (cs%rad_edge*us%L_to_m*cs%rad_max_wind)**cs%Holland_B
574 tmp = ( 0.5*radius_km*absf) * (cs%rho_a*radiusb)
575 holland_axbxdp = holland_a*cs%Holland_B*(cs%pressure_ambient - cs%pressure_central)
576 u10 = ((cs%rad_ambient/(cs%rad_ambient - cs%rad_edge)) - &
577 radius/((cs%rad_ambient - cs%rad_edge)*cs%rad_max_wind)) * &
578 (holland_axbxdp*exp(-holland_a/radiusb) ) / &
579 ( tmp + sqrt(holland_axbxdp*exp(-holland_a/radiusb) * (cs%rho_a*radiusb) + tmp**2) )
580 else
581 u10 = 0.0
582 endif
583 else
584 ! This is mathematically equivalent to the expressions above, but allows for full
585 ! dimensional consistency testing.
586 dp = cs%pressure_ambient - cs%pressure_central
587 if ( (rad_rad_max > 0.001) .and. (rad_rad_max <= cs%rad_edge) ) then
588 rad_max_rad_b = (rad_rad_max)**(-cs%Holland_B)
589 tmpa = (rad_max_rad_b*cs%Holland_B) * dp
590 tmpb = (0.5*radius_km*absf) * cs%rho_a
591 u10 = ( tmpa * exp(-rad_max_rad_b) ) / &
592 ( tmpb + sqrt( (tmpa * cs%rho_a) * exp(-rad_max_rad_b) + tmpb**2) )
593 elseif ( (rad_rad_max > cs%rad_edge) .and. (rad_rad_max < cs%rad_ambient) ) then
594 if (cs%edge_taper_bug) then ! This recreates a bug that was in SCM_idealized_hurricane_wind_forcing.
595 radius = cs%rad_edge * cs%rad_max_wind
596 rad_rad_max = cs%rad_edge
597 endif
598
599 radius_km = cs%rad_edge * cs%rad_max_wind
600 if (cs%BR_Bench) radius_km = radius_km * 0.001
601 rad_max_rad_b = cs%rad_edge**(-cs%Holland_B)
602 tmpa = (rad_max_rad_b*cs%Holland_B) * dp
603 tmpb = (0.5*radius_km*absf) * cs%rho_a
604 u10 = ((cs%rad_ambient - rad_rad_max) * ( tmpa * exp(-rad_max_rad_b) )) / &
605 ((cs%rad_ambient - cs%rad_edge) * &
606 ( tmpb + sqrt((tmpa * cs%rho_a) * exp(-rad_max_rad_b) + tmpb**2) ) )
607 else
608 u10 = 0.0
609 endif
610 endif
611
612 adir = atan2(yy,xx)
613
614 !\
615
616 ! Wind angle model following Zhang and Ulhorn (2012)
617 ! ALPH is inflow angle positive outward.
618 rstr = min(cs%rad_edge, rad_rad_max)
619 if (cs%answer_date < 20240501) then
620 a0 = cs%A0_Rnorm*rstr + cs%A0_speed*cs%max_windspeed + cs%A0_0
621 a1 = -a0*(cs%A1_Rnorm*rstr + cs%A1_speed*cs%hurr_translation_spd + cs%A1_0)
622 p1 = (cs%P1_Rnorm*rstr + cs%P1_speed*cs%hurr_translation_spd + cs%P1_0) * cs%Deg2Rad
623 alph = a0 - a1*cos(cs%hurr_translation_dir-adir-p1)
624 if ( (radius > cs%rad_edge*cs%rad_max_wind) .and. (radius < cs%rad_ambient*cs%rad_max_wind) ) then
625 alph = alph*(cs%rad_ambient - rad_rad_max) / (cs%rad_ambient - cs%rad_edge)
626 elseif (radius > cs%rad_ambient*cs%rad_max_wind) then ! This should be >= to avoid a jump at CS%rad_ambient
627 alph = 0.0
628 endif
629 alph = alph * cs%Deg2Rad
630 else
631 a0 = (cs%A0_Rnorm*rstr + cs%A0_speed*cs%max_windspeed) + cs%A0_0
632 a1 = -a0*((cs%A1_Rnorm*rstr + cs%A1_speed*cs%hurr_translation_spd) + cs%A1_0)
633 p1 = ((cs%P1_Rnorm*rstr + cs%P1_speed*cs%hurr_translation_spd) + cs%P1_0) * cs%Deg2Rad
634 alph = (a0 - a1*cos((cs%hurr_translation_dir- adir) - p1) ) * cs%Deg2Rad
635 if (rad_rad_max > cs%rad_edge) &
636 alph = alph * (max(cs%rad_ambient - rad_rad_max, 0.0) / (cs%rad_ambient - cs%rad_edge))
637 endif
638
639 ! Calculate translation speed components
640 u_ts = cs%hurr_translation_spd * 0.5*cos(cs%hurr_translation_dir)
641 v_ts = cs%hurr_translation_spd * 0.5*sin(cs%hurr_translation_dir)
642
643 ! Set output (relative) winds
644 du = u10*sin(adir-cs%pi-alph) - uocn + u_ts
645 dv = u10*cos(adir-alph) - vocn + v_ts
646
647 ! Use a simple drag coefficient as a function of U10 (from Sullivan et al., 2010)
648 du10 = sqrt((du**2) + (dv**2))
649 cd = simple_wind_scaled_cd(u10, du10, cs)
650
651 ! Compute stress vector
652 tx = us%L_to_Z * cs%rho_a * cd * du10 * du
653 ty = us%L_to_Z * cs%rho_a * cd * du10 * dv
655
656!> This function returns the air-sea drag coefficient using a simple function of the air-sea velocity difference.
657function simple_wind_scaled_cd(u10, du10, CS) result(Cd)
658 real, intent(in) :: u10 !< The 10 m wind speed [L T-1 ~> m s-1]
659 real, intent(in) :: du10 !< The magnitude of the difference between the 10 m wind
660 !! and the ocean flow [L T-1 ~> m s-1]
661 type(idealized_hurricane_cs), pointer :: cs !< Container for SCM parameters
662 real :: cd ! Air-sea drag coefficient [nondim]
663
664 ! Note that these expressions are discontinuous at dU10 = 11 and 20 m s-1.
665 if (du10 < cs%calm_speed) then
666 cd = cs%Cd_calm
667 elseif (du10 < cs%windy_speed) then
668 if (cs%answer_date < 20190101) then
669 cd = (cs%Cd_intercept + cs%dCd_dU10 * u10 )*0.001
670 else
671 cd = (cs%Cd_intercept + cs%dCd_dU10 * du10 )*0.001
672 endif
673 else
674 cd = cs%Cd_windy
675 endif
676
677end function simple_wind_scaled_cd
678
679end module idealized_hurricane