Geology Reference
In-Depth Information
C
IMPLICIT DOUBLE PRECISION(A-H,O-Z)
C Specify value of pi.
PI=3.141592653589793D0
C Set angular frequency of Earth's rotation.
WE=7.292115D-5
C Calculate number of hours in a semi-sidereal day.
HSSD=PI/(WE*3600.D0)
C Set flattening to hydrostatic value at the core-mantle boundary
C for Earth model Cal8.
FLAT=1.D0/394.03D0
C Calculate square of the eccentricity from flattening.
ESQD=1.D0-(1.D0-FLAT)*(1.D0-FLAT)
C Enter azimuthal number M and degree N.
10 CONTINUE
WRITE(6,11)
11 FORMAT(1X,'Enter M and N')
READ(5,*)M,N
C Choose between searching a range of Coriolis frequencies or
C iterating on two trial periods.
12 CONTINUE
WRITE(6,13)
13 FORMAT(1X,'Enter -1 for search a range, 0 for iterate')
READ(5,*)ISW
IF(ISW.GE.0) GO TO 14
C Enter bottom and top of range of Coriolis frequencies
C and number of increments.
21 CONTINUE
WRITE(6,15)
15 FORMAT(1X,'Enter bottom and top of Coriolis frequency range',1X,
1 'and number of increments')
READ(5,*)BOT,TOP,NT
C Find increment in period.
DSIG=(TOP-BOT)/DFLOAT(NT)
NTP1=NT+1
C Find residual ERR1 of secular equation at each of trial periods.
DO 16 I=1,NTP1
C Find Coriolis frequency.
SIG1=(BOT+DFLOAT(I-1)*DSIG)
Z=SIG1*DSQRT((1.D0-ESQD)/(1.D0-ESQD*SIG1*SIG1))
C Find error in satisfying secular equation.
ERR1=ERR(Z,M,N,SIG1)
TS=12.D0/SIG1
T=HSSD/SIG1
WRITE(6,17)SIG1,TS,ERR1,T
16 CONTINUE
17 FORMAT(1X,F7.4,2X,'Sidereal period=',F8.4,2X,'Error=',E10.3,2X,
1 'Period',F7.3,1X,'hours')
C Decide whether to do new search or to iterate.
GO TO 12
C Enter two trial periods for iteration.
14
CONTINUE
WRITE(6,18)
18
FORMAT(1X,'Enter trial periods')
READ(5,*)T1,T2
SIG1=12.D0/T1
Z=SIG1*DSQRT((1.D0-ESQD)/(1.D0-ESQD*SIG1*SIG1))
ERR1=ERR(Z,M,N,SIG1)
SIG2=12.D0/T2
Z=SIG2*DSQRT((1.D0-ESQD)/(1.D0-ESQD*SIG2*SIG2))
ERR2=ERR(Z,M,N,SIG2)
Search WWH ::
Custom Search