Wanderelch
/
ids2020
1
1
Fork 0
Code Issues 1 Pull Requests Releases Wiki Activity
Browse Source

IDS Update 2021-08-04

main
mmetz 11 months ago
parent
30c5452213
commit
5fdd5d955b
7 changed files with 971 additions and 976 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 2
      README.md
  2. 4
      src/idsalloc.f
  3. 192
      src/idsinverse.f
  4. 16
      src/idskernel.f
  5. 937
      src/idsmodel.f
  6. 28
      src/idssmdat.f
  7. 768
      src/idssmgrn.f

2
README.md
Unescape View File

@ -32,7 +32,7 @@ Packaging has been done by Malte Metz.
```
autoreconf -i # only if 'configure' script is missing
./configure
./configure # optional --prefix=/new/bin/directory/
make
sudo make install
```

4
src/idsalloc.f
Unescape View File

@ -24,6 +24,8 @@ c===================================================================
parameter(TPREMIN=10.0d+00,TPSTMIN=40.0d+00)
real*8 HYPOERR
parameter(HYPOERR=1.5d+04)
integer*4 nbpfids
parameter(nbpfids=3)
c===================================================================
c global variables
c===================================================================
@ -52,7 +54,7 @@ c
c IDS parameters
c
integer*4 nt,nr,ng,ntgrn,nf,itmax,nsnap,isfhyp
integer*4 itermax,nbpfg
integer*4 itermax,nbpfg,nearfield
real*8 fcut1,fcut2,fcut1g,fcut2g,swindow,twindow,dt,df
real*8 patchsize,twgrn,dtgrn,vrpeak,vsmean,mw,trise,tsource
logical*2 highpass

192
src/idsinverse.f
Unescape View File

