subroutine vradp(preds,slat,slon,spd,head,vmaxop,ias, + pvmx,ppmin,rm,x,pr34,pr50,pr64) c*********************************************************** c This subroutine was modified April 2005, by CIRA c with the new predictors and coefficients, c developed using global data from 1999-2004, c for estimating the MSW (pvmx), MSLP (ppmin), and wind radii. c Changes are denoted by *** prior to and after changes. c*********************************************************** c This routine estimates the storm maximum wind c and the radii of 34,50 and 64 kt winds from AMSU c analysis parameters. c c Input: c preds - Array of length 18 containing predictors from the c AMSU analysis (see below) c slat - Storm latitude (deg N) c slon - Storm longitude (deg W positive) c spd - storm translational speed c head - storm heading (deg clockwise from north) c vmaxop - Operational estimate of max winds (kt) c ias - Set ias=1 to skip asymmetric wind radii estimation c Set ias=0 to include asymmetric winds c c Output c pvmx - predicted maximum surface winds (kt) c ppmin- predicted minimum sea-level pressure (hPa) c pr34 - radius (nm) of 34 kt winds c pr50 - radius (nm) of 50 kt winds c pr64 - radius (nm) of 64 kt winds c rm - radius of max winds (nm) from vortex model fit c x - Decay exponent from vortex model fit c c Notes: c The 18 input AMSU predictors are as follows: c 1. analyzed pressure (hPa) at r,z=0 c 2. r=600 to r=0 pressure drop (hPa) at z=0 km c 3. r=600 to r=0 pressure drop (hPa) at z=3 km c 4. r=0 Max temperature anomaly (C) c 5. height (km) of r=0 max temp. anomaly c 6. swath spacing (km) c 7. max wind (kt) at z=0 c 8. radius (km) of z=0 max wind c 9. max wind (kt) at z=3 km c 10. radius (km) of z=3 max wind c 11. 0-250 km avg. wind (kt) at z=0 km c 12. 0-250 km avg. wind (kt) at z=3 km c 13. 0-250 km avg. wind (kt) at z=5 km c 14. 250-500 km avg. wind (kt) at z=0 km c 15. 250-500 km avg. wind (kt) at z=3 km c 16. 250-500 km avg. wind (kt) at z=5 km c 17. r=0 to r=100 km avg. CLW c 18. Percent CLW r=0 to r=300 exceeding 0.5 c c vmaxop and lat are included as possible predictors c for wind radii. These are predictors 19. and 20. c c (vmaxop should not be used to predict pvmx or ppmin) c*** c Four additional predictors are added to the pool. c Predictor 21 = tmax^2 (i.e., predictor 4 from above squared). c Predictor 22 = tmax*clwave (i.e., predictors 4*17). c Predictor 23 = clwave^2 (i.e., predictor 17 squared). c Predictor 24 = the pressure at r=600km (i.e., predictors 1+2 c from above), which is not AMSU-derived, rather it is acquired c from the NCEP GFS model and is used to derive pmin and dp0. c*** c c pr34,50,64 are 1-D arrays of dimension 6, where c pr(1) - NE radius c pr(2) - SE radius c pr(3) - SW radius c pr(4) - SW radius c pr(5) - azimuthal mean radius c pr(6) - azimuthal mean radius from vortex fit c dimension pr34(6),pr50(6),pr64(6), dummy(4) c parameter (nco=18) c*** nct is modified to accomodate the 4 additional predictors, 21-24 c parameter (nct=nco+2) parameter (nct=nco+6) c*** c dimension preds(nco),predt(nct) dimension cvmx(0:nct),cpmn(0:nct) dimension cr34(0:nct),cr50(0:nct),cr64(0:nct) c c*** Modified predictor layout to include predictors 21-24 c (i.e., tmax^2, tmax*clwave, clwave^2, and p600, respectively) c Coefficients also were modified to be based on 1999-2004 data. c 0 1 2 3 4 5 c 6 7 8 9 10 c 11 12 13 14 15 c 16 17 18 19 20 c 21 22 23 24 c c Coefficients based upon 1999-2004 global data data cvmx /13.12557, 0.00000,-4.33459, 6.48789, 6.28701, 0.00000, + 0.13380, 0.00000, 0.00000, 0.49635,-0.02713, + 0.00000, 0.00000, 1.72608, 1.85672,-2.48450, + 0.00000,19.84888,-0.26614, 0.00000, 0.00000, + -0.51428, 0.00000, 0.00000, 0.00000/ c data cpmn /15.17819, 0.00000,-0.04260, 0.06316, 0.07395, 0.00000, + 0.00153, 0.00000, 0.00000, 0.00537,-0.00029, + 0.00000, 0.00000, 0.01681, 0.01753,-0.02156, + 0.00000, 0.20820,-0.00209, 0.00000, 0.00000, + -0.00655, 0.00000, 0.00000, -0.01145/ c data cr34 /1734.869,-1.67423, 3.46238, 0.00000,10.41016, 0.00000, + 0.38938,-0.67877, 0.09078, 0.00000, 0.00000, + 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, + 0.00000,-39.6419, 0.30049, 0.00000, 0.00000, + 0.00000,-8.81319,30.62560, 0.00000/ c data cr50 /-5.12785, 0.00000, 2.56301, 0.00000, 6.85088, 0.00000, + 0.00000, 0.89570, 0.06395, 0.00000, 0.00000, + -5.54525,13.59224,-13.0449, 0.00000, 0.00000, + 0.00000,-18.6103, 0.27442, 0.49101, 0.00000, + 0.00000, 0.00000, 0.00000, 0.00000/ c data cr64 /-18.4571, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, + 0.37589, 2.10560, 0.00000,-2.52543, 0.00000, + -4.11899, 6.62236, 1.75455, 0.00000, 0.00000, + -1.94906, 0.00000, 0.15221, 0.00000, 0.00000, + 0.45719, 0.00000, 0.00000, 0.00000/ c*** c Coefficients based upon 1999-2001 data c data cvmx /14.036, 0.000,-2.171, 0.000, 2.917, 0.000, c + 0.190, 0.000, 0.000, 0.000,-0.028, c + 0.000, 0.000, 4.250, 0.581, 0.000, c + 0.000,20.050,-0.206, 0.000, 0.000/ c c data cpmn /450.28, 0.564, 1.808, 0.000,-2.737, 0.000, c + -0.145, 0.000, 0.000, 0.000, 0.020, c + 0.000, 0.000,-2.588,-0.371, 0.000, c + 0.000,-12.77, 0.103, 0.000, 0.000/ c c data cr34 /2406.1,-2.409, 4.482,-10.39, 0.000, 0.000, c + 0.000, 0.000, 0.000, 0.000, 0.000, c + 0.000, 0.000, 0.000, 0.000, 0.000, c + 2.493, 0.000, 0.590, 1.676, 0.891/ c c data cr50 /-63.59, 0.000, 0.000, 0.000,22.037, 0.000, c + 0.000, 0.000, 0.000, 0.000, 0.122, c + 0.000, 8.043,-13.62, 0.000, 0.000, c + 0.000, 0.000, 0.000, 0.000, 0.990/ c c data cr64 /-48.57, 0.000, 0.000, 0.000, 6.286, 0.000, c + 0.427, 0.000, 0.000, 0.000, 0.000, c + 0.000, 0.000, 0.000, 0.000, 0.000, c + 0.000, 0.000, 0.270, 0.000, 0.217/ c c Coefficients based upon 1999-2000 data c data cvmx /26.781, 0.000, 0.000,-2.315, 0.000, 0.000, c + 0.000, 0.000, 0.000, 0.000,-0.038, c + 0.000,-5.158,10.261, 1.169, 0.000, c + -1.477,11.319, 0.000, 0.000, 0.000/ c c data cpmn /1005.5, 0.000, 0.000, 0.000, 0.000, 0.000, c + 0.000, 0.000, 0.000, 0.000, 0.024, c + 0.000, 3.453,-6.190,-0.909, 0.000, c + 1.236,-8.436, 0.000, 0.298, 0.000/ c c data cr34 /-85.57, 0.000, 0.000,-0.658, 0.000, 0.000, c + 0.000, 0.000, 0.000, 0.000, 0.000, c + 0.000, 1.597, 0.000, 0.000, 0.000, c + 0.000,32.184, 0.000, 3.616, 0.873/ c c data cr50 /-72.25, 0.000, 0.000, 3.269, 0.000, 0.000, c + 0.000, 0.000, 0.000, 0.000, 0.000, c + 0.000,-0.890, 0.000, 0.000, 0.000, c + 0.000,19.278, 0.000, 1.609, 0.842/ c c data cr64 /-22.97, 0.000, 0.000, 4.087, 0.000, 0.000, c + 0.000, 0.000, 0.000, 0.000, 0.000, c + 0.000,-1.647, 0.000, 0.000, 0.000, c + 0.000, 8.009, 0.000, 0.231, 0.418/ c c Initialize wind radii to zero do 10 i=1,6 pr34(i) = 0.0 pr50(i) = 0.0 pr64(i) = 0.0 10 continue c c Copy preds array and append it with latitude and vmaxop do 12 i=1,nco predt(i) = preds(i) 12 continue predt(nco+1) = abs(slat) predt(nco+2) = vmaxop c c*** Create predictors 21 (tmax^2), 22 (tmax*clwave), c 23 (clwave^2), and 24 (p600=pmin+dp0) predt(nco+3) = preds(4) * preds(4) predt(nco+4) = preds(4) * preds(17) predt(nco+5) = preds(17) * preds(17) predt(nco+6) = preds(1) + preds(2) c*** c c Predict max wind, minimum pressure, and azimuthally averaged c 34, 50 and 64 kt wind radii pvmx = cvmx(0) ppmin = cpmn(0) rb34 = cr34(0) rb50 = cr50(0) rb64 = cr64(0) c do 15 k=1,nct pvmx = pvmx + predt(k)*cvmx(k) ppmin = ppmin + predt(k)*cpmn(k) rb34 = rb34 + predt(k)*cr34(k) rb50 = rb50 + predt(k)*cr50(k) rb64 = rb64 + predt(k)*cr64(k) 15 continue c c*** Modified prediction of MSLP (ppmin) is of ln(1050-MSLP), c so need to calculate MSLP (ppmin) from that ppmin = 1050. - exp(ppmin) c*** c if (rb34 .lt. 0.0) rb34=0.0 if (rb50 .lt. 0.0) rb50=0.0 if (rb64 .lt. 0.0) rb64=0.0 c pr34(5) = rb34 pr50(5) = rb50 pr64(5) = rb64 c if (ias .eq. 1) return c c Calcuate asymmetric wind radii call wrasym(rb34,rb50,rb64,vmaxop,spd,head,rm,x,pr34,pr50,pr64, . slat) c c Switch assymetry around for southern hemisphere c J. Knaff - January 10, 2003 c if (slat.lt.0.0) then c do i=1,4 c dummy(i)=pr34(i) c enddo c pr34(1)=dummy(4) c pr34(4)=dummy(1) c pr34(2)=dummy(3) c pr34(3)=dummy(2) c do i=1,4 c dummy(i)=pr50(i) c enddo c pr50(1)=dummy(4) c pr50(4)=dummy(1) c pr50(2)=dummy(3) c pr50(3)=dummy(2) c do i=1,4 c dummy(i)=pr64(i) c enddo c pr64(1)=dummy(4) c pr64(4)=dummy(1) c pr64(2)=dummy(3) c pr64(3)=dummy(2) c endif return end subroutine wrasym(rb34,rb50,rb64,vmx,spd,head,rm,x,pr34,pr50,pr64, . slat) c c The routine calculates the wind radii in 4 quadrants relative c to the storm center (NE,SE,SW,NW) given the azimuthally c average wind radii (rb34,rb50,rb64) the max winds (vmx), c and the speed and heading of the storm motion (spd,head). c c An idealized Rankine vortex with a wave number one asymmetry c is fitted to the mean wind radii to give the wind radii in c each qaudrant. The parameters of the Rankine vortex c (maximum wind radius rm and decay exponent x) are also returned. c c Input: rb34 - azimuthally averaged 34 kt wind radius c rb50 - azimuthally averaged 50 kt wind radius c rb64 - azimuthally averaged 64 kt wind radius c vmx - maximum wind (kt) c spd - storm speed of motion (kt) c head - storm heading (deg clockwise from N) c slat - storm latitude (degrees) c c Output: pr34(4) - array with 34 kt wind radii (nm) NE,SE,SW,NW of center c pr50(4) - array with 50 kt wind radii (nm) NE,SE,SW,NW of center c pr64(4) - array with 64 kt wind radii (nm) NE,SE,SW,NW of center c rm - radius of max winds (nm) from vortex model fit c x - Decay exponent from vortex model fit c dimension pr34(6),pr50(6),pr64(6) c c Internal work array dimension cf(0:125,0:125) c c Set weights for climatology penalty terms for x, rm al1 = 0.1 al2 = 0.1 c c Specify angle for adjusting max winds relative to the direction c 90 deg to the right of the direction of motion theta0 = 0.0 c c Specify weight for adjusting the asymmetry factor rasf = 0.6 c c Specify weight for case when wind threshold is too close c to vmx, or set wttc=-1.0 to calculate wttc based upon azimuthal c distance covered by each wind radii. wttc = -1.0 c c Specify search interval for x,rm x0 = 0.01 dx = 0.01 nx = 125 c rm0 = 5.0 dr = 1.0 nr = 125 c c Initialize output variables to zero do 10 k=1,4 pr34(k) = 0.0 pr50(k) = 0.0 pr64(k) = 0.0 10 continue c x = 0.0 rm = 0.0 c c Calculate asymmetry factor from spd if (spd .le. 0.0) then a = 0.0 else a = rasf*1.5*(spd**0.63) endif c vmxa = vmx - a vmx2a = vmx - 2.0*a c c Find azimuth covered by each wind radius if (a .le. 0.0) then ac34 = 360.0 ac50 = 360.0 ac64 = 360.0 else pi = 3.14159 rtd = 180.0/pi c if (vmx2a .ge. 34.0) then ac34 = 360.0 else if (vmx .le. 34.0) then ac34 = 0.0 else ac34 = 2.0*rtd*acos(1.0 - (vmx-34.0)/a) endif endif c if (vmx2a .ge. 50.0) then ac50 = 360.0 else if (vmx .le. 50.0) then ac50 = 0.0 else ac50 = 2.0*rtd*acos(1.0 - (vmx-50.0)/a) endif endif c if (vmx2a .ge. 64.0) then ac64 = 360.0 else if (vmx .le. 64.0) then ac64 = 0.0 else ac64 = 2.0*rtd*acos(1.0 - (vmx-64.0)/a) endif endif endif c c Set wttc variables aclow = 180.0 wtval = 0.1 c if (wttc .lt. 0.0) then if (ac34 .lt. aclow) then wttc34 = wtval elseif (ac34 .ge. 360.0) then wttc34 = 1.0 else wttc34 = wtval + (1.0-wtval)*(ac34-aclow)/(360.0-aclow) endif c if (ac50 .lt. aclow) then wttc50 = wtval elseif (ac50 .ge. 360.0) then wttc50 = 1.0 else wttc50 = wtval + (1.0-wtval)*(ac50-aclow)/(360.0-aclow) endif c if (ac64 .lt. aclow) then wttc64 = wtval elseif (ac64 .ge. 360.0) then wttc64 = 1.0 else wttc64 = wtval + (1.0-wtval)*(ac64-aclow)/(360.0-aclow) endif else wttc34 = wttc wttc50 = wttc wttc64 = wttc endif c c Check maximum wind and set values accordingly. if (vmx .lt. 34.0) then return elseif (vmx .ge. 34.0 .and. vmx .lt. 50.0) then w34 = wttc34 w50 = 0.0 w64 = 0.0 elseif (vmx .gt. 50.0 .and. vmx .lt. 64.0) then w34 = wttc34 w50 = wttc50 w64 = 0.0 else w34 = wttc34 w50 = wttc50 w64 = wttc64 endif c c Calculate climatological rm,x and their standard deviations c from empirical formulas rmc = 54.0 - .27*vmx if (rmc .lt. 18.0) rmc = 18.0 c srm = 33.0 - .21*vmx if (srm .lt. 6.0) srm = 6.0 c xc = .42 + .0025*vmx sx = .10 c c Specify wind radii standard deviations (indep. of vmax) s34 = 43.0 s50 = 30.0 s64 = 22.0 c c write(6,810) rmc,srm,xc,sx,ac34,ac50,ac64,w34,w50,w64,a c 810 format(' rm mean,std: ',f5.1,1x,f5.1,/, c + ' x mean,std: ',f5.3,1x,f5.3,/, c + ' ac34,50,64: ',f5.1,1x,f5.1,1x,f5.1,/, c + ' wt34,50,64: ',f5.3,1x,f5.3,1x,f5.3,/, c + ' asym factor: ',f5.1) c c Start search loop for x,rm do 20 i=0,nx do 20 j=0,nr rmt = rm0 + float(j)*dr xt = x0 + float(i)*dx c c Calculate mean radii for current values of x,rm call rbar(vmx,34.0,a,rmt,xt,tb34) call rbar(vmx,50.0,a,rmt,xt,tb50) call rbar(vmx,64.0,a,rmt,xt,tb64) c c Calculate cost function cf(i,j) = w34*((tb34-rb34)/s34)**2 + + w50*((tb50-rb50)/s50)**2 + + w64*((tb64-rb64)/s64)**2 + + al1*((xt -xc )/sx )**2 + + al2*((rmt -rmc )/srm)**2 20 continue c iprt = 0 if (iprt .eq. 1) then c Print cost function write(6,300) 300 format(/,' COST FUNCTION') c do 30 j=nr,0,-1 write(6,310) j,(cf(i,j),i=0,nx) 310 format(1x,i2,1x,11(f4.1,1x)) 30 continue write(6,320) (i,i=0,nx) 320 format(4x,11(1x,i2,2x)) endif c c Find cost function minimum cmin = 1.0e+10 do 40 j=0,nr do 40 i=0,nx if (cf(i,j) .lt. cmin) then imin = i jmin = j cmin = cf(i,j) endif 40 continue c x = x0 + dx*float(imin) rm= rm0+ dr*float(jmin) c c Calculate best fit mean radii call rbar(vmx,34.0,a,rm,x,fb34) call rbar(vmx,50.0,a,rm,x,fb50) call rbar(vmx,64.0,a,rm,x,fb64) c c write(6,200) fb34,fb50,fb64 c 200 format(' Fit 34,50,64 kt wind radii: ',3(f5.0,1x)) c c Put fit to mean radii in element 6 of pr arrays pr34(6) =fb34 pr50(6) =fb50 pr64(6) =fb64 c c Calculate wind radii in each quadrant pi = 3.14159 dtr= pi/180.0 xi = 1.0/x c hemfac=1.0 if (slat.lt.0.0)then ! a knaff change for SH. hemfac=-1.0 endif do 50 k=1,4 q = 45.0 + 90.0*(float(k-1)) theta = dtr*(head + hemfac*90.0 - q -theta0) c pr34(k) = rm*( (vmx-a)/(34.0-a*cos(theta)) )**xi if (pr34(k) .lt. rm) pr34(k) = 0.0 c pr50(k) = rm*( (vmx-a)/(50.0-a*cos(theta)) )**xi if (pr50(k) .lt. rm) pr50(k) = 0.0 c pr64(k) = rm*( (vmx-a)/(64.0-a*cos(theta)) )**xi if (pr64(k) .lt. rm) pr64(k) = 0.0 50 continue c return end subroutine rbar(vm,v,a,rm,x,rb) c This routine calculates the azimuthally averaged wind radii c (rb) for wind speed v, max wind vm, asymmetry factor a, radius c of max wind rm, and Rankine vortex factor x. c c Check for illegal values of x,rm if (x .le. 0.0 .or. rm .le. 0.0) then write(6,100) 100 format(/,' Illegal values of x or rm input to routine rbar') stop endif c c Check for wind threshold greater than max wind if (v .gt. vm) then rb = 0.0 return endif c pi = 3.14159 dtr= pi/180.0 xi = 1.0/x c nt = 72 dt = 360.0/float(nt) c rb = 0.0 c*** Set up new counter to only average azimuths with wind radii > 0 ncount = 0 c*** do 10 i=1,nt theta = dtr*dt*float(i) fac = (vm-a)/(v-a*cos(theta)) if (fac .lt. 1.0) then fac = 0.0 else fac = fac**xi c*** Increment new counter for azimuths with wind radii > 0 ncount = ncount + 1 c*** endif c rb = rb + fac 10 continue c c*** Check that ncount isn't 0 and modify calculation of rb c so that it's only averaged over azimuths with wind radii > 0 if (ncount .le. 0) then rb = rm*(1.1) else rb = rm*rb/float(ncount) endif c rb = rm*rb/float(nt) c*** c return end