program ocd5
!     This program makes an ocd5 forecast by making a decay-SHIFOR forecast
!     with the CLIPER track. The SHIFOR and CLIPER forecasts are read in from
!     input files. 
!
!     Input files: ocd5.com - contains CARQ lines for the forecast
!                  ocd5.inp - File containing the CLIPER and SHIFOR forecasts
!                             in A-deck format
!
!     Output file: ocd5.dat - OCD5 forecast in ATCF format
!                  ocd5.log - OCD5 log file
!
!     Written by M. DeMaria, June 2013
!     Modified Oct 2016 for Cray  (MD)
!         - global land utilities used
!         - OCD5 intensity set to 15 kt if it is less than 15 kt but
!           track forecast is available
!
!     Modified May 2021 (MD)
!         - Bug fix for EP and CP crossing dateline
!
!     The program uses the following subroutines and related data files:
!        
!        decayo.f writeaidlocal2.f xint.f
!        gdland_table.f90 gfland_table.f (from liblandutils)
!        dataio.f dtgutils.f dataformats.inc dataioparams.inc (from libguidanceio)
!        gdland_table.dat, gfland_table.dat
!
      include 'dataformats.inc'
!
      character *2  bb,bbc,bbs
      character *4  mname
      character *8 strmid
      character *10 aymdh,aymdhc, aymdhs
!
      parameter (imax=11)
      character*1 latNS(imax),lonEW(imax)
      dimension ftime(imax),iftime(imax)
      dimension rlatc(imax),rlonc(imax),vmaxs(imax)
      dimension rlons(imax)
!
!     Variables for ATCF format output
      dimension ifocd5(imax,3)
      dimension iolat(0:imax-1),iolon(0:imax-1),iovmax(0:imax-1)
      character *4 stype(0:imax-1)
!
      dimension vmaxd(imax),dland(imax)
!
      type ( AID_DATA ) comRcd, tauData
!
      data lucom,luinp,ludat,lulog /10,11,12,13/
      data mname /'OCD5'/
!
!     **** Set up activities
!
!     Calculate forecast times
      idelt = 12
      do k=1,imax
         iftime(k) = idelt*(k-1)
         ftime(k)  = float(iftime(k))
      enddo
!
!     Initialize forecast arrays to zero
      do k=1,imax
         vmaxs(k) = 0.0
         rlatc(k) = 0.0
         rlonc(k) = 0.0
      enddo
!
!     **** Open the com, adeck input files and dat,log output files
      open(file='ocd5.log',unit=lulog,form='formatted',status='replace')
!
      write(lulog,499)
  499 format('Begin OCD5 forecast')
!
      open(file='ocd5.com',unit=lucom,form='formatted',status='old', &
      err=900)
!
      open(file='ocd5.inp',unit=luinp,form='formatted',status='old', &
      err=910)
!
      open(file='ocd5.dat',unit=ludat,form='formatted',status='replace')
!
!     **** Read t=0 info from CARQ line
!
      call getARecord (lucom,"CARQ", comRcd, istat)
      if (istat .eq. 0) then
        write (lulog, 970) istat
        stop
      end if
!
      call getSingleTAU ( comRcd, 0, tauData, istat )
      if (istat .ne. 1) then
        write (lulog, 970) istat
        stop
      end if
      aymdh   = tauData%aRecord(1)%DTG
      bb      = tauData%aRecord(1)%basin
      nn      = tauData%aRecord(1)%cyNum
      ivmax00 = tauData%aRecord(1)%vmax
      rlat00  = tauData%aRecord(1)%lat
      rlon00  = tauData%aRecord(1)%lon
      latNS(1)= tauData%aRecord(1)%NS
      lonEW(1)= tauData%aRecord(1)%EW
!
      write(lulog,502) bb,nn,aymdh(1:10)
  502 format(/,'CARQ input: ',/,a2,i2.2,' t=00 date/time: ',a10)
