program shootsr

# Code is shoot copied from tt-resolution March 8, 2000 but imput 
# modified to read output files from dbcrosscr - RMA

# Traces rays through 1d (background) vel model (iaspi91 tables used).
# Tried to get ray bending to work aswell ...not sure of outcome.
# Input: station and event locations, ray type. 
# Output: ray segment lengths associated with nodes -> the matrix 
#         also relocation derivatives 

# Original code: Costas, Mike Deal
# development:   Richard M Allen	Feb 16, 96
# development:   Richard M Allen	Aug 99
# development:   Richard M Allen	May 2000 - read ray input from one file

# Input files:
#  
# file='names'    # names of other files: 1st is input, last 3 output
#   read(88,*)dfile1		  # all the ray param stuff
#                                    ...format output by dbcrosscr[bp2]
#   read(88,*)matfile             # matrix file (unformatted)
#   read(88,*)resfile             # residual file (unformatted)
#   
#
# shoot.params ...the input file, format:
#
#  888 format(i7,1x,f15.5,1x,a5,1x,a8,1x,f9.4,1x,f9.4,1x,f9.4,1x,f9.4,1x,
#		f9.4,1x,f9.4,1x,f8.3,1x,f7.2,1x,f7.2,1x,f8.5,1x,f8.5,1x,f8.5)
#
#  read(31,888,end=999) evday,evtime,phasej,stat1,slatd,slond,stelev,elatd,
#				elond,deps,delta,rdum,iaz,kres,stdev,rav
#
#      INPUTS for each ray to be traced:
#      evday,			event yearjdy for ID only
#      evtime,			event epoch time for ID only
#      phasej,			phase eg P
#      stat1,slatd,slond,stelev,	station name, lat, lon, elev
#      elatd,elond,deps,		event lat lon depth
#      delta,rdum,iaz,		epi dist (deg) <NOTUSED> ev to st azimuth
#      kres,stdev,rav		tt resids, standard deviations, x-corr coef
#
#  
# model.cube - defines grid, as lincomb software
#   read(1,*) ola,olo,oazi,nz,ny,nx,zmin,zmax,ymin,ymax,xmin,xmax,dz,dy,dx
#  
#
# Running SHOOT:
#   get messages to screen about the sucess of each ray tracing
#  
#
# Output:
#  
# matrixfile - contains all matrix info.
#   per ray:
#     write(22) flag,ndvel,ievnt,istind,phasorg,onset,iflgorg
#     write(22) (idvel(j),j=1,ndvel)
#     write(22) (dvel(j),j=1,ndvel)
#     write(22) dtx,dty,dtz
#         ndvel - number of non-zero matrix elements
#         idvel - indices of non-zero matrix elements
#          dvel - non-zero matrix elements
#                  -> km of ray path associated with each node
#         dtx,dty,dtz - hypo center derivatives ie dt/dx, dt/dy ..
#                 - for relocation of evets - units: s/km
#  
# resdiual file - simply writes out travel-time resid as read in
#   pre ray:
#     write(23) robres,flag,ndvel,ievnt,istind,phasorg,onset,iflgorg
#         robres - residual
#         other sutff just for info
#  
# data.stdev file - simply writes out the standard deviation
#         ie certainty of the data for use in covariance matrix
#   per ray:
#     write(24,*)stdev        # stdev of the data
#  
# 


include 'defines3'   	# now in file above
include 'comm.ray3'	# added below

common /der/dvel(NVEL),idvel(NVEL),ndvel
common /gr/n1,n2,n3,x1,x2,y1,y2,z1,z2,dx,dy,dz,tc,phc
common/gri/x1i,x2i,y1i,y2i,z1i,z2i
common/ang/cltc,sltc,clnc,slnc
common/flags/nflag,iflag
common/names/dfile1,dfile2,dfile4,dfile5,matfile,resfile,velmod
#dimension  r(NRAY), rdisc(12), t(NRAY),lr(NRAY)
dimension  r(NRAY),  t(NRAY),lr(NRAY)
dimension mporg(11,65),msorg(11,65),mp(11,65),ms(11,65),kfr(65),klr(65),jsray (65)
dimension input(4320,2160)
dimension iflag(1000)
character*8 statst(200),stat1,kstnm
character*8 name(65)
double precision cltc,sltc,clnc,slnc
real mddelta,dpdd,slat,slon,lat,lon,elat,elon,delta
real robres,robtime,latc,lonc
real x1,x2,y1,y2,z1,z2,dx,dy,dz
real del,baz,trtime,time,elcr,stel,azis,offset,iaz,iazo,kres
character*4 dum4,dum5
character*1 onset
character*8 phasej,phasorg
#character*15 evnt,event(500)
character*30 datafile,statelev,matfile,resfile,velmod
character*30 dfile1,dfile2,dfile4,dfile5,dfile78
integer err,reas,ditim
integer ls,n1,n2,n3,flag,boxx
integer*2  input
logical debug

