[cig-commits] r5114 - in seismo/1D/syndat/trunk: . fdb time utils
barmin at geodynamics.org
barmin at geodynamics.org
Sun Oct 29 00:09:13 PDT 2006
Author: barmin
Date: 2006-10-29 00:09:11 -0700 (Sun, 29 Oct 2006)
New Revision: 5114
Added:
seismo/1D/syndat/trunk/fdb/fdb_eigen.h
seismo/1D/syndat/trunk/fdb/fdb_io.h
seismo/1D/syndat/trunk/fdb/fdb_site.h
seismo/1D/syndat/trunk/fdb/fdb_sitechan.h
seismo/1D/syndat/trunk/fdb/fdb_wfdisc.h
seismo/1D/syndat/trunk/fdb/fdbs.f
seismo/1D/syndat/trunk/fdb/swapn.c
seismo/1D/syndat/trunk/time/
seismo/1D/syndat/trunk/time/time.f
seismo/1D/syndat/trunk/utils/
seismo/1D/syndat/trunk/utils/endi.c
seismo/1D/syndat/trunk/utils/preigen.f
seismo/1D/syndat/trunk/utils/swap.c
seismo/1D/syndat/trunk/utils/wfendi.f
Modified:
seismo/1D/syndat/trunk/Makefile.am
seismo/1D/syndat/trunk/eigcon.f
seismo/1D/syndat/trunk/fdb/fdb_eigen.f
seismo/1D/syndat/trunk/green.f
seismo/1D/syndat/trunk/green.h
seismo/1D/syndat/trunk/minos_bran.f
seismo/1D/syndat/trunk/syndat.f
Log:
The new version without fdb suppori. This version
corresponds CIG workshop in St Louis, Oct 30, 2006.
Modified: seismo/1D/syndat/trunk/Makefile.am
===================================================================
--- seismo/1D/syndat/trunk/Makefile.am 2006-10-28 23:23:24 UTC (rev 5113)
+++ seismo/1D/syndat/trunk/Makefile.am 2006-10-29 07:09:11 UTC (rev 5114)
@@ -20,58 +20,74 @@
# $Id: Makefile.am 2452 2005-10-19 20:03:01Z leif $
bin_PROGRAMS = minos_bran syndat green eigcon \
- ghead glist gtail iutil \
- curse
+ endi wfendi preigen
+# ghead glist gtail iutil \
+# curse
# main programs
minos_bran_SOURCES = minos_bran.f
-syndat_SOURCES = syndat.f $(gfs) fftsubs.f
-green_SOURCES = green.f $(fdb) $(gfs) fftsubs.f rcode.f
-eigcon_SOURCES = eigcon.f $(fdb)
+syndat_SOURCES = syndat.f $(fdb) $(time) fftsubs.f
+green_SOURCES = green.f $(fdb) $(time) fftsubs.f
+eigcon_SOURCES = eigcon.f $(fdb) $(time)
# FDB utils
fdb = \
- fdb/fdb_eigen.f
+ fdb/fdb_eigen.f \
+ fdb/fdbs.f \
+ fdb/fdb_io.h \
+ fdb/fdb_site.h \
+ fdb/fdb_sitechan.h \
+ fdb/fdb_wfdisc.h \
+ fdb/fdb_eigen.h \
+ fdb/swapn.c
+# time utils
+time = \
+ time/time.f
+# Utilities
+endi_SOURCES = utils/endi.c utils/swap.c
+wfendi_SOURCES = utils/wfendi.f $(fdb) $(time)
+preigen_SOURCES = utils/preigen.f $(fdb) $(time)
+
# GFS utils
-ghead_SOURCES = gfs/ghead.f $(gfs)
-glist_SOURCES = gfs/glist.f $(gfs)
-gtail_SOURCES = gfs/gtail.f $(gfs)
-iutil_SOURCES = gfs/iutil.f $(gfs)
-
+# ghead_SOURCES = gfs/ghead.f $(gfs)
+# glist_SOURCES = gfs/glist.f $(gfs)
+# gtail_SOURCES = gfs/gtail.f $(gfs)
+# iutil_SOURCES = gfs/iutil.f $(gfs)
+#
# curse
-curse_SOURCES = gfs/curse.f $(gfs) $(leo)
-curse_LDADD = -L/usr/X11R6/lib -lX11
-
-gfs = \
- gfs/disk2.c \
- gfs/gfs.f \
- gfs/usr.f
-
-leo = \
- gfs/leolib.c \
- gfs/psscreendump.f \
- gfs/scr1.f
-
+# curse_SOURCES = gfs/curse.f $(gfs) $(leo)
+# curse_LDADD = -L/usr/X11R6/lib -lX11
+#
+# gfs = \
+# gfs/disk2.c \
+# gfs/gfs.f \
+# gfs/usr.f
+#
+# leo = \
+# gfs/leolib.c \
+# gfs/psscreendump.f \
+# gfs/scr1.f
+#
# included by leolib.c
-leoinc = \
- gfs/ftoc.c \
- gfs/leodraw.c \
- gfs/leoinp.c \
- gfs/leotext.c \
- $(xvertext)
-
+# leoinc = \
+# gfs/ftoc.c \
+# gfs/leodraw.c \
+# gfs/leoinp.c \
+# gfs/leotext.c \
+# $(xvertext)
+#
# xvertext 5.0 by Alan Richardson
-xvertext = \
- gfs/rotated.c \
- gfs/rotated.h
-
-EXTRA_DIST = \
- LIST.STA \
- locale.cmt.big \
- stations \
- gfs/gfs.doc \
- $(leoinc)
-
+#xvertext = \
+# gfs/rotated.c \
+# gfs/rotated.h
+#
+# EXTRA_DIST = \
+# LIST.STA \
+# locale.cmt.big \
+# stations \
+# gfs/gfs.doc \
+# $(leoinc)
+#
## end of Makefile.am
Modified: seismo/1D/syndat/trunk/eigcon.f
===================================================================
--- seismo/1D/syndat/trunk/eigcon.f 2006-10-28 23:23:24 UTC (rev 5113)
+++ seismo/1D/syndat/trunk/eigcon.f 2006-10-29 07:09:11 UTC (rev 5114)
@@ -1,48 +1,39 @@
program eigcon
c
c=====converts eigenfunction files ( output from minos )
-c into gfs-file containing only info about the
-c upper most 700 km.
-c ---this version for spheroidal and toroidal modes
+c into into .eigen relation containing only info about the
+c upper most dmax km ( 0 < dmax < Rn), where, RO is radius of
+c the free surface in km.
+c ---this version works for all modes: spheroidal, toroidal, and radial.
c
***************************************************************************
* designed: *
* changes: 15.08.91 spheroidal and toroidal modes *
* parameter mnl for max no. layers in model *
* variable irecl for reclength of eigfun file*
-* latest change: 19.08.91 subroutines for sph, rad and tor modes *
+* changes: 19.08.91 subroutines for sph, rad and tor modes *
+* latest change: 09.01.06 output is .eigen relation.
***************************************************************************
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
+ include "fdb/fdb_eigen.h"
c --- other variables
real*8 r(mk),pi
- real*4 rout(mk),buf(6,mk)
+ real*4 param(7),buf(7,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)
+ integer*4 i,ieig,idat,iin,ifanis,ifdeck,k,j,jj,jcom,l,ll,lll
+ integer*4 n,nic,noc,nreg,nn,nlay,nstart,nrad,ni,nrecl
+ integer*4 lnblnk,ierr,iflag,ititle(20)
character*20 str
character*64 dir
character*1 typeo
character*256 fmodel,fflatin,fbinin,fout
+c --
+ common/prm/ nn,ll,ww,qq,rn,vn,accn
+ equivalence (param(1),nn)
c ---
data bigg,tau,rhobar/6.6723e-11,1000.0,5515.0/
pi=4.0d0*datan(1.d0)
@@ -52,30 +43,38 @@
c read control parameters from stdin
c
c read jcom : flag to indicate mode type
+ write(*,*) '============= Program eigcon ===================='
print*,' spheroidals (3) or toroidals (2) or radial (1) or'
print*,' inner core toroidals (4) modes'
- read(*,*)jcom
+ read(*,*) jcom
+ write(*,*) jcom
c read model file name
print*,' enter name of model file'
- read(*,'(a)')fmodel
+ read(*,'(a256)')fmodel
+ write(*,*) fmodel(1:lnblnk(fmodel))
c read max depth for eigenvectors
print *,'enter max depth [km] : '
- read(5,*) dmax
+ read(*,*) dmax
+ write(*,*) dmax
c read minos_bran output text file
write(*,*) ' enter name of minos_bran output text file'
- read(*,'(a)') fflatin
+ read(*,'(a256)') fflatin
+ write(*,*) fflatin(1:lnblnk(fflatin))
c read minos_bran output binary unformatted file
write(*,*) ' minos_bran output binary unformatted file name'
- read(*,'(a)') fbinin
+ read(*,'(a256)') fbinin
+ write(*,*) fbinin(1:lnblnk(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 "path/dbase_name" or "dbase_name", where path is relative or absolute
+c path to existing directory. Under path directory or working
+c directory (no character(s) "/" in input string) it suppose to be created
+c db relation file 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
+ read(*,'(a256)')fout
+ write(*,*) fout(1:lnblnk(fout))
+ write(*,*) '===================================================='
c
c read in radial knots from model
c
@@ -111,8 +110,7 @@
c rn : radius at surface
c wn : frequency normalization
c vn : velocity normalisation
-c gn : acceleration normalisation
-c accn : ??????
+c accn : acceleration normalisation
c n : index of surface grid point
c nstart : index of lowest grid point
c nrad : # of gridpoints of reduced eigenfunctions
@@ -138,10 +136,10 @@
j=0
do i=nstart,n
j=j+1
- rout(j)=r(i)
+ buf(1,j)=r(i)
enddo
nrad=j
- write(*,*) 'n,nstart,nrad = ',n,nstart,nrad
+ write(*,*) 'eigcon: 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
@@ -152,7 +150,7 @@
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'
+ 9 stop 'ERR004:eigcon: 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')
@@ -164,26 +162,19 @@
call null_eigen
eigid_eigen = 1
foff_eigen = 0
- ncol_eigen = 2
- if(jcom.eq.3) ncol_eigen = 6
- npar_eigen = 4
+ ncol_eigen = 3
+ if(jcom.eq.3) ncol_eigen = 7
+ npar_eigen = 7
nraw_eigen = nrad
- nrecl = max0(ncol_eigen*nraw_eigen+npar_eigen,10)*4
- call open_eigen(fout,nrecl,ieig,idat,dir,'w',ierr)
+ nrecl = (ncol_eigen*nraw_eigen+npar_eigen)*4
+ call open_eigen(fout,ieig,idat,nrecl,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)
+ 200 read(8,end=99) nn,ll,ww,qq,qc,((buf(l,lll),lll=1,n),
+ + l=2,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
@@ -205,7 +196,7 @@
typeo_eigen = 'T'
else if(jcom.eq.1) then
if(typeo.ne.'s') iflag = 1
- typeo_eigen = 'R'
+ typeo_eigen = 'S'
grvel_eigen = -1.0
else
write(*,*) 'ERR002: eigcon: Unknown jcom ',jcom
@@ -216,13 +207,13 @@
+ ' does not fit mode ',typeo_eigen
stop
endif
-c additional normalization V, V' or W, W'
+c additional normalization V, V' or W, W' by sqrt(l(l+1))
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
+ do j =2,3
if(typeo_eigen.eq.'T') then
buf(j,i)=buf(j,i)/fl3
else
@@ -232,10 +223,21 @@
enddo
endif
qq=0.5*ww/qq
-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)
+c
+c Form output buffer
+c
+ do k = 2,ncol_eigen
+ j=0
+ do i=nstart,n
+ j=j+1
+ buf(k,j)=buf(k,i)
+ enddo
+ enddo
+ call put_eigen(idat,npar_eigen,param,nraw_eigen,ncol_eigen,buf,
+ + eigid_eigen,ierr)
+ eigid_eigen = eigid_eigen+1
+ foff_eigen = foff_eigen+nrecl
goto 200
99 close(7)
close(8)
Modified: seismo/1D/syndat/trunk/fdb/fdb_eigen.f
===================================================================
--- seismo/1D/syndat/trunk/fdb/fdb_eigen.f 2006-10-28 23:23:24 UTC (rev 5113)
+++ seismo/1D/syndat/trunk/fdb/fdb_eigen.f 2006-10-29 07:09:11 UTC (rev 5114)
@@ -21,45 +21,57 @@
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)
+ subroutine open_eigen(name,ieig,idat,n,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
+ integer*4 i,ieig,idat,l,lnblnk,ierr,n
logical tf
c ---
ierr = 0
l = lnblnk(name)
+c relation itself
namep = name(1:l)//'.eigen'
+ inquire(file=namep,exist=tf)
+ if((.not.tf).and.mode.eq.'r') then
+ write(*,*) 'ERR030: fdb: file ',namep(1:lnblnk(namep)),
+ * ' does not exist.'
+ ierr = 2
+ return
+ endif
+c directory for data
named = name(1:l)//'.eigen.dat'
- inquire(file=named,exist=tf)
+c data file
+ nameb = name(1:l)//'.eigen.dat/eigen'
+c relative directory - dir
+ 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 inquire(file=named,exist=tf)
if(.not.tf) then
c create new dir for data
- cmd = 'mkdir -p '//name(1:l)//'.eigen.dat'
+ cmd = 'mkdir -p '//named
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(ieig,err=98,file=namep,form='formatted',status='unknown')
open(idat,err=99,file=nameb,access='direct',recl=n)
return
+ 98 nameb = namep
99 ierr = 1
+ write(*,*) 'ERR031: fdb: Can not open file: ',
+ * nameb(1:lnblnk(nameb))
return
end
c **********************************************************************
@@ -71,22 +83,8 @@
c **********************************************************************
subroutine write_eigen(ieig,ierr)
implicit none
+ include "fdb_eigen.h"
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,
@@ -101,20 +99,7 @@
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
+ include "fdb_eigen.h"
integer*4 ieig,ierr
c ---
ierr = 0
@@ -132,21 +117,11 @@
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
+ include "fdb_eigen.h"
c---
+ character*17 loctime
+ character*2 endian
+c---
c t4 - SUN IEEE real*4, t4 - VAX/PC IEEE real*4
norder_eigen = -1
lorder_eigen = -1
@@ -159,11 +134,62 @@
nraw_eigen = 0
ncol_eigen = 0
npar_eigen = 0
- datatype_eigen = 't4'
+ datatype_eigen = endian()
foff_eigen = 0
dir_eigen = '-'
dfile_eigen = '-'
commid_eigen = -1
- lddate_eigen = '-'
+ lddate_eigen = loctime()
return
end
+c **********************************************************************
+ subroutine get_eigen(idat,npar,param,nraw,ncol,eig,ieig,ierr)
+ implicit none
+ integer*4 idat,npar,nraw,ncol,ieig,ierr,l,lll
+ real*4 param(npar),eig(7,nraw)
+
+ ierr = 0
+ read(idat,rec=ieig,err=9) (param(l),l=1,npar),
+ + ((eig(l,lll),l=1,ncol),lll=1,nraw)
+ return
+ 9 ierr = 1
+ return
+ end
+c **********************************************************************
+ subroutine put_eigen(idat,npar,param,nraw,ncol,eig,ieig,ierr)
+ implicit none
+ integer*4 idat,npar,nraw,ncol,ieig,ierr,l,lll
+ real*4 param(npar),eig(7,nraw)
+cxx real*4 paramw(npar),eigw(7,nraw)
+cxx character*2 type,endian
+
+ ierr = 0
+cxx type = endian()
+cxx if(type.ne.'t4') then
+cxx call swap(param,paramw,4,npar)
+cxx call swap(eig,eigw,4,7*nraw)
+cxx endif
+ write(idat,rec=ieig,err=9) (param(l),l=1,npar),
+ + ((eig(l,lll),l=1,ncol),lll=1,nraw)
+ return
+ 9 ierr = 1
+ return
+ end
+c **********************************************************************
+c Estimate hardware platform:
+c endian = 't4' - BIG_ENDIAN,
+c endian = 'f4' - LOW_ENDIAN,
+c **********************************************************************
+ character*2 function endian()
+ implicit none
+ integer*4 i,is
+ character*4 s
+ equivalence (is,s)
+c ---
+ s = '+ '
+ i = is/256
+ i = is-i*256
+ endian = 'f4'
+ if(i.eq.32) endian = 't4'
+ return
+ end
Added: seismo/1D/syndat/trunk/fdb/fdb_eigen.h
===================================================================
--- seismo/1D/syndat/trunk/fdb/fdb_eigen.h 2006-10-28 23:23:24 UTC (rev 5113)
+++ seismo/1D/syndat/trunk/fdb/fdb_eigen.h 2006-10-29 07:09:11 UTC (rev 5114)
@@ -0,0 +1,19 @@
+c********************************************************************
+c Header: fdb_eigen.h
+c eigen declarations (3.0)
+c (R) University of Colorado at Boulder,CIEI. Date: 08/15/2006
+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
Added: seismo/1D/syndat/trunk/fdb/fdb_io.h
===================================================================
--- seismo/1D/syndat/trunk/fdb/fdb_io.h 2006-10-28 23:23:24 UTC (rev 5113)
+++ seismo/1D/syndat/trunk/fdb/fdb_io.h 2006-10-29 07:09:11 UTC (rev 5114)
@@ -0,0 +1,9 @@
+c********************************************************************
+c Header: fdb_io.h
+c Data base io declarations (3.0)
+c (R) University of Colorado at Boulder,CIEI. Date: 02/23/2004
+c********************************************************************
+ character*256 dbin, dbout, dbcmt, dbegn, dbwfdisc
+ integer*4 idbin, idbout, idbcmt, idbegn, idbwfdisc
+ common/c_dbio/idbin,idbout,dbin,dbout,dbcmt,idbcmt,dbegn,
+ * idbegn,dbwfdisc, idbwfdisc
Added: seismo/1D/syndat/trunk/fdb/fdb_site.h
===================================================================
--- seismo/1D/syndat/trunk/fdb/fdb_site.h 2006-10-28 23:23:24 UTC (rev 5113)
+++ seismo/1D/syndat/trunk/fdb/fdb_site.h 2006-10-29 07:09:11 UTC (rev 5114)
@@ -0,0 +1,19 @@
+c********************************************************************
+c Header: fdb_site.h
+c site declarations (3.0)
+c (R) University of Colorado at Boulder,CIEI. Date: 02/24/2004
+c********************************************************************
+ integer*4 msite, nrowsite
+ parameter (msite = 10000)
+ character*4 statype_site
+ character*6 sta_site, refsta_site
+ character*17 lddate_site
+ character*50 staname_site
+ integer*4 ondate_site, offdate_site
+ real*4 lat_site,lon_site,elev_site, deast_site, dnorth_site
+ common/c_site/sta_site(msite), ondate_site(msite),
+ * offdate_site(msite),lat_site(msite), lon_site(msite),
+ * elev_site(msite), staname_site(msite),
+ * statype_site(msite), refsta_site(msite),
+ * dnorth_site(msite), deast_site(msite),
+ * lddate_site(msite), nrowsite
Added: seismo/1D/syndat/trunk/fdb/fdb_sitechan.h
===================================================================
--- seismo/1D/syndat/trunk/fdb/fdb_sitechan.h 2006-10-28 23:23:24 UTC (rev 5113)
+++ seismo/1D/syndat/trunk/fdb/fdb_sitechan.h 2006-10-29 07:09:11 UTC (rev 5114)
@@ -0,0 +1,21 @@
+c********************************************************************
+c Header: fdb_sitechan.h
+c sitechan declarations (3.0)
+c (R) University of Colorado at Boulder,CIEI. Date: 02/24/2004
+c********************************************************************
+ integer*4 msitechan, nrowsitechan
+ parameter(msitechan=10000)
+ character*4 ctype_sitechan
+ character*6 sta_sitechan
+ character*8 chan_sitechan
+ character*17 lddate_sitechan
+ character*50 descrip_sitechan
+ integer*4 ondate_sitechan, offdate_sitechan, chanid_sitechan
+ real*4 edepth_sitechan, hang_sitechan, vang_sitechan
+ common/csitechan/sta_sitechan(msitechan),
+ * chan_sitechan(msitechan),
+ * ondate_sitechan(msitechan), chanid_sitechan(msitechan),
+ * offdate_sitechan(msitechan), ctype_sitechan(msitechan),
+ * edepth_sitechan(msitechan), hang_sitechan(msitechan),
+ * vang_sitechan(msitechan), descrip_sitechan(msitechan),
+ * lddate_sitechan(msitechan), nrowsitechan
Added: seismo/1D/syndat/trunk/fdb/fdb_wfdisc.h
===================================================================
--- seismo/1D/syndat/trunk/fdb/fdb_wfdisc.h 2006-10-28 23:23:24 UTC (rev 5113)
+++ seismo/1D/syndat/trunk/fdb/fdb_wfdisc.h 2006-10-29 07:09:11 UTC (rev 5114)
@@ -0,0 +1,28 @@
+c********************************************************************
+c Header: fdb_wfdisc.h
+c wfdisc declarations (3.0)
+c (R) University of Colorado at Boulder,CIEI. Date: 02/24/2004
+c********************************************************************
+ integer*4 mwfdisc, nrowwfdisc
+ parameter(mwfdisc = 10000)
+ character segtype_wfdisc, clip_wfdisc
+ character*2 datatype_wfdisc
+ character*6 sta_wfdisc, instype_wfdisc
+ character*8 chan_wfdisc
+ character*17 lddate_wfdisc
+ character*32 dfile_wfdisc
+ character*64 dir_wfdisc
+ integer*4 wfid_wfdisc, chanid_wfdisc, jdate_wfdisc,
+ * nsamp_wfdisc, foff_wfdisc, commid_wfdisc
+ real*4 samprate_wfdisc, calib_wfdisc, calper_wfdisc
+ real*8 time_wfdisc, endtime_wfdisc
+ common/c_wfdisc/sta_wfdisc(mwfdisc), chan_wfdisc(mwfdisc),
+ * time_wfdisc(mwfdisc), wfid_wfdisc(mwfdisc),
+ * chanid_wfdisc(mwfdisc), jdate_wfdisc(mwfdisc),
+ * endtime_wfdisc(mwfdisc), nsamp_wfdisc(mwfdisc),
+ * samprate_wfdisc(mwfdisc), calib_wfdisc(mwfdisc),
+ * calper_wfdisc(mwfdisc), instype_wfdisc(mwfdisc),
+ * segtype_wfdisc(mwfdisc), datatype_wfdisc(mwfdisc),
+ * clip_wfdisc(mwfdisc), dir_wfdisc(mwfdisc),
+ * dfile_wfdisc(mwfdisc), foff_wfdisc(mwfdisc),
+ * commid_wfdisc(mwfdisc), lddate_wfdisc(mwfdisc), nrowwfdisc
Added: seismo/1D/syndat/trunk/fdb/fdbs.f
===================================================================
--- seismo/1D/syndat/trunk/fdb/fdbs.f 2006-10-28 23:23:24 UTC (rev 5113)
+++ seismo/1D/syndat/trunk/fdb/fdbs.f 2006-10-29 07:09:11 UTC (rev 5114)
@@ -0,0 +1,613 @@
+c*********************************************************************
+c Read an site relation
+c idbin -- database handle -- in fdb_io.h
+c*********************************************************************
+ subroutine read_site
+ implicit none
+ include 'fdb_site.h'
+ include 'fdb_io.h'
+c === Local variables ===
+ integer*4 i, j, n, ind, idbsite, lnblnk
+ character*200 str,str1
+c declarations specific to this routine
+
+cxx print *,'reading site relation.....'
+
+ idbsite=58
+ nrowsite = 0
+ open(idbsite,file=dbin(1:lnblnk(dbin))//'.site',
+ * err = 99, status='old')
+
+ i=1
+ 1 if(i.gt.msite) then
+ print *, 'ERR032: fdb: Dimension of site arrays too small'
+ print *,i,' > ',msite
+ stop
+ endif
+ read(idbsite,2,end=999)
+ * sta_site(i), ondate_site(i), offdate_site(i),
+ * lat_site(i), lon_site(i), elev_site(i), staname_site(i),
+ * statype_site(i), refsta_site(i), dnorth_site(i),
+ * deast_site(i), lddate_site(i)
+c
+c fix IRIS bug in rdseed
+c
+ if(lnblnk(statype_site(i)).eq.0) then
+ read(idbsite,'(a200)',end=999) str
+ n = lnblnk(str)
+ ind = 1
+ do j = 1,n
+ if(str(j:j).ne.' ') then
+ ind = j
+ goto 10
+ endif
+ enddo
+ 10 str1 = str(ind:n)
+ read(str1,3)
+ * statype_site(i), refsta_site(i), dnorth_site(i),
+ * deast_site(i), lddate_site(i)
+ endif
+c
+c fix end
+c
+ 2 format(a6,1x,2(i8,1x),3(f9.4,1x),a50,1x,a4,1x,a6,1x,
+ * 2(f9.4,1x),a17)
+ 3 format(a4,1x,a6,1x,2(f9.4,1x),a17)
+ i=i+1
+ goto 1
+ 999 nrowsite = i - 1
+ close(idbsite)
+ 99 if(nrowsite.le.0) then
+ write(*,*)
+ * 'ERR033: fdb: site relation ',dbin(1:lnblnk(dbin)),
+ * ' appears to be empty, or does not exist'
+ endif
+ return
+ end
+
+c*********************************************************************
+c Read an sitechan relation
+c idbin -- database handle -- in fdb_io.h
+c*********************************************************************
+ subroutine read_sitechan
+ implicit none
+ include 'fdb_sitechan.h'
+ include 'fdb_io.h'
+c === Local variables ===
+ integer*4 i, idbsitechan, lnblnk
+c declarations specific to this routine
+
+cxx print *,'reading sitechan relation.....'
+
+ idbsitechan=58
+ nrowsitechan = 0
+ open(idbsitechan,file=dbin(1:lnblnk(dbin))//'.sitechan',
+ * err = 99, status='old')
+
+ i=1
+ 1 if(i.gt.msitechan) then
+ print *, 'ERR032: Dimension of sitechan arrays too small'
+ print *,i,' > ',msitechan
+ stop
+ endif
+ read(idbsitechan,2,end=999) sta_sitechan(i),chan_sitechan(i),
+ * ondate_sitechan(i),chanid_sitechan(i),offdate_sitechan(i),
+ * ctype_sitechan(i),edepth_sitechan(i),hang_sitechan(i),
+ * vang_sitechan(i),descrip_sitechan(i),lddate_sitechan(i)
+ 2 format(a6,1x,a8,1x,3(i8,1x),a4,1x,f9.4,1x,2(f6.1,1x),a50,1x,a17)
+ i=i+1
+ goto 1
+ 999 nrowsitechan = i - 1
+ close(idbsitechan)
+ 99 if(nrowsitechan.le.0) then
+ write(*,*)
+ * 'ERR033: fdb: sitechan relation ',dbin(1:lnblnk(dbin)),
+ * ' appears to be empty, or does not exist'
+ endif
+ return
+ end
+
+c*********************************************************************
+c Write a site relation
+c*********************************************************************
+ subroutine write_site
+ implicit none
+ include 'fdb_site.h'
+ include 'fdb_io.h'
+c === Local variables ===
+cxx character*24 date
+ integer*4 i, idbsite, lnblnk
+
+cxx print *,'writing site relation.....'
+
+ idbsite=58
+ open(idbsite,file=dbout(1:lnblnk(dbout))//'.site',
+ * access='APPEND')
+cxx call fdate(date)
+cxx lddate_site(1)=date(5:10)//','//date(21:24)
+
+ do i=1,nrowsite
+ write(idbsite,2) sta_site(i), ondate_site(i), offdate_site(i),
+ * lat_site(i), lon_site(i), elev_site(i), staname_site(i),
+ * statype_site(i), refsta_site(i), dnorth_site(i),
+ * deast_site(i), lddate_site(i)
+ 2 format(a6,1x,2(i8,1x),3(f9.4,1x),a50,1x,a4,1x,a6,1x,
+ * 2(f9.4,1x),a17)
+ enddo
+
+ close(idbsite)
+ return
+ end
+
+c*********************************************************************
+c Write a sitechan relation
+c writes sitechan table located in common\csitechan
+c to the specified data base named dbout.
+c idbout -- input in fdb_io.h.
+c WARNING: if JSPC_DBPATH is set, it will override
+c the specified output database dbout.
+c*********************************************************************
+ subroutine write_sitechan
+ implicit none
+ include 'fdb_sitechan.h'
+ include 'fdb_io.h'
+c === Local variables ===
+cxx character*24 date
+ integer*4 i, idbsitechan, lnblnk
+
+cxx print *,'writing sitechan relation.....'
+
+ idbsitechan=58
+ open(idbsitechan,file=dbout(1:lnblnk(dbout))//'.sitechan',
+ * access='APPEND')
+cxx call fdate(date)
+cxx lddate_sitechan(1)=date(5:10)//','//date(21:24)
+
+ do i=1,nrowsitechan
+ write(idbsitechan,2) sta_sitechan(i),chan_sitechan(i),
+ * ondate_sitechan(i),chanid_sitechan(i),offdate_sitechan(i),
+ * ctype_sitechan(i),edepth_sitechan(i),hang_sitechan(i),
+ * vang_sitechan(i),descrip_sitechan(i),lddate_sitechan(i)
+ 2 format(a6,1x,a8,1x,3(i8,1x),a4,1x,f9.4,1x,2(f6.1,1x),a50,
+ * 1x,a17)
+ enddo
+
+ close(idbsitechan)
+ return
+ end
+
+c*********************************************************************
+c sets default values for variables in common c_site, row = ituple
+c*********************************************************************
+ subroutine default_site(ituple)
+ implicit none
+ include 'fdb_site.h'
+ integer*4 ituple
+ character*17 loctime
+c a valid entry for sta is required under the CSS v. 3.0 standard
+ sta_site(ituple) = '-'
+c a valid entry for ondate is required under the CSS v. 3.0 standard
+ ondate_site(ituple) = -1
+ offdate_site(ituple) = -1
+ lat_site(ituple) = -999.0
+ lon_site(ituple) = -999.0
+ elev_site(ituple) = -999.0
+ staname_site(ituple) = '-'
+ statype_site(ituple) = '-'
+ refsta_site(ituple) = '-'
+ dnorth_site(ituple) = 0.0
+ deast_site(ituple) = 0.0
+ lddate_site(ituple) = loctime()
+
+ return
+ end
+
+c*********************************************************************
+c sets default values for variables in common c_sitechan, row = ituple
+c*********************************************************************
+ subroutine default_sitechan(ituple)
+ implicit none
+ include 'fdb_sitechan.h'
+ integer*4 ituple
+ character*17 loctime
+c a valid entry for sta is required under the CSS v. 3.0 standard
+ sta_sitechan(ituple) = '-'
+c a valid entry for chan is required under the CSS v. 3.0 standard
+ chan_sitechan(ituple) = '-'
+c a valid entry for ondate is required under the CSS v. 3.0 standard
+ ondate_sitechan(ituple) = -1
+ chanid_sitechan(ituple) = -1
+ offdate_sitechan(ituple) = -1
+ ctype_sitechan(ituple) = '-'
+ edepth_sitechan(ituple) = -9.9999
+c a valid entry for hang is required under the CSS v. 3.0 standard
+ hang_sitechan(ituple) = -999.0
+c a valid entry for vang is required under the CSS v. 3.0 standard
+ vang_sitechan(ituple) = -999.0
+ descrip_sitechan(ituple) = '-'
+ lddate_sitechan(ituple) = loctime()
+
+ return
+ end
+c*********************************************************************
+c Read an wfdisc relation
+c*********************************************************************
+ subroutine read_wfdisc
+ implicit none
+ include 'fdb_wfdisc.h'
+ include 'fdb_io.h'
+c === Local variables ===
+ integer*4 i, lnblnk
+c declarations specific to this routine
+
+cxx print *,'reading wfdisc relation.....'
+
+ idbwfdisc=58
+ nrowwfdisc = 0
+ open(idbwfdisc,file=dbin(1:lnblnk(dbin))//'.wfdisc',
+ * err = 99, status='old')
+
+ i=1
+ 1 if(i.gt.mwfdisc) then
+ print *, 'ERR032: fdb:Dimension of wfdisc arrays too small'
+ print *,i,' > ',mwfdisc
+ stop
+ endif
+ read(idbwfdisc,2,end=999) sta_wfdisc(i), chan_wfdisc(i),
+ * time_wfdisc(i), wfid_wfdisc(i), chanid_wfdisc(i),
+ * jdate_wfdisc(i),
+ * endtime_wfdisc(i), nsamp_wfdisc(i), samprate_wfdisc(i),
+ * calib_wfdisc(i), calper_wfdisc(i), instype_wfdisc(i),
+ * segtype_wfdisc(i), datatype_wfdisc(i), clip_wfdisc(i),
+ * dir_wfdisc(i), dfile_wfdisc(i), foff_wfdisc(i),
+ * commid_wfdisc(i), lddate_wfdisc(i)
+ 2 format(a6,1x,a8,1x,f17.5,1x,3(i8,1x),f17.5,1x,i8,1x,f11.7,1x,
+ * 2(f16.6,1x),a6,1x,a1,1x,a2,1x,a1,1x,a64,1x,a32,1x,i10,
+ * 1x,i8,1x,a17)
+ i=i+1
+ goto 1
+ 999 nrowwfdisc = i - 1
+ close(idbwfdisc)
+ 99 if(nrowwfdisc.le.0) then
+ write(*,*)
+ * 'ERR033: fdb: wfdisc relation ',dbin(1:lnblnk(dbin)),
+ * ' appears to be empty, or does not exist'
+ endif
+ return
+ end
+
+c*********************************************************************
+c Write a wfdisc relation
+c*********************************************************************
+ subroutine write_wfdisc
+ implicit none
+ include 'fdb_wfdisc.h'
+ include 'fdb_io.h'
+c === Local variables ===
+cxx character*24 date
+ integer*4 i, lnblnk
+
+cxx print *,'writing wfdisc relation.....'
+
+ idbwfdisc=58
+ open(idbwfdisc,file=dbout(1:lnblnk(dbout))//'.wfdisc',
+ * access='APPEND')
+cxx call fdate(date)
+cxx ddate_wfdisc(1)=date(5:10)//','//date(21:24)
+
+ do i=1,nrowwfdisc
+ write(idbwfdisc,2) sta_wfdisc(i), chan_wfdisc(i),
+ * time_wfdisc(i), wfid_wfdisc(i), chanid_wfdisc(i),
+ * jdate_wfdisc(i),
+ * endtime_wfdisc(i), nsamp_wfdisc(i), samprate_wfdisc(i),
+ * calib_wfdisc(i), calper_wfdisc(i), instype_wfdisc(i),
+ * segtype_wfdisc(i), datatype_wfdisc(i), clip_wfdisc(i),
+ * dir_wfdisc(i), dfile_wfdisc(i), foff_wfdisc(i),
+ * commid_wfdisc(i), lddate_wfdisc(i)
+ 2 format(a6,1x,a8,1x,f17.5,1x,3(i8,1x),f17.5,1x,i8,1x,f11.7,1x,
+ * 2(f16.6,1x),a6,1x,a1,1x,a2,1x,a1,1x,a64,1x,a32,1x,i10,
+ * 1x,i8,1x,a17)
+ enddo
+
+ close(idbwfdisc)
+ return
+ end
+
+c***************************************************************
+c select valid sitechan relatons, sort, and split into chn-sta groups
+c***************************************************************
+ subroutine select_sitechan(jdate,nchn,numchan,numsta)
+ implicit none
+ include 'fdb_site.h'
+ include 'fdb_sitechan.h'
+ include 'fdb_io.h'
+ integer*4 jdate,nchn
+ integer*4 numchan(msitechan),numsta(msite)
+c --- local variables
+ integer*4 i,i1,i2,iflag,ioff,ii,n,j
+ character*6 stanew,staold
+ nchn = 0
+c create chanid for .sitechan
+ do i = 1,nrowsitechan
+ chanid_sitechan(i)=i
+ enddo
+c select subset of chanid by ondate
+ nchn = 0
+ do i = 1,nrowsitechan
+ if(ondate_sitechan(i).ne.-1) then
+ ioff = offdate_sitechan(i)
+ if(ioff.eq.-1) ioff = 3001001
+ if(ondate_sitechan(i).le.jdate.and.jdate.le.ioff) then
+ nchn = nchn+1
+ numchan(nchn) = i
+ else
+ chanid_sitechan(i) = -1
+ endif
+ endif
+ enddo
+ if(nchn.le.1) goto 9
+c sort sitechan relation by sta and channel
+ do i = nchn-1,1,-1
+ do j = 1,i
+ i1 = numchan(j)
+ i2 = numchan(j+1)
+ iflag = 0
+ if(sta_sitechan(i2).lt.sta_sitechan(i1)) iflag = 1
+ if(sta_sitechan(i2).eq.sta_sitechan(i1).and.
+ * chan_sitechan(i2).gt.chan_sitechan(i1)) iflag = 1
+ if(iflag.eq.1) then
+ ii = numchan(j)
+ numchan(j) = numchan(j+1)
+ numchan(j+1) = ii
+ endif
+ enddo
+ enddo
+c remove duplicates from sitechan relation
+ n = 1
+ i1 = iabs(numchan(1))
+ do i = 2,nchn
+ i2 = iabs(numchan(i))
+ if(sta_sitechan(i1).ne.sta_sitechan(i2).or.
+ * chan_sitechan(i1).ne.chan_sitechan(i2)) then
+ if(sta_sitechan(i1).ne.sta_sitechan(i2).or.
+ * chan_sitechan(i1)(1:2).ne.chan_sitechan(i2)(1:2)) then
+ numchan(n) = -numchan(n)
+ endif
+ n = n+1
+ numchan(n) = i2
+ i1 = iabs(numchan(n))
+ endif
+ enddo
+ numchan(n) = -numchan(n)
+ nchn = n
+ 9 continue
+ if(nchn.le.1) goto 99
+c remove sitechan for nonexistent stations
+ staold = ' '
+ n = 1
+ do i = 1,nchn
+ i1 = iabs(numchan(i))
+ stanew = sta_sitechan(i1)
+ if(stanew.ne.staold) then
+ iflag = 0
+ do j = 1,nrowsite
+ if(stanew.eq.sta_site(j)) then
+ if(ondate_site(j).ne.-1) then
+ ioff = offdate_site(j)
+ if(ioff.eq.-1) ioff = 3001001
+ if(ondate_site(j).le.jdate.and.jdate.le.ioff) then
+ iflag = j
+ goto 8
+ endif
+ endif
+ endif
+ enddo
+ 8 continue
+ endif
+ if(iflag.ne.0) then
+ numchan(n) = numchan(i)
+ numsta(n) = iflag
+ n = n+1
+ endif
+ staold = stanew
+ enddo
+ nchn = n-1
+ 99 continue
+ return
+ end
+
+c***************************************************************
+c output wfdisc data
+c***************************************************************
+ subroutine put_wfdisc(ituple,n,data,ierr)
+ implicit none
+ include 'fdb_wfdisc.h'
+ include 'fdb_io.h'
+ integer*4 ituple,n,ierr
+ real*4 data(n)
+cxx real*4 data1(n)
+ integer*4 i,i1,iwf,l,n1,nd,lnblnk,factor2
+ character*256 named,cmd
+cxx character*2 endian
+ logical tf
+
+ ierr = 0
+ idbwfdisc = 59
+c check that data directory exists
+ l = lnblnk(dbout)
+ named = dbout(1:l)//'.wfdisc.dat'
+ inquire(file=named,exist=tf)
+ if(.not.tf) then
+c create new dir for data
+ cmd = 'mkdir -p '//named
+ call system(cmd)
+ endif
+c select relative subdurectory
+ l = lnblnk(named)
+ do i = l,1,-1
+ if(named(i:i).eq.'/'.or.named(i:i).eq.' ') then
+ dir_wfdisc(ituple) = named(i+1:l)
+ goto 1
+ endif
+ enddo
+ dir_wfdisc(ituple) = named(1:l)
+c write data
+ 1 dbwfdisc = named(1:l)//'/'//dfile_wfdisc(ituple)
+ n1 = factor2(n)
+ nd = n1*4
+cxx if(datatype_wfdisc(ituple).ne.endian()) then
+cxx call swap(data,data1,4,n)
+cxx endif
+ open(idbwfdisc,err=9,file=dbwfdisc,access='direct',recl=nd)
+ iwf = 0
+ do i = 1,n,n1
+ iwf = iwf+1
+ i1 = i+n1-1
+ write(idbwfdisc,rec=iwf,err=9) (data(l),l=i,i1)
+ enddo
+ close(idbwfdisc)
+ return
+ 9 ierr =1
+ return
+ end
+
+c***************************************************************
+c read wfdisc data
+c***************************************************************
+ subroutine get_wfdisc(ituple,n,data,ierr)
+ implicit none
+ include 'fdb_wfdisc.h'
+ include 'fdb_io.h'
+ integer*4 ituple,n,ierr
+ real*4 data(n)
+ integer*4 i,i1,iwf,l,ll,n1,nd,lnblnk,factor2
+ character*256 dir,named
+ character*2 endian
+ logical tf
+c--
+ ierr = 0
+ idbwfdisc = 59
+c create name to file
+ l = lnblnk(dbin)
+ named = dbin(1:l)//'.wfdisc.dat'
+c path to directory with .wfdisc relation
+ l = lnblnk(named)
+ do i = l,1,-1
+ if(named(i:i).eq.'/'.or.named(i:i).eq.' ') then
+ dir = named(1:i)
+ goto 1
+ endif
+ enddo
+ dir = ''
+ 1 ll = lnblnk(dir_wfdisc(ituple))
+ if(lnblnk(dir).eq.0) then
+ dbwfdisc = dir_wfdisc(ituple)(1:ll)//'/'//dfile_wfdisc(ituple)
+ else
+ dbwfdisc = dir(1:lnblnk(dir))//dir_wfdisc(ituple)(1:ll)//'/'//
+ * dfile_wfdisc(ituple)
+ endif
+ inquire(file=dbwfdisc,exist=tf)
+ if(.not.tf) then
+ l = lnblnk(dbwfdisc)
+ write(*,*) 'ERR030: fdb: file ',dbwfdisc(1:l),' does not exist'
+ goto 9
+ endif
+c read data
+ n1 = factor2(n)
+ nd = n1*4
+ open(idbwfdisc,err=9,file=dbwfdisc,access='direct',recl=nd)
+ iwf = 0
+ do i = 1,n,n1
+ iwf = iwf+1
+ i1 = i+n1-1
+ read(idbwfdisc,rec=iwf,err=9) (data(l),l=i,i1)
+ enddo
+ close(idbwfdisc)
+ if(datatype_wfdisc(ituple).ne.endian()) then
+ call swap1(data,4,n)
+ endif
+ return
+ 9 ierr =1
+ return
+ end
+
+c*********************************************************************
+c sets default values for variables in common c_wfdisc, row = ituple
+c*********************************************************************
+ subroutine default_wfdisc(ituple)
+ implicit none
+ include 'fdb_wfdisc.h'
+ integer*4 ituple
+ character*2 endian
+ character*17 loctime
+c a valid entry for sta is required under the CSS v. 3.0 standard
+ sta_wfdisc(ituple) = '-'
+c a valid entry for chan is required under the CSS v. 3.0 standard
+ chan_wfdisc(ituple) = '-'
+c a valid entry for wfid is required under the CSS v. 3.0 standard
+ wfid_wfdisc(ituple) = 0
+ time_wfdisc(ituple) = -9999999999.0
+ endtime_wfdisc(ituple) = 9999999999.0
+c a valid entry for nsamp is required under the CSS v. 3.0 standard
+ nsamp_wfdisc(ituple) = -1
+c a valid entry for samprate is required under the CSS v. 3.0 standard
+ samprate_wfdisc(ituple) = 0.0
+c a valid entry for calib is required under the CSS v. 3.0 standard
+ calib_wfdisc(ituple) = 0.0
+c a valid entry for calper is required under the CSS v. 3.0 standard
+ calper_wfdisc(ituple) = 0.0
+ instype_wfdisc(ituple) = '-'
+ segtype_wfdisc(ituple) = '-'
+ datatype_wfdisc(ituple) = endian()
+ clip_wfdisc(ituple) = '-'
+c a valid entry for dir is required under the CSS v. 3.0 standard
+ dir_wfdisc(ituple) = '-'
+c a valid entry for dfile is required under the CSS v. 3.0 standard
+ dfile_wfdisc(ituple) = '-'
+c a valid entry for foff is required under the CSS v. 3.0 standard
+ foff_wfdisc(ituple) = 0
+ commid_wfdisc(ituple) = -1
+ lddate_wfdisc(ituple) = loctime()
+
+ return
+ end
+
+c****************************************************************
+c get number of factors 2 3 5 7 close to 2000
+c****************************************************************
+ integer*4 function factor2(n)
+ implicit none
+ integer*4 n
+ integer*4 i,j,n7,nb,ne,icon(4)
+ data icon/2,5,3,7/
+
+ n7 = 0
+ nb = n
+ ne = 1
+c investigate factor 7
+ do i = 1,4
+ if(mod(nb,7).eq.0) then
+ n7 = n7+1
+ nb = nb/7
+ ne = ne*7
+ endif
+ enddo
+ if(n7.ge.4) then
+ factor2 = 7*7*7*7
+ return
+ endif
+c investigate other factors
+ do i =1,11
+ do j = 1,4
+ if(mod(nb,icon(j)).eq.0) then
+ nb = nb/icon(j)
+ ne = ne*icon(j)
+ goto 2
+ endif
+ enddo
+ 2 continue
+ if(ne.ge.1000) goto 1
+ enddo
+ 1 factor2 = ne
+ return
+ end
Added: seismo/1D/syndat/trunk/fdb/swapn.c
===================================================================
--- seismo/1D/syndat/trunk/fdb/swapn.c 2006-10-28 23:23:24 UTC (rev 5113)
+++ seismo/1D/syndat/trunk/fdb/swapn.c 2006-10-29 07:09:11 UTC (rev 5114)
@@ -0,0 +1,35 @@
+/*
+ * Byte swap an array of length n of N byte integer*4 or real*4 elements.
+ * A good compiler should unroll the inner loops. Letting the compiler do it
+ * gives us portability. Note that we might want to isolate the
+ * cases N = 2, 4, 8 (and 16 for long double and perhaps long long)
+ */
+#include <string.h>
+#include "config.h"
+#define SWAP2(a, b) { (a) ^= (b); (b) ^= (a); (a) ^= (b); }
+
+F77_FUNC(swap,SWAP)(c, b, N, n)
+ char *c;
+ char *b;
+ int *N;
+ int *n;
+{
+ int i, j;
+ memcpy((void *)b, (void *)c,(*n)*(*N));
+ for (i = 0; i < (*n)*(*N); i += *N)
+ for (j = 0; j < *N/2; j ++)
+ SWAP2(b[i + j], b[i + *N - j - 1]);
+ return;
+}
+
+F77_FUNC(swap1,SWAP1)(b, N, n)
+ unsigned char *b;
+ int *N;
+ int *n;
+{
+ int i, j;
+ for (i = 0; i < (*n)*(*N); i += *N)
+ for (j = 0; j < *N/2; j ++)
+ SWAP2(b[i + j], b[i + *N - j - 1]);
+ return;
+}
Modified: seismo/1D/syndat/trunk/green.f
===================================================================
--- seismo/1D/syndat/trunk/green.f 2006-10-28 23:23:24 UTC (rev 5113)
+++ seismo/1D/syndat/trunk/green.f 2006-10-29 07:09:11 UTC (rev 5114)
@@ -1,47 +1,93 @@
-
c.....program green
c common blocks :
c dheadX : source info
-c myhedX : record header
c modeX : complex eigenfrequencies
c vecnlX : exitations of modes
c
program green
+ implicit none
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)
+ include "fdb/fdb_io.h"
+ include "fdb/fdb_site.h"
+ include "fdb/fdb_sitechan.h"
+ include "fdb/fdb_wfdisc.h"
+c---
+ real*4 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,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(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(6*mseis)
- character*4 nchn,nsta,ksta
-c (MB) added for consistency
- character*4 ntyp
-
- equivalence (n,scalrs),(hed1,nscan),(alocal,d0)
- data pi,rad,zero,tstart/3.14159265359,57.29578,0.,0./
- data icomp,itype1,ksta/0,1,'xxxx'/
+c--- local variables ---
+ integer*4 lnblnk,numchan(msitechan),numsta(msite)
+ real*8 time,endtime,htoepoch
+ real*4 d0,th,ph,sec,tstart,ss,dt,slat,slon,sdep,
+ * sss, grn, om, a0, omt, a0t, vecnl, vecnlt,
+ * wt,e1,e2,e3,e4,au,av,pi,rad,fmin,fmax,
+ * t0,p0,c0,s0,t1,p1,dp,del,azim
+ integer*4 i,lmax,nmodes,nmodet,j,lp1,indx,
+ * jy,jd,jh,jm,nscan,iy,id,ih,im,jys,jds,jhs,jms,
+ * n,l,len,mchn,icomp,itype1,iscan,ifl,isq
+ integer*4 nchan,inchn,iseq,jdate,ierr,nc(100),ns(100)
+ real*4 ang(3)
+ character*8 ename
+ character*256 efname
+c--- equivalence and data ---
+ equivalence (n,scalrs)
+ data pi,rad,tstart/3.14159265359,57.29578,0./
+ data icomp,itype1/0,1/
c
+c read input parameters
+c
+c....read in dbname with static relations: sta, stachan
+ write(*,*) '============= Program green ===================='
+ write(*,*) 'enter path to db with sta & stachan:'
+ read(*,'(a256)') dbin
+ write(*,*) dbin(1:lnblnk(dbin))
+c....read in name of file containing list of dbnames.
+c....each dbname in list refferes ito db with .eigen relation.
+ write(*,*) 'enter name of file within list of nmodes db:'
+ read(*,'(a256)') fname3
+ write(*,*) fname3(1:lnblnk(fname3))
+c....read in file within event and moment tensor info
+ write(*,*) 'enter input CMT file name:'
+ read(*,'(a256)') fname4
+ write(*,*) fname4(1:lnblnk(fname4))
+c....read in frequency range in mHz
+ write(*,*) 'min and max frequencies to be considered (mHz) : '
+ read(*,*) fmin,fmax
+ write(*,*) fmin,fmax
+c....read in number of samples in Green function
+ write(*,*) 'enter # pts in greens fns .le. ',mseis,' :'
+ read(*,*) iscan
+ write(*,*) iscan
+c.... read in output path to gsf name to store Green functions
+ write(*,*) 'enter Green functions output db file name:'
+ read(*,'(a256)') dbout
+ write(*,*) dbout(1:lnblnk(dbout))
+ write(*,*) '===================================================='
+c
c.....open event file for which you want to compute green's functions
-c.....read first header to get source info
+c.....read first line to get source and moment tensor info
c
- 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
- do i=1,12
- spare(i) =0.0
- enddo
+ open(12,file=fname4,status='old')
+ read(12,*) ename,jys,jds,jhs,jms,sss,slat,slon,sdep,dt
close (12)
+ write(*,1001) ename,jys,jds,jhs,jms,sss,slat,slon
+ write(*,1002) sdep
+ 1001 format(' green: Event: ',a8,2x,i4,1x,i3,1x,i2,':',i2,':',f6.3,1x,
+ * 'lat = ',f8.3,', lon = ',f9.3)
+ 1002 format(' green: source depth =',f5.1,' km')
+ write(*,1003) dt,iscan
+ 1003 format(' green: step = ',f8.3,' sec, nsamples =',i7)
+c convert human time to epoch time
+ time = htoepoch (jys,jds,jhs,jms,dble(sss))
+ endtime = time+dt*(iscan-1)
+ jdate = jys*1000+jds
c--
d0=sdep
th=90.-slat
@@ -54,26 +100,23 @@
c
c.....open output file for green's functions
c
- call gfs_opena(2,'greens function file (output) :',itype1,jret2)
- call gfs_erase(2)
- print*,'min and max frequencies to be considered (mHz) : '
- read*,fmin,fmax
fmin = fmin*pi/500.
fmax = fmax*pi/500.
call source(fmin,fmax,lmax,nmodes,nmodet)
- if(lmax.lt.meig) goto 5
- print*,'max l =',lmax,' must be .le.',meig
- goto 99
- 5 print 105,d0,th,ph
- 105 format('source depth =',f5.0,' km, colat. and long. =',2f7.2)
- print*,'# pts in greens fns ( .le. 9000) :'
- read*,iscan
- iscan = iscan/2
+ if(lmax.le.ml) goto 5
+ write(*,*) 'ERR010: green: max l =',lmax,' must be .le.',ml
+ stop ' '
+cxx goto 99
+ 5 iscan = iscan/2
call factor(iscan)
iscan = 2*iscan
- iscan=min0(iscan,9000)
+ if(iscan.gt.mseis) then
+ write(*,*) 'WARNING: green: # of points in Green ',
+ * 'functions is stripped to ',mseis
+ endif
+ iscan=min0(iscan,mseis)
t0=th/rad
p0=ph/rad
c0=cos(t0)
@@ -82,44 +125,82 @@
c====loop over records
c........read event header
c
- open(12,file='green_in',status='old')
-cxx do 90 ifl=1,jret1
- do 90 ifl=1,3
-cxx call gfs_rwdir(1,hed1,ifl,'r',ierr)
- read(12,*) nscan,nsta,nchn,ntyp,iy,id,ih,im,ss,dt,slat,slon,
- & sdep,jys,jds,jhs,jms,sss
- ierr=0
- write(*,*) nscan,nsta,nchn,ntyp,iy,id,ih,im,ss,dt,slat,slon,
- & sdep,jys,jds,jhs,jms,sss
- if (ierr.ne.0) stop 'could not read event header'
- if (nscan.gt.9000) stop 'Data array too long! '
- if (nsta.eq.ksta) then
- if (nchn(1:4).eq.'MODE') icomp = 0
- if (nchn(1:3).eq.'VHZ') icomp = 0
- told = tstart
- else
- icomp = 0
- told = 0.
+c.....open file with station & channel info ---
+ write(*,*) 'green: Input dbname : ',dbin(1:lnblnk(dbin))
+ call read_site
+ call read_sitechan
+c.....select channel sequence ---
+ call select_sitechan(jdate,nchan,numchan,numsta)
+ open(12,file=fname4,status='old')
+ inchn = 1
+ ifl = 0
+c.....choose new single or triple of station(s)
+ 6 iseq = 1
+cxx
+cxx if(ifl.ge.30) goto 88
+ 66 if(inchn.gt.nchan) goto 88
+ nc(iseq) = iabs(numchan(inchn))
+ ns(iseq) = numsta(inchn)
+ inchn = inchn+1
+ iseq = iseq+1
+ if(numchan(inchn-1).le.0) goto 7
+ goto 66
+ 7 iseq = iseq-1
+c.....print station and channel info ---
+ write(*,1004) sta_site(ns(1)),lat_site(ns(1)),
+ * lon_site(ns(1)),iseq
+ 1004 format(' green: Station: ',a6,1x,f9.4,f10.4,' , Channels: ',i1)
+ do i = 1,iseq
+ write(*,*) 'green: Channel: # ',i,' ',sta_sitechan(nc(i)),
+ * chan_sitechan(nc(i)),hang_sitechan(nc(i)),
+ * vang_sitechan(nc(i))
+ enddo
+c....check channels sequence: Z,N,E ---
+ if(iseq.gt.3) then
+ write(*,*) 'WARNING: green: # of channels is stripped to 3'
+ iseq = 3
endif
- ksta = nsta
- mchn=1
- if(nchn(1:3).eq.'VHN ') mchn=2
- if(nchn(1:3).eq.'VH1 ') mchn=2
- if(nchn(1:3).eq.'VHE ') mchn=3
- if(nchn(1:3).eq.'VH2 ') mchn=3
- if(mchn.eq.2) call newsta(nsta,nchn,iy,id,t1,p1,az1)
- if(mchn.eq.3) call newsta(nsta,nchn,iy,id,t1,p1,az2)
+ if(iseq.eq.2) then
+ write(*,*) 'WARNING: green: # of channels is stripped to 1'
+ iseq = 1
+ endif
+ ang(1) = (vang_sitechan(nc(1))-90.0)/rad
+ if(abs(vang_sitechan(nc(1))).gt.0.5) then
+ write(*,*)
+ * 'WARNING: green: Channel: # ',1,' is not vertical. ',
+ * 'Sequence ignored'
+ goto 6
+ endif
+ do i = 2,iseq
+ ang(i) = hang_sitechan(nc(i))/rad
+ if(abs(vang_sitechan(nc(i))-90.0).gt.0.5) then
+ write(*,*)
+ * 'WARNING: green: Channel: # ',i,' is not horizontal. ',
+ * 'Sequence ignored'
+ goto 6
+ endif
+ enddo
+c
+c compute green functions for selected channels
+c
+ icomp = 0
+ do 90 isq=1,iseq
+c.....extract channel parameters from relation tables
+ ifl = ifl+1
+ mchn = isq
+ iy = jy
+ id = jd
+ ih = jh
+ im = jm
+ ss = sec
+ t1 = (90.0-lat_site(ns(1)))/rad
+ p1 = lon_site(ns(1))
+ if(p1.lt.0.0) p1 = 360.0+p1
+ p1 = p1/rad
- print 110,nsta,mchn
- 110 format(' Station ',a4,' Channel',i2)
-cxx call inresp(iscan)
- tstart=((((iy-jy)*366.+(id-jd))*24.+(ih-jh))*60.+(im-jm))*60.+
- +ss-sec
- if(tstart.ne.told.and.icomp.ne.0) icomp=1
if(icomp.eq.2) goto 500
len=0
if(icomp.eq.1) goto 35
- call allnsta(nsta,t1,p1)
c
c....do some trigonometry
c epicentral distance: co, si
@@ -131,18 +212,18 @@
co=c0*c1+s0*s1*cos(dp)
si=dsqrt(1.d0-co*co)
- del = atan2(si,co)
- del = del /pi*180.
- write(*,'(a,f10.3)')' Epicentral Distance : ',del
+ del = datan2(si,co)
+ del = del/pi*180.
+ write(*,'(a,f10.3)')' green: Epicentral Distance : ',del
sz=s1*sin(dp)/si
cz=(c0*co-c1)/(s0*si)
saz=-s0*sin(dp)/si
caz=(c0-co*c1)/(si*s1)
- azim = atan2(saz,caz)
+ azim = datan2(saz,caz)
azim = azim/pi*180.
- write(*,'(a,f10.3)')' Azimuth of Source : ',azim
+ write(*,'(a,f10.3)')' green: Azimuth of Source : ',azim
c2=2.d0*(cz**2-sz**2)
s2=8.d0*cz*sz
@@ -207,64 +288,51 @@
call prop(omt,a0t,dt,tstart,nmodet,len,iscan,4,prp)
61 continue
call s10t12(grn,iscan,c1,s1,c2,s2)
- call rotate(grn,iscan,caz,saz,az1,az2)
+ call rotate(grn,iscan,caz,saz,ang(2),ang(3))
c
-c....read in data record
-c
-cxx 500 call gfs_rwentry(1,hed1,dat,ifl,'r')
- 500 do ii = 1,nscan
- dat(ii) = 0.5*sin(0.31415926*(ii-1)*dt)
- enddo
-c
c....greens function longer than data ? : add flags to data !
c
- if(iscan.le.nscan) goto 75
- ip1=nscan+1
- do 70 i=ip1,iscan
- 70 dat(i)=5.e15
- 75 nscan=iscan
-c
-c....write out record
-c
- do 77 i=1,9
- 77 hed1(10+i)=alocal(i)
- igr = 2*ifl - 1
- call prhdr(igr,hed1,itype1)
- call gfs_rwentry(2,hed1,dat,igr,'w')
-c
-c....find gaps in data
-c
- call panlsp(dat,nscan,n1,n2,np)
- indx=0
+cxx 500 nscan=iscan
+ 500 indx=0
if(mchn.eq.3) indx=6*iscan
-c
-c....convolve green's function by instrument response
-c
-cxx call grnflt(grn(indx+1),iscan,6)
len=6*iscan
-c if(dcg(mchn).eq.1.0) goto 87
c
-c....multiply green's function by calibration factor
-c
-c do 86 i1=1,len
-c jlong=indx+i1
-c 86 grn(jlong)=grn(jlong)*dcg(mchn)
-c
c....write out green's functions
c
- 87 nscan = 6 * iscan
- igr = 2*ifl
- call prhdr(igr,hed1,itype1)
- call gfs_rwentry(2,hed1,grn(indx+1),igr,'w')
+ nscan = 6 * iscan
+ write(*,1000) ifl,sta_sitechan(nc(1)),chan_sitechan(nc(isq)),
+ * jys,jds,jhs,jms,sss,dt,nscan
+ 1000 format(' green: ',i4,1x,a6,1x,a8,i4,1x,i3,1x,i2,':',i2,':',
+ * f6.3,1x,f8.3,1x,i6)
+c write wfdisc relation with green functions
+ call default_wfdisc(1)
+ nrowwfdisc = 1
+ sta_wfdisc(1) = sta_sitechan(nc(1))
+ chan_wfdisc(1) = chan_sitechan(nc(isq))
+ chanid_wfdisc(1) = nc(isq)
+ time_wfdisc(1) = time
+ wfid_wfdisc(1) = ifl
+ jdate_wfdisc(1) = jdate
+ endtime_wfdisc(1) = time+(nscan-1)/dt
+ nsamp_wfdisc(1) = nscan
+ samprate_wfdisc(1) = 1.0/dt
+ segtype_wfdisc(1) = 'g'
+ foff_wfdisc(1) = 0
+ dir_wfdisc(1) = '.'
+ write(efname,'("g.",i5)') ifl
+ do i = 1,7
+ if(efname(i:i).eq.' ') efname(i:i) = '0'
+ enddo
+ dfile_wfdisc(1) = efname
+ call put_wfdisc(1,nscan,grn(indx+1),ierr)
+ call write_wfdisc
90 continue
+ goto 6
+ 88 continue
close(12)
-cxx 99 call gfs_close(1)
99 continue
- call gfs_close(2)
- stop
end
-
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
@@ -317,6 +385,7 @@
return
end
+
c***************************************************************
c source sub reads eigen functions from flat databases.
c List of dbnames are defined in dbnames.dat file. source read
@@ -325,6 +394,7 @@
subroutine source(fmin,fmax,lmax,nmodes,nmodet)
implicit none
include "green.h"
+ include "fdb/fdb_eigen.h"
integer*4 lmax,nmodes,nmodet
real*4 fmin,fmax
c --- common blocks -------------------------------
@@ -337,37 +407,21 @@
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 pi2,fot,vn,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 npts,nrecl,ieig,idat,ierr,ifl,i,is,j,js
integer*4 ll,lll,m,l,n
+ character*2 tendia,endian
c ---
character*64 dir
character*256 dbname
- real*4 r(mk),buf(6,mk)
+ real*4 fnl(2),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'/
+
pi2 = datan(1.0d0)*8.0d0
nmodes = 0
nmodet = 0
@@ -377,37 +431,21 @@
c===================================================================
c Main loop by dbase names
c===================================================================
- open(7,file='dbnames.dat',status='old')
+ open(7,file=fname3,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 open_eigen(dbname,ieig,idat,nrecl,dir,'r',ierr)
call read_eigen(ieig,ierr)
nrecl = (ncol_eigen*nraw_eigen+npar_eigen)*4
-c
-c....read in parameters and radial grid
-c
- ifl=1
- 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.
- r0=1.-(d0/rs)
- do 5 i=1,npts
- if(r(i).lt.r0)is=i
- 5 continue
- x1=r(is)
- js=is+1
- x2=r(js)
+ call close_eigen(ieig,idat)
c
c....read in spheroidal and toroidal modes
c....in frequency band fmin < f < fmax
c
- call open_eigen(dbname,nrecl,ieig,idat,dir,'r',ierr)
- call read_eigen(ieig,ierr)
+ ifl = 0
+ call open_eigen(dbname,ieig,idat,nrecl,dir,'r',ierr)
10 ifl = ifl+1
call read_eigen(ieig,ierr)
@@ -418,14 +456,45 @@
if (w.lt.fmin) goto 10
if (w.le.fmax) goto 11
goto 10
- 11 read(idat,rec=eigid_eigen) n,l,w,q,
- + ((buf(ll,lll),ll=1,ncol_eigen),lll=1,nraw_eigen)
+ 11 read(idat,rec=eigid_eigen) n,l,w,q,rn,vn,accn,
+ + (r(lll),(buf(ll,lll),ll=1,ncol_eigen-1),lll=1,nraw_eigen)
+c swap bytes if necessary
+ tendia = endian()
+ if(datatype_eigen.ne.tendia) then
+ write(*,*) 'XXX'
+ call swap1(fnl,4,2)
+ call swap1(w,4,1)
+ call swap1(q,4,1)
+ call swap1(rn,4,1)
+ call swap1(vn,4,1)
+ call swap1(accn,4,1)
+ call swap1(r,4,mk)
+ call swap1(buf,4,6*mk)
+ endif
npts = nraw_eigen
- if (typeo_eigen.eq.'S') then
+ wn = vn/rn
c
+c find radius interpolation points
+c
+ rs=rn/1000.
+ r0=1.-(d0/rs)
+ do 5 i=1,npts
+ if(r(i).lt.r0)is=i
+ 5 continue
+ x1=r(is)
+ js=is+1
+ x2=r(js)
+
+ if (typeo_eigen.eq.'S') then
+C
c get source scalars for S mode
c
nmodes=nmodes+1
+ if(nmodes.ge.meig) then
+ write(*,*) 'ERR012: green: # sph. modes in band exceed ',
+ * 'max allowed number ',meig
+ stop ' '
+ endif
om(nmodes)=w
a0(nmodes)=q
lmax = max0(l,lmax)
@@ -466,6 +535,11 @@
c get source scalars for T mode
c
nmodet=nmodet+1
+ if(nmodet.ge.meig) then
+ write(*,*) 'ERR012: green: # tor. modes in band exceed ',
+ * 'max allowed number ',meig
+ stop ' '
+ endif
omt(nmodet)=w
a0t(nmodet)=q
lmax = max0(l,lmax)
@@ -488,77 +562,18 @@
vecnlt(4,nmodet)=-4.*aw*wr
c print 801,n,l,(vecnlt(i,nmodet),i=3,4)
801 format(2i4,6e15.7)
+c 665 continue
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
+ write(*,*) 'green: # sph. modes in band =',nmodes,
+ * ' must be .le. ',meig
+ write(*,*) 'green: # tor. modes in band =',nmodet,
+ * ' must be .le. ',meig
return
end
-
- subroutine grnflt(grn,iscan,ifilt)
-c
-c....multiply green's functions by instrument response
-c
- real*8 s,sh
- complex resp,zed
- common/respX/resp(4501)
- common/workX/s(9002)
- common/myhed/nscan,nsta,nchn,ntyp,iy,id,ih,im,ss,dt,spare(20)
- dimension grn(*)
- nh=iscan/2
- nha=nh+1
- sh=dt
- do 10 j=1,ifilt
- k1=(j-1)*iscan
- do 15 i=1,iscan
- 15 s(i) = grn(k1+i)
- isn = -1
- call fftldp(s,iscan,isn,sh,kerr)
- if (kerr.ne.0) stop 'fft failed'
- i = 1
- do 20 jj=2,nha
- i=i+2
- zed = cmplx(s(i),s(i+1)) * resp(jj)
- s(i) = real(zed)
- 20 s(i+1)= aimag(zed)
- s(1)=0.0d0
- s(2)=0.0d0
- isn=-2
- call fftldp(s,iscan,isn,sh,kerr)
- if (kerr.ne.0) stop 'fft failed'
- do 25 i=1,iscan
- 25 grn(k1+i)=s(i)
- 10 continue
- return
- end
-
- subroutine inresp(iscan)
-c
-c....instrument response
-c
- real*8 df,pi,sh
-cxx (MB) complex evresh,resp,z
- complex evresh,resp
- common/respX/resp(4501)
- common/myhed/nscan,nsta,nchn,ntyp,iy,id,ih,im,ss,dt,spare(20)
- data pi/3.14159265358979d0/
- character*4 ntyp,nsta,nchn
- call pznall(nsta,nchn,ntyp,iy,id,istat)
- nh=iscan/2
- nha=nh+1
- sh=dt
- df=2.d0*pi/(sh*iscan)
- resp(1)=(1.,0.)
- do i=2,nha
- wr = (i-1)*df
- resp(i) = evresh(wr)/cmplx(0.,wr)
- enddo
- return
- end
-
subroutine s4t6(grn,nscan,c1,s1,c2,s2)
implicit real*8(a-h,o-z)
real*4 grn(nscan,*)
@@ -709,98 +724,3 @@
10 elpmm1=elpmm1+1d0
return
end
- subroutine panlsp(a,len,n1,n2,np)
-c find gaps of bad data, previously flagged, zero them and index them.
-c the resulting time series consists of np panels of good data indexed
-c from n1(i) to n2(i), i=1 to np.
- dimension a(*),n1(*),n2(*)
- data zero/0./,badun/3.5e15/
- id=2
- np=0
- do 30 i=1,len
- go to (15,20),id
- 15 if(a(i).lt.badun) goto 30
- a(i)=zero
- n2(np)=i-1
- id=2
- go to 30
- 20 if(a(i).gt.badun) goto 25
- np=np+1
- n1(np)=i
- id=1
- goto 30
- 25 a(i)=zero
- 30 continue
- if(id.eq.2) return
- n2(np)=len
- return
- end
-
-
- integer function lenb(pathname)
-c counts number of non-blank characters in the name and left justifies
-c if there are leading blanks
- character*(*) pathname
- lenb=0
- jchar=len(pathname)
- do 5 i=1,jchar
- k=i
- if(pathname(i:i).ne.' ') goto 10
- 5 continue
- return
- 10 do 20 i=k,jchar
- if(pathname(i:i).eq.' ') goto 30
- 20 lenb=lenb+1
- 30 if(k.eq.1) return
- do 40 i=1,lenb
- j=k+i-1
- 40 pathname(i:i)=pathname(j:j)
- return
- end
-
- subroutine newsta(nsta,nchn,iy,id,th,ph,azim)
- character*90 filename
- character*4 tsta,nchn,nsta
- character*3 tchn
- integer ty1,td1,ty2,td2,nix
- parameter(ix=100000)
- common/staXXX/nix,tsta(ix),tchn(ix),
- & ty1(ix),td1(ix),ty2(ix),td2(ix),
- & tlat(ix),tlon(ix),tele(ix),azi(ix)
- DATA RAD/.0174532/
- data iflag/1/
- if(iflag.eq.1) then
-cxx (MB) filename='/Users/guy/resp.dir/LIST.STA'
- filename='LIST.STA'
- open(99,file=filename(1:lenb2(filename)))
- n=1
- 98 read(99,96,end=97) tsta(n),tchn(n),ty1(n),td1(n),
- & ty2(n),td2(n),tlat(n),tlon(n),tele(n),azi(n)
-cxx (MB) 96 format(a4,1x,a3,2x,i4,1x,i3,1x,i4,1x,i3,1x,f10.5,1x,f10.5,1x,
- 96 format(a4,2x,a3,2x,i4,1x,i3,1x,i4,1x,i3,1x,f10.5,1x,f10.5,1x,
- + f6.0,1x,f10.5)
- n=n+1
- go to 98
- 97 close(99)
- nix=n-1
- if(nix.gt.ix) stop 'resize dimension of newsta'
- iflag=0
- end if
- do 94 i=1,nix
- if(nsta.ne.tsta(i)) goto 94
- if(nchn(1:3).ne.tchn(i)) goto 94
- if(iy.lt.ty1(i)) goto 94
- if(iy.gt.ty2(i)) goto 94
- if(iy.eq.ty1(i).and.id.lt.td1(i)) goto 94
- if(iy.eq.ty2(i).and.id.ge.td2(i)) goto 94
- th=(90.-tlat(i))*rad
- ph=tlon(i)
- if(ph.lt.0.) ph=360+ph
- ph=ph*rad
- azim=azi(i)*rad
- return
- 94 continue
- print *, 'station not found in newsta',nsta
- return
- end
-
Modified: seismo/1D/syndat/trunk/green.h
===================================================================
--- seismo/1D/syndat/trunk/green.h 2006-10-28 23:23:24 UTC (rev 5113)
+++ seismo/1D/syndat/trunk/green.h 2006-10-29 07:09:11 UTC (rev 5114)
@@ -12,7 +12,11 @@
parameter (ml=6000)
c ---
integer*4 meig
- parameter (meig=100000)
+ parameter (meig=200000)
c ---
integer*4 mseis
parameter (mseis=30000)
+c ---
+c coomon block /names/ stores in/out file names
+ character*256 fname1,fname2,fname3,fname4,fname5
+ common/names/fname1,fname2,fname3,fname4,fname5
Modified: seismo/1D/syndat/trunk/minos_bran.f
===================================================================
--- seismo/1D/syndat/trunk/minos_bran.f 2006-10-28 23:23:24 UTC (rev 5113)
+++ seismo/1D/syndat/trunk/minos_bran.f 2006-10-29 07:09:11 UTC (rev 5114)
@@ -93,14 +93,20 @@
print *,'input model file:'
read(5,100) filnam
100 format(a256)
+c MB added one line below
+ print *,filnam(1:lnblnk(filnam))
open(7,file=filnam,status='old',form='formatted',iostat=iret)
print *,'output file:'
read(5,100) filnam
+c MB added one line below
+ print *,filnam(1:lnblnk(filnam))
open(8,file=filnam,form='formatted',iostat=iret)
call model(7,8)
close(7)
print *,'eigenfunction file (output):'
read(5,100) filnam
+c MB added one line below
+ print *,filnam(1:lnblnk(filnam))
ifreq=1
if(filnam(1:4).eq.'none') ifreq=0
open(3,file=filnam,form='unformatted',iostat=iret)
@@ -123,6 +129,8 @@
stepf=1.d0
print *,'enter eps and wgrav'
read(*,*) eps,wgrav
+c MB added one line below
+ print *,eps,wgrav
eps=max(eps,1.d-12)
eps1=eps
eps2=eps
@@ -135,9 +143,13 @@
call steps(eps)
print *,'enter jcom (1=rad;2=tor;3=sph;4=ictor)'
read(*,*) jcom
+c MB added one line below
+ print *,jcom
if(jcom.lt.1.or.jcom.gt.4) jcom=3
print *,'enter lmin,lmax,wmin,wmax,nmin,nmax'
read(*,*) lmin,lmax,wmin,wmax,normin,normax
+c MB added one line below
+ print *,lmin,lmax,wmin,wmax,normin,normax
wmin=max(wmin,0.1d0)
wmin=wmin*cmhz
wmax=wmax*cmhz
Modified: seismo/1D/syndat/trunk/syndat.f
===================================================================
--- seismo/1D/syndat/trunk/syndat.f 2006-10-28 23:23:24 UTC (rev 5113)
+++ seismo/1D/syndat/trunk/syndat.f 2006-10-29 07:09:11 UTC (rev 5114)
@@ -1,168 +1,183 @@
- program syndat
c
c.....syndat makes up synthetic seismograms from green's functions
c
- implicit real*8(a-h,o-z)
- real*4 ss,dt,a,flag,hed1(30),spare,zero
- character*1l2,k2,k3
- common/rcrds/a(108000)
- common/myhed/nscan,nsta,nchn,ntyp,iy,id,ih,im,ss,dt,spare(20)
- common/rhs/rh1(10000)
- dimension sol(6),n1(100),n2(100)
- equivalence (nscan,hed1)
+ program syndat
+ implicit none
+ include 'fdb/fdb_wfdisc.h'
+ include 'fdb/fdb_io.h'
+ include 'green.h'
+c----
+ real*4 a
+ common/rcrds/a(6*mseis)
+ real*8 rh1
+ common/rhs/rh1(mseis)
+ real*4 dt,fnorm
+ integer*4 nscan, itptens, jys, jds, jhs, jms
+ integer*4 l, i, nfw, ifl, ierr, j, ib, nh
+ integer*4 iscan, iscan2, np1, isn, kerr,lnblnk
+ real*8 pi,sss,slat,slon,sdep,tconst,dm0,fscal,ww
+ real*8 si,df,wr,con,arg,sinc2,frp,dati,fip,datr,sol(6)
+ real*8 p1(3),p2(3),un(3),us(3)
+ character*256 efname
+ character*8 ename
+ character*2 endian
+c--
data pi/3.14159265358979d0/
- data itype1,zero,flag/1,0.,5.e15/
- open(7,file='locale.cmt.big')
- 21 write(*,101)
- 101 format('Do you want to use moment tensor in locale ? (y/n) : ',$)
- read(*,'(a1)',err=21)l2
- if(l2.ne.'y'.and.l2.ne.'Y') goto 1
- goto 22
- 1 write(*,102)
- 102 format('Do you want to use strike,dip and slip ? (y/n) : ',$)
- read(*,'(a1)',err=1)k2
- if(k2.ne.'y'.and.k2.ne.'Y') goto 10
- call fault(sol)
- goto 20
- 10 write(*,103)
- 103 format('enter moment tensor elements f(x6) : ',$)
- read(*,*,err=10)sol
- 20 write(*,104)
- 104 format('enter time constant : ',$)
- read(*,*,err=20)tconst
- 22 call gfs_opena(1,'greens function file name : ',itype1,jret1)
- call gfs_opena(2,'output file name : ',itype1,jret2)
- call gfs_erase(2)
- 2 write(*,105)
- 105 format('Do you want synthetics only ? (y/n) : ',$)
- read(*,'(a1)',err=2)k2
- 3 write(*,106)
- 106 format('Do you want the synthetics to be gap free ? (y/n) : ',$)
- read(*,'(a1)',err=3)k3
- ifl=1
- call gfs_rwdir(1,hed1,ifl,'r',ierr)
- if (ierr.ne.0) goto 50
- if(l2.ne.'y'.and.l2.ne.'Y') goto 34
+ data fnorm/1.0e-27/
c
-c.... locate source parameters
+c read input parameters
c
- 31 read(7,100,end=32) jy,jd,jh,jm,sec,t0,p0,d0,fm,
- + (sol(i),i=1,6),tconst,fscal
- 100 format(i4,i4,2i3,f5.1,2f9.2,f4.0,g11.3,7f7.2,g11.3)
- if(iy.eq.jy.and.id.eq.jd.and.jh.eq.ih) goto 33
- goto 31
- 32 print*,'***** event not found in locale! *****'
- close(7)
- l2 = 'n '
- goto 1
- 33 close(7)
- fm=fm*1.e-27
- fscal=fscal*1.e-27
- do i=1,6
- sol(i)=sol(i)*fscal
- enddo
+c....read in file name within event and moment tensor info
+ write(*,*) '============== Program syndat =================='
+ write(*,*) 'enter input CMT file name:'
+ read(*,'(a256)') fname1
+ write(*,*) fname1(1:lnblnk(fname1))
+c.... setup tensor type: moment/Mo,nodal2/Mo,nodal2
+ write(*,*) 'enter tensor type: 0 - moment, ',
+ * '1 - nodal plane 1, 2 - nodal plane 2'
+ read(*,*) itptens
+ write(*,*) itptens
+c.... i/o database name
+ write(*,*) 'enter input dbname'
+ read(*,'(a256)') dbin
+ write(*,*) dbin(1:lnblnk(dbin))
+ write(*,*) 'enter output dbname'
+ read(*,'(a256)') dbout
+ write(*,*) dbout(1:lnblnk(dbout))
c
-c....loop over records
-c
- 34 isy = 0
- do 30 ifl=1,jret1,2
- igr = ifl + 1
+c read source and moment
c
-c....read in record
+ open(12,file=fname1,status='old')
+ read(12,*) ename,jys,jds,jhs,jms,sss,slat,slon,sdep,dt,
+ * tconst,dM0,(sol(l),l=1,6),fscal,(p1(l),l=1,3),
+ * (p2(l),l=1,3)
+ close(12)
+ write(*,1001) ename,jys,jds,jhs,jms,sss,slat,slon,sdep,dt,
+ * tconst,dM0,(fscal*sol(l),l=1,6),(p1(l),l=1,3),
+ * (p2(l),l=1,3)
+ 1001 format('syndat: CMT solution:',/
+ * 'syndat: Event: ',a8,i5,i4,i3,':',i3,':',f5.1,', lat/lon = ',
+ * f8.2,f9.2,',',/
+ * 'syndat: depth = ',f6.1,' km, step = ',f8.3,
+ * ' sec, Duration = ',f5.1, ' sec,',/
+ * 'syndat: M0 = ',e11.3,','/
+ * 'syndat: comp: ',6e11.3,/
+ * 'syndat: plane1: ',3f5.0,' plane2: ',3f5.0)
+ if(itptens.eq.0) then
+ do i = 1,6
+ sol(i) = sol(i)*fscal*fnorm
+ enddo
+ else if(itptens.eq.1) then
+ call udc(p1(1),p1(2),p1(3),un,us,sol)
+ do i = 1,6
+ sol(i) = sol(i)*dM0*fnorm
+ enddo
+ else
+ call udc(p2(1),p2(2),p2(3),un,us,sol)
+ do i = 1,6
+ sol(i) = sol(i)*dM0*fnorm
+ enddo
+ endif
+ write(*,1002) (sol(l)/fnorm,l=1,6)
+ 1002 format('syndat: Selected moment tensor components:',/
+ * 'syndat:',6e11.3,//'syndat: Synthetic waveforms: ')
c
- call gfs_rwentry(1,hed1,a,ifl,'r')
- call panlsp(a,nscan,n1,n2,np)
- if(k2.eq.'y'.or.k2.eq.'Y') goto 60
- call fgapsp(a,nscan,n1,n2,np,flag)
-c
-c....write record to disk
-c
- isy = isy + 1
- call gfs_rwentry(2,hed1,a,isy,'w')
-c
-c....read in green's functions
-c
- 60 nold = nscan
- call gfs_rwentry(1,hed1,a,igr,'r')
- len = nscan
- nscan = nold
- do 71 i=1,nscan
- 71 rh1(i)=0.d0
- do 72 j=1,6
- ib=(j-1)*nscan
- if(k3.ne.'y'.and.k3.ne.'Y')
- * call fgapsp(a(ib+1),nscan,n1,n2,np,zero)
- do 72 i=1,nscan
- 72 rh1(i)=rh1(i)+a(ib+i)*sol(j)
- do 81 i=1,nscan
- 81 a(i)=rh1(i)
-
- if(tconst.le.0.d0) goto 90
-
- si=dt
- nh=(nscan+1)/2
- call facdwn(nh)
- iscan=2*nh
- iscan2=iscan+2
- np1=nscan+1
- do 82 i=np1,iscan2
- 82 a(i)=0.
- isn=-1
- call fftl(a,iscan,isn,kerr)
- df=2.d0*pi/(iscan*si)
- j=0
- do 83 i=3,iscan2,2
- j=j+1
- wr=j*df
- con=0.5d0*tconst*wr
- arg=0.5d0*con
- sinc2=(dsin(arg)/arg)**2
- frp=sinc2*dcos(con)
- fip=-sinc2*dsin(con)
- datr=a(i)*frp-a(i+1)*fip
- dati=a(i)*fip+a(i+1)*frp
- a(i)=datr
- 83 a(i+1)=dati
- isn=-2
- call fftl(a,iscan,isn,kerr)
- 90 if(k3.ne.'y'.and.k3.ne.'Y') call fgapsp(a,nscan,n1,n2,np,flag)
-c
-c....write out synthetics
-c
- isy = isy + 1
-c....force SNR flag of synthetic to be zero
- hed1(30)=0.
- call gfs_rwentry(2,hed1,a,isy,'w')
- call prhdr(isy,hed1,itype1)
+c....loop over records
+c
+c.....read input wfdisc relation
+ call read_wfdisc
+ nfw = nrowwfdisc
+ do 30 ifl=1,nfw
+ nscan =nsamp_wfdisc(ifl)
+ write(*,1000) ifl,sta_wfdisc(ifl),chan_wfdisc(ifl),
+ * 1.0/samprate_wfdisc(ifl), nscan/6
+ 1000 format('syndat: ',i4,1x,a6,1x,a8,1x,f8.3,1x,i7)
+ call get_wfdisc(ifl,nscan,a,ierr)
+ nscan = nscan/6
+ do i=1,nscan
+ rh1(i) = 0.d0
+ enddo
+ do j = 1,6
+ ib=(j-1)*nscan
+ do i = 1,nscan
+ rh1(i)=rh1(i)+a(ib+i)*sol(j)
+ enddo
+ enddo
+ do i=1,nscan
+ a(i)=rh1(i)
+ enddo
+
+ if(tconst.le.0.d0) goto 90
+
+ si=dt
+ nh=(nscan+1)/2
+ call facdwn(nh)
+ iscan=2*nh
+ iscan2=iscan+2
+ np1=nscan+1
+ do i=np1,iscan2
+ a(i)=0.
+ enddo
+ isn=-1
+ call fftl(a,iscan,isn,kerr)
+ df=2.d0*pi/(iscan*si)
+ j=0
+ do i=3,iscan2,2
+ j=j+1
+ wr=j*df
+ con=0.5d0*tconst*wr
+ arg=0.5d0*con
+ sinc2=(dsin(arg)/arg)**2
+ frp=sinc2*dcos(con)
+ fip=-sinc2*dsin(con)
+ datr=a(i)*frp-a(i+1)*fip
+ dati=a(i)*fip+a(i+1)*frp
+ a(i)=datr
+ a(i+1)=dati
+c convert from acceleration to velocity
+ ww = a(i)
+ a(i) = a(i+1)/wr
+ a(i+1) = -ww/wr
+ enddo
+ isn=-2
+ call fftl(a,iscan,isn,kerr)
+ 90 continue
+c write new wfdisc relation
+ nrowwfdisc = 1
+ sta_wfdisc(1) = sta_wfdisc(ifl)
+ chan_wfdisc(1) = chan_wfdisc(ifl)
+ time_wfdisc(1) = time_wfdisc(ifl)
+ wfid_wfdisc(1) = ifl
+ jdate_wfdisc(1) = jdate_wfdisc(ifl)
+ endtime_wfdisc(1) = time_wfdisc(ifl)+(nscan-1)/
+ * samprate_wfdisc(ifl)
+ nsamp_wfdisc(1) = nscan
+ chanid_wfdisc(1) = chanid_wfdisc(ifl)
+ samprate_wfdisc(1) = samprate_wfdisc(ifl)
+ segtype_wfdisc(1) = 'w'
+ foff_wfdisc(1) = 0
+ dir_wfdisc(1) = '.'
+ write(efname,'("w.",i5)') ifl
+ do i = 1,7
+ if(efname(i:i).eq.' ') efname(i:i) = '0'
+ enddo
+ dfile_wfdisc(1) = efname
+ calper_wfdisc(1) = 20.0
+ calib_wfdisc(1) = 1.0
+ datatype_wfdisc(1) = endian()
+ do i = 1,nscan
+ a(i) = a(i)*1.0e9
+ enddo
+ call put_wfdisc(1,nscan,a,ierr)
+ call write_wfdisc
30 continue
- 50 call gfs_close(1)
- call gfs_close(2)
- stop
end
-
- subroutine fault(sol)
- implicit real*8(a-h,o-z)
- common/faultX/un(3),us(3)
- dimension sol(*)
- 1 write(*,100)
- 100 format('enter strike, dip and slip [deg] : ',$)
- read(*,*,err=1)sig,del,gam
- call udc(sig,del,gam,un,us,sol)
- 2 write(*,101)
- 101 format('enter moment : ',$)
- read(*,*,err=2)fmom
- do 5 i=1,6
- 5 sol(i)=fmom*sol(i)
- return
- end
-
-
- subroutine udc(sig,del,gam,un,us,f)
-c
+c*********************************************************************
c udc computes the unit normal, unit slip and unit moment tensor
c given strike(sig), dip(del) and slip(gam) in degrees.
-c
+c*********************************************************************
+ subroutine udc(sig,del,gam,un,us,f)
implicit real*8(a-h,o-z)
dimension un(*),us(*),f(*)
data rad/57.29578/
@@ -189,60 +204,3 @@
f(6)=un(2)*us(3)+us(2)*un(3)
return
end
-
-
-
- subroutine panlsp(a,len,n1,n2,np)
-c find gaps of bad data, previously flagged, zero them and index them.
-c the resulting time series consists of np panels of good data indexed
-c from n1(i) to n2(i), i=1 to np.
- dimension a(*),n1(*),n2(*)
- data zero/0./,badun/3.5e15/
- id=2
- np=0
- do 30 i=1,len
- go to (15,20),id
- 15 if(a(i).lt.badun) goto 30
- a(i)=zero
- n2(np)=i-1
- id=2
- go to 30
- 20 if(a(i).gt.badun) goto 25
- np=np+1
- n1(np)=i
- id=1
- goto 30
- 25 a(i)=zero
- 30 continue
- if(id.eq.2) return
- n2(np)=len
- return
- end
- subroutine fgapsp(a,len,n1,n2,np,flag)
-c fgapsp flags bad gaps in a time series with flag(usually zero
-c or 5.e15). fgapsp presupposes a call to panlsp to put
-c panel structure in n1 and n2.
- dimension a(*),n1(*),n2(*)
- if(np.gt.0) goto 15
- do 10 i=1,len
- 10 a(i) = flag
- return
- 15 if(n1(1).eq.1) go to 25
- mm1=n1(1)-1
- do 20 i=1,mm1
- 20 a(i)=flag
- 25 if(n2(np).eq.len) go to 35
- mp1=n2(np)+1
- do 30 i=mp1,len
- 30 a(i)=flag
- 35 if(np.eq.1) return
- nbad=np-1
- do 45 nb=1,nbad
- m1=n2(nb)+1
- m2=n1(nb+1)-1
- if(m2.lt.m1) go to 45
- do 40 m=m1,m2
- 40 a(m)=flag
- 45 continue
- return
- end
Added: seismo/1D/syndat/trunk/time/time.f
===================================================================
--- seismo/1D/syndat/trunk/time/time.f 2006-10-28 23:23:24 UTC (rev 5113)
+++ seismo/1D/syndat/trunk/time/time.f 2006-10-29 07:09:11 UTC (rev 5114)
@@ -0,0 +1,104 @@
+c
+c convert epoch time to human time
+c
+ subroutine epochtoh(t, year,doy,hour,miin,sec)
+ implicit none
+ real*8 t,sec,fsec
+ integer*4 year,doy,hour,miin
+ integer*4 idate, itime, irest,iysupp,idsupp
+
+c... date parts
+ idate = t/86400.0
+ iysupp = idate/365
+ idsupp = iysupp*365+(iysupp+1)/4
+ if(idate.lt.idsupp) iysupp = iysupp-1
+ idsupp = iysupp*365+(iysupp+1)/4
+c... extract year
+ year = 1970+iysupp
+c... extract doy
+ doy = idate - idsupp+1
+c... time part
+ fsec = t-idate*86400.0d0
+ irest = fsec
+ fsec = fsec - irest
+ itime = irest
+ irest = itime/60
+c... extract seconds
+ sec = itime - irest*60 + fsec
+ itime = irest
+ irest = itime/60
+c... extract mininutes
+ miin = itime - irest*60
+c... extract hours
+ hour = irest
+ return
+ end
+c
+c convert human time to epoch time
+c
+ real*8 function htoepoch (year,doy,hour,miin,sec)
+ implicit none
+ integer*4 year,doy,hour,miin
+ integer*4 vis, ydelay
+ real*8 sec
+ ydelay = year - 1970
+ vis = (ydelay+1)/4
+ htoepoch = (365.0d0*ydelay+vis+doy-1)*86400.0d0+
+ * (hour*60.0+miin)*60.0d0+sec
+ return
+ end
+c
+c... convert month and day to day of year - doy
+c
+ integer*4 function mtodoy(year,mon,day)
+ implicit none
+ integer*4 year,mon,day,i,jday,days(12)
+ data days/31,28,31,30,31,30,31,31,30,31,30,31/
+c ---
+ days(2) = 28
+ if(mod(year,4).eq.0) days(2) = 29
+ jday = 0
+ if(mon.gt.1) then
+ do i = 1,mon-1
+ jday = jday+days(i)
+ enddo
+ endif
+ mtodoy = jday+day
+ return
+ end
+c
+c... convert day of year to month and day
+c
+ subroutine doytom(year,doy,mon,day)
+ implicit none
+ integer*4 year,doy,mon,day,i,days(12)
+ data days/31,28,31,30,31,30,31,31,30,31,30,31/
+c ---
+ days(2) = 28
+ if(mod(year,4).eq.0) days(2) = 29
+ mon = 1
+ day = doy
+ do i = 1,11
+ if(day.le.days(i)) goto 2
+ mon = mon+1
+ day = day-days(i)
+ enddo
+ 2 continue
+ return
+ end
+c
+c get local time in form mm/dd/yy-hh:mm:ss
+c
+ character*17 function loctime()
+ implicit none
+ integer*4 itime,ia(9),time,i
+ itime=time()
+ call ltime(itime,ia)
+ write(loctime,1000) ia(5)+1,ia(4),mod(1900+ia(6),100),
+ + ia(3),ia(2),ia(1)
+ do i = 1,17
+ if(loctime(i:i).eq.' ') loctime(i:i) = '0'
+ enddo
+ return
+ 1000 format(i2,'/',i2,'/',i2,'-',i2,':',i2,':',i2)
+ end
Added: seismo/1D/syndat/trunk/utils/endi.c
===================================================================
--- seismo/1D/syndat/trunk/utils/endi.c 2006-10-28 23:23:24 UTC (rev 5113)
+++ seismo/1D/syndat/trunk/utils/endi.c 2006-10-29 07:09:11 UTC (rev 5114)
@@ -0,0 +1,51 @@
+#include <stdio.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <unistd.h>
+#include <stdlib.h>
+
+/* Function prototype */
+void swapn(unsigned char *b, int N, int n);
+
+int main(int iargc,char *argv[])
+{
+ int i,n,wdt;
+ char *p;
+ unsigned char *b;
+ struct stat st;
+ FILE *in,*out;
+
+ if(iargc < 3) {
+ printf("Usage: endi nwidth file1 [ file2 ... filen]\n");
+ return 1;
+ }
+ for(i = 2; i < iargc; i++) {
+ p = argv[i];
+ if(stat(p,&st) != 0) {
+ printf("File %s does no exists\n",p);
+ continue;
+ }
+ wdt = atoi(argv[1]); /* get swapping width */
+ n = st.st_size/wdt;
+ if((b = (unsigned char *)malloc(st.st_size)) == NULL) {
+ printf(" Can not allocate memory\n");
+ }
+ /* open input file */
+ if((in = fopen(p,"r")) == NULL) {
+ fprintf(stderr,"File %s does nor exist.\n",p);
+ return 1;
+ }
+ fread(b,st.st_size,1,in);
+ fclose(in);
+ swapn(b,wdt,n);
+ /* open output file */
+ if((out = fopen(p,"w")) == NULL) {
+ fprintf(stderr,"File %s does nor exist.\n",p);
+ return 1;
+ }
+ fwrite(b,st.st_size,1,in);
+ fclose(out);
+ free(b);
+ }
+ return 0;
+}
Added: seismo/1D/syndat/trunk/utils/preigen.f
===================================================================
--- seismo/1D/syndat/trunk/utils/preigen.f 2006-10-28 23:23:24 UTC (rev 5113)
+++ seismo/1D/syndat/trunk/utils/preigen.f 2006-10-29 07:09:11 UTC (rev 5114)
@@ -0,0 +1,106 @@
+c read eigen functions
+ 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
+ character*64 dir
+ character*256 fout
+ real*4 rout(mk),buf(6,mk)
+ real*4 pi2,rn,vn,accn,ts,t1,t2
+ real*4 t3,ww,qq
+ integer*4 narg,iargc,nrecl,ieig,idat,ierr
+ integer*4 indreq,nn,ll,lll,i,j,mods
+
+ pi2 = atan(1.0)*8.0
+c get input arguments
+ narg=iargc()
+ if (narg.ne.0) then
+ write(*,*) ' Usage: peigcon '
+ call exit(1)
+ endif
+c ---
+ nrecl = 2000
+ ieig = 9
+ idat = 10
+ write(*,*) 'c Enter eigen database name'
+ read(*,*,end =98) fout
+c get record length of direct eigen file
+ call open_eigen(fout,ieig,idat,nrecl,dir,'r',ierr)
+ call read_eigen(ieig,ierr)
+ nrecl = (ncol_eigen*nraw_eigen+npar_eigen)*4
+ call close_eigen(ieig,idat)
+ 1 call open_eigen(fout,ieig,idat,nrecl,dir,'r',ierr)
+c find record by index indreq and print it out
+ write(*,*) 'c Enter mode and Period:'
+ read(*,*,end=999) mods,ts
+ indreq = 0
+ 4 indreq = indreq+1
+ call read_eigen(ieig,ierr)
+ if(ierr.ne.0) goto 99
+ if(norder_eigen.ne.mods) goto 4
+ t1 = abs(per_eigen-ts)
+ indreq = indreq+1
+ call read_eigen(ieig,ierr)
+ if(ierr.ne.0) goto 3
+ if(norder_eigen.ne.mods) goto 3
+ t2 = abs(per_eigen-ts)
+ if(t2.gt.t1) goto 3
+ 2 indreq = indreq+1
+ call read_eigen(ieig,ierr)
+ if(ierr.ne.0) goto 3
+ if(norder_eigen.ne.mods) goto 3
+ t3 = abs(per_eigen-ts)
+ if((t1-t2)*(t2-t3).le.0.0) goto 3
+ t1 = t2
+ t2 = t3
+ goto 2
+ 3 indreq = indreq-1
+ call close_eigen(ieig,idat)
+ call open_eigen(fout,ieig,idat,nrecl,dir,'r',ierr)
+ 10 call read_eigen(ieig,ierr)
+ if(eigid_eigen.ne.indreq) goto 10
+ read(idat,rec=eigid_eigen) nn,ll,ww,qq,rn,vn,accn,
+ + (rout(lll),(buf(ll,lll),ll=1,ncol_eigen-1),lll=1,nraw_eigen)
+ do i=1,nraw_eigen
+ rout(i) = (1.0-rout(i))*rn/1000.0
+ enddo
+c output data
+ write(*,*) 'c eigid= ',indreq, ', datatype= ',datatype_eigen
+ write(*,*) 'c n order= ',norder_eigen,', mode= ',typeo_eigen,
+ * ', l order= ',lorder_eigen
+ write(*,*) 'c T = ',per_eigen,', Q = ',attn_eigen, ', nrow= ',
+ * nraw_eigen,', ncol= ',ncol_eigen
+ write(*,*) 'c rn = ',rn
+ do i = nraw_eigen,1,-1
+ j = nraw_eigen+1-i
+ if(ncol_eigen.eq.3) then
+ write(*,1000) j,rout(i),(buf(lll,i),lll=1,ncol_eigen-1)
+ else
+ write(*,1001) j,rout(i),(buf(lll,i),lll=1,ncol_eigen-1)
+ endif
+ enddo
+ call close_eigen(ieig,idat)
+ goto 1
+ 98 continue
+ 99 call close_eigen(ieig,idat)
+ goto 1
+ 999 continue
+ call close_eigen(ieig,idat)
+1000 format(i4,f10.3,2e15.6)
+1001 format(i4,f10.3,6e15.6)
+ end
Added: seismo/1D/syndat/trunk/utils/swap.c
===================================================================
--- seismo/1D/syndat/trunk/utils/swap.c 2006-10-28 23:23:24 UTC (rev 5113)
+++ seismo/1D/syndat/trunk/utils/swap.c 2006-10-29 07:09:11 UTC (rev 5114)
@@ -0,0 +1,15 @@
+/*
+ * Byte swap an array of length n of N byte integer*4 or real*4 elements.
+ * A good compiler should unroll the inner loops. Letting the compiler do it
+ * gives us portability. Note that we might want to isolate the
+ * cases N = 2, 4, 8 (and 16 for long double and perhaps long long)
+ */
+#define SWAP2(a, b) { (a) ^= (b); (b) ^= (a); (a) ^= (b); }
+void swapn(unsigned char *b, int N, int n)
+{
+ int i, j;
+ for (i = 0; i < (n)*(N); i += N)
+ for (j = 0; j < N/2; j ++)
+ SWAP2(b[i + j], b[i + N - j - 1]);
+ return;
+}
Added: seismo/1D/syndat/trunk/utils/wfendi.f
===================================================================
--- seismo/1D/syndat/trunk/utils/wfendi.f 2006-10-28 23:23:24 UTC (rev 5113)
+++ seismo/1D/syndat/trunk/utils/wfendi.f 2006-10-29 07:09:11 UTC (rev 5114)
@@ -0,0 +1,50 @@
+c================================================================
+c Reverse order of bytes for whole .wfdisc relation
+c Usage: fwendi ivel ipc INFILE OUTFILE
+c ivel - =1, make multiplication by omega, =0 - No
+c ipc - =1, rdseed for PC, add f4 always, =0 - No
+c================================================================
+ implicit none
+ include "../fdb/fdb_wfdisc.h"
+ include "../fdb/fdb_io.h"
+ integer*4 i,ivel,ierr,ipc,n,narg,iargc,lnblnk
+ real*4 pi2,data(400000)
+ character*256 named,cmd
+ character*2 endian,datatype
+ logical tf
+
+ pi2 = datan(1.0d0)*8.0d0
+ narg=iargc()
+ if (narg.ne.4) STOP ' Usage: wfehdi ivel ipc INFILE OUTFILE'
+ call getarg(1,dbin)
+ read(dbin,*) ivel
+ call getarg(2,dbin)
+ read(dbin,*) ipc
+ call getarg(3,dbin)
+ call getarg(4,dbout)
+c ---
+ named = dbout(1:lnblnk(dbout))//'.wfdisc'
+ call read_wfdisc
+ do i = 1,nrowwfdisc
+ if(ipc.eq.1) datatype = datatype_wfdisc(i)
+ if(ipc.eq.1) datatype_wfdisc(i) = 'f4'
+ n = nsamp_wfdisc(i)
+ call get_wfdisc(i,n,data,ierr)
+ if(ivel.eq.1) then
+ calib_wfdisc(i) = calib_wfdisc(i)*pi2/calper_wfdisc(i)
+ endif
+ datatype_wfdisc(i) = endian()
+ if(datatype_wfdisc(i).eq.'f4') then
+ datatype_wfdisc(i) = 't4'
+ if(ipc.eq.1) datatype_wfdisc(i)=datatype
+ call swap1(data,4,n)
+ call put_wfdisc(i,n,data,ierr)
+ endif
+ enddo
+ inquire(file=named,exist=tf)
+ if(tf) then
+ cmd = 'rm -f '//named
+ call system(cmd)
+ endif
+ call write_wfdisc
+ end
More information about the cig-commits
mailing list