C                    PROGRAMME SOURCES.FOR
C SOURCES.FOR solves the sixth-order spheroidal system for the four singular sources.
C For a given input Earth model, the solutions are obtained by integrating through
C the inner core, fluid outer core, mantle and crust for singular sources located
C at the focal depth in either the mantle or crust. The four source solutions are given,
C as output, in the files source1.dat, source2.dat, source3.dat and source4.dat.
C As required for the calculation of changes in the inertia tensor in epicentral 
C co-ordinates, the radial andtransverse coefficients u and v are integrated 
C over radius for the four source solutions. The results of these integrations are 
C output in the file results.dat.
C 
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION R(100),RHO(100),GZERO(100),RI(300),RHOI(300),
     1  GZEROI(300),ENAME(10),NM(4),NI(4),NK(4),B(98,198),C(100,100),
     1  Y(6,6),YSCAL(6,3),YSCAL2(6,3),YA(6),A(6,6),AS(6,6),CAUG(6,13),
     1  RM1(6,3),RM2(6,3),YP(6,6),YIC(6,3),Y1(6,6),Y2(6,6),DELTA(6,3),
     2  YINFD(6,6),YINFDS(6,6),RAD(1000),FS(6,3,1000),G1(6),G2(6),G3(6),
     2  G4(6),ABC(6,3),PROD(6,3),AA(3,3),AAS(3,3),RHOIC(725),
     2  B1(3),B2(3),B3(3),B4(3),ROC(15),RHOOC(15),YOC(6,15),RM(725),
     3  RHOM(725),YYM(6,6,725),RC(75),RHOC(75),YYC(6,6,65),YF1(6),
     3  YF2(6),YF3(6),YF4(6),UV(6),UV1(6),UV2(6),UV3(6),UV4(6),
     3  CAUG2(3,7),SIC(8,725),SOC(8,15),SM(8,725),SC(8,75),VOC(8,15),
     4  WOC(15),SX(8),RES(8,5),DISC(4,2)
      DOUBLE PRECISION MU(100),LAMBDA(100),MUI(300),LAMBDAI(300),
     1  K1(6,6),MUFD
      CHARACTER*7 ENAME
      CHARACTER*20 EMODEL
C Enter Earth model file name.
      WRITE(6,10)
 10   FORMAT(1X,'Type in Earth model file name.')
      READ(5,11) EMODEL
 11   FORMAT(A20)
C Open Earth model file.
      OPEN(UNIT=1,FILE=EMODEL,STATUS='OLD')
C Open free solutions files.
      OPEN(UNIT=3,FILE='fs1.dat',STATUS='UNKNOWN')
      OPEN(UNIT=4,FILE='fs2.dat',STATUS='UNKNOWN')
      OPEN(UNIT=7,FILE='fs3.dat',STATUS='UNKNOWN')
C Open source solution files.
      OPEN(UNIT=8,FILE='source1.dat',STATUS='UNKNOWN')
      OPEN(UNIT=9,FILE='source2.dat',STATUS='UNKNOWN')
      OPEN(UNIT=10,FILE='source3.dat',STATUS='UNKNOWN')
      OPEN(UNIT=11,FILE='source4.dat',STATUS='UNKNOWN')
      OPEN(UNIT=12,FILE='results.dat',STATUS='UNKNOWN')
C Set angular frequency of Earth's rotation (WGS 84).
      WE=7.292115D-5
C Calculate square of angular frequency of Earth's rotation.
      WES=WE*WE
C Set value of pi.
      PI=3.141592653589793D0
C Set value of universal constant of gravitation (CODATA 2006).
      G=6.67428D-11
C Set maximum relative error tolerance for integrations.
      EPS=1.D-5
C Set minimum starting radius for variable Earth properties.
      RMIN=1200.D0
C Set maximum dimensions for interpolation.
      M1=100
      M2=198
      M3=98
C Read in and write out Earth model.
C Read in and write out headers.
C Earth model name.
      READ(1,12)(ENAME(I),I=1,10)
 12   FORMAT(10A7)
      WRITE(6,13)(ENAME(I),I=1,10)
 13   FORMAT(11X,10A7)
C Read in next line of header.
      READ(1,14)NN
 14   FORMAT(I10)
