Bug Correction

The following programs were corrected:

Two bugs in Program kb2d.

1. The subroutine cova2 used to calculate the 2D covariance values from a given model, uses the input parameters of GSLIB 1.0, not 2.0. Consequently, the exponential and Gaussian variogram models are defined in terms of range ``parameter'', and not in terms of practical range, as is the case with program kt3d, version 2.0. Now the subroutine cova2 has been modified to follow GSLIB 2.0's convention.

2. When using ordinary kriging (OK), there is a mismatch in the estimation variance between kb2d and kt3d. This mismatch is due to the following bug in program kb2d: the augmented OK matrix contains covariance values C(0) instead of ones (for stability purposes), but the contribution of the Lagrange parameter to the OK variance has not been multiplied correspondingly by C(0).

One bug about the multiple grid search in the simulation programs(sgsim, sisim, sisim_gs, sisim_lm)

Line 838 in program sgsim: sim(index) = sim(index) + imult should be
                sim(index) = sim(index) -imult
Line 655 in programs sisim: sim(index) = sim(index) + imult should be
                sim(index) = sim(index) -imult
Line 654 in program sgsim_gs: sim(index) = sim(index) + imult should be
                  sim(index) = sim(index) -imult
Line 599 in program sgsim_lm: sim(index) = sim(index) + imult should be
                  sim(index) = sim(index) -imult
Two bugs were found in program kt3d
The most important one is the new procedure of simple kriging with local means (koption=2).
The problem is that the local means are not taken into account in the program.
The following lines should be modified:
1) line 305: if(ktype.eq.3.or.ktype.eq.2) then
    instead of
    if(ktype.eq.3) then
2) lines 981-984: if(ktype.eq.0) then
    est = est + real(s(j))*(vra(j)-skmean)
    else if(ktype.eq.2) then
    est = est + real(s(j))*(vra(j)-vea(j))
else
    est = est + real(s(j))*vra(j)
endif
instead of:
if(ktype.eq.0.or.ktype.eq.2) then
    est = est + real(s(j))*(vra(j)-skmean)
else
    est = est + real(s(j))*vra(j)
endif
The second problem is that, in the computation of statistics regarding prediction error (cross validation option), unestimated locations are also taken into account, which greatly overestimates the average prediction errors in the case where some locations are not estimated.
The lines 1019-1023:
if(true.ne.UNEST.and.est.ne.UNEST) err=est-true
write(lout,'(7(f12.3,1x))') xloc,yloc,zloc,true, est,estv,err
xkmae = xkmae + abs(err)
xkmse = xkmse + err*err
should be replaced by
if(true.ne.UNEST.and.est.ne.UNEST) then
    err=est-true
    xkmae = xkmae + abs(err)
endif
write(lout,'(7(f12.3,1x))') xloc,yloc,zloc,true, est,estv,err

Correction of the implementation of the conditioning data search in sisim, sisim_lm and sisim_gs


GSLIB program sisim allows relocating the original data to the closest sim ulation grid nodes in order to speed up the conditioning data search. When two d ata are assigned to the same grid node, the previous sisim version retains the closest one, without distinguishing hard data from soft data.

The new sisim version gives priority to the hard datum whenever a hard and a soft datum are close to the same grid node. When data are not relocated to grid nodes, a super block search strategy is used to search the nearest conditioning data. In the previous version of sisim , sisim_lm and sisim_gs, after setting up the super block strategy, the data values do not match the data locations anymore, and the restriction of the number of conditioning soft data usually fails. Those failures are corrected in the new version of the 3 programs affected.

Description - The third angle in all GSLIB programs operates in the opposite direction than specified in the GSLIB book.

Explanation - The books says (pp27) the angle is measured clockwise when looking toward the origin (from the postive principal direction), but it should be counter-clockwise. This is a documentation error. Although rarely used, the correct specification of the third angle is critical if used. For a detailed description of the rotation angles, see a note by Zhanjun Ying.