module gam2_sub
use sep
!
!                Variogram of Data on a 2-D Regular Grid
!                ***************************************
!
! The following Parameters control static dimensioning within gam2:
!
!   MAXX      maximum number of nodes in the X direction
!   MAXY      maximum number of nodes in the Y direction
!   MAXVAR    maximum number of variables
!
!   MAXDIR    maximum number of directions possible at one time
!   MAXLAG    maximum number of lags at one time
!   MXVARG    maximum number of variograms possible at one time
!
!   MXDLV     MAXDIR*MAXLAG*MXVARG - used for dimensioning
!   EPSLON    a small number to avoid dividing by zero
!
!-----------------------------------------------------------------------
!
! User Adjustable:
!
integer, parameter:: MAXX=500,MAXY=500,MAXVAR=3,MAXDIR=4,MAXLAG=200,MXVARG=8 
integer, parameter:: maxlg=maxlag+2,mxdlv=maxdir*maxlg*mxvarg
real, parameter :: EPSLON=0.00001, VERSION=1.202

contains

      subroutine gam2(nlag,nx,ny,ndir,ixd,iyd,vr,tmin,tmax,nvarg,ivtail,&
                     ivhead,ivtype,np,gam,hm,tm,hv,tv)
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!                                                                      %
! Copyright (C) 1992 Stanford Center for Reservoir Forecasting.  All   %
! rights reserved.  Distributed with: C.V. Deutsch and A.G. Journel.   %
! ``GSLIB: Geostatistical Software Library and User's Guide,'' Oxford  %
! University Press, New York, 1992.                                    %
!                                                                      %
! The programs in GSLIB are distributed in the hope that they will be  %
! useful, but WITHOUT ANY WARRANTY.  No author or distributor accepts  %
! responsibility to anyone for the consequences of using them or for   %
! whether they serve any particular purpose or work at all, unless he  %
! says so in writing.  Everyone is granted permission to copy, modify  %
! and redistribute the programs in GSLIB, but only under the condition %
! that this notice and the above copyright notice remain intact.       %
!                                                                      %
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!-----------------------------------------------------------------------
!
!                Variogram of Data on a 2-D Regular Grid
!                ***************************************
!
! This subroutine computes any of eight different measures of spatial
! continuity for regular spaced 2-D data.  Missing values are allowed
! and the grid need not be square.
!
!
!
! INPUT VARIABLES:
!
!   nlag             Number of lags to calculate
!   nx               Number of units in x
!   ny               Number of units in y
!   ndir             Number of directions to consider
!   ixd(ndir)        X indicator of direction - number of X grid columns
!                      that must be shifted to move from a node on the
!                      grid to the next nearest node on the grid which
!                      lies on the directional vector.
!   iyd(ndir)        Y indicator of direction - similar to ixd, number
!                      of grid lines that must be shifted to nearest
!                      node which lies on the directional vector
!   nv               The number of variables
!   vr(nx,ny,nv)     Three dimensional array of data
!   tmin,tmax        Trimming limits
!   nvarg            Number of variograms to compute
!   ivtail(nvarg)    Variable for the tail of each variogram
!   ivhead(nvarg)    Variable for the head of each variogram
!   ivtype(nvarg)    Type of variogram to compute:
!                      1. semivariogram
!                      2. cross-semivariogram
!                      3. covariance
!                      4. correlogram
!                      5. general relative semivariogram
!                      6. pairwise relative semivariogram
!                      7. semivariogram of logarithms
!                      8. rodogram
!                      9. madogram
!                     10. indicator semivariogram: an indicator variable
!                         is constructed in the main program.
!
!
!
! OUTPUT VARIABLES:  The following arrays are ordered by direction,
!                    then variogram, and finally lag, i.e.,
!                      iloc = (id-1)*nvarg*nlag+(iv-1)*nlag+il
!
!   np()             Number of pairs
!   gam()            semivariogram, covariance, correlogram,... value
!   hm()             Mean of the tail data
!   tm()             Mean of the head data
!   hv()             Variance of the tail data
!   tv()             Variance of the head data
!
!
!
! Original:  A.G. Journel                                           1978
! Revisions: B.E. Buxton                                       Apr. 1982
!            C.V. Deutsch                                      July 1990
!-----------------------------------------------------------------------
!      include  'gam2.inc'
      real      vr(MAXX,MAXY,MAXVAR)
      real*8    gam(MXDLV),hm(MXDLV),tm(MXDLV),hv(MXDLV),tv(MXDLV),&
                np(MXDLV)
      integer   ixd(MAXDIR),iyd(MAXDIR),ivtail(MXVARG),ivhead(MXVARG),&
                ivtype(MXVARG)