C Read in number of model points and number of integration steps.
      READ(1,15)(NM(I),NI(I),I=1,4)
 15   FORMAT(8I10)
      WRITE(6,16)
 16   FORMAT(//16X,'Number of model points and number of integration',
     1 ' steps.'/)
      WRITE(6,15)(NM(I),NI(I),I=1,4)
      WRITE(6,17) 
 17   FORMAT(//5X,'Radius',6X,'Rho',4X,'Lambda',7X,'Mu',6X,'Gzero')
      WRITE(6,18)
 18   FORMAT(6X,'(km)',5X,'(gm/cc)',2X,'(kbars)',4X,'(kbars)',2X,
     1 '(cm/sec/sec)'/)
C Read in Earth model, inner core, outer core, mantle and crust.
      K=0
      DO 19 M=1,4
        N1=NM(M)
        READ(1,20)(R(I),RHO(I),LAMBDA(I),MU(I),GZERO(I),I=1,N1)
 20     FORMAT(1X,F10.1,F10.2,F10.1,F10.1,F10.1)
        DO 21 I=1,N1
C Scale Earth model to SI values and store.
          J=K+I
          RI(J)=R(I)*1.D3
          RHOI(J)=RHO(I)*1.D3
          LAMBDAI(J)=LAMBDA(I)*1.D8
          MUI(J)=MU(I)*1.D8
          GZEROI(J)=GZERO(I)*1.D-2
C Write Earth model out to screen.
          WRITE(6,20)R(I),RHO(I),LAMBDA(I),MU(I),GZERO(I)
 21     CONTINUE
        NK(M)=K
        K=K+N1
 19   CONTINUE
C Store surface radius.
      D=RI(NK(4)+NM(4))
C Begin integration through the inner core.
C Set up interpolation for the inner core.
      N1=NM(1)
      N2=2*N1-2
      N3=N1-2
C Put inner core values in active locations.
      DO 22 I=1,N1
        J=NK(1)+I
        R(I)=RI(J)
        RHO(I)=RHOI(J)
        MU(I)=MUI(J)
        LAMBDA(I)=LAMBDAI(J)
        GZERO(I)=GZEROI(J)
 22   CONTINUE
C Store density just inside inner core.
      RHOAM=RHO(N1)
C Calculate gravity gradient at geocentre.
      GAMMA=(4.D0*PI*G*RHO(1)-2.D0*WES)/3.D0
C Set gravity at minimum radius.
      GZERO(1)=GAMMA*RMIN
C Construct interpolation matrix.
      CALL SPMAT(N1,N2,N3,C,R,B,M1,M2,M3)
C Set degree to 2.
      N=2
      AN=DFLOAT(N)
C Set ANGS and COR to zero for dc solutions.
      ANGS=0.D0
      COR=0.D0
C Begin power series expansions for fundamental solutions for N=2.
C Set initial values of fundamental solutions.
C Set initial values of first fundamental solution.
      Y(1,1)=1.D0
      Y(2,1)=2.D0*(AN-1.D0)*MU(1)
      Y(3,1)=1.D0/AN
      Y(4,1)=Y(2,1)/AN
      Y(5,1)=4.D0*PI*G*RHO(1)/AN
      Y(6,1)=0.D0
C Set initial values of second fundamental solution.
      Y(1,2)=0.D0
      Y(2,2)=0.D0
      Y(3,2)=0.D0
      Y(4,2)=0.D0
      Y(5,2)=1.D0/AN
      Y(6,2)=1.D0
C Set initial values of third fundamental solution.
      Y(3,3)=P2(AN,LAMBDA(1),MU(1))/P1(AN,LAMBDA(1),MU(1))
      Y(2,3)=-Q1(AN,LAMBDA(1),MU(1))*Y(3,3)+Q2(AN,LAMBDA(1),MU(1))
      Y(1,3)=-AN*Y(3,3)+1.D0/MU(1)
      Y(4,3)=1.D0
      Y(5,3)=2.D0*PI*G*RHO(1)*((AN+3.D0)*Y(1,3)-AN*(AN+1.D0)*Y(3,3))/
     1 (2.D0*AN+3.D0)
      Y(6,3)=(AN+2.D0)*Y(5,3)-4.D0*PI*G*RHO(1)*Y(1,3)
C Calculate coefficients of second terms in series expansions.
C Calculate coefficients of second terms
C for first fundamental solution.
      YSCAL(3,1)=RHO(1)*((3.D0-AN)*GAMMA+ANGS+2.D0*WES-COR/AN)/
     1 P1(AN,LAMBDA(1),MU(1))
      YSCAL(2,1)=-Q1(AN,LAMBDA(1),MU(1))*YSCAL(3,1)
      YSCAL(1,1)=-AN*YSCAL(3,1)
      YSCAL(4,1)=0.D0
      YSCAL(5,1)=2.D0*PI*G*RHO(1)*((AN+3.D0)*YSCAL(1,1)
     1 -AN*(AN+1.D0)*YSCAL(3,1))/(2.D0*AN+3.D0)
      YSCAL(6,1)=(AN+2.D0)*YSCAL(5,1)-4.D0*PI*G*RHO(1)*YSCAL(1,1)
C Calculate coefficients of second terms
C for second fundamental solution.
      YSCAL(3,2)=RHO(1)/P1(AN,LAMBDA(1),MU(1))
      YSCAL(2,2)=-Q1(AN,LAMBDA(1),MU(1))*YSCAL(3,2)
      YSCAL(1,2)=-AN*YSCAL(3,2)
      YSCAL(4,2)=0.D0
      YSCAL(5,2)=2.D0*PI*G*RHO(1)*((AN+3.D0)*YSCAL(1,2)
     1 -AN*(AN+1.D0)*YSCAL(3,2))/(2.D0*AN+3.D0)
      YSCAL(6,2)=(AN+2.D0)*YSCAL(5,2)-4.D0*PI*G*RHO(1)*YSCAL(1,2)
C Generate matrix for evaluation of second terms in power series
C expansion of third fundamental solution.
      CALL MATRIX(A,RHO,LAMBDA,MU,PI,G,AN,WE)
C Set right-hand side for calculation of second terms
C in third solution.
      YA(1)=0.D0
      YA(2)=RHO(1)*(-(4.D0*GAMMA+ANGS+2.D0*WES)*Y(1,3)
     1 +(AN*(AN+1.D0)*GAMMA+COR)*Y(3,3)-Y(6,3))
      YA(3)=0.D0
      YA(4)=RHO(1)*((GAMMA+COR/(AN*(AN+1.D0)))*Y(1,3)
     1 -(ANGS-COR/(AN*(AN+1.D0)))*Y(3,3)-Y(5,3))
      YA(5)=0.D0
      YA(6)=0.D0
C Store matrix A in AS.
      DO 23 I=1,6
        DO 24 J=1,6
          AS(I,J)=A(I,J)
 24     CONTINUE
 23   CONTINUE
C Solve for second terms in third solution.
      CALL LINSOL(A,YA,6,CAUG,DET,6,13)
      YSCAL(1,3)=YA(1)
      YSCAL(2,3)=YA(2)
      YSCAL(3,3)=YA(3)
      YSCAL(4,3)=YA(4)
      YSCAL(5,3)=YA(5)
      YSCAL(6,3)=YA(6)
C Set right-hand side for calculation of third terms
C in first solution.
      YA(1)=0.D0
      YA(2)=RHO(1)*(-(4.D0*GAMMA+ANGS+2.D0*WES)*YSCAL(1,1)
     1 +(AN*(AN+1.D0)*GAMMA+COR)*YSCAL(3,1)-YSCAL(6,1))
      YA(3)=0.D0
      YA(4)=RHO(1)*((GAMMA+COR/(AN*(AN+1.D0)))*YSCAL(1,1)
     1 -(ANGS-COR/(AN*(AN+1.D0)))*YSCAL(3,1)-YSCAL(5,1))
      YA(5)=0.D0
      YA(6)=0.D0
C Copy stored matrix AS into matrix A.
      DO 25 I=1,6
        DO 26 J=1,6
          A(I,J)=AS(I,J)
 26     CONTINUE
 25   CONTINUE
C Solve for third terms in first solution.
      CALL LINSOL(A,YA,6,CAUG,DET,6,13)
      YSCAL2(1,1)=YA(1)
      YSCAL2(2,1)=YA(2)
      YSCAL2(3,1)=YA(3)
      YSCAL2(4,1)=YA(4)
      YSCAL2(5,1)=YA(5)
      YSCAL2(6,1)=YA(6)
C Set right-hand side for calculation of third terms
C in second solution.
      YA(1)=0.D0
      YA(2)=RHO(1)*(-(4.D0*GAMMA+ANGS+2.D0*WES)*YSCAL(1,2)
     1 +(AN*(AN+1.D0)*GAMMA+COR)*YSCAL(3,2)-YSCAL(6,2))
      YA(3)=0.D0
      YA(4)=RHO(1)*((GAMMA+COR/(AN*(AN+1.D0)))*YSCAL(1,2)
     1 -(ANGS-COR/(AN*(AN+1.D0)))*YSCAL(3,2)-YSCAL(5,2))
      YA(5)=0.D0
      YA(6)=0.D0
C Copy stored matrix AS into matrix A.
      DO 27 I=1,6
        DO 28 J=1,6
          A(I,J)=AS(I,J)
 28     CONTINUE
 27   CONTINUE
C Solve for third terms in second solution.
      CALL LINSOL(A,YA,6,CAUG,DET,6,13)
      YSCAL2(1,2)=YA(1)
      YSCAL2(2,2)=YA(2)
      YSCAL2(3,2)=YA(3)
      YSCAL2(4,2)=YA(4)
      YSCAL2(5,2)=YA(5)
      YSCAL2(6,2)=YA(6)
C Set right-hand side for calculation of third terms
C in third solution.
      YA(1)=0.D0
      YA(2)=RHO(1)*(-(4.D0*GAMMA+ANGS+2.D0*WES)*YSCAL(1,3)
     1 +(AN*(AN+1.D0)*GAMMA+COR)*YSCAL(3,3)-YSCAL(6,3))
      YA(3)=0.D0
      YA(4)=RHO(1)*((GAMMA+COR/(AN*(AN+1.D0)))*YSCAL(1,3)
     1 -(ANGS-COR/(AN*(AN+1.D0)))*YSCAL(3,3)-YSCAL(5,3))
      YA(5)=0.D0
      YA(6)=0.D0
C Copy stored matrix AS into matrix A and add 2.0 to diagonal of A.
      DO 29 I=1,6
        DO 30 J=1,6
          A(I,J)=AS(I,J)
          IF(I.EQ.J)A(I,J)=A(I,J)+2.D0
 30     CONTINUE
 29   CONTINUE
C Solve for third terms in third solution.
      CALL LINSOL(A,YA,6,CAUG,DET,6,13)
      YSCAL2(1,3)=YA(1)
      YSCAL2(2,3)=YA(2)
      YSCAL2(3,3)=YA(3)
      YSCAL2(4,3)=YA(4)
      YSCAL2(5,3)=YA(5)
      YSCAL2(6,3)=YA(6)
C Calculate terms in power series expansions at minimum radius RMIN.
      XV=RMIN
      XVS=XV*XV
      XV3=XV*XVS
      XV4=XVS*XVS
      DO 31 I=1,6
        DO 32 J=1,3
          RM1(I,J)=XVS*YSCAL(I,J)
          RM2(I,J)=XV4*YSCAL2(I,J)
 32     CONTINUE
 31   CONTINUE
C Write out terms in power series expansions.
      DO 33 J=1,3
        WRITE(6,34)J
 34     FORMAT(/5X,'Fundamental solution number ',I3)
        WRITE(6,35)
 35     FORMAT(/6X,'first term',10X,'second term',9X,'third term')
        WRITE(6,36)(Y(I,J),RM1(I,J),RM2(I,J),I=1,6)
 36     FORMAT(5X,D15.8,5X,D15.8,5X,D15.8)
 33   CONTINUE
C Set last three columns of arrays Y(6,6) and YP(6,6) to zero
C for inner core.
      DO 37 I=1,6
        DO 38 J=4,6
          Y(I,J)=0.D0
          YP(I,J)=0.D0
 38     CONTINUE
 37   CONTINUE
C Put inner core solutions at minimum radius in array YIC(6,3) and
C initialize free solutions array in y-variables, FS(6,3,1), at the geocentre.
      DO 39 I=1,6
        DO 40 J=1,3
          YIC(I,J)=Y(I,J)
          FS(I,J,1)=0.D0
          IF(J.EQ.1.AND.(I.EQ.2.OR.I.EQ.4)) FS(I,J,1)=Y(I,J)
 40     CONTINUE
 39   CONTINUE
C Write out starting values of free solutions in y-variables at the geocentre.
      RAD(1)=0.D0
      RHOIC(1)=RHO(1)
      WRITE(3,42)RAD(1),(FS(I,1,1),I=1,6)
      WRITE(4,42)RAD(1),(FS(I,2,1),I=1,6)
      WRITE(7,42)RAD(1),(FS(I,3,1),I=1,6)
C Store minimum radius and density at minimum radius.
      RAD(2)=XV*1.D-3
      CALL INTPL(XV,RHOR,N1,C,R,RHO,M1)
      RHOIC(2)=RHOR
C Convert solutions back to y-variables.
      ANM2=AN-2.D0
      DO 41 J=1,3
        IF(J.EQ.3) ANM2=AN
        CONST=XV**ANM2
        FS(2,J,2)=CONST*(Y(2,J)+RM1(2,J)+RM2(2,J))
        FS(4,J,2)=CONST*(Y(4,J)+RM1(4,J)+RM2(4,J))
        CONST=CONST*XV
        FS(1,J,2)=CONST*(Y(1,J)+RM1(1,J)+RM2(1,J))
        FS(3,J,2)=CONST*(Y(3,J)+RM1(3,J)+RM2(3,J))
        FS(6,J,2)=CONST*(Y(6,J)+RM1(6,J)+RM2(6,J))
        CONST=CONST*XV
        FS(5,J,2)=CONST*(Y(5,J)+RM1(5,J)+RM2(5,J))
 41   CONTINUE
C Write out free solutions at minimum radius, RMIN.
      WRITE(3,42)RAD(2),(FS(I,1,2),I=1,6)
      WRITE(4,42)RAD(2),(FS(I,2,2),I=1,6)
      WRITE(7,42)RAD(2),(FS(I,3,2),I=1,6)
 42   FORMAT(F7.2,6D14.6)
      IRAD=2
C Begin Runge-Kutta integration for the inner core.
C Find maximum ratio of second coefficient to first term
C in power series expansions.
      CALL REL(ERRMAX,YSCAL,YIC,N)
C Set initial stepsize so that fourth-order method would have relative error
C bound EPS.
      H=(EPS/ERRMAX)**(0.2D0)
C Find solutions and derivatives at minimum radius RMIN
C by power series expansions.
      DO 43 I=1,6
        DO 44 J=1,3
C Find solutions.
          Y(I,J)=Y(I,J)+XVS*YSCAL(I,J)+XV4*YSCAL2(I,J)
C Find derivatives.
          YP(I,J)=2.D0*XV*YSCAL(I,J)+4.D0*XV3*YSCAL2(I,J)
 44     CONTINUE
 43   CONTINUE
C Set inner core index.
      IR=1
 45   CONTINUE
C Calculate derivatives at beginning of current step.
      CALL YPRIME(XV,Y,A,K1,N1,C,R,RHO,MU,LAMBDA,GZERO,N,PI,G,
     1 IR,RMIN,ANGS,COR,WES)
C Initialize function values at beginning of integration step.
      DO 46 I=1,6
        DO 47 J=1,6
          Y1(I,J)=Y(I,J)
          Y2(I,J)=Y(I,J)
 47     CONTINUE
 46   CONTINUE
C Integrate solution forward one full step.
      CALL RK4(XV,Y1,A,K1,N1,C,R,RHO,MU,LAMBDA,GZERO,N,PI,G,H,IR,
     1 RMIN,ANGS,COR,WES)
C Integrate solution forward two half-steps.
C Halve step.
      HH=0.5D0*H
C Reset radius.
      XV=XV-H
      CALL RK4(XV,Y2,A,K1,N1,C,R,RHO,MU,LAMBDA,GZERO,N,PI,G,HH,IR,
     1 RMIN,ANGS,COR,WES)
      CALL YPRIME(XV,Y2,A,K1,N1,C,R,RHO,MU,LAMBDA,GZERO,N,PI,G,IR,
     1 RMIN,ANGS,COR,WES)
      CALL RK4(XV,Y2,A,K1,N1,C,R,RHO,MU,LAMBDA,GZERO,N,PI,G,HH,IR,
     1 RMIN,ANGS,COR,WES)
C Find error matrix and maximum relative error.
      DO 48 I=1,6
        DO 49 J=1,3
          DELTA(I,J)=Y1(I,J)-Y2(I,J)
 49     CONTINUE
 48   CONTINUE
C Copy Y2(6,6) into inner core solution array YIC(6,3).
      DO 50 I=1,6
        DO 51 J=1,3
          YIC(I,J)=Y2(I,J)
 51     CONTINUE
 50   CONTINUE
C Find maximum relative error.
      CALL REL(ERRMAX,DELTA,YIC,N)
C Scale maximum relative error by specified relative error tolerance.
      ERRMAX=ERRMAX/EPS
C Test whether step was successful.
      IF(ERRMAX.LT.1.D0) GO TO 52
C Step too large, reduce stepsize and repeat.
      XV=XV-H
      H=0.9D0*H*(ERRMAX**(-0.25D0))
      GO TO 45
 52   CONTINUE
C Step successful, accept solution and improve truncation error.
      DO 53 I=1,6
        DO 54 J=1,3
          Y(I,J)=Y2(I,J)+DELTA(I,J)/15.D0
 54     CONTINUE
 53   CONTINUE
C Prepare to write out free solutions.
      IRAD=IRAD+1
      RAD(IRAD)=XV*1.D-3
C Write out inner core density.
      CALL INTPL(XV,RHOR,N1,C,R,RHO,M1)
      RHOIC(IRAD)=RHOR
C Convert solutions back to y-variables.
      ANM2=AN-2.D0
      DO 55 J=1,3
        IF(J.EQ.3) ANM2=AN
        CONST=XV**ANM2
        FS(2,J,IRAD)=CONST*Y(2,J)
        FS(4,J,IRAD)=CONST*Y(4,J)
        CONST=CONST*XV
        FS(1,J,IRAD)=CONST*Y(1,J)
        FS(3,J,IRAD)=CONST*Y(3,J)
        FS(6,J,IRAD)=CONST*Y(6,J)
        CONST=CONST*XV
        FS(5,J,IRAD)=CONST*Y(5,J)
 55   CONTINUE
      WRITE(3,42)RAD(IRAD),(FS(I,1,IRAD),I=1,6)
      WRITE(4,42)RAD(IRAD),(FS(I,2,IRAD),I=1,6)
      WRITE(7,42)RAD(IRAD),(FS(I,3,IRAD),I=1,6)
C Increase stepsize.
      H=0.9D0*H*(ERRMAX**(-0.20D0))
C Test if solution is finished.
      XVT=XV+H
      IF(XVT.LT.R(N1)) GO TO 45
C Find remaining distance to the ICB.
      H=R(N1)-XV
C Compute fraction of radius.
      RATIO=H/R(N1)
C Continue integration only if remaining distance is more than
C 1.D-9 of radius.
      IF(RATIO.GT.1.D-9) GO TO 45
C Begin integration through the fluid outer core.
C Set up interpolation for the fluid outer core.
      N1=NM(2)
      N2=2*N1-2
      N3=N1-2
C Put fluid outer core values in active locations.
      DO 56 I=1,N1
        J=NK(2)+I
        R(I)=RI(J)
        RHO(I)=RHOI(J)
        MU(I)=MUI(J)
        LAMBDA(I)=LAMBDAI(J)
        GZERO(I)=GZEROI(J)
 56   CONTINUE
C Store density and gravity at bottom of the outer core.
      RHOAP=RHO(1)
      GZEROA=GZERO(1)
C Store density and gravity at top of the outer core.
      RHOBM=RHO(N1)
      GZEROB=GZERO(N1)
C Construct interpolation matrix.
      CALL SPMAT(N1,N2,N3,C,R,B,M1,M2,M3)
C Find starting values of y5 and y6 to within constant factor C.
C Find coefficients cij.
      C11=FS(2,1,IRAD)+RHOAP*FS(5,1,IRAD)-RHOAP*GZEROA*FS(1,1,IRAD)
      C12=FS(2,2,IRAD)+RHOAP*FS(5,2,IRAD)-RHOAP*GZEROA*FS(1,2,IRAD)
      C13=FS(2,3,IRAD)+RHOAP*FS(5,3,IRAD)-RHOAP*GZEROA*FS(1,3,IRAD)
      C21=FS(4,1,IRAD)
      C22=FS(4,2,IRAD)
      C23=FS(4,3,IRAD)
C Find delta (DEL) and epsilon (EPS).
      DEL=(C12*C23-C13*C22)/(C11*C22-C12*C21)
      EPS=(C13*C21-C11*C23)/(C11*C22-C12*C21)
C Find starting values of y5 and y6.
      Y5S=DEL*FS(5,1,IRAD)+EPS*FS(5,2,IRAD)+FS(5,3,IRAD)
      Y6S=DEL*(FS(6,1,IRAD)-4.D0*PI*G*(RHOAP*FS(5,1,IRAD)/GZEROA
     1 -RHOAM*FS(1,1,IRAD)))+EPS*(FS(6,2,IRAD)-4.D0*PI*G*(RHOAP
     2 *FS(5,2,IRAD)/GZEROA-RHOAM*FS(1,2,IRAD)))
     3 +FS(6,3,IRAD)-4.D0*PI*G*(RHOAP*FS(5,3,IRAD)/GZEROA
     4 -RHOAM*FS(1,3,IRAD))
C Initialize propagator matrix for the fluid outer core.
      DO 57 I=1,2
        Y(I,I)=1.D0
        DO 58 J=1,2
          IF(I.NE.J) Y(I,J)=0.D0
 58     CONTINUE
 57   CONTINUE
C Begin Runge-Kutta integration for the fluid outer core.
      M=NI(2)
      AM=DFLOAT(M)
C Find stepsize.
      H=(R(N1)-R(1))/AM
C Set initial radius and density.
      XV=R(1)
      ROC(1)=XV*1.D-3
      RHOOC(1)=RHO(1)
C Set initial value of solution vector.
      YOC(1,1)=Y5S/GZERO(1)
      YOC(2,1)=0.D0
      YOC(3,1)=(2.D0+XV*WES/GZERO(1))*Y5S/(3.D0*GZERO(1))+
     1  XV*Y6S/(6.D0*GZERO(1))
      YOC(4,1)=0.D0
      YOC(5,1)=Y5S
      YOC(6,1)=Y6S
      DO 59 I=1,M
        CALL YPRIMEOC(XV,Y,A,K1,N1,C,R,RHO,GZERO,PI,G,WES)
        CALL RK4OC(XV,Y,A,K1,N1,C,R,RHO,GZERO,PI,G,H,WES)
C Find value GXV of GZERO at radius XV.
        CALL INTPL(XV,GXV,N1,C,R,GZERO,M1)
        IP1=I+1
C Store radius.
        ROC(IP1)=XV*1.D-3
C Write out outer core density.
        CALL INTPL(XV,RHOR,N1,C,R,RHO,M1)
        RHOOC(IP1)=RHOR
C Update solution vector.
        ALPHA=Y(1,1)*Y5S+Y(1,2)*Y6S
        BETA=Y(2,1)*Y5S+Y(2,2)*Y6S
        YOC(1,IP1)=ALPHA/GXV
        YOC(2,IP1)=0.D0
        YOC(3,IP1)=(2.D0+XV*WES/GXV)*ALPHA/(3.D0*GXV)+
     1   XV*BETA/(6.D0*GXV)
        YOC(4,IP1)=0.D0
        YOC(5,IP1)=ALPHA
        YOC(6,IP1)=BETA
 59   CONTINUE
C Input focal depth.
      WRITE(6,60)
 60   FORMAT(1X,'Enter focal depth in kilometres.')
      READ(5,*) FD
C Find radius of focus in kilometres.
      RFD=1.D-3*D-FD
C Find radius at the bottom of the crust.
      J=NK(4)+1
      RC(1)=1.D-3*RI(J)
C Set switch for focus in the mantle.
      IFSW=1
C Set switch for focus in the crust.
      IF(RFD.GT.RC(1)) IFSW=-1
C Begin integration of propagator matrix through the mantle.
C Set up interpolation for the mantle.
      N1=NM(3)
      N2=2*N1-2
      N3=N1-2
C Put mantle values in active locations.
      DO 61 I=1,N1
        J=NK(3)+I
        R(I)=RI(J)
        RHO(I)=RHOI(J)
        MU(I)=MUI(J)
        LAMBDA(I)=LAMBDAI(J)
        GZERO(I)=GZEROI(J)
 61   CONTINUE
C Store density at bottom of the mantle.
      RHOBP=RHO(1)
C Construct interpolation matrix.
      CALL SPMAT(N1,N2,N3,C,R,B,M1,M2,M3)
C Initialize propagator matrix for the mantle and store initial value.
      DO 62 I=1,6
        Y(I,I)=1.D0
        YYM(I,I,1)=1.D0
        DO 63 J=1,6
          IF(I.NE.J) Y(I,J)=0.D0
          IF(I.NE.J) YYM(I,J,1)=0.D0
 63     CONTINUE
 62   CONTINUE
C Begin Runge-Kutta integration for the mantle.
      M=NI(3)
      AM=DFLOAT(M)
C Find stepsize.
      H=(R(N1)-R(1))/AM
C Set initial radius and density.
      XV=R(1)
      RM(1)=XV*1.D-3
      RHOM(1)=RHO(1)
C Set integration region counter.
      IR=3
      DO 64 I=1,M
        CALL YPRIME(XV,Y,A,K1,N1,C,R,RHO,MU,LAMBDA,GZERO,N,PI,G,
     1   IR,RMIN,ANGS,COR,WES)
        CALL RK4(XV,Y,A,K1,N1,C,R,RHO,MU,LAMBDA,GZERO,N,PI,G,H,IR,
     1   RMIN,ANGS,COR,WES)
        IP1=I+1
C Store radius and density.
        RM(IP1)=XV*1.D-3
        CALL INTPL(XV,RHOR,N1,C,R,RHO,M1)
        RHOM(IP1)=RHOR
C Store propagator matrix.
        DO 65 J=1,6
          DO 66 K=1,6 
            YYM(J,K,IP1)=Y(J,K)
 66       CONTINUE
 65     CONTINUE
 64   CONTINUE
C Set exact radius at the top of the mantle.
      RM(IP1)=R(N1)*1.D-3
C For focus in the mantle, calculate coefficient matrix AS(6,6).
      XV=RFD*1.D3
      IF(IFSW.EQ.1) CALL
     1 YPRIME(XV,Y,AS,K1,N1,C,R,RHO,MU,LAMBDA,GZERO,N,PI,G,
     2 IR,RMIN,ANGS,COR,WES)
C If focus is in the mantle, find shear modulus MUFD at focal depth. 
      IF(IFSW.EQ.1) CALL
     1 INTPL(XV,MUFD,N1,C,R,MU,M1)
C If focus is in the mantle, find density at focal depth.
      IF(IFSW.EQ.1) CALL
     1 INTPL(XV,RHOFD,N1,C,R,RHO,M1)
C Begin integration of propagator matrix through the crust.
C Set integration region counter.
      IR=4
C Set number of steps in integration.
      M=NI(4)
      AM=DFLOAT(M)
C Set up interpolation for the crust.
      N1=NM(4)
      N2=2*N1-2
      N3=N1-2
C Put crust values in active locations.
      DO 67 I=1,N1
        J=NK(4)+I
        R(I)=RI(J)
        RHO(I)=RHOI(J)
        MU(I)=MUI(J)
        LAMBDA(I)=LAMBDAI(J)
        GZERO(I)=GZEROI(J)
 67   CONTINUE
C Construct interpolation matrix.
      CALL SPMAT(N1,N2,N3,C,R,B,M1,M2,M3)
C Initialize propagator matrix for the crust.
      DO 68 I=1,6
        DO 69 J=1,6
          Y(I,J)=YYM(I,J,IP1)
          YYC(I,J,1)=Y(I,J)
 69     CONTINUE
 68   CONTINUE
C Find stepsize.
      H=(R(N1)-R(1))/AM
C Set initial radius and density.
      XV=R(1)
      RC(1)=XV*1.D-3
      RHOC(1)=RHO(1)
      DO 70 I=1,M
        CALL YPRIME(XV,Y,A,K1,N1,C,R,RHO,MU,LAMBDA,GZERO,N,PI,G,
     1   IR,RMIN,ANGS,COR,WES)
        CALL RK4(XV,Y,A,K1,N1,C,R,RHO,MU,LAMBDA,GZERO,N,PI,G,H,IR,
     1   RMIN,ANGS,COR,WES)
        IP1=I+1
C Store radius and density.
        RC(IP1)=XV*1.D-3
        CALL INTPL(XV,RHOR,N1,C,R,RHO,M1)
        RHOC(IP1)=RHOR
C Store propagator matrix.
        DO 71 J=1,6
          DO 72 K=1,6
            YYC(J,K,IP1)=Y(J,K)
 72       CONTINUE
 71     CONTINUE
 70   CONTINUE
C If focus is in the crust, calculate coefficient matrix AS(6,6).
      XV=RFD*1.D3
      IF(IFSW.EQ.-1) CALL
     1 YPRIME(XV,Y,AS,K1,N1,C,R,RHO,MU,LAMBDA,GZERO,N,PI,G,
     2 IR,RMIN,ANGS,COR,WES)
C If focus is in the crust, find shear modulus MUFD at focal depth.
      IF(IFSW.EQ.-1) CALL
     1 INTPL(XV,MUFD,N1,C,R,MU,M1)
C If focus is in the crust, find density at focal depth.
      IF(IFSW.EQ.-1) CALL
     1 INTPL(XV,RHOFD,N1,C,R,RHO,M1)
C Determine whether focus is in the mantle or the crust.
      IF(IFSW.EQ.-1) GO TO 73
C Focus is in the mantle.
C Find location of focus in the mantle and select propagator matrix to be
C inverted.
      M=NI(3)
      MP1=M+1
      DO 74 I=2,MP1
        IM1=I-1
        IF(RFD.GT.RM(IM1).AND.RFD.LT.RM(I)) GO TO 75
        IF(RFD.EQ.RM(I)) GO TO 76
        GO TO 74
 75     IEQ=-1
        IL1=I-1
        DO 77 J=1,6
          DO 78 K=1,6
C Take average propagator matrix.
            YINFD(J,K)=(YYM(J,K,I)*(RFD-RM(IM1))+
     1       YYM(J,K,IM1)*(RM(I)-RFD))/(RM(I)-RM(IM1))
 78       CONTINUE
          UV(J)=1.D0
 77     CONTINUE
        GO TO 74
 76     IEQ=1
        IL1=I-1
        DO 79 J=1,6
          DO 80 K=1,6
C Use propagator matrix at this point.
            YINFD(J,K)=YYM(J,K,I)
 80       CONTINUE
          UV(J)=1.D0
 79     CONTINUE
 74   CONTINUE
      GO TO 81
C Focus is in the crust.
C Find location of focus in the crust and select propagator matrix to be
C inverted.
 73   M=NI(4)
      MP1=M+1
      DO 82 I=2,MP1
        IM1=I-1
        IF(RFD.GT.RC(IM1).AND.RFD.LT.RC(I)) GO TO 83
        IF(RFD.EQ.RC(I)) GO TO 84
        GO TO 82
 83     IEQ=-1
        IL1=I-1
        DO 85 J=1,6
          DO 86 K=1,6
C Take average propagator matrix.
            YINFD(J,K)=(YYC(J,K,I)*(RFD-RC(IM1))+
     1       YYC(J,K,IM1)*(RC(I)-RFD))/(RC(I)-RC(IM1))
 86       CONTINUE
          UV(J)=1.D0
 85     CONTINUE
        GO TO 82
 84     IEQ=1
        IL1=I-1
        DO 87 J=1,6
          DO 88 K=1,6
C Use propagator matrix at this point.
            YINFD(J,K)=YYC(J,K,I)
 88       CONTINUE
          UV(J)=1.D0
 87     CONTINUE
 82   CONTINUE
C Store propagator matrix at the focus.
 81   DO 89 I=1,6
        DO 90 J=1,6
          YINFDS(I,J)=YINFD(I,J)
 90     CONTINUE
 89   CONTINUE
C Invert propagator matrix at the focus.
      CALL LINSOL(YINFD,UV,6,CAUG,DET,6,13) 
C Begin solutions for four sources.
C Initialize four source vectors G1,G2,G3,G4.
      XV=RFD*1.D3
C Construct G1.
      FACT=-5.D0*MUFD/(8.D0*PI*XV*XV*XV)
      DO 91 I=1,6
        G1(I)=FACT*YINFD(I,2)
 91   CONTINUE
C Construct G2.
      FACT=FACT*XV
      DO 92 I=1,6
        UV(I)=FACT*AS(I,2)
 92   CONTINUE
      DO 93 I=1,6
        HOLD=0.D0
        DO 94 J=1,6
          HOLD=HOLD+YINFD(I,J)*UV(J)
 94     CONTINUE
        G2(I)=HOLD
 93   CONTINUE
C Construct G3.
      FACT=FACT/XV
      DO 95 I=1,6
        G3(I)=FACT*YINFD(I,4)
 95   CONTINUE
C Construct G4.
      FACT=FACT*XV
      DO 96 I=1,6
        UV(I)=FACT*AS(I,4)
 96   CONTINUE
      DO 97 I=1,6
        HOLD=0.D0
        DO 98 J=1,6
          HOLD=HOLD+YINFD(I,J)*UV(J)
 98     CONTINUE
        G4(I)=HOLD
 97   CONTINUE
C Initialize matrix ABC(6,3) connecting unknowns to the vector y(b).
      ABC(1,1)=1.D0
      ABC(1,2)=0.D0
      ABC(1,3)=0.D0
      ABC(2,1)=RHOBM*GZEROB
      ABC(2,2)=0.D0
      ABC(2,3)=-RHOBM*ALPHA
      ABC(3,1)=0.D0
      ABC(3,2)=1.D0
      ABC(3,3)=0.D0
      ABC(4,1)=0.D0
      ABC(4,2)=0.D0
      ABC(4,3)=0.D0
      ABC(5,1)=0.D0
      ABC(5,2)=0.D0
      ABC(5,3)=ALPHA
      ABC(6,1)=-4.D0*PI*G*RHOBP
      ABC(6,2)=0.D0
      ABC(6,3)=BETA+4.D0*PI*G*RHOBM*ALPHA/GZEROB
C Multiply propagator matrix at Earth's surface into matrix ABC(6,3)
C to get PROD(6,3).
      DO 99 I=1,6
        DO 100 K=1,3
          HOLD=0.D0
          DO 101 J=1,6
            HOLD=HOLD+Y(I,J)*ABC(J,K)
101       CONTINUE
          PROD(I,K)=HOLD
100     CONTINUE
 99   CONTINUE
C Set ABC equal to PROD.
      DO 102 I=1,6
        DO 103 J=1,3
          ABC(I,J)=PROD(I,J)
103     CONTINUE
102   CONTINUE
C Multiply four source vectors by propagator matrix at Earth's surface.
      DO 104 I=1,6
        HOLD1=0.D0
        HOLD2=0.D0
        HOLD3=0.D0
        HOLD4=0.D0
        DO 105 J=1,6
          HOLD1=HOLD1+Y(I,J)*G1(J)
          HOLD2=HOLD2+Y(I,J)*G2(J)
          HOLD3=HOLD3+Y(I,J)*G3(J)
          HOLD4=HOLD4+Y(I,J)*G4(J)
105     CONTINUE
        UV1(I)=HOLD1
        UV2(I)=HOLD2
        UV3(I)=HOLD3
        UV4(I)=HOLD4
104   CONTINUE
C Apply boundary conditions at Earth's surface.
C Construct coefficient matrix AA(3,3).
C Normal stress vanishes.
      DO 106 I=1,3
        AA(1,I)=ABC(2,I)
106   CONTINUE
C Shear stress vanishes.
      DO 107 I=1,3
        AA(2,I)=ABC(4,I)
107   CONTINUE
C Gravity potential is harmonic outside the Earth.
      DO 108 I=1,3
        AA(3,I)=ABC(5,I)+D*ABC(6,I)/3.D0
108   CONTINUE
C Store coefficient matrix AA(3,3) as AAS(3,3).
      DO 109 I=1,3
        DO 110 J=1,3
          AAS(I,J)=AA(I,J)
110     CONTINUE
109   CONTINUE
C Construct constant vectors B1(3), B2(3), B3(3), B4(3) for four sources.
      B1(1)=-UV1(2)
      B1(2)=-UV1(4)
      B1(3)=-UV1(5)-D*UV1(6)/3.D0
C Solve for first source coefficients.
      CALL LINSOL(AA,B1,3,CAUG2,DET,3,7)
      B2(1)=-UV2(2)
      B2(2)=-UV2(4)
      B2(3)=-UV2(5)-D*UV2(6)/3.D0
C Restore matrix AA(3,3).
      DO 111 I=1,3
        DO 112 J=1,3
          AA(I,J)=AAS(I,J)
112     CONTINUE
111   CONTINUE
C Solve for second source coefficients.
      CALL LINSOL(AA,B2,3,CAUG2,DET,3,7)
      B3(1)=-UV3(2)
      B3(2)=-UV3(4)
      B3(3)=-UV3(5)-D*UV3(6)/3.D0
C Restore matrix AA(3,3).
      DO 113 I=1,3
        DO 114 J=1,3
          AA(I,J)=AAS(I,J)
114     CONTINUE
113   CONTINUE
C Solve for third source coefficients.
      CALL LINSOL(AA,B3,3,CAUG2,DET,3,7)
      B4(1)=-UV4(2)
      B4(2)=-UV4(4)
      B4(3)=-UV4(5)-D*UV4(6)/3.D0
C Restore matrix AA(3,3).
      DO 115 I=1,3
        DO 116 J=1,3
          AA(I,J)=AAS(I,J)
116     CONTINUE
115   CONTINUE
C Solve for fourth source coefficients.
      CALL LINSOL(AA,B4,3,CAUG2,DET,3,7)
C Find four source solutions in the inner core.
      DO 117 I=1,IRAD
        DO 118 J=1,6
          UV1(J)=B1(3)*(DEL*FS(J,1,I)+EPS*FS(J,2,I)+FS(J,3,I))
          UV2(J)=B2(3)*(DEL*FS(J,1,I)+EPS*FS(J,2,I)+FS(J,3,I))
          UV3(J)=B3(3)*(DEL*FS(J,1,I)+EPS*FS(J,2,I)+FS(J,3,I))
          UV4(J)=B4(3)*(DEL*FS(J,1,I)+EPS*FS(J,2,I)+FS(J,3,I))
118     CONTINUE
        SIC(1,I)=UV1(1)
        SIC(2,I)=UV1(3)
        SIC(3,I)=UV2(1)
        SIC(4,I)=UV2(3)
        SIC(5,I)=UV3(1)
        SIC(6,I)=UV3(3)
        SIC(7,I)=UV4(1)
        SIC(8,I)=UV4(3)
        WRITE(8,42)RAD(I),(UV1(K),K=1,6)
        WRITE(9,42)RAD(I),(UV2(K),K=1,6)
        WRITE(10,42)RAD(I),(UV3(K),K=1,6)
        WRITE(11,42)RAD(I),(UV4(K),K=1,6)
117   CONTINUE
C Initialize radial displacements at the inner core boundary.
      DISC(1,1)=-UV1(1)
      DISC(2,1)=-UV2(1)
      DISC(3,1)=-UV3(1)
      DISC(4,1)=-UV4(1)
C Find four source solutions in the outer core.
      M=NI(2)
      MP1=M+1
      DO 119 I=1,MP1
        DO 120 J=1,6
          UV1(J)=B1(3)*YOC(J,I)
          UV2(J)=B2(3)*YOC(J,I)
          UV3(J)=B3(3)*YOC(J,I)
          UV4(J)=B4(3)*YOC(J,I)
120     CONTINUE
        SOC(1,I)=UV1(1)
        SOC(2,I)=UV1(3)
        SOC(3,I)=UV2(1)
        SOC(4,I)=UV2(3)
        SOC(5,I)=UV3(1)
        SOC(6,I)=UV3(3)
        SOC(7,I)=UV4(1)
        SOC(8,I)=UV4(3)
C Record discontinuities in radial displacement at the inner core boundary.
        IF(I.EQ.1)DISC(1,1)=DISC(1,1)+UV1(1)
        IF(I.EQ.1)DISC(2,1)=DISC(2,1)+UV2(1)
        IF(I.EQ.1)DISC(3,1)=DISC(3,1)+UV3(1)
        IF(I.EQ.1)DISC(4,1)=DISC(4,1)+UV4(1)
        WRITE(8,42)ROC(I),(UV1(K),K=1,6)
        WRITE(9,42)ROC(I),(UV2(K),K=1,6)
        WRITE(10,42)ROC(I),(UV3(K),K=1,6)
        WRITE(11,42)ROC(I),(UV4(K),K=1,6)
119   CONTINUE
C Initialize displacements at the core-mantle boundary.
      DISC(1,2)=-UV1(1)
      DISC(2,2)=-UV2(1)
      DISC(3,2)=-UV3(1)
      DISC(4,2)=-UV4(1)
C Find four source solutions in the mantle.
C Initialize solution vectors at the base of the mantle.
      UV1(1)=B1(1)
      UV2(1)=B2(1)
      UV3(1)=B3(1)
      UV4(1)=B4(1)
      UV1(2)=-RHOBM*(B1(3)*ALPHA-GZEROB*B1(1))
      UV2(2)=-RHOBM*(B2(3)*ALPHA-GZEROB*B2(1))
      UV3(2)=-RHOBM*(B3(3)*ALPHA-GZEROB*B3(1))
      UV4(2)=-RHOBM*(B4(3)*ALPHA-GZEROB*B4(1))
      UV1(3)=B1(2)
      UV2(3)=B2(2)
      UV3(3)=B3(2)
      UV4(3)=B4(2)
      UV1(4)=0.D0
      UV2(4)=0.D0
      UV3(4)=0.D0
      UV4(4)=0.D0
      UV1(5)=B1(3)*ALPHA
      UV2(5)=B2(3)*ALPHA
      UV3(5)=B3(3)*ALPHA
      UV4(5)=B4(3)*ALPHA
      UV1(6)=B1(3)*BETA+4.D0*PI*G*(RHOBM*B1(3)*ALPHA/GZEROB
     1 -RHOBP*B1(1))
      UV2(6)=B2(3)*BETA+4.D0*PI*G*(RHOBM*B2(3)*ALPHA/GZEROB
     2 -RHOBP*B2(1))
      UV3(6)=B3(3)*BETA+4.D0*PI*G*(RHOBM*B3(3)*ALPHA/GZEROB
     3 -RHOBP*B3(1))
      UV4(6)=B4(3)*BETA+4.D0*PI*G*(RHOBM*B4(3)*ALPHA/GZEROB
     4 -RHOBP*B4(1))
C Solve for four source solutions in the mantle.
      M=NI(3)
      MP1=M+1
C Switch for focus in the mantle.
      IF(IFSW.EQ.1) GO TO 121
      DO 122 I=1,MP1
        DO 123 J=1,6
          HOLD1=0.D0
          HOLD2=0.D0
          HOLD3=0.D0
          HOLD4=0.D0
          DO 124 K=1,6
            HOLD1=HOLD1+YYM(J,K,I)*UV1(K)
            HOLD2=HOLD2+YYM(J,K,I)*UV2(K)
            HOLD3=HOLD3+YYM(J,K,I)*UV3(K)
            HOLD4=HOLD4+YYM(J,K,I)*UV4(K)
124       CONTINUE
          YF1(J)=HOLD1
          YF2(J)=HOLD2
          YF3(J)=HOLD3
          YF4(J)=HOLD4
123     CONTINUE
        SM(1,I)=YF1(1)
        SM(2,I)=YF1(3)
        SM(3,I)=YF2(1)
        SM(4,I)=YF2(3)
        SM(5,I)=YF3(1)
        SM(6,I)=YF3(3)
        SM(7,I)=YF4(1)
        SM(8,I)=YF4(3)
C Record radial displacement discontinuities at the core-mantle boundary.
        IF(I.EQ.1)DISC(1,2)=DISC(1,2)+YF1(1)
        IF(I.EQ.1)DISC(2,2)=DISC(2,2)+YF2(1)
        IF(I.EQ.1)DISC(3,2)=DISC(3,2)+YF3(1)
        IF(I.EQ.1)DISC(4,2)=DISC(4,2)+YF4(1)
        WRITE(8,42)RM(I),(YF1(J),J=1,6)
        WRITE(9,42)RM(I),(YF2(J),J=1,6)
        WRITE(10,42)RM(I),(YF3(J),J=1,6)
        WRITE(11,42)RM(I),(YF4(J),J=1,6)
122   CONTINUE
      GO TO 125
C Focus is in the mantle.
C Output solutions below focal depth.
121   DO 126 I=1,IL1
        DO 127 J=1,6
          HOLD1=0.D0
          HOLD2=0.D0
          HOLD3=0.D0
          HOLD4=0.D0
          DO 128 K=1,6
            HOLD1=HOLD1+YYM(J,K,I)*UV1(K)
            HOLD2=HOLD2+YYM(J,K,I)*UV2(K)
            HOLD3=HOLD3+YYM(J,K,I)*UV3(K)
            HOLD4=HOLD4+YYM(J,K,I)*UV4(K)
128       CONTINUE
          YF1(J)=HOLD1
          YF2(J)=HOLD2
          YF3(J)=HOLD3
          YF4(J)=HOLD4
127     CONTINUE
        SM(1,I)=YF1(1)
        SM(2,I)=YF1(3)
        SM(3,I)=YF2(1)
        SM(4,I)=YF2(3)
        SM(5,I)=YF3(1)
        SM(6,I)=YF3(3)
        SM(7,I)=YF4(1)
        SM(8,I)=YF4(3)
C Record radial displacement discontinuities at the core-mantle boundary.
        IF(I.EQ.1)DISC(1,2)=DISC(1,2)+YF1(1)
        IF(I.EQ.1)DISC(2,2)=DISC(2,2)+YF2(1)
        IF(I.EQ.1)DISC(3,2)=DISC(3,2)+YF3(1)
        IF(I.EQ.1)DISC(4,2)=DISC(4,2)+YF4(1)
        WRITE(8,42)RM(I),(YF1(J),J=1,6)
        WRITE(9,42)RM(I),(YF2(J),J=1,6)
        WRITE(10,42)RM(I),(YF3(J),J=1,6)
        WRITE(11,42)RM(I),(YF4(J),J=1,6)
126   CONTINUE
C Augment vectors UV1, UV2, UV3, UV4.
      DO 129 I=1,6
        UV1(I)=UV1(I)+G1(I)
        UV2(I)=UV2(I)+G2(I)
        UV3(I)=UV3(I)+G3(I)
        UV4(I)=UV4(I)+G4(I)
129   CONTINUE
C  Output solutions at and above focal depth.
      IF(IEQ.EQ.1) GO TO 130
C Find solutions at focus.
      DO 131 I=1,6
        HOLD1=0.D0
        HOLD2=0.D0
        HOLD3=0.D0
        HOLD4=0.D0
        DO 133 J=1,6
          HOLD1=HOLD1+YINFDS(I,J)*UV1(J)
          HOLD2=HOLD2+YINFDS(I,J)*UV2(J)
          HOLD3=HOLD3+YINFDS(I,J)*UV3(J)
          HOLD4=HOLD4+YINFDS(I,J)*UV4(J)
133     CONTINUE
        YF1(I)=HOLD1
        YF2(I)=HOLD2
        YF3(I)=HOLD3
        YF4(I)=HOLD4
131   CONTINUE
      IST=IL1+1
      SM(1,IST)=YF1(1)
      SM(2,IST)=YF1(3)
      SM(3,IST)=YF2(1)
      SM(4,IST)=YF2(3)
      SM(5,IST)=YF3(1)
      SM(6,IST)=YF3(3)
      SM(7,IST)=YF4(1)
      SM(8,IST)=YF4(3)
      WRITE(8,42)RFD,(YF1(I),I=1,6)
      WRITE(9,42)RFD,(YF2(I),I=1,6)
      WRITE(10,42)RFD,(YF3(I),I=1,6)
      WRITE(11,42)RFD,(YF4(I),I=1,6)
C Modify vector RM to include RFD, modify vector RHOM to include RHOFD.
      DO 132 I=IST,MP1
        J=MP1+IST-I+1
        RM(J)=RM(J-1)
        RHOM(J)=RHOM(J-1)
132   CONTINUE
      RM(IST)=RFD
      RHOM(IST)=RHOFD
130   IST=IL1+1
      DO 134 I=IST, MP1
        IP=I
        IF(IEQ.EQ.-1)IP=I+1
        DO 135 J=1,6
          HOLD1=0.D0
          HOLD2=0.D0
          HOLD3=0.D0
          HOLD4=0.D0
          DO 136 K=1,6
            HOLD1=HOLD1+YYM(J,K,I)*UV1(K)
            HOLD2=HOLD2+YYM(J,K,I)*UV2(K)
            HOLD3=HOLD3+YYM(J,K,I)*UV3(K)
            HOLD4=HOLD4+YYM(J,K,I)*UV4(K)
136       CONTINUE
          YF1(J)=HOLD1
          YF2(J)=HOLD2
          YF3(J)=HOLD3
          YF4(J)=HOLD4
135     CONTINUE
        SM(1,IP)=YF1(1)
        SM(2,IP)=YF1(3)
        SM(3,IP)=YF2(1)
        SM(4,IP)=YF2(3)
        SM(5,IP)=YF3(1)
        SM(6,IP)=YF3(3)
        SM(7,IP)=YF4(1)
        SM(8,IP)=YF4(3)
        WRITE(8,42)RM(IP),(YF1(J),J=1,6)
        WRITE(9,42)RM(IP),(YF2(J),J=1,6)
        WRITE(10,42)RM(IP),(YF3(J),J=1,6)
        WRITE(11,42)RM(IP),(YF4(J),J=1,6)
134   CONTINUE
125   CONTINUE
C Find four source solutions in the crust.
      M=NI(4)
      MP1=M+1
C Switch for focus in the crust.
      IF(IFSW.EQ.-1) GO TO 137
      DO 138 I=1,MP1
        DO 139 J=1,6
          HOLD1=0.D0
          HOLD2=0.D0
          HOLD3=0.D0
          HOLD4=0.D0
          DO 140 K=1,6
            HOLD1=HOLD1+YYC(J,K,I)*UV1(K)
            HOLD2=HOLD2+YYC(J,K,I)*UV2(K)
            HOLD3=HOLD3+YYC(J,K,I)*UV3(K)
            HOLD4=HOLD4+YYC(J,K,I)*UV4(K)
140       CONTINUE
          YF1(J)=HOLD1
          YF2(J)=HOLD2
          YF3(J)=HOLD3
          YF4(J)=HOLD4
139     CONTINUE
        SC(1,I)=YF1(1)
        SC(2,I)=YF1(3)
        SC(3,I)=YF2(1)
        SC(4,I)=YF2(3)
        SC(5,I)=YF3(1)
        SC(6,I)=YF3(3)
        SC(7,I)=YF4(1)
        SC(8,I)=YF4(3)
        WRITE(8,42)RC(I),(YF1(J),J=1,6)
        WRITE(9,42)RC(I),(YF2(J),J=1,6)
        WRITE(10,42)RC(I),(YF3(J),J=1,6)
        WRITE(11,42)RC(I),(YF4(J),J=1,6)
138   CONTINUE
      GO TO 141
C Focus is in the crust.
C Output solutions below focal depth.
137   DO 142 I=1,IL1
        DO 143 J=1,6
          HOLD1=0.D0
          HOLD2=0.D0
          HOLD3=0.D0
          HOLD4=0.D0
          DO 144 K=1,6
            HOLD1=HOLD1+YYC(J,K,I)*UV1(K)
            HOLD2=HOLD2+YYC(J,K,I)*UV2(K)
            HOLD3=HOLD3+YYC(J,K,I)*UV3(K)
            HOLD4=HOLD4+YYC(J,K,I)*UV4(K)
144       CONTINUE
          YF1(J)=HOLD1
          YF2(J)=HOLD2
          YF3(J)=HOLD3
          YF4(J)=HOLD4
143     CONTINUE
        SC(1,I)=YF1(1)
        SC(2,I)=YF1(3)
        SC(3,I)=YF2(1)
        SC(4,I)=YF2(3)
        SC(5,I)=YF3(1)
        SC(6,I)=YF3(3)
        SC(7,I)=YF4(1)
        SC(8,I)=YF4(3)
        WRITE(8,42)RC(I),(YF1(J),J=1,6)
        WRITE(9,42)RC(I),(YF2(J),J=1,6)
        WRITE(10,42)RC(I),(YF3(J),J=1,6)
        WRITE(11,42)RC(I),(YF4(J),J=1,6)
142   CONTINUE
C Augment vectors UV1, UV2, UV3, UV4.
      DO 145 I=1,6
        UV1(I)=UV1(I)+G1(I)
        UV2(I)=UV2(I)+G2(I)
        UV3(I)=UV3(I)+G3(I)
        UV4(I)=UV4(I)+G4(I)
145   CONTINUE
C Output solutions at and above focal depth.
      IF(IEQ.EQ.1) GO TO 146
C Find solutions at the focus.
      DO 147 I=1,6
        HOLD1=0.D0
        HOLD2=0.D0
        HOLD3=0.D0
        HOLD4=0.D0
        DO 149 J=1,6
          HOLD1=HOLD1+YINFDS(I,J)*UV1(J)
          HOLD2=HOLD2+YINFDS(I,J)*UV2(J)
          HOLD3=HOLD3+YINFDS(I,J)*UV3(J)
          HOLD4=HOLD4+YINFDS(I,J)*UV4(J)
149     CONTINUE
        YF1(I)=HOLD1
        YF2(I)=HOLD2
        YF3(I)=HOLD3
        YF4(I)=HOLD4
147   CONTINUE
      IST=IL1+1
      SC(1,IST)=YF1(1)
      SC(2,IST)=YF1(3)
      SC(3,IST)=YF2(1)
      SC(4,IST)=YF2(3)
      SC(5,IST)=YF3(1)
      SC(6,IST)=YF3(3)
      SC(7,IST)=YF4(1)
      SC(8,IST)=YF4(3)
      WRITE(8,42)RFD,(YF1(I),I=1,6)
      WRITE(9,42)RFD,(YF2(I),I=1,6)
      WRITE(10,42)RFD,(YF3(I),I=1,6)
      WRITE(11,42)RFD,(YF4(I),I=1,6)
C Modify vector RC to include RFD, modify vector RHOC to include RHOFD.
      DO 148 I=IST,MP1
        J=MP1+IST-I+1
        RC(J)=RC(J-1)
        RHOC(J)=RHOC(J-1)
148   CONTINUE
      RC(IST)=RFD
      RHOC(IST)=RHOFD
146   IST=IL1+1
      DO 150 I=IST,MP1
        IP=I
        IF(IEQ.EQ.-1)IP=I+1
        DO 151 J=1,6
          HOLD1=0.D0
          HOLD2=0.D0
          HOLD3=0.D0
          HOLD4=0.D0
          DO 152 K=1,6
            HOLD1=HOLD1+YYC(J,K,I)*UV1(K)
            HOLD2=HOLD2+YYC(J,K,I)*UV2(K)
            HOLD3=HOLD3+YYC(J,K,I)*UV3(K)
            HOLD4=HOLD4+YYC(J,K,I)*UV4(K)
152       CONTINUE
          YF1(J)=HOLD1
          YF2(J)=HOLD2
          YF3(J)=HOLD3
          YF4(J)=HOLD4
151     CONTINUE
        SC(1,IP)=YF1(1)
        SC(2,IP)=YF1(3)
        SC(3,IP)=YF2(1)
        SC(4,IP)=YF2(3)
        SC(5,IP)=YF3(1)
        SC(6,IP)=YF3(3)
        SC(7,IP)=YF4(1)
        SC(8,IP)=YF4(3)
        WRITE(8,42)RC(IP),(YF1(J),J=1,6)
        WRITE(9,42)RC(IP),(YF2(J),J=1,6)
        WRITE(10,42)RC(IP),(YF3(J),J=1,6)
        WRITE(11,42)RC(IP),(YF4(J),J=1,6)
150   CONTINUE
141   CONTINUE
C Begin integrations of the radial and transverse components of the
C four fundamental solutions over radius.
C Initialize integrals.
      S1IU=0.D0
      S1IV=0.D0
      S2IU=0.D0
      S2IV=0.D0
      S3IU=0.D0
      S3IV=0.D0
      S4IU=0.D0
      S4IV=0.D0
C Begin integration through the inner core using the trapezoidal rule.
C Set number of segments.
      IRADM1=IRAD-1
      DO 153 I=1,IRADM1
        IP1=I+1
C Change radii to metres.
        RAD(I)=RAD(I)*1.D3
        RAD(IP1)=RAD(IP1)*1.D3
C Find cubes of radii.
        RC1=RAD(I)*RAD(I)*RAD(I)
        RC2=RAD(IP1)*RAD(IP1)*RAD(IP1)
C Find increment in radius for segment.
        DR=RAD(IP1)-RAD(I)
C Find densities at segment boundaries.
        RHO1=RHOIC(I)
        RHO2=RHOIC(IP1)
C Add contributions of segment to integrals.
        S1IU=S1IU+DR*(RC1*RHO1*SIC(1,I)+RC2*RHO2*SIC(1,IP1))/2.D0
        S1IV=S1IV+DR*(RC1*RHO1*SIC(2,I)+RC2*RHO2*SIC(2,IP1))/2.D0
        S2IU=S2IU+DR*(RC1*RHO1*SIC(3,I)+RC2*RHO2*SIC(3,IP1))/2.D0
        S2IV=S2IV+DR*(RC1*RHO1*SIC(4,I)+RC2*RHO2*SIC(4,IP1))/2.D0
        S3IU=S3IU+DR*(RC1*RHO1*SIC(5,I)+RC2*RHO2*SIC(5,IP1))/2.D0
        S3IV=S3IV+DR*(RC1*RHO1*SIC(6,I)+RC2*RHO2*SIC(6,IP1))/2.D0
        S4IU=S4IU+DR*(RC1*RHO1*SIC(7,I)+RC2*RHO2*SIC(7,IP1))/2.D0
        S4IV=S4IV+DR*(RC1*RHO1*SIC(8,I)+RC2*RHO2*SIC(8,IP1))/2.D0
153   CONTINUE
C Store results for inner core.
      RES(1,1)=S1IU
      RES(2,1)=S1IV
      RES(3,1)=S2IU
      RES(4,1)=S2IV
      RES(5,1)=S3IU
      RES(6,1)=S3IV
      RES(7,1)=S4IU
      RES(8,1)=S4IV
C Begin integration through the outer core using extended Simpson's rule.
C Find number of model points.
      N1=NM(2)
C Find stepsize.
      M=NI(2)
      AM=DFLOAT(M)
      H=(R(N1)-R(1))/AM
      MP1=M+1
C Establish vector of functions, VOC, at each segment border, and vector,
C WOC, of weights.
      DO 154 I=1,MP1
C Convert radius to metres.
        CUR=ROC(I)*1.D3
C Cube radius.
        CUR=CUR*CUR*CUR
C Construct VOC.
        VOC(1,I)=H*(CUR*RHOOC(I)*SOC(1,I))/3.D0
        VOC(2,I)=H*(CUR*RHOOC(I)*SOC(2,I))/3.D0
        VOC(3,I)=H*(CUR*RHOOC(I)*SOC(3,I))/3.D0
        VOC(4,I)=H*(CUR*RHOOC(I)*SOC(4,I))/3.D0
        VOC(5,I)=H*(CUR*RHOOC(I)*SOC(5,I))/3.D0
        VOC(6,I)=H*(CUR*RHOOC(I)*SOC(6,I))/3.D0
        VOC(7,I)=H*(CUR*RHOOC(I)*SOC(7,I))/3.D0
        VOC(8,I)=H*(CUR*RHOOC(I)*SOC(8,I))/3.D0
154   CONTINUE
C Construct WOC.
      WOC(1)=1.D0
      WOC(2)=4.D0
      WOC(3)=2.D0
      WOC(4)=4.D0
      WOC(5)=2.D0
      WOC(6)=4.D0
      WOC(7)=2.D0
      WOC(8)=4.D0
      WOC(9)=2.D0
      WOC(10)=4.D0
      WOC(11)=1.D0
C Find integrals as the scalar products of the vectors WOC and VOC.
      HOLD1=0.D0
      HOLD2=0.D0
      HOLD3=0.D0
      HOLD4=0.D0
      HOLD5=0.D0
      HOLD6=0.D0
      HOLD7=0.D0
      HOLD8=0.D0
      DO 155 I=1,MP1
        HOLD1=HOLD1+WOC(I)*VOC(1,I)
        HOLD2=HOLD2+WOC(I)*VOC(2,I)
        HOLD3=HOLD3+WOC(I)*VOC(3,I)
        HOLD4=HOLD4+WOC(I)*VOC(4,I)
        HOLD5=HOLD5+WOC(I)*VOC(5,I)
        HOLD6=HOLD6+WOC(I)*VOC(6,I)
        HOLD7=HOLD7+WOC(I)*VOC(7,I)
        HOLD8=HOLD8+WOC(I)*VOC(8,I)
155   CONTINUE
C Update integrals.
      S1IU=S1IU+HOLD1
      S1IV=S1IV+HOLD2
      S2IU=S2IU+HOLD3
      S2IV=S2IV+HOLD4
      S3IU=S3IU+HOLD5
      S3IV=S3IV+HOLD6
      S4IU=S4IU+HOLD7
      S4IV=S4IV+HOLD8
C Store results for outer core.
      RES(1,2)=S1IU-RES(1,1)
      RES(2,2)=S1IV-RES(2,1)
      RES(3,2)=S2IU-RES(3,1)
      RES(4,2)=S2IV-RES(4,1)
      RES(5,2)=S3IU-RES(5,1)
      RES(6,2)=S3IV-RES(6,1)
      RES(7,2)=S4IU-RES(7,1)
      RES(8,2)=S4IV-RES(8,1)
C Begin integration through the mantle using the trapezoidal rule.
C Find number of segments M.
      M=NI(3)
C Switch for focus in the mantle.
      IF(IFSW.EQ.1) GO TO 156
      DO 157 I=1,M
        IP1=I+1
C Change radii to metres.
        RM(I)=RM(I)*1.D3
        RM(IP1)=RM(IP1)*1.D3
C Find cubes of radii.
        RCM1=RM(I)*RM(I)*RM(I)
        RCM2=RM(IP1)*RM(IP1)*RM(IP1)
C Find increment in radius for segment.
        DR=RM(IP1)-RM(I)
C Find densities at segment boundaries.
        RHO1=RHOM(I)
        RHO2=RHOM(IP1)
C Add contributions of segment to integrals.
        S1IU=S1IU+DR*(RCM1*RHO1*SM(1,I)+RCM2*RHO2*SM(1,IP1))/2.D0
        S1IV=S1IV+DR*(RCM1*RHO1*SM(2,I)+RCM2*RHO2*SM(2,IP1))/2.D0
        S2IU=S2IU+DR*(RCM1*RHO1*SM(3,I)+RCM2*RHO2*SM(3,IP1))/2.D0
        S2IV=S2IV+DR*(RCM1*RHO1*SM(4,I)+RCM2*RHO2*SM(4,IP1))/2.D0
        S3IU=S3IU+DR*(RCM1*RHO1*SM(5,I)+RCM2*RHO2*SM(5,IP1))/2.D0
        S3IV=S3IV+DR*(RCM1*RHO1*SM(6,I)+RCM2*RHO2*SM(6,IP1))/2.D0
        S4IU=S4IU+DR*(RCM1*RHO1*SM(7,I)+RCM2*RHO2*SM(7,IP1))/2.D0
        S4IV=S4IV+DR*(RCM1*RHO1*SM(8,I)+RCM2*RHO2*SM(8,IP1))/2.D0
157   CONTINUE
      GO TO 158
C Focus is in the mantle.
C Add contributions of segments to two below focal depth to integrals.
156   IL2=IL1-1
      DO 159 I=1,IL2
        IP1=I+1
C Change radii to metres.
        RM(I)=RM(I)*1.D3
        RM(IP1)=RM(IP1)*1.D3
C Find cubes of radii.
        RCM1=RM(I)*RM(I)*RM(I)
        RCM2=RM(IP1)*RM(IP1)*RM(IP1)
C Find increment in radius for segment.
        DR=RM(IP1)-RM(I)
C Find densities at segment boundaries.
        RHO1=RHOM(I)
        RHO2=RHOM(IP1)
C Add contributions of segment to integrals.
        S1IU=S1IU+DR*(RCM1*RHO1*SM(1,I)+RCM2*RHO2*SM(1,IP1))/2.D0
        S1IV=S1IV+DR*(RCM1*RHO1*SM(2,I)+RCM2*RHO2*SM(2,IP1))/2.D0
        S2IU=S2IU+DR*(RCM1*RHO1*SM(3,I)+RCM2*RHO2*SM(3,IP1))/2.D0
        S2IV=S2IV+DR*(RCM1*RHO1*SM(4,I)+RCM2*RHO2*SM(4,IP1))/2.D0
        S3IU=S3IU+DR*(RCM1*RHO1*SM(5,I)+RCM2*RHO2*SM(5,IP1))/2.D0
        S3IV=S3IV+DR*(RCM1*RHO1*SM(6,I)+RCM2*RHO2*SM(6,IP1))/2.D0
        S4IU=S4IU+DR*(RCM1*RHO1*SM(7,I)+RCM2*RHO2*SM(7,IP1))/2.D0
        S4IV=S4IV+DR*(RCM1*RHO1*SM(8,I)+RCM2*RHO2*SM(8,IP1))/2.D0
159   CONTINUE
C Add contributions at and above focal depth.
C Do linear extrapolation of sources to focus.
      IST=IL1+1
      DO 161 I=1,8
        SX(I)=SM(I,IL2)+(SM(I,IL1)-SM(I,IL2))*(RM(IST)-RM(IL2))/
     1   (RM(IL1)-RM(IL2))
161   CONTINUE
C Add contribution of segment before focus.
C Change radii to metres.
      RM(IL1)=RM(IL1)*1.D3
      RM(IST)=RM(IST)*1.D3
C Find cubes of radii.
      RCM1=RM(IL1)*RM(IL1)*RM(IL1)
      RCM2=RM(IST)*RM(IST)*RM(IST)
C Find increment in radius for segment.
      DR=RM(IST)-RM(IL1)
C Find densities at segment boundaries.
      RHO1=RHOM(IL1)
      RHO2=RHOM(IST)
C Update contributions.
      S1IU=S1IU+DR*(RCM1*RHO1*SM(1,IL1)+RCM2*RHO2*SX(1))/2.D0
      S1IV=S1IV+DR*(RCM1*RHO1*SM(2,IL1)+RCM2*RHO2*SX(2))/2.D0
      S2IU=S2IU+DR*(RCM1*RHO1*SM(3,IL1)+RCM2*RHO2*SX(3))/2.D0
      S2IV=S2IV+DR*(RCM1*RHO1*SM(4,IL1)+RCM2*RHO2*SX(4))/2.D0
      S3IU=S3IU+DR*(RCM1*RHO1*SM(5,IL1)+RCM2*RHO2*SX(5))/2.D0
      S3IV=S3IV+DR*(RCM1*RHO1*SM(6,IL1)+RCM2*RHO2*SX(6))/2.D0
      S4IU=S4IV+DR*(RCM1*RHO1*SM(7,IL1)+RCM2*RHO2*SX(7))/2.D0
      S4IV=S4IV+DR*(RCM1*RHO1*SM(8,IL1)+RCM2*RHO2*SX(8))/2.D0
C Add contributions at focus and above.
      MP1=M+1
      IFIN=MP1
      IF(IEQ.EQ.1) IFIN=M
      DO 160 I=IST,IFIN
        IP1=I+1
C Change radii to metres.
        RM(I)=RM(I)*1.D3
        RM(IP1)=RM(IP1)*1.D3
C Find cubes of radii.
        RCM1=RM(I)*RM(I)*RM(I)
        RCM2=RM(IP1)*RM(IP1)*RM(IP1)
C Find increment in radius for segment.
        DR=RM(IP1)-RM(I)
C Find densities at segment boundaries.
        RHO1=RHOM(I)
        RHO2=RHOM(IP1)
C Add contributions of segment to integrals.
        S1IU=S1IU+DR*(RCM1*RHO1*SM(1,I)+RCM2*RHO2*SM(1,IP1))/2.D0
        S1IV=S1IV+DR*(RCM1*RHO1*SM(2,I)+RCM2*RHO2*SM(2,IP1))/2.D0
        S2IU=S2IU+DR*(RCM1*RHO1*SM(3,I)+RCM2*RHO2*SM(3,IP1))/2.D0
        S2IV=S2IV+DR*(RCM1*RHO1*SM(4,I)+RCM2*RHO2*SM(4,IP1))/2.D0
        S3IU=S3IU+DR*(RCM1*RHO1*SM(5,I)+RCM2*RHO2*SM(5,IP1))/2.D0
        S3IV=S3IV+DR*(RCM1*RHO1*SM(6,I)+RCM2*RHO2*SM(6,IP1))/2.D0
        S4IU=S4IU+DR*(RCM1*RHO1*SM(7,I)+RCM2*RHO2*SM(7,IP1))/2.D0
        S4IV=S4IV+DR*(RCM1*RHO1*SM(8,I)+RCM2*RHO2*SM(8,IP1))/2.D0
160   CONTINUE
158   CONTINUE
C Store results for mantle.
      RES(1,3)=S1IU-RES(1,2)-RES(1,1)
      RES(2,3)=S1IV-RES(2,2)-RES(2,1)
      RES(3,3)=S2IU-RES(3,2)-RES(3,1)
      RES(4,3)=S2IV-RES(4,2)-RES(4,1)
      RES(5,3)=S3IU-RES(5,2)-RES(5,1)
      RES(6,3)=S3IV-RES(6,2)-RES(6,1)
      RES(7,3)=S4IU-RES(7,2)-RES(7,1)
      RES(8,3)=S4IV-RES(8,2)-RES(8,1)
C Begin integration through the crust using the trapezoidal rule.
C Find number of segments M.
      M=NI(4)
C Switch for focus in the crust.
      IF(IFSW.EQ.-1) GO TO 162
      DO 163 I=1,M
        IP1=I+1
C Change radii to metres.
        RC(I)=RC(I)*1.D3
        RC(IP1)=RC(IP1)*1.D3
C Find cubes of radii.
        RCC1=RC(I)*RC(I)*RC(I)
        RCC2=RC(IP1)*RC(IP1)*RC(IP1)
C Find increment in radius for segment.
        DR=RC(IP1)-RC(I)
C Find densities at segment boundaries.
        RHO1=RHOC(I)
        RHO2=RHOC(IP1)
C Add contributions of segment to integrals.
        S1IU=S1IU+DR*(RCC1*RHO1*SC(1,I)+RCC2*RHO2*SC(1,IP1))/2.D0
        S1IV=S1IV+DR*(RCC1*RHO1*SC(2,I)+RCC2*RHO2*SC(2,IP1))/2.D0
        S2IU=S2IU+DR*(RCC1*RHO1*SC(3,I)+RCC2*RHO2*SC(3,IP1))/2.D0
        S2IV=S2IV+DR*(RCC1*RHO1*SC(4,I)+RCC2*RHO2*SC(4,IP1))/2.D0
        S3IU=S3IU+DR*(RCC1*RHO1*SC(5,I)+RCC2*RHO2*SC(5,IP1))/2.D0
        S3IV=S3IV+DR*(RCC1*RHO1*SC(6,I)+RCC2*RHO2*SC(6,IP1))/2.D0
        S4IU=S4IU+DR*(RCC1*RHO1*SC(7,I)+RCC2*RHO2*SC(7,IP1))/2.D0
        S4IV=S4IV+DR*(RCC1*RHO1*SC(8,I)+RCC2*RHO2*SC(8,IP1))/2.D0
163   CONTINUE
      GO TO 164
C Focus is in the crust.
C Add contributions of segments to two below focal depth to integrals.
162   IL2=IL1-1
      DO 165 I=1,IL2
        IP1=I+1
C Change radii to metres.
        RC(I)=RC(I)*1.D3
        RC(IP1)=RC(IP1)*1.D3
C Find cubes of radii.
        RCC1=RC(I)*RC(I)*RC(I)
        RCC2=RC(IP1)*RC(IP1)*RC(IP1)
C Find increment in radius for segment.
        DR=RC(IP1)-RC(I)
C Find densities at segment boundaries.
        RHO1=RHOC(I)
        RHO2=RHOC(IP1)
C Add contributions of segment to integrals.
        S1IU=S1IU+DR*(RCC1*RHO1*SC(1,I)+RCC2*RHO2*SC(1,IP1))/2.D0
        S1IV=S1IV+DR*(RCC1*RHO1*SC(2,I)+RCC2*RHO2*SC(2,IP1))/2.D0
        S2IU=S2IU+DR*(RCC1*RHO1*SC(3,I)+RCC2*RHO2*SC(3,IP1))/2.D0
        S2IV=S2IV+DR*(RCC1*RHO1*SC(4,I)+RCC2*RHO2*SC(4,IP1))/2.D0
        S3IU=S3IU+DR*(RCC1*RHO1*SC(5,I)+RCC2*RHO2*SC(5,IP1))/2.D0
        S3IV=S3IV+DR*(RCC1*RHO1*SC(6,I)+RCC2*RHO2*SC(6,IP1))/2.D0
        S4IU=S4IU+DR*(RCC1*RHO1*SC(7,I)+RCC2*RHO2*SC(7,IP1))/2.D0
        S4IV=S4IV+DR*(RCC1*RHO1*SC(8,I)+RCC2*RHO2*SC(8,IP1))/2.D0
165   CONTINUE
C Add contributions at and above focal depth.
C Do linear extrapolation of sources to focus.
      IST=IL1+1
      DO 166 I=1,8
        SX(I)=SC(I,IL2)+(SC(I,IL1)-SC(I,IL2))*(RC(IST)-RC(IL2))/
     1   (RC(IL1)-RC(IL2))
166   CONTINUE
C Add contribution of segment before focus.
C Change radii to metres.
      RC(IL1)=RC(IL1)*1.D3
      RC(IST)=RC(IST)*1.D3
C Find cubes of radii.
      RCC1=RC(IL1)*RC(IL1)*RC(IL1)
      RCC2=RC(IST)*RC(IST)*RC(IST)
C Find increment in radius for segment.
      DR=RC(IST)-RC(IL1)
C Find densities at segment boundaries.
      RHO1=RHOC(IL1)
      RHO1=RHOC(IST)
C Update contributions.
      S1IU=S1IU+DR*(RCC1*RHO1*SC(1,IL1)+RCC2*RHO2*SX(1))/2.D0
      S1IV=S1IV+DR*(RCC1*RHO1*SC(2,IL1)+RCC2*RHO2*SX(2))/2.D0
      S2IU=S2IU+DR*(RCC1*RHO1*SC(3,IL1)+RCC2*RHO2*SX(3))/2.D0
      S2IV=S2IV+DR*(RCC1*RHO1*SC(4,IL1)+RCC2*RHO2*SX(4))/2.D0
      S3IU=S3IU+DR*(RCC1*RHO1*SC(5,IL1)+RCC2*RHO2*SX(5))/2.D0
      S3IV=S3IV+DR*(RCC1*RHO1*SC(6,IL1)+RCC2*RHO2*SX(6))/2.D0
      S4IU=S4IU+DR*(RCC1*RHO1*SC(7,IL1)+RCC2*RHO2*SX(7))/2.D0
      S4IV=S4IV+DR*(RCC1*RHO1*SC(8,IL1)+RCC2*RHO2*SX(8))/2.D0
C Add contributions at focus and above.
      MP1=M+1
      IFIN=MP1
      IF(IEQ.EQ.1) IFIN=M
      DO 167 I=IST,IFIN
        IP1=I+1
C Change radii to metres.
        RC(I)=RC(I)*1.D3
        RC(IP1)=RC(IP1)*1.D3
C Find cubes of radii.
        RCC1=RC(I)*RC(I)*RC(I)
        RCC2=RC(IP1)*RC(IP1)*RC(IP1)
C Find increment in radius for segment.
        DR=RC(IP1)-RC(I)
C Find densities at segment boundaries.
        RHO1=RHOC(I)
        RHO2=RHOC(IP1)
C Add contributions of segment to integrals.
        S1IU=S1IU+DR*(RCC1*RHO1*SC(1,I)+RCC2*RHO2*SC(1,IP1))/2.D0
        S1IV=S1IV+DR*(RCC1*RHO1*SC(2,I)+RCC2*RHO2*SC(2,IP1))/2.D0
        S2IU=S2IU+DR*(RCC1*RHO1*SC(3,I)+RCC2*RHO2*SC(3,IP1))/2.D0
        S2IV=S2IV+DR*(RCC1*RHO1*SC(4,I)+RCC2*RHO2*SC(4,IP1))/2.D0
        S3IU=S3IU+DR*(RCC1*RHO1*SC(5,I)+RCC2*RHO2*SC(5,IP1))/2.D0
        S3IV=S3IV+DR*(RCC1*RHO1*SC(6,I)+RCC2*RHO2*SC(6,IP1))/2.D0
        S4IU=S4IU+DR*(RCC1*RHO1*SC(7,I)+RCC2*RHO2*SC(7,IP1))/2.D0
        S4IV=S4IV+DR*(RCC1*RHO1*SC(8,I)+RCC2*RHO2*SC(8,IP1))/2.D0
167   CONTINUE
164   CONTINUE
C Store results for crust.
      RES(1,4)=S1IU-RES(1,3)-RES(1,2)-RES(1,1)
      RES(2,4)=S1IV-RES(2,3)-RES(2,2)-RES(2,1)
      RES(3,4)=S2IU-RES(3,3)-RES(3,2)-RES(3,1)
      RES(4,4)=S2IV-RES(4,3)-RES(4,2)-RES(4,1)
      RES(5,4)=S3IU-RES(5,3)-RES(5,2)-RES(5,1)
      RES(6,4)=S3IV-RES(6,3)-RES(6,2)-RES(6,1)
      RES(7,4)=S4IU-RES(7,3)-RES(7,2)-RES(7,1)
      RES(8,4)=S4IV-RES(8,3)-RES(8,2)-RES(8,1)
C Store results for whole Earth.
      RES(1,5)=S1IU
      RES(2,5)=S1IV
      RES(3,5)=S2IU
      RES(4,5)=S2IV
      RES(5,5)=S3IU
      RES(6,5)=S3IV
      RES(7,5)=S4IU
      RES(8,5)=S4IV
C Write out results to screen.
      WRITE(6,168) FD
168   FORMAT(/5X,'Results for focal depth =',1X,F5.1,1X,'km.')
      WRITE(6,169)
169   FORMAT(/11X,'INTEGRALS OVER RADIUS')
      WRITE(6,170)
170   FORMAT(/1X,'source',5X,'integral over u',5X,'integral over v',5X,
     1 'region')
      WRITE(6,171)RES(1,1),RES(2,1)
171   FORMAT(3X,'1'/14X,D11.4,9X,D11.4,7X,'inner core')
      WRITE(6,172)RES(1,2),RES(2,2)
172   FORMAT(14X,D11.4,9X,D11.4,7X,'outer core')
      WRITE(6,173)RES(1,3),RES(2,3)
173   FORMAT(14X,D11.4,9X,D11.4,7X,'mantle')
      WRITE(6,174)RES(1,4),RES(2,4)
174   FORMAT(14X,D11.4,9X,D11.4,7X,'crust')
      WRITE(6,175)RES(1,5),RES(2,5)
175   FORMAT(14X,D11.4,9X,D11.4,7X,'whole Earth')
      WRITE(6,176)RES(3,1),RES(4,1)
176   FORMAT(3X,'2'/14X,D11.4,9X,D11.4,7X,'inner core')
      WRITE(6,177)RES(3,2),RES(4,2)
177   FORMAT(14X,D11.4,9X,D11.4,7X,'outer core')
      WRITE(6,178)RES(3,3),RES(4,3)
178   FORMAT(14X,D11.4,9X,D11.4,7X,'mantle')
      WRITE(6,179)RES(3,4),RES(4,4)
179   FORMAT(14X,D11.4,9X,D11.4,7X,'crust')
      WRITE(6,180)RES(3,5),RES(4,5)
180   FORMAT(14X,D11.4,9X,D11.4,7X,'whole Earth')
      WRITE(6,181)RES(5,1),RES(6,1)
181   FORMAT(3X,'3'/14X,D11.4,9X,D11.4,7X,'inner core')
      WRITE(6,182)RES(5,2),RES(6,2)
182   FORMAT(14X,D11.4,9X,D11.4,7X,'outer core')
      WRITE(6,183)RES(5,3),RES(6,3)
183   FORMAT(14X,D11.4,9X,D11.4,7X,'mantle')
      WRITE(6,184)RES(5,4),RES(6,4)
184   FORMAT(14X,D11.4,9X,D11.4,7X,'crust')
      WRITE(6,185)RES(5,5),RES(6,5)
185   FORMAT(14X,D11.4,9X,D11.4,7X,'whole Earth')
      WRITE(6,186)RES(7,1),RES(8,1)
186   FORMAT(3X,'4'/14X,D11.4,9X,D11.4,7X,'inner core')
      WRITE(6,187)RES(7,2),RES(8,2)
187   FORMAT(14X,D11.4,9X,D11.4,7X,'outer core')
      WRITE(6,188)RES(7,3),RES(8,3)
188   FORMAT(14X,D11.4,9X,D11.4,7X,'mantle')
      WRITE(6,189)RES(7,4),RES(8,4)
189   FORMAT(14X,D11.4,9X,D11.4,7X,'crust')
      WRITE(6,190)RES(7,5),RES(8,5)
190   FORMAT(14X,D11.4,9X,D11.4,7X,'whole Earth')
C Write out discontinuities in radial displacements.
      WRITE(6,191)
191   FORMAT(/11X,'DISCONTINUITIES IN RADIAL DISPLACEMENTS')
      WRITE(6,192)
192   FORMAT(/1X,'source',5X,'discontinuity at ICB',5X,'discontinuity at
     1 CMB')
      WRITE(6,193)DISC(1,1),DISC(1,2)
193   FORMAT(3X,'1'/14X,D11.4,14X,D11.4)
      WRITE(6,194)DISC(2,1),DISC(2,2)
194   FORMAT(3X,'2'/14X,D11.4,14X,D11.4)
      WRITE(6,195)DISC(3,1),DISC(3,2)
195   FORMAT(3X,'3'/14X,D11.4,14X,D11.4)
      WRITE(6,196)DISC(4,1),DISC(4,2)
196   FORMAT(3X,'4'/14X,D11.4,14X,D11.4)
C Write out results to file.
      WRITE(12,197) FD,MUFD
197   FORMAT(1X,2D15.8)
C Write out integrals over radius.
      DO 198 I=1,4
        I2=2*I
        I2M1=I2-1
        WRITE(12,199)(RES(I2M1,J),RES(I2,J),J=1,5)
199     FORMAT(1X,2D20.8)
198   CONTINUE
C Write out discontinuities in radial displacements at the ICB and CMB.
      DO 200 I=1,4
        WRITE(12,199)DISC(I,1),DISC(I,2)
200   CONTINUE
      END
C
      SUBROUTINE SPMAT(N1,N2,N3,C,R,B,M1,M2,M3)
C SPMAT computes the N1 by N1 matrix C connecting
C second derivatives to function values at N1 nodal points R(I).
C B is the N3=N1-2 by N2=2*N1-2 augmented matrix used internally
C in the Gaussian elimination procedure.
C M1, M2, M3 are the respective maximum dimensions for N1, N2, N3
C specified in the main programme.
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION C(M1,M1),R(M1),B(M3,M2)
      N=N3+1
C Construct augmented matrix A.
      DO 11 I=1,N3
        DO 10 J=1,N2
          B(I,J)=0.D0
  10    CONTINUE
        B(I,I)=2.D0*(R(I+2)-R(I))
        B(I,I+N3)=1.D0
        IF(I.EQ.1.OR.I.EQ.N3) GO TO 11
        B(I,I-1)=R(I+1)-R(I)
        B(I,I+1)=R(I+2)-R(I+1)
  11  CONTINUE
C Apply conditions that third derivatives are constant
C in first two and last two intervals.
      B(1,2)=R(3)-R(2)-((R(2)-R(1))**2.D0)/(R(3)-R(2))
      B(N3,N3-1)=R(N3+1)-R(N3)-((R(N1)-R(N1-1))**2.D0)/(R(N1-1)-R(N1-2))
      B(1,1)=B(1,1)+((R(3)-R(1))*(R(2)-R(1)))/(R(3)-R(2))
      B(N3,N3)=B(N3,N3)
     1 +((R(N1)-R(N1-1))*(R(N1)-R(N1-2)))/(R(N1-1)-R(N1-2))
      NM1=N3-1
C Eliminate sub-diagonal.
      DO 13 I=1,NM1
        DIAG=B(I,I)
        ELL=B(I+1,I)
        DO 12 J=1,N2
          B(I,J)=B(I,J)/DIAG
          B(I+1,J)=B(I+1,J)-B(I,J)*ELL
  12    CONTINUE
  13  CONTINUE
C Eliminate super-diagonal.
      DO 15 I=1,NM1
        NO=N-I
        DIAG=B(NO,NO)
        ELL=B(NO-1,NO)
        DO 14 J=1,N2
          B(NO,J)=B(NO,J)/DIAG
          B(NO-1,J)=B(NO-1,J)-B(NO,J)*ELL
  14    CONTINUE
  15  CONTINUE
C Move inverse to position formerly occupied by A.
      DO 17 I=1,N3
        DO 16 J=1,N3
          B(I,J)=B(I,J+N3)
  16    CONTINUE
  17  CONTINUE
C Construct matrix B in N-2 by N augmentation space. 
      DO 19 I=1,N3
        DO 18 J=1,N1
          B(I,J+N3)=0.0D0
  18    CONTINUE
        IN=I+N3
        FACTOR=6.D0/((R(I+2)-R(I+1))*(R(I+1)-R(I)))
        B(I,IN)=FACTOR*(R(I+2)-R(I+1))
        B(I,IN+1)=-FACTOR*(R(I+2)-R(I))
        B(I,IN+2)=FACTOR*(R(I+1)-R(I))
  19  CONTINUE
C Multiply inverse of A into B to get N-2 by N version of C.
      DO 22 I=1,N3
        DO 21 J=1,N1
          SUM=0.D0
          DO 20 K=1,N3
            SUM=SUM+B(I,K)*B(K,J+N3)
  20      CONTINUE
          C(I+1,J)=SUM
  21    CONTINUE
  22  CONTINUE
C Add first and last rows to C.
      DO 23 J=1,N1
        C(1,J)=C(2,J)*(R(3)-R(1))/(R(3)-R(2))
     1 -C(3,J)*(R(2)-R(1))/(R(3)-R(2))
        C(N1,J)=C(N1-1,J)*(R(N1)-R(N1-2))/(R(N1-1)-R(N1-2))
     1 -C(N1-2,J)*(R(N1)-R(N1-1))/(R(N1-1)-R(N1-2))
  23  CONTINUE
      RETURN
      END
C
      SUBROUTINE INTPL(XV,YV,N1,C,R,S,M1)
C INTPL interpolates to find YV at XV, using a cubic spline,
C given the function S(I) at N1 nodal points R(I). C is the
C coefficient matrix, from subroutine SPMAT, relating the
C second derivatives to the function values at the nodal points.
C M1 is the maximum dimension for N1 specified in the main programme.
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION C(M1,M1),R(M1),S(M1)
      NPM1=N1-1
      J=NPM1
C Locate subinterval.
      DO 10 L=1,NPM1
        IF(XV.GE.R(L).AND.XV.LT.R(L+1)) J=L
  10  CONTINUE
C Find second derivative SUM1 at beginning of subinterval
C and second deivative SUM2 at end of subinterval.
      SUM1=0.D0
      SUM2=0.D0
      DO 11 K=1,N1
        SUM1=SUM1+C(J,K)*S(K)
        SUM2=SUM2+C(J+1,K)*S(K)
  11  CONTINUE
C Find interpolate YV at XV.
      YV=S(J)+(S(J+1)-S(J))*(XV-R(J))/(R(J+1)-R(J))
     1 +(XV-R(J))*(XV-R(J+1))*(SUM1*(2.D0-(XV-R(J))/(R(J+1)-R(J)))
     2 +SUM2*(1.D0+(XV-R(J))/(R(J+1)-R(J))))/6.D0
      RETURN
      END
C
      DOUBLE PRECISION FUNCTION Q1(AN,LAMBDA,MU)
C Finds value of function Q1 for series solution at geocentre.
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION LAMBDA,MU
      Q1=2.D0*AN*((AN+2.D0)*LAMBDA+(AN+1.D0)*MU)
      RETURN
      END
C
      DOUBLE PRECISION FUNCTION P1(AN,LAMBDA,MU)
C Finds value of function P1 for series solution at geocentre.
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION LAMBDA,MU
      P1=2.D0*AN*(AN*(AN+2.D0)*LAMBDA+(AN*(AN+2.D0)-1.D0)*MU)
      RETURN
      END
C
      DOUBLE PRECISION FUNCTION P2(AN,LAMBDA,MU)
C Finds value of function P2 for series solution at geocentre.
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION LAMBDA,MU
      P2=AN*(AN+5.D0+(AN+3.D0)*(LAMBDA/MU))
      RETURN
      END
C
      DOUBLE PRECISION FUNCTION Q2(AN,LAMBDA,MU)
C Finds value of function Q2 for series solution at geocentre.
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION LAMBDA,MU
      Q2=2.D0*(AN+1.D0)+(AN+3.D0)*(LAMBDA/MU)
      RETURN
      END
C
      SUBROUTINE MATRIX(A,RHO,LAMBDA,MU,PI,G,AN,WE)
C This subroutine calculates the 6X6 matrix A required to find
C higher terms in the power series expansions of the three
C fundamental solutions at the geocentre.
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(6,6),RHO(100)
      DOUBLE PRECISION LAMBDA(100),MU(100)
      BETA=1.D0/(LAMBDA(1)+2.D0*MU(1))
      DELTA=2.D0*MU(1)*(3.D0*LAMBDA(1)+2.D0*MU(1))*BETA
      GAMMA=(4.D0*PI*G*RHO(1)-2.D0*WE*WE)/3.D0
      AN1=AN*(AN+1.D0)
      EPSLN=4.D0*AN1*MU(1)*(LAMBDA(1)+MU(1))*BETA-2.D0*MU(1)
      A(1,1)=AN+3.D0+2.D0*LAMBDA(1)*BETA
      A(1,2)=-BETA
      A(1,3)=-AN1*LAMBDA(1)*BETA
      A(1,4)=0.D0
      A(1,5)=0.D0
      A(1,6)=0.D0
      A(2,1)=-2.D0*DELTA
      A(2,2)=AN+2.D0+4.D0*MU(1)*BETA
      A(2,3)=AN1*DELTA
      A(2,4)=-AN1
      A(2,5)=0.D0
      A(2,6)=0.D0
      A(3,1)=1.D0
      A(3,2)=0.D0
      A(3,3)=AN+2.D0
      A(3,4)=-1.D0/MU(1)
      A(3,5)=0.D0
      A(3,6)=0.D0
      A(4,1)=DELTA
      A(4,2)=LAMBDA(1)*BETA
      A(4,3)=-EPSLN
      A(4,4)=AN+5.D0
      A(4,5)=0.D0
      A(4,6)=0.D0
      A(5,1)=-3.D0*GAMMA-2.D0*WE*WE
      A(5,2)=0.D0
      A(5,3)=0.D0
      A(5,4)=0.D0
      A(5,5)=AN+4.D0
      A(5,6)=-1.D0
      A(6,1)=0.D0
      A(6,2)=0.D0
      A(6,3)=AN1*(3.D0*GAMMA+2.D0*WE*WE)
      A(6,4)=0.D0
      A(6,5)=-AN1
      A(6,6)=AN+5.D0
      RETURN
      END
C
      SUBROUTINE LINSOL(A,B,N,C,DET,M1,M2)
C LINSOL solves the linear system Ax=b by Gauss-Jordan elimination
C with partial pivoting. A is the N x N coefficient matrix
C and B is the N-length constant vector b. The inverse of A
C is returned in A, the solution vector x is returned in B,
C and the determinant of A is returned in DET. C is the N x 2N+1
C augmented matrix used internally. M1 is the maximum dimension
C for N specified in the main programme and M2=2M1+1 is
C the maximum dimension for 2N+1 used in the main programme.
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(M1,M1),B(M1),C(M1,M2),HOLD(M1,M2),
     1  IR(M1),IC(M1),IDET(M1)
C Construct the augmented matrix C.
      N2=2*N+1
      NP1=N+1
      NM1=N-1
C Insert the matrix A.
      DO 10 I=1,N
        DO 10 J=1,N
          C(I,J)=A(I,J)
 10   CONTINUE
C Insert the unknown vector.
      DO 11 J=1,N
        C(J,NP1)=B(J)
 11   CONTINUE
C Augment with the identity matrix.
      DO 12 I=1,N
        DO 12 J=1,N
          JJ=J+NP1
          C(I,JJ)=0.D0
          IF(I.EQ.J)C(I,JJ)=1.D0
 12   CONTINUE
C Begin Gauss-Jordan elimination.
      DET=1.D0
      DO 13 I=1,N
        IM1=I-1
C Search for largest pivot.
        PIVOT=0.D0
        DO 14 J=1,N
          DO 14 K=1,N
C Scan IR and IC vectors for invalid pivot subscripts.
            IF(I.EQ.1)GO TO 15
            DO 16 ISC=1,IM1
              DO 16 JSC=1,IM1
                IF(J.EQ.IR(ISC))GO TO 14
                IF(K.EQ.IC(JSC))GO TO 14
 16         CONTINUE
 15         IF(DABS(C(J,K)).LE.DABS(PIVOT))GO TO 14
            PIVOT=C(J,K)
            IR(I)=J
            IC(I)=K
 14     CONTINUE
        DET=DET*PIVOT
C Label next row.
        IRP1=IR(I)+1
C Label previous row.
        IRM1=IR(I)-1
C Divide row by pivot.
        DO 17 J=1,N2
             C(IR(I),J)=C(IR(I),J)/PIVOT
 17     CONTINUE
C Reduce elements below pivot to zero.
          IF(IRP1.GT.N)GO TO 18
          DO 19 K=IRP1,N
            ELL=C(K,IC(I))
            DO 20 J=1,N2
              C(K,J)=C(K,J)-C(IR(I),J)*ELL
 20         CONTINUE
 19       CONTINUE
 18       CONTINUE
C Reduce elements above pivot to zero.
          IF(IRM1.LT.1)GO TO 21
          DO 22 K=1,IRM1
            ELL=C(K,IC(I))
            DO 23 J=1,N2
              C(K,J)=C(K,J)-C(IR(I),J)*ELL
 23         CONTINUE
 22       CONTINUE
 21       CONTINUE
 13   CONTINUE     
C Unscramble equations and record row interchanges.
      DO 24 I=1,N
        IDET(IR(I))=IC(I)
        DO 24 J=1,N2
          HOLD(IC(I),J)=C(IR(I),J)
 24   CONTINUE
C Put equations back in original matrix.
      DO 25 I=1,N
        DO 25 J=1,N2
          C(I,J)=HOLD(I,J)
 25   CONTINUE
C Determine sign of determinant.
      ICH=0
      DO 26 I=1,NM1
        IP1=I+1
        DO 26 J=IP1,N
          IF(IDET(J).GE.IDET(I))GO TO 26
          IHOLD=IDET(J)
          IDET(J)=IDET(I)
          IDET(I)=IHOLD
          ICH=ICH+1
 26   CONTINUE
      IF((ICH/2)*2.NE.ICH)DET=-DET
C Put inverse in A and solution in B.
      DO 27 I=1,N
        B(I)=C(I,NP1)
        DO 27 J=1,N
          JJ=J+NP1
          A(I,J)=C(I,JJ)
 27   CONTINUE
      RETURN
      END
C
      SUBROUTINE REL(ERRMAX,DELTA,YIC,N)
C This subroutine is used in the variable step size Runge-Kutta integration
C of the fundamental solutions in the Inner Core. It finds the maximum
C relative error in the elements of the difference array, DELTA(6,3), compared
C to the elements of the fundamental solution array, YIC(6,3),for degree N.
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION DELTA(6,3),YIC(6,3)
C Zero relative error.
      ERRMAX=0.D0
      DO 1 I=1,6
C Branch for degee zero.
        IF(N.EQ.0)GO TO 2
        GO TO 3
C If degree is zero skip y_3 and y_4.
 2      IF(I.EQ.3.OR.I.EQ.4)GO TO 1
 3      DO 4 J=1,3
C Third solution is only fundamental solution for degree zero.
        IF(N.EQ.0.AND.J.NE.3)GO TO 4
        IF(YIC(I,J).EQ.0.D0)GO TO 4
        RATIO=DABS(DELTA(I,J)/YIC(I,J))
        ERRMAX=DMAX1(ERRMAX,RATIO)
 4      CONTINUE
 1    CONTINUE
      RETURN
      END
C
      SUBROUTINE YPRIME(R,Y,A,YP,N1,C,RM,RHOM,MUM,LAMBDAM,GZEROM,N,PI,G,
     1 IR,RMIN,ANGS,COR,WES)
C This subroutine finds the derivatives, YP(6,6), of the six
C fundamental solutions represented by the propagator matrix,
C Y(6,6), in the mantle and crust (shell), at an arbitrary radius R
C for degree N>0. For N=0, the third and fourth solution vectors can
C be ignored. In the inner core, it finds the derivatives of the
C fundamental solutions in terms of the z-variables, which are
C transformations of the y-variables. 
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION RM(100),RHOM(100),GZEROM(100),C(100,100),
     1 Y(6,6),YP(6,6),A(6,6)
      DOUBLE PRECISION MUM(100),LAMBDAM(100),MU,LAMBDA 
C Zero fill derivative matrix YP(I,J) and coefficient matrix A(I,J).
      DO 10 I=1,6
        DO 11 J=1,6
          YP(I,J)=0.D0
          A(I,J)=0.D0
 11     CONTINUE
 10   CONTINUE
C Set maximum number of points for interpolation.
      M1=100
C If R is less than or equal to RMIN, take properties to be those at the
C geocentre.
      IF(R.LE.RMIN) GO TO 12
C Interpolate Earth properties at radius R.
      CALL INTPL(R,RHO,N1,C,RM,RHOM,M1)
      CALL INTPL(R,MU,N1,C,RM,MUM,M1)
      CALL INTPL(R,LAMBDA,N1,C,RM,LAMBDAM,M1)         
      CALL INTPL(R,GZERO,N1,C,RM,GZEROM,M1)
      GO TO 13
C Assign geocentric Earth properties.
 12   RHO=RHOM(1)
      MU=MUM(1)
      LAMBDA=LAMBDAM(1)
      GZERO=(4.D0*PI*G*RHO-2.D0*WES)*R/3.D0
C Construct matrix A.
 13   AN=DFLOAT(N)
C Define factors.
      F1=AN*(AN+1.D0)
      F2=LAMBDA+2.D0*MU
      F3=3.D0*LAMBDA+2.D0*MU
      F4=RHO*GZERO
      F5=4.D0*PI*G*RHO
C Assign values to matrix elements.
      A(1,1)=-2.D0*LAMBDA/(F2*R)
      A(1,2)=1.D0/F2
      A(1,3)=F1*LAMBDA/(F2*R)
      A(1,4)=0.D0
      A(1,5)=0.D0
      A(1,6)=0.D0
      A(2,1)=4.D0*(-F4+MU*F3/(F2*R))/R-2.D0*RHO*WES-ANGS*RHO
      A(2,2)=-4.D0*MU/(F2*R)
      A(2,3)=F1*(F4-2.D0*MU*F3/(F2*R))/R+COR*RHO
      A(2,4)=F1/R
      A(2,5)=0.D0
      A(2,6)=-RHO
      IF(N.EQ.0.AND.IR.EQ.1) GO TO 14
      A(3,1)=-1.D0/R
      A(3,2)=0.D0
      A(3,3)=1.D0/R
      A(3,4)=1.D0/MU
      A(3,5)=0.D0
      A(3,6)=0.D0
      A(4,1)=(F4-2.D0*MU*F3/(F2*R))/R
C If degree is not zero, add Coriolis term.
      IF(N.NE.0) A(4,1)=A(4,1)+COR*RHO/F1
      A(4,2)=-LAMBDA/(F2*R)
      A(4,3)=2.D0*MU*((2.D0*F1-1.D0)*LAMBDA+(2.D0*F1-2.D0)*MU)/(F2*R*R)
     1 -ANGS*RHO
C If degree is not zero, add Coriolis term.
      IF(N.NE.0) A(4,3)=A(4,3)+COR*RHO/F1
      A(4,4)=-3.D0/R
      A(4,5)=-RHO/R
      A(4,6)=0.D0
 14   CONTINUE
      A(5,1)=F5
      A(5,2)=0.D0
      A(5,3)=0.D0
      A(5,4)=0.D0
      A(5,5)=0.D0
      A(5,6)=1.D0
      IF (N.EQ.0.AND.IR.EQ.1) GO TO 15
      A(6,1)=0.D0
      A(6,2)=0.D0
      A(6,3)=-F5*F1/R
      A(6,4)=0.D0
      A(6,5)=F1/(R*R)
 15   CONTINUE
      A(6,6)=-2.D0/R
C Branch to inner core calculation
      IF(IR.EQ.1) GO TO 16
C Multiply A into Y to find derivative matrix for shell.
      DO 17 I=1,6
        DO 18 J=1,6
          DO 19 K=1,6
            YP(I,J)=YP(I,J)+A(I,K)*Y(K,J)
 19       CONTINUE
 18     CONTINUE    
 17   CONTINUE
      GO TO 20
C Modify matrix A for inner core calculation.
 16   A(1,2)=A(1,2)/R
      A(2,1)=A(2,1)*R
      A(2,3)=A(2,3)*R
      A(2,6)=A(2,6)*R
      A(3,4)=A(3,4)/R
      A(4,1)=A(4,1)*R
      A(4,3)=A(4,3)*R
      A(4,5)=A(4,5)*R*R
      A(5,1)=A(5,1)/R
      A(5,6)=A(5,6)/R
      A(6,5)=A(6,5)*R
C Multiply A into Y to find YP for inner core.
      DO 21 I=1,6
        DO 22 J=1,3
          DO 23 K=1,6
            YP(I,J)=YP(I,J)+A(I,K)*Y(K,J)
 23       CONTINUE
 22     CONTINUE    
 21   CONTINUE
C Correct YP to give derivatives of z-variables.
C Set value of alpha plus one for first two solutions.
      ANM1=AN-1.D0
      DO 24 J=1,3
C Set value of alpha plus one for third solution
        IF(J.EQ.3) ANM1=AN+1.D0
        YP(1,J)=YP(1,J)-ANM1*Y(1,J)/R
        YP(2,J)=YP(2,J)-(ANM1-1.D0)*Y(2,J)/R
        YP(3,J)=YP(3,J)-ANM1*Y(3,J)/R
        YP(4,J)=YP(4,J)-(ANM1-1.D0)*Y(4,J)/R
        YP(5,J)=YP(5,J)-(ANM1+1.D0)*Y(5,J)/R
        YP(6,J)=YP(6,J)-ANM1*Y(6,J)/R
 24   CONTINUE
 20   CONTINUE
      RETURN
      END
C
      SUBROUTINE RK4(R,Y,A,K1,N1,C,RM,RHOM,MUM,LAMBDAM,GZEROM,N,
     1 PI,G,H,IR,RMIN,ANGS,COR,WES)
C This subroutine completes the last three steps of a fourth order
C Runge-Kutta integration, given the derivatives K1(6,6) at the starting
C radius R.
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION RM(100),RHOM(100),GZEROM(100),C(100,100),
     1 Y(6,6),YY(6,6),A(6,6)
      DOUBLE PRECISION MUM(100),LAMBDAM(100),
     1 K1(6,6),K2(6,6),K3(6,6),K4(6,6)
C Set half stepsize, one sixth stepsize.
      HH=0.5D0*H
      H6=H/6.D0
C Set number of solutions for shell.
      NS=6
C If region is inner core, reset number of solutions to three.
      IF(IR.EQ.1)NS=3
C Complete remaining steps of Runge-Kutta integration.
      DO 10 I=1,6
        DO 11 J=1,NS
          YY(I,J)=Y(I,J)+HH*K1(I,J)
 11     CONTINUE  
 10   CONTINUE
C Increment radius and solve for new derivatives.
      R=R+HH    
      CALL YPRIME(R,YY,A,K2,N1,C,RM,RHOM,MUM,LAMBDAM,GZEROM,N,PI,G,IR,
     1 RMIN,ANGS,COR,WES)
      DO 12 I=1,6
        DO 13 J=1,NS
          YY(I,J)=Y(I,J)+HH*K2(I,J)
 13     CONTINUE  
 12   CONTINUE
C Find new derivatives.
      CALL YPRIME(R,YY,A,K3,N1,C,RM,RHOM,MUM,LAMBDAM,GZEROM,N,PI,G,IR,
     1 RMIN,ANGS,COR,WES)
      DO 14 I=1,6
        DO 15 J=1,NS
          YY(I,J)=Y(I,J)+H*K3(I,J)
 15     CONTINUE  
 14   CONTINUE
C Increment radius and solve for new derivatives.    
      R=R+HH
      CALL YPRIME(R,YY,A,K4,N1,C,RM,RHOM,MUM,LAMBDAM,GZEROM,N,PI,G,IR,
     1 RMIN,ANGS,COR,WES)
      DO 16 I=1,6
        DO 17 J=1,NS
          Y(I,J)=Y(I,J)+H6*(K1(I,J)+2.D0*(K2(I,J)+K3(I,J))+K4(I,J))
 17     CONTINUE  
 16     CONTINUE
      RETURN
      END    
C
      SUBROUTINE YPRIMEOC(R,Y,A,YP,N1,C,RM,RHOM,GZEROM,PI,G,WES)
C This subroutine finds the derivatives, YP, of the two fundamental solutions in
C the fluid outer core.
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION RM(100),RHOM(100),GZEROM(100),C(100,100),
     1 Y(6,6),YP(6,6),A(6,6)
C Zero fill derivative matrix YP(I,J) and coefficient matrix A(I,J).
      DO 10 I=1,6
        DO 11 J=1,6
          YP(I,J)=0.D0
          A(I,J)=0.D0
 11     CONTINUE
 10   CONTINUE
C Set maximum number of points for interpolation.
      M1=100
C Interpolate Earth properties at radius R.
      CALL INTPL(R,RHO,N1,C,RM,RHOM,M1)
      CALL INTPL(R,GZERO,N1,C,RM,GZEROM,M1)
C Construct matrix A.
      A(1,1)=4.D0*PI*G*RHO/GZERO
      A(1,2)=1.D0
      A(2,1)=-16.D0*PI*G*RHO*(1.D0+R*WES/(2.D0*GZERO))/(GZERO*R)
     1 +6.D0/(R*R)
      A(2,2)=-4.D0*PI*G*RHO/GZERO-2.D0/R
C Multiply A into Y to find derivative matrix YP.
      DO 12 I=1,2
        DO 13 J=1,2
          DO 14 K=1,2
            YP(I,J)=YP(I,J)+A(I,K)*Y(K,J)
 14       CONTINUE
 13     CONTINUE
 12   CONTINUE
      RETURN
      END
C
      SUBROUTINE RK4OC(R,Y,A,K1,N1,C,RM,RHOM,GZEROM,PI,G,H,WES)
C This subroutine completes the last three steps of a fourth order
C Runge-Kutta integration, given the derivatives K1(6,6) at the starting
C radius R.
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION RM(100),RHOM(100),GZEROM(100),C(100,100),
     1 Y(6,6),YY(6,6),A(6,6)
      DOUBLE PRECISION K1(6,6),K2(6,6),K3(6,6),K4(6,6)
C Set half stepsize, one sixth stepsize.
      HH=0.5D0*H
      H6=H/6.D0
C Complete remaining steps of Runge-Kutta integration.
      DO 10 I=1,2
        DO 11 J=1,2
          YY(I,J)=Y(I,J)+HH*K1(I,J)
 11     CONTINUE
 10   CONTINUE
C Increment radius and solve for new derivatives.
      R=R+HH
      CALL YPRIMEOC(R,YY,A,K2,N1,C,RM,RHOM,GZEROM,PI,G,WES)
      DO 12 I=1,2
        DO 13 J=1,2
          YY(I,J)=Y(I,J)+HH*K2(I,J)
 13     CONTINUE
 12   CONTINUE
C Find new derivatives.
      CALL YPRIMEOC(R,YY,A,K3,N1,C,RM,RHOM,GZEROM,PI,G,WES)
      DO 14 I=1,2
        DO 15 J=1,2
          YY(I,J)=Y(I,J)+H*K3(I,J)
 15     CONTINUE
 14   CONTINUE
C Increment radius and solve for new derivatives.
      R=R+HH
      CALL YPRIMEOC(R,YY,A,K4,N1,C,RM,RHOM,GZEROM,PI,G,WES)
      DO 16 I=1,2
        DO 17 J=1,2
          Y(I,J)=Y(I,J)+H6*(K1(I,J)+2.D0*(K2(I,J)+K3(I,J))+K4(I,J))
 17     CONTINUE
 16   CONTINUE
      RETURN
      END
