[cig-commits] r4496 - in seismo/1D/syndat/trunk: . fdb
barmin at geodynamics.org
barmin at geodynamics.org
Fri Sep 8 10:13:42 PDT 2006
Author: barmin
Date: 2006-09-08 10:13:41 -0700 (Fri, 08 Sep 2006)
New Revision: 4496
Added:
seismo/1D/syndat/trunk/fdb/
seismo/1D/syndat/trunk/fdb/fdb_eigen.f
seismo/1D/syndat/trunk/green.h
Modified:
seismo/1D/syndat/trunk/Makefile.am
seismo/1D/syndat/trunk/eigcon.f
seismo/1D/syndat/trunk/green.f
Log:
Upgraded eigcon and green programs. They work with a new format
eigen functions flat dbase.
Modified: seismo/1D/syndat/trunk/Makefile.am
===================================================================
--- seismo/1D/syndat/trunk/Makefile.am 2006-09-08 15:53:55 UTC (rev 4495)
+++ seismo/1D/syndat/trunk/Makefile.am 2006-09-08 17:13:41 UTC (rev 4496)
@@ -26,9 +26,14 @@
# main programs
minos_bran_SOURCES = minos_bran.f
syndat_SOURCES = syndat.f $(gfs) fftsubs.f
-green_SOURCES = green.f $(gfs) fftsubs.f rcode.f
-eigcon_SOURCES = eigcon.f $(gfs)
+green_SOURCES = green.f $(fdb) $(gfs) fftsubs.f rcode.f
+eigcon_SOURCES = eigcon.f $(fdb)
+# FDB utils
+
+fdb = \
+ fdb/fdb_eigen.f
+
# GFS utils
ghead_SOURCES = gfs/ghead.f $(gfs)
glist_SOURCES = gfs/glist.f $(gfs)
Modified: seismo/1D/syndat/trunk/eigcon.f
===================================================================
--- seismo/1D/syndat/trunk/eigcon.f 2006-09-08 15:53:55 UTC (rev 4495)
+++ seismo/1D/syndat/trunk/eigcon.f 2006-09-08 17:13:41 UTC (rev 4496)
@@ -12,30 +12,75 @@
* variable irecl for reclength of eigfun file*
* latest change: 19.08.91 subroutines for sph, rad and tor modes *
***************************************************************************
- dimension ititle(20),hed6(1554)
- real*8 r(300),pi
- character atyp*4, model*8,filen*80
- common/hedX/nhed,model,jcom,nz,atyp,lz,w,q,g,nrad,
- * rn,wn,accn,rout(1540)
- equivalence (nhed,hed6(1))
+ implicit none
+ integer*4 mk
+ parameter (mk=350)
+c --- eigen relation common block
+ real*4 per_eigen,phvel_eigen,grvel_eigen,attn_eigen
+ integer*4 norder_eigen,lorder_eigen,eigid_eigen,
+ + nraw_eigen,ncol_eigen,npar_eigen,foff_eigen,
+ + commid_eigen
+ character*2 datatype_eigen
+ character*64 dir_eigen
+ character*32 dfile_eigen
+ character*17 lddate_eigen
+ character*1 typeo_eigen
+ common/c_eigen/norder_eigen,lorder_eigen,
+ + eigid_eigen,per_eigen,phvel_eigen,grvel_eigen,
+ + attn_eigen,nraw_eigen,ncol_eigen,npar_eigen,
+ + foff_eigen,commid_eigen,typeo_eigen,
+ + datatype_eigen,dir_eigen,dfile_eigen,lddate_eigen
+c --- other variables
+ real*8 r(mk),pi
+ real*4 rout(mk),buf(6,mk)
+ real*4 bigg,tau,rhobar,con,accn
+ real*4 tref,rx,r1,r2,dr,dum,rn,gn,vn,wn
+ real*4 dmax,rmax,ww,qq,qc,tf,fl,fl1,fl3
+ integer*4 i,ieig,idat,iin,ifanis,ifdeck,j,jj,jcom,l,ll,lll
+ integer*4 n,nic,noc,nreg,nn,nlay,nstart,nrad,ni,nrecl,nrest
+ integer*4 ierr,iflag,ititle(20)
+ character*20 str
+ character*64 dir
+ character*1 typeo
+ character*256 fmodel,fflatin,fbinin,fout
+c ---
data bigg,tau,rhobar/6.6723e-11,1000.0,5515.0/
- data rout/1540*0./
- irecl=5600
- iin=11
- pi=4.*atan(1.d0)
+ pi=4.0d0*datan(1.d0)
con=pi*bigg
- print*,'enter name of model : '
- read (*,'(a)') model
-c... hed6 must have dimension of rout+14 (not 13, because of model*8!)
-c jcom : flag to indicate mode type
- print*,' spheroidals (3) or toroidals (2) or radial (1)'
- read(*,*)jcom
+ nstart = 0
c
-c read in radial knots
+c read control parameters from stdin
c
+c read jcom : flag to indicate mode type
+ print*,' spheroidals (3) or toroidals (2) or radial (1) or'
+ print*,' inner core toroidals (4) modes'
+ read(*,*)jcom
+c read model file name
print*,' enter name of model file'
- read(*,'(a)')filen
- open(iin,file=filen,status='old')
+ read(*,'(a)')fmodel
+c read max depth for eigenvectors
+ print *,'enter max depth [km] : '
+ read(5,*) dmax
+c read minos_bran output text file
+ write(*,*) ' enter name of minos_bran output text file'
+ read(*,'(a)') fflatin
+c read minos_bran output binary unformatted file
+ write(*,*) ' minos_bran output binary unformatted file name'
+ read(*,'(a)') fbinin
+c read dbase_name. There might be two form of input string:
+c path/dbase_name or dbase_name, where path is relative or absolute
+c path to existing directory. Under this directory or working
+c directory (no path/ in input string) will be created db relation file
+c dbase_name.eigen, directory dbase_name.eigen.dat, and
+c under dbase_name.eigen.dat binary data file eigen.
+ print*,
+ * ' enter path/dbase_name or dbase_name to store eigenfunctions:'
+ read(*,'(a)')fout
+c
+c read in radial knots from model
+c
+ iin = 7
+ open(iin,file=fmodel,status='old')
read(iin,101) (ititle(i),i=1,20)
101 format(20a4)
read(iin,*) ifanis,tref,ifdeck
@@ -64,8 +109,12 @@
2000 close(iin)
c
c rn : radius at surface
+c wn : frequency normalization
+c vn : velocity normalisation
+c gn : acceleration normalisation
+c accn : ??????
c n : index of surface grid point
-c nstart : " " lowest " "
+c nstart : index of lowest grid point
c nrad : # of gridpoints of reduced eigenfunctions
c
rn=r(n)
@@ -74,207 +123,121 @@
wn=vn/rn
accn=1.e+20/(rhobar*rn**4)
c normalize radius
- do 55 i=1,n
+ do i=1,n
r(i)=r(i)/rn
if(i.gt.1.and.dabs(r(i)-r(i-1)).lt.1.d-7) r(i)=r(i-1)
- 55 continue
- print *,'enter max depth [km] : '
- read(5,*) dmax
+ enddo
+c cut radius knots lower than max depth
rmax=1.-dmax*tau/rn
j=0
- do 25 i=1,n
- if(r(i).gt.rmax) goto 30
- 25 j=j+1
- 30 nstart=j
+ do i=1,n
+ j=j+1
+ if(r(i).ge.rmax) goto 30
+ enddo
+ 30 nstart=max0(j-1,1)
j=0
- do 35 i=nstart,n
+ do i=nstart,n
j=j+1
- 35 rout(j)=r(i)
+ rout(j)=r(i)
+ enddo
nrad=j
- print *,'n,nstart,nrad = ',n,nstart,nrad
-
- ityp = 6
- print*,' enter name of eigenfunction file to be reduced'
- read(*,'(a)')filen
- open(2,file=filen,
- * form='unformatted',iostat=iret)
- call gfs_opena(3,'output gfs-file name :',ityp,ierr)
-c write first header containing radial grid
- kfl = 1
- nz = 0
- atyp = 'grid'
- lz = 0
- w = 0.
- q = 0.
- g = 0.
- nhed = 0
-ccc
- print 40,nhed,model,jcom,nz,atyp,lz,w,q,g,nrad,
- * rn,wn,accn
- 40 format(i6,1x,a8,2i6,1x,a4,i4,3g15.5,i6,3g15.5)
- print 41,(rout(i),i=1,nrad)
- 41 format(5g15.5)
-ccc
- call gfs_rwdir(3,hed6,kfl,'w',ierr)
-c... i/o of eigenfunctions
- if(jcom.eq.3)then
- call conv4(nstart,n,kfl)
- else if(jcom.eq.2) then
- call conv2(nstart,n,kfl)
+ write(*,*) 'n,nstart,nrad = ',n,nstart,nrad
+c open minos_bran plane output file and search mode part
+ open(7,file=fflatin,form='formatted',status='old')
+ 1 read(7,'(a)',end=9) str
+ ni=0
+ do i = 1,20
+ if(str(i:i).eq.'m') ni=i
+ enddo
+ if(ni.eq.0) goto 1
+ if(str(ni:ni+3).eq.'mode') goto 2
+ goto 1
+ 9 stop 'Wrong minos_bran output text file'
+ 2 read(7,'(a)') str
+c open minos_bran binary unformatted file
+ open(8,file=fbinin,form='unformatted',status='old')
+c***************************************************************
+c Create .eigen relations and binary eigenfunctions file
+c***************************************************************
+ ieig = 9
+ idat = 10
+ call null_eigen
+ eigid_eigen = 1
+ foff_eigen = 0
+ ncol_eigen = 2
+ if(jcom.eq.3) ncol_eigen = 6
+ npar_eigen = 4
+ nraw_eigen = nrad
+ nrecl = max0(ncol_eigen*nraw_eigen+npar_eigen,10)*4
+ call open_eigen(fout,nrecl,ieig,idat,dir,'w',ierr)
+ dir_eigen = dir
+ dfile_eigen = 'eigen'
+ call write_eigen(ieig,ierr)
+ nrest=nrecl/4-nrad-9
+ write(idat,rec=eigid_eigen)
+ + jcom,bigg,tau,rhobar,rn,wn,vn,gn,accn,
+ + (rout(l),l=1,nrad),(foff_eigen,lll=1,nrest)
+c
+c Main loop ----
+c
+ 200 eigid_eigen = eigid_eigen+1
+ foff_eigen = foff_eigen+nrecl
+ read(8,end=99) nn,ll,ww,qq,qc,((buf(l,lll),lll=1,n),
+ + l=1,ncol_eigen)
+ read(7,*) norder_eigen,typeo,lorder_eigen,phvel_eigen,
+ + tf,per_eigen,grvel_eigen,attn_eigen
+ if(nn.ne.norder_eigen.or.ll.ne.lorder_eigen) then
+ write(*,*)
+ + 'ERR001: eigcon: Input plane and binary files differ: ',
+ + nn,ll,norder_eigen,lorder_eigen
+ stop
+ endif
+c check that jcom corresponds to minos_bran mode
+ iflag = 0
+ if(jcom.eq.4) then
+ if(typeo.ne.'c') iflag = 1
+ typeo_eigen = 'C'
+ else if(jcom.eq.3)then
+ if(typeo.ne.'s') iflag = 1
+ typeo_eigen = 'S'
+ else if(jcom.eq.2) then
+ if(typeo.ne.'t') iflag = 1
+ typeo_eigen = 'T'
+ else if(jcom.eq.1) then
+ if(typeo.ne.'s') iflag = 1
+ typeo_eigen = 'R'
+ grvel_eigen = -1.0
else
- call conv1(nstart,n,kfl)
+ write(*,*) 'ERR002: eigcon: Unknown jcom ',jcom
+ stop
endif
- close(2)
- call gfs_close(3)
- stop
- end
-
- subroutine conv4(nstart,n,kfl)
- parameter(mnl=323)
- common/bufz/nn,ll,ww,qq,cg,buf(6*mnl)
- common/hedX/nhed,model,jcom,nz,atyp,lz,w,q,g,nrad,
- * rn,wn,accn,x(1540)
- dimension abuf(6*mnl+5),a(4,350),aa(1),hed6(1554)
- character*4 atyp
- character*8 model
- equivalence (nn,abuf),(a,aa),(hed6,nhed)
- atyp ='S '
- nvec = 6*n + 5
- nloop=4
- nlen = nloop*nrad
- nout = nloop*nrad+5
- print *,'nlen (must be .le. 1540 or else change defhdr) : ',nlen,nvec
- if(nlen.gt.1540)stop
- 50 read(2,end=99) (abuf(i),i=1,nvec)
- 5 do 1 i=1,nrad
- do 1 j=1,nloop
- 1 a(j,i)=0.0
- fl=ll
- fl1=fl+1.0
- fl3=fl*fl1
- sfl3=sqrt(fl3)
- print 900,nn,atyp,ll,ww,qq,cg,sfl3,buf(5*n)
- 900 format(i4,a4,i4,5g14.6)
+ if(iflag.eq.1) then
+ write(*,*) 'ERR003: eigcon: jcom=',jcom,
+ + ' does not fit mode ',typeo_eigen
+ stop
+ endif
+c additional normalization V, V' or W, W'
+ if(typeo_eigen.eq.'S'.or.typeo_eigen.eq.'T') then
+ fl=ll
+ fl1=fl+1.0
+ fl3=sqrt(fl*fl1)
+ do i=nstart,n
+ do j =1,2
+ if(typeo_eigen.eq.'T') then
+ buf(j,i)=buf(j,i)/fl3
+ else
+ buf(j+2,i)=buf(j+2,i)/fl3
+ endif
+ enddo
+ enddo
+ endif
qq=0.5*ww/qq
-c*** cg is really the perturbation of the grav potl at the surface!
- cg=buf(5*n)
- if(ll.eq.0) cg=0.
- k=0
- do 10 i=nstart,n
- k=k+1
- a(1,k)=buf(i)
- j=i+n
- a(2,k)=buf(j)
- j=j+n
- a(3,k)=buf(j)/sfl3
- j=j+n
- a(4,k)=buf(j)/sfl3
- 10 continue
- do 20 i=1,nlen
- 20 x(i)=aa(i)
- nz = nn
- lz = ll
- w = ww
- q = qq
- g = cg
- kfl = kfl + 1
- call gfs_rwdir(3,hed6,kfl,'w',ierr)
- goto 50
- 99 return
- end
-
- subroutine conv2(nstart,n,kfl)
- parameter(mnl=323)
- common/bufz/nn,ll,ww,qq,cg,buf(6*mnl)
- common/hedX/nhed,model,jcom,nz,atyp,lz,w,q,g,nrad,
- * rn,wn,accn,x(1540)
- dimension abuf(6*mnl+5),a(2,350),aa(1),hed6(1554)
- character*4 atyp
- character*8 model
- equivalence (nn,abuf(1)),(a(1,1),aa(1)),(hed6(1),nhed)
- atyp ='T '
- nvec = 2*n + 5
- nloop=2
- nlen = nloop*nrad
- nout = nloop*nrad+5
- print *,'nlen (must be .le. 1540 or else change defhdr) : ',nlen
- if(nlen.gt.1540)stop
- 50 read(2,end=99) (abuf(i),i=1,nvec)
- 5 do 1 i=1,nrad
- do 1 j=1,nloop
- 1 a(j,i)=0.0
- fl=ll
- fl1=fl+1.0
- fl3=fl*fl1
- sfl3=sqrt(fl3)
- print 900,nn,atyp,ll,ww,qq,cg,sfl3
- 900 format(i4,a4,i4,5g14.6)
- qq=0.5*ww/qq
- k=0
- do 10 i=nstart,n
- k=k+1
- a(1,k)=buf(i)/sfl3
- a(2,k)=buf(i+n)/sfl3
- 10 continue
- do 20 i=1,nlen
- 20 x(i)=aa(i)
- nz = nn
- lz = ll
- w = ww
- q = qq
- g = 0.
- kfl = kfl + 1
- call gfs_rwdir(3,hed6,kfl,'w',ierr)
- goto 50
- 99 return
- end
-
-
- subroutine conv1(nstart,n,kfl)
-c*** special for radial modes
- parameter(mnl=323)
- common/bufz/nn,ll,ww,qq,cg,buf(6*mnl)
- common/hedX/nhed,model,jcom,nz,atyp,lz,w,q,g,nrad,
- * rn,wn,accn,x(1540)
- dimension abuf(6*mnl+5),a(4,350),aa(1),hed6(1554)
- character*4 atyp
- character*8 model
- equivalence (nn,abuf),(a,aa),(hed6,nhed)
- atyp ='S '
- nvec = 2*n + 5
- nloop=4
- nlen = nloop*nrad
- nout = nloop*nrad+5
- print *,'nlen (must be .le. 1540 or else change defhdr) : ',nlen
- if(nlen.gt.1540)stop
- 50 read(2,end=99) (abuf(i),i=1,nvec)
- 5 do 1 i=1,nrad
- do 1 j=1,nloop
- 1 a(j,i)=0.0
- fl=ll
- fl1=fl+1.0
- fl3=fl*fl1
- sfl3=sqrt(fl3)
- print 900,nn,atyp,ll,ww,qq
- 900 format(i4,a4,i4,5g14.6)
- qq=0.5*ww/qq
- k=0
- do 10 i=nstart,n
- k=k+1
- a(1,k)=buf(i)
- a(2,k)=buf(i+n)
- a(3,k)=0.
- a(4,k)=0.
- 10 continue
- do 20 i=1,nlen
- 20 x(i)=aa(i)
- nz = nn
- lz = ll
- w = ww
- q = qq
- g = 0.
- kfl = kfl + 1
- call gfs_rwdir(3,hed6,kfl,'w',ierr)
- goto 50
- 99 return
- end
+cxx ww = ww/wn
+ call write_eigen(ieig,ierr)
+ write(idat,rec=eigid_eigen) nn,ll,ww,qq,
+ + ((buf(l,lll),l=1,ncol_eigen),lll=nstart,n)
+ goto 200
+ 99 close(7)
+ close(8)
+ call close_eigen(ieig,idat)
+ end
Added: seismo/1D/syndat/trunk/fdb/fdb_eigen.f
===================================================================
--- seismo/1D/syndat/trunk/fdb/fdb_eigen.f 2006-09-08 15:53:55 UTC (rev 4495)
+++ seismo/1D/syndat/trunk/fdb/fdb_eigen.f 2006-09-08 17:13:41 UTC (rev 4496)
@@ -0,0 +1,169 @@
+c ******************************************************************
+c Relation: eigen
+c Description: Eigenfunction and eigenvalue.
+c attribute field storage external character attribure
+c name no type format position description
+c norder 1 i4 i8 1-8 order number n
+c lorder 2 i4 i8 10-17 order number l
+c typeo 3 i4 a1 19-19 type of oscillations
+c eigid 4 i4 i8 21-28 eigen id
+c per 5 f4 f10.5 30-45 eigenvalue, period, sec
+c phvel 6 f4 f10.5 47-62 phase velocity, km/s
+c grvel 7 f4 f10.5 64-79 group velocity, km/s
+c attn 8 f4 f10.5 81-96 attenuation
+c nraw 9 i4 i8 98-105 number of raws
+c ncol 10 i4 i4 107-110 number of columns
+c npar 11 i4 i4 112-115 number of parameters
+c datatype 12 c2 a2 117-118 numeric storage
+c foff 13 i4 i10 120-129 byte offset
+c dir 14 c64 a64 131-194 directory
+c dfile 15 c32 a32 196-227 file name
+c commid 16 i4 i8 229-236 comment id
+c lddate 17 date a17 238-254 load date
+c ******************************************************************
+ subroutine open_eigen(name,n,ieig,idat,dir,mode,ierr)
+ implicit none
+ character mode *(*)
+ character*256 name,namep,named,nameb,cmd
+ character*64 dir
+ integer*4 i,n,ieig,idat,l,lnblnk,ierr
+ logical tf
+c ---
+ ierr = 0
+ l = lnblnk(name)
+ namep = name(1:l)//'.eigen'
+ named = name(1:l)//'.eigen.dat'
+ inquire(file=named,exist=tf)
+ if(.not.tf) then
+c create new dir for data
+ cmd = 'mkdir -p '//name(1:l)//'.eigen.dat'
+ call system(cmd)
+c write(*,*) 'dir created'
+ endif
+ nameb = name(1:l)//'.eigen.dat/eigen'
+ inquire(file=nameb,exist=tf)
+ if(tf.and.mode.eq.'w') then
+c delete file eigen, if exists
+ cmd = 'rm -f '//nameb
+ call system(cmd)
+c write(*,*) 'file deleted'
+ endif
+ l = lnblnk(named)
+ do i = l,1,-1
+ if(named(i:i).eq.'/'.or.named(i:i).eq.' ') then
+ dir = named(i+1:l)
+ goto 1
+ endif
+ enddo
+ dir = named(1:l)
+ 1 open(ieig,err=99,file=namep,form='formatted',status='unknown')
+ open(idat,err=99,file=nameb,access='direct',recl=n)
+ return
+ 99 ierr = 1
+ return
+ end
+c **********************************************************************
+ subroutine close_eigen(ieig,idat)
+ close(idat)
+ close(ieig)
+ return
+ end
+c **********************************************************************
+ subroutine write_eigen(ieig,ierr)
+ implicit none
+c ---
+ real*4 per_eigen,phvel_eigen,grvel_eigen,attn_eigen
+ integer*4 norder_eigen,lorder_eigen,eigid_eigen,
+ + nraw_eigen,ncol_eigen,npar_eigen,foff_eigen,
+ + commid_eigen
+ character*2 datatype_eigen
+ character*64 dir_eigen
+ character*32 dfile_eigen
+ character*17 lddate_eigen
+ character*1 typeo_eigen
+ common/c_eigen/norder_eigen,lorder_eigen,
+ + eigid_eigen,per_eigen,phvel_eigen,grvel_eigen,
+ + attn_eigen,nraw_eigen,ncol_eigen,npar_eigen,
+ + foff_eigen,commid_eigen,typeo_eigen,
+ + datatype_eigen,dir_eigen,dfile_eigen,lddate_eigen
+c ---
+ integer*4 ieig,ierr
+c ---
+ write(ieig,1000) norder_eigen,lorder_eigen,typeo_eigen,
+ + eigid_eigen,per_eigen,phvel_eigen,grvel_eigen,
+ + attn_eigen,nraw_eigen,ncol_eigen,npar_eigen,
+ + datatype_eigen,foff_eigen,dir_eigen,
+ + dfile_eigen,commid_eigen,lddate_eigen
+ return
+ 1000 format(i8,1x,i8,1x,a1,1x,i8,1x,4(f16.5,1x),i8,1x,2(i4,1x),
+ + a2,1x,i10,1x,a64,1x,a32,1x,i8,1x,a17)
+ end
+c **********************************************************************
+ subroutine read_eigen(ieig,ierr)
+ implicit none
+ common/c_eigen/norder_eigen,lorder_eigen,
+ + eigid_eigen,per_eigen,phvel_eigen,grvel_eigen,
+ + attn_eigen,nraw_eigen,ncol_eigen,npar_eigen,
+ + foff_eigen,commid_eigen,typeo_eigen,
+ + datatype_eigen,dir_eigen,dfile_eigen,lddate_eigen
+ real*4 per_eigen,phvel_eigen,grvel_eigen,attn_eigen
+ integer*4 norder_eigen,lorder_eigen,eigid_eigen,
+ + nraw_eigen,ncol_eigen,npar_eigen,foff_eigen,
+ + commid_eigen
+ character*2 datatype_eigen
+ character*64 dir_eigen
+ character*32 dfile_eigen
+ character*17 lddate_eigen
+ character*1 typeo_eigen
+ integer*4 ieig,ierr
+c ---
+ ierr = 0
+ read(ieig,1000,end=9) norder_eigen,lorder_eigen,typeo_eigen,
+ + eigid_eigen,per_eigen,phvel_eigen,grvel_eigen,
+ + attn_eigen,nraw_eigen,ncol_eigen,npar_eigen,
+ + datatype_eigen,foff_eigen,dir_eigen,
+ + dfile_eigen,commid_eigen,lddate_eigen
+ return
+ 9 ierr = 1
+ return
+ 1000 format(i8,1x,i8,1x,a1,1x,i8,1x,4(f16.5,1x),i8,1x,2(i4,1x),
+ + a2,1x,i10,1x,a64,1x,a32,1x,i8,1x,a17)
+ end
+c **********************************************************************
+ subroutine null_eigen
+ implicit none
+ common/c_eigen/norder_eigen,lorder_eigen,
+ + eigid_eigen,per_eigen,phvel_eigen,grvel_eigen,
+ + attn_eigen,nraw_eigen,ncol_eigen,npar_eigen,
+ + foff_eigen,commid_eigen,typeo_eigen,
+ + datatype_eigen,dir_eigen,dfile_eigen,lddate_eigen
+ real*4 per_eigen,phvel_eigen,grvel_eigen,attn_eigen
+ integer*4 norder_eigen,lorder_eigen,eigid_eigen,
+ + nraw_eigen,ncol_eigen,npar_eigen,foff_eigen,
+ + commid_eigen
+ character*2 datatype_eigen
+ character*64 dir_eigen
+ character*32 dfile_eigen
+ character*17 lddate_eigen
+ character*1 typeo_eigen
+c---
+c t4 - SUN IEEE real*4, t4 - VAX/PC IEEE real*4
+ norder_eigen = -1
+ lorder_eigen = -1
+ typeo_eigen = 'P'
+ eigid_eigen = 0
+ per_eigen = -1.0
+ phvel_eigen = -1.0
+ grvel_eigen = -1.0
+ attn_eigen = -1.0
+ nraw_eigen = 0
+ ncol_eigen = 0
+ npar_eigen = 0
+ datatype_eigen = 't4'
+ foff_eigen = 0
+ dir_eigen = '-'
+ dfile_eigen = '-'
+ commid_eigen = -1
+ lddate_eigen = '-'
+ return
+ end
Modified: seismo/1D/syndat/trunk/green.f
===================================================================
--- seismo/1D/syndat/trunk/green.f 2006-09-08 15:53:55 UTC (rev 4495)
+++ seismo/1D/syndat/trunk/green.f 2006-09-08 17:13:41 UTC (rev 4496)
@@ -7,22 +7,23 @@
c vecnlX : exitations of modes
c
program green
+ include "green.h"
integer n1(100),n2(100)
cxx (MB) real*4 hed1(30),alocal(9),scalrs(8),dcg(3)
real*4 hed1(30),alocal(9),scalrs(8)
real*8 co,si,c1,c2,s1,s2,z,zp,cz,sz,caz,saz,prp
common/dheadX/d0,th,ph,jy,jd,jh,jm,sec,tstart
- common/zfXX/z(3,250),zp(3,250)
- common/rcrds/dat(9000)
+ common/zfXX/z(3,ml),zp(3,ml)
+ common/rcrds/dat(mseis)
common/myhed/nscan,nsta,nchn,ntyp,iy,id,ih,im,ss,dt,slat,slon,
& sdep,jys,jds,jhs,jms,sss,spare(12)
- common/grnX/grn(108000)
- common/modeX/om(2760),a0(2760),omt(1650),a0t(1650)
- common/vecnlX/vecnl(8,2760),vecnlt(4,1650)
- common/wts/wt(6,2760)
+ common/grnX/grn(12*mseis)
+ common/modeX/om(meig),a0(meig),omt(meig),a0t(meig)
+ common/vecnlX/vecnl(8,meig),vecnlt(4,meig)
+ common/wts/wt(6,meig)
common/sclrX/n,l,e1,e2,e3,e4,au,av
- common/propX/prp(54000)
+ common/propX/prp(6*mseis)
character*4 nchn,nsta,ksta
c (MB) added for consistency
character*4 ntyp
@@ -32,10 +33,8 @@
data icomp,itype1,ksta/0,1,'xxxx'/
c
c.....open event file for which you want to compute green's functions
+c.....read first header to get source info
c
-cxx call gfs_opena(1,'event file (input) :',itype1,jret1)
-c*** read first header to get source info
-cxx call gfs_rwdir(1,hed1,1,'r',ierr)
open(12,file='green_in',status='old')
read(12,*) nscan,nsta,nchn,ntyp,iy,id,ih,im,ss,dt,slat,slon,
& sdep,jys,jds,jhs,jms,sss
@@ -267,11 +266,12 @@
subroutine prop(om,a0,dt,tst,nmodes,len,npts,nfun,prp)
+ include "green.h"
real*8 ddt,dts,dt0,dt1,prp(nfun,npts),phi,ain
real*8 b0,b1,c0,c1,c2
common/propc/ddt,dts,dt0,dt1,phi,ain,b0,b1,c0,c1,c2
- common/grnX/grn(108000)
- common/wts/wt(6,2760)
+ common/grnX/grn(12*mseis)
+ common/wts/wt(6,meig)
dimension om(*),a0(*)
c
c....stores green's function in multiplexed form
@@ -317,36 +317,84 @@
return
end
-
+c***************************************************************
+c source sub reads eigen functions from flat databases.
+c List of dbnames are defined in dbnames.dat file. source read
+c .eigen relations and select eigen functions for S & T modes
+c****************************************************************
subroutine source(fmin,fmax,lmax,nmodes,nmodet)
- character atyp*4, model*8, name1*30, name2*30
+ implicit none
+ include "green.h"
+ integer*4 mk,lmax,nmodes,nmodet
+ real*4 fmin,fmax
+c --- common blocks -------------------------------
+ real*4 d0,t0,p0,sec,tstart
+ integer*4 jy,jd,jh,jm
common/dheadX/d0,t0,p0,jy,jd,jh,jm,sec,tstart
+ real*4 x1,r0,x2,f
common/sclXXX/x1,r0,x2,f(4,3)
- common/vecnlX/vecnl(8,2760),vecnlt(4,1650)
- common/modeX/om(2760),a0(2760),omt(1650),a0t(1650)
- common/hedX/nhed,model,jcom,n,atyp,l,w,q,p,npts,
- * rn,wn,accn,buf(1540)
- dimension r(300),hed6(1554)
- equivalence (nhed,hed6)
+ real*4 vecnl,vecnlt
+ common/vecnlX/vecnl(8,meig),vecnlt(4,meig)
+ real*4 om,a0,omt,a0t
+ common/modeX/om(meig),a0(meig),omt(meig),a0t(meig)
+c --- eigen relation common block
+ real*4 per_eigen,phvel_eigen,grvel_eigen,attn_eigen
+ integer*4 norder_eigen,lorder_eigen,eigid_eigen,
+ + nraw_eigen,ncol_eigen,npar_eigen,foff_eigen,
+ + commid_eigen
+ character*2 datatype_eigen
+ character*64 dir_eigen
+ character*32 dfile_eigen
+ character*17 lddate_eigen
+ character*1 typeo_eigen
+ common/c_eigen/norder_eigen,lorder_eigen,
+ + eigid_eigen,per_eigen,phvel_eigen,grvel_eigen,
+ + attn_eigen,nraw_eigen,ncol_eigen,npar_eigen,
+ + foff_eigen,commid_eigen,typeo_eigen,
+ + datatype_eigen,dir_eigen,dfile_eigen,lddate_eigen
+c --- other variables
+ real* 4 w,q,p,rn,wn,accn
+ real*4 pi2,fot,bigg,tau,rhobar,vn,gn,rs,fl,fl1,fl3,u,v
+ real*4 wsq,e14,au,av,wr,aw
+ integer*4 jcom,npts,nrecl,ieig,idat,ierr,ifl,i,is,j,js
+ integer*4 ll,lll,m,l,n
+c ---
+ character*64 dir
+ character*256 dbname
+ real*4 r(mk),buf(6,mk)
+c ---
+ real*4 fnl(2)
+ equivalence (fnl(1),n),(fnl(2),l)
data fot/1.33333333333/
cxx (MB) data name1/'/home/guy/lfsyn.dir/sph20m.gfs'/
cxx (MB) data name2/'/home/guy/lfsyn.dir/tor20m.gfs'/
- data name1/'gfs_S'/
- data name2/'gfs_T'/
- data ityp6/6/
-c*** get spheroidal source scalars ***
- call gfs_open(3,name1,ityp6,jret3)
+ pi2 = datan(1.0d0)*8.0d0
+ nmodes = 0
+ nmodet = 0
+ lmax=0
+ ieig = 9
+ idat = 10
+c===================================================================
+c Main loop by dbase names
+c===================================================================
+ open(7,file='dbnames.dat',status='old')
+c*** read dbnames list
+ 8 read(7,'(a256)',end=9) dbname
+ nrecl = 2000
+ call open_eigen(dbname,nrecl,ieig,idat,dir,'r',ierr)
+ call read_eigen(ieig,ierr)
+ nrecl = (ncol_eigen*nraw_eigen+npar_eigen)*4
c
-c....read in radial grid
+c....read in parameters and radial grid
c
ifl=1
- call gfs_rwdir(3,hed6,ifl,'r',ierr)
-c call prhdr(ifl,hed6,ityp6)
- if(atyp.ne.'grid') stop 'grid not found'
- do 4 i=1,npts
- 4 r(i)=buf(i)
+ read(idat,rec=ifl) jcom,bigg,tau,rhobar,rn,wn,vn,gn,accn,
+ + (r(ll),ll=1,nraw_eigen)
+ call close_eigen(ieig,idat)
+ if(typeo_eigen.ne.'P') stop
+ * 'ERR010:eigen: radius pints are not found'
+ npts = nraw_eigen
rs=rn/1000.
- lmax=0
r0=1.-(d0/rs)
do 5 i=1,npts
if(r(i).lt.r0)is=i
@@ -354,22 +402,30 @@
x1=r(is)
js=is+1
x2=r(js)
- iau=4*npts-3
- iav=iau+2
c
-c....read in spheroidal modes in frequency band fmin < f < fmax
+c....read in spheroidal and toroidal modes
+c....in frequency band fmin < f < fmax
c
- nmodes = 0
- do 10 ifl=1,jret3
- call gfs_rwdir(3,hed6,ifl,'r',ierr)
+ call open_eigen(dbname,nrecl,ieig,idat,dir,'r',ierr)
+ call read_eigen(ieig,ierr)
+
+ 10 ifl = ifl+1
+ call read_eigen(ieig,ierr)
+ if(ierr.ne.0) goto 30
+ if(ifl.ne.eigid_eigen) stop
+ * 'ERR011:eigen: flat and bin indices are different.'
+ w = pi2/per_eigen
if (w.lt.fmin) goto 10
if (w.le.fmax) goto 11
-c goto 12
-c*** files may not be sorted
- goto 10
- 11 if (atyp.ne.'S') goto 10
+ goto 10
+ 11 read(idat,rec=eigid_eigen) n,l,w,q,
+ + ((buf(ll,lll),ll=1,ncol_eigen),lll=1,nraw_eigen)
+ npts = nraw_eigen
+ if (typeo_eigen.eq.'S') then
+c
+c get source scalars for S mode
+c
nmodes=nmodes+1
-c call prhdr(ifl,hed6,ityp6)
om(nmodes)=w
a0(nmodes)=q
lmax = max0(l,lmax)
@@ -377,11 +433,12 @@
fl1=fl+1.
fl3=fl*fl1
m=-1
- do 14 i=is,js
+ do i=is,js
m=m+2
- jst=(i-1)*4
- do 14 j=1,4
- 14 f(j,m)=buf(j+jst)
+ do j=1,4
+ f(j,m)=buf(j,i)
+ enddo
+ enddo
call cubic(2)
u=f(1,2)/r0
v=f(3,2)/r0
@@ -389,12 +446,13 @@
e14=f(4,2)+u-v
if(l.eq.0) e14=0.
- au=-((wsq+2.*fot)*buf(iau)+fl1*p)*accn
- av=-(wsq*buf(iav)-buf(iau)*fot-p)*accn
+ p = buf(5,npts)
+ au=-((wsq+2.*fot)*buf(1,npts)+fl1*p)*accn
+ av=-(wsq*buf(3,npts)-buf(1,npts)*fot-p)*accn
if(l.eq.0) av=0.
- vecnl(1,nmodes)=hed6(5)
- vecnl(2,nmodes)=hed6(7)
+ vecnl(1,nmodes)=fnl(1)
+ vecnl(2,nmodes)=fnl(2)
vecnl(3,nmodes)=f(2,2)
vecnl(4,nmodes)=u-.5*fl3*v
vecnl(5,nmodes)=e14
@@ -403,51 +461,39 @@
vecnl(8,nmodes)=av
c print 800,n,l,(vecnl(i,nmodes),i=3,8)
800 format(2i4,6e15.7)
- 10 continue
- 12 call gfs_close(3)
- print *,'# sph. modes in band =',nmodes,' must be .le. 2760'
-c*** get source scalars for toroidal modes ***
- ifl=1
- call gfs_open(3,name2,ityp6,jret3)
- call gfs_rwdir(3,hed6,ifl,'r',ierr)
- iaw=2*npts-1
+ else if (typeo_eigen.eq.'T') then
c
-c....read in toroidal modes in frequency band fmin < f < fmax
+c get source scalars for T mode
c
- nmodet = 0
- do 30 ifl=1,jret3
- call gfs_rwdir(3,hed6,ifl,'r',ierr)
- if (w.lt.fmin) goto 30
- if (w.le.fmax) goto 31
-c goto 32
- goto 30
- 31 if (atyp.ne.'T') goto 30
nmodet=nmodet+1
-c call prhdr(nmodet,hed6,ityp6)
omt(nmodet)=w
a0t(nmodet)=q
lmax = max0(l,lmax)
m=-1
- do 45 i=is,js
+ do i=is,js
m=m+2
- jst=(i-1)*2
- do 45 j=1,2
- 45 f(j,m)=buf(j+jst)
+ do j=1,2
+ f(j,m)=buf(j,i)
+ enddo
+ enddo
call cubic(1)
wr=f(1,2)/r0
wsq=(w/wn)**2
- aw=-buf(iaw)*wsq*accn
+ aw=-buf(1,npts)*wsq*accn
- vecnlt(1,nmodet)=hed6(5)
- vecnlt(2,nmodet)=hed6(7)
+ vecnlt(1,nmodet)=fnl(1)
+ vecnlt(2,nmodet)=fnl(2)
vecnlt(3,nmodet)=aw*(wr-f(2,2))
vecnlt(4,nmodet)=-4.*aw*wr
c print 801,n,l,(vecnlt(i,nmodet),i=3,4)
801 format(2i4,6e15.7)
- 30 continue
- 32 print *,'# tor. modes in band =',nmodet,' must be .le. 1650'
- call gfs_close(3)
+ endif
+ goto 10
+ 30 goto 8
+ 9 close(7)
+ print *,'# sph. modes in band =',nmodes,' must be .le. ',meig
+ print *,'# tor. modes in band =',nmodet,' must be .le. ',meig
return
end
@@ -622,8 +668,9 @@
c b(m,l) is given in G&D (1975) equation (21)
c and X(m,l,theta) is given in G&D (1975) equation (2)
c G&D = Gilbert & Dziewonski
+ include "green.h"
implicit real*8(a-h,o-z)
- common/zfXX/z(3,250),p(3,250)
+ common/zfXX/z(3,ml),p(3,ml)
data pi/3.14159265358979d0/
z(1,1)=1d0/(4.d0*pi)
z(2,1)=0.d0
Added: seismo/1D/syndat/trunk/green.h
===================================================================
--- seismo/1D/syndat/trunk/green.h 2006-09-08 15:53:55 UTC (rev 4495)
+++ seismo/1D/syndat/trunk/green.h 2006-09-08 17:13:41 UTC (rev 4496)
@@ -0,0 +1,18 @@
+c============================================================
+c green program. Parameters for the length of arrays:
+c mk - maximum number knots in 1-D model
+c ml - maximum value of angular mode number l
+c meig - maximum total number of eigen functions involved
+c mseis - maximun length of seismogram
+c============================================================
+ integer*4 mk
+ parameter (mk=350)
+c ---
+ integer*4 ml
+ parameter (ml=6000)
+c ---
+ integer*4 meig
+ parameter (meig=100000)
+c ---
+ integer*4 mseis
+ parameter (mseis=30000)
More information about the cig-commits
mailing list