module krig_base ! that this notice and the above copyright notice remain intact. % ! % !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !----------------------------------------------------------------------- ! ! Ordinary Kriging of a 2-D Rectangular Grid ! ****************************************** ! ! The following Parameters control static dimensioning within okb2d: ! ! MAXX maximum nodes in X ! MAXY maximum nodes in Y ! MAXDAT maximum number of data points ! MAXSAM maximum number of data points to use in one kriging system ! MAXDIS maximum number of discretization points per block ! MAXNST maximum number of nested structures ! ! MAXKRG MAXSAM*MAXSAM - used for dimensioning ! ! UNEST Assigned to unestimated blocks ! !----------------------------------------------------------------------- ! ! User Adjustable: ! parameter(MAXX = 500, MAXY = 500,MAXDAT = 10000,& MAXSAM = 1000, MAXDIS = 80, MAXNST = 4) ! ! Fixed ! parameter(MAXKD=MAXSAM+1,MAXKRG=MAXKD*MAXKD,UNEST=-999.,& EPSLON=0.0000001,VERSION=1.202) end module module krig_it use krig_base contains subroutine okb2d(nd,x,y,vr,nx,ny,xmn,ymn,xsiz,ysiz,nxdis,nydis,& ndmin,ndmax,radius,nst,c0,it,aa,cc,ang,anis,& idbg,ldbg,est,estv) !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! % ! 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. % ! % !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% !----------------------------------------------------------------------- ! ! Ordinary 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 ! ! ! ! INPUT VARIABLES: ! ! nd Number of data (no missing values) ! x(nd),y(nd) X and Y coordinates of the data ! vr(nd) Data values ! nx,ny Number of blocks in X and Y ! xmn Coordinate at the center of the first X Block ! ymn Coordinate at the center of the first Y Block ! all blocks or grid nodes will have X and Y ! coordinates equal to or greater than xmn and ymn ! xsiz,ysiz X and Y spacing of the grid nodes (block size) ! nxdis Number of discretization points/block in X ! nydis Number of discretization points/block in Y ! ndmin Minimum number of samples required for kriging ! ndmax Maximum number of samples to use in kriging. ! radius Maximum search radius (isotropic) ! ! nst Number of variogram structures ! c0 Nugget effect ! it(nst) Variogram type (1=sph,2=exp,3=gaus,4=powr) ! aa(nst) Range parameter except for the power model ("aa" is ! the power parameter) ! cc(nst) Contribution of each nested structure except for ! the power model where "cc" is the slope ! ang(nst) Azimuth angle for anisotropy (measured clockwise ! from the Y axis) for each variogram structure. ! anis(nst) Anisotropy ratio, i.e., the range in the direction ! 90 degrees from "ang" divided by the range given ! by "aa". For example, if the variogram model has ! a range of 20 at an azimuth of N30E and a range ! of 10 in the direction N120E the parameters ! should be entered as: aa=20, ang=30, anis=0.5 ! idbg The debug level: ! 0. None ! 1. Input parameters and summary statistics ! 2. The data used for each kriging along with the ! kriging weights are printed. ! 3. Kriging weights, values and variances / block ! kriging matrices also printed (serious) ! ldbg The output unit for the debugging output ! ! ! ! OUTPUT VARIABLES: ! ! est(nx,ny) the kriged values ! estv(nx,ny) the kriging variances ! ! ! ! EXTERNAL REFERENCES: ! ! cova2 Calculates the covariance given a variogram model ! ksol Linear system solver using Gaussian Elimination ! (double precision) ! ! ! ! Original: A.G. Journel 1978 ! Revisions: B.E. Buxton Apr. 1983 ! C.V. Deutsch July 1990 !----------------------------------------------------------------------- real x(MAXDAT),y(MAXDAT),vr(MAXDAT),aa(MAXNST),cc(MAXNST) real ang(MAXNST),anis(MAXNST),est(MAXX,MAXY),estv(MAXX,MAXY) real xdb(MAXDIS),ydb(MAXDIS),xa(MAXSAM),ya(MAXSAM) real vra(MAXSAM),dist(MAXSAM) real*8 r(MAXSAM+1),rr(MAXSAM+1),s(MAXSAM+1),a(MAXKRG) integer it(MAXNST),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,*) 'OKB2D 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 1 id=1,nd write(ldbg,99) id,x(id),y(id),vr(id) 99 format('Data: ',i5,' at ',2f12.3,' value: ',f12.5) 1 continue 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(*,*) 'ERROR OKB2D: Too many discretization points ' write(*,*) ' Increase MAXDIS or lower n[xy]dis' stop endif xdis = xsiz / amax1(real(nxdis),1.0) ydis = ysiz / amax1(real(nydis),1.0) xloc = -0.5*(xsiz+xdis) i = 0 do 2 ix =1,nxdis xloc = xloc + xdis yloc = -0.5*(ysiz+ydis) do 2 iy=1,nydis yloc = yloc + ydis i = i+1 xdb(i) = xloc ydb(i) = yloc 2 continue ! ! Initialize accumulators: ! uk = 0.0 vk = 0.0 nk = 0 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 3 i=1,ndb do 3 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 3 continue 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: ! 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 5 isam=1,ndmax dist(isam) = 1.0e+20 nums(isam) = 0 5 continue ! ! 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 7 ii=1,n1 k=ii if(h2.lt.dist(ii)) then jk = 0 do 8 jj=k,n1 j = n1-jk jk = jk+1 j1 = j+1 dist(j1) = dist(j) nums(j1) = nums(j) 8 continue dist(k) = h2 nums(k) = id go to 6 endif 7 continue 6 continue ! ! Is there enough samples? ! if(na.lt.ndmin) then if(idbg.ge.2)write(ldbg,*) 'Block ',ix,iy, 'not estimated' est(ix,iy) = UNEST estv(ix,iy) = UNEST go to 4 endif ! ! Put coordinates and values of neighborhood samples into xa,ya,vra: ! do 10 ia=1,na jj = nums(ia) xa(ia) = x(jj) ya(ia) = y(jj) vra(ia) = vr(jj) 10 continue ! ! 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 12 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 12 continue cb = cb / real(ndb) end if est(ix,iy) = vra(1) estv(ix,iy) = cbb - 2.0*cb + cb1 else ! ! Solve the Kriging System with more than one sample: ! neq = na + 1 nn = (neq + 1)*neq/2 ! ! Set up kriging matrices: ! in=0 do 13 j=1,na ! ! Establish Left Hand Side Covariance Matrix: ! do 14 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) ) 14 continue 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 15 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 15 continue cb = cb / real(ndb) end if r(j) = dble(cb) rr(j) = r(j) 13 continue ! ! Set the unbiasedness constraint: ! do 16 i=1,na in = in + 1 a(in) = dble(unbias) 16 continue in = in + 1 a(in) = 0.0 r(neq) = dble(unbias) rr(neq) = r(neq) ! ! Write out the kriging Matrix if Seriously Debugging: ! if(idbg.ge.3) then write(ldbg,101) ix,iy is = 1 do 100 i=1,neq ie = is + i - 1 write(ldbg,102) i,r(i),(a(j),j=is,ie) is = is + i 100 continue 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(*,*) 'WARNING OKB2D: singular matrix' write(*,*) ' for block',ix,iy est(ix,iy) = UNEST estv(ix,iy) = UNEST go to 4 endif ! ! Write the kriging weights and data if requested: ! if(idbg.ge.2) then write(ldbg,*) ' ' write(ldbg,*) 'BLOCK: ',ix,iy write(ldbg,*) ' ' write(ldbg,*) ' Lagrange multiplier: ',s(neq) write(ldbg,*) ' BLOCK EST: x,y,vr,wt ' do 18 i=1,na 18 write(ldbg,'(4f12.3)') xa(i),ya(i),vra(i),s(i) endif ! ! Compute the estimate and the kriging variance: ! ook = 0.0 ookv = cbb-real(s(na+1)) do 19 i=1,na ook = ook + real(s(i))*vra(i) ookv = ookv - real(s(i))*real(rr(i)) 19 continue est(ix,iy) = ook estv(ix,iy) = ookv endif ! ! Accumulate statistics of kriged blocks: ! nk = nk + 1 uk = uk + est(ix,iy) vk = vk + est(ix,iy)*est(ix,iy) ! ! Write more details if requested: ! if(idbg.ge.2) then ! call srite("out",est(ix,iy),4) write(ldbg,*) ' est ',est(ix,iy) write(ldbg,*) ' estv ',estv(ix,iy) write(ldbg,*) ' ' endif ! ! END OF MAIN LOOP OVER ALL THE BLOCKS: ! 4 continue ! ! Write statistics of kriged values: ! if(nk.gt.0.and.idbg.gt.0) then vk = (vk-uk*uk/real(nk))/real(nk) uk = uk/real(nk) write(ldbg,*) ' ' write(ldbg,*) 'Estimated ',nk,' blocks ' write(ldbg,*) ' average ',uk write(ldbg,*) ' variance ',vk endif return end subroutine 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. ! ! ! ! AUTHOR: C.V. Deutsch DATE: July 1990 !----------------------------------------------------------------------- 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 1 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 1 continue 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 2 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(dx1*dx1+dy1*dy1) if(it(is).eq.1) then ! ! Spherical model: ! hr = h/aa(is) if(hr.ge.1.0) go to 2 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 2 continue write(0,*) "called with ",x1,x2,y1,y2,cova2 return end function subroutine ksol(nright,neq,nsb,a,r,s,ising) !----------------------------------------------------------------------- ! ! Solution of a System of Linear Equations ! **************************************** ! ! ! ! INPUT VARIABLES: ! ! nright,nsb number of columns in right hand side matrix. ! for OKB2D: nright=1, nsb=1 ! neq number of equations ! a() upper triangular left hand side matrix (stored ! columnwise) ! r() right hand side matrix (stored columnwise) ! for okb2d, one column per variable ! ! ! ! OUTPUT VARIABLES: ! ! s() solution array, same dimension as r above. ! ising singularity indicator ! 0, no singularity problem ! -1, neq .le. 1 ! k, a null pivot appeared at the kth iteration ! ! ! ! PROGRAM NOTES: ! ! 1. Requires the upper triangular left hand side matrix. ! 2. Pivots are on the diagonal. ! 3. Does not search for max. element for pivot. ! 4. Several right hand side matrices possible. ! 5. USE for ok and sk only, NOT for UK. ! ! !----------------------------------------------------------------------- implicit real*8 (a-h,o-z) real*8 a(*),r(*),s(*) ! ! If there is only one equation then set ising and return: ! if(neq.le.1) then ising = -1 return endif ! ! Initialize: ! tol = 0.1e-06 ising = 0 nn = neq*(neq+1)/2 nm = nsb*neq m1 = neq-1 kk = 0 ! ! Start triangulation: ! do 70 k=1,m1 kk=kk+k ak=a(kk) if(abs(ak).lt.tol) then ising=k return endif km1=k-1 do 60 iv=1,nright nm1=nm*(iv-1) ii=kk+nn*(iv-1) piv=1./a(ii) lp=0 do 50 i=k,m1 ll=ii ii=ii+i ap=a(ii)*piv lp=lp+1 ij=ii-km1 do 30 j=i,m1 ij=ij+j ll=ll+j a(ij)=a(ij)-ap*a(ll) 30 continue do 40 llb=k,nm,neq in=llb+lp+nm1 ll1=llb+nm1 r(in)=r(in)-ap*r(ll1) 40 continue 50 continue 60 continue 70 continue ! ijm=ij-nn*(nright-1) if(abs(a(ijm)).lt.tol) then ising=neq return endif ! ! Finished triangulation, start solving back: ! do 140 iv=1,nright nm1=nm*(iv-1) ij=ijm+nn*(iv-1) piv=1./a(ij) do 100 llb=neq,nm,neq ll1=llb+nm1 s(ll1)=r(ll1)*piv 100 continue i=neq kk=ij do 130 ii=1,m1 kk=kk-i piv=1./a(kk) i=i-1 do 120 llb=i,nm,neq ll1=llb+nm1 in=ll1 ap=r(in) ij=kk do 110 j=i,m1 ij=ij+j in=in+1 ap=ap-a(ij)*s(in) 110 continue s(ll1)=ap*piv 120 continue 130 continue 140 continue ! ! Finished solving back, return: ! return end subroutine end module