[cig-commits] r13405 - in seismo/3D/GRD_CMT3D: cmt3d grid3d scripts
liuqy at geodynamics.org
liuqy at geodynamics.org
Tue Nov 25 15:24:33 PST 2008
Author: liuqy
Date: 2008-11-25 15:24:33 -0800 (Tue, 25 Nov 2008)
New Revision: 13405
Added:
seismo/3D/GRD_CMT3D/scripts/mech_table.pl
Modified:
seismo/3D/GRD_CMT3D/cmt3d/Makefile
seismo/3D/GRD_CMT3D/cmt3d/cmt3d_constants.f90
seismo/3D/GRD_CMT3D/cmt3d/cmt3d_sub.f90
seismo/3D/GRD_CMT3D/cmt3d/cmt3d_sub2.f90
seismo/3D/GRD_CMT3D/cmt3d/cmt3d_sub3.f90
seismo/3D/GRD_CMT3D/cmt3d/cmt3d_variables.f90
seismo/3D/GRD_CMT3D/grid3d/GRID3D.PAR
seismo/3D/GRD_CMT3D/grid3d/grid3d_sub.f90
seismo/3D/GRD_CMT3D/grid3d/grid3d_sub2.f90
seismo/3D/GRD_CMT3D/grid3d/grid3d_sub3.f90
seismo/3D/GRD_CMT3D/grid3d/grid3d_variables.f90
seismo/3D/GRD_CMT3D/scripts/plot_grid.pl
seismo/3D/GRD_CMT3D/scripts/process_data_and_syn.pl
Log:
Update cmt3d and grid3d package
Add more auxilliary scripts
Modified: seismo/3D/GRD_CMT3D/cmt3d/Makefile
===================================================================
--- seismo/3D/GRD_CMT3D/cmt3d/Makefile 2008-11-25 22:57:08 UTC (rev 13404)
+++ seismo/3D/GRD_CMT3D/cmt3d/Makefile 2008-11-25 23:24:33 UTC (rev 13405)
@@ -23,8 +23,7 @@
$(F90) -o $(BIN_DIR)/$@ $(F90_FLAGS) $@.f90 -I $(MOD_DIR) $(CMT_OBJ) $(LIB)
$(CMT_F90_OBJ): $(OBJ_DIR)/%.o : %.f90
- $(F90) -o $@ $(F90_FLAGS) -c $*.f90 -I $(MOD_DIR)
-# @ \mv -f *.mod $(MOD_DIR)
+ $(F90) -o $@ $(F90_FLAGS) -c $*.f90 -I $(MOD_DIR) -J$(MOD_DIR)
$(CMT_F77_OBJ): $(OBJ_DIR)/%.o : %.f
$(F77) -o $@ $(F77_FLAGS) -c $*.f
Modified: seismo/3D/GRD_CMT3D/cmt3d/cmt3d_constants.f90
===================================================================
--- seismo/3D/GRD_CMT3D/cmt3d/cmt3d_constants.f90 2008-11-25 22:57:08 UTC (rev 13404)
+++ seismo/3D/GRD_CMT3D/cmt3d/cmt3d_constants.f90 2008-11-25 23:24:33 UTC (rev 13405)
@@ -46,6 +46,7 @@
! small numbers
real*8, parameter :: EPS2 = 1.0d-2
real*8, parameter :: EPS5 = 1.0d-5
+ real, parameter :: SMALL = -huge(1.0)
! io unit for parameter files
integer, parameter :: IOPAR = 40
@@ -63,7 +64,7 @@
integer, parameter :: NREGIONS = 10
! reference distance for Pnl, Rayleigh and Love wave weighting
- real*8, parameter :: REF_DIST = 100.0
+ real, parameter :: REF_DIST = 100.0
end module cmt3d_constants
Modified: seismo/3D/GRD_CMT3D/cmt3d/cmt3d_sub.f90
===================================================================
--- seismo/3D/GRD_CMT3D/cmt3d/cmt3d_sub.f90 2008-11-25 22:57:08 UTC (rev 13404)
+++ seismo/3D/GRD_CMT3D/cmt3d/cmt3d_sub.f90 2008-11-25 23:24:33 UTC (rev 13405)
@@ -23,8 +23,10 @@
read(IOPAR,'(a)') new_cmt_file
read(IOPAR,*) npar
read(IOPAR,*) ddelta,ddepth,dmoment
- dcmt_par = (/dmoment,dmoment,dmoment, dmoment,dmoment,dmoment,&
- ddepth,ddelta,ddelta, SCALE_CTIME,SCALE_HDUR/) / SCALE_PAR
+ dcmt_par = (/dble(dmoment),dble(dmoment),dble(dmoment), &
+ dble(dmoment),dble(dmoment),dble(dmoment),&
+ dble(ddepth),dble(ddelta),dble(ddelta), &
+ SCALE_CTIME,SCALE_HDUR/) / SCALE_PAR
read(IOPAR,'(a)') flexwin_out_file
@@ -116,7 +118,7 @@
read(IOWIN,'(a)') syn_file
read(IOWIN,*) nwins(i)
if (nwins(i) < 0) stop 'Check nwins(i) '
- print *, trim(data_file), ' ', trim(syn_file)
+! if (DEBUG) print *, trim(data_file), ' ', trim(syn_file)
nwin_total = nwin_total + nwins(i)
do j = 1, nwins(i)
read(IOWIN,*) tstart, tend
@@ -147,7 +149,7 @@
real*8, intent(out) :: A(npar,npar), b(npar)
integer :: ios,nf,nw,nwint,i,j
- real*8 :: tstart, tend, tshift
+ real :: tstart, tend
real*8 :: A1(npar,npar), b1(npar)
open(IOWIN,file=trim(flexwin_out_file),iostat=ios)
@@ -184,7 +186,6 @@
endif
-
end subroutine setup_matrix
!================================================================
@@ -303,8 +304,11 @@
real*8,intent(in) :: dm(npar)
integer :: ios, nf, nwint, nw, is, ie, i, j, npts, ishift, ishift_new
- real*8 :: b, dt, tstart, tend, tshift, cc, dlna, var, &
- tshift_new,cc_new,dlna_new,var_new
+ real :: b, dt, tstart, tend
+ real :: tshift, cc, dlna, v1, v2, tshift_new,cc_new,dlna_new,d1,d2
+ real :: var_all, var_all_new
+ integer :: istart, iend, istart_d, iend_d
+ integer :: istart_n, iend_n, istart_dn, iend_dn
open(IOWIN,file=trim(flexwin_out_file),iostat=ios)
@@ -315,7 +319,7 @@
read(IOWIN,*) nf
write(IOINV,*) nf
- nwint = 0
+ nwint = 0; var_all = 0.; var_all_new = 0.
do i = 1, nf
read(IOWIN,'(a)') data_file
read(IOWIN,'(a)') syn_file
@@ -330,26 +334,37 @@
read(IOWIN,*) tstart, tend
is=max(floor((tstart-b)/dt),1)
ie=min(ceiling((tend-b)/dt),npts)
- call calc_criteria(data,syn,npts,is,ie,dt,tshift,cc,dlna)
- call calc_criteria(data,new_syn,npts,is,ie,dt,tshift_new,cc_new,dlna_new)
+
if (station_correction) then
- ishift=int(tshift/dt)
- ishift_new=int(tshift_new/dt)
- var=sum((syn(is:ie)-data(is+ishift:ie+ishift))**2)/sum(data(is:ie)**2)
- var_new=sum((new_syn(is:ie)-data(is+ishift_new:ie+ishift_new))**2)/sum(data(is:ie)**2)
+ call calc_criteria(data_sngl,syn_sngl,npts,is,ie,ishift,cc,dlna)
+ call calc_criteria(data_sngl,new_syn_sngl,npts,is,ie,ishift_new,cc_new,dlna_new)
+ istart_d=max(1,is+ishift); iend_d=min(npts,ie+ishift)
+ istart_dn=max(1,is+ishift_new); iend_dn=min(npts,ie+ishift_new)
+ istart=istart_d-ishift; iend=iend_d-ishift
+ istart_n=istart_dn-ishift_new; iend_n=iend_dn-ishift_new
else
- var=sum((syn(is:ie)-data(is:ie))**2)/sum(data(is:ie)**2)
- var_new=sum((new_syn(is:ie)-data(is:ie))**2)/sum(data(is:ie)**2)
+ istart_d=is; istart=is; iend_d=ie; iend=ie
+ istart_dn=is; istart_n=is; iend_dn=ie; iend_n=ie
endif
+
+ v1=sum((syn_sngl(istart:iend)-data_sngl(istart_d:iend_d))**2)
+ v2=sum((new_syn_sngl(istart_n:iend_n)-data_sngl(istart_dn:iend_dn))**2)
+ d1=sum(data_sngl(istart_d:iend_d)**2)
+ d2=sum(data_sngl(istart_dn:iend_dn)**2)
+ var_all = var_all + v1*data_weights(nwint)*dt
+ var_all_new = var_all_new + v2*data_weights(nwint)*dt
+
write(IOINV,'(3g15.5)') tstart,tend,data_weights(nwint)
- write(IOINV,'(3x,4g15.5)') tshift, cc, dlna, var
- write(IOINV,'(3x,4g15.5)') tshift_new, cc_new, dlna_new, var_new
+ write(IOINV,'(3x,4g15.5)') ishift*dt, cc, dlna, v1/d1
+ write(IOINV,'(3x,4g15.5)') ishift_new*dt, cc_new, dlna_new, v2/d2
enddo
enddo
close(IOWIN)
close(IOINV)
if (nwint /= nwin_total) stop 'Check nwin_total in variance reduction'
+ print *, 'Total Variance reduced from ', var_all, ' to ', var_all_new, &
+ ' = ', (var_all-var_all_new)/var_all*100, ' %'
end subroutine variance_reduction
Modified: seismo/3D/GRD_CMT3D/cmt3d/cmt3d_sub2.f90
===================================================================
--- seismo/3D/GRD_CMT3D/cmt3d/cmt3d_sub2.f90 2008-11-25 22:57:08 UTC (rev 13404)
+++ seismo/3D/GRD_CMT3D/cmt3d/cmt3d_sub2.f90 2008-11-25 23:24:33 UTC (rev 13405)
@@ -42,7 +42,7 @@
npts = min(npts1, npts2)
! DEBUG
- if (DEBUG) write(*,*) 'DEBUG : b, dt, npts ', b, dt, npts
+! if (DEBUG) write(*,*) 'DEBUG : b, dt, npts ', b, dt, npts
! read event and station header parameters from observation file
call getfhv('evla', evla, nerr)
@@ -60,8 +60,8 @@
stop 'We only deal with Z, R, and T components at the moment'
! CHT: why does this output as: FMP BHR CI(013, 064)?
- if (DEBUG) write(*,'(a, a,a,a)') ' sta, cmp, net ', &
- trim(kstnm), trim(kcmpnm), trim(knetwk)
+! if (DEBUG) write(*,'(a, a,a,a)') ' sta, cmp, net ', &
+! trim(kstnm), trim(kcmpnm), trim(knetwk)
! calculate distances and azimuths (needs testing)
call distaz(evla,evlo,stla,stlo,azimuth,backazimuth,dist_deg,dist_km)
@@ -76,12 +76,12 @@
character*8, dimension(:),intent(in) :: kcmpnm
real, dimension(:), intent(in) :: azimuth, dist_km
- real*8,intent(out) :: data_weights(NWINMAX)
+ real,intent(out) :: data_weights(NWINMAX)
real :: daz
integer :: naz(NREGIONS), nwint, i, j, k
character(len=8) :: comp_name
- real*8 :: cmp_weight(NWINMAX), dist_exp_weight(NWINMAX),max_data_weight
+ real :: cmp_weight(NWINMAX), dist_exp_weight(NWINMAX),max_data_weight
daz = 360./NREGIONS
@@ -147,17 +147,17 @@
subroutine compute_A_b(syn_file,data_file,data_weight,tstart,tend,A1,b1,npar)
character(len=*),intent(in) :: syn_file, data_file
- real*8,intent(in) :: data_weight, tstart, tend
+ real,intent(in) :: data_weight, tstart, tend
real*8, intent(out),dimension(:,:) :: A1
real*8, intent(out),dimension(:) :: b1
integer, intent(in) :: npar
real :: t0, dt, t0_1, dt1
integer :: npts, npts1, npts2, nerr,istart,iend, nshift
- integer :: i, j, istart_d, iend_d
+ integer :: i, j, istart_d, iend_d, istart_s, iend_s
real, dimension(NDATAMAX,NPARMAX) :: dsyn_sngl
character(len=150) :: dsyn_file
- real*8 :: dtt, tshift, dlna, cc
+ real :: dlna, cc
! read in data, syn
@@ -171,17 +171,15 @@
if (istart >= iend) stop 'Check tstart and tend'
if (station_correction) then
- data(1:npts) = dble(data_sngl(1:npts))
- syn(1:npts) = dble(syn_sngl(1:npts))
- dtt = dble(dt)
! matching syn(is:ie) with data(is+it:ie+it)
- call calc_criteria(data,syn,npts,istart,iend,dtt,tshift,cc,dlna)
- nshift = nint(tshift/dt)
- istart_d = istart + nshift
- iend_d = iend + nshift
+ call calc_criteria(data_sngl,syn_sngl,npts,istart,iend,nshift,cc,dlna)
+ istart_d = max(1,istart + nshift)
+ iend_d = min(npts,iend + nshift)
+ istart_s=istart_d-nshift
+ iend_s=iend_d-nshift
else
- istart_d = istart
- iend_d = iend
+ istart_d = istart; iend_d = iend
+ istart_s = istart; iend_s = iend
endif
! read in dsyns
@@ -200,13 +198,13 @@
dsyn_sngl(1:npts,i) = (dsyn_sngl(1:npts,i) - syn_sngl(1:npts)) / dcmt_par(i)
endif
enddo
-
+
! compute A and b taking into account data_weights
do j = 1, npar ! col
do i = 1, j ! row
- A1(i,j) = data_weight * sum(dsyn_sngl(istart:iend,i)*dsyn_sngl(istart:iend,j))*dble(dt)
+ A1(i,j) = data_weight * sum(dsyn_sngl(istart_s:iend_s,i)*dsyn_sngl(istart_s:iend_s,j))*dble(dt)
enddo
- b1(j) = data_weight * sum((data_sngl(istart_d:iend_d)-syn_sngl(istart:iend))*dsyn_sngl(istart:iend,j))*dble(dt)
+ b1(j) = data_weight * sum((data_sngl(istart_d:iend_d)-syn_sngl(istart_s:iend_s))*dsyn_sngl(istart_s:iend_s,j))*dble(dt)
enddo
@@ -281,8 +279,8 @@
character(len=*),intent(in):: data_file, syn_file
integer,intent(out) :: npts
- real*8 ,intent(out) :: b, dt
- real*8,intent(in) :: dm(:)
+ real, intent(out) :: b, dt
+ real*8, intent(in) :: dm(:)
real, dimension(NDATAMAX,NPARMAX) :: dsyn_sngl
real, dimension(NDATAMAX) :: time
@@ -292,20 +290,14 @@
character*8 :: kstnm,knetwk,kcmpnm
! read in data, syn
- call rsac1(data_file,data_sngl,npts1,b1,dt1,NDATAMAX,nerr)
- call rsac1(syn_file,syn_sngl,npts2,b1,dt1,NDATAMAX,nerr)
+ call rsac1(data_file,data_sngl,npts1,b,dt,NDATAMAX,nerr)
+ call rsac1(syn_file,syn_sngl,npts2,b,dt,NDATAMAX,nerr)
npts=min(npts1,npts2)
call getkhv('kstnm', kstnm, nerr)
call getkhv('kcmpnm', kcmpnm, nerr)
call getkhv('knetwk', knetwk, nerr)
- ! prepare b, dt, data(:), syn(:) for calc_criteria
- b=dble(b1)
- dt=dble(dt1)
- data(1:npts) = dble(data_sngl(1:npts))
- syn(1:npts) = dble(syn_sngl(1:npts))
-
! read in dsyns
do i = 1, npar
dsyn_file = trim(syn_file) // '.' // trim(PAR_NAME(i))
@@ -321,11 +313,10 @@
enddo
! update synthetics with linearized relation
- new_syn(1:npts) = syn(1:npts) + matmul(dsyn_sngl(1:npts,1:npar),dm(1:npar))
+ new_syn_sngl(1:npts) = syn_sngl(1:npts) + matmul(dsyn_sngl(1:npts,1:npar),dm(1:npar))
! output new synthetics as sac file
if (write_new_syn) then
- syn_sngl(1:npts) = sngl(new_syn(1:npts))
do i = 1, npts
time(i) = b1 + (i-1)*dt1
enddo
@@ -333,13 +324,13 @@
call setkhv('kstnm',trim(kstnm),nerr)
call setkhv('knetwk',trim(knetwk),nerr)
call setkhv('kcmpnm',trim(kcmpnm),nerr)
- call setfhv('b',sngl(b),nerr)
- call setfhv('delta',sngl(dt),nerr)
+ call setfhv('b',b,nerr)
+ call setfhv('delta',dt,nerr)
call setnhv('npts',npts,nerr)
call setihv('iftype','ITIME',nerr)
call setihv('iztype','IB',nerr)
call setlhv('leven',1,nerr)
- call wsac0(trim(syn_file)//'.new',time(1:npts),syn_sngl(1:npts),nerr)
+ call wsac0(trim(syn_file)//'.new',time(1:npts),new_syn_sngl(1:npts),nerr)
if (nerr .ne. 0) stop 'Error writing new synthetics files'
endif
@@ -347,28 +338,46 @@
! ===========================================================
- subroutine calc_criteria(d,s,npts,i1,i2,dt,tshift,cc_max,dlnA)
+ subroutine calc_criteria(d,s,npts,i1,i2,ishift,cc_max,dlnA)
- double precision, dimension(*), intent(in) :: d, s
+ real, dimension(*), intent(in) :: d, s
integer, intent(in) :: npts,i1,i2
- double precision, intent(in) :: dt
- double precision, intent(out) :: tshift,cc_max,dlnA
+ integer, intent(out) :: ishift
+ real, intent(out) :: cc_max,dlnA
- double precision, dimension(NDATAMAX) :: d_win, s_win
- integer :: ishift
+ real, dimension(NDATAMAX) :: d_win, s_win
+ real :: cc,norm
+ integer :: i,i_left,i_right,id_left,id_right,j,nlen
! do cross-correlation
! CHT: zero the data and synthetics outside the window (see comments in xcorr_calc)
d_win(:) = 0. ; d_win(i1:i2) = d(i1:i2)
s_win(:) = 0. ; s_win(i1:i2) = s(i1:i2)
- call xcorr_calc(d_win,s_win,npts,i1,i2,ishift,cc_max)
+ ishift = 0; cc_max = SMALL
- ! cross-correlation time-shift
- tshift = ishift*dt
+ ! length of window (number of points, including ends)
+ nlen = i2 - i1 + 1
- ! calculate dlnA : definition of Dahlen and Baig (2002), Eq. 3,17,18 : dlnA = Aobs/Asyn - 1
- ! NOTE 1: these records are unshifted
- ! NOTE 2: this measurement will reflect any noise in the data, too
+ i_left = -1*int(nlen/2.0)
+ i_right = int(nlen/2.0)
+
+ do i = i_left, i_right
+ id_left = max(1,i1+i) ! left-most point on the data that will be treated
+ id_right = min(npts,i2+i) ! right-most point on the data that will be treated
+ norm = sqrt(sum(s_win(id_left-i:id_right-i)**2) * sum(d_win(id_left:id_right)**2))
+
+ cc=sum(s_win(id_left-i:id_right-i)*d_win(id_left:id_right))
+
+ cc = cc/norm
+
+ if (cc .gt. cc_max) then
+ cc_max = cc
+ ishift = i
+ endif
+ enddo
+
+ ! calculate dlnA : Dahlen and Baig (2002), 3,17,18 -- dlnA = Aobs/Asyn - 1
+
dlnA = sqrt( ( sum( d_win(i1:i2)*d_win(i1:i2) )) / (sum( s_win(i1:i2)*s_win(i1:i2) )) ) - 1.0
end subroutine calc_criteria
Modified: seismo/3D/GRD_CMT3D/cmt3d/cmt3d_sub3.f90
===================================================================
--- seismo/3D/GRD_CMT3D/cmt3d/cmt3d_sub3.f90 2008-11-25 22:57:08 UTC (rev 13404)
+++ seismo/3D/GRD_CMT3D/cmt3d/cmt3d_sub3.f90 2008-11-25 23:24:33 UTC (rev 13405)
@@ -152,66 +152,6 @@
end subroutine gaussian_elimination
-!==========================================================================
+!===========================================================
- subroutine xcorr_calc(d,s,npts,i1,i2,ishift,cc_max)
-
- ! borrowed from Alessia Maggi's flexwin
- ! inputs:
- ! s(npts) = synthetic
- ! d(npts) = data (or observed)
- ! i1, i2 = start and stop indexes of window within s and d
-
- double precision, dimension(*), intent(in) :: s,d
- integer, intent(in) :: npts, i1, i2
-
- ! outputs:
- ! ishift = index lag (d-s) for max cross correlation
- ! cc_max = maximum of cross correlation (normalised by sqrt(synthetic*data))
- integer, intent(out) :: ishift
- double precision, intent(out) :: cc_max
-
- ! local variables
- integer :: nlen
- integer :: i_left, i_right, i, j, id_left, id_right
- double precision :: cc, norm
-
- ! initialise shift and cross correlation to zero
- ishift = 0
- cc_max = 0.0d0
-
- if (i1.lt.1 .or. i1.gt.i2 .or. i2.gt.npts) then
- write(*,*) 'Error with window limits: i1, i2, npts ', i1, i2, npts
- return
- endif
-
- ! length of window (number of points, including ends)
- nlen = i2 - i1 + 1
-
- ! left and right limits of index (time) shift search
- ! NOTE: This looks OUTSIDE the time window of interest to compute TSHIFT and CC.
- ! THIS IS A VERY IMPORTANT DECISION!
- i_left = -1*int(nlen/2.0)
- i_right = int(nlen/2.0)
-
- ! i -> shift (to be applied to DATA (d) in cc search)
- do i = i_left, i_right
- cc = 0.
- id_left = max(1,i1+i) ! left-most point on the data that will be treated
- id_right = min(npts,i2+i) ! right-most point on the data that will be treated
- norm = sqrt(sum(s(i1:i2)*s(i1:i2)) * sum(d(id_left:id_right)*(d(id_left:id_right))))
- do j = i1, i2
- if((j+i).ge.1 .and. (j+i).le.npts) cc = cc + s(j)*d(j+i)
- enddo
- cc = cc/norm
- if (cc .gt. cc_max) then
- cc_max = cc
- ishift = i
- endif
- enddo
-
-end subroutine xcorr_calc
-
-!------------------------------------------------------------------------
-
end module cmt3d_sub3
Modified: seismo/3D/GRD_CMT3D/cmt3d/cmt3d_variables.f90
===================================================================
--- seismo/3D/GRD_CMT3D/cmt3d/cmt3d_variables.f90 2008-11-25 22:57:08 UTC (rev 13404)
+++ seismo/3D/GRD_CMT3D/cmt3d/cmt3d_variables.f90 2008-11-25 23:24:33 UTC (rev 13405)
@@ -9,12 +9,12 @@
! inversion parameters
character(len=150) :: cmt_file, new_cmt_file
integer :: npar
- real*8 :: ddelta, ddepth, dmoment
+ real :: ddelta, ddepth, dmoment
! data selection
character(len=150) :: flexwin_out_file
logical :: weigh_data_files
- real*8 :: comp_z_weight, comp_r_weight, comp_t_weight,&
+ real :: comp_z_weight, comp_r_weight, comp_t_weight,&
az_exp_weight, &
pnl_dist_weight, rayleigh_dist_weight, love_dist_weight
@@ -32,11 +32,10 @@
! number of files and windows
integer :: nfiles,nwins(NRECMAX),nwin_total
- real*8 :: data_weights(NWINMAX)
+ real :: data_weights(NWINMAX)
! data and syn arrays
- real*8, dimension(NDATAMAX) :: data, syn, new_syn
- real, dimension(NDATAMAX) :: data_sngl, syn_sngl
+ real, dimension(NDATAMAX) :: data_sngl, syn_sngl, new_syn_sngl
character(len=150) :: data_file, syn_file
Modified: seismo/3D/GRD_CMT3D/grid3d/GRID3D.PAR
===================================================================
--- seismo/3D/GRD_CMT3D/grid3d/GRID3D.PAR 2008-11-25 22:57:08 UTC (rev 13404)
+++ seismo/3D/GRD_CMT3D/grid3d/GRID3D.PAR 2008-11-25 23:24:33 UTC (rev 13405)
@@ -3,7 +3,7 @@
flexwin.out
.true. -- weigh_data_files
2 2 1 0 1.15 0.55 0.78 -- weights of data comp, az(exp), dist(exp)
-.true. -- station_correction
+.true. 6.0 -- station_correction
.false. 3 -- global_search, ncalc (=1 if local search)
0 180 10
0 90 10
Modified: seismo/3D/GRD_CMT3D/grid3d/grid3d_sub.f90
===================================================================
--- seismo/3D/GRD_CMT3D/grid3d/grid3d_sub.f90 2008-11-25 22:57:08 UTC (rev 13404)
+++ seismo/3D/GRD_CMT3D/grid3d/grid3d_sub.f90 2008-11-25 23:24:33 UTC (rev 13405)
@@ -27,7 +27,7 @@
az_exp_weight, &
pnl_dist_weight, rayleigh_dist_weight, love_dist_weight
- read(IOPAR,*) station_correction
+ read(IOPAR,*) station_correction, tshift_max
read(IOPAR,*) global_search, ncalc
if (.not. global_search) then
ncalc = 1
@@ -55,7 +55,7 @@
! these are global search control parameters -- can be adjusted
! 19x10x13x5
s_strike = 0; e_strike = 180; d_strike = 10
- s_dip = 0; e_dip = 90; d_dip = 10
+ s_dip = 0; e_dip = 90; d_dip = 10
s_rake = -180; e_rake = 180; d_rake = 30
s_mw = mw * 0.9; e_mw = mw * 1.1; d_mw = mw * 0.05
t_strike=1.5; t_dip=1.5; t_rake=1.5; t_mw=1
@@ -87,6 +87,11 @@
endif
if (station_correction) then
write(*,*) 'shift data to aligh with synthetics before grid search'
+ if (tshift_max > 0) then
+ write(*,'(a,g15.5)') ' with maximum shift allowed', tshift_max
+ else
+ stop 'tshift_max should be > 0'
+ endif
else
write(*,*) 'data are NOT shifted before grid search'
endif
Modified: seismo/3D/GRD_CMT3D/grid3d/grid3d_sub2.f90
===================================================================
--- seismo/3D/GRD_CMT3D/grid3d/grid3d_sub2.f90 2008-11-25 22:57:08 UTC (rev 13404)
+++ seismo/3D/GRD_CMT3D/grid3d/grid3d_sub2.f90 2008-11-25 23:24:33 UTC (rev 13405)
@@ -189,7 +189,7 @@
real,intent(inout) :: misfit(:)
integer :: npts,npts2,nerr,is(nw),ie(nw)
- integer :: iw,it,id,ir,im,i,j,ishift,iss,iee
+ integer :: iw,it,id,ir,im,i,j,ishift,iss,iee,isd,ied,ishift_max
real :: b,dt,b2,dt2,cc,t1(nw),tt
character(len=150) :: dsyn_file
@@ -211,13 +211,19 @@
enddo
! core computation
- ishift=0
+ ishift=0; ishift_max=nint(tshift_max/dt)
do j = 1, n_total
syn(1:npts) = matmul(mij(:,j),dsyn(:,1:npts)) / dmoment
do iw = 1, nw
tt=t1(iw); iss=is(iw);iee=ie(iw)
- if (station_correction) call xcorr_calc(data,syn,npts,iss,iee,ishift,cc)
- misfit(j) = misfit(j)+tt*sum((syn(iss:iee)-data(iss+ishift:iee+ishift))**2)
+ if (station_correction) then
+ call xcorr_calc(data,syn,npts,iss,iee,ishift,cc,ishift_max)
+ isd=max(1,iss+ishift); ied=min(npts,iee+ishift)
+ iss=isd-ishift; iee=ied-ishift
+ else
+ isd=iss; ied=iee
+ endif
+ misfit(j) = misfit(j)+tt*sum((syn(iss:iee)-data(isd:ied))**2)
enddo
enddo
Modified: seismo/3D/GRD_CMT3D/grid3d/grid3d_sub3.f90
===================================================================
--- seismo/3D/GRD_CMT3D/grid3d/grid3d_sub3.f90 2008-11-25 22:57:08 UTC (rev 13404)
+++ seismo/3D/GRD_CMT3D/grid3d/grid3d_sub3.f90 2008-11-25 23:24:33 UTC (rev 13405)
@@ -7,7 +7,7 @@
!==========================================================================
- subroutine xcorr_calc(dd,ss,npts,i1,i2,ishift,cc_max)
+ subroutine xcorr_calc(dd,ss,npts,i1,i2,ishift,cc_max,ishift_max)
! the core of Alessia's version
@@ -15,6 +15,7 @@
integer, intent(in) :: npts, i1, i2
integer, intent(out) :: ishift
real,intent(out) :: cc_max
+ integer,intent(in) :: ishift_max
! local variables
real,dimension(npts) :: d,s
@@ -27,7 +28,7 @@
ishift = 0; cc_max = SMALL
! length of window (number of points, including ends)
- nlen = i2 - i1 + 1
+ nlen = min(i2 - i1 + 1,ishift_max * 2)
i_left = -1*int(nlen/2.0)
i_right = int(nlen/2.0)
Modified: seismo/3D/GRD_CMT3D/grid3d/grid3d_variables.f90
===================================================================
--- seismo/3D/GRD_CMT3D/grid3d/grid3d_variables.f90 2008-11-25 22:57:08 UTC (rev 13404)
+++ seismo/3D/GRD_CMT3D/grid3d/grid3d_variables.f90 2008-11-25 23:24:33 UTC (rev 13405)
@@ -18,7 +18,7 @@
! grid search scheme
logical :: station_correction,global_search
integer :: ncalc
- real :: s_strike,e_strike,d_strike,s_dip,e_dip,d_dip, &
+ real :: tshift_max,s_strike,e_strike,d_strike,s_dip,e_dip,d_dip, &
s_rake,e_rake,d_rake,s_mw,e_mw,d_mw
! misc
Added: seismo/3D/GRD_CMT3D/scripts/mech_table.pl
===================================================================
--- seismo/3D/GRD_CMT3D/scripts/mech_table.pl (rev 0)
+++ seismo/3D/GRD_CMT3D/scripts/mech_table.pl 2008-11-25 23:24:33 UTC (rev 13405)
@@ -0,0 +1,100 @@
+#!/usr/bin/perl
+use lib '/opt/seismo/lib/perl';
+use GMT_PLOT;
+use GMT_PLOT_SC;
+use CMT_TOOLS;
+use POSIX;
+
+
+# parameter files
+$inv_file="INVERSION.PAR.SAVE";
+$outdir="mech"; $psfile="mech.ps";
+
+$weigh_data="2 2 1 0.5 1.15 0.55 0.78";
+$weigh_data_nocomp="1 1 1 0.5 1.15 0.55 0.78";
+$weigh_data_noaz="2 2 1 0 1.15 0.55 0.78";
+
+ at names=("initial","grid_search","7 Par+ZT", "6 Par+ZT","7 Par+no-SC","7 Par+no-AZ","6 Par+dc","7 Par+DC","7 Par+no-Comp");
+ at npars=(6,6, 7, 6,7,7, 6,7,7);
+ at wdatas=($weigh_data, $weigh_data, $weigh_data, $weigh_data,$weigh_data,$weigh_data_noaz, $weigh_data,$weigh_data,$weigh_data_nocomp);
+ at scorrs=("true","true",".true.", ".true.", ".false.", ".true.", ".true.",".true.",".true.");
+ at cons=(".true. true. 0.0", ".true. .true. 0.0", ".true. .false. 0.0",
+ ".true. .false. 0.0", "true. .false. 0.0", "true. .false. 0.0",
+ ".true. .true. 0.0",".true. .true. 0.0",".true. .false. 0.0");
+
+$nbeach=@npars;
+$ncols = 3;
+$nrows = 3;
+# plot size jx/jy, plot origin in xy
+($jx,$jy) = shift_xy("1 1","$ncols $nrows","7 8", \@xy,"0.8 0.8");
+$JX = "-JX$jx/$jy"; $R="-R0/3/0/3";
+$GMT_PLOT::paper_orient="-P";
+
+if (not -f "cmt3d_flexwin") {die("Check if cmt3d_flexwin exists or not\n");}
+
+open(INV,"$inv_file");@all=<INV>;close(INV);
+$cmt=$all[0]; $cmt_new=$all[1]; $npar=$all[2]; $delta=$all[3];
+$flex=$all[4];$lwdata=$all[5]; $wdata=$all[6]; $scorr=$all[7]; $con=$all[8];
+$wns=$all[9];
+chomp($cmt_new);chomp($cmt);
+(undef,undef,$ename)=split(" ",`grep 'event name' $cmt`);
+
+ at cmt=($cmt);
+if (-f "${cmt}_GRD") {@cmt=(@cmt,"${cmt}_GRD"); print BASH "cp -f ${cmt}_GRD $outdir\n";}
+for ($i=2;$i<$nbeach;$i++) {@cmt=(@cmt,"CMTSOLUTION.$i");}
+
+open(BASH,">mech_table.bash");
+print BASH "mkdir -p $outdir\n";
+print BASH "cp -f $cmt $outdir\n";
+if (-f "${cmt}_GRD") {print BASH "cp -f ${cmt}_GRD $outdir\n";}
+
+for ($k=2;$k<$nbeach;$k++) { # cmt3d runs
+ print BASH "echo Running cmt3d_flexwin for INVERSION.$k ...\n";
+ open(INV,">INVERSION.$k");
+ print INV "${cmt}\n${cmt_new}\n$npars[$k]\n$delta${flex}.true.\n$wdatas[$k]\n$scorrs[$k]\n$cons[$k]\n${wns}";
+ close(INV);
+ print BASH "# --- mech $k ----\n";
+ print BASH "cp -f INVERSION.$k INVERSION.PAR\n";
+ print BASH "cmt3d_flexwin > cmt3d.stdout\n";
+ print BASH "if [ \$? != 0 ]; then\n echo \"exiting case $k\";exit\nfi\n";
+ print BASH "mv -f $cmt_new $outdir/$cmt[$k]\n";
+ print BASH "mv -f cmt3d_flexwin.out $outdir/cmt3d_out.$k\n";
+ print BASH "mv -f cmt3d.stdout $outdir/cmt3d_stdout.$k\n";
+ print BASH "mv -f INVERSION.$k $outdir\n";
+}
+close(BASH);
+system("chmod a+x mech_table.bash; mech_table.bash");
+
+for ($k=2;$k<$nbeach;$k++) {
+ (@tmp)=split(" ",`grep Variance $outdir/cmt3d_stdout.$k`);
+ $var[$k]=$tmp[8]; print "$k,$var[$k]\n";}
+#die("here\n");
+
+
+open(BASH,">mech_plot.bash");
+plot_psxy(\*BASH,$psfile,"$JX -K -X0 -Y0 $R","");
+$k=0;
+for ($i=0;$i<$nrows;$i++) {
+ for ($j=0;$j<$ncols;$j++) {
+ print BASH "# ---- Mech $k ------\n";
+ $xy=$xy[$i][$j];($x,$y) = split(/\//,$xy);
+ plot_psxy(\*BASH,$psfile,"$JX -X$x -Y$y $R","");
+ if (-f "$outdir/$cmt[$k]") {
+ plot_psmeca_raw(\*BASH,$psfile,"-JX -R -Sm1.0 -B3/3wesn","1.5 2","$outdir/$cmt[$k]");
+ @output = `cmtsol2faultpar.pl $outdir/$cmt[$k]`;
+ (undef,$tmp1) = split(" ", $output[4]);
+ (undef,$tmp2) = split(" ", $output[6]);
+ ($s,$d,$r,$p,$t) = split(/\//,$tmp2);
+ $tex="1 1 9 0 4 CM $names[$k]\n";
+ $tex.="0.2 0.8 9 0 4 LM Mw/dep/eps/= $tmp1\n0.2 0.6 9 0 4 LM S/D/R= $s/$d/$r";
+ if ($k>=2) {$tex.=sprintf("\n0.2 0.4 9 0 4 LM Var=%6.2f%",$var[$k]);}
+ plot_pstext(\*BASH,$psfile,"-JX -R -B -G0/0/255","$tex");
+ }
+ if ($k==1) {plot_pstext(\*BASH,$psfile,"-JX -R -B -N","1.5 4 12 0 4 CM $ename");}
+ plot_psxy(\*BASH,$psfile,"-JX -X-$x -Y-$y","");
+ $k++;
+ }
+}
+plot_psxy(\*BASH,$psfile,"-JX -O","");
+close(BASH);
+system("chmod a+x mech_plot.bash; mech_plot.bash");
Property changes on: seismo/3D/GRD_CMT3D/scripts/mech_table.pl
___________________________________________________________________
Name: svn:executable
+ *
Modified: seismo/3D/GRD_CMT3D/scripts/plot_grid.pl
===================================================================
--- seismo/3D/GRD_CMT3D/scripts/plot_grid.pl 2008-11-25 22:57:08 UTC (rev 13404)
+++ seismo/3D/GRD_CMT3D/scripts/plot_grid.pl 2008-11-25 23:24:33 UTC (rev 13405)
@@ -7,8 +7,13 @@
use POSIX;
# specify the grid search parameter range
-if (@ARGV == 0) {die(" plot_grid.pl grid-output-files\n");}
+if (@ARGV == 0 or @ARGV < 2) {die(" plot_grid.pl cmt_file grid-output-files \n");}
+$cmt=$ARGV[0];
+if (not -f $cmt) {die("Check $cmt file\n");}
+(undef,undef,$ename)=split(" ",`grep 'event name' $cmt`);
+print "CMT file -- $ename\n";
+
# this controls the paper orientation;
$GMT_PLOT::paper_orient = "-P";
@@ -22,13 +27,13 @@
$Bs[-1]="WeSn"; $Bd[-1]="WeSn"; $Br[-1]="WeSn";
# plot size jx/jy, plot origin in xy
-($jx,$jy) = shift_xy("1 1","$ncols $nrows","7 8", \@xy,"0.84 0.84");
+($jx,$jy) = shift_xy("1 1.5","$ncols $nrows","7 8", \@xy,"0.84 0.84");
$JX = "-JX$jx/$jy";
$GMT_PLOT::paper_orient="-P";
# loop over input files
-for ($k=0;$k<@ARGV;$k++) {
+for ($k=1;$k<@ARGV;$k++) {
$result=$ARGV[$k];
print " == file = $result \n";
@@ -40,6 +45,9 @@
$psfile = "grid_$k.ps";
($ss,$es,$ns,$ts,$sd,$ed,$nd,$td,$sr,$er,$nr,$tr,$min0,$max0, at array) = pick_ds_dd_dr($result,$nrows);
+ # LQY -- hard-wire this maximum color control
+ if ($max0 > 100*$min0) {$max0=100*$min0;}
+
$j=0;
@as=@array[$j..$j+$ns-1]; $j+=$ns;
@ad=@array[$j..$j+$nd-1]; $j+=$nd;
@@ -49,7 +57,7 @@
@pr=@array[$j..$j+$nrows-1]; $j+=$nrows;
print "s,d,r=$ss,$es,$ns,$ts, $sd,$ed,$nd,$td, $sr,$er,$nr,$tr\n";
- print "min/max=$min0,$max0\n";
+ print "min/max=$min0/$max0\n";
print "as=@as\nad=@ad\nar=@ar\n";
print "ps=@ps\npd=@pd\npr=@pr\n";
@@ -72,10 +80,10 @@
# plot position of the text for dip-rake, strike-rake, and strike-dip plot
@pos= ("$ad[1] $ar[-2]","$as[1] $ar[-2]","$as[1] $ad[-2]");
- $arr=$ar[-1]+$ddr*1.1;
+ $arr=$ar[-1]+$ddr*2; # position for title
$pp="$as[$nrows/2] $arr";
- print CSH "makecpt -T$min/$max/$db2 -Crainbow -Z -I> temp.cpt \n";
+ print CSH "makecpt -T$min/$max/$db2 -Cseis -Z -I > temp.cpt \n";
plot_psxy(\*CSH,$psfile,"$JX -K -X0 -Y0","");
@@ -100,12 +108,12 @@
print CSH "grdimage out.grd $R[$i] -JX $B -Ctemp.cpt -K -O -P >> $psfile \n";
print CSH "grdcontour out.grd -R -JX -K -O -P $A>> $psfile\n";
plot_pstext(\*CSH,$psfile,"-JX -W255 ","$pos[$i] 8 0 4 LT $name[$i]=$value");
- if ($i==1 and $j==0) {plot_pstext(\*CSH,$psfile,"-JX -W255 -N","$pp 9 0 4 LT title");}
+ if ($i==1 and $j==0) {plot_pstext(\*CSH,$psfile,"-JX -W255 -N","$pp 14 0 4 LT $ename");}
plot_psxy(\*CSH,$psfile,"-JX -X-$x -Y-$y","");
print CSH "\n";
}
}
-print CSH "psscale -Ctemp.cpt -D4i/0.7i/2.5i/0.17h $db -K -O -V -P >> $psfile\n";
+print CSH "psscale -Ctemp.cpt -D4i/1.0i/2.5i/0.17h $db -K -O -V -P >> $psfile\n";
plot_psxy(\*CSH,$psfile,"$JX -O","");
close(CSH);
@@ -124,6 +132,7 @@
close(FILE);
%nns=(); %nnd=(); %nnr=(); @ps=(); @pd=(); @pr=();
($ss,$es,$sd,$ed,$sr,$er,$min,$max) = split(" ",`minmax -C $file`);
+
for $line (@lines) {
($s,$d,$r,$v) = split(" ",$line);
if ($v == $min) {$mins=$s;$mind=$d;$minr=$r;}
@@ -136,6 +145,8 @@
@ad=sort { $a <=> $b } keys %nnd;
@ar=sort { $a <=> $b } keys %nnr;
$ts=$as[1]-$as[0]; $td=$ad[1]-$ad[0]; $tr=$ar[1]-$ar[0];
+ if ($nrows > @as or $nrows > @ad or $nrows > @ar) {
+ die("Choose a smaller nrows to be in accordance with ns,nd,nr\n");}
$nleft=floor($nrows/2);
$nright=$nrows-$nleft;
@@ -143,12 +154,16 @@
$is=int(($mins-$ss)/$ts); $id=int(($mind-$sd)/$td); $ir=int(($minr-$sr)/$tr);
for ($i=-$nleft;$i<$nright;$i++) {
- @ps=(@ps,$as[$is+$i*$ns1]);
- @pd=(@pd,$ad[$id+$i*$nd1]);
- @pr=(@pr,$ar[$ir+$i*$nr1]);
+ if ($is+$i*$ns1 >= @as) {$iss=$is+$i*$ns1- at as;} else {$iss=$is+$i*$ns1;}
+ if ($id+$i*$nd1 >= @ad) {$idd=$id+$i*$nd1- at ad;} else {$idd=$id+$i*$nd1;}
+ if ($ir+$i*$nr1 >= @ar) {$irr=$ir+$i*$nr1- at ar;} else {$irr=$ir+$i*$nr1;}
+
+ @ps=(@ps,$as[$iss]);
+ @pd=(@pd,$ad[$idd]);
+ @pr=(@pr,$ar[$irr]);
}
-# print "ps=@ps\npd=@pd\npr=@pr\n";
return ($ss,$es,$ns,$ts,$sd,$ed,$nd,$td,$sr,$er,$nr,$tr,$min,$max, at as, at ad, at ar, at ps, at pd, at pr);
+
}
Modified: seismo/3D/GRD_CMT3D/scripts/process_data_and_syn.pl
===================================================================
--- seismo/3D/GRD_CMT3D/scripts/process_data_and_syn.pl 2008-11-25 22:57:08 UTC (rev 13404)
+++ seismo/3D/GRD_CMT3D/scripts/process_data_and_syn.pl 2008-11-25 23:24:33 UTC (rev 13405)
@@ -52,8 +52,10 @@
if ($opt_a) {$sta_file=$opt_a;} else {die("-a option\n");}
if (not -f $cmt_file or not -f $sta_file) {die("Check $cmt_file or $sta_file\n");}
if ($opt_n) {$nsps="-s $opt_n";} else {$nsps = "-s 20";}
- (undef,undef, $hdur) = split(" ",`grep half $cmt_file`);
- if (abs($hdur)<0.1) {die("Check half duration for -h option in syn processing\n");}
+
+# uncomment this if you choose to convolve hdur here
+# (undef,undef, $hdur) = split(" ",`grep half $cmt_file`);
+# if (abs($hdur)<0.1) {die("Check half duration for -h option in syn processing\n");}
}
if ($opt_d) {$dirdat_pp= "${data_dir}_PP";}
More information about the CIG-COMMITS
mailing list