program maps
c-----------------------------------------------------------------------
c
c                  Clean a Categorical Realization
c                  *******************************
c
c
c
c-----------------------------------------------------------------------
      parameter (MAXX  = 100, MAXY  = 100, MAXZ  =  50, MAXRC = 10,
     +           MAXXS =  25, MAXYS =  10, MAXZS =   5, VERSION=1.000)
c
c Needed variable declaration:
c
      character  datafl*40,outfl*40,condfl*40,str*80
      real       wt(MAXXS,MAXYS,MAXZS),tmp(100),prop(MAXRC),
     +           tprop(MAXRC),fac(MAXRC),facact(MAXRC)
      integer*4  cat(MAXRC)
      integer*2  val(MAXX,MAXY,MAXZ),smo(MAXX,MAXY,MAXZ)
      logical    cond(MAXX,MAXY,MAXZ),testfl
      data       lin/1/,lout/2/
c
c Note VERSION number:
c
      write(*,9999) VERSION
 9999 format(/' MAPS 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 (*,'(a20)') str(1:20)
      if(str(1:1).eq.' ') str(1:20) = 'maps.par            '
      inquire(file=str(1:20),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.'maps.par            ') then
                  write(*,*) '        creating a blank parameter file'
                  call makepar
                  write(*,*)
            end if
            stop
      endif
      open(lin,file=str(1:20),status='old')
c
c Find Start of Parameters:
c
 1    read(lin,'(a4)',end=98) str(1:4)
      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,'(a40)',err=98) outfl
      call chknam(outfl,40)
      write(*,*) ' Output file = ',outfl

      read(lin,*,err=98) nx,xmn,xsiz
      write(*,*) ' X grid specification = ',nx,xmn,xsiz
      if(nx.gt.MAXX) stop 'nx too big'

      read(lin,*,err=98) ny,ymn,ysiz
      write(*,*) ' Y grid specification = ',ny,ymn,ysiz
      if(ny.gt.MAXY) stop 'ny too big'

      read(lin,*,err=98) nz,zmn,zsiz
      write(*,*) ' Z grid specification = ',nz,zmn,zsiz
      if(nz.gt.MAXZ) stop 'nz too big'

      read(lin,*,err=98) isim
      write(*,*) ' Realization number = ',isim

      read(lin,*,err=98) ncat
      write(*,*) ' number of categories = ',ncat
      if(ncat.gt.MAXRC) stop 'ncat too big'

      read(lin,*,err=98) (cat(i),i=1,ncat)
      write(*,*) ' categories = ',(cat(i),i=1,ncat)

      read(lin,*,err=98) (tprop(i),i=1,ncat)
      write(*,*) ' target proportions = ',(tprop(i),i=1,ncat)

      read(lin,*,err=98) (fac(i),i=1,ncat)
      write(*,*) ' factors = ',(fac(i),i=1,ncat)

      read(lin,'(a40)',err=98) condfl
      call chknam(condfl,40)
      write(*,*) ' Conditioning data file = ',condfl

      read(lin,*,err=98) icx,icy,icz,iccat
      write(*,*) ' The columns = ',icx,icy,icz,iccat

      read(lin,*,err=98) cfac
      write(*,*) ' Conditioning data factor = ',cfac

      read(lin,*,err=98) niter
      write(*,*) ' Number of iterations = ',niter

      read(lin,*,err=98) nxs,nys,nzs
      write(*,*) ' nxs,nys,nzs = ',nxs,nys,nzs
      if(nxs.gt.MAXXS) stop 'nxs too big'
      if(nys.gt.MAXYS) stop 'nys too big'
      if(nzs.gt.MAXZS) stop 'nzs too big'
      nhxs = int(nxs/2.0)
      nhys = int(nys/2.0)
      nhzs = int(nzs/2.0)
      if(2*nhxs.eq.nxs) stop 'nxs must be odd (not even)'
      if(2*nhys.eq.nys) stop 'nys must be odd (not even)'
      if(2*nhzs.eq.nzs) stop 'nzs must be odd (not even)'

      do iz=1,nzs
            do iy=1,nys
                  read(lin,*,err=98)      (wt(ix,iy,iz),ix=1,nxs)
                  write(*,*) ' weights: ',(wt(ix,iy,iz),ix=1,nxs)
            end do
      end do

      close(lin)