!
! Initialize the summation arrays for each direction, variogram, and lag
!
      nsiz = ndir*nvarg*nlag
      if(nsiz.gt.MXDLV) then
            write(0,*) 'ERROR: available storage in gam2 = ',MXDLV
            write(0,*) '       requested storage         = ',nsiz
            stop
      endif
      do 1 i=1,nsiz
            np(i)  = 0.
            gam(i) = 0.0
            hm(i)  = 0.0
            tm(i)  = 0.0
            hv(i)  = 0.0
            tv(i)  = 0.0
 1    continue
!
! First fix the location of a seed point on the grid (ix,iy):
!
      do 2 ix=1,nx
      do 2 iy=1,ny
!
! For the fixed seed point, loop through all directions:
!
            do 3 id=1,ndir
              ixinc = ixd(id)
              iyinc = iyd(id)
              ix1   = ix
              iy1   = iy
!
! For this direction, loop over all the lags:
!
              do 4 il=1,nlag
!
! Check to be sure that the point being considered is still in the
! grid - if not, then finished with this direction:
!
                ix1 = ix1 + ixinc
                if(ix1.lt.1.or.ix1.gt.nx) go to 3
                iy1 = iy1 + iyinc
                if(iy1.lt.1.or.iy1.gt.ny) go to 3
!
! For this direction and lag, loop over all variograms:
!
                do 5 iv=1,nvarg
                  it = ivtype(iv)
!
! Get the head value, skip this value if missing:
!
                  i     = ivhead(iv)
                  vrt   = vr(ix,iy,i)
                  if(vrt.lt.tmin.or.vrt.ge.tmax) go to 5
!
! Get the tail value, skip this value if missing:
!
                  i     = ivtail(iv)
                  vrh   = vr(ix1,iy1,i)
                  if(vrh.lt.tmin.or.vrh.ge.tmax) go to 5
!
! Need increment for the cross semivariogram only:
!
                  if(it.eq.2) then
                        i     = ivhead(iv)
                        vrtpr = vr(ix1,iy1,i)
                        if(vrtpr.lt.tmin.or.vrtpr.ge.tmax) go to 5
                        i     = ivtail(iv)
                        vrhpr = vr(ix,iy,i)
                        if(vrhpr.lt.tmin.or.vrhpr.ge.tmax) go to 5
                  endif
!
! We have an acceptable pair, therefore accumulate all the statistics
! that are required for the variogram:
!
                  i      = (id-1)*nvarg*nlag+(iv-1)*nlag+il
                  np(i)  = np(i) + 1.
                  tm(i)  = tm(i) + dble(vrt)
                  hm(i)  = hm(i) + dble(vrh)
!
! Choose the correct variogram type and keep relevant sums:
!
                  if(it.le.1.or.it.ge.10) then
                    gam(i) = gam(i) + dble((vrh-vrt)*(vrh-vrt))
                  else if(it.eq.2) then
                    gam(i) = gam(i) + dble((vrhpr-vrh)*(vrt-vrtpr))
                  else if(it.eq.3) then
                    gam(i) = gam(i) +  dble(vrh*vrt)
                  else if(it.eq.4) then
                    gam(i) = gam(i) +  dble(vrh*vrt)
                    hv(i)  = hv(i)  +  dble(vrh*vrh)
                    tv(i)  = tv(i)  +  dble(vrt*vrt)
                  else if(it.eq.5) then
                    gam(i) = gam(i) + dble((vrh-vrt)*(vrh-vrt))
                  else if(it.eq.6) then
                    if((vrt+vrh).lt.EPSLON) then
                        np(i)  = np(i) - 1.
                        tm(i)  = tm(i) - vrt
                        hm(i)  = hm(i) - vrh
                    else
                        gamma  = 2.0*(vrt-vrh)/(vrt+vrh)
                        gam(i) = gam(i) + dble(gamma*gamma)
                    endif
                  else if(it.eq.7) then
                    if(vrt.lt.EPSLON.or.vrh.lt.EPSLON) then
                        np(i)  = np(i) - 1.
                        tm(i)  = tm(i) - dble(vrt)
                        hm(i)  = hm(i) - dble(vrh)
                    else
                        gamma  = alog(vrt)-alog(vrh)
                        gam(i) = gam(i) + dble(gamma*gamma)
                    endif
                  else if(it.eq.8) then
                    gam(i) = gam(i) + dble(sqrt(abs(vrt-vrh)))
                  else if(it.eq.9) then
                    gam(i) = gam(i) + dble(abs(vrt-vrh))
                  endif
 5              continue
 4            continue
 3          continue
 2    continue
