subroutine decay(ftime,rlat,rlon,vmax,vmaxa,dland,lulg) C C This routine adjusts a tropical cyclone intensity C forecast to account for decay over land. this version is C valid for the atlantic basin and was written by M. DeMaria C and J. Kaplan of the Hurricane Research Division, May 1994. C C This version was modified 4/10/97 (MDM) to include the C New England coefficients. The distance inland correction C is disabled in this version (idtlc=0). C C New version created 9/30/2004 that allows for decay proportional C to the fraction of the storm circulation over land. C (Set rcrad > 0.0 to activate this option). The logic of the code c was changed so accomodate this option. The old distance inland c correction was completely eliminated with this modification. C c Modified 12/19/2011 to include west Pacific option and to change c longitude input to deg E positive, deg W negative. The routines xint.f c and fland.f were removed from the source code so those need to be included c in the Makefile or compile script. Common block mparm was removed and save c statements were added. c c Modified Oct 2013 to replace dtland with global gbland routine (MD) c Modified Apr 2020 to update to use new distance to land and fraction of c land tables (SS) c Modified Jan 2023 for NCO specifications (KM) c c Required routines c xint.f c gdland_table.f90 c gfland_table.f (new version with deg E pos, deg W neg and global land routine) c c Required data files c gdland_table.dat c gfland_table.dat c C ********** INPUT ********** C C ftime: the time in hours (for example, 0.,12.,24. ... 72.) C The times need to sequential, but the time interval C does not need to be even. C rlat: The storm latitude (deg N) at the times in array ftime C rlon: The storm longitude (deg E positive, deg W neg) at the times in C array ftime C vmax: The storm maximum wind (kt) at the times in array ftime. C Set vmax=0 for missing forecast times. C lulg: Unit number for write statements C C ********** OUTPUT ********** C C vmaxa: The storm maximum wind (kt) adjusted for decay over land C at the times in array ftime. C C dland: The distance (km) from the storm center (rlat,rlon) to C the nearest major land mass. dland is negative if the C point is storm center is inland. C C ********** METHOD ********* C C The simple exponential decay model developed by M. DeMaria C and J. Kaplan at HRD is used to decay the storm intensity for C the portions of the track over land. C C In this version, the decay rate is proportional to the fraction C of the storm circualtion over land. C C ********** PARAMETER SPECIFICATION ********** C C Specify the maximum number of time values. parameter (imax=29) C C Specify the time interval (hr) for linearly interpolating C the track positions. data dt /1.00/ C C Set interp=1 to print out (to unit lulg) all intermediate C intensity calculations or else set interp=0 for no print data interp /1/ C C Specify decay model parameters C C Coefficients for east/gulf coast (or low latitude systems) data rf1,a1,vb1,rclat1 /0.9,0.095,26.7,36.0/ C C Coefficients for New England (or high latitude systems) data rf2,a2,vb2,rclat2 /0.9,0.183,29.6,40.0/ C C Specify radius of storm circulation (km) for fractional C decay option. Set rcrad to zero to eliminate this option. C data rcrad / 110.0/ c data rcrad / 0.0/ c c common /mparm/ alpha,vb,redfac save rf1,a1,vb1,rclat1 save rf2,a2,vb2,rclat2 save rcrad,dt,interp C C ********** DIMENSION ARRAYS ********** C dimension ftime(imax),rlat(imax),rlon(imax) dimension vmax(imax),vmaxa(imax),dland(imax) C c 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) c C ********** MODEL CODE ********* C C 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 C c Initialize vmaxa array to zero do i=1,imax vmaxa(i) = 0.0 enddo c C Find the number of valid forecast times itimet = 0 do i=1,imax if (vmax(i) .lt. 0.5) exit itimet=i enddo C C There must be at least two valid forecast times if (itimet .lt. 2) return C C Check to make sure times are sequential itime=0 do i=2,itimet if (ftime(i) .le. ftime(i-1)) exit itime=i enddo C if (itime .lt. 2) return c if (interp .gt. 2) then do i=1,itime write(lulg,887) ftime(i),rlat(i),rlon(i),vmax(i) enddo endif c c 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 c c Interpolate the input lat,lon and max winds to the c small time interval c c ++Find vmax on small time grid iflag=1 lflag=0 xi = 0.0 c call xint(ftime,vmax,itime,iflag,lflag,xi,fi,ierr) c iflag=0 do i=1,ntimes call xint(ftime,vmax,itime,iflag,lflag, + ftimes(i),vmaxs(i),ierr) enddo c c ++Find lat on small time grid iflag=1 call xint(ftime,rlat,itime,iflag,lflag,xi,fi,ierr) c iflag=0 do i=1,ntimes call xint(ftime,rlat,itime,iflag,lflag, + ftimes(i),rlats(i),ierr) enddo c c ++Find lon on small time grid iflag=1 call xint(ftime,rlon,itime,iflag,lflag,xi,fi,ierr) c iflag=0 do i=1,ntimes call xint(ftime,rlon,itime,iflag,lflag, + ftimes(i),rlons(i),ierr) enddo c c 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 c c 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 c vmaxas(1) = vmaxs(1) c do i=2,ntimes c At each step in this loop, the decay model is integrated c from t=ftimes(i-1) to t=ftimes(i) c c 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) C redfac = w1*rf1 + w2*rf2 alpha = w1*a1 + w2*a2 vb = w1*vb1 + w2*vb2 endif C vmaxt1 = vmaxas(i-1) c if (dlands(i) .ge. 0.0) then c ++ This is an over-water point c c Check to see if storm just moved over water. c If so, adjust for land/ocean surface roughness differences if (dlands(i-1) .lt. 0.) then vmaxt1 = vmaxt1/redfac endif c vmaxt2 = vmaxt1 + (vmaxs(i)-vmaxs(i-1)) else c ++ This is an over-land point c c Check to see if storm just moved over land. c If so, adjust for ocean/land surface roughness differences if (dlands(i-1) .ge. 0.) then vmaxt1 = redfac*vmaxt1 endif c t = ftimes(i)-ftimes(i-1) fbar = 0.5*(flands(i)+flands(i-1)) vmaxt2 = vb + (vmaxt1-vb)* exp(-fbar*alpha*t) endif c vmaxas(i) = vmaxt2 enddo c if (interp .gt. 2) then do i=1,ntimes write(lulg,887) ftimes(i),vmaxs(i),vmaxas(i), + rlats(i),rlons(i),dlands(i),flands(i) 887 format(8(f6.1,1x)) enddo endif c c Interpolate decay vmaxas back to original forecast times iflag=1 c call xint(ftimes,vmaxas,ntimes,iflag,lflag,xi,fi,ierr) c iflag=0 do i=1,itime call xint(ftimes,vmaxas,ntimes,iflag,lflag, + ftime(i),vmaxa(i),ierr) enddo c c Interpolate dlands back to original forecast times iflag=1 c call xint(ftimes,dlands,ntimes,iflag,lflag,xi,fi,ierr) c iflag=0 do i=1,itime call xint(ftimes,dlands,ntimes,iflag,lflag, + ftime(i),dland(i),ierr) enddo c return end