double precision evtime,event(500)
real eventla(500),eventlo(500),statla(500),statlo(500),statel(500)



# These are the different ray types recognized
# Currently I've turned off all S rays ...RMA: turned back on Aug 2000!
DATA name/'P','PKKP','PKP','PKPPKS','PKPPKP','PKS','PP','PPP',
      'PPS','PS','PSS','PcP','PcPPKP','PcS','PcSPKP','S','SKKS',
      'SKKKS','SKP','SKS','SKSP','SP','SPP','SS','SSP','SSS','ScP',
      'ScS','ScSP','ScSPKP','pP','pPKP','pPKPPKP','sP','sPKP',
      'sPKPPKP','sS','PKIKP','PKKS','(PKP)3','PKPSKS','SKKP','SKPPKP',
      'SKSSKS','pPP','pPcP','sPP','sPcP','sSS','sScS','ScSScS',
      'sScSScS','PKiKP','PKJKP','PPPP','PPPPP','PPPS','PPPPS',
      'PPPSS','SSSS','SSSSS','pPKiKP','pwP',2*'' /

# Ray number 62 is pwP  same multiplicities as pP

DATA kfr/15*3,15*4,3*1,4*2,4*3,3*4,1,1,4*2,4,2,7*3,2*4,1,1,2*0/

data klr/1,1,1,2,1,2,1,1,2,2,  2,1,1,2,1,2,2,2,1,2,
       1,1,1,2,1,2,1,2,1,1,  1,1,1,1,1,1,2,1,2,1,
       2,1,1,2,1,1,1,1,2,2,  2,2,1,1,1,1,  2,2,2,2,2,1,1, 2*0/

data nrays/63/

#       deal: changed MS MP for 11 layer model
#	If you want a background model with more/fewer than 11 layers then
#	you will have to redo the multiplicities accordingly
#       These are the mulitplicities for all rays.  I trace using symmetry
DATA ms/33*0,9*1,13*0,9*1,24*0,9*2,2*0,9*2,2*0,9*4,24*0,9*1,
      2*0,9*1,2*0,9*2,2*0,9*2,2*0,9*2,2*0,9*1,2*0,9*2,2*0,9*2,2*0,
      9*2,2*0,9*2,2*0,9*4,2*0,9*4,2*0,9*6,2*0,9*1,2*0,9*2,2*0,9*2,
      2*0,9*2,2*0,66*0,9*2,2*0,11*0,9*1,13*0,9*2,2*0,9*1,2*0,9*1,2*0,
      9*4,2*0,44*0,9*4,2*0,9*2,2*0,9*4,2*0,9*4,2*0,11*0,10*0,2,
      22*0, 9*2,2*0, 9*2,2*0, 9*4,2*0, 9*8,2*0, 9*10,2*0, 44*0/

DATA mp/9*2,2*0,9*2,4,0,9*2,2,0,9*3,4,0,10*4,0,9*1,2,0,9*4,2*0,
      9*6,2*0,9*4,2*0,9*2,2*0,9*2,2*0,9*2,2*0,9*4,2,0,9*1,2*0,9*3,2,0,
      11*0,9*0,4,0,9*0,6,0,9*1,2,0,9*0,2,0,10*2,0,9*2,2*0,9*4,2*0,9*0,
      2*0,9*2,2*0,11*0,9*1,2*0,11*0,9*1,2*0,10*2,0,9*2,2*0,10*2,0,10*4,0,
      9*2,2*0,10*2,0,10*4,0,11*0,11*2,9*1,4,0,9*6,2*0,9*2,4,0,9*1,4,0,9*3,
      4,0,
      9*0,4,0,9*4,2*0,9*2,2*0,9*4,2*0,9*2,2*0,44*0,10*2,0,10*2,0,
      9*8,0,0,9*10,0,0, 9*6,2*0, 9*8,2*0, 9*6,2*0,11*0,11*0, 10*2,0,9*2,2*0,22*0/

DATA msorg/33*0,9*1,13*0,9*1,24*0,9*2,2*0,9*2,2*0,9*4,24*0,9*1,
	2*0,9*1,2*0,9*2,2*0,9*2,2*0,9*2,2*0,9*1,2*0,9*2,2*0,9*2,2*0,
	9*2,2*0,9*2,2*0,9*4,2*0,9*4,2*0,9*6,2*0,9*1,2*0,9*2,2*0,9*2,
        2*0,9*2,2*0,66*0,9*2,2*0,11*0,9*1,13*0,9*2,2*0,9*1,2*0,9*1,2*0,
        9*4,2*0,44*0,9*4,2*0,9*2,2*0,9*4,2*0,9*4,2*0,11*0,10*0,2,
        22*0, 9*2,2*0, 9*2,2*0, 9*4,2*0, 9*8,2*0, 9*10,2*0, 44*0/

