MOM6
Data Types | Functions/Subroutines | Variables
idealized_hurricane Module Reference

Detailed Description

Forcing for the idealized hurricane and SCM_idealized_hurricane examples.

Data Types

type  idealized_hurricane_cs
 Container for parameters describing idealized wind structure. More...
 

Functions/Subroutines

subroutine, public idealized_hurricane_wind_init (Time, G, US, param_file, CS)
 Initializes wind profile for the SCM idealized hurricane example. More...
 
subroutine, public idealized_hurricane_wind_forcing (sfc_state, forces, day, G, US, CS)
 Computes the surface wind for the idealized hurricane test cases. More...
 
subroutine idealized_hurricane_wind_profile (CS, US, absf, YY, XX, UOCN, VOCN, Tx, Ty)
 Calculate the wind speed at a location as a function of time. More...
 
subroutine, public scm_idealized_hurricane_wind_forcing (sfc_state, forces, day, G, US, CS)
 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. More...
 

Variables

character(len=40) mdl = "idealized_hurricane"
 This module's name.
 

Function/Subroutine Documentation

◆ idealized_hurricane_wind_forcing()

subroutine, public idealized_hurricane::idealized_hurricane_wind_forcing ( type(surface), intent(in)  sfc_state,
type(mech_forcing), intent(inout)  forces,
type(time_type), intent(in)  day,
type(ocean_grid_type), intent(inout)  G,
type(unit_scale_type), intent(in)  US,
type(idealized_hurricane_cs), pointer  CS 
)

Computes the surface wind for the idealized hurricane test cases.

Parameters
[in]sfc_stateSurface state structure
[in,out]forcesA structure with the driving mechanical forces
[in]dayTime in days
[in,out]gGrid structure
[in]usA dimensional unit scaling type
csContainer for idealized hurricane parameters

Compute storm center location

Computes taux

Computes tauy

Get Ustar

Definition at line 208 of file Idealized_Hurricane.F90.

