      program blkavg
c-----------------------------------------------------------------------
c
c            Average a Realization to Coarser Scale
c            **************************************
c
c
c
c
c
c
c
c-----------------------------------------------------------------------
      parameter(MAXXYZ=1000000,EPSLON=0.01,VERSION=1.000)
c
c Dimensioning:
c
      real      var(50),vr(MAXXYZ)
      character datafl*40,outfl*40,str*80
      logical   testfl,inflag
      data      lin/1/,lout/2/
c
c Note VERSION number before anything else:
c
      write(*,9999) VERSION
 9999 format(/' BLKAVG Version: ',f5.3/)
c
c Get the name of the parameter file - try the default name if no input:
c
      write(*,*) 'Which parameter file do you want to use?'
      read (*,'(a40)') str(1:40)
      if(str(1:1).eq.' ')str='blkavg.par                              '
      inquire(file=str(1:40),exist=testfl)
      if(.not.testfl) then
            write(*,*) 'ERROR - the parameter file does not exist,'
            write(*,*) '        check for the file and try again  '
            write(*,*)
            if(str(1:20).eq.'blkavg.par          ') then
                  write(*,*) '        creating a blank parameter file'
                  call makepar
                  write(*,*)
            end if
            stop
      endif
      open(lin,file=str(1:40),status='OLD')
c
c Find Start of Parameters:
c
 1    read(lin,'(a40)',end=98) str(1:40)
      if(str(1:4).ne.'STAR') go to 1
c
c Read Input Parameters:
c
      read(lin,'(a40)',err=98) datafl
      call chknam(datafl,40)
      write(*,*) ' data file = ',datafl

      read(lin,*,err=98) ivr
      write(*,*) ' column for grade = ',ivr

      read(lin,*,err=98) tmin,tmax
      write(*,*) ' trimming limits = ',tmin,tmax

      read(lin,*,err=98) nsim
      write(*,*) ' number of realizations = ',nsim

      read(lin,*,err=98) nx,xmn,xsiz
      write(*,*) ' X grid limits = ',nx,xmn,xsiz

      read(lin,*,err=98) ny,ymn,ysiz
      write(*,*) ' Y grid limits = ',ny,ymn,ysiz

      read(lin,*,err=98) nz,zmn,zsiz
      write(*,*) ' Z grid limits = ',nz,zmn,zsiz

      read(lin,'(a40)',err=98) outfl
      call chknam(outfl,40)
      write(*,*) ' file for output = ',outfl

      read(lin,*,err=98) nxo,xmno,xsizo
      write(*,*) ' X output grid limits = ',nxo,xmno,xsizo

      read(lin,*,err=98) nyo,ymno,ysizo
      write(*,*) ' Y output grid limits = ',nyo,ymno,ysizo

      read(lin,*,err=98) nzo,zmno,zsizo
      write(*,*) ' Z output grid limits = ',nzo,zmno,zsizo

      close(lin)
c
c Make sure the file exists:
c
      inquire(file=datafl,exist=testfl)
      if(.not.testfl) then
            write(*,*) 'ERROR: ',datafl,' does not exist'
            stop
      endif
      if((nx*ny*nz).gt.MAXXYZ) then
            write(*,*) 'ERROR: grid size too large.  have ',MAXXYZ
            stop
      endif
c
c Open the input file:
c
      open(lin,file=datafl, status='OLD')
      read(lin,'(a40)',err=99) str(1:40)
      read(lin,*,err=99)       nvari
      do i=1,nvari
            read(lin,'()',err=99)
      end do
c
c Open the output file:
c
      open(lout,file=outfl, status='UNKNOWN')
      write(lout,100) str(1:40)
 100  format('Block Averaged:',a40,/,'1',/,'value')
c
c Loop over all realizations
c
      do isim=1,nsim
            index = 0
            do iz=1,nz
                  do iy=1,ny
                        do ix=1,nx
                              read(lin,*) (var(i),i=1,nvari)
                              index = index + 1
                              vr(index)  =var(ivr)    
                        end do
                  end do
            end do
c
c Loop over all of the output blocks:
c
            do izo=1,nzo
            do iyo=1,nyo
            do ixo=1,nxo
                  xlo = (xmno + real(ixo-1.5)*xsizo) + EPSLON
                  xhi = (xlo  + xsizo              ) - 2.0*EPSLON
                  call getindx(nx,xmn,xsiz,xlo,ixlo,inflag)
                  call getindx(nx,xmn,xsiz,xhi,ixhi,inflag)
                  ylo = (ymno + real(iyo-1.5)*ysizo) + EPSLON
                  yhi = (ylo  + ysizo              ) - 2.0*EPSLON
                  call getindx(ny,ymn,ysiz,ylo,iylo,inflag)
                  call getindx(ny,ymn,ysiz,yhi,iyhi,inflag)
                  zlo = (zmno + real(izo-1.5)*zsizo) + EPSLON
                  zhi = (zlo  + zsizo              ) - 2.0*EPSLON
                  call getindx(nz,zmn,zsiz,zlo,izlo,inflag)
                  call getindx(nz,zmn,zsiz,zhi,izhi,inflag)

                  navg  = 0
                  vravg = 0
                  do i=ixlo,ixhi
                  do j=iylo,iyhi
                  do k=izlo,izhi
                        index = i + (j-1)*nx + (k-1)*nx*ny
                        if(vr(index).ge.tmin.and.vr(index).lt.tmax) then
                              navg  = navg  + 1
                              vravg = vravg + vr(index)
                        end if
                  end do
                  end do
                  end do

                  if(navg.lt.1) then
                        vravg = -999.0
                  else
                        vravg = vravg / real(navg)
                  end if
                  write(lout,'(f12.5)') vravg
            end do
            end do
            end do