DATA mporg/9*2,2*0,9*2,4,0,9*2,2,0,9*3,4,0,10*4,0,9*1,2,0,9*4,2*0,
      9*6,2*0,9*4,2*0,9*2,2*0,9*2,2*0,9*2,2*0,9*4,2,0,9*1,2*0,9*3,2,0,
      11*0,9*0,4,0,9*0,6,0,9*1,2,0,9*0,2,0,10*2,0,9*2,2*0,9*4,2*0,9*0,
      2*0,9*2,2*0,11*0,9*1,2*0,11*0,9*1,2*0,10*2,0,9*2,2*0,10*2,0,10*4,0,
      9*2,2*0,10*2,0,10*4,0,11*0,11*2,9*1,4,0,9*6,2*0,9*2,4,0,9*1,4,0,9*3,
      4,0,
      9*0,4,0,9*4,2*0,9*2,2*0,9*4,2*0,9*2,2*0,44*0,10*2,0,10*2,0,
      9*8,0,0,9*10,0,0, 9*6,2*0, 9*8,2*0, 9*6,2*0,11*0,11*0, 10*2,0,9*2,2*0,22*0/


debug=.false. 

qwe=3.1415927/180.

input(1,1)=0
7 continue
err=0
numray=0

## read flag file
do i=1,1000 {iflag(i)=0}	# set all flags to zero
iflag(1) = 1			# indicates this is Richards format etc
nflag=1
				# if iflag(2)=1 then ray bending attempted in genxyz
# READ FROM DATA FILE
  open(88,file='names')
  read(88,*)dfile1
  read(88,*)matfile		# matrix file (unformatted)
  read(88,*)resfile		# residual file (unformatted)
  close(88)
  open(31,file=dfile1,status='old')
  dfile78='used.'//dfile1	# for output of used parameters
  open(78,file=dfile78)

open(22,file=matfile,form='unformatted')
open(23,file=resfile,form='unformatted')
open(24,file='data.stdev')
open(25,file='data.stlalo')
open(68,file='data.phslowangles')

221     format(i5,i14,2i5,i3,i8,a4,2i5,i5,i3,a1, i9,a8,i5,i6,i4,i5,i2,i1,4i5)
222 	format(a8,2x,i3,2x,f8.2,2x,f12.6,2x,f12.6,2x,f12.6,2x,f12.6)
224 	format(a8,2x,f10.4,2x,f10.4,3x,i3)
225 	format(i5,1x,f7.1,2x,a4)
230	format(3i5)
231	format(3f12.3)

## read in model cube - oazi not used ...assumed x points E
open(1,file='model.cube',status='old')
rewind 1
#read(1,*) ola,olo,oazi,nz,ny,nx,zmin,zmax,ymin,ymax,xmin,xmax,dz,dy,dx
read(1,*) latc,lonc,rdum,n3,n2,n1,z1,z2,y1,y2,x1,x2,dz,dy,dx
close(1)
write(6,*) 'read model.cube'

x1i=x1+5*dx
x2i=x2-5*dx
y1i=y1+5*dy
y2i=y2-5*dy
z1i=z1+5*dz
z2i=z2
# inner box is 5 grid units smaller than model


#Costas: convert latc & lonc to geocentric colatitude and longitude
#and store them in tc & phc. Also calculate some constants
latc=latc/RAD
latc=atan(.993277*tan(latc))
lonc=lonc/RAD
cltc=dcos(dble(latc))
sltc=dsin(dble(latc))
clnc=dcos(dble(lonc))
slnc=dsin(dble(lonc))
tc=1.57079632679-latc
# tc co-lat of center of model
phc=lonc

# I write to files fort.97-99 for rays kept and not kept, may wiwsh to change sergei
write(99,*)'These are phases not found in tables or not dealt with.'
write(99,*)'                          delta,    ps,        reas'
write(98,*)'Compare and check ray parameters from table values. NB ray accepted'
write(97,*)'These are rays with poor agreement to robs ps'
write(96,*)'These are rays with poor agreement to robs delta' 


qwe=3.1415927/180.

ik=0
irayno=0

### TOP OF LOOP OVER ALL RAYS #########################################
223 continue

# call msmp to reset ms and mp vectors (multiplicities) 
call msmp(mp,ms,mporg,msorg)

## check no exceeding array length
if (ik>=1500000) { go to 999  }

## read in line of data from each file NB everything in radians by end
## new input file format ...output from dbcrosscr: *.sr.* files

#888 format(i7,1x,f16.5,1x,a5,1x,a8,1x,f9.4,1x,f9.4,1x,f9.4,1x,f9.4,1x,f9.4,1x,f9.4,1x,f8.3,1x,f7.2,1x,f7.2,1x,f8.5,1x,f8.5,1x,f8.5) #mei

## XM: changed the 888 format, now consistent with the output format 888 from dbcrosscrbp.r
## !!! very IMPORTANT for pass parameters right

888 format(i7,1x,f16.5,1x,a5,1x,a8,1x,f9.4,1x,f9.4,1x,f9.4,1x,f9.4,1x,f9.4,1x,f9.4,1x,f8.3,1x,f7.2,1x,f7.2,1x,f8.5,1x,f8.5,1x,f8.5)
read(31,888,end=999) evday,evtime,phasej,stat1,slatd,slond,stelev,elatd,elond,deps,delta,rdum,iaz,kres,stdev,rav
#write(*,*) evday,evtime,phasej,stat1