c
c Initialize the conditioning data flags:
c
      do iz=1,nz
            do iy=1,ny
                  do ix=1,nx
                        cond(ix,iy,iz) = .false.
                  end do
            end do
      end do
      inquire(file=condfl,exist=testfl)
      if(testfl) then
            open(lin,file=condfl,status='OLD')
            read(lin,'()',err=99)
            read(lin,*,err=99) nvari
            do i=1,nvari
                  read(lin,'()',err=99)
            end do
            ix = 1
            iy = 1
            iz = 1
 2          read(lin,*,err=99,end=3) (tmp(i),i=1,nvari)
            if(icx.gt.0) call getindx(nx,xmn,xsiz,tmp(icx),ix,testfl)
            if(icy.gt.0) call getindx(ny,ymn,ysiz,tmp(icy),iy,testfl)
            if(icz.gt.0) call getindx(nz,zmn,zsiz,tmp(icz),iz,testfl)
            cond(ix,iy,iz) = .true.
            go to 2
 3          continue
            close(lin)
      endif
c
c Make sure the data file exists:
c
      inquire(file=datafl,exist=testfl)
      if(.not.testfl) then
            write(*,*) 'ERROR - the data file does not exist,'
            write(*,*) '        check for the file and try again  '
            stop
      endif
      open(lin,file=datafl,status='OLD')
c
c Read in the data:
c
      read(lin,'(a20)',err=99) str(1:20)
      read(lin,*,err=99) nvari
      do i=1,nvari
            read(lin,'()',err=99)
      end do
c
c Find the correct realization:
c
      do i=1,isim-1
            do j=1,nx*ny*nz
                  read(lin,*)
            end do
      end do
c
c Read this realization:
c
      ivr = 1
      do iz=1,nz
            do iy=1,ny
                  do ix=1,nx
                        read(lin,*,err=99) (tmp(i),i=1,nvari)
                        kat = int(tmp(ivr)+0.5)
                        iii = 1
                        do ic=1,ncat
                              if(kat.eq.cat(ic)) iii = ic
                        end do
                        smo(ix,iy,iz) = iii
                  end do
            end do
      end do
      close(lin)
c
c Prepare the output file:
c
      open(lout,file=outfl)
      write(lout,998) str(1:20)
 998  format('Cleaned: ',a20,/,'1',/,'category')
c
c MAIN Loop over all of the iterations:
c
      do iter=1,niter
c
c Report the initial / target proportions:
c
      write(*,*)
      do ic=1,ncat
            prop(ic) = 0.0 
      end do
      do iz=1,nz
            do iy=1,ny
                  do ix=1,nx
                        iii = smo(ix,iy,iz)
                        prop(iii) = prop(iii) + 1.0
                  end do
            end do
      end do
      do ic=1,ncat
            prop(ic) = prop(ic) / real(nx*ny*nz)
            write(*,200) ic,prop(ic),tprop(ic)
 200        format(' Category ',i3,' realization = ',f6.4,
     +             ' target = ',f6.4)
            if(prop(ic).le.0.001) prop(ic) = 0.001
            facact(ic) = fac(ic) * ( tprop(ic) / prop(ic) )
      end do
c
c Load "val" array from "smo" array:
c
      do iz=1,nz
            do iy=1,ny
                  do ix=1,nx
                        val(ix,iy,iz) = smo(ix,iy,iz)
                  end do
            end do
      end do
c
c Weighted value:
c
      do iz=1,nz
      do iy=1,ny
      do ix=1,nx
c
c Working on ix,iy,iz:
c
            do ic=1,ncat
                  prop(ic) = 0.0 
            end do
            k = 0
            do jz=iz-nhzs,iz+nhzs
              k = k + 1
              j = 0
              do jy=iy-nhys,iy+nhys
                j = j + 1
                i = 0
                do jx=ix-nhxs,ix+nhxs
                  i = i + 1
                  if(jx.ge.1.and.jx.le.nx.and.
     +               jy.ge.1.and.jy.le.ny.and.
     +               jz.ge.1.and.jz.le.nz) then
                        icode = val(jx,jy,jz)
                        if(cond(jx,jy,jz)) then
                              prop(icode) = prop(icode) + wt(i,j,k)*cfac
                        else
                              prop(icode) = prop(icode) + wt(i,j,k)
                        end if
                  end if
            end do
            end do
            end do
            if(.not.cond(ix,iy,iz)) then
                  pmax = prop(1) * facact(1)
                  smo(ix,iy,iz) = 1
                  do ic=2,ncat
                        prop(ic) = prop(ic) * facact(ic)
                        if(prop(ic).gt.pmax) then
                              pmax = prop(ic)
                              smo(ix,iy,iz) = ic
                        end if
                  end do
            end if
      end do
      end do
      end do
