program rrfcalc 

!----------------------------------------------------------------------
!   RRFCALC 3.0
!----------------------------------------------------------------------

c *** Variable declarations
      implicit none
      
      REAL    :: XCELL, YCELL,rlon0,rlat0,tlat1,tlat2
      REAL    :: matsrrf,days,ppb,bdv,fdv,fdvmax,contrib,bdvmax
      REAL    :: lat,long,dx,xorg,yorg,conc,rcut,temprrf,tempcontrib

      INTEGER :: NCOLS, NROWS, ndays, daylist(10)
      integer    c, r, k, m, n, cells, ii, jj, nsil, nthreshold
      integer :: i,j,ijcell,refcell,date,nxc,nyc,nzc
      real :: usercutoff, userthreshold, usernumber, rrf
      real :: xbase, xfuture, number, cut, value, maxval,abscontrib
      REAL, ALLOCATABLE :: base(:,:,:),future(:,:,:)
      INTEGER, ALLOCATABLE :: today(:),xdayused(:)

      CHARACTER(LEN=9) :: monid,dvtype
      CHARACTER(LEN=16) :: progname,st,cy
      CHARACTER(LEN=16) :: met, outfile
      CHARACTER(LEN=25) :: MESG
      CHARACTER(LEN=150) :: ename,fname
      CHARACTER(LEN=200) :: ifile,ipath
      CHARACTER(LEN=2) :: blank
 
      CHARACTER(LEN=16), ALLOCATABLE :: vnames(:)

      INTEGER :: istatus

      INTEGER :: TRIMLEN
      real       envreal

      EXTERNAL TRIMLEN, envreal



c *** Initialize variables
      progname = 'RRFCALC'
c      outfile  = 'OUTFILE'

c--------Read standard input for grid and filename information

      read(*,'(10x,a)') MESG 
      write(*,*)'Note for output file',MESG

      read(*,'(10x,a)') ifile
      read(ifile,*) usercutoff,userthreshold,usernumber,cells
      write(*,*)'User specified RRF values for cutoff,threshold,N,cell array:',usercutoff,userthreshold,usernumber,cells

      read(*,'(10x,a)') dvtype
      write(*,*) 'RRFs will be applied to user selected design value option:',dvtype

      read(*,'(10x,a)') ifile
      read(ifile,*) ndays
      write(*,'(a,t20,3i10)') 'Number of days on the file',ndays

      read(*,'(10x,a)') ifile
      read(ifile,*) nxc,nyc
      nzc = 1
      write(*,'(a,t20,3i10)') 'Grid size',nxc,nyc,nzc

      read(*,'(10x,a)') ifile
      read(ifile,*) dx,xorg,yorg
      write(*,'(a,t20,3f10.0)') 'Grid spacing and origin',dx,xorg,yorg

      read(*,'(10x,a)') ifile
      read(ifile,*) rlon0,rlat0,tlat1,tlat2
      write(*,'(a,t20,4f10.0)') 'Projection:',rlon0,rlat0,tlat1,tlat2

      read(*,'(10x,a)') ifile
      open(9,file=ifile,form='formatted')
      write(*,*) 'Output file 1: ',ifile

      read(*,'(10x,a)') ifile
      open(19,file=ifile,form='formatted')
      write(*,*) 'Output file 2: ',ifile

      read(*,'(10x,a)') ifile
      open(29,file=ifile,form='formatted')
      write(*,*) 'Output file 3 (unused monitors): ',ifile

      read(*,'(10x,a)') ifile
      open(20,file=ifile,status='old')
      write(*,*) 'Opened input file MATS monitor DV info: ', ifile

      read(*,'(10x,a)') ifile
      open(10,file=ifile,status='old')
      write(*,*) 'Opened input file baseline: ', ifile

      read(*,'(10x,a)') ifile
      open(11,file=ifile,status='old')
      write(*,*) 'Opened input file 2: ', ifile

c-----Allocate variables'

      allocate (base(ndays,nxc,nyc))
      allocate (future(ndays,nxc,nyc))
      allocate (today(ndays))
      allocate (xdayused(ndays))

c-----Read input files

        do n = 1,2  !read over header lines
         read(10,*)
         read(11,*)
        enddo

        do k = 1, ndays
        do r = 1, nxc
        do c = 1, nyc

        read(10,*,end=799)ijcell,blank,lat,long,date,conc

          j = mod(ijcell,1000)
          i = (ijcell-j)/1000

         base(k,i,j)=conc
         today(k) = date

        read(11,*,end=799)ijcell,blank,lat,long,date,conc

          j = mod(ijcell,1000)
          i = (ijcell-j)/1000

         future(k,i,j)=conc

        enddo
        enddo
        enddo

c---- read MATS output file

        do n = 1,1  !read over header lines
         read(20,*)
        enddo

 443    continue

       do k=1,ndays
       xdayused(k)=0
       enddo

        read(20,*,end=444)monid,lat,long,bdv,fdv,bdvmax,fdvmax,ijcell,matsrrf,ppb,days,st,cy

          jj = mod(ijcell,1000)
          ii = (ijcell-jj)/1000

