      program main
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!                                                                      %
! Copyright (C) 1996, The Board of Trustees of the Leland Stanford     %
! Junior University.  All rights reserved.                             %
!                                                                      %
! 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.       %
!                                                                      %
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!-----------------------------------------------------------------------
!
!                 Kriging of a 2-D Rectangular Grid
!                 *********************************
!
! The program is executed with no command line arguments.  The user
! will be prompted for the name of a parameter file.  The parameter
! file is described in the documentation (see the example kb2d.par)
!
!
!
!-----------------------------------------------------------------------
use kbmod
!
! Read the Parameter File and the Data:
!
     call grab_pars()
!
! Call kb2d to krige the grid:
!
      call kb2d
!
! Finished:
!
      write(0,9998) VERSION
 9998 format(/' KB2D Version: ',f5.3, ' Finished'/)
      stop
      end
 
 
      subroutine kb2d
!-----------------------------------------------------------------------
!
!           Ordinary/Simple Kriging of a 2-D Rectangular Grid
!           *************************************************
!
! This subroutine estimates point or block values of one variable by
! ordinary kriging.  All of the samples are rescanned for each block
! estimate; this makes the program simple but inefficient.  The data
! should NOT contain any missing values.  Unestimated points are
! returned as -1.0e21
!
!
!
! Original:  A.G. Journel                                           1978
! Revisions: B.E. Buxton                                       Apr. 1983
!-----------------------------------------------------------------------
     use kbmod
      real      xdb(MAXDIS),ydb(MAXDIS),xa(MAXSAM),ya(MAXSAM),&
                vra(MAXSAM),dist(MAXSAM)
      real*8    r(MAXSAM+1),rr(MAXSAM+1),s(MAXSAM+1),a(MAXKRG)
      integer   nums(MAXSAM)
      logical   first
      data      first/.true./,PMX/9999.0/
!
! Echo the input parameters if debugging flag is >2:
!
      if(idbg.gt.2) then
            write(ldbg,*) 'KB2D Parameters'
            write(ldbg,*)
            write(ldbg,*) 'Variogram Parameters for ',nst,' structures:'
            write(ldbg,*) '  Nugget effect:         ',c0
            write(ldbg,*) '  Types of variograms:   ',(it(i),i=1,nst)
            write(ldbg,*) '  Contribution cc        ',(cc(i),i=1,nst)
            write(ldbg,*) '  Ranges:                ',(aa(i),i=1,nst)
            write(ldbg,*) '  Angle for Continuity:  ',(ang(i),i=1,nst)
            write(ldbg,*) '  Anisotropy Factors:    ',(anis(i),i=1,nst)
            write(ldbg,*) ' '
            write(ldbg,*) 'Grid for Kriging:'
            write(ldbg,*) '  Number of X and Y Blocks:',nx,ny
            write(ldbg,*) '  Origin of X and Y Blocks:',xmn,ymn
            write(ldbg,*) '  Size   of X and Y Blocks:',xsiz,ysiz
            write(ldbg,*) ' '
            write(ldbg,*) 'Discretization of blocks:  ',nxdis,nydis
            write(ldbg,*) 'Search Radius:             ',radius
            write(ldbg,*) 'Minimum number of samples: ',ndmin
            write(ldbg,*) 'Maximum number of samples: ',ndmax
            write(ldbg,*) ' '
      endif
!
! Echo the input data if debugging flag >1:
!
      if(idbg.ge.4) then
            do id=1,nd
                  write(ldbg,99) id,x(id),y(id),vr(id)
 99               format('Data: ',i5,' at ',2f12.3,' value: ',f12.5)
            end do
      endif
