c**************************************************************** c program stormking.yr.f c C TO COMPILE: C f77 stormking.yr.f -o stormking.yr -L/usr/local/lib -lnetopen -luduerr -lnetcdf -ludunits C c This program identifies/tracks cyclones using as a basis c gridded sea level pressure (SLP) analyses. Each cyclone c identified is ascribed a unique number which is maintained c thruought the life history of the system from cyclogenesis c to cyclolysis. The output variables (written to unit 5) c include the position of each cyclone, the cyclone number, the c year, month day and hour of each observation, central pressure, c the local Laplacian of SLP at the cyclone center (a measure of c intensity), pressure tendency (as determined between subsequent c central pressure values) and whether the system represents a c cyclogenesis or cyclolysis event, based on the first and last c observation. The output variables are listed near the bottom c of the program just before the the write/format statements c corresponding to unit 5. All cyclone numbers are resest c at 1 January 0000Z of each year. c The program was originally developed by M. Serreze (Univ. c Colorado, Boulder) for application to 12-hourly NMC fields c in the Octagonal Grid format for the entire Northern Hemishere. c The present version has been modified by F. Lo and M. Serreze c for application to 6-hourly Northern Hemisphere fields from c the NCEP/NCAR reanalysis. With only modest changes, it c can also be applied to SLP fields in other formats (e.g., c ECMWF, GCM or regional model outputs). c The logic of the algorithm is discussed by Serreze [1995] c and Serreze et al. [1997]. Potential users should c contact M. Serreze (serreze@kryos.colorado.edu) for further c details. The two major components of the algorithm are: 1) c detection of cyclones from a series of search patterns, c testing whether a grid-point SLP value is surrounded by grid c point values higher than the central point being tested. The c detection threshold determines how many systems are tracked c (e.g., a threshold of 1 mb will identify more systems than c if a 2 mb threshold is used); 2) system tracking, based c on a "nearest neighbor" analysis of the positions of systems c between time steps with a maximum distance threshold between c candidate pairings, with further checks based on distance c moved in the N/S and W directions and pressure tendency. c These limits are also adjustable. c c The major modifications from the original version are as c follows: c c 1) Prior to identification of cyclone centers, the input c SLP arrays are interpolated to a 250x250 km version of the NSIDC c EASE-grid [Armstrong and Brodzik, 1995]. This is a lower- c resolution form of the same equal-area projection being used at c NSIDC (Boulder, CO) for regridding of passive microwave satellite c data. The interpolation (based on Cressman weights) is necessary for c compatibility with the search logic for identifying system centers c and also promotes flexibility when applying the algorithm to SLP c fields other than the NCEP/NCAR reanalysis. Of course, the c interpoation has the undesreable effect of smoothing the fields. c The user has the option of adjusting the the search radius/number c of points in the interpolation. Future plans are to incorporate c a better interpolation scheme. c c 2) Alteration of the distance and pressure tendency thresholds c for cyclone tracking for use with 6-hourly as opposed to c 12-hourly analyses. The parameter sets as supplied work well c with the NCEP/NCAR analyses, but will need to be adjusted for use c with other SLP fields. Before attempting to do so, contact M. c Serreze. The most important threshold is set by variable c "maxdist". For the NCEP/NCAR data used here, it is set to 800 km, c meaning that the total allowable distance a cyclone can move c between 6-hourly time steps is 800 km (133 km/hr). While seemingly c "too fast" (a speed of 100 km is about the upper limit one c could ever imagine for cyclone motion), it allows for "center jumps" c to be tracked. It is also necessary as since one only has c data at specific grid points, there are only a finite number c of possible distances a cyclone can move (the distances are c "quantized"). "Maxdist" and the other distance thresholds are c adjusted to account for this. c c References: c c Armstrong, R.L. and M.J. Brodzik, 1995: An earth-gridded SSM/I c data set for cryospheric studies and global change monitoring, c Adv. Spave. Res., 16(10), 155-163. c c Serreze, M.C., 1995: Climatological aspects of cyclone development c and decay in the Arctic, Atmos.-Ocean, 33, 1-23. c c Serreze, M.C., F. Carse, R.G. Barry and J.C. Rogers, 1997: c Icelandic Low cyclone activity: Climatological Features, linkages c with the NAO and relationships with recent changes in the c Northern Hemisphere Circulation, J. Climate, 10(3), 453-464. c c c Note: You must set variable iyrlast to the first year that c is being processed. Variable iyrlast is set just below the c opening statement for the EASE grid file. c c********************************************************************** c********************************************************************** c********************************************************************** c Define adjustable parameters c**************************************************************** c c lat/lon positions of the EASE grid to be interpolated to. c parameter(ieqlon=72,jeqlat=72) c c search radius (km)for interpolations of grids. c parameter (guess = 300) c c minimum required number of points used in the interpolation. c parameter (npoints =2) c c pressure threshold (mb) for identifying cyclones. c parameter (pdiff =1) c c maximum search distance (km) for tracking storms between time steps. c parameter (maxdist=800) c c limit on allowable absolute pressure tendency (mb) between time steps. c parameter (pchange =20) c c limit (km) on allowable westward migration between time steps. c parameter (westmax =-600) c c limit (km) on allowable southward migration between time steps. c parameter (southmax =-600) c c limit (km) on allowable northward migration between time steps. c parameter (northmax =600) c c********************************************************************** c********************************************************************** c Define Variables and arrays c integer indexlat, indexlon real lat(ieqlon,jeqlat), long(ieqlon,jeqlat) real pi, r real meters(ieqlon,jeqlat,4) real radlat(ieqlon,jeqlat) real press(ieqlon,jeqlat), *p3(20000,60),zlat2(20000,60),zlong2(20000,60), *dist1(60,60),zmin1(60),zstart(12),zjul(20000), *dlat(60,60),dlat3(60),dp(60,60),tdist(20000,60), *dp2(60),zstart1(12), *psamp(10,ieqlon,jeqlat), *lap(ieqlon,jeqlat),lap1(20000,40), *var3(60),var4(60),var5(60),var6(60),dlong(60,60), *dlong2(60,60),dlat2(60,60),tojo(60,60),tojo1(60,60),dl5(60), *tendency(20000,60),zlong5(8),zlat5(8),p5(8),dlat6(2000,60), *dlong6(2000,60) c real latncept(8000),lonncept(8000), latstop(6000), t, *longstop(6000),pressa(8000),press1(8000),sum1,sum2 c integer tdiff(60),tdiff1(20000,60),iswitch(ieqlon,jeqlat), *ibad(ieqlon,jeqlat), *ncenter(ieqlon,jeqlat),status(20000,60),i2(20000,60),j2(20000,60), *icount(20000),nsys(20000,60),il(60),ih(60),jl(60),jh(60), *nsys1(60,60),nsys3(60),idog(60),idup(60,60),jdup(60,60), *var1(60),var2(60),cyclo(20000,60),lysis(20000,60), *imon1(20000),iday1(20000),ihr1(20000),iyr1(20000),iskip(20000), *npair(8),npair1(8),imonlast,idaylast real sum((ieqlon*jeqlat),500) real isearch((ieqlon*jeqlat), 500),isearch1(ieqlon*jeqlat) integer nchartstart,nchartstop,nchartint character month*3, fnameout*29, filenetcdf*42 c c*************************************************************** c*************************************************************** c**** READ IN RAW DATA IN netCDF FORMAT (4x daily) ************* c*************************************************************** integer fileid, varid, latid, lonid, timid, ilon, jlat, ntime integer vrdim, vrat, vrtyp, lttyp, ltdim, ltat integer tmtyp, tmdim, tmat, lntyp, lndim, lnat real*4 sf, ao character*31 vrnam, tmnam, lnnam, ltnam integer vshp(4), tmshp(4) parameter(ilon=144,jlat=36, ntime=365*4, ntimeleap=366*4) real*8 timencep(ntime) real*8 timencepleap(ntimeleap) real*4 latncep(jlat), lonncep(ilon) integer*2 islp (ilon, jlat, ntime) integer*2 islpleap (ilon, jlat, ntimeleap) integer count(3), start(3) include'/usr/local/include/netcdf.inc' data zstart/0,31,59,90,120,151,181,212,243,273,304,334/ data zstart1/0,31,60,91,121,152,182,213,244,274,305,335/ c**************************************************************** print*, ' What year? [4 digits, e.g. 2001]' read(*,*) iyear print*, ' Up to what month? [3 letters, e.g. dec]' read(*,*) month c**************************************************************** write(filenetcdf,"('sfc_data/slp.',i4.4,'.nc')") iyear ctjp & "('/Datasets/ncep.reanalysis/surface/slp.',i2,'.nc')") iyear c fileid = NCOPN('/Datasets/ncep.reanalysis/surface/slp.97.nc', 0, c & rcode) print*, 'filenetcdf = ',filenetcdf fileid = NCOPN(filenetcdf,0,rcode) varid = NCVID(fileid, 'slp', rcode) latid= NCVID(fileid, 'lat', rcode) lonid= NCVID(fileid, 'lon', rcode) timid= NCVID(fileid, 'time', rcode) slp=0 call ncvinq(fileid,varid,vrnam,vrtyp,vrdim,vshp,vrat,rcode) call ncvinq(fileid,timid,tmnam,tmtyp, tmdim, tmshp, tmat, rcode) call ncvinq(fileid,latid,ltnam,lttyp, ltdim, ltshp, ltat, rcode) call ncvinq(fileid,lonid,lnnam,lntyp, lndim, lnshp, lnat, rcode) call ncvgt(fileid, lonid, 1, ilon, lonncep, rcode) call ncvgt(fileid, latid, 1, jlat, latncep, rcode) if((float(iyear)/4).eq.(AINT(float(iyear)/4))) then c if its a leap year: call ncvgt(fileid, timid, 1, ntimeleap, timencepleap, rcode) else call ncvgt(fileid, timid, 1, ntime, timencep, rcode) endif call ncagt(fileid, varid, 'scale_factor', sf, rcode) call ncagt(fileid, varid, 'add_offset', ao, rcode) start(1)=1 start(2)=1 start(3)=1 count(1)=ilon count(2)=jlat count(3)=ntime if((float(iyear)/4).eq.(AINT(float(iyear)/4))) then c if its a leap year: count(3) = ntimeleap call ncvgt(fileid, varid, start, count, islpleap, rcode) else call ncvgt(fileid, varid, start, count, islp, rcode) endif c *************************************************************** c open input and output units c**************************************************************** c c Output file with cyclone data c c write(fnameout,"('output/',i4.4,'.out')") iyear write(fnameout,"('output/ncepstorms.',i4.4)") iyear open(unit=5,file=fnameout,status='unknown') c**************************************************************** c**************************************************************** c Test outputs c c open(unit=30,name='test1.d',status='unknown') c open(unit=31,name='test2.d',status='unknown') c open(unit=32,file='test60.d',status='unknown') c open(unit=33,file='bad.d',status='unknown') c open(unit=14,file='nmcgridpoints', status='unknown') c open(unit=15,file='pressncept',status='unknown') c open(unit=16,file='interpress',status='unknown') c c*************************************************************** c EASE grid file (the one you want interpolate to) c open(177,file='/home/map/storm/eqarea250.grid',form ='formatted') c*************************************************************** c**************************************************************** c$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ c Set starting year c iyrlast=iyear c c$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ c Read in EASE grid c do j=1,jeqlat do i=1,ieqlon lat(i,j)=999 long(i,j)=999 enddo enddo do k=1,4060 read(177,*) indexlat, indexlon, & lat((indexlon+1),(indexlat+1)), & long((indexlon+1),(indexlat+1)) enddo close(177) do j=1,jeqlat do i=1,ieqlon if (lat(i,j).eq.90) lat(i,j)=89.9 if (long(i,j).lt.0) long(i,j)=360+long(i,j) enddo enddo c c transform two-dimensional array to linear array c do i=1,ieqlon do j=1,jeqlat q=q+1 latstop(q)=lat(i,j) longstop(q)=long(i,j) end do end do cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c equal area grid size is 250km. c figure out length and width of the grids. c pi = 3.1415926 c radius of the earth in km r = 6380 c change lat to radians do j=1,jeqlat do i=1,ieqlon if(lat(i,j).eq.999.)then radlat(i,j) = 999. else radlat(i,j) = lat(i,j)*pi/180. endif enddo enddo do j=2,jeqlat-1 do i=2,ieqlon-1 if ((radlat(i,j).eq.999.).or.(long(i,j).eq.999.).or. & (radlat(i+1,j).eq.999.).or.(long(i+1,j).eq.999.).or. & (radlat(i-1,j).eq.999.).or.(long(i-1,j).eq.999.).or. & (radlat(i,j+1).eq.999.).or.(long(i,j+1).eq.999.).or. & (radlat(i,j-1).eq.999.).or.(long(i,j-1).eq.999.)) then c do k =1,4 c meters(i,j,k)=-999. c enddo goto 111 else c d1 is the lat/lon distance between (i+1,j) and (i,j) d1lon = diflong(long(i+1,j),long(i,j)) & *pi*r*cos(radlat(i,j))/180. d1lat = abs((radlat(i+1,j)-radlat(i,j))*r) c d2 is the lat/lon distance between (i-1,j) and (i,j) d2lon = diflong(long(i,j),long(i-1,j)) & *pi*r*cos(radlat(i,j))/180. d2lat = abs((radlat(i,j)-radlat(i-1,j))*r) c d3 is the lat/lon distance between (i,j-1) and (i,j) d3lon = diflong(long(i,j-1),long(i,j)) & *pi*r*cos(radlat(i,j))/180. d3lat = abs((radlat(i,j-1)-radlat(i,j))*r) c d4 is the lat/lon distance between (i,j+1) and (i,j) d4lon = diflong(long(i,j),long(i,j+1)) & *pi*r*cos(radlat(i,j))/180. d4lat = abs((radlat(i,j)-radlat(i,j+1))*r) c meters(i,j, ) is distance between the points meters(i,j,1) = sqrt( d1lat**2 + d1lon**2) meters(i,j,2) = sqrt( d2lat**2 + d2lon**2) meters(i,j,3) = sqrt( d3lat**2 + d3lon**2) meters(i,j,4) = sqrt( d4lat**2 + d4lon**2) write (98,'(2i3,4f10.3)') i,j,(meters(i,j,k),k=1,4) endif 111 enddo enddo cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c c**************************************************************** c read in NCEP reanalysis pressure array, initialize variables c**************************************************************** c******************************** c THE BIG LOOP: TO LINE 100 * c******************************** c if(month.eq.'jan') then nchartstart = 1 nchartstop = (31*4) +1 imon=1 endif c if its a leap year: if((float(iyear)/4).eq.(AINT(float(iyear)/4))) then if(month.eq.'feb') then nchartstart = 31*4 +1 nchartstop = 60*4 +1 imon=2 elseif(month.eq.'mar') then nchartstart = 60*4 +1 nchartstop = 91*4 +1 imon=3 elseif(month.eq.'apr') then nchartstart = 91*4 +1 nchartstop = 121*4 +1 imon=4 elseif(month.eq.'may') then nchartstart = 121*4 +1 nchartstop = 152*4 +1 imon=5 elseif(month.eq.'jun') then nchartstart = 152*4 +1 nchartstop = 182*4 +1 imon=6 elseif(month.eq.'jul') then nchartstart = 182*4 +1 nchartstop = 213*4 +1 imon=7 elseif(month.eq.'aug') then nchartstart = 213*4 +1 nchartstop = 244*4 +1 imon=8 elseif(month.eq.'sep') then nchartstart = 244*4 +1 nchartstop = 274*4 +1 imon=9 elseif(month.eq.'oct') then nchartstart = 274*4 +1 nchartstop = 305*4 +1 imon=10 elseif(month.eq.'nov') then nchartstart = 305*4 +1 nchartstop = 335*4 +1 imon=11 elseif(month.eq.'dec') then nchartstart = 335*4 +1 nchartstop = 366*4 +1 imon=12 endif else if(month.eq.'feb') then nchartstart = 31*4 +1 nchartstop = 59*4 +1 imon=2 elseif(month.eq.'mar') then nchartstart = 59*4 +1 nchartstop = 90*4 +1 imon=3 elseif(month.eq.'apr') then nchartstart = 90*4 +1 nchartstop = 120*4 +1 imon=4 elseif(month.eq.'may') then nchartstart = 120*4 +1 nchartstop = 151*4 +1 imon=5 elseif(month.eq.'jun') then nchartstart = 151*4 +1 nchartstop = 181*4 +1 imon=6 elseif(month.eq.'jul') then nchartstart = 181*4 +1 nchartstop = 212*4 +1 imon=7 elseif(month.eq.'aug') then nchartstart = 212*4 +1 nchartstop = 243*4 +1 imon=8 elseif(month.eq.'sep') then nchartstart = 243*4 +1 nchartstop = 273*4 +1 imon=9 elseif(month.eq.'oct') then nchartstart = 273*4 +1 nchartstop = 304*4 +1 imon=10 elseif(month.eq.'nov') then nchartstart = 304*4 +1 nchartstop = 334*4 +1 imon=11 elseif(month.eq.'dec') then nchartstart = 334*4 +1 nchartstop = 365*4 +1 imon=12 endif endif nchartstart = 1 nchartstop = 365 * 4 + 1 if(iyear .eq. (iyear/4)*4) nchartstop = 366 * 4 + 1 print*, nchartstart,nchartstop,nchartint nchartint =1 do 100 ncharts=nchartstart,nchartstop,nchartint t=float(ncharts) iyr=iyear ihr = 6 * mod(ncharts-1,4) jday = ((ncharts+3) - (ihr/6))/4 if(iyr .eq. (iyr/4)*4) then do jmon = 1, 12 if(jday .gt. zstart1(jmon)) imon = jmon end do iday = jday - zstart1(imon) else do jmon = 1, 12 if(jday .gt. zstart(jmon)) imon = jmon end do iday = jday - zstart(imon) endif c solve for imon, iday, ihr from the timestep c if((float(iyr)/4).eq.(AINT(float(iyr)/4))) then c if iyr is a leap year c if (imon .eq. 1) iday=AINT(t/4+0.75) c if (imon .eq. 2) iday=AINT(t/4+0.75)-31 c if (imon .eq. 3) iday=AINT(t/4+0.75)-31-29 c if (imon .eq. 4) iday=AINT(t/4+0.75)-31-29-31 c if (imon .eq. 5) iday=AINT(t/4+0.75)-31-29-31-30 c if (imon .eq. 6) iday=AINT(t/4+0.75)-31-29-31-30-31 c if (imon .eq. 7) iday=AINT(t/4+0.75)-31-29-31-30-31-30 c if (imon .eq. 8) iday=AINT(t/4+0.75)-31-29-31-30-31-30-31 c if (imon .eq. 9)iday=AINT(t/4+0.75)-31-29-31-30-31-30-31-31 c if (imon .eq. 10) iday=AINT(t/4+0.75)-244-30 c if (imon .eq. 11) iday=AINT(t/4+0.75)-244-30-31 c if (imon .eq. 12) iday=AINT(t/4+0.75)-244-30-31-30 c else c if (imon .eq. 1) iday=AINT(t/4+0.75) c if (imon .eq. 2) iday=AINT(t/4+0.75)-31 c if (imon .eq. 3) iday=AINT(t/4+0.75)-31-28 c if (imon .eq. 4) iday=AINT(t/4+0.75)-31-28-31 c if (imon .eq. 5) iday=AINT(t/4+0.75)-31-28-31-30 c if (imon .eq. 6) iday=AINT(t/4+0.75)-31-28-31-30-31 c if (imon .eq. 7) iday=AINT(t/4+0.75)-31-28-31-30-31-30 c if (imon .eq. 8) iday=AINT(t/4+0.75)-31-28-31-30-31-30-31 c if (imon .eq. 9)iday=AINT(t/4+0.75)-31-28-31-30-31-30-31-31 c if (imon .eq. 10) iday=AINT(t/4+0.75)-243-30 c if (imon .eq. 11) iday=AINT(t/4+0.75)-243-30-31 c if (imon .eq. 12) iday=AINT(t/4+0.75)-243-30-31-30 c endif c if (((t/4)-AINT(t/4)) .eq. 0.25) then c ihr=0 c elseif (((t/4)-AINT(t/4)) .eq. 0.5) then c ihr=6 c elseif (((t/4)-AINT(t/4)) .eq. 0.75) then c ihr=12 c elseif (((t/4)-AINT(t/4)) .eq. 0) then c ihr=18 c end if PRINT *,'CHART READ IN, NCHARTS =', ncharts PRINT *, 'yr, mon, dy, hr, t', iyr, imon, iday, ihr, t nstorms=0 c c Initialize variables c do k=1,60 tdiff(k)=0 do i=1,ieqlon do j=1,jeqlat iswitch(i,j)=0 ibad(i,j)=0 ncenter(i,j)=0 end do end do end do q=0 k=0 c********************************************************* c********************************************************** c interpolation: reanalysis grid to equal area grid c********************************************************* c*************************************************************** c read in NCEP/NCAR Grid (the original grid) and pressure c data maz=0 do j=1,144 do l=1,36 maz=maz+1 latncept(maz)=latncep(l) lonncept(maz)=lonncep(j) if((float(iyear)/4).eq.(AINT(float(iyear)/4))) then tdash = t if(tdash .gt. real(ntimeleap)) tdash = real(ntimeleap) pressa(maz)=0.01*(ao+islpleap(j,l,tdash)) else tdash = t if(tdash .gt. real(ntime)) tdash = real(ntime) pressa(maz)=0.01*(ao+islp(j,l,tdash)) endif end do end do c open(unit=17,file='press.dat',FORM='UNFORMATTED', c * ACCESS='DIRECT', RECL=144*73*10) c write(17,REC=1) ((islp(i,j,1),i=1,144),j=1,36) c c************************************************************ c reinterpolate reanalysis grid (2.5x2.5) to the equal area c grid using Cressman weights. Interpolate full grid for the first c only and then store distances. c if (ncharts.eq.nchartstart) then do 57 q=1,(ieqlon*jeqlat) ilevs=0 isearch1(q) =0 if ((latstop(q).eq.999).or.(longstop(q).eq.999)) then press1(q)= -999. goto 57 endif do 56 i=1,(144*36) c c************************************************************ c find distance between equal area grid point to be interpolated c to and original point to be interpolated from. dim=(sind(latstop(q))*sind(latncept(i))) dam=(cosd(latstop(q))*cosd(latncept(i))) dam1=longstop(q)-lonncept(i) dam1=cosd(dam1) dist=dim+dam*dam1 dist=acosd(dist)*110.949 c************************************************************* c if distance is .lt. guess then identify those c cases and find the Cressman weights c if(dist .le. guess) then zpiggy=(guess**2)-(dist**2) zpiggy1=(guess**2)+(dist**2) ilevs=ilevs+1 sum(q,ilevs)=zpiggy/zpiggy1 isearch(q,ilevs)=i isearch1(q) = ilevs endif 56 continue if (ilevs.lt.npoints) then press1(q)= -999. goto 57 endif 57 continue endif do 157 q =1,(ieqlon*jeqlat) if (press1(q).ne.-999.) then sum1=0 sum2=0 do 156 i=1,isearch1(q) sum1=sum1+(sum(q,i)*pressa(isearch(q,i))) sum2=sum2+sum(q,i) 156 continue c************************************************************* c find interpolated value c press1(q)=sum1/sum2 endif 157 continue c transform one dimensional reinterpolated c pressure array to two dimensional equal area array c q=0 do i=1,ieqlon do j=1,jeqlat q=q+1 press(i,j)=press1(q) end do end do c************************************************************* c FINISHED INTERPOLATION FOR A CHART c************************************************************* c************************************************************* c************************************************************* c************************************************************* c************************************************************* c c store sample data fields for plotting (will be written into c unit 32) c if(ncharts .gt. 50 .and. ncharts .le. 60) then impy=impy+1 do i=1,ieqlon do j=1,jeqlat psamp(impy,i,j)=press(i,j) end do end do endif c c c skip days with read errors or ridiculous data c do i=1,ieqlon do j=1,jeqlat if(press(i,j) .gt. 1200) goto 100 if(imon .gt. 12) goto 100 if(iday .gt. 32) goto 100 if(press(i,j) .eq. 0.0) then press(i,j)= -999 endif end do end do c c skip is chart out of order c if(ncharts .gt. 1) then if(imon .eq. imonlast .and. iday .lt. idaylast) then print *,'chart out of order',ncharts,iyr,imon,iday,ihr goto 100 endif endif imonlast=imon idaylast=iday c*************************************************************** c$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ c$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ c c CYCLONE DETECTION SECTION c*************************************************************** c*************************************************************** c identify cyclone centers ******************************** c*************************************************************** c*************************************************************** c do 60 i=4,(ieqlon-3) do 50 j=4,(jeqlat-3) itype=0 c c find local laplacian of pressure field c if((press(i,j).eq.-999).or.(press(i+1,j).eq.-999).or. & (press(i-1,j).eq.-999).or.(press(i,j+1).eq.-999).or. & (press(i,j-1).eq.-999).or.(meters(i,j,1).eq.0).or. & (meters(i,j,2).eq.0).or.(meters(i,j,3).eq.0).or. & (meters(i,j,4).eq.0)) then lap(i,j)=-999 else zdiff1=press(i+1,j)-press(i,j) zdiff2=press(i-1,j)-press(i,j) zdiff3=press(i,j-1)-press(i,j) zdiff4=press(i,j+1)-press(i,j) d1 = meters(i,j,1)*(meters(i,j,1)+meters(i,j,2)) d2 = meters(i,j,2)*(meters(i,j,1)+meters(i,j,2)) d3 = meters(i,j,3)*(meters(i,j,3)+meters(i,j,4)) d4 = meters(i,j,4)*(meters(i,j,3)+meters(i,j,4)) lap(i,j)=2*((zdiff1/d1)+(zdiff2/d2)+(zdiff3/d3)+(zdiff4/d4)) lap(i,j)=lap(i,j)*100000 endif c c Identify grids over 7x7 array with pressure .lt. c central value c do ii=i-3,i+3 do jj=j-3,j+3 ibad(ii,jj)=0 ncenter(ii,jj)=0 end do end do c ick=0 t=0 t=press(i,j) do ii=i-3,i+3 do jj=j-3,j+3 if(press(ii,jj)-t .lt. pdiff) then ibad(ii,jj)=1 else ibad(ii,jj)=0 endif end do end do c c Check the 8 surrounding points. If pressure rises at all c 8 points by at least the seclected threshold (pdiff), c identify as cyclone. If pressure falls, skip out of loop and c check next grid. If pressure does not fall, but c is .lt. pdiff, proceeed to additional checks do determine if c central pressure identifies a cyclone center c do ii=i-1,i+1 do jj=j-1,j+1 diff=press(ii,jj)-t if(diff .lt. 0) goto 50 if(diff .lt. 0) print *,diff, i,j,ii,jj end do end do c do ii=i-1,i+1 do jj=j-1,j+1 diff=press(ii,jj)-t if(diff .lt. pdiff) goto 3332 end do end do ncenter(i,j)=2 goto 3333 c c if cyclone not found but pressures at 8 surrounding c points are .lt. pdiff, trace path to see if low is closed. c If not, skip out and check next grid c 3332 ibad(i,j)=0 do mm=i-1,i+1 do nn=j-1,j+1 if(ibad(mm,nn) .eq. 1) then do m1=mm-1,mm+1 do n1=nn-1,nn+1 if(ibad(m1,n1) .eq. 1) then do m2=m1-1,m1+1 do n2=n1-1,n1+1 if(ibad(m2,n2) .eq. 1) then if(m2 .eq. i-3) goto 50 if(m2 .eq. i+3) goto 50 if(n2 .eq. j-3) goto 50 if(n2 .eq. j+3) goto 50 endif end do end do endif end do end do endif end do end do c c PRINT *,'CYCLONE FOUND' c c identify all grid locations within the closed c perimeter c 3333 ncenter(i,j)=1 do mm=i-1,i+1 do nn=j-1,j+1 if(ibad(mm,nn) .eq. 1) then ncenter(mm,nn)=ncenter(mm,nn)+1 do m1=mm-1,mm+1 do n1=nn-1,nn+1 if(ibad(m1,n1) .eq. 1) then ncenter(m1,n1)=ncenter(m1,n1)+1 do m2=m1-1,m1+1 do n2=n1-1,n1+1 if(ibad(m2,n2) .eq. 1) then ncenter(m2,n2)=ncenter(m2,n2)+1 endif end do end do endif end do end do endif end do end do c c sample output to check if cyclones are c properly identified c c if(ncharts .le. 10) then c write(30,6452),ncharts,i,j c do ii=i-3,i+3 c write(30,4517), (press(ii,jj), jj=j-3,j+3) c end do c do ii=i-3,i+3 c write(30,4516), (ncenter(ii,jj),jj=j-3,j+3) c end do c endif c6452 format(3i5) c4517 format(7f7.1) c4516 format(7i7) c c if enclosed perimeter contains multiple grid points c test whether central pressure (press(i,j) is has an adjacent c grid point with the same pressure, which in turn borders c in a grid point with lower pressure. If so, exit and c inspect next grid as cyclone center represented by c press(i,j) is anomalous. c do mm=i-1,i+1 do nn=j-1,j+1 if(press(mm,nn) .eq. press(i,j)) then do m1=mm-1,mm+1 do n1=nn-1,nn+1 if(press(m1,n1) .lt. press(i,j)) then c do mmm=i-3,i+3 c print 456, (press(mmm,nnn), nnn=j-3,j+3) c end do c456 format(7f7.1) goto 50 endif if(press(m1,n1) .eq. press(i,j)) then do m2=m1-1,m1+1 do n2=n1-1,n1+1 if(press(m2,n2) .lt. press(i,j)) goto 50 end do end do endif end do end do endif end do end do c c if i,j grid location already defined as cyclone c center, skip to next location c if(iswitch(i,j) .gt. 0) then goto 50 else k=k+1 c search for duplicate centers (identical adjacent c pressure values) by examining pressures c at grids adjacent to press(i,j). For each cyclone c center, set iswitch(i,j) to the cyclone number (k). c c search for duplicates over 3x3 shell c do i1=i-1,i+1 do j1=j-1,j+1 if(press(i1,j1) .eq. press(i,j)) then iswitch(i1,j1)=k c c if duplicate found, extend search (up to a 5x5 array) c do i4=i1-1,i1+1 do j4=j1-1,j1+1 if(press(i4,j4) .eq. press(i,j)) then iswitch(i4,j4)=k endif end do end do c endif end do end do endif c nstorms=k c c if(nstorms .ge. 39) then c print *,'too many storms, skip out to next chart' c goto 100 c endif c if(nstorms .ge. 70) then print *,'too many storms, skip to next',nstorms, k write(33,5610) ncharts,iyr,imon,iday,ihr,nstorms goto 100 5610 format(6i10) endif c 50 continue 60 continue c print *,ncharts,iyr,imon,iday,ihr,nstorms c c c*************************************************************** c*************************************************************** c now at end of chart ************************************* c*************************************************************** c*************************************************************** c C Retest each cyclone. If duplicate centers found c (identical pressure values), define center as that with c largest (most positive) local laplacian of pressure. c Put variables in new array for subsequent tracking of c cyclones c c if no storms found, assume chart is bogus and skip out c if(nstorms .eq. 0) then print *,'error, no storms found' goto 100 endif c do k=1,nstorms n=0 iii=0 jjj=0 do i=4,(ieqlon-3) do j=4,(jeqlat-3) if(iswitch(i,j) .eq. k) then n=n+1 idup(k,n)=i jdup(k,n)=j tdiff(k)=n endif end do end do c c if no duplicate centers c if(tdiff(k) .eq. 1) then iii=idup(k,1) jjj=jdup(k,1) var1(k)=iii var2(k)=jjj var3(k)=lat(iii,jjj) var4(k)=long(iii,jjj) var5(k)=press(iii,jjj) var6(k)=lap(iii,jjj) else c c if duplicate centers c zit=0 izit=0 jzit=0 do n=1,tdiff(k) if(lap(idup(k,n),jdup(k,n)) .ge. zit) then zit=lap(idup(k,n),jdup(k,n)) izit=idup(k,n) jzit=jdup(k,n) endif c c test output for cases with duplicate cyclone centers c c if(tdiff(k) .ge. 2) then c print *,tdiff(k),n,idup(k,n),jdup(k,n), c * lat(idup(k,n),jdup(k,n)),long(idup(k,n),jdup(k,n)), c * lap(idup(k,n),jdup(k,n)),zit c endif c end do c c if(tdiff(k) .ge. 2) then c do inky=izit-5,izit+5 c print 30671, (press(inky,inky1),inky1= c * jzit-5,jzit+5) c end do c endif c30671 format(11f7.1) c iii=izit jjj=jzit var1(k)=iii var2(k)=jjj var3(k)=lat(iii,jjj) var4(k)=long(iii,jjj) var5(k)=press(iii,jjj) var6(k)=lap(iii,jjj) c c if(tdiff(k) .ge. 2) then c print *,var1(k),var2(k),var3(k),var4(k), c * var5(k),var6(k) c endif endif end do c c m=m+1 icount(m)=nstorms do k=1,nstorms p3(m,k)=var5(k) i2(m,k)=var1(k) j2(m,k)=var2(k) zlat2(m,k)=var3(k) zlong2(m,k)=var4(k) tdiff1(m,k)=tdiff(k) lap1(m,k)=var6(k) end do c c retest to see if multi-lobed systems are c retained c c do k=1,icount(m) c imult=0 c ick=iswitch(i2(m,k),j2(m,k)) c print *,k,ick,i2(m,k),j2(m,k) c do i4=i2(m,k)-2,i2(m,k)+2 c do j4=j2(m,k)-2,j2(m,k)+2 c print *,'were here',i4,j4 c if(iswitch(i4,j4) .ne. 0 .and. c * iswitch(i4,j4) ..ne. ick) then c imult=imult+1 c print *,'were here' c endif c end do c end do c c if(imult .ge. 1) then c print *,'multi-lobed system found' c do i5=i2(m,k)-3,i2(m,k)+3 c print 566, (iswitch(i5,j5),j5=j2(m,k)-3,j2(m,k)+3) c566 format(7i2) c end do c endif c end do c c c if(m .le. 100) then c do i=4,(ieqlon-3) c do j=4,(jeqlat-3) c butt(m,i,j)=press(i,j) c end do c end do c endif c c convert dates to julian days, accounting for leap c years. Put year,month,day and hour into array c do idate=1,12 if(imon .eq. idate) then zjul(m)=zstart(idate)+iday if( iyr.eq.2000 .or. iyr.eq.1996 .or. iyr.eq.1992 * .or. iyr.eq.1988 .or. iyr.eq.1984 * .or. iyr.eq.1980 .or. iyr.eq.1976) * zjul(m)=zstart1(idate)+iday endif end do C if(ihr .eq. 6) zjul(m)=zjul(m)+0.25 if(ihr .eq. 12) zjul(m)=zjul(m)+0.5 if(ihr .eq. 18) zjul(m)=zjul(m)+0.75 iyr1(m)=iyr imon1(m)=imon iday1(m)=iday ihr1(m)=ihr c c c***************************************************************** c***************************************************************** c$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ c c CYCLONE TRACKING SECTION c c$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ c***************************************************************** c***************************************************************** c c set cyclone numbers, pressure tendency, cyclogenesis/ c cyclolysis flags for chart 1 c c do k=1,60 cyclo(m,k)=0 lysis(m,k)=0 tdist(m,k)=-999 tendency(m,k)=-999 end do iskip(m)=0 c if(m .eq. 1) then iskip(m)=1 ilast=0 do k=1,icount(m) nsys(m,k)=0 end do goto 4444 endif c c if first chart in a year, renumber all systems c starting with #1 c if(iyr .ne. iyrlast) then print *, 'new year' iskip(m)=1 do k=1,icount(m) nsys(m,k)=0 end do iyrlast=iyr ilast=0 goto 4444 endif iyrlast=iyr c c if a chart is skipped, do not attempt to c track systems from previous chart c if(zjul(m) .gt. (zjul(m-1)+0.25)) then C PRINT *,'SKIPPED CHART', m iskip(m)=1 goto 4444 endif c c set initial status of systems on previous c day to untracked (status(m-1,k1)=0) c do k1=1,icount(m-1) status(m-1,k1)=0 end do c c********************************************************* c Compare systems on present date (chart m) with c systems on previous day (chart m-1). If system c on chart m-1 falls within a 3x3 shell centered c over system on chart m, system on chart m c is tracked as continuation of the system on chart m-1 c ((status(m-1,k1)=1, nsys(m,k)) is number of chart m-1 c system it was paired with) c********************************************************** c do k=1,icount(m) do k1=1,icount(m-1) if(i2(m-1,k1) .eq. i2(m,k) .and. * j2(m-1,k1) .eq. j2(m,k)) then if(status(m-1,k1) .eq. 0) then nsys(m,k)=nsys(m-1,k1) tendency(m,k)=p3(m,k)-p3(m-1,k1) tdist(m,k)=0 status(m-1,k1)=1 endif endif end do end do c do k=1,icount(m) il(k)=i2(m,k)-1 ih(k)=i2(m,k)+1 jl(k)=j2(m,k)-1 jh(k)=j2(m,k)+1 end do c c c************************************************************* c search for cases in which multiple chart m-1 cyclones c were located in a 3x3 shell centered over a chart m c cyclone c************************************************************* c do k=1,icount(m) mdup=0 if(nsys(m,k) .eq. 0) then do k1=1,icount(m-1) if(i2(m-1,k1) .ge. il(k) .and. i2(m-1,k1) .le. ih(k) * .and. j2(m-1,k1) .ge. jl(k) .and. j2(m-1,k1) * .le. jh(k) .and. status(m-1,k1) .eq. 0) then c c mdup=mdup+1 p5(mdup)=p3(m-1,k1) zlat5(mdup)=zlat2(m-1,k1) zlong5(mdup)=zlong2(m-1,k1) npair(mdup)=nsys(m-1,k1) npair1(mdup)=k1 endif end do c c if potential pair shows the chart m-1 system moving c east, pair with that system (if more than 2 duplicates, c take the one that moved the furthest east) c c test output c c if(mdup .gt. 1) then c print *,'multiple chart m-1 cyclones in' c print *,'chart m shell ',ncharts c print *, m,zlat2(m,k),zlong2(m,k),p3(m,k) c do kk=1,mdup c print *,kk,zlat5(kk),zlong5(kk),p5(kk),npair(kk),npair1(kk) c end do c endif c if(mdup .gt. 1) then zmax=-100 npair2=0 npair3=0 c do kk=1,mdup q2=zlong2(m,k) q1=zlong5(kk) zchange=q2-q1 if((q2 .gt. 270 .and. q2 .le. 360) .and. (q1 .ge. 0 * .and. q1 .lt. 90)) zchange=q2-(q1+360) if((q2 .ge. 0 .and. q2 .lt. 90) .and. (q1 .gt. 270 * .and. q1 .le. 360)) zchange=(q2+360)-q1 c if(zchange .gt. zmax) then zmax=zchange npair2=npair(kk) npair3=npair1(kk) zpair4=zlat5(kk) zpair5=zlong5(kk) zpair6=p5(kk) endif c print *,kk,zchange,zmax end do c nsys(m,k)=npair2 status(m-1,npair3)=1 dim=(sind(zlat2(m,k))*sind(zpair4)) dam=(cosd(zlat2(m,k))*cosd(zpair4)) dam1=zlong2(m,k)-zpair5 dam1=cosd(dam1) dist=dim+dam*dam1 dist=acosd(dist)*110.949 tdist(m,k)=dist tendency(m,k)=p3(m,k)-zpair6 c c test output c c write(31,7306) zjul(m),nsys(m,k) c7306 format(f5.1,i7) c print *,m,zpair4,zpair5,zpair6,nsys(m,k),npair3 c print *,tdist(m,k),tendency(m,k) endif endif end do c c**************************************************************** c search for cases in which a the same chart m-1 cyclone c appeared in 3x3 shells over different chart m cyclones. c**************************************************************** c ndup=0 do k1=1,icount(m-1) ndup=0 if(status(m-1,k1) .eq. 0) then do k=1,icount(m) if(i2(m-1,k1) .ge. il(k) .and. i2(m-1,k1) .le. ih(k) * .and. j2(m-1,k1) .ge. jl(k) .and. j2(m-1,k1) * .le. jh(k) .and. nsys(m,k) .eq. 0) then c ndup=ndup+1 p5(ndup)=p3(m,k) zlat5(ndup)=zlat2(m,k) zlong5(ndup)=zlong2(m,k) c npair(ndup)=nsys(m,k) npair1(ndup)=k endif end do c c test output c c if(ndup .gt. 1) then c print *,'chart m-1 cyclones in' c print *,'multiple chart m shells ',ncharts c print *, zlat2(m-1,k1),zlong2(m-1,k1),p3(m-1,k1),nsys(m-1,k1) c do kk=1,ndup c print *,zlat5(kk),zlong5(kk),p5(kk),npair1(kk) c end do c endif if(ndup .gt. 1) then zmax=-100 npair2=0 npair3=0 c do kk=1,ndup q2=zlong5(kk) q1=zlong2(m-1,k1) zchange=q2-q1 if((q2 .gt. 270 .and. q2 .le. 360) .and. (q1 .ge. 0 * .and. q1 .lt. 90)) zchange=q2-(q1+360) if((q2 .ge. 0 .and. q2 .lt. 90) .and. (q1 .gt. 270 * .and. q1 .le. 360)) zchange=(q2+360)-q1 c if(zchange .gt. zmax) then zmax=zchange c npair2=npair(kk) npair3=npair1(kk) zpair4=zlat5(kk) zpair5=zlong5(kk) zpair6=p5(kk) endif c print *,kk,zchange,zmax end do c nsys(m,npair3)=nsys(m-1,k1) status(m-1,k1)=1 dim=(sind(zlat2(m-1,k1))*sind(zpair4)) dam=(cosd(zlat2(m-1,k1))*cosd(zpair4)) dam1=zlong2(m-1,k1)-zpair5 dam1=cosd(dam1) dist=dim+dam*dam1 dist=acosd(dist)*110.949 tdist(m,npair3)=dist tendency(m,npair3)=zpair6-p3(m-1,k1) c c test output c c write(31,5567) zjul(m),nsys(m,npair3) c5567 format(f5.1,i7) c print *,zpair4,zpair5,zpair6,nsys(m,npair3),npair3 c print *,tdist(m,npair3),tendency(m,npair3) c endif endif end do c c*************************************************************** c retest to make sure above situations dealt with c*************************************************************** c 2626 do k=1,icount(m) mdup=0 if(nsys(m,k) .eq. 0) then do k1=1,icount(m-1) if(i2(m-1,k1) .ge. il(k) .and. i2(m-1,k1) .le. ih(k) * .and. j2(m-1,k1) .ge. jl(k) .and. j2(m-1,k1) * .le. jh(k) .and. status(m-1,k1) .eq. 0) then mdup=mdup+1 endif end do c if(mdup .gt. 1) then print *,'multiple chart m-1 cyclones in' print *,'chart m shell, chart = ',ncharts endif endif end do c do k1=1,icount(m-1) ndup=0 if(status(m-1,k1) .eq. 0) then do k=1,icount(m) if(i2(m-1,k1) .ge. il(k) .and. i2(m-1,k1) .le. ih(k) * .and. j2(m-1,k1) .ge. jl(k) .and. j2(m-1,k1) * .le. jh(k) .and. nsys(m,k) .eq. 0) then ndup=ndup+1 endif end do if(ndup .gt. 1) then print *,'chart m-1 cyclone found in multiple' print *,'chart m shells, chart = ',ncharts endif endif end do c c*************************************************************** c now deal with remaining pairings based on 3x3 c search shell c************************************************************** do k=1,icount(m) if(nsys(m,k) .eq. 0) then do k1=1,icount(m-1) if(i2(m-1,k1) .ge. il(k) .and. i2(m-1,k1) .le. ih(k) * .and. j2(m-1,k1) .ge. jl(k) .and. j2(m-1,k1) * .le. jh(k) .and. i2(m,k) .ne. 0) then if(status(m-1,k1) .eq. 0) then nsys(m,k)=nsys(m-1,k1) tendency(m,k)=p3(m,k)-p3(m-1,k1) status(m-1,k1)=1 c dim=(sind(zlat2(m,k))*sind(zlat2(m-1,k1))) dam=(cosd(zlat2(m,k))*cosd(zlat2(m-1,k1))) dam1=zlong2(m,k)-zlong2(m-1,k1) dam1=cosd(dam1) dist=dim+dam*dam1 dist=acosd(dist)*110.949 tdist(m,k)=dist c endif endif end do endif end do c****************************************************************** c if all systems on previous day accounted for, c move to the next day c****************************************************************** do k1=1,icount(m-1) if(status(m-1,k1) .eq. 0) goto 1178 end do c c PRINT *,'ALL PREVIOUS SYSTEMS ACCOUNTED FOR, c * AT DAY',m c goto 4444 c c*********************************************************** c*********************************************************** c do futher search using nearest neighbor approach c*********************************************************** c*********************************************************** c c Find distances between untracked cyclones on chart m c (nsys(m,k) = 0) and cyclones on chart m-1 not already c paired with chart m system from previous check from the c 3x3 shell (status(k-1,k1)=0). c c1178 print *,'distances' C 1178 do 2200 k=1,icount(m) if(nsys(m,k) .eq. 0) then do 2195 k1=1,icount(m-1) C if(status(m-1,k1) .eq. 0) then dim=(sind(zlat2(m,k))*sind(zlat2(m-1,k1))) dam=(cosd(zlat2(m,k))*cosd(zlat2(m-1,k1))) dam1=zlong2(m,k)-zlong2(m-1,k1) dam1=cosd(dam1) dist=dim+dam*dam1 dist=acosd(dist)*110.949 dist1(k,k1)=dist c c find longitude changes of potential storm pairings c with respect to reference latitude (for pairs c south of 85 deg. north) c if(zlat2(m-1,k1) .le. 85 .and. * zlat2(m,k) .le. 85) then c c average latitude of storm pairs c rotten=(zlat2(m,k)+zlat2(m-1,k1))/2 c c length of deg. longitude at ave. latitude c rotten1=cosd(rotten)*110.949 c c change in longitude of potential storm pairings c q2=zlong2(m,k) q1=zlong2(m-1,k1) rotten2=q2-q1 if((q2 .gt. 270 .and. q2 .le. 360) .and. (q1 .ge. 0 * .and. q1 .lt. 90)) rotten2=q2-(q1+360) if((q2 .ge. 0 .and. q2 .lt. 90) .and. (q1 .gt. 270 * .and. q1 .le. 360)) rotten2=(q2+360)-q1 c c distance in km for longitude change at reference latitude c dlong(k,k1)=rotten1*rotten2 c dlong2(k,k1)=rotten2 dlat2(k,k1)=rotten tojo(k,k1)=q2 tojo1(k,k1)=q1 else dlong(k,k1)=999 dlong2(k,k1)=999 dlat2(k,k1)=999 tojo(k,k1)=999 tojo1(k,k1)=999 endif c c find pressure tendency of potential storm pairings c dp(k,k1)=p3(m,k)-p3(m-1,k1) c c find latitude change of potential storm pairings c dlat(k,k1)=zlat2(m,k)-zlat2(m-1,k1) c c find storm number match of potential pairing c nsys1(k,k1)=nsys(m-1,k1) c c test output c c PRINT 716,L,K,K1,NSYS1(K,K1),DIST1(K,K1),DLAT(K,K1), c * dlong(k,k1),dlong2(k,k1),tojo(k,k1),tojo1(k,k1), c * dlat2(k,k1) c716 FORMAT(4I4,7F8.1) endif 2195 continue endif 2200 continue c c*********************************************************** c Search potential pairings, find minimum distances c of chart m systems with respect to those systems c on chart m-1 not already paired with chart c m system. If chart m system failed other c ad hoc checks, set distance to 6000 km so that it c will never show its minimum distance to that system. c THIS IS THE SECTION WHERE THE PRAMETERS SET AT THE C TOP OF THE PROGRAM ARE USED. c*********************************************************** c print *,'shortest distances' c do k=1,icount(m) if(nsys(m,k) .eq. 0) then not=not+1 idog(not)=k zmin=100000 do k1=1,icount(m-1) if(status(m-1,k1) .eq. 0) then c c limit of total distance moved c if( dist1(k,k1).gt. maxdist) dist1(k,k1)=6000 c c limit of westward motion c if(dlong(k,k1) .le. westmax) dist1(k,k1)=6000 c c limit of southward motion c if((dlat(k,k1).lt.0.).and. & ((dlat(k,k1)*110.949).le.southmax)) dist1(k,k1)=6000 c c limit of northward motion c if((dlat(k,k1).gt.0.).and. & ((dlat(k,k1)*110.949).gt.northmax)) dist1(k,k1)=6000 c c limit of pressure tendency c if(abs(dp(k,k1)) .gt. pchange) dist1(k,k1)=6000 c if(dist1(k,k1) .lt. zmin) then zmin=dist1(k,k1) nsys2=nsys1(k,k1) dlat1=dlat(k,k1) dp1=dp(k,k1) dl4=dlong(k,k1) endif endif end do zmin1(k)=zmin nsys3(k)=nsys2 dlat3(k)=dlat1 dp2(k)=dp1 dl5(k)=dl4 c print 4067,m,m,nsys3(k),zmin1(k),dlat3(k),dp2(K), c * dl5(k) c4067 format(3i4,4x,f8.1) endif end do C C if all systems on previous day tracked, skip C to next day C if(not .eq. 0) goto 4444 C C find minimum and maximum storm numbers on chart m-1 C min=100000 do j=1,not if(nsys3(idog(j)) .le. min) then min=nsys3(idog(j)) endif end do max=0 do j=1,not if(nsys3(idog(j)) .ge. max) then max=nsys3(idog(j)) endif end do c c PRINT *,'NOT,MIN/MAX VALUE',NOT,MIN,MAX c**************************************************************** c If only one system on chart m is untracked, c match with closest system on chart m-1. c*************************************************************** c j=1 if(not .eq. 1) then if(zmin1(idog(j)) .le. maxdist) then c if(abs(dp2(idog(j))) .lt. 20 .and. c * (dlat3(idog(j)) .ge. -4 .and. dlat3(idog(j)) c * .le. 5)) then nsys(m,idog(j))=nsys3(idog(j)) tdist(m,idog(j))=zmin1(idog(j)) dlat6(m,idog(j))=dlat3(idog(j)) dlong6(m,idog(j))=dl5(idog(j)) tendency(m,idog(j))=dp2(idog(j)) c PRINT *,'ONLY 1 UNTRACKED, ASSIGN',L,ZMIN1(IDOG(J)), c * NSYS(m,IDOG(J)),NSYS3(IDOG(J)) c endif endif not=0 c goto 4444 endif c c***************************************************************** c If here, multiple systems on chart m are untracked. c Often, two or more intracked chart m systems find c their minimum distances with respect to the same chart c m-1 system. In these cases, pair the chart m-1 system c with the closest chart m system, using same critera as above. c If does not pass, the chart m system is new (cyclogenesis) and c sucessively further chart m systems are inspected in the c same manner, up to the search km limit. In simpler case c where sinle chart m system has a minimum distance to single c chart m-1 system, pair it with the m-1 system, provided is c passes the criteria set above. c***************************************************************** c istart=min 5001 zcut=maxdist C do 5000 j=1,not if(nsys3(idog(j)) .eq. istart) then if(zmin1(idog(j)) .lt. zcut) then ib=j zcut=zmin1(idog(j)) iobs=iobs+1 endif else goto 5000 endif 5000 continue if(iobs .ge. 1) then nsys(m,idog(ib))=nsys3(idog(ib)) tdist(m,idog(ib))=zmin1(idog(ib)) tendency(m,idog(ib))=dp2(idog(ib)) dlat6(m,idog(ib))=dlat3(idog(ib)) dlong6(m,idog(ib))=dl5(idog(ib)) iobs=0 endif 3334 if(istart .ge. max) goto 4444 istart=istart+1 goto 5001 c c**************************************************************** c Now at end of cyclone search. Ascribe new c cyclone numbers to untracked chart m c systems. c**************************************************************** c 4444 not=0 c c print *,'down to here' c c find the highest numbered system c do k=1,icount(m) if(nsys(m,k) .ge. ilast) then ilast=nsys(m,k) endif end do c c search for remaining untracked chart m systems c (nsys(m,k)=0) if found, define as cyclogenesis events c and ascibe new system numbers (cyclogenesis only if c previous day was not skipped, i.e., iskip(m)=0) do k=1,icount(m) if(nsys(m,k) .eq. 0) then ilast=ilast+1 nsys(m,k)=ilast tdist(m,k)=-999 tendency(m,k)=-999 if(iskip(m) .eq. 0) cyclo(m,k)=1 endif end do c c if systems on previous day were determined c to have filled, define as cyclolysis events. c (only if iskip(m)=0) c if(m .gt. 1 .and. iskip(m) .eq. 0) then do 450 k1=1,icount(m-1) do 445 k=1,icount(m) if(nsys(m,k) .eq. nsys(m-1,k1)) goto 450 445 continue lysis(m-1,k1)=1 450 continue endif c c**************************************************************** c**************************************************************** c c write output into unit 5. Variables are c c 1) julian day c 2) record number (ignore this variable) c 3) year c 4) month c 5) day of month c 6) hour (0000 or 1200 UTC) c 7) total number of systems for that day c 8) number of grid points defining the central pressure c 9) flag to indicate if previous day(s) were skipped c (yes=1, no=0) c 10) cyclone central pressure (mb) c 11) local laplacian of pressure (system intensity) c 12) distance traveled from last observations (-999 c when iskip=1 or is cyclogenesis event) c 13) pressure tendency from last observation (-999 when c iskip=1 or is cycogenesis event) c 14) latitude of system center c 15) longitude of system center (0 to 360 deg.) c 16) NMC grid row c 17) NMC grid column c 18) flag for cyclogenesis event (yes=1,no=0) c 19) flag for cyclolysis event (yes=1, no=0) c 20) system number c 4446 if(m .gt. 1) then c print *,m, icount(m-1) do k1=1,icount(m-1) irecord=(zjul(m-1)*4)-3 write(5,533) zjul(m-1),irecord,iyr1(m-1),imon1(m-1),iday1(m-1), * ihr1(m-1),icount(m-1),tdiff1(m-1,k1),iskip(m-1),p3(m-1,k1), * lap1(m-1,k1),tdist(m-1,k1),tendency(m-1,k1),zlat2(m-1,k1), * zlong2(m-1,k1),i2(m-1,k1),j2(m-1,k1), * cyclo(m-1,k1),lysis(m-1,k1),nsys(m-1,k1) end do endif 533 format(F7.2,8I5,6F8.2,4I3,i7) c c**************************************************************** c**************************************************************** 100 CONTINUE c c sample output of pressure arrays for plotting c c do i=1,ieqlon c do j=1,jeqlat c write(32,1053) i,j,lat(i,j),long(i,j), c * (psamp(n,i,j),n=1,10) c end do c end do c1053 format(2i3,12f8.2) close (5) 99 STOP END C----------------------------------------------------------------------C C This function returns the difference (in degrees longitude) between C the input longitudes. The output is always positive and less than C or equal to 180.0 degrees because the routine is designed to do the C difference the short way around the globe. The input longitudes are C assumed to be constrained to a range of 360 degrees. C----------------------------------------------------------------------C function diflong(xlon1,xlon2) diflong = abs(xlon2 - xlon1) if(abs(xlon1-xlon2).le.180.0) return if(xlon2 .lt. xlon1) then diflong = abs(xlon2 - xlon1 + 360.0) else diflong = abs(xlon2 - xlon1 - 360.0) endif return end C----------------------------------------------------------------------C