subroutine tclip(rlon00,rlat00,rlonm12,rlatm12,cx00,cy00,ipf, + vmx00,vmxm12,iyr,imon,iday,itime,dthr,ndt,ioper, + iftypet,iftypei,flon,flat,fvmax,ierr) c c This routine makes a climatology and persistence forecast by following c a trajectory where the horizontal motion is a weighted sum of c a climatological storm motion vector and persistence. The intensity forecast c is from a version of LGEM with climatological input. c c Input: rlon00 - Longitude at t=0 (0 to 360 deg) c rlat00 - Latitude at t=0 (-90 to 90 deg) c rlonm12 - same as rlon00 at t=-12 hr (not used if ipf=1) c rlatm12 - same as rlat 00 at t=-12 hr (not used if ipf=1) c cx00 - eastward component of storm motion at t=0 (kt) (not used if ipf=0) c cy00 - northward component of storm motion at t=0 (kt) (not used if ipf=0) c ipf - Method for calculation initial motion, c =1 to use cx00,cy00 or c =0 to calculate it from rlon,rlat at t=0 and -12 hr c vmx00 - Maximum wind at t=0 (kt) c vmxm12 - Maximum wind at t=-12 hr (kt) c iyr - 4 digit year of t=0 c imon - 2 digit month of t=0 c iday - 2 digit day of t=0 c itme - 2 digit time of t=0 c dthr - forecast output time step (hr) c ndt - Number of time steps c ioper - Flag for adding path to file open on NCEP operational computers c iftypet - =0 for persistence only, c =1 for climatology only c =2 for climatology and persistence, for track forecast c iftypei - =0 for persistence only, c =1 for climatology only c =2 for climatology and persistence, for intensity forecast c c Output: flon(0:ndt) - Forecast longitudes (0 to 360 deg) c flat(0:ndt) - Forecast latitudes (-90 to 90 deg) c fvmax(0:ndt)- Forecast vmax (kt) c ierr - Error variable (=0 for normal completion) c (=1 if climo files are missing) c Required data files: c oban_WH_1982_2018.dat - climatological motion and growth rate values c clim_rsst_1982_2018.dat - climatological SST files c c Modified Oct 2013 to replace dtland calls with global gbland calls c Modified Oct 2016 (MD) for Cray c Modified Apr 2020 (SS) for new land routine tables c Modified May 2020 (SS/RZ) for updated model configuration and climo files c Modified Jan 2023 (KM) to remove goto statements; create error handling subroutine c c Passed variables dimension flat(0:ndt),flon(0:ndt),fvmax(0:ndt) c c Internal variables parameter (mdt=2400) dimension flons(0:mdt),flats(0:mdt),vmaxs(0:mdt) dimension cx(0:mdt),cy(0:mdt) dimension cxc(0:mdt),cyc(0:mdt) dimension cxp(0:mdt),cyp(0:mdt) dimension wxc(0:mdt),wyc(0:mdt) dimension wxp(0:mdt),wyp(0:mdt) dimension dtlnd(0:mdt),flnd(0:mdt),sstc(0:mdt),vmpi(0:mdt) dimension tf(0:mdt) dimension iyrf(0:mdt),imonf(0:mdt),idayf(0:mdt),itimef(0:mdt) c dimension cappa(0:mdt),cappap(0:mdt),cappac(0:mdt) dimension wcc(0:mdt),wcp(0:mdt) c parameter (mxc=360,myc=181,ntc=12) c dimension uclim(mxc,myc,ntc),vclim(mxc,myc,ntc),cclim(mxc,myc,ntc) dimension sstclim(mxc,myc,ntc) dimension ndmo(ntc) c common /climo/ uclim,vclim,cclim,glon1,glat1,dlon,dlat,nxc,nyc common /climoo/ sstclim,glono1,glato1,dlono,dlato,nxco,nyco common /climot/ ndmo c common /pncon1/ pi,dtr,ckt2ms,deg2m,eradkm,eradm,rmiss common /pncon2/ rf1,rf2,a1,a2,vb1,vb2,rclat1,rclat2,rcrad common /pncon3/ auc,avc,efuc,efvc,aup,avp common /pncon4/ efcp,acc,acp,rnn,beta c data iclimo /1/ c c Set the maximum time step (hr) to use for the integration dtmax = 1.0 c c Set mpired=1 to linearly reduce MPI to vmpizer as SST approaches sstzer, starting c at sstmin mpired=1 sstmin=20.0 sstzer =10.0 vmpizer=1.0 c c Set minimum intensity (kt) vmin = 15.0 c c Set maximum time (hr) over land before dissipation tolmax = 72.0 c c Set mpispd=1 to add a fraction of the translational speed to the mpi mpispd=1 c c Set unit number (lulog) for debug write statements and debug write flag lulog = 6 idbug = 0 c c Initialize forecast to zero do k=0,ndt flat(k) = 0.0 flon(k) = 0.0 fvmax(k)= 0.0 enddo c do i=0,mdt flons(i) = 0.0 flats(i) = 0.0 vmaxs(i) = 0.0 enddo c call pninit c c Initialize the position flon(0) = rlon00 flat(0) = rlat00 flons(0) = rlon00 flats(0) = rlat00 vmaxs(0) = vmx00 fvmax(0) = vmx00 c c Determine which basin the storm starts in call bassel(flat(0),flon(0),ibasin) c ierr = 0 if (iclimo .eq. 1) then c Get climatological u,v,cappa and SST fields call getclimo(ioper,ierrc) if (ierrc .ne. 0) then ierr=ierrc return endif c glon2 = glon1 + dlon*float(nxc-1) glat2 = glat1 + dlat*float(nyc-1) endif c c write(6,750) glon1,glon2,dlon,nxc,glat1,glat2,dlat,nyc c 750 format(' glon1,glon2,dlon,nxc:',3(f6.1,1x),i3, c + ' glat1,glat2,dlat,nyc:',3(f6.1,1x),i3) c c Initialize numerical and physical constants c Calculate initial storm motion vector if necessary and c convert it to m/s if (ipf .eq. 0) then cfac = cos(0.5*dtr*(rlat00+rlatm12)) dlon00 = (rlon00-rlonm12)*cfac dlat00 = (rlat00-rlatm12) c c Check to see if Greenwich meridian was crossed if (dlon00 .gt. 180.0) dlon00 = dlon00 - 360.0 if (dlon00 .lt. -180.0) dlon00 = dlon00 + 360.0 c cxp(0) = (dlon00*deg2m)/(12*3600.0) cyp(0) = (dlat00*deg2m)/(12*3600.0) else cxp(0) = cx00*ckt2ms cyp(0) = cy00*ckt2ms endif c c Compare dtmax and output time step and compute the c interval for saving the output (intsave) and the c total number of time steps to integrate the model if (dthr .lt. dtmax) then intsave = 1 ndts = ndt dts = dthr else intsave = ifix(0.1 + dthr/dtmax) ndts = ndt*intsave dts = dthr/float(intsave) endif dtss = 3600.0*dts c c c Calcuate time (hr) at each time step, c the date/time to the nearest hour at each time step, c and persistence vector at each time step tf(0) = 0.0 iyrf(0) = iyr imonf(0) = imon idayf(0) = iday itimef(0) = itime c do k=1,ndts tf(k) = tf(0) + dts*float(k) c itadd = ifix(tf(k) + 0.49) call tadd(iyr,imon,iday,itime,itadd, + iyrf(k),imonf(k),idayf(k),itimef(k)) cxp(k) = cxp(0) cyp(k) = cyp(0) enddo c if (idbug .eq. 1) then write(lulog,700) cxp(0),cyp(0),dthr,dts,ndts 700 format('Initial motion vector (m/s): ',f7.2,1x,f7.2, + /,'Output time step=',f6.1, + /,'Integration time step=',f6.1, + /,'No. of time steps= ',i6) endif c c *** Forward time step c ++ Track call clsst(flons(0),flats(0),imonf(0),idayf(0),sstc(0)) call uvcal(flons(0),flats(0),vmaxs(0),tf(0), + imonf(0),idayf(0),iftypet, + cxp(0),cyp(0),cxc(0),cyc(0), + wxp(0),wyp(0),wxc(0),wyc(0), + cx(0),cy(0)) c fxm1 = cx(0)/(dtr*eradm*cos(dtr*flats(0))) fym1 = cy(0)/(dtr*eradm) c flons(1) = flons(0) + dtss*fxm1 flats(1) = flats(0) + dtss*fym1 c c ++ Intensity call gdland_table(flons(0),flats(0),dtlnd(0)) call gfland_table(flons(0),flats(0),flnd(0)) call mpicalo(sstc(0),vmpi(0),ibasin) call mpicalo(sstmin,vmpimin,ibasin) if (mpired .eq. 1) call rmpi(sstmin,vmpimin,sstzer,vmpizer, + sstc(0),vmpi(0)) if (mpispd .eq. 1) call ampi(cx(0),cy(0),vmpi(0)) c c Calculate initial growth rate from persistence call capinit(rlon00,rlat00,vmx00,rlonm12,rlatm12,vmxm12, + vmpi(0),iyr,imon,iday,cappa00) c do k=0,ndts cappap(k) = cappa00 cappac(k) = 0.0 cappa(k) = 0.0 wcp(k) = 0.0 wcc(k) = 0.0 enddo c call capcal(flons(0),flats(0),tf(0),imonf(0),idayf(0), + iftypei,cappap(0),cappac(0),wcp(0),wcc(0), + cappa(0)) c fvm0 = cappa(0)*vmaxs(0) - beta*vmaxs(0)*(vmaxs(0)/vmpi(0))**rnn c if (dtlnd(0) .lt. 0.0) then c Use inland decay model for intensity change call dparm(flats(0),redfac,alpha,vb) vmaxs(1) = vmaxs(0) - dts*alpha*flnd(0)*(vmaxs(0)-vb1) else vmaxs(1) = vmaxs(0) + dts*fvm0 endif c if (idbug .eq. 1) then call tsprint(tf(0),flons(0),flats(0),vmaxs(0), + cxp(0),cyp(0),cxc(0),cyc(0), + wxp(0),wyp(0),wxc(0),wyc(0), + cx(0),cy(0),dtlnd(0),sstc(0),vmpi(0), + cappap(0),cappac(0),cappa(0), + wcp(0),wcc(0),lulog) endif c c ++ Prepare for next time step and print if necessary call clsst(flons(1),flats(1),imonf(1),idayf(1),sstc(1)) call uvcal(flons(1),flats(1),vmaxs(1),tf(1), + imonf(1),idayf(1),iftypet, + cxp(1),cyp(1),cxc(1),cyc(1), + wxp(1),wyp(1),wxc(1),wyc(1), + cx(1),cy(1)) c fxm0 = cx(1)/(dtr*eradm*cos(dtr*flats(1))) fym0 = cy(1)/(dtr*eradm) c call gdland_table(flons(1),flats(1),dtlnd(1)) call gfland_table(flons(1),flats(1),flnd(1)) call mpicalo(sstc(1),vmpi(1),ibasin) if (mpired .eq. 1) call rmpi(sstmin,vmpimin,sstzer,vmpizer, + sstc(1),vmpi(1)) if (mpispd .eq. 1) call ampi(cx(1),cy(1),vmpi(1)) c call capcal(flons(1),flats(1),tf(1),imonf(1),idayf(1), + iftypei,cappap(1),cappac(1),wcp(1),wcc(1), + cappa(1)) fvm0 = cappa(1)*vmaxs(1) - beta*vmaxs(1)*(vmaxs(1)/vmpi(1))**rnn c if (idbug .eq. 1) then call tsprint(tf(1),flons(1),flats(1),vmaxs(1), + cxp(1),cyp(1),cxc(1),cyc(1), + wxp(1),wyp(1),wxc(1),wyc(1), + cx(1),cy(1),dtlnd(1),sstc(1),vmpi(1), + cappap(1),cappac(1),cappa(1), + wcp(1),wcc(1),lulog) endif c c *** Adams-Bashforth time steps for track, forward for intensity, to desired time c c Set time over land counter toland = 0.0 c do k=2,ndts c c ++ Track flons(k) = flons(k-1) - 0.5*dtss*fxm1 + 1.5*dtss*fxm0 flats(k) = flats(k-1) - 0.5*dtss*fym1 + 1.5*dtss*fym0 c fxm1 = fxm0 fym1 = fym0 c c ++ Intensity if (dtlnd(k-1) .lt. 0.0) then c Use inland decay model for intensity change call dparm(flats(k-1),redfac,alpha,vb) c c Check to see if storm just moved from water to land. c If so, multiple surface roughness adjustment parameter if (dtlnd(k-1) .lt. 0.0 .and. dtlnd(k-2) .ge. 0.0) then vmaxs(k-1) = vmaxs(k-1)*redfac endif c vmaxs(k) = vmaxs(k-1) - dts*alpha*flnd(k-1)*(vmaxs(k-1)-vb1) c toland = toland + dts if (toland .gt. tolmax) exit else c Check to see if storm just moved from water to land. c If so, multiple surface roughness adjustment parameter if (dtlnd(k-1) .ge. 0.0 .and. dtlnd(k-2) .lt. 0.0) then vmaxs(k-1) = vmaxs(k-1)/redfac endif c c Use climatological version of LGEM for intensity change vmaxs(k) = vmaxs(k-1) + dts*fvm0 c c Reset time of land counter since storm is over water again toland = 0.0 c c Check to see if vmax is below minimum intensity if (vmaxs(k) .lt. vmin) exit endif c c Check to see if storm is out of climo domain if (flons(k) .lt. glon1 .or. flons(k) .gt. glon2) exit if (flats(k) .lt. glat1 .or. flats(k) .gt. glat2) exit c c ++ Prepare for next time step and print if necessary call clsst(flons(k),flats(k),imonf(k),idayf(k),sstc(k)) call uvcal(flons(k),flats(k),vmaxs(k),tf(k), + imonf(k),idayf(k),iftypet, + cxp(k),cyp(k),cxc(k),cyc(k), + wxp(k),wyp(k),wxc(k),wyc(k), + cx(k),cy(k)) c fxm0 = cx(k)/(dtr*eradm*cos(dtr*flats(k))) fym0 = cy(k)/(dtr*eradm) c call gdland_table(flons(k),flats(k),dtlnd(k)) call gfland_table(flons(k),flats(k),flnd(k)) call mpicalo(sstc(k),vmpi(k),ibasin) if (mpired .eq. 1) call rmpi(sstmin,vmpimin,sstzer,vmpizer, + sstc(k),vmpi(k)) if (mpispd .eq. 1) call ampi(cx(k),cy(k),vmpi(k)) c call capcal(flons(k),flats(k),tf(k),imonf(k),idayf(k), + iftypei,cappap(k),cappac(k),wcp(k),wcc(k), + cappa(k)) fvm0 = cappa(k)*vmaxs(k) + -beta*vmaxs(k)*(vmaxs(k)/vmpi(k))**rnn c if (idbug .eq. 1) then call tsprint(tf(k),flons(k),flats(k),vmaxs(k), + cxp(k),cyp(k),cxc(k),cyc(k), + wxp(k),wyp(k),wxc(k),wyc(k), + cx(k),cy(k),dtlnd(k),sstc(k),vmpi(k), + cappap(k),cappac(k),cappa(k), + wcp(k),wcc(k),lulog) endif enddo c c *** Extract forecast at output time steps kk = 0 do k=1,ndts if (mod(k,intsave) .eq. 0) then kk = kk + 1 flon(kk) = flons(k) flat(kk) = flats(k) fvmax(kk)= vmaxs(k) endif enddo c return end subroutine pninit c This routine initializes needed constants c common /pncon1/ pi,dtr,ckt2ms,deg2m,eradkm,eradm,rmiss common /pncon2/ rf1,rf2,a1,a2,vb1,vb2,rclat1,rclat2,rcrad common /pncon3/ aup,avp,efup,efvp,auc,avc common /pncon4/ efcp,acc,acp,rnn,beta c c Numerical constants pi = 3.14159265 dtr = pi/180.0 c c Physical constants eradkm = 6371.0 eradm = eradkm*1000.0 c c Conversion factors deg2m = 2.0*pi*eradm/360.0 ckt2ms = 1.0/1.944 c c Missing data flag rmiss = -999. c c Parameters for inland decay model rf1 = 0.9 rf2 = 0.9 a1 = 0.095 a2 = 0.183 vb1 = 26.7 vb2 = 29.6 rclat1 = 36.0 rclat2 = 40.0 rcrad = 110.0 c c Parameters for weighting climatology and persistence track forecast auc = 1.0 avc = 1.0 aup = 1.0 avp = 1.0 c efup = 1.0/48.0 efvp = 1.0/40.0 c c Parameters for LGEM intensity forecast efcp = 1.0/12.0 acc = 1.0 acp = 1.0 rnn = 2.5 beta = 1.0/24.0 c return end subroutine uvcal(slon,slat,vmax,t,imon,iday,iftype, + cxp,cyp,cxc,cyc, + wxp,wyp,wxc,wyc, + cx,cy) c c This routine estimates the storm motion vector (cx,cy) c as a weighted average of climatology and persistence c c Input: c slon: storm longitude ( 0 to 360) c slat: storm latitude (-90 to 90) c vmax: max wind (kt) c t : forecast time (hr) c imon: Month (1-12) c iday: Day of the month (1-31) c iftype: Method for averaging persistence and climatology c (=0 for persistence, =1 for climo, =2 for time weighted average) c cxp,cyp: t=0 motion vector (m/s) c c Output: c wxp: The weight on the x-component of persistence c wyp: The weight on the y-component of persistence c wxc: The weight on the x-component of climatology c wyc: The weight on the y-component of climatology c cx: The x-component of the motion vector (m/s) c cy: The y-component of the motion vector (m/s) c common /pncon3/ aup,avp,efup,efvp,auc,avc c c Get climatological storm motion vector call cluvc(slon,slat,imon,iday,cxc,cyc,cappa) c c Calcuate the climatology and persistence weights if (iftype .eq. 0) then wxp = 1.0 wyp = 1.0 wxc = 0.0 wyc = 0.0 elseif (iftype .eq. 1) then wxp = 0.0 wyp = 0.0 wxc = 1.0 wyc = 1.0 else wxp = aup*exp(-t*efup) wyp = avp*exp(-t*efvp) wxc = auc*(1.0-wxp) wyc = avc*(1.0-wyp) endif c cx = wxc*cxc + wxp*cxp cy = wyc*cyc + wyp*cyp c return end subroutine tsprint(t,flon,flat,vmax, + cxp,cyp,cxc,cyc, + wxp,wyp,wxc,wyc, + cx,cy,dtlnd,sst,vmpi, + cappap,cappac,cappa, + wcp,wcc,lulog) c c This routine prints the storm info at a single time c tlon = flon if (tlon .gt. 180.0) tlon=tlon-360.0 c write(lulog,100) t,tlon,flat,vmax,cxp,cyp,cxc,cyc, + wxp,wyp,wxc,wyc,cx,cy,dtlnd,sst,vmpi, + 100.0*cappap,100.0*cappac,100.0*cappa, + wcp,wcc 100 format('t=',f7.1,' lon,lat,v=',f7.2,1x,f6.2,1x,f6.1, + ' cpx,y=',f6.2,1x,f6.2,' ccx,y=',f6.2,1x,f6.2, + ' wpx,y=',f5.2,1x,f5.2,' wcx,y=',f5.2,1x,f5.2, + ' cx,y=',f6.2,1x,f6.2,' dtl=',f6.0,' sst=',f5.1, + ' mpi=',f4.0,' kp,kc,k=',3(f5.2,1x), + ' wcp,wcc=',f5.2,1x,f5.2) c return end subroutine dparm(slat,redfac,alpha,vb) c This routine calculates the parameters of the inland decay c model as a function of latitude c common /pncon2/ rf1,rf2,a1,a2,vb1,vb2,rclat1,rclat2,rcrad c if (slat .ge. rclat2) then redfac = rf2 alpha = a2 vb = vb2 elseif (slat .le. rclat1) then redfac = rf1 alpha = a1 vb = vb1 else w1 = (rclat2-slat)/(rclat2-rclat1) w2 = (slat-rclat1)/(rclat2-rclat1) c redfac = w1*rf1 + w2*rf2 alpha = w1*a1 + w2*a2 vb = w1*vb1 + w2*vb2 endif c return end subroutine getclimo(ioper,ierrc) c This routine gets the monthly climatology fields for u,v,cappa and SST c parameter (mxc=360,myc=181,ntc=12) c dimension uclim(mxc,myc,ntc),vclim(mxc,myc,ntc),cclim(mxc,myc,ntc) dimension sstclim(mxc,myc,ntc) dimension ndmo(ntc) c character *256 coef_location,fntemp c common /climo/ uclim,vclim,cclim,glon1,glat1,dlon,dlat,nxc,nyc common /climoo/ sstclim,glono1,glato1,dlono,dlato,nxco,nyco common /climot/ ndmo c dimension txy(nxc,nyc) c character *3 tchar c c Specify the unit number to use lut = 61 c c Scale for cappa cscale = 100.0 c c Specify the number of days in each month ndmo( 1) = 31 ndmo( 2) = 28 ndmo( 3) = 31 ndmo( 4) = 30 ndmo( 5) = 31 ndmo( 6) = 30 ndmo( 7) = 31 ndmo( 8) = 31 ndmo( 9) = 30 ndmo(10) = 31 ndmo(11) = 30 ndmo(12) = 31 c c *** Open and read the u,v,cappa climo file if (ioper .eq. 1) then call getenv("SHIPS_COEF",coef_location ) fntemp = trim( coef_location )//'oban_WH_1982_2018.dat' else fntemp = 'oban_WH_1982_2018.dat' endif c write(6,888) fntemp 888 format(a256) c open(file=fntemp,unit=lut,form='formatted', + status='old',err=901) c c ++ Read the header line and check parameters read(lut,*,err=911) glon1,dlon,nxc,glat1,dlat,nyc c c ++ Read the data do m=1,ntc read(lut,*,err=911) read(lut,*,err=911) do j=nyc,1,-1 read(lut,*,err=911) (uclim(i,j,m),i=1,nxc) enddo c read(lut,*,err=911) read(lut,*,err=911) do j=nyc,1,-1 read(lut,*,err=911) (vclim(i,j,m),i=1,nxc) enddo c read(lut,*,err=911) read(lut,*,err=911) do j=nyc,1,-1 read(lut,*,err=911) (cclim(i,j,m),i=1,nxc) enddo enddo c c Scale cappa do m=1,12 do j=1,nyc do i=1,nxc cclim(i,j,m) = cclim(i,j,m)/cscale enddo enddo enddo c close(lut) c c *** Open and read the SST climo file if (ioper .eq. 1) then call getenv("SHIPS_COEF",coef_location ) fntemp = trim( coef_location )//'clim_rsst_1982_2018.dat' else fntemp = 'clim_rsst_1982_2018.dat' endif c open(file=fntemp,unit=lut,form='formatted', + status='old',err=903) c c ++ Read the header line and check parameters read(lut,*,err=913) iyr1,iyr2,nxco,glono1,dlono,nyco,glato1,dlato c c ++ Read the data do m=1,ntc read(lut,100,err=912) tchar 100 format(a3) read(lut,101,err=912) ptemp 101 format(5x,f6.1) do j=1,nyco read(lut,102,err=912) (sstclim(i,j,m),i=1,nxco) 102 format(1000f8.2) enddo read(lut,100,err=912) tchar enddo close(lut) c ierrc = 0 return c 901 continue ierrc=1 return c 903 continue ierrc=3 return c 911 continue ierrc=11 return c 912 continue ierrc=12 return c 913 continue ierrc=13 return c end subroutine clsst(slon,slat,imon,iday,sst) c This routine interpolates the monthly SST climatology c values to the point (slon,slat) on imon/iday. c common /pncon1/ pi,dtr,ckt2ms,deg2m,eradkm,eradm,rmiss common /climoo/ sstclim,glono1,glato1,dlono,dlato,nxco,nyco common /climot/ ndmo c parameter (mxc=360,myc=181,ntc=12) c dimension sstclim(mxc,myc,ntc) dimension ndmo(ntc) c c Climo only includes 28 days in Feb so adjust Feb 29th idayt = iday if (imon .eq. 2 .and. iday .eq. 29) idayt=28 c c Check for illegal dates if (imon .lt. 1 .or. imon .gt. 12) then sst = rmiss return endif if (idayt .lt. 1 .or. idayt .gt. ndmo(imon)) then sst = rmiss return endif c c Find time weights call tweight(imon,idayt,ndmo,ntc,imon1,imon2,wt1,wt2) c c Find lon,lat array indices of the four points surrounding c the point slon,slat: (i0,j0),(i1,j0),(i0,j1),(i1,j1) and the c weights for linearly interpolating a function to the input point c (w00,w10,w01,w11). call llindx(slon,slat,glono1,glato1,dlono,dlato,nxco,nyco, + i0,i1,j0,j1,w00,w10,w01,w11) c sst1 = w00*sstclim(i0,j0,imon1) + w10*sstclim(i1,j0,imon1) + + w01*sstclim(i0,j1,imon1) + w11*sstclim(i1,j1,imon1) c sst2 = w00*sstclim(i0,j0,imon2) + w10*sstclim(i1,j0,imon2) + + w01*sstclim(i0,j1,imon2) + w11*sstclim(i1,j1,imon2) c sst = wt1*sst1 + wt2*sst2 c return c end subroutine tweight(imon,iday,ndmo,ntc,imon1,imon2,wt1,wt2) c This routine calculates the time weights (wt1,wt2) c for linearly interpolating to a date located between c any two months. The indices imon1 and imon2 of the two c months are also returned. c dimension ndmo(ntc) c c Find indices of 2 months closest to input date imon1 = imon rndmo1 = float(ndmo(imon1)) rmid1 = 0.5*(rndmo1+1.0) rday = float(iday) if (rday .ge. rmid1) then imon2 = imon1 + 1 else imon2 = imon1 - 1 endif c c Impose temporal periodicity if (imon2 .eq. 13) imon2 = 1 if (imon2 .eq. 0) imon2 = 12 c rndmo2 = float(ndmo(imon2)) rmid2 = 0.5*(rndmo2+1.0) c c Calculate time weights dcen1 = abs(rmid1-rday) dedg1 = rmid1-dcen1-1.0 dcen2 = rmid2 d1 = dcen1 d2 = dcen2+dedg1 wt1 = d2/(d1+d2) wt2 = d1/(d1+d2) c c write(6,800) imon1,imon2,rndmo1,rndmo2,d1,d2,wt1,wt2 c 800 format('imon1,imon2,rndmo1,rndmo2: ',2(i4),2(f6.0), c + ' d1,d2: ',f6.1,1x,f6.1,' wt1,wt2: '3(f6.2)) c return end subroutine llindx(slon,slat,glon1,glat1,dlon,dlat,nxc,nyc, + i0,i1,j0,j1,w00,w10,w01,w11) c This routine finds the lon,lat indices (i0,j0) of the point to the c lower left of the input point slon,slat. An evenly spaced lat/lon grid c that is periodic in x is assumed, where the point (1,1) corresponds c to (glon1,glat1) and there are (nx,ny) total grid points. For the x c periodicity it is assumed that f(1,j) = f(nxc+1). c c The incides for defining all four points around the input point (i1,j1) c are also calculated as well as the linear interpolation weights for those c 4 points (w00,w10,w01,w11) c c Find the indices of the point to the lower left of the input point i0 = 1 + ifix( (slon-glon1)/dlon ) j0 = 1 + ifix( (slat-glat1)/dlat ) c c Do not let index be outside the domain if (i0 .gt. nxc) i0 = nxc if (i0 .lt. 1) i0 = 1 if (j0 .gt. nyc) j0 = nyc if (j0 .lt. 1) j0 = 1 c c Find indices of the point to the upper right of the input point i1 = i0 + 1 j1 = j0 + 1 c c Do not let index be outside the domain if (i1 .gt. nxc) i1 = nxc if (i1 .lt. 1) i1 = 1 if (j1 .gt. nyc) j1 = nyc if (j1 .lt. 1) j1 = 1 c c Calculate normalized x,y coordinates in the grid cell x = ((slon-glon1) - float(i0-1)*dlon)/dlon y = ((slat-glat1) - float(j0-1)*dlat)/dlat c c Calculate weights w00 = 1 - x - y + x*y w10 = x - x*y w01 = y - x*y w11 = x*y c c write(6,800) slon,slat,i0,i1,j0,j1,x,y,w00,w10,w01,w11 c 800 format(/,'lon,lat,i0,i1,j0,j1: ',2(f6.1,1x),4(i4), c + ' x,y: ',f6.2,1x,f6.2,' w00,10,01,11: ',4(f7.3,1x)) c return end subroutine ampi(cx,cy,vmpi) c This routine adjusts the MPI by a adding a fraction of the c translational speed. c c Input: cx,cy - the components of the motion vector (m/s) c vmpi - the unadjusted MPI in kt c c Output: vmpi - the adjusted MPI in kt c common /pncon1/ pi,dtr,ckt2ms,deg2m,eradkm,eradm,rmiss c spd = sqrt(cx*cx + cy*cy)/ckt2ms c if (spd .le. 0.0) return vmpi = vmpi + 1.5*(spd**0.63) c return end subroutine cluvc(slon,slat,imon,iday,u,v,cappa) c This routine gets the climatological values of u,v (m/s) c and cappa (1/hr) at the point slon,slat on the date imon/iday c parameter (mxc=360,myc=181,ntc=12) c dimension uclim(mxc,myc,ntc),vclim(mxc,myc,ntc),cclim(mxc,myc,ntc) dimension ndmo(ntc) c common /climo/ uclim,vclim,cclim,glon1,glat1,dlon,dlat,nxc,nyc common /climot/ ndmo c common /pncon1/ pi,dtr,ckt2ms,deg2m,eradkm,eradm,rmiss c c Climo only includes 28 days in Feb so adjust Feb 29th idayt = iday if (imon .eq. 2 .and. idayt .eq. 29) idayt=28 c c Check for illegal dates if (imon .lt. 1 .or. imon .gt. 12) return if (idayt .lt. 1 .or. idayt .gt. ndmo(imon)) return c c Find time weights call tweight(imon,idayt,ndmo,ntc,imon1,imon2,wt1,wt2) c c Find lon,lat array indices of the four points surrounding c the point slon,slat: (i0,j0),(i1,j0),(i0,j1),(i1,j1) and the c weights for linearly interpolating a function to the input point c (w00,w10,w01,w11). call llindx(slon,slat,glon1,glat1,dlon,dlat,nxc,nyc, + i0,i1,j0,j1,w00,w10,w01,w11) c u1 = w00*uclim(i0,j0,imon1) + w10*uclim(i1,j0,imon1) + + w01*uclim(i0,j1,imon1) + w11*uclim(i1,j1,imon1) c u2 = w00*uclim(i0,j0,imon2) + w10*uclim(i1,j0,imon2) + + w01*uclim(i0,j1,imon2) + w11*uclim(i1,j1,imon2) c u = wt1*u1 + wt2*u2 c v1 = w00*vclim(i0,j0,imon1) + w10*vclim(i1,j0,imon1) + + w01*vclim(i0,j1,imon1) + w11*vclim(i1,j1,imon1) c v2 = w00*vclim(i0,j0,imon2) + w10*vclim(i1,j0,imon2) + + w01*vclim(i0,j1,imon2) + w11*vclim(i1,j1,imon2) c v = wt1*v1 + wt2*v2 c cappa1 = w00*cclim(i0,j0,imon1) + w10*cclim(i1,j0,imon1) + + w01*cclim(i0,j1,imon1) + w11*cclim(i1,j1,imon1) c cappa2 = w00*cclim(i0,j0,imon2) + w10*cclim(i1,j0,imon2) + + w01*cclim(i0,j1,imon2) + w11*cclim(i1,j1,imon2) c cappa = wt1*cappa1 + wt2*cappa2 c return c end subroutine capinit(rlon00,rlat00,vmx00,rlonm12,rlatm12,vmxm12, + vmpi,iyr,imon,iday,cappa00) c This routine calculates the initial growth rate (cappa00) c using 12 hr persistence. c common /pncon4/ efcp,acc,acp,rnn,beta c c Check to see if the storm is over land now or 12 hr ago. If so, c skip calculation and use climtological growth rate instead. c call gdland_table(rlon00 ,rlat00 ,dtlnd00) call gdland_table(rlonm12,rlatm12,dtlndm12) c if (dtlnd00 .lt. 0.0 .or. dtlndm12 .lt. 0.0) then cappa00 = 0.01 else cappa00 = beta*(vmx00/vmpi)**rnn + + (1.0/vmx00)*(vmx00-vmxm12)/12.0 endif c return end subroutine capcal(slon,slat,t,imon,iday, + iftype,cappap,cappac,wcp,wcc, + cappa) c This routine calculates the growth rate cappa as a weighted c average between persistence and climatology c c Input: c slon: storm longitude ( 0 to 360) c slat: storm latitude (-90 to 90) c t : forecast time (hr) c imon: Month (1-12) c iday: Day of the month (1-31) c iftype: Method for averaging persistence and climatology c (=0 for persistence, =1 for climo, =2 for time weighted average) c cappap: t=0 growth rate (1/hr) c c Output: c cappac: Climatological growth rate c wcp: Weight on persistence growth rate c wcc: Weight on climo growth ratec c cappa: Combined persistence and climo growth rate c common /pncon4/ efcp,acc,acp,rnn,beta c c Calculate climatological growth rate call cluvc(slon,slat,imon,iday,cxc,cyc,cappac) c c Calcuate the climatology and persistence weights if (iftype .eq. 0) then wcp = 1.0 wcc = 0.0 elseif (iftype .eq. 1) then wcp = 0.0 wcc = 1.0 else wcp = acp*exp(-t*efcp) wcc = acc*(1.0-wcp) endif c cappa = wcp*cappap + wcc*cappac c return end subroutine rmpi(sstmin,vmpimin,sstzer,vmpizer,sst,vmpi) c This routine reduces linearly reduces mpi to vmpizer by c sst=sstzer if sst < sstmin. c if (sst .ge. sstmin) then return elseif (sst .le. sstzer) then vmpi = vmpizer else vmpi = vmpizer + (sst-sstzer)*(vmpimin-vmpizer)/(sstmin-sstzer) endif c return end subroutine tclip_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 -3: TCLIPER not available in basin c -4: error in tclip() 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 opening tcliper.com input file') stop elseif (ierrtype .eq. 2) then write(luerr,970) ierr 970 format(/,'Error reading tcliper.com input file, istat=',i4) stop elseif (ierrtype .eq. 3) then write(luerr,960) 960 format(/,'TCLIPER not available for this basin') stop elseif (ierrtype .eq. 4) then write(luerr,980) ierr 980 format(/,'Error in tclip routine, ierr=',i2) stop else write(luerr,999) 999 format(/,'Unrecognized error in TCLIPER') stop endif c return end subroutine tclip_err_handling