c
c End loop over realizations:
c
      end do
      close(lin)
      close(lout)
c
c Finished:
c
      write(*,9998) VERSION
 9998 format(/' BLKAVG Version: ',f5.3, ' Finished'/)
      stop
 98   stop ' ERROR in parameter file'
 99   stop ' ERROR in data file'
      end



      subroutine chknam(str,len)
c-----------------------------------------------------------------------
c
c                   Check for a Valid File Name
c                   ***************************
c
c This subroutine takes the character string "str" of length "len" and
c removes all leading blanks and blanks out all characters after the
c first blank found in the string (leading blanks are removed first).
c
c
c
c Author: C.V. Deutsch                                Date: January 1993
c-----------------------------------------------------------------------
      parameter (MAXLEN=132)
      character str(MAXLEN)*1
c
c Remove leading blanks:
c
      do i=1,len-1
            if(str(i).ne.' ') then
                  if(i.eq.1) go to 1
                  do j=1,len-i+1
                        k = j + i - 1
                        str(j) = str(k)
                  end do
                  do j=len,len-i+2,-1
                        str(j) = ' '
                  end do
                  go to 1
            end if
      end do
 1    continue
c
c Find first blank and blank out the remaining characters:
c
      do i=1,len-1
            if(str(i).eq.' ') then
                  do j=i+1,len
                        str(j) = ' '
                  end do
                  go to 2
            end if
      end do
 2    continue
c
c Return with modified file name:
c
      return
      end



      subroutine getindx(n,min,siz,loc,index,inflag)
c-----------------------------------------------------------------------
c
c     Gets the coordinate index location of a point within a grid
c     ***********************************************************
c
c
c n       number of "nodes" or "cells" in this coordinate direction
c min     origin at the center of the first cell
c siz     size of the cells
c loc     location of the point being considered
c index   output index within [1,n]
c inflag  true if the location is actually in the grid (false otherwise
c         e.g., if the location is outside then index will be set to
c         nearest boundary
c
c
c
c Author: C.V. Deutsch                               Date: February 1995
c-----------------------------------------------------------------------
      integer   n,index
      real      min,siz,loc
      logical   inflag
c
c Compute the index of "loc":
c
      index = int( (loc-min)/siz + 1.5 )
c
c Check to see if in or out:
c
      if(index.lt.1) then
            index  = 1
            inflag = .false.
      else if(index.gt.n) then
            index  = n
            inflag = .false.
      else
            inflag = .true.
      end if
c
c Return to calling program:
c
      return
      end



      subroutine makepar
c-----------------------------------------------------------------------
c
c                      Write a Parameter File
c                      **********************
c
c
c
c-----------------------------------------------------------------------
      lun = 99
      open(lun,file='blkavg.par')
      write(lun,10)
 10   format('                  Parameters for BLKAVG',/,
     +       '                  *********************',/,/,
     +       'START OF PARAMETERS:')

      write(lun,11)
 11   format('sgsim01.clp               ',
     +       '-file with realization')
      write(lun,12)
 12   format('1                         ',
     +       '-  column for grade')
      write(lun,13)
 13   format('-1.0      1.0e21          ',
     +       '-  trimming limits')
      write(lun,14)
 14   format('1                         ',
     +       '-  number of realizations')
      write(lun,15)
 15   format('160   5202.5    5.0       ',
     +       '-Input size:  nx,xmn,xsiz')
      write(lun,16)
 16   format('180   3602.5    5.0       ',
     +       '-             ny,ymn,ysiz')
      write(lun,17)
 17   format('  1      0.0    1.0       ',
     +       '-             nz,zmn,zsiz')
      write(lun,18)
 18   format('sgsim01.a10               ',
     +       '-file for output')
      write(lun,19)
 19   format(' 80   5205.0   10.0       ',
     +       '-Output size: nx,xmn,xsiz')
      write(lun,20)
 20   format(' 90   3605.0   10.0       ',
     +       '-             ny,ymn,ysiz')
      write(lun,21)
 21   format('  1      0.0    1.0       ',
     +       '-             nz,zmn,zsiz')

      close(lun)
      return
      end
