subroutine decayo(ftime,rlat,rlon,vmax,vmaxa,dland,imax,lulg)
!
!     This routine adjusts a tropical cyclone intensity 
!     forecast to account for decay over land. this version is
!     valid for the atlantic basin and was written by M. DeMaria
!     and J. Kaplan of the Hurricane Research Division, May 1994.
!
!     This version was modified 4/10/97 (MDM) to include the 
!     New England coefficients. The distance inland correction
!     is disabled in this version (idtlc=0). 
!
!     New version created 9/30/2004 that allows for decay proportional
!     to the fraction of the storm circulation over land.
!     (Set rcrad > 0.0 to activate this option). The logic of the code
!     was changed so accomodate this option. The old distance inland
!     correction was completely eliminated with this modification. 
!
!     Modified 12/19/2011 to include west Pacific option and to change
!     longitude input to deg E positive, deg W negative. The routines xint.f
!     and fland.f were removed from the source code so those need to be included
!     in the Makefile or compile script. Common block mparm was removed and save
!     statements were added.
! 
!     Modified 06/02/2013 to make imax an input parameter for use by OCD5
!
!     Modified 10/13/2016 to use global land subroutine instead of dtland (MD)
!     Modified 04/06/2020 to use global land tables (SS)

!     Required routines
!       xint.f
!       gdland_table.f90
!       gfland_table.f
!
!     Required data files
!       gdland_table.dat
!       gfland_table.dat
!
!     ********** INPUT **********
! 
!       ftime: the time in hours (for example, 0.,12.,24. ... 72.)
!              The times need to sequential, but the time interval
!              does not need to be even. 
!       rlat:  The storm latitude (deg N) at the times in array ftime
!       rlon:  The storm longitude (deg E positive, deg W neg) at the times in
!              array ftime
!       vmax:  The storm maximum wind (kt) at the times in array ftime.
!              Set vmax=0 for missing forecast times. 
!       imax:  The maximum dimension of ftime, rlat, rlon, vmax, vmaxa and dland 
!       lulg:  Unit number for write statements
!
!     ********** OUTPUT **********  
!
!       vmaxa: The storm maximum wind (kt) adjusted for decay over land
!              at the times in array ftime.
!
!       dland: The distance (km) from the storm center (rlat,rlon) to 
!              the nearest major land mass. dland is negative if the 
!              point is storm center is inland.
!
!     ********** METHOD *********
!
!     The simple exponential decay model developed by M. DeMaria 
!     and J. Kaplan at HRD is used to decay the storm intensity for
!     the portions of the track over land.
!
!     In this version, the decay rate is proportional to the fraction
!     of the storm circualtion over land. 
!
!     ********** PARAMETER SPECIFICATION **********
!
!     Specify the maximum number of time values 
!     (now an input parameter in this version)
!      parameter (imax=21)
!
!     Specify the time interval (hr) for linearly interpolating 
!     the track positions.
      data dt /1.00/
!
!     Set interp=1 to print out (to unit lulg) all intermediate 
!     intensity calculations or else set interp=0 for no print 
      data interp /1/ 
!
!     Specify decay model parameters
!
!     Coefficients for east/gulf coast (or low latitude systems)
      data rf1,a1,vb1,rclat1 /0.9,0.095,26.7,36.0/
!     
!     Coefficients for New England (or high latitude systems)
      data rf2,a2,vb2,rclat2 /0.9,0.183,29.6,40.0/
!
!     Specify radius of storm circulation (km) for fractional
!     decay option. Set rcrad to zero to eliminate this option. 
!
      data rcrad / 110.0/
!     data rcrad /   0.0/
!
!     common /mparm/ alpha,vb,redfac
      save rf1,a1,vb1,rclat1
      save rf2,a2,vb2,rclat2
      save rcrad,dt,interp
!
!     ********** DIMENSION ARRAYS **********
!
      dimension ftime(imax),rlat(imax),rlon(imax)
      dimension vmax(imax),vmaxa(imax),dland(imax)
!
!     Arrays for small time step
      parameter (imaxs=1000)
      dimension ftimes(imaxs),rlats(imaxs),rlons(imaxs)
      dimension vmaxs(imaxs),vmaxas(imaxs)
      dimension dlands(imaxs),flands(imaxs)
! 
!     ********** MODEL CODE *********
!
!     Write model message to log file
      if (interp .gt. 0) then
         write(lulg,810) rcrad
  810    format(' Decay model with rcrad= ',f5.0,' km')
      endif
!
!     Initialize vmaxa array to zero
      do i=1,imax
	 vmaxa(i) = 0.0
      enddo
! 
!     Find the number of valid forecast times
      itimet = 0
      do i=1,imax
         if (vmax(i) .lt. 0.5) exit
         itimet=i 
      enddo
!
!     There must be at least two valid forecast times
      if (itimet .lt. 2) return