# INPUTS for each ray to be traced:
# evday,			event yearjdy for ID only
# evtime,			event epoch time for ID only
# phasej,			phase eg P
# stat1,slatd,slond,stelev,	station name, lat, lon, elev
# elatd,elond,deps,		event lat lon depth
# delta,rdum,iaz,		epi dist (deg) <NOTUSED> ev to st azimuth
# kres,stdev,rav		tt resids, standard deviations, x-corr coef
#                                                               (not used)

## put in del and az if missing
if (delta == 0. .or. iaz == 0.) {
  call azidl(slatd*qwe,slond*qwe,elatd*qwe,elond*qwe,delta,deldist,azs,aze)
  delta = delta/qwe
  iaz = aze/qwe
}

## correct phase names 
##     - sr files have dbpick terminology truncated by dbcrosscr! ...sorry!
if (phasej == "PKPd") { phasej = "PKIKP" }
if (phasej == "PKIK") { phasej = "PKIKP" }

onset='i'
iflg=0

slat=slatd*qwe; slat=atan(0.993277*tan(slat))	# to rads & geocentric
slat=1.5707963-slat
slon=slond*qwe			
elat=elatd*qwe; 
elat=atan(0.993277*tan(elat))
elat=1.570796327-elat
elon=elond*qwe
iazo=iaz		# store for later
iaz=iaz*qwe
robres=kres
phasorg=phasej

# find/assign a station ID number - istind
if (irayno.eq.0) {
    nstat=1
    statst(1)=stat1
    statla(1)=slatd
    statlo(1)=slond
    statel(1)=stelev
}
do i=1,nstat {
  if(stat1.eq.statst(i)) {
      istind=i
      goto 735
  }
} 
nstat=nstat+1
statst(nstat)=stat1
statla(nstat)=slatd
statlo(nstat)=slond
statel(nstat)=stelev
istind=nstat
735 continue

# find/assign a event ID number - ievnt
if (irayno.eq.0) {
   nevt=1
   event(1)=evtime
   eventla(1)=elatd
   eventlo(1)=elond
}
do i=1,nevt {
   if(evtime.eq.event(i)) {
      ievnt=i
      goto 736
   }
} 
nevt=nevt+1
event(nevt)=evtime
eventla(nevt)=elatd
eventlo(nevt)=elond
ievnt=nevt
736 continue

#write(*,*)'nevt = ',nevt