!
      write(lulog,504) ivmax00,rlat00,latNS(1),rlon00,lonEW(1)
  504 format(5x,'t=0 vmax, lat, lon: ',i3,1x,f5.1,a1,1x,f6.1,a1)
!
!     Initialize t=0 CLIPER,SHIFOR lat,lon,vmax with CARQ info
      if (latNS(1) .eq. 'S') then
         rlatc(1) = -1.0*rlat00
      else
         rlatc(1) =      rlat00
      endif

      if (lonEW(1) .eq. 'E') then
         rlonc(1) = 360.0 - rlon00
      else
         rlonc(1) = rlon00
      endif

      vmaxs(1) = float(ivmax00)
!
!     **** Read CLP5 track. There must be at least a 12  CLP5
!     forecast to make an OCD5 forecast
!
      do nrec=1,1000
        call getARecord (luinp,"CLP5", comRcd, istat)
        if (istat .eq. 0) then
          write (lulog, 980) istat
          stop
        end if
!
        call getSingleTAU( comRcd, 12, tauData, istat)
        if (istat .ne. 1) then 
          write (lulog, 980) istat
          stop
        end if
        aymdhc  = tauData%aRecord(1)%DTG
        bbc     = tauData%aRecord(1)%basin
        nnc     = tauData%aRecord(1)%cyNum
        if (aymdh .eq. aymdhc) exit
      enddo
!
      do k=2,imax
         call getSingleTAU( comRcd, iftime(k), tauData, istat)
         if (istat .ne. 1) exit
         rlatc(k) = tauData%aRecord(1)%lat
         rlonc(k) = tauData%aRecord(1)%lon
         latNS(k) = tauData%aRecord(1)%NS
         lonEW(k) = tauData%aRecord(1)%EW
         if (latNS(k) .eq. 'S') then
            rlatc(k) = -1.0*rlatc(k)
         endif

         if (lonEW(k) .eq. 'E') then
            rlonc(k) = 360.0-rlonc(k)
         endif

      enddo
!
      write(lulog,512) bbc,nnc,aymdhc(1:10)
  512 format(/,'CLP5 input: ',/,a2,i2.2,' t=12 data/time:',a10)
      do k=1,imax
         write(lulog,514) iftime(k),rlatc(k),rlonc(k)
  514    format('t=',i3,' CLP5 lat,lon: ',f5.1,1x,f6.1)
      enddo
!
!     **** Read SHF5 intensity. There must be at least at 12 h SHF5
!     forecast to make an OCD5 forecast
!
      call getARecord (luinp,"SHF5", comRcd, istat)
      if (istat .eq. 0) then
        write (lulog, 980) istat
        stop
      end if
!
      call getSingleTAU( comRcd, 12, tauData, istat)
      if (istat .ne. 1) then
        write (lulog, 980) istat
        stop
      end if
      aymdhs = tauData%aRecord(1)%DTG
      bbs    = tauData%aRecord(1)%basin
      nns    = tauData%aRecord(1)%cyNum
!
      do k=2,imax
         call getSingleTAU( comRcd, iftime(k), tauData, istat)
         if (istat .ne. 1) exit
         vmaxs(k) = float( tauData%aRecord(1)%vmax )
      enddo
 2000 continue
!
      write(lulog,522) bbs,nns,aymdhs(1:10)
  522 format(/,'SHF5 input: ',/,a2,i2.2,' t=12 data/time:',a10)
      do k=1,imax
         write(lulog,524) iftime(k),vmaxs(k)
  524    format('t=',i3,' SHF5 vmax: ',f5.0)
      enddo
!
!     **** Check consistency between CARQ, CLP5 and SHF5 forecasts
      if (aymdh .ne. aymdhc .or. aymdh .ne. aymdhs      ) then
        write (lulog, 990) aymdh, aymdhc, aymdhs, bb, bbc, bbs, nn, nnc, nns
        stop
      end if
