[cig-commits] r19462 - in seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES: . noise_layered noise_layered/masks noise_layered/model_0 noise_layered/model_1 noise_layered/model_2

xie.zhinan at geodynamics.org xie.zhinan at geodynamics.org
Tue Jan 24 11:11:26 PST 2012


Author: xie.zhinan
Date: 2012-01-24 11:11:25 -0800 (Tue, 24 Jan 2012)
New Revision: 19462

Added:
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/adj_mask.f90
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/changemask
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/changerec
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/clean
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/getdirect
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/masks/
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/masks/drillbit
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/masks/uniform
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_0/
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_0/Par_file_fair
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_0/Par_file_fair.before_update_to_r19xxx
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_0/Par_file_good
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_0/Par_file_good.before_update_to_r19xxx
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_0/SOURCE_noise
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_0/interfaces_fair
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_0/interfaces_good
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/Par_file_best
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/Par_file_best.before_update_to_r19xxx
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/Par_file_fair
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/Par_file_fair.before_update_to_r19xxx
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/Par_file_good
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/Par_file_good.before_update_to_r19xxx
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/SOURCE_noise
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/interfaces_best
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/interfaces_fair
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/interfaces_good
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_2/
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_2/Par_file_fair
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_2/Par_file_fair.before_update_to_r19xxx
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_2/Par_file_good
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_2/Par_file_good.before_update_to_r19xxx
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_2/SOURCE_noise
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_2/interfaces_fair
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_2/interfaces_good
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/replicate
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/run_adj
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/run_all
   seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/run_fwd
