c       Program to determine cold frontal passages over South America 
c234567890123456789001234567890123456789012345678901234567890123456789012
      parameter (ndmax=10958,imax=29,jmax=29)
      dimension t925(imax,jmax,ndmax),v925(imax,jmax,ndmax)
      dimension slp(imax,jmax,ndmax),y(imax,jmax),x(imax,jmax)
      dimension cnt(imax,jmax),z(imax,jmax)
      character path*30
      data path/'/cpc/home/wd52vk/samer/fronts/'/
      open (11,file=path//'t925-SA-y1980-2009',
     *access='direct',form='unformatted',recl=imax*jmax*4)
      open (12,file=path//'slp-SA-y1980-2009',
     *access='direct',form='unformatted',recl=imax*jmax*4)
      open (13,file=path//'v925-SA-y1980-2009',
     *access='direct',form='unformatted',recl=imax*jmax*4)
      open (51,file=path//'SA-cf-dates-y1980-2009',
     *access='direct',form='unformatted',recl=28)
      do 999 i=1,imax
      do 990 j=1,jmax
      cnt(i,j)=0.
 990  continue
 999  continue
      do 99 iday=1,ndmax
      read(11,rec=iday)z
      read(12,rec=iday)x
      read(13,rec=iday)y
      do 20 i=1,imax
      do 21 j=1,jmax
      t925(i,j,iday)=z(i,j)
      slp(i,j,iday)=x(i,j)
      v925(i,j,iday)=y(i,j)
 21   continue
 20   continue
 99   continue
c     determine cold frontal passages based on temperature decreases
c       and meridional wind at 925 hPa
      write(6,201)
 201  format(2x,'nday',5x,'mind',9x,'nday',4x,'sum',1x,'v925',1x,'pdif',
     *1x,'cnt')
      nrec=0
      do 1 i=1,imax
      do 2 j=1,jmax
      iday=0
 998  iday=iday+1
      sum=0.
      tot=0.
      mday=iday
      if(iday.gt.ndmax)go to 997
      if(iday+1.gt.ndmax)go to 997
      tdif=t925(i,j,iday+1)-t925(i,j,iday)
      pdif=slp(i,j,iday+1)-slp(i,j,iday)
      v=v925(i,j,iday+1)
      if(i.eq.15.and.j.eq.9)write(6,202)iday,tdif,pdif,v
 202  format(1x,i8,3f6.1)
      if(tdif.lt.0..and.v.gt.0.)go to 30
      go to 998
 30   sum=sum+tdif
      tot=tot+pdif
      iday=iday+1
      if(iday.gt.ndmax)go to 997
      if(iday+1.gt.ndmax)go to 997
      tdif=t925(i,j,iday+1)-t925(i,j,iday)
      pdif=slp(i,j,iday+1)-slp(i,j,iday)
      if(v925(i,j,iday+1).gt.v)v=v925(i,j,iday+1)
      if(i.eq.15.and.j.eq.9)write(6,202)iday,tdif,pdif,v
      if(tdif.lt.0.)go to 30
      nday=iday
ccc   frontal passage criteria is:
ccc      925-hPa Temp decrease of at least 2C
ccc      925-hPa v >= 2 m/s
c     if(sum.le.-2..and.v.ge.2..and.tot.ge.0.)go to 31
      if(sum.le.-2..and.v.ge.2.)go to 31
      go to 998
 31   cnt(i,j)=cnt(i,j)+1.
      nrec=nrec+1
ccc   subroutine xdate computes the date (yyyymmdd)
      call xdate(mday,jdate)
      idate=jdate
      call xdate(nday,jdate)
      ndate=jdate
      write(51,rec=nrec)i,j,idate,ndate,sum,v,tot
      if(i.eq.15.and.j.eq.9)write(6,200)mday,nday,idate,ndate,sum,v,tot
     *,cnt(i,j)
c23456789012345678901234567890123456789012345678901234567890123456789012
 200  format(1x,2i10,2(2x,i10),2x,3f5.1,2x,f5.0)
      go to 998
 997  continue
c     if(sum.le.-2..and.v.ge.2..and.tot.ge.0.)go to 32
      if(sum.le.-2..and.v.ge.2.)go to 32
      go to 2
 32   cnt(i,j)=cnt(i,j)+1.
      nday=iday
      nrec=nrec+1
      call xdate(mday,jdate)
      idate=jdate
      call xdate(nday,jdate)
      ndate=jdate
      write(51,rec=nrec)i,j,idate,ndate,sum,v,tot
      if(i.eq.15.and.j.eq.9)write(6,200)mday,nday,idate,ndate,sum,v,tot
     *,cnt(i,j)
 2    continue
 1    continue
      stop
      end
cccc
ccc   subroutine xdate computes the date (yyyymmdd)
      subroutine xdate(nday,jdate)
      dimension ndmon(12)
      data ndmon/31,28,31,30,31,30,31,31,30,31,30,31/ 
      lpyr=1980
      mday=0
ccc   total number of years 30 (1980-2009)
      do 1 iyr=1,30
ccc   initial year is 1980
      iyear=1979+iyr
      do 2 mon=1,12
      maxd=ndmon(mon)
      if(mon.eq.2.and.iyear.eq.lpyr)maxd=maxd+1
      do 3 iday=1,maxd
      jdate=iyear*10000+mon*100+iday
      mday=mday+1
      if(nday.eq.mday)go to 50
 3    continue
 2    continue 
      if(iyear.eq.lpyr)lpyr=lpyr+4
 1    continue
 50   continue
      return
      end