!
      if (bb .ne. bbc .or. bb .ne. bbs      ) then
        write (lulog, 990) aymdh, aymdhc, aymdhs, bb, bbc, bbs, nn, nnc, nns
        stop
      end if
!
      if (nn .ne. nnc .or. nn .ne. nns      ) then
        write (lulog, 990) aymdh, aymdhc, aymdhs, bb, bbc, bbs, nn, nnc, nns
        stop
      end if
!
!     **** Determine the sign of the longitude and adjust input
!          to decay model accordingly (that used deg E pos, deg W neg)
      if (bb .eq. 'AL' .or. bb .eq. 'EP' .or. bb .eq. 'CP') then
         sgnlon = -1.0
      else
         sgnlon =  1.0
      endif
!
      isgnlon = nint(sgnlon)
!
      do k=1,imax
         rlons(k) = sgnlon*rlonc(k)
      enddo
!
!     **** Call the decay model with CLP5 track, SHF5 intensity
      call decayo(ftime,rlatc,rlons,vmaxs,vmaxd,dland,imax,lulog)
!
!     **** Check the OCD5 forecast for points with no valid intensity
!          but a valid track forecast. Set intensity to 15 kt for those times.
      vmaxdef = 15.0
      do k=1,imax
         if (abs(rlatc(k)) .gt. 0.0 .and. vmaxd(k) .lt. vmaxdef) then
            vmaxd(k) = vmaxdef
         endif
      enddo
!
      write(lulog,550)
  550 format(/,'OCD5 output:', /,'  t   lat    lon vmaxs vmaxd dland')
      do k=1,imax
         write(lulog,552) iftime(k),rlatc(k),rlons(k), &
                          vmaxs(k),vmaxd(k),dland(k)
  552    format(i3,1x,f5.1,1x,f6.1,1x, &
                f5.1,1x,f5.1,1x,f6.0)
      enddo
!
!     **** Write forecast in ATCF format
      do k=1,imax
         ifocd5(k,1) = nint(10.0*rlatc(k))
         ifocd5(k,2) = nint(10.0*rlonc(k))
         ifocd5(k,3) = nint(     vmaxd(k))
!
         iolat(k-1) =         nint(10.0*rlatc(k))
         iolon(k-1) = isgnlon*nint(10.0*rlonc(k))
         iovmax(k-1)=         nint(     vmaxd(k))
!
         stype(k-1) = '    '
      enddo
!
      write(strmid,558) bb,nn,aymdh(1:4)
  558 format(a2,i2.2,a4)
!
!     call newWriteAidRcd(ludat,strmid,aymdh,"OCD5",ifocd5)
!     call writeaidlocal(ludat,strmid,aymdh,mname,
!    +                   iolat,iolon,iovmax,stype,isgnlon,imax-1)
      ndt  = imax-1
      minc = 1
      call writeaidlocal2(ludat,strmid,aymdh,mname, &
                          iolat,iolon,iovmax,stype,iftime,ndt,minc)
!
      write(lulog,599)
  599 format(/,'OCD5 completed normally')
!
      stop
!
!     **** Error processing
  900 write(lulog,950) 
  950 format(/,'Error opening ocd5.com input file')
      stop
!
  910 write(lulog,960) 
  960 format(/,'Error opening ocd5.inp input file')
      stop
!
  920 write(lulog,970) istat
  970 format(/,'Error reading ocd5.com input file, istat=',i4)
      stop
!
  930 write(lulog,980) istat
  980 format(/,'Error reading ocd5.inp input file, istat=',i4)
      stop
!
  940 write(lulog,990) aymdh,aymdhc,aymdhs,bb,bbc,bbs,nn,nnc,nns
  990 format(/,'Date/time or storm ID not consistent:', &
             /,'aymdh,aymdhc,aymds: ',3(1x,a10), &
             /,'bb,      bbc,  bbs: ',3(9x,a2), &
             /,'nn,      nnc,  nns: ',3(7x,i4) )
      stop
!
      end