!
! Set up the discretization points per block.  Figure out how many
! are needed, the spacing, and fill the xdb and ydb arrays with the
! offsets relative to the block center (this only gets done once):
!
      ndb  = nxdis * nydis
      if(ndb.gt.MAXDIS) then
            write(0,*) 'ERROR KB2D: Too many discretization points '
            write(0,*) '            Increase MAXDIS or lower n[xy]dis'
            stop
      endif
      xdis = xsiz  / max(real(nxdis),1.0)
      ydis = ysiz  / max(real(nydis),1.0)
      xloc = -0.5*(xsiz+xdis)
      i    = 0
      do ix =1,nxdis
            xloc = xloc + xdis
            yloc = -0.5*(ysiz+ydis)
            do iy=1,nydis
                  yloc = yloc + ydis
                  i = i+1
                  xdb(i) = xloc
                  ydb(i) = yloc
            end do
      end do
!
! Initialize accumulators:
!
      cbb  = 0.0
      rad2 = radius*radius
!
! Calculate Block Covariance. Check for point kriging.
!
      cov   = cova2(xdb(1),ydb(1),xdb(1),ydb(1),nst,c0,PMX,cc,&
                    aa,it,ang,anis,first)
!
! Keep this value to use for the unbiasedness constraint:
!
      unbias = cov
      first  = .false.
      if (ndb.le.1) then
            cbb = cov
      else
            do i=1,ndb
                  do j=1,ndb
                        cov = cova2(xdb(i),ydb(i),xdb(j),ydb(j),nst,c0,&
                                   PMX,cc,aa,it,ang,anis,first)
                        if(i.eq.j) cov = cov - c0
                        cbb = cbb + cov
                  end do
            end do
            cbb = cbb/real(ndb*ndb)
      endif
      if(idbg.gt.1) then
            write(ldbg,*) ' '
            write(ldbg,*) 'Block Covariance: ',cbb
            write(ldbg,*) ' '
      endif
!
! MAIN LOOP OVER ALL THE BLOCKS IN THE GRID:
!
      nk = 0
      ak = 0.0
      vk = 0.0
      do 4 iy=1,ny
      yloc = ymn + (iy-1)*ysiz
      do 4 ix=1,nx
            xloc = xmn + (ix-1)*xsiz
!
! Find the nearest samples within each octant: First initialize
! the counter arrays:
!
            na = 0
            do isam=1,ndmax
                  dist(isam) = 1.0e+20
                  nums(isam) = 0
            end do
!
! Scan all the samples (this is inefficient and the user with lots of
! data should move to ktb3d):
!
            do 6 id=1,nd
                  dx = x(id) - xloc
                  dy = y(id) - yloc
                  h2 = dx*dx + dy*dy
                  if(h2.gt.rad2) go to 6
!
! Do not consider this sample if there are enough close ones:
!
                  if(na.eq.ndmax.and.h2.gt.dist(na)) go to 6
!
! Consider this sample (it will be added in the correct location):
!
                  if(na.lt.ndmax) na = na + 1
                  nums(na)           = id
                  dist(na)           = h2
                  if(na.eq.1) go to 6
!
! Sort samples found thus far in increasing order of distance:
!
                  n1 = na-1
                  do ii=1,n1
                        k=ii
                        if(h2.lt.dist(ii)) then
                              jk = 0
                              do jj=k,n1
                                    j  = n1-jk
                                    jk = jk+1
                                    j1 = j+1
                                    dist(j1) = dist(j)
                                    nums(j1) = nums(j)
                              end do
                              dist(k) = h2
                              nums(k) = id
                              go to 6
                        endif
                  end do
 6          continue
!
! Is there enough samples?
!
            if(na.lt.ndmin) then
                  if(idbg.ge.2)&
                 write(ldbg,*) 'Block ',ix,iy, 'not estimated'
                  est  = UNEST
                  estv = UNEST
                  go to 1
            endif
!
! Put coordinates and values of neighborhood samples into xa,ya,vra:
!
            do ia=1,na
                  jj      = nums(ia)
                  xa(ia)  = x(jj)
                  ya(ia)  = y(jj)
                  vra(ia) = vr(jj)
            end do
