C
C C ************************************************************** C SUBROUTINE MODEL(LU,LOU,MM,MTEXT) C C APPROXIMATION OF INTERFACES AND VELOCITY DISTRIBUTION C IN INDIVIDUAL LAYERS (ISOVELOCITY DISCONTINUITIES) C CHARACTER*80 MTEXT DIMENSION IPR(101),X(100),FX(100),AUX(12),VEL(6),Y(4) COMMON/MEDIM/V1(19),V2(19),NVS(19),PTOS(19),MAUX COMMON/INTRF/A1(30,20),B1(29,20),C1(30,20),D1(29,20),X1(30,20), 1BRD(2),III(30,20),NPNT(20),NINT COMMON/AUXI/INTR,INT1,IOUT,KRE,IREFR,LAY,ITYPE,NDER,IPRINT,MPRINT, 1NTR,IMET COMMON/DEN/RHO1(19),RHO2(19),NRHO COMMON/ABSR/QP1(19),QP2(19),QP3(19),QS1(19),QS2(19),QS3(19), 1NQP(19),NQS(19),NABS C MI=1 IF(MM.NE.MI)MM=-1 IF(MM.LT.0)RETURN C C INPUT DATA - INTERFACES C READ(LU,*)NINT,(NPNT(I),I=1,NINT) IF(MPRINT.GE.0)WRITE(LOU,102)NINT,(NPNT(I),I=1,NINT) DO 11 I=1,NINT NC=NPNT(I) READ(LU,*)(X(J),FX(J),III(J,I),J=1,NC) IF(MPRINT.GE.0)WRITE(LOU,103)(X(J),FX(J),III(J,I),J=1,NC) C C DETERMINATION OF COEFFICIENTS OF CUBIC PARABOLAS C APPROXIMATING INTERFACES C DO 1 J=1,NC X1(J,I)=X(J) 1 A1(J,I)=FX(J) J1=1 NMIN=1 2 DO 3 J=J1,NC J2=J IF(III(J,I))4,3,6 3 CONTINUE 4 IF(NMIN.EQ.J2)GO TO 5 FX(NMIN)=A1(NMIN,I) CALL SPLIN(X,FX,NMIN,J2) KEY=0 GO TO 8 5 IF(J2.EQ.NC)GO TO 11 J1=J2+1 NMIN=J2 GO TO 2 6 IF(NMIN.EQ.J2)GO TO 7 FX(NMIN)=A1(NMIN,I) CALL SPLIN(X,FX,NMIN,J2) KEY=1 GO TO 8 7 IN=III(J2,I) X1(J2,I)=X1(IN,I-1) A1(J2,I)=A1(IN,I-1) B1(J2,I)=B1(IN,I-1) C1(J2,I)=C1(IN,I-1) D1(J2,I)=D1(IN,I-1) IF(J2.EQ.(NC-1))GO TO 11 J1=J2+1 NMIN=J1 GO TO 2 8 IF((J2-NMIN).EQ.1)GO TO 10 J3=J2-1 DO 9 J=NMIN,J3 H=X(J+1)-X(J) D=(A1(J+1,I)-A1(J,I))/H D1(J,I)=(FX(J+1)-FX(J))/(3.*H) C1(J,I)=FX(J) 9 B1(J,I)=D-.333333*H*(FX(J+1)+2.*FX(J)) IF(KEY)5,5,7 10 D1(NMIN,I)=0. C1(NMIN,I)=0. B1(NMIN,I)=(A1(J2,I)-A1(NMIN,I))/(X(J2)-X(NMIN)) IF(KEY)5,5,7 11 CONTINUE NC=NINT-1 C C INPUT DATA - VELOCITY DISTRIBUTION IN INDIVIDUAL LAYERS WRITE(LOU,117) READ(LU,*)(V1(I),I=1,NC) READ(LU,*)(V2(I),I=1,NC) WRITE(LOU,104)(V1(I),I=1,NC) WRITE(LOU,104)(V2(I),I=1,NC) C C INPUT DATA - DENSITIES AND S VELOCITIES C READ(LU,*)NRO,(NVS(I),I=1,NC),NABS IF(MPRINT.GE.0)WRITE(LOU,102)NRO,(NVS(I),I=1,NC),NABS IF(NRO.EQ.0)GO TO 30 READ(LU,*)(RHO1(I),RHO2(I),I=1,NC) IF(MPRINT.GE.0)WRITE(LOU,104)(RHO1(I),RHO2(I),I=1,NC) 30 IF(NABS.EQ.0)GO TO 34 DO 35 I=1,NC READ(LU,*)NQP(I),NQS(I),QP1(I),QP2(I),QP3(I),QS1(I),QS2(I), 1QS3(I) IF(MPRINT.GE.0)WRITE(LOU,106) 1NQP(I),NQS(I),QP1(I),QP2(I),QP3(I),QS1(I),QS2(I),QS3(I) 35 CONTINUE 34 IF(NRO.NE.0)GO TO 31 DO 32 IRHO=1,NC RHO1(IRHO)=1.7 RHO2(IRHO)=.2 32 CONTINUE 31 CONTINUE READ(LU,*)(PTOS(I),I=1,NC) IF(MPRINT.GE.0)WRITE(LOU,104)(PTOS(I),I=1,NC) DO 33 IRHO=1,NC 33 IF(PTOS(IRHO).LT.1.41421356)PTOS(IRHO)=1.732 IF(PTOS(1).GE.100.)RHO1(1)=1. IF(PTOS(1).GE.100.)RHO2(1)=0. NRHO=0 IF(PTOS(1).GE.100.)NRHO=1 C C INPUT DATA - PRINTER PLOT OF THE VELOCITY DISTRIBUTION C X COORDINATES OF VERTICAL BOUNDARIES OF THE MODEL C READ(LU,*)VMIN,VMAX,BMIN,BMAX,BLEFT,BRIGHT IF(MPRINT.GE.0)WRITE(LOU,104)VMIN,VMAX,BMIN,BMAX,BLEFT,BRIGHT MAUX=0 IF(ABS(BMIN).LT..00001)BMIN=A1(1,1) IF(ABS(BMAX).LT..00001)BMAX=A1(1,NINT) IF(ABS(BLEFT-BRIGHT).GT..001)GO TO 14 NC=NPNT(1) BLEFT=X1(1,1) BRIGHT=X1(NC,1) IF(MPRINT.EQ.0)RETURN C C NUMERICAL FORM OF INTERFACES C 14 AUX1=(BRIGHT-BLEFT)/25. IF(MPRINT.GT.0)WRITE(LOU,101)MTEXT IF(MPRINT.GE.2)WRITE(LOU,105)NINT IF(MPRINT.GE.2)WRITE(LOU,107) DO 19 I=1,NINT NC=NPNT(I)-1 IF(MPRINT.GE.2)WRITE(LOU,112) IF(MPRINT.GE.2)WRITE(LOU,108)I DO 15 J=1,NC IF(MPRINT.GE.2) 1WRITE(LOU,109)D1(J,I),C1(J,I),B1(J,I),A1(J,I),X1(J,I),X1(J+1,I), 2III(J,I) 15 CONTINUE IF(MPRINT.GE.2)WRITE(LOU,110) AUX2=BLEFT AUX3=AUX2 AUX(1)=AUX2 L=1 J=1 16 AUX3=AUX3+AUX1 IF(AUX3.GT.BRIGHT)GO TO 17 L=L+1 AUX(L)=AUX3 IF(L.LT.12)GO TO 16 17 IF(MPRINT.GE.2)WRITE(LOU,111)(AUX(M),M=1,L) K=1 18 IF(AUX(K).GT.X1(J+1,I))J=J+1 IF(AUX(K).GT.X1(J+1,I))GO TO 18 AUX4=AUX(K)-X1(J,I) AUX(K)=((D1(J,I)*AUX4+C1(J,I))*AUX4+B1(J,I))*AUX4+A1(J,I) K=K+1 IF(K.LE.L)GO TO 18 IF(MPRINT.GE.2)WRITE(LOU,111)(AUX(M),M=1,L) IF(MPRINT.GE.2)WRITE(LOU,112) IF((AUX3+AUX1).GT.BRIGHT)GO TO 19 L=0 GO TO 16 19 CONTINUE C C NUMERICAL FORM OF VELOCITY DISTRIBUTION C LAY=1 ITYPE=1 Y(3)=1. VMM=VMAX-VMIN VMMM=VMM/10. IF(MPRINT.GE.1)WRITE(LOU,114)VMIN,VMAX,VMMM DY=(BMAX-BMIN)/50. DX=(BRIGHT-BLEFT)/5. Y(2)=BMIN-DY K1=1 AUX(1)=BLEFT DO 20 I=2,6 20 AUX(I)=AUX(I-1)+DX IF(MPRINT.GE.1)WRITE(LOU,113)(AUX(I),I=1,6) DX=(BRIGHT-BLEFT)/100. NDER=0 DO 29 K=1,51 Y(2)=Y(2)+DY Y(1)=BLEFT-DX DO 28 L=1,101 Y(1)=Y(1)+DX IF(LAY.GE.NINT)GOTO 24 21 LAY1=LAY+1 NC=NPNT(LAY1)-1 DO 22 I=1,NC J=I IF(Y(1).LT.X1(I+1,LAY1))GO TO 23 22 CONTINUE 23 A2=A1(J,LAY1) B2=B1(J,LAY1) C2=C1(J,LAY1) D2=D1(J,LAY1) X2=X1(J,LAY1) AUX1=Y(1)-X2 ZINT=((D2*AUX1+C2)*AUX1+B2)*AUX1+A2 IF(Y(2).GE.ZINT)LAY=LAY+1 IF(LAY.GE.NINT)GO TO 27 IF(Y(2).GT.ZINT)GO TO 21 IF(LAY.LE.0)GO TO 27 24 NC=NPNT(LAY)-1 DO 25 I=1,NC J=I IF(Y(1).LT.X1(I+1,LAY))GO TO 26 25 CONTINUE 26 A2=A1(J,LAY) B2=B1(J,LAY) C2=C1(J,LAY) D2=D1(J,LAY) X2=X1(J,LAY) AUX1=Y(1)-X2 ZINT=((D2*AUX1+C2)*AUX1+B2)*AUX1+A2 IF(Y(2).LT.ZINT)LAY=LAY-1 IF(LAY.LE.0)GO TO 27 IF(Y(2).LT.ZINT)GO TO 24 27 IF(LAY.LE.0.OR.LAY.GE.NINT)IPR(L)=11 IF(LAY.LE.0.OR.LAY.GE.NINT)GO TO 28 CALL VELOC(Y,VEL) AUX1=10.*(VEL(1)-VMIN)/VMM IPR(L)=IFIX(AUX1) IF(AUX1.LT.0..OR.AUX1.GT.10.)IPR(L)=11 ITYPE=1 28 CONTINUE C C PRINTER PLOT OF THE VELOCITY DISTRIBUTION C IF(K1.EQ.K.AND.MPRINT.GE.1)WRITE(LOU,115)Y(2),IPR IF(K1.NE.K.AND.MPRINT.GE.1)WRITE(LOU,116)IPR IF(K1.EQ.K)K1=K1+10 29 CONTINUE C 101 FORMAT(A) 102 FORMAT(16I5) 103 FORMAT(3(2F10.5,I5)) 104 FORMAT(8F10.5) 105 FORMAT(////1X,'M O D E L O F M E D I U M'/2X,'NUMBER OF INTERF 1ACES - ',I3) 106 FORMAT(2I5,6F10.5) 107 FORMAT(///////1X,'I N T E R F A C E S'///1X,'INTERFACES ARE APPROX 1IMATED BY CUBIC PARABOLAS Z=D*(X-X1)**3+C*(X-X1)**2+B*(X-X1)+A B 2ETWEEN X1 AND X2'/1X,'COEFFICIENTS OF PARABOLAS ARE DETERMINED 3BY CUBIC SPLINE INTERPOLATION'///) 108 FORMAT(/1X,'COEFFICIENTS OF PARABOLAS APPROXIMATING INTERFACE',I3/ 1/15X,'D',14X,'C',14X,'B',14X,'A',2X,'IN INTERVAL',5X,'FROM X1', 27X,'TO X2',5X,'INDEX') 109 FORMAT(1X,4E15.5,15X,F10.5,F12.5,I10) 110 FORMAT(////1X,'NUMERICAL FORM OF INTERFACE'/1X,'UPPER ROW - X-COOR 1DINATES OF POINTS OF INTERFACE, LOWER ROW - CORRESPONDING Z-COORDI 2NATES OF POINTS OF INTERFACE'//) 111 FORMAT(1X,13F9.4) 112 FORMAT(/) 113 FORMAT(5X,5(F7.3,13X),F7.3) 114 FORMAT(////1X,'VELOCITY DISTRIBUTION'/1X,'ISOLINES CONSTRUCTED FRO 1M',F10.5,' TO ',F10.5,' WITH INCREMENT',F10.5//) 115 FORMAT(1X,F6.2,2X,101I1) 116 FORMAT(9X,101I1) 117 FORMAT(1X,'ISOVELOCITY DISCONTINUITIES') C RETURN END SUBROUTINE VELOC(Y,VEL) C C DETERMINATION OF VELOCITY AND DERIVATIVES FROM VELOCITY ISOLINES C SAVE AU,BU,CU,DU,XU,AL,BL,CL,DL,XL,JU,JL DIMENSION Y(4),VEL(6) COMMON/RAY/AY(12,400),DS(11,50),KINT(50),MREG,N,IREF,IND,IND1 COMMON/INTRF/A1(30,20),B1(29,20),C1(30,20),D1(29,20),X1(30,20), 1BRD(2),III(30,20),NPNT(20),NINT COMMON/AUXI/INTR,INT1,IOUT,KRE,IREFR,LAY,ITYPE,NDER,IPRINT,MPRINT, 1NTR,IMET COMMON/MEDIM/V1(19),V2(19),NVS(19),PTOS(19),MAUX COMMON/ABSR/QP1(19),QP2(19),QP3(19),QS1(19),QS2(19),QS3(19), 1NQP(19),NQS(19),NABS C JL=0 VU=V1(LAY) VL=V2(LAY) IA=LAY 1 NL=NPNT(IA)-1 DO 2 I=1,NL J=I IF(Y(1).LT.X1(I+1,IA))GO TO 3 2 CONTINUE 3 IF(IA.NE.LAY)GO TO 5 IF(MAUX.EQ.0)GO TO 8 IF(J.EQ.JU.AND.IA1.EQ.LAY)GO TO 11 8 JU=J AU=A1(J,IA) BU=B1(J,IA) CU=C1(J,IA) DU=D1(J,IA) XU=X1(J,IA) MAUX=1 11 AUX=Y(1)-XU AU1=DU*AUX AU2=3.*AU1+CU Z1=((AU1+CU)*AUX+BU)*AUX+AU DZ1=(AU2+CU)*AUX+BU DDZ1=2.*AU2 IA=IA+1 GO TO 1 5 IF(J.EQ.JL.AND.IA1.EQ.LAY)GO TO 16 JL=J AL=A1(J,IA) BL=B1(J,IA) CL=C1(J,IA) DL=D1(J,IA) XL=X1(J,IA) 16 AUX=Y(1)-XL AU1=DL*AUX AU2=3.*AU1+CL Z2=((AU1+CL)*AUX+BL)*AUX+AL DZ2=(AU2+CL)*AUX+BL DDZ2=2.*AU2 AUX=VL-VU IA1=LAY AUX1=Z2-Z1 AUX2=Y(2)-Z1 AUX3=Y(2)-Z2 AUX4=DZ2-DZ1 IF(ABS(AUX1).GT..000001)GO TO 7 13 VEL(1)=VU DO 6 I=2,6 6 VEL(I)=0. GO TO 4 7 AU1=AUX/AUX1 AU2=AU1/AUX1 AU3=DZ1*AUX3 AU4=DZ2*AUX2 VEL(1)=AU1*AUX2+VU VEL(2)=AU2*(AU3-AU4) VEL(3)=AU1 IF(NDER.EQ.0)GO TO 4 VEL(4)=(AU2/AUX1)*(2.*AUX4*(AU4-AU3)+AUX1*(DDZ1*AUX3-DDZ2*AUX2)) VEL(5)=-AU2*AUX4 VEL(6)=0. C 4 IF(ITYPE.GT.0.AND.NVS(LAY).EQ.0)GO TO 10 IF(ITYPE.LT.0.AND.NVS(LAY).EQ.1)GO TO 10 IF(ITYPE.GT.0.AND.NVS(LAY).EQ.1)GO TO 9 IF(PTOS(LAY).GE.100.)GO TO 12 VEL(1)=VEL(1)/PTOS(LAY) VEL(2)=VEL(2)/PTOS(LAY) VEL(3)=VEL(3)/PTOS(LAY) IF(NDER.EQ.0)GO TO 10 VEL(4)=VEL(4)/PTOS(LAY) VEL(5)=VEL(5)/PTOS(LAY) VEL(6)=VEL(6)/PTOS(LAY) GO TO 10 12 VEL(1)=0. VEL(2)=0. VEL(3)=0. IF(NDER.EQ.0)GO TO 10 VEL(4)=0. VEL(5)=0. VEL(6)=0. GO TO 10 9 VEL(1)=PTOS(LAY)*VEL(1) VEL(2)=PTOS(LAY)*VEL(2) VEL(3)=PTOS(LAY)*VEL(3) IF(NDER.EQ.0)GO TO 10 VEL(4)=PTOS(LAY)*VEL(4) VEL(5)=PTOS(LAY)*VEL(5) VEL(6)=PTOS(LAY)*VEL(6) 10 IF(Y(3).LT.100.)RETURN VV=VEL(1) AY(6,N)=VEL(1) AY(7,N)=VEL(2) AY(8,N)=VEL(3) IF(NDER.EQ.0)RETURN AY(9,N)=VEL(4) AY(10,N)=VEL(5) AY(11,N)=VEL(6) AY(12,N)=0. IF(NABS.EQ.0)RETURN IF(ITYPE.LT.0)GO TO 14 Q1=QP1(LAY) Q2=QP2(LAY) Q3=QP3(LAY) NQ=NQP(LAY) GO TO 15 14 Q1=QS1(LAY) Q2=QS2(LAY) Q3=QS3(LAY) NQ=NQS(LAY) 15 IF(NQ.EQ.0)AY(12,N)=Q1+Q2*VV+Q3*VV*VV IF(NQ.EQ.1)AY(12,N)=Q1+Q2*VV+Q3/VV RETURN END C C ************************************************************* C SUBROUTINE SPLIN(X,FX,NMIN,NMAX) C DIMENSION A(100),B(100),H(100),F(100),X(100),FX(100) C C CUBIC SPLINE INTERPOLATION WITH ZERO CURVATURES AT C THE END POINTS C IF((NMAX-NMIN).EQ.1)GO TO 4 NMIN1=NMIN+1 NMAX1=NMAX-1 H(NMIN1)=X(NMIN1)-X(NMIN) D2=(FX(NMIN1)-FX(NMIN))/H(NMIN1) DO 1 I=NMIN1,NMAX1 H(I+1)=X(I+1)-X(I) D1=D2 D2=(FX(I+1)-FX(I))/H(I+1) B(I)=H(I)+H(I+1) FX(I)=3.*(D2-D1)/B(I) A(I)=H(I)/B(I) 1 B(I)=H(I+1)/B(I) 4 FX(NMIN)=0. FX(NMAX)=0. IF((NMAX-NMIN).EQ.1)RETURN H(NMIN)=0. F(NMIN)=0. DO 2 I=NMIN1,NMAX1 XPOM=2.+A(I)*H(I-1) H(I)=-B(I)/XPOM 2 F(I)=(FX(I)-A(I)*F(I-1))/XPOM DO 3 I=NMIN,NMAX1 J=NMAX1-(I-NMIN) 3 FX(J)=H(J)*FX(J+1)+F(J) RETURN END C