!
! Get average values for gam, hm, tm, hv, and tv, then compute
! the correct "variogram" measure:
!
      do 6 id=1,ndir
      do 6 iv=1,nvarg
      do 6 il=1,nlag
            i = (id-1)*nvarg*nlag+(iv-1)*nlag+il
            if(np(i).eq.0.) go to 6
            rnum   = np(i)
            gam(i) = gam(i) / dble(rnum)
            hm(i)  = hm(i)  / dble(rnum)
            tm(i)  = tm(i)  / dble(rnum)
            hv(i)  = hv(i)  / dble(rnum)
            tv(i)  = tv(i)  / dble(rnum)
!
!  1. report the semivariogram rather than variogram
!  2. report the cross-semivariogram rather than variogram
!  3. the covariance requires "centering"
!  4. the correlogram requires centering and normalizing
!  5. general relative requires division by lag mean
!  6. report the semi(pairwise relative variogram)
!  7. report the semi(log variogram)
!  8. report the semi(rodogram)
!  9. report the semi(madogram)
!
            it = ivtype(iv)
            if(it.le.2) then
                  gam(i) = 0.5 * gam(i)
            else if(it.eq.3) then
                  gam(i) = gam(i) - hm(i)*tm(i)
            else if(it.eq.4) then
                  hv(i)  = sqrt(hv(i)-hm(i)*hm(i))
                  tv(i)  = sqrt(tv(i)-tm(i)*tm(i))
                  if((hv(i)*tv(i)).lt.EPSLON) then
                        gam(i) = 0.0
                  else
                        gam(i) =(gam(i)-hm(i)*tm(i))/(hv(i)*tv(i))
                  endif
!
! Square "hv" and "tv" so that we return the variance:
!
                  hv(i)  = hv(i)*hv(i)
                  tv(i)  = tv(i)*tv(i)
            else if(it.eq.5) then
                  htave  = 0.5*real(hm(i)+tm(i))
                  htave  = htave   *   htave
                  if(htave.lt.EPSLON) then
                        gam(i) = 0.0
                  else
                        gam(i) = gam(i)/dble(htave)
                  endif
            else if(it.ge.6) then
                  gam(i) = 0.5 * gam(i)
            endif
 6    continue
      return
  end subroutine

     subroutine readparm(nlag,nx,ny,xsiz,ysiz,ndir,ixd,iyd,vr,tmin,&
                          tmax,nvarg,ivtail,ivhead,ivtype,outfl,names)
!-----------------------------------------------------------------------
!
!                  Initialization and Read Parameters
!                  **********************************
!
! The input parameters and data are read in from their files. Some quick
! error checking is performed and the statistics of all the variables
! being considered are written to standard output.
!
!
!
! Original: C.V. Deutsch                                 Date: July 1990
!-----------------------------------------------------------------------
!      include  'gam2.inc'
      parameter(MV=20)
      real      vr(MAXX,MAXY,MAXVAR),var(MV),cut(MXVARG)
      real*8    avg(MV),ssq(MV)
      real      vrmin(MAXVAR),vrmax(MAXVAR)
      integer   ixd(MAXDIR),iyd(MAXDIR),ivtail(MXVARG),ivhead(MXVARG),&
               ivtype(MXVARG),ivar(MV),num(MV),ivc(MXVARG)
      character datafl*40,outfl*40,str*40,names(MAXVAR)*12
      logical   testfl
      data      lin/1/,ncut/0/
			integer,external :: getch
			integer :: ierr