!
! Handle the situation of only one sample:
!
            if(na.eq.1) then
                  cb1 = cova2(xa(1),ya(1),xa(1),ya(1),nst,c0,&
                              PMX,cc,aa,it,ang,anis,first)
                  xx  = xa(1) - xloc
                  yy  = ya(1) - yloc
!
! Establish Right Hand Side Covariance:
!
                  if(ndb.le.1) then
                        cb = cova2(xx,yy,xdb(1),ydb(1),nst,c0,&
                                   PMX,cc,aa,it,ang,anis,first)
                  else
                        cb  = 0.0
                        do i=1,ndb
                              cb = cb + cova2(xx,yy,xdb(i),ydb(i),nst,&
                                        c0,PMX,cc,aa,it,ang,anis,first)
                              dx = xx - xdb(i)
                              dy = yy - ydb(i)
                              if((dx*dx+dy*dy).lt.EPSLON)&
                              cb = cb - c0
                        end do
                        cb = cb / real(ndb)
                  end if
                  if(ktype.eq.0) then
                        s(1) = cb/cbb
                        est  = s(1)*vra(1) + (1.0-s(1))*skmean
                        estv = cbb - s(1) * cb
                  else
                        est  = vra(1)
                        estv = cbb - 2.0*cb + cb1
                  end if
            else
!
! Solve the Kriging System with more than one sample:
!
                  neq = na + ktype
                  nn  = (neq + 1)*neq/2
!
! Set up kriging matrices:
!
                  in=0
                  do j=1,na
!
! Establish Left Hand Side Covariance Matrix:
!
                        do i=1,j
                              in = in + 1
                              a(in) = dble( cova2(xa(i),ya(i),xa(j),&
                                            ya(j),nst,c0,PMX,cc,aa,it,&
                                            ang,anis,first) )
                        end do
                        xx = xa(j) - xloc
                        yy = ya(j) - yloc
!
! Establish Right Hand Side Covariance:
!
                        if(ndb.le.1) then
                              cb = cova2(xx,yy,xdb(1),ydb(1),nst,c0,&
                                         PMX,cc,aa,it,ang,anis,first)
                        else
                              cb  = 0.0
                              do j1=1,ndb
                                    cb = cb + cova2(xx,yy,xdb(j1),&
                                         ydb(j1),nst,c0,PMX,cc,aa,&
                                         it,ang,anis,first)
                                    dx = xx - xdb(j1)
                                    dy = yy - ydb(j1)
                                    if((dx*dx+dy*dy).lt.EPSLON)&
                                          cb = cb - c0
                              end do
                              cb = cb / real(ndb)
                        end if
                        r(j)  = dble(cb)
                        rr(j) = r(j)
                  end do
!
! Set the unbiasedness constraint:
!
                  if(ktype.eq.1) then
                        do i=1,na
                              in    = in + 1
                              a(in) = dble(unbias)
                        end do
                        in      = in + 1
                        a(in)   = 0.0
                        r(neq)  = dble(unbias)
                        rr(neq) = r(neq)
                  end if
!
! Write out the kriging Matrix if Seriously Debugging:
!
                  if(idbg.ge.3) then
                        write(ldbg,101) ix,iy
                        is = 1
                        do i=1,neq
                              ie = is + i - 1
                              write(ldbg,102) i,r(i),(a(j),j=is,ie)
                              is = is + i
                        end do
 101                    format(/,'Kriging Matrices for Node: ',2i4,&
                                 ' RHS first')
 102                    format('    r(',i2,') =',f7.4,'  a= ',9(10f7.4))
                  endif
!
! Solve the Kriging System:
!
                  call ksol(1,neq,1,a,r,s,ising)
!
! Write a warning if the matrix is singular:
!
                  if(ising.ne.0) then
                        write(0,*) 'WARNING KB2D: singular matrix'
                        write(0,*) '              for block',ix,iy
                        est  = UNEST
                        estv = UNEST
                        go to 1
                  endif
