program ntrack c c This program reads in the file ntrack.inp, and looks for the c official forecast track points. If the points are not at 12-h c intervals, the missing values are filled in using a quadratic c spline interpolation method. Then the forecast is checked to see c if extends to 120 h. If it does not, it is extrapolated up to c an additional 72 h. Finally, the date/time is checked to see c if the official forecast is older than the date and time c of the current forecast (as indicated on the CARQ lines of c the input file). If the ofcl forecast is 6 or 12 hours c old, it is shifted to the current forecast position and time. c c Modified Aug 2006 to fix east longitude problem c include 'dataformats.inc' c character *1 iewl character *2 basin,csnum character *4 tname character *8 strmid character*10 aymdh character*10 fninp,fndat c dimension olat(12),olon(12),otime(12) dimension itime(12) c dimension olat6(23),olon6(23),otime6(23) dimension ifofcl(11,3) c data otime /0.,12.,24.,36.,48.,60.,72.,84.,96.,108.,120.,132./ c data olat /12*0.0/ data olon /12*0.0/ data olat6 /23*0.0/ data olon6 /23*0.0/ c data fninp,fndat /'ntrack.inp','ntrack.dat'/ data luinp,ludat /21,31/ c type ( AID_DATA ) comRcd, tauData c c Set iclip=1 to extend forecast tracks up to 72 extra c hours using CLIPER. If iclip=0, then simple linear extrapolation c is used to extent the forecast from 114 to 120 hr when the c OFCI track is used. iclip = 1 c c Specify maximum and minimum acceptable longitude for including tracks rlongmx=179.0 rlongmn= 1.0 c do k=1,11 itime(k) = ifix(otime(k)) enddo c c call getarg ( 2, fninp ) ! use for HP c c Open the input file open(unit=luinp,file=fninp,status='old',form='formatted') c c **** Get OFCL forecast information **** tname = 'OFCI' call getARecord ( luinp, "OFCL", comRcd, istat ) if (istat .ne. 1) then rewind( luinp ) call getARecord ( luinp, "OFCP", comRcd, istat ) tname = 'OFPI' print *, "OFCP step... ", tname, istat endif if (istat .ne. 1) go to 901 c c t=0 hr call getSingleTAU (comRcd, 00, tauData, istat ) if (istat .ne. 1) go to 901 c aymdh = tauData%aRecord(1)%DTG olat(1) = tauData%aRecord(1)%lat olon(1) = tauData%aRecord(1)%lon iewl = tauData%aRecord(1)%EW if (iewl .eq. 'E') olon(1) = 180.0 + (180.0-olon(1)) c read(aymdh,100) iyro,imono,idayo,itimeo 100 format(i4,3(i2)) c c t=12-120 hrs do k=2,11 call getSingleTAU (comRcd, itime(k), tauData, istat ) if (istat .eq. 1) then olat(k) = tauData%aRecord(1)%lat olon(k) = tauData%aRecord(1)%lon iewl = tauData%aRecord(1)%EW if (iewl .eq. 'E') olon(k) = 180.0 + (180.0-olon(k)) else olat(k) = 0.0 olon(k) = 0.0 endif enddo c c **** Get CARQ information **** c Get storm position/date at forecast initialization time call getARecord ( luinp, "CARQ", comRcd, istat ) if (istat .eq. 0) go to 900 c call getSingleTAU ( comRcd, 0, tauData, istat ) if (istat .ne. 1) go to 900 c basin = tauData%aRecord(1)%basin snum = tauData%aRecord(1)%cyNum aymdh = tauData%aRecord(1)%DTG alat = tauData%aRecord(1)%lat alon = tauData%aRecord(1)%lon iewl = tauData%aRecord(1)%EW if (iewl .eq. 'E') alon = 180.0 + (180.0-alon) c read(aymdh,100) iyra,imona,idaya,itimea isnum = ifix(snum) write(csnum,101) isnum 101 format(i2.2) c c Fill in holes in the official track forecast c do k=1,12 c write(6,800) otime(k),olat(k),olon(k) 800 format(1x,f5.0,1x,f6.2,1x,f6.2) c enddo c call fill(olat,otime,12) call fill(olon,otime,12) c c Fill in holes in the official track forecast c do k=1,12 c write(6,800) otime(k),olat(k),olon(k) c enddo c if (iclip .eq. 1) then c Extend official track for up to 72 hours using CLIPER model call addclip(olat,olon,otime,basin,iyra,imona,idaya,itimea,12) else c c Extrapolate 120 hr forecast to 132 hr if possible if (olat(11) .gt. 0.0 .and. olat(10) .gt. 0.0 + .and. olat(12) .le. 0.0 ) then olat(12) = olat(11) + (olat(11)-olat(10)) olon(12) = olon(11) + (olon(11)-olon(10)) endif endif c c Check for longitude exceed maximum do k=1,12 if (olon(k) .gt. rlongmx) then olon(k) = 0.0 olat(k) = 0.0 endif enddo c c Check for longitude below mininum do k=1,12 if (olon(k) .le. rlongmn) then olon(k) = 0.0 olat(k) = 0.0 endif enddo c c do k=1,12 c write(6,800) otime(k),olat(k),olon(k) c enddo c c Calculate time difference between OFCL track c and forecast time call tdiff(iyra,imona,idaya,itimea, + iyro,imono,idayo,itimeo,idelt) c c write(6,*) 'idelt=',idelt c c Check to see if data can be interpolated if (idelt .gt. 12 .or. idelt .lt. 0) stop c if (idelt .eq. 0) then c OFCL tracks don't need to be time shifted, use them directly do 20 k=1,11 ifofcl(k,3) = 0 if (olat(k) .gt. 0.0) then ifofcl(k,1) = ifix(10.0*olat(k)) ifofcl(k,2) = ifix(10.0*olon(k)) else ifofcl(k,1) = 0 ifofcl(k,2) = 0 endif 20 continue elseif (idelt .eq. 12) then c Move track starting with the 12 h position do 30 k=1,11 ifofcl(k,3) = 0 if (olat(k) .gt. 0.0) then tlat = (alat-olat(2)) + olat(k+1) tlon = (alon-olon(2)) + olon(k+1) ifofcl(k,1) = ifix(10.0*tlat) ifofcl(k,2) = ifix(10.0*tlon) else ifofcl(k,1) = 0 ifofcl(k,2) = 0 endif 30 continue elseif (idelt .eq. 6) then c Move track starting with the 6 h position klast = 0 do 35 k=1,23 otime6(k) = 6.0*float(k-1) if (mod(k+2,2) .eq. 0) go to 35 kk = 1 + k/2 c if (olat(kk) .gt. 0.0) then olat6(k) = olat(kk) olon6(k) = olon(kk) klast = k else olat6(k) = 0.0 olon6(k) = 0.0 endif 35 continue c c do k=1,23 c write(6,800) otime6(k),olat6(k),olon6(k) c enddo c c Interpolate to get the 6 hour positions call fill(olat6,otime6,klast) call fill(olon6,otime6,klast) c c do k=1,klast c write(6,800) otime6(k),olat6(k),olon6(k) c enddo c do 40 k=1,11 ifofcl(k,3) = 0 kk = 2*k if (olat6(kk) .gt. 0.0) then tlat = (alat-olat6(2)) + olat6(kk) tlon = (alon-olon6(2)) + olon6(kk) ifofcl(k,1) = ifix(10.0*tlat) ifofcl(k,2) = ifix(10.0*tlon) else ifofcl(k,1) = 0 ifofcl(k,2) = 0 endif 40 continue endif c c Open the output file open(unit=ludat,file=fndat,form='formatted',status='replace') c strmid(1:2) = basin strmid(3:4) = csnum strmid(5:8) = aymdh(1:4) c call newWriteAidRcd(ludat,strmid,aymdh,tname,ifofcl) c c do k=1,12 c write(6,800) otime(k),olat(k),olon(k) c enddo c stop c 900 write(6,*) ' Error reading CARQ record from input file' stop c 901 write(6,*) ' Error reading OFCL record from input file' stop c end subroutine addclip(olat,olon,otime,basin,iyr,imon,iday,itime,nt) c This routine extends the track in olat,olon for up to 72 hours c using the CLIPER model c dimension olat(nt),olon(nt),otime(nt) c dimension fpclip(5,2),ftclip(6,2) c character *2 basin character *8 strmid c c Find the last time with a valid forecast klast = nt do k=1,nt if (olat(k) .gt. 0.0) klast = k enddo c c Check to see if CLIPER extrapolation is needed if (klast .ge. nt) return c c Check to make sure at least three forecasts are available if (klast .lt. 3) return c c Prepare input for call to CLIPER iymdh = itime + 100*iday + 10000*imon + 1000000*iyr c c default intensity/storm id vmax = 60.0 strmid = ' XXXXXX ' c c Calculate storm speeds alat = olat(klast) alon = olon(klast) ftim = otime(klast) c alatm1 = olat(klast-1) alonm1 = olon(klast-1) ftimm1 = otime(klast-1) c alatm2 = olat(klast-2) alonm2 = olon(klast-2) ftimm2 = otime(klast-2) c c Check for missing positions just to make sure if (alat .le. 0.0 .or. alatm1 .le. 0.0 .or. + alatm2 .le. 0.0 ) return c call distn(-alonm1,alatm1,-alon,alat,dx,dy,rad) u = dx/(ftim-ftimm1) v = dy/(ftim-ftimm1) call ctorh(u,v,spd,dir) c call distn(-alonm2,alatm2,-alonm1,alatm1,dx,dy,rad) u = dx/(ftim-ftimm1) v = dy/(ftim-ftimm1) call ctorh(u,v,spdm1,dirm1) c if (basin .eq. 'AL') then call atclip (iymdh,alat,alon,dir,spd,dirm1,spdm1, + vmax, strmid, fpclip) c c Linear interpolation to get the 60 hour CLIPER position do 10 j=1,2 do k=1,4 ftclip(k,j) = fpclip(k,j) enddo c ftclip(6,j) = fpclip(5,j) ftclip(5,j) = 0.5*(ftclip(6,j)+ftclip(4,j)) 10 continue c elseif (basin .eq. 'EP' .or. basin .eq. 'CP') then call epcl84 ( iymdh, alat, alon, alatm1, alonm1, + alatm2, alonm2, dir, spd, vmax, strmid, + ftclip) else c Unrecognized basin return endif c c Add CLIPER forecasts kk = 0 do k=klast+1,nt kk = kk+1 if (kk .le. 6) then olat(k) = ftclip(kk,1) olon(k) = ftclip(kk,2) endif enddo c return end subroutine fill(ft,t,nt) c This routine uses a quadratic spline to fill in missing values of ft. c If ft=0, then the value of f is assumed to be missing. c dimension ft(nt),t(nt) c parameter (mxt=100) dimension ftem(mxt),ttem(mxt) c if (nt .gt. mxt) then write(6,*) 'mxt too small in routine fill' stop endif c c Copy the data points with valid values to a sequential array kk = 0 do k=1,nt if (ft(k) .gt. 0.0) then kk = kk+1 klast = k c ftem(kk) = ft(k) ttem(kk) = t(k) endif enddo c c At least two valid times must be available for interpolation if (kk .lt. 2) return c c Initialize the interpolation call xint(ttem,ftem,kk,1,1,ttem(1),fi,ierr) c do k=1,klast call xint(ttem,ftem,kk,0,1,t(k),ft(k),ierr) enddo c return end subroutine tdiff(iy2,im2,id2,it2,iy1,im1,id1,it1,idelt) c This routine calculates the number of hours (delt) between c two date/times. c dimension nday(12) c data nday /0,31,59,90,120,151,181,212,243,273,304,334/ c c Calculate reference year iry = iy1-2 if (iy2 .lt. iry) iry=iy2-2 c c Calculate the number of hours from 00 Jan. 1 of the reference year ity1 = 0 do 10 i=iry,iy1-1 if (mod(i,4) .eq. 0) then ity1 = ity1 + 24*366 else ity1 = ity1 + 24*365 endif 10 continue c ity2 = 0 do 15 i=iry,iy2-1 if (mod(i,4) .eq. 0) then ity2 = ity2 + 24*366 else ity2 = ity2 + 24*365 endif 15 continue c ity1 = ity1 + 24*nday(im1) if ((mod(iy1,4) .eq. 0) .and. im1 .gt. 2) ity1=ity1+24 c ity2 = ity2 + 24*nday(im2) if ((mod(iy2,4) .eq. 0) .and. im2 .gt. 2) ity2=ity2+24 c ity1 = ity1 + 24*id1 + it1 ity2 = ity2 + 24*id2 + it2 c idelt = ity2 - ity1 c return end subroutine xint(x,f,n,iflag,lflag,xi,fi,ierr) c This routine applies a quadratic interpolation procedure c to f(x) between x(1) and x(n). f(x) is assumed to be c represented by quadratic polynomials between the points c x(i). The polynomials are chosen so that they equal f(i) c at the points x(i), the first derviatives on either c side of the interior x(i) match at x(i), and the second c derivative of the approximated function integrated c over the domain is minimized. c c This version is for interpolating longitude c c Input: x(1),x(2) ... x(n) The x values (must be sequential) c f(1),f(2) ... f(n) The function values c n The number of x,f pairs c iflag Flag for initialization c =1 for coefficient calculation c =0 to use previous coefficients c lflag Flag for linear interpolation c =0 to perform linear interpolation c =1 to perform quadratic interpolation c xi The x value at which to interpolate f c c Output: fi The interpolated function value c ierr Error flag c =0 Normal return c =1 Parameter nmax is too small c =2 The x values are not sequential c =3 Coefficient iteration did not c converge c =4 Mix-up finding coefficients c =5 if xi .gt. x(n) or .lt. x(1), c xi is set to nearest endpoint c before the interpolation c c Note: fi is set to -99.9 if c ierr=1,2,3 or 4 c parameter (nmax=100) c dimension x(n),f(n) c c Save variables dimension ax(nmax),bx(nmax),cx(nmax) c c Temporary local variables dimension df(nmax),dx(nmax),gm(nmax),ct(nmax) c common /xsave/ ax,bx,cx c c Specify unit number for debug write statements c and debug flag idbug = 0 lutest = 6 c c Initialize error flag ierr = 0 c c Specify minimum reduction in cost function for convergence thresh = 1.0e-10 c c Check to make sure nmax is large enough, and n is .gt. 1 if (n .gt. nmax .or. n .lt. 2) then ierr=1 fi = -99.9 return endif c if (iflag .eq. 1) then c Perform the initialization for later interpolation c c Check to make sure x is sequential do 10 i=1,n-1 if (x(i) .ge. x(i+1)) then ierr=2 fi = -99.9 return endif 10 continue c c Check for special case where n=2. Only linear interpolation c is possible. if (n .eq. 2) then cx(1) = 0.0 bx(1) = (f(2)-f(1))/(x(2)-x(1)) ax(1) = f(1) - bx(1)*x(1) go to 1500 endif c c Calculate x and f differences do 15 i=1,n-1 df(i) = f(i+1)-f(i) dx(i) = x(i+1)-x(i) 15 continue c c Calculate domain size d = x(n) - x(1) c c Check for linearity of input points eps = 1.0e-10 bb = (f(2)-f(1))/(x(2)-x(1)) aa = f(1) - bb*x(1) dev = 0.0 do 12 i=3,n dev = dev + abs(aa + bb*x(i) - f(i)) 12 continue c if (dev .lt. eps .or. lflag .eq. 0) then do 13 i=1,n-1 cx(i) = 0.0 13 continue go to 1000 endif c c Iterate to find the c-coefficients cx(1) = 0.0 nit = 100 slt = 0.01 cfsave = 1.0e+10 c do 20 k=1,nit c Calculate c values do 25 i=2,n-1 cx(i) = -cx(i-1)*dx(i-1)/dx(i) + -df(i-1)/(dx(i)*dx(i-1)) + +df(i )/(dx(i)*dx(i )) 25 continue c c Calculate current value of cost function cf0 = 0.0 do 26 i=1,n-1 cf0 = cf0 + cx(i)*cx(i)*dx(i) 26 continue cf0 = 0.5*cf0/d c if (idbug .ne. 0) then write(lutest,101) cf0 101 format(/,' cf0=',e13.6) endif c c Check for convergence rel = abs(cf0 - cfsave)/abs(cfsave) if (rel .lt. thresh) go to 1000 cfsave = cf0 c c Calculate values of Lagrange multipliers gm(n-1) = cx(n-1)*dx(n-1)/d c if (n .gt. 3) then do 30 i=n-2,2,-1 gm(i) = cx(i)*dx(i)/d - gm(i+1)*dx(i)/dx(i+1) 30 continue endif c c Calculate gradient of cost function with respect to c1 dsdc1 = dx(1)*(cx(1)/d - gm(2)/dx(2)) c c Adjust cx(1) using trial step ct(1) = cx(1) - slt*dsdc1 c c Calculate remaining c values at trial step do 33 i=2,n-1 ct(i) = -ct(i-1)*dx(i-1)/dx(i) + -df(i-1)/(dx(i)*dx(i-1)) + +df(i )/(dx(i)*dx(i )) 33 continue c c Calculate cost function at trial step cft = 0.0 do 31 i=1,n-1 cft = cft + ct(i)*ct(i)*dx(i) 31 continue cft = 0.5*cft/d c c write(6,*) 'dsdc1,cft,cf0',dsdc1,cft,cf0 c Calculate optimal step length and re-adjust cx(1) den = 2.0*((cft-cf0) + slt*dsdc1*dsdc1) if (den .ne. 0.0) then slo = dsdc1*dsdc1*slt*slt/den else slo =0.0 endif c c Adjust slo if desired slo = 1.0*slo c cx(1) = cx(1) - slo*dsdc1 c if (idbug .ne. 0) then write(lutest,100) k,cft,slt,slo 100 format(' Iteration=',i4,' cf1=',e11.4,' slt=',e11.4, + ' slo=',e11.4) c do 99 j=1,n-1 write(lutest,102) j,cx(j) 102 format(' i=',i2,' c=',f8.4) 99 continue endif c c Calculate trial step for next time step slt = 0.5*slo 20 continue c c Iteration did not converge ierr=3 fi=-99.9 return c c Iteration converged 1000 continue c if (idbug .ne. 0) then write(lutest,104) 104 format(/,' Iteration converged') endif c c Calculate b and a coefficients do 40 i=1,n-1 bx(i) = df(i)/dx(i) - cx(i)*(x(i+1) + x(i)) ax(i) = f(i) - bx(i)*x(i) - cx(i)*x(i)*x(i) 40 continue endif c 1500 continue c Interpolate the function c c Check for xi out of bounds if (xi .lt. x(1)) then xi = x(1) ierr = 5 endif c if (xi .gt. x(n)) then xi = x(n) ierr = 5 endif c c Find the interval for the interpolation ii = 1 do 50 i=2,n if (xi .le. x(i)) then ii = i-1 go to 2000 endif 50 continue c fi = -99.9 ierr=4 return c 2000 continue fi = ax(ii) + bx(ii)*xi + cx(ii)*xi*xi c return end subroutine ctorh(x,y,r,theta) c This routine converts from Cartesion coordinates c to radial coordinates, where theta is in c degrees measured clockwise from c the +y-axis (standard meteorological heading). c r = sqrt(x*x + y*y) c if (r .le. 0.0) then theta = 0.0 return endif c rtd = 57.296 theta = rtd*acos(x/r) if (y .lt. 0.0) theta = 360.0 - theta c c Convert theta to heading theta = 90.0 - theta if (theta .lt. 0.0) theta = theta + 360.0 c return end subroutine distn(rlon1,rlat1,rlon2,rlat2,dx,dy,rad) c This routine calculates the distance in nm (rad) between the c points (rlon1,rlat1) and (rlon2,rlat2) using an approximate c formula. The lon and lat are in deg E and N. The east and c north components of the distance (dx,dy) are also calculated. c dtn = 60.0 dtr = 0.0174533 c cfac = cos(0.5*dtr*(rlat1+rlat2)) c dx = dtn*(rlon2-rlon1)*cfac dy = dtn*(rlat2-rlat1) rad = sqrt(dx*dx + dy*dy) c return end subroutine atclip ( iYMDH, Y, X, DNOW, SNOW, & DM12, SM12, WKTS, STRMID, & FPCLIP ) C C** ATLANTIC CLIPER PROGRAM C C** REFER TO NOAA TECHNICAL MEMORANDUM NWS SR-62 (C.J. NEUMANN, C** JAN 1972) FOR A DESCRIPTION OF CLIPER SYSTEM. C dimension FPCLIP(5,2) DIMENSION DISP(10),FULL(10),GULF(10),CNS(14,5),CEW(8,5),ND(12), 1 P(14),XG(5),YG(5),XA(5),YA(5),FP(2,5) C CHARACTER*1 EW CHARACTER*2 FTIME(5) CHARACTER*8 STRMID C C** LIST CONSTANTS FOR MERIDIONAL DISPLACEMENTS. C DATA FTIME /'12','24','36','48','72'/ DATA CNS/7.60553,13.59909,-2.575127,-.0001868914,.00460007, 1.0022623,-.001491833,-.0002678624,-.00006994195,.00004407501, 2-.0002049626,.00007781249,.1430621,-.00008156078,30.30846, 322.91538,-2.484599,.004968631,.009297729,.02511378,-.007839641, 4-.005977559,-.0003504201,.0001557468,-.0009956548,.0004778924, 5.3879478,-.000671019,67.69324,31.94291,-3.697592,.009672852, 6.009539232,.06322068,-.01332191,-.01610885,-.0007292257, 7.0002341562,-.002814067,.001145255,.8940798,-.002181595,120.27143, 838.94701,-4.380877,.01323013,.02292699,.09532129,-.01664333, 9-.03200883,-.001216522,.000315372,-.005450743,.001867606,1.666658, 1-.004353084,263.15653,48.41731,-4.456658,.01126704,.04297187, 2.16962,-.01748416,-.06485724,-.002215515,.0003599106,-.012677, 3.003690242,4.121246,-.01102184/ C C** LIST CONSTANTS FOR ZONAL DISPLACEMENTS. C DATA CEW/-3.52591,13.69309,-2.637347,.8151257,.6867776, 1-.002168753,-.000596625,.1247267,-13.12388,23.30256,-3.215529, 23.584506,3.949364,-.007860247,-.006764091,.5135556,-28.48156, 332.37355,-5.342858,8.073875,9.321241,-.01318171,-.02040824, 41.044624,-44.13759,38.93667,-6.819777,14.10797,16.35476, 5-.01967039,-.03853289,1.698018,-60.23074,46.26022,-8.8089, 629.11625,32.91178,-.02181549,-.08553838,3.291178/ C DATA ND /0,31,59,90,120,151,181,212,243,273,304,334/ C C** WRITE THE INPUT PARAMETER TO THE PRINTER C c WRITE (*,'('' ATLANTIC CLIPER ''/)') c WRITE (*,'('' INPUT PARAMETERS FOR TROPICAL CYCLONE '',A, c & '' ON YMDH = '',I10.10,/)') STRMID,iYMDH c WRITE (*,'('' LATCUR = '',F4.1,''N LONCUR = '',F5.1, c & ''W DIRCUR = '',F4.0,'' DEG SPDCUR = '',F3.0, c & '' KT'',/,'' WNDCUR = '',F4.0,'' KT'',18X, c & ''DIRM12 = '',F4.0,'' DEG SPDM12 = '',F3.0, c & '' KT'',/)') Y, X, DNOW, SNOW, WKTS, DM12, SM12 C C** CHECK FOR STORM LOCATION AND DETERMINE WEIGHTING FACTORS C WA = (X - Y - 61.0)/5.0 WB = (Y - 13.0)/4.0 IF (WA.GT.1.0) WA = 1.0 IF (WA.LT.0.0) WA = 0.0 IF (WB.GT.1.0) WB = 1.0 IF (WB.LT.0.0) WB = 0.0 WGULF = WA*WB WFULL = 1.0 - WGULF C C** COMPUTE DAY NUMBER AND SUBTRACT OFF MEANS C IY = iYMDH/1000000 MO = iYMDH/10000 - IY*100 KDA = iYMDH/100 - MO*100 - IY*10000 C DAN = ND(MO) + KDA DANBR = ND(MO) + KDA - 248 IF (DANBR.LT.-95.0) DANBR = -95.0 IF (DANBR.GT. 82.0) DANBR = 82.0 C WMPH = (WKTS*1.15) - 71.0 ALAT = Y - 24.0 ALON = X - 68.0 C C** CONVERT PRESENT AND PAST MOTION TO U AND V COMPONENTS. C T = 0.0174533 UNOW = SIN(DNOW*T)*SNOW VNOW = COS(DNOW*T)*SNOW UM12 = SIN(DM12*T)*SM12 VM12 = COS(DM12*T)*SM12 C C** CALL THE GULF OF MEXICO CORRECTIONS C CALL CGULF (Y,X,UNOW,VNOW,UM12,VM12,DAN,WKTS,WGULF,GULF) C C** SET UP MERIDIONAL PREDICTORS AND COMPUTE MERIDIONAL DISPLACEMENTS. C P(1) = 1.0 P(2) = VNOW P(3) = VM12 P(4) = VNOW*VM12*VM12 P(5) = WMPH*VM12 P(6) = VNOW*WMPH P(7) = VNOW*VNOW*VM12 P(8) = ALAT*ALAT*VNOW P(9) = DANBR*DANBR*VM12 P(10) = VNOW*DANBR*DANBR P(11) = ALAT*ALAT*DANBR P(12) = WMPH*DANBR*VM12 P(13) = UNOW P(14) = DANBR*DANBR C DO 20 I = 1,5 ZZ = 0.0 C DO 10 J = 1,14 ZZ = ZZ + CNS(J,I)*P(J) 10 CONTINUE C FULL(2*I - 1) = ZZ 20 CONTINUE C C** SET UP ZONAL PREDICTORS AND COMPUTE ZONAL DISPLACEMENTS. C P(2) = UNOW P(3) = UM12 P(4) = ALAT P(5) = VNOW P(6) = VNOW*VNOW*UM12 P(7) = ALAT*VNOW*UM12 P(8) = ALON C DO 40 I = 1,5 ZZ = 0.0 C DO 30 J = 1,8 ZZ = ZZ + CEW(J,I)*P(J) 30 CONTINUE C C** CHANGE SIGN TO DESIGNATE WESTWARD MOTION AS POSITIVE MOTION. C FULL(2*I) = -ZZ 40 CONTINUE C DO 50 I = 1,10 DISP(I) = GULF(I)*WGULF + FULL(I)*WFULL 50 CONTINUE C DO 60 I = 1,5 YG(I) = GULF(2*I - 1)/60.0 + Y XG(I) = GULF(2*I)/COS((YG(I) + Y)*T/2.0)/60.0 + X IF (WGULF.EQ.0.0) YG(I) = -99.9 IF (WGULF.EQ.0.0) XG(I) = -99.9 YA(I) = FULL(2*I - 1)/60.0 + Y XA(I) = FULL(2*I)/COS((YA(I) + Y)*T/2.0)/60.0 + X c FPCLIP(i,1) = DISP(2*I - 1)/60.0 + Y FPCLIP(i,2) = DISP(2*I)/COS((FPCLIP(i,1) + Y)*T/2.0)/60.0 + X c 60 CONTINUE C return END C*************************************************************** SUBROUTINE CGULF (Y,X,U,V,U12,V12,DN,W,WGULF,GULF) C C** THIS SUBROUTINE WAS ADDED JUNE 1980 AND INCLUDES A CLIPER EQUATION C** SET DEVELOPED FOR THE WESTERN GULF OF MEXICO AND NORTHWEST C CARIBBEAN SEA. IT IS CALLED BY SUBROUTINE CLIPER. C DIMENSION GULF(10),X12(37),X24(37),X36(37),X48(37),X72(37), 1 P(37),Y12(37),Y24(37),Y36(37),Y48(37),Y72(37) C DATA Y12/ 0.0888784, 2.0189338, -3.0232588, 5.0112631, 1 -0.4663586, -2.2601204, 1.6083838, -0.0002354, 0.0025722, 2 0.0380934, -0.0006782, -0.0451460, 0.0196852, 0.0053249, 3 -0.0941109, 0.1065220, 0.0371968, -0.0074174, 0.0114476, 4 0.0397118, -0.1690010, -0.1136812, -0.0028529, 0.0703457, 5 -0.0153966, -0.0423261, 0.1190649, 0.0001786, 0.0059967, 6 -0.0224450, -0.0375346, 0.1719109, 0.2308200, -0.1857073, 7 -0.0995993, 0.0085389, 138.45727/ C DATA Y24/ -0.1899584, 8.6911949, -7.7038102, 17.8842992, 1 1.8958812, -26.6646610, 0.8565066, -0.0005621, -0.0008909, 2 0.0850158, 0.0034124, -0.1292738, 0.0426631, 0.0168735, 3 -0.1717200, 0.0842524, 0.0193852, -0.0207373, 0.1633104, 4 0.0636077, -0.5315965, -0.4433265, -0.0009782, 0.2091911, 5 0.1755244, 0.0419906, 0.3086234, -0.0357209, 0.0119171, 6 -0.1367056, -0.0408844, 0.5093530, 0.8909345, -0.5621746, 7 -0.3903772, 0.0554313, 408.02599/ C DATA Y36/ -1.2791723, 18.5437755, -21.4965205, 15.7440568, 1 13.7586086, -56.0136527, -1.0273326, -0.0010307, -0.0173412, 2 -0.0069898, 0.0192933, -0.1414892, 0.0882439, 0.0368659, 3 -0.2220821, 0.1678605, 0.1462510, -0.0429246, -0.0773177, 4 0.1059194, -1.3342283, -0.8653569, 0.0101134, 0.3178639, 5 0.4078366, -0.0984645, 0.7358741, 0.0350543, 0.0221573, 6 0.2828421, -0.1563257, 1.0009459, 1.6636230, -1.0310655, 7 -0.6516668, 0.1635486,1203.12003/ C DATA Y48/ -2.2503651, 29.2124117, -40.3642497, 17.0779631, 1 17.2671773,-104.0852329, 12.0955831, -0.0017919, -0.0235381, 2 -0.0449264, 0.0317965, -0.2155600, 0.1666921, 0.0817288, 3 0.2846074, -0.0195034, 0.2448353, -0.0474278, -0.1532889, 4 0.1664857, -1.9062156, -1.2207714, 0.0097954, -0.1064522, 5 1.0404667, -0.2380823, 0.9551006, -0.0139472, 0.0130756, 6 0.4827187, -0.3369744, 0.9455972, 2.3926721, -0.9962841, 7 -0.9012016, 0.2849380,2280.38586/ C DATA Y72/ -1.4432551, 69.0610138, -95.5025662, -27.2247875, 1 29.0697411,-193.8726068, 19.9328988, -0.0065750, -0.0171360, 2 -0.0367320, 0.0388394, -0.6790639, 0.4813834, 0.2157335, 3 0.8932160, 0.0837914, 0.4136060, -0.0491940, -0.3513719, 4 0.2171922, -2.6074776, -1.9936009, -0.0604859, -1.0243483, 5 2.4907242, -0.6472058, 1.1425050, -0.0513188, -0.0240031, 6 0.6859001, -0.3869234, 0.6265983, 3.8490617, -0.7979159, 7 -1.4598369, 0.3810998,4739.46017/ C DATA X12/ -0.6704965, 1.4162022, -4.6100138, 1.9871832, 1 17.5537836, -2.1569492, -5.3402337, 0.0012716, 0.0030864, 2 0.0278260, 0.0000683, -0.0343753, 0.0304263, 0.0099592, 3 0.0788021, -0.0810208, 0.1417678, 0.0012298, -0.0018823, 4 -0.0419311, 0.0747932, -0.0306750, -0.0050167, -0.0070392, 5 0.0516586, -0.1393882, 0.0512521, 0.0067553, -0.0037287, 6 -0.0017485, 0.0441335, -0.2081854, 0.1950879, 0.0235002, 7 -0.1346830, -0.0285769, 259.47152/ C DATA X24/ -3.2076529, 18.7032497, -51.1117663, 5.6067213, 1 41.0624529, -5.9989359, -10.5029546, 0.0062741, 0.0186705, 2 0.1565185, -0.0014309, -0.3151049, 0.3294738, 0.0497310, 3 0.4094322, -0.3283910, 0.2947786, 0.0020355, 0.0202198, 4 -0.1759863, 0.1977871, -0.0611235, -0.0238707, -0.0192411, 5 0.1737112, -0.2251423, 0.3112377, -0.0422888, -0.0167690, 6 0.0152779, 0.1156717, -0.8251600, 0.5527846, 0.0455937, 7 -0.3788840, -0.1737905,2398.43577/ C DATA X36/ -6.6781488, 45.4614176,-144.0815478, -5.8840358, 1 69.7781262, -2.2740308, -4.6481818, 0.0134484, 0.0379086, 2 0.3815790, -0.0039753, -0.7343225, 0.9044577, 0.1020808, 3 1.2687996, -0.5415504, 0.4282047, 0.0083583, 0.2549853, 4 -0.4604174, 0.0714267, -0.0236414, -0.0567410, -0.3394735, 5 0.3102145, -0.4093667, 0.5577365, -0.0197161, -0.0410260, 6 0.0190018, 0.1124278, -1.3487227, 0.5703272, -0.0312319, 7 -0.3382746, -0.4661488,6604.67254/ C DATA X48/ -11.2833695, 67.7041708,-213.4578700, 7.2101560, 1 103.3901233, -19.5321088, -14.6757722, 0.0229843, 0.0483794, 2 0.4193408, -0.0036961, -0.9820638, 1.3240829, 0.1520701, 3 2.3483309, -1.0601856, 0.4061403, 0.0143723, 0.9105408, 4 -0.9182075, 0.0862826, -0.1668057, -0.0868560, -0.8808472, 5 0.7276123, -0.3159147, 0.5817031, -0.0911473, -0.0663834, 6 -0.4230014, 0.3969854, -2.0183565, 0.7029859, 0.2263893, 7 -0.3186548, -0.7828428,9918.98936/ C DATA X72/ -20.0923662, 140.1359479,-364.9217391, 111.3688979, 1 174.4471782,-107.3250567, -43.3118878, 0.0429588, 0.0710944, 2 0.5375573, -0.0133502, -1.8051834, 2.2847824, 0.2595816, 3 4.8826019, -3.0362177, 0.2986331, 0.0258721, 2.1804674, 4 -1.8921960, 0.2580539, -0.5958524, -0.1663261, -2.3306786, 5 2.2924075, -0.7559446, 0.7681488, 0.0810918, -0.1449324, 6 -1.4833308, 1.2019937, -3.7551080, 1.0364360, 0.7572207, 7 -0.3288727, -1.4009940,16794.94253/ C C** INITIALIZE PREDICTAND ARRAY WITH ZEROS C DO 10 I = 1,10 GULF(I) = 0.0 10 CONTINUE C IF (WGULF.EQ.0.0) RETURN C C** COMPUTE PREDICTORS C P(1) = DN P(2) = Y P(3) = X P(4) = V P(5) = U P(6) = V12 P(7) = U12 C L = 7 DO 20 I = 1,7 K = I DO 20 J = 1,K L = L + 1 P(L) = P(I)*P(J) 20 CONTINUE C P(36) = W P(37) = 1.0 C C** COMPUTE DISPLACEMENTS C DO 30 I = 1,37 GULF(1) = GULF(1) + Y12(I)*P(I) GULF(2) = GULF(2) - X12(I)*P(I) GULF(3) = GULF(3) + Y24(I)*P(I) GULF(4) = GULF(4) - X24(I)*P(I) GULF(5) = GULF(5) + Y36(I)*P(I) GULF(6) = GULF(6) - X36(I)*P(I) GULF(7) = GULF(7) + Y48(I)*P(I) GULF(8) = GULF(8) - X48(I)*P(I) GULF(9) = GULF(9) + Y72(I)*P(I) GULF(10) = GULF(10) - X72(I)*P(I) 30 CONTINUE C RETURN END subroutine epcl84 ( ymdh, la0, lo0, lam12, lom12, lam24, lom24, & dir0, spd0, wnd0, STRMID, fpclip ) c C** Charlie's eastern Pacific CLIPER c c Minor modifications added by MD (4/01) for use with ntrack program. c COMMON /CL84/ CDI(12) COMMON /PP/ P(164) c dimension fpclip( 6, 2 ) DIMENSION PAV(8),INX(88),INDAY(13),cfo(2,6) c integer*4 ymdh real la0,lo0,lam12,lom12,lam24,lom24 DOUBLE PRECISION DISP(2,6),C(88),COF(2,164,6),AICP(2,7) C CHARACTER*8 STRMID character*256 cof_file_name character*256 coef_location CALL getenv( "SHIPS_COEF", coef_location ) c DATA PAV /16.8,120.3,230.6,3.4,-7.5,3.3,-7.8,61.5/ DATA INDAY /0,31,59,90,120,151,181,212,243,273,304,334,365/ c c** Write the input parameter to the printer c c write (*,'('' EASTERN PACIFIC CLIPER'',/)') c write (*,'('' INPUT PARAMETERS FOR TROPICAL CYCLONE '',A, c & '' ON YMDH = '',i10,/)') STRMID, ymdh c write (*,'('' LATCUR = '',f4.1,'' N LONCUR = '',f6.1, c 1 '' W'',/, c & '' DIRCUR = '',f4.0,'' DEG SPDCUR = '',f4.0,'' KT WNDCUR = '', c 2 f4.0,'' KT'',/, c & '' LATM12 = '',f4.1,'' N LONM12 = '',f6.1,'' W'',/, c 3 '' LATM24 = '',f4.1,'' N LONM24 = '',f6.1,'' W'',/)') c 4 la0,lo0,dir0,spd0,wnd0,lam12,lom12,lam24,lom24 c iyear = ymdh/1000000 imonth = ymdh/10000 - iyear*100 iday = ymdh/100 - imonth*100 - iyear*10000 ihour = ymdh - iday*100 - imonth*10000 - iyear*1000000 c DAYN = FLOAT(INDAY(imonth) + iday) + FLOAT(ihour)/24.0 IF (DAYN.LT.135.0) DAYN = 135.0 IF (DAYN.GT.349.0) DAYN = 349.0 DIR = dir0/57.29578 V0 = spd0*COS(DIR) U0 = spd0*SIN(DIR) c CALL STHGPR (lam12,lom12,360.0,1.0,0.0,0.0) CALL LL2XYH (la0,lo0,XF,YF) CALL LL2XYH (lam24,lom24,XP,YP) c V1 = (YF - YP)/24.0 U1 = (XF - XP)/24.0 c P(1) = la0 - PAV(1) P(2) = lo0 - PAV(2) P(3) = DAYN - PAV(3) P(4) = V0 - PAV(4) P(5) = U0 - PAV(5) P(6) = V1 - PAV(6) P(7) = U1 - PAV(7) P(8) = wnd0 - PAV(8) c CALL MULT1 CALL MULT2 CALL MULT3 CALL MULT4 c CALL STHGPR (la0,lo0,360.0,1.0,0.0,0.0) c DO 10 NDIR = 1,2 DO 10 J = 1,164 DO 10 NTIM = 1,6 COF(NDIR,J,NTIM) = 0.0 10 continue c c** Open regression coefficient file c ccc cof_file_name = 'epcliper.cof' cof_file_name = trim( coef_location )//'epcliper.cof' c OPEN ( 22, FILE=cof_file_name, STATUS='OLD', iostat=ios, err=1010) c NTIM = 1 20 READ (22,1,END=40,iostat=ios,err=1020) NDIR,NC, & (INX(I),C(I),I = 1,NC),AICP(NDIR,NTIM) 1 FORMAT (2I5,/,11(4(I3,E15.7),/),E15.7) c cc write (*,1) NDIR,NC,(INX(I),C(I),I = 1,NC),AICP(NDIR,NTIM) cc pause ' Hit to continue.' c DO 30 I = 1,NC J = INX(I) - 2 COF(NDIR,J,NTIM) = C(I) 30 continue c IF (NDIR.EQ.1) GO TO 20 NTIM = NTIM + 1 GO TO 20 c 40 close (22) c DO 50 NDIR = 1,2 DO 50 NTIM = 1,6 DISP(NDIR,NTIM) = AICP(NDIR,NTIM) 50 continue c DO 60 NDIR = 1,2 DO 60 NTIM = 1,6 DO 60 J = 1,164 IF (COF(NDIR,J,NTIM).EQ.0.) GO TO 60 DISP(NDIR,NTIM) = DISP(NDIR,NTIM) + COF(NDIR,J,NTIM)*P(J) NF = 2*(NTIM - 1) + NDIR CDI(NF) = SNGL(DISP(NDIR,NTIM)) 60 CONTINUE c DO 70 NTIM = 1,6 NF = 2*NTIM AX = CDI(NF) AY = CDI(NF - 1) c CALL XY2LLH (AX,AY,CFO(1,Ntim),CFO(2,Ntim)) c 70 continue c do i = 1,6 fpclip(i,1) = cfo(1,i) fpclip(i,2) = cfo(2,i) enddo c c write (*,'('' FORECAST POSITIONS'',/)') c WRITE ( 6, 3 ) la0,lo0,(fpclip(i,1),fpclip(i,2), i = 1, 6) c 3 FORMAT (' 00 HR = ',F5.1,' N ',F6.1,' W', c 1 /,' 12 HR = ',F5.1,' N ',F6.1,' W', c 2 /,' 24 HR = ',F5.1,' N ',F6.1,' W', c 3 /,' 36 HR = ',F5.1,' N ',F6.1,' W', c 4 /,' 48 HR = ',F5.1,' N ',F6.1,' W', c 5 /,' 60 HR = ',F5.1,' N ',F6.1,' W', c 6 /,' 72 HR = ',F5.1,' N ',F6.1,' W',/) c RETURN c c** Error messages c 1010 print *,' Error openning coeficient data file = ',ios stop c 1020 print *,' Error reading coefficient data file = ',ios stop c END c*********************************************************************** SUBROUTINE MULT1 c COMMON /PP/ P(164) c P( 9)=P(1)*P(1) P(10)=P(2)*P(2) P(11)=P(3)*P(3) P(12)=P(4)*P(4) P(13)=P(5)*P(5) P(14)=P(6)*P(6) P(15)=P(7)*P(7) P(16)=P(8)*P(8) P(17)=P(1)*P(2) P(18)=P(1)*P(3) P(19)=P(1)*P(4) P(20)=P(1)*P(5) P(21)=P(1)*P(6) P(22)=P(1)*P(7) P(23)=P(1)*P(8) P(24)=P(2)*P(3) P(25)=P(2)*P(4) P(26)=P(2)*P(5) P(27)=P(2)*P(6) P(28)=P(2)*P(7) P(29)=P(2)*P(8) P(30)=P(3)*P(4) P(31)=P(3)*P(5) P(32)=P(3)*P(6) P(33)=P(3)*P(7) P(34)=P(3)*P(8) P(35)=P(4)*P(5) P(36)=P(4)*P(6) P(37)=P(4)*P(7) P(38)=P(4)*P(8) P(39)=P(5)*P(6) P(40)=P(5)*P(7) P(41)=P(5)*P(8) P(42)=P(6)*P(7) P(43)=P(6)*P(8) P(44)=P(7)*P(8) P(45)=P(1)*P(1)*P(1) P(46)=P(2)*P(2)*P(2) P(47)=P(3)*P(3)*P(3) c RETURN END c***************************************************************** SUBROUTINE MULT2 c COMMON /PP/ P(164) c P(48)=P(4)*P(4)*P(4) P(49)=P(5)*P(5)*P(5) P(50)=P(6)*P(6)*P(6) P(51)=P(7)*P(7)*P(7) P(52)=P(8)*P(8)*P(8) P(53)=P(1)*P(1)*P(2) P(54)=P(1)*P(1)*P(3) P(55)=P(1)*P(1)*P(4) P(56)=P(1)*P(1)*P(5) P(57)=P(1)*P(1)*P(6) P(58)=P(1)*P(1)*P(7) P(59)=P(1)*P(1)*P(8) P(60)=P(2)*P(2)*P(1) P(61)=P(2)*P(2)*P(3) P(62)=P(2)*P(2)*P(4) P(63)=P(2)*P(2)*P(5) P(64)=P(2)*P(2)*P(6) P(65)=P(2)*P(2)*P(7) P(66)=P(2)*P(2)*P(8) P(67)=P(3)*P(3)*P(1) P(68)=P(3)*P(3)*P(2) P(69)=P(3)*P(3)*P(4) P(70)=P(3)*P(3)*P(5) P(71)=P(3)*P(3)*P(6) P(72)=P(3)*P(3)*P(7) P(73)=P(3)*P(3)*P(8) P(74)=P(4)*P(4)*P(1) P(75)=P(4)*P(4)*P(2) P(76)=P(4)*P(4)*P(3) P(77)=P(4)*P(4)*P(5) P(78)=P(4)*P(4)*P(6) P(79)=P(4)*P(4)*P(7) P(80)=P(4)*P(4)*P(8) P(81)=P(5)*P(5)*P(1) P(82)=P(5)*P(5)*P(2) P(83)=P(5)*P(5)*P(3) P(84)=P(5)*P(5)*P(4) P(85)=P(5)*P(5)*P(6) P(86)=P(5)*P(5)*P(7) c RETURN END c************************************************************ SUBROUTINE MULT3 c COMMON /PP/ P(164) c P(87)=P(5)*P(5)*P(8) P(88)=P(6)*P(6)*P(1) P(89)=P(6)*P(6)*P(2) P(90)=P(6)*P(6)*P(3) P(91)=P(6)*P(6)*P(4) P(92)=P(6)*P(6)*P(5) P(93)=P(6)*P(6)*P(7) P(94)=P(6)*P(6)*P(8) P(95)=P(7)*P(7)*P(1) P(96)=P(7)*P(7)*P(2) P(97)=P(7)*P(7)*P(3) P(98)=P(7)*P(7)*P(4) P(99)=P(7)*P(7)*P(5) P(100)=P(7)*P(7)*P(6) P(101)=P(7)*P(7)*P(8) P(102)=P(8)*P(8)*P(1) P(103)=P(8)*P(8)*P(2) P(104)=P(8)*P(8)*P(3) P(105)=P(8)*P(8)*P(4) P(106)=P(8)*P(8)*P(5) P(107)=P(8)*P(8)*P(6) P(108)=P(8)*P(8)*P(7) P(109)=P(1)*P(2)*P(3) P(110)=P(1)*P(2)*P(4) P(111)=P(1)*P(2)*P(5) P(112)=P(1)*P(2)*P(6) P(113)=P(1)*P(2)*P(7) P(114)=P(1)*P(2)*P(8) P(115)=P(1)*P(3)*P(4) P(116)=P(1)*P(3)*P(5) P(117)=P(1)*P(3)*P(6) P(118)=P(1)*P(3)*P(7) P(119)=P(1)*P(3)*P(8) P(120)=P(1)*P(4)*P(5) P(121)=P(1)*P(4)*P(6) P(122)=P(1)*P(4)*P(7) P(123)=P(1)*P(4)*P(8) P(124)=P(1)*P(5)*P(6) P(125)=P(1)*P(5)*P(7) c RETURN END c************************************************************** SUBROUTINE MULT4 c COMMON /PP/ P(164) c P(126)=P(1)*P(5)*P(8) P(127)=P(1)*P(6)*P(7) P(128)=P(1)*P(6)*P(8) P(129)=P(1)*P(7)*P(8) P(130)=P(2)*P(3)*P(4) P(131)=P(2)*P(3)*P(5) P(132)=P(2)*P(3)*P(6) P(133)=P(2)*P(3)*P(7) P(134)=P(2)*P(3)*P(8) P(135)=P(2)*P(4)*P(5) P(136)=P(2)*P(4)*P(6) P(137)=P(2)*P(4)*P(7) P(138)=P(2)*P(4)*P(8) P(139)=P(2)*P(5)*P(6) P(140)=P(2)*P(5)*P(7) P(141)=P(2)*P(5)*P(8) P(142)=P(2)*P(6)*P(7) P(143)=P(2)*P(6)*P(8) P(144)=P(2)*P(7)*P(8) P(145)=P(3)*P(4)*P(5) P(146)=P(3)*P(4)*P(6) P(147)=P(3)*P(4)*P(7) P(148)=P(3)*P(4)*P(8) P(149)=P(3)*P(5)*P(6) P(150)=P(3)*P(5)*P(7) P(151)=P(3)*P(5)*P(8) P(152)=P(3)*P(6)*P(7) P(153)=P(3)*P(6)*P(8) P(154)=P(3)*P(7)*P(8) P(155)=P(4)*P(5)*P(6) P(156)=P(4)*P(5)*P(7) P(157)=P(4)*P(5)*P(8) P(158)=P(4)*P(6)*P(7) P(159)=P(4)*P(6)*P(8) P(160)=P(4)*P(7)*P(8) P(161)=P(5)*P(6)*P(7) P(162)=P(5)*P(6)*P(8) P(163)=P(5)*P(7)*P(8) P(164)=P(6)*P(7)*P(8) c RETURN END C*********************************************************************** BLOCK DATA C** THIS DATASET IS NOW NAMED' NWS.WD80.CJN.SOURCE1(PGM04) C** ALBION D. TAYLOR, MARCH 19, 1982 C** THE HURRICANE GRID IS BASED ON AN OBLIQUE EQUIDISTANT CYLINDRICAL C** MAP PROJECTION ORIENTED ALONG THE TRACK OF THE HURRICANE. C C** THE X (OR I) COORDINATE XI OF A POINT REPRESENTS THE DISTANCE C** FROM THAT POINT TO THE GREAT CIRCLE THROUGH THE HURRICANE, IN C** THE DIRECTION OF MOTION OF THE HURRICANE MOTION. POSITIVE VALUES C** REPRESENT DISTANCES TO THE RIGHT OF THE HURRICANE MOTION, NEGATIVE C** VALUES REPRESENT DISTANCES TO THE LEFT. C** THE Y (OR J) COORDINATE OF THE POINT REPRESENTS THE DISTANCE C** ALONG THE GREAT CIRCLE THROUGH THE HURRICANE TO THE PROJECTION C** OF THE POINT ONTO THAT CIRCLE. POSITIVE VALUES REPRESENT C** DISTANCE IN THE DIRECTION OF HURRICANE MOTION, NEGATIVE VALUES C** REPRESENT DISTANCE IN THE OPPOSITE DIRECTION. C C** SCALE DISTANCES ARE STRICTLY UNIFORM IN THE I-DIRECTION ALWAYS. C** THE SAME SCALE HOLDS IN THE J-DIRECTION ONLY ALONG THE HURRICANE TRAC C** ELSEWHERE, DISTANCES IN THE J-DIRECTION ARE EXAGERATED BY A FACTOR C** INVERSELY PROPORTIONAL TO THE COSINE OF THE ANGULAR DISTANCE FROM C** THE TRACK. THE SCALE IS CORRECT TO 1 PERCENT WITHIN A DISTANCE OF C** 480 NM OF THE STORM TRACK, 5 PERCENT WITHIN 1090 NM, AND C** 10 PERCENT WITHIN 1550 NM. C C** BIAS VALUES ARE ADDED TO THE XI AND YJ COORDINATES FOR CONVENIENCE C** IN INDEXING. C C** A PARTICULAR GRID IS SPECIFIED BY THE USER BY MEANS OF A CALL C** TO SUBROUTINE STHGPR (SET HURRICANE GRID PARAMETERS) C** WITH ARGUMENTS (XLATH,XLONH,BEAR,GRIDSZ,XIO,YJO) C** WHERE c C** XLATH,XLONH = LATITUDE, LONGITUDE OF THE HURRICANE C** BEAR = BEARING OF THE HURRICANE MOTION C** GRIDSZ = SIZE OF GRID ELEMENTS IN NAUTICAL MILES C** XIO, YJO = OFFSETS IN I AND J COORDINATES (OR I AND J C COORDINATES OF HURRICANE) C** AND WHERE c C** LATITUDES, LONGITUDES AND BEARINGS ARE GIVEN IN DEGREES, C** POSITIVE VALUES ARE NORTH AND WEST, NEGATIVE SOUTH AND EAST, C** BEARINGS ARE GIVEN CLOCKWISE FROM NORTH. C C** THE CALL TO STHGPR SHOULD BE MADE ONCE ONLY, AND BEFORE REFERENCE C** TO ANY CALL TO LL2XYH OR XY2LLH. IN DEFAULT, THE SYSTEM C** WILL ASSUME A STORM AT LAT,LONG=0.,0., BEARING DUE NORTH, C** WITH A GRIDSIZE OF 120 NAUTICAL MILES AND OFFSETS OF 0.,0. . C C** TO CONVERT FROM GRID COORDINATES XI AND YJ, USE A CALL TO C** CALL XY2LLH(XI,YJ,XLAT,XLONG) C** THE SUBROUTINE WILL RETURN THE LATITUDE AND LONGITUDE CORRESPONDING C** TO THE GIVEN VALUES OF XI AND YJ. C C** TO CONVERT FROM LATITUDE AND LONGITUDE TO GRID COORDINATES, USE C** CALL LL2XYH(XLAT,XLONG,XI,YJ) C** THE SUBROUTINE WILL RETURN THE I-COORDINATE XI AND Y-COORDINATE C** YJ CORRESPONDING TO THE GIVEN VALUES OF LATITUDE XLAT AND C** LONGITUDE XLONG. c COMMON /HGRPRM/ A(3,3),RADPDG,RRTHNM,DGRIDH,HGRIDX,HGRIDY c DATA A /0.,-1.,0., 1.,0.,0., 0.,0.,1./ DATA RADPDG/1.745 3293 E-2/,RRTHNM /3 440.17/ DATA DGRIDH/120./ DATA HGRIDX,HGRIDY/0.,0./ c END c********************************************************** SUBROUTINE STHGPR (XLATH,XLONH,BEAR,GRIDSZ,XI0,YJ0) c C** ALBION D. TAYLOR, MARCH 19, 1982 c COMMON /HGRPRM/ A(3,3),RADPDG,RRTHNM,DGRIDH,HGRIDX,HGRIDY c CLAT=COS(RADPDG*XLATH) SLAT=SIN(RADPDG*XLATH) SLON=SIN(RADPDG*XLONH) CLON=COS(RADPDG*XLONH) SBEAR=SIN(RADPDG*BEAR) CBEAR=COS(RADPDG*BEAR) c A(1,1)= CLAT*SLON A(1,2)= CLAT*CLON A(1,3)= SLAT A(2,1)= - CLON*CBEAR + SLAT*SLON*SBEAR A(2,2)= SLON*CBEAR + SLAT*CLON*SBEAR A(2,3)= - CLAT* SBEAR A(3,1)= - CLON*SBEAR - SLAT*SLON*CBEAR A(3,2)= SLON*SBEAR - SLAT*CLON*CBEAR A(3,3)= CLAT* CBEAR c DGRIDH=GRIDSZ HGRIDX=XI0 HGRIDY=YJ0 c RETURN END c**************************************************************** SUBROUTINE LL2XYH (XLAT,XLONG,XI,YJ) c C** ALBION D. TAYLOR, MARCH 19, 1982 c COMMON /HGRPRM/ A(3,3),RADPDG,RRTHNM,DGRIDH,HGRIDX,HGRIDY c DIMENSION ZETA(3),ETA(3) c CLAT=COS(RADPDG*XLAT) SLAT=SIN(RADPDG*XLAT) SLON=SIN(RADPDG*XLONG) CLON=COS(RADPDG*XLONG) ZETA(1)=CLAT*SLON ZETA(2)=CLAT*CLON ZETA(3)=SLAT DO 20 I=1,3 ETA(I)=0. DO 20 J=1,3 ETA(I)=ETA(I) + A(I,J)*ZETA(J) 20 CONTINUE R=SQRT(ETA(1)*ETA(1) + ETA(3)*ETA(3)) XI=HGRIDX+RRTHNM*ATAN2(ETA(2),R)/DGRIDH IF(R.LE.0.) GO TO 40 YJ=HGRIDY+RRTHNM*ATAN2(ETA(3),ETA(1))/DGRIDH RETURN c 40 YJ=0. RETURN c END c************************************************************* SUBROUTINE XY2LLH (XI,YJ,XLAT,XLONG) c C** ALBION D. TAYLOR, MARCH 19, 1982 c COMMON /HGRPRM/ A(3,3),RADPDG,RRTHNM,DGRIDH,HGRIDX,HGRIDY c DIMENSION ZETA(3),ETA(3) c CXI=COS(DGRIDH*(XI-HGRIDX)/RRTHNM) SXI=SIN(DGRIDH*(XI-HGRIDX)/RRTHNM) SYJ=SIN(DGRIDH*(YJ-HGRIDY)/RRTHNM) CYJ=COS(DGRIDH*(YJ-HGRIDY)/RRTHNM) ETA(1)=CXI*CYJ ETA(2)=SXI ETA(3)=CXI*SYJ DO 20 I=1,3 ZETA(I)=0. DO 20 J=1,3 ZETA(I)=ZETA(I) + A(J,I)*ETA(J) 20 CONTINUE R=SQRT(ZETA(1)*ZETA(1) + ZETA(2)*ZETA(2)) XLAT=ATAN2(ZETA(3),R)/RADPDG IF(R.LE.0.) GO TO 40 XLONG=ATAN2(ZETA(1),ZETA(2))/RADPDG RETURN c 40 XLONG=0. RETURN c END