Log:
add noise_layered example



Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/adj_mask.f90
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/adj_mask.f90	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/adj_mask.f90	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,440 @@
+program adj_cc
+
+implicit none
+
+!choose whether to differentiate traces 0, 1, or 2 times
+integer, parameter :: differentiate = 1
+
+!choose whether to bandpass filter
+logical, parameter :: use_filtering = .true.
+
+!choose exactly one of the following window options
+logical, parameter :: use_negative_branch = .false.
+logical, parameter :: use_positive_branch = .true.
+logical, parameter :: use_custom_window = .false.
+
+!choose whether to time reverse, carried out subsequent to all other processing
+logical, parameter :: time_reverse = .false.
+
+
+
+! FILTERING PARAMETERS
+real  freq_low,freq_high
+data  freq_low  / 1.d-2 /
+data  freq_high / 1.5d1 /
+
+! CUSTOM WINDOW PARAMETERS
+real :: t_begin, t_end, w_tukey
+data t_begin / 45.d0 /
+data t_end   / 65.d0 /
+data w_tukey / 0.4   /
+!see explanation below
+
+! time variables
+integer :: it, nt, nthalf
+double precision :: dt, off 
+
+! data variables
+double precision, dimension(:), allocatable :: seismo_1, seismo_2, seismo_3, seismo_4, &
+ seismo_adj, t, w
+
+! input/ output
+character(len=512) arg1, arg2
+character(len=512) filename
+integer :: ios
+
+! miscellaneous
+double precision, parameter :: PI = 3.141592653589793
+integer :: it_off, it_wdt, it_begin, it_end, k
+integer :: ifreq, nfreq
+real :: F1,F2,D(8),G,DELT
+real :: alpha, beta
+
+
+! EXPLANATION OF WINDOW PARAMETERS
+
+!To select the desired branch of the cross-correlogram, we employ a Tukey window.  A Tukey taper is just a variant of a cosine taper.
+
+!W_TUKEY is a number between 0 and 1, 0 being pure boxcar and 1 being pure cosine
+
+
+
+
+! Get arguments
+call getarg(1,arg1)
+call getarg(2,arg2)
+filename = arg1
+read(arg2,*) off
+
+! Get file info
+call getlen(filename,nt)
+call getinc(filename,nt,dt)
+nthalf = (nt+1)/2
+
+! Used to mask direct arrivals
+off = abs(off) + 0.30
+
+
+write(*,*) ''
+write(*,*) 'This routine works only for evenly sampled time series.'
+write(*,*) 'Reading from file: '//trim(filename)
+write(*,'(a,i10)')   ' nt: ', nt
+write(*,'(a,f10.3)') ' dt: ', dt
+
+! Allocate, initialize
+allocate(t(nt))
+allocate(w(nt))
+allocate(seismo_1(nt))
+allocate(seismo_2(nt))
+allocate(seismo_3(nt))
+allocate(seismo_4(nt))
+allocate(seismo_adj(nt))
+w(:)          = 0.0d0
+seismo_1(:)   = 0.0d0
+seismo_2(:)   = 0.0d0
+seismo_3(:)   = 0.0d0
+seismo_4(:)   = 0.0d0
+seismo_adj(:) = 0.0d0
+
+
+!!!!!!!!!! READ INPUT !!!!!!!!!!!!!!!!!!!!
+open(unit=1001,file=trim(filename),status='old',action='read')
+do it = 1, nt
+    read(1001,*) t(it), seismo_1(it)
+end do
+close(1001)
+
+
+!!!!!!!!!! DIFFERENTIATE !!!!!!!!!!!!!!!!!!!
+seismo_2(1) = 0.0
+seismo_2(nt) = 0.0
+do it = 2, nt-1
+    if (differentiate == 0) seismo_2(it) = seismo_1(it)
+    if (differentiate == 1) seismo_2(it) = ( seismo_1(it+1) - seismo_1(it-1) ) / (2*dt)
+    if (differentiate == 2) seismo_2(it) = ( seismo_1(it+1) - 2 * seismo_1(it) + seismo_1(it-1) ) / dt**2
+
+end do
+
+
+!!!!!!!!!! FILTER !!!!!!!!!!!!!!!!!!!!
+seismo_3 = seismo_2
+if (use_filtering) then
+! THIS SECTION CALCULATES THE FILTER AND MUST BE CALLED BEFORE          
+! FILTER IS CALLED                                                      
+DELT = 1.0d3 * dt
+F1=freq_low
+F2=freq_high
+call BNDPAS(F1,F2,DELT,D,G,nt)
+!    F1 = LOW FREQUENCY CUTOFF (6 DB DOWN)                             
+!    F2 = HIGH FREQUENCY CUTOFF (6 DB DOWN)                             
+!    DELT = SAMPLE INTERVAL IN MILLISECONDS                             
+!    D = WILL CONTAIN 8 Z DOMAIN COEFICIENTS OF RECURSIVE FILTER        
+!    G = WILL CONTAIN THE GAIN OF THE FILTER,
+call FILTER(seismo_3,nt,D,G,2)
+!    X = DATA VECTOR OF LENGTH N CONTAINING DATA TO BE FILTERED        
+!    D = FILTER COEFFICIENTS CALCULATED BY BNDPAS                      
+!    G = FILTER GAIN                                                   
+!    IG = 1  one pass
+!    IG = 2  two passes
+end if
+
+
+!!!!!!!!!! WINDOW !!!!!!!!!!!!!!!!!!!!
+if (use_custom_window) then
+  it_begin = floor((t_begin - t(1))/dt)
+  it_end   = ceiling((t_end - t(1))/dt)
+  if (it_begin < 1) it_begin = 1
+  if (it_end > nt) it_end = nt
+
+elseif (use_positive_branch) then
+  write(*,*) 'Choosing positive branch'
+  it_begin = nthalf + floor(off/dt)
+  it_end   = nt
+  if (it_begin < nthalf) it_begin = nthalf
+  if (it_end > nt) it_end = nt
+
+elseif (use_negative_branch) then
+  write(*,*) 'Choosing negative branch'
+  it_begin = 1
+  it_end   = nthalf - floor(off/dt)
+  if (it_begin < 1) it_begin = 1
+  if (it_end > nthalf) it_end = nthalf
+
+else
+  write(*,*) 'Must select one of the following: positive_branch, &
+              negative_branch, custom_window.'
+
+endif
+
+write(*,'(a,2f10.3)') ' Time range: ', t(1), t(nt)
+write(*,'(a,2f10.3)') ' Window:     ', t(it_begin), t(it_end)
+write(*,'(a,f10.3,f10.3)') ' Filtering:  ', 1./freq_high, 1./freq_low
+
+!! Tukey taper
+alpha = w_tukey
+k=0
+do it = it_begin,it_end
+  k=k+1
+  beta = real(k-1)/(it_end-it_begin)
+
+  if (beta<alpha/2.) then
+    w(it) = 0.5*(1.+cos(2.*pi/alpha*(beta-alpha/2.)))
+
+  elseif (beta>alpha/2. .and. beta<1.-alpha/2.) then
+    w(it) = 1.0
+
+  else
+    w(it) = 0.5*(1.+cos(2*pi/w_tukey*(beta-1.+alpha/2.)))
+
+  endif
+end do
+seismo_4 = w * seismo_3
+
+
+!!!!!!!!!! NORMALIZE !!!!!!!!!!!!!!!!!!!!
+seismo_adj = - seismo_4/(DOT_PRODUCT(seismo_1,seismo_1)*dt)
+!seismo_adj = w
+
+!!!!!!!!!! WRITE ADJOINT SOURCE !!!!!!!!!!!!!!!!!!!!
+open(unit=1002,file=trim(filename)//'.adj',status='unknown',iostat=ios)
+if (ios /= 0) write(*,*) 'Error opening output file.'
+
+write(*,*) ''
+write(*,*) 'Writing to file: '//trim(filename)//'.adj'
+
+do it = 1,nt
+    if (.not. time_reverse) write(1002,'(f16.12,1pe18.8)'), t(it), seismo_adj(it)
+    if (time_reverse) write(1002,'(f16.12,1pe18.8)'), t(it), seismo_adj(nt+1-it)
+end do
+close(1002)
+
+write(*,*) 'Finished writing to file.'
+write(*,*) ''
+
+
+end program adj_cc
+
+
+
+!=====================================================================
+subroutine getlen(filename,len)
+
+implicit none
+
+!input
+character(len=512) :: filename
+
+!output
+integer :: len
+
+!local
+integer, parameter :: IMAX = 1000000
+integer :: i,ios
+real :: dummy1, dummy2
+
+open(unit=1001,file=trim(filename),status='old',action='read')
+len=0
+do i=1,IMAX
+    read(1001,*,iostat=ios) dummy1, dummy2
+    if (ios==-1) exit
+    len=len+1
+enddo
+close(1001)
+
+end subroutine getlen
+
+
+
+!=====================================================================
+subroutine getinc(filename,len,inc)
+
+implicit none
+
+!input
+character(len=64) :: filename
+integer :: len
+
+!output
+double precision :: inc
+
+!local
+integer :: it
+double precision, dimension(len) :: t
+double precision :: sumdt
+real :: dummy
+
+open(unit=1001,file=trim(filename),status='old',action='read')
+do it=1,len
+    read(1001,*) t(it), dummy
+enddo
+close(1001)
+
+sumdt = 0.0d0
+do it=1,len-1
+    sumdt = sumdt + t(it+1) - t(it)
+enddo
+inc=sumdt/(len-1)
+
+end subroutine getinc
+
+
+!=====================================================================
+SUBROUTINE BNDPAS(F1,F2,DELT,D,G,N)                                         
+! RECURSIVE BUTTERWORTH BAND PASS FILTER (KANASEWICH, TIME SERIES       
+! ANALYSIS IN GEOPHYSICS, UNIVERSITY OF ALBERTA PRESS, 1975; SHANKS,    
+! JOHN L, RECURSION FILTERS FOR DIGITAL PROCESSING, GEOPHYSICS, V32,    
+! FILTER.  THE FILTER WILL HAVE 8 POLES IN THE S PLANE AND IS           
+! APPLIED IN FORWARD AND REVERSE DIRECTIONS SO AS TO HAVE ZERO          
+! PHASE SHIFT.  THE GAIN AT THE TWO FREQUENCIES SPECIFIED AS            
+! CUTOFF FREQUENCIES WILL BE -6DB AND THE ROLLOFF WILL BE ABOUT         
+! THE FILTER TO PREVENT ALIASING PROBLEMS.                              
+    COMPLEX P(4),S(8),Z1,Z2                                           
+    real D(8),XC(3),XD(3),XE(3)                            
+    double precision :: X(N) 
+    DATA ISW/0/,TWOPI/6.2831853/                                      
+! THIS SECTION CALCULATES THE FILTER AND MUST BE CALLED BEFORE          
+! FILTER IS CALLED                                                      
+                                                                       
+!    F1 = LOW FREQUENCY CUTOFF (6 DB DOWN)                             
+!    F2 = HIGH FREQUENCY CUTOFF (6 DB DOWN)                             
+!    DELT = SAMPLE INTERVAL IN MILLISECONDS                             
+!    D = WILL CONTAIN 8 Z DOMAIN COEFICIENTS OF RECURSIVE FILTER        
+!    G = WILL CONTAIN THE GAIN OF THE FILTER,                           
+
+      DT=DELT/1000.0                                                    
+      TDT=2.0/DT                                                        
+      FDT=4.0/DT                                                        
+      ISW=1                                                             
+      P(1)=CMPLX(-.3826834,.9238795)                                    
+      P(2)=CMPLX(-.3826834,-.9238795)                                   
+      P(3)=CMPLX(-.9238795,.3826834)                                    
+      P(4)=CMPLX(-.9238795,-.3826834)                                   
+      W1=TWOPI*F1                                                       
+      W2=TWOPI*F2                                                       
+      W1=TDT*TAN(W1/TDT)                                                
+      W2=TDT*TAN(W2/TDT)                                               
+      HWID=(W2-W1)/2.0                                                  
+      WW=W1*W2                                                          
+      DO 19 I=1,4                                                       
+      Z1=P(I)*HWID                                                      
+      Z2=Z1*Z1-WW                                                       
+      Z2=CSQRT(Z2)                                                      
+      S(I)=Z1+Z2                                                        
+   19 S(I+4)=Z1-Z2                                                      
+      G=.5/HWID                                                         
+      G=G*G                                                             
+      G=G*G                                                             
+      DO 29 I=1,7,2                                                     
+      B=-2.0*REAL(S(I))                                                 
+      Z1=S(I)*S(I+1)                                                    
+      C=REAL(Z1)                                                        
+      A=TDT+B+C/TDT                                                     
+      G=G*A                                                             
+      D(I)=(C*DT-FDT)/A                                                 
+   29 D(I+1)=(A-2.0*B)/A                                                
+      G=G*G                                                           
+    5 FORMAT ('-FILTER GAIN IS ', 9E12.6)                                 
+      RETURN                                                            
+
+      ENTRY FILTER(X,N,D,G,IG)                                          
+                                                                       
+!     X = DATA VECTOR OF LENGTH N CONTAINING DATA TO BE FILTERED        
+!     D = FILTER COEFFICIENTS CALCULATED BY BNDPAS                      
+!     G = FILTER GAIN                                                   
+!     IG = 1  one pass
+!     ig = 2  two passes
+                                                                       
+      IF (ISW.EQ.1) GO TO 31                                            
+      WRITE (6,6)                                                       
+    6 FORMAT ('1BNDPAS MUST BE CALLED BEFORE FILTER')                   
+      return                                                            
+                                                                       
+!     APPLY FILTER IN FORWARD DIRECTION                                 
+                                                                       
+   31 XM2=X(1)                                                          
+      XM1=X(2)                                                          
+      XM=X(3)                                                           
+      XC(1)=XM2                                                         
+      XC(2)=XM1-D(1)*XC(1)                                              
+      XC(3)=XM-XM2-D(1)*XC(2)-D(2)*XC(1)                                
+      XD(1)=XC(1)                                                       
+      XD(2)=XC(2)-D(3)*XD(1)                                            
+      XD(3)=XC(3)-XC(1)-D(3)*XD(2)-D(4)*XD(1)                           
+      XE(1)=XD(1)                                                       
+      XE(2)=XD(2)-D(5)*XE(1)                                            
+      XE(3)=XD(3)-XD(1)-D(5)*XE(2)-D(6)*XE(1)                           
+      X(1)=XE(1)                                                        
+      X(2)=XE(2)-D(7)*X(1)                                              
+      X(3)=XE(3)-XE(1)-D(7)*X(2)-D(8)*X(1)                              
+      DO 39 I=4,N                                                       
+      XM2=XM1                                                           
+      XM1=XM                                                            
+      XM=X(I)                                                           
+      K=I-((I-1)/3)*3                                                   
+      GO TO (34,35,36),K                                                
+   34 M=1                                                               
+      M1=3                                                              
+      M2=2                                                              
+      GO TO 37                                                          
+   35 M=2                                                               
+      M1=1                                                              
+      M2=3                                                              
+      GO TO 37                                                          
+   36 M=3                                                               
+      M1=2                                                              
+      M2=1                                                              
+   37 XC(M)=XM-XM2-D(1)*XC(M1)-D(2)*XC(M2)                              
+      XD(M)=XC(M)-XC(M2)-D(3)*XD(M1)-D(4)*XD(M2)                        
+      XE(M)=XD(M)-XD(M2)-D(5)*XE(M1)-D(6)*XE(M2)                        
+   39 X(I)=XE(M)-XE(M2)-D(7)*X(I-1)-D(8)*X(I-2)                         
+!                                                                       
+!
+      if(ig.eq.1) goto 3333                                             
+      XM2=X(N)                                                          
+      XM1=X(N-1)                                                        
+      XM=X(N-2)                                                         
+      XC(1)=XM2                                                         
+      XC(2)=XM1-D(1)*XC(1)                                              
+      XC(3)=XM-XM2-D(1)*XC(2)-D(2)*XC(1)                                
+      XD(1)=XC(1)                                                       
+      XD(2)=XC(2)-D(3)*XD(1)                                            
+      XD(3)=XC(3)-XC(1)-D(3)*XD(2)-D(4)*XD(1)                           
+      XE(1)=XD(1)                                                       
+      XE(2)=XD(2)-D(5)*XE(1)                                            
+      XE(3)=XD(3)-XD(1)-D(5)*XE(2)-D(6)*XE(1)                           
+      X(N)=XE(1)                                                        
+      X(N-1)=XE(2)-D(7)*X(1)                                            
+      X(N-2)=XE(3)-XE(1)-D(7)*X(2)-D(8)*X(1)                            
+      DO 49 I=4,N                                                       
+      XM2=XM1                                                           
+      XM1=XM                                                            
+      J=N-I+1                                                           
+      XM=X(J)                                                           
+      K=I-((I-1)/3)*3                                                   
+      GO TO (44,45,46),K                                                
+   44 M=1                                                               
+      M1=3                                                              
+      M2=2                                                              
+      GO TO 47                                                          
+   45 M=2                                                               
+      M1=1                                                              
+      M2=3                                                              
+      GO TO 47                                                          
+   46 M=3                                                               
+      M1=2                                                              
+      M2=1                                                              
+   47 XC(M)=XM-XM2-D(1)*XC(M1)-D(2)*XC(M2)                              
+      XD(M)=XC(M)-XC(M2)-D(3)*XD(M1)-D(4)*XD(M2)                        
+      XE(M)=XD(M)-XD(M2)-D(5)*XE(M1)-D(6)*XE(M2)                        
+   49 X(J)=XE(M)-XE(M2)-D(7)*X(J+1)-D(8)*X(J+2)                         
+ 3333 continue
+      if (ig.eq.1) then
+        gg=sqrt(g)   ! if only pass once, modify gain
+      else
+        gg=g
+      endif
+      DO 59 I=1,N                                                       
+   59 X(I)=X(I)/gg                                                      
+      RETURN                                                            
+END                                                               
+

Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/changemask
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/changemask	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/changemask	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,27 @@
+#!/bin/bash
+
+codefile=../../src/specfem2D/noise_tomography.f90
+
+cd masks
+select maskfile in `find -type f`
+do
+case $maskfile in
+    *) echo $maskfile; break;;
+esac
+done
+cd ..
+
+
+line1=`grep -n '^[ \t]*subroutine create_mask_noise' $codefile | cut -d':' -f1`
+line2=`grep -n '^[ \t]*end subroutine create_mask_noise' $codefile | cut -d':' -f1`
+
+if [ `echo $line1 | wc -l` -ne 1 ]; then echo "Error reading noise_tomography.f90"; exit 1; fi
+if [ `echo $line2 | wc -l` -ne 1 ]; then echo "Error reading noise_tomography.f90"; exit 2; fi
+if [ ! -n "$line1" ]; then echo "Error reading noise_tomography.f90"; exit 3; fi
+if [ ! -n "$line2" ]; then echo "Error reading noise_tomography.f90"; exit 4; fi
+
+awk -vline1=$line1 '{if (NR < line1) print $0}' < $codefile > temp1
+awk -vline2=$line2 '{if (NR > line2) print $0}' < $codefile > temp2
+
+cat temp1 masks/$maskfile temp2 > $codefile
+rm temp?


Property changes on: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/changemask
___________________________________________________________________
Name: svn:executable
   + *

Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/changerec
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/changerec	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/changerec	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,26 @@
+#!/bin/bash
+
+if [ $# -ne 2 ]; then echo "Usage:  change_rec  dir  nrec"; exit 1; fi
+if [ ! -d $1 ];  then echo "Directory not found: $1"; exit 1; fi
+if [ $2 -lt 1 ]; then echo "Need at least 1 receiver."; exit 1; fi
+
+DIR=$1
+NREC=$2
+
+cd $DIR
+
+for file in `ls Par_file*`
+do
+
+    echo $file
+
+    Line1=`grep "^nrec\s" $file`
+
+    if [ `echo $Line1 | wc -l` -ne 1 ]; then echo "Error reading file: $file."; exit 1; fi
+
+    sed -i -e "s:${Line1}:nrec                            = ${NREC}:"  $file
+
+done
+
+
+


Property changes on: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/changerec
___________________________________________________________________
Name: svn:executable
   + *

Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/clean
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/clean	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/clean	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,13 @@
+#!/bin/bash
+#
+# runs all three noise tomography simulations
+#
+
+echo "Cleaning up."
+
+rm -rf DATA SEM NOISE_TOMOGRAPHY OUTPUT_FILES OUTPUT_ALL &> /dev/null
+rm xmeshfem2D xspecfem2D &> /dev/null
+rm adj_run &> /dev/null
+rm log.* &> /dev/null
+
+


Property changes on: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/clean
___________________________________________________________________
Name: svn:executable
   + *

Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/getdirect
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/getdirect	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/getdirect	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,10 @@
+#!/bin/sh
+NREC=$1
+ISRC=$2
+IREC=$3
+
+XINC=`echo "scale=3; (2000-300)/($NREC-1)" | bc`
+IOFF=`echo "scale=0; $IREC - $ISRC" | bc`
+XOFF=`echo "scale=3; $IOFF*$XINC" | bc`
+TOFF=`echo "scale=3; ($XOFF)/1300" | bc`
+echo $TOFF


Property changes on: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/getdirect
___________________________________________________________________
Name: svn:executable
   + *

Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/masks/drillbit
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/masks/drillbit	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/masks/drillbit	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,26 @@
+  subroutine create_mask_noise(nglob,coord,mask_noise)
+
+  implicit none
+  include "constants.h"
+
+  !input
+  integer :: nglob
+  real(kind=CUSTOM_REAL), dimension(2,nglob) :: coord
+
+  !output
+  real(kind=CUSTOM_REAL), dimension(nglob) :: mask_noise
+
+  !local
+  integer :: iglob
+  real(kind=CUSTOM_REAL) :: xx,zz
+
+  !specify distribution of noise sources as a function of xx, zz
+  do iglob = 1,nglob
+    xx = coord(1,iglob)
+    zz = coord(2,iglob)
+
+    mask_noise(iglob) = 1.e24 * 30.5577 * exp(-((xx-1250.)**2 + (zz-750)**2)/(2.*150**2) )
+
+  enddo
+
+  end subroutine create_mask_noise

Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/masks/uniform
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/masks/uniform	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/masks/uniform	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,28 @@
+  subroutine create_mask_noise(nglob,coord,mask_noise)
+
+  implicit none
+  include "constants.h"
+
+  !input
+  integer :: nglob
+  real(kind=CUSTOM_REAL), dimension(2,nglob) :: coord
+
+  !output
+  real(kind=CUSTOM_REAL), dimension(nglob) :: mask_noise
+
+  !local
+  integer :: iglob
+  real(kind=CUSTOM_REAL) :: xx,zz
+
+  !specify distribution of noise sources as a function of xx, zz
+  do iglob = 1,nglob
+    xx = coord(1,iglob)
+    zz = coord(2,iglob)
+
+    !below, the noise is assumed to be uniform; users are free to
+    !to change this expression to one involving xx, zz
+    mask_noise(iglob) = 1.
+
+  enddo
+
+  end subroutine create_mask_noise

Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_0/Par_file_fair
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_0/Par_file_fair	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_0/Par_file_fair	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,134 @@
+# title of job
+title                           = Passive Imaging - 0 Interfaces
+
+# forward or adjoint simulation
+SIMULATION_TYPE                 = 1
+NOISE_TOMOGRAPHY                = 1
+SAVE_FORWARD                    = .false.
+
+# parameters concerning partitioning
+nproc                           = 1              # number of processes
+partitioning_method             = 3              # SCOTCH = 3, ascending order (very bad idea) = 1
+PERFORM_CUTHILL_MCKEE           = .true.         # perform inverse Cuthill-McKee (1969) optimization/permutation for mesh numbering
+
+ngnod                           = 9              # number of control nodes per element (4 or 9)
+initialfield                    = .false.        # use a plane wave as source or not
+add_Bielak_conditions           = .false.        # add Bielak conditions or not if initial plane wave
+assign_external_model           = .false.        # define external earth model or not
+READ_EXTERNAL_SEP_FILE          = .false.        # Read external SEP file from DATA/model_velocity.dat_input, or use routine
+ATTENUATION_VISCOELASTIC_SOLID  = .false.        # turn attenuation (viscoelasticity) on or off for non-poroelastic solid parts of the model
+ATTENUATION_PORO_FLUID_PART     = .false.        # turn viscous attenuation on or off for the fluid part of poroelastic parts of the model
+Q0                              =  1             # quality factor for viscous attenuation
+freq0                           =  10            # frequency for viscous attenuation
+p_sv                            = .false.        # set the type of calculation (P-SV or SH/membrane waves)
+
+# time step parameters
+nt                              = 3000           # total number of time steps
+deltat                          = 1.d-3          # duration of a time step
+USER_T0                         = 0.0d0          # use this t0 as earliest starting time rather than the automatically calculated one
+time_stepping_scheme            = 1   # 1 = Newmark (2nd order),     2 = LDDRK4-6 (4th-order 6-stage low storage Runge-Kutta),     3 = classical 4th-order 4-stage Runge-Kutta
+
+# source parameters
+NSOURCES                        = 1              # number of sources [source info read in CMTSOLUTION file]
+force_normal_to_surface         = .false.        # angleforce normal to surface (external mesh and curve file needed)
+
+# constants for attenuation
+N_SLS                           = 2                      # number of standard linear solids for attenuation
+f0_attenuation                  = 5.196152422706633      # (Hz) relevant only if source is a Dirac or a Heaviside, else it is f0
+
+# receiver set parameters for seismograms
+seismotype                      = 1              # record 1=displ 2=veloc 3=accel 4=pressure
+generate_STATIONS               = .true.         # creates a STATION file in ./DATA
+nreceiversets                   = 1              # number of receiver sets
+anglerec                        = 0.d0           # angle to rotate components at receivers
+rec_normal_to_surface           = .false.        # base anglerec normal to surface (external mesh and curve file needed)
+SU_FORMAT                       = .false.        # output seismograms in Seismic Unix format (adjoint traces will be read in the same format)
+
+# first receiver set (repeat these 6 lines and adjust nreceiversets  accordingly)
+nrec                            = 29
+xdeb                            = 150.           # first receiver x in meters
+zdeb                            = 1500.          # first receiver z in meters
+xfin                            = 1850.          # last receiver x in meters (ignored if onlyone receiver)
+zfin                            = 1500.          # last receiver z in meters (ignored if onlyone receiver)
+enreg_surf_same_vertical        = .true.         # receivers inside the medium or at the surface
+
+# display parameters
+NTSTEP_BETWEEN_OUTPUT_INFO      = 250            # display frequency in time steps
+output_postscript_snapshot      = .false.        # output Postscript snapshot of the results
+output_color_image              = .false.        # output color image of the results
+imagetype                       = 1              # display 1=displ 2=veloc 3=accel 4=pressure
+cutsnaps                        = 1.             # minimum amplitude in % for snapshots
+meshvect                        = .true.         # display mesh on vector plots or not
+modelvect                       = .false.        # display velocity model on vector plots
+boundvect                       = .true.         # display boundary conditions on plots
+interpol                        = .true.         # interpolation of the display or not
+pointsdisp                      = 6              # points for interpolation of display (set to 1 for lower-left corner only)
+subsamp_postscript              = 1              # subsampling of color snapshots
+factor_subsample_image          = 1              # factor to subsample color images output by the code (useful for very large models)
+POWER_DISPLAY_COLOR             = 0.30d0         # non linear display to enhance small amplitudes in color images
+DRAW_WATER_CONSTANT_BLUE_IN_JPG = .true.         # display acoustic layers as constant blue in JPEG images, because they likely correspond to water
+sizemax_arrows                  = 1.d0           # maximum size of arrows on vector plots in cm
+US_LETTER                       = .false.        # US letter paper or European A4
+USE_SNAPSHOT_NUMBER_IN_FILENAME = .false.        # use snapshot number in the file name of JPEG color snapshots instead of the time step
+gnuplot                         = .false.        # generate a GNUPLOT file for the grid
+output_grid                     = .false.        # save the grid in a text file or not
+output_energy                   = .false.        # compute and output acoustic and elastic energy (slows down the code significantly)
+output_wavefield_snapshot       = .false.         # output Ux,Uy,Uz text file for each output time (big files)
+
+# velocity and density models
+nbmodels                        = 1              # nb of different models
+# define models as
+# I:   (model_number 1 rho Vp Vs 0 0 Qkappa Qmu 0 0 0 0 0 0) or
+# II:  (model_number 2 rho c11 c13 c15 c33 c35 c55 0 0 0 0 0 0) or
+# III: (model_number 3 rhos rhof phi c kxx kxz kzz Ks Kf Kfr etaf mufr Qmu).
+# For istropic elastic/acoustic material use I and set Vs to zero to make a given model acoustic, for anisotropic elastic use II,
+# and for isotropic poroelastic material use III. The mesh can contain acoustic, elastic, and poroelastic models simultaneously.
+1 1 2500.d0 2700.d0 1600.000d0 0 0 9999 9999 0 0 0 0 0 0
+
+# external mesh or not
+read_external_mesh              = .false.
+
+# absorbing boundary active or not
+absorbing_conditions            = .true.
+
+# for horizontal periodic conditions: detect common points between left and right edges
+ADD_PERIODIC_CONDITIONS         = .false.
+
+# horizontal periodicity distance for periodic conditions
+PERIODIC_horiz_dist             = 0.3597d0
+
+# grid point detection tolerance for periodic conditions
+PERIODIC_DETECT_TOL             = 3.3334d-6
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR EXTERNAL MESHING
+
+# data concerning mesh, when generated using third-party app (more info in README)
+# (see also absorbing_conditions above)
+mesh_file                       = ./DATA/Mesh_canyon/canyon_mesh_file   # file containing the mesh
+nodes_coords_file               = ./DATA/Mesh_canyon/canyon_nodes_coords_file   # file containing the nodes coordinates
+materials_file                  = ./DATA/Mesh_canyon/canyon_materials_file   # file containing the material number for each element
+free_surface_file               = ./DATA/Mesh_canyon/canyon_free_surface_file   # file containing the free surface
+absorbing_surface_file          = ./DATA/Mesh_canyon/canyon_absorbing_surface_file   # file containing the absorbing surface
+tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the curve delimiting the velocity model
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR INTERNAL MESHING
+
+# file containing interfaces for internal mesh
+interfacesfile                  = interfaces.dat
+
+# geometry of the model (origin lower-left corner = 0,0) and mesh description
+xmin                            = 0.d0           # abscissa of left side of the model
+xmax                            = 2000.d0        # abscissa of right side of the model
+nx                              = 40             # number of elements along X
+
+# absorbing boundary parameters (see absorbing_conditions above)
+absorbbottom                    = .true.
+absorbright                     = .true.
+absorbtop                       = .false.
+absorbleft                      = .true.
+
+# define the different regions of the model in the (nx,nz) spectral element mesh
+nbregions                       = 1              # nb of regions and model number for each
+1 40  1 30 1

Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_0/Par_file_fair.before_update_to_r19xxx
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_0/Par_file_fair.before_update_to_r19xxx	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_0/Par_file_fair.before_update_to_r19xxx	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,133 @@
+# title of job
+title                           = Passive Imaging - 0 Interfaces
+
+# forward or adjoint simulation
+SIMULATION_TYPE                 = 1
+NOISE_TOMOGRAPHY                = 1
+SAVE_FORWARD                    = .false.
+
+# parameters concerning partitioning
+nproc                           = 1              # number of processes
+partitioning_method             = 3              # SCOTCH = 3, ascending order (very bad idea) = 1
+PERFORM_CUTHILL_MCKEE           = .true.         # perform inverse Cuthill-McKee (1969) optimization/permutation for mesh numbering
+
+ngnod                           = 9              # number of control nodes per element (4 or 9)
+initialfield                    = .false.        # use a plane wave as source or not
+add_Bielak_conditions           = .false.        # add Bielak conditions or not if initial plane wave
+assign_external_model           = .false.        # define external earth model or not
+READ_EXTERNAL_SEP_FILE          = .false.        # Read external SEP file from DATA/model_velocity.dat_input, or use routine
+ATTENUATION_VISCOELASTIC_SOLID  = .false.        # turn attenuation (viscoelasticity) on or off for non-poroelastic solid parts of the model
+ATTENUATION_PORO_FLUID_PART     = .false.        # turn viscous attenuation on or off for the fluid part of poroelastic parts of the model
+Q0                              =  1             # quality factor for viscous attenuation
+freq0                           =  10            # frequency for viscous attenuation
+p_sv                            = .false.        # set the type of calculation (P-SV or SH/membrane waves)
+
+# time step parameters
+nt                              = 3000           # total number of time steps
+deltat                          = 1.d-3          # duration of a time step
+USER_T0                         = 0.0d0          # use this t0 as earliest starting time rather than the automatically calculated one
+
+# source parameters
+NSOURCES                        = 1              # number of sources [source info read in CMTSOLUTION file]
+force_normal_to_surface         = .false.        # angleforce normal to surface (external mesh and curve file needed)
+
+# constants for attenuation
+N_SLS                           = 2                      # number of standard linear solids for attenuation
+f0_attenuation                  = 5.196152422706633      # (Hz) relevant only if source is a Dirac or a Heaviside, else it is f0
+
+# receiver set parameters for seismograms
+seismotype                      = 1              # record 1=displ 2=veloc 3=accel 4=pressure
+generate_STATIONS               = .true.         # creates a STATION file in ./DATA
+nreceiversets                   = 1              # number of receiver sets
+anglerec                        = 0.d0           # angle to rotate components at receivers
+rec_normal_to_surface           = .false.        # base anglerec normal to surface (external mesh and curve file needed)
+SU_FORMAT                       = .false.        # output seismograms in Seismic Unix format (adjoint traces will be read in the same format)
+
+# first receiver set (repeat these 6 lines and adjust nreceiversets  accordingly)
+nrec                            = 29
+xdeb                            = 150.           # first receiver x in meters
+zdeb                            = 1500.          # first receiver z in meters
+xfin                            = 1850.          # last receiver x in meters (ignored if onlyone receiver)
+zfin                            = 1500.          # last receiver z in meters (ignored if onlyone receiver)
+enreg_surf_same_vertical        = .true.         # receivers inside the medium or at the surface
+
+# display parameters
+NTSTEP_BETWEEN_OUTPUT_INFO      = 250            # display frequency in time steps
+output_postscript_snapshot      = .false.        # output Postscript snapshot of the results
+output_color_image              = .false.        # output color image of the results
+imagetype                       = 1              # display 1=displ 2=veloc 3=accel 4=pressure
+cutsnaps                        = 1.             # minimum amplitude in % for snapshots
+meshvect                        = .true.         # display mesh on vector plots or not
+modelvect                       = .false.        # display velocity model on vector plots
+boundvect                       = .true.         # display boundary conditions on plots
+interpol                        = .true.         # interpolation of the display or not
+pointsdisp                      = 6              # points for interpolation of display (set to 1 for lower-left corner only)
+subsamp_postscript              = 1              # subsampling of color snapshots
+factor_subsample_image          = 1              # factor to subsample color images output by the code (useful for very large models)
+POWER_DISPLAY_COLOR             = 0.30d0         # non linear display to enhance small amplitudes in color images
+DRAW_WATER_CONSTANT_BLUE_IN_JPG = .true.         # display acoustic layers as constant blue in JPEG images, because they likely correspond to water
+sizemax_arrows                  = 1.d0           # maximum size of arrows on vector plots in cm
+US_LETTER                       = .false.        # US letter paper or European A4
+USE_SNAPSHOT_NUMBER_IN_FILENAME = .false.        # use snapshot number in the file name of JPEG color snapshots instead of the time step
+gnuplot                         = .false.        # generate a GNUPLOT file for the grid
+output_grid                     = .false.        # save the grid in a text file or not
+output_energy                   = .false.        # compute and output acoustic and elastic energy (slows down the code significantly)
+output_wavefield_snapshot       = .false.         # output Ux,Uy,Uz text file for each output time (big files)
+
+# velocity and density models
+nbmodels                        = 1              # nb of different models
+# define models as
+# I:   (model_number 1 rho Vp Vs 0 0 Qkappa Qmu 0 0 0 0 0 0) or
+# II:  (model_number 2 rho c11 c13 c15 c33 c35 c55 0 0 0 0 0 0) or
+# III: (model_number 3 rhos rhof phi c kxx kxz kzz Ks Kf Kfr etaf mufr Qmu).
+# For istropic elastic/acoustic material use I and set Vs to zero to make a given model acoustic, for anisotropic elastic use II,
+# and for isotropic poroelastic material use III. The mesh can contain acoustic, elastic, and poroelastic models simultaneously.
+1 1 2500.d0 2700.d0 1600.000d0 0 0 9999 9999 0 0 0 0 0 0
+
+# external mesh or not
+read_external_mesh              = .false.
+
+# absorbing boundary active or not
+absorbing_conditions            = .true.
+
+# for horizontal periodic conditions: detect common points between left and right edges
+ADD_PERIODIC_CONDITIONS         = .false.
+
+# horizontal periodicity distance for periodic conditions
+PERIODIC_horiz_dist             = 0.3597d0
+
+# grid point detection tolerance for periodic conditions
+PERIODIC_DETECT_TOL             = 3.3334d-6
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR EXTERNAL MESHING
+
+# data concerning mesh, when generated using third-party app (more info in README)
+# (see also absorbing_conditions above)
+mesh_file                       = ./DATA/Mesh_canyon/canyon_mesh_file   # file containing the mesh
+nodes_coords_file               = ./DATA/Mesh_canyon/canyon_nodes_coords_file   # file containing the nodes coordinates
+materials_file                  = ./DATA/Mesh_canyon/canyon_materials_file   # file containing the material number for each element
+free_surface_file               = ./DATA/Mesh_canyon/canyon_free_surface_file   # file containing the free surface
+absorbing_surface_file          = ./DATA/Mesh_canyon/canyon_absorbing_surface_file   # file containing the absorbing surface
+tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the curve delimiting the velocity model
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR INTERNAL MESHING
+
+# file containing interfaces for internal mesh
+interfacesfile                  = interfaces.dat
+
+# geometry of the model (origin lower-left corner = 0,0) and mesh description
+xmin                            = 0.d0           # abscissa of left side of the model
+xmax                            = 2000.d0        # abscissa of right side of the model
+nx                              = 40             # number of elements along X
+
+# absorbing boundary parameters (see absorbing_conditions above)
+absorbbottom                    = .true.
+absorbright                     = .true.
+absorbtop                       = .false.
+absorbleft                      = .true.
+
+# define the different regions of the model in the (nx,nz) spectral element mesh
+nbregions                       = 1              # nb of regions and model number for each
+1 40  1 30 1

Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_0/Par_file_good
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_0/Par_file_good	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_0/Par_file_good	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,134 @@
+s# title of job
+4itle                           = Passive Imaging - 0 Interfaces
+
+# forward or adjoint simulation
+SIMULATION_TYPE                 = 1
+NOISE_TOMOGRAPHY                = 1
+SAVE_FORWARD                    = .false.
+
+# parameters concerning partitioning
+nproc                           = 1              # number of processes
+partitioning_method             = 3              # SCOTCH = 3, ascending order (very bad idea) = 1
+PERFORM_CUTHILL_MCKEE           = .true.         # perform inverse Cuthill-McKee (1969) optimization/permutation for mesh numbering
+
+ngnod                           = 9              # number of control nodes per element (4 or 9)
+initialfield                    = .false.        # use a plane wave as source or not
+add_Bielak_conditions           = .false.        # add Bielak conditions or not if initial plane wave
+assign_external_model           = .false.        # define external earth model or not
+READ_EXTERNAL_SEP_FILE          = .false.        # Read external SEP file from DATA/model_velocity.dat_input, or use routine
+ATTENUATION_VISCOELASTIC_SOLID  = .false.        # turn attenuation (viscoelasticity) on or off for non-poroelastic solid parts of the model
+ATTENUATION_PORO_FLUID_PART     = .false.        # turn viscous attenuation on or off for the fluid part of poroelastic parts of the model
+Q0                              =  1             # quality factor for viscous attenuation
+freq0                           =  10            # frequency for viscous attenuation
+p_sv                            = .false.        # set the type of calculation (P-SV or SH/membrane waves)
+
+# time step parameters
+nt                              = 3000           # total number of time steps
+deltat                          = 1.d-3          # duration of a time step
+USER_T0                         = 0.0d0          # use this t0 as earliest starting time rather than the automatically calculated one
+time_stepping_scheme            = 1   # 1 = Newmark (2nd order),     2 = LDDRK4-6 (4th-order 6-stage low storage Runge-Kutta),     3 = classical 4th-order 4-stage Runge-Kutta
+
+# source parameters
+NSOURCES                        = 1              # number of sources [source info read in CMTSOLUTION file]
+force_normal_to_surface         = .false.        # angleforce normal to surface (external mesh and curve file needed)
+
+# constants for attenuation
+N_SLS                           = 2                      # number of standard linear solids for attenuation
+f0_attenuation                  = 5.196152422706633      # (Hz) relevant only if source is a Dirac or a Heaviside, else it is f0
+
+# receiver set parameters for seismograms
+seismotype                      = 1              # record 1=displ 2=veloc 3=accel 4=pressure
+generate_STATIONS               = .true.         # creates a STATION file in ./DATA
+nreceiversets                   = 1              # number of receiver sets
+anglerec                        = 0.d0           # angle to rotate components at receivers
+rec_normal_to_surface           = .false.        # base anglerec normal to surface (external mesh and curve file needed)
+SU_FORMAT                       = .false.        # output seismograms in Seismic Unix format (adjoint traces will be read in the same format)
+
+# first receiver set (repeat these 6 lines and adjust nreceiversets  accordingly)
+nrec                            = 29
+xdeb                            = 400.           # first receiver x in meters
+zdeb                            = 1500.          # first receiver z in meters
+xfin                            = 2100.          # last receiver x in meters (ignored if onlyone receiver)
+zfin                            = 1500.          # last receiver z in meters (ignored if onlyone receiver)
+enreg_surf_same_vertical        = .true.         # receivers inside the medium or at the surface
+
+# display parameters
+NTSTEP_BETWEEN_OUTPUT_INFO      = 250            # display frequency in time steps
+output_postscript_snapshot      = .false.        # output Postscript snapshot of the results
+output_color_image              = .false.        # output color image of the results
+imagetype                       = 1              # display 1=displ 2=veloc 3=accel 4=pressure
+cutsnaps                        = 1.             # minimum amplitude in % for snapshots
+meshvect                        = .true.         # display mesh on vector plots or not
+modelvect                       = .false.        # display velocity model on vector plots
+boundvect                       = .true.         # display boundary conditions on plots
+interpol                        = .true.         # interpolation of the display or not
+pointsdisp                      = 6              # points for interpolation of display (set to 1 for lower-left corner only)
+subsamp_postscript              = 1              # subsampling of color snapshots
+factor_subsample_image          = 1              # factor to subsample color images output by the code (useful for very large models)
+POWER_DISPLAY_COLOR             = 0.30d0         # non linear display to enhance small amplitudes in color images
+DRAW_WATER_CONSTANT_BLUE_IN_JPG = .true.         # display acoustic layers as constant blue in JPEG images, because they likely correspond to water
+sizemax_arrows                  = 1.d0           # maximum size of arrows on vector plots in cm
+US_LETTER                       = .false.        # US letter paper or European A4
+USE_SNAPSHOT_NUMBER_IN_FILENAME = .false.        # use snapshot number in the file name of JPEG color snapshots instead of the time step
+gnuplot                         = .false.        # generate a GNUPLOT file for the grid
+output_grid                     = .false.        # save the grid in a text file or not
+output_energy                   = .false.        # compute and output acoustic and elastic energy (slows down the code significantly)
+output_wavefield_snapshot       = .false.         # output Ux,Uy,Uz text file for each output time (big files)
+
+# velocity and density models
+nbmodels                        = 1              # nb of different models
+# define models as
+# I:   (model_number 1 rho Vp Vs 0 0 Qkappa Qmu 0 0 0 0 0 0) or
+# II:  (model_number 2 rho c11 c13 c15 c33 c35 c55 0 0 0 0 0 0) or
+# III: (model_number 3 rhos rhof phi c kxx kxz kzz Ks Kf Kfr etaf mufr Qmu).
+# For istropic elastic/acoustic material use I and set Vs to zero to make a given model acoustic, for anisotropic elastic use II,
+# and for isotropic poroelastic material use III. The mesh can contain acoustic, elastic, and poroelastic models simultaneously.
+1 1 2500.d0 2700.d0 1600.000d0 0 0 9999 9999 0 0 0 0 0 0
+
+# external mesh or not
+read_external_mesh              = .false.
+
+# absorbing boundary active or not
+absorbing_conditions            = .true.
+
+# for horizontal periodic conditions: detect common points between left and right edges
+ADD_PERIODIC_CONDITIONS         = .false.
+
+# horizontal periodicity distance for periodic conditions
+PERIODIC_horiz_dist             = 0.3597d0
+
+# grid point detection tolerance for periodic conditions
+PERIODIC_DETECT_TOL             = 3.3334d-6
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR EXTERNAL MESHING
+
+# data concerning mesh, when generated using third-party app (more info in README)
+# (see also absorbing_conditions above)
+mesh_file                       = ./DATA/Mesh_canyon/canyon_mesh_file   # file containing the mesh
+nodes_coords_file               = ./DATA/Mesh_canyon/canyon_nodes_coords_file   # file containing the nodes coordinates
+materials_file                  = ./DATA/Mesh_canyon/canyon_materials_file   # file containing the material number for each element
+free_surface_file               = ./DATA/Mesh_canyon/canyon_free_surface_file   # file containing the free surface
+absorbing_surface_file          = ./DATA/Mesh_canyon/canyon_absorbing_surface_file   # file containing the absorbing surface
+tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the curve delimiting the velocity model
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR INTERNAL MESHING
+
+# file containing interfaces for internal mesh
+interfacesfile                  = interfaces.dat
+
+# geometry of the model (origin lower-left corner = 0,0) and mesh description
+xmin                            = 0.d0           # abscissa of left side of the model
+xmax                            = 2500.d0        # abscissa of right side of the model
+nx                              = 50             # number of elements along X
+
+# absorbing boundary parameters (see absorbing_conditions above)
+absorbbottom                    = .true.
+absorbright                     = .true.
+absorbtop                       = .false.
+absorbleft                      = .true.
+
+# define the different regions of the model in the (nx,nz) spectral element mesh
+nbregions                       = 1              # nb of regions and model number for each
+1 50  1 30 1

Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_0/Par_file_good.before_update_to_r19xxx
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_0/Par_file_good.before_update_to_r19xxx	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_0/Par_file_good.before_update_to_r19xxx	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,133 @@
+s# title of job
+4itle                           = Passive Imaging - 0 Interfaces
+
+# forward or adjoint simulation
+SIMULATION_TYPE                 = 1
+NOISE_TOMOGRAPHY                = 1
+SAVE_FORWARD                    = .false.
+
+# parameters concerning partitioning
+nproc                           = 1              # number of processes
+partitioning_method             = 3              # SCOTCH = 3, ascending order (very bad idea) = 1
+PERFORM_CUTHILL_MCKEE           = .true.         # perform inverse Cuthill-McKee (1969) optimization/permutation for mesh numbering
+
+ngnod                           = 9              # number of control nodes per element (4 or 9)
+initialfield                    = .false.        # use a plane wave as source or not
+add_Bielak_conditions           = .false.        # add Bielak conditions or not if initial plane wave
+assign_external_model           = .false.        # define external earth model or not
+READ_EXTERNAL_SEP_FILE          = .false.        # Read external SEP file from DATA/model_velocity.dat_input, or use routine
+ATTENUATION_VISCOELASTIC_SOLID  = .false.        # turn attenuation (viscoelasticity) on or off for non-poroelastic solid parts of the model
+ATTENUATION_PORO_FLUID_PART     = .false.        # turn viscous attenuation on or off for the fluid part of poroelastic parts of the model
+Q0                              =  1             # quality factor for viscous attenuation
+freq0                           =  10            # frequency for viscous attenuation
+p_sv                            = .false.        # set the type of calculation (P-SV or SH/membrane waves)
+
+# time step parameters
+nt                              = 3000           # total number of time steps
+deltat                          = 1.d-3          # duration of a time step
+USER_T0                         = 0.0d0          # use this t0 as earliest starting time rather than the automatically calculated one
+
+# source parameters
+NSOURCES                        = 1              # number of sources [source info read in CMTSOLUTION file]
+force_normal_to_surface         = .false.        # angleforce normal to surface (external mesh and curve file needed)
+
+# constants for attenuation
+N_SLS                           = 2                      # number of standard linear solids for attenuation
+f0_attenuation                  = 5.196152422706633      # (Hz) relevant only if source is a Dirac or a Heaviside, else it is f0
+
+# receiver set parameters for seismograms
+seismotype                      = 1              # record 1=displ 2=veloc 3=accel 4=pressure
+generate_STATIONS               = .true.         # creates a STATION file in ./DATA
+nreceiversets                   = 1              # number of receiver sets
+anglerec                        = 0.d0           # angle to rotate components at receivers
+rec_normal_to_surface           = .false.        # base anglerec normal to surface (external mesh and curve file needed)
+SU_FORMAT                       = .false.        # output seismograms in Seismic Unix format (adjoint traces will be read in the same format)
+
+# first receiver set (repeat these 6 lines and adjust nreceiversets  accordingly)
+nrec                            = 29
+xdeb                            = 400.           # first receiver x in meters
+zdeb                            = 1500.          # first receiver z in meters
+xfin                            = 2100.          # last receiver x in meters (ignored if onlyone receiver)
+zfin                            = 1500.          # last receiver z in meters (ignored if onlyone receiver)
+enreg_surf_same_vertical        = .true.         # receivers inside the medium or at the surface
+
+# display parameters
+NTSTEP_BETWEEN_OUTPUT_INFO      = 250            # display frequency in time steps
+output_postscript_snapshot      = .false.        # output Postscript snapshot of the results
+output_color_image              = .false.        # output color image of the results
+imagetype                       = 1              # display 1=displ 2=veloc 3=accel 4=pressure
+cutsnaps                        = 1.             # minimum amplitude in % for snapshots
+meshvect                        = .true.         # display mesh on vector plots or not
+modelvect                       = .false.        # display velocity model on vector plots
+boundvect                       = .true.         # display boundary conditions on plots
+interpol                        = .true.         # interpolation of the display or not
+pointsdisp                      = 6              # points for interpolation of display (set to 1 for lower-left corner only)
+subsamp_postscript              = 1              # subsampling of color snapshots
+factor_subsample_image          = 1              # factor to subsample color images output by the code (useful for very large models)
+POWER_DISPLAY_COLOR             = 0.30d0         # non linear display to enhance small amplitudes in color images
+DRAW_WATER_CONSTANT_BLUE_IN_JPG = .true.         # display acoustic layers as constant blue in JPEG images, because they likely correspond to water
+sizemax_arrows                  = 1.d0           # maximum size of arrows on vector plots in cm
+US_LETTER                       = .false.        # US letter paper or European A4
+USE_SNAPSHOT_NUMBER_IN_FILENAME = .false.        # use snapshot number in the file name of JPEG color snapshots instead of the time step
+gnuplot                         = .false.        # generate a GNUPLOT file for the grid
+output_grid                     = .false.        # save the grid in a text file or not
+output_energy                   = .false.        # compute and output acoustic and elastic energy (slows down the code significantly)
+output_wavefield_snapshot       = .false.         # output Ux,Uy,Uz text file for each output time (big files)
+
+# velocity and density models
+nbmodels                        = 1              # nb of different models
+# define models as
+# I:   (model_number 1 rho Vp Vs 0 0 Qkappa Qmu 0 0 0 0 0 0) or
+# II:  (model_number 2 rho c11 c13 c15 c33 c35 c55 0 0 0 0 0 0) or
+# III: (model_number 3 rhos rhof phi c kxx kxz kzz Ks Kf Kfr etaf mufr Qmu).
+# For istropic elastic/acoustic material use I and set Vs to zero to make a given model acoustic, for anisotropic elastic use II,
+# and for isotropic poroelastic material use III. The mesh can contain acoustic, elastic, and poroelastic models simultaneously.
+1 1 2500.d0 2700.d0 1600.000d0 0 0 9999 9999 0 0 0 0 0 0
+
+# external mesh or not
+read_external_mesh              = .false.
+
+# absorbing boundary active or not
+absorbing_conditions            = .true.
+
+# for horizontal periodic conditions: detect common points between left and right edges
+ADD_PERIODIC_CONDITIONS         = .false.
+
+# horizontal periodicity distance for periodic conditions
+PERIODIC_horiz_dist             = 0.3597d0
+
+# grid point detection tolerance for periodic conditions
+PERIODIC_DETECT_TOL             = 3.3334d-6
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR EXTERNAL MESHING
+
+# data concerning mesh, when generated using third-party app (more info in README)
+# (see also absorbing_conditions above)
+mesh_file                       = ./DATA/Mesh_canyon/canyon_mesh_file   # file containing the mesh
+nodes_coords_file               = ./DATA/Mesh_canyon/canyon_nodes_coords_file   # file containing the nodes coordinates
+materials_file                  = ./DATA/Mesh_canyon/canyon_materials_file   # file containing the material number for each element
+free_surface_file               = ./DATA/Mesh_canyon/canyon_free_surface_file   # file containing the free surface
+absorbing_surface_file          = ./DATA/Mesh_canyon/canyon_absorbing_surface_file   # file containing the absorbing surface
+tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the curve delimiting the velocity model
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR INTERNAL MESHING
+
+# file containing interfaces for internal mesh
+interfacesfile                  = interfaces.dat
+
+# geometry of the model (origin lower-left corner = 0,0) and mesh description
+xmin                            = 0.d0           # abscissa of left side of the model
+xmax                            = 2500.d0        # abscissa of right side of the model
+nx                              = 50             # number of elements along X
+
+# absorbing boundary parameters (see absorbing_conditions above)
+absorbbottom                    = .true.
+absorbright                     = .true.
+absorbtop                       = .false.
+absorbleft                      = .true.
+
+# define the different regions of the model in the (nx,nz) spectral element mesh
+nbregions                       = 1              # nb of regions and model number for each
+1 50  1 30 1

Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_0/SOURCE_noise
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_0/SOURCE_noise	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_0/SOURCE_noise	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,13 @@
+#source 1.  The components of a moment tensor source must be given in N.m, not in dyne.cm as in the DATA/CMTSOLUTION source file of the 3D version of the code.
+source_surf                     = .false.        # source inside the medium or at the surface
+xs                              = 0.
+zs                              = 0.
+source_type                     = 1              # elastic force or acoustic pressure = 1 or moment tensor = 2
+time_function_type              = 3              # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4, Heaviside = 5
+f0                              = 10.0           # dominant source frequency (Hz) if not Dirac or Heaviside
+tshift                          = 0.0            # time shift when multi sources (if one source, must be zero)
+angleforce                      = 0.             # angle of the source (for a force only)
+Mxx                             = 0.d0           # Mxx component (for a moment tensor source only)
+Mzz                             = 0.d0           # Mzz component (for a moment tensor source only)
+Mxz                             = 0.d0           # Mxz component (for a moment tensor source only)
+factor                          = 0.d0          # amplification factor

Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_0/interfaces_fair
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_0/interfaces_fair	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_0/interfaces_fair	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,26 @@
+#
+# number of interfaces
+#
+ 2
+#
+# for each interface below, we give the number of points and then x,z for each point
+#
+#
+# interface number 1 (bottom of the mesh)
+#
+ 2
+ 0 0
+ 2000 0
+#
+# interface number 2 (topography, top of the mesh)
+#
+ 2
+    0 1500
+ 2000 1500
+#
+# for each layer, we give the number of spectral elements in the vertical direction
+#
+#
+# layer number 1
+#
+ 30

Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_0/interfaces_good
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_0/interfaces_good	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_0/interfaces_good	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,26 @@
+#
+# number of interfaces
+#
+ 2
+#
+# for each interface below, we give the number of points and then x,z for each point
+#
+#
+# interface number 1 (bottom of the mesh)
+#
+ 2
+ 0 0
+ 2500 0
+#
+# interface number 2 (topography, top of the mesh)
+#
+ 2
+    0 1500
+ 2500 1500
+#
+# for each layer, we give the number of spectral elements in the vertical direction
+#
+#
+# layer number 1
+#
+ 30

Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/Par_file_best
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/Par_file_best	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/Par_file_best	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,136 @@
+# title of job
+title                           = PASSIVE IMAGING - 1 INTERFACE
+
+# forward or adjoint simulation
+SIMULATION_TYPE                 = 1
+NOISE_TOMOGRAPHY                = 1
+SAVE_FORWARD                    = .false.
+
+# parameters concerning partitioning
+nproc                           = 1              # number of processes
+partitioning_method             = 3              # SCOTCH = 3, ascending order (very bad idea) = 1
+PERFORM_CUTHILL_MCKEE           = .true.         # perform inverse Cuthill-McKee (1969) optimization/permutation for mesh numbering
+
+ngnod                           = 9              # number of control nodes per element (4 or 9)
+initialfield                    = .false.        # use a plane wave as source or not
+add_Bielak_conditions           = .false.        # add Bielak conditions or not if initial plane wave
+assign_external_model           = .false.        # define external earth model or not
+READ_EXTERNAL_SEP_FILE          = .false.        # Read external SEP file from DATA/model_velocity.dat_input, or use routine
+ATTENUATION_VISCOELASTIC_SOLID  = .false.        # turn attenuation (viscoelasticity) on or off for non-poroelastic solid parts of the model
+ATTENUATION_PORO_FLUID_PART     = .false.        # turn viscous attenuation on or off for the fluid part of poroelastic parts of the model
+Q0                              =  1             # quality factor for viscous attenuation
+freq0                           =  10            # frequency for viscous attenuation
+p_sv                            = .false.        # set the type of calculation (P-SV or SH/membrane waves)
+
+# time step parameters
+nt                              = 3000           # total number of time steps
+deltat                          = 1.d-3          # duration of a time step
+USER_T0                         = 0.0d0          # use this t0 as earliest starting time rather than the automatically calculated one
+time_stepping_scheme            = 1   # 1 = Newmark (2nd order),     2 = LDDRK4-6 (4th-order 6-stage low storage Runge-Kutta),     3 = classical 4th-order 4-stage Runge-Kutta
+
+# source parameters
+NSOURCES                        = 1              # number of sources [source info read in CMTSOLUTION file]
+force_normal_to_surface         = .false.        # angleforce normal to surface (external mesh and curve file needed)
+
+# constants for attenuation
+N_SLS                           = 2                      # number of standard linear solids for attenuation
+f0_attenuation                  = 5.196152422706633      # (Hz) relevant only if source is a Dirac or a Heaviside, else it is f0
+
+# receiver set parameters for seismograms
+seismotype                      = 1              # record 1=displ 2=veloc 3=accel 4=pressure
+generate_STATIONS               = .true.         # creates a STATION file in ./DATA
+nreceiversets                   = 1              # number of receiver sets
+anglerec                        = 0.d0           # angle to rotate components at receivers
+rec_normal_to_surface           = .false.        # base anglerec normal to surface (external mesh and curve file needed)
+SU_FORMAT                       = .false.        # output seismograms in Seismic Unix format (adjoint traces will be read in the same format)
+
+# first receiver set (repeat these 6 lines and adjust nreceiversets  accordingly)
+nrec                            = 29
+xdeb                            = 400.           # first receiver x in meters
+zdeb                            = 1500.          # first receiver z in meters
+xfin                            = 2100.          # last receiver x in meters (ignored if onlyone receiver)
+zfin                            = 1500.          # last receiver z in meters (ignored if onlyone receiver)
+enreg_surf_same_vertical        = .true.         # receivers inside the medium or at the surface
+
+# display parameters
+NTSTEP_BETWEEN_OUTPUT_INFO      = 250            # display frequency in time steps
+output_postscript_snapshot      = .false.        # output Postscript snapshot of the results
+output_color_image              = .false.        # output color image of the results
+imagetype                       = 1              # display 1=displ 2=veloc 3=accel 4=pressure
+cutsnaps                        = 1.             # minimum amplitude in % for snapshots
+meshvect                        = .true.         # display mesh on vector plots or not
+modelvect                       = .false.        # display velocity model on vector plots
+boundvect                       = .true.         # display boundary conditions on plots
+interpol                        = .true.         # interpolation of the display or not
+pointsdisp                      = 6              # points for interpolation of display (set to 1 for lower-left corner only)
+subsamp_postscript              = 1              # subsampling of color snapshots
+factor_subsample_image          = 1              # factor to subsample color images output by the code (useful for very large models)
+POWER_DISPLAY_COLOR             = 0.30d0         # non linear display to enhance small amplitudes in color images
+DRAW_WATER_CONSTANT_BLUE_IN_JPG = .true.         # display acoustic layers as constant blue in JPEG images, because they likely correspond to water
+sizemax_arrows                  = 1.d0           # maximum size of arrows on vector plots in cm
+US_LETTER                       = .false.        # US letter paper or European A4
+USE_SNAPSHOT_NUMBER_IN_FILENAME = .false.        # use snapshot number in the file name of JPEG color snapshots instead of the time step
+gnuplot                         = .false.        # generate a GNUPLOT file for the grid
+output_grid                     = .false.        # save the grid in a text file or not
+output_energy                   = .false.        # compute and output acoustic and elastic energy (slows down the code significantly)
+output_wavefield_snapshot       = .false.         # output Ux,Uy,Uz text file for each output time (big files)
+
+# velocity and density models
+nbmodels                        = 2              # nb of different models
+# define models as
+# I:   (model_number 1 rho Vp Vs 0 0 Qkappa Qmu 0 0 0 0 0 0) or
+# II:  (model_number 2 rho c11 c13 c15 c33 c35 c55 0 0 0 0 0 0) or
+# III: (model_number 3 rhos rhof phi c kxx kxz kzz Ks Kf Kfr etaf mufr Qmu).
+# For istropic elastic/acoustic material use I and set Vs to zero to make a given model acoustic, for anisotropic elastic use II,
+# and for isotropic poroelastic material use III. The mesh can contain acoustic, elastic, and poroelastic models simultaneously.
+1 1 2700.d0 3000.d0 1732.051d0 0 0 9999 9999 0 0 0 0 0 0
+2 1 2500.d0 2700.d0 1600.000d0 0 0 9999 9999 0 0 0 0 0 0
+
+# external mesh or not
+read_external_mesh              = .false.
+
+# absorbing boundary active or not
+absorbing_conditions            = .true.
+
+# for horizontal periodic conditions: detect common points between left and right edges
+ADD_PERIODIC_CONDITIONS         = .false.
+
+# horizontal periodicity distance for periodic conditions
+PERIODIC_horiz_dist             = 0.3597d0
+
+# grid point detection tolerance for periodic conditions
+PERIODIC_DETECT_TOL             = 3.3334d-6
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR EXTERNAL MESHING
+
+# data concerning mesh, when generated using third-party app (more info in README)
+# (see also absorbing_conditions above)
+mesh_file                       = ./DATA/Mesh_canyon/canyon_mesh_file   # file containing the mesh
+nodes_coords_file               = ./DATA/Mesh_canyon/canyon_nodes_coords_file   # file containing the nodes coordinates
+materials_file                  = ./DATA/Mesh_canyon/canyon_materials_file   # file containing the material number for each element
+free_surface_file               = ./DATA/Mesh_canyon/canyon_free_surface_file   # file containing the free surface
+absorbing_surface_file          = ./DATA/Mesh_canyon/canyon_absorbing_surface_file   # file containing the absorbing surface
+tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the curve delimiting the velocity model
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR INTERNAL MESHING
+
+# file containing interfaces for internal mesh
+interfacesfile                  = interfaces.dat
+
+# geometry of the model (origin lower-left corner = 0,0) and mesh description
+xmin                            = 0.d0           # abscissa of left side of the model
+xmax                            = 2500.d0        # abscissa of right side of the model
+nx                              = 50             # number of elements along X
+
+# absorbing boundary parameters (see absorbing_conditions above)
+absorbbottom                    = .true.
+absorbright                     = .true.
+absorbtop                       = .false.
+absorbleft                      = .true.
+
+# define the different regions of the model in the (nx,nz) spectral element mesh
+nbregions                       = 2              # nb of regions and model number for each
+1 50  1 15 1
+1 50 16 40 2

Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/Par_file_best.before_update_to_r19xxx
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/Par_file_best.before_update_to_r19xxx	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/Par_file_best.before_update_to_r19xxx	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,135 @@
+# title of job
+title                           = PASSIVE IMAGING - 1 INTERFACE
+
+# forward or adjoint simulation
+SIMULATION_TYPE                 = 1
+NOISE_TOMOGRAPHY                = 1
+SAVE_FORWARD                    = .false.
+
+# parameters concerning partitioning
+nproc                           = 1              # number of processes
+partitioning_method             = 3              # SCOTCH = 3, ascending order (very bad idea) = 1
+PERFORM_CUTHILL_MCKEE           = .true.         # perform inverse Cuthill-McKee (1969) optimization/permutation for mesh numbering
+
+ngnod                           = 9              # number of control nodes per element (4 or 9)
+initialfield                    = .false.        # use a plane wave as source or not
+add_Bielak_conditions           = .false.        # add Bielak conditions or not if initial plane wave
+assign_external_model           = .false.        # define external earth model or not
+READ_EXTERNAL_SEP_FILE          = .false.        # Read external SEP file from DATA/model_velocity.dat_input, or use routine
+ATTENUATION_VISCOELASTIC_SOLID  = .false.        # turn attenuation (viscoelasticity) on or off for non-poroelastic solid parts of the model
+ATTENUATION_PORO_FLUID_PART     = .false.        # turn viscous attenuation on or off for the fluid part of poroelastic parts of the model
+Q0                              =  1             # quality factor for viscous attenuation
+freq0                           =  10            # frequency for viscous attenuation
+p_sv                            = .false.        # set the type of calculation (P-SV or SH/membrane waves)
+
+# time step parameters
+nt                              = 3000           # total number of time steps
+deltat                          = 1.d-3          # duration of a time step
+USER_T0                         = 0.0d0          # use this t0 as earliest starting time rather than the automatically calculated one
+
+# source parameters
+NSOURCES                        = 1              # number of sources [source info read in CMTSOLUTION file]
+force_normal_to_surface         = .false.        # angleforce normal to surface (external mesh and curve file needed)
+
+# constants for attenuation
+N_SLS                           = 2                      # number of standard linear solids for attenuation
+f0_attenuation                  = 5.196152422706633      # (Hz) relevant only if source is a Dirac or a Heaviside, else it is f0
+
+# receiver set parameters for seismograms
+seismotype                      = 1              # record 1=displ 2=veloc 3=accel 4=pressure
+generate_STATIONS               = .true.         # creates a STATION file in ./DATA
+nreceiversets                   = 1              # number of receiver sets
+anglerec                        = 0.d0           # angle to rotate components at receivers
+rec_normal_to_surface           = .false.        # base anglerec normal to surface (external mesh and curve file needed)
+SU_FORMAT                       = .false.        # output seismograms in Seismic Unix format (adjoint traces will be read in the same format)
+
+# first receiver set (repeat these 6 lines and adjust nreceiversets  accordingly)
+nrec                            = 29
+xdeb                            = 400.           # first receiver x in meters
+zdeb                            = 1500.          # first receiver z in meters
+xfin                            = 2100.          # last receiver x in meters (ignored if onlyone receiver)
+zfin                            = 1500.          # last receiver z in meters (ignored if onlyone receiver)
+enreg_surf_same_vertical        = .true.         # receivers inside the medium or at the surface
+
+# display parameters
+NTSTEP_BETWEEN_OUTPUT_INFO      = 250            # display frequency in time steps
+output_postscript_snapshot      = .false.        # output Postscript snapshot of the results
+output_color_image              = .false.        # output color image of the results
+imagetype                       = 1              # display 1=displ 2=veloc 3=accel 4=pressure
+cutsnaps                        = 1.             # minimum amplitude in % for snapshots
+meshvect                        = .true.         # display mesh on vector plots or not
+modelvect                       = .false.        # display velocity model on vector plots
+boundvect                       = .true.         # display boundary conditions on plots
+interpol                        = .true.         # interpolation of the display or not
+pointsdisp                      = 6              # points for interpolation of display (set to 1 for lower-left corner only)
+subsamp_postscript              = 1              # subsampling of color snapshots
+factor_subsample_image          = 1              # factor to subsample color images output by the code (useful for very large models)
+POWER_DISPLAY_COLOR             = 0.30d0         # non linear display to enhance small amplitudes in color images
+DRAW_WATER_CONSTANT_BLUE_IN_JPG = .true.         # display acoustic layers as constant blue in JPEG images, because they likely correspond to water
+sizemax_arrows                  = 1.d0           # maximum size of arrows on vector plots in cm
+US_LETTER                       = .false.        # US letter paper or European A4
+USE_SNAPSHOT_NUMBER_IN_FILENAME = .false.        # use snapshot number in the file name of JPEG color snapshots instead of the time step
+gnuplot                         = .false.        # generate a GNUPLOT file for the grid
+output_grid                     = .false.        # save the grid in a text file or not
+output_energy                   = .false.        # compute and output acoustic and elastic energy (slows down the code significantly)
+output_wavefield_snapshot       = .false.         # output Ux,Uy,Uz text file for each output time (big files)
+
+# velocity and density models
+nbmodels                        = 2              # nb of different models
+# define models as
+# I:   (model_number 1 rho Vp Vs 0 0 Qkappa Qmu 0 0 0 0 0 0) or
+# II:  (model_number 2 rho c11 c13 c15 c33 c35 c55 0 0 0 0 0 0) or
+# III: (model_number 3 rhos rhof phi c kxx kxz kzz Ks Kf Kfr etaf mufr Qmu).
+# For istropic elastic/acoustic material use I and set Vs to zero to make a given model acoustic, for anisotropic elastic use II,
+# and for isotropic poroelastic material use III. The mesh can contain acoustic, elastic, and poroelastic models simultaneously.
+1 1 2700.d0 3000.d0 1732.051d0 0 0 9999 9999 0 0 0 0 0 0
+2 1 2500.d0 2700.d0 1600.000d0 0 0 9999 9999 0 0 0 0 0 0
+
+# external mesh or not
+read_external_mesh              = .false.
+
+# absorbing boundary active or not
+absorbing_conditions            = .true.
+
+# for horizontal periodic conditions: detect common points between left and right edges
+ADD_PERIODIC_CONDITIONS         = .false.
+
+# horizontal periodicity distance for periodic conditions
+PERIODIC_horiz_dist             = 0.3597d0
+
+# grid point detection tolerance for periodic conditions
+PERIODIC_DETECT_TOL             = 3.3334d-6
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR EXTERNAL MESHING
+
+# data concerning mesh, when generated using third-party app (more info in README)
+# (see also absorbing_conditions above)
+mesh_file                       = ./DATA/Mesh_canyon/canyon_mesh_file   # file containing the mesh
+nodes_coords_file               = ./DATA/Mesh_canyon/canyon_nodes_coords_file   # file containing the nodes coordinates
+materials_file                  = ./DATA/Mesh_canyon/canyon_materials_file   # file containing the material number for each element
+free_surface_file               = ./DATA/Mesh_canyon/canyon_free_surface_file   # file containing the free surface
+absorbing_surface_file          = ./DATA/Mesh_canyon/canyon_absorbing_surface_file   # file containing the absorbing surface
+tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the curve delimiting the velocity model
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR INTERNAL MESHING
+
+# file containing interfaces for internal mesh
+interfacesfile                  = interfaces.dat
+
+# geometry of the model (origin lower-left corner = 0,0) and mesh description
+xmin                            = 0.d0           # abscissa of left side of the model
+xmax                            = 2500.d0        # abscissa of right side of the model
+nx                              = 50             # number of elements along X
+
+# absorbing boundary parameters (see absorbing_conditions above)
+absorbbottom                    = .true.
+absorbright                     = .true.
+absorbtop                       = .false.
+absorbleft                      = .true.
+
+# define the different regions of the model in the (nx,nz) spectral element mesh
+nbregions                       = 2              # nb of regions and model number for each
+1 50  1 15 1
+1 50 16 40 2

Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/Par_file_fair
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/Par_file_fair	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/Par_file_fair	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,136 @@
+# title of job
+title                           = PASSIVE IMAGING - 1 INTERFACE
+
+# forward or adjoint simulation
+SIMULATION_TYPE                 = 1
+NOISE_TOMOGRAPHY                = 1
+SAVE_FORWARD                    = .false.
+
+# parameters concerning partitioning
+nproc                           = 1              # number of processes
+partitioning_method             = 3              # SCOTCH = 3, ascending order (very bad idea) = 1
+PERFORM_CUTHILL_MCKEE           = .true.         # perform inverse Cuthill-McKee (1969) optimization/permutation for mesh numbering
+
+ngnod                           = 9              # number of control nodes per element (4 or 9)
+initialfield                    = .false.        # use a plane wave as source or not
+add_Bielak_conditions           = .false.        # add Bielak conditions or not if initial plane wave
+assign_external_model           = .false.        # define external earth model or not
+READ_EXTERNAL_SEP_FILE          = .false.        # Read external SEP file from DATA/model_velocity.dat_input, or use routine
+ATTENUATION_VISCOELASTIC_SOLID  = .false.        # turn attenuation (viscoelasticity) on or off for non-poroelastic solid parts of the model
+ATTENUATION_PORO_FLUID_PART     = .false.        # turn viscous attenuation on or off for the fluid part of poroelastic parts of the model
+Q0                              =  1             # quality factor for viscous attenuation
+freq0                           =  10            # frequency for viscous attenuation
+p_sv                            = .false.        # set the type of calculation (P-SV or SH/membrane waves)
+
+# time step parameters
+nt                              = 3000           # total number of time steps
+deltat                          = 1.d-3          # duration of a time step
+USER_T0                         = 0.0d0          # use this t0 as earliest starting time rather than the automatically calculated one
+time_stepping_scheme            = 1   # 1 = Newmark (2nd order),     2 = LDDRK4-6 (4th-order 6-stage low storage Runge-Kutta),     3 = classical 4th-order 4-stage Runge-Kutta
+
+# source parameters
+NSOURCES                        = 1              # number of sources [source info read in CMTSOLUTION file]
+force_normal_to_surface         = .false.        # angleforce normal to surface (external mesh and curve file needed)
+
+# constants for attenuation
+N_SLS                           = 2                      # number of standard linear solids for attenuation
+f0_attenuation                  = 5.196152422706633      # (Hz) relevant only if source is a Dirac or a Heaviside, else it is f0
+
+# receiver set parameters for seismograms
+seismotype                      = 1              # record 1=displ 2=veloc 3=accel 4=pressure
+generate_STATIONS               = .true.         # creates a STATION file in ./DATA
+nreceiversets                   = 1              # number of receiver sets
+anglerec                        = 0.d0           # angle to rotate components at receivers
+rec_normal_to_surface           = .false.        # base anglerec normal to surface (external mesh and curve file needed)
+SU_FORMAT                       = .false.        # output seismograms in Seismic Unix format (adjoint traces will be read in the same format)
+
+# first receiver set (repeat these 6 lines and adjust nreceiversets  accordingly)
+nrec                            = 29
+xdeb                            = 150.           # first receiver x in meters
+zdeb                            = 1500.          # first receiver z in meters
+xfin                            = 1850.          # last receiver x in meters (ignored if onlyone receiver)
+zfin                            = 1500.          # last receiver z in meters (ignored if onlyone receiver)
+enreg_surf_same_vertical        = .true.         # receivers inside the medium or at the surface
+
+# display parameters
+NTSTEP_BETWEEN_OUTPUT_INFO      = 250            # display frequency in time steps
+output_postscript_snapshot      = .false.        # output Postscript snapshot of the results
+output_color_image              = .false.        # output color image of the results
+imagetype                       = 1              # display 1=displ 2=veloc 3=accel 4=pressure
+cutsnaps                        = 1.             # minimum amplitude in % for snapshots
+meshvect                        = .true.         # display mesh on vector plots or not
+modelvect                       = .false.        # display velocity model on vector plots
+boundvect                       = .true.         # display boundary conditions on plots
+interpol                        = .true.         # interpolation of the display or not
+pointsdisp                      = 6              # points for interpolation of display (set to 1 for lower-left corner only)
+subsamp_postscript              = 1              # subsampling of color snapshots
+factor_subsample_image          = 1              # factor to subsample color images output by the code (useful for very large models)
+POWER_DISPLAY_COLOR             = 0.30d0         # non linear display to enhance small amplitudes in color images
+DRAW_WATER_CONSTANT_BLUE_IN_JPG = .true.         # display acoustic layers as constant blue in JPEG images, because they likely correspond to water
+sizemax_arrows                  = 1.d0           # maximum size of arrows on vector plots in cm
+US_LETTER                       = .false.        # US letter paper or European A4
+USE_SNAPSHOT_NUMBER_IN_FILENAME = .false.        # use snapshot number in the file name of JPEG color snapshots instead of the time step
+gnuplot                         = .false.        # generate a GNUPLOT file for the grid
+output_grid                     = .false.        # save the grid in a text file or not
+output_energy                   = .false.        # compute and output acoustic and elastic energy (slows down the code significantly)
+output_wavefield_snapshot       = .false.         # output Ux,Uy,Uz text file for each output time (big files)
+
+# velocity and density models
+nbmodels                        = 2              # nb of different models
+# define models as
+# I:   (model_number 1 rho Vp Vs 0 0 Qkappa Qmu 0 0 0 0 0 0) or
+# II:  (model_number 2 rho c11 c13 c15 c33 c35 c55 0 0 0 0 0 0) or
+# III: (model_number 3 rhos rhof phi c kxx kxz kzz Ks Kf Kfr etaf mufr Qmu).
+# For istropic elastic/acoustic material use I and set Vs to zero to make a given model acoustic, for anisotropic elastic use II,
+# and for isotropic poroelastic material use III. The mesh can contain acoustic, elastic, and poroelastic models simultaneously.
+1 1 2700.d0 3000.d0 1732.051d0 0 0 9999 9999 0 0 0 0 0 0
+2 1 2500.d0 2700.d0 1600.000d0 0 0 9999 9999 0 0 0 0 0 0
+
+# external mesh or not
+read_external_mesh              = .false.
+
+# absorbing boundary active or not
+absorbing_conditions            = .true.
+
+# for horizontal periodic conditions: detect common points between left and right edges
+ADD_PERIODIC_CONDITIONS         = .false.
+
+# horizontal periodicity distance for periodic conditions
+PERIODIC_horiz_dist             = 0.3597d0
+
+# grid point detection tolerance for periodic conditions
+PERIODIC_DETECT_TOL             = 3.3334d-6
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR EXTERNAL MESHING
+
+# data concerning mesh, when generated using third-party app (more info in README)
+# (see also absorbing_conditions above)
+mesh_file                       = ./DATA/Mesh_canyon/canyon_mesh_file   # file containing the mesh
+nodes_coords_file               = ./DATA/Mesh_canyon/canyon_nodes_coords_file   # file containing the nodes coordinates
+materials_file                  = ./DATA/Mesh_canyon/canyon_materials_file   # file containing the material number for each element
+free_surface_file               = ./DATA/Mesh_canyon/canyon_free_surface_file   # file containing the free surface
+absorbing_surface_file          = ./DATA/Mesh_canyon/canyon_absorbing_surface_file   # file containing the absorbing surface
+tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the curve delimiting the velocity model
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR INTERNAL MESHING
+
+# file containing interfaces for internal mesh
+interfacesfile                  = interfaces.dat
+
+# geometry of the model (origin lower-left corner = 0,0) and mesh description
+xmin                            = 0.d0           # abscissa of left side of the model
+xmax                            = 2000.d0        # abscissa of right side of the model
+nx                              = 40             # number of elements along X
+
+# absorbing boundary parameters (see absorbing_conditions above)
+absorbbottom                    = .true.
+absorbright                     = .true.
+absorbtop                       = .false.
+absorbleft                      = .true.
+
+# define the different regions of the model in the (nx,nz) spectral element mesh
+nbregions                       = 2              # nb of regions and model number for each
+1 40  1 15 1
+1 40 16 30 2

Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/Par_file_fair.before_update_to_r19xxx
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/Par_file_fair.before_update_to_r19xxx	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/Par_file_fair.before_update_to_r19xxx	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,135 @@
+# title of job
+title                           = PASSIVE IMAGING - 1 INTERFACE
+
+# forward or adjoint simulation
+SIMULATION_TYPE                 = 1
+NOISE_TOMOGRAPHY                = 1
+SAVE_FORWARD                    = .false.
+
+# parameters concerning partitioning
+nproc                           = 1              # number of processes
+partitioning_method             = 3              # SCOTCH = 3, ascending order (very bad idea) = 1
+PERFORM_CUTHILL_MCKEE           = .true.         # perform inverse Cuthill-McKee (1969) optimization/permutation for mesh numbering
+
+ngnod                           = 9              # number of control nodes per element (4 or 9)
+initialfield                    = .false.        # use a plane wave as source or not
+add_Bielak_conditions           = .false.        # add Bielak conditions or not if initial plane wave
+assign_external_model           = .false.        # define external earth model or not
+READ_EXTERNAL_SEP_FILE          = .false.        # Read external SEP file from DATA/model_velocity.dat_input, or use routine
+ATTENUATION_VISCOELASTIC_SOLID  = .false.        # turn attenuation (viscoelasticity) on or off for non-poroelastic solid parts of the model
+ATTENUATION_PORO_FLUID_PART     = .false.        # turn viscous attenuation on or off for the fluid part of poroelastic parts of the model
+Q0                              =  1             # quality factor for viscous attenuation
+freq0                           =  10            # frequency for viscous attenuation
+p_sv                            = .false.        # set the type of calculation (P-SV or SH/membrane waves)
+
+# time step parameters
+nt                              = 3000           # total number of time steps
+deltat                          = 1.d-3          # duration of a time step
+USER_T0                         = 0.0d0          # use this t0 as earliest starting time rather than the automatically calculated one
+
+# source parameters
+NSOURCES                        = 1              # number of sources [source info read in CMTSOLUTION file]
+force_normal_to_surface         = .false.        # angleforce normal to surface (external mesh and curve file needed)
+
+# constants for attenuation
+N_SLS                           = 2                      # number of standard linear solids for attenuation
+f0_attenuation                  = 5.196152422706633      # (Hz) relevant only if source is a Dirac or a Heaviside, else it is f0
+
+# receiver set parameters for seismograms
+seismotype                      = 1              # record 1=displ 2=veloc 3=accel 4=pressure
+generate_STATIONS               = .true.         # creates a STATION file in ./DATA
+nreceiversets                   = 1              # number of receiver sets
+anglerec                        = 0.d0           # angle to rotate components at receivers
+rec_normal_to_surface           = .false.        # base anglerec normal to surface (external mesh and curve file needed)
+SU_FORMAT                       = .false.        # output seismograms in Seismic Unix format (adjoint traces will be read in the same format)
+
+# first receiver set (repeat these 6 lines and adjust nreceiversets  accordingly)
+nrec                            = 29
+xdeb                            = 150.           # first receiver x in meters
+zdeb                            = 1500.          # first receiver z in meters
+xfin                            = 1850.          # last receiver x in meters (ignored if onlyone receiver)
+zfin                            = 1500.          # last receiver z in meters (ignored if onlyone receiver)
+enreg_surf_same_vertical        = .true.         # receivers inside the medium or at the surface
+
+# display parameters
+NTSTEP_BETWEEN_OUTPUT_INFO      = 250            # display frequency in time steps
+output_postscript_snapshot      = .false.        # output Postscript snapshot of the results
+output_color_image              = .false.        # output color image of the results
+imagetype                       = 1              # display 1=displ 2=veloc 3=accel 4=pressure
+cutsnaps                        = 1.             # minimum amplitude in % for snapshots
+meshvect                        = .true.         # display mesh on vector plots or not
+modelvect                       = .false.        # display velocity model on vector plots
+boundvect                       = .true.         # display boundary conditions on plots
+interpol                        = .true.         # interpolation of the display or not
+pointsdisp                      = 6              # points for interpolation of display (set to 1 for lower-left corner only)
+subsamp_postscript              = 1              # subsampling of color snapshots
+factor_subsample_image          = 1              # factor to subsample color images output by the code (useful for very large models)
+POWER_DISPLAY_COLOR             = 0.30d0         # non linear display to enhance small amplitudes in color images
+DRAW_WATER_CONSTANT_BLUE_IN_JPG = .true.         # display acoustic layers as constant blue in JPEG images, because they likely correspond to water
+sizemax_arrows                  = 1.d0           # maximum size of arrows on vector plots in cm
+US_LETTER                       = .false.        # US letter paper or European A4
+USE_SNAPSHOT_NUMBER_IN_FILENAME = .false.        # use snapshot number in the file name of JPEG color snapshots instead of the time step
+gnuplot                         = .false.        # generate a GNUPLOT file for the grid
+output_grid                     = .false.        # save the grid in a text file or not
+output_energy                   = .false.        # compute and output acoustic and elastic energy (slows down the code significantly)
+output_wavefield_snapshot       = .false.         # output Ux,Uy,Uz text file for each output time (big files)
+
+# velocity and density models
+nbmodels                        = 2              # nb of different models
+# define models as
+# I:   (model_number 1 rho Vp Vs 0 0 Qkappa Qmu 0 0 0 0 0 0) or
+# II:  (model_number 2 rho c11 c13 c15 c33 c35 c55 0 0 0 0 0 0) or
+# III: (model_number 3 rhos rhof phi c kxx kxz kzz Ks Kf Kfr etaf mufr Qmu).
+# For istropic elastic/acoustic material use I and set Vs to zero to make a given model acoustic, for anisotropic elastic use II,
+# and for isotropic poroelastic material use III. The mesh can contain acoustic, elastic, and poroelastic models simultaneously.
+1 1 2700.d0 3000.d0 1732.051d0 0 0 9999 9999 0 0 0 0 0 0
+2 1 2500.d0 2700.d0 1600.000d0 0 0 9999 9999 0 0 0 0 0 0
+
+# external mesh or not
+read_external_mesh              = .false.
+
+# absorbing boundary active or not
+absorbing_conditions            = .true.
+
+# for horizontal periodic conditions: detect common points between left and right edges
+ADD_PERIODIC_CONDITIONS         = .false.
+
+# horizontal periodicity distance for periodic conditions
+PERIODIC_horiz_dist             = 0.3597d0
+
+# grid point detection tolerance for periodic conditions
+PERIODIC_DETECT_TOL             = 3.3334d-6
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR EXTERNAL MESHING
+
+# data concerning mesh, when generated using third-party app (more info in README)
+# (see also absorbing_conditions above)
+mesh_file                       = ./DATA/Mesh_canyon/canyon_mesh_file   # file containing the mesh
+nodes_coords_file               = ./DATA/Mesh_canyon/canyon_nodes_coords_file   # file containing the nodes coordinates
+materials_file                  = ./DATA/Mesh_canyon/canyon_materials_file   # file containing the material number for each element
+free_surface_file               = ./DATA/Mesh_canyon/canyon_free_surface_file   # file containing the free surface
+absorbing_surface_file          = ./DATA/Mesh_canyon/canyon_absorbing_surface_file   # file containing the absorbing surface
+tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the curve delimiting the velocity model
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR INTERNAL MESHING
+
+# file containing interfaces for internal mesh
+interfacesfile                  = interfaces.dat
+
+# geometry of the model (origin lower-left corner = 0,0) and mesh description
+xmin                            = 0.d0           # abscissa of left side of the model
+xmax                            = 2000.d0        # abscissa of right side of the model
+nx                              = 40             # number of elements along X
+
+# absorbing boundary parameters (see absorbing_conditions above)
+absorbbottom                    = .true.
+absorbright                     = .true.
+absorbtop                       = .false.
+absorbleft                      = .true.
+
+# define the different regions of the model in the (nx,nz) spectral element mesh
+nbregions                       = 2              # nb of regions and model number for each
+1 40  1 15 1
+1 40 16 30 2

Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/Par_file_good
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/Par_file_good	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/Par_file_good	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,136 @@
+# title of job
+title                           = PASSIVE IMAGING - 1 INTERFACE
+
+# forward or adjoint simulation
+SIMULATION_TYPE                 = 1
+NOISE_TOMOGRAPHY                = 1
+SAVE_FORWARD                    = .false.
+
+# parameters concerning partitioning
+nproc                           = 1              # number of processes
+partitioning_method             = 3              # SCOTCH = 3, ascending order (very bad idea) = 1
+PERFORM_CUTHILL_MCKEE           = .true.         # perform inverse Cuthill-McKee (1969) optimization/permutation for mesh numbering
+
+ngnod                           = 9              # number of control nodes per element (4 or 9)
+initialfield                    = .false.        # use a plane wave as source or not
+add_Bielak_conditions           = .false.        # add Bielak conditions or not if initial plane wave
+assign_external_model           = .false.        # define external earth model or not
+READ_EXTERNAL_SEP_FILE          = .false.        # Read external SEP file from DATA/model_velocity.dat_input, or use routine
+ATTENUATION_VISCOELASTIC_SOLID  = .false.        # turn attenuation (viscoelasticity) on or off for non-poroelastic solid parts of the model
+ATTENUATION_PORO_FLUID_PART     = .false.        # turn viscous attenuation on or off for the fluid part of poroelastic parts of the model
+Q0                              =  1             # quality factor for viscous attenuation
+freq0                           =  10            # frequency for viscous attenuation
+p_sv                            = .false.        # set the type of calculation (P-SV or SH/membrane waves)
+
+# time step parameters
+nt                              = 3000           # total number of time steps
+deltat                          = 1.d-3          # duration of a time step
+USER_T0                         = 0.0d0          # use this t0 as earliest starting time rather than the automatically calculated one
+time_stepping_scheme            = 1   # 1 = Newmark (2nd order),     2 = LDDRK4-6 (4th-order 6-stage low storage Runge-Kutta),     3 = classical 4th-order 4-stage Runge-Kutta
+
+# source parameters
+NSOURCES                        = 1              # number of sources [source info read in CMTSOLUTION file]
+force_normal_to_surface         = .false.        # angleforce normal to surface (external mesh and curve file needed)
+
+# constants for attenuation
+N_SLS                           = 2                      # number of standard linear solids for attenuation
+f0_attenuation                  = 5.196152422706633      # (Hz) relevant only if source is a Dirac or a Heaviside, else it is f0
+
+# receiver set parameters for seismograms
+seismotype                      = 1              # record 1=displ 2=veloc 3=accel 4=pressure
+generate_STATIONS               = .true.         # creates a STATION file in ./DATA
+nreceiversets                   = 1              # number of receiver sets
+anglerec                        = 0.d0           # angle to rotate components at receivers
+rec_normal_to_surface           = .false.        # base anglerec normal to surface (external mesh and curve file needed)
+SU_FORMAT                       = .false.        # output seismograms in Seismic Unix format (adjoint traces will be read in the same format)
+
+# first receiver set (repeat these 6 lines and adjust nreceiversets  accordingly)
+nrec                            = 29
+xdeb                            = 400.           # first receiver x in meters
+zdeb                            = 1500.          # first receiver z in meters
+xfin                            = 2100.          # last receiver x in meters (ignored if onlyone receiver)
+zfin                            = 1500.          # last receiver z in meters (ignored if onlyone receiver)
+enreg_surf_same_vertical        = .true.         # receivers inside the medium or at the surface
+
+# display parameters
+NTSTEP_BETWEEN_OUTPUT_INFO      = 250            # display frequency in time steps
+output_postscript_snapshot      = .false.        # output Postscript snapshot of the results
+output_color_image              = .false.        # output color image of the results
+imagetype                       = 1              # display 1=displ 2=veloc 3=accel 4=pressure
+cutsnaps                        = 1.             # minimum amplitude in % for snapshots
+meshvect                        = .true.         # display mesh on vector plots or not
+modelvect                       = .false.        # display velocity model on vector plots
+boundvect                       = .true.         # display boundary conditions on plots
+interpol                        = .true.         # interpolation of the display or not
+pointsdisp                      = 6              # points for interpolation of display (set to 1 for lower-left corner only)
+subsamp_postscript              = 1              # subsampling of color snapshots
+factor_subsample_image          = 1              # factor to subsample color images output by the code (useful for very large models)
+POWER_DISPLAY_COLOR             = 0.30d0         # non linear display to enhance small amplitudes in color images
+DRAW_WATER_CONSTANT_BLUE_IN_JPG = .true.         # display acoustic layers as constant blue in JPEG images, because they likely correspond to water
+sizemax_arrows                  = 1.d0           # maximum size of arrows on vector plots in cm
+US_LETTER                       = .false.        # US letter paper or European A4
+USE_SNAPSHOT_NUMBER_IN_FILENAME = .false.        # use snapshot number in the file name of JPEG color snapshots instead of the time step
+gnuplot                         = .false.        # generate a GNUPLOT file for the grid
+output_grid                     = .false.        # save the grid in a text file or not
+output_energy                   = .false.        # compute and output acoustic and elastic energy (slows down the code significantly)
+output_wavefield_snapshot       = .false.         # output Ux,Uy,Uz text file for each output time (big files)
+
+# velocity and density models
+nbmodels                        = 2              # nb of different models
+# define models as
+# I:   (model_number 1 rho Vp Vs 0 0 Qkappa Qmu 0 0 0 0 0 0) or
+# II:  (model_number 2 rho c11 c13 c15 c33 c35 c55 0 0 0 0 0 0) or
+# III: (model_number 3 rhos rhof phi c kxx kxz kzz Ks Kf Kfr etaf mufr Qmu).
+# For istropic elastic/acoustic material use I and set Vs to zero to make a given model acoustic, for anisotropic elastic use II,
+# and for isotropic poroelastic material use III. The mesh can contain acoustic, elastic, and poroelastic models simultaneously.
+1 1 2700.d0 3000.d0 1732.051d0 0 0 9999 9999 0 0 0 0 0 0
+2 1 2500.d0 2700.d0 1600.000d0 0 0 9999 9999 0 0 0 0 0 0
+
+# external mesh or not
+read_external_mesh              = .false.
+
+# absorbing boundary active or not
+absorbing_conditions            = .true.
+
+# for horizontal periodic conditions: detect common points between left and right edges
+ADD_PERIODIC_CONDITIONS         = .false.
+
+# horizontal periodicity distance for periodic conditions
+PERIODIC_horiz_dist             = 0.3597d0
+
+# grid point detection tolerance for periodic conditions
+PERIODIC_DETECT_TOL             = 3.3334d-6
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR EXTERNAL MESHING
+
+# data concerning mesh, when generated using third-party app (more info in README)
+# (see also absorbing_conditions above)
+mesh_file                       = ./DATA/Mesh_canyon/canyon_mesh_file   # file containing the mesh
+nodes_coords_file               = ./DATA/Mesh_canyon/canyon_nodes_coords_file   # file containing the nodes coordinates
+materials_file                  = ./DATA/Mesh_canyon/canyon_materials_file   # file containing the material number for each element
+free_surface_file               = ./DATA/Mesh_canyon/canyon_free_surface_file   # file containing the free surface
+absorbing_surface_file          = ./DATA/Mesh_canyon/canyon_absorbing_surface_file   # file containing the absorbing surface
+tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the curve delimiting the velocity model
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR INTERNAL MESHING
+
+# file containing interfaces for internal mesh
+interfacesfile                  = interfaces.dat
+
+# geometry of the model (origin lower-left corner = 0,0) and mesh description
+xmin                            = 0.d0           # abscissa of left side of the model
+xmax                            = 2500.d0        # abscissa of right side of the model
+nx                              = 50             # number of elements along X
+
+# absorbing boundary parameters (see absorbing_conditions above)
+absorbbottom                    = .true.
+absorbright                     = .true.
+absorbtop                       = .false.
+absorbleft                      = .true.
+
+# define the different regions of the model in the (nx,nz) spectral element mesh
+nbregions                       = 2              # nb of regions and model number for each
+1 50  1 15 1
+1 50 16 30 2

Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/Par_file_good.before_update_to_r19xxx
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/Par_file_good.before_update_to_r19xxx	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/Par_file_good.before_update_to_r19xxx	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,135 @@
+# title of job
+title                           = PASSIVE IMAGING - 1 INTERFACE
+
+# forward or adjoint simulation
+SIMULATION_TYPE                 = 1
+NOISE_TOMOGRAPHY                = 1
+SAVE_FORWARD                    = .false.
+
+# parameters concerning partitioning
+nproc                           = 1              # number of processes
+partitioning_method             = 3              # SCOTCH = 3, ascending order (very bad idea) = 1
+PERFORM_CUTHILL_MCKEE           = .true.         # perform inverse Cuthill-McKee (1969) optimization/permutation for mesh numbering
+
+ngnod                           = 9              # number of control nodes per element (4 or 9)
+initialfield                    = .false.        # use a plane wave as source or not
+add_Bielak_conditions           = .false.        # add Bielak conditions or not if initial plane wave
+assign_external_model           = .false.        # define external earth model or not
+READ_EXTERNAL_SEP_FILE          = .false.        # Read external SEP file from DATA/model_velocity.dat_input, or use routine
+ATTENUATION_VISCOELASTIC_SOLID  = .false.        # turn attenuation (viscoelasticity) on or off for non-poroelastic solid parts of the model
+ATTENUATION_PORO_FLUID_PART     = .false.        # turn viscous attenuation on or off for the fluid part of poroelastic parts of the model
+Q0                              =  1             # quality factor for viscous attenuation
+freq0                           =  10            # frequency for viscous attenuation
+p_sv                            = .false.        # set the type of calculation (P-SV or SH/membrane waves)
+
+# time step parameters
+nt                              = 3000           # total number of time steps
+deltat                          = 1.d-3          # duration of a time step
+USER_T0                         = 0.0d0          # use this t0 as earliest starting time rather than the automatically calculated one
+
+# source parameters
+NSOURCES                        = 1              # number of sources [source info read in CMTSOLUTION file]
+force_normal_to_surface         = .false.        # angleforce normal to surface (external mesh and curve file needed)
+
+# constants for attenuation
+N_SLS                           = 2                      # number of standard linear solids for attenuation
+f0_attenuation                  = 5.196152422706633      # (Hz) relevant only if source is a Dirac or a Heaviside, else it is f0
+
+# receiver set parameters for seismograms
+seismotype                      = 1              # record 1=displ 2=veloc 3=accel 4=pressure
+generate_STATIONS               = .true.         # creates a STATION file in ./DATA
+nreceiversets                   = 1              # number of receiver sets
+anglerec                        = 0.d0           # angle to rotate components at receivers
+rec_normal_to_surface           = .false.        # base anglerec normal to surface (external mesh and curve file needed)
+SU_FORMAT                       = .false.        # output seismograms in Seismic Unix format (adjoint traces will be read in the same format)
+
+# first receiver set (repeat these 6 lines and adjust nreceiversets  accordingly)
+nrec                            = 29
+xdeb                            = 400.           # first receiver x in meters
+zdeb                            = 1500.          # first receiver z in meters
+xfin                            = 2100.          # last receiver x in meters (ignored if onlyone receiver)
+zfin                            = 1500.          # last receiver z in meters (ignored if onlyone receiver)
+enreg_surf_same_vertical        = .true.         # receivers inside the medium or at the surface
+
+# display parameters
+NTSTEP_BETWEEN_OUTPUT_INFO      = 250            # display frequency in time steps
+output_postscript_snapshot      = .false.        # output Postscript snapshot of the results
+output_color_image              = .false.        # output color image of the results
+imagetype                       = 1              # display 1=displ 2=veloc 3=accel 4=pressure
+cutsnaps                        = 1.             # minimum amplitude in % for snapshots
+meshvect                        = .true.         # display mesh on vector plots or not
+modelvect                       = .false.        # display velocity model on vector plots
+boundvect                       = .true.         # display boundary conditions on plots
+interpol                        = .true.         # interpolation of the display or not
+pointsdisp                      = 6              # points for interpolation of display (set to 1 for lower-left corner only)
+subsamp_postscript              = 1              # subsampling of color snapshots
+factor_subsample_image          = 1              # factor to subsample color images output by the code (useful for very large models)
+POWER_DISPLAY_COLOR             = 0.30d0         # non linear display to enhance small amplitudes in color images
+DRAW_WATER_CONSTANT_BLUE_IN_JPG = .true.         # display acoustic layers as constant blue in JPEG images, because they likely correspond to water
+sizemax_arrows                  = 1.d0           # maximum size of arrows on vector plots in cm
+US_LETTER                       = .false.        # US letter paper or European A4
+USE_SNAPSHOT_NUMBER_IN_FILENAME = .false.        # use snapshot number in the file name of JPEG color snapshots instead of the time step
+gnuplot                         = .false.        # generate a GNUPLOT file for the grid
+output_grid                     = .false.        # save the grid in a text file or not
+output_energy                   = .false.        # compute and output acoustic and elastic energy (slows down the code significantly)
+output_wavefield_snapshot       = .false.         # output Ux,Uy,Uz text file for each output time (big files)
+
+# velocity and density models
+nbmodels                        = 2              # nb of different models
+# define models as
+# I:   (model_number 1 rho Vp Vs 0 0 Qkappa Qmu 0 0 0 0 0 0) or
+# II:  (model_number 2 rho c11 c13 c15 c33 c35 c55 0 0 0 0 0 0) or
+# III: (model_number 3 rhos rhof phi c kxx kxz kzz Ks Kf Kfr etaf mufr Qmu).
+# For istropic elastic/acoustic material use I and set Vs to zero to make a given model acoustic, for anisotropic elastic use II,
+# and for isotropic poroelastic material use III. The mesh can contain acoustic, elastic, and poroelastic models simultaneously.
+1 1 2700.d0 3000.d0 1732.051d0 0 0 9999 9999 0 0 0 0 0 0
+2 1 2500.d0 2700.d0 1600.000d0 0 0 9999 9999 0 0 0 0 0 0
+
+# external mesh or not
+read_external_mesh              = .false.
+
+# absorbing boundary active or not
+absorbing_conditions            = .true.
+
+# for horizontal periodic conditions: detect common points between left and right edges
+ADD_PERIODIC_CONDITIONS         = .false.
+
+# horizontal periodicity distance for periodic conditions
+PERIODIC_horiz_dist             = 0.3597d0
+
+# grid point detection tolerance for periodic conditions
+PERIODIC_DETECT_TOL             = 3.3334d-6
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR EXTERNAL MESHING
+
+# data concerning mesh, when generated using third-party app (more info in README)
+# (see also absorbing_conditions above)
+mesh_file                       = ./DATA/Mesh_canyon/canyon_mesh_file   # file containing the mesh
+nodes_coords_file               = ./DATA/Mesh_canyon/canyon_nodes_coords_file   # file containing the nodes coordinates
+materials_file                  = ./DATA/Mesh_canyon/canyon_materials_file   # file containing the material number for each element
+free_surface_file               = ./DATA/Mesh_canyon/canyon_free_surface_file   # file containing the free surface
+absorbing_surface_file          = ./DATA/Mesh_canyon/canyon_absorbing_surface_file   # file containing the absorbing surface
+tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the curve delimiting the velocity model
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR INTERNAL MESHING
+
+# file containing interfaces for internal mesh
+interfacesfile                  = interfaces.dat
+
+# geometry of the model (origin lower-left corner = 0,0) and mesh description
+xmin                            = 0.d0           # abscissa of left side of the model
+xmax                            = 2500.d0        # abscissa of right side of the model
+nx                              = 50             # number of elements along X
+
+# absorbing boundary parameters (see absorbing_conditions above)
+absorbbottom                    = .true.
+absorbright                     = .true.
+absorbtop                       = .false.
+absorbleft                      = .true.
+
+# define the different regions of the model in the (nx,nz) spectral element mesh
+nbregions                       = 2              # nb of regions and model number for each
+1 50  1 15 1
+1 50 16 30 2

Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/SOURCE_noise
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/SOURCE_noise	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/SOURCE_noise	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,13 @@
+#source 1.  The components of a moment tensor source must be given in N.m, not in dyne.cm as in the DATA/CMTSOLUTION source file of the 3D version of the code.
+source_surf                     = .false.        # source inside the medium or at the surface
+xs                              = 0.
+zs                              = 0.
+source_type                     = 1              # elastic force or acoustic pressure = 1 or moment tensor = 2
+time_function_type              = 3              # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4, Heaviside = 5
+f0                              = 10.0           # dominant source frequency (Hz) if not Dirac or Heaviside
+tshift                          = 0.0            # time shift when multi sources (if one source, must be zero)
+angleforce                      = 0.             # angle of the source (for a force only)
+Mxx                             = 0.d0           # Mxx component (for a moment tensor source only)
+Mzz                             = 0.d0           # Mzz component (for a moment tensor source only)
+Mxz                             = 0.d0           # Mxz component (for a moment tensor source only)
+factor                          = 0.d0          # amplification factor

Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/interfaces_best
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/interfaces_best	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/interfaces_best	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,36 @@
+#
+# number of interfaces
+#
+ 3
+#
+# for each interface below, we give the number of points and then x,z for each point
+#
+#
+# interface number 1 (bottom of the mesh)
+#
+ 2
+ 0 0
+ 2500 0
+#
+# interface number 2
+#
+ 2
+    0  750
+ 2500  750
+#
+# interface number 3 (topography, top of the mesh)
+#
+ 2
+    0 2000
+ 2500 2000
+#
+# for each layer, we give the number of spectral elements in the vertical direction
+#
+#
+# layer number 1 (bottom layer)
+#
+ 15
+#
+# layer number 2 (top layer)
+#
+ 25

Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/interfaces_fair
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/interfaces_fair	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/interfaces_fair	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,36 @@
+#
+# number of interfaces
+#
+ 3
+#
+# for each interface below, we give the number of points and then x,z for each point
+#
+#
+# interface number 1 (bottom of the mesh)
+#
+ 2
+ 0 0
+ 2000 0
+#
+# interface number 2
+#
+ 2
+    0  750
+ 2000  750
+#
+# interface number 3 (topography, top of the mesh)
+#
+ 2
+    0 1500
+ 2000 1500
+#
+# for each layer, we give the number of spectral elements in the vertical direction
+#
+#
+# layer number 1 (bottom layer)
+#
+ 15
+#
+# layer number 2 (top layer)
+#
+ 15

Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/interfaces_good
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/interfaces_good	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_1/interfaces_good	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,36 @@
+#
+# number of interfaces
+#
+ 3
+#
+# for each interface below, we give the number of points and then x,z for each point
+#
+#
+# interface number 1 (bottom of the mesh)
+#
+ 2
+ 0 0
+ 2500 0
+#
+# interface number 2
+#
+ 2
+    0  750
+ 2500  750
+#
+# interface number 3 (topography, top of the mesh)
+#
+ 2
+    0 1500
+ 2500 1500
+#
+# for each layer, we give the number of spectral elements in the vertical direction
+#
+#
+# layer number 1 (bottom layer)
+#
+ 15
+#
+# layer number 2 (top layer)
+#
+ 15

Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_2/Par_file_fair
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_2/Par_file_fair	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_2/Par_file_fair	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,138 @@
+# title of job
+title                           = PASSIVE IMAGING - 2 INTERFACES
+
+# forward or adjoint simulation
+SIMULATION_TYPE                 = 1
+NOISE_TOMOGRAPHY                = 1
+SAVE_FORWARD                    = .false.
+
+# parameters concerning partitioning
+nproc                           = 1              # number of processes
+partitioning_method             = 3              # SCOTCH = 3, ascending order (very bad idea) = 1
+PERFORM_CUTHILL_MCKEE           = .true.         # perform inverse Cuthill-McKee (1969) optimization/permutation for mesh numbering
+
+ngnod                           = 9              # number of control nodes per element (4 or 9)
+initialfield                    = .false.        # use a plane wave as source or not
+add_Bielak_conditions           = .false.        # add Bielak conditions or not if initial plane wave
+assign_external_model           = .false.        # define external earth model or not
+READ_EXTERNAL_SEP_FILE          = .false.        # Read external SEP file from DATA/model_velocity.dat_input, or use routine
+ATTENUATION_VISCOELASTIC_SOLID  = .false.        # turn attenuation (viscoelasticity) on or off for non-poroelastic solid parts of the model
+ATTENUATION_PORO_FLUID_PART     = .false.        # turn viscous attenuation on or off for the fluid part of poroelastic parts of the model
+Q0                              =  1             # quality factor for viscous attenuation
+freq0                           =  10            # frequency for viscous attenuation
+p_sv                            = .false.        # set the type of calculation (P-SV or SH/membrane waves)
+
+# time step parameters
+nt                              = 3000           # total number of time steps
+deltat                          = 1.d-3          # duration of a time step
+USER_T0                         = 0.0d0          # use this t0 as earliest starting time rather than the automatically calculated one
+time_stepping_scheme            = 1   # 1 = Newmark (2nd order),     2 = LDDRK4-6 (4th-order 6-stage low storage Runge-Kutta),     3 = classical 4th-order 4-stage Runge-Kutta
+
+# source parameters
+NSOURCES                        = 1              # number of sources [source info read in CMTSOLUTION file]
+force_normal_to_surface         = .false.        # angleforce normal to surface (external mesh and curve file needed)
+
+# constants for attenuation
+N_SLS                           = 2                      # number of standard linear solids for attenuation
+f0_attenuation                  = 5.196152422706633      # (Hz) relevant only if source is a Dirac or a Heaviside, else it is f0
+
+# receiver set parameters for seismograms
+seismotype                      = 1              # record 1=displ 2=veloc 3=accel 4=pressure
+generate_STATIONS               = .true.         # creates a STATION file in ./DATA
+nreceiversets                   = 1              # number of receiver sets
+anglerec                        = 0.d0           # angle to rotate components at receivers
+rec_normal_to_surface           = .false.        # base anglerec normal to surface (external mesh and curve file needed)
+SU_FORMAT                       = .false.        # output seismograms in Seismic Unix format (adjoint traces will be read in the same format)
+
+# first receiver set (repeat these 6 lines and adjust nreceiversets  accordingly)
+nrec                            = 29             # number of receivers
+xdeb                            = 150.           # first receiver x in meters
+zdeb                            = 1500.          # first receiver z in meters
+xfin                            = 1850.          # last receiver x in meters (ignored if onlyone receiver)
+zfin                            = 1500.          # last receiver z in meters (ignored if onlyone receiver)
+enreg_surf_same_vertical        = .true.         # receivers inside the medium or at the surface
+
+# display parameters
+NTSTEP_BETWEEN_OUTPUT_INFO      = 250            # display frequency in time steps
+output_postscript_snapshot      = .false.        # output Postscript snapshot of the results
+output_color_image              = .false.        # output color image of the results
+imagetype                       = 1              # display 1=displ 2=veloc 3=accel 4=pressure
+cutsnaps                        = 1.             # minimum amplitude in % for snapshots
+meshvect                        = .true.         # display mesh on vector plots or not
+modelvect                       = .false.        # display velocity model on vector plots
+boundvect                       = .true.         # display boundary conditions on plots
+interpol                        = .true.         # interpolation of the display or not
+pointsdisp                      = 6              # points for interpolation of display (set to 1 for lower-left corner only)
+subsamp_postscript              = 1              # subsampling of color snapshots
+factor_subsample_image          = 1              # factor to subsample color images output by the code (useful for very large models)
+POWER_DISPLAY_COLOR             = 0.30d0         # non linear display to enhance small amplitudes in color images
+DRAW_WATER_CONSTANT_BLUE_IN_JPG = .true.         # display acoustic layers as constant blue in JPEG images, because they likely correspond to water
+sizemax_arrows                  = 1.d0           # maximum size of arrows on vector plots in cm
+US_LETTER                       = .false.        # US letter paper or European A4
+USE_SNAPSHOT_NUMBER_IN_FILENAME = .false.        # use snapshot number in the file name of JPEG color snapshots instead of the time step
+gnuplot                         = .false.        # generate a GNUPLOT file for the grid
+output_grid                     = .false.        # save the grid in a text file or not
+output_energy                   = .false.        # compute and output acoustic and elastic energy (slows down the code significantly)
+output_wavefield_snapshot       = .false.         # output Ux,Uy,Uz text file for each output time (big files)
+
+# velocity and density models
+nbmodels                        = 3              # nb of different models
+# define models as
+# I:   (model_number 1 rho Vp Vs 0 0 Qkappa Qmu 0 0 0 0 0 0) or
+# II:  (model_number 2 rho c11 c13 c15 c33 c35 c55 0 0 0 0 0 0) or
+# III: (model_number 3 rhos rhof phi c kxx kxz kzz Ks Kf Kfr etaf mufr Qmu).
+# For istropic elastic/acoustic material use I and set Vs to zero to make a given model acoustic, for anisotropic elastic use II,
+# and for isotropic poroelastic material use III. The mesh can contain acoustic, elastic, and poroelastic models simultaneously.
+1 1 2700.d0 3000.d0 1732.051d0 0 0 9999 9999 0 0 0 0 0 0
+2 1 2500.d0 2700.d0 1600.000d0 0 0 9999 9999 0 0 0 0 0 0
+3 1 2200.d0 2500.d0 1443.375d0 0 0 9999 9999 0 0 0 0 0 0
+
+# external mesh or not
+read_external_mesh              = .false.
+
+# absorbing boundary active or not
+absorbing_conditions            = .true.
+
+# for horizontal periodic conditions: detect common points between left and right edges
+ADD_PERIODIC_CONDITIONS         = .false.
+
+# horizontal periodicity distance for periodic conditions
+PERIODIC_horiz_dist             = 0.3597d0
+
+# grid point detection tolerance for periodic conditions
+PERIODIC_DETECT_TOL             = 3.3334d-6
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR EXTERNAL MESHING
+
+# data concerning mesh, when generated using third-party app (more info in README)
+# (see also absorbing_conditions above)
+mesh_file                       = ./DATA/Mesh_canyon/canyon_mesh_file   # file containing the mesh
+nodes_coords_file               = ./DATA/Mesh_canyon/canyon_nodes_coords_file   # file containing the nodes coordinates
+materials_file                  = ./DATA/Mesh_canyon/canyon_materials_file   # file containing the material number for each element
+free_surface_file               = ./DATA/Mesh_canyon/canyon_free_surface_file   # file containing the free surface
+absorbing_surface_file          = ./DATA/Mesh_canyon/canyon_absorbing_surface_file   # file containing the absorbing surface
+tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the curve delimiting the velocity model
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR INTERNAL MESHING
+
+# file containing interfaces for internal mesh
+interfacesfile                  = interfaces.dat
+
+# geometry of the model (origin lower-left corner = 0,0) and mesh description
+xmin                            = 0.d0           # abscissa of left side of the model
+xmax                            = 2000.d0        # abscissa of right side of the model
+nx                              = 40             # number of elements along X
+
+# absorbing boundary parameters (see absorbing_conditions above)
+absorbbottom                    = .true.
+absorbright                     = .true.
+absorbtop                       = .false.
+absorbleft                      = .true.
+
+# define the different regions of the model in the (nx,nz) spectral element mesh
+nbregions                       = 3              # nb of regions and model number for each
+1 40  1 10 1
+1 40 11 20 2
+1 40 21 30 3

Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_2/Par_file_fair.before_update_to_r19xxx
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_2/Par_file_fair.before_update_to_r19xxx	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_2/Par_file_fair.before_update_to_r19xxx	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,137 @@
+# title of job
+title                           = PASSIVE IMAGING - 2 INTERFACES
+
+# forward or adjoint simulation
+SIMULATION_TYPE                 = 1
+NOISE_TOMOGRAPHY                = 1
+SAVE_FORWARD                    = .false.
+
+# parameters concerning partitioning
+nproc                           = 1              # number of processes
+partitioning_method             = 3              # SCOTCH = 3, ascending order (very bad idea) = 1
+PERFORM_CUTHILL_MCKEE           = .true.         # perform inverse Cuthill-McKee (1969) optimization/permutation for mesh numbering
+
+ngnod                           = 9              # number of control nodes per element (4 or 9)
+initialfield                    = .false.        # use a plane wave as source or not
+add_Bielak_conditions           = .false.        # add Bielak conditions or not if initial plane wave
+assign_external_model           = .false.        # define external earth model or not
+READ_EXTERNAL_SEP_FILE          = .false.        # Read external SEP file from DATA/model_velocity.dat_input, or use routine
+ATTENUATION_VISCOELASTIC_SOLID  = .false.        # turn attenuation (viscoelasticity) on or off for non-poroelastic solid parts of the model
+ATTENUATION_PORO_FLUID_PART     = .false.        # turn viscous attenuation on or off for the fluid part of poroelastic parts of the model
+Q0                              =  1             # quality factor for viscous attenuation
+freq0                           =  10            # frequency for viscous attenuation
+p_sv                            = .false.        # set the type of calculation (P-SV or SH/membrane waves)
+
+# time step parameters
+nt                              = 3000           # total number of time steps
+deltat                          = 1.d-3          # duration of a time step
+USER_T0                         = 0.0d0          # use this t0 as earliest starting time rather than the automatically calculated one
+
+# source parameters
+NSOURCES                        = 1              # number of sources [source info read in CMTSOLUTION file]
+force_normal_to_surface         = .false.        # angleforce normal to surface (external mesh and curve file needed)
+
+# constants for attenuation
+N_SLS                           = 2                      # number of standard linear solids for attenuation
+f0_attenuation                  = 5.196152422706633      # (Hz) relevant only if source is a Dirac or a Heaviside, else it is f0
+
+# receiver set parameters for seismograms
+seismotype                      = 1              # record 1=displ 2=veloc 3=accel 4=pressure
+generate_STATIONS               = .true.         # creates a STATION file in ./DATA
+nreceiversets                   = 1              # number of receiver sets
+anglerec                        = 0.d0           # angle to rotate components at receivers
+rec_normal_to_surface           = .false.        # base anglerec normal to surface (external mesh and curve file needed)
+SU_FORMAT                       = .false.        # output seismograms in Seismic Unix format (adjoint traces will be read in the same format)
+
+# first receiver set (repeat these 6 lines and adjust nreceiversets  accordingly)
+nrec                            = 29             # number of receivers
+xdeb                            = 150.           # first receiver x in meters
+zdeb                            = 1500.          # first receiver z in meters
+xfin                            = 1850.          # last receiver x in meters (ignored if onlyone receiver)
+zfin                            = 1500.          # last receiver z in meters (ignored if onlyone receiver)
+enreg_surf_same_vertical        = .true.         # receivers inside the medium or at the surface
+
+# display parameters
+NTSTEP_BETWEEN_OUTPUT_INFO      = 250            # display frequency in time steps
+output_postscript_snapshot      = .false.        # output Postscript snapshot of the results
+output_color_image              = .false.        # output color image of the results
+imagetype                       = 1              # display 1=displ 2=veloc 3=accel 4=pressure
+cutsnaps                        = 1.             # minimum amplitude in % for snapshots
+meshvect                        = .true.         # display mesh on vector plots or not
+modelvect                       = .false.        # display velocity model on vector plots
+boundvect                       = .true.         # display boundary conditions on plots
+interpol                        = .true.         # interpolation of the display or not
+pointsdisp                      = 6              # points for interpolation of display (set to 1 for lower-left corner only)
+subsamp_postscript              = 1              # subsampling of color snapshots
+factor_subsample_image          = 1              # factor to subsample color images output by the code (useful for very large models)
+POWER_DISPLAY_COLOR             = 0.30d0         # non linear display to enhance small amplitudes in color images
+DRAW_WATER_CONSTANT_BLUE_IN_JPG = .true.         # display acoustic layers as constant blue in JPEG images, because they likely correspond to water
+sizemax_arrows                  = 1.d0           # maximum size of arrows on vector plots in cm
+US_LETTER                       = .false.        # US letter paper or European A4
+USE_SNAPSHOT_NUMBER_IN_FILENAME = .false.        # use snapshot number in the file name of JPEG color snapshots instead of the time step
+gnuplot                         = .false.        # generate a GNUPLOT file for the grid
+output_grid                     = .false.        # save the grid in a text file or not
+output_energy                   = .false.        # compute and output acoustic and elastic energy (slows down the code significantly)
+output_wavefield_snapshot       = .false.         # output Ux,Uy,Uz text file for each output time (big files)
+
+# velocity and density models
+nbmodels                        = 3              # nb of different models
+# define models as
+# I:   (model_number 1 rho Vp Vs 0 0 Qkappa Qmu 0 0 0 0 0 0) or
+# II:  (model_number 2 rho c11 c13 c15 c33 c35 c55 0 0 0 0 0 0) or
+# III: (model_number 3 rhos rhof phi c kxx kxz kzz Ks Kf Kfr etaf mufr Qmu).
+# For istropic elastic/acoustic material use I and set Vs to zero to make a given model acoustic, for anisotropic elastic use II,
+# and for isotropic poroelastic material use III. The mesh can contain acoustic, elastic, and poroelastic models simultaneously.
+1 1 2700.d0 3000.d0 1732.051d0 0 0 9999 9999 0 0 0 0 0 0
+2 1 2500.d0 2700.d0 1600.000d0 0 0 9999 9999 0 0 0 0 0 0
+3 1 2200.d0 2500.d0 1443.375d0 0 0 9999 9999 0 0 0 0 0 0
+
+# external mesh or not
+read_external_mesh              = .false.
+
+# absorbing boundary active or not
+absorbing_conditions            = .true.
+
+# for horizontal periodic conditions: detect common points between left and right edges
+ADD_PERIODIC_CONDITIONS         = .false.
+
+# horizontal periodicity distance for periodic conditions
+PERIODIC_horiz_dist             = 0.3597d0
+
+# grid point detection tolerance for periodic conditions
+PERIODIC_DETECT_TOL             = 3.3334d-6
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR EXTERNAL MESHING
+
+# data concerning mesh, when generated using third-party app (more info in README)
+# (see also absorbing_conditions above)
+mesh_file                       = ./DATA/Mesh_canyon/canyon_mesh_file   # file containing the mesh
+nodes_coords_file               = ./DATA/Mesh_canyon/canyon_nodes_coords_file   # file containing the nodes coordinates
+materials_file                  = ./DATA/Mesh_canyon/canyon_materials_file   # file containing the material number for each element
+free_surface_file               = ./DATA/Mesh_canyon/canyon_free_surface_file   # file containing the free surface
+absorbing_surface_file          = ./DATA/Mesh_canyon/canyon_absorbing_surface_file   # file containing the absorbing surface
+tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the curve delimiting the velocity model
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR INTERNAL MESHING
+
+# file containing interfaces for internal mesh
+interfacesfile                  = interfaces.dat
+
+# geometry of the model (origin lower-left corner = 0,0) and mesh description
+xmin                            = 0.d0           # abscissa of left side of the model
+xmax                            = 2000.d0        # abscissa of right side of the model
+nx                              = 40             # number of elements along X
+
+# absorbing boundary parameters (see absorbing_conditions above)
+absorbbottom                    = .true.
+absorbright                     = .true.
+absorbtop                       = .false.
+absorbleft                      = .true.
+
+# define the different regions of the model in the (nx,nz) spectral element mesh
+nbregions                       = 3              # nb of regions and model number for each
+1 40  1 10 1
+1 40 11 20 2
+1 40 21 30 3

Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_2/Par_file_good
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_2/Par_file_good	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_2/Par_file_good	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,138 @@
+# title of job
+title                           = PASSIVE IMAGING - 2 INTERFACES
+
+# forward or adjoint simulation
+SIMULATION_TYPE                 = 1
+NOISE_TOMOGRAPHY                = 1
+SAVE_FORWARD                    = .false.
+
+# parameters concerning partitioning
+nproc                           = 1              # number of processes
+partitioning_method             = 3              # SCOTCH = 3, ascending order (very bad idea) = 1
+PERFORM_CUTHILL_MCKEE           = .true.         # perform inverse Cuthill-McKee (1969) optimization/permutation for mesh numbering
+
+ngnod                           = 9              # number of control nodes per element (4 or 9)
+initialfield                    = .false.        # use a plane wave as source or not
+add_Bielak_conditions           = .false.        # add Bielak conditions or not if initial plane wave
+assign_external_model           = .false.        # define external earth model or not
+READ_EXTERNAL_SEP_FILE          = .false.        # Read external SEP file from DATA/model_velocity.dat_input, or use routine
+ATTENUATION_VISCOELASTIC_SOLID  = .false.        # turn attenuation (viscoelasticity) on or off for non-poroelastic solid parts of the model
+ATTENUATION_PORO_FLUID_PART     = .false.        # turn viscous attenuation on or off for the fluid part of poroelastic parts of the model
+Q0                              =  1             # quality factor for viscous attenuation
+freq0                           =  10            # frequency for viscous attenuation
+p_sv                            = .false.        # set the type of calculation (P-SV or SH/membrane waves)
+
+# time step parameters
+nt                              = 3000           # total number of time steps
+deltat                          = 1.d-3          # duration of a time step
+USER_T0                         = 0.0d0          # use this t0 as earliest starting time rather than the automatically calculated one
+time_stepping_scheme            = 1   # 1 = Newmark (2nd order),     2 = LDDRK4-6 (4th-order 6-stage low storage Runge-Kutta),     3 = classical 4th-order 4-stage Runge-Kutta
+
+# source parameters
+NSOURCES                        = 1              # number of sources [source info read in CMTSOLUTION file]
+force_normal_to_surface         = .false.        # angleforce normal to surface (external mesh and curve file needed)
+
+# constants for attenuation
+N_SLS                           = 2                      # number of standard linear solids for attenuation
+f0_attenuation                  = 5.196152422706633      # (Hz) relevant only if source is a Dirac or a Heaviside, else it is f0
+
+# receiver set parameters for seismograms
+seismotype                      = 1              # record 1=displ 2=veloc 3=accel 4=pressure
+generate_STATIONS               = .true.         # creates a STATION file in ./DATA
+nreceiversets                   = 1              # number of receiver sets
+anglerec                        = 0.d0           # angle to rotate components at receivers
+rec_normal_to_surface           = .false.        # base anglerec normal to surface (external mesh and curve file needed)
+SU_FORMAT                       = .false.        # output seismograms in Seismic Unix format (adjoint traces will be read in the same format)
+
+# first receiver set (repeat these 6 lines and adjust nreceiversets  accordingly)
+nrec                            = 29             # number of receivers
+xdeb                            = 400.           # first receiver x in meters
+zdeb                            = 1500.          # first receiver z in meters
+xfin                            = 2100.          # last receiver x in meters (ignored if onlyone receiver)
+zfin                            = 1500.          # last receiver z in meters (ignored if onlyone receiver)
+enreg_surf_same_vertical        = .true.         # receivers inside the medium or at the surface
+
+# display parameters
+NTSTEP_BETWEEN_OUTPUT_INFO      = 250            # display frequency in time steps
+output_postscript_snapshot      = .false.        # output Postscript snapshot of the results
+output_color_image              = .false.        # output color image of the results
+imagetype                       = 1              # display 1=displ 2=veloc 3=accel 4=pressure
+cutsnaps                        = 1.             # minimum amplitude in % for snapshots
+meshvect                        = .true.         # display mesh on vector plots or not
+modelvect                       = .false.        # display velocity model on vector plots
+boundvect                       = .true.         # display boundary conditions on plots
+interpol                        = .true.         # interpolation of the display or not
+pointsdisp                      = 6              # points for interpolation of display (set to 1 for lower-left corner only)
+subsamp_postscript              = 1              # subsampling of color snapshots
+factor_subsample_image          = 1              # factor to subsample color images output by the code (useful for very large models)
+POWER_DISPLAY_COLOR             = 0.30d0         # non linear display to enhance small amplitudes in color images
+DRAW_WATER_CONSTANT_BLUE_IN_JPG = .true.         # display acoustic layers as constant blue in JPEG images, because they likely correspond to water
+sizemax_arrows                  = 1.d0           # maximum size of arrows on vector plots in cm
+US_LETTER                       = .false.        # US letter paper or European A4
+USE_SNAPSHOT_NUMBER_IN_FILENAME = .false.        # use snapshot number in the file name of JPEG color snapshots instead of the time step
+gnuplot                         = .false.        # generate a GNUPLOT file for the grid
+output_grid                     = .false.        # save the grid in a text file or not
+output_energy                   = .false.        # compute and output acoustic and elastic energy (slows down the code significantly)
+output_wavefield_snapshot       = .false.         # output Ux,Uy,Uz text file for each output time (big files)
+
+# velocity and density models
+nbmodels                        = 3              # nb of different models
+# define models as
+# I:   (model_number 1 rho Vp Vs 0 0 Qkappa Qmu 0 0 0 0 0 0) or
+# II:  (model_number 2 rho c11 c13 c15 c33 c35 c55 0 0 0 0 0 0) or
+# III: (model_number 3 rhos rhof phi c kxx kxz kzz Ks Kf Kfr etaf mufr Qmu).
+# For istropic elastic/acoustic material use I and set Vs to zero to make a given model acoustic, for anisotropic elastic use II,
+# and for isotropic poroelastic material use III. The mesh can contain acoustic, elastic, and poroelastic models simultaneously.
+1 1 2700.d0 3000.d0 1732.051d0 0 0 9999 9999 0 0 0 0 0 0
+2 1 2500.d0 2700.d0 1600.000d0 0 0 9999 9999 0 0 0 0 0 0
+3 1 2200.d0 2500.d0 1443.375d0 0 0 9999 9999 0 0 0 0 0 0
+
+# external mesh or not
+read_external_mesh              = .false.
+
+# absorbing boundary active or not
+absorbing_conditions            = .true.
+
+# for horizontal periodic conditions: detect common points between left and right edges
+ADD_PERIODIC_CONDITIONS         = .false.
+
+# horizontal periodicity distance for periodic conditions
+PERIODIC_horiz_dist             = 0.3597d0
+
+# grid point detection tolerance for periodic conditions
+PERIODIC_DETECT_TOL             = 3.3334d-6
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR EXTERNAL MESHING
+
+# data concerning mesh, when generated using third-party app (more info in README)
+# (see also absorbing_conditions above)
+mesh_file                       = ./DATA/Mesh_canyon/canyon_mesh_file   # file containing the mesh
+nodes_coords_file               = ./DATA/Mesh_canyon/canyon_nodes_coords_file   # file containing the nodes coordinates
+materials_file                  = ./DATA/Mesh_canyon/canyon_materials_file   # file containing the material number for each element
+free_surface_file               = ./DATA/Mesh_canyon/canyon_free_surface_file   # file containing the free surface
+absorbing_surface_file          = ./DATA/Mesh_canyon/canyon_absorbing_surface_file   # file containing the absorbing surface
+tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the curve delimiting the velocity model
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR INTERNAL MESHING
+
+# file containing interfaces for internal mesh
+interfacesfile                  = interfaces.dat
+
+# geometry of the model (origin lower-left corner = 0,0) and mesh description
+xmin                            = 0.d0           # abscissa of left side of the model
+xmax                            = 2500.d0        # abscissa of right side of the model
+nx                              = 50             # number of elements along X
+
+# absorbing boundary parameters (see absorbing_conditions above)
+absorbbottom                    = .true.
+absorbright                     = .true.
+absorbtop                       = .false.
+absorbleft                      = .true.
+
+# define the different regions of the model in the (nx,nz) spectral element mesh
+nbregions                       = 3              # nb of regions and model number for each
+1 50  1 10 1
+1 50 11 20 2
+1 50 21 30 3

Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_2/Par_file_good.before_update_to_r19xxx
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_2/Par_file_good.before_update_to_r19xxx	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_2/Par_file_good.before_update_to_r19xxx	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,137 @@
+# title of job
+title                           = PASSIVE IMAGING - 2 INTERFACES
+
+# forward or adjoint simulation
+SIMULATION_TYPE                 = 1
+NOISE_TOMOGRAPHY                = 1
+SAVE_FORWARD                    = .false.
+
+# parameters concerning partitioning
+nproc                           = 1              # number of processes
+partitioning_method             = 3              # SCOTCH = 3, ascending order (very bad idea) = 1
+PERFORM_CUTHILL_MCKEE           = .true.         # perform inverse Cuthill-McKee (1969) optimization/permutation for mesh numbering
+
+ngnod                           = 9              # number of control nodes per element (4 or 9)
+initialfield                    = .false.        # use a plane wave as source or not
+add_Bielak_conditions           = .false.        # add Bielak conditions or not if initial plane wave
+assign_external_model           = .false.        # define external earth model or not
+READ_EXTERNAL_SEP_FILE          = .false.        # Read external SEP file from DATA/model_velocity.dat_input, or use routine
+ATTENUATION_VISCOELASTIC_SOLID  = .false.        # turn attenuation (viscoelasticity) on or off for non-poroelastic solid parts of the model
+ATTENUATION_PORO_FLUID_PART     = .false.        # turn viscous attenuation on or off for the fluid part of poroelastic parts of the model
+Q0                              =  1             # quality factor for viscous attenuation
+freq0                           =  10            # frequency for viscous attenuation
+p_sv                            = .false.        # set the type of calculation (P-SV or SH/membrane waves)
+
+# time step parameters
+nt                              = 3000           # total number of time steps
+deltat                          = 1.d-3          # duration of a time step
+USER_T0                         = 0.0d0          # use this t0 as earliest starting time rather than the automatically calculated one
+
+# source parameters
+NSOURCES                        = 1              # number of sources [source info read in CMTSOLUTION file]
+force_normal_to_surface         = .false.        # angleforce normal to surface (external mesh and curve file needed)
+
+# constants for attenuation
+N_SLS                           = 2                      # number of standard linear solids for attenuation
+f0_attenuation                  = 5.196152422706633      # (Hz) relevant only if source is a Dirac or a Heaviside, else it is f0
+
+# receiver set parameters for seismograms
+seismotype                      = 1              # record 1=displ 2=veloc 3=accel 4=pressure
+generate_STATIONS               = .true.         # creates a STATION file in ./DATA
+nreceiversets                   = 1              # number of receiver sets
+anglerec                        = 0.d0           # angle to rotate components at receivers
+rec_normal_to_surface           = .false.        # base anglerec normal to surface (external mesh and curve file needed)
+SU_FORMAT                       = .false.        # output seismograms in Seismic Unix format (adjoint traces will be read in the same format)
+
+# first receiver set (repeat these 6 lines and adjust nreceiversets  accordingly)
+nrec                            = 29             # number of receivers
+xdeb                            = 400.           # first receiver x in meters
+zdeb                            = 1500.          # first receiver z in meters
+xfin                            = 2100.          # last receiver x in meters (ignored if onlyone receiver)
+zfin                            = 1500.          # last receiver z in meters (ignored if onlyone receiver)
+enreg_surf_same_vertical        = .true.         # receivers inside the medium or at the surface
+
+# display parameters
+NTSTEP_BETWEEN_OUTPUT_INFO      = 250            # display frequency in time steps
+output_postscript_snapshot      = .false.        # output Postscript snapshot of the results
+output_color_image              = .false.        # output color image of the results
+imagetype                       = 1              # display 1=displ 2=veloc 3=accel 4=pressure
+cutsnaps                        = 1.             # minimum amplitude in % for snapshots
+meshvect                        = .true.         # display mesh on vector plots or not
+modelvect                       = .false.        # display velocity model on vector plots
+boundvect                       = .true.         # display boundary conditions on plots
+interpol                        = .true.         # interpolation of the display or not
+pointsdisp                      = 6              # points for interpolation of display (set to 1 for lower-left corner only)
+subsamp_postscript              = 1              # subsampling of color snapshots
+factor_subsample_image          = 1              # factor to subsample color images output by the code (useful for very large models)
+POWER_DISPLAY_COLOR             = 0.30d0         # non linear display to enhance small amplitudes in color images
+DRAW_WATER_CONSTANT_BLUE_IN_JPG = .true.         # display acoustic layers as constant blue in JPEG images, because they likely correspond to water
+sizemax_arrows                  = 1.d0           # maximum size of arrows on vector plots in cm
+US_LETTER                       = .false.        # US letter paper or European A4
+USE_SNAPSHOT_NUMBER_IN_FILENAME = .false.        # use snapshot number in the file name of JPEG color snapshots instead of the time step
+gnuplot                         = .false.        # generate a GNUPLOT file for the grid
+output_grid                     = .false.        # save the grid in a text file or not
+output_energy                   = .false.        # compute and output acoustic and elastic energy (slows down the code significantly)
+output_wavefield_snapshot       = .false.         # output Ux,Uy,Uz text file for each output time (big files)
+
+# velocity and density models
+nbmodels                        = 3              # nb of different models
+# define models as
+# I:   (model_number 1 rho Vp Vs 0 0 Qkappa Qmu 0 0 0 0 0 0) or
+# II:  (model_number 2 rho c11 c13 c15 c33 c35 c55 0 0 0 0 0 0) or
+# III: (model_number 3 rhos rhof phi c kxx kxz kzz Ks Kf Kfr etaf mufr Qmu).
+# For istropic elastic/acoustic material use I and set Vs to zero to make a given model acoustic, for anisotropic elastic use II,
+# and for isotropic poroelastic material use III. The mesh can contain acoustic, elastic, and poroelastic models simultaneously.
+1 1 2700.d0 3000.d0 1732.051d0 0 0 9999 9999 0 0 0 0 0 0
+2 1 2500.d0 2700.d0 1600.000d0 0 0 9999 9999 0 0 0 0 0 0
+3 1 2200.d0 2500.d0 1443.375d0 0 0 9999 9999 0 0 0 0 0 0
+
+# external mesh or not
+read_external_mesh              = .false.
+
+# absorbing boundary active or not
+absorbing_conditions            = .true.
+
+# for horizontal periodic conditions: detect common points between left and right edges
+ADD_PERIODIC_CONDITIONS         = .false.
+
+# horizontal periodicity distance for periodic conditions
+PERIODIC_horiz_dist             = 0.3597d0
+
+# grid point detection tolerance for periodic conditions
+PERIODIC_DETECT_TOL             = 3.3334d-6
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR EXTERNAL MESHING
+
+# data concerning mesh, when generated using third-party app (more info in README)
+# (see also absorbing_conditions above)
+mesh_file                       = ./DATA/Mesh_canyon/canyon_mesh_file   # file containing the mesh
+nodes_coords_file               = ./DATA/Mesh_canyon/canyon_nodes_coords_file   # file containing the nodes coordinates
+materials_file                  = ./DATA/Mesh_canyon/canyon_materials_file   # file containing the material number for each element
+free_surface_file               = ./DATA/Mesh_canyon/canyon_free_surface_file   # file containing the free surface
+absorbing_surface_file          = ./DATA/Mesh_canyon/canyon_absorbing_surface_file   # file containing the absorbing surface
+tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the curve delimiting the velocity model
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR INTERNAL MESHING
+
+# file containing interfaces for internal mesh
+interfacesfile                  = interfaces.dat
+
+# geometry of the model (origin lower-left corner = 0,0) and mesh description
+xmin                            = 0.d0           # abscissa of left side of the model
+xmax                            = 2500.d0        # abscissa of right side of the model
+nx                              = 50             # number of elements along X
+
+# absorbing boundary parameters (see absorbing_conditions above)
+absorbbottom                    = .true.
+absorbright                     = .true.
+absorbtop                       = .false.
+absorbleft                      = .true.
+
+# define the different regions of the model in the (nx,nz) spectral element mesh
+nbregions                       = 3              # nb of regions and model number for each
+1 50  1 10 1
+1 50 11 20 2
+1 50 21 30 3

Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_2/SOURCE_noise
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_2/SOURCE_noise	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_2/SOURCE_noise	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,13 @@
+#source 1.  The components of a moment tensor source must be given in N.m, not in dyne.cm as in the DATA/CMTSOLUTION source file of the 3D version of the code.
+source_surf                     = .false.        # source inside the medium or at the surface
+xs                              = 0.
+zs                              = 0.
+source_type                     = 1              # elastic force or acoustic pressure = 1 or moment tensor = 2
+time_function_type              = 3              # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4, Heaviside = 5
+f0                              = 10.0           # dominant source frequency (Hz) if not Dirac or Heaviside
+tshift                          = 0.0            # time shift when multi sources (if one source, must be zero)
+angleforce                      = 0.             # angle of the source (for a force only)
+Mxx                             = 0.d0           # Mxx component (for a moment tensor source only)
+Mzz                             = 0.d0           # Mzz component (for a moment tensor source only)
+Mxz                             = 0.d0           # Mxz component (for a moment tensor source only)
+factor                          = 0.d0          # amplification factor

Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_2/interfaces_fair
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_2/interfaces_fair	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_2/interfaces_fair	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,46 @@
+#
+# number of interfaces
+#
+ 4
+#
+# for each interface below, we give the number of points and then x,z for each point
+#
+#
+# interface number 1 (bottom of the mesh)
+#
+ 2
+ 0 0
+ 2000 0
+#
+# interface number 2
+#
+ 2
+    0  500
+ 2000  500
+#
+# interface number 3
+#
+ 2
+    0 1000
+ 2000 1000
+#
+# interface number 4 (topography, top of the mesh)
+#
+ 2
+    0 1500
+ 2000 1500
+#
+# for each layer, we give the number of spectral elements in the vertical direction
+#
+#
+# layer number 1 (bottom layer)
+#
+ 10
+#
+# layer number 2
+#
+ 10
+#
+# layer number 3 (top layer)
+#
+ 10

Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_2/interfaces_good
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_2/interfaces_good	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/model_2/interfaces_good	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,46 @@
+#
+# number of interfaces
+#
+ 4
+#
+# for each interface below, we give the number of points and then x,z for each point
+#
+#
+# interface number 1 (bottom of the mesh)
+#
+ 2
+ 0 0
+ 2500 0
+#
+# interface number 2
+#
+ 2
+    0  500
+ 2500  500
+#
+# interface number 3
+#
+ 2
+    0 1000
+ 2500 1000
+#
+# interface number 4 (topography, top of the mesh)
+#
+ 2
+    0 1500
+ 2500 1500
+#
+# for each layer, we give the number of spectral elements in the vertical direction
+#
+#
+# layer number 1 (bottom layer)
+#
+ 10
+#
+# layer number 2
+#
+ 10
+#
+# layer number 3 (top layer)
+#
+ 10

Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/replicate
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/replicate	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/replicate	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,41 @@
+#!/bin/bash
+
+if [ $# -ne 2 ]; then echo "Usage:  replicate  par_file  dir"; exit 1; fi
+
+if [ ! -e $1 ];  then echo "File not found: $1"; exit 1; fi
+if [ ! -d $2 ];  then echo "Directory not found: $2"; exit 1; fi
+
+FILE=$1
+DIR=$2
+
+FILE_BAK=$DIR/replicate.bak
+cp $FILE $FILE_BAK
+cd $DIR
+
+Line1=`grep ^SIMULATION_TYPE $FILE_BAK`
+Line2=`grep ^NOISE_TOMOGRAPHY $FILE_BAK`
+Line3=`grep ^SAVE_FORWARD $FILE_BAK`
+
+if [ `echo $Line1 | wc -l` -ne 1 ]; then echo "Error reading SIMULATION_TYPE."; exit 1; fi
+if [ `echo $Line2 | wc -l` -ne 1 ]; then echo "Error reading NOISE_TOMOGRAPHY."; exit 1; fi
+if [ `echo $Line3 | wc -l` -ne 1 ]; then echo "Error reading SAVE_FORWARD."; exit 1; fi
+
+
+#PAR_FILE_NOISE_1
+cp $FILE_BAK Par_file_noise_1
+sed -i -e "s:${Line1}:SIMULATION_TYPE                 = 1:"  Par_file_noise_1
+sed -i -e "s:${Line2}:NOISE_TOMOGRAPHY                = 1:"  Par_file_noise_1
+sed -i -e "s:${Line3}:SAVE_FORWARD                    = .false.:"  Par_file_noise_1
+
+#PAR_FILE_NOISE_2
+cp $FILE_BAK Par_file_noise_2
+sed -i -e "s:${Line1}:SIMULATION_TYPE                 = 1:"  Par_file_noise_2
+sed -i -e "s:${Line2}:NOISE_TOMOGRAPHY                = 2:"  Par_file_noise_2
+sed -i -e "s:${Line3}:SAVE_FORWARD                    = .true.:"  Par_file_noise_2
+
+#PAR_FILE_NOISE_3
+cp $FILE_BAK Par_file_noise_3
+sed -i -e "s:${Line1}:SIMULATION_TYPE                 = 2:"  Par_file_noise_3
+sed -i -e "s:${Line2}:NOISE_TOMOGRAPHY                = 3:"  Par_file_noise_3
+sed -i -e "s:${Line3}:SAVE_FORWARD                    = .false.:"  Par_file_noise_3
+


Property changes on: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/replicate
___________________________________________________________________
Name: svn:executable
   + *

Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/run_adj
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/run_adj	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/run_adj	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,58 @@
+#!/bin/bash
+#
+# noise tomography simulation 3
+#
+
+if [ $# -ne 2 ]; then echo "Usage:  run_adj  nrec  irec"; exit 1; fi
+
+
+NREC=$1
+IREC=$2
+ISRC=$IREC
+
+rm log.*
+rm DATA/Par_file
+#rm OUTPUT_FILES/image*
+rm OUTPUT_FILES/*.semd
+
+
+#create adjoint sources
+CODE_ADJ=../adj_mask.f90
+gfortran $CODE_ADJ -o adj_run
+
+cd SEM
+for ((ii=1; ii<=$NREC; ii++))
+do
+
+    WF_DIR=`printf ../../WF_DIR/noise_cc_%02d/step_2/ $IREC`
+    YTRACE=`printf "S%04d.AA.BXY" $ii`
+    cp ../OUTPUT_ALL/step_2/$YTRACE.semd .
+    TOFF=`../../getdirect $NREC $ISRC $ii` 
+    ../adj_run $YTRACE.semd $TOFF >> ../log.adj
+    rename '.semd' '' $YTRACE.semd.adj
+
+    XTRACE=`printf "S%04d.AA.BXX" $ii`
+    awk '{printf(" %20.10f %20.10f\n",$1,0.)}' < $YTRACE.adj > $XTRACE.adj
+
+    ZTRACE=`printf "S%04d.AA.BXZ" $ii`
+    awk '{printf(" %20.10f %20.10f\n",$1,0.)}' < $YTRACE.adj > $ZTRACE.adj
+
+done
+cd ..
+
+
+########### SIMULATION 3 ##########
+echo "Simulation 3"
+cp Par_file_noise_3  DATA/Par_file
+./xmeshfem2D > log.mesh3
+./xspecfem2D > log.spec3
+
+
+#save results
+mkdir OUTPUT_ALL/step_3
+cp log.*                  OUTPUT_ALL/step_3
+cp SEM/*Y.adj             OUTPUT_ALL/step_3
+cp DATA/Par_file          OUTPUT_ALL/step_3
+cp OUTPUT_FILES/proc*     OUTPUT_ALL/step_3
+#cp OUTPUT_FILES/image*    OUTPUT_ALL/step_3
+cp OUTPUT_FILES/snapshot* OUTPUT_ALL/


Property changes on: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/run_adj
___________________________________________________________________
Name: svn:executable
   + *

Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/run_all
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/run_all	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/run_all	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,21 @@
+#!/bin/bash
+
+nrec=29
+
+dir=`basename $PWD`
+
+../clean;
+
+for ((irec=1; irec<=$nrec; irec=irec+2))
+do
+
+    ii=`printf '%02d' $irec`
+    echo -e "\nRECEIVER $ii of $nrec"
+
+    ../run_fwd $nrec $irec
+    ../run_adj $nrec $irec
+    rm -rf ../${dir}_${ii}; mv OUTPUT_ALL ../${dir}_${ii}
+    ../clean;
+
+done
+


Property changes on: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/run_all
___________________________________________________________________
Name: svn:executable
   + *

Added: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/run_fwd
===================================================================
--- seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/run_fwd	                        (rev 0)
+++ seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/run_fwd	2012-01-24 19:11:25 UTC (rev 19462)
@@ -0,0 +1,72 @@
+#!/bin/bash
+#
+# runs first and second noise tomography simulations
+#
+
+if [ $# -ne 2 ]; then echo "Usage:  run_fwd  nrec  irec"; exit 1; fi
+
+
+NREC=$1
+IREC=$2
+
+RUNDIR=$PWD
+CODE_MASK=$PWD/`basename $PWD`.f90
+
+
+# prepare directories
+rm -rf   DATA SEM OUTPUT_FILES OUTPUT_ALL
+mkdir -p DATA SEM OUTPUT_FILES OUTPUT_ALL DATA/NOISE_TOMOGRAPHY OUTPUT_FILES/NOISE_TOMOGRAPHY
+
+
+# prepare files
+../changerec . $NREC &> /dev/null
+
+cp SOURCE_noise DATA/SOURCE
+cp interfaces*  DATA/
+echo $IREC >    DATA/NOISE_TOMOGRAPHY/irec_master
+if [ -f model* ]; then cp model* DATA/; fi
+if [ -f S_squared ]; then cp S_squared DATA/NOISE_TOMOGRAPHY; fi
+
+
+# compile
+echo "Compiling."
+cd ../../..
+cp $RUNDIR/Par_file_noise_1  DATA/Par_file
+if [ -f $CODE_MASK ]; then cp $CODE_MASK src/specfem2D/noise_tomography.f90; fi
+make > $RUNDIR/log.make
+ln -s $PWD/bin/xmeshfem2D -t $RUNDIR
+ln -s $PWD/bin/xspecfem2D -t $RUNDIR
+cd $RUNDIR
+
+
+
+########## SIMULATION 1 ##########
+echo "Simulation 1"
+cp Par_file_noise_1  DATA/Par_file
+./xmeshfem2D > log.mesh1
+./xspecfem2D > log.spec1
+
+# save results
+mkdir OUTPUT_ALL/step_1
+mv log.*                  OUTPUT_ALL/step_1
+mv DATA/Par_file          OUTPUT_ALL/step_1
+#mv OUTPUT_FILES/image*    OUTPUT_ALL/step_1
+mv OUTPUT_FILES/*Y.semd   OUTPUT_ALL/step_1
+mv OUTPUT_FILES/mask_*    OUTPUT_ALL/
+mv OUTPUT_FILES/mesh_???? OUTPUT_ALL/
+
+
+
+########## NOISE SIMULATION 2 ##########
+echo "Simulation 2"
+cp Par_file_noise_2  DATA/Par_file
+./xmeshfem2D > log.mesh2
+./xspecfem2D > log.spec2
+
+# save results
+mkdir OUTPUT_ALL/step_2
+cp log.*                  OUTPUT_ALL/step_2
+cp DATA/Par_file          OUTPUT_ALL/step_2
+#cp OUTPUT_FILES/image*    OUTPUT_ALL/step_2
+cp OUTPUT_FILES/*Y.semd   OUTPUT_ALL/step_2
+


Property changes on: seismo/2D/SPECFEM2D/branches/new_branch_for_Xie_Zhinan/trunk/EXAMPLES/noise_layered/run_fwd
___________________________________________________________________
Name: svn:executable
   + *



More information about the CIG-COMMITS mailing list