c
c Finish main loop:
c
      end do
c
c Write out the results:
c
      do ic=1,ncat
            prop(ic) = 0.0 
      end do
      do iz=1,nz
            do iy=1,ny
                  do ix=1,nx
                        iii = smo(ix,iy,iz)
                        prop(iii) = prop(iii) + 1.0
                        kat = cat(iii)
                        write(lout,'(i2)') kat
                  end do
            end do
      end do
      close(lout)
c
c Compute number in "smo"
c
      write(*,*)
      do ic=1,ncat
            prop(ic) = prop(ic) / real(nx*ny*nz)
            write(*,200) ic,prop(ic),tprop(ic)
            if(prop(ic).le.0.001) prop(ic) = 0.001
            facact(ic) = fac(ic) * tprop(ic) / prop(ic)
      end do
c
c Finished:
c
      write(*,9998) VERSION
 9998 format(/' MAPS Version: ',f5.3, ' Finished'/)
      stop
 98   stop 'ERROR in parameters somewhere'
 99   stop 'ERROR in data somewhere'
      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-----------------------------------------------------------------------
      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 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-----------------------------------------------------------------------
      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 makepar
c-----------------------------------------------------------------------
c
c                      Write a Parameter File
c                      **********************
c
c
c
c-----------------------------------------------------------------------
      lun = 99
      open(lun,file='maps.par',status='UNKNOWN')
      write(lun,10)
 10   format('                 Parameters for MAPS',/,
     +       '                 *******************',/,/,
     +       'START OF PARAMETERS:')

      write(lun,11)
 11   format('sisim.out                    ',
     +       '-file with input initial realization')
      write(lun,12)
 12   format('maps.out                     ',
     +       '-file for output cleaned realization')
      write(lun,13)
 13   format('100  0.5  1.0                ',
     +       '-nx, xmn, xsiz')
      write(lun,14)
 14   format('100  0.5  1.0                ',
     +       '-ny, ymn, ysiz')
      write(lun,15)
 15   format('  1  0.5  1.0                ',
     +       '-nz, zmn, zsiz')
      write(lun,16)
 16   format('1                            ',
     +       '-realization number')
      write(lun,17)
 17   format('4                            ',
     +       '-number of categories')
      write(lun,18)
 18   format('1    2    3    4             ',
     +       '-   categories')
      write(lun,19)
 19   format('0.05 0.20 0.30 0.45          ',
     +       '-   target global proportions')
      write(lun,20)
 20   format('1.0  1.0  1.0  1.0           ',
     +       '-   factors for target proportions')
      write(lun,21)
 21   format('data.dat                     ',
     +       '-file with conditioning data')
      write(lun,22)
 22   format('1 2 3 4                      ',
     +       '-   columns for X, Y, Z, and category')
      write(lun,23)
 23   format('4.0                          ',
     +       '-   C = conditioning data weight')
      write(lun,24)
 24   format('1                            ',
     +       '-number of iterations')
      write(lun,25)
 25   format('5  5  1                      ',
     +       '-size of smoothing window W')
      write(lun,26)
 26   format(' 1.0 1.0 1.0 1.0 1.0         ',
     +       '-   weights iy=u+2')
      write(lun,27)
 27   format(' 1.0 2.0 3.0 2.0 1.0         ',
     +       '-   weights iy=u+1')
      write(lun,28)
 28   format(' 1.0 3.0 5.0 3.0 1.0         ',
     +       '-   weights iy=u')
      write(lun,29)
 29   format(' 1.0 2.0 3.0 2.0 1.0         ',
     +       '-   weights iy=u-1')
      write(lun,30)
 30   format(' 1.0 1.0 1.0 1.0 1.0         ',
     +       '-   weights iy=u-2')
	
      close(lun)
      return
      end