      PROGRAM onset_demise
      IMPLICIT NONE
      INTEGER nlon,nlat,nday,nyear,x,y,t,tt,yy,dd,total,lrec
      INTEGER i,j,k,irec,undef,leap,ly,tot,tot1,ll,l,t1,ll1,ll2
      PARAMETER (x=329,y=347,t=366,total=3652,yy=10,dd=181,t1=275)
      INTEGER iu, k1, m, n
      REAL intmin(yy),intmax(yy), rawp(t), smtp(t)
      REAL aa(x,y,total),cc(x,y,total),MX,MN
      REAL var1(x,y,t,yy),sp_av(t,yy)
      REAL store
      REAL anmn_grd(x,y,yy),anmn_spav(yy),climo_grd(x,y),climo_spav
      REAL danom_grd(x,y,t,yy),canom_grd(x,y,t,yy)
      REAL canom_spav(t,yy),danom_spav(t,yy)
      REAL onset_grd(x,y,yy),peak_grd(x,y,yy),demise_grd(x,y,yy)
      REAL sl_grd(x,y,yy),mask(x,y),acc(x,y,yy)
      REAL po_grd(x,y,yy),pd_grd(x,y,yy),pod_grd(x,y)
      REAL ao_grd(x,y,yy),ad_grd(x,y,yy)
      CHARACTER*80 name(yy), yr(yy),nameA(yy)
      REAL a1(x,y),a2(x,y),oslp(x,y)
      REAL b1(x,y),b2(x,y),dslp(x,y)
      REAL c1(x,y),c2(x,y),sslp(x,y)
      REAL clm_ons(x,y),clm_dms(x,y),clm_sl(x,y)
      
!!!!!!!!!!!!!!!!!!!!!!!!!!!RETREIVING VARIABLE!!!!!!!!!!!!!!!!!!!!!!!!

!      In shell(not program), create ordered file listing all
!      necessary datasets(each "nyear") with ls -v *.grd > list_all
!      Open file(1001 identity is arbitrary)

       Open(1122,file='list-MM10',status='old')
       Open(1121,file='list-aMM10',status='old')

       Do i = 1 , yy
         READ(1122,*) name(i)
         READ(1121,*) nameA(i)
!         print*, name(i)
       Enddo 

       open(3113,file='accES.grd',form='unformatted',
     & access='direct',status='unknown',recl=x*y)
       open(3114,file='slgthES.grd',form='unformatted',
     & access='direct',status='unknown',recl=x*y)
       open(3115,file='onsetES.grd',form='unformatted',
     & access='direct',status='unknown',recl=x*y)
       open(3116,file='dmiseES.grd',form='unformatted',
     & access='direct',status='unknown',recl=x*y)

       open(9991,file='ConsES.grd',form='unformatted',
     & access='direct',status='unknown',recl=x*y)
       open(9992,file='CdmsES.grd',form='unformatted',
     & access='direct',status='unknown',recl=x*y)
       open(9993,file='CslES.grd',form='unformatted',
     & access='direct',status='unknown',recl=x*y)

!!!!!!!!!!!!!!!!MM10 DATA 1986 to 1995 DATA!!!!!!!!!!!!!!!!
         ly = 0
         k1 = 1
        Do iu = 1 , yy
         OPEN(iu,file=name(iu),form='unformatted',
     &   access='direct',status='old',recl=x*y)
         print*, name(iu)
         
         IF(iu.EQ.1+ly*4) then
            print*,ly,iu
           ly = ly + 1
           tt = 366
         ELSE
           tt = 365
         ENDIF

         irec = 1
        DO k = 1 , tt-1
!         print*,irec
         READ(iu,rec=irec)((aa(i,j,k),i=1,x),j=1,y)
         irec = irec + 1
         ENDDO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!  Akhilesh Mishra !!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!  11 point smoothing !!!!!!!!!!!!!!!!!!!!!!!!!
!!! call subroutine to smooth daily pcp  !!!!!!!
          Do nlat = 1 , y
           Do nlon = 1 , x
             Do nday = 1 , tt
             rawp(nday) = aa(nlon,nlat,nday)
             ENDDO
            call MovingAverage(rawp,smtp,tt,10)
             Do nday = 1 , tt
             aa(nlon,nlat,nday) = smtp(nday)
             ENDDO
           ENDDO
          ENDDO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
          DO k = 1, tt
          cc(:,:,k1) = aa(:,:,k)
          k1 = k1 + 1
          ENDDO

       Do k = 1 , total
        Do j = 1 , y
         Do i = 1 , x
           IF(abs(cc(i,j,k)).GT.500.0) cc(i,j,k) = 0.0
         ENDDO
        ENDDO
       ENDDO
      ENDDO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ly = 0   !1986 start year 
        ll = 1   ! will count yy*days in year 
        Do nyear = 1 , yy

