C                             PROGRAMME SHIFT.FOR
C SHIFT.FOR calculates the secular polar shift for a given earthquake.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION RES(8,5),DISC(4,2),TIDS(3,3),TISS(3,3)
      DOUBLE PRECISION L,MUFD
C Define pi.
      PI= 3.141592653589793D0
C Open results.dat file, containing integrals of radial and transverse coefficients
C u and v, alongwith discontinuities in radial displacements at the ICB and the CMB.
      OPEN(UNIT=1,FILE='results.dat',STATUS='OLD')
C Read in focal depth and shear modulus there.
      READ(1,1)FD,MUFD
 1    FORMAT(1X,2D15.8)
C Read in integrals of radial and transverse coefficients u and v.
      DO 2 I=1,4
           I2=2*I
           I2M1=I2-1
           READ(1,3)(RES(I2M1,J),RES(I2,J),J=1,5)
 3         FORMAT(1X,2D20.8)
 2    CONTINUE
C Read in discontinuites in radial displacements at the ICB and the CMB.
      DO 4 I=1,4
           READ(1,3)DISC(I,1),DISC(I,2)
 4         CONTINUE
C Write results to screen.
      WRITE(6,5)FD
 5    FORMAT(/5X,'Results for focal depth =',1X,F5.1,1X,'km.')
      WRITE(6,6)
 6    FORMAT(/11X,'INTEGRALS OVER RADIUS')
      WRITE(6,7)
 7    FORMAT(/1X,'source',5X,'integral over u',5X,'integral over v',5X,
     1     'region')
      WRITE(6,8)RES(1,1),RES(2,1)
 8    FORMAT(3X,'1'/14X,D11.4,9X,D11.4,7X,'inner core')
      WRITE(6,9)RES(1,2),RES(2,2)
 9    FORMAT(14X,D11.4,9X,D11.4,7X,'outer core')
      WRITE(6,10)RES(1,3),RES(2,3)
 10   FORMAT(14X,D11.4,9X,D11.4,7X,'mantle')
      WRITE(6,11)RES(1,4),RES(2,4)
 11   FORMAT(14X,D11.4,9X,D11.4,7X,'crust')
      WRITE(6,12)RES(1,5),RES(2,5)
 12   FORMAT(14X,D11.4,9X,D11.4,7X,'whole Earth')
      WRITE(6,13)RES(3,1),RES(4,1)
 13   FORMAT(3X,'2',/,14X,D11.4,9X,D11.4,7X,'inner core')
      WRITE(6,14)RES(3,2),RES(4,2)
 14   FORMAT(14X,D11.4,9X,D11.4,7X,'outer core')
      WRITE(6,15)RES(3,3),RES(4,3)
 15   FORMAT(14X,D11.4,9X,D11.4,7X,'mantle')
      WRITE(6,16)RES(3,4),RES(4,4)
 16   FORMAT(14X,D11.4,9X,D11.4,7X,'crust')
      WRITE(6,17)RES(3,5),RES(4,5)
 17   FORMAT(14X,D11.4,9X,D11.4,7X,'whole Earth')
      WRITE(6,18)RES(5,1),RES(6,1)
 18   FORMAT(3X,'3',/,14X,D11.4,9X,D11.4,7X,'inner core')
      WRITE(6,19)RES(5,2),RES(6,2)
 19   FORMAT(14X,D11.4,9X,D11.4,7X,'outer core')
      WRITE(6,20)RES(5,3),RES(6,3)
 20   FORMAT(14X,D11.4,9X,D11.4,7X,'mantle')
      WRITE(6,21)RES(5,4),RES(6,4)
 21   FORMAT(14X,D11.4,9X,D11.4,7X,'crust')
      WRITE(6,22)RES(5,5),RES(6,5)
 22   FORMAT(14X,D11.4,9X,D11.4,7X,'whole Earth')
      WRITE(6,23)RES(7,1),RES(8,1)
 23   FORMAT(3X,'4',/,14X,D11.4,9X,D11.4,7X,'inner core')
      WRITE(6,24)RES(7,2),RES(8,2)
 24   FORMAT(14X,D11.4,9X,D11.4,7X,'outer core')
      WRITE(6,25)RES(7,3),RES(8,3)
 25   FORMAT(14X,D11.4,9X,D11.4,7X,'mantle')
      WRITE(6,26)RES(7,4),RES(8,4)
 26   FORMAT(14X,D11.4,9X,D11.4,7X,'crust')
      WRITE(6,27)RES(7,5),RES(8,5)
 27   FORMAT(14X,D11.4,9X,D11.4,7X,'whole Earth')
      WRITE (6,28)
 28   FORMAT(/11X,'DISCONTINUITIES IN RADIAL DISPLACEMENTS')
      WRITE(6,29)
 29   FORMAT(/1X,'source',5X,'discontinuity at ICB',5X,'discontinuity at
     1 CMB')
      WRITE(6,30)DISC(1,1),DISC(1,2)
 30   FORMAT(3X,'1',/14X,D11.4,14X,D11.4)
      WRITE(6,31)DISC(2,1),DISC(2,2)
 31   FORMAT(3X,'2',/14X,D11.4,14X,D11.4)
      WRITE(6,32)DISC(3,1),DISC(3,2)
 32   FORMAT(3X,'3',/14X,D11.4,14X,D11.4)
      WRITE(6,33)DISC(4,1),DISC(4,2)
 33   FORMAT(3X,'4',/14X,D11.4,14X,D11.4)
C Enter fault parameters.
      WRITE(6,34)
 34   FORMAT(1X,'Enter dip angle in degrees.')
      READ(5,*)DIP
      WRITE(6,35)
 35   FORMAT(1X,'Enter slip in metres.')
      READ(5,*)SLIP
      WRITE(6,36)
 36   FORMAT(1X,'Enter rake angle in degrees.')
      READ(5,*)RAKE
      WRITE(6,37)
 37   FORMAT(1X,'Enter upper edge of fault, down- dip in kilometres.')
      READ(5,*)DU
      DU=DU*1.D3
      WRITE(6,38)
 38   FORMAT(1X,'Enter lower edge of fault, down-dip in kilometres.')
      READ(5,*)DL
      DL=DL*1.D3
      WRITE(6,39)
 39   FORMAT(1X,'Enter fault half-length in kilometres.')
      READ(5,*)L
      L=L*1.D3
C Convert dip and rake to radians.
      DIP=DIP*PI/180.D0
      RAKE=RAKE*PI/180.D0
C Find double dip angle.
      DDIP=DIP+DIP
C Find dip-slip and strike-slip components.
      SLIPD=SLIP*DCOS(RAKE)
      SLIPS=SLIP*DSIN(RAKE)
C Find fault area.
      AREA=2.D0*L*(DL-DU)
C Find actual focal depth in kilometres.
      FDA=(DU+DL)/2.D0
      FDA=(DU+DL)*(DSIN(DIP)*1.D-3)/2.D0
C Write out actual focal depth.
      WRITE(6,42)FDA
 42   FORMAT(/5X,'Actual focal depth =',1X,F5.2,1X,'km.')
C Determine if it is necessary to rerun sources.for for actual focal depth.
      RATIO=(FD-FDA)/FDA
      IF(RATIO.GT.1.D-4) GO TO 40
C Find seismic moment.
      SM=MUFD*SLIP*AREA
C Write out seismic moment.
      WRITE(6,44)SM
 44   FORMAT(/5X,'Siesmic moment =',1X,D10.3,1X,'N m.')
c      write(6,100)fda,sm
c 100  format(1x,2d15.8)
C Construct dip-slip inertia tensor.
      TIDS(1,1)=4.D0*(RES(1,5)+RES(3,5))*DSIN(DDIP)/15.D0
      TIDS(1,1)=TIDS(1,1)+4.D0*(RES(2,5)+RES(4,5))*DSIN(DDIP)/5.D0
      TIDS(2,2)=4.D0*(RES(1,5)+RES(3,5)-RES(5,5))*DSIN(DDIP)/15.D0
      TIDS(2,2)=TIDS(2,2)+4.D0*(RES(2,5)+RES(4,5)-RES(6,5))*DSIN(DDIP)
     1/15.D0
      TIDS(3,3)=4.D0*(RES(5,5)-2.D0*RES(1,5)-2.D0*RES(3,5))*DSIN(DDIP)
     1/5.D0
      TIDS(1,2)=0.0D0
      TIDS(1,3)=0.0D0
      TIDS(2,3)=2.D0*(6.D0*RES(1,5)-RES(5,5)-RES(7,5))*DCOS(DDIP)/15.D0
      TIDS(2,3)=TIDS(2,3)+2.D0*(6.D0*RES(2,5)-RES(6,5)-RES(8,5))*
     1 DCOS(DDIP)/5.D0
C Multiply by 4pi x slip x fault area.
      DO 45 I=1,3
          DO 46 J=1,3
               TIDS(I,J)=4.D0*PI*SLIPD*AREA*TIDS(I,J)
 46         CONTINUE
 45   CONTINUE
C Construct strike-slip inertia tensor.
      TISS(1,1)=0.0D0
      TISS(2,2)=0.0D0
      TISS(3,3)=0.0D0
      TISS(1,2)=4.D0*RES(5,5)*DSIN(DIP)/15.D0+4.D0*RES(6,5)*
     1 DSIN(DIP)/5.D0
      TISS(1,3)=2.D0*(RES(5,5)+RES(7,5)-6.D0*RES(1,5))*DCOS(DIP)/15.D0
      TISS(2,3)=0.0D0
C Multiply by 4pi x slip X fault area.
      DO 47 I=1,3
          DO 48 J=1,3
               TISS(I,J)=4.D0*PI*SLIPS*AREA*TISS(I,J)
 48         CONTINUE
 47   CONTINUE  
      GO TO 43
 40   CONTINUE
      WRITE(6,41)
 41   FORMAT(1X,'Rerun sources.for for actual focal depth.')
 43   CONTINUE
      END
 
