C PEDERSEN-GORDON "N-LAYER" NORMAL-MODE PROGRAM C Double Precision Version (Real*8) C REFERENCE: N.O.S.C. TECH REPORT 393 (1979) BY D F GORDON. C THIS IS THE MAIN PROGRAM C LINK STATEMENT IS in NLAYER.MAK C Parameter Running Index Dimension Particular Value C C Layers L LDN LX C Modes M MDN MX C Receivers J JDN JX C Ranges K KDN KX C KEYS: K1 = 1, Output individual Depth Function values; C K6: 2, small quantity of mode data output; 1, medium quantity C 0, large quantity C K8 = 1: read value for TLIM (default is 25) PROGRAM NLAYER include 'param.inc' real*8 F(JDN), FI(JDN) real*8 GAMMAI(LDN), BLPK(LDN) real*8 DPK(LDN), CI(LDN), Thick(LDN) real*8 C(LDN), GAMMA(LDN) real*8 RECVRS(JDN), source real ANORM(JDN, MDN), RNORM(JDN, MDN) character*40 caption /'NLAYER '/ common / functions / wuu(jdn, mdn) INCLUDE 'COMMON.inc' write (*,*) ' Enter the Filename (up to 8 characters):' READ (*,'(a)') CAPTION(1:8) OPEN (JIN, file = CAPTION(1:8)//'.dat', STATUS = 'OLD') C READ IN PARAMETERS READ (JIN, '(9i4)') K1, K2, K3, K4, K5, K6, K7, K8, K9 write (*, '(9i4)') K1, K2, K3, K4, K5, K6, K7, K8, K9 if (k8 .eq. 1) then read (jin,*) tlim else TLIM = 25 end if SLIM = - TLIM**(2.D0 /3) EXPONT = EXP((TLIM + TLIM) / 3.) READ (JIN,'(i2, f10.1)') LX, FREQ write (*,'(i2, f10.1)') LX, FREQ IF (LX .EQ. 0) STOP 130 OPEN (JOUT, file = CAPTION(1:8)//'.out') WRITE (JOUT,135) CAPTION 135 FORMAT (' Filename: ', A) WRITE (JOUT,165) K1, K2, K3, K4, K5, K6, K7, K8, K9 165 FORMAT(' KEYS:', 9I8) WRITE (CAPTION(11:15), 105) NINT (FREQ) 105 FORMAT (I5) WRITE (JOUT,125) LX, FREQ 125 FORMAT (' No. of layers:', I3, ' Freq:', F10.1, 3H Hz) READ (JIN, '(8f10.4)') (Z(L), L=1,LX) WRITE (JOUT,305) (Z(L), L=1,LX) WRITE (*,305) (Z(L), L=1,LX) 305 FORMAT ('~', f13.5, 8F14.5) READ (JIN, '(8f10.4)') (C(L), L = 1,LX) WRITE (JOUT,305) (C(L), L = 1,LX) WRITE (*,305) (C(L), L = 1,LX) READ (JIN, '(8f10.4)') (CB(L), L = 1,LX) WRITE (JOUT,305) (CB(L), L = 1,LX) READ (JIN, '(8f10.4)') (GAMMA(L), L = 1,LX) WRITE (JOUT,305) (GAMMA(L), L = 1,LX) READ (JIN, '(8f10.4)') (DPK(L), L = 1,LX) WRITE (JOUT,305) (DPK(L), L = 1,LX) READ (JIN, '(8f10.4)') (BLPK(L), L = 1,LX) WRITE (JOUT,305) (BLPK(L), L = 1,LX) READ (JIN, '(8f10.4)') (RHO(L), L = 1,LX) WRITE (JOUT,305) (RHO(L), L = 1,LX) C ABSORPTION: FK = FREQ /1000. fsq = fk**2 atten = 0.1 * fsq /(1 + fsq) + 40 * fsq /(4100 + fsq) atten = atten * 1.0936 ! convert from dB/kyd to dB/km WRITE (JOUT,205) atten 205 FORMAT (' Assumed rate of absorption:', 1P e9.2, ' dB/Km'/) GAMMAI(LX) = 0. C COMPLETE PROFILE GAMMAI(LX) = 0 DO 33 L = 1, LX IF (L .LT. LX) Thick(L) = Z(L+1) - Z(L) IF (RHO(L) .EQ. 0.) RHO(L) = 1.02 IF (CB(L) .NE. 0.) GO TO 31 IF (L .LT. LX) CB(L) = C(L+1) 31 IF (C(L) .NE. 0.) GO TO 32 C(L) = CB(L-1) 32 IF (DPK(L) .GE. 0.) GO TO 34 CI(L) = CBI(L-1) GO TO 36 34 CI(L) = 0. IF (DPK(L) .EQ. 0.) GO TO 36 T = 27287.52708 / DPK(L) CI(L) = T - SQRT((T - C(L)) * (T + C(L))) 36 CBI(L) = 0. IF (GAMMA(L) .NE. 0.) GO TO 38 C **BOTH SOUND SPEEDS GIVEN IF (BLPK(L) .LE. 0.) GO TO 37 T = 27287.52708 / BLPK(L) CBI(L) = T - SQRT ((T - CB(L)) * (T + CB(L))) 37 T = C(L) * (C(L)**2 - 3. * CI(L)**2) TI = CI(L) * (3. * C(L)**2 - CI(L)**2) IF (BLPK(L) .LT. 0.) GO TO 39 TEMP = CB(L)**2 - CBI(L)**2 TEMP1 = 2. * CB(L) * CBI(L) DENOM = TEMP**2 + TEMP1**2 TEMP = TEMP / DENOM TEMP1 = -TEMP1 / DENOM GAMMA(L) = 0.5 * (C(L) - (T * TEMP - TI * TEMP1)) / Thick(L) GAMMAI(L) = 0.5 * (CI(L) - (T * TEMP1 + TI * TEMP)) / Thick(L) IF (L .EQ. LX) GO TO 27 GO TO 33 C **SPECIAL CASE, GRADIENT REAL NUMBER 39 IF (CI(L) .EQ. 0.) GO TO 42 TEMP = CB(L)**2 TEM = TEMP**2 TEMP1 = CI(L) COEF1 = CI(L) COEF2 = 2. * TEMP * CI(L) + TI COEF3 = 2. * T * CB(L) COEF4 = TEM * CI(L) - TEMP * TI OLDFN = 1.D20 DO 41 N = 1,10 FN= (((COEF1 * TEMP1) + COEF2) * TEMP1 + COEF3) * TEMP1 + COEF4 FP = ((4. * COEF1 * TEMP1) + 2. * COEF2) * TEMP1 + COEF3 TEMP1 = TEMP1 - FN / FP IF (FN .GE. OLDFN) GO TO 43 OLDFN = FN 41 CONTINUE 43 CBI(L) = TEMP1 GAMMA(L) = .5 *(.5 * (CI(L) *(TEMP - CBI(L)**2) - TI) / * (CB(L) * TEMP1) + C(L)) / Thick(L) GO TO 28 42 GAMMA(L) = C(L) - T / (CB(L)**2 * Thick(L) * 2.) GO TO 33 C **SOUND SPEED AND GRADIENT GIVEN 38 IF (L .EQ. LX) GO TO 33 T = C(L) * (C(L)**2 - 3. * CI(L)**2) TI = CI(L) * (3. * C(L)**2 - CI(L)**2) IF (BLPK(L) .EQ. 0.) GO TO 29 IF (BLPK(L) .LT. 0.) GO TO 28 TEMP = (BLPK(L) / 54575.05416)**2 / Thick(L) * 0.5 T = T * TEMP TI = TI * TEMP TEMP = .5 * C(L) / Thick(L) - GAMMA(L) + T T = -(TI - SQRT( TI * TI + T * TEMP)) / T CB(L) = 54575.05416 * T / BLPK(L) / (1. + T * T) CBI(L) = CB(L) / T GO TO 37 C **SPECIAL CASE, GRADIENT REAL NUMBER 28 TEMP = C(L) - 2. * GAMMA(L) * Thick(L) TEM = TEMP**2 + CI(L)**2 XRE = (T * TEMP + TI * CI(L)) / TEM XIM = -(T * CI(L) - TI * TEMP) / TEM TEM = XRE**2 + XIM**2 CB(L) = SQRT((XRE + SQRT(TEM)) * .5) CBI(L) = .5 * XIM / CB(L) GAMMAI(L) = 0. GO TO 33 29 TEMP = C(L) - 2. * GAMMA(L) * Thick(L) CB(L) = SQRT (T / TEMP) GAMMAI(L) = .5 * (CI(L) - TI / CB(L)**2) / Thick(L) GO TO 33 27 LX = LX - 1 33 CONTINUE C COMPUTE USEFUL QUANTITIES WRITE (JOUT, 58) 58 FORMAT (7X,6H RE G ,8X,6H IM G ,9X,5H L/KM,8X,6H RE C ,8X, * 6H IM C ,5X,12H RE C BOTTOM,4X,12H IM C BOTTOM,10X,9H GRADIENT ) OMEGA = 6.283185307D0 * FREQ DO 40 L = 1, LX TEMP = C(L)**2 + CI(L)**2 CAY(L) = OMEGA * C(L) / TEMP CAYI(L) = -OMEGA * CI(L) / TEMP CAY SQ(L) = CAY(L)**2 - CAYI(L)**2 CAY SQI(L) = 2.D0 * CAY(L) * CAYI(L) TEMDR = -2. * (GAMMA(L) * CAY SQ(L) - GAMMAI(L) * CAY SQI(L)) TEMDI = -2. * (GAMMA(L) * CAY SQI(L) + GAMMAI(L) * CAY SQ(L)) G CU(L) = (TEMDR * C(L) + TEMDI * CI(L)) / TEMP G CUI(L) = (TEMDI * C(L) - TEMDR * CI(L)) / TEMP TEM1 = SQRT( GAMMA(L)**2 + GAMMAI(L)**2) * 2.*OMEGA**2 TEM1 = -TEM1**(1.D0 /3.) TEM1I = ATAN (ABS(GAMMAI(L) /GAMMA(L))) CRTG = TEM1 * COS(TEM1I /3) CRTGI = TEM1 * SIN(TEM1I /3) IF (GAMMA(L) .LT. 0.) CRTG = -CRTG IF (GAMMAI(L).LT. 0.) CRTGI = -CRTGI G(L) = (C(L) * CRTG + CI(L) * CRTGI) / TEMP GI(L) = (C(L) * CRTGI - CI(L) * CRTG) / TEMP CON(L) = G(L) * C(L) - GI(L) * CI(L) CON(L) = OMEGA**2 / CON(L)**2 DPK(L) = -8686.D0 * CAYI(L) WRITE (JOUT, 305) G(L), GI(L), DPK(L), C(L), CI(L), CB(L), : CBI(L), GAMMA(L), GAMMAI(L) G SQI(L) = 2. * G(L) * GI(L) 40 G SQ(L) = G(L)**2 - GI(L)**2 C FIND MODES NXTRA=0 IJ FLAG=0 CALL MODES (NXTRA, IJFLAG) write (*,*) 'Mx:', MX, ' Nxtra:', NXTRA, ' Ijflag:', IJFLAG WRITE (CAPTION(16:20), 105) MX C READ SOURCE AND RECEIVERS DEPTHS READ (JIN,*) SOURCE IF (SOURCE .LT. 0.) THEN WRITE (JOUT,*) ' Source depth is negative. Run halted.' STOP END IF WRITE (CAPTION(21:25), 105) NINT (SOURCE) READ (JIN,'(3f10.4)') HPMIN, HPMAX, DELHP c Data files for the 1979 version contain a blank line to tell c the program that there are no further receiver depths to read: read (jin,*) dummy ! C JX is the number of receiver depths if (delhp .eq. 0.) then jx = 1 else JX = (HPMAX - HPMIN) /DELHP + 1 end if C HPMIN = recvrs(1); ... C Hpmax = recvrs(JX) IF (JX .GE. JDN) JX = JDN - 1 DO 310 J = 1, JX 310 RECVRS(J) = HPMIN + (J - 1) * DELHP WRITE (CAPTION(26:30), 105) NINT(RECVRS(1)) WRITE (CAPTION(31:35), 105) NINT(RECVRS(JX)) WRITE (CAPTION(36:40), 105) NINT(DELHP) WRITE (JOUT,*) WRITE (JOUT,135) CAPTION WRITE (JOUT, 303) 303 FORMAT (/21H SOURCE AND RECEIVERS ) WRITE (JOUT, 21) SOURCE, (RECVRS(J), J = 1, JX) 21 FORMAT ('~', 8F10.1) C COMPUTE DEPTH FUNCTIONS DO 500 M = 1, MX CALL EIGENFUNCTION (M, SOURCE, FS, FIS) DO 390 J = 1, JX CALL EIGENFUNCTION (M, RECVRS(J), F(J), FI(J)) 390 CONTINUE TEMDEN = D(M)*D(M) + DI(M)*DI(M) TEMRE = FS*D(M) + FIS*DI(M) FD = TEMRE/TEMDEN FDI = (D(M) * FIS - DI(M) * FS) / TEMDEN DO 380 J = 1, JX SIZENORM = SQRT((F(J)**2 + FI(J)**2) /SQRT(TEMDEN)) IF (SIZENORM .NE. 0.) THEN ARGNORM = ATAN2(FI(J), F(J)) - ATAN2(DI(M), D(M)) /2 RNORM (J,M) = SIZENORM * COS (ARGNORM) ANORM (J,M) = SIZENORM * SIN (ARGNORM) ELSE RNORM (J,M) = 0. ANORM (J,M) = 0. END IF 380 CONTINUE DO 400 J = 1, JX FF = FD * F (J) - FDI * FI(J) FFI = FD * FI(J) + FDI * F(J) WUU(J,M) = CMPLX(FF, FFI) 400 CONTINUE 500 CONTINUE IF (K1 .eq. 1) THEN WRITE (JOUT, 705) (NINT(RECVRS(J)), J = 1, JX) 705 FORMAT (' Products of magnitudes of source and receiver depth', : ' functions'/ ' Mode', i8, 9i9) DO 700 M = 1, MX 700 WRITE (JOUT, 715) M, (ABS(WUU(J,M)), J = 1, JX) 715 FORMAT ('~', I3, 1X, 1P 10E9.1) END IF WRITE (JOUT, 925) MX 925 FORMAT (' No. of MODES IN SUM:', i8/) C CALCULATE ATTENUATION AND READ IN RANGES READ (JIN,'(3f10.4)') RMIN, RMAX, DELR write (*,*) nint(rmin), nint(rmax), nint(delr) IF (RMIN .LE. 0.) THEN WRITE (JOUT,*) ' Initial Range <= 0. Run halted.' STOP END IF call losses(rmin, rmax, delr, jx, atten) stop end