209  type(surface), intent(in) :: sfc_state !< Surface state structure
210  type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces
211  type(time_type), intent(in) :: day !< Time in days
212  type(ocean_grid_type), intent(inout) :: G !< Grid structure
213  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
214  type(idealized_hurricane_CS), pointer :: CS !< Container for idealized hurricane parameters
215 
216  ! Local variables
217  integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq
218  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
219 
220  real :: TX, TY !< wind stress components [R L Z T-2 ~> Pa]
221  real :: Uocn, Vocn !< Surface ocean velocity components [L T-1 ~> m s-1]
222  real :: YY, XX !< storm relative position [L ~> m]
223  real :: XC, YC !< Storm center location [L ~> m]
224  real :: f_local !< Local Coriolis parameter [T-1 ~> s-1]
225  real :: fbench !< The benchmark 'f' value [T-1 ~> s-1]
226  real :: fbench_fac !< A factor that is set to 0 to use the
227  !! benchmark 'f' value [nondim]
228  real :: rel_tau_fac !< A factor that is set to 0 to disable
229  !! current relative stress calculation [nondim]
230 
231  ! Bounds for loops and memory allocation
232  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
233  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
234  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
235  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
236 
237  ! Allocate the forcing arrays, if necessary.
238  call allocate_mech_forcing(g, forces, stress=.true., ustar=.true.)
239 
240  if (cs%relative_tau) then
241  rel_tau_fac = 1.
242  else
243  rel_tau_fac = 0. !Multiplied to 0 surface current
244  endif
245 
246  !> Compute storm center location
247  xc = cs%Hurr_cen_X0 + (time_type_to_real(day)*us%s_to_T * cs%hurr_translation_spd * &
248  cos(cs%hurr_translation_dir))
249  yc = cs%Hurr_cen_Y0 + (time_type_to_real(day)*us%s_to_T * cs%hurr_translation_spd * &
250  sin(cs%hurr_translation_dir))
251 
252 
253  if (cs%BR_Bench) then
254  ! f reset to value used in generated wind for benchmark test
255  fbench = 5.5659e-05 * us%T_to_s
256  fbench_fac = 0.0
257  else
258  fbench = 0.0
259  fbench_fac = 1.0
260  endif
261 
262  !> Computes taux
263  do j=js,je
264  do i=is-1,ieq
265  uocn = sfc_state%u(i,j) * rel_tau_fac
266  if (cs%answers_2018) then
267  vocn = 0.25*(sfc_state%v(i,j)+sfc_state%v(i+1,j-1)&
268  +sfc_state%v(i+1,j)+sfc_state%v(i,j-1))*rel_tau_fac
269  else
270  vocn =0.25*((sfc_state%v(i,j)+sfc_state%v(i+1,j-1)) +&
271  (sfc_state%v(i+1,j)+sfc_state%v(i,j-1))) * rel_tau_fac
272  endif
273  f_local = abs(0.5*(g%CoriolisBu(i,j)+g%CoriolisBu(i,j-1)))*fbench_fac + fbench
274  ! Calculate position as a function of time.
275  if (cs%SCM_mode) then
276  yy = yc + cs%dy_from_center
277  xx = xc
278  else
279  yy = g%geoLatCu(i,j)*1000.*us%m_to_L - yc
280  xx = g%geoLonCu(i,j)*1000.*us%m_to_L - xc
281  endif
282  call idealized_hurricane_wind_profile(cs, us, f_local, yy, xx, uocn, vocn, tx, ty)
283  forces%taux(i,j) = g%mask2dCu(i,j) * tx
284  enddo
285  enddo
286  !> Computes tauy
287  do j=js-1,jeq
288  do i=is,ie
289  if (cs%answers_2018) then
290  uocn = 0.25*(sfc_state%u(i,j)+sfc_state%u(i-1,j+1) + &
291  sfc_state%u(i-1,j)+sfc_state%u(i,j+1))*rel_tau_fac
292  else
293  uocn = 0.25*((sfc_state%u(i,j)+sfc_state%u(i-1,j+1)) + &
294  (sfc_state%u(i-1,j)+sfc_state%u(i,j+1))) * rel_tau_fac
295  endif
296  vocn = sfc_state%v(i,j) * rel_tau_fac
297  f_local = abs(0.5*(g%CoriolisBu(i-1,j)+g%CoriolisBu(i,j)))*fbench_fac + fbench
298  ! Calculate position as a function of time.
299  if (cs%SCM_mode) then
300  yy = yc + cs%dy_from_center
301  xx = xc
302  else
303  yy = g%geoLatCv(i,j)*1000.*us%m_to_L - yc
304  xx = g%geoLonCv(i,j)*1000.*us%m_to_L - xc
305  endif
306  call idealized_hurricane_wind_profile(cs, us, f_local, yy, xx, uocn, vocn, tx, ty)
307  forces%tauy(i,j) = g%mask2dCv(i,j) * ty
308  enddo
309  enddo
310 
311  !> Get Ustar
312  do j=js,je
313  do i=is,ie
314  ! This expression can be changed if desired, but need not be.
315  forces%ustar(i,j) = g%mask2dT(i,j) * sqrt(us%L_to_Z * (cs%gustiness/cs%Rho0 + &
316  sqrt(0.5*(forces%taux(i-1,j)**2 + forces%taux(i,j)**2) + &
317  0.5*(forces%tauy(i,j-1)**2 + forces%tauy(i,j)**2))/cs%Rho0))
318  enddo
319  enddo
320 
321  return

◆ idealized_hurricane_wind_init()

subroutine, public idealized_hurricane::idealized_hurricane_wind_init ( type(time_type), intent(in)  Time,
type(ocean_grid_type), intent(in)  G,
type(unit_scale_type), intent(in)  US,
type(param_file_type), intent(in)  param_file,
type(idealized_hurricane_cs), pointer  CS 
)

Initializes wind profile for the SCM idealized hurricane example.

Parameters
[in]timeModel time
[in]gGrid structure
[in]usA dimensional unit scaling type
[in]param_fileInput parameter structure
csParameter container for this module

Definition at line 95 of file Idealized_Hurricane.F90.

