bruno
this has probably been fixed, try attached version
jens
Den 2020-11-23 18:54 skreiv Bruno Faria:
Dear All,
Despite there is about 20 years that I'm working with Seisan, I'm
facing some difficulties to install the version 11 on Centos 8 (64
bits).
Uncompressing Seisan file and and configurating (as I did so many
times) the programs (e.g. mulplt) don't work. Trying to compile with
gfortran (ver. 8.3.1) the routines in both folders, PRO and LIB the
compiling process is finding a problem with the subroutine
hyposub6.for:73:26:
/ call ttcal(n_layer,v,z,d,vr,grad,z_hypo,pcrit,phi,dcrt,tcrt,
1
Error: Actual argument contains too few elements for dummy argument
‘v’ (150/200) at (1)/
I'm missing something. I'm aware that the gfortran is new, so probably
there is and option that I'll have to add to Makefile in order to
compile SEISAN.
May I ask, please, what is missing me?
Thank you very much in advance for your help.
Best regards,
Bruno
_______________________________________________
seisan mailing list
[email protected]
https://mailman.uib.no/listinfo/seisan
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
[email protected]
https://mailman.uib.no/listinfo/seisan
