#! /bin/sh # This is a shell archive. Remove anything before this line, then unpack # it by saving it into a file and typing "sh file". To overwrite existing # files, type "sh file -c". You can also feed this as standard input via # unshar, or by typing "sh 'dno/README' <<'END_OF_FILE' Xc************************* README **************************************c Xc c Xc Copyright (C) 1989, 1990 William W. Symes c Xc c Xc This file is part of the Down-n-Out Traveltime Computation c Xc Project ("DNO"). c Xc c Xc The DNO codes are distributed in the hope that they will be c Xc useful, but WITHOUT ANY WARRANTY, without even the implied c Xc warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. c Xc c Xc You are free to modify this file provided that you leave this c Xc header intact. I ask that you not redistribute this file c Xc without first notifying me. This isn't because I'm trying c Xc to be nasty here; it's simply that I would like to keep track c Xc of who has the program. c Xc c Xc These codes were developed by William W. Symes with support from c Xc the National Science Foundation, the Office of Naval Research, c Xc and the Air Force Office of Scientific Research. Thanks are due c Xc to Jim Carazzone, Dave Hale, and Jos van Trier for helpful c Xc suggestions. c Xc c Xc William W. Symes c Xc Department of Mathematical Sciences c Xc Rice University c Xc Houston, Texas 77005 c Xc USA c Xc c Xc e-mail: symes@rice.edu c Xc c Xc March 1990 c Xc c Xc***********************************************************************c Xc Xc For some introductory reading matter, see the paper Xc "Upwind finite difference calculation of traveltimes", Xc by J. van Trier and W. W. Symes, available from the Xc authors (Mar. '90); also to appear in GEOPHYSICS sometime Xc in the next year or so. Xc Xc================================================================ Xc Xc PACKAGING: Xc Xc You should have: Xc dno.f...........main program module Xc slow.f..........generates example slowness files Xc init.f..........sets up initial time patches Xc datafile.dat....a collection of sample default and job Xc files - MUST be picked apart and Xc placed in the active directory. Xc Xc The ".f" files are unitary source files, which can be Xc compiled on their own. UNIX users will want to "fsplit" Xc them and set up an appropriate makefile. Xc Xc================================================================ Xc Xc COMPATIBILITY Xc Xc This code has been used recently with success on: Xc Xc 1. a 386 clone with Weitek fpu, DOS, compiled with NDP Xc 386Fortran; Xc 2. Sun3 (68020) and Sun4/Sparc machines, Sun UNIX, f77; Xc 3. Ardent Titan, their version of System V UNIX, X11, Xc fc @ O3; Xc 4. Apollo 4500, Domain OS, f77; Xc 5. Cray 1S, COS, CFT. Xc Xc Good luck. Xc Xc================================================================ Xc Xc USING IT: Xc Xc INPUTS: Xc Xc Defaults - consisting of unit numbers, various parameters Xc unlikely to be changed, and the vector/scalar loop option, Xc flag - are read by SETDEF from the default Xc file "dnodef.dat", which must be present (a sample is Xc included in the "datafile.dat"). Xc Xc Two options for basic input: Xc Xc 1. Interactive: Xc Inputs are given on "stdin" (unit 5), and various levels Xc of output on "stdout" (unit 6) (set in "dnodef.dat"); Xc Xc 2. batch: Xc Inputs are read form the file "dno.job" (sample included), Xc signaled by input unit >< 5 AND output unit >< 6 in Xc "dnodef.dat". Xc Xc Job Name - character*20 Xc will be the name of the output files Xc containing the computed travel times - Xc i.e. [jobname].t etc. (see below for Xc a fuller description of output files). Xc Xc dstep, lstep, rstep - integer Xc the calculation proceeds in a cycle, Xc first dstep steps down, the lstep steps Xc to the left, then rstep steps to the Xc right. thus these numbers determine Xc the rate of expansion and shape of the Xc computational front. Xc Xc chkshot - integer Xc the number of constant-offset ("checkshot") Xc profiles, which give Xc the travel-time on vertical lines ("wells"); Xc none if set = 0 Xc Xc xchk(*) - real array Xc offset for each "well" Xc Xc horizons - integer Xc the number of constant-depth profiles Xc Xc zhor(*) - real array Xc depth for each horizon Xc Xc OUTPUT: Xc Xc on "stdout" (interactive option) or output unit as specified Xc in "dnodef.dat" (batch option): a record of the calculation, at varying Xc levels of verbosity according to the parameter "idump", Xc set before each subroutine call in MAIN: Xc idump < 0: a dour, silent option; merely Xc remarks the entry into each subroutine, Xc except for input queries and Xc complaints about errors. Xc idump = 0: notes successful completion of each Xc step, amongst other things. Xc idump > 0: also spews vastly detailed records Xc of intermediate computational steps - Xc really useful only for debugging. Xc Xc For the batch option, the external filename for this output is Xc "dno.out". Xc Xc on [jobname].t: an array containing the computed travel Xc times. Xc Xc on [jobname].xn, n=1,chkshot: constant-offset profiles Xc Xc on [jobname].zn, n=1,horizons: constant-depth profiles Xc Xc Xc OTHER NECESSARY FILES: Xc Xc 1. the default file "dnodef.dat" Xc Xc 2. for the batch option, the job file "dno.job" Xc Xc 3. Arrays containing the slowness field and initial time data Xc are read from files. The file names are read in INIT. Xc Xc [slowfile]:array of slowness values. Implicitly Xc defines maximal computational domain. Xc Xc [timefile]:array of times. Initially set on a small Xc subrectangle of the computational domain Xc containing the source location, with at least two rows Xc and columns. For accuracy reasons, a larger Xc rectangle around a point source gives better Xc results than a small one (since then the wavefront Xc curvature is less). Typically the initial time Xc rectangle will be a small subset of the slowness Xc rectangle. Xc Xc Formats set up for these files are as specified in the headers Xc of the i/o routines SREADMAT and SWRITMAT. These are simple "ascii flat Xc file" formats, chosen for compatibility with Matlab (which I use Xc extensively for data preparation and graphics). Xc Xc The programs SLOW and INIT compute examples of Xc these files, with filenames "slow.dat" and "time0.dat". Xc Xc NOTE THAT SLOW AND INIT HAVE THEIR OWN DEFAULT FILES "slowdef.dat" Xc and "initdef.dat" respectively. THESE MUST BE PRESENT - examples Xc are included in "datafile.dat". Xc Xc====================== program dno ============================== Xc Xc GLOSSARY of Variables Xc Xc nxmin,nxmax maximal x-limits for all arrays, Xc set as parameters in MAIN Xc nzmin,nzmax maximal z-limits for all arrays, Xc set as parameters in MAIN Xc ier error flag - if > 0 on return Xc exit immediately Xc idump dump switch - see above Xc ipr output unit number (usually 6) Xc irpt unit number for time output files, Xc set as default in INIT Xc jminus,jplus x-limits of computational rectangle, Xc initialized to those of the initial Xc time data rectangle, subsequently adjusted Xc in SOLVE at the completion of each Xc successful block of steps (down, left, Xc of right). Xc kminus,kplus z-limits of the computational rectangle, Xc initialized to those of the initial Xc time data rectangle, subsequently adjusted Xc in SOLVE at the completion of each Xc successful block of steps (down, left, Xc of right). Xc dzswh fixed/variable step switch, set as default Xc in INIT. Xc dzswh = 1: fixed steps (dangerous!) Xc else: adaptive steps Xc jstop switch for stopping at a corner inflow Xc condition. In order to proceed, the Xc flow must be OUT at each end of the faces Xc of the rectangular compuational front, Xc else an upwind step cannot be taken there. Xc jstop = 0: drop the endpoint in question, Xc continue computing on smaller Xc face ( this option only works Xc in the present implementation Xc for marching in ONE DIRECTION Xc ONLY, e.g. down). Prevents com- Xc putation of full rectangle of Xc times. Xc jstop > 0: stop computation in current Xc direction when inflow is Xc detected. This option is Xc consistent with the design of Xc the algorithm. Xc jcut flux cutoff switch: Xc jcut = 0: step computed adaptively Xc jcut ><0: current flow variable (x or Xc z gradient component) truncated Xc to limit signal velocity to 95% Xc of the maximum permissible. Xc This option Xc ensures that the calculation Xc will go to completion even when Xc it shouldn't. Xc vector character*3 switch to control main finite Xc difference loop coding: Xc vector='no ': scalar register variables used Xc vector='yes': direct array references used Xc The first option should be more efficient in Xc architecture/compiler combinations which make use Xc of scalar registers, whereas the second is NECESSARY Xc to achieve vectorization under Cray CFT. Xc Set in INIT. Xc maxstep max. number if partial steps allowed - set Xc in INIT. Xc jxmin,jxmax x-limits of computational domain, set during Xc reading of slowness array. Xc kzmax max. z-limit of computational domain, set Xc during reading of slowness array. Note that Xc the minimum z-limit remains at its initial Xc value, i.e. kminus - down and out, but not Xc up! Xc maxsweep max. number of "down, left, right" cycles Xc permitted Xc dstep, lstep, rstep Xc down, left, and right steps per cycle (see Xc above). Xc xsam0, zsam0 x and z steps as read from slowness file, and Xc checked against initial time file. Xc cfl Courant-Friedrichs-Lewy factor, used in Xc setting adaptive step. Set as default in Xc INIT. Xc time array for initial and computed time values Xc slow array for input slowness values Xc Xc *** the next set of variables provide storage for Xc the time field and its gradient components on Xc the three parts (bottom, left, and right) of Xc the computational front (no top component - Xc again, "down and out"!). For reasons we won't Xc go into, w is used to denote the x-derivative Xc and u the z-derivative, to which are added "b", Xc "l", and "r" as appropriate. For some arrays we Xc also keep a "**temp" array for updated values, Xc to avoid vectorization-destroying dependencies. Xc Xc wb, wbtemp, ub, tb Xc ul, ultemp, wl, tl Xc ur, urtemp, wr, tr Xc Xc sb, slr temporary storage for interpolated Xc slowness values on bottom and Xc left/right fronts, respectively. Xc Xc jobname character*20 for job name Xc jobsize length of [jobname] Xc maxchk maximum number of constant-offset profiles Xc maxhor maximum number of constant-depth profiles Xc chkshot number of constant-offset profiles requested Xc horizons number of constant-depth profiles requested Xc xchk(*) offsets Xc zhor(*) depths Xc Xc Xc================================================================ Xc X X X END_OF_FILE if test 14029 -ne `wc -c <'dno/README'`; then echo shar: \"'dno/README'\" unpacked with wrong size! fi # end of 'dno/README' fi if test -f 'dno/datafile.dat' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'dno/datafile.dat'\" else echo shar: Extracting \"'dno/datafile.dat'\" \(2056 characters\) sed "s/^X//" >'dno/datafile.dat' <<'END_OF_FILE' X// dnodef.dat - begin ================================================ X// DNODEF.DAT - sets defaults for dno X5 / input unit number X6 / output (report) unit number X3 / output (file) unit number X0.98e+00 / CFL number X0 / adjust (0) or not (1) step (def.: 0) X20 / maximum partial steps per step X1 / stop (1) or not (0) at inflow condition (def.: 1) X0 / truncate slope (1) or don't (0) (def.: 0) Xyes / use "vector" (yes) or "scalar" (no ) main loop X//dnodef.dat - end =================================================== X//dno.job -begin ===================================================== Xtest Xslow.dat Xtime0.dat X10 X10 X10 X1 X.2 X1 X.5 X// dno.job - end ==================================================== X// slowdef.dat - begin ============================================== X5 / input unit number X6 / output (report) unit number X3 / output (file) unit number X0.98e+00 / CFL number X0 / adjust (0) or not (1) step (def.: 0) X20 / maximum partial steps per step X1 / stop (1) or not (0) at inflow condition (def.: 1) X0 / truncate slope (1) or don't (0) (def.: 0) Xno / use "vector" (yes) or "scalar" (no ) main loop X// slowdef.dat - end ================================================ X// initdef.dat - begin ============================================== X5 / input unit number X6 / output (report) unit number X3 / output (file) unit number X0.98e+00 / CFL number X0 / adjust (0) or not (1) step (def.: 0) X20 / maximum partial steps per step X1 / stop (1) or not (0) at inflow condition (def.: 1) X0 / truncate slope (1) or don't (0) (def.: 0) Xno / use "vector" (yes) or "scalar" (no ) main loop X// initdef.dat - end ================================================= X END_OF_FILE if test 2056 -ne `wc -c <'dno/datafile.dat'`; then echo shar: \"'dno/datafile.dat'\" unpacked with wrong size! fi # end of 'dno/datafile.dat' fi if test -f 'dno/dno.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'dno/dno.f'\" else echo shar: Extracting \"'dno/dno.f'\" \(68102 characters\) sed "s/^X//" >'dno/dno.f' <<'END_OF_FILE' Xc***********************************************************************c Xc c Xc Copyright (C) 1989, 1990 William W. Symes c Xc c Xc This file is part of the Down-n-Out Traveltime Computation c Xc Project ("DNO"). c Xc c Xc The DNO codes are distributed in the hope that they will be c Xc useful, but WITHOUT ANY WARRANTY, without even the implied c Xc warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. c Xc c Xc You are free to modify this file provided that you leave this c Xc header intact. I ask that you not redistribute this file c Xc without first notifying me. This isn't because I'm trying c Xc to be nasty here; it's simply that I would like to keep track c Xc of who has the program. c Xc c Xc These codes were developed by William W. Symes with support from c Xc the National Science Foundation, the Office of Naval Research, c Xc and the Air Force Office of Scientific Research. Thanks are due c Xc to Jim Carazzone, Dave Hale, and Jos van Trier for helpful c Xc suggestions. c Xc c Xc William W. Symes c Xc Department of Mathematical Sciences c Xc Rice University c Xc Houston, Texas 77005 c Xc USA c Xc c Xc e-mail: symes@rice.edu c Xc c Xc March 1990 c Xc c Xc***********************************************************************c Xc Xc-------------------------------------------------------------- Xc X program dno Xc Xc parameters Xc X integer*4 nxmin,nxmax,nzmin,nzmax,maxstep X parameter (nxmin=-200,nxmax=200) X parameter (nzmin=0,nzmax=200) X parameter (maxstep=25) Xc X integer*4 ier,idump,ipr,irpt,jminus,jplus, X & kminus,kplus, X & dzswh,jstop,jcut,naxstep, X & kzmax,jxmin,jxmax,maxsweep, X & dstep,lstep,rstep Xc X real*4 xsam0,zsam0,cfl Xc X real*4 time(nxmin:nxmax,nzmin:nzmax), X & slow(nxmin:nxmax,nzmin:nzmax) Xc X real*4 wb(nxmin:nxmax),wbtemp(nxmin:nxmax), X & ul(nzmin:nzmax),ultemp(nzmin:nzmax), X & ur(nzmin:nzmax),urtemp(nzmin:nzmax), X & sb(nxmin:nxmax),slr(nzmin:nzmax), X & tb(nxmin:nxmax), X & tl(nzmin:nzmax), tr(nzmin:nzmax), X & wl(nzmin:nzmax), wr(nzmin:nzmax), X & ub(nxmin:nxmax) Xc X integer jobsize X character*20 jobname X character*3 vector X integer chkshot,horizons X integer maxchk,maxhor X parameter (maxchk=10,maxhor=10) X real xchk(maxchk),zhor(maxhor) Xc Xc initialize work arrays Xc X jobname=' ' X idump=-1 X call init(ier,idump,ipr,irpt,dzswh,jstop,jcut, X & vector, X & nxmin,nxmax,nzmin,nzmax,xsam0,zsam0,cfl, X & naxstep,jminus,jplus,kminus,kplus, X & kzmax,jxmin,jxmax,maxsweep,dstep,lstep,rstep, X & time,slow,jobname,jobsize, X & maxchk,maxhor,chkshot,horizons,xchk,zhor) X if (ier.gt.0) go to 1000 X if(naxstep.gt.maxstep) naxstep=maxstep Xc Xc propagate in depth Xc X idump=-1 X call solve(ier,idump,ipr,dzswh,jstop,jcut, X & vector, X & nxmin,nxmax,nzmin,nzmax,xsam0,zsam0, X & maxstep,jminus,jplus,kminus,kplus, X & kzmax,jxmin,jxmax,maxsweep, X & dstep,lstep,rstep, X & cfl,time,slow,sb,slr,wb,wbtemp,tb, X & ul,ultemp,tl,ur,urtemp,tr,wl,wr,ub) X if (ier.gt.0) go to 1000 Xc Xc write results to files Xc X idump=-1 X call output(ier,idump,ipr,irpt, X & nxmin,nxmax,nzmin,nzmax, X & jminus,jplus,kminus,kplus,xsam0,zsam0, X & time,wb,ul,jobname,jobsize, X & maxchk,maxhor,chkshot,horizons,xchk,zhor) Xc X1000 continue X stop X end Xc=============== subroutine init ======================== Xc X subroutine init(ier,idump,ipr,ips,dzswh,jstop, X & jcut,vector, X & nxmin,nxmax,nzmin,nzmax,xsam0,zsam0,cfl, X & maxstep,jminus,jplus,kminus,kplus, X & kzmax,jxmin,jxmax,maxsweep, X & dstep,lstep,rstep, X & time,slow, X & jobname,jobsize, X & maxchk,maxhor,chkshot,horizons,xchk,zhor) Xc Xc *** parameters (described in header of MAIN) Xc X integer ier,idump,ipr,ips,dzswh,nxmin,nxmax, X & nzmin,nzmax,maxstep,jminus,jplus,kminus,kplus, X & kzmax,jxmin,jxmax,jstop,jcut, X & maxsweep,dstep,lstep,rstep X real*4 xsam0,zsam0,cfl X real*4 time(nxmin:nxmax,nzmin:nzmax), X & slow(nxmin:nxmax,nzmin:nzmax) X integer jobsize X character*20 jobname X character*3 vector X integer chkshot,horizons X integer maxchk,maxhor X real xchk(maxchk),zhor(maxhor) Xc Xc *** workspace Xc Xc ipin input unit number (this is the only procedure Xc in the program which receives input) Xc i,j,k loop counters Xc jmin,jmax,kmin,kmax Xc x- and z- limits of initial time array, as Xc read (ascii). Xc kzmin lower z-limit of slowness array. Replaced by Xc kminus on return to MAIN. Xc test,tol workspace used in check on sample rate Xc compatibility of the two input files. Xc blank ' ' Xc junk workspace for caption lines Xc slowname name of slowness file Xc slowsize length of slowname Xc timename name of initial time file Xc timesize length of timename Xc defname name of defaults file Xc defsize length of defname Xc X integer i,j,k,ipin,kmin,kmax,jmin,jmax, X & kzmin,slowsize,timesize,defsize X real*4 test,tol,xsam1,zsam1 X character*1 blank X character*20 blank20 X character*80 slowname,timename,junk,defname Xc Xc *** twenty blanks Xc X blank20=' ' Xc Xc *** set unit numbers and defaults Xc X defsize=10 X defname(1:10)='dnodef.dat' X call setdef(ipin,ipr,ips,cfl,dzswh,maxstep,jstop, X . jcut,vector,ier,defname,defsize) X if (ier.ne.0) go to 1000 X blank=' ' Xc Xc *** if input unit is not stdin (i.e. ipin.ne.5) then Xc open job file, also set up output comment file Xc X jobname=blank20 X if ((ipr.eq.6).and.(ipin.eq.5)) then X write(6,*)' Enter job name' X read(5,2800)jobname X else if ((ipr.ne.6).and.(ipin.ne.5)) then X open(unit=ipin,file='dno.job',status='old', X . iostat=ier) X if (ier.ne.0) then X write(ipr,*)' INIT: failed to open job file' X go to 1000 X end if X open(unit=ipr,file='dno.out',status='unknown', X . iostat=ier) X if (ier.ne.0) then X write(ipr,*)' INIT: failed to open output file' X go to 1000 X end if X read(ipin,2800)jobname X else X write(ipr,*)' ERROR in INIT: impossible unit numbers' X ier=888 X go to 1000 X end if X2800 format(a20) X jobsize=index(jobname,blank) X if (jobsize.le.1) then X write(ipr,*)' ERROR in INIT: impossible job name' X ier=889 X go to 1000 X end if X jobsize=jobsize-1 Xc Xc *** say hello, Gracie... Xc X write(ipr,3562) X3562 format(' DOWNNOUT: Vidale-Enquist-Osher traveltime solver', X & /) X write(ipr,*)' INIT...' X write(ipr,*)' job:' X write(ipr,*)jobname(1:jobsize) Xc Xc get names of input slowness file Xc input source times file Xc X slowname=blank20 X junk=blank20//blank20//blank20 X junk=' Enter filename for slowness' X call filnam(slowname,slowsize,junk,ipin,ipr,ier) X if (ier.gt.0) then X go to 1000 X end if Xc X timename=blank20 X junk=blank20//blank20//blank20 X junk=' Enter filename for initial time patch' X call filnam(timename,timesize,junk,ipin,ipr,ier) X if (ier.gt.0) then X go to 1000 X end if Xc X ier=0 Xc Xc *** get job name Xc Xc X write(ipr,*)' Enter number of steps per sweep:' X write(ipr,*)' down:' X read(ipin,*)dstep X write(ipr,*)dstep X write(ipr,*)' left:' X read(ipin,*)lstep X write(ipr,*)lstep X write(ipr,*)' right:' X read(ipin,*)rstep X write(ipr,*)rstep Xc Xc *** inquire about desired checkshot profiles Xc X43 continue X write(ipr,*)' Enter number of constant-offset profiles' X write(ipr,*)' desired: > 0, < ',maxchk X read(ipin,*)chkshot X if ((chkshot.gt.maxchk).or.(chkshot.lt.0)) then X write(ipr,*)' Try again, and this time...' X write(ipr,*)' FOLLOW THE INSTRUCTIONS!!!' X go to 43 X end if X write(ipr,*)chkshot X do 44 i=1,chkshot X write(ipr,*)' Enter offset for profile ', i X read(ipin,*)xchk(i) X write(ipr,*)xchk(i) X44 continue X45 continue X write(ipr,*)' Enter number of constant-depth profiles' X write(ipr,*)' desired: > 0, < ',maxhor X read(ipin,*)horizons X if ((horizons.gt.maxhor).or.(horizons.lt.0)) then X write(ipr,*)' Try again, and this time...' X write(ipr,*)' FOLLOW THE INSTRUCTIONS!!!' X go to 45 X end if X write(ipr,*)horizons X do 46 i=1,horizons X write(ipr,*)' Enter depth for profile ', i X read(ipin,*)zhor(i) X write(ipr,*)zhor(i) X46 continue Xc Xc *** initialize time and slowness arrays Xc X if((nxmax-nxmin).lt.1) then X write(ipr,*)' ERROR: not enough points in x' X ier=3 X go to 1000 X end if X if((nzmax-nzmin).lt.1) then X write(ipr,*)' ERROR: not enough points in z' X ier=3 X go to 1000 X end if Xc X do 30 k=nzmin,nzmax X do 20 j=nxmin,nxmax X time(j,k)=0.0e+00 X slow(j,k)=0.0e+00 X20 continue X30 continue Xc Xc *** read slowness file Xc X write(ipr,*)' INIT: reading slowness file...' X write(ipr,*)' Filename = ',slowname(1:slowsize) Xc X call sreadmat(nxmin,nxmax,nzmin,nzmax, X . jxmin,jxmax,kzmin,kzmax, X . xsam0,zsam0, X . slow,slowname,slowsize,ips,ier) Xc X if (ier.gt.0) then X write(ipr,*)' INIT: error on read of slowness file' X write(ipr,*)' ier = ',ier X go to 1000 X end if Xc Xc *** read initial time data Xc X write(ipr,*)' INIT: reading initial time file...' X write(ipr,*)' Filename = ',timename(1:timesize) Xc Xc X call sreadmat(nxmin,nxmax,nzmin,nzmax, X . jmin,jmax,kmin,kmax, X . xsam1,zsam1, X . time,timename,timesize,ips,ier) Xc X if (ier.gt.0) then X write(ipr,*)' INIT: error on read of initial time file' X write(ipr,*)' ier = ',ier X go to 1000 X end if Xc Xc *** test for compatibility of sample rates Xc X test=abs(zsam0-zsam1) X tol=zsam0*(1.0e-04) X if (test.gt.tol) then X write(ipr,*)' ERROR: incompatible z sample rates' X ier=1 X go to 1000 X end if X test=abs(xsam0-xsam1) X tol=xsam0*(1.0e-04) X if (test.gt.tol) then X write(ipr,*)' ERROR: incompatible x sample rates' X ier=1 X go to 1000 X end if Xc Xc *** ascii read formats Xc X8000 format(a80) X8900 format(4(2x,e15.7)) Xc Xc *** dump Xc X if (idump.gt.0) then X write(ipr,*)' INIT: time0 data files' X write(ipr,*)' xsam0, zsam0' X write(ipr,*)xsam0,zsam0 X write(ipr,*)' jmin,jmax,kmin,kmax' X write(ipr,*)jmin,jmax,kmin,kmax X do 5500 k=kmin,kmax X write(ipr,*)' k= ',k X write(ipr,*)' time:' X write(ipr,8900)(time(j,k),j=jmin,jmax) X write(ipr,*)'---------------------------------------' X5500 continue X end if Xc Xc *** set-up for down-and-out continuation Xc X jminus=max(jmin,jxmin) X jplus =min(jmax,jxmax) X kminus= max(kmin,kzmin) X kplus= min(kmax,kzmax) X maxsweep=max(abs(jxmax),abs(jxmin)) X maxsweep=max(maxsweep,abs(kzmin)) X maxsweep=max(maxsweep,abs(kzmax)) X101 continue X if (idump.gt.0) then X write(ipr,*)' INIT: xsam0,zsam0,jminus,jplus,kminus,kplus' X write(ipr,*) xsam0,zsam0,jminus,jplus,kminus,kplus X end if Xc X1000 continue X return X end Xc Xc============================================================== Xc X subroutine setdef(ipin,ipr,ips,cfl,dzswh,maxstep,jstop, X . jcut,vector,ier,name,size) Xc------------------------------------------------------------- Xc Xc reads unit numbers and flags from file "dnodef.dat" Xc Xc------------------------------------------------------------- Xc X integer ipin,ipr,ips,dzswh,maxstep,jstop,jcut,ier,size X real cfl X character*3 vector X character*80 name X open(unit=2,file=name(1:size),status='old', X . iostat=ier,err=999) X read(2,*)ipin X read(2,*)ipr X read(2,*)ips X read(2,*)cfl X read(2,*)dzswh X read(2,*)maxstep X read(2,*)jstop X read(2,*)jcut X read(2,3000)vector X close(unit=2) X999 continue X3000 format(a3) X return X end Xc Xc================================================================== Xc X subroutine sreadmat(ld1min, ld1max, ld2min, ld2max, X . n1min, n1max, n2min, n2max, d1, d2, X . a, nameout, sizeout, iounit, X . ier) Xc Xc READMAT - A Fortran subroutine that reads a single-precision Xc matrix of uniformly spaced samples of a real function of Xc one or two real variables from an ASCII flat file Xc Xc Xc Description of arguments: Xc Xc ld1min| Xc ld1max|______> declared index limits EXACTLY as in Xc ld2min| calling program Xc ld2max| Xc Xc n1min| Xc n1max|_______> ranges of indices to be read in; Xc n2min| stored as array words 1,2,3,4 Xc n2max| Xc Xc d1|__________> sample increments (words 5,6) Xc d2| Xc Xc a - input array of real*4 or real*8 numbers, Xc regarded internally as a 2-d array (words 7,...) Xc Xc nameout - name for input file Xc sizeout - length of name Xc Xc iounit - logical unit number of input file Xc ier - error flag Xc Xc ASCII Flat File: Xc Xc All entries are stored as single-precision reals in Fortran Xc fomatted i/o. The first six entries are n1min, n1max, n2min, n2max, Xc d1 and d2 in the notation explained above. These are read and the Xc appropriate type conversions performed. The file is then re-wound Xc and the rest of the array read into memory in the given index limits. Xc Xc This routine was written by W. W. Symes, 2-25-90. Xc Standard disclaimer is hereby declared. Xc Xc================================================================ Xc Xc *** arguments: Xc Xc integer arguments: X integer*4 ld1min,ld1max,ld2min,ld2max, X . n1min,n1max,n2min,n2max, X . sizeout, iounit, ier Xc Xc real arguments: X real*4 d1,d2,a(ld1min:ld1max,ld2min:ld2max) Xc Xc character arguments: X character*80 nameout Xc Xc loop counters X integer*4 i, k Xc Xc real storage for ascii read X real*4 xs1,xs2,xs3,xs4,xs5,xs6 Xc Xc Xc write(6,*)' open file for input:' Xc X open(unit=iounit,file=nameout(1:sizeout), X . status='unknown', X . access='sequential',form='formatted', X . iostat=ier,err=992) X if (ier.ne.0) go to 992 X read(iounit,*)xs1,xs2,xs3,xs4,xs5,xs6 Xc write(6,*)' in sreadmat...' Xc write(6,*)xs1,xs2,xs3,xs4,xs5,xs6 X n1min=nint(xs1) X n1max=nint(xs2) X n2min=nint(xs3) X n2max=nint(xs4) X d1=xs5 X d2=xs6 X if ((n1min.lt.ld1min).or.(n1min.gt.ld1max)) go to 997 X if ((n2min.lt.ld2min).or.(n2min.gt.ld2max)) go to 997 X rewind iounit X read(iounit,*)xs1,xs2,xs3,xs4,xs5,xs6, X . ((a(i,k),i=n1min,n1max),k=n2min,n2max) X X close(unit=iounit) X return Xc Xc error during file opening X992 ier=8 X close(unit=iounit) X return Xc Xc impossible dimension X997 ier = 3 Xc write(6,*)'ld1min,ld1max,ld2min,ld2max' Xc write(6,*)ld1min,ld1max,ld2min,ld2max Xc write(6,*)'n1min,n1max,n2min,n2max' Xc write(6,*)n1min,n1max,n2min,n2max X close(unit=iounit) X return Xc X end Xc Xc=============================================================== Xc X subroutine swritmat(ld1min, ld1max, ld2min, ld2max, X . n1min, n1max, n2min, n2max, d1, d2, X . a, nameout, sizeout, iounit, X . ier) Xc Xc SWRITMAT - A Fortran subroutine that writes a single-precision Xc matrix of uniformly spaced samples of a real function Xc of one or two real variables to an ASCII flat file. Xc Xc Xc Description of arguments: Xc Xc ld1min| Xc ld1max|______> declared index limits EXACTLY as in Xc ld2min| calling program Xc ld2max| Xc Xc n1min| Xc n1max|_______> ranges of indices to be written out; Xc n2min| stored as array words 1,2,3,4 Xc n2max| Xc Xc d1|__________> sample increments (words 5,6) Xc d2| Xc Xc a - output array of real*4 numbers, Xc regarded internally as a 2-d array (words 7,...) Xc Xc nameout - name for output file Xc sizeout - length of name Xc Xc iounit - logical unit number of output file Xc ier - error flag Xc Xc ASCII Flat File: Xc Xc All entries are stored as single-precision reals in Fortran Xc fomatted i/o. The first six entries are n1min, n1max, n2min, n2max, Xc d1 and d2 in the notation explained above. Xc Xc This routine was written by W. W. Symes, 2-25-90. Xc Standard disclaimer is hereby declared. Xc Xc Xc================================================================ Xc Xc *** arguments: Xc Xc integer arguments: X integer*4 ld1min,ld1max,ld2min,ld2max, X . n1min,n1max,n2min,n2max, X . sizeout, iounit, ier Xc Xc real arguments: X real*4 d1,d2,a(ld1min:ld1max,ld2min:ld2max) Xc Xc character arguments: X character*80 nameout Xc Xc loop counters, limit X integer i, k, nwrit Xc Xc workspace for type conversion X real*4 xs1, xs2, xs3, xs4 Xc Xc padding for ascii records X real*4 zeroes(1:5) Xc Xc *** open file for output: Xc X open(unit=iounit,file=nameout(1:sizeout), X . status='unknown', X . access='sequential',form='formatted', X . iostat=ier,err=992) X if (ier.ne.0) go to 992 Xc Xc type conversions, padding Xc X xs1=real(n1min) X xs2=real(n1max) X xs3=real(n2min) X xs4=real(n2max) X do 800 i=1,5 X zeroes(i)=0.0e+00 X800 continue Xc Xc *** write stuff Xc X nwrit=6+(n1max-n1min+1)*(n2max-n2min+1) X nwrit=5*((nwrit/5)+1)-nwrit X write(iounit,1000)xs1,xs2,xs3,xs4,d1,d2, X . ((a(i,k),i=n1min,n1max),k=n2min,n2max), X . (zeroes(i),i=1,nwrit) X1000 format(5(2x,e12.5)) Xc Xc good write X close(unit=iounit) X return Xc Xc error during file opening X992 ier=8 X close(unit=iounit) X return Xc X end Xc Xc================= subroutine solve ======================= Xc X subroutine solve(ier,idump,ipr,dzswh,jstop,jcut, X & vector, X & nxmin,nxmax,nzmin,nzmax,xsam0,zsam0, X & maxstep,jminus,jplus,kminus,kplus, X & kzmax,jxmin,jxmax,maxsweep, X & dstep,lstep,rstep, X & cfl,time,slow,sb,slr,wb,wbtemp,tb, X & ul,ultemp,tl,ur,urtemp,tr,wl,wr,ub) Xc Xc X integer*4 ier,idump,ipr,nxmin,nxmax,nzmin,nzmax, X & jminus,jplus,kminus,kplus,kzmax,jxmin,jxmax, X & maxstep,dzswh,jstop,jcut,maxsweep, X & dstep,lstep,rstep X real*4 xsam0,zsam0,cfl X real*4 time(nxmin:nxmax,nzmin:nzmax), X & slow(nxmin:nxmax,nzmin:nzmax) X real*4 wb(nxmin:nxmax), wbtemp(nxmin:nxmax), X & wl(nzmin:nzmax),wr(nzmin:nzmax), X & ul(nzmin:nzmax),ultemp(nzmin:nzmax), X & ur(nzmin:nzmax),urtemp(nzmin:nzmax), X & ub(nxmin:nxmax), X & sb(nxmin:nxmax),slr(nzmin:nzmax), X & tb(nxmin:nxmax), X & tl(nzmin:nzmax), tr(nzmin:nzmax) X character*3 vector Xc Xc Xc GLOSSARY of Variables Xc Xc nxmin,nxmax maximal x-limits for all arrays, Xc set as parameters in MAIN Xc nzmin,nzmax maximal z-limits for all arrays, Xc set as parameters in MAIN Xc ier error flag - if > 0 on return Xc exit immediately Xc idump dump switch - see above Xc ipr output unit number (usually 6) Xc jminus,jplus x-limits of computational rectangle, Xc initialized to those of the initial Xc time data rectangle, subsequently adjusted Xc in SOLVE at the completion of each Xc successful block of steps (down, left, Xc of right). Xc kminus,kplus z-limits of the computational rectangle, Xc initialized to those of the initial Xc time data rectangle, subsequently adjusted Xc in SOLVE at the completion of each Xc successful block of steps (down, left, Xc of right). Xc dzswh fixed/variable step switch, set as default Xc in INIT. Xc dzswh = 1: fixed steps (dangerous!) Xc else: adaptive steps Xc jstop switch for stopping at a corner inflow Xc condition. In order to proceed, the Xc flow must be OUT at each end of the faces Xc of the rectangular compuational front, Xc else an upwind step cannot be taken there. Xc jstop = 0: drop the endpoint in question, Xc continue computing on smaller Xc face ( this option only works Xc in the present implementation Xc for marching in ONE DIRECTION Xc ONLY, e.g. down). Prevents com- Xc putation of full rectangle of Xc times. Xc jstop > 0: stop computation in current Xc direction when inflow is Xc detected. This option is Xc consistent with the design of Xc the algorithm. Xc jcut flux cutoff switch: Xc jcut = 0: step computed adaptively Xc jcut ><0: current flow variable (x or Xc z gradient component) truncated Xc to limit signal velocity to 95% Xc of the maximum permissible. Xc This option Xc ensures that the calculation Xc will go to completion even when Xc it shouldn't. Xc maxstep max. number if partial steps allowed - set as Xc default in INIT Xc jxmin,jxmax x-limits of computational domain, set during Xc reading of slowness array. Xc kzmax max. z-limit of computational domain, set Xc during reading of slowness array. Note that Xc the minimum z-limit remains at its initial Xc value, i.e. kminus - down and out, but not Xc up! Xc maxsweep max. number of "down, left, right" cycles Xc permitted Xc dstep, lstep, rstep Xc down, left, and right steps per cycle (see Xc above). Xc xsam0, zsam0 x and z steps as read from "slow.dat", and Xc checked against "time0.dat". Xc cfl Courant-Friedrichs-Lewy factor, used in Xc setting adaptive step. Set as default in Xc INIT. Xc time array for computed time values Xc slow array for input slowness values Xc Xc *** the next set of variables provide storage for Xc the time field and its gradient components on Xc the three parts (bottom, left, and right) of Xc the computational front (no top component - Xc again, "down and out"!). For reasons we won't Xc go into, w is used to denote the x-derivative Xc and u the z-derivative, to which are added "b", Xc "l", and "r" as appropriate. For some arrays we Xc also keep a "**temp" array for updated values, Xc to avoid vectorization-destroying dependencies. Xc Xc wb, wbtemp, ub, tb Xc ul, ultemp, wl, tl Xc ur, urtemp, wr, tr Xc Xc sb, slr temporary storage for interpolated Xc slowness values on bottom and Xc left/right fronts, respectively. Xc Xc Xc---------------------------------------------------------------------------- Xc Xc Subroutine SOLVE: Xc Xc Vectorized version of numerical solution of the Eikonal equation of Xc geometric optics: Xc Xc ( DT(x,z)/Dx )**2 + ( DT(x,z)/Dz )**2 = ( s(x,z) )**2 Xc Xc where x is the horizonatal coordinate in a box, Xc and z is the vertical coordinate in a box, Xc and T is the ray traveltime, Xc and s is the slowness ( 1/velocity ). Xc Xc Let w(x,z) = DT(x,z)/Dx Xc Xc Let u(x,z) = DT(x,z)/Dz Xc Xc Then, w**2 + u**2 = s**2 Xc Xc if rays are moving downward, then Xc Xc u = + sqrt( s**2 - w**2 ) Xc Xc apply D/Dx to both sides and use Du/Dx = Dw/Dz Xc Xc---------------------------------------------------------------------------- Xc Xc Depth equation: rays going down. Xc Xc Dw/Dz + Df(w)/Dx = 0 Xc Xc f(w) = - sqrt( s**2 - w**2 ) Xc Xc---------------------------------------------------------------------------- Xc Xc Right equation: rays going to right: Xc Xc Du/dx + Df(u)/dz = 0 Xc Xc f(u) = - sqrt( s**2 - u**2 ) Xc Xc---------------------------------------------------------------------------- Xc Xc left equation: rays going to left: Xc Xc - Du/dx + Df(u)/dz = 0 Xc Xc f(u) = - sqrt( s**2 - u**2 ) Xc Xc---------------------------------------------------------------------------- Xc Xc Xc *** workspace Xc Xc j,jx,k,kz,isweep loop counters Xc *minus1 (*=j,k) *minus+1 Xc *plus1 *plus-1 Xc idun* (*=b,l,r) completion flag ("I-done-with-...") Xc for the three computational directions Xc jnext,knext temporary limits for x, z steps in each Xc cycle Xc xsam,zsam partial steps, determined adaptively Xc slope rate of change of primary conservation Xc law variable (w on bottom, u on sides) Xc Xc *** register variables Xc Xc flux* (*=l,m,r) conservation law flux (i.e. the other Xc gradient component) at point to left, Xc current point, point to right Xc fmin* (*=l,m,r) flux(min(u or w,0)) at left, middle, Xc or right point - construct in EO scheme Xc fmax* (*=l,m,r) similar Xc s* (*=l,m,r) similar, for slowness field Xc hs workspace for Heaviside function Xc one what do you think? Xc zero one guess Xc fluxeol, fluxeor left and right Enquist-Osher numerical Xc fluxes Xc %#* (% = w, u register variables for nearby values Xc # = b,l,r of the gradient components on the Xc * = l,m,r) three computational fronts Xc uc,up,wc,wp storage for use in update loops Xc theta1,theta2 piecewise-linear interpolation factors Xc xc,zc coordinates for current step Xc xp,zp coordinates for next step Xc x,z coordintes for partial step Xc umin,tolstep tolerances used to avoid Xc continuing computation when Xc ray is detected parallel to the Xc computational front. Xc umax,wmax cutoffs used to avoid constructing rays Xc parallel to the computational front Xc tol,test workspace for testing the adherence Xc of the initial gradient components Xc to the Eikonal equation. Xc cmax workspace used in CFL adaptive step Xc computation (as is test) Xc X integer j,k,jx,kz,jminus1,jplus1,isweep,idunb, X & idunl,idunr,jnext,knext,kplus1,kminus1 X real*4 zsam,xsam,slope,fluxl,fluxm,fluxr,fminl,fmaxl, X & fminm,fmaxm,fminr,fmaxr,fluxeol,fluxeor, X & z,zc,zp,x,xc,xp,umin,umax,wmax,tol,tolstep, X & wbr,wbm,wbl,ull,ulm,ulr, url,urm,urr, X & sm,sl,test,cmax,hs,one,zero, X & uc,up,wc,wp,theta1,theta2 Xc Xc *** say hello, Gracie... Xc X write(ipr,*)' SOLVE...' Xc Xc *** set defaults Xc X ier=0 X one=1.0e+00 X zero=0.0e+00 X umin=1.0e-08 X tolstep=1.0e-05 X wmax=0.97 X umax=0.97 X Xc Xc *** initialize work arrays: Xc load time, Xc create approximate gradient, check for gross errors Xc i. e. tachyonic rays Xc Xc bottom: initialize tb as times on bottom of rectange Xc ub as DT/Dz Xc Xc jminus, jplus are current limits in x Xc kminus, kplus are current limits in z Xc Xc use upwind depth z-derivative Xc X do 6000 j=jminus,jplus X tb(j)=time(j,kplus) X ub(j)=(time(j,kplus)-time(j,kplus-1))/zsam0 X6000 continue Xc Xc initialize wb as DT/Dx on bottom Xc use one-side numerical derivatives on the edges Xc X wb(jminus)=(time(jminus+1,kplus) X & - time(jminus,kplus))/xsam0 X wb(jplus)=(time(jplus,kplus)-time(jplus-1,kplus))/xsam0 Xc Xc use centered x-difference inside: Xc X do 6010 j=jminus+1,jplus-1 X wb(j)=0.5*(time(j+1,kplus)-time(j-1,kplus))/xsam0 X6010 continue Xc Xc *** test for bad rays Xc X do 6020 j=jminus,jplus X test=(slow(j,kplus)-ub(j))*(slow(j,kplus)+ub(j)) X tol=-xsam0*slow(j,kplus)*slow(j,kplus) X if (test.lt.tol) then X write(ipr,*)' ERROR in INIT: ' X write(ipr,*)' initial gradient inadmissible' X write(ipr,*)' On bottom, j = ',j X write(ipr,*)' gradz = ', ub(j) X ier=999 X go to 9999 X end if X test=(slow(j,kplus)-wb(j))*(slow(j,kplus)+wb(j)) X tol=-zsam0*slow(j,kplus)*slow(j,kplus) X if (test.lt.tol) then X write(ipr,*)' ERROR in INIT: ' X write(ipr,*)' initial gradient inadmissible' X write(ipr,*)' On bottom, j = ',j X write(ipr,*)' gradx = ',wb(j) X ier=999 X go to 9999 X end if X6020 continue Xc Xc left: initialize tl as times on left side of rectangle Xc wl as DT/dx Xc Xc use upwind left x-derivative Xc X do 6100 k=kminus,kplus X tl(k)=time(jminus,k) X wl(k)=(time(jminus+1,k)-time(jminus,k))/xsam0 X6100 continue Xc Xc initialize ul as DT/Dz on left Xc use one-sided numerical derivatives on the edges Xc X ul(kminus)=(time(jminus,kminus+1) X & - time(jminus,kminus))/zsam0 X ul(kplus)=(time(jminus,kplus)-time(jminus,kplus-1))/zsam0 Xc Xc use centered-difference z-derivatives inside Xc X do 6110 k=kminus+1,kplus-1 X ul(k)=0.5*(time(jminus,k+1)-time(jminus,k-1))/zsam0 X6110 continue Xc Xc *** test for bad rays Xc X do 6120 k=kminus,kplus X test=(slow(jminus,k)-ul(k))*(slow(jminus,k)+ul(k)) X tol=-zsam0*slow(jminus,k)*slow(jminus,k) X if (test.lt.tol) then X write(ipr,*)' ERROR in INIT: ' X write(ipr,*)' initial gradient inadmissible' X write(ipr,*)' On left, k = ',k X write(ipr,*)' gradz = ', ul(k) X ier=999 X go to 9999 X end if X test=(slow(jminus,k)-wl(k))*(slow(jminus,k)+wl(k)) X tol=-zsam0*slow(jminus,k)*slow(jminus,k) X if (test.lt.tol) then X write(ipr,*)' ERROR in INIT: ' X write(ipr,*)' initial gradient inadmissible' X write(ipr,*)' On left, k = ',k X write(ipr,*)' gradx = ',wl(k) X ier=999 X go to 9999 X end if X6120 continue Xc Xc right: initialize tr as times on right side of rectangle Xc wr as DT/Dx Xc Xc use upwind right x-derivative Xc X do 6200 k=kminus,kplus X tr(k)=time(jplus,k) X wr(k)=(time(jplus,k)-time(jplus-1,k))/xsam0 X6200 continue Xc Xc initialize ur as DT/dz on right Xc use one-sided numerical derivatives on edges Xc X ur(kminus)=(time(jplus,kminus+1) X & - time(jplus,kminus))/zsam0 X ur(kplus)=(time(jplus,kplus)-time(jplus,kplus-1))/zsam0 Xc Xc use centered difference inside: Xc X do 6210 k=kminus+1,kplus-1 X ur(k)=0.5*(time(jplus,k+1)-time(jplus,k-1))/zsam0 X6210 continue Xc Xc *** test for bad rays Xc X do 6220 k=kminus,kplus X test=(slow(jplus,k)-ur(k))*(slow(jplus,k)+ur(k)) X tol=-zsam0*slow(jplus,k)*slow(jplus,k) X if (test.lt.tol) then X write(ipr,*)' ERROR in INIT: ' X write(ipr,*)' initial gradient inadmissible' X write(ipr,*)' On right, k = ',k X write(ipr,*)' gradz = ', ur(k) X ier=999 X go to 9999 X end if X test=(slow(jplus,k)-wr(k))*(slow(jplus,k)+wr(k)) X tol=-zsam0*slow(jplus,k)*slow(jplus,k) X if (test.lt.tol) then X write(ipr,*)' ERROR in INIT: ' X write(ipr,*)' initial gradient inadmissible' X write(ipr,*)' On right, k = ',k X write(ipr,*)' gradx = ',wr(k) X ier=999 X go to 9999 X end if X6220 continue Xc Xc *** main loop: Xc Xc will try to evolve outward from initial source-box Xc taking pre-specified numbers of down steps, left steps, Xc right steps. Xc X idunb=0 X idunl=0 X idunr=0 Xc X do 5000 isweep = 1,maxsweep Xc Xc**************** 5000 loop ************************ Xc* * Xc* *** if we haven't reached the bottom, * Xc* attemp dstep steps down: * Xc* * Xc**************** start of down loop *************** Xc X if (dstep.le.0) then X idunb=1 X go to 1999 X end if Xc Xc *** conditional entrance to down loop Xc X if ((kplus.lt.kzmax).and.(idunb.eq.0)) then Xc X knext=min(kzmax,kplus+dstep) Xc X do 1000 kz=kplus+1,knext Xc Xc 0. set current depth Xc X zc = float(kz-1)*zsam0 X zp = float(kz)*zsam0 X z = zc Xc Xc *** set slowness vector Xc X do 20 j=jminus,jplus X sb(j)=slow(j,kz-1) X20 continue Xc Xc *** take as many intermediate steps as necessary Xc X do 100 k=1,maxstep Xc Xc *** quit if we've made it to zp Xc X if (((zp-z)/zsam0).lt.tolstep) then X if (idump.ge.0) then X write(ipr,*)' SOLVE: down step ',kz,' completed' X end if X go to 800 X end if Xc Xc *** set local depth step according to CFL Xc Note: if jcut >< 0, this yields the largest step Xc consistent with stability, which is always at Xc least 1/maxstep of the total step Xc Xc for depth stepping CFL checks Max( w/u ) * Dz/Dx < 1 Xc (depth) Xc Xc will put floor under ub of umin. Xc floor under cmax of 1.0e-3 Xc X cmax=1.0e-3 X do 15 j=jminus,jplus X wbm=wb(j) X sm=sb(j) X fluxm=ub(j) X fluxm=max(umin,fluxm) X test= abs( wbm/fluxm ) X cmax= max(cmax,test) X 15 continue Xc X test= cfl*xsam0/cmax X if (dzswh.eq.1) then X zsam=zsam0 X if (zsam.gt.test) then X write(ipr,*)' SOLVE: fixed step too large, step',kz X ier=1 X go to 9999 X end if X else X zsam = min((zp-z),test) X end if Xc Xc *** step impossible, error exit Xc X if (zsam.lt.(0.5*tolstep*zsam0)) then X write(ipr,*)' SOLVE: down step ',kz X write(ipr,*)' internal step ',k X write(ipr,*)' zsam < tolstep*zsam0' X ier=0 X idunb=1 X kplus=kz-1 X go to 1998 X end if Xc X j =jminus+2 X if (j.gt.jplus) then X write(ipr,*)' SOLVE: run out of room' X ier=0 X idunb=1 X kplus=kz-1 X go to 1998 X end if Xc Xc *** step wb in z Xc Xc 1. check endpoint slopes, update end values if possible Xc X jminus1=jminus+1 X jplus1=jplus-1 Xc Xc check for inflow at right, wb(xmax) < 0 ? Xc X wbr=wb(jplus) X if (wbr.lt.zero) then X if (jstop.gt.0) then X write(ipr,*)' SOLVE: inflow at right' X write(ipr,*)' kz = ',kz,' k = ',k X ier=0 X idunr=1 X idunb=1 X kplus=kz-1 X go to 1998 X else X wbtemp(jplus)=zero X jplus=jplus-1 X jplus1=jplus-1 X end if X else Xc Xc no inflow - step w at xmax forward in z: Xc X fluxr=ub(jplus) X fluxm=ub(jplus1) X slope = (fluxr-fluxm)/xsam0 X wbtemp(jplus)=wbr+zsam*slope X end if X if (idump.gt.0) then X write(ipr,*)' SOLVE: slope @ jplus = ', slope X end if Xc X wbl=wb(jminus) X if (wbl.gt.zero) then X if (jstop.gt.0) then X write(ipr,*)' SOLVE: inflow at left' X write(ipr,*)' kz = ',kz,' k = ',k X idunl=1 X idunb=1 X ier=0 X kplus=kz-1 X go to 1998 X else X wbtemp(jminus)=zero X jminus=jminus+1 X jminus1=jminus+1 X wbl=wb(jminus) X end if X else Xc Xc no inflow at left - step w(xmin) forward in z: Xc X fluxm=ub(jminus1) X fluxl=ub(jminus) X slope=(fluxm-fluxl)/xsam0 X wbtemp(jminus)=wbl+zsam*slope X end if X if (idump.gt.0) then X write(ipr,*)' SOLVE: slope @ jminus = ', slope X end if Xc Xc 2. step interior values... Xc Believe it or not, this is an Xc implementation of the E/O scheme. Note that the only Xc value of slowness used in the definitions is the Xc value at the center point. This "frozen coefficient" Xc scheme is first-order accurate. More straightforward Xc application of the E/O formulas actually gives an Xc inconsistent scheme. Xc Xc=============================================================== Xc Xc Some equations, courtesy of J. J. Carazzone: Xc Xc Dw/Dz - DFlux(w)/Dx = 0 Xc Xc Flux(w) = sqrt( s**2 - w**2 ) = u(w) Xc Xc Flux(0) = s Xc Xc E-O Scheme: Xc Xc Xc f(w) = sqrt( s**2 - w**2 ) Xc Xc Df(w)/Dx = ( s * s' - w * w' )/f(w) Xc Xc Assume s' = Ds/Dx <<<< w' = Dw/Dx Xc Xc so w=0 makes Df(w)/Dx = 0 (approximately) Xc Xc w=0 is the One and Only stagnation point! Xc Xc w(j,k+1) = Xc w(j,k) + zsam/xsam * ( D+( Flux-( w(j,k) ) + D-( Flux+( w(j,k) ) ) ) Xc Xc Xc D+( f ) = f(j+1) - f(j) Xc Xc D-( f ) = f(j) - f(j-1) Xc Xc Flux+( w ) = Flux( max(w,0) ) Xc Xc Flux-( w ) = Flux( min(w,0) ) Xc Xc= w(j,k) + zsam/xsam * ( Flux( min(w(j+1,k),0) ) - Flux( min(w(j,k),0) ) Xc + Flux( max(w(j,k),0) ) - Flux( max(w(j-1,k),0) )) Xc Xc= w(j,k)+zsam/xsam * ( [ Flux( min(w(j+1,k),0) ) + Flux( max(w(j,k),0) ) ] Xc - [ Flux( min(w(j,k),0) ) + Flux( max(w(j-1,k),0) ) ]) Xc Xc= w(j,k)+zsam/xsam * ( fluxeor - fluxeol ) Xc Xc=============================================================== Xc Xc Scalar vs. Vector loops: Xc Xc We provide two ways of computing the main loop, toggled by the Xc parameter 'vector': Xc scalar - uses scalar register variables, should be more efficient Xc on ix86, motorola, and similar machines (vector='no ') Xc vector - uses vector array references directly, for machines with Xc vector registers and compilers like CFT (vector='yes') Xc Xc *** scalar form Xc X if (vector(1:3).eq.'no ') then Xc X wbl=wb(jminus) X fluxl=ub(jminus) X wbm=wb(jminus1) X fluxm=ub(jminus1) X do 52 j=jminus1,jplus1 X wbr=wb(j+1) X fluxr=ub(j+1) X sm=sb(j) X hs=0.5*(1.0+sign(one,wbl)) X fminl=hs*sm+(1.0-hs)*fluxl X fmaxl=hs*fluxl+(1.0-hs)*sm X hs=0.5*(1.0+sign(one,wbm)) X fminm=hs*sm+(1.0-hs)*fluxm X fmaxm=hs*fluxm+(1.0-hs)*sm X hs=0.5*(1.0+sign(one,wbr)) X fminr=hs*sm+(1.0-hs)*fluxr X fmaxr=hs*fluxr+(1.0-hs)*sm X fluxeol=fminm+fmaxl X fluxeor=fminr+fmaxm X slope=(fluxeor-fluxeol)/xsam0 X wbtemp(j)=wbm+zsam*slope X wbl=wbm X fluxl=fluxm X wbm=wbr X fluxm=fluxr X52 continue Xc Xc *** vector form Xc X else if (vector(1:3).eq.'yes') then Xc X do 152 j=jminus1,jplus1 X hs=0.5*(1.0+sign(one,wb(j-1))) X fminl=hs*sb(j)+(1.0-hs)*ub(j-1) X fmaxl=hs*ub(j-1)+(1.0-hs)*sb(j) X hs=0.5*(1.0+sign(one,wb(j))) X fminm=hs*sb(j)+(1.0-hs)*ub(j) X fmaxm=hs*ub(j)+(1.0-hs)*sb(j) X hs=0.5*(1.0+sign(one,wb(j+1))) X fminr=hs*sb(j)+(1.0-hs)*ub(j+1) X fmaxr=hs*ub(j+1)+(1.0-hs)*sb(j) X fluxeol=fminm+fmaxl X fluxeor=fminr+fmaxm X slope=(fluxeor-fluxeol)/xsam0 X wbtemp(j)=wb(j)+zsam*slope X152 continue Xc Xc *** nonexistant option Xc X else Xc X ier=88 X write(ipr,*)' SOLVE: improper vector option, 52 loop' X return Xc X end if Xc Xc Xc *** update depth, tb, wb, ub at left endpoint Xc also update slowness for next depth step Xc X z=z+zsam X theta1=(z-zc)/zsam0 X theta2=(zp-z)/zsam0 X do 55 j=jminus,jplus X sb(j)=theta1*slow(j,kz) + theta2*slow(j,kz-1) X55 continue Xc Xc *** update wb, ub, integrate ub=DT/Dz to get tb using Simpson's rule Xc X test = zero X if (jcut.eq.0) then X do 60 j=jminus,jplus X sl=sb(j) X wbl=wbtemp(j) X up=(sl-wbl)*(sl+wbl) X test=min(up,test) X up=sqrt(abs(up)) X uc=ub(j) X ub(j)=up X tb(j)=tb(j)+0.5*zsam*(up+uc) X wb(j)=wbl X60 continue X else X do 61 j=jminus,jplus X sl=sb(j) X wbl=min(wbtemp(j),wmax*sl) X wbl=max(wbl,-wmax*sl) X up=(sl-wbl)*(sl+wbl) X test=min(up,test) X up=sqrt(abs(up)) X uc=ub(j) X ub(j)=up X tb(j)=tb(j)+0.5*zsam*(up+uc) X wb(j)=wbl X61 continue X end if Xc Xc *** debug dump Xc X if (idump.gt.0) then X write(ipr,*)' SOLVE: down step ',kz X write(ipr,*)' SOLVE: internal step ',k X write(ipr,*)' xsam0, zsam0, zsam, z, zc, zp:' X write(ipr,*) xsam0, zsam0, zsam, z, zc, zp X write(ipr,*)'jminus,jplus,jminus1,jplus1:' X write(ipr,*)jminus,jplus,jminus1,jplus1 X write(ipr,*)'(ub(j),j=jminus,jplus) at depth ',z X write(ipr,*)' ' X write(ipr,4998)jminus,jplus,xsam0 X write(ipr,4999)(ub(j),j=jminus,jplus) X write(ipr,*)'(wb(j),j=jminus,jplus) at depth ',z X write(ipr,*)' ' X write(ipr,4998)jminus,jplus,xsam0 X write(ipr,4999)(wb(j),j=jminus,jplus) X write(ipr,*)'(tb(j),j=jminus,jplus) at depth ',z X write(ipr,*)' ' X write(ipr,4998)jminus,jplus,xsam0 X write(ipr,4999)(tb(j),j=jminus,jplus) X end if Xc Xc *** horizontal ray condition Xc X if (test.lt.zero) then X write(ipr,*)' SOLVE: detects horizontal ray' X write(ipr,*)' down step ',kz X write(ipr,*)' internal step ',k X write(ipr,*)' depth ',z X ier=0 X idunb=1 X kplus=kz-1 X go to 1998 X end if Xc X100 continue Xc Xc *** didn't make the step Xc X write(ipr,*)' SOLVE: failed to complete down step ' X write(ipr,*)' step ',kz X kplus=kz-1 X idunb=1 X go to 1998 Xc X800 continue Xc Xc 3. update time on original grid Xc X do 900 j=jminus,jplus X time(j,kz)=tb(j) X900 continue Xc Xc *** pick off new values for left, right arrays Xc X wl(kz)=wb(jminus) X ul(kz)=ub(jminus) X tl(kz)=tb(jminus) X wr(kz)=wb(jplus) X ur(kz)=ub(jplus) X tr(kz)=tb(jplus) Xc X1000 continue Xc Xc *** reset count of successfully completed down steps Xc X kplus=knext Xc X1998 continue Xc X end if Xc X1999 continue Xc Xc**************** end of depth loop ***************************** Xc* * Xc* *** if we haven't reached the left side, attemp * Xc* lstep steps left: * Xc* * Xc**************** start of left loop **************************** Xc X if (lstep.le.0) then X idunl=1 X go to 2999 X end if Xc Xc *** conditional entrance to left loop: Xc X if ((jminus.gt.jxmin).and.(idunl.eq.0)) then Xc X jnext=max(jxmin,jminus-lstep) Xc X do 2000 jx=jminus-1,jnext,-1 Xc Xc 0. set current left boundary Xc X xc = float(jx+1)*xsam0 X xp = float(jx)*xsam0 X x = xc Xc Xc *** set slowness vector Xc X do 1020 k=kminus,kplus X slr(k)=slow(jx+1,k) X1020 continue Xc Xc *** take as many intermediate steps as necessary Xc X do 1100 j=1,maxstep Xc Xc *** quit if we've made it to xp Xc X if (((xp-x)/xsam0).gt.-tolstep) then X if (idump.ge.0) then X write(ipr,*)' SOLVE: left step ',jx,' completed' X end if X go to 1800 X end if Xc Xc *** set local depth step according to CFL Xc Note: if jcut >< 0, this yields the largest step Xc consistent with stability, which is always at Xc least 1/maxstep of the total step Xc X cmax=1.0e-3 X do 1015 k=kminus,kplus X ulm=ul(k) X fluxm=wl(k) X fluxm=max(umin,abs(fluxm)) X test= abs( ulm/fluxm ) X cmax= max(cmax,test) X1015 continue Xc X test= cfl*zsam0/cmax X if (dzswh.eq.1) then X xsam=xsam0 X if (xsam.gt.test) then X write(ipr,*)' SOLVE: fixed left step too large' X write(ipr,*)' step',jx X ier=1 X go to 9999 X end if X else X xsam = min((x-xp),test) X end if Xc Xc *** step impossible, error exit Xc X if (xsam.lt.(0.5*tolstep*xsam0)) then X write(ipr,*)' SOLVE: left step ',jx X write(ipr,*)' internal step ',j X write(ipr,*)' xsam < tolstep*xsam0' X ier=0 X idunl=1 X jminus=jx+1 X go to 2998 X end if Xc X k =kminus+2 X if (k.gt.kplus) then X write(ipr,*)' SOLVE: run out of room on left' X ier=0 X go to 2998 X end if Xc Xc *** step ul in x Xc Xc 1. check endpoint slopes, update end values if possible Xc X kminus1=kminus+1 X kplus1=kplus-1 Xc X ull=ul(kplus) X if (ull.lt.zero) then X if (jstop.gt.0) then X write(ipr,*)' SOLVE: inflow at bottom left' X write(ipr,*)' jx = ',jx,' j = ',j X idunl=1 X idunb=1 X jminus=jx+1 X go to 2998 X else X ultemp(kplus)=zero X kplus=kplus-1 X kplus1=kplus-1 X end if X else X fluxl=wl(kplus) X fluxm=wl(kplus1) X slope = (fluxm-fluxl)/zsam0 X ultemp(kplus)=ull+xsam*slope X end if X if (idump.gt.0) then X write(ipr,*)' SOLVE: slope @ kplus = ', slope X write(ipr,*)' updated ul value = ',ultemp(kplus) X end if Xc X ulr=ul(kminus) X if (ulr.gt.zero) then X if (jstop.gt.0) then X write(ipr,*)' SOLVE: inflow at top left' X write(ipr,*)' jx = ',jx,' j = ',j X idunl=1 X jminus=jx+1 X go to 2998 X else X ultemp(kminus)=zero X kminus=kminus+1 X kminus1=kminus+1 X ulr=ul(kminus) X end if X else X fluxm=wl(kminus1) X fluxr=wl(kminus) X slope=(fluxr-fluxm)/zsam0 X ultemp(kminus)=ulr+xsam*slope X end if X if (idump.gt.0) then X write(ipr,*)' SOLVE: slope @ kminus = ', slope X write(ipr,*)' updated ul value = ', ul(kminus) X end if Xc Xc 2. step interior values Xc Xc *** scalar form Xc X if (vector(1:3).eq.'no ') then Xc X ulr=ul(kminus) X fluxr=wl(kminus) X ulm=ul(kminus1) X fluxm=wl(kminus1) X do 1052 k=kminus1,kplus1 X ull=ul(k+1) X fluxl=wl(k+1) X sm=slr(k) X hs=0.5*(1.0+sign(one,ulr)) X fminr=-hs*sm+(1.0-hs)*fluxr X fmaxr=hs*fluxr-(1.0-hs)*sm X hs=0.5*(1.0+sign(one,ulm)) X fminm=-hs*sm+(1.0-hs)*fluxm X fmaxm=hs*fluxm-(1.0-hs)*sm X hs=0.5*(1.0+sign(one,ull)) X fminl=-hs*sm+(1.0-hs)*fluxl X fmaxl=hs*fluxl-(1.0-hs)*sm X fluxeol=fminl+fmaxm X fluxeor=fminm+fmaxr X slope=(fluxeor-fluxeol)/zsam0 X ultemp(k)=ulm+xsam*slope X fluxr=fluxm X fluxm=fluxl X ulr=ulm X ulm=ull X1052 continue Xc Xc *** vector form Xc X else if (vector(1:3).eq.'yes') then Xc X do 1152 k=kminus1,kplus1 X hs=0.5*(1.0+sign(one,ul(k-1))) X fminr=-hs*slr(k)+(1.0-hs)*wl(k-1) X fmaxr=hs*wl(k-1)-(1.0-hs)*slr(k) X hs=0.5*(1.0+sign(one,ul(k))) X fminm=-hs*slr(k)+(1.0-hs)*wl(k) X fmaxm=hs*wl(k)-(1.0-hs)*slr(k) X hs=0.5*(1.0+sign(one,ul(k+1))) X fminl=-hs*slr(k)+(1.0-hs)*wl(k+1) X fmaxl=hs*wl(k+1)-(1.0-hs)*slr(k) X fluxeol=fminl+fmaxm X fluxeor=fminm+fmaxr X slope=(fluxeor-fluxeol)/zsam0 X ultemp(k)=ul(k)+xsam*slope X1152 continue Xc Xc *** nonexistant vector option Xc X else Xc X ier=89 X write(ipr,*)' SOLVE: improper vector option in 1052 loop' X return Xc X end if Xc Xc *** update offset, slowness Xc X x=x-xsam X theta1=-(x-xc)/xsam0 X theta2=-(xp-x)/xsam0 X do 1055 k=kminus,kplus X slr(k)=theta1*slow(jx,k) + theta2*slow(jx+1,k) X1055 continue Xc Xc *** integrate ultemp in z to get time, update ul, wl Xc X test = zero X if (jcut.eq.0) then X do 1060 k=kminus,kplus X sl=slr(k) X ull=ultemp(k) X wp=(sl-ull)*(sl+ull) X test=min(wp,test) X wp=-sqrt(abs(wp)) X wc=wl(k) X wl(k)=wp X tl(k)=tl(k)-0.5*xsam*(wp+wc) X ul(k)=ull X1060 continue X else X do 1061 k=kminus,kplus X sl=slr(k) X ull=min(ultemp(k),umax*sl) X ull=max(ull,-umax*sl) X wp=(sl-ull)*(sl+ull) X test=min(wp,test) X wp=-sqrt(abs(wp)) X wc=wl(k) X wl(k)=wp X tl(k)=tl(k)-0.5*xsam*(wp+wc) X ul(k)=ull X1061 continue X end if Xc Xc *** debug dump Xc X if (idump.gt.0) then X write(ipr,*)' SOLVE: left step ',jx X write(ipr,*)' SOLVE: internal step ',j X write(ipr,*)' xsam0, zsam0, xsam, x, xc, xp:' X write(ipr,*) xsam0, zsam0, xsam, x, xc, xp X write(ipr,*)'kminus,kplus,kminus1,kplus1:' X write(ipr,*)kminus,kplus,kminus1,kplus1 X write(ipr,*)'(ul(k),k=kminus,kplus) at offset ',x X write(ipr,*)' ' X write(ipr,4998)kminus,kplus,zsam0 X write(ipr,4999)(ul(k),k=kminus,kplus) X write(ipr,*)'(wl(k),k=kminus,kplus) at offset ',x X write(ipr,*)' ' X write(ipr,4998)kminus,kplus,zsam0 X write(ipr,4999)(wl(k),k=kminus,kplus) X write(ipr,*)'(tl(k),k=kminus,kplus) at offset ',x X write(ipr,*)' ' X write(ipr,4998)kminus,kplus,zsam0 X write(ipr,4999)(tl(k),k=kminus,kplus) X end if Xc Xc *** vertical ray condition Xc X if (test.lt.zero) then X write(ipr,*)' SOLVE: detects vertical ray at left' X write(ipr,*)' left step ',jx X write(ipr,*)' internal step ',j X write(ipr,*)' offset ',x X ier=0 X idunl=1 X jminus=jx+1 X go to 2998 X end if Xc X1100 continue Xc X write(ipr,*)' SOLVE: failure to complete left step' X write(ipr,*)' step = ',jx X idunl=1 X jminus=jx+1 X go to 2998 Xc X1800 continue Xc Xc 3. update time, gradient Xc X do 1900 k=kminus,kplus X time(jx,k)=tl(k) X1900 continue Xc Xc *** pick off new values for bottom arrays Xc X wb(jx)=wl(kplus) X ub(jx)=ul(kplus) X tb(jx)=tl(kplus) Xc X2000 continue Xc Xc *** reset count of successfully completed left steps Xc X jminus=jnext Xc X2998 continue Xc X end if Xc X2999 continue Xc Xc*************** end of left step ********************* Xc* * Xc* *** if possible take a step to the right * Xc* * Xc*************** start of right step ****************** Xc Xc X if (rstep.le.0) then X idunr=1 X go to 3999 X end if Xc Xc *** conditional entrance to right loop: Xc X if ((jplus.lt.jxmax).and.(idunr.eq.0)) then Xc X jnext=min(jxmax,jplus+rstep) Xc X do 3000 jx=jplus+1,jnext Xc Xc 0. set current right boundary Xc X xc = float(jx-1)*xsam0 X xp = float(jx)*xsam0 X x = xc Xc Xc *** set slowness vector Xc X do 2020 k=kminus,kplus X slr(k)=slow(jx-1,k) X2020 continue Xc Xc *** take as many intermediate steps as necessary Xc Note: if jcut >< 0, this yields the largest step Xc consistent with stability, which is always at Xc least 1/maxstep of the total step Xc X do 2100 j=1,maxstep Xc Xc *** quit if we've made it to xp Xc X if (((xp-x)/xsam0).lt.tolstep) then X if (idump.ge.0) then X write(ipr,*)' SOLVE: right step ',jx,' completed' X end if X go to 2800 X end if Xc Xc *** set local depth step according to CFL Xc X cmax=1.0e-3 X do 2015 k=kminus,kplus X urm=ur(k) X fluxm=wr(k) X fluxm=max(umin,abs(fluxm)) X test= abs( urm/fluxm ) X cmax= max(cmax,test) X2015 continue Xc X test= cfl*zsam0/cmax X if (dzswh.eq.1) then X xsam=xsam0 X if (xsam.gt.test) then X write(ipr,*)' SOLVE: fixed right step too large' X write(ipr,*)' step',jx X ier=1 X go to 9999 X end if X else X xsam = min((xp-x),test) X end if Xc Xc *** step impossible, error exit Xc X if (xsam.lt.(0.5*tolstep*xsam0)) then X write(ipr,*)' SOLVE: right step ',jx X write(ipr,*)' internal step ',j X write(ipr,*)' xsam < tolstep*xsam0' X ier=0 X idunr=1 X jplus=jx-1 X go to 3998 X end if Xc X k =kminus+2 X if (k.gt.kplus) then X write(ipr,*)' SOLVE: run out of room on right' X ier=0 X go to 3998 X end if Xc Xc *** step ur in x Xc Xc 1. check endpoint slopes, update end values if possible Xc X kminus1=kminus+1 X kplus1=kplus-1 Xc X urr=ur(kplus) X if (urr.lt.zero) then X if (jstop.gt.0) then X write(ipr,*)' SOLVE: inflow at bottom right' X write(ipr,*)' jx = ',jx,' j = ',j X idunr=1 X idunb=1 X jplus=jx-1 X go to 3998 X else X urtemp(kplus)=zero X kplus=kplus-1 X kplus1=kplus-1 X end if X else X fluxm=wr(kplus) X fluxl=wr(kplus1) X slope = (fluxm-fluxl)/zsam0 X urtemp(kplus)=urr+xsam*slope X end if X if (idump.gt.0) then X write(ipr,*)' SOLVE: slope @ kplus = ', slope X write(ipr,*)' updated ur value = ',urtemp(kplus) X end if Xc X url=ur(kminus) X if (url.gt.zero) then X if (jstop.gt.0) then X write(ipr,*)' SOLVE: inflow at top right' X write(ipr,*)' jx = ',jx,' j = ',j X idunr=1 X jplus=jx-1 X go to 3998 X else X urtemp(kminus)=zero X kminus=kminus+1 X kminus1=kminus+1 X url=ur(kminus) X end if X else X fluxm=wr(kminus) X fluxr=wr(kminus1) X slope=(fluxr-fluxm)/zsam0 X urtemp(kminus)=url+xsam*slope X end if X if (idump.gt.0) then X write(ipr,*)' SOLVE: slope @ kminus = ', slope X write(ipr,*)' updated ur value = ', ur(kminus) X end if Xc Xc 2. step interior values Xc Xc *** scalar form Xc X if (vector(1:3).eq.'no ') then Xc X url=ur(kminus) X urm=ur(kminus1) X fluxl=wr(kminus) X fluxm=wr(kminus1) X do 2052 k=kminus1,kplus1 X urr=ur(k+1) X fluxr=wr(k+1) X slr(k)=slr(k) X hs=0.5*(1.0+sign(one,url)) X fminl=hs*slr(k)+(1.0-hs)*fluxl X fmaxl=hs*fluxl+(1.0-hs)*slr(k) X hs=0.5*(1.0+sign(one,urm)) X fminm=hs*slr(k)+(1.0-hs)*fluxm X fmaxm=hs*fluxm+(1.0-hs)*slr(k) X hs=0.5*(1.0+sign(one,urr)) X fminr=hs*slr(k)+(1.0-hs)*fluxr X fmaxr=hs*fluxr+(1.0-hs)*slr(k) X fluxeol=fmaxl+fminm X fluxeor=fmaxm+fminr X slope=(fluxeor-fluxeol)/zsam0 X urtemp(k)=urm+xsam*slope X url=urm X fluxl=fluxm X urm=urr X fluxm=fluxr X2052 continue Xc Xc *** vector form Xc X else if (vector(1:3).eq.'yes') then Xc X do 2152 k=kminus1,kplus1 X hs=0.5*(1.0+sign(one,ur(k-1))) X fminl=hs*slr(k)+(1.0-hs)*wr(k-1) X fmaxl=hs*wr(k-1)+(1.0-hs)*slr(k) X hs=0.5*(1.0+sign(one,ur(k))) X fminm=hs*slr(k)+(1.0-hs)*wr(k) X fmaxm=hs*wr(k)+(1.0-hs)*slr(k) X hs=0.5*(1.0+sign(one,ur(k+1))) X fminr=hs*slr(k)+(1.0-hs)*wr(k+1) X fmaxr=hs*wr(k+1)+(1.0-hs)*slr(k) X fluxeol=fmaxl+fminm X fluxeor=fmaxm+fminr X slope=(fluxeor-fluxeol)/zsam0 X urtemp(k)=ur(k)+xsam*slope X2152 continue Xc Xc *** option screwed up Xc X else Xc X ier=90 X write(ipr,*)' SOLVE: improper vector option in 2052 loop' X return Xc X end if Xc Xc *** update offset, slowness Xc X x=x+xsam X theta1=(x-xc)/xsam0 X theta2=(xp-x)/xsam0 X do 2055 k=kminus,kplus X slr(k)=theta1*slow(jx,k) + theta2*slow(jx-1,k) X2055 continue Xc Xc *** integrate urtemp in z to get time, update ur, wr Xc X test = zero X if (jcut.eq.0) then X do 2060 k=kminus,kplus X sl=slr(k) X url=urtemp(k) X wp=(sl-url)*(sl+url) X test=min(wp,test) X wp=sqrt(abs(wp)) X wc=wr(k) X wr(k)=wp X tr(k)=tr(k)+0.5*xsam*(wc+wp) X ur(k)=url X2060 continue X else X do 2061 k=kminus,kplus X sl=slr(k) X url=min(urtemp(k),umax*sl) X url=max(url,-umax*sl) X wp=(sl-url)*(sl+url) X test=min(wp,test) X wp=sqrt(abs(wp)) X wc=wr(k) X wr(k)=wp X tr(k)=tr(k)+0.5*xsam*(wc+wp) X ur(k)=url X2061 continue X end if Xc Xc *** debug dump Xc X if (idump.gt.0) then X write(ipr,*)' SOLVE: right step ',jx X write(ipr,*)' SOLVE: internal step ',j X write(ipr,*)' xsam0, zsam0, xsam, x, xc, xp:' X write(ipr,*) xsam0, zsam0, xsam, x, xc, xp X write(ipr,*)'kminus,kplus,kminus1,kplus1:' X write(ipr,*)kminus,kplus,kminus1,kplus1 X write(ipr,*)'(ur(k),k=kminus,kplus) at offset ',x X write(ipr,*)' ' X write(ipr,4998)kminus,kplus,zsam0 X write(ipr,4999)(ur(k),k=kminus,kplus) X write(ipr,*)'(wr(k),k=kminus,kplus) at offset ',x X write(ipr,*)' ' X write(ipr,4998)kminus,kplus,zsam0 X write(ipr,4999)(wr(k),k=kminus,kplus) X write(ipr,*)'(tr(k),k=kminus,kplus) at offset ',x X write(ipr,*)' ' X write(ipr,4998)kminus,kplus,zsam0 X write(ipr,4999)(tr(k),k=kminus,kplus) X end if Xc Xc *** vertical ray condition Xc X if (test.lt.zero) then X write(ipr,*)' SOLVE: detects vertical ray at right' X write(ipr,*)' right step ',jx X write(ipr,*)' internal step ',j X write(ipr,*)' offset ',x X ier=0 X idunr=1 X jplus=jx-1 X go to 3998 X end if Xc X2100 continue Xc X write(ipr,*)' SOLVE: failure to complete right step' X write(ipr,*)' step = ',jx X idunr=1 X jplus=jx-1 X go to 3998 Xc X2800 continue Xc Xc 3. update time, gradient Xc X do 2900 k=kminus,kplus X time(jx,k)=tr(k) X2900 continue Xc Xc *** pick off new values for bottom arrays Xc X wb(jx)=wr(kplus) X ub(jx)=ur(kplus) X tb(jx)=tr(kplus) Xc X3000 continue Xc Xc *** reset count of successfully completed left steps Xc X jplus=jnext Xc X3998 continue Xc X end if Xc X3999 continue Xc Xc Xc*************** end of right step ******************** Xc* * Xc* *** continue loop by returning to down step, * Xc* or end * Xc* * Xc*************** 5000 loop **************************** Xc X if (idunl*idunb*idunr.ne.0) then X write(ipr,*)' SOLVE: done with steps' X go to 9999 X end if Xc X5000 continue Xc Xc*************** end of it all ************************ Xc X9999 continue X write(ipr,*)' SOLVE: exit' X if (idump.ge.0) then X write(ipr,*)' jminus = ',jminus X write(ipr,*)' jplus = ',jplus X write(ipr,*)' kplus = ',kplus X end if Xc X4998 format(i3,', ',i3,', ',e10.4) X4999 format(6(2x,e10.4)) Xc X return X end Xc===================== subroutine output ======================= Xc X subroutine output(ier,idump,ipr,irpt, X & nxmin,nxmax,nzmin,nzmax,jminus0,jplus0, X & kminus,kstep,xsam0,zsam0, X & time,wb,ul,jobname,jobsize, X & maxchk,maxhor,chkshot,horizons,xchk,zhor) Xc Xc *** arguments Xc X integer ier,idump,ipr,irpt, X & nxmin,nxmax,nzmin,nzmax, X & jminus0,jplus0,kminus,kstep X real*4 xsam0,zsam0 X real*4 time(nxmin:nxmax,nzmin:nzmax) X real*4 wb(nxmin:nxmax),ul(nzmin:nzmax) Xc X integer jobsize,outsize X character*20 jobname,blk20 X character*80 outname X integer chkshot,horizons X integer maxchk,maxhor X real xchk(maxchk),zhor(maxhor) Xc Xc *** internal workspace Xc X integer i,j,k,ijnk X character*1 ichar Xc Xc *** scratch unit number Xc X ijnk=8 Xc Xc *** twenty blanks Xc X blk20=' ' Xc Xc *** say hello, Gracie... Xc X write(ipr,*)' OUTPUT...' Xc Xc *** write out time field Xc X outname=blk20//blk20//blk20//blk20 X outsize=jobsize X outname(1:jobsize)=jobname(1:jobsize) X outname(jobsize+1:jobsize+4)= '.dat' X outsize= jobsize+4 X write(ipr,*)' OUTPUT: time field on ',outname(1:outsize) X call swritmat(nxmin,nxmax,nzmin,nzmax, X . jminus0,jplus0,kminus,kstep,xsam0,zsam0, X . time,outname,outsize,irpt,ier) X if (ier.ne.0) then X write(ipr,*)' Error: OUTPUT from SWRITMAT' X write(ipr,*)' ier = ',ier X go to 999 X end if Xc Xc *** write out checkshot profiles Xc X open(unit=ijnk,status='scratch') X do 300 i=1,chkshot X j=int(xchk(i)/xsam0) X if (j.lt.jminus0) go to 300 X if (j.gt.jplus0) go to 300 X rewind ijnk X write(ijnk,5000)i X rewind ijnk X read(ijnk,5001)ichar X outsize=jobsize+2 X outname(1:outsize)=jobname(1:jobsize)//'x'//ichar Xc X outname(outsize+1:outsize+4)='.dat' X outsize=outsize+4 X do 275 k=kminus,kstep X ul(k)=time(j,k) X275 continue X call swritmat(nzmin,nzmax,1,1,kminus,kstep,1,1, X . zsam0,xsam0,ul,outname,outsize,irpt, X . ier) X if (ier.ne.0) then X write(ipr,*)' OUTPUT: error on write in .MAT format' X write(ipr,*)' ier = ',ier X go to 999 X end if Xc X300 continue X do 310 i=1,horizons X k=int(zhor(i)/zsam0) X if (k.lt.kminus) go to 310 X if (k.gt.kstep) go to 310 X rewind ijnk X write(ijnk,5000)i X rewind ijnk X read(ijnk,5001)ichar X outsize=jobsize+2 X outname(1:outsize)=jobname(1:jobsize)//'z'//ichar Xc X outname(outsize+1:outsize+4)='.dat' X outsize=outsize+4 X do 280 j=jminus0,jplus0 X wb(j)=time(j,k) X280 continue X call swritmat(nxmin,nxmax,1,1,jminus0,jplus0,1,1, X . xsam0,zsam0,wb,outname,outsize,irpt, X . ier) X if (ier.ne.0) then X write(ipr,*)' OUTPUT: error on write in .MAT format' X write(ipr,*)' ier = ',ier X go to 999 X end if Xc X310 continue X close(unit=ijnk) Xc X999 continue Xc X5000 format(i1) X5001 format(a1) X return X end Xc Xc======================================================== Xc X subroutine filnam(name,size,line,ipin,ipout,ierr) X integer j,size,ipin,ipout,ierr X character*80 name,line,fname X character blank X character*10 blank10 X data blank10/' '/ X data blank /' '/ Xc Xc prompt for name Xc X write(ipout,*)line Xc Xc read file name Xc X fname(1:10)= blank10 X fname(11:20)= blank10 X fname(21:30)= blank10 X fname(31:40)= blank10 X fname(41:50)= blank10 X fname(51:60)= blank10 X fname(61:70)= blank10 X fname(71:80)= blank10 Xc X read(ipin,1111) fname X1111 format(a80) Xc Xc get size of file name Xc X name(1:10)= blank10 X name(11:20)= blank10 X name(21:30)= blank10 X name(31:40)= blank10 X name(41:50)= blank10 X name(51:60)= blank10 X name(61:70)= blank10 X name(71:80)= blank10 Xc X j = index(fname,blank) X if( j .gt. 1 ) then X size = j-1 X name(1:size) = fname(1:size) X else X write(ipout,*) ' Bad file name' X ierr=98 X return X endif X write(ipout,*) 'file name = ', fname(1:size) Xc X return X end END_OF_FILE if test 68102 -ne `wc -c <'dno/dno.f'`; then echo shar: \"'dno/dno.f'\" unpacked with wrong size! fi # end of 'dno/dno.f' fi if test -f 'dno/init.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'dno/init.f'\" else echo shar: Extracting \"'dno/init.f'\" \(9357 characters\) sed "s/^X//" >'dno/init.f' <<'END_OF_FILE' Xc***********************************************************************c Xc c Xc Copyright (C) 1989, 1990 William W. Symes c Xc c Xc This file is part of the Down-n-Out Traveltime Computation c Xc Project ("DNO"). c Xc c Xc The DNO codes are distributed in the hope that they will be c Xc useful, but WITHOUT ANY WARRANTY, without even the implied c Xc warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. c Xc c Xc You are free to modify this file provided that you leave this c Xc header intact. I ask that you not redistribute this file c Xc without first notifying me. This isn't because I'm trying c Xc to be nasty here; it's simply that I would like to keep track c Xc of who has the program. c Xc c Xc These codes were developed by William W. Symes with support from c Xc the National Science Foundation, the Office of Naval Research, c Xc and the Air Force Office of Scientific Research. Thanks are due c Xc to Jim Carazzone, Dave Hale, and Jos van Trier for helpful c Xc suggestions. c Xc c Xc William W. Symes c Xc Department of Mathematical Sciences c Xc Rice University c Xc Houston, Texas 77005 c Xc USA c Xc c Xc e-mail: symes@rice.edu c Xc c Xc March 1990 c Xc c Xc***********************************************************************c Xc Xc produces initial time file for DOWN N' OUT - use and Xc structure self-explanatory, so no need for comments, Xc RIGHT? Xc X integer nxmin,nxmax,nzmin,nzmax X parameter (nxmin=-100,nxmax=100,nzmin=0,nzmax=200) X integer i,j,k,ipout,hdimx,ipr,ipin,ndat, X & im,ip,jm,jp,ier X real x,z,xsam,zsam,xm,xp,zm,zp,xs,zs, X & test,tol,slow X real time(nxmin:nxmax,nzmin:nzmax) X character*80 name Xc Xc *** spurious defaults Xc X real cfl X integer dzswh,maxstep,jstop,jcut X character*3 vector Xc Xc Xc *** set unit numbers and defaults Xc X name(1:11)='initdef.dat' X call setdef(ipin,ipr,ipout,cfl,dzswh,maxstep,jstop, X . jcut,vector,ier,name,11) X if (ier.ne.0) then X write(6,*)' Error: SETDEF, ier = ', ier X go to 999 X end if Xc X write(ipr,1000) X1000 format(' Initial data generator for DNO',//, X & ' Generates point source time field in a',/, X & ' prescribed box, which is to be continued',/, X & ' into a larger domain by DNO',/, X & ' NOTE: it is assumed that the velocity or slow-',/, X & ' ness field is HOMOGENEOUS throughout the box ',/, X & ' and in a neighboring layer of gridpoints.',/, X & ' If this provision is violated, DNO will will',/, X & ' encounter an error condition.',/) Xc X12 continue X write(ipr,*)' Enter x coordinates for ends of data box' X read(ipin,*)xm,xp X13 continue X write(ipr,*)' Enter x step ' X read(ipin,*)xsam X tol=max(abs(xm),abs(xp)) X test=xsam+tol X if (test.le.tol) then X write(ipr,*)' xsam < or = 0; try again' X go to 13 X end if X im=int(xm/xsam) X if (im.lt.nxmin) then X write(ipr,*)' too far left' X go to 12 X end if X ip=int(xp/xsam) X if (ip.gt.nxmax) then X write(ipr,*)' too far right' X go to 12 X end if X14 continue X write(ipr,*)' Enter z coordinates for ends of data box' X read(ipin,*)zm,zp X15 continue X write(ipr,*)' Enter z step' X read(ipin,*)zsam X tol=max(abs(zm),abs(zp)) X test=zsam+tol X if (test.le.tol) then X write(ipr,*)' zsam < or = 0; try again' X go to 15 X end if X jm=int(zm/zsam) X if (jm.lt.nzmin) then X write(ipr,*)' too high' X go to 14 X end if X jp=int(zp/zsam) X if (jp.gt.nzmax) then X write(ipr,*)' too low' X go to 14 X end if X write(ipr,*)' Enter x and z coordinates for source point' X read(ipin,*)xs,zs X write(ipr,*)' Enter slowness in box' X read(ipin,*)slow X do 40 j=jm,jp X do 30 i=im,ip X x=i*xsam-xs X z=j*zsam-zs X time(i,j)=slow*sqrt(x*x+z*z) X30 continue X40 continue Xc X name(1:9)='time0.dat' Xc X call swritmat(nxmin,nxmax,nzmin,nzmax, X . im,ip,jm,jp,xsam,zsam, X . time,name,9,ipout,ier) X if (ier.ne.0) then X write(ipr,*)' INIT: error on swritmat, ier = ',ier X end if Xc X999 continue X stop X end X subroutine setdef(ipin,ipr,ips,cfl,dzswh,maxstep,jstop, X . jcut,vector,ier,name,size) Xc------------------------------------------------------------- Xc Xc reads unit numbers and flags from file "dnodef.dat" Xc Xc------------------------------------------------------------- Xc X integer ipin,ipr,ips,dzswh,maxstep,jstop,jcut,ier,size X real cfl X character*3 vector X character*80 name X open(unit=2,file=name(1:size),status='old', X . iostat=ier,err=999) X read(2,1000)ipin X read(2,1000)ipr X read(2,1000)ips X read(2,2000)cfl X read(2,1000)dzswh X read(2,1000)maxstep X read(2,1000)jstop X read(2,1000)jcut X read(2,3000)vector X close(unit=2) X 999 continue X 1000 format(i2) X 2000 format(e10.4) X 3000 format(a3) X return X end Xc Xc=============================================================== Xc X subroutine swritmat(ld1min, ld1max, ld2min, ld2max, X . n1min, n1max, n2min, n2max, d1, d2, X . a, nameout, sizeout, iounit, X . ier) Xc Xc SWRITMAT - A Fortran subroutine that writes a single-precision Xc matrix of uniformly spaced samples of a real function Xc of one or two real variables to an ASCII flat file. Xc Xc Xc Description of arguments: Xc Xc ld1min| Xc ld1max|______> declared index limits EXACTLY as in Xc ld2min| calling program Xc ld2max| Xc Xc n1min| Xc n1max|_______> ranges of indices to be written out; Xc n2min| stored as array words 1,2,3,4 Xc n2max| Xc Xc d1|__________> sample increments (words 5,6) Xc d2| Xc Xc a - output array of real*4 numbers, Xc regarded internally as a 2-d array (words 7,...) Xc Xc nameout - name for output file Xc sizeout - length of name Xc Xc iounit - logical unit number of output file Xc ier - error flag Xc Xc ASCII Flat File: Xc Xc All entries are stored as single-precision reals in Fortran Xc fomatted i/o. The first six entries are n1min, n1max, n2min, n2max, Xc d1 and d2 in the notation explained above. Xc Xc This routine was written by W. W. Symes, 2-25-90. Xc Standard disclaimer is hereby declared. Xc Xc Xc================================================================ Xc Xc *** arguments: Xc Xc integer arguments: X integer*4 ld1min,ld1max,ld2min,ld2max, X . n1min,n1max,n2min,n2max, X . sizeout, iounit, ier Xc Xc real arguments: X real*4 d1,d2,a(ld1min:ld1max,ld2min:ld2max) Xc Xc character arguments: X character*80 nameout Xc Xc loop counters, limit X integer i, k, nwrit Xc Xc workspace for type conversion X real*4 xs1, xs2, xs3, xs4 Xc Xc padding for ascii records X real*4 zeroes(1:5) Xc Xc *** open file for output: Xc X open(unit=iounit,file=nameout(1:sizeout), X . status='unknown', X . access='sequential',form='formatted', X . iostat=ier,err=992) X if (ier.ne.0) go to 992 Xc Xc type conversions, padding Xc X xs1=real(n1min) X xs2=real(n1max) X xs3=real(n2min) X xs4=real(n2max) X do 800 i=1,5 X zeroes(i)=0.0e+00 X800 continue Xc Xc *** write stuff Xc X nwrit=6+(n1max-n1min+1)*(n2max-n2min+1) X nwrit=5*((nwrit/5)+1)-nwrit X write(iounit,1000)xs1,xs2,xs3,xs4,d1,d2, X . ((a(i,k),i=n1min,n1max),k=n2min,n2max), X . (zeroes(i),i=1,nwrit) X1000 format(5(2x,e12.5)) Xc Xc good write X close(unit=iounit) X return Xc Xc error during file opening X992 ier=8 X close(unit=iounit) X return Xc X end END_OF_FILE if test 9357 -ne `wc -c <'dno/init.f'`; then echo shar: \"'dno/init.f'\" unpacked with wrong size! fi # end of 'dno/init.f' fi if test -f 'dno/slow.f' -a "${1}" != "-c" ; then echo shar: Will not clobber existing file \"'dno/slow.f'\" else echo shar: Extracting \"'dno/slow.f'\" \(13945 characters\) sed "s/^X//" >'dno/slow.f' <<'END_OF_FILE' Xc***********************************************************************c Xc c Xc Copyright (C) 1989, 1990 William W. Symes c Xc c Xc This file is part of the Down-n-Out Traveltime Computation c Xc Project ("DNO"). c Xc c Xc The DNO codes are distributed in the hope that they will be c Xc useful, but WITHOUT ANY WARRANTY, without even the implied c Xc warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. c Xc c Xc You are free to modify this file provided that you leave this c Xc header intact. I ask that you not redistribute this file c Xc without first notifying me. This isn't because I'm trying c Xc to be nasty here; it's simply that I would like to keep track c Xc of who has the program. c Xc c Xc These codes were developed by William W. Symes with support from c Xc the National Science Foundation, the Office of Naval Research, c Xc and the Air Force Office of Scientific Research. Thanks are due c Xc to Jim Carazzone, Dave Hale, and Jos van Trier for helpful c Xc suggestions. c Xc c Xc William W. Symes c Xc Department of Mathematical Sciences c Xc Rice University c Xc Houston, Texas 77005 c Xc USA c Xc c Xc e-mail: symes@rice.edu c Xc c Xc March 1990 c Xc c Xc***********************************************************************c Xc Xc Xc slowness field generator for DOWNNOUT Xc Xc inputs: size and grid spacing of slowness array Xc type of example: Xc 1. constant (== 1.0) Xc 2. two half-spaces with linear transition Xc 3. cosine lens Xc parameters for each type Xc Xc variable definitions Xc mx,mz max horizontal array limit Xc s(-mx:mx) Xc nxmin, horizontal array limits Xc nxmax Xc nzmin, vertical array limits Xc nzmax Xc xsam horizontal sample rate Xc zsam vertical sample rate Xc ips output unit number Xc s slowness array Xc name output file name Xc Xc plus a bunch of workspace Xc Xc============================================================ Xc Xc *** declarations Xc X integer mx X parameter (mx=400) X integer nxmin,nxmax,nzmin,nzmax,i,j,ipr,ips,ipin X real xsam,zsam,xmin,xmax,zmin,zmax,av,buf,s0, X & slow,s(-mx:mx,-mx:mx),prof(-mx:mx),temp(-mx:mx) X integer itype,idump,ism,n,nsm,ier X real x,z,xleft,xright,ztop,zbot,scum,test,tol X character yn X character*80 name Xc Xc *** spurious defaults Xc X real cfl X integer dzswh,maxstep,jstop,jcut X character*3 vector Xc Xc *** set unit numbers and defaults Xc X ier=0 X name(1:11)='slowdef.dat' X call setdef(ipin,ipr,ips,cfl,dzswh,maxstep,jstop, X . jcut,vector,ier,name,11) X if (ier.ne.0) go to 999 Xc Xc *** dump switch for slowness function Xc X idump=0 Xc Xc *** prompt for discretization information Xc X write(ipr,3287) X3287 format(///,' SLOW - slowness in a box',//) Xc X12 continue X write(ipr,*)' Enter x-limits of box' X read(ipin,*)xmin,xmax X14 continue X write(ipr,*)' Enter x sample interval' X read(ipin,*)xsam X tol=max(abs(xmin),abs(xmax)) X test=tol+xsam X if (test.le.tol) then X write(ipr,*)' xsam < or = 0, fathead; try again' X go to 14 X end if X nxmin=int(xmin/xsam) X nxmax=int(xmax/xsam) X if (nxmin.lt.-mx) then X write(ipr,*)' too far left' X go to 12 X end if X if (nxmax.gt.mx) then X write(ipr,*)' too far right' X go to 12 X end if X15 continue X write(ipr,*)' Enter z-limits of box' X read(ipin,*)zmin,zmax X16 continue X write(ipr,*)' Enter z sample interval' X read(ipin,*)zsam X tol=max(abs(zmin),abs(zmax)) X test=zsam+tol X if (test.le.tol) then X write(ipr,*)' zsam < or = 0, fathead; try again' X go to 16 X end if X nzmin=int(zmin/zsam) X nzmax=int(zmax/zsam) X if (nzmin.lt.-mx) then X write(ipr,*)' too far up' X go to 15 X end if X if (nzmax.gt.mx) then X write(ipr,*)' too far down' X go to 15 X end if X write(ipr,*)' x-limits:' X write(ipr,*)' ',nxmin,nxmax X write(ipr,*)' z-limits:' X write(ipr,*)' ',nzmin,nzmax X write(ipr,*)' proceed? (y/n)' X read(ipin,1111)yn X1111 format(a1) X if (yn.ne.'y') go to 12 Xc Xc *** choose example type and parameters Xc X18 continue X write(ipr,*)' Choose:' X write(ipr,*)' 1. constant slowness (== 1.0)' X write(ipr,*)' 2. two layers of relative slowness 1, s' X write(ipr,*)' connected by linear transition' X write(ipr,*)' 3. cosine lens with specified relative' X write(ipr,*)' slowness at center' X write(ipr,*)' NOTE: all choices of slowness are relative' X write(ipr,*)' to slowness in vicinity of source, for the' X write(ipr,*)' value of which you will be prompted!' X read(ipin,*)itype X if ((itype.lt.1).or.(itype.gt.3)) then X write(ipr,*)' invalid choice, fathead!' X go to 18 X end if X if (itype.eq.2) then X write(ipr,*)' enter top, bottom layer depths' X read(ipin,*)ztop,zbot X write(ipr,*)' enter relative slowness in bottom layer' X read(ipin,*)scum X x=0.0 X do 34 j=nzmin,nzmax X z=j*zsam X prof(j)=slow(x,z,xleft,xright,ztop,zbot,scum, X & itype,idump,ipr) X34 continue X write(ipr,*)' enter width of moving average, # pts' X read(ipin,*)ism X av=float(2*ism+1) X av=1.0/av X write(ipr,*)' enter number of m. a. iterations' X read(ipin,*)nsm X do 39 n=1,nsm X do 38 i=nzmin,nzmax X temp(i)=prof(i) X38 continue X do 36 i=nzmin+ism,nzmax-ism X buf =0. X do 35 j=i-ism,i+ism X buf=buf+prof(j) X35 continue X temp(i)=buf*av X36 continue X do 37 i=nzmin,nzmax X prof(i)=temp(i) X37 continue X39 continue X else if (itype.eq.3) then X write(ipr,*)' enter depths for top and bottom of lens' X read(ipin,*)ztop,zbot X write(ipr,*)' enter x-values for left and right ' X write(ipr,*)' sides of lens' X read(ipin,*)xleft,xright X write(ipr,*)' enter relative slowness at center of lens' X read(ipin,*)scum X end if Xc Xc *** reference slowness Xc X write(ipr,*)' enter (reference) slowness near source ' X read(ipin,*)s0 Xc Xc *** compute slowness Xc X do 200 j=nzmin,nzmax X z=j*zsam X do 100 i=nxmin,nxmax X x=i*xsam X if (itype.ne.2) then X s(i,j)=s0*slow(x,z,xleft,xright,ztop,zbot,scum, X & itype,idump,ipr) X else X s(i,j)=s0*prof(j) X end if X100 continue X200 continue Xc Xc *** output Xc X name(1:8)='slow.dat' Xc X ier=0 X call swritmat(-mx, mx, -mx, mx, X . nxmin, nxmax, nzmin, nzmax, xsam, zsam, X . s, name, 8, ips, ier) Xc X999 continue X stop X end Xc Xc==================== function slow ====================== Xc X real function slow(x,z,left,right,top,bottom, X & smax,itype,idbg2,ipr) Xc Xc this function computes the slowness field Xc Xc no error-checking - presumes good input Xc X real x, z X integer itype,idbg2,ipr X real top,bottom,left,right,smax Xc X real pi, xn, zn Xc X pi=4.0e+00*atan(1.0e+00) Xc X if (itype.eq.1) then X slow = 1.0e+00 X if (idbg2.gt.0) then X write(ipr,*)' RHS: itype = ', itype,' case = 1' X end if X else if (itype.eq.2) then X if (idbg2.gt.0) then X write(ipr,*)' RHS: itype = ', itype,' case = 2' X end if X if (z.lt.top) then X slow =1.0e+00 X else if (z.lt.bottom) then X slow = 1. + (z-top)*(smax-1.0)/(bottom-top) X else X slow = smax X end if X else if (itype.eq.3) then X if (idbg2.gt.0) then X write(ipr,*)' RHS: itype = ', itype,' case = 3' X end if X if (z.lt.top) then X slow = 1.0e+00 X else if (z.lt.bottom) then X if ((x.lt.left).or.(x.gt.right)) then X slow = 1.0e+00 X else X xn = 2.0*pi*(x-left)/(right-left) X zn = 2.0*pi*(z-top)/(bottom-top) X slow = 1.0 + 0.25*(smax-1.0)*((1.0 - cos(zn)) X & *(1.0 - cos(xn))) X end if X else X slow = 1.0e+00 X end if X end if X if (idbg2.gt.0) then X write(ipr,*)' RHS:' X write(ipr,*)' left, right, top, bottom, pi, smax' X write(ipr,*) left, right, top, bottom, pi, smax X write(ipr,*)' z, x, zn, xn, slow' X write(ipr,*) z, x, zn, xn, slow X end if Xc X return X end X X subroutine setdef(ipin,ipr,ips,cfl,dzswh,maxstep,jstop, X . jcut,vector,ier,name,size) Xc------------------------------------------------------------- Xc Xc reads unit numbers and flags from file "dnodef.dat" Xc Xc------------------------------------------------------------- Xc X integer ipin,ipr,ips,dzswh,maxstep,jstop,jcut,ier,size X real cfl X character*3 vector X character*80 name X open(unit=2,file=name(1:size),status='old', X . iostat=ier,err=999) X read(2,*)ipin X read(2,*)ipr X read(2,*)ips X read(2,*)cfl X read(2,*)dzswh X read(2,*)maxstep X read(2,*)jstop X read(2,*)jcut X read(2,3000)vector X close(unit=2) X999 continue X3000 format(a3) X return X end Xc Xc=============================================================== Xc X subroutine swritmat(ld1min, ld1max, ld2min, ld2max, X . n1min, n1max, n2min, n2max, d1, d2, X . a, nameout, sizeout, iounit, X . ier) Xc Xc SWRITMAT - A Fortran subroutine that writes a single-precision Xc matrix of uniformly spaced samples of a real function Xc of one or two real variables to an ASCII flat file. Xc Xc Xc Description of arguments: Xc Xc ld1min| Xc ld1max|______> declared index limits EXACTLY as in Xc ld2min| calling program Xc ld2max| Xc Xc n1min| Xc n1max|_______> ranges of indices to be written out; Xc n2min| stored as array words 1,2,3,4 Xc n2max| Xc Xc d1|__________> sample increments (words 5,6) Xc d2| Xc Xc a - output array of real*4 numbers, Xc regarded internally as a 2-d array (words 7,...) Xc Xc nameout - name for output file Xc sizeout - length of name Xc Xc iounit - logical unit number of output file Xc ier - error flag Xc Xc ASCII Flat File: Xc Xc All entries are stored as single-precision reals in Fortran Xc fomatted i/o. The first six entries are n1min, n1max, n2min, n2max, Xc d1 and d2 in the notation explained above. Xc Xc This routine was written by W. W. Symes, 2-25-90. Xc Standard disclaimer is hereby declared. Xc Xc Xc================================================================ Xc Xc *** arguments: Xc Xc integer arguments: X integer*4 ld1min,ld1max,ld2min,ld2max, X . n1min,n1max,n2min,n2max, X . sizeout, iounit, ier Xc Xc real arguments: X real*4 d1,d2,a(ld1min:ld1max,ld2min:ld2max) Xc Xc character arguments: X character*80 nameout Xc Xc loop counters, limit X integer i, k, nwrit Xc Xc workspace for type conversion X real*4 xs1, xs2, xs3, xs4 Xc Xc padding for ascii records X real*4 zeroes(1:5) Xc Xc *** open file for output: Xc X open(unit=iounit,file=nameout(1:sizeout), X . status='unknown', X . access='sequential',form='formatted', X . iostat=ier,err=992) X if (ier.ne.0) go to 992 Xc Xc type conversions, padding Xc X xs1=real(n1min) X xs2=real(n1max) X xs3=real(n2min) X xs4=real(n2max) X do 800 i=1,5 X zeroes(i)=0.0e+00 X800 continue Xc Xc *** write stuff Xc X nwrit=6+(n1max-n1min+1)*(n2max-n2min+1) X nwrit=5*((nwrit/5)+1)-nwrit X write(iounit,1000)xs1,xs2,xs3,xs4,d1,d2, X . ((a(i,k),i=n1min,n1max),k=n2min,n2max), X . (zeroes(i),i=1,nwrit) X1000 format(5(2x,e12.5)) Xc Xc good write X close(unit=iounit) X return Xc Xc error during file opening X992 ier=8 X close(unit=iounit) X return Xc X end END_OF_FILE if test 13945 -ne `wc -c <'dno/slow.f'`; then echo shar: \"'dno/slow.f'\" unpacked with wrong size! fi # end of 'dno/slow.f' fi echo shar: End of archive 1 \(of 1\). cp /dev/null ark1isdone MISSING="" for I in 1 ; do if test ! -f ark${I}isdone ; then MISSING="${MISSING} ${I}" fi done if test "${MISSING}" = "" ; then echo You have the archive. rm -f ark[1-9]isdone else echo You still need to unpack the following archives: echo " " ${MISSING} fi ## End of shell archive. exit 0