Geology Reference
In-Depth Information
CALL QUAD(INT,EX,Q,NIS,W,X,N1,C,R,M,M1)
C Store starting value for next region.
INTI(2,MP1)=INT(N1,M)
C Form first estimate of ETA and store in ETAI.
DO 31 I=1,N1
IF(I.EQ.1.AND.M.EQ.1)GO TO 32
ETA(I)=5.D0*INT(I,M)/(RHOB(I)*(R(I)**5))
ETA(I)=ETA(I)*ETA(I)-1.D0
GO TO 33
32 ETA(I)=0.D0
C Store ETA.
33 J=NK(M)+I
ETAI(J)=ETA(I)
31 CONTINUE
C Perform first estimate of figure parameter e.
C Set exponent in integral.
EX=-1.D0
C Initialize starting value of integral INT.
INT(1,M)=INTI(3,M)
C Form integrand function Q.
DO 34 I=1,N1
Q(I)=ETA(I)
34 CONTINUE
C Perform quadrature.
CALL QUAD(INT,EX,Q,NIS,W,X,N1,C,R,M,M1)
C Store starting value for next region.
INTI(3,MP1)=INT(N1,M)
C Form first estimate of e(r)/e(0) and store in EI.
DO 35 I=1,N1
IF(I.EQ.1.AND.M.EQ.1)GO TO 36
E(I)=DEXP(INT(I,M))
GO TO 37
36 E(I)=1.D0
C Store e(r)/e(0) in EI.
37 J=NK(M)+I
EI(J)=E(I)
35 CONTINUE
22 CONTINUE
C Find total number of model points.
NT=NK(4)+NM(4)
C Find surface value of Froude number or dimensionless parameter m(d).
FRD=3.D0*WE*WE/(4.D0*PI*G*RHOBI(NT))
C Find surface value e(d).
C Generate coefficients of quadratic.
BB=-(14.D0+7.D0*ETAI(NT)+6.D0*FRD)/4.D0
CC=FRD*(35.D0/8.D0+5.D0*FRD/6.D0)
ED=-(BB+DSQRT(BB*BB-4.D0*CC))/2.D0
C Find e(0).
E0=ED/EI(NT)
C Find e(r).
DO 38 I=1,NT
EI(I)=E0*EI(I)
38 CONTINUE
C Begin second pass on integration of the Clairaut equation
C by the Radau transformation.
DO 39 M=1,4
MP1=M+1
C Fix number of points for region M.
N1=NM(M)
C Restore interpolation matrix for region M.
DO 40 I=1,N1
Search WWH ::
Custom Search