96  type(time_type), intent(in) :: Time !< Model time
97  type(ocean_grid_type), intent(in) :: G !< Grid structure
98  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
99  type(param_file_type), intent(in) :: param_file !< Input parameter structure
100  type(idealized_hurricane_CS), pointer :: CS !< Parameter container for this module
101 
102  ! Local variables
103  real :: dP ! The pressure difference across the hurricane [R L2 T-2 ~> Pa]
104  real :: C
105  logical :: default_2018_answers ! The default setting for the various 2018_ANSWERS flags.
106 
107  ! This include declares and sets the variable "version".
108 # include "version_variable.h"
109 
110  if (associated(cs)) then
111  call mom_error(fatal, "idealized_hurricane_wind_init called "// &
112  "with an associated control structure.")
113  return
114  endif
115 
116  allocate(cs)
117 
118  cs%pi = 4.0*atan(1.0)
119  cs%Deg2Rad = cs%pi/180.
120 
121  ! Read all relevant parameters and write them to the model log.
122  call log_version(param_file, mdl, version, "")
123 
124  ! Parameters for computing a wind profile
125  call get_param(param_file, mdl, "IDL_HURR_RHO_AIR", cs%rho_a, &
126  "Air density used to compute the idealized hurricane wind profile.", &
127  units='kg/m3', default=1.2, scale=us%kg_m3_to_R)
128  call get_param(param_file, mdl, "IDL_HURR_AMBIENT_PRESSURE", cs%pressure_ambient, &
129  "Ambient pressure used in the idealized hurricane wind profile.", &
130  units='Pa', default=101200., scale=us%m_s_to_L_T**2*us%kg_m3_to_R)
131  call get_param(param_file, mdl, "IDL_HURR_CENTRAL_PRESSURE", cs%pressure_central, &
132  "Central pressure used in the idealized hurricane wind profile.", &
133  units='Pa', default=96800., scale=us%m_s_to_L_T**2*us%kg_m3_to_R)
134  call get_param(param_file, mdl, "IDL_HURR_RAD_MAX_WIND", &
135  cs%rad_max_wind, "Radius of maximum winds used in the "//&
136  "idealized hurricane wind profile.", units='m', &
137  default=50.e3, scale=us%m_to_L)
138  call get_param(param_file, mdl, "IDL_HURR_MAX_WIND", cs%max_windspeed, &
139  "Maximum wind speed used in the idealized hurricane"// &
140  "wind profile.", units='m/s', default=65., scale=us%m_s_to_L_T)
141  call get_param(param_file, mdl, "IDL_HURR_TRAN_SPEED", cs%hurr_translation_spd, &
142  "Translation speed of hurricane used in the idealized "//&
143  "hurricane wind profile.", units='m/s', default=5.0, scale=us%m_s_to_L_T)
144  call get_param(param_file, mdl, "IDL_HURR_TRAN_DIR", cs%hurr_translation_dir, &
145  "Translation direction (towards) of hurricane used in the "//&
146  "idealized hurricane wind profile.", units='degrees', &
147  default=180.0, scale=cs%Deg2Rad)
148  call get_param(param_file, mdl, "IDL_HURR_X0", cs%Hurr_cen_X0, &
149  "Idealized Hurricane initial X position", &
150  units='m', default=0., scale=us%m_to_L)
151  call get_param(param_file, mdl, "IDL_HURR_Y0", cs%Hurr_cen_Y0, &
152  "Idealized Hurricane initial Y position", &
153  units='m', default=0., scale=us%m_to_L)
154  call get_param(param_file, mdl, "IDL_HURR_TAU_CURR_REL", cs%relative_tau, &
155  "Current relative stress switch "//&
156  "used in the idealized hurricane wind profile.", &
157  units='', default=.false.)
158 
159  ! Parameters for SCM mode
160  call get_param(param_file, mdl, "IDL_HURR_SCM_BR_BENCH", cs%BR_BENCH, &
161  "Single column mode benchmark case switch, which is "// &
162  "invoking a modification (bug) in the wind profile meant to "//&
163  "reproduce a previous implementation.", units='', default=.false.)
164  call get_param(param_file, mdl, "IDL_HURR_SCM", cs%SCM_MODE, &
165  "Single Column mode switch "//&
166  "used in the SCM idealized hurricane wind profile.", &
167  units='', default=.false.)
168  call get_param(param_file, mdl, "IDL_HURR_SCM_LOCY", cs%dy_from_center, &
169  "Y distance of station used in the SCM idealized hurricane "//&
170  "wind profile.", units='m', default=50.e3, scale=us%m_to_L)
171  call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, &
172  "This sets the default value for the various _2018_ANSWERS parameters.", &
173  default=.false.)
174  call get_param(param_file, mdl, "IDL_HURR_2018_ANSWERS", cs%answers_2018, &
175  "If true, use expressions driving the idealized hurricane test case that recover "//&
176  "the answers from the end of 2018. Otherwise use expressions that are rescalable "//&
177  "and respect rotational symmetry.", default=default_2018_answers)
178 
179  ! The following parameters are model run-time parameters which are used
180  ! and logged elsewhere and so should not be logged here. The default
181  ! value should be consistent with the rest of the model.
182  call get_param(param_file, mdl, "RHO_0", cs%Rho0, &
183  "The mean ocean density used with BOUSSINESQ true to "//&
184  "calculate accelerations and the mass for conservation "//&
185  "properties, or with BOUSSINSEQ false to convert some "//&
186  "parameters from vertical units of m to kg m-2.", &
187  units="kg m-3", default=1035.0, scale=us%kg_m3_to_R, do_not_log=.true.)
188  call get_param(param_file, mdl, "GUST_CONST", cs%gustiness, &
189  "The background gustiness in the winds.", units="Pa", &
190  default=0.0, scale=us%kg_m3_to_R*us%m_s_to_L_T**2*us%L_to_Z, do_not_log=.true.)
191 
192  if (cs%BR_BENCH) then
193  cs%rho_a = 1.2*us%kg_m3_to_R
194  endif
195  dp = cs%pressure_ambient - cs%pressure_central
196  if (cs%answers_2018) then
197  c = cs%max_windspeed / sqrt( us%R_to_kg_m3 * dp )
198  cs%Holland_B = c**2 * us%R_to_kg_m3*cs%rho_a * exp(1.0)
199  else
200  cs%Holland_B = cs%max_windspeed**2 * cs%rho_a * exp(1.0) / dp
201  endif
202  cs%Holland_A = (us%L_to_m*cs%rad_max_wind)**cs%Holland_B
203  cs%Holland_AxBxDP = cs%Holland_A*cs%Holland_B*dp
204 

