C *****************************************************************
C PROGRAM: monthly.f
C
C LANGUAGE: Fortran 77
C
C MACHINE: VM
C ----------------------------------------------------------------*
C PURPOSE: Generate 600-month (50-year) daily precipitation data ends on
C last month
C ----------------------------------------------------------------*
C USAGE:
C    monthly.x
C ----------------------------------------------------------------*
C INPUT FILES:
C    work/spi/spi_date - target date in yyyymmdd format
C    work/spi/p.conus.cur - up to date daily precipitation analysis
C    work/spi/archive - 30 and 60-day averaged precipitation
C    input/spi/mask.p.xie - land mask
C OUTPUT FILES:
C    work/spi/monthly.600m - 600-month daily precipitation data
C ----------------------------------------------------------------*
C SUBROUTINES USED:
C FUNCTIONS USED:
C ----------------------------------------------------------------*
C INPUT VARIABLES:
C LOCAL VARIABLES:
C ----------------------------------------------------------------*
C AUTHOR(S): Kingtse Mo
C ----------------------------------------------------------------*
C DATE 2000-01-01
C ----------------------------------------------------------------*
C MODIFICATIONS:
C      2018-03-01 Hai-Tien Lee R2O standard compliance
C *****************************************************************

      parameter (jm=60,im=150)
      dimension a(150,60),p(150,60),pd(150,60),iday(24)
      dimension fg(150,60),isum(12),pmonth(150,60)
      dimension  vmask(150,60)
      integer endp
      open(10,file='work/spi/spi_date')
      open(11,file='input/spi/p.conus.cur',
     c   access='direct',recl=150*60*4,form='unformatted')
      open(12,file='work/spi/archive',
     c   access='direct',recl=150*60*4,form='unformatted')
      open(15,file='input/spi/mask.p.xie',
     c   form='unformatted')
      open(20,file='work/spi/monthly.600m',
     c   access='direct',recl=150*60*4,form='unformatted')
      data iday/31,28,31,30,31,30,31,31,30,31,30
     c   ,31,31,28,31,30,31,30,31,31,30,31,30,31/
      data isum/0,31,59,90,120,151,181,212,243,273,304,334/

c read in the ending date
c convert to month and day

c the time series up to the end of 2008 is 21550
c read in the last day needed to calculate monthly means
      read(15) vmask
      read(10,*) jdate
      one=1.
      mdate=jdate/100
      print *,mdate
      myr=mdate/10000
      mmo=(mdate-myr*10000)/100
      mm1=(mdate-myr*10000)/100-1
      if(mm1.eq.0) then
         mm1=12
      endif

c mdd is the ending date for monthly means
      mdd=(mdate-myr*10000-100*mmo)
      print *,myr,mmo,mm1,mdd
      mystart=myr-50
      print *,mystart
      kkstart=mystart-1
      jul=isum(mmo)+mdd
      print *,' jul date ',jul, ' the end year ',kkstart

c calculate the starting date
      irec=0
      do kk=1950,kkstart
         irec=irec+365
         jq=kk/4
         ddm=float(kk)/4
         if(jq.eq.ddm) irec=irec+1
      enddo
      print *,'  days from years ',irec
      irec=irec+jul
      mstart=mmo
      if(mstart.le.0.) mstart=mstart+12
      mend=mstart+11
      print *,jul,irec,mstart,mend
      iout=0
      kyr=myr-1
      do iy=mystart,kyr
         do mm=mstart,mend
            limit=iday(mm)
            print *,'  iyr ',iy,' mm ',mm, 'irec', irec
            print *,'  limit ', limit
            if((mm.eq.2).or.(mm.eq.14)) then
               jq=iy/4
               ddm=float(iy)/4
               if(jq.eq.ddm) limit=29
               print *,iy,jq,ddm
            endif
            print *,'  limit ', limit

            do j=1,jm
               do i=1,im
                  pmonth(i,j)=0.
                  fg(i,j)=0.
               enddo
            enddo

            do nday=1,limit
               irec=irec+1
               read(11,rec=irec) p
               do j=1,jm
                  do i=1,im
                     if(abs(p(i,j)).le.100.) then
                        if(vmask(i,j).eq.1)  then
                           fg(i,j)=fg(i,j)+one
                           pmonth(i,j)=pmonth(i,j)+p(i,j)
                        endif
                     endif
                  enddo
               enddo
            enddo

            do j=1,jm
               do i=1,im
                  pd(i,j)=-999.
                  if(fg(i,j).ne.0.) pd(i,j)=pmonth(i,j)/fg(i,j)
               enddo
            enddo
            iout=iout+1
            write(20,rec=iout) pd
            print *,iy,mm,iout,fg(60,30),pd(60,30),irec
         enddo
      enddo

      do kk=1,2
         read(12,rec=kk) p
         do j=1,jm
            do i=1,im
            pd(i,j)=-999.
            if(abs(p(i,j)).le.100.) then
               if(vmask(i,j).eq.1) then
                  pd(i,j)=p(i,j)
               endif
            endif
            enddo
         enddo
         iout=iout+1
         write(20,rec=iout) pd
      enddo
      stop
      end