SUBROUTINE SCATT C C MODIFIED 26-MAR-1986 TO N-Kirchhoff scattering for free surface C MODIFIED 4-DEC-1986 TO N-Kirchhoff scattering for fluid/fluid. C MODIFIED 13-JAN-1987 TO N-Kirchhoff scattering for fluid/SOLID. C MODIFIED 7-MAY-1987 TO N-Kirchhoff scattering for all interfaces c PARAMETER (SQ2PI=2.506628,ONEO2PI=1.591549E-1) INCLUDE 'compar.f' INCLUDE 'comnla.f' INCLUDE 'comnp.f' C C LOCAL INTERMEDIATE MATRICES DEFINED C COMPLEX CC COMPLEX PCO(16),PRHS(32),PCOIN(16),PQ(32),PRES(32),BRES(32) COMPLEX DRES(32),RRES(32,isrowmax),RDRES(32) COMPLEX SCATI1(16),SCATI2(16),SCATI3(16),SCATI4(16) COMPLEX RPCO(16),RBRES(16) COMPLEX ALF0,ALF1,BET0,BET1 COMPLEX CETA,ALF,BET,CSUM,KSI,TSSKS,SSP LOGICAL NONKIRCH COMPLEX DIFK C C STATEMENT FUNCTIONS FOR N-K SCATTERING C c P(CETA,IND)=SQ2PI*CLEN(IND)*EXP(-0.5*CLEN(IND)*CLEN(IND) c 1 *CETA*CETA) ALF(CETA,IND)=SQRT(CETA*CETA-AK2(IND,1)) BET(CETA,IND)=SQRT(CETA*CETA-AK2(IND,2)) TSSKS(CETA,IND)=CON1(IND)*(2E0*CETA*CETA-AK2(IND,2)) DO 10 IN=1,NUMI IN1=IN+1 IF (ROUGH2(IN1).LT.1E-10) GO TO 10 SCCOIM = DLWNK(IN1) IF (V(IN,2).GT.1E-10) GO TO 100 C C VACUUM UPPER HALFSPACE C ********************** IF (V(IN1,3).GT.1E-10) GO TO 50 C NEXT LAYER OR HALFSPACE IS LIQUID C --------------------------------- C SCATTERING INTEGRALS IF (IMX(IN1).GE.1) THEN CSUM=CMPLX(0E0,0E0) C IMX(IN1)=5E0/(CLEN(IN1)*DLWNK(IN1)) DO 30 IK=-IMX(IN1),IMX(IN1) DDR=IK*DLWNK(IN1) DDI=SIGN(SCCOIM,REAL(WVNO)+DDR) DIFK=CMPLX(DDR,DDI) KSI=WVNO+DIFK 30 CSUM=CSUM+P(DIFK,IN1)*ALF(KSI,IN1) CSUM=ONEO2PI*ROUGH2(IN1)*DLWNK(IN1)*ALFA(IN1)*CSUM ELSE CSUM=ROUGH2(IN1)*ALFA(IN1)*ALFA(IN1) END IF C CORRECT LOCAL COEFFICIENT MATRICES AUP(IN1,3,1)=AUP(IN1,3,1) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN1)*ALFA(IN1)-CSUM) AUP(IN1,3,3)=AUP(IN1,3,3) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN1)*ALFA(IN1)+CSUM) C SOURCE TERMS SIMMILARLY (ONLY PROP TOWARDS INTERFACE) do is=1,nsrow RUIN(IN,3,IS)=RUIN(IN,3,IS) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN1)*ALFA(IN1)+CSUM) end do GO TO 10 50 CONTINUE C NEXT LAYER OR HALFSPACE IS SOLID C -------------------------------- C SA1=0.5+SIGN(0.5,REAL(AK2(IN1,1))-REAL(S2)) C SB1=0.5+SIGN(0.5,REAL(AK2(IN1,2))-REAL(S2)) SA1=1E0 SB1=1E0 NPEQ=2 C CLEAR SCATTERING TERMS CALL VCLR(PRES,1,16) CALL VCLR(DRES,1,16) do is=1,nsrow CALL VCLR(RRES(1,IS),1,8) end do CALL VCLR(RDRES,1,8) IF (IMX(IN1).GE.1) THEN NONKIRCH=.TRUE. CSUM=CMPLX(0E0,0E0) CALL VCLR(SCATI1,1,8) CALL VCLR(SCATI2,1,8) CALL VCLR(SCATI3,1,8) CALL VCLR(SCATI4,1,8) DO 501 IK=-IMX(IN1),IMX(IN1) DDR=IK*DLWNK(IN1) DDI=SIGN(SCCOIM,REAL(WVNO)+DDR) DIFK=CMPLX(DDR,DDI) KSI=WVNO+DIFK C VERTICAL WAVENUMBERS ALF1=ALF(KSI,IN1) BET1=BET(KSI,IN1) PCO(1)=-TSSKS(KSI,IN1) PCO(2)=-CON1(IN1)*2E0*KSI*ALF1 PCO(3)=CON1(IN1)*2E0*KSI*BET1 PCO(4)=-PCO(1) C COEFFICIENT MATRIX FOR ROTATION RPCO(1)=-2E0*AI*PCO(2) RPCO(2)=2E0*(-AI)*CON1(IN1)* & (AK2(IN1,1)-2E0*KSI*KSI) RPCO(3)=-2E0*AI*PCO(4) RPCO(4)=2E0*(-AI)*PCO(3) C C INVERT COEFFICIENT MATRIX C CALL CMATIN(NPEQ,PCO,PCOIN) C DERIVATIVES OF P AND Q (MULT BY VERT WAVENO) DO 502 JJ=1,NPEQ PQ(NPEQ*(JJ-1)+1)=PCOIN(NPEQ*(JJ-1)+1)*(-ALF1) PQ(NPEQ*(JJ-1)+2)=PCOIN(NPEQ*(JJ-1)+2)*(-BET1) 502 CONTINUE C SCATTERING INTEGRALS CALL CMMUL(PCO,NPEQ,NPEQ,PQ,NPEQ,BRES) CALL CMMUL(RPCO,NPEQ,NPEQ,PCOIN,NPEQ,RBRES) SSP=P(DIFK,IN1) DO 501 JJ=1,4 SCATI4(JJ)=SCATI4(JJ)+AI*DIFK*SSP*RBRES(JJ) SCATI3(JJ)=SCATI3(JJ)-AI*DIFK*SSP*BRES(JJ) SCATI2(JJ)=SCATI2(JJ)+DIFK*DIFK*SSP*RBRES(JJ) 501 SCATI1(JJ)=SCATI1(JJ)+SSP*BRES(JJ) DO 503 JJ=1,4 SCATI4(JJ)=ONEO2PI*DLWNK(IN1)*SCATI4(JJ) SCATI3(JJ)=ONEO2PI*DLWNK(IN1)*SCATI3(JJ) SCATI2(JJ)=ONEO2PI*DLWNK(IN1)*SCATI2(JJ) 503 SCATI1(JJ)=ONEO2PI*DLWNK(IN1)*SCATI1(JJ) ELSE NONKIRCH=.FALSE. C KIRCHHOFF CASE PCO(1)=AUP(IN1,3,1) PCO(2)=AUP(IN1,4,1) PCO(3)=AUP(IN1,3,2) PCO(4)=AUP(IN1,4,2) C C INVERT COEFFICIENT MATRIX C CALL CMATIN(NPEQ,PCO,PCOIN) C DERIVATIVES OF P AND Q (MULT BY VERT WAVENO) DO 53 JJ=1,NPEQ PQ((JJ-1)*NPEQ+1)=PCOIN((JJ-1)*NPEQ+1)*(-ALFA(IN1)) PQ(JJ*NPEQ)=PCOIN(JJ*NPEQ)*(-BETA(IN1)) 53 CONTINUE C CORRECTIONS FOR FIELD PARAMETERS CALL CMMUL(PCO,NPEQ,NPEQ,PQ,NPEQ,SCATI1) END IF C HOMEGENEOUS PART NRHS=4 DO 52 JJ=1,2 PRHS(JJ)=-ROUGH2(IN1)*(-ALFA(IN1))*AUP(IN1,JJ+2,1) PRHS(JJ+2)=-ROUGH2(IN1)*(-BETA(IN1))*AUP(IN1,JJ+2,2) PRHS(JJ+4)=-ROUGH2(IN1)*SA1*ALFA(IN1)*AUP(IN1,JJ+2,3) PRHS(JJ+6)=-ROUGH2(IN1)*SB1*BETA(IN1)*AUP(IN1,JJ+2,4) 52 CONTINUE C MULTIPLY WITH INVERTED MATRIX TO FIND P,Q CALL CMMUL(SCATI1,NPEQ,NPEQ,PRHS,NRHS,PRES) C IN KIRCHHOFF CASE ADD OTHER TERM IF (NONKIRCH) THEN CALL CMMUL(SCATI4,NPEQ,NPEQ,PRHS,NRHS,DRES) CALL VADD(PRES,1,DRES,1,PRES,1,2*NRHS*NPEQ) END IF C SOURCE TERMS SIMMILARLY (ONLY PROP TOWARDS INTERFACE) NRHS=1 do is=1,nsrow DO 56 JJ=1,NPEQ PRHS(JJ)=-ROUGH2(IN1)*SA1*ALFA(IN1)*RUIN(IN,JJ+2,is) 56 CONTINUE CALL CMMUL(SCATI1,NPEQ,NPEQ,PRHS,NRHS,RRES(1,IS)) C IN KIRCHHOFF CASE ADD OTHER TERM IF (NONKIRCH) THEN CALL CMMUL(SCATI4,NPEQ,NPEQ,PRHS,NRHS,RDRES) CALL VADD(RRES(1,IS),1,RDRES,1,RRES(1,IS),1,2*NRHS*NPEQ) END IF end do IF (NONKIRCH) THEN c c rotation terms only in n-k case c CALL VCLR(PRHS,1,16) NRHS=4 C rotation correction for normal stress PRHS(1)=-ROUGH2(IN1)*(-2*AI*AUP(IN1,4,1)) PRHS(3)=-ROUGH2(IN1)*(-2*AI*AUP(IN1,4,2)) PRHS(5)=-ROUGH2(IN1)*(-2*AI*AUP(IN1,4,3)) PRHS(7)=-ROUGH2(IN1)*(-2*AI*AUP(IN1,4,4)) C rotation correction for shear stress PRHS(2)=-ROUGH2(IN1)*2E0*(-AI)*CON1(IN1)* & (-ALFA(IN1)*AUP(IN1,1,1)-WVNO*AUP(IN1,2,1)) PRHS(4)=-ROUGH2(IN1)*2E0*(-AI)*CON1(IN1)* & (-BETA(IN1)*AUP(IN1,1,2)-WVNO*AUP(IN1,2,2)) PRHS(6)=-ROUGH2(IN1)*2E0*(-AI)*CON1(IN1)* & (ALFA(IN1)*AUP(IN1,1,3)-WVNO*AUP(IN1,2,3)) PRHS(8)=-ROUGH2(IN1)*2E0*(-AI)*CON1(IN1)* & (BETA(IN1)*AUP(IN1,1,4)-WVNO*AUP(IN1,2,4)) CALL CMMUL(SCATI2,NPEQ,NPEQ,PRHS,NRHS,DRES) CALL VADD(PRES,1,DRES,1,PRES,1,2*NPEQ*NRHS) CALL CMMUL(SCATI3,NPEQ,NPEQ,PRHS,NRHS,DRES) CALL VADD(PRES,1,DRES,1,PRES,1,2*NPEQ*NRHS) C SOURCE TERMS CALL VCLR(PRHS,1,12) NRHS=1 do is=1,nsrow PRHS(1)=-ROUGH2(IN1)*(-2*AI*RUIN(IN,4,is)) PRHS(2)=-ROUGH2(IN1)*2*(-AI)*CON1(IN1)* & (ALFA(IN1)*RUIN(IN,1,is)-WVNO*RUIN(IN,2,is)) CALL CMMUL(SCATI2,NPEQ,NPEQ,PRHS,NRHS,RDRES) CALL VADD(RRES(1,IS),1,RDRES,1,RRES(1,IS),1,2*NPEQ*NRHS) CALL CMMUL(SCATI3,NPEQ,NPEQ,PRHS,NRHS,RDRES) CALL VADD(RRES(1,IS),1,RDRES,1,RRES(1,IS),1,2*NPEQ*NRHS) end do END IF C CORRECT LOCAL COEFFICIENT MATRICES C DO 54 JJ=1,NRHS-1,2 DO 55 II=3,4 JJ=1 AUP(IN1,II,JJ)=AUP(IN1,II,JJ) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN1)*ALFA(IN1)) 2 +PRES(II-2+(JJ-1)*NPEQ) AUP(IN1,II,JJ+1)=AUP(IN1,II,JJ+1) 1 *(1E0+0.5*ROUGH2(IN1)*BETA(IN1)*BETA(IN1)) 2 +PRES(II-2+JJ*NPEQ) JJ=3 AUP(IN1,II,JJ)=AUP(IN1,II,JJ) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN1)*ALFA(IN1)) 2 +PRES(II-2+(JJ-1)*NPEQ) AUP(IN1,II,JJ+1)=AUP(IN1,II,JJ+1) 1 *(1E0+0.5*ROUGH2(IN1)*BETA(IN1)*BETA(IN1)) 2 +PRES(II-2+JJ*NPEQ) 55 CONTINUE 54 CONTINUE do is=1,nsrow DO 57 JJ=3,4 RUIN(IN,JJ,is)=RUIN(IN,JJ,is) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN1)*ALFA(IN1)) 1 +RRES(JJ-2,is) 57 CONTINUE end do GO TO 10 100 IF (V(IN,3).GT.1E-10) GO TO 200 C C LAYER OR HALFSPACE OVER INTERFACE IS LIQUID C ******************************************* IF (V(IN1,2).GT.1E-10) GO TO 125 C LOWER HALFSPACE IS VACUUM C ------------------------- C CORRECT LOCAL COEFFICIENT MATRICES C SCATTERING INTEGRALS IF (IMX(IN1).GE.1) THEN CSUM=CMPLX(0E0,0E0) C IMX(IN1)=5E0/(CLEN(IN1)*DLWNK(IN1)) DO 110 IK=-IMX(IN1),IMX(IN1) DDR=IK*DLWNK(IN1) DDI=SIGN(SCCOIM,REAL(WVNO)+DDR) DIFK=CMPLX(DDR,DDI) KSI=WVNO+DIFK 110 CSUM=CSUM+P(DIFK,IN1)*ALF(KSI,IN) CSUM=ONEO2PI*ROUGH2(IN1)*DLWNK(IN1)*ALFA(IN)*CSUM ELSE CSUM=ROUGH2(IN1)*ALFA(IN)*ALFA(IN) END IF ALO(IN,3,1)=ALO(IN,3,1) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN)*ALFA(IN)+CSUM) ALO(IN,3,3)=ALO(IN,3,3) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN)*ALFA(IN)-CSUM) C SOURCE TERMS SIMMILARLY (ONLY PROP TOWARDS INTERFACE) do is=1,nsrow ROIN(IN,3,is)=ROIN(IN,3,is) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN)*ALFA(IN)+CSUM) end do GO TO 10 125 IF (V(IN1,3).GT.1E-10) GO TO 150 C NEXT LAYER OR HALFSPACE IS LIQUID C --------------------------------- NPEQ=2 C SCATTERING INTEGRALS C SA=0.5+SIGN(0.5,REAL(AK2(IN,1))-REAL(S2)) C SA1=0.5+SIGN(0.5,REAL(AK2(IN1,1))-REAL(S2)) SA=1E0 SA1=1E0 CALL VCLR(PRES,1,16) CALL VCLR(DRES,1,16) do is=1,nsrow CALL VCLR(RRES(1,IS),1,8) end do CALL VCLR(RDRES,1,8) IF (IMX(IN1).GE.1) THEN NONKIRCH=.TRUE. CSUM=CMPLX(0E0,0E0) C IMX(IN1)=5E0/(CLEN(IN1)*DLWNK(IN1)) CALL VCLR(SCATI1,1,8) CALL VCLR(SCATI2,1,8) CALL VCLR(SCATI3,1,8) CALL VCLR(SCATI4,1,8) DO 130 IK=-IMX(IN1),IMX(IN1) DDR=IK*DLWNK(IN1) DDI=SIGN(SCCOIM,REAL(WVNO)+DDR) DIFK=CMPLX(DDR,DDI) KSI=WVNO+DIFK PCO(1)=ALF(KSI,IN) PCO(2)=-RCON1(IN) PCO(3)=ALF(KSI,IN1) PCO(4)=RCON1(IN1) C COEFFICIENT MATRIX FOR ROTATION RPCO(1)=-(-AI*KSI) RPCO(2)=0 RPCO(3)=-AI*KSI RPCO(4)=0 C C INVERT COEFFICIENT MATRIX C CALL CMATIN(NPEQ,PCO,PCOIN) C DERIVATIVES OF P AND Q (MULT BY VERT WAVENO) DO 1302 JJ=1,NPEQ PQ(NPEQ*JJ-1)=PCOIN(NPEQ*JJ-1)*PCO(1) PQ(NPEQ*JJ)=PCOIN(NPEQ*JJ)*(-PCO(3)) 1302 CONTINUE C CORRECTIONS FOR FIELD PARAMETERS CALL CMMUL(PCO,NPEQ,NPEQ,PQ,NPEQ,BRES) CALL CMMUL(RPCO,NPEQ,NPEQ,PCOIN,NPEQ,RBRES) SSP=P(DIFK,IN1) DO 130 JJ=1,4 SCATI4(JJ)=SCATI4(JJ)+AI*DIFK*SSP*RBRES(JJ) SCATI3(JJ)=SCATI3(JJ)-AI*DIFK*SSP*BRES(JJ) SCATI2(JJ)=SCATI2(JJ)+DIFK*DIFK*SSP*RBRES(JJ) 130 SCATI1(JJ)=SCATI1(JJ)+SSP*BRES(JJ) DO 1303 JJ=1,4 SCATI4(JJ)=ONEO2PI*DLWNK(IN1)*SCATI4(JJ) SCATI3(JJ)=ONEO2PI*DLWNK(IN1)*SCATI3(JJ) SCATI2(JJ)=ONEO2PI*DLWNK(IN1)*SCATI2(JJ) 1303 SCATI1(JJ)=ONEO2PI*DLWNK(IN1)*SCATI1(JJ) ELSE NONKIRCH=.FALSE. PCO(1)=ALFA(IN) PCO(2)=-RCON1(IN) PCO(3)=AUP(IN1,1,1) PCO(4)=AUP(IN1,3,1) C C INVERT COEFFICIENT MATRIX C CALL CMATIN(NPEQ,PCO,PCOIN) DO 1252 JJ=1,NPEQ PQ(NPEQ*JJ-1)=PCOIN(NPEQ*JJ-1)*ALFA(IN) PQ(NPEQ*JJ)=PCOIN(NPEQ*JJ)*(-ALFA(IN1)) 1252 CONTINUE C SCATTERING INTEGRAL FOR KIRCHHOFF CASE CALL CMMUL(PCO,NPEQ,NPEQ,PQ,NPEQ,SCATI1) CALL VCLR(SCATI2,1,8) END IF C HOMEGENEOUS PART NRHS=4 DO 1251 JJ=1,2 PRHS(JJ)= -ROUGH2(IN1)*(-SA*ALFA(IN))*ALO(IN,JJ*2-1,1) PRHS(JJ+2)= -ROUGH2(IN1)*ALFA(IN)*ALO(IN,JJ*2-1,3) PRHS(JJ+4)= -ROUGH2(IN1)*(-ALFA(IN1))*AUP(IN1,JJ*2-1,1) PRHS(JJ+6)= -ROUGH2(IN1)*SA1*ALFA(IN1)*AUP(IN1,JJ*2-1,3) 1251 CONTINUE C CREATE SCATTERING TERMS CALL CMMUL(SCATI1,NPEQ,NPEQ,PRHS,NRHS,PRES) C IN KIRCHOFF CASE ADD OTHER TERM IF (NONKIRCH) THEN CALL CMMUL(SCATI4,NPEQ,NPEQ,PRHS,NRHS,DRES) CALL VADD(PRES,1,DRES,1,PRES,1,2*NPEQ*NRHS) END IF C SOURCE TERMS SIMMILARLY (ONLY PROP TOWARDS INTERFACE) NRHS=2 do is=1,nsrow DO 1254 JJ=1,2 PRHS(JJ)=-ROUGH2(IN1)*(-SA*ALFA(IN))*ROIN(IN,2*JJ-1,is) PRHS(JJ+2)=-ROUGH2(IN1)*SA1*ALFA(IN1)*RUIN(IN,2*JJ-1,is) 1254 CONTINUE CALL CMMUL(SCATI1,NPEQ,NPEQ,PRHS,NRHS,RRES(1,IS)) C IN KIRCHOFF CASE ADD OTHER TERM IF (NONKIRCH) THEN CALL CMMUL(SCATI4,NPEQ,NPEQ,PRHS,NRHS,RDRES) CALL VADD(RRES(1,IS),1,RDRES,1,RRES(1,IS),1,2*NPEQ*NRHS) END IF end do IF (NONKIRCH) THEN c c rotation terms only in n-k case c CALL VCLR(PRHS,1,16) NRHS=4 PRHS(1)=-ROUGH2(IN1)*(-AI*ALO(IN,2,1)) PRHS(3)=-ROUGH2(IN1)*(-AI*ALO(IN,2,3)) PRHS(5)=-ROUGH2(IN1)*(-AI*AUP(IN1,2,1)) PRHS(7)=-ROUGH2(IN1)*(-AI*AUP(IN1,2,3)) CALL CMMUL(SCATI2,NPEQ,NPEQ,PRHS,NRHS,DRES) CALL VADD(PRES,1,DRES,1,PRES,1,2*NPEQ*NRHS) CALL CMMUL(SCATI3,NPEQ,NPEQ,PRHS,NRHS,DRES) CALL VADD(PRES,1,DRES,1,PRES,1,2*NPEQ*NRHS) C SOURCE TERMS CALL VCLR(PRHS,1,8) NRHS=2 do is=1,nsrow PRHS(1)=-ROUGH2(IN1)*(-AI*ROIN(IN,2,is)) PRHS(3)=-ROUGH2(IN1)*(-AI*RUIN(IN,2,is)) CALL CMMUL(SCATI2,NPEQ,NPEQ,PRHS,NRHS,RDRES) CALL VADD(RRES(1,IS),1,RDRES,1,RRES(1,IS),1,2*NPEQ*NRHS) CALL CMMUL(SCATI3,NPEQ,NPEQ,PRHS,NRHS,RDRES) CALL VADD(RRES(1,IS),1,RDRES,1,RRES(1,IS),1,2*NPEQ*NRHS) end do END IF C CORRECT LOCAL COEFFICIENT MATRICES DO 1253 JJ=1,NPEQ ALO(IN,2*JJ-1,1)=ALO(IN,2*JJ-1,1) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN)*ALFA(IN)) 2 +PRES(JJ) ALO(IN,2*JJ-1,3)=ALO(IN,2*JJ-1,3) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN)*ALFA(IN)) 2 +PRES(JJ+2) AUP(IN1,2*JJ-1,1)=AUP(IN1,2*JJ-1,1) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN1)*ALFA(IN1)) 2 +PRES(JJ+4) AUP(IN1,2*JJ-1,3)=AUP(IN1,2*JJ-1,3) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN1)*ALFA(IN1)) 2 +PRES(JJ+6) 1253 CONTINUE do is=1,nsrow do JJ=1,2 ROIN(IN,2*JJ-1,is)=ROIN(IN,2*JJ-1,is) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN)*ALFA(IN)) 2 +RRES(JJ,is) RUIN(IN,2*JJ-1,is)=RUIN(IN,2*JJ-1,is) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN1)*ALFA(IN1)) 2 +RRES(JJ+2,is) end do end do GO TO 10 C C NEXT LAYER OR HALFSPACE IS SOLID C -------------------------------- 150 CONTINUE NPEQ=3 C SA=0.5+SIGN(0.5,REAL(AK2(IN,1))-REAL(S2)) C SA1=0.5+SIGN(0.5,REAL(AK2(IN1,1))-REAL(S2)) C SB1=0.5+SIGN(0.5,REAL(AK2(IN1,2))-REAL(S2)) SA=1E0 SA1=1E0 SB1=1E0 C CLEAR SCATTERING TERMS CALL VCLR(PRES,1,36) CALL VCLR(DRES,1,36) do is=1,nsrow CALL VCLR(RRES(1,IS),1,18) end do CALL VCLR(RDRES,1,18) IF (IMX(IN1).GE.1) THEN NONKIRCH=.TRUE. CSUM=CMPLX(0E0,0E0) C IMX(IN1)=5E0/(CLEN(IN1)*DLWNK(IN1)) CALL VCLR(SCATI1,1,18) CALL VCLR(SCATI2,1,18) CALL VCLR(SCATI3,1,18) CALL VCLR(SCATI4,1,18) DO 1501 IK=-IMX(IN1),IMX(IN1) DDR=IK*DLWNK(IN1) DDI=SIGN(SCCOIM,REAL(WVNO)+DDR) DIFK=CMPLX(DDR,DDI) KSI=WVNO+DIFK C VERTICAL WAVENUMBERS ALF0=ALF(KSI,IN) ALF1=ALF(KSI,IN1) BET1=BET(KSI,IN1) PCO(1)=ALF0 PCO(2)=-RCON1(IN) PCO(3)=0E0 PCO(4)=ALF1 PCO(5)=-TSSKS(KSI,IN1) PCO(6)=-CON1(IN1)*2E0*KSI*ALF1 PCO(7)=-KSI PCO(8)=CON1(IN1)*2E0*KSI*BET1 PCO(9)=-PCO(5) C COEFFICIENT MATRIX FOR ROTATION RPCO(1)=-(-AI*KSI) RPCO(2)=0 RPCO(3)=0 RPCO(4)=-AI*KSI RPCO(5)=-2E0*AI*PCO(6) RPCO(6)=2E0*(-AI)*CON1(IN1)* & (AK2(IN1,1)-2E0*KSI*KSI) RPCO(7)=AI*BET1 RPCO(8)=-2E0*AI*PCO(9) RPCO(9)=2E0*(-AI)*PCO(8) C C INVERT COEFFICIENT MATRIX C CALL CMATIN(NPEQ,PCO,PCOIN) C DERIVATIVES OF P AND Q (MULT BY VERT WAVENO) DO 1502 JJ=1,NPEQ PQ(NPEQ*(JJ-1)+1)=PCOIN(NPEQ*(JJ-1)+1)*ALF0 PQ(NPEQ*(JJ-1)+2)=PCOIN(NPEQ*(JJ-1)+2)*(-ALF1) PQ(NPEQ*(JJ-1)+3)=PCOIN(NPEQ*(JJ-1)+3)*(-BET1) 1502 CONTINUE C SCATTERING INTEGRALS CALL CMMUL(PCO,NPEQ,NPEQ,PQ,NPEQ,BRES) CALL CMMUL(RPCO,NPEQ,NPEQ,PCOIN,NPEQ,RBRES) SSP=P(DIFK,IN1) DO 1501 JJ=1,9 SCATI4(JJ)=SCATI4(JJ)+AI*DIFK*SSP*RBRES(JJ) SCATI3(JJ)=SCATI3(JJ)-AI*DIFK*SSP*BRES(JJ) SCATI2(JJ)=SCATI2(JJ)+DIFK*DIFK*SSP*RBRES(JJ) SCATI1(JJ)=SCATI1(JJ)+SSP*BRES(JJ) 1501 continue DO 1503 JJ=1,9 SCATI4(JJ)=ONEO2PI*DLWNK(IN1)*SCATI4(JJ) SCATI3(JJ)=ONEO2PI*DLWNK(IN1)*SCATI3(JJ) SCATI2(JJ)=ONEO2PI*DLWNK(IN1)*SCATI2(JJ) 1503 SCATI1(JJ)=ONEO2PI*DLWNK(IN1)*SCATI1(JJ) ELSE NONKIRCH=.FALSE. C KIRCHHOFF CASE PCO(1)=ALFA(IN) PCO(2)=-RCON1(IN) PCO(3)=CMPLX(0E0,0E0) PCO(4)=AUP(IN1,1,1) PCO(7)=AUP(IN1,1,2) DO 151 JJ=2,3 PCO(JJ+3)=AUP(IN1,JJ+1,1) PCO(JJ+6)=AUP(IN1,JJ+1,2) 151 CONTINUE C C INVERT COEFFICIENT MATRIX C CALL CMATIN(NPEQ,PCO,PCOIN) C DERIVATIVES OF P AND Q (MULT BY VERT WAVENO) DO 153 JJ=1,NPEQ PQ(NPEQ*(JJ-1)+1)=PCOIN(NPEQ*(JJ-1)+1)*ALFA(IN) PQ(NPEQ*(JJ-1)+2)=PCOIN(NPEQ*(JJ-1)+2)*(-ALFA(IN1)) PQ(NPEQ*(JJ-1)+3)=PCOIN(NPEQ*(JJ-1)+3)*(-BETA(IN1)) 153 CONTINUE C SCATTERING INTEGRALS FOR KIRCHHOFF CASE CALL CMMUL(PCO,NPEQ,NPEQ,PQ,NPEQ,SCATI1) CALL VCLR(SCATI2,1,18) END IF C HOMEGENEOUS PART NRHS=6 PRHS(1)= -ROUGH2(IN1)*(-SA*ALFA(IN))*ALO(IN,1,1) PRHS(4)= -ROUGH2(IN1)*ALFA(IN)*ALO(IN,1,3) PRHS(7)=-ROUGH2(IN1)*(-ALFA(IN1))*AUP(IN1,1,1) PRHS(10)=-ROUGH2(IN1)*(-BETA(IN1))*AUP(IN1,1,2) PRHS(13)=-ROUGH2(IN1)*SA1*ALFA(IN1)*AUP(IN1,1,3) PRHS(16)=-ROUGH2(IN1)*SB1*BETA(IN1)*AUP(IN1,1,4) DO 152 JJ=2,3 PRHS(JJ)= -ROUGH2(IN1)*(-SA*ALFA(IN))*ALO(IN,JJ+1,1) PRHS(JJ+3)= -ROUGH2(IN1)*ALFA(IN)*ALO(IN,JJ+1,3) PRHS(JJ+6)=-ROUGH2(IN1)*(-ALFA(IN1))*AUP(IN1,JJ+1,1) PRHS(JJ+9)=-ROUGH2(IN1)*(-BETA(IN1))*AUP(IN1,JJ+1,2) PRHS(JJ+12)=-ROUGH2(IN1)*SA1*ALFA(IN1)*AUP(IN1,JJ+1,3) PRHS(JJ+15)=-ROUGH2(IN1)*SB1*BETA(IN1)*AUP(IN1,JJ+1,4) 152 CONTINUE C CREATE SCATTERING TERMS CALL CMMUL(SCATI1,NPEQ,NPEQ,PRHS,NRHS,PRES) C IN KIRCHHOFF CASE ADD OTHER TERM IF (NONKIRCH) THEN CALL CMMUL(SCATI4,NPEQ,NPEQ,PRHS,NRHS,DRES) CALL VADD(PRES,1,DRES,1,PRES,1,2*NRHS*NPEQ) END IF C SOURCE TERMS SIMMILARLY (ONLY PROP TOWARDS INTERFACE) NRHS=2 do is=1,nsrow PRHS(1)=-ROUGH2(IN1)*(-SA*ALFA(IN))*ROIN(IN,1,is) PRHS(4)=-ROUGH2(IN1)*SA1*ALFA(IN1)*RUIN(IN,1,is) DO 157 JJ=2,3 PRHS(JJ)=-ROUGH2(IN1)*(-SA*ALFA(IN))*ROIN(IN,JJ+1,is) PRHS(JJ+NPEQ)=-ROUGH2(IN1)*SA1*ALFA(IN1)*RUIN(IN,JJ+1,is) 157 CONTINUE CALL CMMUL(SCATI1,NPEQ,NPEQ,PRHS,NRHS,RRES(1,IS)) C IN KIRCHHOFF CASE ADD OTHER TERM IF (NONKIRCH) THEN CALL CMMUL(SCATI4,NPEQ,NPEQ,PRHS,NRHS,RDRES) CALL VADD(RRES(1,IS),1,RDRES,1,RRES(1,IS),1,2*NRHS*NPEQ) END IF end do IF (NONKIRCH) THEN c c rotation terms only in n-k case c CALL VCLR(PRHS,1,36) NRHS=6 C rotation correction for vertical displacement PRHS(1)=-ROUGH2(IN1)*(-AI*ALO(IN,2,1)) PRHS(4)=-ROUGH2(IN1)*(-AI*ALO(IN,2,3)) PRHS(7)=-ROUGH2(IN1)*(-AI*AUP(IN1,2,1)) PRHS(10)=-ROUGH2(IN1)*(-AI*AUP(IN1,2,2)) PRHS(13)=-ROUGH2(IN1)*(-AI*AUP(IN1,2,3)) PRHS(16)=-ROUGH2(IN1)*(-AI*AUP(IN1,2,4)) C rotation correction for normal stress PRHS(2)=0 PRHS(5)=0 PRHS(8)=-ROUGH2(IN1)*(-2*AI*AUP(IN1,4,1)) PRHS(11)=-ROUGH2(IN1)*(-2*AI*AUP(IN1,4,2)) PRHS(14)=-ROUGH2(IN1)*(-2*AI*AUP(IN1,4,3)) PRHS(17)=-ROUGH2(IN1)*(-2*AI*AUP(IN1,4,4)) C rotation correction for shear stress PRHS(3)=0 PRHS(6)=0 C PRHS(9)=-ROUGH2(IN1)*2E0*(-AI)*CON1(IN1)* C & (-ALFA(IN1)*AUP(IN1,1,1)-(-AI*WVNO)*AI*AUP(IN1,2,1)) PRHS(9)=-ROUGH2(IN1)*2E0*(-AI)*CON1(IN1)* & (-ALFA(IN1)*AUP(IN1,1,1)-WVNO*AUP(IN1,2,1)) PRHS(12)=-ROUGH2(IN1)*2E0*(-AI)*CON1(IN1)* & (-BETA(IN1)*AUP(IN1,1,2)-WVNO*AUP(IN1,2,2)) PRHS(15)=-ROUGH2(IN1)*2E0*(-AI)*CON1(IN1)* & (ALFA(IN1)*AUP(IN1,1,3)-WVNO*AUP(IN1,2,3)) PRHS(18)=-ROUGH2(IN1)*2E0*(-AI)*CON1(IN1)* & (BETA(IN1)*AUP(IN1,1,4)-WVNO*AUP(IN1,2,4)) CALL CMMUL(SCATI2,NPEQ,NPEQ,PRHS,NRHS,DRES) CALL VADD(PRES,1,DRES,1,PRES,1,2*NPEQ*NRHS) CALL CMMUL(SCATI3,NPEQ,NPEQ,PRHS,NRHS,DRES) CALL VADD(PRES,1,DRES,1,PRES,1,2*NPEQ*NRHS) C SOURCE TERMS CALL VCLR(PRHS,1,12) NRHS=2 do is=1,nsrow PRHS(1)=-ROUGH2(IN1)*(-AI*ROIN(IN,2,is)) PRHS(4)=-ROUGH2(IN1)*(-AI*RUIN(IN,2,is)) PRHS(5)=-ROUGH2(IN1)*(-2*AI*RUIN(IN,4,is)) PRHS(6)=-ROUGH2(IN1)*2*(-AI)*CON1(IN1)* & (ALFA(IN1)*RUIN(IN,1,is)-WVNO*RUIN(IN,2,is)) CALL CMMUL(SCATI2,NPEQ,NPEQ,PRHS,NRHS,RDRES) CALL VADD(RRES(1,IS),1,RDRES,1,RRES(1,IS),1,2*NPEQ*NRHS) CALL CMMUL(SCATI3,NPEQ,NPEQ,PRHS,NRHS,RDRES) CALL VADD(RRES(1,IS),1,RDRES,1,RRES(1,IS),1,2*NPEQ*NRHS) end do END IF C CORRECT LOCAL COEFFICIENT MATRICES III=1 II=1 ALO(IN,1,II)=ALO(IN,1,II) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN)*ALFA(IN)) 2 +PRES(1+(III-1)*NPEQ) AUP(IN1,1,II)=AUP(IN1,1,II) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN1)*ALFA(IN1)) 2 +PRES(1+(II-1)*NPEQ+6) AUP(IN1,1,II+1)=AUP(IN1,1,II+1) 1 *(1E0+0.5*ROUGH2(IN1)*BETA(IN1)*BETA(IN1)) 2 +PRES(1+II*NPEQ+6) III=2 II=3 ALO(IN,1,II)=ALO(IN,1,II) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN)*ALFA(IN)) 2 +PRES(1+(III-1)*NPEQ) AUP(IN1,1,II)=AUP(IN1,1,II) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN1)*ALFA(IN1)) 2 +PRES(1+(II-1)*NPEQ+6) AUP(IN1,1,II+1)=AUP(IN1,1,II+1) 1 *(1E0+0.5*ROUGH2(IN1)*BETA(IN1)*BETA(IN1)) 2 +PRES(1+II*NPEQ+6) DO 155 JJ=2,3 III=1 II=1 ALO(IN,JJ+1,II)=ALO(IN,JJ+1,II) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN)*ALFA(IN)) 2 +PRES(JJ+(III-1)*NPEQ) AUP(IN1,JJ+1,II)=AUP(IN1,JJ+1,II) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN1)*ALFA(IN1)) 2 +PRES(JJ+(II-1)*NPEQ+6) AUP(IN1,JJ+1,II+1)=AUP(IN1,JJ+1,II+1) 1 *(1E0+0.5*ROUGH2(IN1)*BETA(IN1)*BETA(IN1)) 2 +PRES(JJ+II*NPEQ+6) III=2 II=3 ALO(IN,JJ+1,II)=ALO(IN,JJ+1,II) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN)*ALFA(IN)) 2 +PRES(JJ+(III-1)*NPEQ) AUP(IN1,JJ+1,II)=AUP(IN1,JJ+1,II) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN1)*ALFA(IN1)) 2 +PRES(JJ+(II-1)*NPEQ+6) AUP(IN1,JJ+1,II+1)=AUP(IN1,JJ+1,II+1) 1 *(1E0+0.5*ROUGH2(IN1)*BETA(IN1)*BETA(IN1)) 2 +PRES(JJ+II*NPEQ+6) 155 CONTINUE do is=1,nsrow ROIN(IN,1,is)=ROIN(IN,1,is) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN)*ALFA(IN)) 2 +RRES(1,is) RUIN(IN,1,is)=RUIN(IN,1,is) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN1)*ALFA(IN1)) 2 +RRES(4,is) DO 159 JJ=2,3 ROIN(IN,JJ+1,is)=ROIN(IN,JJ+1,is) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN)*ALFA(IN)) 2 +RRES(JJ,is) RUIN(IN,JJ+1,is)=RUIN(IN,JJ+1,is) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN1)*ALFA(IN1)) 2 +RRES(JJ+3,is) 159 CONTINUE end do GO TO 10 200 CONTINUE C LAYER OR HALFSPACE OVER INTERFACE IS SOLID C ******************************************* C LOWER HALFSPACE VACUUM C ---------------------- IF (V(IN1,2).GT.1E-10) GO TO 225 C SA=0.5+SIGN(0.5,REAL(AK2(IN,1))-REAL(S2)) C SB=0.5+SIGN(0.5,REAL(AK2(IN,2))-REAL(S2)) SA=1E0 SB=1E0 NPEQ=2 C CLEAR SCATTERING TERMS CALL VCLR(PRES,1,16) CALL VCLR(DRES,1,16) do is=1,nsrow CALL VCLR(RRES(1,IS),1,8) end do CALL VCLR(RDRES,1,8) C ICHECK=ICHECK+1 IF (IMX(IN1).GE.1) THEN NONKIRCH=.TRUE. CSUM=CMPLX(0E0,0E0) CALL VCLR(SCATI1,1,8) CALL VCLR(SCATI2,1,8) CALL VCLR(SCATI3,1,8) CALL VCLR(SCATI4,1,8) DO 2001 IK=-IMX(IN1),IMX(IN1) DDR=IK*DLWNK(IN1) DDI=SIGN(SCCOIM,REAL(WVNO)+DDR) DIFK=CMPLX(DDR,DDI) KSI=WVNO+DIFK C VERTICAL WAVENUMBERS ALF0=ALF(KSI,IN) BET0=BET(KSI,IN) PCO(1)=TSSKS(KSI,IN) PCO(2)=-CON1(IN)*2E0*KSI*ALF0 PCO(3)=CON1(IN)*2E0*KSI*BET0 PCO(4)=-PCO(1) C IF (ICHECK.EQ.10.AND.IK.EQ.0) THEN C WRITE(49,'(4E15.6)') (PCO(JJ),JJ=1,16) C END IF C COEFFICIENT MATRIX FOR ROTATION RPCO(1)=-2*AI*PCO(2) RPCO(2)=2E0*(-AI)*CON1(IN)* & (2E0*KSI*KSI-AK2(IN,1)) RPCO(3)=-2*AI*PCO(4) RPCO(4)=-2*AI*PCO(3) C C INVERT COEFFICIENT MATRIX C CALL CMATIN(NPEQ,PCO,PCOIN) C DERIVATIVES OF P AND Q (MULT BY VERT WAVENO) DO 2002 JJ=1,NPEQ PQ(NPEQ*(JJ-1)+1)=PCOIN(NPEQ*(JJ-1)+1)*ALF0 PQ(NPEQ*(JJ-1)+2)=PCOIN(NPEQ*(JJ-1)+2)*BET0 2002 CONTINUE C SCATTERING INTEGRALS CALL CMMUL(PCO,NPEQ,NPEQ,PQ,NPEQ,BRES) CALL CMMUL(RPCO,NPEQ,NPEQ,PCOIN,NPEQ,RBRES) SSP=P(DIFK,IN1) DO 2001 JJ=1,NPEQ*NPEQ SCATI4(JJ)=SCATI4(JJ)+AI*DIFK*SSP*RBRES(JJ) SCATI3(JJ)=SCATI3(JJ)-AI*DIFK*SSP*BRES(JJ) SCATI2(JJ)=SCATI2(JJ)+DIFK*DIFK*SSP*RBRES(JJ) 2001 SCATI1(JJ)=SCATI1(JJ)+SSP*BRES(JJ) DO 2003 JJ=1,NPEQ*NPEQ SCATI4(JJ)=ONEO2PI*DLWNK(IN1)*SCATI4(JJ) SCATI3(JJ)=ONEO2PI*DLWNK(IN1)*SCATI3(JJ) SCATI2(JJ)=ONEO2PI*DLWNK(IN1)*SCATI2(JJ) 2003 SCATI1(JJ)=ONEO2PI*DLWNK(IN1)*SCATI1(JJ) ELSE NONKIRCH=.FALSE. C KIRCHHOFF CASE PCO(1)=ALO(IN,3,3)*EZALFM(IN) PCO(2)=ALO(IN,4,3)*EZALFM(IN) PCO(3)=ALO(IN,3,4)*EZBETM(IN) PCO(4)=ALO(IN,4,4)*EZBETM(IN) C PCO(1)=CON2(IN) C PCO(2)=-CON3(IN) C PCO(3)=CON4(IN) C PCO(4)=-CON2(IN) C C INVERT COEFFICIENT MATRIX C CALL CMATIN(NPEQ,PCO,PCOIN) C DERIVATIVES OF P AND Q (MULT BY VERT WAVENO) DO 203 JJ=1,NPEQ PQ((JJ-1)*NPEQ+1)=PCOIN((JJ-1)*NPEQ+1)*ALFA(IN) PQ(JJ*NPEQ)=PCOIN(JJ*NPEQ)*BETA(IN) 203 CONTINUE C KIRCHHOFF SCATTERING INTEGRAL CALL CMMUL(PCO,NPEQ,NPEQ,PQ,NPEQ,SCATI1) CALL VCLR(SCATI2,1,8) END IF C HOMEGENEOUS PART NRHS=4 DO 202 JJ=1,2 PRHS(JJ)=-ROUGH2(IN1)*(-SA*ALFA(IN))*ALO(IN,JJ+2,1) PRHS(JJ+2)=-ROUGH2(IN1)*(-SB*BETA(IN))*ALO(IN,JJ+2,2) PRHS(JJ+4)=-ROUGH2(IN1)*ALFA(IN)*ALO(IN,JJ+2,3) PRHS(JJ+6)=-ROUGH2(IN1)*BETA(IN)*ALO(IN,JJ+2,4) 202 CONTINUE C MULTIPLY BY SCATTERING INTEGRAL CALL CMMUL(SCATI1,NPEQ,NPEQ,PRHS,NRHS,PRES) C IN NON-KIRCHHOFF CASE ADD OTHER TERM IF (NONKIRCH) THEN CALL CMMUL(SCATI4,NPEQ,NPEQ,PRHS,NRHS,DRES) CALL VADD(PRES,1,DRES,1,PRES,1,2*NRHS*NPEQ) END IF C SOURCE TERMS SIMMILARLY (ONLY PROP TOWARDS INTERFACE) NRHS=1 do is=1,nsrow DO 206 JJ=1,NPEQ PRHS(JJ)=-ROUGH2(IN1)*(-SA*ALFA(IN))*ROIN(IN,JJ+2,is) 206 CONTINUE CALL CMMUL(SCATI1,NPEQ,NPEQ,PRHS,NRHS,RRES(1,IS)) C IN KIRCHHOFF CASE ADD OTHER TERM IF (NONKIRCH) THEN CALL CMMUL(SCATI4,NPEQ,NPEQ,PRHS,NRHS,RDRES) CALL VADD(RRES(1,IS),1,RDRES,1,RRES(1,IS),1,2*NRHS*NPEQ) END IF end do IF (NONKIRCH) THEN c c rotation terms only in n-k case c CALL VCLR(PRHS,1,16) NRHS=4 C rotation correction for normal stress PRHS(1)=-ROUGH2(IN1)*(-2E0*AI*ALO(IN,4,1)) PRHS(3)=-ROUGH2(IN1)*(-2E0*AI*ALO(IN,4,2)) PRHS(5)=-ROUGH2(IN1)*(-2E0*AI*ALO(IN,4,3)) PRHS(7)=-ROUGH2(IN1)*(-2E0*AI*ALO(IN,4,4)) C rotation correction for shear stress PRHS(2)=-ROUGH2(IN1)*2E0*(-AI)*CON1(IN)* & (-ALFA(IN)*ALO(IN,1,1)-WVNO*ALO(IN,2,1)) PRHS(4)=-ROUGH2(IN1)*2E0*(-AI)*CON1(IN)* & (-BETA(IN)*ALO(IN,1,2)-WVNO*ALO(IN,2,2)) PRHS(6)=-ROUGH2(IN1)*2E0*(-AI)*CON1(IN)* & (ALFA(IN)*ALO(IN,1,3)-WVNO*ALO(IN,2,3)) PRHS(8)=-ROUGH2(IN1)*2E0*(-AI)*CON1(IN)* & (BETA(IN)*ALO(IN,1,4)-WVNO*ALO(IN,2,4)) CALL CMMUL(SCATI2,NPEQ,NPEQ,PRHS,NRHS,DRES) CALL VADD(PRES,1,DRES,1,PRES,1,2*NPEQ*NRHS) CALL CMMUL(SCATI3,NPEQ,NPEQ,PRHS,NRHS,DRES) CALL VADD(PRES,1,DRES,1,PRES,1,2*NPEQ*NRHS) C SOURCE TERMS CALL VCLR(PRHS,1,4) NRHS=1 do is=1,nsrow PRHS(1)=-ROUGH2(IN1)*(-2E0*AI*ROIN(IN,4,is)) PRHS(2)=-ROUGH2(IN1)*2*(-AI)*CON1(IN)* & (-ALFA(IN)*ROIN(IN,1,is)-WVNO*ROIN(IN,2,is)) CALL CMMUL(SCATI2,NPEQ,NPEQ,PRHS,NRHS,RDRES) CALL VADD(RRES(1,IS),1,RDRES,1,RRES(1,IS),1,2*NPEQ*NRHS) CALL CMMUL(SCATI3,NPEQ,NPEQ,PRHS,NRHS,RDRES) CALL VADD(RRES(1,IS),1,RDRES,1,RRES(1,IS),1,2*NPEQ*NRHS) end do END IF C CORRECT LOCAL COEFFICIENT MATRICES C DO 204 JJ=1,NRHS-1,2 DO 205 II=3,4 JJ=1 ALO(IN,II,JJ)=ALO(IN,II,JJ) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN)*ALFA(IN)) 2 +PRES(II-2+(JJ-1)*NPEQ) ALO(IN,II,JJ+1)=ALO(IN,II,JJ+1) 1 *(1E0+0.5*ROUGH2(IN1)*BETA(IN)*BETA(IN)) 2 +PRES(II-2+JJ*NPEQ) JJ=3 ALO(IN,II,JJ)=ALO(IN,II,JJ) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN)*ALFA(IN)) 2 +PRES(II-2+(JJ-1)*NPEQ) ALO(IN,II,JJ+1)=ALO(IN,II,JJ+1) 1 *(1E0+0.5*ROUGH2(IN1)*BETA(IN)*BETA(IN)) 2 +PRES(II-2+JJ*NPEQ) 205 CONTINUE do is=1,nsrow DO 207 JJ=3,4 ROIN(IN,JJ,is)=ROIN(IN,JJ,is) & *(1E0+0.5*ROUGH2(IN1)*ALFA(IN)*ALFA(IN))+RRES(JJ-2,is) 207 CONTINUE end do GO TO 10 225 IF (V(IN1,3).GT.1E-10) GO TO 250 C NEXT LAYER OR HALFSPACE IS LIQUID C --------------------------------- C SA=0.5+SIGN(0.5,REAL(AK2(IN,1))-REAL(S2)) C SB=0.5+SIGN(0.5,REAL(AK2(IN,2))-REAL(S2)) C SA1=0.5+SIGN(0.5,REAL(AK2(IN1,1))-REAL(S2)) SA=1E0 SB=1E0 SA1=1E0 NPEQ=3 C CLEAR SCATTERING TERMS CALL VCLR(PRES,1,36) CALL VCLR(DRES,1,36) do is=1,nsrow CALL VCLR(RRES(1,IS),1,18) end do CALL VCLR(RDRES,1,18) C ICHECK=ICHECK+1 IF (IMX(IN1).GE.1) THEN NONKIRCH=.TRUE. CSUM=CMPLX(0E0,0E0) CALL VCLR(SCATI1,1,18) CALL VCLR(SCATI2,1,18) CALL VCLR(SCATI3,1,18) CALL VCLR(SCATI4,1,18) DO 2201 IK=-IMX(IN1),IMX(IN1) DDR=IK*DLWNK(IN1) DDI=SIGN(SCCOIM,REAL(WVNO)+DDR) DIFK=CMPLX(DDR,DDI) KSI=WVNO+DIFK C VERTICAL WAVENUMBERS ALF0=ALF(KSI,IN) BET0=BET(KSI,IN) ALF1=ALF(KSI,IN1) PCO(1)=ALF0 PCO(2)=TSSKS(KSI,IN) PCO(3)=-CON1(IN)*2E0*KSI*ALF0 PCO(4)=KSI PCO(5)=CON1(IN)*2E0*KSI*BET0 PCO(6)=-PCO(2) PCO(7)=ALF1 PCO(8)=RCON1(IN1) PCO(9)=0 C IF (ICHECK.EQ.10.AND.IK.EQ.0) THEN C WRITE(49,'(4E15.6)') (PCO(JJ),JJ=1,9) C END IF C COEFFICIENT MATRIX FOR ROTATION RPCO(1)=-AI*(-KSI) RPCO(2)=-2*AI*PCO(3) RPCO(3)=2E0*(-AI)*CON1(IN)* & (2E0*KSI*KSI-AK2(IN,1)) RPCO(4)=-AI*(-BET0) RPCO(5)=-2*AI*PCO(6) RPCO(6)=-2*AI*PCO(5) RPCO(7)=-AI*KSI RPCO(8)=0 RPCO(9)=0 C C INVERT COEFFICIENT MATRIX C CALL CMATIN(NPEQ,PCO,PCOIN) C DERIVATIVES OF P AND Q (MULT BY VERT WAVENO) DO 2202 JJ=1,NPEQ PQ(NPEQ*(JJ-1)+1)=PCOIN(NPEQ*(JJ-1)+1)*ALF0 PQ(NPEQ*(JJ-1)+2)=PCOIN(NPEQ*(JJ-1)+2)*BET0 PQ(NPEQ*(JJ-1)+3)=PCOIN(NPEQ*(JJ-1)+3)*(-ALF1) 2202 CONTINUE C SCATTERING INTEGRALS CALL CMMUL(PCO,NPEQ,NPEQ,PQ,NPEQ,BRES) CALL CMMUL(RPCO,NPEQ,NPEQ,PCOIN,NPEQ,RBRES) SSP=P(DIFK,IN1) DO 2201 JJ=1,9 SCATI4(JJ)=SCATI4(JJ)+AI*DIFK*SSP*RBRES(JJ) SCATI3(JJ)=SCATI3(JJ)-AI*DIFK*SSP*BRES(JJ) SCATI2(JJ)=SCATI2(JJ)+DIFK*DIFK*SSP*RBRES(JJ) 2201 SCATI1(JJ)=SCATI1(JJ)+SSP*BRES(JJ) DO 2203 JJ=1,9 SCATI4(JJ)=ONEO2PI*DLWNK(IN1)*SCATI4(JJ) SCATI3(JJ)=ONEO2PI*DLWNK(IN1)*SCATI3(JJ) SCATI2(JJ)=ONEO2PI*DLWNK(IN1)*SCATI2(JJ) 2203 SCATI1(JJ)=ONEO2PI*DLWNK(IN1)*SCATI1(JJ) ELSE NONKIRCH=.FALSE. C KIRCHHOFF CASE C PCO(1)=ALFA(IN) C PCO(2)=CON2(IN) C PCO(3)=-CON3(IN) C PCO(4)=WVNO C PCO(5)=CON4(IN) C PCO(6)=-CON2(IN) PCO(1)=ALO(IN,1,3)*EZALFM(IN) PCO(2)=ALO(IN,3,3)*EZALFM(IN) PCO(3)=ALO(IN,4,3)*EZALFM(IN) PCO(4)=ALO(IN,1,4)*EZBETM(IN) PCO(5)=ALO(IN,3,4)*EZBETM(IN) PCO(6)=ALO(IN,4,4)*EZBETM(IN) PCO(7)=AUP(IN1,1,1) PCO(8)=AUP(IN1,3,1) PCO(9)=AUP(IN1,4,1) C C INVERT COEFFICIENT MATRIX C CALL CMATIN(NPEQ,PCO,PCOIN) C DERIVATIVES OF P AND Q (MULT BY VERT WAVENO) DO 2253 JJ=1,NPEQ PQ(NPEQ*(JJ-1)+1)=PCOIN(NPEQ*(JJ-1)+1)*ALFA(IN) PQ(NPEQ*(JJ-1)+2)=PCOIN(NPEQ*(JJ-1)+2)*BETA(IN) PQ(NPEQ*(JJ-1)+3)=PCOIN(NPEQ*(JJ-1)+3)*(-ALFA(IN1)) 2253 CONTINUE C KIRCHHOFF SCATTERING INTEGRAL CALL CMMUL(PCO,NPEQ,NPEQ,PQ,NPEQ,SCATI1) CALL VCLR(SCATI2,1,18) END IF C HOMEGENEOUS PART NRHS=6 PRHS(1)=-ROUGH2(IN1)*(-SA*ALFA(IN))*ALO(IN,1,1) PRHS(4)=-ROUGH2(IN1)*(-SB*BETA(IN))*ALO(IN,1,2) PRHS(7)=-ROUGH2(IN1)*ALFA(IN)*ALO(IN,1,3) PRHS(10)=-ROUGH2(IN1)*BETA(IN)*ALO(IN,1,4) PRHS(13)=-ROUGH2(IN1)*(-ALFA(IN1))*AUP(IN1,1,1) PRHS(16)=-ROUGH2(IN1)*SA1*ALFA(IN1)*AUP(IN1,1,3) DO 2252 JJ=2,3 PRHS(JJ)=-ROUGH2(IN1)*(-SA*ALFA(IN))*ALO(IN,JJ+1,1) PRHS(JJ+3)=-ROUGH2(IN1)*(-SB*BETA(IN))*ALO(IN,JJ+1,2) PRHS(JJ+6)=-ROUGH2(IN1)*ALFA(IN)*ALO(IN,JJ+1,3) PRHS(JJ+9)=-ROUGH2(IN1)*BETA(IN)*ALO(IN,JJ+1,4) PRHS(JJ+12)=-ROUGH2(IN1)*(-ALFA(IN1))*AUP(IN1,JJ+1,1) PRHS(JJ+15)=-ROUGH2(IN1)*SA1*ALFA(IN1)*AUP(IN1,JJ+1,3) 2252 CONTINUE C MULTIPLY WITH SCATTERING INTEGRALS CALL CMMUL(SCATI1,NPEQ,NPEQ,PRHS,NRHS,PRES) C IN NON-KIRCHHOFF CASE ADD OTHER TERM IF (NONKIRCH) THEN CALL CMMUL(SCATI4,NPEQ,NPEQ,PRHS,NRHS,DRES) CALL VADD(PRES,1,DRES,1,PRES,1,2*NRHS*NPEQ) END IF C SOURCE TERMS SIMMILARLY (ONLY PROP TOWARDS INTERFACE) NRHS=2 do is=1,nsrow PRHS(1)=-ROUGH2(IN1)*(-SA*ALFA(IN))*ROIN(IN,1,is) PRHS(4)=-ROUGH2(IN1)*SA1*ALFA(IN1)*RUIN(IN,1,is) DO 2257 JJ=2,3 PRHS(JJ)=-ROUGH2(IN1)*(-SA*ALFA(IN))*ROIN(IN,JJ+1,is) PRHS(JJ+NPEQ)=-ROUGH2(IN1)*SA1*ALFA(IN1)*RUIN(IN,JJ+1,is) 2257 CONTINUE CALL CMMUL(SCATI1,NPEQ,NPEQ,PRHS,NRHS,RRES(1,IS)) C IN NON-KIRCHHOFF CASE ADD OTHER TERM IF (NONKIRCH) THEN CALL CMMUL(SCATI4,NPEQ,NPEQ,PRHS,NRHS,RDRES) CALL VADD(RRES(1,IS),1,RDRES,1,RRES(1,IS),1,2*NRHS*NPEQ) END IF end do IF (NONKIRCH) THEN c c rotation terms only in n-k case c CALL VCLR(PRHS,1,36) NRHS=6 C rotation correction for vertical displacements PRHS(1)=-ROUGH2(IN1)*(-AI*ALO(IN,2,1)) PRHS(4)=-ROUGH2(IN1)*(-AI*ALO(IN,2,2)) PRHS(7)=-ROUGH2(IN1)*(-AI*ALO(IN,2,3)) PRHS(10)=-ROUGH2(IN1)*(-AI*ALO(IN,2,4)) PRHS(13)=-ROUGH2(IN1)*(-AI*AUP(IN1,2,1)) PRHS(16)=-ROUGH2(IN1)*(-AI*AUP(IN1,2,3)) C rotation correction for normal stress PRHS(2)=-ROUGH2(IN1)*(-2*AI*ALO(IN,4,1)) PRHS(5)=-ROUGH2(IN1)*(-2*AI*ALO(IN,4,2)) PRHS(8)=-ROUGH2(IN1)*(-2*AI*ALO(IN,4,3)) PRHS(11)=-ROUGH2(IN1)*(-2*AI*ALO(IN,4,4)) PRHS(14)=0 PRHS(17)=0 C rotation correction for shear stress PRHS(3)=-ROUGH2(IN1)*2E0*(-AI)*CON1(IN)* & (-ALFA(IN)*ALO(IN,1,1)-WVNO*ALO(IN,2,1)) PRHS(6)=-ROUGH2(IN1)*2E0*(-AI)*CON1(IN)* & (-BETA(IN)*ALO(IN,1,2)-WVNO*ALO(IN,2,2)) PRHS(9)=-ROUGH2(IN1)*2E0*(-AI)*CON1(IN)* & (ALFA(IN)*ALO(IN,1,3)-WVNO*ALO(IN,2,3)) PRHS(12)=-ROUGH2(IN1)*2E0*(-AI)*CON1(IN)* & (BETA(IN)*ALO(IN,1,4)-WVNO*ALO(IN,2,4)) PRHS(15)=0 PRHS(18)=0 CALL CMMUL(SCATI2,NPEQ,NPEQ,PRHS,NRHS,DRES) CALL VADD(PRES,1,DRES,1,PRES,1,2*NPEQ*NRHS) CALL CMMUL(SCATI3,NPEQ,NPEQ,PRHS,NRHS,DRES) CALL VADD(PRES,1,DRES,1,PRES,1,2*NPEQ*NRHS) C SOURCE TERMS CALL VCLR(PRHS,1,12) NRHS=2 do is=1,nsrow PRHS(1)=-ROUGH2(IN1)*(-AI*ROIN(IN,2,is)) PRHS(2)=-ROUGH2(IN1)*(-2E0*AI*ROIN(IN,4,is)) PRHS(3)=-ROUGH2(IN1)*2*(-AI)*CON1(IN)* & (-ALFA(IN)*ROIN(IN,1,is)-WVNO*ROIN(IN,2,is)) PRHS(4)=-ROUGH2(IN1)*(-AI*RUIN(IN,2,is)) CALL CMMUL(SCATI2,NPEQ,NPEQ,PRHS,NRHS,RDRES) CALL VADD(RRES(1,IS),1,RDRES,1,RRES(1,IS),1,2*NPEQ*NRHS) CALL CMMUL(SCATI3,NPEQ,NPEQ,PRHS,NRHS,RDRES) CALL VADD(RRES(1,IS),1,RDRES,1,RRES(1,IS),1,2*NPEQ*NRHS) end do END IF C CORRECT LOCAL COEFFICIENT MATRICES C DO 2254 III=1,2 C II=(III-1)*2+1 III=1 II=1 ALO(IN,1,II)=ALO(IN,1,II) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN)*ALFA(IN)) 2 +PRES(1+(II-1)*NPEQ) ALO(IN,1,II+1)=ALO(IN,1,II+1) 1 *(1E0+0.5*ROUGH2(IN1)*BETA(IN)*BETA(IN)) 2 +PRES(1+II*NPEQ) AUP(IN1,1,II)=AUP(IN1,1,II) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN1)*ALFA(IN1)) 2 +PRES(1+(III-1)*NPEQ+12) III=2 II=3 ALO(IN,1,II)=ALO(IN,1,II) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN)*ALFA(IN)) 2 +PRES(1+(II-1)*NPEQ) ALO(IN,1,II+1)=ALO(IN,1,II+1) 1 *(1E0+0.5*ROUGH2(IN1)*BETA(IN)*BETA(IN)) 2 +PRES(1+II*NPEQ) AUP(IN1,1,II)=AUP(IN1,1,II) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN1)*ALFA(IN1)) 2 +PRES(1+(III-1)*NPEQ+12) C2254 CONTINUE DO 2255 JJ=2,3 C DO 2256 III=1,2 C II=(III-1)*2+1 III=1 II=1 ALO(IN,JJ+1,II)=ALO(IN,JJ+1,II) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN)*ALFA(IN)) 2 +PRES(JJ+(II-1)*NPEQ) ALO(IN,JJ+1,II+1)=ALO(IN,JJ+1,II+1) 1 *(1E0+0.5*ROUGH2(IN1)*BETA(IN)*BETA(IN)) 2 +PRES(JJ+II*NPEQ) AUP(IN1,JJ+1,II)=AUP(IN1,JJ+1,II) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN1)*ALFA(IN1)) 2 +PRES(JJ+(III-1)*NPEQ+12) III=2 II=3 ALO(IN,JJ+1,II)=ALO(IN,JJ+1,II) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN)*ALFA(IN)) 2 +PRES(JJ+(II-1)*NPEQ) ALO(IN,JJ+1,II+1)=ALO(IN,JJ+1,II+1) 1 *(1E0+0.5*ROUGH2(IN1)*BETA(IN)*BETA(IN)) 2 +PRES(JJ+II*NPEQ) AUP(IN1,JJ+1,II)=AUP(IN1,JJ+1,II) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN1)*ALFA(IN1)) 2 +PRES(JJ+(III-1)*NPEQ+12) C2256 CONTINUE 2255 CONTINUE do is=1,nsrow ROIN(IN,1,is)=ROIN(IN,1,is) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN)*ALFA(IN)) 2 +RRES(1,is) RUIN(IN,1,is)=RUIN(IN,1,is) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN1)*ALFA(IN1)) 2 +RRES(4,is) DO 2259 JJ=2,3 ROIN(IN,JJ+1,is)=ROIN(IN,JJ+1,is) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN)*ALFA(IN)) 2 +RRES(JJ,is) RUIN(IN,JJ+1,is)=RUIN(IN,JJ+1,is) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN1)*ALFA(IN1)) 2 +RRES(JJ+3,is) 2259 CONTINUE end do GO TO 10 250 CONTINUE C NEXT LAYER OR HALFSPACE IS SOLID C -------------------------------- C SA=0.5+SIGN(0.5,REAL(AK2(IN,1))-REAL(S2)) C SB=0.5+SIGN(0.5,REAL(AK2(IN,2))-REAL(S2)) C SA1=0.5+SIGN(0.5,REAL(AK2(IN1,1))-REAL(S2)) C SB1=0.5+SIGN(0.5,REAL(AK2(IN1,2))-REAL(S2)) SA=1E0 SB=1E0 SA1=1E0 SB1=1E0 NPEQ=4 C CLEAR SCATTERING TERMS CALL VCLR(PRES,1,64) CALL VCLR(DRES,1,64) do is=1,nsrow CALL VCLR(RRES(1,IS),1,32) end do CALL VCLR(RDRES,1,32) C ICHECK=ICHECK+1 IF (IMX(IN1).GE.1) THEN NONKIRCH=.TRUE. CSUM=CMPLX(0E0,0E0) CALL VCLR(SCATI1,1,32) CALL VCLR(SCATI2,1,32) CALL VCLR(SCATI3,1,32) CALL VCLR(SCATI4,1,32) DO 2501 IK=-IMX(IN1),IMX(IN1) DDR=IK*DLWNK(IN1) DDI=SIGN(SCCOIM,REAL(WVNO)+DDR) DIFK=CMPLX(DDR,DDI) KSI=WVNO+DIFK C VERTICAL WAVENUMBERS ALF0=ALF(KSI,IN) BET0=BET(KSI,IN) ALF1=ALF(KSI,IN1) BET1=BET(KSI,IN1) PCO(1)=ALF0 PCO(2)=-KSI PCO(3)=TSSKS(KSI,IN) PCO(4)=-CON1(IN)*2E0*KSI*ALF0 PCO(5)=KSI PCO(6)=-BET0 PCO(7)=CON1(IN)*2E0*KSI*BET0 PCO(8)=-PCO(3) PCO(9)=ALF1 PCO(10)=KSI PCO(11)=-TSSKS(KSI,IN1) PCO(12)=-CON1(IN1)*2E0*KSI*ALF1 PCO(13)=-KSI PCO(14)=-BET1 PCO(15)=CON1(IN1)*2E0*KSI*BET1 PCO(16)=-PCO(11) C IF (ICHECK.EQ.10.AND.IK.EQ.0) THEN C WRITE(49,'(4E15.6)') (PCO(JJ),JJ=1,16) C END IF C COEFFICIENT MATRIX FOR ROTATION RPCO(1)=-AI*PCO(2) RPCO(2)=-AI*PCO(1) RPCO(3)=-2*AI*PCO(4) RPCO(4)=2E0*(-AI)*CON1(IN)* & (2E0*KSI*KSI-AK2(IN,1)) RPCO(5)=-AI*PCO(6) RPCO(6)=-AI*PCO(5) RPCO(7)=-2*AI*PCO(8) RPCO(8)=-2*AI*PCO(7) RPCO(9)=-AI*PCO(10) RPCO(10)=-AI*PCO(9) RPCO(11)=-2E0*AI*PCO(12) RPCO(12)=2E0*(-AI)*CON1(IN1)* & (AK2(IN1,1)-2E0*KSI*KSI) RPCO(13)=-AI*PCO(14) RPCO(14)=-AI*PCO(13) RPCO(15)=-2E0*AI*PCO(16) RPCO(16)=2E0*(-AI)*PCO(15) C C INVERT COEFFICIENT MATRIX C CALL CMATIN(NPEQ,PCO,PCOIN) C DERIVATIVES OF P AND Q (MULT BY VERT WAVENO) DO 2502 JJ=1,NPEQ PQ(NPEQ*(JJ-1)+1)=PCOIN(NPEQ*(JJ-1)+1)*ALF0 PQ(NPEQ*(JJ-1)+2)=PCOIN(NPEQ*(JJ-1)+2)*BET0 PQ(NPEQ*(JJ-1)+3)=PCOIN(NPEQ*(JJ-1)+3)*(-ALF1) PQ(NPEQ*(JJ-1)+4)=PCOIN(NPEQ*(JJ-1)+4)*(-BET1) 2502 CONTINUE C SCATTERING INTEGRALS CALL CMMUL(PCO,NPEQ,NPEQ,PQ,NPEQ,BRES) CALL CMMUL(RPCO,NPEQ,NPEQ,PCOIN,NPEQ,RBRES) SSP=P(DIFK,IN1) DO 2501 JJ=1,16 SCATI4(JJ)=SCATI4(JJ)+AI*DIFK*SSP*RBRES(JJ) SCATI3(JJ)=SCATI3(JJ)-AI*DIFK*SSP*BRES(JJ) SCATI2(JJ)=SCATI2(JJ)+DIFK*DIFK*SSP*RBRES(JJ) 2501 SCATI1(JJ)=SCATI1(JJ)+SSP*BRES(JJ) DO 2503 JJ=1,16 SCATI4(JJ)=ONEO2PI*DLWNK(IN1)*SCATI4(JJ) SCATI3(JJ)=ONEO2PI*DLWNK(IN1)*SCATI3(JJ) SCATI2(JJ)=ONEO2PI*DLWNK(IN1)*SCATI2(JJ) 2503 SCATI1(JJ)=ONEO2PI*DLWNK(IN1)*SCATI1(JJ) ELSE NONKIRCH=.FALSE. C KIRCHHOFF CASE C PCO(1)=ALFA(IN) C PCO(2)=-WVNO C PCO(3)=CON2(IN) C PCO(4)=-CON3(IN) C PCO(5)=WVNO C PCO(6)=-BETA(IN) C PCO(7)=CON4(IN) C PCO(8)=-CON2(IN) DO 251 JJ=1,4 PCO(JJ)=ALO(IN,JJ,3)*EZALFM(IN) PCO(JJ+4)=ALO(IN,JJ,4)*EZBETM(IN) PCO(JJ+8)=AUP(IN1,JJ,1) PCO(JJ+12)=AUP(IN1,JJ,2) 251 CONTINUE C IF (ICHECK.EQ.10) THEN C WRITE(49,'(4E15.6)') (PCO(JJ),JJ=1,16) C END IF C C INVERT COEFFICIENT MATRIX C CALL CMATIN(NPEQ,PCO,PCOIN) C DERIVATIVES OF P AND Q (MULT BY VERT WAVENO) DO 253 JJ=1,NPEQ PQ(NPEQ*(JJ-1)+1)=PCOIN(NPEQ*(JJ-1)+1)*ALFA(IN) PQ(NPEQ*(JJ-1)+2)=PCOIN(NPEQ*(JJ-1)+2)*BETA(IN) PQ(NPEQ*(JJ-1)+3)=PCOIN(NPEQ*(JJ-1)+3)*(-ALFA(IN1)) PQ(NPEQ*(JJ-1)+4)=PCOIN(NPEQ*(JJ-1)+4)*(-BETA(IN1)) 253 CONTINUE C KIRCHHOFF SCATTERING INTEGRAL CALL CMMUL(PCO,NPEQ,NPEQ,PQ,NPEQ,SCATI1) END IF C HOMEGENEOUS PART NRHS=8 DO 252 JJ=1,4 PRHS(JJ)= -ROUGH2(IN1)*(-SA*ALFA(IN))*ALO(IN,JJ,1) PRHS(JJ+4)= -ROUGH2(IN1)*(-SB*BETA(IN))*ALO(IN,JJ,2) PRHS(JJ+8)= -ROUGH2(IN1)*ALFA(IN)*ALO(IN,JJ,3) PRHS(JJ+12)=-ROUGH2(IN1)*BETA(IN)*ALO(IN,JJ,4) PRHS(JJ+16)=-ROUGH2(IN1)*(-ALFA(IN1))*AUP(IN1,JJ,1) PRHS(JJ+20)=-ROUGH2(IN1)*(-BETA(IN1))*AUP(IN1,JJ,2) PRHS(JJ+24)=-ROUGH2(IN1)*SA1*ALFA(IN1)*AUP(IN1,JJ,3) PRHS(JJ+28)=-ROUGH2(IN1)*SB1*BETA(IN1)*AUP(IN1,JJ,4) 252 CONTINUE C MULTIPLY BY SCATTERING INTEGRAL CALL CMMUL(SCATI1,NPEQ,NPEQ,PRHS,NRHS,PRES) C IN KIRCHHOFF CASE ADD OTHER TERM IF (NONKIRCH) THEN CALL CMMUL(SCATI4,NPEQ,NPEQ,PRHS,NRHS,DRES) CALL VADD(PRES,1,DRES,1,PRES,1,2*NRHS*NPEQ) END IF C SOURCE TERMS SIMMILARLY (ONLY PROP TOWARDS INTERFACE) NRHS=2 do is=1,nsrow DO 256 JJ=1,NPEQ PRHS(JJ)=-ROUGH2(IN1)*(-SA*ALFA(IN))*ROIN(IN,JJ,is) PRHS(JJ+NPEQ)=-ROUGH2(IN1)*SA1*ALFA(IN1)*RUIN(IN,JJ,is) 256 CONTINUE CALL CMMUL(SCATI1,NPEQ,NPEQ,PRHS,NRHS,RRES(1,IS)) C IN NON-KIRCHHOFF CASE ADD OTHER TERM IF (NONKIRCH) THEN CALL CMMUL(SCATI4,NPEQ,NPEQ,PRHS,NRHS,RDRES) CALL VADD(RRES(1,IS),1,RDRES,1,RRES(1,IS),1,2*NRHS*NPEQ) END IF end do IF (NONKIRCH) THEN c c rotation terms only in n-k case c CALL VCLR(PRHS,1,64) NRHS=8 C rotation correction for displacements and normal stress DO 2504 JJ=1,NRHS/2 PRHS(NPEQ*(JJ-1)+1)=-ROUGH2(IN1)*(-AI*ALO(IN,2,JJ)) PRHS(NPEQ*(JJ-1)+2)=-ROUGH2(IN1)*(-AI*ALO(IN,1,JJ)) PRHS(NPEQ*(JJ-1)+3)=-ROUGH2(IN1)*(-2E0*AI*ALO(IN,4,JJ)) PRHS(NPEQ*(JJ-1+NPEQ)+1)=-ROUGH2(IN1)*(-AI*AUP(IN1,2,JJ)) PRHS(NPEQ*(JJ-1+NPEQ)+2)=-ROUGH2(IN1)*(-AI*AUP(IN1,1,JJ)) PRHS(NPEQ*(JJ-1+NPEQ)+3)=-ROUGH2(IN1)*(-2E0*AI*AUP(IN1,4,JJ)) 2504 CONTINUE C rotation correction for shear stress PRHS(4)=-ROUGH2(IN1)*2E0*(-AI)*CON1(IN)* & (-ALFA(IN)*ALO(IN,1,1)-WVNO*ALO(IN,2,1)) PRHS(8)=-ROUGH2(IN1)*2E0*(-AI)*CON1(IN)* & (-BETA(IN)*ALO(IN,1,2)-WVNO*ALO(IN,2,2)) PRHS(12)=-ROUGH2(IN1)*2E0*(-AI)*CON1(IN)* & (ALFA(IN)*ALO(IN,1,3)-WVNO*ALO(IN,2,3)) PRHS(16)=-ROUGH2(IN1)*2E0*(-AI)*CON1(IN)* & (BETA(IN)*ALO(IN,1,4)-WVNO*ALO(IN,2,4)) PRHS(20)=-ROUGH2(IN1)*2E0*(-AI)*CON1(IN1)* & (-ALFA(IN1)*AUP(IN1,1,1)-WVNO*AUP(IN1,2,1)) PRHS(24)=-ROUGH2(IN1)*2E0*(-AI)*CON1(IN1)* & (-BETA(IN1)*AUP(IN1,1,2)-WVNO*AUP(IN1,2,2)) PRHS(28)=-ROUGH2(IN1)*2E0*(-AI)*CON1(IN1)* & (ALFA(IN1)*AUP(IN1,1,3)-WVNO*AUP(IN1,2,3)) PRHS(32)=-ROUGH2(IN1)*2E0*(-AI)*CON1(IN1)* & (BETA(IN1)*AUP(IN1,1,4)-WVNO*AUP(IN1,2,4)) CALL CMMUL(SCATI2,NPEQ,NPEQ,PRHS,NRHS,DRES) CALL VADD(PRES,1,DRES,1,PRES,1,2*NPEQ*NRHS) CALL CMMUL(SCATI3,NPEQ,NPEQ,PRHS,NRHS,DRES) CALL VADD(PRES,1,DRES,1,PRES,1,2*NPEQ*NRHS) C SOURCE TERMS CALL VCLR(PRHS,1,16) NRHS=2 do is=1,nsrow PRHS(1)=-ROUGH2(IN1)*(-AI*ROIN(IN,2,is)) PRHS(2)=-ROUGH2(IN1)*(-AI*ROIN(IN,1,is)) PRHS(3)=-ROUGH2(IN1)*(-2E0*AI*ROIN(IN,4,is)) PRHS(4)=-ROUGH2(IN1)*2*(-AI)*CON1(IN)* & (-ALFA(IN)*ROIN(IN,1,is)-WVNO*ROIN(IN,2,is)) PRHS(5)=-ROUGH2(IN1)*(-AI*RUIN(IN,2,is)) PRHS(6)=-ROUGH2(IN1)*(-AI*RUIN(IN,1,is)) PRHS(7)=-ROUGH2(IN1)*(-2*AI*RUIN(IN,4,is)) PRHS(8)=-ROUGH2(IN1)*2*(-AI)*CON1(IN1)* & (ALFA(IN1)*RUIN(IN,1,is)-WVNO*RUIN(IN,2,is)) CALL CMMUL(SCATI2,NPEQ,NPEQ,PRHS,NRHS,RDRES) CALL VADD(RRES(1,IS),1,RDRES,1,RRES(1,IS),1,2*NPEQ*NRHS) CALL CMMUL(SCATI3,NPEQ,NPEQ,PRHS,NRHS,RDRES) CALL VADD(RRES(1,IS),1,RDRES,1,RRES(1,IS),1,2*NPEQ*NRHS) end do END IF C CORRECT LOCAL COEFFICIENT MATRICES DO 254 JJ=1,NPEQ DO 255 II=1,3,2 ALO(IN,JJ,II)=ALO(IN,JJ,II) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN)*ALFA(IN)) 2 +PRES(JJ+(II-1)*NPEQ) ALO(IN,JJ,II+1)=ALO(IN,JJ,II+1) 1 *(1E0+0.5*ROUGH2(IN1)*BETA(IN)*BETA(IN)) 2 +PRES(JJ+II*NPEQ) AUP(IN1,JJ,II)=AUP(IN1,JJ,II) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN1)*ALFA(IN1)) 2 +PRES(JJ+(II-1)*NPEQ+16) AUP(IN1,JJ,II+1)=AUP(IN1,JJ,II+1) 1 *(1E0+0.5*ROUGH2(IN1)*BETA(IN1)*BETA(IN1)) 2 +PRES(JJ+II*NPEQ+16) 255 CONTINUE 254 CONTINUE do is=1,nsrow DO 258 JJ=1,NPEQ ROIN(IN,JJ,is)=ROIN(IN,JJ,is) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN)*ALFA(IN))+RRES(JJ,is) RUIN(IN,JJ,is)=RUIN(IN,JJ,is) 1 *(1E0+0.5*ROUGH2(IN1)*ALFA(IN1)*ALFA(IN1))+RRES(JJ+4,is) 258 CONTINUE end do 10 CONTINUE RETURN END