Geology Reference
In-Depth Information
C Set switch IC to zero for particular integral.
IC=0
C Set initial value of YYP vector at beginning of region.
IF(M.EQ.1)GO TO 59
YYP(1,1)=INITYP(1,M)
YYP(2,1)=INITYP(2,M)
GO TO 60
59 YYP(1,1)=0.D0
CALL DERIV(R(1),RHOP,N1,C,R,RHO,M1)
YYP(2,1)=-RHOP*E(1)*E(1)/(64.D0*RHO(1))
60 CONTINUE
C Begin Runge-Kutta extension of the particular integral through region.
C Construct right-hand side.
DO 61 I=1,N1
IF(I.EQ.1.AND.M.EQ.1)GO TO 62
RHS(I)=3.D0*E(I)*E(I)*(1.D0-RHO(I)/RHOB(I))*
1 (4.D0+6.D0*ETA(I)+3.D0*ETA(I)*ETA(I))/(4.D0*R(I)*R(I))
2 -E(I)*E(I)*(7.D0+5.D0*ETA(I))*ETA(I)/(2.D0*R(I)*R(I))
GO TO 61
62 RHS(1)=0.D0
61 CONTINUE
C Do Runge-Kutta extension.
CALL RK(YYP,IC,RHS,RHO,RHOB,E,NIS,N1,C,R,M1)
C Store starting value for next region.
INITYP(1,MP1)=YYP(1,N1)
INITYP(2,MP1)=YYP(2,N1)
C Store YYP in YYPI.
DO 63 I=1,N1
J=NK(M)+I
YYPI(1,J)=YYP(1,I)
YYPI(2,J)=YYP(2,I)
63 CONTINUE
52 CONTINUE
C Apply boundary condition to YYCI and YYPI.
C determine amplitude DD of YYCI at the surface.
DD=-RI(NT)*YYPI(2,NT)-4.D0*YYPI(1,NT)
1 -5.D0*FRD*(EI(NT)-5.D0*FRD/4.D0)/4.D0
DD=DD/(RI(NT)*YYCI(2,NT)+4.D0*YYCI(1,NT))
C Find YYI and flattening FI.
DO 64 I=1,NT
YYI(1,I)=DD*YYCI(1,I)+YYPI(1,I)
YYI(2,I)=DD*YYCI(2,I)+YYPI(2,I)
FI(I)=EI(I)+5.D0*EI(I)*EI(I)/42.D0-4.D0*YYI(1,I)/7.D0
64 CONTINUE
C Write out inner core data file, innerc.dat.
C Write out number of model points.
WRITE(3,65)NM(1)
65 FORMAT(1X,I9)
C Write out density at the bottom of the outer core.
NM1P1=NM(1)+1
WRITE(3,66)RHOI(NM1P1)
66 FORMAT(1X,D15.8)
WRITE(3,67)(RI(I),RHOI(I),RHOBI(I),FI(I),I=1,NM(1))
67 FORMAT(1X,4D15.8)
C Write out radius, density, mean density, 1/f and kappa.
C Write out headings.
WRITE(2,68)
WRITE(6,68)
68
FORMAT(//10X,'Radius',8X,'Density',6X,'Mean Density',6X,'1/f',
1 11X,'kappa' )
WRITE(2,69)
Search WWH ::




Custom Search