program spi_obs

c This code computes observed SPI for July 1, 2022 (doy=182) using NLDAS-2 daily precip
c of most recent 41 years (excluding Feb 29 of leap years) as the base. 

c To compile: gfortran -O spi.sample.f spi_daily.sub.f

      real imiss
      integer yrbegin,yrend,yyyy,ydiff,doy
      parameter(yrbegin=1980,yrend=2020,maxyrs=yrend+1-yrbegin)
      parameter(yyyy=2022,ydiff=yyyy-yrend,doy=182)
      parameter(limit=maxyrs*365,limit0=limit+365*ydiff)
      parameter(ratio=24000.0)  !convert Noah P output to mm/day

c im, jm: # of grids in x and y directions
C nlen:   # of SPIs to compute
      parameter(im=117,jm=49,imiss=-999.0,nlen=5)

      character*99 fname
      integer dstart,dtotal
      dimension len(nlen)

c len: timescales of SPI in terms of number of days
c 31,91,182,365,730 respectively correspond to 1, 3, 6, 12 and 24 months
      data len/31,91,182,365,730/

      real aa(im,jm),p(im,jm,limit),pp(im,jm,limit0)
      real mask(im,jm)

      real prec(limit),beta(365),gamm(365),pzero(365)
      real spi(limit,nlen)
      real b(im,jm,limit),r(im,jm)

c input: land mask (the SPI will only be computed for land grids)
      open(10,file='/cpc/home/hawang/bin/noahmp.USmask3439',
     *  form='unformatted',access='direct',recl=im*jm*4)
      read(10,rec=1) mask      
      close(10)

c output: SPI output for time scales of 31,91,182,365,730 days
c         (i.e., 1, 3, 6, 12 and 24 months)
      open(21,file='@varout1361224.noah.opn.@year@mon@day',
     *  form='unformatted',access='direct',recl=im*jm*4)

c loop: year
        irec=0
        do iyr=yrbegin,yyyy

c input: daily precip for each year (excluding Feb 29)
           write(fname,'("@var_daily.noah.opn.",i4)') iyr

           open(1,file=
     *"/cpc/drought/fcst/Noah/conus_0p5/out/opn/"//fname,
     *form='unformatted',access='direct',recl=im*jm*4)
           
           ndays=365
           if(iyr.eq.yyyy) ndays=doy

           do iday=1,ndays
              irec=irec+1
              read(1,rec=iday) aa

              do j=1,jm 
              do i=1,im
                 pp(i,j,irec)=aa(i,j)
              enddo 
              enddo 

           enddo        !iday

           close(1)
        enddo           !iyr

c Only use the most recent daily data of maxyrs years (which in this code 
c is 365*41) for the SPI calculation
c For convenience, daily data for Feb 29 of leap years are not used
      dtotal=(yyyy-yrbegin)*365+doy
      dstart=dtotal-(maxyrs*365)+1

        print*,dtotal,dstart

      irec=0
      do iday=dstart,dtotal
         irec=irec+1
         do j=1,jm
         do i=1,im
            p(i,j,irec)=pp(i,j,iday)
         enddo
         enddo
      enddo

      do it=1,limit
         do j=1,jm
         do i=1,im
            b(i,j,it)=0.0
         enddo
         enddo
      enddo

c compute SPIs 
      do 10 ilen=1,nlen

         print*, 'compute SPI',ilen

         do 6000 jj=1,jm
         do 6000 ii=1,im
            if(mask(ii,jj).ge.0.5) then
               do it=1,limit 
                  prec(it)=p(ii,jj,it)
               enddo

               call spigam(len(ilen),prec,beta,gamm,pzero,spi(1,ilen))

c skip leading missings
               ifirst = 0
20             continue
               ifirst = ifirst + 1
               if(spi(ifirst,ilen).eq.imiss) goto 20

c skip trailing missings
               last = maxyrs*365 + 1
30             continue
               last = last - 1
               if(spi(last,ilen).eq.imiss) goto 30

c output SPIs
               do 40 i=ifirst,limit
                  if(spi(i,ilen).eq.imiss) goto 50
                  b(ii,jj,i)=spi(i,ilen)
40             continue
50             continue
            endif
6000     continue
                
         do j=1,jm
         do i=1,im
            r(i,j)=imiss
            if(mask(i,j).ge.0.5) then
            if(b(i,j,limit).ge.-6.0.and.b(i,j,limit).le.6.0) then
               r(i,j)=b(i,j,limit)
            endif
            endif
         enddo
         enddo
         write(21,rec=ilen) r
10    continue

      stop
      end