!
! Note VERSION number:
!
!      write(0,9999) VERSION
 9999 format(/' GAM2M Version: ',f5.3/)
!
! Get the name of the parameter file - try the default name if no input:
!
!      write(0,*) 'Which parameter file do you want to use?'
!      read (0,'(a40)') str
!      if(str(1:1).eq.' ')str='gam2.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  '
!            stop
!      endif
!      open(lin,file=str,status='OLD')
!
! Find Start of Parameters:
!
! 1    read(lin,'(a4)',end=98) str(1:4)
!      if(str(1:4).ne.'STAR') go to 1
!
! Read Input Parameters:
!
!      read(lin,'(a40)',err=98) datafl
!      write(*,*) ' data file:           ',datafl
!      read(lin,*,err=98)       nvar
	
		
			call from_history("n3",nvar,1)
      write(0,*) ' number of variables: ',nvar
!      backspace lin
			 
			ierr=getch("vars","d",ivar)	 
!      read(lin,*,err=98)       j,(ivar(i),i=1,nvar)
      write(0,*) ' columns:             ',(ivar(i),i=1,nvar)
			call from_param("tmin",tmin,-1.e21)
			call from_param("tmax",tmax,1.e21)
!      read(lin,*,err=98)       tmin,tmax
      write(0,*) ' trimming limits:     ',tmin,tmax
!      read(lin,'(a40)',err=98) outfl
!      write(*,*) ' output file:         ',outfl
			call from_history("n1",nx)
			call from_history("n2",ny)
!      read(lin,*,err=98)       nx,ny,isim
			isim=1
      write(0,*) ' nx,ny,isim:          ',nx,ny,isim
			call from_history("d1",xsiz)
			call from_history("d2",ysiz)
!      read(lin,*,err=98)       xsiz,ysiz
      write(0,*) ' xsiz,ysiz:           ',xsiz,ysiz
			call from_param("ndir",ndir)
			call from_param("nlag",nlag)
!      read(lin,*,err=98)       ndir,nlag
      write(0,*) ' ndir,nlag:           ',ndir,nlag
			if(ndir.ne.getch("xshift","d",ixd))&
				call erexit("xshift must have ndir elements")
			if(ndir.ne.getch("yshift","d",iyd))&
				call erexit("yshift must have ndir elements")
      do 2 i=1,ndir
!            read(lin,*,err=98)       ixd(i),iyd(i)
            write(0,*) ' direction:           ',ixd(i),iyd(i)
 2    continue
	 nvarg=getch("itail","d",ivtail)
	 if(getch("ihead","d",ivhead).ne.nvarg) &
		call erexit("number of head elemenst not equal number of tail elements")
	 if(getch("itype","d",ivtype).ne.nvarg) &
		call erexit("number of head elemenst not equal number of tail elements")
!   read(lin,*,err=98)       nvarg
      write(0,*) ' number of variograms:',nvarg
      do 3 i=1,nvarg
!            read(lin,*,err=98) ivtail(i),ivhead(i),ivtype(i)
            write(0,*) ' tail,head,type:    ',&
                        ivtail(i),ivhead(i),ivtype(i)
            if(ivtype(i).eq.10) then
                   ncut = ncut + 1
                   if((nvar+ncut).gt.MAXVAR) then
                         write(0,*) ' Too many indicator cutoffs!'
                         write(0,*) ' use fewer or increase MAXVAR'
                         stop
                   endif
                   backspace lin
                   read(lin,*,err=98) ii,jj,kk,cut(ncut)
                   ivc(ncut) = ivtail(i)
                   ivtail(i) = nvar + ncut
                   ivhead(i) = nvar + ncut
                   write(names(nvar+ncut),140) ncut
 140               format('Indicator ',i2)
                   write(0,*) ' indicator cutoff:  ',cut(ncut)
            endif
 3    continue