!
! Write the kriging weights and data if requested:
!
                  if(idbg.ge.2) then
                        write(ldbg,*) '       '
                        write(ldbg,*) 'BLOCK: ',ix,iy
                        write(ldbg,*) '       '
                        if(ktype.eq.1) write(ldbg,*) &
                        '  Lagrange multiplier: ',s(neq)*unbias
                        write(ldbg,*) '  BLOCK EST: x,y,vr,wt '
                        do i=1,na
                        write(ldbg,'(4f12.3)') xa(i),ya(i),vra(i),s(i)
                        end do
                  endif
!
! Compute the estimate and the kriging variance:
!
                  est  = 0.0
                  estv = cbb
                  sumw = 0.0
                  if(ktype.eq.1) estv = estv - real(s(na+1))
                  do i=1,na
                        sumw = sumw + real(s(i))
                        est  = est  + real(s(i))*vra(i)
                        estv = estv - real(s(i)*rr(i))
                  end do
                  if(ktype.eq.0) est = est + (1.0-sumw)*skmean
            endif
            if(idbg.ge.2) then
                  write(ldbg,*) '  est  ',est
                  write(ldbg,*) '  estv ',estv
                  write(ldbg,*) ' '
            endif
!
! Write the result to the output file:
!
 1          write(lout,'(f8.3,1x,f8.3)') est,estv
            call srite("out",est,4)
            if(est.gt.UNEST) then
                  nk = nk + 1
                  ak = ak + est
                  vk = vk + est*est
            end if
!
! END OF MAIN LOOP OVER ALL THE BLOCKS:
!
 4    continue
      if(nk.ge.1) then
            ak = ak / real(nk)
            vk = vk/real(nk) - ak*ak
            write(ldbg,105) nk,ak,vk
            write(0,   105) nk,ak,vk
 105        format(/,'Estimated   ',i8,' blocks ',/,&
                     '  average   ',f9.4,/,'  variance  ',f9.4,/)
      end if
      return
      end
 
 
 
      real function cova2(x1,y1,x2,y2,nst,c0,PMX,cc,aa,it,&
                          ang,anis,first)
!-----------------------------------------------------------------------
!
!              Covariance Between Two Points (2-D Version)
!              *******************************************
!
! This function returns the covariance associated with a variogram model
! that is specified by a nugget effect and possibly four different
! nested varigoram structures.  The anisotropy definition can be
! different for each of the nested structures (spherical, exponential,
! gaussian, or power).
!
!
!
! INPUT VARIABLES:
!
!   x1,y1            Coordinates of first point
!   x2,y2            Coordinates of second point
!   nst              Number of nested structures (max. 4).
!   c0               Nugget constant (isotropic).
!   PMX              Maximum variogram value needed for kriging when
!                      using power model.  A unique value of PMX is
!                      used for all nested structures which use the
!                      power model.  therefore, PMX should be chosen
!                      large enough to account for the largest single
!                      structure which uses the power model.
!   cc(nst)          Multiplicative factor of each nested structure.
!   aa(nst)          Parameter "a" of each nested structure.
!   it(nst)          Type of each nested structure:
!                      1. spherical model of range a;
!                      2. exponential model of parameter a;
!                           i.e. practical range is 3a
!                      3. gaussian model of parameter a;
!                           i.e. practical range is a*sqrt(3)
!                      4. power model of power a (a must be gt. 0  and
!                           lt. 2).  if linear model, a=1,c=slope.
!   ang(nst)         Azimuth angle for the principal direction of
!                      continuity (measured clockwise in degrees from Y)
!   anis(nst)        Anisotropy (radius in minor direction at 90 degrees
!                      from "ang" divided by the principal radius in 
!                      direction "ang")
!   first            A logical variable which is set to true if the
!                      direction specifications have changed - causes
!                      the rotation matrices to be recomputed.
!
!
!
! OUTPUT VARIABLES: returns "cova2" the covariance obtained from the
!                   variogram model.
!
!
!
!-----------------------------------------------------------------------
      parameter(DTOR=3.14159265/180.0,EPSLON=0.0000001)
      real      aa(*),cc(*),ang(*),anis(*),rotmat(4,4),maxcov
      integer   it(*)
      logical   first
      save      rotmat,maxcov
