28 use mom_time_manager,
only : time_type,
operator(+),
operator(/), time_type_to_real
33 implicit none ;
private
35 #include <MOM_memory.h>
37 public idealized_hurricane_wind_init
39 public idealized_hurricane_wind_forcing
41 public scm_idealized_hurricane_wind_forcing
49 real :: pressure_ambient
50 real :: pressure_central
53 real :: hurr_translation_spd
54 real :: hurr_translation_dir
65 real :: holland_axbxdp
67 logical :: relative_tau
77 real :: dy_from_center
86 #include "version_variable.h"
88 character(len=40) :: mdl =
"idealized_hurricane"
93 subroutine idealized_hurricane_wind_init(Time, G, param_file, CS)
99 intent(in) :: param_file
106 #include "version_variable.h"
108 if (
associated(cs))
then
109 call mom_error(fatal,
"idealized_hurricane_wind_init called "// &
110 "with an associated control structure.")
116 cs%pi = 4.0*atan(1.0)
117 cs%Deg2Rad = cs%pi/180.
123 call get_param(param_file, mdl,
"IDL_HURR_RHO_AIR", cs%rho_a, &
124 "Air density used to compute the idealized hurricane "//&
125 "wind profile.", units=
'kg/m3', default=1.2)
126 call get_param(param_file, mdl,
"IDL_HURR_AMBIENT_PRESSURE", &
127 cs%pressure_ambient,
"Ambient pressure used in the "//&
128 "idealized hurricane wind profile.", units=
'Pa', &
130 call get_param(param_file, mdl,
"IDL_HURR_CENTRAL_PRESSURE", &
131 cs%pressure_central,
"Central pressure used in the "//&
132 "idealized hurricane wind profile.", units=
'Pa', &
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', &
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.)
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)
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', &
148 cs%hurr_translation_dir = cs%hurr_translation_dir * cs%Deg2Rad
149 call get_param(param_file, mdl,
"IDL_HURR_X0", cs%Hurr_cen_X0, &
150 "Idealized Hurricane initial X position", &
151 units=
'm', default=0.)
152 call get_param(param_file, mdl,
"IDL_HURR_Y0", cs%Hurr_cen_Y0, &
153 "Idealized Hurricane initial Y position", &
154 units=
'm', default=0.)
155 call get_param(param_file, mdl,
"IDL_HURR_TAU_CURR_REL", cs%relative_tau, &
156 "Current relative stress switch "//&
157 "used in the idealized hurricane wind profile.", &
158 units=
'', default=.false.)
161 call get_param(param_file, mdl,
"IDL_HURR_SCM_BR_BENCH", cs%BR_BENCH, &
162 "Single column mode benchmark case switch, which is "// &
163 "invoking a modification (bug) in the wind profile meant to "//&
164 "reproduce a previous implementation.", units=
'', default=.false.)
165 call get_param(param_file, mdl,
"IDL_HURR_SCM", cs%SCM_MODE, &
166 "Single Column mode switch "//&
167 "used in the SCM idealized hurricane wind profile.", &
168 units=
'', default=.false.)
169 call get_param(param_file, mdl,
"IDL_HURR_SCM_LOCY", cs%DY_from_center, &
170 "Y distance of station used in the SCM idealized hurricane "//&
171 "wind profile.", units=
'm', default=50.e3)
176 call get_param(param_file, mdl,
"RHO_0", cs%Rho0, &
177 "The mean ocean density used with BOUSSINESQ true to "//&
178 "calculate accelerations and the mass for conservation "//&
179 "properties, or with BOUSSINSEQ false to convert some "//&
180 "parameters from vertical units of m to kg m-2.", &
181 units=
"kg m-3", default=1035.0, do_not_log=.true.)
182 call get_param(param_file, mdl,
"GUST_CONST", cs%gustiness, &
183 "The background gustiness in the winds.", units=
"Pa", &
184 default=0.00, do_not_log=.true.)
187 if (cs%BR_BENCH)
then
190 dp = cs%pressure_ambient - cs%pressure_central
191 c = cs%max_windspeed / sqrt( dp )
192 cs%Holland_B = c**2 * cs%rho_a * exp(1.0)
193 cs%Holland_A = (cs%rad_max_wind)**cs%Holland_B
194 cs%Holland_AxBxDP = cs%Holland_A*cs%Holland_B*dp
197 end subroutine idealized_hurricane_wind_init
200 subroutine idealized_hurricane_wind_forcing(state, forces, day, G, US, CS)
203 type(time_type),
intent(in) :: day
209 integer :: i, j, is, ie, js, je, isq, ieq, jsq, jeq
210 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
225 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
226 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
227 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
228 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
231 call allocate_mech_forcing(g, forces, stress=.true., ustar=.true.)
233 if (cs%relative_tau)
then
240 xc = cs%Hurr_cen_X0 + (time_type_to_real(day)*cs%hurr_translation_spd*&
241 cos(cs%hurr_translation_dir))
242 yc = cs%Hurr_cen_Y0 + (time_type_to_real(day)*cs%hurr_translation_spd*&
243 sin(cs%hurr_translation_dir))
246 if (cs%BR_Bench)
then
258 uocn = state%u(i,j)*rel_tau_fac
259 vocn = 0.25*(state%v(i,j)+state%v(i+1,j-1)&
260 +state%v(i+1,j)+state%v(i,j-1))*rel_tau_fac
261 f = abs(0.5*us%s_to_T*(g%CoriolisBu(i,j)+g%CoriolisBu(i,j-1)))*fbench_fac + fbench
263 if (cs%SCM_mode)
then
264 yy = yc + cs%dy_from_center
267 lat = g%geoLatCu(i,j)*1000.
268 lon = g%geoLonCu(i,j)*1000.
272 call idealized_hurricane_wind_profile(&
273 cs,f,yy,xx,uocn,vocn,tx,ty)
274 forces%taux(i,j) = g%mask2dCu(i,j) * tx
280 uocn = 0.25*(state%u(i,j)+state%u(i-1,j+1)&
281 +state%u(i-1,j)+state%u(i,j+1))*rel_tau_fac
282 vocn = state%v(i,j)*rel_tau_fac
283 f = abs(0.5*us%s_to_T*(g%CoriolisBu(i-1,j)+g%CoriolisBu(i,j)))*fbench_fac + fbench
285 if (cs%SCM_mode)
then
286 yy = yc + cs%dy_from_center
289 lat = g%geoLatCv(i,j)*1000.
290 lon = g%geoLonCv(i,j)*1000.
294 call idealized_hurricane_wind_profile(cs, f, yy, xx, uocn, vocn, tx, ty)
295 forces%tauy(i,j) = g%mask2dCv(i,j) * ty
303 forces%ustar(i,j) = us%m_to_Z*us%T_to_s * g%mask2dT(i,j) * sqrt(cs%gustiness/cs%Rho0 + &
304 sqrt(0.5*(forces%taux(i-1,j)**2 + forces%taux(i,j)**2) + &
305 0.5*(forces%tauy(i,j-1)**2 + forces%tauy(i,j)**2))/cs%Rho0)
310 end subroutine idealized_hurricane_wind_forcing
313 subroutine idealized_hurricane_wind_profile(CS, absf, YY, XX, UOCN, VOCN, Tx, Ty)
321 real,
intent(in) :: absf
322 real,
intent(in) :: YY
323 real,
intent(in) :: XX
324 real,
intent(in) :: UOCN
325 real,
intent(in) :: VOCN
326 real,
intent(out) :: Tx
327 real,
intent(out) :: Ty
355 radius = sqrt(xx**2 + yy**2)
363 if (cs%BR_Bench)
then
364 radius_km = radius/1000.
369 radiusb = (radius)**cs%Holland_B
374 if ( (radius/cs%rad_max_wind .gt. 0.001) .and. &
375 (radius/cs%rad_max_wind .lt. 10.) )
then
376 u10 = sqrt(cs%Holland_AxBxDP*exp(-cs%Holland_A/radiusb)/(cs%rho_A*radiusb)&
377 +0.25*(radius_km*absf)**2) - 0.5*radius_km*absf
378 elseif ( (radius/cs%rad_max_wind .gt. 10.) .and. &
379 (radius/cs%rad_max_wind .lt. 15.) )
then
381 radius10 = cs%rad_max_wind*10.
383 if (cs%BR_Bench)
then
384 radius_km = radius10/1000.
388 radiusb=radius10**cs%Holland_B
390 u10 = (sqrt(cs%Holland_AxBxDp*exp(-cs%Holland_A/radiusb)/(cs%rho_A*radiusb)&
391 +0.25*(radius_km*absf)**2)-0.5*radius_km*absf) &
392 * (15.-radius/cs%rad_max_wind)/5.
401 rstr = min(10.,radius / cs%rad_max_wind)
402 a0 = -0.9*rstr - 0.09*cs%max_windspeed - 14.33
403 a1 = -a0*(0.04*rstr + 0.05*cs%Hurr_translation_spd + 0.14)
404 p1 = (6.88*rstr - 9.60*cs%Hurr_translation_spd + 85.31) * cs%Deg2Rad
405 alph = a0 - a1*cos(cs%hurr_translation_dir-adir-p1)
406 if ( (radius/cs%rad_max_wind.gt.10.) .and.&
407 (radius/cs%rad_max_wind.lt.15.) )
then
408 alph = alph*(15.0-radius/cs%rad_max_wind)/5.
409 elseif (radius/cs%rad_max_wind.gt.15.)
then
412 alph = alph * cs%Deg2Rad
415 u_ts = cs%hurr_translation_spd/2.*cos(cs%hurr_translation_dir)
416 v_ts = cs%hurr_translation_spd/2.*sin(cs%hurr_translation_dir)
419 du = u10*sin(adir-cs%Pi-alph) - uocn + u_ts
420 dv = u10*cos(adir-alph) - vocn + v_ts
423 du10 = sqrt(du**2+dv**2)
424 if (du10.lt.11.)
then
426 elseif (du10.lt.20.0)
then
427 cd = (0.49 + 0.065*u10)*1.e-3
433 tx = cs%rho_A * cd * sqrt(du**2 + dv**2) * du
434 ty = cs%rho_A * cd * sqrt(du**2 + dv**2) * dv
437 end subroutine idealized_hurricane_wind_profile
443 subroutine scm_idealized_hurricane_wind_forcing(state, forces, day, G, US, CS)
446 type(time_type),
intent(in) :: day
451 integer :: i, j, is, ie, js, je, isq, ieq, jsq, jeq
452 integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
454 real :: u10, a, b, c, r, f, du10, rkm
461 real :: alph,rstr, a0, a1, p1, adir, transdir, v_ts, u_ts
464 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
465 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
466 isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
467 isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
471 call allocate_mech_forcing(g, forces, stress=.true., ustar=.true.)
472 pie = 4.0*atan(1.0) ; deg2rad = pie/180.
480 dp = cs%pressure_ambient - cs%pressure_central
481 c = cs%max_windspeed / sqrt( dp )
482 b = c**2 * cs%rho_a * exp(1.0)
485 b = c**2 * 1.2 * exp(1.0)
487 a = (cs%rad_max_wind/1000.)**b
488 f = us%s_to_T*g%CoriolisBu(is,js)
495 xx = ( t0 - time_type_to_real(day)) * cs%hurr_translation_spd * cos(transdir)
496 r = sqrt(xx**2 + cs%DY_from_center**2)
515 if (r/cs%rad_max_wind > 0.001 .AND. r/cs%rad_max_wind < 10.)
then
516 u10 = sqrt( a*b*dp*exp(-a/rb)/(1.2*rb) + 0.25*(rkm*f)**2 ) - 0.5*rkm*f
517 elseif (r/cs%rad_max_wind > 10. .AND. r/cs%rad_max_wind < 12.)
then
518 r=cs%rad_max_wind*10.
526 u10 = ( sqrt( a*b*dp*exp(-a/rb)/(1.2*rb) + 0.25*(rkm*f)**2 ) - 0.5*rkm*f) &
527 * (12. - r/cs%rad_max_wind)/2.
531 adir = atan2(cs%DY_from_center,xx)
536 rstr = min(10.,r / cs%rad_max_wind)
537 a0 = -0.9*rstr -0.09*cs%max_windspeed - 14.33
538 a1 = -a0 *(0.04*rstr +0.05*cs%hurr_translation_spd+0.14)
539 p1 = (6.88*rstr -9.60*cs%hurr_translation_spd+85.31)*pie/180.
540 alph = a0 - a1*cos( (transdir - adir ) - p1)
541 if (r/cs%rad_max_wind > 10. .AND. r/cs%rad_max_wind < 12.)
then
542 alph = alph* (12. - r/cs%rad_max_wind)/2.
543 elseif (r/cs%rad_max_wind > 12.)
then
546 alph = alph * deg2rad
551 u_ts = cs%hurr_translation_spd/2.*cos(transdir)
552 v_ts = cs%hurr_translation_spd/2.*sin(transdir)
558 do j=js,je ;
do i=is-1,ieq
568 du = u10*sin(adir-pie-alph) - uocn + u_ts
569 dv = u10*cos(adir-alph) - vocn + v_ts
574 du10=sqrt(du**2+dv**2)
577 elseif (du10 < 20.)
then
578 cd = (0.49 + 0.065 * u10 )*0.001
582 forces%taux(i,j) = cs%rho_a * g%mask2dCu(i,j) * cd*sqrt(du**2+dv**2)*du
586 do j=js-1,jeq ;
do i=is,ie
590 du = u10*sin(adir-pie-alph) - uocn + u_ts
591 dv = u10*cos(adir-alph) - vocn + v_ts
592 du10=sqrt(du**2+dv**2)
595 elseif (du10 < 20.)
then
596 cd = (0.49 + 0.065 * u10 )*0.001
600 forces%tauy(i,j) = cs%rho_a * g%mask2dCv(i,j) * cd*du10*dv
603 do j=js,je ;
do i=is,ie
605 forces%ustar(i,j) = us%m_to_Z*us%T_to_s * g%mask2dT(i,j) * sqrt(cs%gustiness/cs%Rho0 + &
606 sqrt(0.5*(forces%taux(i-1,j)**2 + forces%taux(i,j)**2) + &
607 0.5*(forces%tauy(i,j-1)**2 + forces%tauy(i,j)**2))/cs%Rho0)
610 end subroutine scm_idealized_hurricane_wind_forcing