! FILE FOR3D.F VERSION 1.1 - 09/04/92 ! ****************************************************************** ! PROGRAM FOR3D PROGRAM main ! ****************************************************************** ! * THIS IS AN IMPROVED VERSION OF THE RESEARCH CODE FOR3D. * ! ****************************************************************** ! ****************************************************************** ! * MATHEMATICAL MODEL BY: * ! * * ! * DR. D. LEE , CODE 3122 * ! * ENVIRONMENTAL AND TACTICAL SUPPORT SYSTEMS DEPARTMENT * ! * NAVAL UNDERWATER SYSTEMS CENTER * ! * NEW LONDON, CT 06320, U.S.A. * ! * * ! * PROF. M. H. SCHULTZ * ! * DEPARTMENT OF COMPUTER SCIENCE * ! * YALE UNIVERSITY * ! * NEW HAVEN, CT 06520, U.S.A. * ! * * ! * PROF. YOUCEF SAAD * ! * DEPARTMENT OF COMPUTER SCIENCE * ! * UNIVERSITY OF MINNESOTA * ! * MINNEAPOLIS, MN 55455-0159, U.S.A. * ! * * ! * PROF. W. L. SIEGMANN * ! * DEPARTMENT OF MATHEMATICAL SCIENCES * ! * RENSSELAER POLYTECHNIC INSTITUTE * ! * TROY, NY 12180-3590, U.S.A. * ! ****************************************************************** ! * ORIGINAL VERSION OF MATHEMATICAL MODEL REPORTED IN: * ! * AN EFFICIENT METHOD FOR SOLVING THE THREE-DIMENSIONAL WIDE * ! * ANGLE WAVE EQUATION, RESEARCH REPORT YALEU/DCS/RR-643, * ! * YALE UNIVERSITY, NEW HAVEN, CT., OCT 1986. * ! ****************************************************************** ! ****************************************************************** ! * COMPUTER MODEL BY: * ! * G. BOTSEAS , CODE 3122 * ! * ENVIRONMENTAL AND TACTICAL SUPPORT SYSTEMS DEPARTMENT * ! * NAVAL UNDERWATER SYSTEMS CENTER * ! * NEW LONDON, CONNECTICUT 06320, U.S.A. * ! ****************************************************************** ! * ORIGINAL VERSION OF COMPUTER MODEL REPORTED IN: * ! * FOR3D: A COMPUTER MODEL FOR SOLVING THE LSS THREE-DIMENSIONAL * ! * WIDE ANGLE WAVE EQUATION, NUSC TR 7943, 14AUG87, * ! * G. BOTSEAS, D. LEE, AND D. KING * ! ****************************************************************** ! ****************************************************************** ! * PROJECT SPONSORED BY: * ! * DNL - PROGRAM MANAGER, DR. G.W. MORTON * ! * NUSC - IN-HOUSE LABORATORY INDEPENDENT RESEARCH * ! * PROGRAM MANAGERS, DRS. W.A. VONWINKLE AND K. LIMA * ! * ONR - PROGRAM MANAGERS, DRS. R.L. LAU, R. OBROCHTA, M. ORR, * ! * AND R.BAER * ! * NORDA - PROGRAM MANAGER, R.W. MCGIRR * ! ****************************************************************** ! ****************************************************************** ! * CRAY VERSION - FORTRAN * ! * CRAY UNICOS JOB SCRIPT IS IN FILE FOR3D.JOB. * ! ****************************************************************** ! * MAIN PROGRAM IS FOR3D. * ! ****************************************************************** ! * COMMON BLOCK IS IN FILE FOR3D.CMN * ! * DIMENSIONAL PARAMETERS ARE IN FILE FOR3D.PAR * !ljh * Now commons.mod and params.mod ! ****************************************************************** ! * ALPHABETICAL LIST OF SUBROUTINES FOLLOWS * ! ****************************************************************** ! * AMIFD3 - COMPUTES DIAGONALS FOR MATRIX A. CALLED BY FOR3D. * ! * BCON3D - SUPPLIES BOTTOM CONDITION. CALLED BY RHS. * ! - IF ITYPEB = 0, BOTX AND BOTY ARE SET TO 0. * ! - IF ITYPEB = 1, UBCON3D IS CALLED. * ! - IF ITYPEB = 3, AN ABSORBING LAYER IS INTRODUCED. * ! * BMIFD3 - COMPUTES DIAGONALS FOR MATRIX B. CALLED BY FOR3D. * ! * INDX3D - COMPUTES INDEX OF REFRACTION ARRAYS XN1 AND XN2. * ! - CALLED BY FOR3D. * ! * PORT2D - GENERATES 2D SOLUTION AT PORT SIDEWALL IF ITYPPW=2 * ! * PORT3D - SUPPLIES BOUNDARY CONDITIONS AT PORT SIDEWALL. * ! - IF ITYPPW = 0, THE FIELD ALONG THE WALL IS SET TO 0.* ! - IF ITYPPW = 1, UPORT3D IS CALLED. * ! - IF ITYPPW = 2, PORT3D PREPARES FOR A 2D SOLUTION * ! ALONG THE PORT WALL. * ! - CALLED BY RHS IF NDIM = 3. * ! * PRINTP - PRINTS SELECTED INPUT PARAMETERS. CALLED BY FOR3D. * ! * RHS - COMPUTES RIGHT HAND SIDE OF SYSTEM OF EQUATIONS. * ! - CALLED BY FOR3D. * ! * SCON3D - SUPPLIES SURFACE CONDITION. CALLED BY RHS. * ! - IF ITYPES = 0, SURX AND SURY ARE SET TO 0. * ! - IF ITYPEB = 1, USCON3D IS CALLED. * ! * SFLD3D - GENERATES GAUSSIAN STARTING FIELD IF ISF = 0. * ! - CALLS USFLD3D IF ISF = 1. * ! - GENERATES GREEN'S STARTING FIELD IF ISF = 2. * ! - CALLED BY FOR3D. * ! * STBD2D - GENERATES 2D SOLUTION AT STBD SIDEWALL IF ITYPSW=2 * ! - AND NDIM = 3. CALLED BY TWOSTEP. * ! * STBD3D - SUPPLIES BOUNDARY CONDITIONS AT STARBOARD SIDEWALL. * ! - IF ITYPSW = 0, THE FIELD ALONG THE WALL IS SET TO 0.* ! - IF ITYPSW = 1, USTBD3D IS CALLED. * ! - IF ITYPSW = 2, STBD3D PREPARES FOR A 2D SOLUTION. * ! ALONG THE STARBOARD WALL. * ! - CALLED BY RHS IF NDIM = 3. * ! * SVP3D - READS SOUND SPEED PROFILES IN INPUT RUNSTREAM. * ! - CALLED BY FOR3D IF KSVP = 0. * ! * TRID3D - SOLVES A TRIDIAGONAL SYSTEM OF EQUATIONS. * ! - CALLED BY PORT2D, STBD2D, AND TWOSTEP. * ! * TWOSTEP - MARCHES THE ACOUSTIC FIELD FORWARD. * ! - FIRST STEP SOLVES IN RANGE. * ! - SECOND STEP SOLVES IN AZIMUTH. * ! - CALLED BY FOR3D. * ! - IF NDIM = 3 AND ITYPPW = 2, TWOSTEP CALLS PORT2D. * ! - IF NDIM = 3 AND ITYPSW = 2, TWOSTEP CALLS STBD2D. * ! - IF NDIM = 1 OR 2, TRID3D IS CALLED ONCE. * ! - IF NDIM = 3, TRID3D IS CALLED TWICE. * ! ****************************************************************** ! * USER SUBROUTINES * ! ****************************************************************** ! * UBCON3D - USER SUPPLIED BOTTOM CONDITIONS. CALLED BY BCON3D. * ! * UBOTTOM - USER SUPPLIED BOTTOM PROFILE INPUT ROUTINE. * ! - CALLED BY FOR3D IF KBOT NE 0. * ! * UEXACT - EXACT SOLUTION. USER WRITTEN. FOR TEST PURPOSES. * ! - CALLED BY FOR3D. * ! * UPORT3D - USER SUPPLIED PORT WALL CONDITIONS. CALLED BY PORT3D* ! * USCON3D - USER SUPPLIED SURFACE CONDITIONS. CALLED BY SCON3D. * ! * USFLD3D - USER SUPPLIED STARTING FIELD. * ! - CALLED BY SFLD3D IF ISF = 1. * ! * USTBD3D - USER SUPPLIED STBD WALL CONDITIONS. CALLED BY STBD3D* ! * USVP3D - USER SOUND SPEED PROFILE INPUT ROUTINE. * ! - CALLED BY FOR3D IF KSVP NOT 0. * ! ****************************************************************** ! * ALPHABETICAL LIST OF SOME PROGRAM VARIABLES FOLLOWS: * ! ****************************************************************** ! * AL - ARRAY - LOWER DIAGONAL OF MATRIX A - COMPLEX ! * ALFA - SET TO 1/12 IF DOUGLAS METHOD REQUESTED ! * ALPHA - VOLUME ATTENUATION - DB/METER ! * AM - ARRAY - MAIN DIAGONAL OF MATRIX A - COMPLEX ! * ANG - RELATIVE BEARING - DEGREES ! * AU - ARRAY - UPPER DIAGONAL OF MATRIX A - COMPLEX ! * ATT - ATTENUATION COEFFICIENT FOR ARTIFICIAL ABSORBING LAYER ! * BETA - ARRAY - ATTENUATION IN LAYERS - DB/WAVELENGTH ! * BL - ARRAY - LOWER DIAGONAL OF MATRIX B - COMPLEX ! * BM - ARRAY - MAIN DIAGONAL OF MATRIX B - COMPLEX ! * BOTX - ARRAY - COMPLEX PRESSURE AT BOTTOM AT ADVANCED RANGE RA+DR ! * BOTY - ARRAY - COMPLEX PRESSURE AT BOTTOM AT PRESENT RANGE RA ! * BTA - ARRAY - PARTIAL SOLUTION OF SYSTEM OF EQUATIONS ! * BU - ARRAY - UPPER DIAGONAL OF MATRIX B - COMPLEX ! * C0 - REFERENCE SPEED OF SOUND - METERS/SEC ! * CSVP - ARRAY - SOUND VELOCITY - METERS/SEC ! * D - ARRAY - RIGHT HAND SIDE OF SYSTEM OF EQUATIONS - COMPLEX ! * DEG - CONVERSION FACTOR - DEGREES/RADIAN ! * DOUGRA - 0 - DON'T USE DOUGLAS METHOD ! NOT 0 - RANGE AT WHICH TO BEGIN DOUGLAS METHOD - METERS ! * DR - RANGE STEP - METERS ! * DTH - WIDTH OF SECTOR - DEGREES ! * DZ - DEPTH INCREMENT OF SOLUTION - METERS ! * DZZ - DEPTH INCREMENT FOR ADJUSTING LAYER DEPTHS IN SLOPING BOTTOM ! * FLDW - WIDTH OF FIELD - DEGREES ! * FRQ - FREQUENCY - HZ ! * HNK - HANKEL FUNCTION H0(1) ! * HNKL - EXTERNAL FUNCTION - COMPUTES HANKEL FUNCTION H0(1) ! * IBOT - BOTTOM DEPTH PRINT FLAG ! 0 - DO NOT PRINT BOTTOM DEPTHS ! 1 - PRINT BOTTOM DEPTHS ! * IHNK - HANKEL FUNCTION FLAG ! 0 - HANKEL FUNCTION NOT USED. 10*LOG(R) ADDED TO SOLUTION. ! 1 - STARTING FIELD DIVIDED BY HANKEL FUNCTION. ! SOLUTION MULTIPLIED BY HANKEL FUNCTION BEFORE ! COMPUTING PROPAGATION LOSS. ! * ILYR - INDEX FOR ARRAYS BETA, ZLYR, AND RHO ! * IPZ - EVERY IPZ'TH VALUE IN DEPTH IS PRINTED ! * ISF - STARTING FIELD FLAG ! 0 - PROGRAM GENERATES GAUSSIAN STARTING FIELD ! AT RANGE = 0.0. SEE SUBROUTINE SFLD3D. ! 1 - USER SUPPLIES STARTING FIELD. SEE SUBROUTINE USFLD3D. ! 2 - GREENS WIDE ANGLE STARTER ! 3 - AVAILABLE ! * ISFLD - STARTING FIELD PRINT FLAG ! 0 - DO NOT PRINT STARTING FIELD ! 1 - PRINT STARTING FIELD ! * ISVP - SVP PRINT FLAG ! 0 - DO NOT PRINT SOUND VELOCITY PROFILE ! 1 - PRINT SOUND VELOCITY PROFILE ! * ITEMP - TEMPORARY VARIABLE ! * ITRK - ARRAY - INDICES FOR ARRAY TRACK ! * ITYPEB - TYPE OF BOTTOM ! 0 - BCON3D SETS BOTY AND BOTX = 0.0 ! 1 - USER SUPPLIES BOTTOM CONDITION. SEE SUBROUTINE UBCON3D. ! 2 - SPARE. ! 3 - ABSORBING LAYER INTRODUCED - BOTTOM OF LAYER IS FLAT ! 4 - SPARE. ! * ITYPES - TYPE OF SURFACE ! 0 - PRESSURE RELEASE. SCON3D SETS SURY AND SURX = 0.0 ! 1 - USER SUPPLIES SURFACE CONDITION. SEE SUBROUTINE USCON3D. ! 2 - SPARE ! * ITYPPW - TYPE OF PORT SIDEWALL BOUNDARY CONDITION. ! 0 - FIELD ALONG PORT SIDEWALL IS SET TO 0.0. ! 1 - USER SUPPLIED. SEE SUBROUTINE UPORT3D. ! 2 - MODEL GENERATES 2D SOLUTION IF NDIM = 3. ! * ITYPSW - TYPE OF STBD SIDEWALL BOUNDARY CONDITION. ! 0 - FIELD ALONG STBD SIDEWALL IS SET TO 0.0. ! 1 - USER SUPPLIED. SEE SUBROUTINE USTBD3D. ! 2 - MODEL GENERATES 2D SOLUTION IF NDIM = 3. ! * IWZ - INDEX INCREMENT OF RECEIVER SOLUTIONS TO BE WRITTEN ON DISK. ! * IXSVP - ARRAY OF POINTERS WHICH POINT TO ENTRIES IN CSVP AND ZSVP. ! - IXSVP(1,L) POINTS TO BOTTOM DEPTH AND SPEED IN LAYER 1 ! ON SECTOR BOUNDARY L. ! - IXSVP(2,L) POINTS TO BOTTOM DEPTH AND SPEED IN LAYER 2 ! ON SECTOR BOUNDARY L. ! ETC. ! * KSVP - SVP PROFILE FLAG. ! 0 - PROFILE IS IN RUNSTREAM - SEE SUBROUTINE SVP3D. ! NOT 0 - USER ROUTINE GENERATES PROFILES - SEE SUBROUTINE ! USVP3D. ! - KSVP MAY BE USED AS AN INDEX IN A GOTO STATEMENT IN USVP3D. ! * KBOT - BOTTOM PROFILE FLAG. ! 0 - PROFILE IS IN RUNSTREAM . ! NOT 0 - USER ROUTINE GENERATES PROFILES - SEE SUBROUTINE ! UBOTTOM. ! * MLYR - TEMPORARY - NUMBER OF LAYERS INVOLVED IN SPECIFIC CALCULATION. ! * MM - INDEX - MM+1 POINTS TO FIRST VALUE OF ARTIFICIAL ABSORBING ! LAYER IN ARRAY U. ! * MXLYR - PARAMETER - MAXIMUM NUMBER OF LAYERS. ! - MAX DIMENSION OF ARRAYS BETA,RHO,ZLYR,IXSVP. ! * MXN - PARAMETER - MAXIMUM DIMENSION OF C, X, Y, AND U ARRAYS. ! * MXSVP - PARAMETER - MAXIMUM DIMENSION OF ARRAYS CSVP AND ZSVP. ! * MXTRK - PARAMETER - MAXIMUM DIMENSION OF ARRAY TRACK. ! * N - NUMBER OF EQUI-SPACED POINTS IN U. ! - INCLUDES BOTTOM POINT - DOES NOT INCLUDE SURFACE POINT. ! * N1 - DIAGONAL ELEMENTS N1 THRU N WILL BE COMPUTED - NOT USED. ! * NA - NUMBER OF POINTS IN ABSORBING LAYER. ! * NDIV - IF DOUGLAS METHOD IS REQUESTED, DIVIDE N BY NDIV. ! * NIU - PARAMETER - INPUT UNIT NUMBER. ! * NLYR - NUMBER OF LAYERS. ! * NOLD - NUMBER OF RECEIVER DEPTHS ON ENTRY TO ROUTINE CRNK. ! * NOU - PARAMETER - OUTPUT UNIT NUMBER. ! * NPU - PARAMETER - PRINTER UNIT NUMBER. ! * NSEC - NUMBER OF SECTORS. ! * NSOL - NUMBER OF SOLUTIONS (NSEC+1). ! * NSVP - NUMBER OF POINTS IN CSVP AND ZSVP. ! * NTOT - TOTAL POINTS IN SOLUTION (N*NSOL). ! * NWSVP - NUMBER OF POINTS IN LAYER 1 SVP. ! * NZ - NUMBER OF SOLUTION DEPTHS TO BE WRITTEN ON DISK. ! * OLDR - RANGE INCREMENT AT START OF PROBLEM. ! * PDR - RANGE INCREMENT AT WHICH SOLUTION IS PRINTED - METERS. ! * PDTH - AZIMUTHAL INCREMENT AT WHICH SOLUTION IS PRINTED - DEGREES. ! * PDZ - DEPTH INCREMENT AT WHICH SOLUTION IS PRINTED - METERS. ! * PHI - WIDTH OF SECTOR IN RADIANS. ! * PI - 3.141592654 ! * PL - PROPAGATION LOSS - DB. ! * PORTX - ARRAY - COMPLEX FIELD AT PORT SIDEWALL - ADVANCED RANGE. ! * PORTY - ARRAY - COMPLEX FIELD AT PORT SIDEWALL - PRESENT RANGE. ! * R0 - INITIAL RANGE - METERS. ! * R1 - RANGE AT WHICH BOTTOM DEPTH IS AVAILABLE - METERS. ! * R2 - NEXT RANGE AT WHICH BOTTOM DEPTH IS AVAILABLE - METERS. ! * RA - HORIZONTAL RANGE OF STARTING FIELD FROM SOURCE - METERS. ! - RA IS SET TO 0.0 IF STARTING FIELD IS GAUSSIAN. RA IS ! - INCREMENTED BY DR AS SOLUTION IS MARCHED OUT IN RANGE. ! * RA+DR - RANGE TO WHICH SOLUTION IS TO BE ADVANCED - METERS. ! * RHO - ARRAY - DENSITY IN LAYER . ! * RHOG - ARRAY - DENSITY GRADIENT IN LAYER. ! * RMAX - MAXIMUM RANGE OF SOLUTION - METERS. ! * RSVP - NEXT RANGE AT WHICH SVP IS AVAILABLE - METERS. ! * STBDX - ARRAY - COMPLEX FIELD AT STBD SIDEWALL - ADVANCED RANGE. ! * STBDY - ARRAY - COMPLEX FIELD AT STBD SIDEWALL - PRESENT RANGE. ! * SURX - ARRAY - COMPLEX PRESSURE AT SURFACE AT ADVANCED RANGE RA+DR. ! * SURY - ARRAY - COMPLEX PRESSURE AT SURFACE AT PRESENT RANGE RA. ! * TEMP - TEMPORARY VARIABLE - COMPLEX. ! * THETA - ARRAY - SLOPE OF BOTTOM ALONG SECTOR BOUNDARIES - RADIANS. ! .EQ.0 -- FLAT BOTTOM . ! .GT.0 -- SHALLOW TO DEEP. ! .LT.0 -- DEEP TO SHALLOW. ! * TITLE - 80 CHARACTER TITLE. ! * TM - ARRAY - TIME OF DAY. ! * TRACK - 3 DIM. ARRAY - RANGE AND DEPTH OF WATER - METERS. ! - ASSUMES FLAT BOTTOM AT PRESENT TIME. ! * U - ARRAY - COMPLEX ACOUSTIC PRESSURE FIELD. ! * U1-U12 - USER MAY USE THESE VARIABLES TO INPUT DATA REQUIRED BY USER ! - ROUTINES - REAL - IN COMMON BLOCK. ! * UA-UE - USER MAY USE THESE VARIABLES - COMPLEX - IN COMMON BLOCK ! * V - ARRAY - PARTIAL SOLUTION OF COMPLEX FIELD ! * WDR - RANGE STEP AT WHICH SOLUTION IS WRITTEN ON DISK - METERS ! * WDTH - AZIMUTHAL INCREMENT AT WHICH SOLUTION IS WRITTEN ON DISK - DEG ! * WDZ - DEPTH INCREMENT AT WHICH SOLUTION IS WRITTEN ON DISK - METERS ! - WDZ SHOULD BE SELECTED SO THAT PLOT PROGRAM DOES NOT ! - INTERPOLATE BETWEEN WIDELY SPACED RECEIVERS. ! * WZ1 - FIRST RECEIVER DEPTH AT WHICH SOLUTION IS WRITTEN ON DISK ! * WZ2 - LAST RECEIVER DEPTH AT WHICH SOLUTION IS WRITTEN ON DISK ! - IN OTHER WORDS, WRITE WZ1 TO WZ2 BY WDZ - METERS ! * XK0 - REFERENCE WAVE NUMBER ! * XN2 - ARRAY - INDEX OF REFRACTION SQUARED - COMPLEX ! * XPR - RANGE AT WHICH SOLUTION IS PRINTED - METERS ! * XWR - RANGE AT WHICH SOLUTION IS WRITTEN ON DISK - METERS ! * Z1 - DEPTH OF WATER AT RANGE R1 - METERS ! * Z2 - DEPTH OF WATER AT RANGE R2 - METERS ! * ZA - DEPTH OF FIELD AT RANGE RA - METERS ! - INITIAL DEPTH OF STARTING FIELD AT RANGE RA IS ! - AS FOLLOWS: ! - IF ITYPEB = 0 OR 1, ZA IS MAXIMUM DEPTH OF ! - BOTTOM-MOST SEDIMENT LAYER AT INITIAL RANGE OF ! - STARTING FIELD. IF ITYPEB = 3, ZA IS MAXIMUM ! - DEPTH OF ARTIFICIAL ABSORBING LAYER AT INITIAL ! - RANGE OF STARTING FIELD. PROGRAM INSERTS LAYER. ! - RHO AND BETA ARE OBTAINED FROM LAYER ABOVE. ! - SPEED IS BOTTOM-MOST SPEED FROM LAYER ABOVE. ! - BOTTOM OF ABSORBING LAYER REMAINS FLAT. ! * ZI - DEPTH OF RECEIVER 'I' - METERS. ! * ZLYR - ARRAY - DEPTH TO BOTTOM OF LAYER - METERS. ! * ZS - SOURCE DEPTH - METERS. ! * ZSVP - ARRAY - DEPTH OF SOUND VELOCITY - METERS. ! ****************************************************************** ! ****************************************************************** ! * INPUT UNIT NUMBER = NIU ! * INPUT FILE NAME = for3d.in - SEE CRAY UNICOS JOB SCRIPT ! * CONTENTS: INPUT RUNSTREAM IN FREE FORMAT ! ! LINE INPUT DATA ! ! 1 TITLE ! ! 2 NDIM ! ! 3 FRQ,ZS,C0,ISF,RA,ZA,N,IHNK,ITYPES,ITYPEB,ITYPPW,ITYPSW,FLDW,NSEC ! ! 4 RMAX,DR,WDR,WZ1,WZ2,WDZ,WDTH,PDR,PDZ,PDTH,ISFLD,ISVP,IBOT ! ! 5 DOUGRA,NDIV ! ! 6 U1,U2,U3,U4,U5,U6,U7,U8,U9,U10,U11,U12 ! ! 7 KBOT ! ! 8 R1,Z1 ** IF KBOT = 0, BOTTOM PROFILE IS IN RUNSTREAM. ! 9 R2,Z2 * RANGE, WATER DEPTH (METERS). ! . ** -1,-1 MARKS END OF THIS PROFILE. ! . *** ENTER PROFILE FOR EACH SECTOR BOUNDARY. ! . ** PORT BOUNDARY FIRST. ! N . * IF KBOT NE 0, UBOTTOM IS CALLED. OMIT ! N+1 -1,-1 ** LINES 8 THRU N+1. ! ! ************************************************ ! N+2 RSVP * ! N+3 KSVP * ! ************************************** * ! N+4 NLYRS(L) * * ! **************************** * ** ! N+5 ZLYR(I,L),RHO(I,L),RHOG(I,L), * ** *** NOTE ! BETA(I,L),BETAG(I,L) ** *** NOTE *** 1 ! N+6 ZSVP(1,L),CSVP(1,L) *** NOTE *** 2 ** ! N+7 ZSVP(2,L),CSVP(2,L) ** 3 ** * ! . . . * * * ! N+M ZSVP(J,L),CSVP(J,L) ************ * * ! ********************** * ! ******************************** ! ! ****************************************************************** ! *** QUICK REFERENCE AND NOTES FOR INPUT RUNSTREAM ! *** UNITS: METERS AND METERS/SEC EXCEPT AS NOTED ! ****************************************************************** ! LINE INPUT DATA ! ! 1 TITLE = USER NOTE. 80 CHARACTERS MAXIMUM. ! ! 2 NDIM = NUMBER OF DIMENSIONS. ! 1 - 2D SOLUTION. ! 2 - N x 2D SOLUTIONS. ! 3 - 3D SOLUTIONS. ! ! 3 FRQ = FREQUENCY (HZ) ! ! ZS = SOURCE DEPTH ! ! C0 = REFERENCE SOUND SPEED. IF C0 = 0.0, C0 IS SET TO AVERAGE ! SPEED IN FIRST LAYER. ! ! ISF = STARTING FIELD FLAG. ! 0 = GAUSSIAN. ! 1 = USER FIELD. ! 2 = GREENS STARTER. ! ! RA = HORIZONTAL RANGE FROM SOURCE TO STARTING FIELD. ! RA IS SET TO 0.0 IF ISF = 0. ! ! ZA = DEPTH OF STARTING FIELD AT RANGE RA. IF ZA = 0.0, ZA IS SET ! TO MAX DEPTH OF BOTTOM LAYER IN FIRST PROFILE. IF ITYPEB = ! 2 OR 3 AND ZA = 0.0, ZA IS SET TO (4/3)*MAX DEPTH OF BOTTOM ! LAYER. IF ITYPEB = 3 AND ZA NOT ZERO, THE ARTIFICIAL ! BOTTOM LAYER IS EXTENDED TO ZA METERS PROVIDED THAT ZA ! IS GREATER THAN OR EQUAL TO MAX DEPTH OF BOTTOM LAYER ! IN FIRST PROFILE. ! ! N = NUMBER OF EQUISPACED RECEIVERS IN STARTING FIELD. IF N = 0, ! N IS SET SO THAT THE RECEIVER DEPTH INCREMENT IS ! EQUAL TO 1 METER. IF N IS GREATER THAN MXN ! N IS SET TO MXN. ! ! IHNK = HANKEL FUNCTION FLAG. IHNK = 0, DON'T USE HANKEL FUNCTION. ! IHNK = 1, DIVIDE STARTING FIELD BY HANKEL FUNCTION, THEN ! MULTIPLY THE SOLUTION FIELD BY HANKEL FUNCTION BEFORE ! COMPUTING PROPAGATION LOSS. IF STARTING FIELD IS GAUSSIAN, ! IHNK SHOULD BE SET TO 0. IF STARTING FIELD IS ELLIPTIC, ! IHNK SHOULD BE SET TO 1. ! ! ITYPES = TYPE OF SURFACE ! 0 - PRESSURE RELEASE. SCON3D SETS SURY AND SURX = 0.0 ! 1 - USER SUPPLIES SURFACE CONDITION. SEE SUBROUTINE USCON3D. ! 2 - SPARE. ! ! ITYPEB = TYPE OF BOTTOM ! 0 - PRESSURE RELEASE. BCON3D SETS BOTY AND BOTX = 0.0 ! 1 - USER SUPPLIES BOTTOM CONDITION. SEE SUBROUTINE UBCON3D. ! 2 - SPARE. ! 3 - ABSORBING LAYER INTRODUCED - FLAT BOTTOM ! 4 - SPARE ! ! ITYPPW = TYPE OF PORT SIDEWALL BOUNDARY CONDITION ! 0 - FIELD ALONG PORT SIDEWALL IS SET TO 0.0 ! 1 - USER SUPPLIED. SEE SUBROUTINE UPORT3D. ! 2 - MODEL GENERATES 2D SOLUTION IF NDIM = 3. ! ! ITYPSW = TYPE OF STARBOARD SIDEWALL BOUNDARY CONDITION ! 0 - FIELD ALONG STARBOARD SIDEWALL IS SET TO 0.0 ! 1 - USER SUPPLIED. SEE SUBROUTINE USTBD3D. ! 2 - MODEL GENERATES 2D SOLUTION IF NDIM = 3. ! ! FLDW = WIDTH OF FIELD IN DEGREES. IGNORED IF NDIM = 1. ! ! NSEC = NUMBER OF SECTORS IN FIELD. IGNORED IF NDIM = 1. ! NUMBER OF SOLUTIONS = NSOL = NSEC+1. ! ! 4 RMAX = MAXIMUM RANGE OF SOLUTION ! ! DR = RANGE STEP. IF DR = 0, DR IS SET TO 1 METER. ! IF BOTTOM OF PROBLEM IS NOT FLAT, DR IS RECOMPUTED ! SO THAT MAX DEPTH IS EITHER INCREMENTED OR DECREMENTED ! BY DZ. SOLUTION IS COMPUTED EVERY DR METERS. ! ! WDR = RANGE STEP AT WHICH SOLUTION IS WRITTEN ON DISK. ! IF WDR NOT 0, AN OUTPUT DISK FILE IS ASSIGNED. ! WDR IS ROUNDED TO NEAREST DR. ! ! WZ1 = FIRST RECEIVER DEPTH AT WHICH SOLUTION IS WRITTEN ON DISK. ! ! WZ2 = LAST RECEIVER DEPTH AT WHICH SOLUTION IS WRITTEN ON DISK. ! ! WDZ = DEPTH INCREMENT AT WHICH SOLUTION IS WRITTEN ON DISK. ! ROUNDED TO NEAREST DZ. ! ! WDTH = AZIMUTHAL INCREMENT AT WHICH SOLUTION IS WRITTEN ON DISK. ! ROUNDED TO NEAREST DTH. ! ! PDR = RANGE STEP AT WHICH SOLUTION IS PRINTED. ! ROUNDED TO NEAREST DR. ! ! PDZ = DEPTH INCREMENT AT WHICH SOLUTION IS PRINTED. ! ROUNDED TO NEAREST DZ. ! ! PDTH = AZIMUTHAL INCREMENT AT WHICH SOLUTION IS PRINTED. ! ROUNDED TO NEAREST DTH. ! ! ISFLD = 0 - DON'T PRINT STARTING FIELD. ! = 1 - PRINT STARTING FIELD. ! ! ISVP = 0 - DON'T PRINT SOUND VELOCITY PROFILE. ! = 1 - PRINT SOUND VELOCITY PROFILE. ! ! IBOT = 0 - DON'T PRINT BOTTOM DEPTHS. ! = 1 - PRINT BOTTOM DEPTHS. ! ! 5 DOUGRA = RANGE AT WHICH TO SWITCH FROM CRANK-NICOLSON METHOD TO ! DOUGLAS METHOD. IF DOUGRA = 0, USE CRANK-NICOLSON METHOD ! ONLY. SUGGESTED VALUE IS 5000 METERS. ! ! NDIV = IF DOUGRA NE 0, N IS DIVIDED BY NDIV RESULTING IN ONLY ! N/NDIV SOLUTIONS IN DEPTH. SPEED UP IS NDIV TIMES. ! SUGGESTED VALUE IS 5. ! ! 6 U1-U12 = USER VARIABLES - REAL, SINGLE PRECISION. ! SEE HARVARD SUBROUTINE HARVARD.F FOR EXAMPLE. ! ! 7 KBOT = BOTTOM PROFILE FLAG. ! = 0 - BOTTOM PROFILE IN INPUT RUNSTREAM. ! = NOT ZERO. PROFILE IS SUPPLIED BY SUBROUTINE UBOTTOM. ! USER WRITES UBOTTOM. IF NSOL IS LARGE, UBOTTOM IS ! PREFERRED METHOD OF INPUT. ! ! BOTTOM PROFILE AT LEFTMOST SECTOR BOUNDARY. ! 8 R1,Z1 = RANGE AND DEPTH OF WATER. ! 9 R2,Z2 = ETC. ! . ! N+1 -1,-1 = MARKS THE END OF THIS BOTTOM PROFILE. ! ! NEXT PROFILE. ! 8 R1,Z1 = RANGE AND DEPTH OF WATER. ! 9 R2,Z2 = ETC. ! . ! N+1 -1,-1 = MARKS THE END OF THIS BOTTOM PROFILE. ! ! ETC. REPEAT UNTIL ALL PROFILES ENTERED. ! ! N+2 RSVP = RANGE OF SOUND SPEED PROFILES. ! SEE NOTE 1. ! ! N+3 KSVP = SVP FLAG. ! = 0 - SOUND SPEED PROFILES IN INPUT RUNSTREAM. ! = NOT ZERO. PROFILE (LINES N+4 THRU N+M) IS SUPPLIED BY ! USER. USER WRITES SUBROUTINE USVP3D. KSVP MAY BE USED IN ! COMPUTED GOTO STATEMENT TO TRANSFER CONTROL IN USVP3D. ! IF NSOL IS LARGE, USVP3D IS PREFERRED METHOD OF INPUT. ! ! N+4 NLYRS(L) = NUMBER OF LAYERS. IF ITYPEB = 3, PROGRAM INSERTS ! AN ARTIFICIAL LAYER AND INCREMENTS NLYRS(L) BY 1. ! SEE NOTE 2. ! ! N+5 ZLYR(I,L) = MAX DEPTH OF LAYER I IN PROFILE. ! ! RHO(I,L) = DENSITY IN LAYER I (G/CM**3). ! ! RHOG(I,L) = DENSITY GRADIENT IN LAYER I (G/CM**3/M). ! ! BETA(I,L) = ATTENUATION IN LAYER I (DB/WAVELENGTH). ! IF BETA(I,L) IS NEGATIVE, ATTENUATION IS COMPUTED. ! SEE NOTE 3. ! ! BETAG(I,L)= ATTENUATION GRADIENT IN LAYER I (DB/WAVELENGTH/M). ! ! N+6 ZSVP(1,L) = DEPTH TO TOP OF LAYER I ! ! CSVP(1,L) = SPEED OF SOUND AT TOP OF LAYER I ! . ! . ! N+M ZSVP(J,L) = DEPTH TO BOTTOM OF LAYER I ! ! CSVP(J,L) = SPEED OF SOUND AT BOTTOM OF LAYER I ! IF ONLY ONE SVP INPUTTED, IT IS USED THRU ENTIRE PROBLEM. ! IF MORE THAN ONE SVP INPUTTED, LAST SVP IS USED THRU ! REMAINDER OF PROBLEM. ! ! NOTE 1 REPEAT LINES N+2 THRU N+M FOR EACH SET OF PROFILES TO BE ENTERED. ! A SET OF PROFILES CONSISTS OF NSOL (NSEC+1) PROFILES LOCATED ON ! SECTOR BOUNDARIES AND EQUISPACED FROM PORT TO STARBOARD IN AZIMUTH. ! NSEC IS THE NUMBER OF SECTORS IN THE FIELD OF INTEREST. NSEC ! ADJACENT SECTORS HAVE NSOL BOUNDARIES. ! ! NOTE 2 REPEAT LINES N+4 THRU N+M FOR EACH PROFILE IN THE SET. AT RANGE ZERO ! ALL PROFILES ARE THE SAME. BEYOND RANGE ZERO, L PROFILES ARE ! ENTERED FROM PORT TO STARBOARD AS VIEWED FROM RANGE R0. L IS ! VARIED FROM 1 THRU NSOL. ENTERING -1,-1 AFTER PROFILE L ! WILL CAUSE THAT PROFILE TO BE ENTERED REPEATEDLY UNTIL ! NSOL PROFILES HAVE BEEN ENTERED. ! ! NOTE 3 REPEAT LINES N+5 THRU N+M FOR EACH LAYER I IN THE PROFILE. ! J IS THE NUMBER OF POINTS IN LAYER I. ! ! ****************************************************************** ! *** OUTPUT ! ****************************************************************** ! OUTPUT UNIT NUMBER = NOU ! OUTPUT FILE NAME = for3d.out - SEE CRAY UNICOS JOB SCRIPT ! CONTENTS: AS FOLLOWS - UNFORMATTED ! ! THE FOLLOWING VARIABLES ARE WRITTEN ONCE: ! ! WRITE(NOU)NDIM,FRQ,ZS,C0,ISF,RA,ZA,N,IHNK,ITYPES,ITYPEB,ITYPPW,ITYPSW,FLDW, ! NSEC,NSOL,RMAX,DR,WDR,WZ1,WZ2,WDZ,DZ,DOUGRA,NDIV, ! U1,U2,U3,U4,U5,U6,U7,U8,U9,U10,U11,U12 ! ! THE VALUES OF THE VARIABLES AND ARRAYS LISTED IN THE FOLLOWING ! WRITE STATEMENT ARE WRITTEN AT R0 AND AT EACH WDR (WRITE RANGE INCREMENT) ! THEREAFTER. ! ! WRITE(NOU)ANG,NZ,RA,WZ1,WDZ,(U(M+I),I=IWZ1,IWZ2,IWZ) ! ! where M = (J-1)*N ! J = SOLUTION INDEX. J VARIES FROM 1 TO NSOL. PORT TO STARBOARD. ! N = NUMBER OF RECEIVERS IN VERTICAL COLUMN. ! ANG = RELATIVE BEARING OF SOLUTION IN DEGREES. PORT IS NEGATIVE. ! NZ = NUMBER OF RECEIVERS WRITTEN ON DISK. ! RA = RANGE IN METERS. ! WZ1 = DEPTH OF FIRST RECEIVER. ! WDZ = DEPTH INCREMENT OF RECEIVERS WRITTEN ON DISK IN METERS. ! U = SOLUTION ARRAY. COMPLEX. ! ! PRINTED OUTPUT IS DIRECTED TO PRINTER UNIT NUMBER NPU. ! ****************************************************************** ! ****************************************************************** ! ****************************************************************** ! ! *** INCLUDE COMMON BLOCK !ljh INCLUDE 'FOR3D.CMN' USE commons IMPLICIT REAL (KIND=8) (A-H,O-Z) IMPLICIT INTEGER (KIND=4) (I-N) CHARACTER*80 TITLE ! cljh COMPLEX HNKL ! cljh ! ! *** TIME OF DAY AT START OF RUN - CRAY ROUTINE !CC CALL SECOND(CPT1) ! ! *** READ INPUT PARAMETERS. ! *** FILES ARE LINKED TO UNIT NUMBERS IN CRAY JOB SCRIPT. ! OPEN(NIU) ! CALL ASSIGN(NIU,'FOR3D.IN') !ljh OPEN(UNIT=NIU,FILE='FOR3D.IN') ! READ (NIU, 15 ) TITLE WRITE (NPU, 15 ) TITLE 15 FORMAT(A80) ! ! *** READ SOLUTION DIMENSION READ(NIU,*)NDIM IF(NDIM.LT.1.OR.NDIM.GT.3)THEN WRITE(NPU,*)'INPUT ERROR. CHECK NDIM VALUE.' STOP ELSE ENDIF ! READ(NIU,*)FRQ,ZS,C0,ISF,RA,ZA,N,IHNK,ITYPES,ITYPEB, & ITYPPW,ITYPSW,FLDW,NSEC ! READ(NIU,*) RMAX,DR,WDR,WZ1,WZ2,WDZ,WDTH,PDR,PDZ,PDTH, & ISFLD,ISVP,IBOT ! READ(NIU,*) DOUGRA,NDIV ! READ(NIU,*) U1,U2,U3,U4,U5,U6,U7,U8,U9,U10,U11,U12 ! ! SET DOUGLAS COEFFICIENT TO 0 ALFA=0.0 ! ! *** IF C0 NE 0.0 COMPUTE REFERENCE WAVE NUMBER IF(C0.NE.0.0)XK0=2.0*PI*FRQ/C0 ! ! *** COMPUTE NUMBER OF SOLUTIONS. ONE FOR EACH SECTOR BOUNDARY. ! IF(NDIM.EQ.1)THEN ! *** USER REQUESTED A 2D SOLUTION. ONLY ONE SOLUTION IS GENERATED. FLDW=0.0 NSEC=0 NSOL=1 NSOLW=1 NSOLP=1 DTH=0 WDTH=0.0 PDTH=0.0 ANG=0.0 PHI=0.0 ! ELSE ! *** MORE THAN 1 SOLUTION. NSOL=NSEC+1 ! *** COMPUTE SECTOR WIDTH IN DEGREES. DTH=FLDW/NSEC ! *** COMPUTE SECTOR WIDTH IN RADIANS. PHI=PI*FLDW/(NSEC*180) ! *** COMPUTE NUMBER OF SOLUTIONS IN AZIMUTH TO WRITE TO DISK. ! *** AVOID CRAY ROUND UP PROBLEMS RNSOLW=(WDTH/DTH+.5) NSOLW=INT(RNSOLW) IF(NSOLW.EQ.0)NSOLW=1 ! *** COMPUTE NUMBER OF SOLUTIONS IN AZIMUTH TO PRINT. RNSOLP=(PDTH/DTH+.5) NSOLP=INT(RNSOLP) IF(NSOLP.EQ.0)NSOLP=1 ENDIF ! IF(NSOL.GT.MXSOL) THEN WRITE(NPU,*)'TOO MANY SOLUTIONS REQUESTED.' WRITE(NPU,*)'INCREASE PARAMETER MXSEC AND RECOMPILE.' STOP ELSE ENDIF ! ! *** IF GAUSSIAN STARTING FIELD, RA MUST BE 0. IF(ISF.EQ.0) RA=0.0 R0=RA ! ! *** READ BOTTOM PROFILE FLAG READ(NIU,*) KBOT ! ! *** IS PROFILE IN RUNSTREAM? IF(KBOT.EQ.0) THEN ! *** YES ! *** READ BOTTOM PROFILE - RANGE,DEPTH ZBOT=0.0 DO 24 L=1,NSOL DO 22 I=1,MXTRK+1 READ(NIU,*) TRACK(I,1,L),TRACK(I,2,L) ! *** END OF PROFILE? IF(TRACK(I,1,L).LT.0.0) GO TO 23 ! *** NO NTRK=I IF(ZBOT.LT.TRACK(I,2,L))ZBOT=TRACK(I,2,L) 22 CONTINUE 23 CONTINUE ! ! *** ERROR DETECTED? IF(NTRK.EQ.1.OR.NTRK.GT.MXTRK) THEN ! *** YES. ITEMP=MXTRK WRITE(NPU,29) ITEMP STOP ELSE ! *** NO ERROR DETECTED. ! *** EXTEND LAST DEPTH BEYOND MAX RANGE. TRACK(NTRK+1,1,L)=1.0E+38 TRACK(NTRK+1,2,L)=TRACK(NTRK,2,L) R2(L)=TRACK(1,1,L) Z2(L)=TRACK(1,2,L) ENDIF 24 CONTINUE ! ELSE ! *** USER WRITTEN SUBROUTINE UBOTTOM SUPPLIES PROFILE. CALL UBOTTOM ENDIF ! ! *** FIND BOTTOM SEGMENT WHICH CONTAINS STARTING RANGE DO 31 L=1,NSOL ITRK(L)=1 25 R1(L)=R2(L) Z1(L)=Z2(L) ITRK(L)=ITRK(L)+1 ! IF(ITRK(L).GT.MXTRK) THEN ITEMP=MXTRK WRITE(NPU,29) ITEMP 29 FORMAT(1X,'ERROR. BOTTOM PROFILE MISSING OR ARRAY OVERFLOW. ', & 'MAX IS ',I5, ',. CHECK PARAMETER MXTRK.') STOP ELSE ENDIF ! R2(L)=TRACK(ITRK(L),1,L) Z2(L)=TRACK(ITRK(L),2,L) ! ! *** ADVANCE TRACK IF NECESSARY SO THAT R1.LE.RA.LT.R2 IF(RA.GE.R2(L)) GO TO 25 ! ! *** BE SURE THAT R1 LE RA IF(R1(L).GT.RA) THEN WRITE(NPU,*)'DEPTH OF BOTTOM AT INITIAL RANGE MISSING,' WRITE(NPU,*)'OR PARAMETER MXTRK SHOULD BE INCREASED' STOP ELSE ENDIF ! ! *** COMPUTE SLOPE OF BOTTOM IN RANGE THETA(L)=ATAN2(Z2(L)-Z1(L),R2(L)-R1(L)) ! ! *** END OF DO LOOP 31 CONTINUE ! ! *** READ RANGE OF SVP READ(NIU,*,END=33) RSVP ! ! *** SVP BEYOND START RANGE? IF(RSVP.GT.RA) GO TO 33 ! ! *** NO. DETERMINE IF SVP IN RUNSTREAM OR SUPPLIED BY SUBROUTINE USVP3D. 32 READ(NIU,*,END=33) KSVP IF(KSVP.EQ.0) CALL SVP3D IF(KSVP.NE.0) CALL USVP3D ! ! *** ERROR IN SVP? IF(NSVP.NE.0) GO TO 35 ! *** YES. 33 WRITE(NPU,34) RA 34 FORMAT(1X,'SVP MISSING OR INPUT ERROR. RANGE = ',F8.1,' M.') STOP ! 35 CONTINUE ! *** NO ERROR DETECTED. ! *** IF C0 NOT SPECIFIED, SET C0 TO AVERAGE SPEED OF FIRST PROFILES ! *** IN FIRST LAYERS AT INITIAL RANGE ! IF(C0.EQ.0.0) THEN ZWSVP=0.0 DO 40 L=1,NSOL NWSVP=IXSVP(1,L) ZWSVP=ZWSVP+ZSVP(NWSVP,L) DO 40 I=2,NWSVP C0=C0+(ZSVP(I,L)-ZSVP(I-1,L))*(CSVP(I-1,L)+.5*(CSVP(I,L)- & CSVP(I-1,L))) 40 CONTINUE C0=C0/ZWSVP ELSE ENDIF ! ! *** COMPUTE REFERENCE WAVE NUMBER XK0=2.0*PI*FRQ/C0 ! ! *** COMPUTE ATTENUATION - SACLANT MEMO SM-121 (JENSEN + FERLA) ! *** MODIFIED AS FOLLOWS: ! *** IF INPUTTED BETA IS LT 0.0, ALPHA IS COMPUTED IN DB/METER ! *** AND USED TO COMPUTE BETA IN DB/WAVELENGTH ROUTINE INDX3D. ! ALPHA=FRQ*FRQ*(.007+(.155*1.7)/(1.7*1.7+FRQ*FRQ*.000001))*1.0E-09 ! ! *** ADJUST LAYER DEPTHS IN CASE BOTTOM SLOPES AND RA.GT.RSVP. ! *** ASSUMES LAYERS FOLLOW BOTTOM CONTOUR. ! DO 50 L=1,NSOL NLYR=NLYRS(L) DO 50 ILYR=1,NLYR ZLYR(ILYR,L)=ZLYR(ILYR,L)+(RA-RSVP)*TAN(THETA(L)) 50 CONTINUE ! ! *** GET RANGE OF NEXT SVP ! *** IF END, NO MORE PROFILES FOLLOW. ! 53 READ(NIU,*,END=55) RSVP ! ! *** TEST FOR 3D INPUT RUNSTREAM - USED DURING DEBUG PHASE. ! *** WHEN COMPARING 2D SOLUTION WITH AZIMUTH INDEPENDENT NX2D SOLUTIONS. ! *** SEE NOTE 2 REGARDING -1,-1 IN DESCRIPTION OF INPUT PARAMETERS. ! IF(NDIM.EQ.1.AND.RSVP.LT.0.0)GO TO 53 GO TO 56 ! ! *** ONLY ONE SVP - SET RSVP LARGE SO SAME SVP USED FOR ENTIRE PROBLEM 55 RSVP=1.0E+38 ! 56 CONTINUE ! ! *** IF STARTING FIELD IS BEYOND NEXT SVP, GO BACK AND GET NEXT SVP IF(RSVP.LE.RA) GO TO 32 ! ! *** ABSORBING LAYER REQUESTED? IF(ITYPEB.EQ.3) THEN ! *** YES. ! *** EXTEND BOTTOM 4/3 MAX DEPTH IF ZA = ZERO IF(ZA.EQ.0.0) ZA=4.0*ZBOT/3.0 ! *** ABSORBING LAYER DEEP ENOUGH? IF(ZA.LT.ZBOT) THEN ! *** NO. EXTEND BOTTOM 4/3 MAX DEPTH. WRITE(NPU,57) 57 FORMAT(1X,'ERROR. ZA RESET TO MAX DEPTH OF BOTTOM LAYER.') ZA=4.0*ZBOT/3.0 ELSE ENDIF ! ! *** INSERT PARAMETERS FOR EXTENDED BOTTOM IN APPROPRIATE ARRAYS ZMAX=ZA DO 59 L=1,NSOL NLYR=NLYRS(L) ZG=ZLYR(NLYR,L) NSVP=IXSVP(NLYR,L) IF(NLYR.GT.1)ZG=ZG-ZLYR(NLYR-1,L) NLYRS(L)=NLYRS(L)+1 NLYR=NLYRS(L) ! *** STORE DEPTH OF ARTIFICIAL ABSORBING LAYER ZLYR(NLYR,L)=ZA ! *** USE BOTTOM DENSITY AND ATTENUATION OF ABOVE LAYER RHO(NLYR,L)=RHO(NLYR-1,L)+ZG*RHOG(NLYR-1,L) RHOG(NLYR,L)=0.0 BETA(NLYR,L)=BETA(NLYR-1,L)+ZG*BETAG(NLYR-1,L) BETAG(NLYR,L)=0.0 ! *** USE BOTTOM SPEED OF ABOVE LAYER IXSVP(NLYR,L)=NSVP+2 ZSVP(NSVP+1,L)=ZSVP(NSVP,L) CSVP(NSVP+1,L)=CSVP(NSVP,L) ZSVP(NSVP+2,L)=ZLYR(NLYR,L) CSVP(NSVP+2,L)=CSVP(NSVP+1,L) 59 CONTINUE ELSE ! *** IF BOTTOM NOT EXTENDED AND ZA=0, SET ZA TO MAX DEPTH INPUTTED ! *** WITH FIRST SVP. ASSUMES DEPTH IS CONSTANT IN AZIMUTH. IF(ZA.EQ.0.0) ZA=ZLYR(NLYR,1) ENDIF ! ! *** IF N NOT SPECIFIED, SET N SO THAT DZ .LE. 1/10 WAVELENGTH IF(N.EQ.0) THEN FN=((10*ZA*FRQ/C0)+1) N=INT(FN) ELSE ENDIF ! IF(N.GT.MXNZ) THEN WRITE(NPU,*)'ERROR. N TOO LARGE.' WRITE(NPU,*)'INCREASE PARAMETER MXNZ AND RECOMPILE.' STOP ELSE ENDIF ! ! *** COMPUTE RECEIVER DEPTH INCREMENT. DZ MAY BE SUCH THAT RECEIVERS DO ! *** NOT LIE EXACTLY ON LAYER INTERFACES DZ=ZA/N ! ! ADJUST N SO THAT DEPTH (N+1)*DZ IS BOTTOM IF(ITYPEB.EQ.1.OR.ITYPEB.EQ.3)N=N-1 ! ! *** COMPUTE TOTAL SOLUTION DEPTHS NTOT=N*NSOL ! ! *** IF RANGE STEP NOT SPECIFIED, SET IT EQUAL TO 1/10 WAVELENGTH. ! *** THIS STEP SIZE IS ARBITRARY. ! IF(DR.EQ.0.0) DR=.1*C0/FRQ ! ! *** GET STARTING FIELD CALL SFLD3D ! ! *** COMPUTE DEPTH WRITE INCREMENT TO NEAREST DZ IF(WZ1.LT.DZ)WZ1=DZ RIWZ=(WZ1/DZ+.5) IWZ1=INT(RIWZ) ! IF(WZ2.LT.DZ)WZ2=DZ IF(ITYPEB.EQ.1.OR.ITYPEB.EQ.3.AND.WZ2.GT.N*DZ)WZ2=N*DZ RIWZ=(WZ2/DZ+.5) IWZ2=INT(RIWZ) ! IF(WDZ.LT.DZ)WDZ=WZ1 RIWZ=(WDZ/DZ+.5) IWZ=INT(RIWZ) WDZ=IWZ*DZ ! NZ=(IWZ2-IWZ1)/IWZ+1 ! ! *** COMPUTE DEPTH PRINT INCREMENT TO NEAREST DZ RIPZ=(PDZ/DZ+.5) IPZ=INT(RIPZ) IF(PDZ.GT.0.0.AND.IPZ.EQ.0) IPZ=1 PDZ=IPZ*DZ ! ! *** PRINT INPUT PARAMETERS CALL PRINTP ! 84 FORMAT(3X,'LAYER ',I3,3X,'DEPTH = ',F8.1,3X, & 'DENSITY = ',F6.2,3X,'GRADIENT = ',F7.4,/,3X,'BETA = ',F8.3, & 3X,'GRADIENT = ',F7.4) 93 FORMAT(//1X,'SOUND VELOCITY PROFILE AT RANGE ',F10.1,' METERS',/) 94 FORMAT(1X,I4,3X,F8.1,3X,F8.1) ! ! *** OUTPUT REQUESTED? IF(WDR.NE.0.0) THEN ! *** YES. WRITE SELECTED PARAMETERS TO OUTPUT FILE. ! CALL ASSIGN(NOU,'FOR3D.OUT') !ljh OPEN(UNIT=NOU,FILE='FOR3D.OUT') WRITE(NOU,994)NDIM,FRQ,ZS,C0,ISF,RA,ZA,N,IHNK,ITYPES,ITYPEB, & ITYPPW,ITYPSW,FLDW,NSEC,NSOL,RMAX,DR,WDR,WZ1,WZ2,WDZ,DZ, & DOUGRA,NDIV,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10,U11,U12 994 FORMAT(1X,I1,3F10.2,I3,2F8.0,I5,5I2,F9.3,/,2I5,F10.1,F5.1, & 5F7.1,F10.1,I4,/,(5F15.7)) ! ! *** WRITE STARTING FIELD TO OUTPUT FILE. DO 97 J=1,NSOL,NSOLW M=(J-1)*N IF(NDIM.NE.1)ANG=-FLDW/2.0+((J-1)*FLDW/NSEC) WRITE(NOU,995) ANG,NZ,RA,WZ1,WDZ,(U(M+I),I=IWZ1,IWZ2,IWZ) 995 FORMAT(1X,F10.5,I6,F8.0,F8.1,F8.1,/,2(E15.7,E15.7)) 97 CONTINUE ELSE ! *** NO OUTPUT REQUESTED. ENDIF !cljh ! ! ! Write SHDFil header !ljh SHDFileName = 'for3d.shd' Ntheta = NSOL do it = 1, Ntheta thetaSHD(it) = -FLDW/2.0+((it-1)*FLDW/NSEC) end do sx = 0 sy = 0 sd = ZS Nsx = 1 Nsy = 1 Nsy = 1 Nsd = 1 Nrd = (ZA /WDZ) + 1 Nr = (((RMAX - R0) + 1) /WDR) + 1 rd(1) = WDZ do id = 2, Nrd rd(id) = rd(id-1) + WDZ end do rr(1) = R0 do ir = 2, Nr rr(ir) = rr(ir-1) + WDR end do freq = FRQ atten = ALPHA !ljh PlotType = ' ' write(6,*) 'sd, Ntheta, Nr, Nrd =', sd, Ntheta, Nr, Nrd write(6,*) 'thetas' write(6,*) thetaSHD(1:Ntheta) write(6,*) 'Receiver Depths' write(6,*) rd(1:Nrd) write(6,*) 'Receiver Ranges' write(6,*) rr(1:Nr) CALL WriteHeader( 'for3d.shd', TITLE, thetaSHD, & Ntheta, sx, Nsx, sy, Nsy, sd, Nsd, rd, Nrd, rr, & Nr, freq, atten, ' ' ) ALLOCATE( P( Ntheta, Nr, Nrd ) ) ! itheta = 0 ir = 0 ! ! !cljh ! ! *** INITIALIZE RANGE VARIABLE AT WHICH SOLUTION IS TO BE PRINTED XPR=RA+PDR ! ! *** INITIALIZE RANGE VARIABLE AT WHICH SOLUTION IS TO BE RECORDED XWR=RA+WDR IF(WDR.EQ.0.0) XWR=RA+RMAX+1.0 ! ! *** SAVE RANGE STEP OLDR=DR ! N1=1 ! ! ****************************************************************** ! ****************************************************************** ! *** MAIN LOOP STARTS HERE *** ! *** SOLUTION WILL BE ADVANCED FROM RANGE RA TO RANGE RA+DR *** ! ****************************************************************** ! ****************************************************************** ! 100 CONTINUE DRLAST=DR ! ! *** COMPUTE INDEX OF REFRACTION TABLE AT PRESENT RANGE CALL INDX3D ! ! *** THIS SECTION OF CODE DETERMINES WHETHER OR NOT TO RESTORE STEP ! *** SIZE TO THE ORIGINAL DR. ! !ljh Well, this looks important! How far does "Next Section" extend??? ! *** COMMENT NEXT SECTION IF INPUT DR IS TO BE USED DO 101 L=1,NSOL IF(THETA(L).NE.0.0) GO TO 102 101 CONTINUE ! ! *** FLAT BOTTOM - RESTORE DR IF IT HAS BEEN CHANGED DR=OLDR ! 102 CONTINUE ! DO 129 L=1,NSOL DRMIN(L)=DR ! *** DOES SVP CHANGE BEFORE BOTTOM PROFILE? IF(RSVP.GE.R2(L)) GO TO 105 ! *** YES. DOES SVP CHANGE BEFORE NEXT SOLUTION RANGE? IF(RA+DRMIN(L).LE.RSVP) GO TO 126 ! *** YES. ADJUST DR SO THAT SOLUTION ADVANCES TO RSVP DRMIN(L)=RSVP-RA GO TO 126 ! *** BOTTOM PROFILE CHANGES BEFORE OR AT SAME RANGE AS SVP 105 CONTINUE ! ! *** DETERMINE IF BOTTOM DEPTHS ARE TO BE UPDATED ! *** FIRST PASS - RA WILL BE .LT. R2 ! ! *** UPDATE BOTTOM DEPTHS? IF(RA.GE.R2(L)) GO TO 110 ! ! *** NO. ! *** DETERMINE IF RANGE STEP TOO LARGE . IF(RA+DRMIN(L).LE.R2(L)) GO TO 126 ! ! *** RANGE STEP IS TOO LARGE - RESET DR - ADV SOLUTION TO R2 DRMIN(L)=R2(L)-RA GO TO 126 ! ! *** UPDATE BOTTOM DEPTHS 110 CONTINUE R1(L)=R2(L) Z1(L)=Z2(L) ITRK(L)=ITRK(L)+1 R2(L)=TRACK(ITRK(L),1,L) Z2(L)=TRACK(ITRK(L),2,L) ! ! *** TWO DEPTHS AT SAME RANGE INDICATE VERTICAL DISCONTINUITY. ! *** ADVANCE TRACK FORWARD. IF(R1(L).EQ.R2(L)) GO TO 110 ! ! *** RESTORE DR IN CASE BOTTOM IS FLAT DRMIN(L)=OLDR ! ! *** COMPUTE SLOPE OF BOTTOM THETA(L)=ATAN2(Z2(L)-Z1(L),R2(L)-R1(L)) ! ! *** IF BOTTOM IS NOT FLAT, COMPUTE NEW RANGE STEP IF(THETA(L).EQ.0) GO TO 112 ! ! *** COMPUTE NEW RANGE STEP. DR MAY BECOME VERY SMALL. DRMIN(L)=ABS(DZ*COS(THETA(L))/SIN(THETA(L))) ! ! *** IF RANGE STEP TOO LARGE, PRINT WARNING MESSAGE. RECOMPUTE DR. 112 IF(RA+DRMIN(L).LE.R2(L)) GO TO 120 ! WRITE(NPU,115)RA,L,DRMIN(L) !ljh 115 FORMAT(1X,'AT RANGE ',F7.0,', ISOL ',I3,', RANGE STEP ',F5.1, !ljh +' TOO LARGE FOR BOTTOM CONTOUR') ! STOP !..................................... ! ! *** ADJUST DR AND CONTINUE. DRMIN(L)=R2(L)-RA ! WRITE(NPU,*)'DR ADJUSTED TO ',DRMIN(L) ! 120 CONTINUE IF(RA+DRMIN(L).GT.RSVP) DRMIN(L)=RSVP-RA ! ! *** PRINT BOTTOM DEPTHS IF(IBOT.EQ.0) GO TO 126 WRITE(NPU,93) RA IF(NDIM.NE.1)ANG=-FLDW/2.0+((L-1)*FLDW/NSEC) ! WRITE(NPU,*)' ' WRITE(NPU,*)' RELATIVE BEARING = ',ANG,' DEGREES' ! WRITE(NPU,*)' ' TH=THETA(L)*DEG WRITE(NPU,123) R1(L),Z1(L),R2(L),Z2(L),TH 123 FORMAT(//1X,'BOTTOM PROFILE ',//,1X, & 'R1 = ',F10.1,' M',/,1X, & 'Z1 = ',F8.1,' M',/,1X, & 'R2 = ',F10.1,' M',/,1X, & 'Z2 = ',F8.1,' M',/,1X, & 'SLOPE = ',F8.1,' DEG') 126 CONTINUE ! 129 CONTINUE ! *** END OF LOOP ! ! *** SELECT SMALLEST DR WITHIN FOLLOWING ARBITRARY LIMITS. ! *** MAXIMUM ALLOWED IS INPUT DR. ! *** MINIMUM ALLOWED IS 0.5*(DZ*DZ). DR=DRMIN(1) DO 128 L=2,NSOL IF(DR.GT.DRMIN(L))DR=DRMIN(L) 128 CONTINUE IF(DR.LT.0.5*(DZ*DZ))THEN DR=.5*DZ*DZ ! WRITE(NPU,*)'AT RANGE ',RA, ' DR RESET TO ',DR ELSE ENDIF IF(DR.GT.OLDR)THEN DR=OLDR ! WRITE(NPU,*)'AT RANGE ',RA, ' DR RESET TO ',DR ELSE ENDIF ! ! IF(DR.NE.DRLAST)WRITE(NPU,*)'AT RANGE ',RA,' DR WAS CHANGED FROM & ! ',DRLAST,' TO ',DR IF(DR.NE.DRLAST)WRITE(NPU,*)'AT RANGE ',RA,' DR WAS CHANGED FROM ',DRLAST,' TO ',DR ! ! *** COMPUTE MATRIX A AT PRESENT RANGE CALL AMIFD3 ! ! *** IF 3D RUN, COMPUTE B MATRIX. IF(NDIM.EQ.3) CALL BMIFD3 ! ! *** ADVANCE RANGE ONE STEP RA=RA+DR NOLD=N ! ! *** COMPUTE RIGHT HAND SIDE CALL RHS ! ! *** MARCH SOLUTION 1 DR FORWARD CALL TWOSTEP ! ! *** APPLY ABSORPTION IF ITYPEB = 3 IF(ITYPEB.EQ.3) THEN DO 180 J=1,NSOL NLYR=NLYRS(J) RMM=(ZLYR(NLYR-1,J)/DZ) MM=INT(RMM) NA=N-MM IF(NA.GT.0) GO TO 175 IF(NA.GE.-1) GO TO 180 WRITE(NPU,174) RA 174 FORMAT(1X,'ERR IN THICKNESS OF ABSORBING LAYER AT ',F8.1,' M') STOP ! 175 CONTINUE M=(J-1)*N DO 179 I=1,NA ATT=EXP(-.01*DR*EXP(-((I-NA)/(NA/3.0))**2)) U(M+MM+I)=U(M+MM+I)*ATT IF(J.GT.1)GO TO 176 IF(NDIM.EQ.3.AND.ITYPPW.EQ.2)PORTY(MM+I)=PORTY(MM+I)*ATT 176 IF(J.LT.NSOL)GO TO 179 IF(NDIM.EQ.3.AND.ITYPSW.EQ.2)STBDY(MM+I)=STBDY(MM+I)*ATT 179 CONTINUE 180 CONTINUE ELSE ENDIF ! ! *** IF SOLUTION IS TO BE WRITTEN ON DISK, ! *** WRITE SELECTED PRESSURE FIELD AT RANGE RA IF(XWR.GT.RA) GO TO 260 ! IF(WZ1.LT.DZ)WZ1=DZ RIWZ=(WZ1/DZ+.5) IWZ1=INT(RIWZ) ! IF(WZ2.LT.DZ)WZ2=DZ IF(ITYPEB.EQ.1.OR.ITYPEB.EQ.3.AND.WZ2.GT.N*DZ)WZ2=N*DZ RIWZ=(WZ2/DZ+.5) IWZ2=INT(RIWZ) IF(WDZ.LT.DZ)WDZ=WZ1 ! RIWZ=(WDZ/DZ+.5) IWZ=INT(RIWZ) WDZ=IWZ*DZ ! NZ=(IWZ2-IWZ1)/IWZ+1 ! ! !cljh !cljh Save pressure U in variable P ! ir = ir + 1 ! The big do-loop is by range step do itheta = 1, Ntheta ! Ntheta = NSOL M=(itheta-1)*N ! See lines 1230-1232 ird = 0 do I = M+IWZ1,M+IWZ2,IWZ ird = ird + 1 P(itheta, ir, ird) = U(I) end do end do ! ! !cljh !cljh DO 259 J=1,NSOL,NSOLW M=(J-1)*N IF(NDIM.NE.1)ANG=-FLDW/2.0+((J-1)*FLDW/NSEC) WRITE(NOU,995) ANG,NZ,RA,WZ1,WDZ,(U(M+I),I=IWZ1,IWZ2,IWZ) 259 CONTINUE ! *** DETERMINE NEXT RANGE AT WHICH TO WRITE SOLUTION ON DISK XWR=XWR+WDR 260 CONTINUE ! ! *** DETERMINE IF SOLUTION IS TO BE PRINTED IF(XPR.GT.RA.OR.IPZ.EQ.0.OR.PDR.EQ.0.0) GO TO 350 ! ! *** PRINT RANGE RTEMP=RA/1000.0 ! OPEN(UNIT=NOU,FILE='FOR3D.OUT') ! *** COMPUTE AND PRINT PROPAGATION LOSS AT EACH IPZ'TH DEPTH 275 FORMAT(2X,'RECEIVER',4X,'DEPTH(M)',4X,'LOSS(DB)',14X,'U(I)') JSOL=NSOL/2 DO 300 J=1,NSOL,NSOLP M=(J-1)*N IF(NDIM.NE.1)ANG=-FLDW*PI*(.5*NSEC-(J-1))/(180.0*NSEC) WRITE(NPU,277)J,RA,ANG*DEG 277 FORMAT(/,2X,'SOLUTION',I4,', RANGE = ',F8.1,' M', & ', RELATIVE BEARING = ',F6.1, & ' DEGREES',/) WRITE(NPU,275) UA=0 DO 295 I=IPZ,N,IPZ ZI=I*DZ UA=U(M+I) IF(IHNK.EQ.1)UA=UA*HNKL(XK0*RA) !ljh PL=CABS(UA) PL=ABS(UA) ! cabs no longer in use IF(PL.LE.0.0) GO TO 295 !ljh PL=-20.0*ALOG10(PL) PL=-20.0*DLOG10(PL) ! dlog10(real*8) IF(IHNK.EQ.0) PL=PL+10.0*ALOG10(RA) GO TO 289 PL=999.9 289 CONTINUE WRITE(NPU,294) I,ZI,PL,UA 294 FORMAT(2X,I5,(3X,F10.2,3X,F10.3),3X,'(',E12.5,2X,E12.5,' )') !ljh Another important bit of code to change, hidden deep in the program! ! *** UNCOMMENT IF EXACT SOLUTION SUBROUTINE AVAILABLE. ! CALL UEXACT ! WRITE(NPU,*)' ' 295 CONTINUE 300 CONTINUE ! ! *** DETERMINE NEXT RANGE AT WHICH TO PRINT SOLUTION XPR=XPR+PDR ! 350 CONTINUE ! ! *** SWITCH TO DOUGLAS METHOD? IF(DOUGRA.NE.0.0.AND.RA.GE.DOUGRA)THEN ! *** YES. OLDZ=DZ NOLD=N NTOLD=NTOT IF(ITYPEB.EQ.1)N=N+1 IF(ITYPEB.EQ.3)N=N+1 ! ! *** DIVIDE N BY NDIV AND MAKE NECESSARY ADJUSTMENTS. N=INT(N/NDIV+.5) DZ=ZA/N IF(ITYPEB.EQ.1)N=N-1 IF(ITYPEB.EQ.3)N=N-1 NTOT=N*NSOL RIPZ=(PDZ/DZ+.5) IPZ=INT(RIPZ) ! ! *** ADJUST SOLUTION FIELD. DO 360 J=1,NSOL MOLD=(J-1)*NOLD M=(J-1)*N DO 359 I=1,N ZI=I*DZ ! *** AVOID ROUNDOFF PROBLEM. RK=ZI/OLDZ+.5 K=INT(RK) U(M+I)=U(MOLD+K)+(ZI-K*OLDZ)*(U(MOLD+K+1)-U(MOLD+K))/OLDZ 359 CONTINUE 360 CONTINUE ! ! *** ADJUST WALL CONDITIONS. DO 365 I=1,N ZI=I*DZ ! AVOID ROUNDOFF PROBLEM. RK=ZI/OLDZ+.5 K=INT(RK) PORTY(I)=PORTY(K) STBDY(I)=STBDY(K) 365 CONTINUE ! ! *** SET ALFA FOR DOUGLAS METHOD. ALFA=1.0/12.0 ! *** SET DOUGRA LARGE SO THAT WE DON'T ENTER THIS CODE AGAIN. DOUGRA=1.0E+38 ! ! WRITE(NPU,*)' ' WRITE(NPU,*)'AT RANGE ',RA,' METERS. SWITCH FROM CRANK-NICOLSON TO DOUGLAS METHOD.' WRITE(NPU,*)'CRANK-NICOLSON N ,DZ, NTOT ',NOLD,OLDZ,NTOLD WRITE(NPU,*)'DOUGLAS N ,DZ, NTOT ',N, DZ,NTOT !CCC CALL SECOND(CPT2) !ljh WRITE(NPU,370)CPT2-CPT1 !ljh 370 FORMAT(1X,'ELAPSED TIME = ',F12.2,' CRAY XMP 2/8 CPU SECONDS.') ! WRITE(NPU,*)' ' ! ! *** RECOMPUTE DEPTH WRITE INCREMENT TO NEAREST DZ IF(WZ1.LT.DZ)WZ1=DZ RIWZ=(WZ1/DZ+.5) IWZ1=INT(RIWZ) ! IF(WZ2.LT.DZ)WZ2=DZ RIWZ=(WZ2/DZ+.5) IWZ2=INT(RIWZ) ! IF(WDZ.LT.DZ)WDZ=WZ1 RIWZ=(WDZ/DZ+.5) IWZ=INT(RIWZ) WDZ=IWZ*DZ ! NZ=(IWZ2-IWZ1)/IWZ+1 ! ! *** RECOMPUTE DEPTH PRINT INCREMENT TO NEAREST DZ RIPZ=(PDZ/DZ+.5) IPZ=INT(RIPZ) IF(PDZ.GT.0.0.AND.IPZ.EQ.0) IPZ=1 PDZ=IPZ*DZ ! ELSE ENDIF ! End of If SWITCH TO DOUGLAS METHOD?, line 1211 ! ! *** TERMINATE RUN IF SOLUTION AT MAXIMUM RANGE HAS BEEN OBTAINED IF(RA.LT.RMAX)THEN ! ! *** GET ENVIRONMENT FOR NEXT RANGE ! *** READ NEW SVP PROFILE FLAG IF(RA.GE.RSVP) READ(NIU,*,END=33) KSVP ! ! *** KSVP MAY HAVE BEEN SET IN USER ROUTINE USVP3D. ! *** IF KSVP NOT 0, SUBROUTINE USVP3D IS CALLED. USER IS ! *** RESPONSIBLE FOR ADJUSTING DEPTHS OF LAYERS. IF(KSVP.EQ.0) GO TO 534 CALL USVP3D GO TO 548 ! ! *** NEW SVP AT ADVANCED RANGE? 534 IF(RA.GE.RSVP) GO TO 545 ! !ljh 535 CONTINUE ! ! *** UPDATE DEPTHS OF LAYERS ! *** ASSUMES LAYERS ARE PARALLEL AND FOLLOW BOTTOM CONTOUR ! *** IF ITYPEB = 3, BOTTOM OF ARTIFICIAL LAYER REMAINS FLAT. DO 537 L=1,NSOL DZZ=DR*TAN(THETA(L)) ! *** IS BOTTOM FLAT? IF(THETA(L).EQ.0.0) GO TO 537 NLYR=NLYRS(L) MLYR=NLYRS(L) IF(ITYPEB.EQ.3) MLYR=NLYRS(L)-1 DO 536 ILYR=1,MLYR ZLYR(ILYR,L)=ZLYR(ILYR,L)+DZZ IF(ITYPEB.EQ.3.AND.ZLYR(NLYR,L).LT.ZLYR(MLYR,L))THEN WRITE(NPU,*)'ERROR. ARTIFICIAL LAYER NOT DEEP ENOUGH.' STOP ELSE ENDIF 536 CONTINUE 537 CONTINUE ! ! *** ADJUST DEPTHS OF PROFILES IN SLOPING LAYERS ! *** ASSUMES SAME SVP IN LAYERS DO 540 L=1,NSOL DZZ=DR*TAN(THETA(L)) NLYR=NLYRS(L) NWSVP=IXSVP(1,L) NSVP=IXSVP(NLYR,L) DO 540 I=NWSVP+1,NSVP ZSVP(I,L)=ZSVP(I,L)+DZZ 540 CONTINUE GO TO 570 ! ! *** GET NEXT SVP FROM INPUT RUNSTREAM. 545 CONTINUE CALL SVP3D ! 548 CONTINUE ! ! *** ERROR DETECTED? IF(NSVP.EQ.0) GO TO 33 ! ! *** NO. READ RANGE OF NEXT PROFILE? IF(RA.GE.RSVP) THEN READ(NIU,*,END=549) RSVP ! *** CHECK FOR 3D RUNSTREAM IF(NDIM.EQ.1.AND.RSVP.LT.0.0)READ(NIU,*,END=549) RSVP ELSE ENDIF GO TO 550 ! 549 CONTINUE ! *** SET RSVP LARGE SO THAT LAST PROFILE IS USED FOR REMAINDER OF PROBLEM RSVP=1.0E+38 ! 550 CONTINUE ! *** ARTIFICIAL ABSORBING LAYER? IF(ITYPEB.EQ.3) THEN ! *** YES. UPDATE DENSITY, ATTEN AND SPEED. DO 551 L=1,NSOL NLYR=NLYRS(L) ZG=ZLYR(NLYR,L) NSVP=IXSVP(NLYR,L) IF(ZA.LT.ZLYR(NLYR,L)) GO TO 551 IF(NLYR.GT.1)ZG=ZLYR(NLYR,L)-ZLYR(NLYR-1,L) NLYRS(L)=NLYRS(L)+1 NLYR=NLYRS(L) ZLYR(NLYR,L)=ZA RHO(NLYR,L)=RHO(NLYR-1,L)+ZG*RHOG(NLYR-1,L) RHOG(NLYR,L)=0.0 BETA(NLYR,L)=BETA(NLYR-1,L)+ZG*BETAG(NLYR-1,L) BETAG(NLYR,L)=0.0 IXSVP(NLYR,L)=NSVP+2 ZSVP(NSVP+1,L)=ZSVP(NSVP,L) CSVP(NSVP+1,L)=CSVP(NSVP,L) ZSVP(NSVP+2,L)=ZLYR(NLYR,L) CSVP(NSVP+2,L)=CSVP(NSVP+1,L) 551 CONTINUE ELSE ENDIF ! ! *** PRINT SVP IF(ISVP.NE.0) THEN WRITE(NPU,93) RA DO 562 L=1,NSOL IF(NDIM.NE.1)ANG=-FLDW/2.0+((L-1)*FLDW/NSEC) ! WRITE(NPU,*)' ' WRITE(NPU,*)' RELATIVE BEARING = ',ANG,' DEGREES' ! WRITE(NPU,*)' ' K=0 NLYR=NLYRS(L) DO 561 ILYR=1,NLYR WRITE(NPU,84)ILYR,ZLYR(ILYR,L),RHO(ILYR,L),RHOG(ILYR,L), & BETA(ILYR,L),BETAG(ILYR,L) ! WRITE(NPU,*)' ' J=K+1 K=IXSVP(ILYR,L) DO 560 I=J,K WRITE(NPU,94) I,ZSVP(I,L),CSVP(I,L) 560 CONTINUE ! WRITE(NPU,*)' ' 561 CONTINUE 562 CONTINUE ELSE ENDIF 570 CONTINUE ! GO TO 100 ELSE !cljh go to 588 !cljh ENDIF ! IF RA .LT. RMAX !cljh ! ! Write out the field 588 DO ird = 1, Nrd DO itheta = 1, Ntheta iRec = 9 + ( itheta - 1 ) * Nrd + ird WRITE( SHDFile, REC = iRec ) P( itheta, 1 : Nr, ird ) END DO END DO ! !cljh ! *** TERMINATE RUN !ljh 600 CONTINUE !ljh IF(WDR.NE.0.0)CALL CLOSE(NOU) IF(WDR.NE.0.0) CLOSE(NOU) !CC CALL SECOND(CPT2) !CC WRITE(NPU,601)CPT2-CPT1 !ljh 601 FORMAT(1X,'END OF RUN. CRAY XMP 2/8 CPU TIME = ',F12.2, !ljh +' SECONDS.') STOP END ! ********** ! * AMIFD3 * ! ********** SUBROUTINE AMIFD3 ! ****************************************************************** ! ****************************************************************** ! * COMPUTE DIAGONALS FOR 'A' MATRIX * ! ****************************************************************** !ljh INCLUDE 'FOR3D.CMN' USE commons IMPLICIT REAL (KIND=8) (A-H,O-Z) IMPLICIT INTEGER (KIND=4) (I-N) COMPLEX T1,T3,XN12,XNSUM ! *** RA HAS BEEN INCREMENTED BY DR IN MAIN PROGRAM T1=CMPLX(.25,-.25*XK0*DR) T2=1.0/(XK0*XK0*DZ*DZ) T3=T1*T2 DO 100 L=1,NSOL M=(L-1)*N DO 50 I=1,N XN12=T1*((XN1(M+I)-1)/RHO1(M+I)+(XN2(M+I)-1)/RHO2(M+I)) R12=1.0/RHO1(M+I)+1.0/RHO2(M+I) RSUM=RHO1(M+I)+RHO2(M+I) XNSUM=T1*(XN1(M+I)+XN2(M+I)-2.0) !** UPPER DIAGONAL AU(M+I)=(T3+ALFA+ALFA*T1*(XN2(M+I)-1))/RHO2(M+I) IF(I.EQ.N) AU(M+I)=0.0 !** MAIN DIAGONAL AM(M+I)=(2.0+XNSUM)/RSUM-(ALFA+T3)*R12-ALFA*XN12 !** LOWER DIAGONAL AL(M+I)=(T3+ALFA+ALFA*T1*(XN1(M+I)-1))/RHO1(M+I) IF(I.EQ.1)AL(M+I)=0.0 50 CONTINUE 100 CONTINUE RETURN END ! ********** ! * BCON3D * ! ********** SUBROUTINE BCON3D ! ***************************************************************** ! *** BOTTOM BOUNDARY CONDITION SUBROUTINE ! ***************************************************************** !ljh INCLUDE 'FOR3D.CMN' USE commons IMPLICIT REAL (KIND=8) (A-H,O-Z) IMPLICIT INTEGER (KIND=4) (I-N) GO TO (100,200,300,400,500),ITYPEB+1 WRITE(NPU,50)ITYPEB 50 FORMAT(1X,'ERROR - BOTTOM BOUNDARY CONDITION ',I2, & ' NOT IMPLEMENTED') IERR=1 RETURN ! 100 CONTINUE ! *** ITYPEB=0 ! *** PRESSURE RELEASE BOTTOM DO 150 L=1,NSOL+2 BOTY(L)=0.0 BOTX(L)=0.0 150 CONTINUE RETURN ! 200 CONTINUE ! *** ITYPEB=1 ! *** USER BOTTOM BOUNDARY CONDITION - USER WRITES SUBROUTINE UBCON3D CALL UBCON3D RETURN ! 300 CONTINUE ! *** ITYPEB=2 IERR=1 WRITE(NPU,50)ITYPEB RETURN ! 400 CONTINUE ! *** ITYPEB=3 ! *** ARTIFICIAL ABSORBING BOTTOM - FLAT DO 450 L=1,NSOL+2 BOTY(L)=0.0 BOTX(L)=0.0 450 CONTINUE RETURN ! 500 CONTINUE ! *** ITYPEB=4 IERR=1 WRITE(NPU,50)ITYPEB RETURN END ! ********** ! * BMIFD3 * ! ********** SUBROUTINE BMIFD3 ! ****************************************************************** ! * COMPUTE DIAGONALS FOR B MATRIX * ! ****************************************************************** !ljh INCLUDE 'FOR3D.CMN' USE commons IMPLICIT REAL (KIND=8) (A-H,O-Z) IMPLICIT INTEGER (KIND=4) (I-N) !ljh COMPLEX T1,T2 COMPLEX T1 T1=-CMPLX(0.0,XK0*DR/(4.0*XK0*XK0*(RA+.5*DR)*(RA+.5*DR)*PHI*PHI)) DO 100 I=1,N M=(I-1)*NSOL DO 50 L=1,NSOL LM1=L-1 IF(LM1.EQ.0)LM1=1 LP1=L IF(LP1.GT.NSOL)LP1=NSOL J=(LM1-1)*N RH3=RHO1(J+I) J=(LP1-1)*N RH4=RHO1(J+I) R34=(1.0/RH3+1.0/RH4) ! *** UPPER DIAGONAL BU(M+L)=.5*(rh3+rh4)*(T1)/RH4 IF(L.EQ.NSOL)BU(M+L)=0.0 ! *** MAIN DIAGONAL BM(M+L)=1.0+(-T1)*.5*(rh3+rh4)*R34 ! *** LOWER DIAGONAL BL(M+L)=.5*(rh3+rh4)*(T1)/RH3 IF(L.EQ.1)BL(M+L)=0.0 50 CONTINUE 100 CONTINUE RETURN END ! ********** ! * HNKL * ! ********** COMPLEX FUNCTION HNKL(X) ! ****************************************************************** ! * HANKEL FUNCTION H0(1) - POLYNOMIAL APPROXIMATION * ! * HANDBOOK OF MATH FUNCTIONS - N.B.S. - NOV 1967 * ! ****************************************************************** ! REAL J0 DATA PI/3.141592654/ ! IF(X.GT.3.) GO TO 10 ! ! *** (-3.0.LE.X.LE.3.0) Y=X*X/9.0 J0=1.+Y*(-2.2499997+Y*(+1.2656208+Y*(-0.3163866+Y*(+0.0444479 & +Y*(-0.0039444+Y*(+0.0002100)))))) ! ! *** (0.0.LT.X.LE.3.0) Y0=2.0*LOG(0.5*X)*J0/PI+0.36746691 & +Y*(+0.60559366+Y*(-0.74350384+Y*(+0.25300117+Y*(-0.04261214 & +Y*(+0.00427916+Y*(-0.00024846)))))) HNKL=CMPLX(J0,Y0) RETURN ! ! *** (3.0.LE.X.LT.INFINITY) 10 Y=3.0/X F0= 0.79788456+Y*(-0.00000077+Y*(-0.00552740+Y*(-0.00009512 & +Y*(+0.00137237+Y*(-0.00072805+Y*(+0.00014476)))))) T0=X-0.78539816+Y*(-0.04166397+Y*(-0.00003954+Y*(+0.00262573 & +Y*(-0.00054125+Y*(-0.00029333+Y*(+0.00013558)))))) HNKL=F0*CEXP(CMPLX(0.0,T0))/SQRT(X) RETURN END ! ********** ! * INDX3D * ! ********** SUBROUTINE INDX3D ! ****************************************************************** ! * COMPUTE INDEX OF REFRACTION TABLE ! ****************************************************************** !ljh INCLUDE 'FOR3D.CMN' USE commons IMPLICIT REAL (KIND=8) (A-H,O-Z) IMPLICIT INTEGER (KIND=4) (I-N) !ljh DIMENSION COEF(4,MXSVP),WORK(MXSVP+4),WORK2(MXSVP+4) DO 200 ISOL=1,NSOL MM=(ISOL-1)*N ! *** GET INDICES, DENSITY, AND ATTENUATION FOR FIRST LAYER ILYR=1 L=1 M=IXSVP(ILYR,ISOL) RH1=RHO(ILYR,ISOL) RH1G=RHOG(ILYR,ISOL) BETA1=BETA(ILYR,ISOL) BETA1G=BETAG(ILYR,ISOL) ZG1=0.0 NLYR=NLYRS(ISOL) DO 100 I=1,N ! *** ZI IS RECEIVER DEPTH ZI=I*DZ ! *** IS RECEIVER IN THIS LAYER? 20 IF(ZI.LE.ZLYR(ILYR,ISOL)) GO TO 49 ! *** NO. ALL LAYERS CHECKED? IF(ILYR.EQ.NLYR) GO TO 52 ZG1=ZLYR(ILYR,ISOL) ! *** NO. SET INDICES, DENSITY, AND ATTENUATION FOR NEXT LAYER. ILYR=ILYR+1 L=M+1 M=IXSVP(ILYR,ISOL) RH1=RHO(ILYR,ISOL) RH1G=RHOG(ILYR,ISOL) BETA1=BETA(ILYR,ISOL) BETA1G=BETAG(ILYR,ISOL) GO TO 20 ! ! *** DEPTH ZI IS IN THIS LAYER. 49 CONTINUE ! *** DETERMINE SOUND SPEED CI AT DEPTH ZI DO 50 J=L,M IF(ZI.GT.ZSVP(J,ISOL)) GO TO 50 L=J ! *** INTERPOLATE CI=CSVP(J-1,ISOL)+(CSVP(J,ISOL)-CSVP(J-1,ISOL))* & (ZI-ZSVP(J-1,ISOL) & )/(ZSVP(J,ISOL)-ZSVP(J-1,ISOL)) GO TO 60 50 CONTINUE ! ! *** EXTRAPOLATE 52 CONTINUE IF(ZSVP(M,ISOL).NE.ZSVP(M-1,ISOL))GO TO 53 CI=CSVP(M,ISOL) GO TO 54 53 CI=CSVP(M-1,ISOL)+(CSVP(M,ISOL)-CSVP(M-1,ISOL))* & (ZI-ZSVP(M-1,ISOL))/(ZSVP(M,ISOL)-ZSVP(M-1,ISOL)) 54 CONTINUE IF(IEX.EQ.1) GO TO 60 WRITE(NPU,56) 56 FORMAT(1X,'WARNING. EXTRAPOLATION OF SVP PERFORMED.') IEX=1 60 CONTINUE C1=CI ! *** IS RECEIVER ON OR WITHIN 1 DZ OF INTERFACE? IF(ZLYR(ILYR,ISOL)-ZI.GE.DZ) GO TO 65 ! *** YES. IS THIS BOTTOM LAYER? IF(ILYR.EQ.NLYR) GO TO 65 ! *** NO. GET PARAMS FOR MEDIUM 2. ZG2=ZLYR(ILYR,ISOL) RH2=RHO(ILYR+1,ISOL) RH2G=RHOG(ILYR+1,ISOL) BETA2=BETA(ILYR+1,ISOL) BETA2G=BETAG(ILYR+1,ISOL) C2=CSVP(IXSVP(ILYR,ISOL)+1,ISOL) GO TO 70 65 CONTINUE ! *** NOT AN INTERFACE. C2=C1 RH2=RH1 RH2G=RH1G BETA2=BETA1 BETA2G=BETA1G ZG2=ZG1 70 CONTINUE RHO1(MM+I)=RH1+(ZI-ZG1)*RH1G RHO2(MM+I)=RH2+(ZI+DZ-ZG2)*RH2G BETA1=BETA1+(ZI-ZG1)*BETA1G BETA2=BETA2+(ZI+DZ-ZG2)*BETA2G XN=C0/C1 IF(BETA1.LT.0.0) BETA1=ALPHA*C1/FRQ XN1(MM+I)=CMPLX(XN*XN,XN*XN*BETA1/27.287527) XN=C0/C2 IF(BETA2.LT.0.0)BETA2=ALPHA*C2/FRQ XN2(MM+I)=CMPLX(XN*XN,XN*XN*BETA2/27.287527) 100 CONTINUE 200 CONTINUE RETURN END ! ********** ! * PORT2D * ! ********** SUBROUTINE PORT2D ! ****************************************************************** ! * SOLVE FOR FIELD AT PORT WALL BOUNDARY ! * FIELD AT ADVANCED RANGE IN PORTX ! ****************************************************************** !ljh INCLUDE 'FOR3D.CMN' USE commons IMPLICIT REAL (KIND=8) (A-H,O-Z) IMPLICIT INTEGER (KIND=4) (I-N) ! *** COMPUTE BTA(N1) THROUGH BTA(N) ! *** ARRAY BTA CONTAINS PARTIAL SOLUTION OF SYSTEM OF EQUATIONS BTA(1)=AM(1) M=2 DO 20 I=M,N BTA(I)=AM(I)-AL(I)*AU(I-1)/BTA(I-1) 20 CONTINUE DO 21 I=1,N BTA(I)=1.0/BTA(I) 21 CONTINUE !ljh TA=0.0 !ljh TB=0.0 ! *** STEP 1 - CALL TRIDIAGONAL SOLVER !ljh CALL TRID3D(AL(1),AM(1),AU(1),PORTX(1),DP(1),BTA(1),N,N, !ljh +TA,TB) CALL TRID3D(AL(1),AU(1),PORTX(1),DP(1),BTA(1),N,N) ! *** SOLUTION IS IN ARRAY PORTX RETURN END ! ********** ! * PORT3D * ! ********** SUBROUTINE PORT3D ! ****************************************************************** ! * PORT SIDEWALL BOUNDARY CONDITION ! ****************************************************************** !ljh INCLUDE 'FOR3D.CMN' USE commons IMPLICIT REAL (KIND=8) (A-H,O-Z) IMPLICIT INTEGER (KIND=4) (I-N) COMPLEX T1,T3,T4,TAL,TAM,TAU,XN12,XNSUM ! GO TO (100,200,300),ITYPPW+1 IERR=1 WRITE(NPU,50) ITYPPW 50 FORMAT(1X,'ERROR - PORT SIDEWALL BOUNDARY CONDITION ',I2, & ' NOT IMPLEMENTED') RETURN ! 100 CONTINUE ! *** ITYPPW = 0 ! *** SET FIELD AT PORT SIDEWALL TO 0.0 DO 150 I=1,N PORTX(I)=0.0 PORTY(I)=0.0 150 CONTINUE RETURN ! 200 CONTINUE ! *** ITYPPW = 1 ; USER SUPPLIES SUBROUTINE UPORT3D CALL UPORT3D RETURN ! 300 CONTINUE ! *** ITYPPW = 2 ! *** ASSUMPTION IS THAT ENVIRONMENTAL CONDITIONS AT PORT SIDEWALL ! *** BOUNDARY ARE EQUAL TO CONDITIONS AT LEFTMOST SOLUTION BOUNDARY. ! *** SOLUTION AT BOUNDARY AT ADVANCED RANGE IS SOLVED IN 2D. DELTA=XK0*DR T1=CMPLX(.25,.25*DELTA) T2=1.0/(XK0*XK0*DZ*DZ) T3=T1*T2 T4=CMPLX(.25,-.25*DELTA) RH1=RHO1(1) RH2=RHO2(1) R12=1.0/RH1+1.0/RH2 XN12=T1*((XN1(1)-1)/RH1+(XN2(1)-1)/RH2) RSUM=RH1+RH2 XNSUM=T1*(XN1(1)+XN2(1)-2.0) TAU=(T3+ALFA+ALFA*T1*(XN2(1)-1))/RH2 TAL=(T3+ALFA+ALFA*T1*(XN1(1)-1))/RH1 TAM=(2.0+XNSUM)/RSUM-(ALFA+T3)*R12-ALFA*XN12 ! *** COMPUTE RIGHT HAND SIDE D DP(1)=TAL*SURY(1)+TAM*PORTY(1)+TAU*PORTY(2)- & (T4*T2+ALFA+ALFA*T4*(XN1(1)-1))/RH1*SURX(1) DO 350 I=2,N-1 RH1=RHO1(I) RH2=RHO2(I) R12=1.0/RH1+1.0/RH2 XN12=T1*((XN1(I)-1)/RH1+(XN2(I)-1)/RH2) RSUM=RH1+RH2 XNSUM=T1*(XN1(I)+XN2(I)-2.0) TAU=(T3+ALFA+ALFA*T1*(XN2(I)-1))/RH2 TAL=(T3+ALFA+ALFA*T1*(XN1(I)-1))/RH1 TAM=(2.0+XNSUM)/RSUM-(ALFA+T3)*R12-ALFA*XN12 ! *** COMPUTE RIGHT HAND SIDE D DP(I)=TAL*PORTY(I-1)+TAM*PORTY(I)+TAU*PORTY(I+1) 350 CONTINUE RH1=RHO1(N) RH2=RHO2(N) R12=1.0/RH1+1.0/RH2 XN12=T1*((XN1(N)-1)/RH1+(XN2(N)-1)/RH2) RSUM=RH1+RH2 XNSUM=T1*(XN1(N)+XN2(N)-2.0) TAU=(T3+ALFA+ALFA*T1*(XN2(N)-1))/RH2 TAL=(T3+ALFA+ALFA*T1*(XN1(N)-1))/RH1 TAM=(2.0+XNSUM)/RSUM-(ALFA+T3)*R12-ALFA*XN12 ! *** COMPUTE RIGHT HAND SIDE D DP(N)=TAL*PORTY(N-1)+TAM*PORTY(N)+TAU*BOTY(1)- & (T4*T2+ALFA+ALFA*T4*(XN1(N)-1))/RH2*BOTX(1) RETURN END ! ********** ! * PRINTP * ! ********** SUBROUTINE PRINTP ! ***************************************************************** ! *** PRINT PROBLEM PARAMETERS ! ***************************************************************** !ljh INCLUDE 'FOR3D.CMN' USE commons IMPLICIT REAL (KIND=8) (A-H,O-Z) IMPLICIT INTEGER (KIND=4) (I-N) CHARACTER*80 TITLE ! cljh !ljh WRITE(NPU,70) TITLE !ljh 70 FORMAT(//1X,80A1) !ljh 70 FORMAT(//1X,A80) WRITE(NPU,72) 72 FORMAT(/1X,'NUSC FOR3D MODEL') IF(NDIM.EQ.1)THEN WRITE(NPU,*)'2D SOLUTION' ELSE IF(NDIM.EQ.2)THEN WRITE(NPU,*)'N X 2D SOLUTION' ELSE WRITE(NPU,*)'3D SOLUTION' ENDIF ENDIF IF(ISF.EQ.0) WRITE(NPU,73) 73 FORMAT(1X,'GAUSSIAN STARTING FIELD') IF(ISF.EQ.1) WRITE(NPU,74) 74 FORMAT(1X,'USER STARTING FIELD') IF(ISF.EQ.2) WRITE(NPU,71) 71 FORMAT(1X,'GREENS STARTING FIELD') IF(IHNK.NE.0) WRITE(NPU,75) 75 FORMAT(1X,'STARTING FIELD DIVIDED BY HANKEL FUNCTION') IF(ITYPES.EQ.1) WRITE(NPU,76) 76 FORMAT(1X,'USER SURFACE CONDITION') IF(ITYPEB.EQ.1) WRITE(NPU,77) 77 FORMAT(1X,'USER BOTTOM CONDITION') IF(ITYPPW.EQ.1) WRITE(NPU,78) IF(NDIM.EQ.3.AND.ITYPPW.EQ.2) & WRITE(NPU,*)'2D SOLUTION AT PORT SIDEWALL BOUNDARY' 78 FORMAT(1X,'USER PORT SIDEWALL BOUNDARY CONDITION') IF(ITYPSW.EQ.1) WRITE(NPU,79) IF(NDIM.EQ.3.AND.ITYPSW.EQ.2) & WRITE(NPU,*)'2D SOLUTION AT STBD SIDEWALL BOUNDARY' 79 FORMAT(1X,'USER STARBOARD SIDEWALL BOUNDARY CONDITION') IF(NDIM.NE.1)ANG=FLDW/NSEC WRITE(NPU,80) FRQ,ZS,C0,ISF,RA,ZA,N,IHNK,ITYPES,ITYPEB, & ITYPPW,ITYPSW,FLDW,NSEC,ANG,NSOL,NTOT 80 FORMAT(1X,'FRQ = ',F8.2,' HZ',/,1X, & 'ZS = ',F8.2,' M',/,1X, & 'C0 = ',F8.2,' M/SEC',/,1X, & 'ISF = ',I5,/,1X, & 'RA = ',F8.2,' M',/,1X, & 'ZA = ',F8.2,' M',/,1X, & 'N = ',I5,/,1X, & 'IHNK = ',I5,/,1X, & 'ITYPES = ',I5,/,1X, & 'ITYPEB = ',I5,/,1X, & 'ITYPPW = ',I5,/,1X, & 'ITYPSW = ',I5,/,1X, & 'FLDW = ',F8.3,' (WIDTH OF FIELD IN DEGREES)',/,1X, & 'NSEC = ',I5,' SECTORS',/,1X, & 'PHI = ',F8.3,' DEGREES/SECTOR',/,1X, & 'NSOL = ',I5,' SOLUTIONS',/,1X, & 'NTOT = ',I5) WRITE(NPU,81) DR,DZ,ANG,WDR,WZ1,WZ2,WDZ,WDTH,PDR,PDZ,PDTH,RMAX, & DOUGRA,NDIV,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10,U11,U12 81 FORMAT(1X,'DR = ',F8.2,' M',/,1X, & 'DZ = ',F8.2,' M',/,1X, & 'DTH = ',F8.2,' DEGREES',/,1X, & 'WDR = ',F8.2,' M',/,1X, & 'WZ1 = ',F8.2,' M',/,1X, & 'WZ2 = ',F8.2,' M',/,1X, & 'WDZ = ',F8.2,' M',/,1X, & 'WDTH = ',F8.2,' DEGREES',/,1X, & 'PDR = ',F8.2,' M',/,1X, & 'PDZ = ',F8.2,' M',/,1X, & 'PDTH = ',F8.2,' DEGREES',/,1X, & 'RMAX =',F9.1,' M',/,1X, & 'DOUGRA =',F9.1,' M',/,1X, & 'NDIV = ',I5,/,1X, & 'USER 1 = ',F12.4,/,1X, & 'USER 2 = ',F12.4,/,1X, & 'USER 3 = ',F12.4,/,1X, & 'USER 4 = ',F12.4,/,1X, & 'USER 5 = ',F12.4,/,1X, & 'USER 6 = ',F12.4,/,1X, & 'USER 7 = ',F12.4,/,1X, & 'USER 8 = ',F12.4,/,1X, & 'USER 9 = ',F12.4,/,1X, & 'USER 10= ',F12.4,/,1X, & 'USER 11= ',F12.4,/,1X, & 'USER 12= ',F12.4,/,/) 84 FORMAT(3X,'LAYER ',I3,3X,'DEPTH = ',F8.1,3X, & 'DENSITY = ',F6.2,3X,'GRADIENT = ',F7.4,/, & 3X,'BETA = ',F8.3,3X,'GRADIENT = ',F7.4) ! ! *** PRINT BOTTOM DEPTHS IF REQUESTED IF(IBOT.EQ.0) GO TO 86 DO 85 L=1,NSOL TH=THETA(L)*DEG IF(NDIM.NE.1)ANG=-FLDW/2.0+((L-1)*FLDW/NSEC) ! WRITE(NPU,*)' ' WRITE(NPU,*)'RELATIVE BEARING = ',ANG,' DEGREES' WRITE(NPU,123) R1(L),Z1(L),R2(L),Z2(L),TH 85 CONTINUE 86 IF(ISFLD.EQ.0) GO TO 92 ! ! *** PRINT STARTING FIELD WRITE(NPU,87) 87 FORMAT(/,1X,'STARTING FIELD') DO 90 J=1,NSOL IF(NDIM.NE.1)ANG=-FLDW/2.0+((J-1)*FLDW/NSEC) ! WRITE(NPU,*)' ' WRITE(NPU,*)'RELATIVE BEARING = ',ANG,' DEGREES' ! WRITE(NPU,*)' ' M=(J-1)*N DO 90 I=1,N ZI=I*DZ IF(U(M+I).NE.0.0)WRITE(NPU,89) I,ZI,U(M+I) 89 FORMAT(1X,I4,3X,F10.2,3X,'(',E12.5,2X,E12.5,' )') 90 CONTINUE 92 CONTINUE ! IF(ISVP.EQ.0) GO TO 110 ! *** PRINT SVP WRITE(NPU,93) RA 93 FORMAT(//1X,'SOUND VELOCITY PROFILE AT RANGE ',F10.1,' METERS',/) DO 102 L=1,NSOL NLYR=NLYRS(L) IF(NDIM.NE.1)ANG=-FLDW/2.0+((L-1)*FLDW/NSEC) ! WRITE(NPU,*)' ' WRITE(NPU,*)'RELATIVE BEARING = ',ANG,' DEGREES' ! WRITE(NPU,*)' ' K=0 DO 101 ILYR=1,NLYR WRITE(NPU,84)ILYR,ZLYR(ILYR,L),RHO(ILYR,L),RHOG(ILYR,L), & BETA(ILYR,L),BETAG(ILYR,L) ! WRITE(NPU,*)' ' J=K+1 K=IXSVP(ILYR,L) DO 100 I=J,K WRITE(NPU,95) I,ZSVP(I,L),CSVP(I,L) 95 FORMAT(1X,I4,3X,F8.1,3X,F8.1) 100 CONTINUE ! WRITE(NPU,*)' ' 101 CONTINUE 102 CONTINUE 110 CONTINUE RETURN 123 FORMAT(1X,'BOTTOM DEPTHS ',/,1X, & 'R1 = ',F10.1,' M',/,1X, & 'Z1 = ',F8.1,' M',/,1X, & 'R2 = ',F10.1,' M',/,1X, & 'Z2 = ',F8.1,' M',/,1X, & 'THETA = ',F8.1,' DEG') END ! ********** ! * RHS * ! ********** SUBROUTINE RHS ! *** COMPUTE RIGHT HAND SIDE. ! *** RA IS NEXT SOLUTION RANGE. ! *** CAN SAVE TIME BY COMPUTING 2D AND 3D SOLUTIONS IN SEPARATE ! *** SUBROUTINES i.e. RHS2 AND RHS3 ETC. !ljh INCLUDE 'FOR3D.CMN' USE commons IMPLICIT REAL (KIND=8) (A-H,O-Z) IMPLICIT INTEGER (KIND=4) (I-N) COMPLEX T1,T3,T5,VM,VMM1,VMP1,TBL,TBM,TBU,TBL2,TBM2,TBU2 COMPLEX TAL,TAM,TAU,XN12,XNSUM ! *** GET SURFACE CONDITION CALL SCON3D ! *** GET BOTTOM CONDITION CALL BCON3D ! *** GET PORT WALL CONDITION IF(NDIM.EQ.3) CALL PORT3D ! *** GET STARBOARD WALL CONDITION IF(NDIM.EQ.3)CALL STBD3D T1=CMPLX(.25,.25*XK0*DR) T2=1.0/(XK0*XK0*DZ*DZ) T3=T1*T2 IF(NDIM.EQ.1.OR.NDIM.EQ.2.OR.RA.EQ.0.0)THEN T5=0.0 ELSE T5=CMPLX(0.0,XK0*DR/(4.0*XK0*XK0*(RA-.5*DR)*(RA-.5*DR)*PHI*PHI)) ENDIF DO 300 L=1,NSOL LP1=L+1 LP2=L+2 LM1TN=(L-1)*N LTN=L*N LM2TN=(L-2)*N J=(L-2)*N IF(J.LT.0)J=0 K=(L-1)*N DO 200 I=1,N RH3=RHO1(J+I) RH4=RHO1(K+I) R34=(1.0/RH3+1.0/RH4) TBU=.5*(rh3+rh4)*(T5)/RH4 TBM=1.0+(-T5)*.5*(rh3+rh4)*R34 TBL=.5*(rh3+rh4)*(T5)/RH3 IF(I.LT.N)THEN RH5=RHO1(J+I+1) RH6=RHO1(K+I+1) R56=(1.0/RH5+1.0/RH6) TBU2=.5*(rh5+rh6)*(T5)/RH6 TBM2=1.0+(-T5)*.5*(rh5+rh6)*R56 TBL2=.5*(rh5+rh6)*(T5)/RH5 ELSE ENDIF IF(NDIM.EQ.1.OR.NDIM.EQ.2) THEN TBM=1.0 TBM2=1.0 TBL=0.0 TBU=0.0 TBL2=0.0 TBU2=0.0 ELSE ENDIF IM1=I-1 RH1=RHO1(LM1TN+I) RH2=RHO2(LM1TN+I) R12=1.0/RH1+1.0/RH2 IF(L.GT.1) GO TO 70 IF(I.GT.1) GO TO 50 ! ! *** I=1; L=1 VMM1=TBL*SURY(L)+TBM*SURY(LP1)+TBU*SURY(LP2) VM =TBL*PORTY(I)+TBM*U(LM1TN+I)+TBU*U(LTN+I) VMP1=TBL2*PORTY(I+1)+TBM2*U(LM1TN+I+1)+TBU2*U(LTN+I+1) GO TO 130 50 CONTINUE IF(I.EQ.N) GO TO 60 ! ! *** I=2 TO N-1; L=1 VMM1=VM VM =VMP1 VMP1=TBL2*PORTY(I+1)+TBM2*U(LM1TN+I+1)+TBU2*U(LTN+I+1) GO TO 130 60 CONTINUE ! ! *** I=N; L=1 VMM1=VM VM =VMP1 VMP1=TBL*BOTY(L)+TBM*BOTY(LP1)+TBU*BOTY(LP2) GO TO 130 70 CONTINUE IF(L.EQ.NSOL) GO TO 100 IF(I.GT.1) GO TO 80 ! ! *** I=1; L=2 TO NSOL-1 VMM1=TBL*SURY(L)+TBM*SURY(LP1)+TBU*SURY(LP2) VM =TBL*U(LM2TN+I)+TBM*U(LM1TN+I)+TBU*U(LTN+I) VMP1=TBL2*U(LM2TN+I+1)+TBM2*U(LM1TN+I+1)+TBU2*U(LTN+I+1) GO TO 130 80 CONTINUE IF(I.EQ.N) GO TO 90 ! ! *** I=2 TO N-1; L=2 TO NSOL-1 VMM1=VM VM =VMP1 VMP1=TBL2*U(LM2TN+I+1)+TBM2*U(LM1TN+I+1)+TBU2*U(LTN+I+1) GO TO 130 90 CONTINUE ! ! *** I=N; L=2 TO NSOL-1 VMM1=VM VM =VMP1 VMP1=TBL*BOTY(L)+TBM*BOTY(LP1)+TBU*BOTY(LP2) GO TO 130 100 CONTINUE IF(I.GT.1) GO TO 110 ! ! *** I=1; L=NSOL VMM1=TBL*SURY(L)+TBM*SURY(LP1)+TBU*SURY(LP2) VM =TBL*U(LM2TN+I)+TBM*U(LM1TN+I)+TBU*STBDY(I) VMP1=TBL2*U(LM2TN+I+1)+TBM2*U(LM1TN+I+1)+TBU2*STBDY(I+1) GO TO 130 110 CONTINUE IF(I.EQ.N) GO TO 120 ! ! *** I=2 TO N-1; L=NSOL VMM1=VM VM =VMP1 VMP1=TBL2*U(LM2TN+I+1)+TBM2*U(LM1TN+I+1)+TBU2*STBDY(I+1) GO TO 130 120 CONTINUE ! ! *** I=N; L=NSOL VMM1=VM VM =VMP1 VMP1=TBL*BOTY(L)+TBM*BOTY(LP1)+TBU*BOTY(LP2) 130 CONTINUE XN12=T1*((XN1(LM1TN+I)-1)/RH1+(XN2(LM1TN+I)-1)/RH2) RSUM=RH1+RH2 XNSUM=T1*(XN1(LM1TN+I)+XN2(LM1TN+I)-2.0) TAU=(T3+ALFA+ALFA*T1*(XN2(LM1TN+I)-1))/RH2 TAL=(T3+ALFA+ALFA*T1*(XN1(LM1TN+I)-1))/RH1 TAM=(2.0+XNSUM)/RSUM-(ALFA+T3)*R12-ALFA*XN12 ! *** COMPUTE RIGHT HAND SIDE D D(LM1TN+I)=TAL*VMM1+TAM*VM+TAU*VMP1 200 CONTINUE 300 CONTINUE RETURN END ! ********** ! * SCON3D * ! ********** SUBROUTINE SCON3D ! ***************************************************************** ! *** SURFACE BOUNDARY CONDITION SUBROUTINE ! ***************************************************************** !ljh INCLUDE 'FOR3D.CMN' USE commons IMPLICIT REAL (KIND=8) (A-H,O-Z) IMPLICIT INTEGER (KIND=4) (I-N) GO TO (100,200,300),ITYPES+1 WRITE(NPU,50)ITYPES 50 FORMAT(1X,'ERROR - SURFACE BOUNDARY CONDITION ',I2, & ' NOT IMPLEMENTED') IERR=1 RETURN ! 100 CONTINUE ! *** ITYPES=0 ! *** PRESSURE RELEASE SURFACE DO 150 L=1,NSOL+2 SURY(L)=0.0 SURX(L)=0.0 150 CONTINUE RETURN ! 200 CONTINUE ! *** ITYPES=1 ! *** USER SURFACE CONDITION CALL USCON3D RETURN ! 300 CONTINUE ! *** ITYPES = 2 WRITE(NPU,50)ITYPES IERR=1 RETURN END ! ********** ! * SFLD3D * ! ********** SUBROUTINE SFLD3D ! *************************************************************** ! *** GENERATE 3D STARTING FIELD ! *************************************************************** !ljh INCLUDE 'FOR3D.CMN' USE commons IMPLICIT REAL (KIND=8) (A-H,O-Z) IMPLICIT INTEGER (KIND=4) (I-N) COMPLEX HNK,HNKL ! cljh GO TO (100,200,300,400),ISF+1 WRITE(NPU,50)ISF 50 FORMAT(1X,'ERROR - STARTING FIELD ISF ',I2,' NOT IMPLEMENTED') IERR=1 RETURN ! 100 CONTINUE ! *** ISF=0 ! *** GAUSSIAN STARTING FIELD. GW=C0/(PI*FRQ) GA=1.0/SQRT(GW) DO 120 I=1,N ZM=I*DZ PR=GA*(EXP(-((ZM-ZS)/GW)**2)-EXP(-((-ZM-ZS)/GW)**2)) U(I)=CMPLX(PR,0.0) IF(NDIM.EQ.1)GO TO 120 DO 110 L=2,NSOL M=(L-1)*N+I U(M)=U(I) 110 CONTINUE 120 CONTINUE IF(NDIM.EQ.1)GO TO 900 IF(ITYPPW.NE.2) GO TO 135 DO 130 I=1,N 130 PORTY(I)=U(I) 135 IF(ITYPSW.NE.2) GO TO 900 DO 140 I=1,N 140 STBDY(I)=U(I) GO TO 900 ! 200 CONTINUE ! *** ISF=1 ! *** USER STARTING FIELD CALL USFLD3D GO TO 900 ! 300 CONTINUE ! *** ISF=2 ! *** GREENS WIDE ANGLE STARTING FIELD XK0=2.0*PI*FRQ/C0 DO 320 I=1,N ZM=I*DZ PR=SQRT(XK0)*(1.45-0.42*XK0*XK0*(ZM-ZS)**2.0)* & EXP(-(XK0*XK0*(ZM-ZS)**2)/3.0512) U(I)=CMPLX(PR,0.0) IF(NDIM.EQ.1)GO TO 320 DO 310 L=2,NSOL M=(L-1)*N+I U(M)=U(I) 310 CONTINUE 320 CONTINUE IF(NDIM.EQ.1)GO TO 900 IF(ITYPPW.NE.2) GO TO 335 DO 330 I=1,N 330 PORTY(I)=U(I) 335 IF(ITYPSW.NE.2) GO TO 900 DO 340 I=1,N 340 STBDY(I)=U(I) GO TO 900 ! 400 CONTINUE ! *** ISF=3 WRITE(NPU,50) ISF IERR=1 RETURN ! 900 CONTINUE ! *** DIVIDE STARTING FIELD BY HANKEL FUNCTION IF REQUESTED BY USER IF(IHNK.EQ.0) GO TO 950 HNK=HNKL(XK0*RA) DO 930 J=1,NSOL M=(J-1)*N DO 930 I=1,N U(M+I)=U(M+I)/HNK 930 CONTINUE DO 935 I=1,N 935 PORTY(I)=PORTY(I)/HNK DO 945 I=1,N 945 STBDY(I)=STBDY(I)/HNK 950 CONTINUE RETURN END ! ********** ! * STBD2D * ! ********** SUBROUTINE STBD2D ! ****************************************************************** ! * SOLVE FOR FIELD AT STBD WALL BOUNDARY ! * FIELD AT ADVANCED RANGE IN STBDX ! ****************************************************************** !ljh INCLUDE 'FOR3D.CMN' USE commons IMPLICIT REAL (KIND=8) (A-H,O-Z) IMPLICIT INTEGER (KIND=4) (I-N) ! *** COMPUTE BTA(N1) THROUGH BTA(N) ! *** ARRAY BTA CONTAINS PARTIAL SOLUTION OF SYSTEM OF EQUATIONS J=(NSOL-1)*N BTA(1)=AM(J+1) M=2 DO 20 I=M,N BTA(I)=AM(J+I)-AL(J+I)*AU(J+I-1)/BTA(I-1) 20 CONTINUE DO 21 I=1,N BTA(I)=1.0/BTA(I) 21 CONTINUE !ljh TA=0.0 !ljh TB=0.0 ! *** STEP 1 - CALL TRIDIAGONAL SOLVER !ljh CALL TRID3D(AL(J+1),AM(J+1),AU(J+1),STBDX(1),DS(1),BTA(1),N,N, !ljh +TA,TB) CALL TRID3D(AL(J+1),AU(J+1),STBDX(1),DS(1),BTA(1),N,N) ! *** SOLUTION IS IN ARRAY STBDX RETURN END ! ********** ! * STBD3D * ! ********** SUBROUTINE STBD3D ! ****************************************************************** ! * STARBOARD SIDEWALL BOUNDARY CONDITION ! ****************************************************************** !ljh INCLUDE 'FOR3D.CMN' USE commons IMPLICIT REAL (KIND=8) (A-H,O-Z) IMPLICIT INTEGER (KIND=4) (I-N) COMPLEX T1,T3,T4,TAL,TAU,TAM,XN12,XNSUM ! GO TO (100,200,300),ITYPSW+1 IERR=1 WRITE(NPU,50) ITYPSW 50 FORMAT(1X,'ERROR - STARBOARD SIDEWALL BOUNDARY CONDITION ',I2, & ' NOT IMPLEMENTED') RETURN ! 100 CONTINUE ! *** ITYPSW = 0 ! *** SET FIELD AT STBD SIDEWALL TO 0.0 DO 150 I=1,N STBDX(I)=0.0 STBDY(I)=0.0 150 CONTINUE RETURN ! 200 CONTINUE ! *** ITYPSW = 1 ; USER SUPPLIES SUBROUTINE USTBD3D CALL USTBD3D RETURN ! 300 CONTINUE ! *** ITYPSW = 2 ! *** ASSUMPTION IS THAT ENVIRONMENTAL CONDITIONS AT STBD SIDEWALL ! *** BOUNDARY ARE EQUAL TO CONDITIONS AT RIGHTMOST SOLUTION BOUNDARY. ! *** SOLUTION AT BOUNDARY AT ADVANCED RANGE IS SOLVED IN 2D. DELTA=XK0*DR T1=CMPLX(.25,.25*DELTA) T2=1.0/(XK0*XK0*DZ*DZ) T3=T1*T2 T4=CMPLX(.25,-.25*DELTA) M=(NSOL-1)*N RH1=RHO1(M+1) RH2=RHO2(M+1) R12=1.0/RH1+1.0/RH2 XN12=T1*((XN1(M+1)-1)/RH1+(XN2(M+1)-1)/RH2) RSUM=RH1+RH2 XNSUM=T1*(XN1(M+1)+XN2(M+1)-2.0) TAU=(T3+ALFA+ALFA*T1*(XN2(M+1)-1))/RH2 TAL=(T3+ALFA+ALFA*T1*(XN1(M+1)-1))/RH1 TAM=(2.0+XNSUM)/RSUM-(ALFA+T3)*R12-ALFA*XN12 ! *** COMPUTE RIGHT HAND SIDE D DS(1)=TAL*SURY(NSOL+2)+TAM*STBDY(1)+TAU*STBDY(2)- & (T4*T2+ALFA+ALFA*T4*(XN1(M+1)-1))/RH1*SURX(NSOL+2) DO 350 I=2,N-1 RH1=RHO1(M+I) RH2=RHO2(M+I) R12=1.0/RH1+1.0/RH2 XN12=T1*((XN1(M+I)-1)/RH1+(XN2(M+I)-1)/RH2) RSUM=RH1+RH2 XNSUM=T1*(XN1(M+I)+XN2(M+I)-2.0) TAU=(T3+ALFA+ALFA*T1*(XN2(M+I)-1))/RH2 TAL=(T3+ALFA+ALFA*T1*(XN1(M+I)-1))/RH1 TAM=(2.0+XNSUM)/RSUM-(ALFA+T3)*R12-ALFA*XN12 ! *** COMPUTE RIGHT HAND SIDE D DS(I)=TAL*STBDY(I-1)+TAM*STBDY(I)+TAU*STBDY(I+1) 350 CONTINUE RH1=RHO1(M+N) RH2=RHO2(M+N) R12=1.0/RH1+1.0/RH2 XN12=T1*((XN1(M+N)-1)/RH1+(XN2(M+N)-1)/RH2) RSUM=RH1+RH2 XNSUM=T1*(XN1(M+N)+XN2(M+N)-2.0) TAU=(T3+ALFA+ALFA*T1*(XN2(M+N)-1))/RH2 TAL=(T3+ALFA+ALFA*T1*(XN1(M+N)-1))/RH1 TAM=(2.0+XNSUM)/RSUM-(ALFA+T3)*R12-ALFA*XN12 ! *** COMPUTE RIGHT HAND SIDE D DS(N)=TAL*STBDY(N-1)+TAM*STBDY(N)+TAU*BOTY(NSOL+2)- & (T4*T2+ALFA+ALFA*T4*(XN1(M+N)-1))/RH2*BOTX(NSOL+2) RETURN END ! ********** ! * SVP3D * ! ********** SUBROUTINE SVP3D ! ****************************************************************** ! *** SOUND VELOCITY PROFILE SUBROUTINE ! ****************************************************************** ! *** CALLING PROGRAM SUPPLIES: NOTHING ! *** SVP SUBROUTINE RETURNS: ! NLYRS - NUMBER OF LAYERS. LAYER 1 IS WATER. OTHERS ARE SEDIMENT ! ZLYR - ARRAY - DEPTH OF EACH LAYER. FIRST IS DEPTH OF WATER. ! RHO - ARRAY - DENSITY OF EACH LAYER. GRAMS/CUBIC CM ! RHOG - ARRAY - DENSITY GRADIENT. ! BETA - ARRAY - ATTENUATION IN EACH LAYER. DB/WAVELENGTH ! BETAG- ARRAY - ATTENUATION GRADIENT. ! IXSVP- ARRAY - CONTAINS POINTERS. POINTS TO LAST VALUE OF SVP ! IN CORRESPONDING LAYER. SVP IS STORED IN ARRAYS ZSVP ! AND CSVP. IXSVP(1,L) POINTS TO LAST SVP POINT IN WATER ! ALONG SECTOR BOUNDARY L. L=1,NSOL. ! IXSVP(NLYR,L) POINTS TO LAST SVP POINT IN BOTTOM-MOST LAYER. ! NSVP - NUMBER OF POINTS IN ZSVP AND CSVP. ZSVP AND CSVP ! CONTAIN THE PROFILES FOR ALL LAYERS. ! ZSVP - ARRAY - SVP DEPTHS - METERS ! CSVP - ARRAY - SOUND SPEED - METERS/SEC ! ****************************************************************** ! !ljh INCLUDE 'FOR3D.CMN' USE commons IMPLICIT REAL (KIND=8) (A-H,O-Z) IMPLICIT INTEGER (KIND=4) (I-N) ! *** READ NUMBER OF LAYERS DO 75 L=1,NSOL NSVP=0 READ(NIU,*,END=100) NLYRS(L) NLYR=NLYRS(L) ! ! *** IF NLYR IS NEGATIVE, REPEAT LAST PROFILE UNTIL NSOL PROFILES ! *** HAVE BEEN ENTERED. IF(NLYR.LT.0) GO TO 200 LNLYR=NLYR LASTL=L ! *** FIRST LAYER IS WATER. OTHERS ARE SEDIMENT. DO 55 I=1,NLYR ! *** READ DEPTH OF LAYER, DENSITY AND ATTENUATION. READ(NIU,*,END=100) ZLYR(I,L),RHO(I,L),RHOG(I,L),BETA(I,L), & BETAG(I,L) ! *** READ PROFILE. 50 READ(NIU,*,END=100) ZV,CV NSVP=NSVP+1 ZSVP(NSVP,L)=ZV CSVP(NSVP,L)=CV IF(ZV.LT.ZLYR(I,L)) GO TO 50 IF(ZV.GT.ZLYR(I,L)) GO TO 100 IXSVP(I,L)=NSVP 55 CONTINUE 75 CONTINUE RETURN ! ! *** ERROR EXIT 100 NSVP=0 RETURN ! ! *** REPEAT LAST PROFILE UNTIL NSOL PROFILES ENTERED. 200 CONTINUE NLYR=LNLYR NSVP=IXSVP(NLYR,LASTL) DO 230 L=LASTL+1,NSOL NLYRS(L)=NLYR DO 210 I=1,NLYR ZLYR(I,L)=ZLYR(I,LASTL) RHO(I,L)=RHO(I,LASTL) RHOG(I,L)=RHOG(I,LASTL) BETA(I,L)=BETA(I,LASTL) BETAG(I,L)=BETAG(I,LASTL) IXSVP(I,L)=IXSVP(I,LASTL) 210 CONTINUE DO 220 J=1,NSVP ZSVP(J,L)=ZSVP(J,LASTL) CSVP(J,L)=CSVP(J,LASTL) 220 CONTINUE 230 CONTINUE RETURN END ! ********** ! * TRID3D * ! ********** SUBROUTINE TRID3D(XL,XU,U,D,BTA,L,M) !ljh SUBROUTINE TRID3D(XL,XM,XU,U,D,BTA,L,M,TA,TB) ! ****************************************************************** ! *** SPECIALIZED VERSION OF METHOD PRESENTED ON PG 442 OF CARNAHAN ! *** ET AL FOR SOLVING A TRIDIAGONAL MATRIX. ! *** TRID RETURNS THE SOLUTION FIELD IN U ! *** BTA IS PARTIAL SOLUTION COMPUTED ELSEWHERE ! *** D CONTAINS R.H.S. ! *** GAMMA - TEMPORARY STORAGE ! *** TA - TEMPORARY VARIABLE USED IN MANIPULATION OF BOTTOM POINT ! *** TB - TEMPORARY VARIABLE USED IN MANIPULATION OF BOTTOM POINT ! ****************************************************************** ! !ljh INCLUDE 'FOR3D.PAR' USE params IMPLICIT REAL (KIND=8) (A-H,O-Z) IMPLICIT INTEGER (KIND=4) (I-N) !ljh COMPLEX XL(MXN),XM(MXN),XU(MXN),U(MXN),GAMMA(MXN),D(MXN), !ljh + BTA(MXN),TA,TB COMPLEX XL(MXN),XU(MXN),U(MXN),GAMMA(MXN),D(MXN),BTA(MXN) ! *** TA AND TB MAY BE USED AT A LATER DATE. GAMMA(1)=D(1)*BTA(1) DO 10 I=2,L GAMMA(I)=(D(I)-XL(I)*GAMMA(I-1))*BTA(I) 10 CONTINUE ! *** IF L IS LESS THAN M, U(M) IS KNOWN. IF(L.EQ.M) U(M)=GAMMA(M) DO 70 I=M,2,-1 U(I-1)=GAMMA(I-1)-XU(I-1)*U(I)*BTA(I-1) 70 CONTINUE RETURN END ! ********** ! * TWOSTEP* ! ********** SUBROUTINE TWOSTEP ! ****************************************************************** ! * SOLVE V ! * SOLVE U ! * FIELD AT ADVANCED RANGE IN U ! ****************************************************************** !ljh INCLUDE 'FOR3D.CMN' USE commons IMPLICIT REAL (KIND=8) (A-H,O-Z) IMPLICIT INTEGER (KIND=4) (I-N) COMPLEX T1,T2,T3,T4,T5,TBL,TBM,TBU,TAL,TAU !ljh COMPLEX T1,T2,T3,T4,T5,VM,VMM1,VMP1,TBL,TBM,TBU,TAL,TAU !ljh COMPLEX VBOTX,VBOTY,VSURX,VSURY !ljh COMPLEX VBOTX,VSURX,VSURY ! Got complaints for those not used. COMPLEX VBOTX,VSURX IF(NDIM.EQ.3.AND.ITYPPW.EQ.2) CALL PORT2D IF(NDIM.EQ.3.AND.ITYPSW.EQ.2) CALL STBD2D T1=CMPLX(.25,.25*XK0*DR) T2=CMPLX(.25,-.25*XK0*DR) T4=T1/(XK0*XK0*DZ*DZ) IF(NDIM.EQ.1.OR.NDIM.EQ.2.OR.RA.EQ.0.0)THEN T5=0.0 ELSE T5=-CMPLX(0.0,XK0*DR/(4.0*XK0*XK0*(RA-.5*DR)*(RA-.5*DR)*PHI*PHI)) ENDIF DO 10 L=1,NSOL M=(L-1)*N J=(L-2)*N IF(J.LT.0)J=0 RH3=RHO1(J+1) K=(L-1)*N RH4=RHO1(K+1) R34=(1.0/RH3+1.0/RH4) TBU=.5*(rh3+rh4)*(T5)/RH4 TBM=1.0+(-T5)*.5*(rh3+rh4)*R34 TBL=.5*(rh3+rh4)*(T5)/RH3 IF(NDIM.EQ.1.OR.NDIM.EQ.2)THEN TBL=0.0 TBM=1.0 TBU=0.0 ELSE ENDIF VSURX=TBL*SURX(L)+TBM*SURX(L+1)+TBU*SURX(L+2) T3=CMPLX(.25,-.25*XK0*DR)/(XK0*XK0*DZ*DZ) TAL=(T3+ALFA+ALFA*T2*(XN1(M+1)-1))/RHO1(M+1) D(M+1)=D(M+1)-TAL*VSURX RH3=RHO1(J+N) RH4=RHO1(K+N) R34=(1.0/RH3+1.0/RH4) TBU=.5*(rh3+rh4)*(T5)/RH4 TBM=1.0+(-T5)*.5*(rh3+rh4)*R34 TBL=.5*(rh3+rh4)*(T5)/RH3 IF(NDIM.EQ.1.OR.NDIM.EQ.2)THEN TBL=0.0 TBM=1.0 TBU=0.0 ELSE ENDIF VBOTX=TBL*BOTX(L)+TBM*BOTX(L+1)+TBU*BOTX(L+2) TAU=(T3+ALFA+ALFA*T2*(XN2(M+N)-1))/RHO2(M+N) D(M+N)=D(M+N)-TAU*VBOTX 10 CONTINUE ! *** COMPUTE BTA(N1) THROUGH BTA(N) ! *** ARRAY BTA CONTAINS PARTIAL SOLUTION OF SYSTEM OF EQUATIONS BTA(1)=AM(1) M=2 DO 20 I=M,NTOT BTA(I)=AM(I)-AL(I)*AU(I-1)/BTA(I-1) 20 CONTINUE DO 21 I=1,NTOT BTA(I)=1.0/BTA(I) 21 CONTINUE !ljh TA=0.0 !ljh TB=0.0 ! *** STEP 1 - CALL TRIDIAGONAL SOLVER !ljh CALL TRID3D(AL,AM,AU,U,D,BTA,NTOT,NTOT,TA,TB) CALL TRID3D(AL,AU,U,D,BTA,NTOT,NTOT) ! *** SOLUTION IS IN ARRAY U IF(NDIM.EQ.1.OR.NDIM.EQ.2) RETURN J=0 DO 50 I=1,N DO 50 L=1,NSOL J=J+1 M=(L-1)*N+I D(J)=U(M) 50 CONTINUE J=(NSOL-1)*N DO 60 I=1,N M=(I-1)*NSOL RH3=RHO1(I) RH4=RHO1(J+I) ! TBL=.5*(rh3+rh3)*(T5)/RH3 TBL=(T5) D(M+1)=D(M+1)-PORTX(I)*TBL ! TBU=.5*(rh4+rh4)*(T5)/RH4 TBU=(T5) D(M+NSOL)=D(M+NSOL)-STBDX(I)*TBU 60 CONTINUE BTA(1)=BM(1) M=2 DO 110 I=M,NTOT BTA(I)=BM(I)-BL(I)*BU(I-1)/BTA(I-1) 110 CONTINUE DO 111 I=1,NTOT BTA(I)=1.0/BTA(I) 111 CONTINUE !ljh TA=0.0 !ljh TB=0.0 ! *** STEP 2 - CALL TRIDIAGONAL SOLVER !ljh CALL TRID3D(BL,BM,BU,V,D,BTA,NTOT,NTOT,TA,TB) CALL TRID3D(BL,BU,V,D,BTA,NTOT,NTOT) ! *** SOLUTION IN V BUT MUST BE REORDERED. ! *** REORDER SOLUTION J=0 DO 150 I=1,N IF(ITYPPW.EQ.2)PORTY(I)=PORTX(I) IF(ITYPSW.EQ.2)STBDY(I)=STBDX(I) DO 150 L=1,NSOL J=J+1 M=(L-1)*N+I U(M)=V(J) 150 CONTINUE ! *** SOLUTION FIELD AT ADVANCED RANGE IS IN ARRAY U. RETURN END