!!!!!!!LEAP YEAR DAY !!!!!!!!!!!!!!!
         IF(nyear.EQ.1+ly*4) then
            ly = ly + 1
            print*,ly,nyear,1985+nyear
           tt = 366
         ELSE
           tt = 365
         ENDIF
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         Do nday = 1 , tt
            ll1 = 0
          Do nlat = 1 , y
            Do nlon = 1 , x
            var1(nlon,nlat,nday,nyear) = cc(nlon,nlat,ll)
       IF(cc(nlon,nlat,ll).EQ.0.0) then 
          var1(nlon,nlat,nday,nyear) = 0.0
            ll1 = ll1 + 1 
       ENDIF
            ENDDO
          ENDDO
             ll = ll + 1
         ENDDO
        ENDDO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!      Compute annual mean variable at each grid point and spatially 
!      averaged
       ly = 0
       anmn_grd = 0.0
       DO nyear = 1 , yy 

        IF(nyear.EQ.1+ly*4) THEN
          leap = nyear
          ly = ly + 1
         ENDIF
         IF(nyear .EQ. leap) THEN
           tt=366
         ELSE
           tt=365
         ENDIF

         DO nlat = 1, y
           DO nlon = 1, x  
             DO nday = 1, tt
               anmn_grd(nlon,nlat,nyear) = anmn_grd(nlon,nlat,nyear) + 
     &         var1(nlon,nlat,nday,nyear)
             ENDDO
             anmn_grd(nlon,nlat,nyear) = 
     &       anmn_grd(nlon,nlat,nyear)/FLOAT(tt)
           ENDDO
         ENDDO
       ENDDO
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!      Compute climatological annual mean variable at each grid point 
!      and spatially averaged
       climo_grd = 0.0
       DO nlat = 1, y
         DO nlon = 1, x
           DO nyear = 1, yy
              climo_grd(nlon,nlat) = climo_grd(nlon,nlat) +
     &        anmn_grd(nlon,nlat,nyear)
           ENDDO
           climo_grd(nlon,nlat) = climo_grd(nlon,nlat)/FLOAT(nyear)
         ENDDO
       ENDDO

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!      Compute annual daily and cumulative variable anomalies at each 
!      grid point and spatially averaged (USEFUL FOR ANALYSIS OF ONSET 
!      AND DEMISE DATES)
        ly = 0
       DO nyear = 1, yy

         IF(nyear.EQ.1+ly*4) then
           leap = nyear
           ly = ly + 1
         ENDIF

         IF(nyear .EQ. leap) THEN
           tt=366
         ELSE
           tt=365
         ENDIF

         DO nlat = 1, y
           DO nlon = 1, x
             store = 0.0
             DO nday = 1, tt
               danom_grd(nlon,nlat,nday,nyear) = 
     &         var1(nlon,nlat,nday,nyear) -  climo_grd(nlon,nlat)
               store = store + danom_grd(nlon,nlat,nday,nyear)
               canom_grd(nlon,nlat,nday,nyear) = store
             ENDDO
           ENDDO
         ENDDO
       ENDDO

!! Write daily anomaly 

           ly = 0
       DO nyear = 1, yy
         OPEN(nyear,file=nameA(nyear),form='unformatted',
     &   access='direct',status='unknown',recl=x*y)
!
         IF(nyear.EQ.1+ly*4) then
           leap = nyear
           ly = ly + 1
         ENDIF

         IF(nyear .EQ. leap) THEN
           tt=366
         ELSE
           tt=365
         ENDIF

           lrec = 0
         DO nday = 1, tt
         lrec = lrec + 1
            pod_grd(:,:) = canom_grd(:,:,nday, nyear)
         WRITE(nyear,rec=lrec)((pod_grd(i,j),i=1,x),j=1,y)
         ENDDO
         close(nyear)
       ENDDO

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!ONSET AND DEMISE!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!      Determine annual onset and demise dates at each grid point
!!      print*,'t-92=', t-92
       DO nyear = 1, yy
         DO nlat = 1, y
           DO nlon = 1, x

!   Get the demise date
                intmax(nyear)=-9.99e+08
!             DO nday = onset_grd(nlon,nlat,nyear), t-1
             DO nday = 275, t-1
               IF(canom_grd(nlon,nlat,nday,nyear).GT.intmax(nyear)) THEN
                intmax(nyear) = canom_grd(nlon,nlat,nday,nyear)
                demise_grd(nlon,nlat,nyear) = nday + 1
               ENDIF
             ENDDO

! Onset should be minimum 60 days before demise
               
                intmin(nyear)=+9.99e+08
