Sarma D.D. Geostatistics with Applications in Earth Sciences

Подождите немного. Документ загружается.

156 Geostatistics with Applications in Earth Sciences

283 T(1)=5 .

284

IF(T

(2))285,28

,28

285

IF(

-T

(2)-5

.)28

7,286,286

286 T

(2)=-5

289

287

IF(T

(2) -5 .

)289,288,288

288 T(2)=5 .

289

1,2

R(L)=T(L)

IF(R

(L))802,803

,803

802 R(L)=-R(L )

803 X=R(L)

)

=l.

804

II=2

, 31

III

=II-

804 P

(II

)

=(P(III)

(II)

)

*X *X

P (32 )=X*.3989 4228

805

UM+P

(JJ)*D

805 D=-D

15 PP(L )=.5 +S

*P(32)

IF(

T (

1))20,2

1,2

IF(

(2))22,22,23

22 Q(I )=PP(1 )-PP (2)

23 Q

(I)

=PP(1 )+PP (2)

-1.

IF(T

(2))23

,24 ,24

24 Q

(I)

=PP(2 )-PP (1)

50 E

(I)

(I)*TOT

UM+E(I

)

IF(E

(I ) .EQ. O.O)GO

1007

1008

CHI

HI+(FR

(I)-E

(I)

)* (FR(I ) - E( I )

E(I)

301

RMA

T(3X, ' FREQ. ' ,2X, '

%AGE'

,3X, '

RMAL

PROB

.LI

MITS

' ,

1 6X, 'Q (I ) , , 2X, 'EXP.

FREQ'

1007 WRI TE(2, 998 )FR(I ) ,

PCR

(I) , T(1) , T(2) , Q

(I)

, E

(I)

998

RMAT(F6

.0,

lX,F8

1,F8

.2,

'TO' ,F8

.2,3X,F8

.3,F8

.1 )

107

NTI

NUE

=TO

TAL

WRITE(2, 1111 )

WRITE( 2, 300)

300

RMA

(I/

20X , 'TOT

NO.

POINTS

'F5

ZAR

=ZBAR+

/2 .

=EXP

(ZAR

)-C

FACT

1=1 .+

/2 .+ (

*(TN-1 . ) )

/(8

.0 *(T

N+

1 . ) )+

1 ((V2**3) * ( (

-l.

) **2) ) / (48 . 0* (

+l.

) * (T

N+3

. ) )

Computer Software 157

FACT2 = ( (V2 * *4 ) * ( (TN

-1.

) * *3 ) ) / (3 84 . * (TN+

) * (TN+3 . ) * (TN+5 . ) ) +

(V2**5)*(

(TN

-1

.)**4))/(3840

.*(TN+1

.)*(TN+3

.)*(TN+5.)*(TN+7.))

=EXP(ZBAR)*

(FACT1+FACT2)

-C

=((T2**2)*(EXP(V2/TN)

-1.))

=EXP(ZBAR)

=U3**2/V2**

=U4/V2**2

WRI TE (2 , 8 2 ) ZBAR

FORMAT

(5X , ' A. M. BASED ON FREQUENCY

DISTN

. , (ZBAR)

= ' F 8 . 2 1)

WRI TE ( 2 , 8 3 ) SD2

FORMAT

(5X , ' STD. DEVIATION

FROM

ZBAR

. , = ' F 8 . 2 1)

WRI TE (2 , 92 ) ESUM

FORMAT

(5X,

' S

OF THE EXPECTED FREQUENCIES = ' F 8 . 2 1)

WRI TE ( 2 , 84 ) CHI

FORMAT

(5X , ' OBSERVED CHI -SQUARE VALUE. ,

...

•

•..

• • • • •

F8 . 2 1)

WRI TE (2 , 8 5 ) T1

FORMAT

(5X , ' KRI GE ' S ESTIMATOR •

.••

.•••

.••

. . .

F8 . 2 1)

WRI TE (2 , 8 9 ) T2

FORMAT

(5X , ' S I EHEL ' S

ESTIMATOR..

•

• •

•..

•

•...

F8.

21)

WRI TE (2 , 8 6 ) GM

FORMAT

(5X , ' GEOMETRI C

MEAN

F8 . 2 1)

WRITE(2,88)VE

FORMAT

(5X , ' VARI ANCE OF T

-ESTIMAT

••

. . . • • . . • • • • • = ' F 8 . 2 1)

WRI TE (2 , 90 ) CS

FORMAT

(5X , ' COEFF I CI ENT OF SKEWNESS

••

. . . • • . . • • • • • = ' F 8 . 2 1)

WRI TE (2 , 91 ) CK

FORMAT

(5X , ' COEFF I CI ENT OF KURTOSIS = ' F8 . 2 1)

GO TO 5

END

************************************************************

C SUBROUTINE

HISTO

. FOR FOR HISTOGRAM

************************************************************

SUBROUTINE

HISTO

(N1 , FR)

CHARACTER *1 ICHAR1

,ICHAR2

DIMENSION

B(104,100)

,FR(104)

DATA

ICHAR1,

ICHAR2 / ' H' " , /

DO 4 I =1

,N1

IFR

=FR(I)

DO 3 J

=1,1

IF(J

-IFR)1,1,2

1 B (I , J )= I CHAR1

GO TO 3

2 B (I , J )= I CHAR2

3 CONTINUE

4 CONTINUE

WRI TE (2 , 5 )

FORMAT

(7X,

' 0 '

,9X,

'10

,8X,

' 2 0 '

,8X,

30'

, 8X,

'40

,8X,

'50

8X,'60

',8X,'7

',8X

'80')

158 Geostatistics with Applications in Earth Sciences

WRI TE(2, 6)

FORMAT(lX

,80(lH-))

9 I =1,N1

WRI TE(2, 8)

(B(I,J)

=1,80)

WRI TE

(2,7)

FR(I) , (B

(I,

=l,

80)

7 FORMAT(lX,F6.1,2X,100A1)

WRI TE(2, 8)

(B(I,J)

=1,80)

FORMAT(9X,100A1)

CONTINUE

RETURN

END

PROGRAM:

LN.FOR

************************************************************

FITTING

LOG-NORMAL

DISTRIBUTION

AND

TESTING

FOR

SUBROUTINE

REQUIRED

AND

APPENDED

I S HI

STO

FOR

************************************************************

CHARACTER

PNE

DIMENSION

A(31) ,AN(104) ,FR(100) ,AV(100) ,FL(100) ,EL(100)

DIMENSION

V(100)

,Sl(100)

,S2(100)

,S3(100)

,S4(100) ,PCR(100)

DIMENSION

R(2)

,T(2)

,P(32)

,PP(2) ,Q(52) ,E(52)

,ALEVEL(6)

************************************************************

THE

NORMAL

PROBABILITIES

ARE

CALCULATED

THE

BASIS

DISCRETE

APPROXIMATION

************************************************************

DATA

A/1.

,6 .0 ,6

.6666,8

,10

.2857 ,12

.2222,14.1818,16.1538

1 18.1333,20 .1176,22 .1052,24.0952,26.0869,28.08,30.0740,32 .0689,

2 34.0645,36 .0606,38 .0571,40.0540,42.0512,44.0487,46.0465,48 .0444,

3 50.0425,52.0408,54 .

0392

,56.

0377

,58 .0363,6 0.

035

0,62.03/

************************************************************

PNE

=CHAR

(12 )

OPEN(UNIT

=6,FILE

='LN.DAT

' , STATUS=' OLD' )

OPEN(UNIT

=2,FILE ='LN.RES' , STATUS='

' )

************************************************************

GENERATING

CLASS

LIMITS

(I)

FOR

ACCUMULATION

************************************************************

READ(6,51)

(ALEVEL(I)

,1 =1,6)

FORMAT

(6A4)

READ(6,*)KEY,N1

READ(6,697)AN(1) ,AN(2),W

697

FORMAT(3F10

.0)

=N1*2

100 K=3,N2 ,2

AN(K)

=AN(K-1)

AN(K+1)

=AN(K)+

100

CONTINUE

************************************************************

STANDS

FOR

THE

NUMBER

CLASS

INTERVALS

6 96

000

001

444

445

446

555

666

Computer Software 159

A CONSTANT TO BE ADDED I N CASE THE

DISTRIBUTION

A 3 -PARAMETER LOG

NORMAL

DISTRIBUTION.

C CAN BE 0 ALSO)

************************************************************

READ(6,696)C

FORMAT(F1 0 . 0)

SM1=0 . 0

SM2=0 . 0

SM3=0 . 0

SM4

=0.0

ESUM

=O.O

CHI

=O.O

************************************************************

) STANDS FOR THE FREQUE

NCIES

************************************************************

DO 8 I =1

,N1

FR(I)

=O.

N=O

READ(6

01)XC

,YC

,Zl,Z2

,Z3

ORMAT

(5F1

0 )

IF(XC.EQ

. 9 9

.99

)GO TO

1005

(KEY . EQ

.1)

GO TO

444

IF(KEY.EQ

.2)

GO TO 44 5

IF(KEY.EQ

GO TO

446

=Zl

GO TO

Z=Z2

GO TO

555

Z=Z3

K=l

Z=Z+C

PROD=Z

************************************************************

DISTRI

BUTES THE OBSERVATIONS I NTO VARI OUS CLASS INTERVALS

************************************************************

DO 9 I =1

,N1

IF(PROD

.GE

.AN(K)

. AND. PROD. L

AN(K+1 ) ) FR (I

FR (I ) +1

=K+2

N=N+1

GO TO 1 00 0

TOTAL=N

************************************************************

COMPUTES THE

MEAN

OF THE LOGS OF Z+C OBSERVATIONS AND

VARIANCE THERE ON .

THIS

BASED ON FREQUENCY D

ISTRIBUTIO

N .

************************************************************

WRI TE ( 2 ,

791)

791

ORMAT

(4X , ' FI TTI NG LOGN

ORMAL

DISTRIBUTIO

N '

WRI TE (2 , 7 9 2 )

(ALEVEL(I)

=1,6)

160 Geostatistics with Applications in Earth Sciences

792

FORMAT(4X

,6A4/)

WRITE

(2,178)

178

FORMAT

(4X

'CLAS

LIMIT

S' , 4X, ' AVG.OF

' ,3 X,

'LOGARI

THMIC

LMTS

' ,

1 7X, 'AVG' , 4X, 'OBS.FREQ' )

DO 14 IP =1,N2 ,2

AV(K)

=(AN(IP)+AN(IP+1))/2.

FLOG

=ALO

G(AN(IP) )

ELOG

=ALOG(AN(IP+1))

FL(K)

=FLOG

EL(K)

=ELOG

AAV(K)

=(FLOG+ELOG)/2

************************************************************

FREQUENCY

DISTRIBUT

ION

AND

OTHER

RELEVANT

STATISTICS

************************************************************

WRI TE(2, 778)

(I P) ,

(I P+1)

,AV(K)

,FLO

,ELOG,AAV(K)

,FR(K)

778

FORMAT(F7

. 0, ' TO' , F7. 0, F10 .1, F10 . 1, ' TO' ,F10. 1 , 2F10. 1 )

14 K=K+1

DO 140 I =1,N1

Sl(I)

=FR(I)*AV(I)

S2(I)

=FR(I)*AAV(I)

S3(I)

=FR(I)*AV(I)**2

S4(I)

=FR(I)*AAV(I)**2

SM1

=SM1+S1(I)

SM2

=SM2+S2(I)

SM3

=SM

3+S

3(I)

140

SM4

=SM4+S4(I)

************************************************************

C X

BAR

AND

VAR1

ARE

MEAN

AND

STANDARD

DEVIATION

BASED

UNTRANSFORMED

OBSERVATIONS

BUT

COMPUTED

FREQUENCY

DISTRIBUTION

APPROACH.

C Z

BAR

AND

VAR2

ARE

MEAN

AND

STANDARD

DEVIATION

BASED

TRANSFORMED

BSERVATIONS

BUT

COMPUTED

FREQUENCY

DISTRIBUTION

APPROACH.

************************************************************

BAR

=SM1/TOTAL-C

BAR

=SM2/TOTAL

VAR1

=SM3/TOTAL-XBAR**2

VAR2

=SM4/T

OTAL

-ZBAR**2

SD1=

RT (VAR1 )

SD2

=SQRT(VAR

************************************************************

CALLING

SUBROUTINE

HISTO

FOR FOR

HISTOGRAM

************************************************************

WRI TE(2, 1111 )PN

IF(KEY

.EQ.1)

3000

IF(KEY

.EQ.3)

3001

WRI TE(2, 2002)

Computer Software 161

************************************************************

SERVED

EXPECTED

FREQ

NCIES

AND

HER

STATI STI

RMA

/ / 25X, 'HI S

TOG

RAM

FOR

GRADE'

/ /

TO 4000

WRITE(2, 2000)

RMA

/ / 25X, ' HIS

TOG

RAM

FOR

ACC

ULA

TION' / /

TO 4000

WRITE(2, 2001)

RMA

/ / 25X, ' HIS

TOG

RAM

FOR

THI

CKN

ESS' / /

ALL

HI STO

(N1

, FR)

2002

3000

2000

3001

2001

4000

1009

978

WRITE (2, 1111)

WRITE(2, 1009)

RMAT

15X,

'OBSER

VED

AND

EXPECTED

FREQ

NCIES'

WRITE(2, 978)

RMAT

(5X,

'OBS

.FREQ'

,3X, '

%AGE'

, 4X,

'PROBABILITY

MTS'

,5X, 'PROW,

3X, 'EXP.

FREQ'

)

107 I =1, N1

PCR(I )=(FR

(I)

/TOTAL

)*100 .

T(1)=(FL(I )

-ZBAR)

/SD2

T(2) =(EL(I )

-ZBAR)

/SD2

IF(T

(1))280,282

,282

280

IF(

-T (

1)-5

.)282,28

1,28

281 T(1 )

=-5.

284

282

IF(

T (1) -5 . )28

4,283,283

283 T(1) =5.

284

IF(T

(2))285,287

,287

285

IF(

-T (2)

-5.)287,286

,286

286 T2=

-5.

289

287

IF(

T (2) -5 .

)289,288,288

288 T(2)=5 .

289

15 L=1,2

R(L)=T(L )

IF(R

(L))802,803

,803

802 R(L)=-R(L )

803 X=R(L)

P (

l)=l.

804

II=2,

III

=II-1

804

P(II

)=(P(III)

(II

) )*

X*X

P (32 )=

.39894228

805

=l , 31

UM+P

(JJ)

805 D=-D

162 Geostatistics with Applications in Earth Sciences

15 PP(L )= .5+S

UM*P(32)

IF(T

(1))20

,21 ,21

IF(T

(2))22,22,23

22 Q

(I)=PP(

1)-PP (2)

23 Q

(I)=PP(

1 )+PP (2 )- 1 .

IF(T

(2))23,2

4,2

24 Q(I )=PP(2 )-PP (1)

50 E(I )=Q(I)*TOTAL

UM=E

UM+E(I

)

IF(E

(I ) .EQ. O.O)GO

1007

1008

CHI=C

HI+(FR

(I)-E

(I)

)* (FR(I ) - E( I )

E( I)

1007 WRITE(2, 998 )FR(I ) ,

PCR

(I ) , T(1) , T(2) , Q(I ) , E(I )

998

AT(5X,F6.0, lX,F8 .

1,F8.

1,'

',F8

1,F8

.3,2X,F8.

1 )

107

NTI

NUE

****

*** **** **

***********

******* *

***********

******** ** *

******

C ESTI

MATIO

RAMETERS

****

*** **** **

***

****

******* *

****

***

******** ** *

****

AR=ZBAR+

VAR2

/2.

EST=EXP(ZAR)

-C

POP=

( (EST**2)* (E

(

VAR2

)-1

.) )

/TOTAL

1=EXP(ZBAR

VAR2

)

WRITE( 2,7 9)TOTAL

WRITE(2, 80) XBAR

WRITE( 2, 81) S

WRITE(2, 82) ZBAR

WRITE( 2, 83) SD2

WRITE(2, 84) E

SUM

WRITE( 2, 85)

WRITE( 2, 86) EST

WRITE( 2, 88) POP

RMAT

(/ / 4X, ' TOT

NUM

BER

OBSER

VATIO

=' F8. 2/)

RMAT

(4X, 'M

EAN

BASED

FREQ

NCY

DISTRIB

UTIO

N =' F8. 2/)

RMAT

(4X, ' S

TAN

ATIO

BASED

FREQ

.DIST

. . =' F8. 2/)

RMAT

(4X, 'M

EAN

BASED

AVE

LOGS

CLS

=' F8. 2/)

RMAT

(4X, ' S

TAN

ATIO

BASED

ABO

EAN

=' F8. 2/)

RMAT

(4X, '

SUM

EXPECTED

FREQ

NCIES

='F 8.2/)

FORMAT

(4X, 'CO

MPUTED

CHI-SQ

UARE

=' F8. 2/ )

FORMAT

(4X,

'KRIGES

STIMATE

=' F8. 2/ )

RMAT

(4X, ' PO

LATI

VAR

IAN

CE =' F8. 2/ )

1111

RMAT

(A)

************************************************************

C S

UBRO

UTI

HISTO

FOR FOR

HISTOGRAM

************************************************************

UBRO

UTI

HISTO

(

FR)

Computer Software 163

CHARACTER*1 ICHAR1

,ICHAR2,

B(104

,100)

DIMENSION

FR(104)

DATA

ICHAR1,

ICHAR2

" , /

DO 4 I =1

,N1

IFR

=FR(I)

DO 3 J =1

0 0

IF(J

-IFR)1

,1,2

B(I,J)

=ICHAR1

GO TO 3

2 B (I , J )= I CHAR2

3 CONTINUE

4 CONTINUE

WRI TE (2 , 5 )

FORMAT

(7X , ' 0 ' , 9X ,

'10

,8X

, ' 2 0 ' , 8X,

'30'

,8X

'40

' , 8X,

'50

' ,

8X,

'60'

,8X

, 8X, ' 8

)

WRI TE ( 2 , 6 )

6 FORMAT(lX

,80(lH

-))

DO 9 I =1

,N1

WRI TE (2 , 8 )

(B(I

,J)

, J =

8 0 )

WRI TE

(2,7)

(I)

, (B

, J =l ,

80)

7 FORMAT(lX

,F6

,2X,100A1)

WRI TE (2 , 8 )

(B(I

,J)

,80)

8 FORMAT( 9X

0A1)

9 CONTINUE

RETURN

END

PROGRAM

AR.FOR

************************************************************

C ESTIMATION BY AUTO

-REGRESSIVE

METHOD.

C THE AR

COEFFICIENTS

ARE OBTAINED BY SOLVING THE YULE-

C WALKER SCHEME.

THIS

PROGRAM

WORKS FOR ANY ORDER.

C NP=LLL

-1

GIVES

THE

DESIRED

ORDER.

************************************************************

THIS

PRO

GRAM

COMPUTES THE AUTO-CORRELATION C

OEFFICIENTS

ALSO

************************************************************

C SUBROUTINE S REQUIRED ARE YWS. FOR (WHICH INCLUDES

CRP.

FOR

C CDOT•FOR) AND SRL • FOR)

************************************************************

C XX

(I)

ARE THE LOG TRANSFORMED

INPUT

VALUES WHICH ARE

C DETRENDED SUBSEQUENTLY BY REMOVING THE

MEAN

C ACV(N+

ARE THE AUTOCOVARIANCE

COEFFICIENTS

FFPE

(N+

ARE THE

FINAL

PREDICTION

ERROR

C C

OEFFICIENTS

USED TO DETERMINE THE

MAXIMUM

C LENGTH OF THE OPERATOR ACCORDING TO AKAIKE ' S

CRITERION

C ALPHA

(I)

=- GG

(I)

ARE THE

PREDICTI

ON ERROR

COEFFICIENT

************************************************************

C READS ONE RECORD PER CARD I .

, XE, YN,

(ACCUMULATION)

164

Geostatistics with Applications in Earth Sciences

C Z2(GRADE) ,Z3(THICKNESS)

KEY

= 1

MEANS

LOG

TRANSFORMATION

AND

DEVIATIONS

TAKEN

FROM

THE

MEAN

KEY

MEANS

NO. LOG.

TRANSFORMATION

BUT

DEVIATIONS

TAKEN

FROM

THE

MEAN

DIMENSION WP (5 0 )

,R(50)

DIMENSION

XX(1200)

,AC

V(1200)

,FFPE(1200)

,GG(1200)

,XXC(1200)

DIMENSION XM(1200) , XXT (12 00 ) , ALPHA (1 0 )

DIMENSION ALEVEL(10)

COMMON

NX,

************************************************************

OPEN(UNIT =6

,FILE

='AR.DAT ' , STATU

S='

OLD' )

OPEN(UNIT =2

,FILE

='AR.RES

' , STATU

S='

' )

************************************************************

PNEW

=CHAR

(12

)

3 WRI TE (2 , 999 9 ) PN

9999

FORMAT

(A)

************************************************************

READ(6,22)

(ALEVEL(I)

=1,10)

FORMAT

(10A4)

WRI TE (2 , 33 ) (ALEVEL(I) , 1= 1 , 10 )

33 FORMAT(lX,10A4)

************************************************************

KEY

MEANS

COMPUTATIONS

UNTRANSFORMED

OBSERVATIONS

OTHER

WISE ON

LOG

-TRANSFORMED

OBSERVATIONS

KKK

THE UPPER LIMIT OF LLL . EG:

LLL=5

THEN

ORDER

C THE PROCESS

************************************************************

READ(6, *) KEY,

KKK

NX=O

SUM

=O.O

************************************************************

READS

DATA

PERFORMS

LOGARITHMIC

TRANSFORMATIO

NEEDED

AND

COMPUTES

MEAN

************************************************************

5 READ(6

,40)XE,YN

,Zl,Z2,Z3

40 FORMAT(5F10

.2)

IF(XE

.EQ

.99.99)GO

TO 66

NX=NX+1

XX(NX)

=Zl

XM(NX)

=XX(NX)

IF(KEY .EQ .O)GO TO

223

XX(NX) =ALOG(XM(NX))

223

SUM

=SUM+XX(NX)

TO 5

66 AN

=NX

=SUM/AN

C498

772

776

497

555

777

778

Computer Software 165

*******************************************************

TRENDS

DATA

BY SUBTRACTING

MEAN

FROM

EACH

DATA

POINT

************************************************************

I=1,

XX(I)=XX(I)

-AM

************************************************************

(I)

ARE THE DETRE

NDED

VALUES :

CONTL=1 THEN XX

(I)

ARE LOG.

TRANSFORMED

AND

TRENDED

VALUES WHI LE

XM(I) ARE THE ORIGINAL VALUES .

OTHER

WISE XM(I) ARE THE

ORIGINAL

VALUES

AND

(I)

ARE THE

DETRENDED

VALUES .

************************************************************

(KEY. EQ. 1 )

TO 77 6

WRI TE( 2 , 498 )

FORMAT

(1X, ' PRI NTS

ACCUMULATIO

DATA

WRI TE (2 , 5 55 ) (XM(I)

=1,10)

WRI TE (2 ,

772

)

FORMAT

(/14X,

'MEAN

DATA

='F8.2)

778

WRI TE

(2,497)

FORMAT

(/1X,

'PRINTS

DATA

LOG

.TRANS.DEV.TAKEN

DATA'I)

WRI TE (2 , 555 ) (XM(I)

,XX(I)

, 1= 1 , 10 )

FORMAT(4X

,10F8.3)

WRI TE

(2,777)AM

FORMAT

(/2X,

'MEAN

TRANSFORMED

DATA

F8 . 2)

************************************************************

COMPUTES

AUTO-CORRELATION

THROUGH

SUB-ROUTINE

'SRL.

FOR'

************************************************************

WRI TE (2 , 41 )

FORMAT

(/2X

, 'AUTO.CORRELATION COEFFICIENT FOR

ZERO

LAG'I)

************************************************************

IF(KEY.EQ.O.O)GO

444

************************************************************

CALL SRL (XX, CPO, RZERO,WP , R)

************************************************************

445

************************************************************

444

CALL SRL (XM,CPO,RZERO, WP,R)

************************************************************

445

WRI TE (2 , 42 )CPO, RZERO

42 FORMAT(3X

,F10

,4X,F8.2

WRI TE (2 , 43 )

FORMAT

(2X, ' AUTOCORRELATI ON COEFFICIENTS FOR

LAGS

1 TO

M')

IF(KEY

.EQ .O .O)GO TO

446

************************************************************

CALL SRL (XX, CPO, RZERO,WP , R)

************************************************************

447

************************************************************

446

CALL SRL(XM,CPO,RZERO, WP,R)