@ -1,97 +1,97 @@
subroutine idsinverse(ierr)
use idsalloc
implicit none
integer*4 ierr
c
integer*4 i,ism,isf,jsf
real*8 alfa,beta,sigma,summa,ttp,tkfp,slwp,tts,tkfs,slws
real*8 re,rn,dis,dis0,azi,azi0,plg,plg0,x,x0,y,y0,z,z0
c
c define neighboring patches
c
do isf=1,nsf
rsmth(isf)=dmax1(0.25d0*vsmean/fcut2,1.5d0*patchsize)
100 i=0
do jsf=1,nsf
if(jsf.ne.isf.and.dis3dsf(isf,jsf).le.rsmth(isf))then
i=i+1
isfnb(i,isf)=jsf
endif
enddo
if(i.le.0)then
stop ' Error in idsinverse: isolated fault patch!'
else if(i.lt.4)then
rsmth(isf)=rsmth(isf)+patchsize
goto 100
endif
nsfnb(isf)=i
enddo
c
allocate(smvar(3,nsm),stat=ierr)
if(ierr.ne.0)stop ' Error in idsinverse: smvar not allocated!'
allocate(dsmvar(3,nsm),stat=ierr)
if(ierr.ne.0)stop ' Error in idsinverse: dsmvar not allocated!'
c
if(itermax.le.0)then
call idsforward(ierr)
else
write(*,'(a,i4,a)')' initial estimate of rise time <= ',
& idnint(trise),' s'
write(*,'(a,i4,a)')' signal window of data used <= ',
& idnint(tsource),' s'
do isf=1,nsf
dis0=dis3dsf(isfhyp,isf)
azi0=azisf2hyp(isf)*DEG2RAD
plg0=plgsf2hyp(isf)*DEG2RAD
vrpeak=0.d0
summa=0.d0
do jsf=1,nsf
dis=dis3dsf(isfhyp,jsf)
if(dis.le.dis0)then
x0=dis*dsin(plg0)*dcos(azi0)
y0=dis*dsin(plg0)*dsin(azi0)
z0=dis*dcos(plg0)
azi=azisf2hyp(jsf)*DEG2RAD
plg=plgsf2hyp(jsf)*DEG2RAD
x=dis*dsin(plg)*dcos(azi)
y=dis*dsin(plg)*dsin(azi)
z=dis*dcos(plg)
if(dsqrt((x-x0)**2+(y-y0)**2+(z-z0)**2).le.patchsize)then
vrpeak=vrpeak+sfvs(jsf)*sflen(jsf)*sfwid(jsf)
summa=summa+sflen(jsf)*sfwid(jsf)
endif
endif
enddo
c
vrpeak=vrpeak/summa
c
alfa=dis3dsf(isfhyp,isf)/(VR2VSLW*vrpeak)
beta=dis3dsf(isfhyp,isf)/(VR2VSUP*vrpeak)
sigma=dmax1(HYPOERR/sfvp(isfhyp),alfa-beta)
it1sf(isf)=1+idint(dmax1(0.d0,alfa-sigma)/dt)
it2sf(isf)=1+idint(dmin1(alfa+trise,tsource)/dt)
c vrpeak=vrpeak/summa
c
c alfa=dis3dsf(isfhyp,isf)/(0.5d0*VRVSRATIO*vrpeak)
c it1sf(isf)=1+idint(dmax1(0.d0,alfa-trise)/dt)
c it2sf(isf)=1+idint(dmin1(alfa+trise,tsource)/dt)
c
do ism=1,nsm
call disazi(REARTH,sflat(isf),sflon(isf),
& smlat(ism),smlon(ism),rn,re)
dis=dsqrt(re*re+rn*rn)
call idspstt(dis,sfdep(isf),ttp,tkfp,slwp,tts,tkfs,slws)
it1sf(isf)=max0(it1sf(isf),1+idint((tpsm(ism)-ttp)/dt))
it2sf(isf)=min0(it2sf(isf),
& 1+idint((tsource+tpsm(ism)-ttp)/dt))
if(dis3dsf(isfhyp,isf).le.HYPOERR)then
it1sf(isf)=max0(1,it1sf(isf)-
& idint(HYPOERR/sfvp(isfhyp)/dt))
endif
enddo
enddo
call idskernel(ierr)
endif
c
return
subroutine idsinverse(ierr)
use idsalloc
implicit none
integer*4 ierr
c
integer*4 i,ism,isf,jsf
real*8 alfa,beta,sigma,summa,ttp,tkfp,slwp,tts,tkfs,slws
real*8 re,rn,dis,dis0,azi,azi0,plg,plg0,x,x0,y,y0,z,z0
c
c define neighboring patches
c
do isf=1,nsf
rsmth(isf)=dmax1(0.25d0*vsmean/fcut2,1.5d0*patchsize)
100 i=0
do jsf=1,nsf
if(jsf.ne.isf.and.dis3dsf(isf,jsf).le.rsmth(isf))then
i=i+1
isfnb(i,isf)=jsf
endif
enddo
if(i.le.0)then
stop ' Error in idsinverse: isolated fault patch!'
else if(i.lt.4)then
rsmth(isf)=rsmth(isf)+patchsize
goto 100
endif
nsfnb(isf)=i
enddo
c
allocate(smvar(3,nsm),stat=ierr)
if(ierr.ne.0)stop ' Error in idsinverse: smvar not allocated!'
allocate(dsmvar(3,nsm),stat=ierr)
if(ierr.ne.0)stop ' Error in idsinverse: dsmvar not allocated!'
c
if(itermax.le.0)then
call idsforward(ierr)
else
write(*,'(a,i4,a)')' initial estimate of rise time <= ',
& idnint(trise),' s'
write(*,'(a,i4,a)')' signal window of data used <= ',
& idnint(tsource),' s'
do isf=1,nsf
dis0=dis3dsf(isfhyp,isf)
azi0=azisf2hyp(isf)*DEG2RAD
plg0=plgsf2hyp(isf)*DEG2RAD
vrpeak=0.d0
summa=0.d0
do jsf=1,nsf
dis=dis3dsf(isfhyp,jsf)
if(dis.le.dis0)then
x0=dis*dsin(plg0)*dcos(azi0)
y0=dis*dsin(plg0)*dsin(azi0)
z0=dis*dcos(plg0)
azi=azisf2hyp(jsf)*DEG2RAD
plg=plgsf2hyp(jsf)*DEG2RAD
x=dis*dsin(plg)*dcos(azi)
y=dis*dsin(plg)*dsin(azi)
z=dis*dcos(plg)
if(dsqrt((x-x0)**2+(y-y0)**2+(z-z0)**2).le.patchsize)then
vrpeak=vrpeak+sfvs(jsf)*sflen(jsf)*sfwid(jsf)
summa=summa+sflen(jsf)*sfwid(jsf)
endif
endif
enddo
c
vrpeak=vrpeak/summa
c
alfa=dis3dsf(isfhyp,isf)/(VR2VSLW*vrpeak)
beta=dis3dsf(isfhyp,isf)/(VR2VSUP*vrpeak)
if(nearfield.gt.0)then
sigma=dmax1(HYPOERR/sfvp(isfhyp),alfa-beta)
it1sf(isf)=1+idint(dmax1(0.d0,alfa-sigma)/dt)
it2sf(isf)=1+idint(dmin1(alfa+trise,tsource)/dt)
else
it1sf(isf)=1+idint(dmax1(0.d0,beta)/dt)
it2sf(isf)=1+idint(dmin1(beta+trise,tsource)/dt)
endif
c
do ism=1,nsm
call disazi(REARTH,sflat(isf),sflon(isf),
& smlat(ism),smlon(ism),rn,re)
dis=dsqrt(re*re+rn*rn)
call idspstt(dis,sfdep(isf),ttp,tkfp,slwp,tts,tkfs,slws)
it1sf(isf)=max0(it1sf(isf),1+idint((tpsm(ism)-ttp)/dt))
it2sf(isf)=min0(it2sf(isf),
& 1+idint((tsource+tpsm(ism)-ttp)/dt))
if(dis3dsf(isfhyp,isf).le.HYPOERR)then
it1sf(isf)=max0(1,it1sf(isf)-
& idint(HYPOERR/sfvp(isfhyp)/dt))
endif
enddo
enddo
call idskernel(ierr)
endif
c
return
end

16
src/idskernel.f
Unescape View File

@ -685,6 +685,14 @@ c
do isf=1,nsf
clust(isf)=dmin1(1.d0,clust(isf)/alpha)
enddo
c
if(iter.eq.1)then
do isf=1,nsf
do i=1,nf
sstf(i,isf)=sstf(i,isf)*dcmplx(clust(isf),0.d0)
enddo
enddo
endif
c
c suppress peak slips
c
@ -699,14 +707,6 @@ c
endif
enddo
endif
c
if(iter.eq.1)then
do isf=1,nsf
do i=1,nf
sstf(i,isf)=sstf(i,isf)*dcmplx(clust(isf),0.d0)
enddo
enddo
endif
c
c re-scaling
c

937
src/idsmodel.f
Unescape View File