irayno=irayno+1


   # calc epicentral dist
   d=acos( cos(slat)*cos(elat)+sin(slat)*sin(elat)*cos(slon-elon) )
   mddelta=d*RAD
   # mddelta in degrees!
   #write(*,*)'nev, mddelta',nev,mddelta,delta 
   #write(*,*)'st ',slat,slon,elat,elon 
   #write(*,*)
   #write(*,*)'   delta,mddelta: ',delta,mddelta


    # Get rayparamter from tables

    243 reas=0
    # Subroutine pss returns ray parameter from interpolation through 
    # ray tables created for iaspi91
    # phasej is phase name, mddelta is dist to station, rayparam is returned, 
    # reas is reason for failure      
    # dpdd is the derivative of p with respect to dist, used later for 
    # fine tuning to make sure ray
    # reaches exactly to the station
    if(debug)write(*,*)' phase (pre pss): ',phasej(1:6)
    call pss(phasej,deps,mddelta,rayparam,ftime,reas,iflg,dpdd)
    if(debug)write(*,*)'From TABLE ps= ',rayparam,ftime,' phase: ',phasej(1:6)

    ## some info for later
    vel1 = 2.25
    vel2 = 4.10
    vel3 = 3.50
    vel4 = 8.00
    rayangle1 = (180./3.14159)*asin(vel1*rayparam/111.11)
    rayangle2 = (180./3.14159)*asin(vel2*rayparam/111.11)
    rayangle3 = (180./3.14159)*asin(vel3*rayparam/111.11)
    rayangle4 = (180./3.14159)*asin(vel4*rayparam/111.11)
    write(68,*)phasej, rayparam,rayangle1,rayangle2,rayangle3,rayangle4

    # GET RAY TYPES
    do i=1,nrays
    jsray(i)=0
    do  j=1,nrays {
	   if(phasej.eq.name(j)) {
	     jsray(j)=1
	     idray=j
	     go to 232
	   }
    }
    write(*,*) 'RAY TYPE NOT MATCHED!!',phasej,nrays,nev,dum4
    write(99,*)'RAY TYPE NOT MATCHED!!',phasej,nrays,nev,dum4
    go to 223

    232 continue
    if(debug)write(*,*)'ray type matched'

    # DONT USE pwP at this TIME!
    if (idray==63) { 
      print *,'write to 99, idray=63 (pwP) ',phasej,delta,rayparam,reas
      write(99,224)phasej,delta,rayparam,reas 
      go to 223
    }

    ## RMA Aug 2000: turned S back on by commenting this out
    # # It will not yet work for PKS type waves
    # # becuase there is no symmetry
    # sum=0.0
    # do i=1,11
    #    sum=sum+ms(i,idray)
    # if (sum>0) {
    #    write(*,*)'write to 99, S component: ',phasej,delta,rayparam,reas
    #    write(99,224)phasej,delta,rayparam,reas
    #    go to 223
    # }


    # if (phasej=='pwP') {phasej='pP' }

    if (reas.ne.0) { 
	    print *,'write to 99, reas= ',reas
	    write(99,224)phasej,delta,rayparam,reas
	    go to 223 
    }

    # CREATE IASPI91 DISCONTINUITIES rd(1-12)
    call mdiaspd

    # GET EVENT DEPTH
    rsource=rdisc(1)-deps
    # Find source layer # (LS) and velocities in source level
    # Subrtouine mdiaspi returns vps and vss (velocities) in layer ls 
    # at depth rsouce
    ls=0
    call mdiaspi(rsource,ls,vps,vss,dummy)
    if(debug)write(*,*)'source layer=',ls


    # IDRAY is ray id number 1-61
    # CORRECT MULTIPLICITIES FOR pP, sS and P,S
    #write(*,*)'kfr,ls  ',kfr(idray),ls
    icps=0
    301 call msmp(mp,ms,mporg,msorg)
    call corm(kfr,ms,mp,ls,nrays,idray)

    ## ds is the step size for ray3d - choose vertical or ray step in ray3d
    #ds=min(25.,mddelta*10.)	
    ds=min(10.,mddelta*10.)
    #ds=min(5.,mddelta*10.)
    #ds=min(2.,mddelta*10.)

    # CALL RAY3D TO TRACE RAY, ALSO SEND KRF,MP,MS
    # Rich: statement call surfacecor (line 917 of raylib.r) commented out 
    #       as don't want to read earth topo
    # ray3d returns ray path as (r,theta) with station at north pole 
    # (doen't matter since 1d)
    # ray goes from station to event.
    time=0.0 # time will be sum of tcor values only for now! (Corrections only)
    call ray3d(rayparam,rsource,ds,itype,r,t,rdisc,lr,n,time,kfr,klr,mp,ms,idray,slat,slon,elat,elon,mddelta,input,err)

    if (err.ne.0) {
      write(*,*)'err.ne.0  err = ',err
      stop
    }

    if(debug)write(*,*)'EpiDist error: ',mddelta-t(n)/qwe
    # Check how close we are to event location, do we need to fine tune?
    # if more the .60 degree away not worth trying to fine tune, discard!
    # if more than .2 sec/deg difference in ps then discard 
    #if (abs(mddelta-t(n)/qwe)>=.60) { 
    if (abs(mddelta-t(n)/qwe)>=.90) { 
      write(96,222)phasej,idray,deps,mddelta,t(n)/qwe,mddelta-t(n)/qwe 
      write(*,*)'ray discarded: large EpiDist err ',phasej,' err: ',mddelta-t(n)/qwe
      goto 223  # ...next ray
    }

    #  if ray dist is more than .005 degrees from station event distance 
    #  then fine tune
    if (abs(mddelta-t(n)/qwe)>=.005) {   
            if(debug)write(*,*)' Fine tuning the tracing...'
	    rayparam=rayparam+ dpdd*(mddelta-t(n)/qwe)
            #	Here we are slightly modifying ps but not the time
            #	For accurate calculated time will have to change this
	    icps=icps+1
	    if (icps>3) { 
	      write(*,*)'ray discarded: EpiDist err (tried tuning) ',phasej,' err: ',mddelta-t(n)/qwe
	      write(96,222)phasej,idray,deps,mddelta,t(n)/qwe,mddelta-t(n)/qwe
	    goto 223
	    }
	    go to 301
    #	write(*,*)'new ps=',rayparam
    }

    ###############
    # RAY ACCEPTED! keep a list in this file fort.98

    numray=numray+1
    # write out all params for used rays
    write(78,888) evday,evtime,phasej,stat1,slatd,slond,stelev,elatd,elond,deps,delta,rdum,iazo,kres,stdev,rav

    ## Do station elevation correction - not needed if using station corrections in inversion
    ## I only correct for station elevation and bounce points above sea level.
    ## Rob has already corrected the residuals
    ## for boiunce points below sea level.  Sergei, change!
    #call elevcor(stelev,rayparam,tcor,idray)
    ## uses elevation on 5 min. intervals
    #time=time+tcor  # because of station elevation
    ## my time is not travel time but correction to robs residual, sergei
    #if (err==1) {write(49,*)'nev,station,phase=',phasej; write(49,*); go to 223}
    #if (err>1) {write(50,*)'nev,station,phase,idray=',phasej,idray; write(50,*); go to 223}

    # to check ray3Ds work:
    #write(*,*)'done with ray3d'
    #write(*,*)'n= ',n
    #do i=1,n
    #write(*,*)'r,t = ',r(i),t(i),t(i)*360./2.0/3.141593

    #call genxyz(rec lat (geogr degress),re long,source lat and long,r,t(from ray3d),lr,n)
    lat=elat; lon=elon

    #Costas: Zero number of derivatives for this ray
    ndvel=0
    
    # genxyz, convert (r,theta) to x,y,z and generate row of matrix
    call genxyz(slat,slon,lat,lon,r,t,lr,n,dtx,dty,dtz)
    
    ## some successful ray infp to screen
    #write(*,*)'sucessful ray number: ',irayno,ievnt,istind,'  ',phasorg,ndvel
    write(*,*)'success: ray# #-non-0 phase: ',irayno,ndvel,'  ',phasorg

    ### sum ray length - for info to screen only - RMA
    #dvelsum=0.0
    #do iii=1,ndvel { dvelsum=dvelsum+dvel(iii) }
    #write(*,*)'         ray length (in model box): ',dvelsum

    ## Write row of derivative matrix ie the matrix!
    ##    phasorg is the phase (ie. P, PKPab, etc. character*8)
    ##    onset is i,e char*1
    ##    iflgorg =0 first arriving P phase used;
    ##            =1 first arriving P not used;
    ##            =2 later phase used;
    ##            =3 later phase not used
    ##    flg 0,1,2, num derivatives, num of event, sta num from stat_elev
    ##    (idvel(j),j=1,ndvel) - indices
    ##    (dvel(j),j=1,ndvel)  - derivatives
    ##    dtx,dty,dtz - hypo center derivatives even if outside
    write(22) flag,ndvel,ievnt,istind,phasorg,onset,iflgorg
    write(22) (idvel(j),j=1,ndvel)
    write(22) (dvel(j),j=1,ndvel)
    write(22) dtx,dty,dtz

    ## Write travel-time residual to resfile
    ##     robres is residual at this point corrected!!
    ##      + other variables for stats!
    write(23) robres,flag,ndvel,ievnt,istind,phasorg,onset,iflgorg

    ## Write data certainty info for use in covariance matrix
    write(24,*)stdev	# stdev of the data 

    ## Write station info
    write(25,*)stat1,slatd,slond


    #write(*,*)'pause;'
    #if (lp==0) {
    #read(*,*)lp
    #}


    go to 223    # NEXT RAY!

    999 continue

    # output station ids if Richards data
    if(iflag(1).eq.1) { 
      open(67,file='stat.ids')
      do i = 1,nstat {
	write(67,'(i4,x,a8,3f10.4)')i,statst(i),statla(i),statlo(i),statel(i)
      }
      close(67)
      open(68,file='evnt.ids')
      do i = 1,nevt {
	write(68,'(i4,x,f16.5,2f10.4)')i,event(i),eventla(i),eventlo(i)
  }
  close(68)

}
write(*,*)'done, number of rays used: ',numray
end