!      close(lin)
!
! Perform some quick error checking:
!
      if(nx.lt.1)         stop 'nx is too small: check parameters'
      if(nx.gt.MAXX)      stop 'nx is too big: check gam2.inc'
      if(ny.lt.1)         stop 'ny is too small: check parameters'
      if(ny.gt.MAXY)      stop 'ny is too big: check gam2.inc'
      if(ndir.lt.1)       stop 'ndir is too small: check parameters'
      if(ndir.gt.MAXDIR)  stop 'ndir is too big: check gam2.inc'
      if(nvar.lt.1)       stop 'nvar is too small: check parameters'
      if(nvar.gt.MAXVAR)  stop 'nvar is too big: check gam2.inc'
    if(nlag.lt.1)       stop 'nlag is too small: check parameters'
      if(nlag.gt.MAXLAG)  stop 'nlag is too big: check gam2.inc'
      if(nvarg.lt.1)      stop 'nvarg is too small: check parameters'
      if(nvarg.gt.MXVARG) stop 'nvarg is too big: check gam2.inc'
      if(ncut.ge.1) then
            if(tmin.gt.0.0) stop 'tmin interferes with indicators!'
            if(tmax.le.1.0) stop 'tmax interferes with indicators!'
      endif
!
! Check to make sure the data file exists, then either read in the
! data or write an error message and stop:
!
!      inquire(file=datafl,exist=testfl)
!      if(testfl) then
!
! The data file exists so open the file and read in the header
! information. Initialize the storage that will be used to summarize
! the data found in the file:
!
!            open(lin,file=datafl,status='OLD')
!            read(lin,'(a40)',err=99) str
!            read(lin,*,err=99)       nvari
names="none"
            do 6 i=1,nvar
!                  read(lin,'(a40)',err=99) str
                  do 60 iv=1,nvar
                        j=ivar(iv)
                        if(i.eq.j) names(iv) = str(1:12)
 60               continue
                  num(i) = 0
                  avg(i) = 0.0
                  ssq(i) = 0.0
 6          continue
!
! Read the regualar grid information row wise (x cycles fastest):


call sreed("in",vr(1:nx,1:ny,1:nvar),4*nx*ny*nvar)
!
                  do 7 iy=1,ny
                  do 7 ix=1,nx
!                              read(lin,*,err=99) (var(i),i=1,nvari)
                              do 8 iv=1,nvar
!                                    i=ivar(iv)
!                                    vr(ix,iy,iv) = var(i)
                                    if(vr(ix,iy,iv).lt.tmin.or.vr(ix,iy,iv).ge.tmax)&
                                      go to 8
                                    num(iv) = num(iv) + 1
                                    avg(iv) = avg(iv) + dble(var(i))
                                    ssq(iv) = ssq(iv) + &
                                              dble(var(i)*var(i))
 8                            continue
 7                continue
!            close(lin)
!
! Compute the averages and variances as an error check for the user:
!
            do 9 iv=1,nvar
                  if(num(iv).gt.0) then
                  avg(iv) = avg(iv) / dble(num(iv))
                  ssq(iv) =(ssq(iv) / dble(num(iv)))&
                         - avg(iv) * avg(iv)
                  write(0,*) 'Variable number ',iv
                  write(0,*) '  Number   = ',num(iv)
                  write(0,*) '  Average  = ',real(avg(iv))
                  write(0,*) '  Variance = ',real(ssq(iv))
                  endif
 9          continue
!      else
!            write(0,*) 'ERROR data file ',datafl,' does not exist!'
!            stop
!      endif
!
! Construct Indicator Variables if necessary:
!
      do 10 ic=1,ncut
            iv = ivc(ic)
            jv = nvar + ic
            do 11 ix=1,nx
            do 11 iy=1,ny
                  if(vr(ix,iy,iv).lt.tmin.or.vr(ix,iy,iv).ge.tmax) then
                        vr(ix,iy,jv) = tmin - EPSLON
                  else if(vr(ix,iy,iv).lt.cut(ic)) then
                        vr(ix,iy,jv) = 0.0
                  else
                        vr(ix,iy,jv) = 1.0
                  endif
 11         continue
 10   continue
!
! Establish minimums and maximums:
!
      do 50 i=1,MAXVAR
            vrmin(i) =  1.0e21
            vrmax(i) = -1.0e21
 50   continue
      do 51 ix=1,nx
      do 51 iy=1,ny
      do 51 iv=1,nvar+ncut
            if(vr(ix,iy,iv).lt.tmin.or.&
               vr(ix,iy,iv).ge.tmax)      go to 51
            if(vr(ix,iy,iv).lt.vrmin(iv)) vrmin(iv) = vr(ix,iy,iv)
            if(vr(ix,iy,iv).gt.vrmax(iv)) vrmax(iv) = vr(ix,iy,iv)
 51   continue
