PB87-170353
Disjunctive Kriging Program for
Two Dimensions
(U.S.) Robert S. Kerr Environmental
Research Lab., Ada, OK
1986
J
-------
TECHNICAL REPORT DATA
(Pteott rtaa Instructions on the rtvtnt before completii
[1. REPORT NO.
EPA/600/J-86/247
2.
3.1
*fid7-170J53
.TITLE AND SUBTITLE
A DISJUNCTIVE KRIGING PROGRAM FOR TWO DIMENSIONS
(Journal Version)
5. REPORT DATE
1986
.PERFORMING ORGANIZATION CODE
7. AUTHOR(S)
8. PERFORMING ORGANIZATION REPORT NO.
S. R. Yates,1 A. W. Warrick2 and D. E. Myers2
9. PERFORMING ORGANIZATION NAME AND ADDRESS
^.S. Kerr Environmental Research Laboratory,
U.S. EPA, Ada, OK 74820
2The University of Arizona, Tucson, AZ 85721
10. PROGRAM ELEMENT NO.
ABWD1A
11. CONTRACT/GRANT NO.
In-House
12. SPONSORING AGENCY NAME AND ADDRESS
Robert S. Kerr Environmental Research Lab. - Ada, OK
U.S. Environmental Protection Agency
Post Office Box 1198
Ada, OK 74820
13. TVPE OF REPORT AND PERIOD COVERED
Journal Article
14. SPONSORING AGENCV CODE
EPA/600/15
is. SUPPLEMENTARY NOTES
Published in COMPUTERS & GEOSCIENCES, 12(3):281-313, 1986.
16. ABSTRACT
This paper describes the disjunctive kriging (OK) method for two dimensional
spatially variable properties. A brief mathematical description is given which
includes all pertinent equations to obtain an estimated value, disjunctive
kriging variance, and conditional probability that the value of a property at a
location is above a known cutoff level. This is followed by a description of the
steps which are necessary to implement DK and an example illustrating the method.
In order to use DK, a series of complex calculations must be carried out. To
facilitate this two FORTRAN programs were developed. The programs, example
input and output files as well as information necessary to use the programs is
included.
17.
KEY WORDS AND DOCUMENT ANALYSIS
DESCRIPTORS
b.IDENTIFIERS/OPEN ENDED TERMS C. COSATI Field/Croup
REPRODUCED BY
U.S. DEPARTMENT OF COMMERCE
HATCNM.TECHMCAL
SPflNGFELD.VA 22181
18. DISTRIBUTION STATEMENT
RELEASE TO THE PUBLIC.
tPA fotm 2220-J (R«v. 4-77) »«lviou§ «OITIOM i* OBSOLCTC
19. SECURITY CLASS iTMl Report)
UNCLASSIFIED.
21. NO. OF PAGES
20. SECURITY CLASS iThil pott)
UNCLASSIFIED.
22. PRICE
-------
EPA/600/J-86/247
JOURNAL ARTICLE
uitn A Gnurimm Vol. 12. No. >. pp. MI-JI3. 1986
Printed in Grai Briuin
OWX-MXM/M S3.00 t 0.00
Pergamon Journal* Ltd
A DISJUNCTIVE KRIGING PROGRAM FOR TWO DIMENSIONS
S.R. YATO
R. S. Kerr Environmental Research Uboratory. P.O. Box 1198. Ada. OK 74820. U.S.A.
A.W. WARRICK
Department of Soils, Water and Engineering. The University of Arizona. Tucson. AZ 85721. U.S.A.
and
D.E. MYERS
Department of Mathematics. The University of Arizona. Tucson. AZ 85721. U.S.A.
(Rfcrirfd 27 July 1985: retard 20 January 1986)
AbstractThis paper describes the disjunctive kriging (DK) method for two dimensional spatially variable
properties. A brief mathematical description is given which includes all pertinent equations to obtain an
cstimati-d value, disjunctive kriging variance, and conditional probability that the value of a property at
a location is above a known cutoff level. This is followed by a description of the steps which are necessary
to implement DK and an example illustrating the method. In order to use DK. a series of complex
calculations must be carried out. To facilitate this two FORTRAN programs were developed. The
programs, example input and output files as well as information necessary to use the programs is included.
Key tt'ordf. Conditional probability. Disjunctive kriging. Joint density. Kriging. Nonlinear estimate.
Regionalized random variable. Spatial variability.
INTRODUCTION
Recently, much attention has been given lo linear
kriging methods in the analysis of spatially dependent
phenomena. These include simple, ordinary, and uni-
versal kriging either on punctual or block support.
Many examples occur in the literature. A few of these
include Burgess and Webster (1980a) who used punc-
tual kriging in the analysis of the spatial distribution
of the sodium content, cover loam, and stone content.
Burgess and Webster (1980b) and Webster and Bur-
gess (1980) also used block and universal kriging in
subsequent analyses. Warrick. Myers, and Nielsen
(1985) give an example of kriging the electrical con-
ductivity for a southwestern Arizona soil (soil type:
Typic Haplargid). Vieria, Nielsen, and Biggar (1981)
investigated the correlation between krigcd and actual
values based on different number of sample values
used in the estimation process in an attempt lo deter-
mine the minimum number of samples necessary to
obtain a given level of information. Vauclin and
others (1983) used ordinary kriging and cokriging to
estimate the spatial distribution of the available water
content and sand content of a Tunisian field. They
determined that cokriging provided improved esti-
mates of the available water content because of ad-
ditional information added to the problem by the
sand content. Journel and Huijbregts (1978) give
many examples of kriging which have mining applica-
tions.
These studies represent but a small sample of the
work which uses linear kriging estimators in the
analysis of the spatial variability of a physical
property. Generally, the linear kriging estimator is
not the best possible estimator. The minimum vari-
ance unbiased estimator of a random variable Y in
terms of random variables AT,. X} X. is me con-
ditional expectation of Y given X,, X}...., X,. For
the most general situation, computing this conditional
expectation requires knowledge of the joint density of
Y. Xt. X: X. although this information is dif-
ficult to obtain in practice. This nonlinear estimator
has the form
x.)
(i)
It also is possible to relax the requirement that the
joint density of (n + I) variables be known and define
another nonlinear estimator
- £/(*,)
(2)
where each/ is a nonlinear function of one X variable
only. This is the disjunctive kriging (DK) estimator
and requires that only the bivariate densities be
known.
The linear kriging estimator is a special form of
Equation (2) where each / is a linear function and
therefore only the constants need to be determined
i,X,.
(3)
In terms of estimating the value of a random variable
at an unsampled location. Equation (2) generally is a
NTIS a wtkorbri to repntfKi Md ut tkte
report. Ptnataici lor turtta npradticttM
aura te obltiittf Inn ttu upr>K)M cwur.
281
-------
282
S. R. YATCS. A. W. WARRICK. and D. E. MYERS
belter estimator than Equation (3) in the sense of
reduced kriging variance. Moreover, if an estimate of
the conditional probability distribution is required a
nonlinear estimator is essential.
The purpose of this paper is to present the disjunc-
tive kriging programs and to describe the basic
algorithm used to obtain a nonlinear estimate of a
regionalized random variable (at an unsampled loca-
tion), the estimation variance and the conditional
probability that the estimate is above a given cutoff
level. For brevity, the DK equations will be given
without derivation. They are given in full in Yates,
Warrick, and Myers (1985a).
THEORY
Consider Z(x) to be an isotropic second-order
stationary random function which is sampled on a
point support at p locations: .x,. .v2. . . . , xf with x a
vector in two-dimensional space. For a stationary
system, E\Z(x)] is constant, the variance is defined
and independent of location, and the spatial covari-
ance function, C(.x, xt), is defined and depends
solely on the separation distance, .t, xt. Also,
assume that a transform function <£[>'(*)] exists which
will transform Z(x) with an arbitrary frequency distri-
bution into Y(x) which has a standard normal
gaussian distribution. The transformation operates
such that the uni- r.'.d bivariate distributions are nor-
mal and jointly normal, respectively.
Disjunctive kriging
Summarizing the results of Yates, Warrick, and
Myers (1985a), the disjunctive kriging estimator is
given as
(4)
C,
I -v;/2]
»-o 1-1
)]. (5)
Z(x.) = lY(x,)}
The coefficients. Ct, are determined by orthogonality
and numerical integration as follows
(6)
where J is the total number of abscissas and weights
factors, v, and I*-,, respectively, used in the Hcrmite
integration (Abramowitz and Stegun, 1965).
Requiring an unbiased estimator which has mini-
mum variance of errors produces the following dis-
junctive kriging system of equations (see Journel and
Huijbregts, 1978; Rendu, 1980 and Yates. Warrick.
and Myers. I985a for further details)
(7)
where the series has been trunctated to K terms and
//4[K(.v.)] in Equation (7) is estimated by
(8)
where/« from Equation (4) is written as C^b* in
Equations (7) and (8). The iu's (the weights analo-
gous to simple kriging) are determined from
where .t is a vector in 2-D space which represents the
coordinates of the sample location, n is the number of
transformed sample values. Y(x,), used in the estima-
tion process, / is a function to be determined and is
expressed on the right hand side of Equation (4) as a
series of Hcrmite polynomials of order k (i.e. //*[arg])
where the /it's are the coefficients of the Hermite
expansion.
The disjunctive kriging process is based on a stan-
dard normal random function, K(.r,), which is related
to Z(x,) through a transform. The transform relation-
ship.
-------
Disjunctive kriging program Tor two dimensions
283
probability is (Journal and Huijbregts. 1978; Yates.
Warrick, and Myers. 1985a)
/"W-O1 = I - GO-,) + g(yt)
where Hk[Y(x.)} is estimated using Equations (8) and
(9). In Equation (12), g(yr) and C(yt) are the
gaussian density and cumulative distribution func-
tions, respectively.
The conditional probability density function,
Pdf*(u), is determined by taking the derivative of
Equation (12) with respect to >,.
PJf(u)
(13)
Two examples using approximately lognormal and
normal data and the disjunctive kriging programs
described in this paper are given by Yatcs, Warrick,
and Myers (1985b).
THE PROGRAMS: DISCALC.FTN AND
DISJUNCJTN
The programs DISCALC.FTN and
DISJUNC.FTN (given in Appendix A) were written
in FORTRAN 77 and run on a DEC PDF 1 1/70. The
compiled version of each program requires less than
64 K. memory for the source code, program constants
and many arrays. Program DISJUNC.FTN also
requires virtual memory for storage of a large array
(20 by 220). For machines which can access more than
64 K directly (i.e. do not need or have virtual memory)
the "virtual" statements can be replaced with dimen-
sion statements.
The program allows up to 220 samples, 25 coef-
ficients (C/s), 20 samples to be used in the estimation
process (n) and up to a 20th order polynomial fit of
the sample distribution [i.e. between Z(x) and Y(x)].
In order to improve the program's computational
efficiency, an (virtual) array (name HH) 20 by 220 was
used to store the Hk[Y(x,)] values for
* = I, 2 ..... 20 and / - I, 2 ..... 220. In this
way, only one calculation of the Hermite polynomials
associated with the samples was necessary. Also, the
program was written such that the kriging matrix was
formed only once for each estimate (for k I), was
saved and then a working matrix (raised to the appro-
priate power) was used to determine the b*'i. Using
this strategy requires only one determination of the
nearest neighbors per estimate.
Solving the DK equations (program
'DISJUNC.FTN) requires approximately (No. of
C,'s - 1) times more computer time than ordinary
kriging (with a Lagrange multiplier), holding every-
thing else constant. This does not include the
preparatory steps (Steps 1 and 2) outlined next. The
extra computer time necessary to solve the DK equa-
tions is because of the more complex calculations such
as taking powers and calculating Hermite poly-
nomials as well as having to solve the simple kriging
equations K times. Although DK requires more com-
puter time, which for a microcomputer can be in-
consequential, information about the conditional
probability distributions at the estimation site is avail-
able.
The following steps give a brief description of the
execution of the DK programs. From these programs
one can obtain an estimate of a random variable at an
unsamplcd location, the disjunctive kriging variance
and the conditional probability for up to 18 cutoff
levels using the disjunctive kriging method outlined.
Steps 1 and 2 are preliminary to the actual estimation
process and the calculations described in Step I are
made in the program DISCALC.FTN whereas those
described in Steps 2 to 7 arc made in DISJUNC.FTN.
Step I. The first step in the DK process is to
determine the coefficients. Ct, which define the trans-
form, 0(X,) in Equation (5).
(a) The calculation begins by sorting the data from
the smallest to largest value, while keeping track of
the associated .t's, (subroutine SORT). Next, an
empirical cumulative frequency distribution is
generated, from which an estimate of the probability
is obtained. The trial transformed value, K,, is ob-
tained by inverting the probability function (Function
PRBI). The estimate of the probability used in
DISCALC.FTN is
P[Z < .-.] * (i - O.S)lp
(14)
where t is the total number of Z(x.) less than or equal
to : and p is the total number of samples.
(b) Once the pairs (Z(.v,). Y(x,)] have been deter-
mined, the values of
-------
Disjunctive kriging program for two dimensions
283
probability is (Journel and Huijbregts, 1978; Yates,
Warrick, and Myers. 1985a)
J^M*.)] = i - GO-,) + s(y<)
Ht(Y(x.)yk'.
(12)
where //»[}'(*.)! is estimated using Equations (8) and
(9). In Equation (12), g(y,) and G(yt) are the
gaussian density and cumulative distribution func-
tions, respectively.
The conditional probability density function,
Pdf*(u), is determined by taking the derivative of
Equation (12) with respect to y,.
(13)
Two examples using approximately lognormal and
normal data and the disjunctive kriging programs
described in this paper are given by Yates, Warrick,
and Myers (1985b).
THE PROGRAMS: DISCALC.FTN AND
D1SJUNCFTN
The programs DISCALC.FTN and
DISJUNC.FTN (given in Appendix A) were written
in FORTRAN 77 and run on a DEC PDF 1 1/70. The
compiled version of each program requires less than
64 K memory for the source code, program constants
and many arrays. Program DISJUNC.FTN also
requires virtual memory for storage of a large array
(20 by 220). For machines which can access more than
64 K. directly (i.e. do not need or have virtual memory)
the "virtual" statements can be replaced with dimen-
sion statements.
The program allows up to 220 samples, 25 coef-
ficients (C/s). 20 samples to be used in the estimation
process (n) and up to a 20th order polynomial fit of
the sample distribution [i.e. between Z(x) and Y\x)].
In order to improve the program's computational
efficiency, an (virtual) array (name HH) 20 by 220 was
used to store the fft[Y(x,)\ values for
k = 1,2 ..... 20 and / - 1,2 ..... 220. In this
way, only one calculation of the Hermite polynomials
associated with the samples was necessary. Also, the
program was written such that the kriging matrix was
formed only once for each estimate (for k I), was
saved and then a working matrix (raised to the appro-
priate power) was used to determine the b^'t. Using
this strategy requires only one determination of the
nearest neighbors per estimate.
Solving the DK equations (program.
DISJUNC.FTN) requires approximately (No. of
Ct's 1) limes more computer time than ordinary
kriging (with a Lagrangc multiplier), holding every-
thing else constant. This does not include the
preparatory steps (Steps I and 2) outlined next. The
extra computer lime necessary to solve the DK equa-
tions is because of the more complex calculations such
as taking powers and calculating Hermite poly-
nomials as well as having to solve the simple kriging
equations K times. Although DK requires more com-
puter time, which for a microcomputer can be in-
consequential, information about the conditional
probability distributions at the estimation site is avail-
able.
The following steps give a brief description of the
execution of the DK programs. From these programs
one can obtain an estimate ofa random variable at an
unsamplcd location, the disjunctive kriging variance
and the conditional probability for up to 18 cutoff
levels using the disjunctive kriging method outlined.
Steps I and 2 are preliminary to the actual estimation
process and the calculations described in Step I are
made in the program DISCALC.FTN whereas those
described in Steps 2 to 7 are made in DISJUNC.FTN.
Step I. The first step in the DK process is to
determine the coefficients. C,. which define the trans-
form. (Y,) in Equation (5).
(a) The calculation begins by sorting the data from
the smallest to largest value, while keeping track of
the associated Vs. (subroutine SORT). Next, an
empirical cumulative frequency distribution is
generated, from which an estimate of the probability
is obtained. The trial transformed value, Y,, is ob-
tained by inverting the probability function (Function
PRBI). The estimate of the probability used in
DISCALC.FTN is
P[Z
(14)
where / is the total number of Z(.x,) less than or equal
to :. and p is the total number of samples.
(b) Once the pairs [Z(.t,). X(.r,)J have been deter-
mined, the values of (V,) in Equation (6) must be
determined. The method used by DISCALC.FTN is
to fit an nth order polynomial to the data pairs and
determine the coefficients of the polynomial via a
least-squares fit (subroutine LEAST). Function PHI
evaluates the nth order polynomial for each abscissa
in Equation (6). Other possibilities exist for deter-
mining this relationship, but the nth order polynomial
is adequate and better than linear interpolation.
(c) After the relationship is determined, the coef-
ficients, C\, are calculated by Hermite integration
(Abramowitz and Stegun, 1965) in subroutine CKI.
From the coefficients the mean and variance of the
original data can be determined (exactly as A: inf)
and arc calculated in the MAIN program and printed
out for a comparison with the actual values for the
data set.
(d) In general, a truncated series is used for the
^I^Ui)! relationship which represents an approxima-
tion in terms of transforming Y(x,) into Z(.x,).
Therefore, to make the inversion exact, new values of
Y(x,) are cakulated given the C»'s just determined.
The subroutine that inverts Equation (5) is HYZ1NV,
-------
S. R. YATES. A. W. WARRICK. and D. E. MVKRS
which uses Newton's iterative method. If the iterative
method diverges, an alternate method of bisection is
used. The values of Y(x,) (array Y7. in program) are
printed out along with other necessary information
and saved in a file which is used as an input file for the
second program DISJUNC.FTN.
Step 2. Once the transform relationship is defined,
the next step is to calculate a table of Hermite poly-
nomials for Jt = 0. I K and for all sample val-
ues Y(\,). This K x p matrix stores the values of
/MV'l.v,)) rather than recalculating them each time a
sample is used in the estimation process. The sub-
routine (in DISJUNC.FTN) that sets up the array is
HCALC and the subroutine that looks up the value in
the matrix for use in the program is H.
The following steps arc required to obtain an esti-
mate, the kriging variance and the conditional proba-
bility at an unsampled location. The sequence is
repealed for each estimate. The program is written
such that estimates are produced on a grid system
beginning at .V = XMINand Y = YMIN and incre-
menting by DX and DY. (NX)-(NY) limes.
Step 3. The first step in obtaining &n estimate is to
generate the DK matrix (subroutine MATRX). This
is done only once per estimate for k *> I by sorting
through the data and selecting the "n" nearest
samples. The spatial correlation between the points is
calculated (subroutine VARIO) and stored in array
AA.
Step 4. In order to obtain the weights, 6*, in
Equation (9). the AA matrix is raised to the tth power
by subroutine ATOK and solved by subroutine
MATINV. Because MATINV modifies the matrix
during the solution, a working matrix AWRK is used
by MATINV. For Jt = 0 the weights are equal to l//i.
Therefore the matrix solving step is skipped and ex-
ecution transfers to step 5.
Step 5. Various intermediate quantities are cal-
culated in subroutine SOLN. These include Hf[ Y(x.)]
in Equation (8) and the right-most series in Equations
(I I) and (12). After returning to the MAIN program.
these intermediate quantities (which are summed on
the data used in the estimation or on the number of
cutoff values desired) are summed on the current
value of k.
Step 6. If k < K, then k is incremented and execu-
tion is returned to step 4.
Step 7. If A- = K. the estimated value, the estima-
tion variance and the conditional probability (for
each cutoff level) arc printed out along with the spatial
coordinates.
EXAMPLE
An example using the DK. programs is given here.
Ninety-two bare soil-surface temperatures were rec-
orded at random locations in a I ha field at the Uni-
versity of Arizona's Maricopa Agricultural Center.
The sample mean, variance, skew, and kurtosis for
the data were: 63.89. 3.08.0.61. and 2.64, respectively.
The data also were tested for type of distribution by
the Kolomogorov-Smirnov (KS) test (Rao and
others, 1979; Rohlf and Sokal. 1981). The KS lest
statistic calculated from the original and log-
transformed data arc 0.121 and 0.116. respectively.
The KS critical value for 92 points at the O.I prob-
ability level is 0.084. Based on this test it was con-
cluded that the data set were neither normally nor
lognormally distributed.
The sample covariancc function (sec Journcl and
Huijbrcgts, 1978. csp. p. 194) was calculated and then
validated using (he jackknifc procedure (Vauclin,
Nielsen and Biggar. 1983). A spherical model was
fitted to the sample covariancc function with par-
ameters: 0.7''C for the nugget. S.T'C5 for the sill and
23.0m for the range. The autocorrelation function,
which is used in Equation (9), was calculated by using
Q(X, - x,) - C(.v. - .v,)/qO).
The coefficients. Ci% used to define the transform
given in Equation (5) for k = 0 to K arc: 63.891.
1.709.0.1832. -0.06818. -0.04437. -3.676 x 10°
and 5.849 x 10"'. respectively. From these Ct's,
estimates of the mean and variance for the sample
distribution can be obtained and are 63.89 C and
3.09°C; whereas the actual values for the. data set are
63.89 C and 3.08 C;. respectively.
In order to compare the DK method with ordinary
kriging (OK) both methods were applied to the data
set. From the results, a comparison based on the
average kriging variance and the sum of squares
deviation (SSQ) between the actual and estimated
values can be made.
The kriging variance was calculated by each
method at 231 locations located on & 6 x 8-m grid
system superimposed over the field. The average krig-
ing variance was determined from these values. From
the OK method, the average kriging variance was
2.202 C: and for DK was 2.098JCJ. This represents a
4.7% improvement over OK.
The SSQ between the actual and estimated value
was determined by making an estimate a: each sample
location (92 total) but not using the sample value in
the estimation (i.e. jackkniting). From the two values
the SSQ can be calculated. The results indicate that
there is an overall improvement (for these data) of
2.6% when DK (SSQ of 189.8 C:) is used instead of
OK (SSQ of 194.9 C:). Improvements of approx. 6%
have been demonstrated for another random variable
by Yates. Warrick. and Myers (I985b).
Figure la gives the estimated value of soil tem-
perature at 9 locutions using the DK and OK
methods. The number above each location (circle) is
(he estimated value using DK and the number below
is from OK. From this son of information, ihe spatial
pattern of temperatures in the field can be determined.
although one would usually use a finer grid system.
Figure Ib gives the associated estimation variance
where the number above and below each location
(circles) arc for DK and OK. respectively. Using this
information, a level of confidence in the estimated
-------
Disjunctive kriging program Tor two dimensions
285
zoc
Y 100 -
o -
(a)
661 622
66
64
4
64
6'L
.B 624
8 636
i
.7 63.5
..6 64.0
6J.4 6J./
(b)
62.6 16 13 1.
e;
6!
6!
_6!
!4 1.
1.2 2.
k i
i.5 2.
i.S 1.
S 1.J 1
« 1.9 2.
5 1.6 2
> 2.5 ?.
63.4 2.0 2.6 2
(c)
1.00 034 049
092 0.7) 06}
0.76 0.80 072
IB
J4
18
Figure I. Results of ordinary and disjunctive kriging of soil surface temperature. Numbers above and below
loealions (circles) indicate estimates by disjunctive and ordinary kriging. respectively. In a. b, and c are
temperature estimates, kriging variance, and conditional probability that temperature is above 6J.J"C.
respectively.
value can be obtained where large variance is asso-
ciated with small confidence.
Because OK produces a nonlinear estimator, val-
ues of the conditional probability that the actual but
unknown value of the temperature is above a specified
cutoff level can be calculated. Figure Ic gives an
example of this where the numbers indicate the prob-
ability that the temperature is above 62.5°C. This
information is unavailable typically when using linear
kriging methods.
Disjunctive kriging also allows one to estimate the
conditional probability density function and (he
cumulative probability distribution. These are shown
in Figure 2 for three of the points given in Figure 1.
The solid, dashed and dotted lines arc for the points:
(2. 200), (34, 200). and (18, 0), respectively. In Figure
2a the probability distribution is with respect to the
cutoff value whereas in Figure 2b the probability
density is plotted as a function of the transformed
value Y. From (his figure it is possible (o obtain an
indication of how the samples combine together to
form the estimate as well as obtaining the probability
level of an occurrence given a specified cutoff value.
CONCLUSIONS
The disjunctive kriging programs,
DISCALC.FTN and DISJUNC.FTN. produce a
nonlinear estimate of a spatially variable property as
well as the conditional probability that the value is
above an arbitrary cutoff level. Generally, the DK
estimator is better in the sense of reduced kriging
variance compared to linear kriging estimators but is
not as good as the conditional expectation (unless the
random variable is multivariate normal). Disjunctive
kriging has the advantage over the conditional expec-
tation in that only the bivariale joint probability dis-
tributions need be known.
The calculations required to estimate a random
i o
O 5
00
I 0
0 S
0 0
(0)
0
r.
(b)
0
Y
Figure 2. Conditional probability distributions. In a is cumu-
lative probability distribution as function of cutoff value. \\.
In b is probability density with respect to transformed vari-
able, IT.
variable arc basically the same as the simple kriging
method but must be performed K times per estimate.
The other differences include: a transform is necessary
and defined by the coefficients, Ct, that are appro-
priate for the data set and Hermitc polynomials up to
order AT must be calculated for each sample used in the
estimation. These steps are done at the beginning of
the problem. A third difference is that the interpola-
tion is with respect to Hk[Y(x)\ for DK whereas it is
based on Z(x) for ordinary (linear) kriging methods.
This result also is used in calculation of the con-
ditional probability. For other examples using the
DK programs see Yates, Warrick. and Myers (I985b).
-------
286
S. R. YATCS. A. W. WAWUCK. and D. E. MYEM
AcknowledgmentsAlthough part of the research described
in this article has been supported by the United Slates En-
vironmental Protection Agency through contract No.
68-7176 to Computer Sciences Corporation, it has not been
subjected to Agency review and therefore does not necess-
arily reflect the views of the Agency and no official endorse-
ment should be inferred.
REFERENCES
Abramowitz. M., and Stcgun, A.. 1965. Handbook of math-
ematical functions: Dover Publications Inc.. New York,
1046 p.
Burgess. T. M.. and Webster. R.. I980a. Optimal interp-
olation and isarilhmic mapping of soil properties. I. The
semivariogram and punctual kriging: Jour. Soil Sci., v.
31. no. 2. p. 315-331.
Burgess. T. M.. and Webster. R.. 1980b. Optimal interp-
olation and isarithmic mapping of soil properties. II.
Block kriging: Jour. Soil Sci.. v. 31. no. 2. p. 333-341.
Journel. A. G.. and Huijbregts. Ch. J.. 1978, Mining geo-
statistics: Academic Press, New York. 600 p.
Kim. Y. C. Myers. D. E.. and Knudscn. H. P.. 1977. Ad-
vanced geostatistics in ore reserve estimation and mine
planning (practitioner's guide): Report to the U.S.
Energy Research and Development Administration.
Subcontract No. 76-003-E. Phase II. 154 p.
Mathcron. G.. 1976. A simple substitute for conditional
expectation: the disjunctive kriging: Proceedings NATO
ASI "Geostal 75". D. Reidel Publishing Co.. Dordrecht.
Netherlands, p. 221-236.
Rao. P. V., Rao. P. S. C. Davidson. J. M.. and Hammond.
L. C.. 1979. Use of goodness-of-fit tests for charac-
terizing the spatial variability of soil properties: Soil Sci.
Soc. Am. Jour. v. 43. p. 274-278.
Rendu. J. M., 1980, Disjunctive kriging: A simplified theory:
Jour. Math. Geology, v. 12. no. 4. p. 306-321.
Rohlf. F. J., and Sokal. R. R.. 1981. Statistical tables (2nd
ed.): W. H. Freeman and Co.. New York. 219 p.
Vauclin. M., Vieira. S. R.. Vauchaud. G., and Nielsen.
D. R.. 1983. The use of cokriging with limited field soil
observations: Soil Sci. Soc. Am. Jour. v. 47, p. 175-184.
Vieira. S. R., Nielsen. D. R.. and Biggar, J. W.. 1981.
Spatial variability of field-measured infiltration rate: Soil
Sci. Soc. Am. Jour. v. 45. p. 1040-1048.
Warrick. A. W., Myers, D. E.. and Nielsen. D. R.. 1985.
Geostatistical methods applied to soil science: to be in
Methods of Soil Analysis. ASA.
Webster, R., and Burgess. T. M.. 1980. Optimal interp-
olation and isarithmic mapping of soil properties. III.
Changing drift and universal kriging: Jour. Soil Sci.. v.
31. no. 3. p. 505-524.
Yatcs. S. R.. Warrick. A. W.. and Myers. D. E.. I985a.
Disjunctive kriging. I. Overview of estimation and con-
ditional probability: Water Resour. Res., v. 22. no. 5. p.
615-621.
Yates. S. R.. Warrick, A. W.. and Myers. D. E.. 1985b.
Disjunctive kriging. II. Examples: Water Resour. Res., v.
22. no. 5. p. 623-630.
APPENDIX 1
Disjunctive Kriging Program} DISCALC.FTN and DISJUNC.FTN
C
C**
C"
C"
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
**»****** ************** ***************************************************
* DISJUNCTIVE KRIGING PROGRAMS FOR TWO-DIMENSIONAL,
***** cnATVftiftWttAnTAnt^nn/^n^nYf^P *******^*.
***** SPATIALLY-VARIABLE PROPERTIES **»*»
ft***************************************.************************************
Solves the disjunctive kriging equations and calculates an estimate
of a random variable, the associated estimation variance and the
conditional probability (for up to 13 cutoff levels). In order to use
these programs, values of a random variable Of Interest along with the
spatial coordinates must be obtained. Also, the correlation structure
written In terns of the varlogram must be determined.
The programs use the following devices:
Disk Input file -- currently set as unit 11
Disk output file -- currently set as unit 12
Disk plot file -- currently set as unit »3
Terminal -- currently set as unit 15
Line printer -- currently set as unit »6
These values are declared at the beginning of each propram.
<
To reduce problem execution time and storage requirements, two programs
are used for disjunctive kriging. The first program calculates the
Hermlte coefficients (I.e. CK's) and then converts the data (ZD) and
cutoff values (ZCUT) Into transformed data (YZ) and transformed cutoff
values (YCUT). Approaching the solution In this manner has the advantage
that the coefficients and transformed values only need to be calculated
once for a given data set. To run the programs, a data file must be
00001
00002
f\f\f\/\ ^
00003
f\A/\/% M
00004
D0005
00006
f\f\t\ *\ ^
00007
00008
h AA A A
00009
00010
00011
00012
D0013
00014
00015
000 IS
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
-------
Disjunctive kriging program Tor two dimensions
287
c
c
c
c
c
c
c
c
c
c
c
c
c
c
C*'
C*'
C*'
C"
C*'
C*1
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c'reated according to the "Data File Input Instructions" olven below.
The "Interactive Data Instructions For Program One" describes the run-
time information the program requires such as instructions for choosing
input/output devices and file names.
The second program uses the output from the first program directly. No
modification on the data file is necessary. The proaram asks for the
name of the file interactively, therefore there aren't any additional
"Data File Input Instructions" for the second program. To run the
program use the "Interactive Data Instructions For Program Two" given
below.
r* *»***»****************»***«*»**»***.****»»*,*****»»**.**«*»*
r******************»********************»********»***»***********»**t********
1***** ********
DATA FILE INPUT INSTRUCTIONS
****** ********
it***************************************************************************
RECORDS 1-3: FORMAT (A76/A76/A76)
TITLE(3) Three lines of title.
RECORD 4: FORMAT(5I5.2F10.3)
NZ Maximum number of data to be included in
calculating an estimate using disjunctive
kriging. that is, the number of nearest
neighbors used in estimation process.
Default: (i.e. zero or blank) is NZ 7.
NHK Number of terms in the Hertnlte expansion.
There will be NHK Hermite coefficients
calculated bv the program.
Default: NHK 7.
NPOLY Number of terms to be used in the least
squares fit to the (7.0, YZ] pairs. The
least squares polynomial is used only for
calculating the abscissa points for Hermite
integration -'to determine the CK's).
Default: NPOlY 7.
NTRM Number of terms to be used in Hermite
intearation. NTRH may be 5 or 10.
Default: NTRM 5.
NZCUT Number of conditional probabilities
to be calculated using disjunctive kriging.
For each conditional probability calculated
a ZCUT value must be provided. The ZCUT
values, if any, are input in RECORD 7.
NZCUT must he less than 19. i
/
RKAX Maximum radial search distance for a data
point to be included in the estimation
process.
Default: RMAX 10000. If linear model
for the spatial correlation function Is
used Default: RMAX RANGE (see Record 8).
The program OISJUNC.FTN will not allow a
value of RMAX > RANGE for linear model.
CUTOFF Data values larger than CUTOFF will not
be included in the analysis.
Default: CUTOFF * 10000.
D0037
00038
00039
00040
00041
D0042
D0043
D0044
0004 S
00046
00047
00048
00049
onoso
DOOS1
00052
OOOS3
00054
00055
D0056
r\f\t\f ^
D0057
00058
00059
00060
00061
00062
D0063
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074
00075
00076
00077
00078
00079
00080
0008 1
00082
D0083
00084
00085
00086
00087
00088
00089
00090
D0091
00092
00093
00094
00095
00096
00097
00098
00099
00100
00101
00102
00103
00104
DO 105
00106
00107
DO 108
-------
288
S. R. YATO, A. W. WAJUUOC, mod D. E. Mras
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
r
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
RECORD
FORMAT(3I5,3F10.3)
IBLK
NXBLK
NYBLK
umx
WIOY
BLKCOV
If: 0 Punctual disjunctive krtglng
If: 1 -- Block disjunctive krlgtno
Number of block subdivisions In the X
direction. NXBLK'NYBLK 1$ the total
number of subblocks used In discrete
approximations to the Integrals of
the correlation (or covarlance) over
the block. If IRLK equals 1 then
Default: NXBLK 5.
Number of block subdivisions In the V
direction. If IBLK equals 1 then
Default: NYBU 5.
Length of the block In the X direction.
Total block area Is UIDX*UIDY.
Length of the block In the Y direction.
The Inner block covarlance (I.e. the
covarlance value resulttna from the double
integral of the covartance over the block.
where each vector Independently describes
the Interior of the block). If a value Is
given BLKCOV will NOT be calculated.
Default: RLKCOV will be calculated.
RECORD
FORHAT(2I5,4F10.3)
NX
NY
XHIN
YHIN
DX
OY
OPTIONAL RECORDS 7a.7b:
(INCLUDE IF NZCUT > 0)
ZCUT(I)
Number of X coordinates In the grid system
of estimates. There will be NX«NY estimates
total. Default: NX 10.
Number of Y coordinates In the grid system
of estimates. There will be NX'NY estimates
total. Default: NY 10.
The minimum value of X on the orid system.
Default: XHIN » 0.
The minimum value of Y on the grid system.
Default: YHIN -0. " *
The distance between estimates in the X
direction on the grid system. The length
of the grid system In the X direction
is (NX-l)'DX.
Default: DX (Max. X in data file )/9.
The distance between estimates in the Y
direction on the grid system. The length
of the .grid system in the Y direction Is
(NY-1)«DY
Default: DY » (Hax. Y in data file 1/9.
FORHAT(8F10.3)
The cutoff values (in terms of the original
variable). A conditional probability will
be calculated for each cutoff value. A
of 18 values are allowed.
00109
00110
D0111
00112
00113
00114
00115
00116
00117
D0118
00119
00120
D01?l
00122
00123
00124
DO 125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
D0137
D0138
DO 139
00140
DOH1
0014?
00143
00144
00145
00146
00147
00148
00140
00150
00151
00152
00153
00154
001S5
00157
00153
D0159
00160
00161
00162
00163
00164
00165
D0166
00167
00168
D0169
00170
00171
0017?
00173
00174
00175
00176
00177
00178
00179
00180
00131
-------
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
r
r
C
C
C
C
f
%
C
C
C
C
^
I
C
C
C
C
C
C
C
C
C
C
C
C
C
RECORD 8:
Disjunctive kriging program for two dimensions
FORMAT(I5,3F10.3)
28V
NODE
CO)
C(3)
RECORD
FORHAT(315)
I OX
IDY
IDZ
The varloqram model number where:
-1
0
1
exponential:
[C(D * C(2)«E!(P(-R/C(3)]
linear (with a sill):
[C(l) * C(2)*R/C(3)]
spherical:
* C(2)|1.5'R/C(3) -
The correlation function is calculated by
correlation I. - varloqran/slll
Note: In two dimensions a linear rodel with a
«« sill Is an Invalid type of variooram
(I.e. it Is not of the conditionally
positive definite type). Therefore, thj
prooram will not allow RMAX to be greater
than RANGE If a linear model Is chosen.
NUGGET value of the varioqram.
SILL value minus the NUGGET value of the
varioaram.
RANGE of the varioqram. For the exponential
variogran this is the constant that the
distance Is divided by.
00182
D01»3
00184
00185
00186
00187
D01S8
00189
00190
00191
00192
00193
00194
D0195
00196
00197
00198
00199
00200
00201
P020?
00203
00204
00?05
00206
00207
The colunn In the data file that contains
the X coordinate data. Note: an Index (or
sample number) Is always assumed to be
in the first column of the data file and
the Index column Is not counted.
Therefore, IDX 1, 2, or 3.
The column In the data file that contains
the Y coordinate data. Note: an 1nde» (or
sample number) Is always assumed to be
In the first column of the data file and
the Index column is not counted.
Therefore, IDY 1, 2. or 3.
The colu.im In the data file that contains
the property to be kriged. Note: an index
(or sample number) is always assumed to be
in the first column of the data file and
the Index column is not counted.
Therefore, IOZ 1, 2, or 3.
RECORD 10: FORHAT(4X.A6)
A na«e which describes the oroperty to be
kriqed.
11:
FKT
Format specification for the data.
should be of the form;
This
00209
D0210
00211
00212
00213
00214
00215
00216
00217
00213
0021"
00220
00221
00222
00223
00224
00225
0022R
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00233
00239
00240
00241
00242
002 43
00244
00245
00246
00247
0024S
00249
00250
00251
00252
OP253
00254
-------
290
S. R. YATES, A. W. WAKRICK. and D. E. Mrau
c
r
r
r
%
c
C 3£CO»0'
c
r
r
r
C
c
c
c
c
c
c
00255
(IS,2F10.3,5»,F10.3) D0256
00257
DO? 58
00259
12- : FORMAT IS GIVEN BY FMT (RECORD 11) 00260
D0?61
DATA VALUES The data set should have a location nt.nber 0026?
(or lnde«). » coordinates, Y coordinates P0263
and the property to be Vrlgeri. Data Is 00254
read until an end-of-ftle marker is found. 00265
Therefore, there should not. be any blank 00266
lines In this part of the data file. 00267
It Is assumed that the location number Is 00268
in the first column. D0269
00270
00271
C****************************************************************************** 00? 7?
£******************«************************************* 00? 7 3
£*****
c *
£***
00274
£******************************* 00?/8
C RECORD
C
C
C
C
C
C
C RECORD
C
C
C
c
c
c
c
c
c
C RECORD
C
C
C
C
C
C
C
C
C RECORD
r
r
c
c
C
C RECORD
C
C
C
C
C
C
1: 281
FILE Input file name. File nane must contain 002?2
less than 31 characters. This Is t»>e 00283
n»e«e of the data Hie created usino 0028*
DATA FILE INPUT INSTRUCTIONS' Qtveo above. 00285
00285
00287
2: FILE Output file name. File name i»ui» contain D0283
less than 31 characters. This file is the 00239
input file for DISJUNC.FTN. Only those D0290
data and steering pjrarr-eters needed by the 00291
program OISJUNC.FTN are printed out to this 00292
file. This file is used 'as is' by the 00293
proorao -OISJU1C.FTK-. D0?9«
00295
00296
00297
3: 4NSW Output device nurter where the outout/ 00298
sumary will be sent. Two choices are D0299
allowed: 00300
M301
P -- Output sent to line printer 0030?
T -- output sent to terminal D0303
P030*
00305
******** M307
>«* DP310
P03U
1: 00315
FlLf Input 'ile name. File name must Contain 00316
less than 31 characters. This is t»e P0317
OUTPUT filo Created using OlSCALC.fTN. D0319
D0319
00320
2a: *NSW Output device designation. Choices 00321
allowed: 00322
00323
0 - disk file 00324
T - terminal 00325
P - 1 ine printer 00325
00327
-------
Disjunctive kriginj program for Iwo dimensions
291
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C"
C
C
C"
C"
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
C
OPTIONAL RECORD 2b: (Include only If ANSW.EQ.'O' In RECORD 2a)
FILE Output file name. File name must contain
less than 31 characters.
RECORD 3: FILE Plot file name. Produces a file that
can be used to make contour maps. etc.
The file contains x, y, estimates,
estimation variance and conditional
probabilities. File name must contain
less than 31 characters. If a plot file
is NOT wanted type a return.
A**************************************************************************
k . ********
MISCELLANEOUS INFORMATION «»«...*
***** ***
A***************************************************************************
A*************************************************************************
IMPORTANT ARRAYS:
Contains the following Information:
INX(I) - Index or sairple number for the Ith data record.
XD(I) - x coordinate datum for the Ith data record.
YO(I) - Y coordinate datum for the Ith data record.
ZD(I) - Property of Interest that Is to be estimated.
YZ(1) - Transform of ZD used in disjunctive kriging equations.
PR(I) - Probability of the Ith datum.
CK(J) - Coefficients of the Hermite expansion for the transform.
ZCUT(K) - Cutoff values in terms of the sampled variable (I.e. ZD),
YCUT(IC) - Cutoff values in terns of the transform var. (I.e. YZ).
PROB(r) - Conditional probability for cutoff level number X.
CS(L) - Coefficients for the least squares relationship.
C(M) - Variogram model coefficients.
AA(N»1,N+1) . Coefficient matrix of the kriqtng equations.
Right hand side of equations is in AA(i.N»l), i
BB(N»1) - Solution vector (i.e. the krtolng weights!.
ILOC(X) - The location numbers for the nearest nelohbors.
RLOC(N) - The radial distance between each nearest neighbor and
the estimation site.
GOX(N) - The correlation between each nearest neighbor and
the estimation site.
1....N*!
WERE:
I
J
K
L
M
N
1,2,....*INO(number of Samples, 220)
1.2.....MIHO(number of Heraite polynomials. 25)
l,2,...,MINO(number of cutoff levels, 18)
1,2,...,MlNO(number of least souares polynomials. 20)
1.'.3
1.2 MINO(number of nearest neighbors, 20)
STEPS NECESSARY TO USE DISJUNCTIVE KRIGING PROGRAMS:
1) Obtain the data set. Must Include an X. Y and Z (property of
interest)
2) Determine the varlogram. (Not Included In these programs).
3) Create data file using "Data File Input Instructions"
4) Run OISCALC.FTN
5) Answer Questions prompted by program using
"Interactive data instructions for program one", aenei'her that
output from this program Is Inout for OISJU'iC.FTN
6) Run OISJUNC.FTN
7) Answer questions pron'pted by arogram usino
"Interactive data instructions for program two". 3ene<-ber to
use the output froi" OISCALC.FTH as Input to this
00328
00329
D0330
00331
D0332
00333
00334
00335
003 36
D0337
00338
00339
00340
DO 341
00342
00343
00344
00345
00346
D0347
0034?
D0349
00350
00351
03352
D0353
00354
D0355
00356
00357
D0358
D0359
00360
00361
00362
00363
00364
H0365
00366
00367
00368
D0369
00370
"0371
00372
00373
00374
00375
P0376
00377
00378
00379
00380
DO 381
00382
00383
00384
00385
00386
D0387
00388
D0389
00390
00391
00392
00393
0"394
00395
P0396
00397
00398
003<>9
00400
00401
00402
-------
292 S. R. YATQ. A. W. WAMUCK. and D. E MYEU
C
< «»»»»»»«»»«»»*«««**«»»» HA 00 3
£*».»«* «**«««» PA004
C«*"«" PROGRAM ONE: OISCALC.FTN ........ MA005
£»..*.» «.....« MA006
£********************************************** PA007
C KA008
C PURPOSE: W009
C Calculates Hernlte coefficients (I.e. Ck's). MA010
C Converts the sample (» ' MA058
REAO(IT,700) FMT «059
OPENtUNlT'IN.FILE'FKT.STATUS.'OLD1 .READONLY) MA060
WRITE(IT,701)' GIVE OUTPUT FILE NAHE (input for krigina) >» ' W061
REAO(IT,700) FHT MA062
OPEN' UNIT- 10. FILE-FMT, STATUS' 'UNKNOWN', CARR I AGECONTROL' 'LIST') MA063
1 WRITEf IT, 850) >-A064
READ(IT,R55) ANSW fA065
IF(ANSW.EO.'P') THEN VA067
KOUT'IL ^068
ELSE IF(ANSW.EO.'T') THEN MA069
KOUTMT MA070
ELSE
WRITE(IT,865)
GOTO 1 MA073
END IF MAO 74
-------
Disjunctive kriging program for Iwo dimensions
293
700 FORWU(ASO)
701 FORMT(////,A50,S)
MA075
MA076
300 FORMAT(A76)
305 FORMAT(4(4X,A6)1
310 FORMAT(5I5,6F10.3)
315 FORMAT(315,6F10.3)
S20 FORMAT(2I5,7F10.3)
825 FORKAT(IS,7F10.3)
S30 FORMAr(8F10.6)
335 FORMAT(A30)
3^0 FORMAT(1P5E15.7)
345 FORMAT(I5,?F10.3,F10.5,F1?.S)
850 FPRMAT(////' GIVE OUTPUT DEVICE (FOR THE OUTPUT/SUMMARY) \/,T20
#,'P . send to line printer1,/,T?0,'T --- send to terminal'.T50
«.' >» ',S)
355 FORMAT(AJ)
?60 FOKMATf///// ERROR: Array overflow. Execution will cease.',/
*.A5,' must he less than or eaual to ',13.'.'.///)
P65 FORMATf////1 ERROR: NOT AN ALLOWED INPUT. (Try aqalnV.//)
370 FORMAT(F10.4)
read input parameters
REAO(IN.SOO) (TITLE(n,I»1.3)
BEAD(IN,810) NZ.NHK.NPOLY.NTfJM.NZaiT.RMAX,CUTOFF
READ(IN.815) IBLK.NXRIK.NYRLK.WIDX.WIDY.BIKCOV
REAP(IN,820) NX.NY.XHIN.YMIM.nx.DY
default parameters (steering)
IF(NZ.IE.O) NZ « 7
IFJNHK.LE.O) NHK ' 7
IF(NPOLY.LE.O) NPOLY » 7
IF(KTRM.LE.O) NTR." » 5
I^RMAX.LE.O.) RHAX « 10000.
iF{CL'TOFF.LE.O.) CUTOFF 10000.
IFflZCUT.GT.NB.OR.NZCUT.LT.O) GPTP 10
IF(S2CUT.GT.O) REAO(IN.?30) (ZCUT( I ) ,1«1.NZCUT)
5EAO(IN.S25) W>E.(Ctn.!'1.3)
RE4PCN.P15) IXO.IYD.IZP
F"T
«r'AX«-99Q99.
y«»4<». 99999.
;r-(ar.-P( i^n). GT. CUTOFF) GOTO 5
IXC"
-------
294 S. R. YATO. A. W. WAUUCK. ind D. E. MYHU
IF(NPOLY.GT.NE) WRITE( IT,860) 'NPPLY1.IFLAG
IF(NZCUT.GT.NB) IFLAG-1
IF(NZCUT.GT.NR) WRITE(IT.P60) 'NZCUT'.KB HA150
IF(IFlAf..GT.O) STOP M*151
C HA152
C «A1S3
C calculate the sample mean and variance MA154
NOB'I MA155
MEAN « MEAN/FLOAT(NOB) MAISS
VAR * 0. MA157
00 13 I»1,NOR HA158
13 VAR VAR + (ZO(l)-HEAN)«(ZD(l)-MEAN) PA159
VAR « VAR/FLOAT(NOB) HA160
C MA161
C MA162
C determine the Hermite coefficients HA163
CALL SORT(NA,N08,INX.XO,YO.ZD,YZ.PR} MA164
CALL LEAST(NOB.NA,NC,NE.ZO.PR.AA.NPOLY.CS) HA165
CALL CK1(NE.NTRM.NHK.CK.NPOLY.CS) HA165
C HA167
C MA168
C Invert PHI(Y) relationship for the data HA169
CALL HYZINVfNOB.INX.ZO.YZ.NHK.CK) MA170
C MA171
C MA172
C invert PHI(Y) relationshio for the Ycut values ."A173
IF(NZCUT.EO.O) GOTO 20 HA174
00 15 I-l.NZCUT MA175
15 YCUT(I) « (ZCUT(I)-MEAN)/VAR HA176
CALL HYZINI(NZCUT.INX.ZCUT.YCUT.NHK.CK) HA177
C MAI 78
C MAI 79
C calculate the sample variance fron- the coefficients HA180
20 SUH = CK(2)*CK(2) MA181
FACT « 1. MA132
00 25 I«3.tlHX MA183
FACT ' FACT*FLOAT(M) MA134
25 SU.H » SUM + FACT«CMI)*CK(I) MA185
C MA186
C MA187
C default parameters for matrix NA138
IF(OX.Lc.O) DX XMAX/9. MA189
IF(DY.LE.O) PY « YWX/9. ' MA190
IF(NX.LE.O) NX » 10 MA191
IF(NY.Lf.O) IIY 10 >'A192
IF(1BLK.E0.1.AMO.»XBLK.LE.O) NXBLK»5 MA193
IF(IRLK.FO.I.ANO.HYBLK.LE.O) NYBLK»S «Ai9«
C MA195
C fA196
C write data file for disjunctive Vrioino MA197
WRITEtIO.800) (TITLE(I).1-1.3) MA198
W«ITE(10.815) NZ.NHK.NZCUT.RKAX.CUTOFF MAigg
WRITE( 10.815) IBLK.NXBLK.NYBLK.UinX.WIOY.RLKCOY >'A200
ws[TE(IO.P20) NX.NY.X-IM.YMIK.OX.OY MA201
IF(NZCUT.flE.O)WS|TE(!0,ft30) (ZCUT( D.I-l.HZCUT) M«202
IF(NZCUr.ME.O)WR|IE;iO,830) (YCUT(I) .1-1.NZCUT) MA2Q3
WRITE!I".P<1) (CK(!),I»1.NHK) Ma204
W»ITF(10.8?5) MOOE.(C(I).I-1.3) MA205
WRITE(10,805) NAM MA20*
C MA207
00 30 I-1.N08 MA208
30 WRITE(10,8«5) lNX(I).XO(I).Yn(I),ZO(I).YZ(I) MA209
C MA210
C MA211
C print out summary of calculation MA21?
WRITE(KOUT,900) HA213
00 35 1-1.3 KA214
35 URITE(KO:lT,905) TITLE(I) MA215
VRITECKO-jT.giO) MA216
C MA217
NLIN 31 HA218
URITE(ltOUT.915) NAM,NOB.NZ.NTRM.RHAX.CUTOFF MA219
IF(IBLK.GT.O) THEN MA220
-------
Disjunctive kriging program for two dimensions
295
BI ' CALCULATE'
IF(BUCOV.GT.O.) WRITE(8I,870) 8LKCOV
VRITE(KOUT,920) BI.NXBLK.NYRLK.WIDX.VIOY
NLIN NLIN+4
ELSE
END IF
IF(MOOE.LT.O) BI ' EXPONENT*
IF(MOOE.EO.O) BI « LINEAR-
IF(MCOE.GT.O) BI ' SPHERICAL1
VRITE(KOUT.925) BI.CdKC(2).C{3)
WRITE(iCOllT,930) NAM,NAM
print data
00 40 I'l.NOS/2+1
NLIN NLIN + 1
IF(NLIN.GE.58) WRITE(KOUT,855) 'I1
IF(NLIN.G£.5P) NLIN 0
II » NOB/2 * I * 1
IF(INX(1).EO.O) INX(I) - I
IF(INX(II).EO.O) INX(II) II
IF(H.GT.NOB) INX(II) « 0
40 WRITE(KOUT,935) l!tX(I),XD(I) ,YO( I).ZO(I) ,YZ(I) ,INX( II) ,XO( II)
.YD(II).ZD(II),YZdn
NLIN » NLIN + 7
IF(NLIN.GE.SO) W31TE(ICOUT,955) 'I1
IF(NLIN.GE.SO) NLIM * 0
WPITEOCOUT.940) M£AJi.CK(l),YAI»,SUM
orint coefficients
ITW'MAXOfNHK.NPOLY)
"LIN » NLIN » 7
IFfXLIN.GE.SP) WSfK(iCOi|T,
IF(MN.LE.NPOLY.ANQ.MN.LE.NHK) WRITE(Kn«lT.950)
IF(MN.LE.NPOLY.ftHO.^.GT.NHK) V«UTE(
-------
2% S. R. YAIB, A. W. WAJWiric. and D. E. Mraa
945 FORMAT(/,11X,'COEFFICIENTS:'/.11X,13(1H«).//.11X,'LEAST SQUARES:'. ' f'A294
&18X.-HERHITE POLYNOMIALS:'./.HX.'tcun. distriu)',/) "A295"
950 FORMAT(11X.F14.6,22X.F14.6) HA296
955 FORHAT(11X.F14.6) W7
960 FORMAT(47X,F14.6) HA298
C MA299
END HA300
C H8001
ft********************************************************** »»*»*»»****** MBOC2
ft******************************************************************** ^8003
£*«.. *»«« M8004
C*"««* PROGRAM TWO: OISJUfC.FTN * KRC05
£*»»*« »»»* M8006
Q**4*********4*******«*********«******4*******************»«tA.**************** M-B007
c HBOOS
C PURPOSE: H8009
C Locates the NZ nearest neighbors to be used in the estimation process. MB010
C Calculates the spatial correlations and constructs coefficient matrix HB011
C Solves matrix equations for disjunctive krtging weights. PB012
C Calculates estimates, error variances and conditional probabilities MB013
C Writes out results Into a user specified device. MBOH
C Creates a plot file if specified hy user. MB015
C ."8016
C N8017
C DEVICES REQUIRED: KBOlfl
C Unit '1 -- input file disk device (Output from PISCALC.FT'I) PE019
C Unit *2 -- output file disk device (Used only 1f disk output specified) KR020
C Unit »3 -- olot file disk device (Used only if plot file specified) M«021
C Unit «5 -- terminal W022
C Unit «6 -- line printer (Required only If user specifies surwnary "3023
C Of results to be sent to line printer. MBO?3
C This is an Interactive input). "8025
C K8026
C MB027
C SUBROUTINE AND FUNCTION CALLS: W23
C 1) SOLN 2) HJTRX 3) ATOK 4) HATlNV ''8029
C 5) AVECCR 6) BLKPTS 7) VARIO 3) HCALC M8030
C 9) H J, HI (entry point) 10) HK 11) PRR KB031
C MB032
C MB033
C M8034
C WITTEN BY: MB035
C Scott R. Yates, Computer Sciences Coro. M8036
C R.S. Kerr Environ. Res. Lab., Ada, OX, 74820 MB037
C V3033
£*******»»********»*****»*******«**»***««***»****» MB039
^********************»«*****ft*ft**4**»******«*******************«*****«* ^gQ^Q
CHARACTER T!TLE(3)*76,FMT»30,MI(3)*n,SI*6,.NAH«fi,/iNSW«l,AZCUT«15, V8042
< AZEXP(18)*4 Mpo43
DOUBLE PRECISION AA(21,21) s'3044
REAL XO{220),YP(220).Zn(??0).YZ(220).C(3).CM?V..PROB{lS). M8045
-------
Disjunctive kriging program for cwo dimensions
297
C
c
C
open files on POP 11/70
WRITE(IT,701)' RIVE INPUT FILE NAME »> '
REAO(IT,700) FM.T
OPEN(UNIT»IN,FILE'FMT.STATUS-'OLO',READONLY)
1 WJITEUT.703)
REAOUT.702) »NSW
IF(ANSW.EQ.'T' .OR. ANSW.EO.'t') THEN
IO-IT
ELSE IF(Af;SU.EO.'P' .OR. ANSW.EQ.'p1) THEN
IO-IL
ELSE IF(ANSW.EO.'0' .OR. ANSW.EO.'d1) THEN
WRITE(IT,70l)' GIVE OUTPUT FILE NAME >» '
READ(IT,700) FMT
OPENJUNIT-ID.FILE'FMT,STATUS-'UNKNOWN')
IO»IO
ELSE
WRITE(IT,«)' INCORRECT VALUE GIVEN MUST BE [T], [P] OR [0]'
GOTO 1
END IF
WRITE(IT,701)' GIVE PLOT FILE NAME (return-none) »> '
READ(IT.700) FHT
IF(FMT.EO.' ') GOTO 5
OPEN(UNIT-IP,FILE-FMT,STATUS-'UNKNOWN')
700 FORMAT(A30)
701 FORKAT(////.A50,S)
702 FORMAT(Al)
703 FPRKAT(////.' GIVE THE OUTPUT nEVICE: './.T20,
*' 0 -- create a disk file './.TZO,' P -- send to line printer '
«./,T20.' T -- send to terminal',T46.'»> '$)
BOO FORMAT(A76)
805 FORMAT(3I5,7F10.3)
810 FORMAT(2I5,7F10.3)
815 FORMAT(8F10.3)
820 FORMAT{1P5E15.8)
P25 FORMAT(I5,7F10.3)
830 FORMAT(4X,A6)
835 FORMAT(15,3F10.3.F1?.8)
840 FORMATf////,' ERROR: array overflow -- ',A5,' Is oreater than '
* ,I3.//)
845 FORMAT(F6.3)
read steering parameters
5 REAO(IN,800) (TITLE(I).1-1,3)
REAOfIN.805) NZ,NHK.NZCUT.RMAX.CUTOFF
REACI(IN,805) IBLK.NXeLK.NYBLK.UlDX.UIDY.BLKCOV
REAO(IN,8IO) NX.NY.XMIH.YMIN.OX.OY
IF(NZCUT.Kr.O) REAO(IN.815) (ZCUT(I).l»l.NZCUT)
IF(NZCUT.NE.O) REAn(IN,815) (YCUT( D.I-1.N7CUT)
REAO(IN.P?0) (CK(I).I'l.HHK)
REAO(IN,825) HOOE,(C(I).1-1.3)
IF(>«OD£.EO.O.AND.RMAX.GT.C(3)) RMAX C(3)
REAO(IK,830) NAM
transforn varfogram into the correlation function
CVARK . C(U»C(?)
C(1)'C(1)/CVARK
C(2)'C(?)/CVARK
read data from file
1*0
MEAN - n.
15 I'I+1
20 RE«D( IN,835.^0-25) INX( I ).XO( I),Yf)( I) ,ZO( I) ,YZ( I)
MB066
M8067
MB068
MB069
MB070
M8071
MR072
MB073
M8074
MRQ75
MB076
MR077
MB07«
M8079
MB080
M8081
MB082
MB083
MB084
MB085
K8086
MB087
MB088
MB089
MB090
MR091
MB092
MB093
M8094
MB095
MB096
MB097
M8098
MB099
M3100
M8101
MB 102
MB103
MH104
MP105
,«8106
^6107
MB 108
MB 109
MB 110
M8111
MB112
MB113
MB114
MB11S
M8116
M8117
MB118
MB 119
MB120
M8121
MB122
M8123
MB124
MR125
MB 126
MB127
MB128
M8129
MB 130
MB131
MB132
MB 133
KB 134
MB 135
MB 136
MB137
-------
298
S. R. YATTS. A. W. WMUUCK. and D. E. Mvws
C
C
C
IF{ZD(I).GT.CUTOFF) GOTO 20
MEAN MEAN + ZD(l)
GO TO 15
25 HOB-l-1
check for array overflow
IFtAC,«0
IF(NOB.GT.NA) IFLAG-1
IF(NOR.GT.NA) WRITE(IT,840)' NOB ',NA
IF(NHK.GT.ND) IFLAG'l
IF(NWC.GT.NO) WRITE(IT.840)' NHK '.NO
IF(NZ.GT.NE) IFLAG-1
IF(HZ.GT.NE) WRITE(IT,840)' NZ ',NE
IF(NZCUT.GT.NB)IFLAG«1
IF(NZC.UT.GT.NB)URITE(1T.840)'NZCUT',NB
IF(IFLAG.EO.l) STOP
calculate the sample statistics
MEAN . MEAN/FLOAT(N08)
VAR « 0.
00 27 I-1.N08
27 VAR V6R + (ZD(I)-MEAN)'(ZD(I)-MFAN)
VAR - VAR/FLOAT(NOB)
SUM . CK(2)«CK(2)
FACT 1.
00 30 I-3.NHK
FACT « FACT«FLOAT(I-1)
30 SUM « SUM + FACT««( I )«(!)
calculate inner block covaHanre. Function SECNOS 1s
a system routine which returns the number of seconds
from start of day and used to "randomUe" the seed
IF(IBLK.EO.O) BUCOV » SUM
IF(IBUC.EQ.O) GOTO 35
ISEED IF!X(SECNOS(0.0))
!F(BUCOV.GT.O.O) GOTO 35
CALL AVECOR(ND.8LKCOV.C.MOOE.NHX,CK)
35 CONTINUE
----- write out summary of input data -----
WRITE (10, 900)
IF(IBLK.FO.O) WRITEU0.905)
IF(tBLK.NE.O) WRITE(IO,°10)
00 40 1-1.3
40 WfiITE(I0.915) TITLE(I)
UR[TE( 10.920)
NUN 31
URITE(I0.925) NAM. NQB. NI.RMAX. CUTOFF
IF(IBLK.C-T.O) WRITE{10.930) BLKCOW .NXRLK .NYBLK .W10X
I«~(;BLK.GT.O) KLIN - m.iN*4
URITE(I0.935) MI(MOOE»2).C(1).C(2),C(3)
WH|TE(I0.940) NAU.NAM
----- print out data .....
00 45 I«1.N08/2»1 ,
NLIN NLIN * 1
IF(NLIN.GE.58) WR1TE(I0.702) '!'
IF(NLIN.GE.58) NLIN 0
II HOB/2 » I + 1
IF(INX(I).EO.O) IHX(I) 1
IF(INX(Il).EO.O)INXUI) II
IF(II.GT.NOB) INX(II)-0
45 WRITEII0.945) 1NX{ I),XO( I),YO( I).ZO(1),Y2( I) ,INX( I I).XO{ 1 1)
I .YO{II),ZO(II).YZ(I1)
NLIN NLIN » 7
IF(NLIN.r.E.SO) WRITE(I0.702) 'I1
IF(NLIN.GE.SO) NLIN - 0
URITEU0.950) MEAN, C<(1). VAR. SUM
MB138
MB139-
MB140
M3141
M8142
MB144
M8145
MB146
M8147
H8148
M8149
MB150
K-B151
MB152
MB153
MB 1 54
M8155
MB156
MB15/
PB153
H8159
MP160
MB161
MB162
MB163
H8164
M8165
MB166
M8167
MB 168
MB 1 69
MB 170
M8171
MB 172
M8173
MB174
MB175
MB 176
MB177
MB 178
MB179
MB 180
MB181
MB182
MB183
MB184
MB185
MBlSfi
MB187
MB 188
MB189
MB190
MB191
M8193
MB 194
M8195
MB196
M8197
MB 198
KB 199
MB200
MB201
M8202
M8203
M8204
MB 20 5
MB206
M3207
MB208
MP209
MB? 10
MB211
-------
Disjunctive kriging program for two dimensions
299
50
55
----- print out coefficients
NLIN « KLIN * 5
IF(NLIN.GE.SO) URITE(I0.702) 'I1
IF(*LIN.r,E.50) NLIN 0
VRITE{10,955)
00 50 I-l.NHK/2+1
II NHK/? +1+1
K - 1-1
KK II-l
K'.IN « NLIN + 1
WRITE(IO,960) K,«(I),r.K.CK(II)
NLIN NLIN + 7
IF(NLIN.GE.S8) WRITE(IO,702) 'I1
IF(NLIN.GE.58) NLIN 0
81 « ' Zcut '
00 55 I'l.NZCUT
WRITE(AZCUT,820) ZCUT(I)
REAO(AZCUT(1:11),845) ZCUT(I)
AZEXP(I) - AZCUT(12:15)
CONTINUE
IV-MINO(9. NZCUT)
JV«HAXO(NZCUT-IV.l)
KV-MAXO(l,(IV8-50)/2)
LV»MAXO(l.(IY*8-16)/2)
WRITE( 10,965} (Bl, 1-1. HINDU .NZCUT))
WRITEU0.970) (BI.l.l,HINO(l,NZCUT))
V
-------
300
S. R. YATH. A. W. WAXWCX. and D. E. MYERS
C these values are set to either 0 or 1 M8285
IF(NZCUT.EQ.O) GOTO 95 M8286
00 90 M'l.NZCUT M8287
TMP » YCUT(K) MB288
PROB(M) « 1.0 - PRB(TMP) + EXP(-TMP«TMP/2.)*PROB(H)/?I2 H8289
IF(PPOB{M).GT.1.0) PROB(M) . 1.0 M8290
90 IF(PROB(M.).LT.O.O) PROB(M) « 0.0 MB291
C MB292
95 VARK RLKCOV-VARK MB293
C M8?94
C orint out results MB295
IF(NLIN.GE.S8) WRITE(I0.702) 'I1 MB296
IF(NLIN.GE.58) NLIN « 0 MR297
NUN » NLIN » 1 MR298
IF(IP.E0.3)WR|TE(I°,990) XO.YO.EST,VARK.(PR08(H),M»1.HZCUT) M8299
VRITE(IO,995) N1.XO,YC,EST,VARK,(PROB(M),M.1.N/.C'JT) MB300
65 YO " YO + OY MB301
60 XO " XO + DX MB302
STOP MB303
C MB 304
C end of problem MP305
900 FPRr'AT(lHl,10X,32(lH*)/ll*,lH*.SOX.1H«/11X,IH«,22X.'DISJUNCTIVE KR MB306
IGING ON A GRIO MATRIX'.22X.1H*./.11X.1H«.18X.'FOR ESTIMATING SPAT M8307
IALLLY DEPENDENT PROPERTIES',16X.1H«/11X.1H«.BOX,1H-) MB308
905 FORMAT(11X.1H«,32X.'POINT ESTIMATION',3?X,1H«,/.I IX,1H«,SOX,1H«) M8309
910 FORMATJ1IX,1H-.32X.'BLOCK ESTIMATION'.32X.1H',/,!IX,1H«,SOX,1H*) MB310
915 FORMAT(11X.1H*.2X.A76,2X.1H«) M8311
920 FORW(llX,lH«,2X.76X,2X,lH«./.llX,fl2(lH«)) HB312
925 FORW(//11X.'INPUT PARAMETERS'/11X.16(1H-)/ M8313
*11X,'NUMBER OF '.A6."'s INPUT '.HO./. MB314
»11X,'MAXIMUM NUMPER OF NEAREST NEIGHBORS '.HO,/. M.B315
«1 IX.'MAXIMUM RADIUS '.F10.4./, HB316
*11X,'MAXIMUM ALLOWED DATA VALUE '.Fl!).4./) M8317
930 FORMAT( MB318
11X.-DISPERSION COVARIANCE WITHIN BLOCK '.F10.4./ M8319
*11X.'NUMBER OF CELLS',34(1H.),5X,- NX « '.HO.SX.'HY . '.HO./ MB320
«11X,'BLOCK SIZE'.5X,34{1H.),5X,' DX » ',F10.3.5X,'OY '.F10.3,/) MB3?1
935 FORMAT{11X,'CORRELATION FUNCTION MODEL TYPE '.All MR322
«. /.I IX. 'NUGGET '.F10.4,/ MR3?3
.IIX.'SILL MINUS NUGGET '.F10.4./ MB324
C.I IX,'RANGE '.F10.4) MB325
940 FOP.MAT(//11X.-OBSERVED DATA- ,/HX,13( 1H»)/.7X.2( 'OBS. NO.'.7X,lHX. MB326
«HX.lHY,6X.A6,8X,2HYz,15X)) M8327
945 FORMAT(2(I12,2F12.3,F12.3.F12.6,5X)) MBJ28
950 FOR"AT(//,11X,'OATA SET:.23X.-HERNITE POLYNOMIALS: './.15X, MB329
'MEAN -,F12.5,9X,-MEAN . -.F12.5./.15X. M8330
'VARIANCE '.F12.5.9X.'VARIANCE '.F12.5.//) MB331
955 FORyATf/.llX.'HERMITE POLYNOMIAL COEFFICIENTS:'/.11X.32(IH«).//. MB332
960 FORMAT(11X J5.1PE15!6.6X.I5.1PE15.6) MB334
965 FORMAT(///.11X.-DISJUNCTIVE KRIGING ESTIMATES' ,19X,: ,(' '),50(1H»),//.59X.( '), MB337
«'Value(s) of ,A6) MB338
975 FORMAT(/.59X,:,(:,'(',F6.3,')'),/,59X,<1V>(:,'( '.A4,')'),: MB339
,/.5'JX,(:,'(' ,F6.3,')'},/, 59X,(:,'( I,A4,')')) ' MB340
930 FnRMATdOX.'l'.SX.'X'./X.'Y'.SX.'EST'.fiX.-VARK'.J) 1*3341
985 FORMAT('»',15X,('» *')) MB342
990 FORI1AT(2F8.2.2F9.3.4X,2(9(1X,F6.3.1X),/.38X)) MB343
995 FORMAT(9X.I2.?FS.2,?F9.3.14X,2(9(1X.F6.3.1X)./.59X)) M8344
C M8345
END «B346
C S0001
C*************************<************»*****«*******»*t****»*«»***«»*ftt»**»*«* SOOO2
(^..^t. jQoo3
C««*««*» «.»*««« SQ004
c....... SUBROUTINES AND FUNCTIONS FOR DISJUNCTIVE KRIGING PROGRAMS * SOQ05
C**"*** *.«. sob06
C"«" **».**« S0007
^*********************************************** ****»****»******»**** sooos
C S0009
C S0010
-------
Disjunctive kriging program Tor two dimensions
301
A******************************************************************
********
SUBROUTINE SORT: SORTS ZD{I) DATA INTO INCREASING ORDER
......
...... ........
(USE WITH PROGRAM ONE) »
*******..*.**..*.*.*.***.*.*************************»*«****
SUBROUTINE SORT(NA,NOB,INX,XO,YD,ZO,YZ,PR)
DIMENSION 1NX(NA),ID(NA),YD(NA),2D(NA),YZ(NA).PR(NA)
00 10 I'l.NOB
X.«MN » 99999999.
DO 15 J«I. MOP
IF(XMIN.UT.ZO(J)) GOTO 15
XMIN Zn(J)
IJ « J
15 CONTINUE
TX»XD(U
TY»YD(I)
Zn(I)-ZO(IJ)
YOd)-YO(IJ)
INX(I).INXdJ)
ZD(IJ)»TZ
XO(U)"TX
YOdJj-TY
INX(1J)«LOC
10 CONTINUE
LCNT « 0
00 20 I'l.NOP
IF(ZO(I).NE.Zn(I»l).AHD.LCNT.EO.O) GOTO 35
lF(ZDd).HE.ZD(I*l).ANO.lCNT.NE.O) GOTO 25
----- duplicate ZO value -----
LCNT » LCKT » 1
GOTO ?0
----- no more duplicates add and calculate results
25 T? , (FlOATd)-.5)/FlOAT(NOB)
LCNT « LCNT » 1
T"P » PRBI(l.O-TZ)
DO 30 KJ«1.LCNT
PR(I-KJ»1) TZ
YZ(I-KJ*1) » THP
30 CONTINUE
LCNT 0
GOTO 20
..... no duplicates (LCNT 0) .....
35 CONTINUE
PH(I) (FLOAT(t)-.5)/FLOAT(NOB)
YZ(I) PRBKl.O-PR(I))
20 CONTINUE
RETURN
EMO
SUBROUTINE CK1:
..*.*........*...*....**.**..*.*.....
CALCULATED THE COEFFICIENTS «" BY
NUMERICAL HERMITE INTEGRATION
U's W(i)*EXP(YU(i)«YW(i)/2.)
***
********
*****
*******
SUBROUTINE CIU(NE,NTRH.NHK,CX,NPOIY,CS)
DIMENSION W(15).YW(15).CK(NHK),CS(NE)
NWC « Maximum order for the Herwite poly **
********
(USE WITH PROGRAM ONE) ****
I**************************************************
soon
S0012
soon
S0014
S0015
S0016
S0017
S0013
S0019
S0020
S00?l
S0022
S0023
50024
S0025
S0026
S0027
S0028
S00?9
50030
S0031
5003?
S0033
S0034
S0035
S0036
$0037
S003*
S0039
S0040
S0041
S0042
S0043
S0041
S0045
S0046
S0047
S0048
S0049
S0050
SOOil
sons?
SOOS3
S0054
S0055
S0056
S0057
S0058
S0059
S0060
S0061
S0062
S0063
S0064
S0065
S0066
S006/
S006G
CKOOl
CK002
DC 003
CK004
CK005
CK006
CK007
CK008
CK009
CK010
won
CK01Z
CK013
CK014
0(015
-------
302 S. R. YATTS. A. W. WAIIIUCK. and D. E. MYHU
C
C
C
C
c
..... VALUFS FOR INTEGRATION: .....
..... YW(I) 8, U(l); 1-1.2 ..... 10 are for 10 term .....
..... YU(I) & V(l); 1-11, 12, ...15 are for 5 term .....
DATA NJ/15/,YW(1)/0.245340713E»00/
fc, YV(2)/0.73747372<»E*00/,U(2)/0.376261601E»00/
&. YW(3)/0.123407622E*01/,U(3)/0.233452289r»00/
*. YV(4)/0.1738S3771£»01/,V(4)/0. 1124517761 *00/
ft, YW(5)/0.22S497400E+01/,U(5)/0.412310259E-01/
*. YW(6)/0.278880606E+01/.U(6)/0.111539518E-01/
&, YU{7)/0.33478S457E*01/.V(7)/0.211861101E-02/
&. YU(8)/0.394476404E*01/.U(8)/0.259968734E-03/
J. YW(9)/0.460368245E+01/,U(9)/0.176028288E-04/
&. YU(10)/0.53a748089E»01/,U(10)/0.447583714E-06/
t. YW(11)/0.34290133/ .«(!!)/. 647852317/,
HYW(12)/1.0366108/.W(12)/.410960S74/,YV(13)/1.75668365/
t.W(13)/.i58479939/,YH(14)/2.53273167/.W(14)/.03320S690/.
&YW(15)/ 3.43615912/,W(15)/.00279908848/
DATA S02PI/2.50662829/
IF(NTRM.EQ.IO) NJ«10
IF(NTRM.EQ.IO) NTRM-0
IF(NTRH.EQ.S) NTRM-10
FACT 1.
00 5 IM.NHK
K (I-U
CK(I) « 0.
IFJK.GE.2) FACT » FACT«FLOAT(K) .
DO 10 J-1»NTRM.NJ
YH « -YW(J)
YP « YW(J)
PHIH « PHI(NE.YH,NPOLY.CS)
PHIP « PHI(NE.YP.NPOa.CS)
CMM ' «(I) * U(J)«PHIfWHK(K,YC)
CK(1) CK(I) + W(J)'PHIP*HK(K.YP)
10 CONTINUE
CK(l) « CK(1)/(FACT«SQ2PI)
£NC
A********
*t
c
C"" **
c
C" SUBROUTINE LEAST: CALCULATES THE LE*ST SP'JARES REGRESSION
<; LINE OF THE FORM:
C
C Y Cl * C2M » C3*X'X »
c
C««"«» (USE WITH PROGRAM ONE)
£***»************************4**********************************ft***«*
C
SUBROOTINF LEAST(Nn8.NA,NC.NE.ZO.PR.AA.NPOLY.CS)
DOUBLE PRECISION AA.OBLE
01KENSION ZO(NA),PR(HA),CS(NE).AA(NC.NC)
(O.IT
***
********
* * *****
NOOU «
NCOL NRQU » 1
DO 1 L'l.NROW
CKTJ16
CX017
C<018
CK019
CK020
CK021
CK02?
CK023
CK024
CK025
CK026
CK027
CK028
CK029
CK030
CK031
CK03?
CK033
CK034
CK035
CK036
CK037
CK033
CK039
CK040
CK041
CK042
CK043
CK044
CK045
CK046
CK047
CIC048
CK049
CK050
CK051
CK052
CKOS3
CK054
CK055
CK056
CK057
LE001
LE002
LE003
LE004
LE005
LE006
LE007
LE008
LE009
LE010
LE011
LE012
LE013
LE014
LEOI5
LE016
LE017
LE019
LE020
LEO?1
oo 5 ««L.NOOW
i^(L.v) '0.000
oo 10 (I.XOB
10
5
00 IS I'l.tOB
15 iA(L.NCOL) A
1 CONTINUE
n8LE(PR(I)"FLOAT(L-l)*PR(I)**FLnAT(H-l))
OBLE(PR(I)"FLOAT(L-1)*ZO(I))
LEO?3
LEO?4
LE025
LEO?6
LCO?7
LE02S
LEO?9
LE030
LE031
LE032
-------
Disjunctive kriging program for two dimensions
303
c '
c
C--..- ^nl v*» thp nxtriv . .*..
» ^ v i v\r winr wo L i i A » »
CALL MAT1NV(NROU, NCOL.NC.AA.CS)
C
RETURN
END
C
£«>.«>>.»>.«...»»».><«>««.<«<>..««>...»«»».>»<'<..>.>>
£******
c. ...... SUBROUTINE HYZINV: CALCULATES THE NORMALIZED YZ FROM THE Z
c******* DATA BY INVERTING THE EQUATION
£«*****
C....... 2 . c(0) , Cd)HX(l.Y) + C(2)HK(2.Y) * C(N)HK(N.Y)
p....***
£....... (ujf uffn PROGRAM ONE)
C
SUBROUTINE HYZINV(NOB.INX,ZO,YZ.NHK.CK)
DIMENSION CK(NHK) .INX(NCB) .ZO(NOB) .YZ(NOB)
CHARACTER'S NAME
DATA TOL/l.OE-6/. NIT/20/
COMrKJN 10. IT
C
NAME ' DATA '
5 DO 10 I'l.NOB
YM1 . YZ(I)-.005*YZ(I) - .0005
Y « YZdKOOS'YZd) » .0005
N 0
C
C ______ Hpnf n c.prc(. o(i i oi/t''U*i|nfliiU'i iw
C..... Npwt nn *c Intpr^tlwP mothrwl -. ...
» r»tfw i un » intcralivc i"t tnuu
15 FM1 » 0.0
F « 0.0
PO ?0 J'l.NHK
ic (j-n
FMl FMl + CK(J)*HK(K.YM1)
20 F » F + «{J)*HK(IC.Y)
C
FMl « FMl-Zn(l)
F « F -ZD{I)
SLOPE - (F-FM1)/(Y-YM1)
IF(ABS(SLOPE).LT.1.0E-10) GOTO 25
IF(ABS(SLOPE).GT.1.0E*10) GOTO 25
YP1 ' Y - F/SLOPE
YM1 Y
Y ' YP1
X » N«l
IF{ABS(Y).GT.10.0) GOTO 25
IF(ASS(Y-YM1).LT.TOL) GOTO 4$
IF(N.GT.NIT) GOTO 25
GOTO 15
C
C» - ^oliition riitf(*rnino - - trw Jtltornjito HIA t hnH .
^ >UIULIV uivtrruiiiu »'jr QI^W'MQVC nnf L**vQ »...
25 IFLAf.-O
N<0
Y YZ(I)
YL Y - .1*A8S(Y) . 1.
YP Y * .l'ABS(Y) * 1.
30 Zl'0.0
00 35 J'l.NHK
K'J-1
35 Zl Zl + CK(J)«HK(K.Y)
IF(ABS(Z1-ZO(I)).LT.. 00001) GOTO 45
IF(N.GT.SO) GOTO 40
C
IF(Zl.GT.ZOd)) YR.Y
IF(Z1.LT.ZD(I)) YL'Y
N'N«1
Y-(YR»YL)/?.
GOTO 30
LE033
LE034
LE035
LE036
LE037
LE038
LE039
HYOOl
......... HYno2
........ HYQ03
........ HY004
........ HY005
« HY006
........ HY007
........ HYOOS
AAAAA**A uw n A A
**" HT009
HYOll
HY012
HY013
HY014
HYQ15
HY016
HY017
HY018
HY019
HY020
HY021
HY022
HY023
Hvn?d
n t \j£ H
UVflOC
~ " Uc 3
HY026
HYQ27
HY028
HY029
"Y030
HYQ31
HY032
HY033
MY034
HY035
HY036
HY037
HY038
HY039
HY040
HY041
HY042
HY043
HY044
HY045
HY046
HY047
HYQ4R
HY049
HY050
HY051
HY052
HY053
HY054
HY055
HY056
HY057
HY058
HY059
HY060
HYQ61
HY062
HY063
HY064
-------
304 S. R. YATTS. A. W. WAJUUCK. and D. E. Mvws
c
c
c
c
c
c
c
c
c
c
c
c
c
c
C1
C'
C '
c
40 IF(IFLAG.E0.21 GOTO 50
YL -3.
YR 3.
Y 0.
N « 0
IFLAG IFLAG » 1
IF(IFLAG.EO.l) GOTO 30
YL -4.
YR 4.
Y .51434556
COTO 30
45 IF(ARS(Y).GT.10.0) POTO 25
YZ(U Y
GOTO 10
____ ^nliit inn rlld nnt rnnvprop KB~»
50 II INX(I)
IF(NAME.EO.' ZCUT ') II'I
WUTE(IT.POO) tl.KA»*,YZ(n.Y.YL.YR,ZO(n.n
10 CONTINUE
RETURN
ENTRY HYZINI(NOI.INX.ZO.YZ.NHK.CIC)
HAKE ' ZCUT '
GOTO 5
800 rORWT(' SOLUTION DID NOT CONVERGE FOR ',A6.' tLOC '.I5.5X
*/.' YZ '.FU.S.SX,' Y .FIZ.S./.' YL '.F12.5.3X.
»' YU '.F12.5./.' 7.0 '.F12.5.3X.' Z '.F12.5./7/)
END
** ftlUfTtrtU Out . OCTllOMt TUP 1t\ uAl Itf fQf\U f tt It W rtDHCO
!* rUHCI ION "HI ; Rt lURKi Int tU V*»LUt * KOr r-n Kin 0**UtR
POLYNOMIAL FIT TO THE NORMALIZED YZ VALUES
» lltfC LJITU BO/V*fiAU f\ttC \
i...... (USE W1IH POOGRAH ONt )
REAL FUNCTION PHI(NE.Y.NPOLY.CS)
01HEMSIOK CS(NE)
P PRB(Y)
SUM . CS(1)
DO 1 I-2.NPOL.'
1 SU- SUM . CS(;)'P"FLOAT(I-1)
PHI SUH
RETURN
ENO
«» FUNCTION PR8I: CALCULATES THE INVERSE OF THE PROBABILITY
>...... FUNCTION GIVES THE X A550CIaTlD WITH D"pO.
« (USE WITH PROGRAM ONE)
REAL FUNCTION PRBI(X)
T AKlNl(X.l.-X)
T « SORT(-2.'ALOG(T*1.E-20))
PRPI T-((.010328'T». 8028531*1*2. 515517)/
& (((.OC130a-T».lP9269)*T»1.4327P8)*T»l.)
IF(I.LT.O.S) RETURN
PRBI » -PPB|
RETURN
ENO
HT065
HY066
HY067
HY068
HY069
HY070
HY071
HY072
HY073
HV074
HY075
HYO'6
HY077
HY078
HYP79
HYOflO
HY081
MYOP3
MYO«4
MYOP-5
HV086
HYPP7
HY033
HY099
MY 090
MY091
HY093
HYQ9«
HY095
MY096
HY097
HY098
HY099
MY 100
HYlOl
PH001
PH009
PH010
PHOII
PH012
PH013
PH014
PH015
PH015
PH017
PM018
PH019
PHP20
PtOOl
........ p|OQ4
........ 9j 006
........ 9[Q07
PI 009
P1010
»IP13
0(014
» 10 15
P[ni6
PI017
P1013
-------
Disjunctive kriging program for two dimensions
305
*********
C*....** «**»«*
c....... SUBROUTINE SOLN: CALCULATES THE DISJUNCTIVE KRIGING ...
C«*"«" SOLUTION: EST, VARK AND PROS ........
r******* ********
C*****" (USE WITH PROGRAK TWO) ........
£*************************************************************************
C
SUBROUTINE SOLN(NB.NC.NO.NE,K.N1.EST.VAPK,ILOC.GOX,CK,PROB.YCUT
» .NZCUT.BB)
REAL GOX{N£).«(NO).PR08(NR).BB(NC),YCUT(NB)
INTEGER ILOC(NF)
COHHON /H/ HHA
C
J K+l
VAR1 0.
HKV ' 0.
IF(K.LE.l) FACT 1.
IF(K.GE.Z) FACT FACT'FLOAT(K)
IF(K.GT.O) r-OTO 10
for K«0; 8B(I) 1/Nl and GOX(I) « 1.0
00 5 I'l.NI
VAR1 . VAB1 + 1.0/FLOAT(N1)
CALL H1(J.ILOC(I))
5 WV » HKV + HHA/FLOAT(H1)
GOTO 30
for K>0; BB(I) from matrix solution
10 00 15 I-l.Nl
CALL Hl(J.ILOCU))
HKV ' HKV + BB(I)«HHA
15 VAR1- VAR1* B8(I)*GOX(I)
IF(NZCUT.EO.O) GOTO 25
00 20 M«1,NZCUT
20 PROB(f) > PROB(M) * HK(K-1,YCUT(M))«HKV/FACT
25 VARK ' VARK + CK(J)*CK(J)«FACT«VAR1
30 EST « EST * CK(J)-HKV
RETURN
EHO
C******************************************************************************
£* ********
c.****** SUBROUTINE MATRX: CREATES THE DISJUNCTIVE KRIGING MATRIX ****«
C******* EQUATION. ********
r******* ********
C**«**«* (USE WITH PROGRAM TWO) ********
£************************************»********************************
C
SUBROUTINE MATRX(NA,NC,NE,NZ,Nl,N08,XO,YO,XD,YO,ILOC,RLOC,RfAX
« .r-OX.NCOL.NROW.fOOE.C.IBLK.AA.IFLG)
DOUBLE PRECISION AA(NC.NC).OBLE
REAL XD(NA).YO(NA).RLOC(NE),c(3),GOX(NE)
INTEGfR ILOC(NE)
locate NZ nearest neighbors
DO 1 J-l.NZ
ILOC(J) 0
1 RLOC(J) 9.9E5
00 10 I'l.NOB
OX XO(I)-XO
DY YD(I)-YO
R « SORT{nX»DX+OY«DY)
IF(R.GT.RMAX) GOTO 10
SN001
SN002
SN003
SN004
SN005
SNOO*1
SN007
SN008
SN009
SN010
SN011
SN01?
SN013
SN014
SN015
SN016
SN017
SN018
SN019
SH020
SN021
SN022
SN023
SN024
SN025
SW26
SN027
SN028
SM029
SN030
SN031
SN03?
SN033
SN034
SK035
SN036
SM037
SN03S
SN039
SN040
SN041
SH042
SN043
SN044
SN045
SN046
SH047
MX 001
MX002
MX003
VX004
HX005
MX006
MX007
MX008
MX 009
MX010
MX011
MX012
MX013
HX014
MX015
HX016
MX017
HX018
MX019
HX020
MX021
MX022
MX023
MX024
MX025
-------
306 S. R. YATO, A. W. WAWUOC. ind D. E. Mvws
C HX026
C ..... check If 8 goes In list ..... K*0?7
ZMX « -9.9E5 HX028
DO 5 O'l.NZ HX029
IF(ZMX.LT.RIOC(J)) XI J W30
IF(ZMX.LT.RLOC{J)) ZHX » KLOC(J) MX031
5 CONTINUE MX032
C MX033
IF(R.GT.ZMX) GOTO 10 WI034
RlOC(Kl) R PX035
ILOC(Kl) I MX036
10 CONTINUE MX037
C MX038
Nl«0 MX039
DO 15 I'l.NZ HX040
IFdlOCdJ.LE.O) GOTO 15 PX041
Nl Nl + 1 MX042
15 CONTINUE MX043
IF(Nl.EO.O) IFLG'l MX044
IF(IFLG.EO.l) RETURN HX045
NCOL NU1 HX046
NROU Nl HX047
C HX048
C ----- fill In main diaaonal and rhs of matrix ----- MX049
00 25 I'l.Nl MX050
AA(I,I) l.ODOO KX051
C MX052
C ..... block disjunctive kriglnq ..... MX053
IF(IBLK.EQ.O) GOTO 20 HX054
1F(IBLK.NF.O)AA(I.N1+1) MX055
08LE{8UPTS(XO.rO.XD(KOCn)).YO{ROCU)),C.MOOE)) MI(056
GOTO 25 MX057
20 R « RLOC(I) HX059
AA(I.N1*1) » DBLE( ViR.'OIS.-.HODE) ) "K060
25 GOX(I) AA(T,NU1) MX061
C ^X062
C ..... fill in off-diagonals (calculate upoer-half ..... «063
C ..... and assign to lower) ..... MX064
DO 35 I*1,N1-1 MX065
00 30 J-I+1.N1 MX066
DX « XO(ILOC(JJ) - XOHLOCd)) MX067
Or YOULOC(J)) - YO(1LOC(D) HX068
R ' SOST(DX*OX+OY*OY) MX069
C ..... upper half ..... M070
AA(I.J) OBLE( VARIO(R.C.KODE) ) HX071
C ..... lower half ..... " MX072
AA(J.I) AA(I.J) HX073
30 CONTINUE HX074
35 CONTINUE HX075
C MX076
RETURN MX077
END ^X078
C AV001
r*********a********«*******«****»********************ft********»**************** ^ y QQ 2
£****** ******** AV003
C*.....* SUBROUTINE AVECO«: CALCULATES THE INNEP BLOCK COVARJANCE » AV004
C******* ******** AV005
C******* (USE WITH PROGRAH TWO) *** AV006
Q*****************************»*******************************************ft**** AV007
C AV008
SUBROUTINE AVECOR(ND,BUCOV,C,MOOE,NHK,CK) 4V009
DIMENSION QC(HO),C(3) avOlO
INTFGER'4 1SF.ED AV011
COMMOfl /BLK/ NX.NY.WIOX.WIDY.ISEEO AVOU
C AW 3
SUM=0.0 AV014
OX » WIOX/NX AVOI5
OY » WIDY/NY A\016
XO « (DX-WIDX)/2. AV017
C AV018
DO 10 I'l.NX AV019
XO « XO » OX AV020
-------
Disjunctive kriging program for two dimensions
307
YO (OY-WIOY)/2.
00 10 J'l.NY
YO « YO + OY
XI « (DX-WIDXJ/2.
00 10 K'l.NX
XI « XI * DX
Yl « (DY-WIDY)/2.
DO 10 L-l.NY
Yl Yl * OY
move a random deviation from center of sub-block
random number assumed to be [0,1]
Xll XI « OX*(.5-RAN(ISEEO))
Yll » Yl * DY*(.5-RAN(1SEEO))
R SQRT({XO-X11)*(XO-X11) + (YO-Y11)*(YO-Y11))
SUM SUM + VARIO(R.C.MOOE)
CNT CNT + 1.0
10 CONTINUE
BLKCOV « SUH/CNT
FACT - 1.0
SUM « 0.0
DO 20 I-l.NHK
IF(I.GE.2) FACT - FLOAT(I)«FACT
SUM SUM * FACT«CK(I+l)«CK(l+l)*BLKCOV**FLOAT(I)
20 CONTINUE
BLKCOV « SUM
RETURN
END
ft*********************************
A**************************
£***** ********
c....... SUBROUTINE 8UPTS: CALCULATES THE BLOCK VARIANCE BETWUEEN **«
£»«** SAMPLE POINTS AND BLOCK *..*..**
£»*«« ********
c.****** (USE WITH PROGRAM TWO) .»*.**..
£*«*******************«****»**«**********
C
REAL FUNCTION 6LKPTS(X3,YB.X,Y,C.HOOE)
DIMENSION C(3)
INTEGER** ISEEO
COKMON /BLK/ NX.NY.UIDX.UIDY.ISEED
C
SUM-0.0
CUT»0.0
OX WIOX/NX
DY « WIOY/NY
C
xo XB » (ox-wiox)/?.
no 5 I-I.NX
YO YB » (OY-UIDY)/?.
DO 10 J'l.NY
C
C
C
C
10
5
XOO
YOO
R
SUM
CNT
YO
XO
8LKPTS
RETURN
END
- move a random deviation from center of sub-block
random number assumed to be [0.1]
XO + OX«(.5-RAN(ISEEO))
YO + DY«(.5-RANMSEEO))
SQRT((XOO-X)«(XOO-X) * (YOO-Y)-(YOO-Y))
SUM + VARIO(R.C.WOE)
CNT » 1.0
YO + OY
XO + OX
SUM/CNT
PA*****************************************************************************
£** ********
c.* SUBROUTINE ATOK: RAISES THE "AA" MATRIX TO THE K/(K-1) *...*
c******* POWER so THAT THE -AA- MATRIX NEED NOT
C*«**«»» BE RECOMPUTED FOR EACH K .......
AV021
AV022
AV023
AV024
AV025
AV026
AV027
AV028
AV029
AV030
AV031
AV032
AV033
AV034'
AV035
AV036
AV037
AV038
AV039
AV040
AV041
AV042
AV043
AV044
AV045
AV046
AV047
AV048
AV049
AV050
BL001
BL002
8L003
BL004
BL005
BL006
BL007
BL008
BL009
8L010
BL011
BL012
BL013
6L014
BL015
BL016
BL017
8L018
BL019
BL020
BL021
BL022
BL023
BL024
BL025
BL026
BL027
8L028
BL029
BL030
BL031
8L032
BL033
BL034
BL035
BL036
BL037
AT001
AT002
AT003
AT004
AT005
AT006
-------
308
S. R. YATCS, A. W. WAMUCX. and D. E. MYWS
£
c
V
20
10
C
C
p
r
1
f
10
r
^
?0
21
22
C
C
30
31
C
C
1000
800
C
c....-
C****'
c
VUit Ml in rKlHjKAr iwuj ""
SUBROUTINE ATOMNC.NE.K.NROW.NCOL.GOX.AA)
RE*L GOX(NE)
DOUBLE PRECISION AA(NC.NC) .KK.OBLE
KK « OBLE(FLOAT(K)/FLOAT(IC-1))
DO 10 I'l.NROW
00 20 J-l.NCOL
IF(AA(I,J).EO. 1.00*00) GOTO 20
IF(AA(I,J).LT.1.0D-10) GOTO 20
AA(I.J) « AA(I.J)**KK
CONTINUE
GOX(I) - SNGL(AA(I.NCOD)
CONTINUE
RETVRN
END
**
'* FUNCTION VARIO: RETURNS SILL - VARIOGRM TO CALLING PROGRAM
'* (USE WITH PROGRAM TWO) »*
SEAL FUNCTION VARIO(H .C.WJDE)
DIMENSION C(3)
IF(R.GT.O.O) GOTO 1
VA^IO « cm + cm
3 E TURN
GOTO(10.20,30), MOOE+2
viRIO » C(i) » C(2)*(1.0 - EXP(-R/C(3)))
GOTO 1000
:rf.n).LE.n.o) GOTO ?i
IF(R.GE.C(3)) GOTO 22
VARIO « C(l) * C(2)*R/C(3)
GOTO 1000
VARIO C(l) + C(2)
GOTO 1000
IF(R.GE.C(3)) GOTO 31
TMP « R/C(3)
VARIO - C(l) * C(2)*(1.5*TMP..5»(TMP*TMP*THP))
GOTO 1000
VARIO - C{1) * C(2)
VARIO c(i) * c(2) - VARIO'
FORMA T(A50)
RETURN
END
>»
** SUBROUTINE HCALC: CALCULATES THE HERMITE POLYNOMIAL FOR ***«
» ORDERS 0 TO NHK FOR EACH NORMALIZED ****
» DATA POINT. THIS SPEEDS UP COMPUTATIONS. ****
(USE WITH PROGRAM TWO) ***'
SUBROUTINE HCALC(NA.NOB.NHK.YZ)
DIMENSION YZ(NA)
COMMON /H/ HHA
AT007
AT010
ATOM
AT012
AT013
ATOM
AT015
AT016
AT017
AT018
AT019
AT020
AT021
AT022
AT023
AT024
AT025
VA001
VAQ08
VA009
VAQ10
VA011
VA012
VA014
VA015
VA01K
VA017
V»013
VA019
VAO?0
VAQ?1
VA02?
VAO?3
VA024
VA025
VA026
VA027
VA028
VA029
VA030
VA031
VA032
VA033
VA034
VA03S
VA036
VA037
VAQ38
VA039
VA040
VA041
HC001
HC010
HC011
HC012
HC013
-------
Disjunctive kriging program for (wo dimensions 309
C ' HCOH
DO 1 J'l.NHK HC015
K « (J-l) HC016
DO 5 I»1,NOB HC017
IF(K.EQ.O) ROTO 10 HC018
IF(K.EQ.l) GOTO 15 HC019
C HC020
C FOR K > 1 CALCULATE THE HERHITE POLYNOMIAL HCO?1
HHA « HX(K,YZ(I)) HC022
GOTO 5 HC023
C HC024
C FOR K'O HERHITE POLYNOMIAL « 1.0 HC025
10 HHA » 1.0 ' HC026
GOTO 5 HC027
C HC02S
f FOR KM HER^ITE POLYNOMIAL » v; HC029
15 HHA « YZ(I) HC030
5 CALL H(J.I) K031
1 CONTINUE HC03?
C «C033
RETURN MC034
FNd HC035
C HP001
£.****.********«*******.**********************.*******************.******** HP002
C*.....* ........ HP003
C**...*» SUBROUTINE H: GIVES THE VALUE OF THE HERMITE POLYNOMIAL «* HP004
C**»«**» FO1? ORDER K AND LOCATION (OR DATA) I. *« HP005
(;......« >...«.. HP006
C*"**«* (USE WITH PROGRAM TWO) * HP007
Q***************************»*******ft******************************«**********4 HP008
C HP009
SUBROUTINE H(K,I) HP010
COMMON /H/ HHA HP011
VIRTUAL HH(?0.220) HP012
C HP013
C store hermite polynomial in table HPOU
HH(K.I) HHA HP015
RETURN HP016
C HP017
C look up hermite polynomial in table HP018
ENTRY Hl(K.I) HP019
20 HHA HH(K,I) HP020
RETURN HP021
END HP022
C W001
£**************************************************************«***** MVQQ2
£ »« MVQ03
c....... SUBROUTINE MATINV -- SOLVES THE MATRIX EQUATION HV004
Q******> ******** PVQQ5
C"**"* (USE UITH BOTH PROGRAMS) ........ HV006
(^.......*.......»......********.*********»*#***»********».***............» HV007
C MV008
MV009
SUBROUTINE MATINV(NROW.NCOL,NC,AA,BB) MV010
DOUBLE PRECISION AA(NC.NC),AURK(21,21).DBLE.SUM MV011
REAL BB(NROW) MV012
COMHON 10,IT MV013
C MV014
C MV015
C put a a Into working matrix MV016
DO 100 I-l.NROU MV017
00 101 J'l.NCOL MV018
101 AWRK(I.J) « AA(l.J) MV019
100 CONTINUE MV020
C MV021
C MV022
C partial pivoting Is not necessary MV023
C since maximum value 1s in main MV024
C diagonal of krlgtng matrix MV025
00 1 I'2,NCOL ' MV026
IF(AWRK(1.1).NE.0.0000) GOTO 1 MV027
VH[TE(IT,800) MV023
VRITEU0.800) MV029
STOP *V030
-------
310 S. R. YATO. A. W. WARRICK. and D. E. Mras
1 AWRK(l.I) « AWRK(l.I)/AVRK(l,l) MV031
C MV032
00 5 L-2.NROV MV033
DO 10 1'L.NROW HV034
SUM « 0.0000 MV035
C HV036
00 15 K-l.L-1 MY037
15 SUM SUM * AWRK(I.K)*AWRK(K,L) W038
10 AWRK(I.L) « AVRK(I,L)-SUM MV039
C MV040
00 20 J-L+l.NCOL MV041
SUM*O.ODOO MV042
C MV043
00 25 K"1,L-1 MV044
25 SUM » SUM + AWRK(L,K)*AWRK(K,J) ^045
IF(AWRK(L.L).NE.O.ODOO) GOTO 20 MV046
WRITE(IT,800) MV047
WRITE(I0.800) MV048
STOP KV049
20 AWRK(L.J) « (AWRK(L,J)-SUH)/AWRK(L,L) HV050
5 CONTINUE MV051
C MV052
C solution MV053
88(NRPU) « SNr,L(AWRK(NROW.NCOL)} MV054
00 30 H=1,NROU-1 MV055
I - NROW-M KV056
SUM - 0.0000 MV057
C MV058
00 35 JM + l.NPOW MV059
35 SUM « SUM + AWRK{I,J)*OBLE{BB(J)) MV060
30 BB(I) ' SNGL(AWRK{I,NCOL) - SUM) MV061
C . MV062
800 FORMATf////' ERROR: Zero value in main diagonal of matrix.',/, MV063
*SX.'Execution will cease.'///) MV064
RETURN MV065
END MV066
C P8001
eft***************************************************************************** P800?
^******* ******** PB003
c******* FUNCTION PRB: CALCULATES THE PROBABILITY INTEGRAL ******** PB004
£******* ******** pBOOS
C******* (USE WITH BOTH PROGRAMS) ****** PB006
C****************************************************************************** P3007
C PBOOS
REAL FUNCTION PRB(X) PB009
Z APS(X) PB010
IF(Z.LT.5.) GO TO 1 PB011
PRB .5*(1.»Z/X) PB012
RETURN PB01J
1 PRB » .5*(1.+(((.019527*Z+.000344)*Z+.115194)«Z+.196854)*Z) PB014
*"(-4) PB015
IF{X.LT.O.) RETURN PB016
PRB » l.-PRB P8017
RETURN PB018
END PB019
C HK001
Q**A*************************************************************************** HK002
£< tttttttt HK003
c....... FUNCTION HK: CALCULATES THE HERMITE POLYNOMIAL OF ORDER ***** HK004
C«*«"" K USING A RECURSIVE RELATIONSHIP ** HK005
(^....... ..«** HK006
C*"**" (USE WITH BOTH PROGRAMS) ** HK007
Q***A***************************************************************tt*«******** HK008
C HK009
REAL FUNCTION HK(K,Y) HK010
C HK011
HK » 1. HK012
IF(K.EO.O) GOTO 5 HK013
HKM1 ' HK HK014
HK « Y . HK015
IF(K.EO.l) GOTO 5 HK016
C HK017
C HK013
-------
Disjunctive kriging program for two dimensions 311
C arbitrary K HKO.">
00 1 l-l.K-1 HK020
HKP1 « Y«HK - FLOAT(I)*HKrn HK021
HKM1 HK HK022
HK » HKP1 HK023
1 CONTINUE HK024
5 RETURN HK02S
END HK026
C HK027
APPENDIX 2A
Example input data for DISCALC.FTN. Data file H-OJ created using "Tne Data File Input Instructions" from
program documentation. For brevity, only the first and last ten entries of the I. X. Y. and Z data are reported
OA001
BARE SOIL TEMPERATURE DATA -- EXAMPLE OA002
OA003
57755 23.000 100.000 OA004
000 0.000 0.000 0.000 OA005
3 3 2.000 0.000 16.000 100.000 DA006
62.50 64.00 65.0000 66.00 67.00 OAOO7
1 0.700 2.400 23.000 OA008
1234 OA009
TEMP OA010
(I5.2F10.3.F10.5) OA011
108 17.000 200.000 60.95000 OA012
94 33.000 243.000 61.07000 OA013
75 27.000 198.000 61.19000 OA014
95 33.000 181.000 61.35000 OA015
100 34.000 188.000 61.44000 OA016
87 31.000 37.000 61.45000 OA017
99 33.000 147.000 61.65000 OA01S
104 17.000 93.000 61.71000 DA019
60 23.000 176.000 61.78000 OA020
96 33.000 94.000 61.81000 DA021
OA022
DA023
OA024
OA025
5 2.000 198.000 66.47000 DA026
20 7.000 149.000 66.74000 DA027
34 13.000 183.000 67.05000 (U02S
7 3.000 206.000 67.33000 0*02"
15 4.000 187.000 67.43^00 OA030
27 11.000 172.000 67.43000 OA031
12 4.000 55.000 67.48000 OA032
17 5.000 149.000 67.61000 PA033
3 1.000 190.000 67.83000 OA034
18 5.000 252.000 68.25000 OA035
APPENDIX 2B
Example of an output file from DISCALC.FTN. This file is used, without modification, as an input file for
D1SJUNC.FTN. For brevity, only the first and last ten entries of the I. X. Y. Z. and YZ data are reported
OB001
BARE SOIL TEMPERATURE DATA -- EXAMPLE OB002
08003
5 7 5 23.000 100.000 DB004
000 0.000 0.000 0.000 OB005
3 3 2.000 0.000 16.000 100.000 08006
62.500000 64.000000 65.000000 66.000000 67.000000 OB007
-0.682487 0.178525 0.674979 1.127083 1.545975 DB008
6.3890934E+01 1.7091980E+00 1.8315849E-01 -6.8178676E-02 -4.4367205E-C2 OB009
-3.6764962E-03 -5.8487770E-03 OBOiO
1 0.700 2.400 23.000 OB011
TEMP 08012
108 17.000 200.000 60.95000 -2.89787889 08013
94 33.000 243.000 61.07000 -2.81422019 OB014
-------
312
S. R. YATES, A W. WAJUICK. and D. E. MYEJU
75 27.000 198.000 61.19000 -2.69484258
95 33.000 181.000 61.35000 -2.21075654
100 34.000 18P.OOO 61.44000 -1.76944888
87 31.000 37.000 61.45000 -1.74350667
99 33.000 147.000 61.65000 -1.40478790
104 17.000 93.000 61.71000 -1.33265126
60 23.000 176.000 61.78000 -1.25641036
96 33.000 94.000 61.81000 -1.22578073
5
20
34
7
27
15
12
17
3
18
2.000
7.000
13.000
3.000
11.000
4.000
4.000
5.000
1.000
5.000
198.000
149.000
183.000
206.000
172.000
187.000
ss.noo
149.000
190.000
252.000
66.47000
66.74000
67.05000
67.33000
67.43000
67.43000
67.48000
67.61000
67.83000
68.25000
1.32675445
1.43894422
1.56650829
1.6816266*
1.72296810
1.72296810
1.74372768
1.79804289
1.89167035
2. OS 196902
OB015
OB016
OB017
06013
06019
08020
08021
08022
08023
08024
08025
06026
08027
08028
OB029
OR030
D8031
08032
D"033
06034
08035
08036
APPENDIX 2C
Example input file for DISJUNC.FTN. This file was created using the data file given in Appendix 2B (note:
only a part of the complete file is given in Appendix 2B). The data in this file wer* used to obtain the
estimates at the points: (2. 200). (IS. 200) and (34. 200) in Figure 2
BARE SOIL TEMPERATURE DATA -- EXAMPLE
5 7 5 23.000 100.000
000 0.000 0.000 0.000
3 1 2.000 200.000 16.000 0.000
62.500000 64.000000 65.000000 66.000000 67.000000
-0.632487 0.178525 0.674979 1.127083 1.545975
6.3890934E»01 1.7091980E+00 1.8315849E-01 -6.8178676E-02 -4.4367205E-02
-3.6764962E-03 -5.8487770E-n3
1
TEMP
0.700
2.400 23.000
108
75
91
72
54
71
64
56
19
5
7
15
3
17
27
32
26
20
26
24
21
6
2
3
4
1
.000
.000
.000
.000
.000
.000
.000
.000
.000
.000
.000
.000
.000
200
198
201
203
207
201
199
199
188
198
206
187
190
.000
.000
.000
.000
.000
.000
.000
.000
.000
.000
.000
.000
.000
60
61
62
62
63
63
63
64
65
66
67
67
67
.95000
.19000
.43000
.91000
.10000
.60000
.65000
.59000
.54000
.47000
.33000
.43000
.83000
-2
-2
-0
-0
-0
-0
-0
0
0
1
1
1
1
.89787889
.6948425?
.72949600
.42498398
.31316853
.03369232
.00667833
.47700810
.92427117
.32675445
.68162668
.72296810
.89167035
OC001
OC002
OC003
OC004
OC005
PC006
OC007
ncoos
OC009
OC010
ocon
OC012
OC013
DC014
OC015
OC016
OC017
OC013
DC019
OC020
OC021
OC022
OC023
OC024
OC025
APPENDIX 2D
Example of output from DISJUNC.FTN
OISJU«CII»1 1«IGI«G CM » C1IO NMIII
rcw (StiMiiKC y»M*uir omnociir
CSIIMIIM
8"[ SOIL TfWEBATUBC OAU .. llttfil
00001
00002
0000)
0000*
00005
t>000«
00007
00008
00004
00010
00011
0001?
00013
DOOM
-------
Disjunctive kriging program for (wo dimensions
313
wiKVi* *j>-8t» OF HE ABES' MICH
CO«El»tlO>l FUHCIIW MOOfL TtPf
VS3EI
SILL "ItlS MXiSET
SJSGE
OB'.f'vtO 0»TA
13
80«S S
? 3. 0000
100.0000
S?H(>ICAl
0.??S8
0.774?
73.0000
:. %?. i i TEW it oes. no. i t it»w
i;« 17. KO 709.000 60.950 -?. 897879 56 ?1.000 199.000 64.590
". 77.000 198.000 (1.190 -7. (91841 19 (.000 188.000 6S.S40
11 1?.'AO 701.000 (?.4SO -0.7?949( 5 7.000 198.000 66.470
72 }'..vn 703.000 67.910 -0.4?4984 7 3.000 ?06.000 67.330
'.4 K.IVt ?07.000 63:100 -0.3131(9 IS 4.000 187.000 (7.430
1 ?(.t.Vn 199.00R 61.650 -0.006(76 0 0.000 0.000 0.000
Ci'i tf: MtBHIIf POiTliONIAlS:
-fin 63.8941)1 xfM « 63.89093
09,-j-iCF. 3.07(19 ««a|«HCl 3.0A9S4
MEBH1H POL'KWIAi COCFF1CIEK1S:
k Ck
0 6.3990-m-Ol
1 1.70919SE-00
3 -6!9I7863E-0?
OISJUXCUtE (>IGI«G (SIIM1ES
k Ck
4 -4.4]«7?1E-0?
S -3.676496E-03
6 -S.848777E-03
7 O.OOOOOOC-01
P»0&A«ILIIt THAI tMf ESIIKAIE IS GBfAlEU IHAH /CUT
11
0.477001
0.9?4?71
1.3?(7S4
l.6816?7
1.7??968
1.891670
0.000000
I I < ES' *(
5 ?.00 ?On.OO ((.?70 l.S'»
5 13.00 ?00.00 6?.187 1.316
5 W.OO JOO.OO (?.61S 1.754
7. cut
( i.?SO)( (.400)1 6.5fX3){ (.600)1 6.700)
( C*01)( I-OIK f'ODI (01)( [-01)
......§-...---...-.............»
1.000 0.999 0.777 0.500 0.?63
0.335 0.?46 0.081 0.001 O.OCO
0.494 0.1?? 0.0?0 0.000 0.000
000 IS
00016
00017
00018
00019
00070
00071
0007?
00073
00074
0007!
00076
00077
00078
00079
00030
00031
0001?
00033
00034
P003S
00036
00017
00018
00039
00040
00041
0004?
00043
PD044
00045
00046
00047
00049
00050
00051
COM?
OOOS3
00054
0005 S
00056
00057
00058
000(0
000(1
0006?
000(3
000(4
Yllue^ m*en entire dit« ^rt u\e4 *n c«Uul«tion
000(6
ooo6;
000(3
000(9
00070
OPO/1
0007?
------- |