subroutine sphere c Response Q(n,p) for potential fields c Q(n) is the ratio of internal to external SHA expansion c coefficient to degree n and period T C use coef_module,only: nm,mm, | qnnmm ! correlates external to internal coeff. from cm3e ! implicit real*8 (a-h,o-z) save parameter(nmax=nm,nm1=2*nmax) COMPLEX*16 CK(0:nm1),CG(0:nm1),CF(0:nm1),ZT, 1 ZU,BETAM,C0,C1,C,G1,G2,RG(nmax),Q(nmax),CN,c01 REAL*8 RA(nm1) real*8 r(20),dd(20) integer im PI=4*ATAN(1.d0) A=1.d0 ! WRITE (*,'(a)') 'Period in hours ? ' ! READ (*,*) TP ! do im = 1,mm ! loop over mm tp = 24./float(im) tp=tp*3600.d0 FR=1.d0/TP nn=nmax ! maximum degree of Q(n) C c conductivity model OLSEN (1998), based on Sq/Dst in Europa C N5=4 ! number of layers dd(1)=520.d0 ! thickness of layer 1 (in km) dd(2)=126.d0 dd(3)=145 RA(1)=6371.d0 ! conversion to radius do n=2,n5 ra(n)=ra(n-1)-dd(n-1) enddo R(1)=76.d0 ! resistivity of layer 1 (in Omega*m) R(2)=4.5d0 R(3)=5.9d0 r(4)=.662 ! resistivity of layer 4 (inner sphere) DO 30 J=1,N5 20 CK(J)=(1.d0,1.d0)*SQRT(.4d0*PI*PI*FR/R(J)) 30 CONTINUE ZT=CK(N5)*RA(N5) TZ=ABS(ZT) 50 DO 200 N=1,NN T=N IF (A.EQ.0.d0) GO TO 100 IF (N.GT.5) GO TO 60 IF (TZ.GE.1.d0) GO TO 70 IF (TZ.LE.0.2d0) GO TO 80 IF (N.LE.3) GO TO 70 60 CALL FUCH (N,ZT,G1) GO TO 90 70 CALL STOK (N,ZT,G1) GO TO 90 80 CALL ZEMI (N,ZT,G1) 90 CG(N5)=G1 CF(N5)=G1 C=CG(N5)/CK(N5) 100 BETAM=C*CK(N5-1) C1=BETAM IF (N5.EQ.1) RG(N)=CG(N5) DO 190 L=1,2 BETAM=C1 IF ((R(N5).LE.0.d0).OR.(L.EQ.1)) GO TO 110 C0=CK(N5-1)/CK(N5) 110 IF (N5.EQ.1) GO TO 190 DO 180 J=2,N5 M=N5+1-J ZT=CK(M)*RA(M) TZ=ABS(ZT) ZU=CK(M)*RA(M+1) UZ=ABS(ZU) RO=RA(M+1)/RA(M) IF (N.GT.5) GO TO 120 IF (UZ.GT.1.d0) GO TO 130 IF (TZ.LE.0.2d0) GO TO 140 IF (N.LE.3) GO TO 130 120 CALL FUCHUR (N,ZT,ZU,BETAM,G1,G2) GO TO 150 130 CALL STOKES (N,ZT,ZU,BETAM,G1,G2) GO TO 150 140 CALL ZETMIN (N,ZT,RO,BETAM,G1,G2) 150 IF (L-1) 160,160,170 160 CG(M)=G1 C=CG(M)/CK(M) BETAM=C*CK(M-1) IF (M.EQ.1) RG(N)=CG(M) GO TO 180 170 CF(M)=G1 C0=CK(M-1)/CK(M) BETAM=1.d0/(C0/CF(M)+(T/ZT)*(1.d0/C0-C0)) 180 CONTINUE 190 CONTINUE 200 CONTINUE ! write(6,'(a3,i4)') 'm= ',im ! print*, 'Flat Earth Response C_0' call wait(r,dd,n5,tp,c01) ! Model response plane Earth ! WRITE (*,'(1x,i5,2f8.0,3f8.4,i6,f8.4)') 0,c01 ! C_0 for plane Earth ! print* ! print*, 'Response of layered sphere' ! print*, ' n C [km] Re{Q} Im{Q} |Q| arg{Q}' ! print* DO 210 N=1,NN Q(N)=N/(N+1.d0)*(1.d0-RG(N)*(2*N+1)/(RA(1)*CK(1))) CN=6371.d0/(N+1)*(1-(N+1.d0)/N*Q(N))/(1.d0+Q(N)) ! WRITE (*,'(1x,i5,2f8.0,3f8.4,i6)') N,cn,Q(N) ! * ,abs(q(n)),nint(carg(q(n))) qnnmm(n,im) = q(n) 210 CONTINUE qnnmm(0,:) = 0. qnnmm(:,0) = 0. enddo ! loop over mm c print* c write(*,301) (real(q(n)),n=1,nn) c print* c write(*,301) (imag(q(n)),n=1,nn) c 301 format(40(5x,1h*,4x,4(e12.7,1h,)/)) END !----------------------------------------------------------------------- c g1 einer homogenen kugel mit wunderformel SUBROUTINE FUCH (N,X,G1) implicit real*8 (a-h,o-z) save COMPLEX*16 R1,R2,D,P,Z COMPLEX*16 X,G1 Z=X A=N+0.5D0 R1=SQRT(Z*Z+A*A) R2=SQRT(Z*Z+(A-1)*(A-1)) D=SQRT(R2/R1) D=D*((A-1+R2)/(A+R1))**(A-1) D=D*(Z/(A+R1)) D=D*(P(A,R1)/P(A-1,R2)) D=D*EXP(R1-R2) G1=D RETURN END SUBROUTINE FUCHUR (N,XT,XU,BETA,G1,G2) implicit real*8 (a-h,o-z) save COMPLEX*16 G11,G22,R,D1,D2,D3,B(5),GM,GN,D4,P,FN,FM,RN,RM,BETAM, 1 ZT,ZU,X1(5),X2(5),X COMPLEX*16 XU,XT,G1,G2,BETA INTEGER L1(5),L2(5) BETAM=BETA ZT=XT ZU=XU L1(1)=N L1(2)=N L1(3)=N-1 L1(4)=N L1(5)=N-1 L2(1)=N-1 L2(2)=N L2(3)=N-1 L2(4)=N-1 L2(5)=N X1(1)=ZT X1(2)=ZU X1(3)=ZT X1(4)=ZU X1(5)=ZT X2(1)=ZU X2(2)=ZT X2(3)=ZU X2(4)=ZT X2(5)=ZT DO 50 K=1,5 A=L1(K)+0.5D0 X=X1(K) R=SQRT(X**2+A**2) D1=-0.5d0*(LOG(X)+LOG(R)) D2=A*(LOG(A+R)-LOG(X)) A=-A D4=1.D0*P(A,R) D3=LOG(D4) A=-A GN=D1+D2+D3 D4=1.D0*P(A,R) D3=LOG(D4) FN=LOG(0.5d0)+D1-D2+D3 RN=R A=L2(K)+0.5d0 X=X2(K) R=SQRT(A**2+X**2) D1=-0.5d0*(LOG(X)+LOG(R)) D2=A*(LOG(A+R)-LOG(X)) A=-A D3=LOG(P(A,R)) A=-A GM=D1+D2+D3 D3=LOG(P(A,R)) FM=LOG(0.5d0)+D1-D2+D3 RM=R GN=GN+FM+RM-RN GM=GM+FN+RN-RM GN=EXP(GN) GM=EXP(GM) IF (L1(K)-2*(L1(K)/2)) 10,20,10 10 GN=-GN 20 IF (L2(K)-2*(L2(K)/2)) 30,40,30 30 GM=-GM 40 B(K)=GN-GM 50 CONTINUE B(1)=BETAM*B(1) B(3)=BETAM*B(3) G11=(B(1)+B(2))/(B(3)+B(4)) G22=(-B(4)+B(2)/G11)/B(5) G1=G11 G2=G22 RETURN END SUBROUTINE STOK (N,X,G1) c g1 einer homogenen kugel mit hyp.fu implicit real*8 (a-h,o-z) save COMPLEX*16 D0,D1,D2,G2,S1(2),S2(2),Z,CT COMPLEX*16 X,G1 INTEGER Y(2) Z=X Y(1)=-1 Y(2)=1 DO 30 L=1,2 K=Y(L) D0=(0.D0,0.D0) D1=(1.D0,0.D0) M=N IF (K.EQ.1) M=M-1 IF (M.EQ.0) GO TO 20 DO 10 I=1,M D2=D0-(2.d0*I+K)*D1/Z D0=D1 10 D1=D2 20 S2(L)=D1 S1(L)=D0 30 CONTINUE CT=(EXP(Z)-EXP(-Z))/(EXP(Z)+EXP(-Z)) G2=S2(1)*CT+S2(2) G1=G2/(S1(1)*CT+S1(2)) RETURN END SUBROUTINE STOKES (N,XT,XU,BETA,G1,G2) implicit real*8 (a-h,o-z) save COMPLEX*16 T,C,D0,D1,D2,G(2,2),H(2,2),B(5),G11,ZT,ZU,BETAM,X(2),Z COMPLEX*16 XU,XT,BETA,G1,G2 INTEGER Y(2) BETAM=BETA ZT=XT ZU=XU X(1)=ZU X(2)=ZT Y(1)=-1 Y(2)=1 DO 40 J=1,2 DO 30 L=1,2 K=Y(L) D0=(0.D0,0.D0) D1=(1.D0,0.D0) M=N IF (K.EQ.1) M=M-1 IF (M.EQ.0) GO TO 20 DO 10 I=1,M D2=D0-(2*I+K)*D1/X(J) D0=D1 10 D1=D2 20 G(L,J)=D0 H(L,J)=D1 30 CONTINUE 40 CONTINUE Z=ZT-ZU T=(EXP(Z)-EXP(-Z))/(EXP(Z)+EXP(-Z)) C=(EXP(Z)+EXP(-Z))/2.d0 B(1)=T*(H(2,2)*G(2,1)-H(1,2)*G(1,1)) B(1)=B(1)+H(1,2)*G(2,1)-H(2,2)*G(1,1) B(1)=BETAM*B(1) B(2)=T*(H(1,2)*H(1,1)-H(2,2)*H(2,1)) B(2)=B(2)-H(2,1)*H(1,2)+H(1,1)*H(2,2) B(3)=T*(-G(1,2)*G(1,1)+G(2,2)*G(2,1)) B(3)=B(3)+G(1,2)*G(2,1)-G(2,2)*G(1,1) B(3)=BETAM*B(3) B(4)=T*(G(1,2)*H(1,1)-G(2,2)*H(2,1)) B(4)=B(4)+G(2,2)*H(1,1)-G(1,2)*H(2,1) B(5)=H(2,2)*G(1,2)-H(1,2)*G(2,2) G11=(B(1)+B(2))/(B(3)+B(4)) G2=(ZT/ZU)*C*(-B(4)+B(2)/G11)/B(5) G1=G11 RETURN END c g1 einer homogenen kugel fuer betrag z kleiner 0.2d0 SUBROUTINE ZEMI (N,Z,G1) implicit real*8 (a-h,o-z) save COMPLEX*16 Z,G1,Y A=1.d0/(2.d0*N+1.d0) B=1.d0/(2.d0*N+3.d0) C=1.d0/(2.d0*N+5.d0) Y=.5d0*Z*Z G1=Z*A*(1.d0+Y*B*(1.d0+Y*.5d0*C))/(1.d0+Y*A*(1.d0+Y*.5d0*B)) RETURN END SUBROUTINE ZETMIN (N,ZT,RO,BETAM,G1,G2) C g1 einer geschichteten kugel fuer kleine argumente implicit real*8 (a-h,o-z) save COMPLEX*16 ZU,ZT,BETAM,G1,G2,D(5),G(2,2),H(2,2),X(2) REAL*8 RO,A,B,C INTEGER Y(2) Y(1)=N-1 Y(2)=N ZU=ZT*RO X(1)=ZU X(2)=ZT DO 20 I=1,2 DO 10 K=1,2 A=2.d0*Y(K) G(I,K)=1.d0+.5d0*X(I)**2/(A+3.d0)+(.5d0*X(I)**2)**2/ * (2.d0*(A+3.d0)*(A+5)) 10 H(I,K)=1.d0+.5d0*X(I)**2/(1-A)+(.5d0*X(I)**2)**2/ * (2.d0*(1.d0-A)*(3.d0-A)) 20 CONTINUE A=2.d0*N+1.d0 B=2.d0*N-1.d0 C=2.d0*N D(1)=(G(1,1)*RO**B*H(2,2)/ZT**2+G(2,2)*H(1,1)/(A*B))*BETAM D(2)=(G(2,2)*H(1,2)/RO-RO**C*G(1,2)*H(2,2))/(A*ZT) D(3)=(G(2,1)*H(1,1)-RO**B*G(1,1)*H(2,1))*BETAM/(B*ZT) D(4)=G(2,1)*H(1,2)/(RO*ZT**2)+RO**C*G(1,2)*H(2,1)/(A*B) D(5)=-RO**N*(H(2,1)*G(2,2)/(A*B)+H(2,2)*G(2,1)/ZT**2) G1=(D(1)+D(2))/(D(3)+D(4)) G2=(D(2)/G1-D(4))/D(5) RETURN END C POLYNOM NACH DOVER 9.3d0.9d0.d0 fuer komplexes Z COMPLEX*16 FUNCTION P(A,R) implicit real*8 (a-h,o-z) save COMPLEX*16 R,V,U1,U2,U3,U4,UU,UG,F REAL*8 A,D,E(6),O(8) E(1)=.12500000000000D000 E(2)=-.20833333333333D000 E(3)=.73242187500000D-001 E(4)=-.89121093750000D000 E(5)=.18464626736111D+001 E(6)=-.10258125964506D+001 O(1)=.70312500000000D-001 O(2)=-.40104166666667D000 O(3)=.33420138888889D000 O(4)=.11215209960938D000 O(5)=-.23640869140625D001 O(6)=.87891235351562D001 O(7)=-.11207002616223D002 O(8)=.46695844234262D001 V=1.D0/(R*R) D=A*A F=D*V U1=E(1)*V*R+E(2)*V*D/R U3=E(6)*F+E(5) DO 10 I=4,3,-1 10 U3=U3*F+E(I) U3=U3*V/R UU=U1+U3 U2=O(1)*V+O(2)*V*F+O(3)*V*F*F U4=O(8)*F+O(7) DO 20 I=6,4,-1 20 U4=U4*F+O(I) U4=U4*V*V UG=U2+U4 IF (A) 30,50,40 30 P=1+UG-UU GO TO 50 40 P=1+UG+UU 50 RETURN END subroutine wait(rho,d,mmax,t,c0) c c Calculates C-response (in km) for a layered Earth-model c with WAIT' s formula. c c Input: mmax: number of layers c rho(m): resistivities [ohm*m] c d(m): layer thickness [km] c T: Period [s] c Output C: C-Response [km] c c May 1994, Nils Olsen. c parameter(nmax=20) implicit real*8 (a-h,o-z) save real*8 rho(nmax),d(nmax) real*8 t complex*16 c0 complex*16 k(nmax),g(nmax),ct,ctanh,z ctanh(z)=(exp(z)-exp(-z))/(exp(z)+exp(-z)) if(mmax.gt.nmax) stop 'MMAX > NMAX in UP WAIT' do 10 m=1,mmax 10 k(m)=sqrt((0.d0,7.8956835d0)/t/rho(m)) ! in [km] g(mmax)=1.d0 do 20 m=mmax-1,1,-1 ct=ctanh(k(m)*d(m)) 20 g(m)=(k(m)*g(m+1)+k(m+1)*ct)/(k(m+1)+k(m)*g(m+1)*ct) c0=g(1)/k(1) ! C-Response in [km] return end FUNCTION CARG (Z) implicit real*8 (a-h,o-z) save COMPLEX*16 Z if(abs(z).eq.0) then carg=0.d0 return endif CARG=180.d0/3.1415926d0*ATAN2(IMAG(Z),DBLE(Z)) IF (CARG.LT.0.d0) CARG=CARG+360.d0 END