program ntrack c c This program reads in the file ntrack.com, and looks for the c official NHC or JTWC 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. c c 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 To handle global storms, the longitude convention (east or west) is determined c by the initial stor position. c c For NHC/CPHC AOR (western hemisphere) the longitude convention is deg W negative. c For JTWC AOR (eastern hemisphere) the longitude convection is deg E postive. c c ivmax set to default value for all cases where valid track point exists. c This is needed by routine writeaidlocal2. c c Modified 9/15/2013 for global tropical cyclones. This includes changing array c indices so that t=0 is element 0 rather than 1. c Also, the old CLIPER routines were removed. c c Modified 10/31/2016 (MD) for WCOSS Cray c Input file name changed to ntrack.com for consistency with other c NHC guidance suite programs. c Subroutines tdiff, xint, ctorh, distn removed. The are now c linked via libraries c include 'dataformats.inc' c character *1 latNS,lonEW character *2 basin,csnum character *4 agenid character *8 strmid character*10 aymdh,aymdh00 character*10 fninp,fndat c parameter (mft = 20) parameter (mft2 = mft*2) dimension olat(0:mft),olon(0:mft),otime(0:mft) dimension itime(0:mft) dimension ilat(0:mft),ilon(0:mft),ivmax(0:mft) character *4 sclass(0:mft) c dimension olat6(0:mft2),olon6(0:mft2),otime6(0:mft2) c data fninp,fndat /'ntrack.com','ntrack.dat'/ data luinp,ludat /21,31/ c type ( AID_DATA ) comRcd, tauData c c Set default value for ivmax when track is valid ivmaxd = 25 c c Calculate times and initialize lat/lon variables to zero delt06 = 6.0 delt12 = 12.0 c do k=0,mft otime(k) = delt12*float(k) itime(k) = nint(otime(k)) c olat(k) = 0.0 olon(k) = 0.0 ivmax(k) = 0 sclass(k) = ' ' enddo c do k=0,mft2 otime6(k) = delt06*float(k) c olat6(k) = 0.0 olon6(k) = 0.0 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 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 latNS = tauData%aRecord(1)%NS lonEW = tauData%aRecord(1)%EW c c Save CARQ date/time group aymdh00 = aymdh c c Check for SH if (latNS .eq. 'S') alat = -1.0*abs(alat) c c Calculate initial longitude in 0-360 deg convention alon360 = alon if (lonEW .eq. 'W') alon360 = 360.0-alon c c Determine if this is an NHC/CPHC or JTWC case if (alat .lt. 0.0) then iloncon = 2 else if (alon360 .ge. 180.0 .or. alon360 .lt. 15.0) then iloncon = 1 else iloncon = 2 endif endif c c Set the analysis lat/lon by longitude convention alata = alat if (iloncon .eq. 1) then if (alon360 .gt. 15.0) then alona = alon360 - 360.0 else alona = alon360 endif else alona = alon360 endif c read(aymdh,100) iyra,imona,idaya,itimea isnum = nint(snum) write(csnum,101) isnum 101 format(i2.2) c c **** Get OFCL or JTWI forecast information if it is available **** if (iloncon .eq. 2) then rewind(luinp) call getARecord ( luinp, "JTWI", comRcd, istat ) else iprov = 0 rewind( luinp ) call getARecord ( luinp, "OFCL", comRcd, istat ) if (istat .ne. 1) then rewind( luinp ) call getArecord( luinp, "OFCP", comRcd, istat ) iprov = 1 endif if (istat .ne. 1) then rewind( luinp ) call getArecord( luinp, "OFCI", comRcd, istat ) iprov = -1 endif endif if (istat .ne. 1) go to 901 c c t=0 hr call getSingleTAU (comRcd, 0, tauData, istat ) if (istat .ne. 1) go to 901 c aymdh = tauData%aRecord(1)%DTG olat(0) = tauData%aRecord(1)%lat olon(0) = tauData%aRecord(1)%lon latNS = tauData%aRecord(1)%NS lonEW = tauData%aRecord(1)%EW if (latNS .eq. 'S') olat(0) = -1.0*abs(olat(0)) if (lonEW .eq. 'W') olon(0) = 360.0-olon(0) c read(aymdh,100) iyro,imono,idayo,itimeo 100 format(i4,3(i2)) c c t=12,24, ... hr do k=1,mft call getSingleTAU (comRcd, itime(k), tauData, istat ) if (istat .eq. 1) then olat(k) = tauData%aRecord(1)%lat olon(k) = tauData%aRecord(1)%lon latNS = tauData%aRecord(1)%NS lonEW = tauData%aRecord(1)%EW if (latNS .eq. 'S') olat(k) = -1.0*abs(olat(k)) if (lonEW .eq. 'W') olon(k) = 360.0-olon(k) else olat(k) = 0.0 olon(k) = 0.0 endif enddo c if (iloncon .eq. 1) then c Use deg W negative convention except for a small region east of the c Greenwich Meridian do k=0,mft iolatt = nint(olat(k)) if (iolatt .ne. 0) then if (olon(k) .gt. 15.0) then olon(k) = olon(k) - 360.0 endif endif enddo endif 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 800 format(1x,f5.0,1x,f6.2,1x,f6.2) c enddo c nextrap = 0 call fill(olat,otime,mft,nextrap) call fill(olon,otime,mft,nextrap) 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 do not need to be time shifted, use them directly nshift = 0 do 20 k=0,mft iolatt = nint(olat(k)) if (iolatt .ne. 0) then ilat(k) = nint(10.0*olat(k)) ilon(k) = nint(10.0*olon(k)) else ilat(k) = 0 ilon(k) = 0 endif 20 continue else c Shift track 6 or 12 hr for OFCI forecast c c Make sure the time difference is either 6 or 12 hr if (idelt .ne. 6 .and. idelt .ne. 12) stop c nshift = idelt/6 c c Copy the 12 h interval points to the 6 hour arrays do k=0,mft kk = 2*k olat6(kk) = olat(k) olon6(kk) = olon(k) enddo c c Fill in missing 6 hr values and extrapolate is desired nextrap = 2 call fill(olat6,otime6,mft2,nextrap) call fill(olon6,otime6,mft2,nextrap) c c do m=0,mft2 c write(6,*) otime6(m),olat6(m),olon6(m) c enddo c c Calculate the lat and lon offset values iolatt = nint(olat(nshift)) if (iolatt .eq. 0) stop c oslat = alata - olat6(nshift) oslon = alona - olon6(nshift) c c write(6,*) 'oslat,oslon: ',oslat,oslon c Extract 12 h interval postions from 6 hr positions and shift accordingly do 30 k=0,mft kk = 2*k + nshift if (kk .lt. mft2) then iolatt = olat6(kk) if (iolatt .ne. 0) then ilat(k) = nint( 10.0*(olat6(kk) + oslat) ) ilon(k) = nint( 10.0*(olon6(kk) + oslon) ) else ilat(k) = 0 ilon(k) = 0 endif endif 30 continue endif c c Set ivmax to default value for good track points do k=0,mft if (ilat(k) .ne. 0 .and. ilon(k) .ne. 0) ivmax(k) = ivmaxd enddo 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 if (iloncon .eq. 2) then agenid = 'JTWI' else if (nshift .eq. 0) then if (iprov .eq. 0) then agenid = 'OFCL' elseif (iprov .eq. -1) then agenid = 'OFCI' else agenid = 'OFCP' endif else if (iprov .eq. 0) then agenid = 'OFCI' else agenid = 'OFPI' endif endif endif c minc = 1 call writeaidlocal2(ludat,strmid,aymdh00,agenid, + ilat,ilon,ivmax,sclass,itime,mft,minc) c stop c 900 write(6,*) ' Error reading CARQ record from input file' stop c 901 write(6,*) ' OFCL, OFCP, OFCI or JTWI not found in .com file. ', + ' Use track from track.dat file instead.' stop c end subroutine fill(ft,t,nt,nextrap) 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. Set nextrap c to the number of time periods to extrapolate the data if it is missing c dimension ft(0:nt),t(0: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=0,nt iftt = nint(ft(k)) if (iftt .ne. 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 if (nextrap .gt. 0) then c Linearly extrapolate nextrap time intervals ftlastm0 = ft(klast) ftlastm1 = ft(klast-1) tlastm0 = t(klast) tlastm1 = t(klast-1) c if (tlastm0 .eq. tlastm1) return dfdt = (ftlastm0-ftlastm1)/(tlastm0-tlastm1) c klaste = klast + nextrap do k=klast+1,klaste if (k .le. nt) then ft(k) = ftlastm0 + dfdt*(t(k)-tlastm0) endif enddo endif 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 c 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