!================================================================= !================================================================= !================================================================= ! ===== ===== ! ===== MODULE vortex ===== ! ===== ===== !================================================================= !================================================================= !================================================================= !================================================================= ! This module allows one to create an asymmetric hurricane vortex ! object from National Hurricane Center forecast advisories, then ! compute wind and pressure fields suitable for forcing a storm ! surge forecast model. The shape and intensity of the vortex can ! be continuously adjusted during the forecast period, as it moves ! through the model domain. ! ! Revision history: ! Date Programmer Description of change ! ---- ---------- --------------------- ! 05/23/06 Craig Mattocks, UNC-CEP Wrote original code ! 06/30/06 Cristina Forbes, UNC-CEP Tested in ADCIRC model ! 08/25/06 Craig Mattocks, UNC-CEP Applied wind reduction: ! top of SFC layer -> SFC ! in subroutine uvp. ! 09/12/06 Craig Mattocks, UNC-CEP Subtracted translational ! wind speed from Vmax in ! nws9get. ! 05/12/09 Robert Weaver, UNC-IMS ! Modified damping of translational ! velocity based on ratio of V/Vmax ! 05/19/09 Cristina Forbes, UNC-IMS ! Implemented old bug fix: units conversion ! 05/22/09 Cristina Forbes, UNC-IMS ! Reverted damping of translational ! velocity back to the original formulation ! developed by Craig Mattocks ! 05/26/09 Cristina Forbes, UNC-IMS ! Fixed implementation of Rick Luettich ! V/Vmax tapering formulation (as in NWS=8) ! for future experimentation - not activated ! 06/2009 Robert Weaver, UNC-IMS ! changed the gradient wind formula to ! use Vmax instead of pressure difference ! to enable easier manipulation of storm ! characteristics ! 07/2009 Robert Weaver, UNC-IMS ! using modules from this code, ! wrote a preprocessor to compute Rmax and B ! prior to running adcirc, and write out an ! input file like the ATCF file with these ! variables to make it easier to manipulate ! the storm parameters. Changed this code ! to read in that input file and use the values ! found there ! 02/2013 Jie Gao, UNC-IMS ! wrote subroutines and functions for ! the Generalized Asymmetric Holland Model (GAHM) ! !================================================================= MODULE vortex USE SIZES, ONLY : sz USE GLOBAL, ONLY : sphericalDistance, deg2rad, rad2deg, & omega, mb2pa, rhoAir, kt2ms, nm2m, m2nm, & ms2kt, Rearth, windReduction, one2ten IMPLICIT NONE SAVE PUBLIC :: newVortex, calcRmaxes, fitRmaxes, & Rmw, uvp, uvtrans, fang, & latlon2xy, xy2latlon, & getShapeParameter , & setShapeParameter , & getLatestRmax, & getLatestAngle , & getUseQuadrantVr , getRmaxes, & setUseQuadrantVr , setRmaxes, & setIsotachWindSpeed , & getShapeParameters , getPhiFactors, & setIsotachWindSpeeds, setIsotachRadii, & setVortex, spInterp,interpR, & setUseVmaxesBL, & getVmaxesBL, setVmaxesBL, & fitRmaxes4,calcRmaxesFull, & uvpr,newVortexFull, & Rmaxes4, QuadFlag4, QuadIr4, Bs4, VmBL4 !jgfdebug PRIVATE INTEGER, PARAMETER :: nQuads = 4 ! Number of quadrants for ! which wind radii are ! provided INTEGER, PARAMETER :: nPoints = nQuads+2 ! Number of (theta, Rmax) ! points for curve fit REAL(sz), DIMENSION(nPoints) :: Rmaxes ! Radius of maximum winds REAL(sz), DIMENSION(nPoints,4) :: Rmaxes4 ! (nautical miles) REAL(sz) :: Pn ! Ambient surface pressure (mb) REAL(sz) :: Pc ! Surface pressure at center of ! storm (mb) REAL(sz) :: cLat ! Latitude of storm center ! (degrees north) REAL(sz) :: cLon ! Longitude of storm center ! (degrees east ) REAL(sz) :: Vmax ! Max sustained wind velocity ! in storm (knots) REAL(sz) :: B ! Exponential shape parameter REAL(sz) :: corio ! Coriolis force (1/s) REAL(sz) :: Vr ! Velocity @ wind radii (knots) REAL(sz) :: Phi REAL(sz), DIMENSION(nPoints) :: Bs REAL(sz), DIMENSION(nPoints) :: Phis ! Correction factor to B and Vh REAL(sz), DIMENSION(nPoints) :: VmBL REAL(sz), DIMENSION(nPoints,4) :: Bs4 REAL(sz), DIMENSION(nPoints,4) :: Phis4 ! Correction factor to B and Vh REAL(sz), DIMENSION(nPoints,4) :: VmBL4 INTEGER, DIMENSION(nPoints,4) :: QuadFlag4 REAL(sz), DIMENSION(nPoints,4) :: QuadIr4 REAL(sz), DIMENSION(nQuads) :: VrQuadrant REAL(sz), DIMENSION(nQuads) :: radius ! Wind radii - the distance ! winds of velocity Vr extend ! outward from center of storm ! (nautical miles) INTEGER :: quad ! Quadrant counter REAL(sz) :: latestRmax ! most recently calculated ! value of fitted rmax REAL(sz) :: latestAngle ! angle of the most recently ! calculated node w.r.t. the ! storm location LOGICAL :: useQuadrantVr LOGICAL :: useVmaxesBL CONTAINS !================= ! Vortex functions !================= !================================================================= ! Create a new Vortex object. ! ! On input: ! Pn Ambient surface pressure (mb) ! Pc Surface pressure at center of storm (mb) ! cLat Latitude of storm center (degrees north) ! cLon Longitude of storm center (degrees east ) ! Vmax Max sustained wind velocity in storm (knots) ! ! On output: ! A new vortex is created with essential parameters calculated. !================================================================= SUBROUTINE newVortex(pinf, p0, lat, lon, vm) REAL(sz), INTENT(IN) :: Pinf REAL(sz), INTENT(IN) :: P0 REAL(sz), INTENT(IN) :: lat REAL(sz), INTENT(IN) :: lon REAL(sz), INTENT(IN) :: vm ! set instance variables Pn = pinf Pc = p0 cLat = lat cLon = lon Vmax = vm ! evaluate basic physical params corio = 2.0d0 * omega * SIN(deg2rad*cLat) B = (Vmax*kt2ms)**2 * RhoAir * EXP(1.d0) & / ((Pn - Pc) * mb2pa) B = MAX( MIN(B,2.5d0), 1.0d0) ! limit B to range 1.0->2.5 ! added for compatibility of calcRmaxes to use with simplified nws20 Bs(1:6)=B VmBL(1:6)=Vmax END SUBROUTINE newVortex ! A new vortex is created for the full gradient wind balance SUBROUTINE newVortexFull(pinf, p0, lat, lon, vm) REAL(sz), INTENT(IN) :: Pinf REAL(sz), INTENT(IN) :: P0 REAL(sz), INTENT(IN) :: lat REAL(sz), INTENT(IN) :: lon REAL(sz), INTENT(IN) :: vm ! set instance variables Pn = pinf Pc = p0 cLat = lat cLon = lon Vmax = vm ! evaluate basic physical params corio = 2.0d0 * omega * SIN(deg2rad*cLat) B = (Vmax*kt2ms)**2 * RhoAir * EXP(1.d0) & / ((Pn - Pc) * mb2pa) Phi=1.d0 Bs(1:6)=B Phis(1:6)=Phi VmBL(1:6)=Vmax ! Jie 2013.01 ! B = MAX( MIN(B,2.5d0), 1.0d0) ! limit B to range 1.0->2.5 END SUBROUTINE newVortexFull !================================================================= ! Set basic parameter for a new Vortex object. ! ! On input: ! Pn Ambient surface pressure (mb) ! Pc Surface pressure at center of storm (mb) ! cLat Latitude of storm center (degrees north) ! cLon Longitude of storm center (degrees east ) ! ! On output: ! Aim is to define Pn, Pc, and corio !================================================================= SUBROUTINE setVortex(pinf, p0, lat, lon) REAL(sz), INTENT(IN) :: Pinf REAL(sz), INTENT(IN) :: P0 REAL(sz), INTENT(IN) :: lat REAL(sz), INTENT(IN) :: lon ! set instance variables Pn = pinf Pc = p0 cLat = lat cLon = lon ! evaluate basic physical params corio = 2.0d0 * omega * SIN(deg2rad*cLat) END SUBROUTINE setVortex !=============================================================== ! Calculate the radius of maximum winds for all storm quadrants. ! ! On input: ! none ! ! On output: ! Rmax radius of maximum winds (nm) in all quadrants, plus ! 2 extra values to tie down circular periodicity ! Jie 2014.07 Modified with quadrant-varying VmBL, which not only ! works for nws19 but for the simplified nws20 !=============================================================== SUBROUTINE calcRmaxes() ! REAL(sz), EXTERNAL :: func REAL(sz) :: root ! Radius of maximum winds REAL(sz), PARAMETER :: innerRadius = 1.d0 REAL(sz), PARAMETER :: outerRadius = 400.d0 REAL(sz), PARAMETER :: accuracy = 0.0001d0 REAL(sz), PARAMETER :: zoom = 0.01d0 INTEGER , PARAMETER :: itmax = 3 REAL(sz) :: r1,r2, r3,r4, dr INTEGER :: n, iter, i REAL(sz) :: vicinity !----------------------------- ! Loop over quadrants of storm !----------------------------- DO n = 1, nQuads ! set B and Vmax values for each quadrant ! for nws19, B and Vmax are constant ! for simplified nws20, B is constant, while Vmax is not B = Bs(n+1) Vmax = VmBL(n+1) quad = n root = -1.d0 r1 = innerRadius r2 = outerRadius dr = 1.d0 DO iter = 1, itmax root = findRoot(VhWithCori, r1,r2, dr, r3,r4) r1 = r3 r2 = r4 dr = dr * zoom END DO ! determine if Rmax is actually in the vicinity of the ! isotach radius that we are using to solve for Rmax, ! and if so, take another shot at finding the ! Rmax using the gradient wind balance that neglects ! coriolis (and is appropriate in the vicinity of Rmax) vicinity = abs(root-radius(quad))/root if ( (root.lt.0.d0).or.(vicinity.le.0.1d0)) then r1 = innerRadius r2 = outerRadius dr = 1.d0 DO iter = 1, itmax root = findRoot(VhNoCori, r1,r2, dr, r3,r4) r1 = r3 r2 = r4 dr = dr * zoom END DO endif Rmaxes(n+1) = root END DO END SUBROUTINE calcRmaxes !=============================================================== ! Calculate the radius of maximum winds for all storm quadrants. ! Solving the full gradient wind equation without the assumption ! of cyclostrohpic balance. ! ! On input: ! none ! ! On output: ! Rmax radius of maximum winds (nm) in all quadrants, plus ! 2 extra values to tie down circular periodicity ! ! Jie 2013.02 added looping procedures to calculate Bs and Phis !=============================================================== SUBROUTINE calcRmaxesFull() ! REAL(sz), EXTERNAL :: func REAL(sz) :: root ! Radius of maximum winds REAL(sz), PARAMETER :: innerRadius = 1.d0 REAL(sz), PARAMETER :: outerRadius = 500.d0 REAL(sz), PARAMETER :: accuracy = 0.0001d0 REAL(sz), PARAMETER :: zoom = 0.01d0 INTEGER , PARAMETER :: itmax = 3 REAL(sz) :: r1,r2, r3,r4, dr INTEGER :: n, iter, i, norootflg REAL(sz) :: vicinity REAL(sz) :: B_new, B_new1 REAL(sz) :: Phi_new INTEGER , PARAMETER :: cont = 400 ! Max # of iterations INTEGER :: icont,ibcont ! iteration counter 211 FORMAT(a7,x,i2,x,a38) !----------------------------- ! Loop over quadrants of storm !----------------------------- DO n = 1, nQuads norootflg=0 ! initialize B and Phi values for each quadrant B = Bs(n+1) Phi = Phis(n+1) Vmax = VmBL(n+1) ! Loop the root-solving process to converge B, for in the ! new wind formulation, B is a function of Rmax, Vmax, f, and Phi DO icont = 1, cont ! logical expre. is at the end to exit the loop norootflg=0 quad = n root = -1.d0 r1 = innerRadius r2 = outerRadius dr = 1.d0 DO iter = 1, itmax root = findRoot(VhWithCoriFull, r1,r2, dr, r3,r4) r1 = r3 r2 = r4 dr = dr * zoom END DO ! avoid invalid B value when root is not found if (root.lt.0.d0) then ! r1 = innerRadius ! r2 = outerRadius ! dr = 1.d0 ! DO iter = 1, itmax ! root = findRoot(VhNoCori, r1,r2, dr, r3,r4) ! r1 = r3 ! r2 = r4 ! dr = dr * zoom ! END DO root=1.0*radius(quad) norootflg=1 endif Rmaxes(n+1) = root ! Jie 2013.02 ! determine if B converges, if yes, break loop and assign ! values to Rmaxes, if not, continue the loop to re-calculate ! root and re-evaluate Bs Phi_new = 1+Vmax*kt2ms*root*nm2m*corio/ & (B*((Vmax*kt2ms)**2+Vmax*kt2ms*root*nm2m*corio)) B_new = ((Vmax*kt2ms)**2+Vmax*kt2ms*root*nm2m*corio) & *RhoAir*EXP(Phi_new)/(Phi_new*(Pn-Pc)* mb2pa) do ibcont = 1, cont B_new1 = B_new Phi_new = 1+Vmax*kt2ms*root*nm2m*corio/ & (B_new*((Vmax*kt2ms)**2+Vmax*kt2ms*root*nm2m*corio)) B_new = ((Vmax*kt2ms)**2+Vmax*kt2ms*root*nm2m*corio) & *RhoAir*EXP(Phi_new)/(Phi_new*(Pn-Pc)* mb2pa) if (abs(B_new-B_new1).le.0.01d0) then EXIT endif end do ! debug with aswip ! if (ibcont.ge.cont) then ! write(1111,211) ! & "iquad=",n, ! & "B_new did not fully converge, procede" ! endif ! end debug with aswip if (abs(B-B_new).le.0.01d0) then EXIT endif ! update B and Phi for next iteration ! warning: modifications made here also affect other subroutines B=B_new Phi=Phi_new END DO !icont = 1, cont ! update to the latest values for aswip output Bs(n+1)=B_new Phis(n+1)=Phi_new ! debug with aswip ! if (icont.ge.cont) then ! write(1111,211) ! & "iquad=",n, ! & "B did not fully converge, procede" ! endif ! end debug with aswip ! determine if Rmax is actually in the vicinity of the ! isotach radius that we are using to solve for Rmax, ! and if so, take another shot at finding the ! Rmax using the gradient wind equation that neglects ! coriolis (and is appropriate in the vicinity of Rmax) ! Jie 2013.01 !vicinity = abs(root-radius(quad))/root if (norootflg==1) then write(*,*) & "iquad=",n, & "No root found, return dist. to Isotach" endif END DO !n = 1, nQuads END SUBROUTINE calcRmaxesFull !==================================================== ! External function f(x) = 0 for which a root is ! sought using Brent's root-finding method. ! ! On input: ! x iterative values which converge to root ! ! On output: ! func f(x) ! ! Internal parameters: ! vortex instance variables via accessor functions ! ! Jie 2013.02 Modified to use the full gradient wind eq. !==================================================== FUNCTION VhWithCoriFull(testRmax) REAL(sz), INTENT(IN) :: testRmax REAL(sz) :: VhWithCoriFull REAL(sz) :: thisVr ! the radial wind speed we've been given REAL(sz) :: Vh !------------------------- ! func(x = Rmax) = Vh - Vr !------------------------- if ( getUseQuadrantVr().eqv..true. ) then thisVr = VrQuadrant(quad) else thisVr = Vr endif ! Jie 2013.02 Vh = ms2kt * ( SQRT( & ((Vmax*kt2ms)**2.d0+Vmax*kt2ms*testRmax*nm2m*corio) & * (testRmax/radius(quad))**B & * EXP(Phi*(1.0d0-(testRmax/radius(quad))**B)) & + (nm2m*radius(quad)*corio/2.d0)**2.d0) & - nm2m*radius(quad)*corio/2.d0 ) VhWithCoriFull = Vh - thisVr END FUNCTION VhWithCoriFull !=============================================================== !==================================================== ! External function f(x) = 0 for which a root is ! sought using Brent's root-finding method. ! ! On input: ! x iterative values which converge to root ! ! On output: ! func f(x) ! ! Internal parameters: ! vortex instance variables via accessor functions !==================================================== FUNCTION VhWithCori(testRmax) REAL(sz), INTENT(IN) :: testRmax REAL(sz) :: VhWithCori REAL(sz) :: thisVr ! the radial wind speed we've been given REAL(sz) :: Vh !------------------------- ! func(x = Rmax) = Vh - Vr !------------------------- if ( getUseQuadrantVr().eqv..true. ) then thisVr = VrQuadrant(quad) else thisVr = Vr endif Vh = ms2kt * ( SQRT( & (Vmax*kt2ms)**2.d0 * (testRmax/radius(quad))**B & * EXP(1.0d0-(testRmax/radius(quad))**B) & + (nm2m*radius(quad)*corio/2.d0)**2.d0) & - nm2m*radius(quad)*corio/2.d0 ) VhWithCori = Vh - thisVr END FUNCTION VhWithCori !=============================================================== !=============================================================== REAL(sz) FUNCTION VhNoCori(testRmax) REAL(sz), INTENT(IN) :: testRmax REAL(sz) :: thisVr ! the radial wind speed we've been given if ( getUseQuadrantVr().eqv..true. ) then thisVr = VrQuadrant(quad) else thisVr = Vr endif VhNoCori = ABS(ms2kt * & SQRT( (Vmax*kt2ms)**2.0d0 * (testRmax/radius(quad))**B & * EXP(1 - (testRmax/radius(quad))**B))) & - thisVr !=============================================================== END FUNCTION VhNoCori !=============================================================== !============================================================== ! Use brute-force marching to find a root the interval [x1,x2]. ! ! On input: ! func function f(x)=0 for which root is sough ! x1 left side of interval ! x2 right side of interval ! dx x increment for march ! ! On output: ! a left side of interval that brackets the root ! b right side of interval that brackets the root ! findRoot root returned !============================================================== FUNCTION findRoot(func, x1, x2, dx, a,b) REAL(sz), EXTERNAL :: func REAL(sz), INTENT(IN) :: x1, x2 ! Search interval [x1,x2] REAL(sz), INTENT(IN) :: dx ! Marching increment REAL(sz), INTENT(OUT) :: a, b ! x values that bracket root REAL(sz) :: findRoot ! The root found INTEGER , PARAMETER :: itmax = 400 ! Max # of iterations INTEGER :: iter ! iteration counter REAL(sz) :: fa,fb ! function values f(x) ! ! Initialize left side of interval a = x1 fa = func(a) ! ! March along interval until root is found ! or solution diverges. findRoot = a DO iter = 1,itmax b = x1 + iter * dx fb = func(b) ! Check progress IF ((fa*fb < 0.0d0) .OR. (ABS(fb) > ABS(fa))) THEN ! Assign root IF (ABS(fb) > ABS(fa)) THEN findRoot = a ELSE findRoot = b END IF EXIT END IF ! ! Move right search interval values to left side ! for next iteration. a = b fa = fb END DO IF (iter >= itmax) THEN PRINT *, & "FUNCTION findRoot: exceeded max # of iterations" findRoot = -99999.0 END IF END FUNCTION findRoot !================================================================= ! Calculate (u,v) wind components and surface pressure from an ! asymmetric hurricane wind model. ! ! On input: ! lat Latitude of nodal point (degrees north) ! lon Longitude of nodal point (degrees east ) ! uTrans x component of translational velocity (kts) ! vTrans y component of translational velocity (kts) ! ! On output: ! u x component of wind velocity at nodal point (m/s) ! v y component of wind velocity at nodal point (m/s) ! p Surface pressure at nodal point (Pa) ! ! Internal parameters: ! dampRadii How far out (# of Rmax radii) to begin damping ! out the translational velocity ! ! Note: ! Subroutine directly accesses global class instance variables !================================================================= SUBROUTINE uvp(lat,lon, uTrans,vTrans, u,v, p) REAL(sz), INTENT(IN) :: lat REAL(sz), INTENT(IN) :: lon REAL(sz), INTENT(IN) :: uTrans REAL(sz), INTENT(IN) :: vTrans REAL(sz), INTENT(OUT) :: u REAL(sz), INTENT(OUT) :: v REAL(sz), INTENT(OUT) :: p REAL(sz) :: TransSpdX !NWS8-style translation speed REAL(sz) :: TransSpdY !NWS8-style translation speed REAL(sz) :: avgLat REAL(sz) :: dx REAL(sz) :: dy REAL(sz) :: dist REAL(sz) :: rmx REAL(sz) :: angle REAL(sz) :: speed REAL(sz) :: maxspeed REAL(sz) :: uf REAL(sz) :: vf REAL(sz) :: damp REAL(sz) :: percentCoriolis REAL(sz) :: speedAtRmax REAL(sz) :: vmaxFactor INTEGER :: i !------------------------------------------------------ ! Calculate distance and angle between eye of hurricane ! and input nodal point !------------------------------------------------------ dx = deg2rad * Rearth * (lon - cLon) * COS(deg2rad*cLat) dy = deg2rad * Rearth * (lat - cLat) dist = SQRT(dx*dx + dy*dy) !---------------------------------------- ! Handle special case at eye of hurricane !---------------------------------------- ! in eye velocity is zero not translational velocity !---------------------------------------- IF (dist < 1.d0) THEN u = 0.0d0 v = 0.0d0 p = Pc * mb2pa RETURN END IF dist = m2nm * dist angle = 360.0d0 + rad2deg * ATAN2(dx,dy) IF (angle > 360.d0) angle = angle - 360.d0 latestAngle = angle rmx = Rmw(angle) latestRmax = rmx !--------------------------------------------------- ! Compute (u,v) wind velocity components from the ! asymmetric hurricane vortex. ! ! Note: the vortex winds are valid at the top of the ! surface layer, so reduce the winds to the surface. ! Also convert the winds from max sustained 1-minute ! averages to 10-minute averages for the storm surge ! model. !--------------------------------------------------- percentCoriolis=1.0d0 speed = SQRT( & (Vmax*kt2ms)**2.d0 * (rmx/dist)**B & * EXP(1.d0-(rmx/dist)**B) & + (nm2m*dist*percentCoriolis & *corio/2.d0)**2.d0) & - nm2m*dist*percentCoriolis*corio/2.d0 ! calculate the wind speed (m/s) at Rmax, using ! equation that includes full coriolis speedAtRmax = SQRT( & (Vmax*kt2ms)**2.d0 * & * EXP(0.d0) & + (nm2m*dist*percentCoriolis & *corio/2.d0)**2.d0) & - nm2m*dist*percentCoriolis*corio/2.d0 ! calculate a factor to place the velocity profile so that ! it hits vmax vmaxFactor = Vmax*kt2ms / speedAtRmax ! jgf20111007: Calculate NWS8-like translation speed TransSpdX = (abs(speed/speedAtRmax))*uTrans*kt2ms TransSpdY = (abs(speed/speedAtRmax))*vTrans*kt2ms speed = speed * vmaxFactor ! now reduce the wind speed to the surface speed = speed * windReduction u = -speed * COS(deg2rad*angle) v = speed * SIN(deg2rad*angle) ! ! Alter wind direction by adding a frictional inflow angle CALL rotate(u,v, fang(dist,rmx), cLat, uf,vf) u = uf v = vf ! ! jgf20111007: Add in the translation velocity u = u + TransSpdX v = v + TransSpdY ! ! convert from 1 minute averaged winds to 10 minute averaged ! winds for use in ADCIRC u = u * one2ten v = v * one2ten ! Compute surface pressure from asymmetric hurricane vortex p = mb2pa * (Pc + (Pn - Pc) * EXP(-(rmx/dist)**B)) ! cut off the vortex field after 401nm ! ! TODO: 401nm should be replaced with something less ! arbitrary ... and find a better way to blend this !if ( dist.gt.401.d0 ) then ! u = 0.0d0 ! v = 0.0d0 ! p = mb2pa * Pn !endif END SUBROUTINE uvp !================================================================= ! Calculate (u,v) wind components and surface pressure from an ! asymmetric hurricane wind model. ! ! On input: ! pinf hurricane ambient pressure ! p0 hurricane central pressure ! idist dist to hurricane center in nautical mile ! irmx RMW ! iangle Azimuth Angle ! iB Holland B parameter ! iVm vortex maximum velocity at upper boundary ! iPhi vortex correction factor ! uTrans x component of translational velocity (kts) ! vTrans y component of translational velocity (kts) ! ! On output: ! u x component of wind velocity at nodal point (m/s) ! v y component of wind velocity at nodal point (m/s) ! p Surface pressure at nodal point (Pa) ! ! Internal parameters: ! dampRadii How far out (# of Rmax radii) to begin damping ! out the translational velocity ! ! Note: ! Subroutine directly accesses global class instance variables !================================================================= SUBROUTINE uvpr(idist,iangle,irmx,irmx_true,iB,iVm,iPhi, & uTrans,vTrans,geof, u,v,p) REAL(sz), INTENT(IN) :: idist REAL(sz), INTENT(IN) :: iangle REAL(sz), INTENT(IN) :: irmx REAL(sz), INTENT(IN) :: irmx_true REAL(sz), INTENT(IN) :: iB REAL(sz), INTENT(IN) :: iVm REAL(sz), INTENT(IN) :: iPhi REAL(sz), INTENT(IN) :: uTrans REAL(sz), INTENT(IN) :: vTrans INTEGER , INTENT(IN) :: geof REAL(sz), INTENT(OUT) :: u REAL(sz), INTENT(OUT) :: v REAL(sz), INTENT(OUT) :: p REAL(sz) :: TransSpdX !NWS8-style translation speed REAL(sz) :: TransSpdY !NWS8-style translation speed REAL(sz) :: rmx REAL(sz) :: avgLat REAL(sz) :: speed REAL(sz) :: uf REAL(sz) :: vf REAL(sz) :: damp REAL(sz) :: percentCoriolis INTEGER :: i rmx=irmx B=iB Vmax=iVm Phi=iPhi !---------------------------------------- ! Handle special case at eye of hurricane !---------------------------------------- ! in eye velocity is zero not translational velocity !---------------------------------------- IF (idist < 1.d0) THEN u = 0.0d0 v = 0.0d0 p = Pc * mb2pa RETURN END IF !--------------------------------------------------- ! Compute (u,v) wind velocity components from the ! asymmetric hurricane vortex. ! ! Note: the vortex winds are valid at the top of the ! surface layer, so reduce the winds to the surface. ! Also convert the winds from max sustained 1-minute ! averages to 10-minute averages for the storm surge ! model. !--------------------------------------------------- percentCoriolis=1.0d0 ! Jie 2014.07 IF (geof == 1) THEN speed = SQRT( & ((Vmax*kt2ms)**2.d0+Vmax*kt2ms*rmx*nm2m & *percentCoriolis*corio) & * (rmx/idist)**B & * EXP(Phi*(1.d0-(rmx/idist)**B)) & + (nm2m*idist*percentCoriolis & *corio/2.d0)**2.d0) & - nm2m*idist*percentCoriolis*corio/2.d0 ELSE speed = SQRT( & (Vmax*kt2ms)**2.d0 * (rmx/idist)**B & * EXP(1.d0-(rmx/idist)**B) & + (nm2m*idist*percentCoriolis & *corio/2.d0)**2.d0) & - nm2m*idist*percentCoriolis*corio/2.d0 ENDIF ! jgf20111007: Calculate NWS8-like translation speed TransSpdX = (abs(speed/(Vmax*kt2ms)))*uTrans*kt2ms TransSpdY = (abs(speed/(Vmax*kt2ms)))*vTrans*kt2ms ! now reduce the wind speed to the surface speed = speed * windReduction u = -speed * COS(deg2rad*iangle) v = speed * SIN(deg2rad*iangle) ! ! Alter wind direction by adding a frictional inflow angle CALL rotate(u,v, fang(idist,irmx_true), cLat, uf,vf) u = uf v = vf ! ! jgf20111007: Add in the translation velocity u = u + TransSpdX v = v + TransSpdY ! ! convert from 1 minute averaged winds to 10 minute averaged ! winds for use in ADCIRC u = u * one2ten v = v * one2ten ! Compute surface pressure from asymmetric hurricane vortex IF (geof == 1) THEN p = mb2pa * (Pc + (Pn - Pc) * EXP(-Phi*(rmx/idist)**B)) ELSE p = mb2pa * (Pc + (Pn - Pc) * EXP(-(rmx/idist)**B)) ENDIF ! cut off the vortex field after 401nm ! ! TODO: 401nm should be replaced with something less ! arbitrary ... and find a better way to blend this !if ( dist.gt.401.d0 ) then ! u = 0.0d0 ! v = 0.0d0 ! p = mb2pa * Pn !endif END SUBROUTINE uvpr !================================================================= ! Spatial Interpolation function based on angle and r. ! ! On input: ! angle Azimuthal angle (degrees) ! r Distnace to storm Center (nm) ! ! On output: ! interpolated value for Rmax/Vmax/B ! ! INTEGER valid_Isot is used as a marker to indicate how many isotachs ! are available in a certain quadrant ! select case(valid_Isot) ! case(1): 1 situation ! case(2): 3 situations ! case(3): 4 situations ! case(4): 5 situations !================================================================= REAL(sz) FUNCTION spInterp(angle,dist,opt) REAL(sz), INTENT(IN) :: angle,dist INTEGER, INTENT(IN) :: opt REAL(sz), DIMENSION(nPoints,4) :: param REAL(sz) :: temp1, temp2 REAL(sz) :: delta_angle INTEGER :: iquad if (opt == 1) then param = Rmaxes4 else if (opt == 2) then param = Bs4 else if (opt == 3) then param = VmBL4 end if if (angle.le.45.0d0) then iquad = 5 delta_angle = 45.d0 + angle else if (angle.le.135.d0) then iquad = 2 delta_angle = angle - 45.d0 else if (angle.le.225.d0) then iquad = 3 delta_angle = angle - 135.d0 else if (angle.le.315.d0) then iquad = 4 delta_angle = angle - 225.d0 else if (angle.gt.315.d0) then iquad = 5 delta_angle = angle - 315.d0 endif ! nearest neighbor weighted interpolation if ( delta_angle.lt.1.d0 ) then spInterp=interpR(param, iquad, dist) else if (delta_angle.gt.89.d0) then spInterp=interpR(param, iquad+1, dist) else temp1= interpR(param, iquad, dist) temp2= interpR(param, iquad+1, dist) spInterp = ( temp1/delta_angle**2.d0 & + temp2/(90.0-delta_angle)**2.d0 ) & / ( 1.d0/delta_angle**2.d0 & + 1.d0/(90.d0-delta_angle)**2.d0 ) end if END FUNCTION spInterp REAL(sz) FUNCTION interpR(val,quad,r) REAL(sz), DIMENSION(nPoints,4),INTENT(IN) :: val INTEGER, INTENT(IN) :: quad REAL(sz), INTENT(IN) :: r REAL(sz) :: fac INTEGER :: total_Isot total_Isot=sum(QuadFlag4(quad,:)) select case(total_Isot) case(1) interpR = val(quad,maxloc(QuadFlag4(quad,:),1)) case(2) if (r > QuadIr4(quad,1)) then interpR = val(quad,1) else if (r > QuadIr4(quad,2)) then fac = (r-QuadIr4(quad,2))/ & (QuadIr4(quad,1)-QuadIr4(quad,2)) interpR = val(quad,1)* fac+ & val(quad,2)* (1-fac) else interpR = val(quad,2) end if case(3) if (r > QuadIr4(quad,1)) then interpR = val(quad,1) else if (r > QuadIr4(quad,2)) then fac = (r-QuadIr4(quad,2))/ & (QuadIr4(quad,1)-QuadIr4(quad,2)) interpR = val(quad,1)* fac+ & val(quad,2)* (1-fac) else if (r > QuadIr4(quad,3)) then fac = (r-QuadIr4(quad,3))/ & (QuadIr4(quad,2)-QuadIr4(quad,3)) interpR = val(quad,2)* fac+ & val(quad,3)* (1-fac) else interpR = val(quad,3) end if case(4) if (r > QuadIr4(quad,1)) then interpR = val(quad,1) else if (r > QuadIr4(quad,2)) then fac = (r-QuadIr4(quad,2))/ & (QuadIr4(quad,1)-QuadIr4(quad,2)) interpR = val(quad,1)* fac+ & val(quad,2)* (1-fac) else if (r > QuadIr4(quad,3)) then fac = (r-QuadIr4(quad,3))/ & (QuadIr4(quad,2)-QuadIr4(quad,3)) interpR = val(quad,2)* fac+ & val(quad,3)* (1-fac) else if (r > QuadIr4(quad,4)) then fac = (r-QuadIr4(quad,4))/ & (QuadIr4(quad,3)-QuadIr4(quad,4)) interpR = val(quad,3)* fac+ & val(quad,4)* (1-fac) else interpR = val(quad,4) end if case default ! For whatever reason if our algorithm fails, add the following ! line to avoid run-time errors interpR = val(quad,maxloc(QuadFlag4(quad,:),1)) write(*,*) "ERROR: interpR failed in nws20get." end select END function interpR !================================================================= ! calculate the radius of maximum winds. ! ! On input: ! angle Azimuthal angle (degrees) ! ! On output: ! Rmw Radius of maximum winds (meters) from curve fit ! I DO NOT BELIEVE IT IS IN METERS rjw !================================================================= REAL(sz) FUNCTION Rmw(angle) REAL(sz), INTENT(IN) :: angle INTEGER :: base_quadrant REAL(sz) :: delta_angle if (angle.le.45.0d0) then base_quadrant = 5 delta_angle = 45.d0 + angle else if (angle.le.135.d0) then base_quadrant = 2 delta_angle = angle - 45.d0 else if (angle.le.225.d0) then base_quadrant = 3 delta_angle = angle - 135.d0 else if (angle.le.315.d0) then base_quadrant = 4 delta_angle = angle - 225.d0 else if (angle.gt.315.d0) then base_quadrant = 5 delta_angle = angle - 315.d0 endif ! nearest neighbor weighted interpolation if ( delta_angle.lt.1.d0 ) then Rmw = Rmaxes(base_quadrant) ! avoid div by zero else if ( delta_angle.gt.89.d0 ) then Rmw = Rmaxes(base_quadrant+1) ! avoid div by zero else Rmw = ( Rmaxes(base_quadrant)/delta_angle**2.d0 & + Rmaxes(base_quadrant+1)/(90.0-delta_angle)**2.d0 ) & / ( 1.d0/delta_angle**2.d0 & + 1.d0/(90.d0-delta_angle)**2.d0 ) endif ! linearly interpolate C Rmw = (delta_angle/90.0d0) C & *(Rmaxes(base_quadrant+1)-Rmaxes(base_quadrant)) C & + Rmaxes(base_quadrant) END FUNCTION Rmw !======================================================== ! Rotate a 2D vector (x,y) by an angle. ! ! On input: ! x x component of vector ! y y component of vector ! angle angle to rotate vector (degrees) ! whichWay direction of rotation: ! - = clockwise, + = counter-clockwise ! ! On output: ! xr x component of rotated vector ! yr y component of rotated vector !======================================================== SUBROUTINE rotate(x,y, angle, whichWay, xr,yr) REAL(sz), INTENT(IN) :: x REAL(sz), INTENT(IN) :: y REAL(sz), INTENT(IN) :: angle REAL(sz), INTENT(IN) :: whichWay REAL(sz), INTENT(OUT) :: xr REAL(sz), INTENT(OUT) :: yr REAL(sz) :: A, cosA, sinA A = SIGN(1.d0, whichWay) * deg2rad * angle cosA = COS(A) sinA = SIN(A) xr = x * cosA - y * sinA yr = x * sinA + y * cosA END SUBROUTINE rotate !======================================================= ! Compute a wind angle to parameterize frictional inflow ! across isobars. ! ! On input: ! r distance from center of storm ! rmx radius of maximum winds ! ! On output: ! fang frictional inflow angle (degrees) !======================================================= REAL(sz) FUNCTION fang(r, rmx) REAL(sz), INTENT(IN) :: r REAL(sz), INTENT(IN) :: rmx IF (0.0d0 <= r .AND. r < rmx) THEN fang = 10.0d0 * r/rmx ELSE IF (rmx <= r .AND. r < 1.2d0*rmx) THEN fang = 10.0d0 + 75.0d0 * (r/rmx - 1.d0) ELSE IF (r >= 1.2d0*rmx) THEN fang = 25.0d0 ELSE fang = 0.0d0 END IF END FUNCTION fang !============================================================ ! Calculate the translational velocity of a moving hurricane. ! ! On input: ! latOld Previous latitude of center (degrees north) ! lonOld Previous longitude of center (degrees east ) ! latNew Current latitude of center (degrees north) ! lonNew Current longitude of center (degrees east ) ! tOld Previous time (seconds) ! tNew Current time (seconds) ! ! On output: ! uTrans x component of translational velocity (m/s) ! vTrans y component of translational velocity (m/s) !============================================================ SUBROUTINE uvtrans(latOld,lonOld, latNew,lonNew, tOld,tNew, & uTrans,vTrans) REAL(sz), INTENT(IN) :: latOld REAL(sz), INTENT(IN) :: lonOld REAL(sz), INTENT(IN) :: latNew REAL(sz), INTENT(IN) :: lonNew REAL(sz), INTENT(IN) :: tOld REAL(sz), INTENT(IN) :: tNew REAL(sz), INTENT(OUT) :: uTrans REAL(sz), INTENT(OUT) :: vTrans C REAL(sz) :: dx REAL(sz) :: dy REAL(sz) :: dt REAL(sz) :: avgLat avgLat = (latOld + latNew) / 2.d0 dx = sphericalDistance( ((lonNew-lonOld)*deg2rad), & 0.0d0, avgLat, avgLat) dy = sphericalDistance(0.0d0,((latNew-latOld)*deg2rad), & latOld,latNew) dt = tNew - tOld uTrans = sign(dx/dt,(lonNew-lonOld)) vTrans = sign(dy/dt,(latNew-latOld)) END SUBROUTINE uvtrans !================================================================= ! Transform (lat,lon) --> (x,y) coordinates. ! ! On input: ! lat Latitude (degrees north) ! lon Longitude (degrees east ) ! lat0 Latitude where projection is true (degrees north) ! lon0 Longitude where projection is true (degrees east ) ! ! On output: ! x x (meters) ! y y (meters) !================================================================= SUBROUTINE latlon2xy(lat,lon, lat0,lon0, x,y) REAL(sz), INTENT(IN) :: lat ,lon REAL(sz), INTENT(IN) :: lat0,lon0 REAL(sz), INTENT(OUT) :: x,y x = deg2rad * Rearth * (lon - lon0) * COS(deg2rad*lat0) y = deg2rad * Rearth * lat END SUBROUTINE latlon2xy !================================================================= ! Transform (x,y) --> (lat,lon) coordinates. ! ! On input: ! x x (meters) ! y y (meters) ! lat0 Latitude where projection is true (degrees north) ! lon0 Longitude where projection is true (degrees east ) ! ! On output: ! lat Latitude (degrees north) ! lon Longitude (degrees east ) !================================================================= SUBROUTINE xy2latlon(x,y, lat0,lon0, lat,lon) REAL(sz), INTENT(IN) :: x,y REAL(sz), INTENT(IN) :: lat0,lon0 REAL(sz), INTENT(OUT) :: lat ,lon lat = y / (deg2rad * Rearth) lon = lon0 + x / (deg2rad * Rearth * COS(deg2rad*lat0)) END SUBROUTINE xy2latlon !=============================================================== ! RJW 07 - 2009 ! Calculate the coefficients that fit the given ! radius of maximum winds for all storm quadrants. ! ! On input: ! Rmax in all 4 quadrants plus 2 extra values to tie down circular periodicity ! ! On output: ! Rmax radius of maximum winds (nm) in all quadrants, plus ! 2 extra values to tie down circular periodicity !=============================================================== SUBROUTINE fitRmaxes() REAL(sz) :: root ! Radius of maximum winds INTEGER :: n, iter,i ! ! Generate 2 additional (theta,Rmax) points for curve-fit Rmaxes(1) = Rmaxes(5) Rmaxes(6) = Rmaxes(2) END SUBROUTINE fitRmaxes SUBROUTINE fitRmaxes4() REAL(sz) :: root ! Radius of maximum winds INTEGER :: n, iter,i ! ! Generate 2 additional points for curve-fit QuadFlag4(1,1:4) = QuadFlag4(5,1:4) QuadFlag4(6,1:4) = QuadFlag4(2,1:4) QuadIr4(1,1:4) = QuadIr4(5,1:4) QuadIr4(6,1:4) = QuadIr4(2,1:4) Rmaxes4(1,1:4) = Rmaxes4(5,1:4) Rmaxes4(6,1:4) = Rmaxes4(2,1:4) Bs4(1,1:4) = Bs4(5,1:4) Bs4(6,1:4) = Bs4(2,1:4) Phis4(1,1:4) = Phis4(5,1:4) Phis4(6,1:4) = Phis4(2,1:4) VmBL4(1,1:4) = VmBL4(5,1:4) VmBL4(6,1:4) = VmBL4(2,1:4) END SUBROUTINE fitRmaxes4 SUBROUTINE setUseQuadrantVr(u) LOGICAL, INTENT(IN) :: u useQuadrantVr = u END SUBROUTINE setUseQuadrantVr SUBROUTINE setIsotachWindSpeeds(vrq) REAL(sz), DIMENSION(4), INTENT(IN) :: vrq VrQuadrant(:) = vrq(:) END SUBROUTINE setIsotachWindSpeeds SUBROUTINE setIsotachWindSpeed(sp) REAL(sz), INTENT(IN) :: sp Vr = sp END SUBROUTINE setIsotachWindSpeed SUBROUTINE setIsotachRadii(ir) REAL(sz), DIMENSION(4), INTENT(IN) :: ir radius(:) = ir(:) END SUBROUTINE setIsotachRadii REAL(sz) FUNCTION getShapeParameter() getShapeParameter = B END FUNCTION getShapeParameter SUBROUTINE setShapeParameter(param) REAL(sz) :: param B = param END SUBROUTINE setShapeParameter SUBROUTINE getRmaxes(rmaxw) REAL(sz), DIMENSION(4), INTENT(OUT) :: rmaxw INTEGER :: i do i=1,4 rmaxw(i) = Rmaxes(i+1) end do END SUBROUTINE getRmaxes SUBROUTINE setRmaxes(rmaxw) REAL(sz), DIMENSION(4), INTENT(IN) :: rmaxw INTEGER :: i do i=1,4 Rmaxes(i+1) = rmaxw(i) end do END SUBROUTINE setRmaxes REAL(sz) FUNCTION getLatestRmax() getLatestRmax = latestRmax END FUNCTION getLatestRmax REAL(sz) FUNCTION getLatestAngle() getLatestAngle = latestAngle END FUNCTION getLatestAngle LOGICAL FUNCTION getUseQuadrantVr() getUseQuadrantVr = useQuadrantVr END FUNCTION getUseQuadrantVr !Jie 2013.02 SUBROUTINE setUseVmaxesBL(u) LOGICAL, INTENT(IN) :: u useVmaxesBL = u END SUBROUTINE setUseVmaxesBL SUBROUTINE getVmaxesBL(vmaxw) REAL(sz), DIMENSION(4), INTENT(OUT) :: vmaxw INTEGER :: i do i=1,4 vmaxw(i) = VmBL(i+1) end do END SUBROUTINE getVmaxesBL SUBROUTINE setVmaxesBL(vmaxw) REAL(sz), DIMENSION(4), INTENT(IN) :: vmaxw INTEGER :: i do i=1,4 VmBL(i+1) = vmaxw(i) end do END SUBROUTINE setVmaxesBL FUNCTION getShapeParameters() REAL(sz),DIMENSION(4) :: getShapeParameters INTEGER :: i do i=1,4 getShapeParameters(i) = Bs(i+1) end do END FUNCTION getShapeParameters FUNCTION getPhiFactors() REAL(sz),DIMENSION(4) :: getPhiFactors INTEGER :: i do i=1,4 getPhiFactors(i) = Phis(i+1) end do END FUNCTION getPhiFactors END MODULE vortex