!
! The first time around, re-initialize the cosine matrix for the
! variogram structures:
!
      if(first) then
            maxcov = c0
            do is=1,nst
                  azmuth       = (90.0-ang(is))*DTOR
                  rotmat(1,is) =  cos(azmuth)
                  rotmat(2,is) =  sin(azmuth)
                  rotmat(3,is) = -sin(azmuth)
                  rotmat(4,is) =  cos(azmuth)
                  if(it(is).eq.4) then
                        maxcov = maxcov + PMX
                  else
                        maxcov = maxcov + cc(is)
                  endif
            end do
      endif
!
! Check for very small distance:
!
      dx = x2-x1
      dy = y2-y1
      if((dx*dx+dy*dy).lt.EPSLON) then
            cova2 = maxcov
            return
      endif
!
! Non-zero distance, loop over all the structures:
!
      cova2 = 0.0
      do is=1,nst
!
! Compute the appropriate structural distance:
!
            dx1 = (dx*rotmat(1,is) + dy*rotmat(2,is))
            dy1 = (dx*rotmat(3,is) + dy*rotmat(4,is))/anis(is)
            h   = sqrt(max((dx1*dx1+dy1*dy1),0.0))
            if(it(is).eq.1) then
!
! Spherical model:
!
                  hr = h/aa(is)
                  if(hr.lt.1.0) cova2 = cova2 &
                                      + cc(is)*(1.-hr*(1.5-.5*hr*hr))
            else if(it(is).eq.2) then
!
! Exponential model:
!
                  cova2 = cova2 +cc(is)*exp(-h/aa(is))
            else if(it(is).eq. 3) then
!
! Gaussian model:
!
                  hh=-(h*h)/(aa(is)*aa(is))
                  cova2 = cova2 +cc(is)*exp(hh)
            else
!
! Power model:
!
                  cov1  = PMX - cc(is)*(h**aa(is))
                  cova2 = cova2 + cov1
            endif
      end do
      return
      end



      subroutine makepar
!-----------------------------------------------------------------------
!
!                      Write a Parameter File
!                      **********************
!
!
!
!-----------------------------------------------------------------------
      lun = 99
      open(lun,file='kb2d.par',status='UNKNOWN')
      write(lun,10)
 10   format('                  Parameters for KB2D',/,&
             '                  *******************',/,/,&
             'START OF PARAMETERS:')

      write(lun,11)
 11   format('../data/cluster.dat          ','-file with data')
      write(lun,12)
 12   format('1   2   3                    ', '-   columns for X, Y, and variable')
      write(lun,13)
 13   format('-1.0e21   1.0e21             ',&
             '-   trimming limits')
      write(lun,14)
 14   format('3                            ',&
             '-debugging level: 0,1,2,3')
      write(lun,15)
 15   format('kb2d.dbg                     ',&
             '-file for debugging output')
      write(lun,16)
 16   format('kb2d.out                     ',&
             '-file for kriged output')
      write(lun,17)
 17   format('5    5.0  10.0               ',&
             '-nx,xmn,xsiz')
      write(lun,18)
 18   format('5    5.0  10.0               ',&
             '-ny,ymn,ysiz')
      write(lun,19)
 19   format('1    1                       ',&
             '-x and y block discretization')
      write(lun,20)
 20   format('4    8                       ',&
             '-min and max data for kriging')
      write(lun,21)
 21   format('20.0                         ',&
             '-maximum search radius')
      write(lun,22)
 22   format('1    2.302                   ',&
             '-0=SK, 1=OK,  (mean if SK)')
      write(lun,23)
 23   format('1   2.0                      ',&
             '-nst, nugget effect')
      write(lun,24)
 24   format('1   8.0  0.0  10.0  10.0     ',&
             '-it, c, azm, a_max, a_min')

      close(lun)
      return
      end