#---------------------------------------------------
# Deal: background discontinuity levels for iasp
# from iasp91 GJI 105,429-469, 1991
# sets up discontinuities in 1d model

subroutine mdiaspd
dimension rd(12)
include 'comm.ray3'
data rd/6371.,6351.,6336.,6251.,6161.,5961.,5711.,5611.,3631.,3482.,1217.1,0.0/
ndisc=12
do i=1,ndisc {
  rdisc(i)=rd(i)
  rdisc2(i)=rd(i)**2
}
return
end



#-----------------------------------------------


#-------------------------------------------------------
#  This subroutine is not used any more
       subroutine getray(name,jsray,nrays)

#      reads rays from input
       dimension jsray(65)
       character*8 name(65),r

       pp1=1.0e30
       pp2=0.
10     continue
       write(6,*) 'Give ray name or <stop>:'
       read(5,fmt='(a8)') r
       if(r.eq.'stop'.or.r.eq.'xxxx') return
#
       if(r.eq.'help') {
         write(6,fmt='(5(a8,2x))') (name(n),n=1,65)
         goto 10
       }
       if(r.eq.'all') {
               do  j=1,nrays
                 jsray(j)=1
               goto 10
       }

       if(r.eq.'core'.or.r.eq.'mantle') {
               jsm=0
               jsc=1
               if(r.eq.'mantle') {
                       jsm=1
                       jsc=0
               }

               do  j=1,nrays
                 jsray(j)=jsc

               jsray(1)=jsm
               jsray(7)=jsm
               jsray(8)=jsm
               jsray(9)=jsm
               jsray(10)=jsm
               jsray(11)=jsm
               jsray(16)=jsm
               jsray(22)=jsm
               jsray(23)=jsm
               jsray(24)=jsm
               jsray(25)=jsm
               jsray(26)=jsm
               jsray(31)=jsm
               jsray(34)=jsm
               jsray(37)=jsm
               jsray(45)=jsm
               jsray(47)=jsm
               jsray(49)=jsm
               do  i=55,61
                 jsray(i)=jsm
               pp1=min(p1,253.841)
               pp2=max(p2,1846.667)
               goto 10
       }
       jj=0
       do  j=1,nrays {
       if(r.eq.name(j)) {
         jsray(j)=1
         jj=1
       }
	}
     
       if(jj.eq.0) {write(6,*) 'Ray type not present in the table'; go to 10 }
	return
       end