◆ idealized_hurricane_wind_profile()

subroutine idealized_hurricane::idealized_hurricane_wind_profile ( type(idealized_hurricane_cs), pointer  CS,
type(unit_scale_type), intent(in)  US,
real, intent(in)  absf,
real, intent(in)  YY,
real, intent(in)  XX,
real, intent(in)  UOCN,
real, intent(in)  VOCN,
real, intent(out)  Tx,
real, intent(out)  Ty 
)
private

Calculate the wind speed at a location as a function of time.

Parameters
csContainer for idealized hurricane parameters
[in]usA dimensional unit scaling type
[in]absfInput Coriolis magnitude [T-1 ~> s-1]
[in]yyLocation in m relative to center y [L ~> m]
[in]xxLocation in m relative to center x [L ~> m]
[in]uocnX surface current [L T-1 ~> m s-1]
[in]vocnY surface current [L T-1 ~> m s-1]
[out]txX stress [R L Z T-2 ~> Pa]
[out]tyY stress [R L Z T-2 ~> Pa]

Definition at line 325 of file Idealized_Hurricane.F90.

326  ! Author: Brandon Reichl
327  ! Date: Nov-20-2014
328  ! Aug-14-2018 Generalized for non-SCM configuration
329 
330  ! Input parameters
331  type(idealized_hurricane_CS), pointer :: CS !< Container for idealized hurricane parameters
332  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
333  real, intent(in) :: absf !< Input Coriolis magnitude [T-1 ~> s-1]
334  real, intent(in) :: YY !< Location in m relative to center y [L ~> m]
335  real, intent(in) :: XX !< Location in m relative to center x [L ~> m]
336  real, intent(in) :: UOCN !< X surface current [L T-1 ~> m s-1]
337  real, intent(in) :: VOCN !< Y surface current [L T-1 ~> m s-1]
338  real, intent(out) :: Tx !< X stress [R L Z T-2 ~> Pa]
339  real, intent(out) :: Ty !< Y stress [R L Z T-2 ~> Pa]
340 
341  ! Local variables
342 
343  ! Wind profile terms
344  real :: U10 ! The 10 m wind speed [L T-1 ~> m s-1]
345  real :: radius ! The distance from the hurricane center [L ~> m]
346  real :: radius10 ! 10 times the distance from the hurricane center [L ~> m]
347  real :: radius_km ! The distance from the hurricane center, perhaps in km [L ~> m] or [1000 L ~> km]
348  real :: radiusB
349  real :: tmp ! A temporary variable [R L T-1 ~> kg m-2 s-1]
350  real :: du10 ! The magnitude of the difference between the 10 m wind and the ocean flow [L T-1 ~> m s-1]
351  real :: du ! The difference between the zonal 10 m wind and the zonal ocean flow [L T-1 ~> m s-1]
352  real :: dv ! The difference between the meridional 10 m wind and the zonal ocean flow [L T-1 ~> m s-1]
353  real :: CD
354 
355  !Wind angle variables
356  real :: Alph !< The resulting inflow angle (positive outward)
357  real :: Rstr
358  real :: A0
359  real :: A1
360  real :: P1
361  real :: Adir
362  real :: V_TS ! Meridional hurricane translation speed [L T-1 ~> m s-1]
363  real :: U_TS ! Zonal hurricane translation speed [L T-1 ~> m s-1]
364 
365  ! Implementing Holland (1980) parameteric wind profile
366 
367  radius = sqrt(xx**2 + yy**2)
368 
369  !/ BGR
370  ! rkm - r converted to km for Holland prof.
371  ! used in km due to error, correct implementation should
372  ! not need rkm, but to match winds w/ experiment this must
373  ! be maintained. Causes winds far from storm center to be a
374  ! couple of m/s higher than the correct Holland prof.
375  if (cs%BR_Bench) then
376  radius_km = radius/1000.
377  else
378  ! if not comparing to benchmark, then use correct Holland prof.
379  radius_km = radius
380  endif
381  radiusb = (us%L_to_m*radius)**cs%Holland_B
382 
383  !/
384  ! Calculate U10 in the interior (inside of 10x radius of maximum wind),
385  ! while adjusting U10 to 0 outside of 12x radius of maximum wind.
386  if (cs%answers_2018) then
387  if ( (radius > 0.001*cs%rad_max_wind) .and. (radius < 10.*cs%rad_max_wind) ) then
388  u10 = sqrt(cs%Holland_AxBxDP*exp(-cs%Holland_A/radiusb) / (cs%rho_a*radiusb) + &
389  0.25*(radius_km*absf)**2) - 0.5*radius_km*absf
390  elseif ( (radius > 10.*cs%rad_max_wind) .and. (radius < 15.*cs%rad_max_wind) ) then
391  radius10 = cs%rad_max_wind*10.
392  if (cs%BR_Bench) then
393  radius_km = radius10/1000.
394  else
395  radius_km = radius10
396  endif
397  radiusb = (us%L_to_m*radius10)**cs%Holland_B
398 
399  u10 = (sqrt(cs%Holland_AxBxDp*exp(-cs%Holland_A/radiusb) / (cs%rho_a*radiusb) + &
400  0.25*(radius_km*absf)**2) - 0.5*radius_km*absf) &
401  * (15. - radius/cs%rad_max_wind)/5.
402  else
403  u10 = 0.
404  endif
405  else ! This is mathematically equivalent to that is above but more accurate.
406  if ( (radius > 0.001*cs%rad_max_wind) .and. (radius < 10.*cs%rad_max_wind) ) then
407  tmp = ( 0.5*radius_km*absf) * (cs%rho_a*radiusb)
408  u10 = (cs%Holland_AxBxDP * exp(-cs%Holland_A/radiusb)) / &
409  ( tmp + sqrt(cs%Holland_AxBxDP*exp(-cs%Holland_A/radiusb) * (cs%rho_a*radiusb) + tmp**2) )
410  elseif ( (radius > 10.*cs%rad_max_wind) .and. (radius < 15.*cs%rad_max_wind) ) then
411  radius_km = 10.0 * cs%rad_max_wind
412  if (cs%BR_Bench) radius_km = radius_km/1000.
413  radiusb = (10.0*us%L_to_m*cs%rad_max_wind)**cs%Holland_B
414  tmp = ( 0.5*radius_km*absf) * (cs%rho_a*radiusb)
415  u10 = (3.0 - radius/(5.0*cs%rad_max_wind)) * (cs%Holland_AxBxDp*exp(-cs%Holland_A/radiusb) ) / &
416  ( tmp + sqrt(cs%Holland_AxBxDp*exp(-cs%Holland_A/radiusb) * (cs%rho_a*radiusb) + tmp**2) )
417  else
418  u10 = 0.0
419  endif
420  endif
421 
422  adir = atan2(yy,xx)
423 
424  !\
425 
426  ! Wind angle model following Zhang and Ulhorn (2012)
427  ! ALPH is inflow angle positive outward.
428  rstr = min(10., radius / cs%rad_max_wind)
429  a0 = -0.9*rstr - 0.09*us%L_T_to_m_s*cs%max_windspeed - 14.33
430  a1 = -a0*(0.04*rstr + 0.05*us%L_T_to_m_s*cs%hurr_translation_spd + 0.14)
431  p1 = (6.88*rstr - 9.60*us%L_T_to_m_s*cs%hurr_translation_spd + 85.31) * cs%Deg2Rad
432  alph = a0 - a1*cos(cs%hurr_translation_dir-adir-p1)
433  if ( (radius > 10.*cs%rad_max_wind) .and.&
434  (radius < 15.*cs%rad_max_wind) ) then
435  alph = alph*(15.0 - radius/cs%rad_max_wind)/5.
436  elseif (radius > 15.*cs%rad_max_wind) then
437  alph = 0.0
438  endif
439  alph = alph * cs%Deg2Rad
440 
441  ! Calculate translation speed components
442  u_ts = cs%hurr_translation_spd * 0.5*cos(cs%hurr_translation_dir)
443  v_ts = cs%hurr_translation_spd * 0.5*sin(cs%hurr_translation_dir)
444 
445  ! Set output (relative) winds
446  du = u10*sin(adir-cs%Pi-alph) - uocn + u_ts
447  dv = u10*cos(adir-alph) - vocn + v_ts
448 
449  ! Use a simple drag coefficient as a function of U10 (from Sullivan et al., 2010)
450  du10 = sqrt(du**2+dv**2)
451  if (du10 < 11.0*us%m_s_to_L_T) then
452  cd = 1.2e-3
453  elseif (du10 < 20.0*us%m_s_to_L_T) then
454  if (cs%answers_2018) then
455  cd = (0.49 + 0.065*us%L_T_to_m_s*u10)*1.e-3
456  else
457  cd = (0.49 + 0.065*us%L_T_to_m_s*du10)*1.e-3
458  endif
459  else
460  cd = 1.8e-3
461  endif
462 
463  ! Compute stress vector
464  tx = us%L_to_Z * cs%rho_a * cd * sqrt(du**2 + dv**2) * du
465  ty = us%L_to_Z * cs%rho_a * cd * sqrt(du**2 + dv**2) * dv
466 