!
!     Check to make sure times are sequential
      itime=0
      do i=2,itimet
         if (ftime(i) .le. ftime(i-1)) exit
         itime=i
      enddo
!
      if (itime .lt. 2) return
!         
      if (interp .gt. 2) then
         do i=1,itime
            write(lulg,887) ftime(i),rlat(i),rlon(i),vmax(i)
         enddo
      endif
!
!     Calcuate the time values at the small time interval points
      ntimes = 1 + (ftime(itime)-ftime(1))/dt
      do i=1,ntimes
	 ftimes(i) = ftime(1) + dt*float(i-1)
      enddo
!
!     Interpolate the input lat,lon and max winds to the 
!     small time interval
!
!       ++Find vmax on small time grid
      iflag=1
      lflag=0
      xi = 0.0
! 
      call xint(ftime,vmax,itime,iflag,lflag,xi,fi,ierr)
!
      iflag=0
      do i=1,ntimes
         call xint(ftime,vmax,itime,iflag,lflag, &
                        ftimes(i),vmaxs(i),ierr)
      enddo
!
!       ++Find lat on small time grid
      iflag=1
      call xint(ftime,rlat,itime,iflag,lflag,xi,fi,ierr)
!
      iflag=0
      do i=1,ntimes
         call xint(ftime,rlat,itime,iflag,lflag, &
                        ftimes(i),rlats(i),ierr)
      enddo
!
!       ++Find lon on small time grid
      iflag=1
      call xint(ftime,rlon,itime,iflag,lflag,xi,fi,ierr)
!
      iflag=0
      do i=1,ntimes
         call xint(ftime,rlon,itime,iflag,lflag, &
                        ftimes(i),rlons(i),ierr)
      enddo
!
!     Calcuate distance to land and fractional land at small time points
      do i=1,ntimes
         call gdland_table(rlons(i),rlats(i),dlands(i))
      enddo
!
!     Integrate the decay model over the small time points
      do i=1,ntimes
	 if (rcrad .gt. 0.0) then
	    call gfland_table(rlons(i),rlats(i),flands(i))
         else
	    flands(i) = 1.0
         endif
      enddo
!
      vmaxas(1) = vmaxs(1)
!
      do i=2,ntimes
!        At each step in this loop, the decay model is integrated
!        from t=ftimes(i-1) to t=ftimes(i)
!
!        Calculate decay model parameters at current latitude
         rlatt = abs(rlats(i-1))
         if     (rlatt .ge. rclat2) then
             redfac = rf2
             alpha  = a2
             vb     = vb2
         elseif (rlatt .le. rclat1) then 
             redfac = rf1
             alpha  = a1
             vb     = vb1
         else
             w1 = (rclat2-rlatt)/(rclat2-rclat1)
             w2 = (rlatt-rclat1)/(rclat2-rclat1)
!
             redfac = w1*rf1 + w2*rf2
             alpha  = w1*a1  + w2*a2
             vb     = w1*vb1 + w2*vb2
         endif
!
        vmaxt1 = vmaxas(i-1)
!
	if (dlands(i) .ge. 0.0) then
!          ++ This is an over-water point
!
!          Check to see if storm just moved over water.
!          If so, adjust for land/ocean surface roughness differences
           if (dlands(i-1) .lt. 0.) then
	      vmaxt1 = vmaxt1/redfac
           endif
!
	   vmaxt2 = vmaxt1 + (vmaxs(i)-vmaxs(i-1))
         else
!          ++ This is an over-land point
!
!          Check to see if storm just moved over land.
!          If so, adjust for ocean/land surface roughness differences
           if (dlands(i-1) .ge. 0.) then
	      vmaxt1 = redfac*vmaxt1
           endif
!
           t      = ftimes(i)-ftimes(i-1)
           fbar   = 0.5*(flands(i)+flands(i-1))
	   vmaxt2 = vb + (vmaxt1-vb)* exp(-fbar*alpha*t)
	 endif
!
	 vmaxas(i) = vmaxt2
      enddo
!
      if (interp .gt. 2) then
         do i=1,ntimes
	    write(6,887) ftimes(i),vmaxs(i),vmaxas(i),rlats(i),rlons(i), &
                         dlands(i),flands(i)
  887       format(8(f6.1,1x))
         enddo
      endif
!
!     Interpolate decay vmaxas back to original forecast times
      iflag=1
! 
      call xint(ftimes,vmaxas,ntimes,iflag,lflag,xi,fi,ierr)
!
      iflag=0
      do i=1,itime
         call xint(ftimes,vmaxas,ntimes,iflag,lflag, &
                        ftime(i),vmaxa(i),ierr)
      enddo
!
!     Interpolate dlands back to original forecast times
      iflag=1
! 
      call xint(ftimes,dlands,ntimes,iflag,lflag,xi,fi,ierr)
!
      iflag=0
      do i=1,itime
         call xint(ftimes,dlands,ntimes,iflag,lflag, &
                        ftime(i),dland(i),ierr)
      enddo
!
      return
      end