!             DO nday = 90, t-92
             DO nday = 90, demise_grd(nlon,nlat,nyear)-60
               IF(canom_grd(nlon,nlat,nday,nyear).LT.intmin(nyear)) then
                intmin(nyear)= canom_grd(nlon,nlat,nday,nyear)
                onset_grd(nlon,nlat,nyear) = nday + 1
               ENDIF
             ENDDO
 
 
                intmax(nyear)=-9.99e+08
        DO nday = onset_grd(nlon,nlat,nyear),demise_grd(nlon,nlat,nyear)
               IF(danom_grd(nlon,nlat,nday,nyear).GT.intmax(nyear)) THEN
                intmax(nyear) = danom_grd(nlon,nlat,nday,nyear)
                peak_grd(nlon,nlat,nyear) = nday
               ENDIF
        ENDDO
             
             IF (onset_grd(nlon,nlat,nyear) .GE. 365) THEN 
               onset_grd(nlon,nlat,nyear) = 0.0
            ENDIF

             IF (peak_grd(nlon,nlat,nyear) .GE. 365) THEN 
               peak_grd(nlon,nlat,nyear) = 0.0
             ENDIF

!             IF (demise_grd(nlon,nlat,nyear) .GE. 365) THEN 
!               demise_grd(nlon,nlat,nyear) = 0.0
!             ENDIF

           ENDDO
         ENDDO
       ENDDO 
     
            sl_grd = demise_grd - onset_grd

         DO nlat = 1, y
           DO nlon = 1, x
            clm_ons(nlon,nlat)  = 0.0
            clm_dms (nlon,nlat) = 0.0
            clm_sl (nlon,nlat) = 0.0
            DO nyear = 1, yy
         clm_ons(nlon,nlat)=clm_ons(nlon,nlat)+onset_grd(nlon,nlat,yy)
         clm_dms(nlon,nlat)=clm_dms(nlon,nlat)+demise_grd(nlon,nlat,yy)
         clm_sl(nlon,nlat)=clm_sl(nlon,nlat)+sl_grd(nlon,nlat,yy)
            ENDDO
         clm_ons(nlon,nlat) = clm_ons(nlon,nlat) / real(yy)
         clm_dms(nlon,nlat) = clm_dms(nlon,nlat) / real(yy)
         clm_sl(nlon,nlat) = clm_sl(nlon,nlat) / real(yy)
           ENDDO
         ENDDO
          
!!!!!!!!!!!!!!!!!!!Accumulated rainfall !!!!!!!!!!!!!!!!!!!
        acc = 0.0
       DO nyear = 1, yy
         DO nlat = 1, y
          DO nlon = 1, x
        DO nday = onset_grd(nlon,nlat,nyear),demise_grd(nlon,nlat,nyear)
           acc(nlon,nlat,nyear) = acc(nlon,nlat,nyear) +
     &        var1(nlon,nlat,nday,nyear)
         ENDDO

         ENDDO
        ENDDO
       ENDDO

       irec = 1
       DO nyear = 1, yy
         WRITE(3113,rec=irec)((acc(i,j,nyear),i=1,x),j=1,y)
         WRITE(3114,rec=irec)((sl_grd(i,j,nyear),i=1,x),j=1,y)
         WRITE(3115,rec=irec)((onset_grd(i,j,nyear),i=1,x),j=1,y)
         WRITE(3116,rec=irec)((demise_grd(i,j,nyear),i=1,x),j=1,y)
         irec = irec + 1
       ENDDO

         WRITE(9991,rec=1)((clm_ons(i,j),i=1,x),j=1,y)
         WRITE(9992,rec=1)((clm_dms(i,j),i=1,x),j=1,y)
         WRITE(9993,rec=1)((clm_sl(i,j),i=1,x),j=1,y)
       
      STOP
      END
      Subroutine  MovingAverage(X,Y,MAX_SIZE,WINDOW)
      IMPLICIT NONE
      INTEGER            :: MAX_SIZE       ! array size
      REAL, DIMENSION(1:MAX_SIZE) :: X, Y     ! arrays
      REAL               :: Sum                 ! for computation use
      INTEGER            :: WINDOW                   ! window size
C      INTEGER            :: Size                ! actual array size
      INTEGER            :: i, j                ! indices
      DO i = WINDOW, MAX_SIZE-WINDOW+1       ! for each xi
       Sum = 0.0                        ! compute the moving average
      DO j = i-WINDOW+1, i+WINDOW-1                   ! of xi, x(i+1), ...,
         Sum = Sum + X(j)                    ! x(i+Window-1).
      END DO
        Y(i) = Sum / (1+(2*WINDOW))                  ! save the result
      END DO
      DO i = 1, WINDOW 
          Y(i) = X(i)
      ENDDO
      DO i = MAX_SIZE-WINDOW+1,  MAX_SIZE
          Y(i) = X(i)
      ENDDO

C ---      WRITE(*,*)  (Y(i), i = 1, MAX_SIZE)
!      WRITE(*,*)  (Y(11))
!      WRITE(*,*)  (X(i), i = 1, 21)
      return
      end subroutine MovingAverage