#----------------------------------------------------------

#----------------------------------------------------------------------
      subroutine corm (kfr,ms,mp,ls,nrays,idray)

#     Correct multiplicities for source depth in the standard ray
# Correct for symmetry arguments if ray at depth or bounces up
      dimension kfr(65),ms(11,65),mp(11,65)
	i=idray
#      do  i=1,nrays {
 	if (kfr(i)==1) go to 10
	if (kfr(i)==2) go to 30
 	if (kfr(i)==3) go to 50
 	if (kfr(i)==4) go to 70
#      GO TO (10,30,50,70), KFR(I)
#     p. Add 1 for every layer above the source.
   10 l=ls
   20 l=l-1
      if (l.LE.0) GO TO 100
      mp(l,i)=mp(l,i)+1
      go to 20
#     s. Add 1
   30 l=ls
   40 l=l-1
      if (l.EQ.0) GO TO 100
      ms(l,i)=ms(l,i)+1
      go to 40
#     P. Substract source layer, and above it.
   50 l=ls+1
   60 l=l-1
#	write(*,*)'in corm l=',l
      if (l.LE.0) go to 100
      mp(l,i)=mp(l,i)-1
      go to 60
#     S. Substract.
   70 l=ls+1
   80 l=l-1
      if (l.EQ.0) go to 100
      ms(l,i)=ms(l,i)-1
      GO TO 80
  100  continue
#    }
      END
#-----------------------------------------------------------------


subroutine  msmp(mp,ms,mporg,msorg)
dimension mp(11,65),ms(11,65),mporg(11,65),msorg(11,65)

# Used to reset multiplicities

do i=1,65
do j=1,11
ms(j,i)=msorg(j,i)
do i=1,65
do j=1,11
mp(j,i)=mporg(j,i)

return
end

#-----------------------------------------------------------------

	Subroutine deriv(s10,s10o,i1n,i2n,i3n,f1,f2,f3)
#c Function to compute velocity derivatives at a point, using the velocity
#c nodes surrounding it and the fractional distances from these nodes
#c Takes account of twodimensional degeneracy.
#c s10,s10o = Lengths of the ray segments for which velocity derivatives
#c	will be estimated, with respect to the point which connects these
#c	segments.
#c 
	common/der/dvel(NVEL),idvel(NVEL),ndvel
	common/gr/n1,n2,n3,x1,x2,y1,y2,z1,z2,dx,dy,dz,tc,phc
	dimension f(8),nod(8)
#c	The node convension is the same with P ploting program
#c	the first node is at x1,y2,z2, the x-index varies first,
#c	the y-index second & the z-index last.
	pro=0.5*(s10+s10o)
	i11=min(n1,i1n+1)
	i21=n2-min(n2,i2n+1)
	i31=n3-min(n3,i3n+1)
	i2=n2-i2n
	i3=n3-i3n
	i1=i1n
	f4=f3*f2
	f5=f3*f1
	f6=f2*f1
	f7=f6*f3
	f(8)=f7
	f(7)=f4-f7
	f(6)=f5-f7
	f(5)=f3-f5-f(7)
	f(4)=f6-f7
	f(3)=f2-f4-f(4)
	f(2)=f1-f5-f(4)
	f(1)=1.-f2-f3+f4-f(2)
	k1=n2*i3
	k2=n2*i31
	k3=(k1+i2)*n1
	k4=(k1+i21)*n1
	k5=(k2+i2)*n1
	k6=(k2+i21)*n1
	nod(1)=i1+k3
	nod(2)=i11+k3
	nod(3)=i1+k4
	nod(4)=i11+k4
	nod(5)=i1+k5
	nod(6)=i11+k5
	nod(7)=i1+k6
	nod(8)=i11+k6
	do i=1,8 {
	if (f(i).le. 0.) next
	ibr=0
	  do j=1,ndvel {
	   if (nod(i).eq.idvel(j)) {
	   dvel(j)=dvel(j)+pro*f(i)
	   ibr=1
	   }
	  if (ibr == 1) break
	  }
	if (ibr == 1) next
	ndvel=ndvel+1
	idvel(ndvel)=nod(i)
	dvel(ndvel)=pro*f(i)
	}
	return
	end

#-----------------------------------------------------------------
	Subroutine cell(xx,yy,zz,i1,i2,i3,f1,f2,f3)
