      program main
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C                                                                      %
C Copyright (C) 1996, The Board of Trustees of the Leland Stanford     %
C Junior University.  All rights reserved.                             %
C                                                                      %
C The programs in GSLIB are distributed in the hope that they will be  %
C useful, but WITHOUT ANY WARRANTY.  No author or distributor accepts  %
C responsibility to anyone for the consequences of using them or for   %
C whether they serve any particular purpose or work at all, unless he  %
C says so in writing.  Everyone is granted permission to copy, modify  %
C and redistribute the programs in GSLIB, but only under the condition %
C that this notice and the above copyright notice remain intact.       %
C                                                                      %
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C                                                                      %
C flagsamp.F - Program to flag a composite or sample if it falls       %
C              inside a block that has a value greater or equal than   %      
C              a threshold                                             %
C                                                                      %
C                                                                      %
C by Julian Ortiz C. August 2001                                       %
C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
      use msflib
      parameter(EPSLON=1.e-12,UNEST=-999.,VERSION=1.000,MAXLEN=512)
      character infile*60,infile2*60,outfl*60,str*60,strlin*512 
      real*8    var(30)
      integer*4 ix,iy,iz,test
      logical   testfl,testind
	data      lin/1/,lin2/2/,lpar/3/, lout/4/
 
     
c
c Dynamic allocation of the arrays.
c
      real, allocatable :: x(:),y(:),z(:),sim(:),simtmp(:)
	+						,in(:)						 
	integer, allocatable :: flag(:),order(:)
c
c Note VERSION number:
c
      write(*,9999) VERSION
 9999 format(/' flagsamp Version: ',f5.3/)
c					  
c Get the name of the parameter file - try the default name if no inputmat:
c
      str(1:1) =' '
      call getarg(1,str)
      if(str(1:1).eq.' ')then
            write(*,*) 'Which parameter file do you want to use?'
            read (*,'(a20)') str(1:20)
      end if
      if(str(1:1).eq.' ')str='flagsamp.par                       '
      inquire(file=str,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.'flagsamp.par         ') then
                  write(*,*) '        creating a blank parameter file'
                  call makepar
                  write(*,*)
            end if
            stop
      endif
 
      open(lpar,file=str,status='OLD')

c
c Find Start of Parameters:
c
 1    read(lpar,'(a4)',end=97) str(1:4)
      if(str(1:4).ne.'STAR') go to 1
c
c Read inputmat Parameters:
c
    
 
      read(lpar,*,err=97) infile
      write(*,*) '  data file to flag = ',infile

	read(lpar,*,err=97) ixl,iyl,izl, ivrl
      write(*,*) ' column numbers = ',ixl,iyl,izl,ivrl
     	     
      read(lpar,*,err=97) infile2,icol
      write(*,*) ' gridded data file and column = ',infile2,icol
     
      read(lpar,*,err=97) nreal
      write(*,*) ' nreal = ',nreal

      read(lpar,*,err=97) nx,xmn,xsiz
      write(*,*) ' nx, xmn, xsiz = ',nx,xmn,xsiz

      read(lpar,*,err=97) ny,ymn,ysiz
      write(*,*) ' ny, ymn, ysiz = ',ny,ymn,ysiz

      read(lpar,*,err=97) nz,zmn,zsiz
      write(*,*) ' nz, zmn, zsiz = ',nz,zmn,zsiz

      read(lpar,*,err=97) outfl
      write(*,*) '  outputfile = ',outfl

      read(lpar,*,err=97) cutoff
      write(*,*) '  threshold = ',cutoff

     
      close(lpar)

c
c Read in the data (if the file exists): 
c
	write(*,*) 'Reading Data...'
      inquire(file=infile,exist=testfl)
      if(.not.testfl) then
            write(*,*) 'ERROR - the data file does not exist,'
            write(*,*) '        check for the file and try again  '
            stop
      endif

c
c From Locmap Code - The first data file exists so open the file and read in the header
c information. Find MAXDAT and allocate the storage that will be used 
c to summarize the data found in the file:
c
 
c
c  Open outfile and write the header
c

      open(lout,file=outfl,status='UNKNOWN')

      open(lin,file=infile,status='OLD')

      read(lin,'(a40)',err=99) str
      write(lout,'(a40)') str
      read(lin,*,err=99) nvari
      write(lout,'(i3)') nvari+1
      do i=1,nvari
            read(lin,'(a40)',err=99) str
            write(lout,'(a40)')      str
      end do
      write(lout,101)
 101  format('Flag')

      maxdat = 0
 213  read(lin,*,end=400,err=97) (var(j),j=1,nvari)
      maxdat = maxdat + 1
      go to 213
 400  continue
      
	write(*,*) 'Allocating Memory...'

      allocate (x(maxdat), stat = test)
            if (test.ne.0) then
                  write(*,*) 'Error: Allocation failed due to ',
     +                 'insufficient memory!', test
                  stop
            end if
c            
c  
      allocate (y(maxdat), stat = test)
            if (test.ne.0) then
                  write(*,*) 'Error: Allocation failed due to ',
     +                 'insufficient memory!', test
                  stop
            end if
c
c
      allocate (z(maxdat), stat = test)
            if (test.ne.0) then
                  write(*,*) 'Error: Allocation failed due to ',
     +                 'insufficient memory!', test
                  stop
            end if
c            
      allocate (sim(maxdat), stat = test)
            if (test.ne.0) then
                  write(*,*) 'Error: Allocation failed due to ',
     +                 'insufficient memory!', test
                  stop
            end if
