subroutine lgeim(nft,ftime,rlat,rlon,vmpi,cappa,v0, + rmm,rnn,beta,vscale,lulg,vmax,dland,ierr) c c This code is for an intensity forecast model based upon c a simple two term logistical growth equation method. A separate procedure is c used if the storm is over land. c c Modified Oct 2013 for global distance to land function (MD) c and SH cases c c Modified Apr 2020 for global distance to land (gdland_table) and c fraction of land (gfland_table) tables c c Modified Jan 2023 for NCO specifications (KM) c c ********** INPUT ********** c c nft: The number of forecast times c ftime(nft ): The time in hours (for example, 0.,6.,12. ... 120.) c The times need to sequential, but the time interval c does not need to be even. c rlat(nft): The storm latitude (deg N) at the times in array ftime c rlon(nft): The storm longitude (deg W negative or deg E positive) c at the times in array fitme c vmpi(nft): The maximum potential intensity (kt) versus time c cappa(nft): The coefficient of the synoptic term (1/hr) in the simple two term model c versus time c rmm: The order of the growth term (usually 1.1 > rmm > 0) c rnn: The order of the SST term (usually 2 or 4) c beta: The time scale (1/hr) of the SST term c vscale: The velocity scale (kt) to make growth term non-dimensional c v0: The initial maximum winds (kt) c lulg: Unit number for write statements c c ********** OUTPUT ********** c c vmax(nft): The forecast maximum wind (kt) versus time c dland(nft): 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 ierr: Error flag (=0 for normal return, =1 for error) c c ********** PARAMETER SPECIFICATION ********** c c Specify the time interval (hr) for numerical integration data dt /1.00/ c c Set interp=1,2 or 3 to print out (to unit lulg) intermediate c intensity calculations or else set interp=0 for no print data interp /0/ c c Inland decay model coefficients for east/gulf coast data rf1,a1,vb1,rclat1 /0.9,0.095,26.7,36.0/ c c Inland decay model coefficients for New England 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 c data rcrad / 0.0/ data rcrad / 110.0/ c data rcrad / 110.0/ c c ********** Passed arrays ********** c dimension ftime(nft),rlat(nft),rlon(nft) dimension vmpi(nft),cappa(nft),dland(nft) dimension vmax(nft) c c ********** Arrays for numerical integration ********** parameter (imaxs=1000) dimension ftimes(imaxs),rlats(imaxs),rlons(imaxs) dimension vmpis(imaxs),cappas(imaxs) dimension vmaxs(imaxs) dimension redfacs(imaxs),vbs(imaxs),alphas(imaxs) dimension dlands(imaxs),flands(imaxs) c c ********** MODEL CODE ********* c c Initial intensity forecast vmax(1) = v0 do i=2,nft vmax(i) = 0.0 enddo c c Find the number of valid forecast times itimet = 0 do i=1,nft if (abs(rlat(i)) .lt. 0.5) exit if (abs(rlat(i)) .gt. 90.0) exit if (vmpi(i) .lt. 0.5) exit if (abs(cappa(i)) .gt. 5.0) exit itimet=i enddo c c There must be at least two valid forecast times if (itimet .lt. 2) then ierr=1 return endif 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) then ierr=1 return endif c if (interp .gt. 2) then do i=1,itime write(6,887) ftime(i),rlat(i),rlon(i),cappa(i),vmpi(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,vmpi,cappa to the c small time interval c c ++ lat lflag=0 iflag=1 call xint(ftime,rlat,itime,iflag,lflag,xi,fi,ierrx) c iflag=0 do i=1,ntimes call xint(ftime,rlat,itime,iflag,lflag, + ftimes(i),rlats(i),ierrx) enddo c c ++ lon iflag=1 call xint(ftime,rlon,itime,iflag,lflag,xi,fi,ierrx) c iflag=0 do i=1,ntimes call xint(ftime,rlon,itime,iflag,lflag, + ftimes(i),rlons(i),ierrx) enddo c c ++ vmpi iflag=1 lflag=0 xi = 0.0 c call xint(ftime,vmpi,itime,iflag,lflag,xi,fi,ierrx) c iflag=0 do i=1,ntimes call xint(ftime,vmpi,itime,iflag,lflag, + ftimes(i),vmpis(i),ierrx) enddo c c ++ cappa iflag=1 lflag=0 xi = 0.0 c call xint(ftime,cappa,itime,iflag,lflag,xi,fi,ierrx) c iflag=0 do i=1,ntimes call xint(ftime,cappa,itime,iflag,lflag, + ftimes(i),cappas(i),ierrx) 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)) call gfland_table(rlons(i),rlats(i),flands(i)) enddo c c Calculate decay model parameters at small time points do i=1,ntimes rlatt = abs(rlats(i)) if (rlatt .ge. rclat2) then redfacs(i) = rf2 alphas(i) = a2 vbs(i) = vb2 elseif (rlatt .le. rclat1) then redfacs(i) = rf1 alphas(i) = a1 vbs(i) = vb1 else w1 = (rclat2-rlatt)/(rclat2-rclat1) w2 = (rlatt-rclat1)/(rclat2-rclat1) c redfacs(i) = w1*rf1 + w2*rf2 alphas(i) = w1*a1 + w2*a2 vbs(i) = w1*vb1 + w2*vb2 endif enddo c c Perform time integration vmaxs(1) = v0 c do 99 i=1,ntimes-1 c vnow = vmaxs(i) c c Check to see if the storm is over water or land if (dlands(i) .lt. 0.0) then c Over land case c if (i .gt. 1) then c Check to see if the storm just moved over land c If so, apply the sea/land reduction factor if (dlands(i-1) .ge. 0.0) then vnow = vnow*redfacs(i) endif endif c c Forward time step vmaxs(i+1) = vnow - alphas(i)*flands(i)*(vnow-vbs(i))*dt else c c Over water case if (i .gt. 1) then c Check to see if the storm just moved over water c If so, apply the inverse sea/land reduction factor if (dlands(i-1) .lt. 0.0) then vnow = vnow/redfacs(i) endif endif c c Forward time step c c Calculate coefficient of MPI term cmpi = beta*(vnow/vmpis(i))**rnn vgtm = (vnow/vscale)**rmm vmaxs(i+1) = vnow + ( cappas(i)*vgtm - cmpi*vnow )*dt endif 99 continue c if (interp .gt. 2) then do i=1,ntimes write(6,887) ftimes(i),rlats(i),rlons(i),vmpis(i),cappas(i), + redfacs(i),vbs(i),alphas(i),dlands(i), + flands(i),vmaxs(i) 887 format(4(f6.1,1x),1x,f8.4,1x,f6.1,1x,f6.1,1x, + f8.4,1x,f6.0,1x,f8.4,1x,f6.1) enddo endif c c Interpolate decay vmaxs back to original forecast times iflag=1 c call xint(ftimes,vmaxs,ntimes,iflag,lflag,xi,fi,ierrx) c iflag=0 do i=2,itime call xint(ftimes,vmax,ntimes,iflag,lflag, + ftime(i),vmax(i),ierrx) 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),ierrx) enddo c return end