PROGRAM TRIM82 INTEGER RAYFIL, INPFIL, PRTFIL, INTFIL, DSTFIL, IGNFIL PARAMETER ( RAYFIL = 21, INPFIL = 5, PRTFIL = 6, INTFIL = 7, & DSTFIL = 37, IGNFIL = 41 ) DIMENSION Z(2,50),V(2,50) COMMON /MIRROR/ FSRL, BRLT(200), BPST(200) DATA BRLT/200*1./, BPST/200*0./ COMMON /INFO/ R0, R, OMEGA, ATT,ICOH,IRAN,IT2,IT3,IGNRAY,IFLMR COMMON /PROF/ RONE, N1, Z1(50), V1(50), 1 RTWO, N2, Z2(50), V2(50) COMMON /LOUDNS/ NRCD, RCD(100), XINT(400) COMMON /TLPLOT/ RKM(1000), XTL(1000,8), ITYP(1000), IRCD(1000) DATA ITYP/1000*-1./ COMMON /PRO/ IPROP, ICUR, REARTH, R1KM, DRKM, R2KM COMMON /PATTRN/ SD, ITBP, DATD, SORLEV COMMON /PIDEF/ PI, DEGREE, TWO PI COMMON /TLE/ ITITLE (18) COMMON/ABC/ NBRS, NSRS, ALIM COMMON /RANST/ RSKM, IREC, IFSK, SDS DATA XLG/1000./, YHGT/6./, XL/1./, YL/6.4/, MAXPLOT/1000/ DATA IPROP/1/, ICUR/1/, IFLMR/1/, REARTH/6371221.3/ C***** VAX CONVERSION OF TRIMAIN--CF. N R L REPORT # 7827, 12/16/74. C E B WRIGHT CODE 5160 NAVAL RESEARCH LAB, WASHINGTON D C 20375. C***** TRY TO OBSERVE VARIABLE-NAME CONVENTION: LENGTHS WITHOUT 'KM' C OR SOME OTHER UNIT QUALIFIER ARE ALWAYS IN METERS. PI = 3.141592654 TWO PI = PI + PI DEGREE = 180./PI READ( INPFIL, 900 ) ITITLE 900 FORMAT(20A4) 901 FORMAT(16I5) 902 FORMAT(10F8.0) 903 FORMAT(16F5.0) 907 FORMAT(' RECEIVER RANGES R1KM, DRKM, R2KM = ', 3F10.3) 908 FORMAT(' RECEIVER DEPTHS(M) = ', 10F10.3) WRITE( PRTFIL, 900) ITITLE READ( INPFIL, 902 ) SD, FHZ, SORLEV, DATD, SLDB, XATT, XBOT IFLMR = 1 IF(SD.GT.200. .OR. SD*FHZ.GT.15000.) IFLMR = 0 ITBP = 0 IF(DATD.NE.0.) ITBP = 1 IF(DATD.GT.90.) DATD = 0. READ( INPFIL, 901 ) ICOH, IRAN, IT2, IT3, IGNRAY, IRD, IRF, 1 IPL, NSRS, NBRS, LIMDB INT = IRAN IF(ICOH.LE.0 .AND. IGNRAY.LE.0) GO TO 15 IF(ICOH.GT.0 .AND. IRAN.LE.0) IRAN = ICOH IF(IGNRAY.GT.0 .AND. IRAN.LE.0) IRAN = IGNRAY IF(ICOH.GT.0) INT = MIN0(ICOH,IRAN) IF(IGNRAY.GT.0) INT = MIN0(IGNRAY,INT) IF(IGNRAY.GT.0 .AND. INT.EQ.0) INT = IGNRAY 15 IF(IT2.GT.0 .AND. INT.LE.0) INT = IT2 IF(IT3.GT.0 .AND. INT.LE.0) INT = IT3 FSRL = 10.**(-.1*SLDB) OMEGA = TWO PI*FHZ FSQ = FHZ*FHZ FKHZSQ = FSQ*1.E-6 WRITE( PRTFIL, 902) SD, SORLEV, FHZ, DATD SDS = SD IF (ICUR.NE.0) SD = SD*(1.+.5*SD/REARTH) ATT=0. C***** RELAXATION FREQUENCY AT 4.4 C IS ABOUT 70 KHZ. ATT IN DB/KM. C THE FIRST FORMULA IS FOR THORP ATTENUATION. C THE SECOND IS CYBULSKI/GOODMAN (?). IF(XATT.NE.0.) ATT = 0.10936*FKHZSQ/(1.+FKHZSQ) + 1 50.*FKHZSQ/(5000. + FKHZSQ) C IF(XATT.EQ.2.) ATT = 0.00044*FKHZSQ + 50.*FKHZSQ/(5000. + FKHZSQ) WRITE( PRTFIL, 915) ATT 915 FORMAT(' VOLUME ATTENUATION IN DB/KM = ', E10.3) IF(NBRS.LE.0) NBRS = 9999 IF(NSRS.LE.0) NSRS = 9999 IF(LIMDB.LE.0) LIMDB=200 WRITE( PRTFIL, 901) ICOH, IRAN, IT2, IT3, NSRS, NBRS, LIMDB ALIM = 10.**(-.1*LIMDB) NRC = 7 NRCD = 0 READ( INPFIL, 902 ) R1KM, DRKM, R2KM, (RCD(I),I=1,7) IF(R1KM.EQ.DRKM) R1KM = 0. DELR12KM = R2KM - R1KM R1 = 1000.*R1KM DR = 1000.*DRKM R2 = 1000.*R2KM IF(RCD(7).LE.0.) GO TO 20 NRC = 17 READ( INPFIL, 902 ) (RCD(I),I=8,17) IF(RCD(17).LE.0.) GO TO 20 NRC = 27 READ( INPFIL, 902 ) (RCD(I),I=18,27) 20 DO 30 I=1,NRC IF(RCD(I).LE.0.) GO TO 30 NRCD = I IF(ICUR.NE.0) RCD(I) = RCD(I)*(1.+.5*RCD(I)/REARTH) 30 CONTINUE 31 WRITE( PRTFIL, 907) R1KM, DRKM, R2KM IF(NRCD.GT.0 .OR. INT.EQ.0) GO TO 32 NRCD = 1 RCD(1) = 1. 32 WRITE( PRTFIL, 908) (RCD(I),I=1,NRCD) CALL INITRY RBKM = 50000. IF(XBOT.GE.0.) GO TO 36 WRITE( PRTFIL, *) 'BOTTOM IS PERFECTLY ABSORBING.' DO 35 I=1,200 BRLT(I) = 0. 35 CONTINUE RBLKM = 50000. GO TO 40 36 IF(XBOT.GT.0.) GO TO 38 WRITE( PRTFIL, *) 'BOTTOM IS PERFECTLY REFLECTING.' RBLKM = 50000. GO TO 40 38 WRITE( PRTFIL, *) 'BOTTOM LOSS CLASS TO BE INPUT.' XBOT = ICOH CALL BOTLOSS(XBOT) RBLKM = XBOT 40 RBL = 1000.*RBLKM CALL NEWPRO(ZM) CALL NEWPRO(ZM) CALL CONNCT IF(ZMAX.EQ.0.) ZMAX=ZM R0 = 0. R = 0. IF(R1.NE.0.) GO TO 50 IF(IRF.GT.0) CALL RAYOUT GO TO 95 50 IF(RTWO.GE.R1 .AND. RBL.GE.R1) GO TO 75 R = AMIN1(RBL,RTWO) CALL ADVNCE R0 = R IF(R.LT.RBL) GO TO 55 CALL NEWBLT(RBLKM) RBL = 1000.*RBLKM 55 IF(R.LT.RTWO) GO TO 50 CALL NEWPRO(ZM) CALL CONNCT GO TO 50 75 R = R1 CALL ADVNCE R0 = R IF(INT.GT.0) CALL INTSTY(0) IF(IRD.GT.0) CALL ZDIST IF(IRF.GT.0) CALL RAYOUT IF(R.NE.RBL) GO TO 90 CALL NEWBLT(RBLKM) 90 IF(R.NE.RTWO) GO TO 95 CALL NEWPRO(ZM) CALL CONNCT 95 NRS = (R2 - R1)/DR + .5 NPLOTS = 0 IF(IPL.EQ.0 .OR. NRS.LE.MAXPLOT) GO TO 105 WRITE( PRTFIL, 986) NPLOTS 986 FORMAT(' TOO MANY RANGE STEPS FOR NPLOTS=', I6) STOP 105 NEXTI = INT NEXTRD = IRD NEXTRF = IRF NEXTBL = (RBLKM-R1KM)/DRKM + .5 NEXTP = (RTWO-R1)/DR + .5 DO 1000 NR=1,NRS R = R + DR CALL ADVNCE R0 = R IF(NR.NE.NEXTI) GO TO 200 CALL INTSTY(NEXTI) NEXTI = NEXTI + INT 200 IF(NR.NE.NEXTRD) GO TO 300 CALL ZDIST NEXTRD = NEXTRD + IRD 300 IF(NR.NE.NEXTRF) GO TO 400 CALL RAYOUT NEXTRF = NEXTRF + IRF 400 IF(NR.NE.NEXTP) GO TO 500 CALL NEWPRO(ZM) CALL CONNCT NEXTP = (RTWO-R1)/DR + .5 500 IF(NR.NE.NEXTBL) GO TO 800 CALL NEWBLT(RBLKM) NEXTBL = (RBLKM-R1KM)/DRKM + .5 800 IF(INT.GT.0) CALL INTPRT 1000 CONTINUE IF(IPL.EQ.0) GO TO 73 NPLOTS = (ICOH+IRAN+IT2+IT3)*NRCD XLG = 8. IF(NRS.GT.500) XLG = 12. YHGT = 6. DELRKM = 100. IF(DELR12KM.LT.675.) DELRKM = 50. IF(DELR12KM.LT.375.) DELRKM = 20. IF(DELR12KM.LT.175.) DELRKM = 10. IF(DELR12KM.LT.65.) DELRKM = 5. IF(DELR12KM.LT.33.) DELRKM = 2. IF(DELR12KM.LT.15.) DELRKM = 1. 73 STOP END SUBROUTINE BOTLOSS(RBKM) INTEGER RAYFIL, INPFIL, PRTFIL, INTFIL, DSTFIL, IGNFIL PARAMETER ( RAYFIL = 21, INPFIL = 5, PRTFIL = 6, INTFIL = 7, & DSTFIL = 37, IGNFIL = 41 ) DATA IENT/0/ COMMON /PIDEF/ PI, DEGREE, TWO PI COMMON /INFO/ ROLD, RANGE, OMEGA, IDUMM(7) COMMON /MIRROR/ FSRL,BRLT(200),BPST(200) DIMENSION DB(50),DG(50),DTSTR(600),BR(91,4),BP(91,4),ISW(4), * RBLKM(50),ICL(50), DTH(600) COMMON/ABC/NBRS,NSRS,ALIM EQUIVALENCE (DB,DTSTR(551)),(DG,DTSTR(501)),(ICL,DTSTR(451)), *(RBLKM,DTSTR(401)),(ISW,DTSTR(397)),(BR,DTSTR),(ICOH,DTSTR(396)), * (IE,DTSTR(395)),(IB,DTSTR(394)),(K,DTSTR(393)), * (IL,DTSTR(392)),(IH,DTSTR(391)),(IC,DTSTR(390)) EQUIVALENCE (N,DTSTR(389)),(T,DTSTR(388)),(F,DTSTR(387)), * (PHI,DTSTR(386)),(DCB,DTSTR(385)),(RKM,DTSTR(384)), * (RNKM,DTSTR(383)) EQUIVALENCE (BP,DTSTR(183)), (DTSTR,DTH1) COMMON/BLDATA/DTH1(100),DTH2(100),DTH3(100),DTH4(100),DTH5(100), 1 DTH6(100) DATA DTH1 /7.5, 11., 13.5, 15., 17.5, 20., 25., 27.5, 2*90., 2 7.5, 10., 11.5, 14., 15., 17.5, 20., 25., 29.5, 90., 3 7.5, 12., 15., 18., 20., 22.5, 25., 27., 2*90., 4 7.5, 10., 14.5, 17.5, 20., 22.5, 25., 28., 2*90., 5 2.5, 5.5, 8.0, 10., 15., 17.5, 20., 23., 2*90., 1 6.5, 8.6, 9.8, 10.8, 11., 11.7, 12.7, 3*13., 2 7.0, 8.8, 10., 10.6, 11.1, 11.9, 12.5, 13.4, 2*14., 3 8.0, 11., 12.8, 14., 14.8, 15.4, 15.8, 3*16., 4 9.0, 11., 13.5, 14.8, 15.6, 16.2, 16.6, 3*17., 5 10., 13., 15., 16., 17.8, 18.4, 18.7, 3*19./ DATA DTH2 /2.5, 10., 12.5, 15., 20., 25., 30., 35., 40., 90., 2 2.5, 7.5, 12.5, 20., 25., 32.5, 37.5, 42.5, 2*90., 3 2.5, 5.0, 7.5, 12.5, 17.5, 22.5, 27.5, 32.5, 2*90., 4 2.5, 7.5, 10., 11.5, 13., 15., 20., 24., 2*90., 5 2.5, 5.0, 8.5, 10., 12.5, 15., 20., 22.5, 27.5, 90., 1 3.5, 5.4, 6.0, 6.6, 7.7, 8.6, 9.4, 9.8, 2*10., 2 6.5, 8.0, 9.5, 11., 11.9, 13., 13.4, 3*13.5, 3 7.0, 9.4, 11.5, 14., 14.7, 15.1, 15.4, 3*15.5, 4 9.0, 13.5, 15.2, 16., 16.6, 17., 17.4, 3*17.5, 5 11., 13.6, 16., 17., 18.1, 18.8, 19.6, 19.8, 2*20./ DATA DTH3 /18., 23., 30., 40., 50., 55., 4*90., 2 17., 20., 25., 30., 35., 45., 55., 3*90., 3 13., 25., 30., 35., 42.5, 45., 50., 3*90., 4 9.0, 12.5, 20., 22.5, 27., 5*90., 5 2.5, 4.5, 7.5, 10., 14.5, 17.5, 20., 23., 2*90., 1 0., 2.0, 4.0, 6.0, 7.4, 5*8., 2 0., 1.7, 3.9, 5.6, 6.8, 8.8, 4*10., 3 3.0, 6.6, 8.0, 9.1, 10.4, 10.7, 4*11., 4 7.0, 8.0, 12., 13.1, 6*14., 5 9.0, 10., 12., 13.9, 16., 17.1, 17.7, 3*18./ DATA DTH4 /19., 25., 35., 45., 51., 5*90., 2 18., 25., 35., 45., 55., 5*90., 3 14., 20., 30., 40., 50., 53., 4*90., 4 8.0, 11., 20., 22.5, 26., 5*90., 5 2.5, 4.5, 7.5, 10., 14.5, 17.5, 20., 23., 2*90., 1 0., 1.5, 3.5, 5.1, 6*6., 2 0., 2.8, 5.6, 7.6, 6*9., 3 3.0, 5.0, 7.5, 9.4, 10.8, 5*11., 4 5.0, 6.1, 12.2, 13.2, 6*14., 5 9.0, 10., 12., 13.9, 16., 17.1, 17.7, 3*18./ DATA DTH5 /18., 30., 40., 50., 55., 5*90., 2 17., 20., 30., 40., 50., 55., 4*90., 3 13., 20., 35., 45., 52.5, 5*90., 4 7.5, 11., 20., 25., 30., 34., 4*90., 5 2.5, 5.0, 7.5, 15., 17.5, 20., 22., 3*90., 1 0., 2.6, 4.4, 5.6, 6*6., 2 0., 1.3, 4.6, 7.0, 8.5, 5*9., 3 3.0, 5.3, 8.7, 10.3, 6*11., 4 3.0, 4.0, 10.4, 12.7, 13.7, 5*14., 5 8.0, 8.7, 10., 14.1, 15., 15.8, 4*16./ DATA DTH6 /19., 35., 40., 52., 6*90., 2 19., 25., 35., 45., 55., 5*90., 3 11., 20., 25., 35., 45., 56., 4*90., 4 11., 20., 30., 40., 55., 5*90., 5 5., 7.5, 16.5, 17.5, 19., 20., 21., 22.5, 2*90., 1 0., 3.7, 4.5, 7*6., 2 0., 2.3, 4.6, 6.5, 6*8., 3 0., 3.0, 4.4, 6.7, 8.5, 5*10., 4 2., 5.2, 7.7, 9.8, 6*12., 5 4., 5.2, 12., 12.8, 13.4, 13.7, 13.9, 3*14./ IF(IENT.EQ.1) GO TO 2 DO 1 J=1,600 1 DTH(J) = DTSTR(J) IENT = 1 GO TO 3 2 DO 4 J=1,600 4 DTSTR(J) = DTH(J) IB=0 3 FR=OMEGA/TWOPI IF(FR.LT.300.) GO TO 25 IF(FR.GT.1500.) GO TO 10 IF(FR.LT.750.) GO TO 5 IQ=301 GO TO 20 5 IQ=401 GO TO 20 10 IF(FR.LT.2750.) GO TO 15 IQ=101 IF(FR.GT.6000.) IQ=1 GO TO 20 15 IQ=201 20 DO 24 I=501,600 DTSTR(I)=DTH1(IQ) 24 IQ=IQ+1 25 DO 30 I=1,500 30 DTSTR(I)=0. ICOH=RBKM 40 IB=IB+1 READ( INPFIL, 900 ) RBLKM(IB), ICL(IB) IF(RBLKM(IB)) 41,41,42 41 RBLKM(IB)= 50000. 42 IF(ICL(IB)) 43,43,44 43 ICL(IB)=0 GO TO 60 44 IF(ICL(IB)-9) 46,46,45 45 ICL(IB)=ICL(IB)/10 46 IF(ICL(IB)-5) 60,60,47 47 K=ICL(IB)-5 IF(ISW(K)) 48,48,60 48 ISW(K)=1 IL=1 49 IH=MIN0(91,IL+15) READ( INPFIL, 901 ) (BR(I,K), I=IL,IH) IL=IH + 1 IF(BR(IH,K)) 51,51,49 50 IF(BR(IH,K)) 51,51,53 51 IH=IH-1 IF(IH) 52,52,50 52 IH=1 53 DO 54 I=IH,90 54 BR(I+1,K)=BR(IH,K) IL=1 IF(ICOH) 60,60,55 55 IH=MIN0(91,IL+15) READ( INPFIL, 901 ) (BP(I,K),I=IL,IH) IL=IH + 1 IF(BP(IH,K)) 55,60,55 60 IF(RBLKM(IB).LT.40000.) GO TO 40 DO 65 K=6,9 IF(ISW(K-5)) 65,65,61 61 WRITE( PRTFIL, 902) K WRITE( PRTFIL, 903) (I,BR(I+1,K-5),BP(I+1,K-5),I=IE,90) 65 CONTINUE WRITE( PRTFIL, 900) IF(IB.LT.2) GO TO 71 DO 70 I=2,IB RNKM =RBLKM(I-1) WRITE( PRTFIL, 904) RKM,RNKM,ICL(I-1) 70 RKM =RNKM 71 WRITE( PRTFIL, 905) RKM,ICL(IB) IE=1 75 IC=ICL(IE) IF(IC) 76,76,80 76 DO 77 I=1,200 BRLT(I)=1. 77 BPST(I)=0. GO TO 100 80 DO 99 I=1,200 T=DEGREE*ATAN(.01*(I-1)) K=IC-5 IF(K) 81,81,95 81 IL=10*IC IF(DG(IL-9).GE.T) GO TO 93 82 IF(DG(IL-8).GE.T) GO TO 92 IL=IL+1 GO TO 82 92 DCB=DB(IL-9)+(T-DG(IL-9))*(DB(IL-8)-DB(IL-9))/(DG(IL-8)-DG(IL-9)) GO TO 94 93 DCB=DB(IL-9) 94 PHI=0. GO TO 98 95 N=T F=T-N DCB=BR(N+1,K)+F*(BR(N+2,K)-BR(N+1,K)) PHI=0. IF(ICOH.NE.0) PHI=BP(N+1,K)+F*(BP(N+2,K)-BP(N+1,K)) 98 BPST(I)=PHI 99 BRLT(I) = 10.**(-.1*DCB) 100 RBKM=RBLKM(IE) RETURN ENTRY NEWBLT(RBKM) IE=IE+1 IF(IC-ICL(IE)) 75,100,75 900 FORMAT(F8.4,I2) 901 FORMAT(16F5.2) 902 FORMAT(38H0USER SUPPLIED BOTTOM LOSS TABLE CLASS ,I5/, 1 5(7X,'GA DB PHI ')) 903 FORMAT(5(I9,F9.3,F8.3)) 904 FORMAT(5H FROM,F10.2,3H TO,F10.2,16H KM BOTTOM CLASS,I5) 905 FORMAT(5H FROM,F10.2,3H TO,3X,23HEND OF RUN BOTTOM CLASS,I5) END SUBROUTINE NEWPRO(ZMAX) INTEGER RAYFIL, INPFIL, PRTFIL, INTFIL, DSTFIL, IGNFIL PARAMETER ( RAYFIL = 21, INPFIL = 5, PRTFIL = 6, INTFIL = 7, & DSTFIL = 37, IGNFIL = 41 ) COMMON /INFO/ ROLD, RANGE, OMEGA,ATT,ICOH,IRAN,IT2,IT3,IGNRAY,IFLMR COMMON /PROF/ R1,N1,Z1(50),V1(50),R2,N2,Z2(50),V2(50) COMMON/INPUT/ RBOTKM(250),ZBOT(250),ZIL(2,50),VIL(2,50), * NN(2), ICON(2,50) COMMON /TLE/ ITITLE (18) COMMON /PRO/ IPROP, ICUR, REARTH, R1KM, DRKM, R2KM DIMENSION TITLE(10),PM(3),LIN (101) COMMON/ABC/NBRS,NSRS,ALIM COMMON /RANST/ RSKM,IDUM,IDUMM,SDS DATA NBPT/0/ DATA IFEOF /0/ DATA LIN/101*1H / DATA IENT/0/ DATA IBLANK /' '/ , IVEE/'V '/ , ISTAR/'* '/ DATA RAKM/0./, RBKM/0./ 952 IF (IENT.EQ.1) GO TO 950 NBPT=0 IFEOF=0 IENT=1 N2=0 R2 = 0. N1=0 R1 = 0. DO 951 J=1,101 951 LIN(J)=IBLANK 950 IF(NBPT.GT.0) GO TO 20 C FIRST ENTRY. READ IN BOTTOM TRACK IL=1 2 IH=IL+4 READ( INPFIL, 900 ) (RBOTKM(I), ZBOT(I), I=IL,IH) IF(RBOTKM(IH).LE.0.) GO TO 3 IL=IL+5 IF(IH.LT.250) GO TO 2 IH = 251 3 IH=IH-1 IF(RBOTKM(IH).LE.0.) GO TO 3 NBPT=IH IF (ICUR.NE.0) GO TO 1 DO 5 I=1,NBPT 5 ZBOT(I)=ZBOT(I)*(1.+.5*ZBOT(I)/REARTH) 1 CONTINUE 611 FORMAT (20A4,/,(10F8.3)) 911 FORMAT (10F8.3) ZMAX=0. DO 4 I=1,NBPT XI = (RBOTKM(I)-R1KM)/DRKM II = XI + .5 IF(XI.GE.0.) RBOTKM(I) = R1KM + II*DRKM IF(ZBOT(I).GT.ZMAX) ZMAX=ZBOT(I) 4 CONTINUE WRITE( PRTFIL, 901)(I,RBOTKM(I),ZBOT(I),I=1,NBPT) IBPT=0 DO 55 I=1,2 DO 55 J=1,50 ZIL(I,J) = 0. VIL(I,J) = 0. 55 ICON(I,J)=0 NN(1)=0 NIB=0 DO 56 J=1,50 Z1(J) = 0. V1(J) = 0. Z2(J) = 0. 56 V2(J) = 0. ASSIGN 6 TO IRET GO TO 70 C AT 70 IS THE READIN AND CONNECT ROUTINE. RETURN IS TO IRET 6 ASSIGN 20 TO IRET GO TO 70 20 IBPT=IBPT+1 IF(IBPT.LE.NBPT) GO TO 21 RBTKM=.001*R2+100. GO TO 22 21 RBTKM=RBOTKM(IBPT) ZBT=ZBOT(IBPT) 22 IF(RBTKM.LE.RBKM) GO TO 25 C GET NEW INPUT PROFILE ASSIGN 23 TO IRET GO TO 70 23 IF (IBPT.GT.NBPT .OR. RBOTKM(IBPT-1).EQ.RAKM) GO TO 25 IBPT=IBPT-1 ZBT=ZBT+(RAKM-RBTKM)*(ZBT-ZBOT(IBPT))/(RBTKM-RBOTKM(IBPT)) C THIS IS THE BOTTOM DEPTH AT PROFILE RANGE RBTKM=RAKM 25 R1=R2 N1=N2 C MOVE PROFILE 2 TO PROFILE 1 IF(N2.LT.0 .OR. N2.GT.50) GO TO 27 DO 26 I=1,N2 Z1(I)=Z2(I) 26 V1(I)=V2(I) 27 R2=1000.*RBTKM FA=(RBKM-RBTKM)/(RBKM-RAKM) FB=1.-FA N2=1 29 IF (ICON(1,N2).GT.NN(1)) ICON (1,N2)=NN(1) IF (ICON(2,N2).GT.NN(2)) ICON (2,N2)=NN(2) Z=FA*ZIL(1,ICON(1,N2))+FB*ZIL(2,ICON(2,N2)) V=FA*VIL(1,ICON(1,N2))+FB*VIL(2,ICON(2,N2)) IF (N2.LE.1 .OR. Z.GE.Z2(N2-1)) GO TO 126 125 Z=Z2(N2-1) V=V2(N2-1) 126 IF(Z.GT.ZBT-1.) GO TO 35 Z2(N2)=Z V2(N2)=V N2=N2+1 IF(N2.LE.NIB) GO TO 29 C EXTRAPOLATION TO BOTTOM Z2(N2)=ZBT M=N2-1 31 M=M-1 Z=Z2(M) IF(ABS(Z-Z2(N2-1)).LT.Z*.00001) GO TO 31 V=V2(M) V2(N2)=V2(M) +(ZBT-Z)*(V2(N2-1)-V)/(Z2(N2-1)-Z) GO TO 40 35 Z2(N2)=ZBT V2(N2)=V2(N2-1)+(ZBT-Z2(N2-1))*(V-V2(N2-1))/(Z-Z2(N2-1)) C REMOVE DUPLICATE POINTS 40 M=2 IDC=0 42 IF (Z2(M).GT.Z2(M-1)+1.) GO TO 45 IDC=1 N2=N2-1 DO 44 I=M,N2 Z2(I)=Z2(I+1) 44 V2(I)=V2(I+1) 45 M=M+1 IF(M.LE.N2) GO TO 42 IF (IDC.GT.0) GO TO 40 IF (Z2(N2+1).EQ.0.) GO TO 46 N5=N2+1 DO 47 I=N5,50 Z2(I) = 0. 47 V2(I) = 0. 46 IF (IPROP.EQ.0.OR.IPROP.EQ.1) GO TO 54 WRITE( PRTFIL, 902) R2 WRITE( PRTFIL, 914) DO 7 J=1,N2 NV=V2(J)-1455.5 IF (NV.LT.4.OR.NV.GT.101) GO TO 48 LIN(NV)= ISTAR WRITE( PRTFIL, 915) Z2(J),V2(J),LIN LIN(NV)= IBLANK GO TO 7 48 WRITE( PRTFIL, 915) Z2(J),V2(J) 7 CONTINUE C54 IF (IPFL.EQ.2) CALL PRPLIT (RBKM,N3,INPF) 54 IF (RSKM.GE.RBOTKM(2).OR.IFEOF.EQ.1) RETURN WRITE( PRTFIL, 275) STOP 275 FORMAT (5X,'PROGRAM ABORTED, TWO SOUND SPEED PROFILES INPUT BEFORE * SECOND BOTTOM POINT') 70 IF(NN(2).EQ.0) GO TO 75 NB=NN(2) DO 71 I=1,NB ZIL(1,I)=ZIL(2,I) 71 VIL(1,I)=VIL(2,I) RAKM=RBKM 75 NN(1)=NN(2) IF(IFEOF.EQ.0) GO TO 76 RBKM=RBKM+1.E6 GO TO 82 76 READ( INPFIL, 903 ) NCUR, RBKM, TITLE IF(NCUR.GE.0) GO TO 78 77 IFEOF=1 RBKM=RAKM+1.E6 GO TO 82 78 IL=1 RSKM=RBKM 79 IH=IL+4 READ( INPFIL, 900 ) (ZIL(2,I), VIL(2,I), I=IL,IH) IF(VIL(2,IH).LE.0.) GO TO 81 IL=IL+5 IF(IH.LT.50) GO TO 79 IH = 51 81 IH=IH-1 IF(VIL(2,IH).LE.0.) GO TO 81 NN(2)=IH N3=IH INPF=0 IF (IPROP.EQ.0) GO TO 74 WRITE( PRTFIL, 904) RBKM,TITLE WRITE( PRTFIL, 914) 74 DO 28 I=1,IH IF (VIL(2,I).GE.1400. .AND. VIL(2,I).LE.10000.) GO TO 57 58 WRITE( PRTFIL, 920) RBKM,ZIL(2,I),VIL(2,I) 920 FORMAT(1H1,5X,'INPUT ERROR IN PROFILE AT RANGE',F10.4,5X,'VALUE AT * DEPTH',F10.3,'=',F10.3,' M/SEC') STOP 57 IF(ZIL(2,I).NE.SDS) GO TO 257 ZIL(2,I)=ZIL(2,I)+0.1 WRITE( PRTFIL, 258) SDS,RBKM, ZIL(2,I) 258 FORMAT (5X,'PROGRAM WARNING, INPUT PROFILE DEPTH',F10.3,1X,'=SOURC *E DEPTH, AT RANGE',F10.3,' KM',5X,'NEW PROFILE DEPTH=',F10.3,/) 257 IF (I.EQ.1) GO TO 157 C***** ISOVELOCITY LAYERS ELIMINATED 9/80 BECAUSE THEY CAN LEAD TO C ROTATION ANGLE OF 90 DEGREES FOR A TRIANGLE AND DIVISION BY THE C COSINE OF ROTATION ANGLE=CT IS USED IN RAY INTENSITY IN ADVNCE. C***** THE ASC PRECISION FOR REAL*4 NUMBERS BETWEEN 1024. AND 2048. IS C 2048/16**6 = .000122. USE .0004 AS A DISTINGUISHABLE DIFFERENCE. IF(VIL(2,I).EQ.VIL(2,I-1)) VIL(2,I) = VIL(2,I) + .0004 IF (ZIL(2,I).LT.ZIL (2,I-1)) GO TO 58 157 IF (IPROP.EQ.0) GO TO 28 NV=VIL(2,I)-1455.5 IF (NV.LT.4.OR.NV.GT.101) GO TO 38 LIN(NV)= IVEE WRITE( PRTFIL, 915) ZIL(2,I),VIL(2,I),LIN LIN (NV)= IBLANK GO TO 28 38 WRITE( PRTFIL, 915) ZIL(2,I),VIL(2,I) 28 CONTINUE IF (NCUR.NE.0) GO TO 89 WRITE( PRTFIL, *) 'CURVED-EARTH CORRECTION APPLIED TO PROFILE.' DO 80 I=1,IH VIL(2,I)=VIL(2,I)*(1.+ZIL(2,I)/REARTH) 80 ZIL(2,I)=ZIL(2,I)*(1.+.5*ZIL(2,I)/REARTH) C89 IF (IPFL.EQ.2) CALL PRPLIT (RBKM,N3,INPF) 89 IF(NN(1).EQ.0) GO TO IRET,(6,20,23) 82 ICON(1,1)=1 ICON(2,1)=1 NA=2 NB=2 NIB=2 83 IF (ZIL(2,NIB).EQ.ZIL(1,NIB)) GO TO 92 IF(NB.LE.NN(2)) GO TO 84 IP=2 GO TO 88 84 IF(NA.LE.NN(1)) GO TO 85 IP=1 GO TO 88 85 DVMA=VIL(1,NA)-VIL(1,NA-1) DVLA=DVMA DVHA=DVMA IF(NA.GT.2) DVLA=VIL(1,NA-1)-VIL(1,NA-2) IF(NA.LT.NN(1)) DVHA=VIL(1,NA+1)-VIL(1,NA) DVMB=VIL(2,NB)-VIL(2,NB-1) DVLB=DVMB DVHB=DVMB IF(NB.GT.2) DVLB=VIL(2,NB-1)-VIL(2,NB-2) IF(NB.LT.NN(2)) DVHB=VIL(2,NB+1)-VIL(2,NB) PM(1)=ABS(ZIL(1,NA-1)-ZIL(2,NB)) IF(DVLA*DVMB.GE.0..AND.DVMA*DVHB.GE.0.) PM(1)=PM(1)-250. PM(2)=ABS(ZIL(1,NA)-ZIL(2,NB-1)) IF(DVMA*DVLB.GE.0..AND.DVHA*DVMB.GE.0.) PM(2)=PM(2)-250. PM(3)=ABS(ZIL(1,NA)-ZIL(2,NB)) IF(DVMA*DVMB.GE.0..AND.DVHA*DVHB.GE.0.) PM(3)=PM(3)-250. BP=PM(1) IP=1 DO 87 I=2,3 IF(PM(I).GE.BP) GO TO 87 IP=I BP=PM(I) 87 CONTINUE 88 GO TO (90,91,92),IP 90 ICON(1,NIB)=NA-1 ICON(2,NIB)=NB NB=NB+1 GO TO 93 91 ICON(1,NIB)=NA ICON(2,NIB)=NB-1 NA=NA+1 GO TO 93 92 ICON(1,NIB)=NA ICON(2,NIB)=NB NA=NA+1 NB=NB+1 93 IF(NA.GT.NN(1).AND.NB.GT.NN(2)) GO TO IRET,(6,20,23) NIB=NIB+1 IF (NIB.LE.50) GO TO 83 WRITE( PRTFIL, 905) STOP 900 FORMAT(10F8.4) 901 FORMAT('0LISTING OF BOTTOM TRACK'//5(9X, 'R(KM) Z(M)') 1 //(5(I6, F8.2, F6.0))) 902 FORMAT(30H0INTERPOLATED PROFILE AT RANGE,F10.0,2H M//) 903 FORMAT(I8,F8.0,10A4) 904 FORMAT('1INPUT PROFILE AT RANGE' ,F11.3,3H KM,5X,10A4//) 914 FORMAT (3X,'DEPTH (M)',2X,'C(M/SEC)','1460',6X,'1470',6X,'1480', *6X,'1490',6X,'1500',6X,'1510',6X,'1520',6X,'1530',6X,'1540',6X, *'1550',6X,'1560',//) 915 FORMAT (1X,F11.4,F10.4,101A1) 905 FORMAT ( '0TOO MANY POINTS INTERPOLATED' ) END SUBROUTINE CHANNL(ZR,T,ZTO,ZTU) INTEGER RAYFIL, INPFIL, PRTFIL, INTFIL, DSTFIL, IGNFIL PARAMETER ( RAYFIL = 21, INPFIL = 5, PRTFIL = 6, INTFIL = 7, & DSTFIL = 37, IGNFIL = 41 ) COMMON /VELPRF/ N,Z(100),C(100) REAL*8 ARG,V,DR,CM,TG1,TG2,F,RCYCLE CALL VELCAL 20 DO 25 I=2,N IF(ZR.GT.Z(I)) GO TO 25 IZ=I F=(ZR-Z(I-1))/(Z(I)-Z(I-1)) CM=DSQRT(1.D0+T*T)*(F*C(I)+(1.D0-F)*C(I-1)) GO TO 30 25 CONTINUE ZTU=-1. ZTO=-1. RETURN 30 IN=IZ-1 31 IF(C(IN).GT.CM) GO TO 35 IN=IN-1 IF(IN.NE.0) GO TO 31 ZTO=0. GO TO 40 35 F=(CM-C(IN+1))/(C(IN)-C(IN+1)) ZTO=F*Z(IN)+(1.-F)*Z(IN+1) 40 IN=IZ 41 IF(C(IN).GT.CM) GO TO 45 IN=IN+1 IF(IN.LE.N) GO TO 41 ZTU=Z(N) RETURN 45 F=(CM-C(IN-1))/(C(IN)-C(IN-1)) ZTU=F*Z(IN)+(1.-F)*Z(IN-1) RETURN ENTRY RCALC(ZR,T,ZTO,ZTU) IF(ZTO.GT.-.5) GO TO 50 RCYCLE=1. D30 RETURN 50 IN=1 RCYCLE=0. TG2=0.D0 DO 60 I=1,N IF(Z(I).GT.ZTO.AND.Z(I).LT.ZTU) GO TO 55 IF(IN.EQ.1) GO TO 60 RCYCLE = RCYCLE + 2.D0*(ZTU-Z(I-1))/TG1 RETURN 55 ARG = CM/C(I) ARG = ARG*ARG-1.D0 TG1 = 0.D0 IF(ARG.GT.0.D0) TG1 = DSQRT(ARG) IF(I.EQ.1) GO TO 59 IF(IN.EQ.1) GO TO 56 DR = 2.D0*(Z(I)-Z(I-1))/(TG1+TG2) GO TO 57 56 DR = 1. D30 IF(TG1.NE.0.D0) DR = 2.D0*(Z(I)-ZTO)/TG1 57 RCYCLE=RCYCLE+DR 59 TG2=TG1 IN=2 60 CONTINUE RETURN ENTRY WDENS(ZR,T) IF(RCYCLE.GT.1.D10) GO TO 71 IN=IZ 65 IF(Z(IN).GT.ZR) GO TO 70 IF(Z(IN+1).GT.ZR) GO TO 75 IN=IN+1 IF(IN.GE.N) GO TO 71 GO TO 65 70 IF(Z(IN-1).LE.ZR) GO TO 74 IN=IN-1 IF(IN.GE.1) GO TO 65 71 T=0. RETURN 74 IN=IN-1 75 V=C(IN)+(C(IN+1)-C(IN))*(ZR-Z(IN))/(Z(IN+1)-Z(IN)) ARG = CM/V ARG = ARG*ARG - 1.D0 T = 1. E+10 IF(ARG.GT.0.D0) T = 1.D0/(RCYCLE*DSQRT(ARG)) RETURN END SUBROUTINE ZDIST INTEGER RAYFIL, INPFIL, PRTFIL, INTFIL, DSTFIL, IGNFIL PARAMETER ( RAYFIL = 21, INPFIL = 5, PRTFIL = 6, INTFIL = 7, & DSTFIL = 37, IGNFIL = 41 ) COMMON/RAYS/NRAY,TH0(1000),TGAM(1000),ZZ(1000),SS(1000),DSUM(1000) 1 ,TIME(1000),PHASE(1000),NBR(1000),NTU(1000),NSR(1000),NTO(1000) COMMON /START/ CT0(1000) COMMON /INFO/ ROLD, RANGE, OMEGA,ATT,ICOH,IRAN,IT2,IT3,IGNRAY,IFLMR COMMON /SCRACH/ DDD(182) COMMON /TLE/ ITITLE (18) DIMENSION LINE(51) COMMON /PIDEF/ PI, DEGREE, TWO PI EQUIVALENCE (LINE,DDD) COMMON /ABC/ NBRS, NSRS, ALIM DIMENSION IDEPTH(6) DATA IBLANK/' '/, IPLUS/'+'/, IMINUS/'-'/, ISTAR/'*'/ DATA IBEE/'B'/, IESS/'S'/, IDOT/'.'/ 3 WRITE( DSTFIL, 2) ITITLE 1 CALL CHANNL(1.,1000.,ZS,ZB) DO 5 I=1,51 5 LINE(I) = IBLANK DO 51 I=1,51,10 LINE(I) = IDOT 51 CONTINUE RKM = .001*RANGE SCALE = .2*ZB DO 75 I=1,6 IDEPTH(I) = (I-1)*SCALE + .5 75 CONTINUE WRITE( DSTFIL, 900) RKM, IDEPTH, LINE DO 8 I=1,51 8 LINE(I) = IBLANK DO 40 I=1,NRAY 35 IF(SS(I).GT.ALIM .AND.NBR(I).LE.NBRS .AND.NSR(I).LE.NSRS) GO TO 12 15 WRITE( DSTFIL, 901) I, TH0(I) 900 FORMAT('0RAY DEPTH DISTRIBUTION AT', F9.2,' KM.', 51X, 'DEPTH (M)' 1 , /, 60X, 5I10, I9, /, 2 ' NRAY THETA0 NBR NTU NSR NTO DEPTH THETA LEVEL TIME', 3 ' CMEAN ', 51A1) 901 FORMAT(1X, I5, F7.3, 10X, 'TERMINATED.') GO TO 40 12 T=TIME(I) CMEAN = RANGE/T Z=ZZ(I) LINE(1)= IESS LINE(51)= IBEE CALL CHANNL(Z,TGAM(I),ZTO,ZTU) IZP = 50.*ZTU/ZB + 1.5 IF(IZP.LT.1 .OR. IZP.GT.51) GO TO 45 LINE(IZP)= IPLUS 45 IZM = 50.*ZTO/ZB + 1.5 IF(IZM.LT.1 .OR. IZM.GT.51) GO TO 55 LINE(IZM)= IMINUS 55 IZR = 50.*Z/ZB + 1.5 IF(IZR.LT.1 .OR. IZR.GT.51) GO TO 60 LINE(IZR)= ISTAR 60 TGR = TGAM(I) STR = TGR/SQRT(1.+TGR*TGR) TH = DEGREE*ATAN(TGR) RTL = SS(I) IF(DSUM(I).NE.0) RTL = RTL*CT0(I)/(RANGE*STR*DSUM(I)) RTL = 10.*ALOG10(ABS(RTL)) 25 WRITE( DSTFIL, 902) I, TH0(I), NBR(I), NTU(I), NSR(I), NTO(I), 1 Z, TH, RTL, T, CMEAN, LINE 902 FORMAT(1X, I5, F7.3, 4I4, F6.0, 2F7.1, F9.3, F7.1, 3X, 51A1) 30 IF (IZP.LE.51.AND.IZP.GE.1) LINE(IZP)= IBLANK IF (IZM.LE.51.AND.IZM.GE.1) LINE(IZM)= IBLANK IF (IZR.LE.51.AND.IZR.GE.1) LINE(IZR)= IBLANK 40 CONTINUE 2 FORMAT ('1',20A4,/) RETURN END SUBROUTINE INITRY INTEGER RAYFIL, INPFIL, PRTFIL, INTFIL, DSTFIL, IGNFIL PARAMETER ( RAYFIL = 21, INPFIL = 5, PRTFIL = 6, INTFIL = 7, & DSTFIL = 37, IGNFIL = 41 ) DIMENSION VBP(91) COMMON /PATTRN/ SD,ITBP,DATD,SORLEV COMMON /PIDEF/ PI, DEGREE, TWO PI COMMON/RAYS/NRAY,TH0(1000),TGAM(1000),ZZ(1000),SS(1000),DSUM(1000) 1 ,TIME(1000),PHASE(1000),NBR(1000),NTU(1000),NSR(1000),NTO(1000) COMMON/START/CT0(1000) NRAY=0 OG=1000. NFAN = 0 WRITE( PRTFIL, 905) 905 FORMAT('0 THE SPECIFIED ANGULAR APERTURES FOLLOW. WITHIN EACH FAN' 1 ,' THE RAYS TRACED ARE CENTERED IN THE ANGULAR INCREMENTS.',/, 2 ' FAN DOWNLIMIT UPLIMIT INCREMENT NRAYS LEVEL PHASE') 10 NFAN = NFAN + 1 READ( INPFIL, 900 ) GAMDD, DGAMD, GAMUD, SL, PH, XCONT NRAYS = (GAMUD - GAMDD)/DGAMD + .1 WRITE( PRTFIL, 906) NFAN, GAMDD, GAMUD, DGAMD, NRAYS, SL, PH 906 FORMAT(1H0, I4, 3F10.4, I5, 2F10.2) SL = 10.**(.1*(SL+SORLEV)) C IFF=1 IF(ABS(OG-GAMDD).LT..001) GO TO 20 IF(NRAY.EQ.0) GO TO 15 C***** THE RAY WITH ZERO INTENSITY GOES BETWEEN NON-CONTIGUOUS FANS. NRAY=NRAY+1 ZZ(NRAY) = SD TIME(NRAY) = 0. NBR(NRAY) = 0 NTU(NRAY) = 0 NSR(NRAY) = 0 NTO(NRAY) = 0 PHASE(NRAY) = PH/DEGREE TH0(NRAY) = OG S1 = SIN(OG/DEGREE) TGAM(NRAY) = S1/SQRT(1. - S1*S1) SS(NRAY)=0. CT0(NRAY) = 0. 15 OG=GAMDD SOG=SIN(OG/DEGREE) C IFF=0 20 G=OG+DGAMD S=SIN(G/DEGREE) NRAY=NRAY+1 IF(NRAY.GT.1000) GO TO 101 TH0(NRAY) = G - .5*DGAMD S1 = SIN(TH0(NRAY)/DEGREE) CT0(NRAY) = SQRT(1. - S1*S1) TGAM(NRAY) = S1/CT0(NRAY) ZZ(NRAY)=SD TIME(NRAY)= 0. DSUM(NRAY) = 0. NBR(NRAY) = 0 NTU(NRAY) = 0 NSR(NRAY) = 0 NTO(NRAY) = 0 PHASE(NRAY)=PH/DEGREE SS(NRAY) = SL*ABS(S-SOG) CT0(NRAY) = CT0(NRAY)/SS(NRAY) C C IFF=1 OG=G SOG=S IF ( (G + .5*DGAMD).LT. GAMUD ) GO TO 20 IF(XCONT.NE.0.) GO TO 10 23 IF(ITBP.EQ.0) GO TO 35 IL=1 25 IH=IL+15 IF(IH.GT.91) IH = 91 READ( INPFIL, 901 ) (VBP(I),I=IL,IH) WRITE( PRTFIL, 907) (VBP(I),I=IL,IH) 907 FORMAT('0VBP ',16F5.1) IF(VBP(IH).EQ.0.) GO TO 27 IL=IL+16 GO TO 25 27 IH=IH-1 IF(IH.EQ.1) GO TO 28 IF(VBP(IH).EQ.0.) GO TO 27 28 DO 29 I=IH,90 29 VBP(I+1)=VBP(IH) DO 31 I=1,NRAY D=ABS(DEGREE*ATAN(TGAM(I))+DATD) N=D D=D-N D=(1.-D)*VBP(N+1)+D*VBP(N+2) 31 SS(I) = SS(I)*10.**(-.1*D) WRITE( PRTFIL, 902) WRITE( PRTFIL, 903) (I, TH0(I), SS(I), PHASE(I), I=1,NRAY) 35 RETURN 101 NRAY=1000 WRITE( PRTFIL, 908) 100 IF(XCONT.EQ.0.) GO TO 23 READ( INPFIL, 900 ) XCONT GO TO 100 900 FORMAT (10F8.0) 901 FORMAT(16F5.1) 902 FORMAT('0 INITIAL RAY ANGLE, SIGNAL LEVEL, AND PHASE'//) 903 FORMAT(5(I6, F8.3, F7.4, F5.2)) 908 FORMAT ('0TOO MANY RAYS INPUT, FIRST 1000 RAYS WILL BE USED') END SUBROUTINE INTSTY(INT) INTEGER RAYFIL, INPFIL, PRTFIL, INTFIL, DSTFIL, IGNFIL PARAMETER ( RAYFIL = 21, INPFIL = 5, PRTFIL = 6, INTFIL = 7, & DSTFIL = 37, IGNFIL = 41 ) COMMON /LOUDNS/ NRCD,RCD(100),TYPE1C(100),TYPE1R(100), * TYPE2(100),TYPE3(100) COMMON /TLE/ ITITLE (18) COMMON/RAYS/NRAY,TH0(1000),TGAM(1000),ZZ(1000),SS(1000),DSUM(1000) 1 ,TIME(1000),PHASE(1000),NBR(1000),NTU(1000),NSR(1000),NTO(1000) COMMON /INFO/ROLD,RANGE,OMEGA,ATT,ICOH,IRAN,IT2,IT3,IGNRAY,IFLMR COMMON /PIDEF/ PI, DEGREE, TWO PI COMMON /VELPRF/NN,Z(100),C(100) COMMON /SCRACH/ DDD(182) DIMENSION ID(1000), IHP(100), SINT(100,4), SSL(100), QUAD(100), 1 QUS(100), TYPS(100) EQUIVALENCE (QUAD,DDD) COMMON/ABC/NBRS,NSRS,ALIM DATA IBLANK/' '/, ISTAR/'* '/, ISTARS/'** '/ DATA IHP /100*0/, SINT /400*0./, ID /1000*0/, IENT/0/ SEC(T)=FVA*SQRT(1.+T*T) 3 FVA=1. IF(ATT.NE.0.) FVA = 10.**(-.0001*ATT*RANGE) CALL VELCAL WAVEL=TWOPI*C(1)/OMEGA AK=TWOPI/WAVEL FLM=1. DO 5 I=1,NRCD TYPE1C(I) = 0. TYPE1R(I) = 0. TYPE2(I) = 0. TYPE3(I) = 0. 5 QUAD(I) = 0. IF(IRAN.EQ.0.AND.IT2.EQ.0) GO TO 60 CALL CHANNL(1.,1000.,ZS,ZB) IF(IRAN.EQ.0) GO TO 35 DO 20 I=2,NRAY NDTUO = IABS(NTO(I)-NTO(I-1)) + IABS(NTU(I)-NTU(I-1)) NDSBR = IABS(NSR(I)-NSR(I-1)) + IABS(NBR(I)-NBR(I-1)) ID(I) = NDSBR + NDTUO IF(ID(I).GT.3) ID(I)=3 W=AMIN1(SS(I),SS(I-1)) IF(W.LE.ALIM .OR. NBR(I).GT.NBRS .OR. NSR(I).GT.NSRS) ID(I)=3 20 CONTINUE DO 21 I=1,NRCD SSL(I) = 0. QUS(I) = 0. TYPS(I) = 0. IHP(I) = 0. 21 CONTINUE IPR=0 INW = ID(2) DO 30 I=2,NRAY IF(I+1.LE.NRAY) INXT = ID(I+1) IF(I+2.LE.NRAY) INXTT = ID(I+2) IF(INW.GE.2) GO TO 29 IQ1=0 IQ2=0 IQ3=0 Z1=ZZ(I-1) S1=SEC(TGAM(I-1))*SS(I-1) DS=SEC(TGAM(I))*SS(I)-S1 T1=TIME(I-1) DT=TIME(I)-T1 TG1=TGAM(I-1) DTG=TGAM(I)-TG1 DZ=ZZ(I)-Z1 C********** TRY THIS FOR THAT 1.-IN-1.E6 CHANCE *********** IF(DZ.EQ.0.) GO TO 29 C************************************************************ DO 25 J=1,NRCD ZR=RCD(J) IF(ZR.GT.ZB) GO TO 25 F=(ZR-Z1)/DZ C DZ1=ZZ(I+1)-ZZ(I) C F1=(ZR-ZZ(I))/DZ1 IF(F.LT.-.5 .OR. F.GT.1.5) GO TO 25 IF(INW.NE.0 .OR. F.GE.1. .OR. F.LE.0.) GO TO 24 C IF (INW.EQ.0 .AND. F.LE.1. .AND. F.GE.0. .AND. IHP(J).EQ.0) 23,24 C***** AT ST. 23, RAYS I & I-1 BRACKET THE RECEIVER AND HAVE SAME C TURNAROUND HISTORY, SO WE HAVE A QUALITY 1 EIGENRAY. 23 IQ1=1 GO TO 26 24 IF(F.LE.1. .OR. INW.NE.0 .OR. INXT.EQ.INW .OR. IHP(J).NE.0) 1 GO TO 32 C***** AT ST. 31 RAYS I & I-1 HAVE SAME TURNAROUND HISTORY BUT DO NOT C BRACKET THE RECEIVER; BUT THEY COME CLOSE AND NO QUALITY 1 EIGENRAY C HAS BEEN FOUND BETWEEN THEM OR WILL BE FOUND BETWEEN THE NEXT PAIR. C WE ACCEPT THIS AS A QUALITY 2 EIGENRAY. NOTE THAT FOR QUALITY 2 WE C HAVE 1. < (ZR - Z(I-1))/(Z(I)-Z(I-1)) < 1.5, SO THAT WE ARE LOOKING C FOR A RECEIVER LOCATION ALMOST REACHED BY THE "ADVANCING" RAY FAN. 31 IQ2=1 GO TO 26 32 IF(F.GE.0. .OR. IPR.EQ.INW .OR. INW.NE.0) GO TO 34 C***** AT ST. 33 RAYS I & I-1 HAVE SAME TURNAROUND HISTORY AND COME C CLOSE TO BRACKETING RECEIVER; NO QUALITY 1 OR 2 EIGENRAY HAS BEEN C FOUND BETWEEN THEM AND NO QUALITY 1 EIGENRAY WAS FOUND BETWEEN THE C PREVIOUS PAIR. WE ACCEPT THIS AS A QUALITY 3 EIGENRAY. NOTE THAT C THIS RECEIVER LOCATION HAD -.5 < (ZR-Z(I-1))/(Z(I)-Z(I-1)) < 0. , C SO THAT IT WAS MISSED IN THE WAKE OF THE RAY FAN. 33 IQ3=1 GO TO 26 C***** AT ST. 34 NO EIGENRAY HAS BEEN FOUND BECAUSE RAYS I & I-1 HAVE C DIFFERENT TURNAROUND HISTORIES. 34 IQ1=0 IQ2=0 IQ3=0 GO TO 25 C***** AT ST. 26 AN EIGENRAY HAS BEEN FOUND; INTERPOLATE OR EXTRAPOLATE C IN DEPTH ITS TRAJECTORY PARAMETERS FROM THOSE OF RAYS I & I-1. 26 SL=(S1+F*DS)/(RANGE*ABS(DZ)) IF (SL.LE.0.) GO TO 25 C***** THE NEXT STATEMENT IS TO IGNORE THIS EIGENRAY IF THE RAYS ARE C VERY CLOSE TOGETHER NEAR THIS RECEIVER; I.E., IF AT A CAUSTIC. IF (ABS(DZ).LT..001) GO TO 25 TG=TG1+F*DTG IF (IQ1.EQ.1) IHP(J)=1 TH=DEGREE*ATAN(TG) IF(INW.GE.1.AND.INXT.GE.1) GO TO 25 IF(IFLMR.NE.0) FLM=2.*SIN(ZR*AK*TG/SQRT(1.+TG*TG))**2 T=T1+F*DT IF (IQ1.EQ.1 .AND. SINT(J,1).EQ.0.) SINT(J,4)=FLM IF (IQ1.EQ.1 .AND. SINT(J,1).EQ.0.) SINT(J,1)=SL IF (IQ2.EQ.1 .AND. SINT(J,1).EQ.0. .AND. 1 SINT(J,2).EQ.0. .AND. SINT(J,3).EQ.0.) SINT(J,4)=FLM IF (IQ2.EQ.1 .AND. SINT(J,1).EQ.0. .AND. 1 SINT(J,2).EQ.0. .AND. SINT(J,3).EQ.0.) SINT(J,2)=SL IF (IQ3.EQ.1 .AND. SINT(J,1).EQ.0. .AND. 1 SINT(J,2).EQ.0. .AND. SINT(J,3).EQ.0.) SINT(J,4)=FLM IF (IQ3.EQ.1 .AND. SINT(J,1).EQ.0. .AND. 1 SINT(J,2).EQ.0. .AND. SINT(J,3).EQ.0.) SINT(J,3)=SL IF(IGNRAY.EQ.0) GO TO 22 IF (IENT.NE.0) GO TO 41 WRITE( IGNFIL, 1899) ITITLE 1899 FORMAT ('1 EIGENRAY SET ',20A4) WRITE( IGNFIL, 899) IENT = 1 41 IF(SL.LT.ALIM) GO TO 22 PR=10.*ALOG10(SL) IF (IQ3.EQ.1) IQUAL=3 IF (IQ2.EQ.1) IQUAL=2 IF (IQ1.EQ.1) IQUAL=1 IREM=IBLANK IF (IQ3.EQ.1 .AND. SINT(J,3).NE.0.) IREM=ISTAR IF (SINT(J,3).NE.0. .AND. IQ1.EQ.1) IREM=ISTARS RMKM = .001*RANGE WRITE( IGNFIL, 900) I, TH0(I), NBR(I), NTU(I), NSR(I), NTO(I), 1 RMKM, ZR, TH, T, PR, IQUAL, IREM 899 FORMAT('0NRAY THETA0 NBR NTU NSR NTO RANGE(km) DEPTH(m)' 1 ,' THETA(deg) TIME(sec)', 4X,'SL(dB)',7X,'IQUAL',4X,'REM') 900 FORMAT(I5, F7.3, 4I5, F15.3, 2F10.2, F13.5, F10.1, 6X, I5, 5X, A2) 1878 FORMAT(I5, F9.3, 4I5, 3F8.2, F12.5, F10.1) 22 FLM=SINT(J,4) IF (SINT(J,1).NE.0.) SL=SINT(J,1) IF (SINT(J,2).NE.0. .AND. SINT(J,1).EQ.0.) SL=SINT(J,2) IF (SINT(J,3).NE.0. .AND. SINT(J,1).EQ.0.) SL=SINT(J,3) IF (SINT(J,3).EQ.0.) GO TO 27 IF (IQ1.NE.1 .OR. SINT(J,1).EQ.0.) GO TO 72 71 TYPE1R(J)=SSL(J) IF (ICOH.EQ.0) GO TO 27 TYPE1C(J)=TYPS(J) QUAD(J)=QUS(J) GO TO 27 72 SSL(J)=TYPE1R(J) 27 TYPE1R(J)=TYPE1R(J)+SL*FLM IF(ICOH.EQ.0) GO TO 25 PR=SQRT(SL) P=OMEGA*T+PHASE(I) IF (IQ3.EQ.1 .AND. IQ1.EQ.0) QUS(J)=QUAD(J) QUAD(J)=QUAD(J)+PR*SIN(P) IF (IQ3.EQ.1 .AND. IQ1.EQ.0) TYPS(J)=TYPE1C(J) TYPE1C(J)=TYPE1C(J)+PR*COS(P) 25 IF (INXT.GE.1 .OR. IPR.NE.INW .AND. INW.NE.INXT)IHP(J)=0 29 IPR=INW INW=INXT C IF (INW.EQ.0 .OR. INW.LT.1 .AND. INXTT.LT.1) GO TO 30 IF(INW.EQ.0) GO TO 30 DO 28 LL=1,NRCD SINT(LL,1)=0. SINT(LL,2)=0. SINT(LL,3)=0. SINT(LL,4)=0. 28 CONTINUE 30 CONTINUE 35 DO 37 I=1,NRCD TYPE1C(I)=TYPE1C(I)*TYPE1C(I) + QUAD(I)*QUAD(I) 37 CONTINUE IF(IT2.EQ.0) GO TO 60 SW = 0. SDZ = 0. DO 38 I=2,NRAY W = AMIN1(SS(I),SS(I-1)) SDZ = SDZ + W*ABS(ZZ(I) - ZZ(I-1)) SW = SW + W 38 CONTINUE DZBAR = WAVEL IF(SW.GT.0.) DZBAR = AMAX1(DZBAR,SDZ/SW) DZM = ZB/SQRT(FLOAT(NRAY)) IF(DZBAR.GT.DZM) DZBAR = DZM C***** ASC WORD FORMAT LIMITS MAGNITUDE OF EXPONENTIAL ARGUMENT. DZBAR = AMAX1(DZBAR,.016*ZB) DO 385 J=1,NRCD IF(RCD(J).LE.ZB) QUAD(J) = EXP(-RCD(J)/DZBAR) 385 CONTINUE F=2.*RANGE*DZBAR EB=EXP(-ZB/DZBAR) DO 50 I=1,NRAY IF(ABS(SS(I)).LT.ALIM) GO TO 50 EZ=EXP(-ZZ(I)/DZBAR) SL=SEC(TGAM(I))*SS(I)/F SGK=TGAM(I)*AK/SQRT(1.+TGAM(I)*TGAM(I)) SL=SL/(1.-.5*(EB/EZ+EZ)) DO 40 J=1,NRCD IF(RCD(J).GT.ZB) GO TO 50 ER=QUAD(J) IF(IFLMR.NE.0) FLM=2.*SIN(RCD(J)*SGK)**2 IF(ER.GT.EZ) GO TO 39 S=SL*ER/EZ GO TO 40 39 S=SL*EZ/ER 40 TYPE2(J)=TYPE2(J)+S*FLM 50 CONTINUE 60 IF(IT3.EQ.0) RETURN DO 70 I=1,NRAY CALL CHANNL(ZZ(I),TGAM(I),ZTO,ZTU) CALL RCALC(ZZ(I),TGAM(I),ZTO,ZTU) IF(IFLMR.NE.0) CALL WDENS(ZZ(I),CTRR) SL=FVA*SS(I)/RANGE DO 65 J=1,NRCD IF(RCD(J).GT.ZTU) GO TO 70 IF(RCD(J).LT.ZTO) GO TO 65 CALL WDENS(RCD(J),S) IF(IFLMR.EQ.0) GO TO 64 TG=TGAM(I)*CTRR/S FLM=2.*SIN(AK*RCD(J)*TG/SQRT(1.+TG*TG))**2 C************************** VAX FLP OVERFLOW W/TG=4.5E19. S=S*FLM 64 TYPE3(J)=TYPE3(J)+S*SL 65 CONTINUE 70 CONTINUE RETURN END DOUBLE PRECISION FUNCTION DELB(DR, TGRP, TGRPN) REAL*8 DR, T0, T1, TGRP, TGRPN, T0T1 T0 = TGRP T1 = TGRPN T0T1 = T0*T1 IF(T0T1.NE.0.D0) DELB = -(1. + 1./T0T1) IF(T0.EQ.0.D0) DELB = 1./(T1*T1) - 1. IF(T1.EQ.0.D0) DELB = 1./(T0*T0) - 1. DELB = DR*DELB RETURN END SUBROUTINE CONNCT IMPLICIT REAL*8(A-H,O-Z) COMMON /PROF/ R1,N1,Z1(50),V1(50),R2,N2,Z2(50),V2(50) COMMON /TRIANG/ AP(100,2),BP(100,2),AL(100),BL(100),ZZERO(100), * RZERO(100),AA(100),BB(100),SST(100),CCT(100),NTRI REAL*4 R1, Z1, V1, R2, Z2, V2 NA=1 NB=1 NTRI=0 IADD=0 10 NTRI=NTRI+1 IF(NB.EQ.N2) GO TO 12 IF (IADD.EQ.1.OR.NA.EQ.N1) GO TO 30 12 RZERO(NTRI)=R2 ZZERO(NTRI)=Z2(NB) AL(NTRI)=Z1(NA+1) BL(NTRI)=(Z2(NB)-Z1(NA+1))/(R2-R1) AA(NTRI)=1.D0/V2(NB)**2 BZ=(1.D0/V1(NA)**2-1.D0/V1(NA+1)**2)/(Z1(NA)-Z1(NA+1)) BR=(1.D0/V2(NB)**2-1.D0/V1(NA)**2-BZ*(Z2(NB)-Z1(NA)))/(R2-R1) NA=NA+1 IADD=1 GO TO 40 30 RZERO(NTRI)=R1 ZZERO(NTRI)=Z1(NA) AL(NTRI)=Z1(NA) BL(NTRI)=(Z1(NA)-Z2(NB+1))/(R1-R2) AA(NTRI)=1.D0/V1(NA)**2 BZ=(1.D0/V2(NB)**2-1.D0/V2(NB+1)**2)/(Z2(NB)-Z2(NB+1)) BR=(1.D0/V1(NA)**2-1.D0/V2(NB)**2-BZ*(Z1(NA)-Z2(NB)))/(R1-R2) NB=NB+1 IADD=0 40 BB(NTRI) = DSIGN(DSQRT(BZ*BZ + BR*BR), BZ) C IF(BB(NTRI).EQ.0.) GO TO 41 IF(DABS(BB(NTRI)).LT.1.D-20) GO TO 41 SST(NTRI)=BR/BB(NTRI) CCT(NTRI)=BZ/BB(NTRI) GO TO 42 41 CCT(NTRI)=1.D0 SST(NTRI)=0. 42 AL(NTRI)=AL(NTRI)-BL(NTRI)*R1 TANT=SST(NTRI)/DMAX1(CCT(NTRI),1.D-30) BP(NTRI,2)=(BL(NTRI)+TANT)/(1.-BL(NTRI)*TANT) B=BL(NTRI-1) IF(NTRI.EQ.1) B=0 BP(NTRI,1)=(B+TANT)/(1.D0-B*TANT) AP(NTRI,1) = 0. AP(NTRI,2) = 0. IF(NA.LT.N1.OR.NB.LT.N2) GO TO 10 RETURN END SUBROUTINE RAYOUT INTEGER RAYFIL, INPFIL, PRTFIL, INTFIL, DSTFIL, IGNFIL PARAMETER ( RAYFIL = 21, INPFIL = 5, PRTFIL = 6, INTFIL = 7, & DSTFIL = 37, IGNFIL = 41 ) COMMON/RAYS/NRAY,TH0(1000),TGAM(1000),ZZ(1000),SS(1000),DSUM(1000) 1 ,TIME(1000),PHASE(1000),NBR(1000),NTU(1000),NSR(1000),NTO(1000) COMMON /INFO/ROLD,RANGE,OMEGA,ATT,ICOH,IRAN,IT2,IT3,IGNRAY,IFLMR COMMON /TRIANG/ AP(100,2),BP(100,2),AL(100),BL(100),ZZERO(100), * RZERO(100),AA(100),BB(100),SST(100),CCT(100),NTRI REAL*8 AP,BP,AL,BL,ZZERO,RZERO,AA,BB,SST,CCT COMMON /TLE/ ITITLE (18) COMMON/ABC/NBRS, NSRS, ALIM COMMON /RANST/ RDUM,IREC,IFSK,SDUM DATA IFSK/-1/ DATA IENT/0/ ZB = AL(NTRI) + BL(NTRI)*RANGE C IF (IENT.EQ.1) GO TO 1 C WRITE( RAYFIL )ITITLE, ALIM, NBRS, NSRS C IENT=1 C IREC=1 C IFSK=IFSK+1 1 WRITE( RAYFIL ) NRAY, RANGE, ZB, & (TGAM(I),I=1,NRAY), (ZZ(I),I=1,NRAY), 1 (SS(I),I=1,NRAY), (TIME(I),I=1,NRAY), (PHASE(I),I=1,NRAY), 2 (NBR(I),I=1,NRAY), (NTU(I),I=1,NRAY), (NSR(I),I=1,NRAY), 3 (NTO(I),I=1,NRAY) IREC=IREC+1 RETURN END SUBROUTINE VELCAL COMMON /INFO/ RDUM, RANGE, OMEGA, IDUMM(7) COMMON /VELPRF/ N,Z(100),C(100) COMMON /TRIANG/ AP(100,2),BP(100,2),AL(100),BL(100),ZZERO(100), * RZERO(100),AA(100),BB(100),SST(100),CCT(100),NTRI REAL*8 AP,BP,AL,BL,ZZERO,RZERO,AA,BB,SST,CCT COMMON/ABC/NBRS,NSRS,ALIM DATA ROLD/-1. E+20/ 3 IF(ABS(RANGE-ROLD).LT.1.) GO TO 20 ROLD=RANGE N=NTRI+1 DO 10 I=1,N IF(I.EQ.1) GO TO 1 J=I-1 ZZ=AL(J)+BL(J)*RANGE GO TO 2 1 J=1 ZZ=0. 2 ZP=CCT(J)*(ZZ-ZZERO(J))+SST(J)*(RANGE-RZERO(J)) C(I)=1.D0/DSQRT(AA(J)+BB(J)*ZP) 10 Z(I)=ZZ 20 RETURN END SUBROUTINE ADVNCE IMPLICIT REAL*8(A-H,O-Z) INTEGER RAYFIL, INPFIL, PRTFIL, INTFIL, DSTFIL, IGNFIL PARAMETER ( RAYFIL = 21, INPFIL = 5, PRTFIL = 6, INTFIL = 7, & DSTFIL = 37, IGNFIL = 41 ) COMMON /TRIANG/ AP(100,2),BP(100,2),AL(100),BL(100),ZZERO(100), * RZERO(100), AA(100), BB(100), SST(100), CCT(100), NTRI COMMON/RAYS/NRAY,TH0(1000),TGAM(1000),ZZ(1000),SS(1000),DSUM(1000) 1 ,TIME(1000),PHASE(1000),NBR(1000),NTU(1000),NSR(1000),NTO(1000) COMMON /START/ CT0(1000) COMMON /INFO/RSTART,RMAX,OMEGA,ATT,ICOH,IRAN,IT2,IT3,IGNRAY,IFLMR COMMON /PIDEF/ PI, DEGREE, TWO PI COMMON /MIRROR/ SL, BRLT(200), BPST(200) COMMON/ABC/ NBRS, NSRS, ALIM DIMENSION RPNEW(4), RNEW(4), ZNEW(4), ZPNEW(4), D(7), RTNEXT(1000) REAL*4 TH0, TGAM, ZZ, SS, TIME, PHASE, RSTART, RMAX, OMEGA, ATT REAL*4 DEGREE, PI, TWO PI, BRLT, BPST, ALIM, DSUM, CT0, SL DATA ONE MM /.001D0/, FIVE MM /.005D0/, RTNEXT /1000*1.D0/ DATA DLIMIT /1.D-30/ LOGICAL ONUP,ONLW DELT(DR,T,U) = DABS(DR)*DSQRT(CMIS)*(1.D0+(T*(T+U)+U*U)/3.D0) C LOOP OVER ALL THE RAYS DO 100 IRAY=1,NRAY IHIST = 0 IF(NBRS.EQ.100000 .AND. NSRS.EQ.100000) GO TO 10 IF(NBR(IRAY).GT.NBRS .OR. NSR(IRAY).GT.NSRS) SS(IRAY) = .1*ALIM 10 IF(ABS(SS(IRAY)).LT.ALIM) GO TO 100 ZR = ZZ(IRAY) RR = RSTART TGR = TGAM(IRAY) TRAY = TIME(IRAY) DRDTH0 = DSUM(IRAY) DO 20 I=1,NTRI C FIND CORRECT LAYER P = ZR-AL(I)-BL(I)*RR IF(P.GT.ONE MM) GO TO 20 NTRR = I IF(DABS(P).GT.ONE MM .OR. TGR.GT.-BL(I)) GO TO 30 20 CONTINUE GO TO 86 30 CT = CCT(NTRR) C IF(ITLC.NE.9) GO TO 301 C STR = TGR/DSQRT(1.+TGR*TGR) C RTL = 0. C IF(RR.GT.0. .AND. DRDTH0.NE.0.) C 1 RTL = 10.*DLOG10(DABS(CT0(IRAY)*SS(IRAY)/(RR*DRDTH0*STR))) C WRITE( PRTFIL, 901) IRAY, NTRR, IHIST, RR, ZR, TGR, TRAY, C 1 RTL, PHASE(IRAY), NBR(IRAY), NTU(IRAY), NSR(IRAY), NTO(IRAY) 301 ST = SST(NTRR) Z0 = ZZERO(NTRR) C TRANSFORM TO PRIMED COORDINATES R0 = RZERO(NTRR) ZRP = CT*(ZR-Z0)+ST*(RR-R0) RRP = CT*(RR-R0)-ST*(ZR-Z0) CIS = AA(NTRR)+BB(NTRR)*ZRP TANT = ST/DMAX1(CT,DLIMIT) TGRP = (TANT - TGR)/(1. + TANT*TGR) CMIS = CIS/(1.+TGRP*TGRP) C CALCULATE PARABOLIC RAY PATH. ALPHA IS THE CURVATURE ALPHA = .25*BB(NTRR)/CMIS TWO ALP = 2.*ALPHA C FIND INTERSECTIONS WITH UPPER AND LOWER LAYER BOUNDARIES DO 40 I=1,2 C = AP(NTRR,I) + BP(NTRR,I)*RRP - ZRP P = BP(NTRR,I) - TGRP IF(DABS(ALPHA).GT.1.D-25) GO TO 303 C RAY PATH IS LINEAR FOR THESE STATEMENTS DRP = -C/P K = 2*I-1 RPNEW(K) = RRP+DRP ZPNEW(K) = ZRP + TGRP*DRP ZNEW(K) = CT*ZPNEW(K) - ST*RPNEW(K) + Z0 RNEW(K) = CT*RPNEW(K) + ST*ZPNEW(K) + R0 ZNEW(K+1) = 1.D+30 ZPNEW(K+1) = 1.D+30 RNEW(K+1) = 1.D+30 RPNEW(K+1) = 1.D+30 GO TO 40 C SOLVE QUADRATIC EQUATION FOR PARABOLIC RAY PATH 303 DISC = P*P+4.*ALPHA*C IF(DISC.LT.0.) GO TO 32 DISC = DSQRT(DISC) DO 31 J=1,2 K = 2*I+J-2 S = 2*J-3 F = P+S*DISC DRP = F/TWO ALP C ITERATION TO IMPROVE ACCURACY FOR SMALL CURVATURE RAYS IF(DABS(F).LT..1*DABS(P)) DRP=(ALPHA*((ALPHA*DRP*DRP-C)/P)**2-C)/P RPNEW(K) = RRP+DRP ZPNEW(K) = ZRP+DRP*(TGRP+ALPHA*DRP) ZNEW(K) = CT*ZPNEW(K) - ST*RPNEW(K) + Z0 31 RNEW(K) = CT*RPNEW(K) + ST*ZPNEW(K) + R0 GO TO 40 32 DO 33 J=1,2 K = 2*I + J-2 33 RNEW(K) = -1.D+30 40 CONTINUE ONLW= DABS(ZR-AL(NTRR)-BL(NTRR)*RR).LT.FIVEMM.AND.TGR.GT.-BL(NTRR) C SELECT CORRECT INTERSECTION AS NEXT POSITION IF(NTRR.EQ.1) GO TO 405 ONUP = DABS(AL(NTRR-1)+BL(NTRR-1)*RR-ZR).LT.FIVE MM * .AND. TGR.LT.-BL(NTRR-1) GO TO 41 405 ONUP = DABS(ZR).LT.FIVE MM .AND. TGR.LT.0. 41 IUP = -1 IDN = -1 IF(.NOT.ONUP) GO TO 411 IUP = 1 IF(DABS(RR-RNEW(1)).LT.DABS(RR-RNEW(2))) IUP = 2 IF(RNEW(IUP).LT.RR .OR. RNEW(IUP).GT.RMAX) IUP = 0 411 IF(.NOT.ONLW) GO TO 412 IDN = 3 IF(DABS(RR-RNEW(3)).LT.DABS(RR-RNEW(4))) IDN = 4 IF(RNEW(IDN).LT.RR .OR. RNEW(IDN).GT.RMAX) IDN = 0 412 IF(IUP.GE.0) GO TO 414 IUP = 0 RTRY = RMAX DO 413 I=1,2 IF(RNEW(I).LT.RR .OR. RNEW(I).GT.RTRY) GO TO 413 RTRY = RNEW(I) IUP = I 413 CONTINUE 414 IF(IDN.GE.0) GO TO 416 IDN = 0 RTRY = RMAX DO 415 I=3,4 IF(RNEW(I).LT.RR .OR. RNEW(I).GT.RTRY) GO TO 415 IDN = I RTRY = RNEW(I) 415 CONTINUE 416 IF(IDN.EQ.0) GO TO 420 IF(IUP.NE.0) GO TO 418 417 ITRY = IDN RTRY = RNEW(IDN) GO TO 429 418 IF(RNEW(IUP).GE.RNEW(IDN)) GO TO 417 419 ITRY = IUP RTRY = RNEW(IUP) GO TO 429 420 IF(IUP.EQ.0) GO TO 50 GO TO 419 429 DELR = RPNEW(ITRY) - RRP TGRPN = TGRP + TWO ALP*DELR TRAY = TRAY + DELT(DELR, TGRP, TGRPN) DRDTH0 = DRDTH0 + DELB(DELR, TGRP, TGRPN)*TGRP*RTNEXT(IRAY)/CT TRP = TGRP IF(DABS(TGRP).LT.DLIMIT) TRP = DSIGN(DLIMIT,TGRP) TRPN = TGRPN IF(DABS(TGRPN).LT.DLIMIT) TRPN = DSIGN(DLIMIT,TGRPN) RTNEXT(IRAY) = RTNEXT(IRAY)*TRP/TRPN TGRN = (TANT - TGRPN)/(1. + TGRPN*TANT) RR = RTRY ZR = ZNEW(ITRY) C CHECK FOR SURFACE AND BOTTOM HITS, TURNOVERS AND TURNUNDERS IF(TGR*TGRN.GE.0.) GO TO 44 IF(TGR.GT.0.) GO TO 43 NTU(IRAY) = NTU(IRAY) + 1 GO TO 44 43 NTO(IRAY) = NTO(IRAY) + 1 44 IF(ITRY.GT.2) GO TO 45 NTRR = NTRR-1 IF(NTRR.GT.0) GO TO 46 NTRR = 1 SS(IRAY) = SS(IRAY)*SL PHASE(IRAY) = PHASE(IRAY) + PI NSR(IRAY) = NSR(IRAY) + 1 TGR = -TGRN GO TO 30 45 NTRR = NTRR + 1 IF(NTRR.LE.NTRI) GO TO 46 NTRR = NTRI NBR(IRAY) = NBR(IRAY) + 1 TGRAZE = (TGRN + BL(NTRI))/(1.D0 - TGRN*BL(NTRI)) S = 100.D0*DABS(TGRAZE) + 1.D0 N = S S = S-N IF(N.GT.199) N = 199 SS(IRAY) = SS(IRAY)*((1.-S)*BRLT(N) + S*BRLT(N + 1)) PHASE(IRAY) = PHASE(IRAY) + ((1.-S)*BPST(N) + S*BPST(N + 1)) TGR = -(TGRAZE + BL(NTRI))/(1. - TGRAZE*BL(NTRI)) IF(TGR.LT.-BL(NTRI)) GO TO 86 GO TO 30 46 TGR = TGRN GO TO 30 C CALCULATE RAY INTERSECTION WITH VERTICAL BOUNDARY, R = RMAX 50 IF(DABS(ST).LT.2.D-4) GO TO 60 BV = -CT/ST AV = (RMAX-R0)/ST C = AV + RRP*BV - ZRP P = BV-TGRP IF(DABS(ALPHA).GT.1.D-25) GO TO 501 C LINEAR RAYS RPNEW(1) = RRP-C/P RPNEW(2) = 1.D+30 GO TO 505 C PARABOLIC RAYS 501 DISC = P*P+4.D0*ALPHA*C IF(DISC.LT.0.) GO TO 86 DISC = DSQRT(DISC) S = 1. DO 503 I=1,2 F = P + S*DISC DRP = F/TWO ALP IF(DABS(F).LT..1*DABS(P)) DRP=(ALPHA*((ALPHA*DRP*DRP-C)/P)**2-C)/P RPNEW(I) = RRP+DRP 503 S = -S 505 DO 51 I=1,2 ZPNEW(I) = ZRP + (RPNEW(I)-RRP)*(TGRP + ALPHA*(RPNEW(I)-RRP)) 51 ZNEW(I) = CT*ZPNEW(I) - ST*RPNEW(I) + Z0 I = 2 C SELECT CORRECT INTERSECTION IF(DABS(ZNEW(1)-ZR).LT.DABS(ZNEW(2)-ZR)) I = 1 DELR = RPNEW(I) - RRP TGRPN = TGRP + TWO ALP*DELR TRAY = TRAY + DELT(DELR, TGRP, TGRPN) DRDTH0 = DRDTH0 + DELB(DELR, TGRP, TGRPN)*TGRP*RTNEXT(IRAY)/CT TRP = TGRP IF(DABS(TGRP).LT.DLIMIT) TRP = DSIGN(DLIMIT,TGRP) TRPN = TGRPN IF(DABS(TGRPN).LT.DLIMIT) TRPN = DSIGN(DLIMIT,TGRPN) RTNEXT(IRAY) = RTNEXT(IRAY)*TRP/TRPN TGRN = (TANT - TGRPN)/(1. + TGRPN*TANT) C CHECK FOR TURNOVERS AND TURNUNDERS IF(TGR*TGRN.GE.0.) GO TO 54 IF(TGR.LT.0.) GO TO 53 NTO(IRAY) = NTO(IRAY) + 1 GO TO 54 53 NTU(IRAY) = NTU(IRAY) + 1 54 TGAM(IRAY) = TGRN IF(NTRR.LE.1) GO TO 545 C CHECK THAT RAY IS WITHIN THE PROPER LIMITS IF(ZNEW(I).LT.AL(NTRR-1) + BL(NTRR-1)*RMAX-.1D0) GO TO 86 545 IF(ZNEW(I).GT.AL(NTRR) + BL(NTRR)*RMAX + .1D0) GO TO 86 C THROW IN VOLUME ATTENUATION IF(ATT.NE.0) 1 SS(IRAY)=SS(IRAY)*10.**(-.0001*ATT*(TRAY-TIME(IRAY))/DSQRT(CIS)) 55 ZZ(IRAY) = ZNEW(I) TIME(IRAY) = TRAY DSUM(IRAY) = DRDTH0 GO TO 100 60 IHIST = 60 C WHEN SMALL OR ZERO HORIZONTAL GRADIENTS, EQUATION FOR SECTION C 50 IS SINGULAR. THIS IS AN ITERATIVE SOLUTION FOR SUCH A CASE C GUESS VALUE FOR ITERATION ZNW = ZR-TGR*(RMAX-RR) 66 DO 70 I=1,7 RPNW = CT*(RMAX-R0)-ST*(ZNW-Z0) DRP = RPNW-RRP ZPNW = ZRP + DRP*(TGRP + ALPHA*DRP) ZNW2 = CT*ZPNW-ST*RPNW+Z0 D(I) = DABS(ZNW2-ZNW) ZNW = ZNW2 70 CONTINUE C CHECK FOR CONVERGENCE IF(D(7).LT..01D0) GO TO 68 IF(D(7).LT..5*D(3)) GO TO 66 IHIST = 66 GO TO 86 68 RPNEW(1) = RPNW RPNEW(2) = RPNW GO TO 505 C PRINT ERROR MESSAGE 86 SS(IRAY) = 0. WRITE( PRTFIL, 986) IRAY, IHIST, RR, ZR, TGR, RMAX, ZNEW(I), TGRN GO TO 55 100 CONTINUE RETURN 901 FORMAT(' RAY', I4, ' IN LAYER', I4, ' IHIST=', I3, 1 ' R=', F10.2, ' Z=', F8.2, ' TG=', F10.5, 2 ' T=', F10.4, ' L=', F6.1, ' P=', F8.3, ' NCTR=', 4I4) 986 FORMAT(' RAY', I5, ' TERMINATED; IHIST=', I2, ' RR, ZR, TGR =', 1 3D11.3, ' NEW R, Z, TG =', 3D11.3) END SUBROUTINE INTPRT INTEGER RAYFIL, INPFIL, PRTFIL, INTFIL, DSTFIL, IGNFIL PARAMETER ( RAYFIL = 21, INPFIL = 5, PRTFIL = 6, INTFIL = 7, & DSTFIL = 37, IGNFIL = 41 ) COMMON /INFO/ RSTART,RMAX,OMEGA,ATT,ICOH,IRAN,IT2,IT3,IGNRAY,IFLMR COMMON/LIMITS/ R1 , DR , R2 COMMON/ABC/NBRS,NSRS,ALIM COMMON /TLE/ ITITLE (18) DATA IENT/0/ DIMENSION IS(4), ITT(4) EQUIVALENCE (ICOH,ITT) COMMON /LOUDNS/ NRCD,RCD(100),XINT(100,4) COMMON /TLPLOT/ RKM(1000), XTL(1000,8), ITYP(1000), IRCD(1000) DATA IS/1HC, 1HR, 1H2, 1H3/, NR/0/ DATA IBLANK/' '/ 40 IF(IENT.EQ.1) GO TO 41 WRITE( PRTFIL, 900) ITITLE 900 FORMAT ('1',20A4) IENT=1 WRITE( PRTFIL, 905) (RCD(I),I=1,NRCD) WRITE( PRTFIL, 908) 905 FORMAT('0RECEIVED LEVEL VS RANGE',19X,'AT DEPTHS'/ * '0 R(KM) TYPE',13F9.3,/,15X,13F9.3) 908 FORMAT(1X) 41 NR = NR + 1 RKM(NR) = .001*RMAX NPLOT=0 NPR = 0 DO 50 I=1,4 IF(ITT(I).EQ.0) GO TO 50 DO 45 J=1,NRCD XINT(J,I) = 10.*ALOG10(AMAX1(ALIM ,XINT(J,I))) NPLOT = NPLOT + 1 XTL(NR,NPLOT) = ABS(XINT(J,I)) ITYP(NPLOT) = I-1 IRCD(NPLOT) = RCD(J) + .5 45 CONTINUE IF(NPR.EQ.0) GO TO 46 C WRITE(INTFIL, 909) RKM(NR), (XINT(J,I),J=1,NRCD) WRITE( PRTFIL, 906) IS(I),(XINT(J,I),J=1,NRCD) 906 FORMAT(13X,A1,1X,13F9.1,/,15X,13F9.1) GO TO 50 46 NPR = 1 C WRITE(INTFIL, 909) RKM(NR), (XINT(J,I),J=1,NRCD) WRITE( PRTFIL, 907) RKM(NR),IS(I),(XINT(J,I),J=1,NRCD) 907 FORMAT(F9.3,4X,A1,1X,13F9.1,/,15X,13F9.1) 909 FORMAT(10F8.2) 50 CONTINUE RETURN END