|
Hi Rod, we fixed some of the problems which show up when compiling on mac last year, can you copy the attached file to your LIB folder and try again? I think you will get some other errors, please let me know what these are and I may have to send another few files. This should fix the first problem. The second should be resolved when LIB all compiles and you get the seisan.a. Cheers, Lars
On 26/08/2019 17.16, Roderick Stewart
wrote:
Haa anyone had success compiling seisan on Macos Mojave? -- ------------------------------------------------------------------------- Lars Ottemöller Department of Earth Science, University of Bergen, Norway web http://folk.uib.no/lot081 email lars.ottemol...@geo.uib.no phone +47-5558 2616 ------------------------------------------------------------------------- |
c**************************************************************
c subroutine to find travel times for a given offset using a
c gradient model. Inputs are the same as for dtdx2
c Presently, the routine only calculates the minimum
c travel time, regardless of ptype
c Refraction can also only occur in the bottom (half-space)
c layer
c no phase ID is returned in phsid
c
c Barry Lienert June 1998
c
c changes
c
c sep 98 by jh : ----------- seisan version 7.0 check ----------------
c no changes
c
c oct 28 99 bmt : linux changed, save and *
c nov 19 99 jh : nl and parm in common block
c may 15 00 jh : fix bugs cause by above
c
subroutine dtdxg(xh,x0,ptype,nmoho,nconrad,iulst,tmin,
&delta,ann,iflag,phsid)
save
include 'hypparm.inc'
character*1 ptype
character*8 phsid
c dimension dcrit(200),z(200),grad(200),vr(200)
dimension dcrit(nlayer),z(nlayer),grad(nlayer),vr(nlayer) ! lo 16/9/2018
c data pi2/1.57079/
data pi2/1.57079633/
c get epicentral distance from event to station
call delaz(x0(2),x0(1),delta,dedeg,azz0,xh(2),xh(1))
delta=abs(delta)
c calculate needed parameters
nn=nl*2-1
n_layer=nl
c n_layer=(nn+1)/2 ! nov 19 jh
n_iter=20 !# iterations to find delta
phi_incr=.001 !slowness increment in iterations
sum=0.0
do i=1,n_layer
cx v(i)=parm(i) !layer velocity
cx d(i)=parm(nl+i) !layer thickness
z(i)=sum !layer depth
sum=sum+parm(n_layer+i)
end do
z_hypo=xh(3)
c find n_hypo, the layer that the hypocenter is in
n_hypo=n_layer
do i=1,n_layer-1
d(i)=z(i+1)-z(i)
vr(i)=v(i+1)/v(i)
grad(i)=(v(i+1)-v(i))/d(i)
if(z(i).le.z_hypo)n_hypo=i
end do
c now find v_hypo
if(n_hypo.lt.n_layer)then
v_hypo=v(n_hypo)+grad(n_hypo)*(z_hypo-z(n_hypo))
else
v_hypo=v(n_layer)
endif
c find v_max, max vel above hypocenter
v_max=0.0
do i=1,n_hypo
if(v(i).gt.v_max)v_max=v(i)
end do
if(delta.gt.0.0)then
c find the critical distance, dcrt & time, tcrt
pcrit=1./v(n_layer)
phi=2.*pi2-asin(pcrit*v(n_hypo))
call ttcal(n_layer,v,z,d,vr,grad,z_hypo,pcrit,phi,dcrt,tcrt,
& n_hypo,v_hypo,v_max,iflagd)
c extend the travel time from dcrt to delta using v(n_layer)
if(delta.ge.dcrt)then
tmin=tcrt+(delta-dcrt)/v(n_layer)
iflagd=1
return
endif
c find p for the ray that exits the hypocenter layer at a
c horizontal offset delta (see Lee & Stewart, 1981, p 93)
eta=z_hypo-z(n_hypo)
xc=0.5*(delta*delta-2.*eta*v(n_hypo)/grad(n_hypo)-eta*eta)/delta
c zc=v(1)*d(1)/(v(2)-v(1))
zc=v(n_hypo)/grad(n_hypo)
c phi is the ray direction measured clockwise from the upward vertical
phi=atan2(-zc-eta,-xc)
p=abs(sin(phi))/v(n_hypo)
c phi=atan(2.*zc/delta)
else
c zero delta case
phi=2.*pi2
p=0.0
if(z_hypo.eq.0.0)then
tmin=0.0
return
endif
endif
if(z_hypo.ne.0.0)then
if(delta.eq.0.0)then
call ttcal(n_layer,v,z,d,vr,grad,z_hypo,p,phi,dist,tmin,
& n_hypo,v_hypo,v_max,iflagd)
c p=0: don't need to iterate here
return
else
c try a ray at the p calculated to make the ray leave the
c hypocenter layer at distance delta
phi0=phi
p0=p
call ttcal(n_layer,v,z,d,vr,grad,z_hypo,p,phi,dist,tmin,
& n_hypo,v_hypo,v_max,iflagd)
c ray did not come out: try smaller phi's
if(iflagd.eq.0)then
itr=0
phi=phi0
do while(iflagd.eq.0.and.itr.lt.n_iter)
phi=phi-phi_incr
p=abs(sin(phi))/v(n_hypo)
call ttcal(n_layer,v,z,d,vr,grad,z_hypo,p,phi,dist,tmin,
& n_hypo,v_hypo,v_max,iflagd)
itr=itr+1
end do
if(iflagd.eq.0)then !try other direction
itr=0
phi=phi0
do while(iflagd.eq.0.and.itr.lt.niter)
phi=phi+phi_incr
p=abs(sin(phi))/v(n_hypo)
call ttcal(n_layer,v,z,d,vr,grad,z_hypo,p,phi,dist,tmin,
& n_hypo,v_hypo,v_max,iflagd)
itr=itr+1
end do
endif
else
c dist0=dist
phi=phi0
do while(iflagd.eq.0.and.itr.lt.niter)
phi=phi-phi_incr
p=abs(sin(phi))/v(n_hypo)
call ttcal(n_layer,v,z,d,vr,grad,z_hypo,p,phi,dist,tmin,
& n_hypo,v_hypo,v_max,iflagd)
c write(iulst,'(3f12.4)')delta,dist0,dist
end do
if(iflagd.eq.0)then
do while(iflagd.eq.0.and.itr.lt.niter)
phi=phi+phi_incr
p=abs(sin(phi))/v(n_hypo)
call ttcal(n_layer,v,z,d,vr,grad,z_hypo,p,phi,dist,tmin,
& n_hypo,v_hypo,v_max,iflagd)
end do
endif
itr=0
do while (abs(dist-delta).gt.0.01.and.itr.lt.niter)
if(dist.ne.dist0)then
slope=(phi-phi0)/(dist-dist0)
if(slope.gt..1)slope=.1
if(slope.lt.-.1)slope=-.1
endif
phi0=phi
p0=p
phi=phi0+slope*(delta-dist)
p=abs(sin(phi))/v(n_hypo)
dist0=dist
call ttcal(n_layer,v,z,d,vr,grad,z_hypo,p,phi,dist,
& tmin,n_hypo,v_hypo,v_max,iflagd)
c write(iulst,'(i3,4f10.4)')itr,p,phi,dist0,dist
itr=itr+1
end do
endif
endif
else
call ttcal(n_layer,v,z,d,vr,grad,z_hypo,p,phi,dist,tmin,n_hypo,
& v_hypo,v_max,iflagd)
endif
c write(iulst,*)p,delta,dist,itr
return
end
c********************************************************************
c calculates distance and travel time as a function of ray parameter
c for a gradient model
c modified from Fred Klein's TTCAL routine by Barry Lienert, June 98
c********************************************************************
subroutine ttcal(n_layer,v,z,d,vr,grad,z_hypo,p,phi,dist,tp,
& n_hypo,v_hypo,v_max,iflagd)
c Inputs: n_layer=number of layers in velocity model
c v(i)=velocities in km/s
c z(i)=layer top positions in km
c z_hypo=hypocentral depth in km
c p=ray parameter in s/km
c
c Outputs: dist=delta in km
c tp=travel time in s
c iflagd=1 for valid arrival, zero otherwise
c
save
c parameter (maxlayer=200)
parameter (maxlayer=150) ! lo match hypparm.inc maxlayer 16/9/2018
real cp(maxlayer),vr(maxlayer)
real v(maxlayer),z(maxlayer),d(maxlayer),grad(maxlayer)
c data pi2/1.57079/
data pi2/1.57079633/
c--case of a ray travelling straight up (p=0)
if(p.eq.0.0)then
phid=0.
dist=0.
tp=(z_hypo-z(n_hypo))/v(n_hypo)
if (grad(n_hypo).ne.0.)then
tp=alog(v_hypo/v(n_hypo))/grad(n_hypo)
endif
if (n_hypo.gt.1)then
do i=1,n_hypo-1
temp=d(i)/v(i)
if (grad(i).ne.0.) temp=alog(vr(i))/grad(i)
tp=tp+temp
end do
endif
tr=tp
iflagd=1
return
endif
c sph is sin(phi) where phi is angle at which ray leaves the hypocenter
sph=sin(phi)
cph=cos(phi)
p=sph/v_hypo
c case of no upward ray
if(p.ge.1./v_hypo)then
tp=-1.
dist=-1.
iflagd=0
return
endif
c vb is the bottoming velocity
vb=1./p
c--n_bottom is the layer in which downgoing ray bottoms
n_bottom=1
do while (v(n_bottom+1).lt. vb.and.n_bottom.lt.n_layer)
n_bottom=n_bottom+1
end do
c--calc depth at ray bottom, zb
zb=z(n_bottom)
if (grad(n_bottom).ne.0.) zb=zb+(vb-v(n_bottom))/grad(n_bottom)
c++++++++++++ ray lost to waveguide ++++++++++++
c--now have case of ray going into waveguide formed by a lvz
c--flag the nonexistent time and dist for this case
if (vb.le.v_max .and. v_hypo.lt.v_max)then
tp=-1.
dp=-1.
tr=-1.
n_bottom=0
iflagd=0
return
endif
c++++++++++++ calculate time and distance ++++++++++++
c--calc t and d since it is valid to do so. consider all 3 cases:
c 1 ray nearly horiz in layer of zero gradient
c 2 upgoing ray
c 3 downgoing ray
c--for each term must consider zero and non-zero gradients as sep cases
c--case 1. nearly horiz non-curving ray+++++++++++
if (grad(n_hypo).eq.0..and.phi.le.0.001)then
c--set dist to be some arbitrary but large no
dist=100000.
tp=100000./v(n_hypo)
iflagd=1
return
endif
c--case 2. upgoing ray+++++++++++++++++++++++++++
c if (phi.le.pi2+.0005)then
if (phi.le.pi2)then
c--calc cosines of emergence angles at each interface
ss=1.0
if(phi.gt.pi2.and.phi.le.2.*pi2)ss=-1.
if(phi.gt.2.*pi2.and.phi.le.3.*pi2)ss=-1.
if(phi.gt.3.*pi2.and.phi.le.4.*pi2)ss=1.
do i=1,n_hypo
cp(i)=ss*sqrt(abs(1.-(v(i)*p)**2))
end do
c--contribution to tp and dist of layer in which hypo lies
if (grad(n_hypo).gt.0.)then
dist=(cp(n_hypo)-cph)/(p*grad(n_hypo))
if(cph.gt.-1.0)then
tp=alog(v_hypo*(1.+cp(n_hypo))/(v(n_hypo)
& *(1.+cph)))/grad(n_hypo)
else
tmin=-1.0
dist=-1.0
iflagd=0
return
endif
else
dist=(z_hypo-z(n_hypo))*sph/cph
tp=(z_hypo-z(n_hypo))/(cph*v_hypo)
endif
c--add contributions from layers above hypo
if (n_hypo.gt.1)then
do i=1,n_hypo-1
if (grad(i).gt.0.)then
dist=dist+(cp(i)-cp(i+1))/(p*grad(i))
tp=tp+alog(vr(i)*(1.+cp(i))/(1.+cp(i+1)))/grad(i)
else
dist=dist+d(i)*p*v(i)/cp(i)
tp=tp+d(i)/(cp(i)*v(i))
endif
end do
endif
iflagd=1
return
else
c case 3 downgoing ray +++++++++++++++++++++++++++++
c--calc cosines of emergence angles at layer interfaces
do i=1,n_bottom
cp(i)=sqrt(1.-(v(i)*p)**2)
end do
c--contribution to t and d from layer in which ray bottoms
if (grad(n_bottom).gt.0.)then
dist=2.*cp(n_bottom)/(p*grad(n_bottom))
tp=2.*alog((1.+cp(n_bottom))/(v(n_bottom)*p))/grad(n_bottom)
else
c--a ray can't bottom in a homogeneous layer
c--if this case should arise, reflect ray from top of the layer
dist=0.
tp=0.
endif
c--subtract contribution from path immed above hypo in its layer
if (grad(n_hypo).gt.0.)then
dist=dist-(cp(n_hypo)-cph)/(p*grad(n_hypo))
tp=tp-alog(v_hypo*(1.+cp(n_hypo))/(v(n_hypo)
& *(1.+cph)))/grad(n_hypo)
else
dist=dist-(z_hypo-z(n_hypo))*sph/cph
tp=tp-(z_hypo-z(n_hypo))/(cph*v_hypo)
endif
c--sum terms over layers above hypo
if (n_hypo.gt.1)then
do i=1,n_hypo-1
if (grad(i).gt.0.)then
dist=dist+(cp(i)-cp(i+1))/(p*grad(i))
tp=tp+alog(vr(i)*(1.+cp(i))/(1.+cp(i+1)))/grad(i)
else
dist=dist+d(i)*p*v(i)/cp(i)
tp=tp+d(i)/(cp(i)*v(i))
endif
end do
endif
c--sum terms over layers between hypo and ray bottom (if any)
c--include layer in which hypo lies
if (n_hypo.ne.n_bottom)then
do i=n_hypo,n_bottom-1
if (grad(i).ne.0.)then
dist=dist+2.*(cp(i)-cp(i+1))/(p*grad(i))
tp=tp+2.*alog(vr(i)*(1.+cp(i))/(1.+cp(i+1)))/grad(i)
else
dist=dist+2.*d(i)*p*v(i)/cp(i)
tp=tp+2.*d(i)/(cp(i)*v(i))
endif
end do
endif
iflagd=1
return
endif
end
_______________________________________________ seisan mailing list seisan@geo.uib.no https://mailman.uib.no/listinfo/seisan