!
! Check on the variogams that were requested:
!
      call check(nvarg,ivtail,ivhead,ivtype,vrmin,vrmax,names,MAXVAR)
!
! Return:
!
      return
!
! Error in an Input File Somewhere:
!
 98   stop 'ERROR in parameter file!'
 99   stop 'ERROR in data file!'
end subroutine



      subroutine writeout(nlag,ndir,ixd,iyd,xsiz,ysiz,nvarg,ivtail,&
                         ivhead,ivtype,np,gam,hm,tm,hv,tv,outfl,names)
!-----------------------------------------------------------------------
!
!                  Write Out the Results of GAM2
!                  *****************************
!
! An output file will be written which contains each directional
! variogram ordered by direction and then variogram (the directions
! cycle fastest then the variogram number).  For each variogram there
! will be a one line description and then "nlag" lines with:
!
!        a) lag number (increasing from 1 to nlag)
!        b) separation distance
!        c) the "variogram" value
!        d) the number of pairs for the lag
!        e) the mean of the data contributing to the tail
!        f) the mean of the data contributing to the head
!        g) IF the correlogram - variance of tail values
!        h) IF the correlogram - variance of head values
!
!
!
!
!
! Original: C.V. Deutsch                                 Date: July 1990
!-----------------------------------------------------------------------
!      include  'gam2.inc'
      real*8    gam(MXDLV),hm(MXDLV),tm(MXDLV),hv(MXDLV),tv(MXDLV),&
                np(MXDLV)
      integer   ixd(MAXDIR),iyd(MAXDIR),ivtail(MXVARG),ivhead(MXVARG),&
                ivtype(MXVARG)
      character outfl*40,title*80,names(MAXVAR)*12
      data      lout/1/
!
! Loop over all the variograms that have been computed:
!
!      open(lout,file=outfl,status='NEW',carriagecontrol='LIST')
	call to_history("n3",nvarg)
	call to_history("n1",nlag)
  call to_history("n2",ndir)
	call to_history("o1",1.)
	call to_history("d1",1.)
	call to_history("esize",8)
		
!	outfl="dff.var"
!      open(lout,file=outfl,status='NEW')
      do 1 iv=1,nvarg
!
! Construct a title that reflects the variogram type and the variables
! that were used to calculate the variogram:
!
      it = ivtype(iv)
      if(it.eq. 1) title(1:24) = 'Semivariogram           '
      if(it.eq. 2) title(1:24) = 'Cross Semivariogram     '
      if(it.eq. 3) title(1:24) = 'Covariance              '
      if(it.eq. 4) title(1:24) = 'Correlogram             '
      if(it.eq. 5) title(1:24) = 'General Relative        '
      if(it.eq. 6) title(1:24) = 'Pairwise Relative       '
      if(it.eq. 7) title(1:24) = 'Variogram of Logarithms '
      if(it.eq. 8) title(1:24) = 'Rodogram (square root)  '
      if(it.eq. 9) title(1:24) = 'Madogram (absolute)     '
      if(it.eq.10) title(1:24) = 'Indicator 1/2 Variogram '
!      write(title(25:62),100) names(ivtail(iv)),names(ivhead(iv))
! 100  format('tail:',a12,' head:',a12)
!
! Loop over all the directions (note the direction in the title):
!
      do 1 id=1,ndir
!            write(title(62:74),101) id
 101        format('direction ',i2)
!            write(lout,'(a74)') title(1:74)
!
! Compute the unit lag distance along the directional vector:
!
            dis = sqrt((ixd(id)*xsiz)**2+(iyd(id)*ysiz)**2)
!
! Write out all the lags:
!
            do 1 il=1,nlag
                  i = (id-1)*nvarg*nlag+(iv-1)*nlag+il
                  disl = real(il)*dis
                  nump = int(np(i))
									call srite("out",cmplx(disl,gam(i)),8)
!                  if(it.eq.4) then
!                        write(lout,102) il,disl,gam(i),nump,&
!                                       hm(i),tm(i),hv(i),tv(i)
!                  else
!                        write(lout,102) il,disl,gam(i),nump,&
!                                        hm(i),tm(i)
!                  endif
 102              format(1x,i3,1x,f12.3,1x,f12.5,1x,i8,4(1x,f14.5))
 1    continue