c---estimate rrf

         rcut = usercutoff

 212     continue

         number = 0.
         xbase = 0.
         xfuture = 0.
         maxval = 0.
         nsil = 0
         nthreshold = 0

          do k = 1,ndays

          if(cells.gt.0) then
         
           value = future(k,ii,jj)
           i=ii
           j=jj
           do m = (i - cells), (i + cells) 
            do n = (j - cells), (j + cells)
            if (future(k,m,n).gt.value) then
            value=future(k,m,n)
            i=m
            j=n
            else
            endif
            enddo
           enddo

          else ! stay at cell where monitor is located

          i=ii
          j=jj

          endif

         if(base(k,i,j).gt.rcut) then
          number = number + 1.0
          xbase = xbase + base(k,i,j)
          xfuture = xfuture + future(k,i,j)
          xdayused(k) = 1

c           if(dvtype.eq.'FDV') then
c           contrib = fdv - (fdv*(future(k,i,j)/base(k,i,j)))
           contrib = base(k,i,j) - future(k,i,j)
c           elseif(dvtype.eq.'FDMAX') then
c           contrib = fdvmax - (fdvmax*(future(k,i,j)/base(k,i,j)))
c           else
c           contrib = bdv - (bdv*(future(k,i,j)/base(k,i,j)))
c           endif

           if(contrib.ge.1.0) then
            nsil = nsil + 1
           endif
           if(contrib.ge.0.7) then
            nthreshold = nthreshold + 1
           endif
           if(contrib.ge.maxval) then
            maxval = contrib
           endif

         endif

          enddo !end loop over days

         if(number.lt.usernumber.and.rcut.gt.userthreshold) then
          rcut = rcut - 0.01
          goto 212
         else !criteria satisfied
         endif

c---average contrib calc

         if(number.ge.5) then !calc rrf if 5 days or more
          xbase = xbase/number
          xfuture = xfuture/number
          rrf =  xfuture / xbase
          cut = rcut

          if(dvtype.eq.'FDV') then
          contrib = fdv - (fdv * rrf)
          elseif(dvtype.eq.'FDVMAX') then
          contrib = fdvmax - (fdvmax * rrf)
          else
          contrib = bdv - (bdv * rrf)
          endif

         else !not enough days

          xbase = -9
          xfuture = -9
          rrf =  -9
          cut = rcut

          if(dvtype.eq.'FDV') then
          contrib = -9
          elseif(dvtype.eq.'FDVMAX') then
          contrib = -9
          else
          contrib = -9
          endif

         endif

c---write out aggregate contrib file

        if(bdv.gt.0.and.fdv.gt.0.and.number.ge.5) then
        write(19,304)monid,contrib,lat,long,bdv,bdvmax,fdv,fdvmax,rrf,number,
     &             cut,i,j,maxval,nsil,nthreshold,MESG
        else
        write(29,304)monid,contrib,lat,long,bdv,bdvmax,fdv,fdvmax,rrf,number,
     &             cut,i,j,maxval,nsil,nthreshold,MESG
        endif

c---write out daily contrib file

        if(bdv.gt.0.and.number.ge.5) then

         do k = 1, ndays

           if(cells.gt.0) then

           value = base(k,ii,jj)
           i=ii
           j=jj
            do m = (i - cells), (i + cells)
            do n = (j - cells), (j + cells)
            if (base(k,m,n).gt.value) then
            value=base(k,m,n)
            i=m
            j=n
            else
            endif
            enddo
            enddo

           else ! stay at cell where monitor is located

           i=ii
           j=jj

           endif !cells greater than 0 condition ends
         
         if(bdv.gt.0.and.fdv.gt.0.and.xdayused(k).eq.1) then
         temprrf = future(k,i,j) / base(k,i,j)
         tempcontrib = fdv - (fdv*temprrf)
         abscontrib = base(k,i,j) - future(k,i,j)
         write(9,305)monid,ii,jj,i,j,bdv,bdvmax,fdv,fdvmax,today(k),temprrf,base(k,i,j),future(k,i,j),abscontrib,tempcontrib,MESG
         endif !end base model estimate greater than threshold

        enddo !loop over days

        endif !end monitor output condition

c---misc formatting and loop back over mats file

        goto 443
 444    continue

 304    format(a9,1x,f8.5,1x,f12.6,1x,f12.6,4(1x,f6.2),1x,f8.5,1x,2(1x,f5.1),2(1x,i4),1x,f5.2,2(1x,i3),1x,a25)
 305    format(a9,4(1x,i3),4(1x,f6.2),1x,i8,1x,f8.5,4(1x,f8.4),1x,a25)

c--- end of calculations 

        goto 800
 799     continue
        write(*,*)'Early end of file!'
        write(*,*)ijcell,blank,lat,long,date,conc
 800     continue
      END

      REAL FUNCTION round (X, PLACES)

      REAL X
      INTEGER PLACES
      REAL TEMPX, POWER

      if(places.ge.0) then
       power = 10.0**places
       tempx = x*power
       round = nint(tempx)/power
      else
       round = x
      endif

      RETURN
      END