@ -1,471 +1,468 @@
subroutine idsmodel(ierr)
use idsalloc
implicit none
integer*4 ierr
c
integer*4 i,j,ind,ig,ir,isf,jsf,iztp,nzmod,nsfred
real*8 pe,pn,sma,smb,smc,bga,bgc,re,rn,d2,d2min
real*8 st,di,t1,t2,dis1,dis2,farea,fsize,pswid
real*8 x,y,x0,y0,dx,dy
real*8 rstk(3),rddp(3)
logical*2 head,diff
character*180 text
c
do ind=120,1,-1
if(smgrndir(ind:ind).ne.' ')goto 100
enddo
100 continue
if(ind.ge.120)then
stop ' Error in idsmodel: too long folder name of smgrndir!'
endif
if(smgrndir(ind:ind).ne.'/'.or.
& smgrndir(ind:ind).ne.'\')then
smgrndir=smgrndir(1:ind)//'/'
ind=ind+1
endif
open(10,file=smgrndir(1:ind)//'GreenInfo.dat',status='old')
call skipdoc(10)
read(10,*)twgrn,dtgrn,ntgrn
if(ntgrn.le.0)stop ' Error in idsmodel: wrong value for ntgrn!'
c
call skipdoc(10)
read(10,*)nr
if(nr.lt.2)then
stop ' Error in idsmodel: GF distance coverage insufficient!'
endif
c
allocate (r(nr),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: r not allocated!'
c
read(10,*)(r(i),i=1,nr)
do i=1,nr
r(i)=r(i)*KM2M
enddo
c
call skipdoc(10)
read(10,*)ng
if(ng.lt.2)then
stop ' Error in idsmodel: GF depth coverage insufficient!'
endif
c
allocate (deps(ng),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: deps not allocated!'
allocate (smgrnfile(ng),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: smgrnfile not allocated!'
c
read(10,*)(deps(ig),smgrnfile(ig),ig=1,ng)
do ig=1,ng
smgrnfile(ig)=smgrndir(1:ind)//smgrnfile(ig)
deps(ig)=deps(ig)*KM2M
enddo
c
c Layered earth model
c
call skipdoc(10)
read(10,*)nlayer
c
allocate(zpmod(nlayer),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: zpmod not allocated!'
allocate(hpmod(nlayer),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: hpmod not allocated!'
allocate(vpmod(nlayer),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: vpmod not allocated!'
allocate(vsmod(nlayer),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: vsmod not allocated!'
allocate(romod(nlayer),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: romod not allocated!'
allocate(mumod(nlayer),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: mumod not allocated!'
c
call skipdoc(10)
do i=1,nlayer
read(10,*)j,hpmod(i),vpmod(i),vsmod(i),romod(i)
hpmod(i)=hpmod(i)*KM2M
vpmod(i)=vpmod(i)*KM2M
vsmod(i)=vsmod(i)*KM2M
romod(i)=romod(i)*KM2M
mumod(i)=romod(i)*vsmod(i)**2
enddo
zpmod(1)=0.d0
nzmod=1
do i=2,nlayer
zpmod(i)=zpmod(i-1)+hpmod(i-1)
if(zpmod(i).le.hypdep+stdismax)then
nzmod=i
endif
enddo
c
call skipdoc(10)
read(10,*)tptable
tptable=smgrndir(1:ind)//tptable
call skipdoc(10)
read(10,*)tstable
tstable=smgrndir(1:ind)//tstable
c
close(10)
c
allocate (tp(nr,ng),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: tp not allocated!'
allocate(tkftp(nr,ng),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: tkftp not allocated!'
allocate(slwtp(nr,ng),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: slwtp not allocated!'
allocate (ts(nr,ng),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: ts not allocated!'
allocate(tkfts(nr,ng),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: tkfts not allocated!'
allocate(slwts(nr,ng),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: slwts not allocated!'
c
open(11,file=tptable,form='unformatted',status='old')
read(11)((tp(ir,ig),tkftp(ir,ig),slwtp(ir,ig),ir=1,nr),ig=1,ng)
close(11)
open(12,file=tstable,form='unformatted',status='old')
read(12)((ts(ir,ig),tkfts(ir,ig),slwts(ir,ig),ir=1,nr),ig=1,ng)
close(12)
c
c empirical scaling laws
c
trise=10.d0**(0.5d0*mwini-2.32d0)
fsize=10.d0**(0.5d0*mwini+1.16d0)
c
fcut2=dmax1(FCUTUP,dmin1(0.75d0*fcut2g,5.0d0/trise))
fcut1=FCUTLW
c
tsource=2.0d0*trise
tpre=dmax1(TPREMIN,fsize/VPREF)
tpst=dmax1(TPSTMIN,fsize/VSREF)
ttap=dmax1(2.5d0/fcut2g,0.25d0*trise)
c
swindow=tsource+dmax1(trise,tpst)
twindow=swindow+tpre+tpst
c
if(idisc.ne.0)then
reflat=hyplat
reflon=hyplon
refdep=hypdep
c
farea=4.d0*fsize**2
patchsize=dmax1(VSREF/(4.d0*fcut2),dsqrt(farea/dble(NPATCH)))
c
if(hypdep/dsin(dip*DEG2RAD).gt.fsize)then
wid1=dmin1(WIDTHMAX,fsize)
else
wid1=hypdep/dsin(dip*DEG2RAD)
endif
wid2=dmin1(WIDTHMAX,fsize,(deps(ng)-hypdep)/dsin(dip*DEG2RAD))
c
len1=dmin1(LENGTHMAX,0.5d0*farea/(wid1+wid2))
len2=len1
c
length=len1+len2
width=wid1+wid2
nlen=max0(1,idnint(length/patchsize))
nwid=max0(1,idnint(width/patchsize))
nsf=nlen*nwid
endif
c
c allocation of sub-fault parameters
c
allocate(sflat(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sflat not allocated!'
allocate(sflon(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sflon not allocated!'
allocate(sfdep(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sfdep not allocated!'
allocate(sflen(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sflen not allocated!'
allocate(sfwid(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sfwid not allocated!'
allocate(sfxln(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sfxln not allocated!'
allocate(sfywd(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sfywd not allocated!'
allocate(sfmue(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sfmue not allocated!'
allocate(sfro(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sfro not allocated!'
allocate(sfvp(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sfvp not allocated!'
allocate(sfvs(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sfvs not allocated!'
allocate(sfstk(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sfstk not allocated!'
allocate(sfdip(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sfdip not allocated!'
allocate(sfrak(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sfrak not allocated!'
c
if(idisc.eq.0)then
open(20,file=finitefault,status='old')
read(20,'(a1)')text
patchsize=0.d0
isf=0
do i=1,nsf
isf=isf+1
read(20,*)sflat(isf),sflon(isf),sfdep(i),sflen(isf),
& sfwid(isf),sfstk(isf),sfdip(isf),sfrak(isf)
sfdep(isf)=sfdep(isf)*KM2M
sflen(isf)=sflen(isf)*KM2M
sfwid(isf)=sfwid(isf)*KM2M
if(iref.ge.1.and.iref.le.4)then
st=sfstk(isf)*DEG2RAD
di=sfdip(isf)*DEG2RAD
if(iref.eq.1)then
sfdep(isf)=sfdep(isf)+0.5d0*sfwid(isf)*dsin(di)
pn=0.5d0*sflen(isf)*dcos(st)
& -0.5d0*sfwid(isf)*dcos(di)*dsin(st)
pe=0.5d0*sflen(isf)*dsin(st)
& +0.5d0*sfwid(isf)*dcos(di)*dcos(st)
else if(iref.eq.2)then
sfdep(isf)=sfdep(isf)+0.5d0*sfwid(isf)*dsin(di)
pn=-0.5d0*sflen(isf)*dcos(st)
& -0.5d0*sfwid(isf)*dcos(di)*dsin(st)
pe=-0.5d0*sflen(isf)*dsin(st)
& +0.5d0*sfwid(isf)*dcos(di)*dcos(st)
else if(iref.eq.3)then
sfdep(isf)=sfdep(isf)-0.5d0*sfwid(isf)*dsin(di)
pn=0.5d0*sflen(isf)*dcos(st)
& +0.5d0*sfwid(isf)*dcos(di)*dsin(st)
pe=0.5d0*sflen(isf)*dsin(st)
& -0.5d0*sfwid(isf)*dcos(di)*dcos(st)
else if(iref.eq.4)then
sfdep(isf)=sfdep(isf)-0.5d0*sfwid(isf)*dsin(di)
pn=-0.5d0*sflen(isf)*dcos(st)
& +0.5d0*sfwid(isf)*dcos(di)*dsin(st)
pe=-0.5d0*sflen(isf)*dsin(st)
& -0.5d0*sfwid(isf)*dcos(di)*dcos(st)
endif
c
c determine central point of the subfault
c
c spherical triangle:
c A = pole, B = source position, C = reference position
c
sma=dsqrt(pn**2+pe**2)/REARTH
smb=0.5d0*PI-sflat(isf)*DEG2RAD
bgc=datan2(pe,pn)
smc=dacos(dcos(sma)*dcos(smb)
& +dsin(sma)*dsin(smb)*dcos(bgc))
bga=dasin(dsin(sma)*dsin(bgc)/dsin(smc))
c
c geographic coordinate of the equivalent point source
c
sflat(isf)=90.d0-smc/DEG2RAD
sflon(isf)=dmod(sflon(isf)+bga/DEG2RAD,360.d0)
endif
c
call disazi(REARTH,hyplat,hyplon,
& sflat(isf),sflon(isf),rn,re)
d2=dsqrt(rn**2+re**2+(sfdep(isf)-hypdep)**2)
patchsize=patchsize+sflen(isf)*sfwid(isf)
enddo
nsf=isf
patchsize=dsqrt(patchsize/dble(nsf))
close(20)
else
st=strike*DEG2RAD
di=dip*DEG2RAD
x0=-len1
dx=length/dble(nlen)
y0=-wid1
dy=width/dble(nwid)
isf=0
do i=1,nlen
x=x0+(dble(i-1)+0.5d0)*dx
do j=1,nwid
y=y0+(dble(j-1)+0.5d0)*dy
isf=isf+1
sflen(isf)=dx
sfwid(isf)=dy
sfstk(isf)=strike
sfdip(isf)=dip
sfrak(isf)=rake
sfdep(isf)=refdep+y*dsin(di)
c
c determine central point of the subfault
c
c spherical triangle:
c A = pole, B = source position, C = reference position
c
pn=x*dcos(st)-y*dcos(di)*dsin(st)
pe=x*dsin(st)+y*dcos(di)*dcos(st)
sma=dsqrt(pn**2+pe**2)/REARTH
smb=0.5d0*PI-reflat*DEG2RAD
bgc=datan2(pe,pn)
smc=dacos(dcos(sma)*dcos(smb)
& +dsin(sma)*dsin(smb)*dcos(bgc))
bga=dasin(dsin(sma)*dsin(bgc)/dsin(smc))
c
c geographic coordinate of the equivalent point source
c
sflat(isf)=90.d0-smc/DEG2RAD
sflon(isf)=dmod(reflon+bga/DEG2RAD,360.d0)
enddo
enddo
nsf=isf
endif
c
isfhyp=0
d2min=0.d0
do isf=1,nsf
call disazi(REARTH,hyplat,hyplon,
& sflat(isf),sflon(isf),rn,re)
d2=re*re+rn*rn+(sfdep(isf)-hypdep)**2
if(isf.eq.1)then
isfhyp=1
d2min=d2
else if(d2min.gt.d2)then
isfhyp=isf
d2min=d2
endif
enddo
write(*,'(a,i10)')' total number of discrete fault patches: ',nsf
write(*,'(a,i10,4(a,f8.2))')' hypocentre-nearest patch number: ',
& isfhyp,' at (',sflat(isfhyp),' deg_N, ',sflon(isfhyp),
& ' deg_E, ',sfdep(isfhyp)/KM2M,' km)'
c
c
c local coordinates of patches
c
st=sfstk(isfhyp)*DEG2RAD
di=sfdip(isfhyp)*DEG2RAD
rstk(1)=dcos(st)
rstk(2)=dsin(st)
rstk(3)=0.d0
rddp(1)=-dcos(di)*dsin(st)
rddp(2)=dcos(di)*dcos(st)
rddp(3)=dsin(di)
do isf=1,nsf
call disazi(REARTH,hyplat,hyplon,
& sflat(isf),sflon(isf),rn,re)
sfxln(isf)=rstk(1)*rn+rstk(2)*re
sfywd(isf)=rddp(1)*rn+rddp(2)*re+rddp(3)*(sfdep(isf)-hypdep)
enddo
c
allocate(it1sf(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: it1sf not allocated!'
allocate(it2sf(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: it2sf not allocated!'
allocate (dis3dsf(nsf,nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: dis3dsf not allocated!'
allocate (azisf2hyp(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: azisf2hyp not allocated!'
allocate (plgsf2hyp(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: plgsf2hyp not allocated!'
c
c Inter-fault-patch 3D distances
c
do isf=1,nsf
do jsf=1,isf-1
call disazi(REARTH,sflat(isf),sflon(isf),
& sflat(jsf),sflon(jsf),rn,re)
dis3dsf(isf,jsf)=dsqrt(rn*rn+re*re
& +(sfdep(jsf)-sfdep(isf))**2)
enddo
dis3dsf(isf,isf)=0.d0
enddo
do isf=1,nsf
do jsf=isf+1,nsf
dis3dsf(isf,jsf)=dis3dsf(jsf,isf)
enddo
enddo
c
do isf=1,nsf
if(isf.ne.isfhyp)then
call disazi(REARTH,sflat(isfhyp),sflon(isfhyp),
& sflat(isf),sflon(isf),rn,re)
azisf2hyp(isf)=dmod(360.d0+datan2(re,rn)/DEG2RAD,360.d0)
plgsf2hyp(isf)=dmod(180.d0+datan2(dsqrt(rn*rn+re*re),
& sfdep(isf)-sfdep(isfhyp))/DEG2RAD,180.d0)
else
azisf2hyp(isf)=0.d0
plgsf2hyp(isf)=0.d0
endif
enddo
c
do isf=1,nsf
do i=2,nzmod
if(sfdep(isf).le.zpmod(i))goto 300
enddo
300 i=i-1
sfvp(isf)=vpmod(i)
sfvs(isf)=vsmod(i)
sfro(isf)=romod(i)
sfmue(isf)=mumod(i)
enddo
c
vsmean=0.d0
do isf=1,nsf
vsmean=vsmean+sflen(isf)*sfwid(isf)*sfvs(isf)
enddo
vsmean=vsmean/(dble(nsf)*patchsize**2)
c
dt=dmax1(dtgrn,twindow/dble(NTMAX-1))
c
nt=1
400 nt=2*nt
if(dble(nt-1)*dt.lt.twindow)goto 400
if(nt.lt.NTMIN)then
nt=NTMIN
else if(nt.gt.NTMAX)then
nt=NTMAX
endif
dt=twindow/dble(nt-1)
nf=nt/2
df=1.d0/(dble(nt)*dt)
c
allocate(lpfs(nf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: lpfs not allocated!'
allocate(omi(nf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: omi not allocated!'
c
c allocate working space
c
allocate (cfct(2*nf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: cfct not allocated!'
allocate (dfct(4*nf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: dfct not allocated!'
c
allocate (sstf(nf,nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sstf not allocated!'
allocate (dstf(nf,nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: dstf not allocated!'
allocate (stfswap(nf,nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: stfswap not allocated!'
c
allocate(estf(2*nf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: estf not allocated!'
allocate(pstf(2*nf,nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: pstf not allocated!'
allocate (sftr(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sftr not allocated!'
allocate (sfslp(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sfslp not allocated!'
allocate (sfswap(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sfswap not allocated!'
allocate (subf(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: subf not allocated!'
allocate (isfnb(nsf,nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: isfnb not allocated!'
allocate (nsfnb(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: nsfnb not allocated!'
c
allocate (slpswap(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: slpswap not allocated!'
allocate (rsmth(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: rsmooth not allocated!'
c
c Butterworth bandpass filter
c
do i=1,nf
omi(i)=PI2*dble(i-1)*df
enddo
c
call butterworth(6,sfvs(isfhyp)/patchsize,df,nf,lpfs)
do i=1,nf
lpfs(i)=lpfs(i)*dconjg(lpfs(i))
enddo
c
return
subroutine idsmodel(ierr)
use idsalloc
implicit none
integer*4 ierr
c
integer*4 i,j,ind,ig,ir,isf,jsf,iztp,nzmod,nsfred
real*8 pe,pn,sma,smb,smc,bga,bgc,re,rn,d2,d2min
real*8 st,di,t1,t2,dis1,dis2,farea,fsize,pswid
real*8 x,y,x0,y0,dx,dy
real*8 rstk(3),rddp(3)
logical*2 head,diff
character*180 text
c
do ind=120,1,-1
if(smgrndir(ind:ind).ne.' ')goto 100
enddo
100 continue
if(ind.ge.120)then
stop ' Error in idsmodel: too long folder name of smgrndir!'
endif
if(smgrndir(ind:ind).ne.'/'.or.
& smgrndir(ind:ind).ne.'\')then
smgrndir=smgrndir(1:ind)//'/'
ind=ind+1
endif
open(10,file=smgrndir(1:ind)//'GreenInfo.dat',status='old')
call skipdoc(10)
read(10,*)twgrn,dtgrn,ntgrn
if(ntgrn.le.0)stop ' Error in idsmodel: wrong value for ntgrn!'
c
call skipdoc(10)
read(10,*)nr
if(nr.lt.2)then
stop ' Error in idsmodel: GF distance coverage insufficient!'
endif
c
allocate (r(nr),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: r not allocated!'
c
read(10,*)(r(i),i=1,nr)
do i=1,nr
r(i)=r(i)*KM2M
enddo
c
call skipdoc(10)
read(10,*)ng
if(ng.lt.2)then
stop ' Error in idsmodel: GF depth coverage insufficient!'
endif
c
allocate (deps(ng),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: deps not allocated!'
allocate (smgrnfile(ng),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: smgrnfile not allocated!'
c
read(10,*)(deps(ig),smgrnfile(ig),ig=1,ng)
do ig=1,ng
smgrnfile(ig)=smgrndir(1:ind)//smgrnfile(ig)
deps(ig)=deps(ig)*KM2M
enddo
c
c Layered earth model
c
call skipdoc(10)
read(10,*)nlayer
c
allocate(zpmod(nlayer),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: zpmod not allocated!'
allocate(hpmod(nlayer),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: hpmod not allocated!'
allocate(vpmod(nlayer),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: vpmod not allocated!'
allocate(vsmod(nlayer),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: vsmod not allocated!'
allocate(romod(nlayer),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: romod not allocated!'
allocate(mumod(nlayer),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: mumod not allocated!'
c
call skipdoc(10)
do i=1,nlayer
read(10,*)j,hpmod(i),vpmod(i),vsmod(i),romod(i)
hpmod(i)=hpmod(i)*KM2M
vpmod(i)=vpmod(i)*KM2M
vsmod(i)=vsmod(i)*KM2M
romod(i)=romod(i)*KM2M
mumod(i)=romod(i)*vsmod(i)**2
enddo
zpmod(1)=0.d0
nzmod=1
do i=2,nlayer
zpmod(i)=zpmod(i-1)+hpmod(i-1)
if(zpmod(i).le.hypdep+stdismax)then
nzmod=i
endif
enddo
c
call skipdoc(10)
read(10,*)tptable
tptable=smgrndir(1:ind)//tptable
call skipdoc(10)
read(10,*)tstable
tstable=smgrndir(1:ind)//tstable
c
close(10)
c
allocate (tp(nr,ng),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: tp not allocated!'
allocate(tkftp(nr,ng),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: tkftp not allocated!'
allocate(slwtp(nr,ng),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: slwtp not allocated!'
allocate (ts(nr,ng),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: ts not allocated!'
allocate(tkfts(nr,ng),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: tkfts not allocated!'
allocate(slwts(nr,ng),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: slwts not allocated!'
c
open(11,file=tptable,form='unformatted',status='old')
read(11)((tp(ir,ig),tkftp(ir,ig),slwtp(ir,ig),ir=1,nr),ig=1,ng)
close(11)
open(12,file=tstable,form='unformatted',status='old')
read(12)((ts(ir,ig),tkfts(ir,ig),slwts(ir,ig),ir=1,nr),ig=1,ng)
close(12)
c
c empirical scaling laws
c
trise=10.d0**(0.5d0*mwini-2.32d0)
fsize=10.d0**(0.5d0*mwini+1.16d0)
c
fcut2=dmax1(FCUTUP,dmin1(0.75d0*fcut2g,5.0d0/trise))
fcut1=dmax1(FCUTLW,0.5d0/trise)
c
tsource=2.0d0*trise
tpre=dmax1(TPREMIN,fsize/VPREF)
tpst=dmax1(TPSTMIN,fsize/VSREF)
ttap=dmax1(2.5d0/fcut2g,0.25d0*trise)
c
swindow=tsource+dmax1(trise,tpst)
twindow=swindow+tpre+tpst
c
if(idisc.ne.0)then
reflat=hyplat
reflon=hyplon
refdep=hypdep
c
farea=4.d0*fsize**2
patchsize=dmax1(VSREF/(4.d0*fcut2),dsqrt(farea/dble(NPATCH)))
c
if(hypdep/dsin(dip*DEG2RAD).gt.fsize)then
wid1=dmin1(WIDTHMAX,fsize)
else
wid1=hypdep/dsin(dip*DEG2RAD)
endif
wid2=dmin1(WIDTHMAX,fsize,(deps(ng)-hypdep)/dsin(dip*DEG2RAD))
c
len1=dmin1(LENGTHMAX,0.5d0*farea/(wid1+wid2))
len2=len1
c
length=len1+len2
width=wid1+wid2
nlen=max0(1,idnint(length/patchsize))
nwid=max0(1,idnint(width/patchsize))
nsf=nlen*nwid
endif
c
c allocation of sub-fault parameters
c
allocate(sflat(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sflat not allocated!'
allocate(sflon(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sflon not allocated!'
allocate(sfdep(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sfdep not allocated!'
allocate(sflen(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sflen not allocated!'
allocate(sfwid(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sfwid not allocated!'
allocate(sfxln(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sfxln not allocated!'
allocate(sfywd(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sfywd not allocated!'
allocate(sfmue(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sfmue not allocated!'
allocate(sfro(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sfro not allocated!'
allocate(sfvp(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sfvp not allocated!'
allocate(sfvs(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sfvs not allocated!'
allocate(sfstk(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sfstk not allocated!'
allocate(sfdip(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sfdip not allocated!'
allocate(sfrak(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sfrak not allocated!'
c
if(idisc.eq.0)then
open(20,file=finitefault,status='old')
read(20,'(a1)')text
patchsize=0.d0
isf=0
do i=1,nsf
isf=isf+1
read(20,*)sflat(isf),sflon(isf),sfdep(i),sflen(isf),
& sfwid(isf),sfstk(isf),sfdip(isf),sfrak(isf)
sfdep(isf)=sfdep(isf)*KM2M
sflen(isf)=sflen(isf)*KM2M
sfwid(isf)=sfwid(isf)*KM2M
if(iref.ge.1.and.iref.le.4)then
st=sfstk(isf)*DEG2RAD
di=sfdip(isf)*DEG2RAD
if(iref.eq.1)then
sfdep(isf)=sfdep(isf)+0.5d0*sfwid(isf)*dsin(di)
pn=0.5d0*sflen(isf)*dcos(st)
& -0.5d0*sfwid(isf)*dcos(di)*dsin(st)
pe=0.5d0*sflen(isf)*dsin(st)
& +0.5d0*sfwid(isf)*dcos(di)*dcos(st)
else if(iref.eq.2)then
sfdep(isf)=sfdep(isf)+0.5d0*sfwid(isf)*dsin(di)
pn=-0.5d0*sflen(isf)*dcos(st)
& -0.5d0*sfwid(isf)*dcos(di)*dsin(st)
pe=-0.5d0*sflen(isf)*dsin(st)
& +0.5d0*sfwid(isf)*dcos(di)*dcos(st)
else if(iref.eq.3)then
sfdep(isf)=sfdep(isf)-0.5d0*sfwid(isf)*dsin(di)
pn=0.5d0*sflen(isf)*dcos(st)
& +0.5d0*sfwid(isf)*dcos(di)*dsin(st)
pe=0.5d0*sflen(isf)*dsin(st)
& -0.5d0*sfwid(isf)*dcos(di)*dcos(st)
else if(iref.eq.4)then
sfdep(isf)=sfdep(isf)-0.5d0*sfwid(isf)*dsin(di)
pn=-0.5d0*sflen(isf)*dcos(st)
& +0.5d0*sfwid(isf)*dcos(di)*dsin(st)
pe=-0.5d0*sflen(isf)*dsin(st)
& -0.5d0*sfwid(isf)*dcos(di)*dcos(st)
endif
c
c determine central point of the subfault
c
c spherical triangle:
c A = pole, B = source position, C = reference position
c
sma=dsqrt(pn**2+pe**2)/REARTH
smb=0.5d0*PI-sflat(isf)*DEG2RAD
bgc=datan2(pe,pn)
smc=dacos(dcos(sma)*dcos(smb)
& +dsin(sma)*dsin(smb)*dcos(bgc))
bga=dasin(dsin(sma)*dsin(bgc)/dsin(smc))
c
c geographic coordinate of the equivalent point source
c
sflat(isf)=90.d0-smc/DEG2RAD
sflon(isf)=dmod(sflon(isf)+bga/DEG2RAD,360.d0)
endif
c
call disazi(REARTH,hyplat,hyplon,
& sflat(isf),sflon(isf),rn,re)
d2=dsqrt(rn**2+re**2+(sfdep(isf)-hypdep)**2)
patchsize=patchsize+sflen(isf)*sfwid(isf)
enddo
nsf=isf
patchsize=dsqrt(patchsize/dble(nsf))
close(20)
else
st=strike*DEG2RAD
di=dip*DEG2RAD
x0=-len1
dx=length/dble(nlen)
y0=-wid1
dy=width/dble(nwid)
isf=0
do i=1,nlen
x=x0+(dble(i-1)+0.5d0)*dx
do j=1,nwid
y=y0+(dble(j-1)+0.5d0)*dy
isf=isf+1
sflen(isf)=dx
sfwid(isf)=dy
sfstk(isf)=strike
sfdip(isf)=dip
sfrak(isf)=rake
sfdep(isf)=refdep+y*dsin(di)
c
c determine central point of the subfault
c
c spherical triangle:
c A = pole, B = source position, C = reference position
c
pn=x*dcos(st)-y*dcos(di)*dsin(st)
pe=x*dsin(st)+y*dcos(di)*dcos(st)
sma=dsqrt(pn**2+pe**2)/REARTH
smb=0.5d0*PI-reflat*DEG2RAD
bgc=datan2(pe,pn)
smc=dacos(dcos(sma)*dcos(smb)
& +dsin(sma)*dsin(smb)*dcos(bgc))
bga=dasin(dsin(sma)*dsin(bgc)/dsin(smc))
c
c geographic coordinate of the equivalent point source
c
sflat(isf)=90.d0-smc/DEG2RAD
sflon(isf)=dmod(reflon+bga/DEG2RAD,360.d0)
enddo
enddo
nsf=isf
endif
c
isfhyp=0
d2min=0.d0
do isf=1,nsf
call disazi(REARTH,hyplat,hyplon,
& sflat(isf),sflon(isf),rn,re)
d2=re*re+rn*rn+(sfdep(isf)-hypdep)**2
if(isf.eq.1)then
isfhyp=1
d2min=d2
else if(d2min.gt.d2)then
isfhyp=isf
d2min=d2
endif
enddo
write(*,'(a,i10)')' total number of discrete fault patches: ',nsf
write(*,'(a,i10,4(a,f8.2))')' hypocentre-nearest patch number: ',
& isfhyp,' at (',sflat(isfhyp),' deg_N, ',sflon(isfhyp),
& ' deg_E, ',sfdep(isfhyp)/KM2M,' km)'
c
c
c local coordinates of patches
c
st=sfstk(isfhyp)*DEG2RAD
di=sfdip(isfhyp)*DEG2RAD
rstk(1)=dcos(st)
rstk(2)=dsin(st)
rstk(3)=0.d0
rddp(1)=-dcos(di)*dsin(st)
rddp(2)=dcos(di)*dcos(st)
rddp(3)=dsin(di)
do isf=1,nsf
call disazi(REARTH,hyplat,hyplon,
& sflat(isf),sflon(isf),rn,re)
sfxln(isf)=rstk(1)*rn+rstk(2)*re
sfywd(isf)=rddp(1)*rn+rddp(2)*re+rddp(3)*(sfdep(isf)-hypdep)
enddo
c
allocate(it1sf(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: it1sf not allocated!'
allocate(it2sf(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: it2sf not allocated!'
allocate (dis3dsf(nsf,nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: dis3dsf not allocated!'
allocate (azisf2hyp(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: azisf2hyp not allocated!'
allocate (plgsf2hyp(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: plgsf2hyp not allocated!'
c
c Inter-fault-patch 3D distances
c
do isf=1,nsf
do jsf=1,isf-1
call disazi(REARTH,sflat(isf),sflon(isf),
& sflat(jsf),sflon(jsf),rn,re)
dis3dsf(isf,jsf)=dsqrt(rn*rn+re*re
& +(sfdep(jsf)-sfdep(isf))**2)
enddo
dis3dsf(isf,isf)=0.d0
enddo
do isf=1,nsf
do jsf=isf+1,nsf
dis3dsf(isf,jsf)=dis3dsf(jsf,isf)
enddo
enddo
c
do isf=1,nsf
if(isf.ne.isfhyp)then
call disazi(REARTH,sflat(isfhyp),sflon(isfhyp),
& sflat(isf),sflon(isf),rn,re)
azisf2hyp(isf)=dmod(360.d0+datan2(re,rn)/DEG2RAD,360.d0)
plgsf2hyp(isf)=dmod(180.d0+datan2(dsqrt(rn*rn+re*re),
& sfdep(isf)-sfdep(isfhyp))/DEG2RAD,180.d0)
else
azisf2hyp(isf)=0.d0
plgsf2hyp(isf)=0.d0
endif
enddo
c
do isf=1,nsf
do i=2,nzmod
if(sfdep(isf).le.zpmod(i))goto 300
enddo
300 i=i-1
sfvp(isf)=vpmod(i)
sfvs(isf)=vsmod(i)
sfro(isf)=romod(i)
sfmue(isf)=mumod(i)
enddo
c
vsmean=0.d0
do isf=1,nsf
vsmean=vsmean+sflen(isf)*sfwid(isf)*sfvs(isf)
enddo
vsmean=vsmean/(dble(nsf)*patchsize**2)
c
dt=dmax1(dtgrn,twindow/dble(NTMAX-1))
c
nt=1
400 nt=2*nt
if(dble(nt-1)*dt.lt.twindow)goto 400
if(nt.lt.NTMIN)then
nt=NTMIN
else if(nt.gt.NTMAX)then
nt=NTMAX
endif
dt=twindow/dble(nt-1)
nf=nt/2
df=1.d0/(dble(nt)*dt)
c
allocate(lpfs(nf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: lpfs not allocated!'
allocate(omi(nf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: omi not allocated!'
c
c allocate working space
c
allocate (cfct(2*nf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: cfct not allocated!'
allocate (dfct(4*nf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: dfct not allocated!'
c
allocate (sstf(nf,nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sstf not allocated!'
allocate (dstf(nf,nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: dstf not allocated!'
allocate (stfswap(nf,nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: stfswap not allocated!'
c
allocate(estf(2*nf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: estf not allocated!'
allocate(pstf(2*nf,nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: pstf not allocated!'
allocate (sftr(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sftr not allocated!'
allocate (sfslp(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sfslp not allocated!'
allocate (sfswap(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: sfswap not allocated!'
allocate (subf(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: subf not allocated!'
allocate (isfnb(nsf,nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: isfnb not allocated!'
allocate (nsfnb(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: nsfnb not allocated!'
c
allocate (slpswap(nsf),stat=ierr)
if(ierr.ne.0)stop ' Error in idsmodel: slpswap not allocated!'