SUBROUTINE MODEL (MTEXT,LU) C C ROUTINE CONTROLLING APPROXIMATION OF THE MEDIUM C DIMENSION A66U(6,6),A66L(6,6),ANGU(3),ANGL(3) 1 ,AUX(12),DEP(6),Y(2),IPR(101) CHARACTER*4 MTEXT(20) COMMON /AUXI/ IANI(20),INTR,INT1,IOUT,KRE,IREFR,LAY,NDER,IPRINT, 1 MPRINT,NTR,ISQRT,NAUX,ISOUR,MAUX,MREG,MDIM,IPOL,MSCON,LOU, 2 IAMP,MTRNS,ICOEF,IAD,IRHO,ISHEAR,IAC,IRT,mori COMMON /AUXX/ MMX(20),MMY(20),MMXY(20) COMMON /EPAR/ E66U(6,6,20),E66L(6,6,20) COMMON /DENS/ RHO(20) COMMON /INTRF/ Z(1000),SX(350),SY(350),NX(20),NY(20),BRD(4),NINT, 1 XINTA C READ(LU,101) MPRINT,NINT IF(MPRINT.GE.0)WRITE(LOU,101) MPRINT,NINT NLAY=NINT-1 C C INPUT FOR 3D INTERFACES C MX2=0 MY2=0 MXY2=0 DO 13 I=1,NINT MX1=MX2+1 MY1=MY2+1 MXY1=MXY2+1 READ(LU,101) MX,MY IF(MPRINT.GE.0) WRITE(LOU,101) MX,MY NX(I)=MX NY(I)=MY MX2=MX1+MX-1 MY2=MY1+MY-1 MXY2=MXY1+MX*MY-1 READ(LU,102)(SX(J),J=MX1,MX2) READ(LU,102)(SY(J),J=MY1,MY2) IF(MPRINT.GE.0) WRITE(LOU,102)(SX(J),J=MX1,MX2) IF(MPRINT.GE.0) WRITE(LOU,102)(SY(J),J=MY1,MY2) M1=MXY1 DO 12 L=1,MX M2=M1+MY-1 READ(LU,102)(Z(J),J=M1,M2) IF(MPRINT.GE.0) WRITE(LOU,102)(Z(J),J=M1,M2) 12 M1=M2+1 C C DETERMINATION OF COEFFICIENTS OF BICUBIC SPLINES C APPROXIMATING INTERFACES C CALL BIAP(MX1,MX,MY1,MY,MXY1) MMX(I)=MX1 MMY(I)=MY1 MMXY(I)=MXY1 13 CONTINUE NB1=MMX(1) NB2=MMX(2)-1 BRD(1)=SX(NB1) BRD(2)=SX(NB2) NB1=MMY(1) NB2=MMY(2)-1 BRD(3)=SY(NB1) BRD(4)=SY(NB2) C C INPUT DATA FOR PRINTER PLOT OF 3-D INTERFACES C DO 50 INTR=1,NINT IF(MPRINT.GE.1) WRITE(LOU,109) INTR READ(LU,102) ZMIN,ZMAX IF(MPRINT.GE.0) WRITE(LOU,102) ZMIN,ZMAX IF(MPRINT.EQ.0)GO TO 50 C C NUMERICAL FORM OF 3-D INTERFACES C ZMM=ZMAX-ZMIN ZMMM=ZMM/10. IF(MPRINT.GE.0) WRITE(LOU,102) ZMIN,ZMAX,ZMMM DY=(BRD(4)-BRD(3))/50. DX=(BRD(2)-BRD(1))/5. Y(2)=BRD(3)-DY K1=1 AUX(1)=BRD(1) DO 20 I=2,6 AUX(I)=AUX(I-1)+DX 20 CONTINUE IF(MPRINT.GE.1) WRITE(LOU,103)(AUX(I),I=1,6) DX=(BRD(2)-BRD(1))/100. MAUX=0 NDER=0 DO 29 K=1,51 Y(2)=Y(2)+DY Y(1)=BRD(1)-DX DO 28 L=1,101 Y(1)=Y(1)+DX CALL DISC(Y,DEP) AUX1=10.*(DEP(1)-ZMIN)/ZMM IPR(L)=IFIX(AUX1) IF(AUX1.LT.0.0.OR.AUX1.GT.10) IPR(L)=11 28 CONTINUE C C PRINTER PLOT OF 3-D INTERFACES C IF(K1.EQ.K.AND.MPRINT.GE.1) WRITE(LOU,104)Y(2),IPR IF(K1.NE.K.AND.MPRINT.GE.1) WRITE(LOU,105)IPR IF(K1.EQ.K)K1=K1+10 29 CONTINUE C C END OF LOOP OVER ALL INTERFACES C 50 CONTINUE C C INPUT OF ELASTIC PARAMETERS C READ(LU,101)ISQRT,IRHO IF(MPRINT.EQ.0)WRITE(LOU,101)ISQRT,IRHO isqrt=0 IF(ISQRT.EQ.0.AND.MPRINT.GT.0) WRITE(LOU,114) IF(ISQRT.NE.0.AND.MPRINT.GT.0) WRITE(LOU,113) if(irho.ne.0)then READ(LU,102)(RHO(L),L=1,NLAY) IF(MPRINT.GE.0)WRITE(LOU,102)(RHO(L),L=1,NLAY) DO 51 L=1,NLAY IF(ABS(RHO(L)).LT..00001)RHO(L)=1. 51 CONTINUE end if C C INPUT OF MATRIX OF ELASTIC PARAMETERS IN COMPRESSED FORM. C ELEMENTS OF THE MATRIX HAVE TO BE GIVEN IN (KM**2/SEC**2) C WHICH CORRESPONDS TO ELASTIC PARAMETERS DIVIDED BY DENSITY C C IF THE CRYSTAL IS GIVEN IN OTHER COORDINATE SYSTEM THAN C THE COORDINATE SYSTEM IN WHICH RAY TRACING IS PERFORMED, C ROUTINE TRFMAT MAKES THE CORRESPONDING TRANSFORMATION C DO 30 L=1,NLAY READ(LU,'(I10,6F10.4)')IANI(L),ANGU(1),ANGU(2),ANGU(3), 1ANGL(1),ANGL(2),ANGL(3) IF(MPRINT.EQ.0) 1WRITE(LOU,'(I10,6F10.4)')IANI(L),ANGU(1),ANGU(2),ANGU(3), 2ANGL(1),ANGL(2),ANGL(3) IROT1=1 IROT2=1 IF(ABS(ANGU(1)+ANGU(2)+ANGU(3)).LT.0.001) IROT1=0 IF(ABS(ANGL(1)+ANGL(2)+ANGL(3)).LT.0.001) IROT2=0 IF(IANI(L).NE.0) THEN READ(LU,111)((A66U(J,K),J=1,6),K=1,6) IF(MPRINT.EQ.0) 1 WRITE(LOU,111)((A66U(J,K),J=1,6),K=1,6) DO 55 K=1,6 DO 55 J=1,6 A66U(K,J)=A66U(J,K) 55 CONTINUE IF(MPRINT.GT.1) THEN WRITE(LOU,115) L WRITE(LOU,111)((A66U(J,K),J=1,6),K=1,6) END IF IF(IROT1.NE.0) THEN CALL TRFMAT(A66U,ANGU) IF(MPRINT.GE.1)WRITE(LOU,116)(ANGU(K),K=1,3) IF(MPRINT.GT.1) WRITE(LOU,111)((A66U(J,K),J=1,6),K=1,6) END IF READ(LU,111)((A66L(J,K),J=1,6),K=1,6) IF(MPRINT.EQ.0) 1 WRITE(LOU,111)((A66L(J,K),J=1,6),K=1,6) DO 60 K=1,6 DO 60 J=1,6 A66L(K,J)=A66L(J,K) 60 CONTINUE IF(MPRINT.GT.1) THEN WRITE(LOU,117) WRITE(LOU,111)((A66L(J,K),J=1,6),K=1,6) END IF IF(IROT2.NE.0) THEN CALL TRFMAT(A66L,ANGL) IF(MPRINT.GE.1)WRITE(LOU,116)(ANGL(K),K=1,3) IF(MPRINT.GT.1) WRITE(LOU,111)((A66L(J,K),J=1,6),K=1,6) END IF ELSE IF(MPRINT.GE.1)WRITE(LOU,118) L READ(LU,111) A66U(1,1),A66U(4,4) READ(LU,111) A66L(1,1),A66L(4,4) IF(MPRINT.EQ.0)WRITE(LOU,102) A66U(1,1),A66U(4,4) IF(MPRINT.EQ.0)WRITE(LOU,102) A66L(1,1),A66L(4,4) IF(MPRINT.GT.1) 1 WRITE(LOU,119)A66U(1,1),A66U(4,4),A66L(1,1),A66L(4,4) END IF DO 40 J=1,6 DO 40 K=1,6 E66U(K,J,L)=A66U(K,J) E66L(K,J,L)=A66L(K,J) 40 CONTINUE 30 CONTINUE WRITE(LOU,107)MTEXT C C FORMATS C 101 FORMAT(14I5) 102 FORMAT(8F10.5) 103 FORMAT(5X,5(F7.3,13X),F7.3) 104 FORMAT(F7.3,1X,101I1) 105 FORMAT(8X,101I1) 107 FORMAT(////,3X,20A///) 109 FORMAT(///,' INTERFACE NUMBER ',I5) 111 FORMAT(6F10.5) 113 FORMAT(//' INTERPOLATION IN SQUARE ROOTS OF ELASTIC PARAMETERES'/) 114 FORMAT(//' INTERPOLATION IN VALUES OF ELASTIC PARAMETERES'/) 115 FORMAT(/' LAYER',I4,' IS ANISOTROPIC ',/,' MATRIX OF 1 ELASTIC PARAMETERS GIVEN IN (KM/SEC)**2: 2 IT CONTAINS ELASTIC PARAMETERS/DENSITY', 3/,' MATRIX IS SPECIFIED IMMEDIATELY BELOW THE UPPER 4 BOUNDARY OF THE LAYER'/) 116 FORMAT(/' ROTATION APPLIED AROUND X1 WITH FI1=',F10.5,/,18X,'ARO 2UND X2 WITH FI2=',F10.5,/,18X,'AROUND X3 WITH FI3=',F10.5,/) 117 FORMAT(/' MATRIX OF ELASTIC PARAMETERS GIVEN IN (KM/SEC)**2: 1IT CONTAINS ELASTIC PARAMETERS/DENSITY', 2/,' MATRIX IS SPECIFIED IMMEDIATELY 3 ABOVE THE LOWER BOUNDARY OF THE LAYER'/) 118 FORMAT(/' LAYER',I4,' IS ISOTROPIC'/) 119 FORMAT(' VP**2=',F10.5,' (KM/SEC)**2 VS=',F10.5,' (KM/SEC)**2 1 IMMEDIATELY BELOW THE UPPER BOUNDARY OF THE LAYER,',/, 2' VP**2=',F10.5,' (KM/SEC)**2 VS=',F10.5,' (KM/SEC)*' 3,'*2 IMMEDIATELY ABOVE THE LOWER BOUNDARY OF THE LAYER') RETURN END