# Subroutine to calculate node indices and distance ratios
# used in derivative interpolation for a point at xx,yy,zz
# Does not takes account of twodimensional degeneracy
# Does not check if point is within grid limits
common/gr/n1,n2,n3,x1,x2,y1,y2,z1,z2,dx,dy,dz,tc,phc
#****	i1, i2, i3 = Minimum corner Indices (closest to origin)
#	             of the cell containing point xx,yy,zz
#	f1, f2, f3 = Interpolation coefficients
#
x=xx-x1
y=yy-y1
z=zz-z1
i1=int(x/dx)+1
f1=mod(x,dx)/dx
i2=int(y/dy)+1
f2=mod(y,dy)/dy
i3=int(z/dz)+1
f3=mod(z,dz)/dz
return
end

#-----------------------------------------------------------------

Function boxx(lat,lon,h)
common/gr/n1,n2,n3,x1,x2,y1,y2,z1,z2,dx,dy,dz,tc,phc
common/gri/x1i,x2i,y1i,y2i,z1i,z2i
real lat,lon,h,x,y,z,slat
integer boxx
slat=1.57079632679-lat
call sphe2cart(slat,lon,h,x,y,z)
#write(*,*)x,y,z,x1i,x2i,y1i,y2i,z1i,z2i
boxx=0
if (x > x1 && x < x2 && y > y1 && y < y2 && z > z1 && z < z2) { boxx=1 } 
if (x > x1i && x < x2i && y > y1i && y < y2i && z > z1i && z < z2i) { boxx=2 }
return
end

#-----------------------------------------------------------------

#	Calculate cartesian coordinates from lat & lon

	subroutine sphe2cart(slat,slon,h,x,y,z)
	common/gr/n1,n2,n3,x1,x2,y1,y2,z1,z2,dx,dy,dz,tc,phc
	common/ang/cltc,sltc,clnc,slnc
	double precision clt,slt,cln,sln,cltc,sltc,clnc,slnc,sla,slo
	double precision cdif,sdif,t1,t2,t3,t4,t5,azi
#***********************************************************************
#
#	sla, slo = Input Latitude and Longitude in radians
#	h = Elevetion from sea surface in Km (positive up)
#	x, y & z = Output catresian coordinates in Km
#
#***********************************************************************
	sla = Dble(slat)
	slo = Dble(slon)
	clt = Dcos(sla)
	slt = Dsin(sla)
	cln = Dcos(slo)
	sln = Dsin(slo)
	cdif = Dcos(slo-dble(phc))
	sdif = Dsin(slo-dble(phc))
	t1 = clt*cln*cltc*clnc
	t2 = clt*sln*cltc*slnc
	t3 = slt*sltc
	sla = Dasin(t1+t2+t3)
	t4 = slt*cltc-clt*sltc*cdif
	t5 = clt*sdif
	if (t4 != 0.d0) 
	  azi = Datan2(t5,t4)
	else
	  azi = 3.141592654d0-dsign(1.57079632679d0,t5)
	slo = 1.57079632679d0-azi
	radi = 6371. + h
	x = radi * Dcos(sla) * Dcos(slo)
	y = radi * Dcos(sla) * Dsin(slo)
	z = radi * Dsin(sla)
	return
	end


#########################################################################
# OLD COMMENTS ON OLD INPUT FILE FORMATS 
#                         ....for code shoot (rather than shootsr)
#
# Input files:
#  
# file='flagfile'  
#         1       # number of flags
#         1       # 1st flag: RMAs format datafiles
#  
# file='names'    # names of other files: 1st 4 input, last 3 output
#   read(88,*)dfile1
#   read(88,*)dfile2
#   read(88,*)dfile4
#   read(88,*)dfile5
#   read(88,*)matfile             # matrix file (unformatted)
#   read(88,*)resfile             # residual file (unformatted)
#   read(88,*)velmod              # residual file (unformatted)
#  
# file='inputfile'        # model grid
#   read(12,230) n1,n2,n3         # Number of nodes in each direction
#   read(12,231) x1,y1,z1         # Minimum distance from center
#   read(12,231) x2,y2,z2         # Maximum distance from center
#   read(12,231) dx,dy,dz         # Step size in each direction
#   read(12,231) latc,lonc   # geographic latitude and long
#   NB the z-coord is RADIUS so for a box from 600km depth to
#      surface z1=5771. z2=6371.
#  
# also a series of data files ...very messy, should clean up!:
#  
# datafile1:
#   read(31,*,end=999)nev,phasej,evnt
#                 nev - event ref #
#                 phasej - phase eg P
#                 evnt - event name - not used
#  
# datafile2:
#   read(32,47)stat1,slat,slon,stelev,elat,elon,deps,delta,dum2,dum3,iaz        
# 	stat,stla,stlo,elev,evla,evlo,evdp,delta,delkm,azis,azie
#         everything in deg, some things not used
#  
# datafile3:
#  read(34,fmt='(a4,4f10.2,f10.3)') kstnm,dtest,pred,tstar,stdev,rav
#        only use stdev - certaintly of data (just written out)
#  
# datafile4:
#   read(35,*) offset,kres
#         offset (not used) and travel-time residual (just written out)
#