c            
      allocate (simtmp(maxdat), stat = test)
            if (test.ne.0) then
                  write(*,*) 'Error: Allocation failed due to ',
     +                 'insufficient memory!', test
                  stop
            end if
c            
      allocate (in(maxdat), stat = test)
            if (test.ne.0) then
                  write(*,*) 'Error: Allocation failed due to ',
     +                 'insufficient memory!', test
                  stop
            end if
c            
      allocate (order(maxdat), stat = test)
            if (test.ne.0) then
                  write(*,*) 'Error: Allocation failed due to ',
     +                 'insufficient memory!', test
                  stop
            end if
c            
      allocate (flag(maxdat), stat = test)
            if (test.ne.0) then
                  write(*,*) 'Error: Allocation failed due to ',
     +                 'insufficient memory!', test
                  stop
            end if
c            
      rewind(lin)
      nd = 0
      read(lin,'(a)',err=99)str
      read(lin,*,err=99)     nvari
      do i=1,nvari
            read(lin,'()',err=99)
      end do
212   read(lin,*,end=3) (var(j),j=1,nvari) 
            nd     = nd + 1
		  order(nd)=nd
            x(nd)  = var(ixl)
            y(nd)  = var(iyl)
            z(nd)  = var(izl)
      go to 212
 3    continue

c
c Read in the data (if the file exists): 
c
      write(*,*) 'Opening Gridded File...'
      inquire(file=infile2,exist=testfl)
      if(.not.testfl) then
            write(*,*) 'ERROR - the gridded data file does not exist,'
            write(*,*) '        check for the file and try again  '
            stop
      endif
c
c Open file with simulated values
c
	open(lin2,file=infile2,status='OLD')
      read(lin2,'(a40)',err=99) str
      read(lin2,*,err=99) nvari2
      do i=1,nvari2
            read(lin2,'(a40)',err=99) str
      end do
c
c Find index of data into the simulated grid
c
      write(*,*) 'Finding Indexes for Data...'
      do i=1,maxdat
	      call getindx(nx,xmn,xsiz,x(i),ix,testind)
		  if(testind) then
	      call getindx(ny,ymn,ysiz,y(i),iy,testind)
		  if(testind) then
	      call getindx(nz,zmn,zsiz,z(i),iz,testind)
            if(testind) then
                  in(i)=ix+(iy-1)*nx+(iz-1)*nx*ny
	      else
	            in(i)=-99.0
		  end if
	      else
	            in(i)=-99.0
		  end if 
		  else  
	            in(i)=-99.0
            end if
	end do
c
c Sort data based on their indexes
c
      call sortem(1,maxdat,in,1,order,c,d,e,f,g,h)
c
c Find the realization...
c
      write(*,*) ' Finding realization...'
	do i=1,nreal-1
	      do j=1,nx
	      do k=1,ny
	      do l=1,nz
	            read(lin2,*)
            end do
            end do
            end do
	end do
c
c Read the corresponding node in the simulated file
c
      write(*,*) ' Number of Data ',maxdat
	i=0
2     i=i+1
      write(*,*) 'Working on Data ',i
	if(in(i).lt.0) then
		flag(i)=0
		goto 2
	end if
	if(i.eq.1) then
		delta=in(i)
	else
	if (in(i-1).lt.0) then
            delta=in(i)
	else
		  delta=in(i)-in(i-1)
	end if
	end if
	if (delta.eq.0) then
		flag(i)=flag(i-1)
	else
	      do j=1,delta-1
	            read(lin2,*)
	      end do
            read(lin2,*) (var(j),j=1,nvari2)
		  sim(i)=var(icol)
		  if(sim(i).ge.cutoff) then
				flag(i)=1
			else
				flag(i)=0
		  end if
	end if
      if (i.lt.maxdat) goto 2
c
c Close gridded file
c
	close(lin2)
c
c Sort data to their original order
c
      call sortem(1,maxdat,order,2,flag,in,d,e,f,g,h)

      rewind(lin)
      read(lin,'(a40)',err=99)
      read(lin,*,err=99)
      do i=1,nvari
            read(lin,'(a40)',err=99)
      end do
      do i = 1, maxdat
		  read(lin,'(a)') strlin
		  call strlen(strlin,MAXLEN,lostr)
		  if(flag(i).eq.1) then
			write(lout,'(a,1x,i3)') strlin(1:lostr),flag(i)
		  end if
      end do

	close(lin)
	close(lout)

c
c Finished:
c

      write(*,9998) VERSION
9998  format(/' flagsamp Version: ',f5.3, ' Finished'/)
      stop
 
99	stop 'ERROR in the data file'
97    stop 'ERROR in parameter file!'
      end



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

      write(lun,11)
 11   format('data.dat                                 ',
     +       '- data file to flag')

      write(lun,12)
 12   format('1  2  3                                  ',
     +       '- x, y, z columns')

      write(lun,13)
 13   format('ik3d.out    1                            ',
     +       '- gridded file and column')
     
      write(lun,135)
 135   format('1                                       ',
     +       '- realization number')

      write(lun,14)
 14   format('400 0.5 1.0                              ',
     +       '- nx xmn xsiz')
	  
      write(lun,15)
 15   format('600 0.5 1.0                              ',
     +       '- ny ymn ysiz')
	  
      write(lun,16)
 16   format('200 0.5 1.0                              ',
     +       '- nz, zmn, zsiz')

      write(lun,17)
 17   format('flagsamp.out                             ',
     +       '- outputfile')

      write(lun,18)
 18   format('0.001                                    ',
     +       '- flag if value greater than threshold')
	  
      close(lun)
      return
      end