!      close(lout)
      return
      end subroutine



      subroutine check(nvarg,ivtail,ivhead,ivtype,vrmin,vrmax,&
                       names,MAXVAR)
!-----------------------------------------------------------------------
!
!                Error Check and Note Variogram Types
!                ************************************
!
! Go through each variogram type and note the type to the screen and
! report any possible errors.
!
!
!
!
!
! Original: C.V. Deutsch                             Date: February 1992
!-----------------------------------------------------------------------
      real      vrmin(*),vrmax(*)
      integer   ivtail(*),ivhead(*),ivtype(*)
      character title*80,names(MAXVAR)*12
!
! Loop over all the variograms to be computed:
!
      do 1 iv=1,nvarg
!
! Note the variogram type and the variables being used:
!
      it = ivtype(iv)
      if(it.eq. 1) title(1:24) = 'Semivariogram          :'
      if(it.eq. 2) title(1:24) = 'Cross Semivariogram    :'
      if(it.eq. 3) title(1:24) = 'Covariance             :'
      if(it.eq. 4) title(1:24) = 'Correlogram            :'
      if(it.eq. 5) title(1:24) = 'General Relative       :'
      if(it.eq. 6) title(1:24) = 'Pairwise Relative      :'
      if(it.eq. 7) title(1:24) = 'Variogram of Logarithms:'
      if(it.eq. 8) title(1:24) = 'Rodogram (square root) :'
      if(it.eq. 9) title(1:24) = 'Madogram (absolute)    :'
      if(it.eq.10) title(1:24) = 'Indicator 1/2 Variogram:'
      write(title(25:64),100) names(ivtail(iv)),names(ivhead(iv))
 100  format('  tail=',a12,' head=',a12)
      write(0,101) iv,title(1:64)
 101  format(' Variogram ',i2,1x,a64)
!
! Check for possible errors or inconsistencies:
!
      if(it.eq.2) then
            if(ivtail(iv).eq.ivhead(iv)) write(0,201)
 201        format('  WARNING: cross variogram with the same variable!')
      else if(it.eq.5) then
            if(ivtail(iv).ne.ivhead(iv)) write(0,501)
            if(vrmin(ivtail(iv)).lt.0.0.and.vrmax(ivtail(iv)).gt.0.0)&
                 write(0,502)
            if(vrmin(ivhead(iv)).lt.0.0.and.vrmax(ivhead(iv)).gt.0.0)&
                  write(0,502)
 501        format('  WARNING: cross general relative variogram are',&
                   ' difficult to interpret!')
 502        format('  WARNING: there are both positive and negative',&
                   ' values - lag mean could be zero!')
      else if(it.eq.6) then
            if(ivtail(iv).ne.ivhead(iv)) write(0,601)
            if(vrmin(ivtail(iv)).lt.0.0.and.vrmax(ivtail(iv)).gt.0.0)&
                  write(0,602)
            if(vrmin(ivhead(iv)).lt.0.0.and.vrmax(ivhead(iv)).gt.0.0)&
                  write(0,602)
 601        format('  WARNING: cross pairwise relative variogram are',&
                   ' difficult to interpret!')
 602        format('  WARNING: there are both positive and negative',&
                   ' values - pair means could be zero!')
      else if(it.eq.7) then
            if(ivtail(iv).ne.ivhead(iv)) write(0,701)
            if(vrmin(ivtail(iv)).lt.0.0.or.vrmin(ivhead(iv)).lt.0.0)&
            write(0,702)
 701        format('  WARNING: cross logarithmic variograms may be',&
                   ' difficult to interpret!')
 702        format('  WARNING: there are zero or negative',&
                   ' values - logarithm undefined!')
      else if(it.eq.8) then
            if(ivtail(iv).ne.ivhead(iv)) write(0,801)
 801        format('  WARNING: cross rodograms may be difficult to',&
                   ' interpret!')
      else if(it.eq.9) then
            if(ivtail(iv).ne.ivhead(iv)) write(0,901)
 901        format('  WARNING: cross madograms may be difficult to',&
                   ' interpret!')
      endif
!
! Loop over all variograms:
!
 1    continue
      return
      end subroutine
end module

