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 (ntrack.inp) 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 c Modified 02/02/2023 (KM) for NCO specifications 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 = 28) 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) then call ntrack_err_handling(1,6,istat) endif c call getSingleTAU ( comRcd, 0, tauData, istat ) if (istat .ne. 1) then call ntrack_err_handling(1,6,istat) endif 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) then call ntrack_err_handling(2,6,istat) endif c c t=0 hr call getSingleTAU (comRcd, 0, tauData, istat ) if (istat .ne. 1) then call ntrack_err_handling(2,6,istat) endif 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 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 enddo 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 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 enddo 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 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 ntrack_err_handling(ierrtype,luerr,ierr) c This routine handles error processing for tcliper. c Types of errors handled (ierrtype=): c -1: input file failed to open c -2: error reading input file c Input c ierrtype - error type to process c luerr - unit number to write to (assumes file is already open) c ierr - error code to report c IMPLICIT NONE c c ++ Passed variables integer, intent(in) :: ierrtype, luerr, ierr c if (ierrtype .eq. 1) then write(luerr,950) 950 format(/,' Error reading CARQ record from input file') stop elseif (ierrtype .eq. 2) then write(luerr,970) 970 format(/,' OFCL, OFCP, OFCI or JTWI not found in .com file. ', + ' Use track from track.dat file instead.') stop else write(luerr,999) 999 format(/,'Unrecognized error in ntrack') stop endif c return end subroutine ntrack_err_handling