◆ scm_idealized_hurricane_wind_forcing()

subroutine, public idealized_hurricane::scm_idealized_hurricane_wind_forcing ( type(surface), intent(in)  sfc_state,
type(mech_forcing), intent(inout)  forces,
type(time_type), intent(in)  day,
type(ocean_grid_type), intent(inout)  G,
type(unit_scale_type), intent(in)  US,
type(idealized_hurricane_cs), pointer  CS 
)

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.

Parameters
[in]sfc_stateSurface state structure
[in,out]forcesA structure with the driving mechanical forces
[in]dayTime in days
[in,out]gGrid structure
[in]usA dimensional unit scaling type
csContainer for SCM parameters

Definition at line 473 of file Idealized_Hurricane.F90.

474  type(surface), intent(in) :: sfc_state !< Surface state structure
475  type(mech_forcing), intent(inout) :: forces !< A structure with the driving mechanical forces
476  type(time_type), intent(in) :: day !< Time in days
477  type(ocean_grid_type), intent(inout) :: G !< Grid structure
478  type(unit_scale_type), intent(in) :: US !< A dimensional unit scaling type
479  type(idealized_hurricane_CS), pointer :: CS !< Container for SCM parameters
480  ! Local variables
481  integer :: i, j, is, ie, js, je, Isq, Ieq, Jsq, Jeq
482  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB
483  real :: pie, Deg2Rad
484  real :: du10 ! The magnitude of the difference between the 10 m wind and the ocean flow [L T-1 ~> m s-1]
485  real :: U10 ! The 10 m wind speed [L T-1 ~> m s-1]
486  real :: A, B, C ! For wind profile expression
487  real :: rad ! The distance from the hurricane center [L ~> m]
488  real :: rkm ! The distance from the hurricane center, sometimes scaled to km [L ~> m] or [1000 L ~> km]
489  real :: f_local ! The local Coriolis parameter [T-1 ~> s-1]
490  real :: xx ! x-position [L ~> m]
491  real :: t0 !for location
492  real :: dP ! The pressure difference across the hurricane [R L2 T-2 ~> Pa]
493  real :: rB
494  real :: Cd ! Air-sea drag coefficient
495  real :: Uocn, Vocn ! Surface ocean velocity components [L T-1 ~> m s-1]
496  real :: dU, dV ! Air-sea differential motion [L T-1 ~> m s-1]
497  !Wind angle variables
498  real :: Alph,Rstr, A0, A1, P1, Adir, transdir
499  real :: V_TS, U_TS ! Components of the translation speed [L T-1 ~> m s-1]
500  logical :: BR_Bench
501  ! Bounds for loops and memory allocation
502  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
503  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
504  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
505  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
506 
507  ! Allocate the forcing arrays, if necessary.
508 
509  call allocate_mech_forcing(g, forces, stress=.true., ustar=.true.)
510  pie = 4.0*atan(1.0) ; deg2rad = pie/180.
511  !/ BR
512  ! Implementing Holland (1980) parameteric wind profile
513  !------------------------------------------------------|
514  br_bench = .true. !true if comparing to LES runs |
515  t0 = 129600. !TC 'eye' crosses (0,0) at 36 hours|
516  transdir = pie !translation direction (-x) |
517  !------------------------------------------------------|
518  dp = cs%pressure_ambient - cs%pressure_central
519  if (cs%answers_2018) then
520  c = cs%max_windspeed / sqrt( us%R_to_kg_m3*dp )
521  b = c**2 * us%R_to_kg_m3*cs%rho_a * exp(1.0)
522  if (br_bench) then ! rho_a reset to value used in generated wind for benchmark test
523  b = c**2 * 1.2 * exp(1.0)
524  endif
525  elseif (br_bench) then ! rho_a reset to value used in generated wind for benchmark test
526  b = (cs%max_windspeed**2 / dp ) * 1.2*us%kg_m3_to_R * exp(1.0)
527  else
528  b = (cs%max_windspeed**2 /dp ) * cs%rho_a * exp(1.0)
529  endif
530 
531  a = (us%L_to_m*cs%rad_max_wind / 1000.)**b
532  f_local = g%CoriolisBu(is,js) ! f=f(x,y) but in the SCM is constant
533  if (br_bench) then
534  ! f reset to value used in generated wind for benchmark test
535  f_local = 5.5659e-05*us%T_to_s
536  endif
537  !/ BR
538  ! Calculate x position as a function of time.
539  xx = us%s_to_T*( t0 - time_type_to_real(day)) * cs%hurr_translation_spd * cos(transdir)
540  rad = sqrt(xx**2 + cs%dy_from_center**2)
541  !/ BR
542  ! rkm - rad converted to km for Holland prof.
543  ! used in km due to error, correct implementation should
544  ! not need rkm, but to match winds w/ experiment this must
545  ! be maintained. Causes winds far from storm center to be a
546  ! couple of m/s higher than the correct Holland prof.
547  if (br_bench) then
548  rkm = rad/1000.
549  rb = (us%L_to_m*rkm)**b
550  else
551  ! if not comparing to benchmark, then use correct Holland prof.
552  rkm = rad
553  rb = (us%L_to_m*rad)**b
554  endif
555  !/ BR
556  ! Calculate U10 in the interior (inside of 10x radius of maximum wind),
557  ! while adjusting U10 to 0 outside of 12x radius of maximum wind.
558  ! Note that rho_a is set to 1.2 following generated wind for experiment
559  if (rad > 0.001*cs%rad_max_wind .AND. rad < 10.*cs%rad_max_wind) then
560  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
561  elseif (rad > 10.*cs%rad_max_wind .AND. rad < 12.*cs%rad_max_wind) then
562  rad=(cs%rad_max_wind)*10.
563  if (br_bench) then
564  rkm = rad/1000.
565  rb = (us%L_to_m*rkm)**b
566  else
567  rkm = rad
568  rb = (us%L_to_m*rad)**b
569  endif
570  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) &
571  * (12. - rad/cs%rad_max_wind)/2.
572  else
573  u10 = 0.
574  endif
575  adir = atan2(cs%dy_from_center,xx)
576 
577  !/ BR
578  ! Wind angle model following Zhang and Ulhorn (2012)
579  ! ALPH is inflow angle positive outward.
580  rstr = min(10., rad / cs%rad_max_wind)
581  a0 = -0.9*rstr - 0.09*us%L_T_to_m_s*cs%max_windspeed - 14.33
582  a1 = -a0 *(0.04*rstr + 0.05*us%L_T_to_m_s*cs%hurr_translation_spd + 0.14)
583  p1 = (6.88*rstr - 9.60*us%L_T_to_m_s*cs%hurr_translation_spd + 85.31)*pie/180.
584  alph = a0 - a1*cos( (transdir - adir ) - p1)
585  if (rad > 10.*cs%rad_max_wind .AND. rad < 12.*cs%rad_max_wind) then
586  alph = alph* (12. - rad/cs%rad_max_wind)/2.
587  elseif (rad > 12.*cs%rad_max_wind) then
588  alph = 0.0
589  endif
590  alph = alph * deg2rad
591  !/BR
592  ! Prepare for wind calculation
593  ! X_TS is component of translation speed added to wind vector
594  ! due to background steering wind.
595  u_ts = cs%hurr_translation_spd*0.5*cos(transdir)
596  v_ts = cs%hurr_translation_spd*0.5*sin(transdir)
597 
598  ! Set the surface wind stresses, in [Pa]. A positive taux
599  ! accelerates the ocean to the (pseudo-)east.
600  ! The i-loop extends to is-1 so that taux can be used later in the
601  ! calculation of ustar - otherwise the lower bound would be Isq.
602  do j=js,je ; do i=is-1,ieq
603  !/BR
604  ! Turn off surface current for stress calculation to be
605  ! consistent with test case.
606  uocn = 0. ! sfc_state%u(I,j)
607  vocn = 0. ! 0.25*( (sfc_state%v(i,J) + sfc_state%v(i+1,J-1)) + &
608  ! (sfc_state%v(i+1,J) + sfc_state%v(i,J-1)) )
609  !/BR
610  ! Wind vector calculated from location/direction (sin/cos flipped b/c
611  ! cyclonic wind is 90 deg. phase shifted from position angle).
612  du = u10*sin(adir-pie-alph) - uocn + u_ts
613  dv = u10*cos(adir-alph) - vocn + v_ts
614  !/----------------------------------------------------|
615  !BR
616  ! Add a simple drag coefficient as a function of U10 |
617  !/----------------------------------------------------|
618  du10 = sqrt(du**2+dv**2)
619  if (du10 < 11.0*us%m_s_to_L_T) then
620  cd = 1.2e-3
621  elseif (du10 < 20.0*us%m_s_to_L_T) then
622  if (cs%answers_2018) then
623  cd = (0.49 + 0.065 * us%L_T_to_m_s*u10 )*0.001
624  else
625  cd = (0.49 + 0.065 * us%L_T_to_m_s*du10 )*0.001
626  endif
627  else
628  cd = 0.0018
629  endif
630  forces%taux(i,j) = cs%rho_a * us%L_to_Z * g%mask2dCu(i,j) * cd*du10*du
631  enddo ; enddo
632  !/BR
633  ! See notes above
634  do j=js-1,jeq ; do i=is,ie
635  uocn = 0. ! 0.25*( (sfc_state%u(I,j) + sfc_state%u(I-1,j+1)) + &
636  ! (sfc_state%u(I-1,j) + sfc_state%u(I,j+1)) )
637  vocn = 0. ! sfc_state%v(i,J)
638  du = u10*sin(adir-pie-alph) - uocn + u_ts
639  dv = u10*cos(adir-alph) - vocn + v_ts
640  du10=sqrt(du**2+dv**2)
641  if (du10 < 11.0*us%m_s_to_L_T) then
642  cd = 1.2e-3
643  elseif (du10 < 20.0*us%m_s_to_L_T) then
644  if (cs%answers_2018) then
645  cd = (0.49 + 0.065 * us%L_T_to_m_s*u10 )*0.001
646  else
647  cd = (0.49 + 0.065 * us%L_T_to_m_s*du10 )*0.001
648  endif
649  else
650  cd = 0.0018
651  endif
652  forces%tauy(i,j) = cs%rho_a * us%L_to_Z * g%mask2dCv(i,j) * cd*du10*dv
653  enddo ; enddo
654  ! Set the surface friction velocity [Z T-1 ~> m s-1]. ustar is always positive.
655  do j=js,je ; do i=is,ie
656  ! This expression can be changed if desired, but need not be.
657  forces%ustar(i,j) = g%mask2dT(i,j) * sqrt(us%L_to_Z * (cs%gustiness/cs%Rho0 + &
658  sqrt(0.5*(forces%taux(i-1,j)**2 + forces%taux(i,j)**2) + &
659  0.5*(forces%tauy(i,j-1)**2 + forces%tauy(i,j)**2))/cs%Rho0))
660  enddo ; enddo
661