[cig-commits] r19642 - seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Thu Feb 16 07:18:09 PST 2012
Author: xie.zhinan
Date: 2012-02-16 07:18:08 -0800 (Thu, 16 Feb 2012)
New Revision: 19642
Modified:
seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/compute_arrays_source.f90
seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/compute_forces_viscoelastic.f90
seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/gmat01.f90
seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/prepare_absorb.f90
seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/read_databases.f90
seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/specfem2D.F90
seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/write_seismograms.F90
Log:
add adjoint implementation for viscoelastic simulation
Modified: seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/compute_arrays_source.f90
===================================================================
--- seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/compute_arrays_source.f90 2012-02-16 05:45:45 UTC (rev 19641)
+++ seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/compute_arrays_source.f90 2012-02-16 15:18:08 UTC (rev 19642)
@@ -129,26 +129,23 @@
! ------------------------------------------------------------------------------------------------------
-! subroutine compute_arrays_adj_source(adj_source_file,xi_receiver,gamma_receiver,adj_sourcearray, & !old
-! xigll,zigll,NSTEP) !old
- subroutine compute_arrays_adj_source(adj_source_file,xi_receiver,gamma_receiver,adj_sourcearray, & !xiezhinan
- xigll,zigll,NSTEP,stage_time_scheme) !xiezhinan
+ subroutine compute_arrays_adj_source(adj_source_file,xi_receiver,gamma_receiver,adj_sourcearray, &
+ xigll,zigll,NSTEP)
+
implicit none
include 'constants.h'
! input
integer NSTEP
- integer stage_time_scheme,i_stage !xiezhinan
double precision xi_receiver, gamma_receiver
character(len=*) adj_source_file
! output
-! real(kind=CUSTOM_REAL), dimension(NSTEP,3,NGLLX,NGLLZ) :: adj_sourcearray !old
- real(kind=CUSTOM_REAL), dimension(NSTEP,stage_time_scheme,3,NGLLX,NGLLZ) :: adj_sourcearray !xiezhinan
+ real(kind=CUSTOM_REAL), dimension(NSTEP,3,NGLLX,NGLLZ) :: adj_sourcearray
! Gauss-Lobatto-Legendre points of integration and weights
double precision, dimension(NGLLX) :: xigll
@@ -156,8 +153,7 @@
double precision :: hxir(NGLLX), hpxir(NGLLX), hgammar(NGLLZ), hpgammar(NGLLZ)
-! real(kind=CUSTOM_REAL) :: adj_src_s(NSTEP,3) !old
- real(kind=CUSTOM_REAL) :: adj_src_s(NSTEP,stage_time_scheme,3) !xiezhinan
+ real(kind=CUSTOM_REAL) :: adj_src_s(NSTEP,3)
integer icomp, itime, i, k, ios
double precision :: junk
@@ -167,8 +163,7 @@
call lagrange_any(xi_receiver,NGLLX,xigll,hxir,hpxir)
call lagrange_any(gamma_receiver,NGLLZ,zigll,hgammar,hpgammar)
-! adj_sourcearray(:,:,:,:) = 0. !old
- adj_sourcearray(:,:,:,:,:) = 0. !xiezhinan
+ adj_sourcearray(:,:,:,:) = 0.
comp = (/"BXX","BXY","BXZ"/)
@@ -179,10 +174,7 @@
if (ios /= 0) call exit_MPI(' file '//trim(filename)//'does not exist')
do itime = 1, NSTEP
- do i_stage = 1,stage_time_scheme !xiezhinan
-! read(IIN,*) junk, adj_src_s(itime,icomp) !old
- read(IIN,*) junk, adj_src_s(itime,i_stage,icomp) !xiezhinan
- enddo
+ read(IIN,*) junk, adj_src_s(itime,icomp)
enddo
close(IIN)
@@ -190,8 +182,7 @@
do k = 1, NGLLZ
do i = 1, NGLLX
-! adj_sourcearray(:,:,i,k) = hxir(i) * hgammar(k) * adj_src_s(:,:) !old
- adj_sourcearray(:,:,:,i,k) = hxir(i) * hgammar(k) * adj_src_s(:,:,:) !xiezhinan
+ adj_sourcearray(:,:,i,k) = hxir(i) * hgammar(k) * adj_src_s(:,:)
enddo
enddo
Modified: seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/compute_forces_viscoelastic.f90
===================================================================
--- seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/compute_forces_viscoelastic.f90 2012-02-16 05:45:45 UTC (rev 19641)
+++ seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/compute_forces_viscoelastic.f90 2012-02-16 15:18:08 UTC (rev 19642)
@@ -63,10 +63,10 @@
e1_LDDRK,e11_LDDRK,e13_LDDRK,alpha_LDDRK,beta_LDDRK, &
e1_initial_rk,e11_initial_rk,e13_initial_rk,e1_force_RK, e11_force_RK, e13_force_RK, &
stage_time_scheme,i_stage,&
- b_veloc_elastic,b_e1,b_e11,b_e13,&
- b_deltat,b_deltatcube,b_deltatfourth,b_twelvedeltat,b_fourdeltatsquare,&
- b_dux_dxl_n,b_duz_dzl_n,b_duz_dxl_n,b_dux_dzl_n,& !xiezhinan
- b_dux_dxl_np1,b_duz_dzl_np1,b_duz_dxl_np1,b_dux_dzl_np1,b_e1_LDDRK,b_e11_LDDRK,b_e13_LDDRK)
+ b_veloc_elastic,b_e1,b_e11,b_e13,&
+ b_deltat,b_deltatcube,b_deltatfourth,b_twelvedeltat,b_fourdeltatsquare,&
+ b_dux_dxl_n,b_duz_dzl_n,b_duz_dxl_n,b_dux_dzl_n,&
+ b_dux_dxl_np1,b_duz_dzl_np1,b_duz_dxl_np1,b_dux_dzl_np1)
! compute forces for the elastic elements
@@ -94,7 +94,9 @@
logical :: SAVE_FORWARD
double precision :: deltatcube,deltatfourth,twelvedeltat,fourdeltatsquare
- double precision :: b_deltat,b_deltatcube,b_deltatfourth,b_twelvedeltat,b_fourdeltatsquare !xiezhinan
+!xiezhinan
+ double precision :: b_deltat,b_deltatcube,b_deltatfourth,b_twelvedeltat,b_fourdeltatsquare
+!xiezhinan
double precision, dimension(NSOURCES) :: angleforce
integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
@@ -116,41 +118,41 @@
real(kind=CUSTOM_REAL), dimension(NSOURCES,NSTEP,stage_time_scheme) :: source_time_function
real(kind=CUSTOM_REAL), dimension(NSOURCES,NDIM,NGLLX,NGLLZ) :: sourcearray
-! real(kind=CUSTOM_REAL), dimension(3,nglob) :: b_accel_elastic,b_displ_elastic !old
- real(kind=CUSTOM_REAL), dimension(3,nglob) :: b_accel_elastic,b_veloc_elastic,b_displ_elastic !xiezhinan
- real(kind=CUSTOM_REAL), dimension(3,nglob) :: b_temp_displ_elastic !xiezhinan
-
-! real(kind=CUSTOM_REAL), dimension(nrec,NSTEP,3,NGLLX,NGLLZ) :: adj_sourcearrays !old
- real(kind=CUSTOM_REAL), dimension(nrec,NSTEP,stage_time_scheme,3,NGLLX,NGLLZ) :: adj_sourcearrays !xiezhinan
+! real(kind=CUSTOM_REAL), dimension(3,nglob) :: b_accel_elastic,b_displ_elastic
+!xiezhinan
+ real(kind=CUSTOM_REAL), dimension(3,nglob) :: b_accel_elastic,b_veloc_elastic,b_displ_elastic
+ real(kind=CUSTOM_REAL), dimension(3,nglob) :: b_temp_displ_elastic
+!xiezhinan
+ real(kind=CUSTOM_REAL), dimension(nrec,NSTEP,3,NGLLX,NGLLZ) :: adj_sourcearrays
real(kind=CUSTOM_REAL), dimension(nglob) :: mu_k,kappa_k
-! real(kind=CUSTOM_REAL), dimension(3,NGLLZ,nspec_left,NSTEP) :: b_absorb_elastic_left !old
-! real(kind=CUSTOM_REAL), dimension(3,NGLLZ,nspec_right,NSTEP) :: b_absorb_elastic_right !old
-! real(kind=CUSTOM_REAL), dimension(3,NGLLX,nspec_top,NSTEP) :: b_absorb_elastic_top !old
-! real(kind=CUSTOM_REAL), dimension(3,NGLLX,nspec_bottom,NSTEP) :: b_absorb_elastic_bottom !old
+ real(kind=CUSTOM_REAL), dimension(3,NGLLZ,nspec_left,NSTEP) :: b_absorb_elastic_left
+ real(kind=CUSTOM_REAL), dimension(3,NGLLZ,nspec_right,NSTEP) :: b_absorb_elastic_right
+ real(kind=CUSTOM_REAL), dimension(3,NGLLX,nspec_top,NSTEP) :: b_absorb_elastic_top
+ real(kind=CUSTOM_REAL), dimension(3,NGLLX,nspec_bottom,NSTEP) :: b_absorb_elastic_bottom
- real(kind=CUSTOM_REAL), dimension(3,NGLLZ,nspec_left,NSTEP,stage_time_scheme) :: b_absorb_elastic_left !xiezhinan
- real(kind=CUSTOM_REAL), dimension(3,NGLLZ,nspec_right,NSTEP,stage_time_scheme) :: b_absorb_elastic_right !xiezhinan
- real(kind=CUSTOM_REAL), dimension(3,NGLLX,nspec_top,NSTEP,stage_time_scheme) :: b_absorb_elastic_top !xiezhinan
- real(kind=CUSTOM_REAL), dimension(3,NGLLX,nspec_bottom,NSTEP,stage_time_scheme) :: b_absorb_elastic_bottom !xiezhinan
-
integer :: N_SLS
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e1,e11,e13
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: b_e1,b_e11,b_e13 !xiezhinan
+!xiezhinan
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: b_e1,b_e11,b_e13
+!xiezhinan
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e1_LDDRK,e11_LDDRK,e13_LDDRK
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: b_e1_LDDRK,b_e11_LDDRK,b_e13_LDDRK
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e1_initial_rk,e11_initial_rk,e13_initial_rk
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS,stage_time_scheme) :: e1_force_RK, e11_force_RK, e13_force_RK
double precision, dimension(NGLLX,NGLLZ,nspec,N_SLS) :: inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2
double precision, dimension(NGLLX,NGLLZ,nspec) :: Mu_nu1,Mu_nu2
real(kind=CUSTOM_REAL) :: e1_sum,e11_sum,e13_sum
- real(kind=CUSTOM_REAL) :: b_e1_sum,b_e11_sum,b_e13_sum !xiezhinan
+!xiezhinan
+ real(kind=CUSTOM_REAL) :: b_e1_sum,b_e11_sum,b_e13_sum
+!xiezhinan
integer :: i_sls
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: &
dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n,dux_dxl_np1,duz_dzl_np1,duz_dxl_np1,dux_dzl_np1
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: & !xiezhinan
- b_dux_dxl_n,b_duz_dzl_n,b_duz_dxl_n,b_dux_dzl_n, & !xiezhinan
- b_dux_dxl_np1,b_duz_dzl_np1,b_duz_dxl_np1,b_dux_dzl_np1 !xiezhinan
+!xiezhinan
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: &
+ b_dux_dxl_n,b_duz_dzl_n,b_duz_dxl_n,b_dux_dzl_n,&
+ b_dux_dxl_np1,b_duz_dzl_np1,b_duz_dxl_np1,b_dux_dzl_np1
+!xiezhinan
! derivatives of Lagrange polynomials
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx,hprimewgll_xx
@@ -196,7 +198,9 @@
! for attenuation
real(kind=CUSTOM_REAL) :: Un,Unp1,tauinv,Sn,Snp1,theta_n,theta_np1,tauinvsquare,tauinvcube,tauinvUn
+!xiezhinan
real(kind=CUSTOM_REAL) :: b_theta_n,b_theta_np1
+!xiezhinan
! for anisotropy
double precision :: c11,c15,c13,c33,c35,c55
@@ -226,15 +230,14 @@
! compute Grad(displ_elastic) at time step n for attenuation
if(ATTENUATION_VISCOELASTIC_SOLID) then
- if(SIMULATION_TYPE == 1)then
- temp_displ_elastic = displ_elastic + deltat * veloc_elastic + 0.5 * deltat**2 * accel_elastic
- call compute_gradient_attenuation(displ_elastic,dux_dxl_n,duz_dxl_n, &
- dux_dzl_n,duz_dzl_n,xix,xiz,gammax,gammaz,ibool,elastic,hprime_xx,hprime_zz,nspec,nglob)
- else
- b_temp_displ_elastic = b_displ_elastic + b_deltat * b_veloc_elastic + 0.5 * b_deltat**2 * b_accel_elastic
+ temp_displ_elastic = displ_elastic + deltat * veloc_elastic
+ call compute_gradient_attenuation(displ_elastic,dux_dxl_n,duz_dxl_n, &
+ dux_dzl_n,duz_dzl_n,xix,xiz,gammax,gammaz,ibool,elastic,hprime_xx,hprime_zz,nspec,nglob)
+!xiezhinan
+ b_temp_displ_elastic = b_displ_elastic + b_deltat * b_veloc_elastic
call compute_gradient_attenuation(b_displ_elastic,b_dux_dxl_n,b_duz_dxl_n, &
b_dux_dzl_n,b_duz_dzl_n,xix,xiz,gammax,gammaz,ibool,elastic,hprime_xx,hprime_zz,nspec,nglob)
- endif
+!xiezhinan
endif
ifirstelem = 1
@@ -373,7 +376,6 @@
lambdalplus2mul_relaxed_viscoel = lambdal_relaxed_viscoelastic + TWO*mul_relaxed_viscoelastic
! compute the stress using the unrelaxed Lame parameters (Carcione 2007 page 125)
- if(SIMULATION_TYPE == 1) then
sigma_xx = lambdaplus2mu_unrelaxed_elastic*dux_dxl + lambdal_unrelaxed_elastic*duz_dzl
sigma_xz = mul_unrelaxed_elastic*(duz_dxl + dux_dzl)
sigma_zz = lambdaplus2mu_unrelaxed_elastic*duz_dzl + lambdal_unrelaxed_elastic*dux_dxl
@@ -395,10 +397,9 @@
sigma_xz = sigma_xz + mul_relaxed_viscoelastic * e13_sum
sigma_zz = sigma_zz + (lambdal_relaxed_viscoelastic + mul_relaxed_viscoelastic) * e1_sum &
- TWO * mul_relaxed_viscoelastic * e11_sum
- endif
+ if(SIMULATION_TYPE == 2)then
- ! compute the stress using the unrelaxed Lame parameters (Carcione 2007 page 125) adjoint
- if(SIMULATION_TYPE == 2) then
+ ! compute the stress using the unrelaxed Lame parameters (Carcione 2007 page 125)
b_sigma_xx = lambdaplus2mu_unrelaxed_elastic*b_dux_dxl + lambdal_unrelaxed_elastic*b_duz_dzl
b_sigma_xz = mul_unrelaxed_elastic*(b_duz_dxl + b_dux_dzl)
b_sigma_zz = lambdaplus2mu_unrelaxed_elastic*b_duz_dzl + lambdal_unrelaxed_elastic*b_dux_dxl
@@ -413,14 +414,14 @@
b_e1_sum = b_e1_sum + b_e1(i,j,ispec,i_sls)
b_e11_sum = b_e11_sum + b_e11(i,j,ispec,i_sls)
b_e13_sum = b_e13_sum + b_e13(i,j,ispec,i_sls)
- enddo
+ enddo
b_sigma_xx = b_sigma_xx + (lambdal_relaxed_viscoelastic + mul_relaxed_viscoelastic) * b_e1_sum &
+ TWO * mul_relaxed_viscoelastic * b_e11_sum
b_sigma_xz = b_sigma_xz + mul_relaxed_viscoelastic * b_e13_sum
b_sigma_zz = b_sigma_zz + (lambdal_relaxed_viscoelastic + mul_relaxed_viscoelastic) * b_e1_sum &
- TWO * mul_relaxed_viscoelastic * b_e11_sum
- endif
+ endif
else
@@ -633,33 +634,20 @@
if(SAVE_FORWARD .and. SIMULATION_TYPE ==1) then
if(p_sv)then !P-SV waves
-! b_absorb_elastic_left(1,j,ib_left(ispecabs),it) = (tx + traction_x_t0)*weight !old
-! b_absorb_elastic_left(3,j,ib_left(ispecabs),it) = (tz + traction_z_t0)*weight !old
-
- b_absorb_elastic_left(1,j,ib_left(ispecabs),it,i_stage) = (tx + traction_x_t0)*weight !xiezhinan
- b_absorb_elastic_left(3,j,ib_left(ispecabs),it,i_stage) = (tz + traction_z_t0)*weight !xiezhinan
+ b_absorb_elastic_left(1,j,ib_left(ispecabs),it) = (tx + traction_x_t0)*weight
+ b_absorb_elastic_left(3,j,ib_left(ispecabs),it) = (tz + traction_z_t0)*weight
else !SH (membrane) waves
-! b_absorb_elastic_left(2,j,ib_left(ispecabs),it) = ty*weight !old
- b_absorb_elastic_left(2,j,ib_left(ispecabs),it,i_stage) = ty*weight
+ b_absorb_elastic_left(2,j,ib_left(ispecabs),it) = ty*weight
endif
elseif(SIMULATION_TYPE == 2) then
if(p_sv)then !P-SV waves
-! b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - & !old
-! b_absorb_elastic_left(1,j,ib_left(ispecabs),NSTEP-it+1) !old
-! b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) - & !old
-! b_absorb_elastic_left(3,j,ib_left(ispecabs),NSTEP-it+1) !old
-
- b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - & !xiezhinan
- b_absorb_elastic_left(1,j,ib_left(ispecabs),NSTEP-it+1,stage_time_scheme-i_stage+1) !xiezhinan
- b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) - & !xiezhinan
- b_absorb_elastic_left(3,j,ib_left(ispecabs),NSTEP-it+1,stage_time_scheme-i_stage+1) !xiezhinan
+ b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - &
+ b_absorb_elastic_left(1,j,ib_left(ispecabs),NSTEP-it+1)
+ b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) - &
+ b_absorb_elastic_left(3,j,ib_left(ispecabs),NSTEP-it+1)
else !SH (membrane) waves
-! b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - & !old
-! b_absorb_elastic_left(2,j,ib_left(ispecabs),NSTEP-it+1) !old
-
- b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - & !xiezhinan
- b_absorb_elastic_left(2,j,ib_left(ispecabs),NSTEP-it+1,stage_time_scheme-i_stage+1) !xiezhinan
-
+ b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - &
+ b_absorb_elastic_left(2,j,ib_left(ispecabs),NSTEP-it+1)
endif
endif
@@ -737,35 +725,20 @@
if(SAVE_FORWARD .and. SIMULATION_TYPE ==1) then
if(p_sv)then !P-SV waves
-! b_absorb_elastic_right(1,j,ib_right(ispecabs),it) = (tx - traction_x_t0)*weight !old
-! b_absorb_elastic_right(3,j,ib_right(ispecabs),it) = (tz - traction_z_t0)*weight !old
-
- b_absorb_elastic_right(1,j,ib_right(ispecabs),it,i_stage) = (tx - traction_x_t0)*weight !xiezhinan
- b_absorb_elastic_right(3,j,ib_right(ispecabs),it,i_stage) = (tz - traction_z_t0)*weight !xiezhinan
-
+ b_absorb_elastic_right(1,j,ib_right(ispecabs),it) = (tx - traction_x_t0)*weight
+ b_absorb_elastic_right(3,j,ib_right(ispecabs),it) = (tz - traction_z_t0)*weight
else! SH (membrane) waves
-! b_absorb_elastic_right(2,j,ib_right(ispecabs),it) = ty*weight !old
- b_absorb_elastic_right(2,j,ib_right(ispecabs),it,i_stage) = ty*weight !xiezhinan
+ b_absorb_elastic_right(2,j,ib_right(ispecabs),it) = ty*weight
endif
elseif(SIMULATION_TYPE == 2) then
if(p_sv)then !P-SV waves
-! b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - & !old
-! b_absorb_elastic_right(1,j,ib_right(ispecabs),NSTEP-it+1) !old
-! b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) - & !old
-! b_absorb_elastic_right(3,j,ib_right(ispecabs),NSTEP-it+1) !old
-
- b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - & !xiezhinan
- b_absorb_elastic_right(1,j,ib_right(ispecabs),NSTEP-it+1,stage_time_scheme-i_stage+1) !xiezhinan
- b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) - & !xiezhinan
- b_absorb_elastic_right(3,j,ib_right(ispecabs),NSTEP-it+1,stage_time_scheme-i_stage+1) !xiezhinan
-
+ b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - &
+ b_absorb_elastic_right(1,j,ib_right(ispecabs),NSTEP-it+1)
+ b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) - &
+ b_absorb_elastic_right(3,j,ib_right(ispecabs),NSTEP-it+1)
else! SH (membrane) waves
-! b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - & !old
-! b_absorb_elastic_right(2,j,ib_right(ispecabs),NSTEP-it+1) !old
-
- b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - & !xiezhinan
- b_absorb_elastic_right(2,j,ib_right(ispecabs),NSTEP-it+1,stage_time_scheme-i_stage+1) !xiezhinan
-
+ b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - &
+ b_absorb_elastic_right(2,j,ib_right(ispecabs),NSTEP-it+1)
endif
endif
@@ -858,34 +831,20 @@
if(SAVE_FORWARD .and. SIMULATION_TYPE ==1) then
if(p_sv)then !P-SV waves
-! b_absorb_elastic_bottom(1,i,ib_bottom(ispecabs),it) = (tx + traction_x_t0)*weight !old
-! b_absorb_elastic_bottom(3,i,ib_bottom(ispecabs),it) = (tz + traction_z_t0)*weight !old
-
- b_absorb_elastic_bottom(1,i,ib_bottom(ispecabs),it,i_stage) = (tx + traction_x_t0)*weight !xiezhinan
- b_absorb_elastic_bottom(3,i,ib_bottom(ispecabs),it,i_stage) = (tz + traction_z_t0)*weight !xiezhinan
-
+ b_absorb_elastic_bottom(1,i,ib_bottom(ispecabs),it) = (tx + traction_x_t0)*weight
+ b_absorb_elastic_bottom(3,i,ib_bottom(ispecabs),it) = (tz + traction_z_t0)*weight
else!SH (membrane) waves
-! b_absorb_elastic_bottom(2,i,ib_bottom(ispecabs),it) = ty*weight !old
- b_absorb_elastic_bottom(2,i,ib_bottom(ispecabs),it,i_stage) = ty*weight !xiezhinan
+ b_absorb_elastic_bottom(2,i,ib_bottom(ispecabs),it) = ty*weight
endif
elseif(SIMULATION_TYPE == 2) then
if(p_sv)then !P-SV waves
-! b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - & !old
-! b_absorb_elastic_bottom(1,i,ib_bottom(ispecabs),NSTEP-it+1) !old
-! b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) - & !old
-! b_absorb_elastic_bottom(3,i,ib_bottom(ispecabs),NSTEP-it+1) !old
-
- b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - & !xiezhinan
- b_absorb_elastic_bottom(1,i,ib_bottom(ispecabs),NSTEP-it+1,stage_time_scheme-i_stage+1) !xiezhinan
- b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) - & !xiezhinan
- b_absorb_elastic_bottom(3,i,ib_bottom(ispecabs),NSTEP-it+1,stage_time_scheme-i_stage+1) !xiezhinan
-
+ b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - &
+ b_absorb_elastic_bottom(1,i,ib_bottom(ispecabs),NSTEP-it+1)
+ b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) - &
+ b_absorb_elastic_bottom(3,i,ib_bottom(ispecabs),NSTEP-it+1)
else!SH (membrane) waves
-! b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - & !old
-! b_absorb_elastic_bottom(2,i,ib_bottom(ispecabs),NSTEP-it+1) !old
-
- b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - & !xiezhinan
- b_absorb_elastic_bottom(2,i,ib_bottom(ispecabs),NSTEP-it+1,stage_time_scheme-i_stage+1) !xiezhinan
+ b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - &
+ b_absorb_elastic_bottom(2,i,ib_bottom(ispecabs),NSTEP-it+1)
endif
endif
@@ -970,30 +929,17 @@
if(SAVE_FORWARD .and. SIMULATION_TYPE ==1) then
if(p_sv)then !P-SV waves
-! b_absorb_elastic_top(1,i,ib_top(ispecabs),it) = (tx- traction_x_t0)*weight !old
-! b_absorb_elastic_top(3,i,ib_top(ispecabs),it) = (tz- traction_z_t0)*weight !old
-
- b_absorb_elastic_top(1,i,ib_top(ispecabs),it,i_stage) = (tx- traction_x_t0)*weight !xiezhinan
- b_absorb_elastic_top(3,i,ib_top(ispecabs),it,i_stage) = (tz- traction_z_t0)*weight !xiezhinan
-
+ b_absorb_elastic_top(1,i,ib_top(ispecabs),it) = (tx- traction_x_t0)*weight
+ b_absorb_elastic_top(3,i,ib_top(ispecabs),it) = (tz- traction_z_t0)*weight
else!SH (membrane) waves
-! b_absorb_elastic_top(2,i,ib_top(ispecabs),it) = ty*weight !old
- b_absorb_elastic_top(2,i,ib_top(ispecabs),it,i_stage) = ty*weight !xiezhinan
+ b_absorb_elastic_top(2,i,ib_top(ispecabs),it) = ty*weight
endif
elseif(SIMULATION_TYPE == 2) then
if(p_sv)then !P-SV waves
-! b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - b_absorb_elastic_top(1,i,ib_top(ispecabs),NSTEP-it+1) !old
-! b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) - b_absorb_elastic_top(3,i,ib_top(ispecabs),NSTEP-it+1) !old
-
- b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - &
- b_absorb_elastic_top(1,i,ib_top(ispecabs),NSTEP-it+1,stage_time_scheme-i_stage+1) !xiezhinan
- b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) - &
- b_absorb_elastic_top(3,i,ib_top(ispecabs),NSTEP-it+1,stage_time_scheme-i_stage+1) !xiezhinan
-
+ b_accel_elastic(1,iglob) = b_accel_elastic(1,iglob) - b_absorb_elastic_top(1,i,ib_top(ispecabs),NSTEP-it+1)
+ b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) - b_absorb_elastic_top(3,i,ib_top(ispecabs),NSTEP-it+1)
else!SH (membrane) waves
-! b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - b_absorb_elastic_top(2,i,ib_top(ispecabs),NSTEP-it+1) !old
- b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - &
- b_absorb_elastic_top(2,i,ib_top(ispecabs),NSTEP-it+1,stage_time_scheme-i_stage+1) !xiezhinan
+ b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - b_absorb_elastic_top(2,i,ib_top(ispecabs),NSTEP-it+1)
endif
endif
@@ -1038,7 +984,6 @@
sourcearray(i_source,1,i,j)*source_time_function(i_source,NSTEP-it+1,stage_time_scheme-i_stage+1)
b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) + &
sourcearray(i_source,2,i,j)*source_time_function(i_source,NSTEP-it+1,stage_time_scheme-i_stage+1)
-
enddo
enddo
endif !endif SIMULATION_TYPE == 1
@@ -1061,18 +1006,11 @@
do j=1,NGLLZ
do i=1,NGLLX
iglob = ibool(i,j,ispec_selected_rec(irec))
- if(p_sv)then !P-SV waves
-! accel_elastic(1,iglob) = accel_elastic(1,iglob) + adj_sourcearrays(irec_local,NSTEP-it+1,1,i,j) !old
-! accel_elastic(3,iglob) = accel_elastic(3,iglob) + adj_sourcearrays(irec_local,NSTEP-it+1,3,i,j) !old
-
- accel_elastic(1,iglob) = accel_elastic(1,iglob) + & !xiezhinan
- adj_sourcearrays(irec_local,NSTEP-it+1,stage_time_scheme-i_stage+1,1,i,j) !xiezhinan
- accel_elastic(3,iglob) = accel_elastic(3,iglob) + & !xiezhinan
- adj_sourcearrays(irec_local,NSTEP-it+1,stage_time_scheme-i_stage+1,3,i,j) !xiezhinan
+ if(p_sv)then !P-SH waves
+ accel_elastic(1,iglob) = accel_elastic(1,iglob) + adj_sourcearrays(irec_local,NSTEP-it+1,1,i,j)
+ accel_elastic(3,iglob) = accel_elastic(3,iglob) + adj_sourcearrays(irec_local,NSTEP-it+1,3,i,j)
else !SH (membrane) waves
-! accel_elastic(2,iglob) = accel_elastic(2,iglob) + adj_sourcearrays(irec_local,NSTEP-it+1,2,i,j) !old
- accel_elastic(2,iglob) = accel_elastic(2,iglob) + & !xiezhinan
- adj_sourcearrays(irec_local,NSTEP-it+1,stage_time_scheme-i_stage+1,2,i,j) !xiezhinan
+ accel_elastic(2,iglob) = accel_elastic(2,iglob) + adj_sourcearrays(irec_local,NSTEP-it+1,2,i,j)
endif
enddo
enddo
@@ -1088,7 +1026,6 @@
! implement attenuation
if(ATTENUATION_VISCOELASTIC_SOLID) then
- if(SIMULATION_TYPE == 1) then
! compute Grad(displ_elastic) at time step n+1 for attenuation
call compute_gradient_attenuation(temp_displ_elastic,dux_dxl_np1,duz_dxl_np1, &
dux_dzl_np1,duz_dzl_np1,xix,xiz,gammax,gammaz,ibool,elastic,hprime_xx,hprime_zz,nspec,nglob)
@@ -1256,10 +1193,9 @@
enddo
enddo
enddo
- endif
if(SIMULATION_TYPE == 2) then
- ! compute Grad(displ_elastic) at time step n+1 for attenuation
+ ! compute Grad(displ_elastic) at time step n-1 for attenuation
call compute_gradient_attenuation(b_temp_displ_elastic,b_dux_dxl_np1,b_duz_dxl_np1, &
b_dux_dzl_np1,b_duz_dzl_np1,xix,xiz,gammax,gammaz,ibool,elastic,hprime_xx,hprime_zz,nspec,nglob)
@@ -1292,43 +1228,7 @@
b_e1(i,j,ispec,i_sls) = Unp1
endif
- if(stage_time_scheme == 6) then
- tauinv = inv_tau_sigma_nu1(i,j,ispec,i_sls)
- Un = b_e1(i,j,ispec,i_sls)
- temp_time_scheme = phi_nu1(i,j,ispec,i_sls)
- b_e1_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * b_e1_LDDRK(i,j,ispec,i_sls) + &
- b_deltat * (b_theta_n * temp_time_scheme - Un * tauinv)
- b_e1(i,j,ispec,i_sls) = Un + beta_LDDRK(i_stage) * b_e1_LDDRK(i,j,ispec,i_sls)
- endif
-
- if(stage_time_scheme == 4) then
-
- tauinv = inv_tau_sigma_nu1(i,j,ispec,i_sls)
- Un = e1(i,j,ispec,i_sls)
- temp_time_scheme = phi_nu1(i,j,ispec,i_sls)
- e1_force_RK(i,j,ispec,i_sls,i_stage) = deltat * (theta_n * temp_time_scheme - Un * tauinv)
-
- if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then
- if(i_stage == 1)weight_rk = 0.5d0
- if(i_stage == 2)weight_rk = 0.5d0
- if(i_stage == 3)weight_rk = 1.0d0
-
- if(i_stage==1)then
- e1_initial_rk(i,j,ispec,i_sls) = e1(i,j,ispec,i_sls)
- endif
-
- e1(i,j,ispec,i_sls) = e1_initial_rk(i,j,ispec,i_sls) + weight_rk * e1_force_RK(i,j,ispec,i_sls,i_stage)
-
- elseif(i_stage==4)then
-
- e1(i,j,ispec,i_sls) = e1_initial_rk(i,j,ispec,i_sls) + 1.0d0 / 6.0d0 * &
- (e1_force_RK(i,j,ispec,i_sls,1) + 2.0d0 * e1_force_RK(i,j,ispec,i_sls,2) + &
- 2.0d0 * e1_force_RK(i,j,ispec,i_sls,3) + e1_force_RK(i,j,ispec,i_sls,4))
- endif
- endif
-
-
! evolution e11
if(stage_time_scheme == 1) then
Un = b_e11(i,j,ispec,i_sls)
@@ -1345,37 +1245,8 @@
b_e11(i,j,ispec,i_sls) = Unp1
endif
- if(stage_time_scheme == 6) then
- temp_time_scheme = b_dux_dxl_n(i,j,ispec)- b_theta_n/TWO
- temper_time_scheme = phi_nu2(i,j,ispec,i_sls)
- b_e11_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * b_e11_LDDRK(i,j,ispec,i_sls) &
- + b_deltat * (temp_time_scheme * temper_time_scheme)- &
- b_deltat * (b_e11(i,j,ispec,i_sls) * inv_tau_sigma_nu2(i,j,ispec,i_sls))
- b_e11(i,j,ispec,i_sls) = b_e11(i,j,ispec,i_sls)+beta_LDDRK(i_stage)*b_e11_LDDRK(i,j,ispec,i_sls)
- endif
- if(stage_time_scheme == 4) then
- temp_time_scheme = dux_dxl_n(i,j,ispec)-theta_n/TWO
- temper_time_scheme = phi_nu2(i,j,ispec,i_sls)
- e11_force_RK(i,j,ispec,i_sls,i_stage) = deltat * (temp_time_scheme * temper_time_scheme- &
- e11(i,j,ispec,i_sls) * inv_tau_sigma_nu2(i,j,ispec,i_sls))
-
- if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then
- if(i_stage == 1)weight_rk = 0.5d0
- if(i_stage == 2)weight_rk = 0.5d0
- if(i_stage == 3)weight_rk = 1.0d0
- if(i_stage==1)then
- e11_initial_rk(i,j,ispec,i_sls) = e11(i,j,ispec,i_sls)
- endif
- e11(i,j,ispec,i_sls) = e11_initial_rk(i,j,ispec,i_sls) + weight_rk * e11_force_RK(i,j,ispec,i_sls,i_stage)
- elseif(i_stage==4)then
- e11(i,j,ispec,i_sls) = e11_initial_rk(i,j,ispec,i_sls) + 1.0d0 / 6.0d0 * &
- (e11_force_RK(i,j,ispec,i_sls,1) + 2.0d0 * e11_force_RK(i,j,ispec,i_sls,2) + &
- 2.0d0 * e11_force_RK(i,j,ispec,i_sls,3) + e11_force_RK(i,j,ispec,i_sls,4))
- endif
- endif
-
! evolution e13
if(stage_time_scheme == 1) then
Un = b_e13(i,j,ispec,i_sls)
@@ -1392,35 +1263,6 @@
b_e13(i,j,ispec,i_sls) = Unp1
endif
- if(stage_time_scheme == 6) then
- temp_time_scheme=b_dux_dzl_n(i,j,ispec) + b_duz_dxl_n(i,j,ispec)
- temper_time_scheme=phi_nu2(i,j,ispec,i_sls)
- b_e13_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * b_e13_LDDRK(i,j,ispec,i_sls)+&
- b_deltat * (temp_time_scheme*temper_time_scheme)- &
- b_deltat * (b_e13(i,j,ispec,i_sls) * inv_tau_sigma_nu2(i,j,ispec,i_sls))
- b_e13(i,j,ispec,i_sls) = b_e13(i,j,ispec,i_sls)+beta_LDDRK(i_stage) * b_e13_LDDRK(i,j,ispec,i_sls)
- endif
-
- if(stage_time_scheme == 4) then
- temp_time_scheme=dux_dzl_n(i,j,ispec) + duz_dxl_n(i,j,ispec)
- temper_time_scheme=phi_nu2(i,j,ispec,i_sls)
- e13_force_RK(i,j,ispec,i_sls,i_stage) = deltat * (temp_time_scheme*temper_time_scheme- &
- e13(i,j,ispec,i_sls) * inv_tau_sigma_nu2(i,j,ispec,i_sls))
- if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then
- if(i_stage == 1)weight_rk = 0.5d0
- if(i_stage == 2)weight_rk = 0.5d0
- if(i_stage == 3)weight_rk = 1.0d0
- if(i_stage==1)then
- e13_initial_rk(i,j,ispec,i_sls) = e13(i,j,ispec,i_sls)
- endif
- e13(i,j,ispec,i_sls) = e13_initial_rk(i,j,ispec,i_sls) + weight_rk * e13_force_RK(i,j,ispec,i_sls,i_stage)
- elseif(i_stage==4)then
- e13(i,j,ispec,i_sls) = e13_initial_rk(i,j,ispec,i_sls) + 1.0d0 / 6.0d0 * &
- (e13_force_RK(i,j,ispec,i_sls,1) + 2.0d0 * e13_force_RK(i,j,ispec,i_sls,2) + &
- 2.0d0 * e13_force_RK(i,j,ispec,i_sls,3) + e13_force_RK(i,j,ispec,i_sls,4))
- endif
- endif
-
enddo
enddo
Modified: seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/gmat01.f90
===================================================================
--- seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/gmat01.f90 2012-02-16 05:45:45 UTC (rev 19641)
+++ seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/gmat01.f90 2012-02-16 15:18:08 UTC (rev 19642)
@@ -281,6 +281,14 @@
! material can be acoustic (fluid) or elastic (solid)
if(poroelastcoef(2,1,n) > TINYVAL) then ! elastic
write(IOUT,200) n,cp,cs,density(1),poisson,lambda,mu,kappa,young,QKappa,Qmu
+ if(poisson < 0.d0) then
+ write(IOUT,*)
+ write(IOUT,*) 'Materials with a negative Poisson''s ratio can exist,'
+ write(IOUT,*) 'see e.g. R. Lakes, "Science" vol. 235, p. 1038-1040 (1987),'
+ write(IOUT,*) 'but are extremely rare.'
+ write(IOUT,*) 'Hope you know what you are doing...'
+ write(IOUT,*)
+ endif
else ! acoustic
write(IOUT,300) n,cp,density(1),kappa,QKappa,Qmu
endif
@@ -290,6 +298,14 @@
! material is poroelastic (solid/fluid)
write(iout,500) n,sqrt(cpIsquare),sqrt(cpIIsquare),sqrt(cssquare)
write(iout,600) density(1),poisson_s,lambda_s,mu_s,kappa_s,young_s
+ if(poisson_s < 0.d0) then
+ write(IOUT,*)
+ write(IOUT,*) 'Materials with a negative Poisson''s ratio can exist,'
+ write(IOUT,*) 'see e.g. R. Lakes, "Science" vol. 235, p. 1038-1040 (1987),'
+ write(IOUT,*) 'but are extremely rare.'
+ write(IOUT,*) 'Hope you know what you are doing...'
+ write(IOUT,*)
+ endif
write(iout,700) density(2),kappa_f,eta_f
write(iout,800) lambda_fr,mu_fr,kappa_fr,porosity_array(n),tortuosity_array(n),&
permeability(1,n),permeability(2,n),permeability(3,n),Qmu
@@ -313,11 +329,11 @@
'P-wave velocity. . . . . . . . . . . (cp) =',1pe15.8,/5x, &
'S-wave velocity. . . . . . . . . . . (cs) =',1pe15.8,/5x, &
'Mass density. . . . . . . . . . (density) =',1pe15.8,/5x, &
- 'Poisson''s ratio . . . . . . . .(poisson) =',1pe15.8,/5x, &
+ 'Poisson''s ratio. . . . . . . . .(poisson) =',1pe15.8,/5x, &
'First Lame parameter Lambda. . . (lambda) =',1pe15.8,/5x, &
'Second Lame parameter Mu. . . . . . .(mu) =',1pe15.8,/5x, &
'Bulk modulus Kappa . . . . . . . .(kappa) =',1pe15.8,/5x, &
- 'Young''s modulus E . . . . . . . .(young) =',1pe15.8,/5x, &
+ 'Young''s modulus E. . . . . . . . .(young) =',1pe15.8,/5x, &
'QKappa_attenuation . . . . . . .(QKappa) =',1pe15.8,/5x, &
'Qmu_attenuation . . . . . . . . . . (Qmu) =',1pe15.8)
@@ -354,11 +370,11 @@
600 format(//5x,'-------------------------------',/5x, &
'-- Solid phase properties --',/5x, &
'Mass density. . . . . . . . . . (density_s) =',1pe15.8,/5x, &
- 'Poisson''s ratio. . . . . . . . (poisson_s) =',1pe15.8,/5x, &
+ 'Poisson''s ratio . . . . . . . . (poisson_s) =',1pe15.8,/5x, &
'First Lame parameter Lambda. . . (lambda_s) =',1pe15.8,/5x, &
'Second Lame parameter Mu. . . . . . .(mu_s) =',1pe15.8,/5x, &
'Solid bulk modulus Kappa . . . . .(kappa_s) =',1pe15.8,/5x, &
- 'Young''s modulus E. . . . . . . .(young_s) =',1pe15.8)
+ 'Young''s modulus E . . . . . . . .(young_s) =',1pe15.8)
700 format(//5x,'-------------------------------',/5x, &
'-- Fluid phase properties --',/5x, &
@@ -385,7 +401,7 @@
'H. . . . . . . . . .=',1pe15.8,/5x, &
'C. . . . . . . . . .=',1pe15.8,/5x, &
'M. . . . . . . . . .=',1pe15.8,/5x, &
- 'characteristic freq =',1pe15.8)
+ 'Characteristic freq =',1pe15.8)
end subroutine gmat01
Modified: seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/prepare_absorb.f90
===================================================================
--- seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/prepare_absorb.f90 2012-02-16 05:45:45 UTC (rev 19641)
+++ seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/prepare_absorb.f90 2012-02-16 15:18:08 UTC (rev 19642)
@@ -254,67 +254,45 @@
!-------------------------------------------------------------------------------------------------
!
-! subroutine prepare_absorb_elastic(NSTEP,p_sv, & !old
-! nspec_left,nspec_right,nspec_bottom,nspec_top, & !old
-! b_absorb_elastic_left,b_absorb_elastic_right, & !old
-! b_absorb_elastic_bottom,b_absorb_elastic_top) !old
+ subroutine prepare_absorb_elastic(NSTEP,p_sv, &
+ nspec_left,nspec_right,nspec_bottom,nspec_top, &
+ b_absorb_elastic_left,b_absorb_elastic_right, &
+ b_absorb_elastic_bottom,b_absorb_elastic_top)
- subroutine prepare_absorb_elastic(NSTEP,p_sv, & !xiezhinan
- nspec_left,nspec_right,nspec_bottom,nspec_top, & !xiezhinan
- b_absorb_elastic_left,b_absorb_elastic_right, & !xiezhinan
- b_absorb_elastic_bottom,b_absorb_elastic_top,stage_time_scheme)!xiezhinan
-
implicit none
include "constants.h"
logical :: p_sv
integer :: nspec_left,nspec_right,nspec_bottom,nspec_top
integer :: NSTEP
- integer :: stage_time_scheme,i_stage !xiezhinan
+ real(kind=CUSTOM_REAL) :: b_absorb_elastic_left(3,NGLLZ,nspec_left,NSTEP)
+ real(kind=CUSTOM_REAL) :: b_absorb_elastic_right(3,NGLLZ,nspec_right,NSTEP)
+ real(kind=CUSTOM_REAL) :: b_absorb_elastic_bottom(3,NGLLX,nspec_bottom,NSTEP)
+ real(kind=CUSTOM_REAL) :: b_absorb_elastic_top(3,NGLLX,nspec_top,NSTEP)
-! real(kind=CUSTOM_REAL) :: b_absorb_elastic_left(3,NGLLZ,nspec_left,NSTEP) !old
-! real(kind=CUSTOM_REAL) :: b_absorb_elastic_right(3,NGLLZ,nspec_right,NSTEP) !old
-! real(kind=CUSTOM_REAL) :: b_absorb_elastic_bottom(3,NGLLX,nspec_bottom,NSTEP) !old
-! real(kind=CUSTOM_REAL) :: b_absorb_elastic_top(3,NGLLX,nspec_top,NSTEP) !old
-
- real(kind=CUSTOM_REAL) :: b_absorb_elastic_left(3,NGLLZ,nspec_left,NSTEP,stage_time_scheme) !xiezhinan
- real(kind=CUSTOM_REAL) :: b_absorb_elastic_right(3,NGLLZ,nspec_right,NSTEP,stage_time_scheme) !xiezhinan
- real(kind=CUSTOM_REAL) :: b_absorb_elastic_bottom(3,NGLLX,nspec_bottom,NSTEP,stage_time_scheme)!xiezhinan
- real(kind=CUSTOM_REAL) :: b_absorb_elastic_top(3,NGLLX,nspec_top,NSTEP,stage_time_scheme) !xiezhinan
-
! local parameters
- integer :: ispec,i,it
+ integer :: ispec,i,it
do it =1, NSTEP
- do i_stage =1, stage_time_scheme !xiezhinan
-
!--- left absorbing boundary
if(nspec_left >0) then
do ispec = 1,nspec_left
if(p_sv)then!P-SV waves
do i=1,NGLLZ
-! read(35) b_absorb_elastic_left(1,i,ispec,it) !old
- read(35) b_absorb_elastic_left(1,i,ispec,it,i_stage) !xiezhinan
+ read(35) b_absorb_elastic_left(1,i,ispec,it)
enddo
do i=1,NGLLZ
-! read(35) b_absorb_elastic_left(3,i,ispec,it) !old
- read(35) b_absorb_elastic_left(3,i,ispec,it,i_stage) !xiezhinan
+ read(35) b_absorb_elastic_left(3,i,ispec,it)
enddo
-! b_absorb_elastic_left(2,:,ispec,it) = ZERO !old
- b_absorb_elastic_left(2,:,ispec,it,i_stage) = ZERO !xiezhinan
+ b_absorb_elastic_left(2,:,ispec,it) = ZERO
else!SH (membrane) waves
do i=1,NGLLZ
-! read(35) b_absorb_elastic_left(2,i,ispec,it) !old
- read(35) b_absorb_elastic_left(2,i,ispec,it,i_stage) !xiezhinan
+ read(35) b_absorb_elastic_left(2,i,ispec,it)
enddo
-! b_absorb_elastic_left(1,:,ispec,it) = ZERO !old
-! b_absorb_elastic_left(3,:,ispec,it) = ZERO !old
-
- b_absorb_elastic_left(1,:,ispec,it,i_stage) = ZERO !xiezhinan
- b_absorb_elastic_left(3,:,ispec,it,i_stage) = ZERO !xiezhinan
-
+ b_absorb_elastic_left(1,:,ispec,it) = ZERO
+ b_absorb_elastic_left(3,:,ispec,it) = ZERO
endif
enddo
@@ -326,26 +304,18 @@
if(p_sv)then!P-SV waves
do i=1,NGLLZ
-! read(36) b_absorb_elastic_right(1,i,ispec,it) !old
- read(36) b_absorb_elastic_right(1,i,ispec,it,i_stage) !xiezhinan
+ read(36) b_absorb_elastic_right(1,i,ispec,it)
enddo
do i=1,NGLLZ
-! read(36) b_absorb_elastic_right(3,i,ispec,it) !old
- read(36) b_absorb_elastic_right(3,i,ispec,it,i_stage) !xiezhinan
+ read(36) b_absorb_elastic_right(3,i,ispec,it)
enddo
-! b_absorb_elastic_right(2,:,ispec,it) = ZERO !old
- b_absorb_elastic_right(2,:,ispec,it,i_stage) = ZERO !xiezhinan
+ b_absorb_elastic_right(2,:,ispec,it) = ZERO
else!SH (membrane) waves
do i=1,NGLLZ
-! read(36) b_absorb_elastic_right(2,i,ispec,it) !old
- read(36) b_absorb_elastic_right(2,i,ispec,it,i_stage) !xiezhinan
+ read(36) b_absorb_elastic_right(2,i,ispec,it)
enddo
-! b_absorb_elastic_right(1,:,ispec,it) = ZERO !old
-! b_absorb_elastic_right(3,:,ispec,it) = ZERO !old
-
- b_absorb_elastic_right(1,:,ispec,it,i_stage) = ZERO !xiezhinan
- b_absorb_elastic_right(3,:,ispec,it,i_stage) = ZERO !xiezhinan
-
+ b_absorb_elastic_right(1,:,ispec,it) = ZERO
+ b_absorb_elastic_right(3,:,ispec,it) = ZERO
endif
enddo
@@ -357,24 +327,18 @@
if(p_sv)then!P-SV waves
do i=1,NGLLX
-! read(37) b_absorb_elastic_bottom(1,i,ispec,it) !old
- read(37) b_absorb_elastic_bottom(1,i,ispec,it,i_stage) !xiezhinan
+ read(37) b_absorb_elastic_bottom(1,i,ispec,it)
enddo
do i=1,NGLLX
-! read(37) b_absorb_elastic_bottom(3,i,ispec,it) !old
- read(37) b_absorb_elastic_bottom(3,i,ispec,it,i_stage) !xiezhinan
+ read(37) b_absorb_elastic_bottom(3,i,ispec,it)
enddo
-! b_absorb_elastic_bottom(2,:,ispec,it) = ZERO !old
- b_absorb_elastic_bottom(2,:,ispec,it,i_stage) = ZERO !xiezhinan
+ b_absorb_elastic_bottom(2,:,ispec,it) = ZERO
else!SH (membrane) waves
do i=1,NGLLZ
-! read(37) b_absorb_elastic_bottom(2,i,ispec,it) !old
- read(37) b_absorb_elastic_bottom(2,i,ispec,it,i_stage) !xiezhinan
+ read(37) b_absorb_elastic_bottom(2,i,ispec,it)
enddo
-! b_absorb_elastic_bottom(1,:,ispec,it) = ZERO !old
-! b_absorb_elastic_bottom(3,:,ispec,it) = ZERO !old
- b_absorb_elastic_bottom(1,:,ispec,it,i_stage) = ZERO !xiezhinan
- b_absorb_elastic_bottom(3,:,ispec,it,i_stage) = ZERO !xiezhinan
+ b_absorb_elastic_bottom(1,:,ispec,it) = ZERO
+ b_absorb_elastic_bottom(3,:,ispec,it) = ZERO
endif
enddo
@@ -386,31 +350,23 @@
if(p_sv)then!P-SV waves
do i=1,NGLLX
-! read(38) b_absorb_elastic_top(1,i,ispec,it) !old
- read(38) b_absorb_elastic_top(1,i,ispec,it,i_stage) !xiezhinan
+ read(38) b_absorb_elastic_top(1,i,ispec,it)
enddo
do i=1,NGLLX
-! read(38) b_absorb_elastic_top(3,i,ispec,it) !old
- read(38) b_absorb_elastic_top(3,i,ispec,it,i_stage) !xiezhinan
+ read(38) b_absorb_elastic_top(3,i,ispec,it)
enddo
-! b_absorb_elastic_top(2,:,ispec,it) = ZERO !old
- b_absorb_elastic_top(2,:,ispec,it,i_stage) = ZERO !xiezhinan
+ b_absorb_elastic_top(2,:,ispec,it) = ZERO
else!SH (membrane) waves
do i=1,NGLLZ
-! read(38) b_absorb_elastic_top(2,i,ispec,it) !old
- read(38) b_absorb_elastic_top(2,i,ispec,it,i_stage) !xiezhinan
+ read(38) b_absorb_elastic_top(2,i,ispec,it)
enddo
-! b_absorb_elastic_top(1,:,ispec,it) = ZERO !old
-! b_absorb_elastic_top(3,:,ispec,it) = ZERO !old
-
- b_absorb_elastic_top(1,:,ispec,it,i_stage) = ZERO !xiezhinan
- b_absorb_elastic_top(3,:,ispec,it,i_stage) = ZERO !xiezhinan
-
+ b_absorb_elastic_top(1,:,ispec,it) = ZERO
+ b_absorb_elastic_top(3,:,ispec,it) = ZERO
endif
enddo
endif
- enddo !xiezhinan
+
enddo
end subroutine prepare_absorb_elastic
Modified: seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/read_databases.f90
===================================================================
--- seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/read_databases.f90 2012-02-16 05:45:45 UTC (rev 19641)
+++ seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/read_databases.f90 2012-02-16 15:18:08 UTC (rev 19642)
@@ -262,12 +262,12 @@
read(IIN,*) NSTEP,deltat
if (myrank == 0 .and. ipass == 1) write(IOUT,703) NSTEP,deltat,NSTEP*deltat
- if( SIMULATION_TYPE == 1 .and. SAVE_FORWARD .and. &
- (ATTENUATION_VISCOELASTIC_SOLID .or. ATTENUATION_PORO_FLUID_PART) ) then
- ! print*, '*************** WARNING ***************' !xiezhinan
- ! print*, 'Anisotropy & Attenuation & Viscous damping are not presently implemented for adjoint calculations'
- ! stop
- endif
+ ! if( SIMULATION_TYPE == 1 .and. SAVE_FORWARD .and. &
+ ! (ATTENUATION_VISCOELASTIC_SOLID .or. ATTENUATION_PORO_FLUID_PART) ) then
+ ! print*, '*************** WARNING ***************'
+ ! print*, 'Anisotropy & Attenuation & Viscous damping are not presently implemented for adjoint calculations'
+ ! stop
+ ! endif
NTSTEP_BETWEEN_OUTPUT_SEISMO = min(NSTEP,NTSTEP_BETWEEN_OUTPUT_INFO)
@@ -300,7 +300,7 @@
'Add Bielak conditions . . . .(add_Bielak_conditions) = ',l6/5x, &
'Assign external model . . . .(assign_external_model) = ',l6/5x, &
'Read external SEP file . . .(READ_EXTERNAL_SEP_FILE) = ',l6/5x, &
- 'Turn attenuation on or off. . .(ATTENUATION_VISCOELASTIC_SOLID) = ',l6/5x, &
+ 'Attenuation on/off .(ATTENUATION_VISCOELASTIC_SOLID) = ',l6/5x, &
'Save grid in external file or not. . . (output_grid) = ',l6/5x, &
'Save a file with total energy or not.(output_energy) = ',l6)
Modified: seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/specfem2D.F90 2012-02-16 05:45:45 UTC (rev 19641)
+++ seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/specfem2D.F90 2012-02-16 15:18:08 UTC (rev 19642)
@@ -387,8 +387,7 @@
character(len=150) dummystring
! for seismograms
- double precision, dimension(:,:,:), allocatable :: sisux,sisuz,siscurl !xiezhinan
-! double precision, dimension(:,:), allocatable :: sisux,sisuz,siscurl !old
+ double precision, dimension(:,:), allocatable :: sisux,sisuz,siscurl
integer :: seismo_offset, seismo_current
! vector field in an element
@@ -543,11 +542,15 @@
double precision :: f0_attenuation
integer nspec_allocate
double precision :: deltatsquare,deltatcube,deltatfourth,twelvedeltat,fourdeltatsquare
+!xiezhinan
double precision :: b_deltatsquare,b_deltatcube,b_deltatfourth,b_twelvedeltat,b_fourdeltatsquare
+!xiezhinan
+
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: e1,e11,e13
+!xiezhinan
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: b_e1,b_e11,b_e13
+!xiezhinan
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: e1_LDDRK,e11_LDDRK,e13_LDDRK
- real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: b_e1_LDDRK,b_e11_LDDRK,b_e13_LDDRK !xiezhinan
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: e1_initial_rk,e11_initial_rk,e13_initial_rk
real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: e1_force_rk,e11_force_rk,e13_force_rk
double precision, dimension(:,:,:,:), allocatable :: inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2
@@ -557,10 +560,11 @@
real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: &
dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n,dux_dxl_np1,duz_dzl_np1,duz_dxl_np1,dux_dzl_np1
-
+!xiezhinan
real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: &
b_dux_dxl_n,b_duz_dzl_n,b_duz_dxl_n,b_dux_dzl_n,&
- b_dux_dxl_np1,b_duz_dzl_np1,b_duz_dxl_np1,b_dux_dzl_np1
+ b_dux_dxl_np1,b_duz_dzl_np1,b_duz_dxl_np1,b_dux_dzl_np1
+!xiezhinan
! for viscous attenuation
double precision, dimension(:,:,:), allocatable :: &
@@ -640,8 +644,6 @@
b_accelw_poroelastic,b_velocw_poroelastic,b_displw_poroelastic
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: &
b_accel_elastic,b_veloc_elastic,b_displ_elastic
- real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: & !xiezhinan
- b_veloc_elastic_LDDRK,b_displ_elastic_LDDRK !xiezhinan
real(kind=CUSTOM_REAL), dimension(:), allocatable :: &
b_potential_dot_dot_acoustic,b_potential_dot_acoustic,b_potential_acoustic
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: accel_ac,b_displ_ac,b_accel_ac
@@ -664,29 +666,16 @@
character(len=150) :: adj_source_file
integer :: irec_local,nadj_rec_local
double precision :: xx,zz,rholb,tempx1l,tempx2l,b_tempx1l,b_tempx2l,bb_tempx1l,bb_tempx2l
-! real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: adj_sourcearray !old
- real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: adj_sourcearray !xiezhinan
-! real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: adj_sourcearrays !old
- real(kind=CUSTOM_REAL), dimension(:,:,:,:,:,:), allocatable :: adj_sourcearrays !xiezhinan
-
-
-! real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: & !old
-! b_absorb_elastic_left,b_absorb_poro_s_left,b_absorb_poro_w_left !old
-! real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: & !old
-! b_absorb_elastic_right,b_absorb_poro_s_right,b_absorb_poro_w_right !old
-! real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: & !old
-! b_absorb_elastic_bottom,b_absorb_poro_s_bottom,b_absorb_poro_w_bottom !old
-! real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: & !old
-! b_absorb_elastic_top,b_absorb_poro_s_top,b_absorb_poro_w_top !old
-
- real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: & !xiezhinan
- b_absorb_elastic_left,b_absorb_elastic_right,b_absorb_elastic_bottom,b_absorb_elastic_top !xiezhinan
-
- real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: & !xiezhinan
- b_absorb_poro_s_left,b_absorb_poro_w_left,b_absorb_poro_s_right,b_absorb_poro_w_right, & !xiezhinan
- b_absorb_poro_s_bottom,b_absorb_poro_w_bottom,b_absorb_poro_s_top,b_absorb_poro_w_top !xiezhinan
-
-
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: adj_sourcearray
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: adj_sourcearrays
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: &
+ b_absorb_elastic_left,b_absorb_poro_s_left,b_absorb_poro_w_left
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: &
+ b_absorb_elastic_right,b_absorb_poro_s_right,b_absorb_poro_w_right
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: &
+ b_absorb_elastic_bottom,b_absorb_poro_s_bottom,b_absorb_poro_w_bottom
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: &
+ b_absorb_elastic_top,b_absorb_poro_s_top,b_absorb_poro_w_top
real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: &
b_absorb_acoustic_left,b_absorb_acoustic_right,&
b_absorb_acoustic_bottom, b_absorb_acoustic_top
@@ -885,8 +874,7 @@
!<SU_FORMAT
integer(kind=4) :: r4head(60)
character(len=512) :: filename
-! real(kind=4),dimension(:,:),allocatable :: adj_src_s !old
- real(kind=4),dimension(:,:,:),allocatable :: adj_src_s !xiezhinan
+ real(kind=4),dimension(:,:),allocatable :: adj_src_s
integer(kind=2) :: header2(2)
!>SU_FORMAT
@@ -992,8 +980,8 @@
POWER_DISPLAY_COLOR,PERFORM_CUTHILL_MCKEE,SU_FORMAT,USER_T0, time_stepping_scheme, &
ADD_PERIODIC_CONDITIONS,PERIODIC_horiz_dist,PERIODIC_DETECT_TOL)
if(nproc_read_from_database < 1) stop 'should have nproc_read_from_database >= 1'
-! if(SIMULATION_TYPE == 2 .and.(time_stepping_scheme == 2 .or. time_stepping_scheme == 3)) &
-! stop 'RK and LDDRK time scheme not supported for adjoint inversion'
+ if(SIMULATION_TYPE == 2 .and.(time_stepping_scheme == 2 .or. time_stepping_scheme == 3)) &
+ stop 'RK and LDDRK time scheme not supported for adjoint inversion'
if(nproc /= nproc_read_from_database) stop 'must always have nproc == nproc_read_from_database'
!! DK DK Dec 2011: add a crack manually
npgeo_ori = npgeo
@@ -1257,42 +1245,27 @@
e1(:,:,:,:) = 0._CUSTOM_REAL
e11(:,:,:,:) = 0._CUSTOM_REAL
e13(:,:,:,:) = 0._CUSTOM_REAL
- if(SIMULATION_TYPE == 2)then
+!xiezhinan
allocate(b_e1(NGLLX,NGLLZ,nspec_allocate,N_SLS))
allocate(b_e11(NGLLX,NGLLZ,nspec_allocate,N_SLS))
allocate(b_e13(NGLLX,NGLLZ,nspec_allocate,N_SLS))
b_e1(:,:,:,:) = 0._CUSTOM_REAL
b_e11(:,:,:,:) = 0._CUSTOM_REAL
b_e13(:,:,:,:) = 0._CUSTOM_REAL
- endif
+!xiezhinan
if(time_stepping_scheme == 2)then
allocate(e1_LDDRK(NGLLX,NGLLZ,nspec_allocate,N_SLS))
allocate(e11_LDDRK(NGLLX,NGLLZ,nspec_allocate,N_SLS))
allocate(e13_LDDRK(NGLLX,NGLLZ,nspec_allocate,N_SLS))
- if(SIMULATION_TYPE == 2)then
- allocate(b_e1_LDDRK(NGLLX,NGLLZ,nspec_allocate,N_SLS))
- allocate(b_e11_LDDRK(NGLLX,NGLLZ,nspec_allocate,N_SLS))
- allocate(b_e13_LDDRK(NGLLX,NGLLZ,nspec_allocate,N_SLS))
- endif
else
allocate(e1_LDDRK(1,1,1,1))
allocate(e11_LDDRK(1,1,1,1))
allocate(e13_LDDRK(1,1,1,1))
- if(SIMULATION_TYPE == 2)then
- allocate(b_e1_LDDRK(1,1,1,1))
- allocate(b_e11_LDDRK(1,1,1,1))
- allocate(b_e13_LDDRK(1,1,1,1))
endif
- endif
e1_LDDRK(:,:,:,:) = 0._CUSTOM_REAL
e11_LDDRK(:,:,:,:) = 0._CUSTOM_REAL
e13_LDDRK(:,:,:,:) = 0._CUSTOM_REAL
- if(SIMULATION_TYPE == 2)then
- b_e1_LDDRK(:,:,:,:) = 0._CUSTOM_REAL
- b_e11_LDDRK(:,:,:,:) = 0._CUSTOM_REAL
- b_e13_LDDRK(:,:,:,:) = 0._CUSTOM_REAL
- endif
if(time_stepping_scheme == 3)then
allocate(e1_initial_rk(NGLLX,NGLLZ,nspec_allocate,N_SLS))
@@ -1324,7 +1297,7 @@
allocate(duz_dzl_np1(NGLLX,NGLLZ,nspec_allocate))
allocate(duz_dxl_np1(NGLLX,NGLLZ,nspec_allocate))
allocate(dux_dzl_np1(NGLLX,NGLLZ,nspec_allocate))
- if(SIMULATION_TYPE == 2)then
+!xiezhinan
allocate(b_dux_dxl_n(NGLLX,NGLLZ,nspec_allocate))
allocate(b_duz_dzl_n(NGLLX,NGLLZ,nspec_allocate))
allocate(b_duz_dxl_n(NGLLX,NGLLZ,nspec_allocate))
@@ -1333,8 +1306,7 @@
allocate(b_duz_dzl_np1(NGLLX,NGLLZ,nspec_allocate))
allocate(b_duz_dxl_np1(NGLLX,NGLLZ,nspec_allocate))
allocate(b_dux_dzl_np1(NGLLX,NGLLZ,nspec_allocate))
- endif
-
+!xiezhinan
allocate(Mu_nu1(NGLLX,NGLLZ,nspec))
allocate(Mu_nu2(NGLLX,NGLLZ,nspec))
endif
@@ -1482,27 +1454,15 @@
! Files to save absorbed waves needed to reconstruct backward wavefield for adjoint method
if(ipass == 1) then
if(any_elastic .and. (SAVE_FORWARD .or. SIMULATION_TYPE == 2)) then
-! allocate(b_absorb_elastic_left(3,NGLLZ,nspec_left,NSTEP)) !old
-! allocate(b_absorb_elastic_right(3,NGLLZ,nspec_right,NSTEP)) !old
-! allocate(b_absorb_elastic_bottom(3,NGLLX,nspec_bottom,NSTEP)) !old
-! allocate(b_absorb_elastic_top(3,NGLLX,nspec_top,NSTEP)) !old
-
- allocate(b_absorb_elastic_left(3,NGLLZ,nspec_left,NSTEP,stage_time_scheme)) !xiezhinan
- allocate(b_absorb_elastic_right(3,NGLLZ,nspec_right,NSTEP,stage_time_scheme)) !xiezhinan
- allocate(b_absorb_elastic_bottom(3,NGLLX,nspec_bottom,NSTEP,stage_time_scheme)) !xiezhinan
- allocate(b_absorb_elastic_top(3,NGLLX,nspec_top,NSTEP,stage_time_scheme)) !xiezhinan
-
+ allocate(b_absorb_elastic_left(3,NGLLZ,nspec_left,NSTEP))
+ allocate(b_absorb_elastic_right(3,NGLLZ,nspec_right,NSTEP))
+ allocate(b_absorb_elastic_bottom(3,NGLLX,nspec_bottom,NSTEP))
+ allocate(b_absorb_elastic_top(3,NGLLX,nspec_top,NSTEP))
else
-! allocate(b_absorb_elastic_left(1,1,1,1)) !old
-! allocate(b_absorb_elastic_right(1,1,1,1)) !old
-! allocate(b_absorb_elastic_bottom(1,1,1,1)) !old
-! allocate(b_absorb_elastic_top(1,1,1,1)) !old
-
- allocate(b_absorb_elastic_left(1,1,1,1,1)) !xiezhinan
- allocate(b_absorb_elastic_right(1,1,1,1,1)) !xiezhinan
- allocate(b_absorb_elastic_bottom(1,1,1,1,1)) !xiezhinan
- allocate(b_absorb_elastic_top(1,1,1,1,1)) !xiezhinan
-
+ allocate(b_absorb_elastic_left(1,1,1,1))
+ allocate(b_absorb_elastic_right(1,1,1,1))
+ allocate(b_absorb_elastic_bottom(1,1,1,1))
+ allocate(b_absorb_elastic_top(1,1,1,1))
endif
if(any_poroelastic .and. (SAVE_FORWARD .or. SIMULATION_TYPE == 2)) then
allocate(b_absorb_poro_s_left(NDIM,NGLLZ,nspec_left,NSTEP))
@@ -1539,15 +1499,10 @@
else
if(.not. allocated(b_absorb_elastic_left)) then
-! allocate(b_absorb_elastic_left(1,1,1,1)) !old
-! allocate(b_absorb_elastic_right(1,1,1,1)) !old
-! allocate(b_absorb_elastic_bottom(1,1,1,1)) !old
-! allocate(b_absorb_elastic_top(1,1,1,1)) !old
-
- allocate(b_absorb_elastic_left(1,1,1,1,1)) !xiezhinan
- allocate(b_absorb_elastic_right(1,1,1,1,1)) !xiezhinan
- allocate(b_absorb_elastic_bottom(1,1,1,1,1)) !xiezhinan
- allocate(b_absorb_elastic_top(1,1,1,1,1)) !xiezhinan
+ allocate(b_absorb_elastic_left(1,1,1,1))
+ allocate(b_absorb_elastic_right(1,1,1,1))
+ allocate(b_absorb_elastic_bottom(1,1,1,1))
+ allocate(b_absorb_elastic_top(1,1,1,1))
endif
if(.not. allocated(b_absorb_poro_s_left)) then
@@ -2160,14 +2115,11 @@
nadj_rec_local = nadj_rec_local + 1
endif
enddo
-! if(ipass == 1) allocate(adj_sourcearray(NSTEP,3,NGLLX,NGLLZ)) !old
- if(ipass == 1) allocate(adj_sourcearray(NSTEP,stage_time_scheme,3,NGLLX,NGLLZ)) !xiezhinan
+ if(ipass == 1) allocate(adj_sourcearray(NSTEP,3,NGLLX,NGLLZ))
if (nadj_rec_local > 0 .and. ipass == 1) then
-! allocate(adj_sourcearrays(nadj_rec_local,NSTEP,3,NGLLX,NGLLZ)) !old
- allocate(adj_sourcearrays(nadj_rec_local,NSTEP,stage_time_scheme,3,NGLLX,NGLLZ)) !xiezhinan
+ allocate(adj_sourcearrays(nadj_rec_local,NSTEP,3,NGLLX,NGLLZ))
else if (ipass == 1) then
-! allocate(adj_sourcearrays(1,1,1,1,1)) !old
- allocate(adj_sourcearrays(1,1,1,1,1,1)) !xiezhinan
+ allocate(adj_sourcearrays(1,1,1,1,1))
endif
if (.not. SU_FORMAT) then
@@ -2177,42 +2129,30 @@
if(myrank == which_proc_receiver(irec))then
irec_local = irec_local + 1
adj_source_file = trim(station_name(irec))//'.'//trim(network_name(irec))
-! call compute_arrays_adj_source(adj_source_file, & !old
-! xi_receiver(irec), gamma_receiver(irec), & !old
-! adj_sourcearray, xigll,zigll,NSTEP) !old
-
- call compute_arrays_adj_source(adj_source_file, & !xiezhinan
- xi_receiver(irec), gamma_receiver(irec), & !xiezhinan
- adj_sourcearray, xigll,zigll,NSTEP,stage_time_scheme) !xiezhinan
-
-! adj_sourcearrays(irec_local,:,:,:,:) = adj_sourcearray(:,:,:,:) !old
- adj_sourcearrays(irec_local,:,:,:,:,:) = adj_sourcearray(:,:,:,:,:) !xiezhinan
+ call compute_arrays_adj_source(adj_source_file, &
+ xi_receiver(irec), gamma_receiver(irec), &
+ adj_sourcearray, xigll,zigll,NSTEP)
+ adj_sourcearrays(irec_local,:,:,:,:) = adj_sourcearray(:,:,:,:)
endif
enddo
else
irec_local = 0
write(filename, "('./SEM/Ux_file_single.bin.adj')")
-! open(111,file=trim(filename),access='direct',recl=240+4*NSTEP,iostat = ios) !old
- open(111,file=trim(filename),access='direct',recl=240+4*NSTEP*stage_time_scheme,iostat = ios) !xiezhinan
+ open(111,file=trim(filename),access='direct',recl=240+4*NSTEP,iostat = ios)
if (ios /= 0) call exit_MPI(' file '//trim(filename)//'does not exist')
write(filename, "('./SEM/Uz_file_single.bin.adj')")
-! open(113,file=trim(filename),access='direct',recl=240+4*NSTEP,iostat = ios) !old
- open(113,file=trim(filename),access='direct',recl=240+4*NSTEP*stage_time_scheme,iostat = ios) !xiezhinan
+ open(113,file=trim(filename),access='direct',recl=240+4*NSTEP,iostat = ios)
if (ios /= 0) call exit_MPI(' file '//trim(filename)//'does not exist')
-! allocate(adj_src_s(NSTEP,3)) !old
- allocate(adj_src_s(NSTEP,stage_time_scheme,3)) !xiezhinan
+ allocate(adj_src_s(NSTEP,3))
do irec = 1, nrec
if(myrank == which_proc_receiver(irec))then
irec_local = irec_local + 1
-! adj_sourcearray(:,:,:,:) = 0.0 !old
- adj_sourcearray(:,:,:,:,:) = 0.0 !xiezhinan
-! read(111,rec=irec,iostat=ios) r4head, adj_src_s(:,1) !old
- read(111,rec=irec,iostat=ios) r4head, adj_src_s(:,:,1)!xiezhinan
+ adj_sourcearray(:,:,:,:) = 0.0
+ read(111,rec=irec,iostat=ios) r4head, adj_src_s(:,1)
if (ios /= 0) call exit_MPI(' file '//trim(filename)//' read error')
-! read(113,rec=irec,iostat=ios) r4head, adj_src_s(:,3) !old
- read(113,rec=irec,iostat=ios) r4head, adj_src_s(:,:,3)!xiezhinan
+ read(113,rec=irec,iostat=ios) r4head, adj_src_s(:,3)
if (ios /= 0) call exit_MPI(' file '//trim(filename)//' read error')
header2=r4head(29)
if (irec==1) print*, r4head(1),r4head(19),r4head(20),r4head(21),r4head(22),header2(2)
@@ -2220,12 +2160,10 @@
call lagrange_any(gamma_receiver(irec),NGLLZ,zigll,hgammar,hpgammar)
do k = 1, NGLLZ
do i = 1, NGLLX
-! adj_sourcearray(:,:,i,k) = hxir(i) * hgammar(k) * adj_src_s(:,:) !old
- adj_sourcearray(:,:,:,i,k) = hxir(i) * hgammar(k) * adj_src_s(:,:,:) !xiezhinan
+ adj_sourcearray(:,:,i,k) = hxir(i) * hgammar(k) * adj_src_s(:,:)
enddo
enddo
-! adj_sourcearrays(irec_local,:,:,:,:) = adj_sourcearray(:,:,:,:) !old
- adj_sourcearrays(irec_local,:,:,:,:,:) = adj_sourcearray(:,:,:,:,:) !xiezhinan
+ adj_sourcearrays(irec_local,:,:,:,:) = adj_sourcearray(:,:,:,:)
endif ! if(myrank == which_proc_receiver(irec))
enddo ! irec
close(111)
@@ -2233,8 +2171,7 @@
deallocate(adj_src_s)
endif
else if (ipass == 1) then
-! allocate(adj_sourcearrays(1,1,1,1,1)) !old
- allocate(adj_sourcearrays(1,1,1,1,1,1)) !xiezhinan
+ allocate(adj_sourcearrays(1,1,1,1,1))
endif
if (ipass == 1) then
@@ -2458,13 +2395,9 @@
! allocate seismogram arrays
if(ipass == 1) then
- allocate(sisux(NTSTEP_BETWEEN_OUTPUT_SEISMO,stage_time_scheme,nrecloc)) !xiezhinan
- allocate(sisuz(NTSTEP_BETWEEN_OUTPUT_SEISMO,stage_time_scheme,nrecloc)) !xiezhinan
- allocate(siscurl(NTSTEP_BETWEEN_OUTPUT_SEISMO,stage_time_scheme,nrecloc)) !xiezhinan
-
-! allocate(sisux(NTSTEP_BETWEEN_OUTPUT_SEISMO,nrecloc)) !old
-! allocate(sisuz(NTSTEP_BETWEEN_OUTPUT_SEISMO,nrecloc)) !old
-! allocate(siscurl(NTSTEP_BETWEEN_OUTPUT_SEISMO,nrecloc)) !old
+ allocate(sisux(NTSTEP_BETWEEN_OUTPUT_SEISMO,nrecloc))
+ allocate(sisuz(NTSTEP_BETWEEN_OUTPUT_SEISMO,nrecloc))
+ allocate(siscurl(NTSTEP_BETWEEN_OUTPUT_SEISMO,nrecloc))
endif
! check if acoustic receiver is exactly on the free surface because pressure is zero there
@@ -2556,12 +2489,6 @@
allocate(b_displ_elastic(3,nglob))
allocate(b_veloc_elastic(3,nglob))
allocate(b_accel_elastic(3,nglob))
-
- if(time_stepping_scheme == 2)then !xiezhinan
- allocate(b_veloc_elastic_LDDRK(3,nglob)) !xiezhinan
- allocate(b_displ_elastic_LDDRK(3,nglob)) !xiezhinan
- endif !xiezhinan
-
allocate(rho_kl(NGLLX,NGLLZ,nspec))
allocate(rho_k(nglob))
allocate(rhol_global(nglob))
@@ -2582,12 +2509,6 @@
allocate(b_displ_elastic(1,1))
allocate(b_veloc_elastic(1,1))
allocate(b_accel_elastic(1,1))
-
- if(time_stepping_scheme == 2)then !xiezhinan
- allocate(b_veloc_elastic_LDDRK(1,1)) !xiezhinan
- allocate(b_displ_elastic_LDDRK(1,1)) !xiezhinan
- endif !xiezhinan
-
allocate(rho_kl(1,1,1))
allocate(rho_k(1))
allocate(rhol_global(1))
@@ -3269,16 +3190,11 @@
! reads in absorbing boundary data
if(any_elastic) then
-! call prepare_absorb_elastic(NSTEP,p_sv, & !old
-! nspec_left,nspec_right,nspec_bottom,nspec_top, & !old
-! b_absorb_elastic_left,b_absorb_elastic_right, & !old
-! b_absorb_elastic_bottom,b_absorb_elastic_top) !old
+ call prepare_absorb_elastic(NSTEP,p_sv, &
+ nspec_left,nspec_right,nspec_bottom,nspec_top, &
+ b_absorb_elastic_left,b_absorb_elastic_right, &
+ b_absorb_elastic_bottom,b_absorb_elastic_top)
- call prepare_absorb_elastic(NSTEP,p_sv, & !xiezhinan
- nspec_left,nspec_right,nspec_bottom,nspec_top, & !xiezhinan
- b_absorb_elastic_left,b_absorb_elastic_right, & !xiezhinan
- b_absorb_elastic_bottom,b_absorb_elastic_top,stage_time_scheme)!xiezhinan
-
endif
if(any_poroelastic) then
call prepare_absorb_poroelastic(NSTEP, &
@@ -3474,23 +3390,21 @@
endif ! initialfield
- if(SIMULATION_TYPE == 1)then
deltatsquare = deltat * deltat
deltatcube = deltatsquare * deltat
deltatfourth = deltatsquare * deltatsquare
twelvedeltat = 12.d0 * deltat
fourdeltatsquare = 4.d0 * deltatsquare
- endif
- if(SIMULATION_TYPE == 2)then
+!xiezhinan
b_deltatsquare = b_deltat * b_deltat
b_deltatcube = b_deltatsquare * b_deltat
b_deltatfourth = b_deltatsquare * b_deltatsquare
b_twelvedeltat = 12.d0 * b_deltat
b_fourdeltatsquare = 4.d0 * b_deltatsquare
- endif
+!xiezhinan
! compute the source time function and store it in a text file
if(.not. initialfield) then
@@ -3885,7 +3799,7 @@
if(ATTENUATION_PORO_FLUID_PART .and. any_poroelastic .and. &
(time_stepping_scheme == 2.or. time_stepping_scheme == 3)) &
stop 'RK and LDDRK time scheme not supported poroelastic simulations with attenuation'
-
+
if(coupled_elastic_poro) then
if(ATTENUATION_VISCOELASTIC_SOLID .or. ATTENUATION_PORO_FLUID_PART) &
@@ -4463,13 +4377,10 @@
seismo_current = seismo_current + 1
! compute current time
-! time = (it-1)*deltat
+ time = (it-1)*deltat
do i_stage=1, stage_time_scheme
- if(time_stepping_scheme == 1)time = (it-1)*deltat
- if(time_stepping_scheme == 2)time = (it-1)*deltat + c_LDDRK(i_stage)
-
! update displacement using finite-difference time scheme (Newmark)
if(any_elastic) then
@@ -4493,19 +4404,12 @@
endif
if(SIMULATION_TYPE == 2) then ! Adjoint calculation
- if(time_stepping_scheme==1)then
b_displ_elastic = b_displ_elastic &
+ b_deltat*b_veloc_elastic &
+ b_deltatsquareover2*b_accel_elastic
b_veloc_elastic = b_veloc_elastic + b_deltatover2*b_accel_elastic
b_accel_elastic = ZERO
endif
-
- if(time_stepping_scheme==2)then
- b_accel_elastic = ZERO
- endif
-
- endif
endif
if(any_poroelastic) then
@@ -5070,10 +4974,8 @@
do j=1,NGLLZ
do i=1,NGLLX
iglob = ibool(i,j,ispec_selected_rec(irec))
-! potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) & !old
-! + adj_sourcearrays(irec_local,NSTEP-it+1,1,i,j) !old
- potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) & !xiezhinan
- + adj_sourcearrays(irec_local,NSTEP-it+1,stage_time_scheme-i_stage+1,1,i,j) !xiezhinan
+ potential_dot_dot_acoustic(iglob) = potential_dot_dot_acoustic(iglob) &
+ + adj_sourcearrays(irec_local,NSTEP-it+1,1,i,j)
enddo
enddo
endif ! if element acoustic
@@ -5239,10 +5141,10 @@
e1_LDDRK,e11_LDDRK,e13_LDDRK,alpha_LDDRK,beta_LDDRK, &
e1_initial_rk,e11_initial_rk,e13_initial_rk,e1_force_rk, e11_force_rk, e13_force_rk, &
stage_time_scheme,i_stage,&
- b_veloc_elastic,b_e1,b_e11,b_e13,&
- b_deltat,b_deltatcube,b_deltatfourth,b_twelvedeltat,b_fourdeltatsquare,&
- b_dux_dxl_n,b_duz_dzl_n,b_duz_dxl_n,b_dux_dzl_n,& !xiezhinan
- b_dux_dxl_np1,b_duz_dzl_np1,b_duz_dxl_np1,b_dux_dzl_np1,b_e1_LDDRK,b_e11_LDDRK,b_e13_LDDRK)
+ b_veloc_elastic,b_e1,b_e11,b_e13,&
+ b_deltat,b_deltatcube,b_deltatfourth,b_twelvedeltat,b_fourdeltatsquare,&
+ b_dux_dxl_n,b_duz_dzl_n,b_duz_dxl_n,b_dux_dzl_n,&
+ b_dux_dxl_np1,b_duz_dzl_np1,b_duz_dxl_np1,b_dux_dzl_np1)
if(anyabs .and. SAVE_FORWARD .and. SIMULATION_TYPE == 1) then
!--- left absorbing boundary
@@ -5250,17 +5152,14 @@
do ispec = 1,nspec_left
if(p_sv)then!P-SV waves
do i=1,NGLLZ
-! write(35) b_absorb_elastic_left(1,i,ispec,it) !old
- write(35) b_absorb_elastic_left(1,i,ispec,it,i_stage) !xiezhinan
+ write(35) b_absorb_elastic_left(1,i,ispec,it)
enddo
do i=1,NGLLZ
-! write(35) b_absorb_elastic_left(3,i,ispec,it) !old
- write(35) b_absorb_elastic_left(3,i,ispec,it,i_stage) !xiezhinan
+ write(35) b_absorb_elastic_left(3,i,ispec,it)
enddo
else!SH (membrane) waves
do i=1,NGLLZ
-! write(35) b_absorb_elastic_left(2,i,ispec,it) !old
- write(35) b_absorb_elastic_left(2,i,ispec,it,i_stage) !xiezhinan
+ write(35) b_absorb_elastic_left(2,i,ispec,it)
enddo
endif
enddo
@@ -5271,17 +5170,14 @@
do ispec = 1,nspec_right
if(p_sv)then!P-SV waves
do i=1,NGLLZ
-! write(36) b_absorb_elastic_right(1,i,ispec,it) !old
- write(36) b_absorb_elastic_right(1,i,ispec,it,i_stage) !xiezhinan
+ write(36) b_absorb_elastic_right(1,i,ispec,it)
enddo
do i=1,NGLLZ
-! write(36) b_absorb_elastic_right(3,i,ispec,it) !old
- write(36) b_absorb_elastic_right(3,i,ispec,it,i_stage) !xiezhinan
+ write(36) b_absorb_elastic_right(3,i,ispec,it)
enddo
else!SH (membrane) waves
do i=1,NGLLZ
-! write(36) b_absorb_elastic_right(2,i,ispec,it) !old
- write(36) b_absorb_elastic_right(2,i,ispec,it,i_stage) !xiezhinan
+ write(36) b_absorb_elastic_right(2,i,ispec,it)
enddo
endif
enddo
@@ -5292,18 +5188,14 @@
do ispec = 1,nspec_bottom
if(p_sv)then!P-SV waves
do i=1,NGLLX
-! write(37) b_absorb_elastic_bottom(1,i,ispec,it) !old
- write(37) b_absorb_elastic_bottom(1,i,ispec,it,i_stage) !xiezhinan
+ write(37) b_absorb_elastic_bottom(1,i,ispec,it)
enddo
do i=1,NGLLX
-! write(37) b_absorb_elastic_bottom(3,i,ispec,it) !old
- write(37) b_absorb_elastic_bottom(3,i,ispec,it,i_stage) !xiezhinan
-
+ write(37) b_absorb_elastic_bottom(3,i,ispec,it)
enddo
else!SH (membrane) waves
do i=1,NGLLX
-! write(37) b_absorb_elastic_bottom(2,i,ispec,it) !old
- write(37) b_absorb_elastic_bottom(2,i,ispec,it,i_stage) !xiezhinan
+ write(37) b_absorb_elastic_bottom(2,i,ispec,it)
enddo
endif
enddo
@@ -5314,17 +5206,14 @@
do ispec = 1,nspec_top
if(p_sv)then!P-SV waves
do i=1,NGLLX
-! write(38) b_absorb_elastic_top(1,i,ispec,it) !old
- write(38) b_absorb_elastic_top(1,i,ispec,it,i_stage)!xiezhinan
+ write(38) b_absorb_elastic_top(1,i,ispec,it)
enddo
do i=1,NGLLX
-! write(38) b_absorb_elastic_top(3,i,ispec,it) !old
- write(38) b_absorb_elastic_top(3,i,ispec,it,i_stage)!xiezhinan
+ write(38) b_absorb_elastic_top(3,i,ispec,it)
enddo
else!SH (membrane) waves
do i=1,NGLLX
-! write(38) b_absorb_elastic_top(2,i,ispec,it) !old
- write(38) b_absorb_elastic_top(2,i,ispec,it,i_stage) !xiezhinan
+ write(38) b_absorb_elastic_top(2,i,ispec,it)
enddo
endif
enddo
@@ -5931,8 +5820,6 @@
endif
if(SIMULATION_TYPE == 2) then
-
- if(time_stepping_scheme == 1)then
b_accel_elastic(1,:) = b_accel_elastic(1,:) * rmass_inverse_elastic(:)
b_accel_elastic(2,:) = b_accel_elastic(2,:) * rmass_inverse_elastic(:)
b_accel_elastic(3,:) = b_accel_elastic(3,:) * rmass_inverse_elastic(:)
@@ -5940,21 +5827,6 @@
b_veloc_elastic = b_veloc_elastic + b_deltatover2*b_accel_elastic
endif
- if(time_stepping_scheme == 2)then
- b_accel_elastic(1,:) = b_accel_elastic(1,:) * rmass_inverse_elastic
- b_accel_elastic(2,:) = b_accel_elastic(2,:) * rmass_inverse_elastic
- b_accel_elastic(3,:) = b_accel_elastic(3,:) * rmass_inverse_elastic
-
- b_veloc_elastic_LDDRK = alpha_LDDRK(i_stage) * b_veloc_elastic_LDDRK - deltat * b_accel_elastic
- b_displ_elastic_LDDRK = alpha_LDDRK(i_stage) * b_displ_elastic_LDDRK - deltat * b_veloc_elastic
-
- b_veloc_elastic = b_veloc_elastic + beta_LDDRK(i_stage) * b_veloc_elastic_LDDRK
- b_displ_elastic = b_displ_elastic + beta_LDDRK(i_stage) * b_displ_elastic_LDDRK
-
- endif
-
- endif
-
endif !if(any_elastic)
@@ -6735,6 +6607,42 @@
!*******************************************************************************
! assembling the displacements on the elastic-poro boundaries
!*******************************************************************************
+
+!
+! Explanation of the code below, from Christina Morency and Yang Luo, January 2012:
+!
+! Coupled elastic-poroelastic simulations imply continuity of traction and
+! displacement at the interface.
+! For the traction we pass on both sides n*(T + Te)/2 , that is, the average
+! between the total stress (from the poroelastic part) and the elastic stress.
+! For the displacement, we enforce its continuity in the assembling stage,
+! realizing that continuity of displacement correspond to the continuity of
+! the acceleration we have:
+!
+! accel_elastic = rmass_inverse_elastic * force_elastic
+! accels_poroelastic = rmass_s_inverse_poroelastic * force_poroelastic
+!
+! Therefore, continuity of acceleration gives
+!
+! accel = (force_elastic + force_poroelastic)/
+! (1/rmass_inverse_elastic + 1/rmass_inverse_poroelastic)
+!
+! Then
+!
+! accel_elastic = accel
+! accels_poroelastic = accel
+! accelw_poroelastic = 0
+!
+! From there, the velocity and displacement are updated.
+! Note that force_elastic and force_poroelastic are the right hand sides of
+! the equations we solve, that is, the acceleration terms before the
+! division by the inverse of the mass matrices. This is why in the code below
+! we first need to recover the accelerations (which are then
+! the right hand sides forces) and the velocities before the update.
+!
+! This implementation highly helped stability especially with unstructured meshes.
+!
+
if(coupled_elastic_poro) then
icount(:)=ZERO
@@ -7042,7 +6950,7 @@
enddo
endif
- ! enddo !LDDRK or RK !initially set the end of LDDRK
+ enddo !LDDRK or RK
! ********************************************************************************************
! reading lastframe for adjoint/kernels calculation
@@ -7141,8 +7049,68 @@
endif
+
+
!>NOISE_TOMOGRAPHY
+! ********************************************************************************************
+! kernels calculation
+! ********************************************************************************************
+ if(any_elastic .and. SIMULATION_TYPE == 2) then ! kernels calculation
+ do iglob = 1,nglob
+ rho_k(iglob) = accel_elastic(1,iglob)*b_displ_elastic(1,iglob) +&
+ accel_elastic(2,iglob)*b_displ_elastic(2,iglob) +&
+ accel_elastic(3,iglob)*b_displ_elastic(3,iglob)
+ rhorho_el_hessian_temp1(iglob) = accel_elastic(1,iglob)*accel_elastic(1,iglob) +&
+ accel_elastic(2,iglob)*accel_elastic(2,iglob) +&
+ accel_elastic(3,iglob)*accel_elastic(3,iglob)
+ rhorho_el_hessian_temp2(iglob) = accel_elastic(1,iglob)*b_accel_elastic(1,iglob) +&
+ accel_elastic(2,iglob)*b_accel_elastic(2,iglob) +&
+ accel_elastic(3,iglob)*b_accel_elastic(3,iglob)
+ enddo
+ endif
+
+ if(any_poroelastic .and. SIMULATION_TYPE ==2) then
+ do iglob =1,nglob
+ rhot_k(iglob) = accels_poroelastic(1,iglob) * b_displs_poroelastic(1,iglob) + &
+ accels_poroelastic(2,iglob) * b_displs_poroelastic(2,iglob)
+ rhof_k(iglob) = accelw_poroelastic(1,iglob) * b_displs_poroelastic(1,iglob) + &
+ accelw_poroelastic(2,iglob) * b_displs_poroelastic(2,iglob) + &
+ accels_poroelastic(1,iglob) * b_displw_poroelastic(1,iglob) + &
+ accels_poroelastic(2,iglob) * b_displw_poroelastic(2,iglob)
+ sm_k(iglob) = accelw_poroelastic(1,iglob) * b_displw_poroelastic(1,iglob) + &
+ accelw_poroelastic(2,iglob) * b_displw_poroelastic(2,iglob)
+ eta_k(iglob) = velocw_poroelastic(1,iglob) * b_displw_poroelastic(1,iglob) + &
+ velocw_poroelastic(2,iglob) * b_displw_poroelastic(2,iglob)
+ enddo
+ endif
+
+!---- compute kinetic and potential energy
+ if(output_energy) &
+ call compute_energy(displ_elastic,veloc_elastic, &
+ displs_poroelastic,velocs_poroelastic, &
+ displw_poroelastic,velocw_poroelastic, &
+ xix,xiz,gammax,gammaz,jacobian,ibool,elastic,poroelastic,hprime_xx,hprime_zz, &
+ nspec,nglob_acoustic,nglob_elastic,nglob_poroelastic, &
+ assign_external_model,it,deltat,t0,kmato,poroelastcoef,density, &
+ porosity,tortuosity, &
+ vpext,vsext,rhoext,c11ext,c13ext,c15ext,c33ext,c35ext,c55ext, &
+ anisotropic,anisotropy,wxgll,wzgll,numat, &
+ pressure_element,vector_field_element,e1,e11, &
+ potential_dot_acoustic,potential_dot_dot_acoustic, &
+ ATTENUATION_VISCOELASTIC_SOLID,Mu_nu1,Mu_nu2,N_SLS,p_sv)
+
+!---- display time step and max of norm of displacement
+ if(mod(it,NTSTEP_BETWEEN_OUTPUT_INFO) == 0 .or. it == 5 .or. it == NSTEP) then
+ call check_stability(myrank,time,it,NSTEP,NOISE_TOMOGRAPHY, &
+ nglob_acoustic,nglob_elastic,nglob_poroelastic, &
+ any_elastic_glob,any_elastic,displ_elastic, &
+ any_poroelastic_glob,any_poroelastic, &
+ displs_poroelastic,displw_poroelastic, &
+ any_acoustic_glob,any_acoustic,potential_acoustic, &
+ timestamp_seconds_start)
+ endif
+
!---- loop on all the receivers to compute and store the seismograms
do irecloc = 1,nrecloc
@@ -7278,80 +7246,21 @@
! rotate seismogram components if needed, except if recording pressure, which is a scalar
if(seismotype /= 4 .and. seismotype /= 6) then
if(p_sv) then
- sisux(seismo_current,i_stage,irecloc) = cosrot_irec(irecloc)*valux + sinrot_irec(irecloc)*valuz !xiezhinan
- sisuz(seismo_current,i_stage,irecloc) = - sinrot_irec(irecloc)*valux + cosrot_irec(irecloc)*valuz !xiezhinan
+ sisux(seismo_current,irecloc) = cosrot_irec(irecloc)*valux + sinrot_irec(irecloc)*valuz
+ sisuz(seismo_current,irecloc) = - sinrot_irec(irecloc)*valux + cosrot_irec(irecloc)*valuz
else
- sisux(seismo_current,i_stage,irecloc) = valuy !xiezhinan
- sisuz(seismo_current,i_stage,irecloc) = ZERO !xiezhinan
+ sisux(seismo_current,irecloc) = valuy
+ sisuz(seismo_current,irecloc) = ZERO
endif
else
- sisux(seismo_current,i_stage,irecloc) = valux !xiezhinan
- sisuz(seismo_current,i_stage,irecloc) = ZERO !xiezhinan
+ sisux(seismo_current,irecloc) = valux
+ sisuz(seismo_current,irecloc) = ZERO
endif
- siscurl(seismo_current,i_stage,irecloc) = valcurl !xiezhinan
+ siscurl(seismo_current,irecloc) = valcurl
enddo
-enddo !LDDRK or RK
-! ********************************************************************************************
-! kernels calculation
-! ********************************************************************************************
- if(any_elastic .and. SIMULATION_TYPE == 2) then ! kernels calculation
- do iglob = 1,nglob
- rho_k(iglob) = accel_elastic(1,iglob)*b_displ_elastic(1,iglob) +&
- accel_elastic(2,iglob)*b_displ_elastic(2,iglob) +&
- accel_elastic(3,iglob)*b_displ_elastic(3,iglob)
- rhorho_el_hessian_temp1(iglob) = accel_elastic(1,iglob)*accel_elastic(1,iglob) +&
- accel_elastic(2,iglob)*accel_elastic(2,iglob) +&
- accel_elastic(3,iglob)*accel_elastic(3,iglob)
- rhorho_el_hessian_temp2(iglob) = accel_elastic(1,iglob)*b_accel_elastic(1,iglob) +&
- accel_elastic(2,iglob)*b_accel_elastic(2,iglob) +&
- accel_elastic(3,iglob)*b_accel_elastic(3,iglob)
- enddo
- endif
-
- if(any_poroelastic .and. SIMULATION_TYPE ==2) then
- do iglob =1,nglob
- rhot_k(iglob) = accels_poroelastic(1,iglob) * b_displs_poroelastic(1,iglob) + &
- accels_poroelastic(2,iglob) * b_displs_poroelastic(2,iglob)
- rhof_k(iglob) = accelw_poroelastic(1,iglob) * b_displs_poroelastic(1,iglob) + &
- accelw_poroelastic(2,iglob) * b_displs_poroelastic(2,iglob) + &
- accels_poroelastic(1,iglob) * b_displw_poroelastic(1,iglob) + &
- accels_poroelastic(2,iglob) * b_displw_poroelastic(2,iglob)
- sm_k(iglob) = accelw_poroelastic(1,iglob) * b_displw_poroelastic(1,iglob) + &
- accelw_poroelastic(2,iglob) * b_displw_poroelastic(2,iglob)
- eta_k(iglob) = velocw_poroelastic(1,iglob) * b_displw_poroelastic(1,iglob) + &
- velocw_poroelastic(2,iglob) * b_displw_poroelastic(2,iglob)
- enddo
- endif
-
-!---- compute kinetic and potential energy
- if(output_energy) &
- call compute_energy(displ_elastic,veloc_elastic, &
- displs_poroelastic,velocs_poroelastic, &
- displw_poroelastic,velocw_poroelastic, &
- xix,xiz,gammax,gammaz,jacobian,ibool,elastic,poroelastic,hprime_xx,hprime_zz, &
- nspec,nglob_acoustic,nglob_elastic,nglob_poroelastic, &
- assign_external_model,it,deltat,t0,kmato,poroelastcoef,density, &
- porosity,tortuosity, &
- vpext,vsext,rhoext,c11ext,c13ext,c15ext,c33ext,c35ext,c55ext, &
- anisotropic,anisotropy,wxgll,wzgll,numat, &
- pressure_element,vector_field_element,e1,e11, &
- potential_dot_acoustic,potential_dot_dot_acoustic, &
- ATTENUATION_VISCOELASTIC_SOLID,Mu_nu1,Mu_nu2,N_SLS,p_sv)
-
-!---- display time step and max of norm of displacement
- if(mod(it,NTSTEP_BETWEEN_OUTPUT_INFO) == 0 .or. it == 5 .or. it == NSTEP) then
- call check_stability(myrank,time,it,NSTEP,NOISE_TOMOGRAPHY, &
- nglob_acoustic,nglob_elastic,nglob_poroelastic, &
- any_elastic_glob,any_elastic,displ_elastic, &
- any_poroelastic_glob,any_poroelastic, &
- displs_poroelastic,displw_poroelastic, &
- any_acoustic_glob,any_acoustic,potential_acoustic, &
- timestamp_seconds_start)
- endif
-
!----- writing the kernels
! kernels output
if(SIMULATION_TYPE == 2) then
@@ -8069,13 +7978,8 @@
call write_seismograms(sisux,sisuz,siscurl,station_name,network_name,NSTEP, &
nrecloc,which_proc_receiver,nrec,myrank,deltat,seismotype,st_xval,t0, &
NTSTEP_BETWEEN_OUTPUT_SEISMO,seismo_offset,seismo_current,p_sv,&
- st_zval,x_source(1),z_source(1),SU_FORMAT,stage_time_scheme,SIMULATION_TYPE,c_LDDRK) !xiezhinan
+ st_zval,x_source(1),z_source(1),SU_FORMAT)
-! call write_seismograms(sisux,sisuz,siscurl,station_name,network_name,NSTEP, &
-! nrecloc,which_proc_receiver,nrec,myrank,deltat,seismotype,st_xval,t0, &
-! NTSTEP_BETWEEN_OUTPUT_SEISMO,seismo_offset,seismo_current,p_sv,&
-! st_zval,x_source(1),z_source(1),SU_FORMAT) !old
-
seismo_offset = seismo_offset + seismo_current
seismo_current = 0
Modified: seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/write_seismograms.F90
===================================================================
--- seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/write_seismograms.F90 2012-02-16 05:45:45 UTC (rev 19641)
+++ seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/write_seismograms.F90 2012-02-16 15:18:08 UTC (rev 19642)
@@ -44,15 +44,10 @@
! write seismograms to text files
-! subroutine write_seismograms(sisux,sisuz,siscurl,station_name,network_name, &
-! NSTEP,nrecloc,which_proc_receiver,nrec,myrank,deltat,seismotype,st_xval,t0, &
-! NTSTEP_BETWEEN_OUTPUT_SEISMO,seismo_offset,seismo_current,p_sv, &
-! st_zval,x_source,z_source,SU_FORMAT) !old
-
subroutine write_seismograms(sisux,sisuz,siscurl,station_name,network_name, &
NSTEP,nrecloc,which_proc_receiver,nrec,myrank,deltat,seismotype,st_xval,t0, &
NTSTEP_BETWEEN_OUTPUT_SEISMO,seismo_offset,seismo_current,p_sv, &
- st_zval,x_source,z_source,SU_FORMAT,stage_time_scheme,SIMULATION_TYPE,c_LDDRK) !xiezhinan
+ st_zval,x_source,z_source,SU_FORMAT)
implicit none
@@ -63,9 +58,6 @@
integer :: nrec,NSTEP,seismotype
integer :: NTSTEP_BETWEEN_OUTPUT_SEISMO,seismo_offset,seismo_current
- integer :: stage_time_scheme,SIMULATION_TYPE !xiezhinan
- integer :: i_stage !xiezhinan
- double precision, dimension(Nstages), intent(in) :: c_LDDRK !xiezhinan
double precision :: t0,deltat
! output seismograms in Seismic Unix format (adjoint traces will be read in the same format)
@@ -76,10 +68,8 @@
integer, intent(in) :: nrecloc,myrank
integer, dimension(nrec),intent(in) :: which_proc_receiver
-! double precision, dimension(NTSTEP_BETWEEN_OUTPUT_SEISMO,nrecloc), intent(in) :: sisux,sisuz,siscurl !old
+ double precision, dimension(NTSTEP_BETWEEN_OUTPUT_SEISMO,nrecloc), intent(in) :: sisux,sisuz,siscurl
- double precision, dimension(NTSTEP_BETWEEN_OUTPUT_SEISMO,stage_time_scheme,nrecloc), intent(in) :: sisux,sisuz,siscurl !xiezhinan
-
double precision :: st_xval(nrec)
character(len=MAX_LENGTH_STATION_NAME), dimension(nrec) :: station_name
@@ -92,8 +82,7 @@
character(len=150) sisname
! to write seismograms in single precision SEP and double precision binary format
-! double precision, dimension(:,:), allocatable :: buffer_binary !old
- double precision, dimension(:,:,:), allocatable :: buffer_binary
+ double precision, dimension(:,:), allocatable :: buffer_binary
! scaling factor for Seismic Unix xsu dislay
double precision, parameter :: FACTORXSU = 1.d0
@@ -140,8 +129,7 @@
number_of_components = NDIM
endif
-! allocate(buffer_binary(NTSTEP_BETWEEN_OUTPUT_SEISMO,number_of_components)) !old
- allocate(buffer_binary(NTSTEP_BETWEEN_OUTPUT_SEISMO,stage_time_scheme,number_of_components)) !xiezhinan
+ allocate(buffer_binary(NTSTEP_BETWEEN_OUTPUT_SEISMO,number_of_components))
if ( myrank == 0 .and. seismo_offset == 0 ) then
@@ -216,17 +204,12 @@
if ( which_proc_receiver(irec) == myrank ) then
irecloc = irecloc + 1
-! buffer_binary(:,1) = sisux(:,irecloc) !old
- buffer_binary(:,:,1) = sisux(:,:,irecloc) !xiezhinan
+ buffer_binary(:,1) = sisux(:,irecloc)
if ( number_of_components == 2 ) then
-! buffer_binary(:,2) = sisuz(:,irecloc) !old
- buffer_binary(:,:,2) = sisuz(:,:,irecloc) !xiezhinan
+ buffer_binary(:,2) = sisuz(:,irecloc)
else if ( number_of_components == 3 ) then
-! buffer_binary(:,2) = sisuz(:,irecloc) !old
-! buffer_binary(:,3) = siscurl(:,irecloc) !old
-
- buffer_binary(:,:,2) = sisuz(:,:,irecloc) !xiezhinan
- buffer_binary(:,:,3) = siscurl(:,:,irecloc) !xiezhinan
+ buffer_binary(:,2) = sisuz(:,irecloc)
+ buffer_binary(:,3) = siscurl(:,irecloc)
end if
#ifdef USE_MPI
@@ -243,8 +226,6 @@
call MPI_RECV(buffer_binary(1,3),NTSTEP_BETWEEN_OUTPUT_SEISMO,MPI_DOUBLE_PRECISION,&
which_proc_receiver(irec),irec,MPI_COMM_WORLD,status,ierror)
end if
-
-
#endif
end if
@@ -297,179 +278,71 @@
! make sure we never write more than the maximum number of time steps
! subtract offset of the source to make sure travel time is correct
do isample = 1,seismo_current
- do i_stage = 1, stage_time_scheme !xiezhinan
if(iorientation == 1) then
-! write(11,*) sngl(dble(seismo_offset+isample-1)*deltat - t0),' ', & !old
-! sngl(buffer_binary(isample,iorientation)) !old
-
- if(stage_time_scheme == 1 )then !xiezhinan
- write(11,*) sngl(dble(seismo_offset+isample-1)*deltat - t0),' ', & !xiezhinan
- sngl(buffer_binary(isample,i_stage,iorientation)) !xiezhinan
- endif
-
- if(stage_time_scheme == 6 .and. SIMULATION_TYPE == 1)then !xiezhinan
- write(11,*) sngl(dble(seismo_offset+isample-1)*deltat + c_LDDRK(i_stage)*deltat - t0),' ', & !xiezhinan
- sngl(buffer_binary(isample,i_stage,iorientation)) !xiezhinan
- endif
-
- if(stage_time_scheme == 6 .and. SIMULATION_TYPE == 2)then !xiezhinan
- write(11,*) sngl(dble(seismo_offset+isample-1)*deltat + c_LDDRK(i_stage)*deltat - t0),' ', & !xiezhinan
- sngl(buffer_binary(isample,i_stage,iorientation)) !xiezhinan
- endif
-
+ write(11,*) sngl(dble(seismo_offset+isample-1)*deltat - t0),' ', &
+ sngl(buffer_binary(isample,iorientation))
else
-! write(11,*) sngl(dble(seismo_offset+isample-1)*deltat - t0),' ', & !old
-! sngl(buffer_binary(isample,iorientation)) !old
-
- if(stage_time_scheme == 1 )then !xiezhinan
- write(11,*) sngl(dble(seismo_offset+isample-1)*deltat - t0),' ', & !xiezhinan
- sngl(buffer_binary(isample,i_stage,iorientation)) !xiezhinan
- endif
-
- if(stage_time_scheme == 6 .and. SIMULATION_TYPE == 1)then !xiezhinan
- write(11,*) sngl(dble(seismo_offset+isample-1)*deltat + c_LDDRK(i_stage)*deltat - t0),' ', & !xiezhinan
- sngl(buffer_binary(isample,i_stage,iorientation)) !xiezhinan
- endif
-
- if(stage_time_scheme == 6 .and. SIMULATION_TYPE == 2)then !xiezhinan
- write(11,*) sngl(dble(seismo_offset+isample-1)*deltat + c_LDDRK(i_stage)*deltat - t0),' ', & !xiezhinan
- sngl(buffer_binary(isample,i_stage,iorientation)) !xiezhinan
- endif
-
+ write(11,*) sngl(dble(seismo_offset+isample-1)*deltat - t0),' ', &
+ sngl(buffer_binary(isample,iorientation))
endif
- enddo
enddo
close(11)
end do
! write binary seismogram
do isample = 1, seismo_current
- do i_stage = 1, stage_time_scheme !xiezhinan
-! write(12,rec=(irec-1)*NSTEP+seismo_offset+isample) sngl(buffer_binary(isample,1)) !old
-! write(13,rec=(irec-1)*NSTEP+seismo_offset+isample) buffer_binary(isample,1) !old
-
- write(12,rec=(irec-1)*NSTEP*stage_time_scheme+seismo_offset+isample) sngl(buffer_binary(isample,i_stage,1)) !xiezhinan
- write(13,rec=(irec-1)*NSTEP*stage_time_scheme+seismo_offset+isample) buffer_binary(isample,i_stage,1) !xiezhinan
-
+ write(12,rec=(irec-1)*NSTEP+seismo_offset+isample) sngl(buffer_binary(isample,1))
+ write(13,rec=(irec-1)*NSTEP+seismo_offset+isample) buffer_binary(isample,1)
if ( seismotype /= 4 .and. seismotype /= 6 .and. p_sv) then
-! write(14,rec=(irec-1)*NSTEP+seismo_offset+isample) sngl(buffer_binary(isample,2)) !old
-! write(15,rec=(irec-1)*NSTEP+seismo_offset+isample) buffer_binary(isample,2) !old
-
- write(14,rec=(irec-1)*NSTEP*stage_time_scheme+seismo_offset+isample) &
- sngl(buffer_binary(isample,i_stage,2)) !xiezhinan
- write(15,rec=(irec-1)*NSTEP*stage_time_scheme+seismo_offset+isample) &
- buffer_binary(isample,i_stage,2) !xiezhinan
-
+ write(14,rec=(irec-1)*NSTEP+seismo_offset+isample) sngl(buffer_binary(isample,2))
+ write(15,rec=(irec-1)*NSTEP+seismo_offset+isample) buffer_binary(isample,2)
end if
-
if ( seismotype == 5 ) then
-! write(16,rec=(irec-1)*NSTEP+seismo_offset+isample) sngl(buffer_binary(isample,3)) !old
-! write(17,rec=(irec-1)*NSTEP+seismo_offset+isample) buffer_binary(isample,3) !old
-
- write(16,rec=(irec-1)*NSTEP+seismo_offset+isample) sngl(buffer_binary(isample,i_stage,3)) !xiezhinan
- write(17,rec=(irec-1)*NSTEP+seismo_offset+isample) buffer_binary(isample,i_stage,3) !xiezhinan
-
+ write(16,rec=(irec-1)*NSTEP+seismo_offset+isample) sngl(buffer_binary(isample,3))
+ write(17,rec=(irec-1)*NSTEP+seismo_offset+isample) buffer_binary(isample,3)
end if
-
- enddo !xiezhinan
-
enddo
else
if (seismo_offset==0) then
! write SU headers (refer to Seismic Unix for details)
-! write(12,rec=(irec-1)*60+(irec-1)*NSTEP+1) irec ! receiver ID !old
-! write(12,rec=(irec-1)*60+(irec-1)*NSTEP+10) NINT(st_xval(irec)-x_source) ! offset !old
-! write(12,rec=(irec-1)*60+(irec-1)*NSTEP+19) NINT(x_source) ! source location xs !old
-! write(12,rec=(irec-1)*60+(irec-1)*NSTEP+20) NINT(z_source) ! source location zs !old
-! write(12,rec=(irec-1)*60+(irec-1)*NSTEP+21) NINT(st_xval(irec)) ! receiver location xr !old
-! write(12,rec=(irec-1)*60+(irec-1)*NSTEP+22) NINT(st_zval(irec)) ! receiver location zr !old
-! if (nrec>1) write(12,rec=(irec-1)*60+(irec-1)*NSTEP+48) SNGL(st_xval(2)-st_xval(1)) ! receiver interval!old
-
- ! write SU headers (refer to Seismic Unix for details)
- write(12,rec=(irec-1)*60+(irec-1)*NSTEP*stage_time_scheme+1) &
- irec ! receiver ID !xiezhinan
- write(12,rec=(irec-1)*60+(irec-1)*NSTEP*stage_time_scheme+10) &
- NINT(st_xval(irec)-x_source) ! offset !xiezhinan
- write(12,rec=(irec-1)*60+(irec-1)*NSTEP*stage_time_scheme+19) &
- NINT(x_source) ! source location xs !xiezhinan
- write(12,rec=(irec-1)*60+(irec-1)*NSTEP*stage_time_scheme+20) &
- NINT(z_source) ! source location zs !xiezhinan
- write(12,rec=(irec-1)*60+(irec-1)*NSTEP*stage_time_scheme+21) &
- NINT(st_xval(irec)) ! receiver location xr !xiezhinan
- write(12,rec=(irec-1)*60+(irec-1)*NSTEP*stage_time_scheme+22) &
- NINT(st_zval(irec)) ! receiver location zr !xiezhinan
- if (nrec>1) write(12,rec=(irec-1)*60+(irec-1)*NSTEP*stage_time_scheme+48) &
- SNGL(st_xval(2)-st_xval(1)) ! receiver interval !xiezhinan
-
+ write(12,rec=(irec-1)*60+(irec-1)*NSTEP+1) irec ! receiver ID
+ write(12,rec=(irec-1)*60+(irec-1)*NSTEP+10) NINT(st_xval(irec)-x_source) ! offset
+ write(12,rec=(irec-1)*60+(irec-1)*NSTEP+19) NINT(x_source) ! source location xs
+ write(12,rec=(irec-1)*60+(irec-1)*NSTEP+20) NINT(z_source) ! source location zs
+ write(12,rec=(irec-1)*60+(irec-1)*NSTEP+21) NINT(st_xval(irec)) ! receiver location xr
+ write(12,rec=(irec-1)*60+(irec-1)*NSTEP+22) NINT(st_zval(irec)) ! receiver location zr
+ if (nrec>1) write(12,rec=(irec-1)*60+(irec-1)*NSTEP+48) SNGL(st_xval(2)-st_xval(1)) ! receiver interval
header2(1)=0 ! dummy
-! header2(2)=NSTEP !old
-! write(12,rec=(irec-1)*60+(irec-1)*NSTEP+29) header2 !old
-! header2(1)=NINT(deltat*1.0d6) ! deltat (unit: 10^{-6} second) !old
-! header2(2)=0 ! dummy !old
-
- header2(2)=NSTEP*stage_time_scheme !xiezhinan
- write(12,rec=(irec-1)*60+(irec-1)*NSTEP*stage_time_scheme+29) header2 !xiezhinan
- header2(1)=NINT(deltat*1.0d6) ! deltat (unit: 10^{-6} second) !xiezhinan
- header2(2)=0 ! dummy !xiezhinan
-
-! write(12,rec=(irec-1)*60+(irec-1)*NSTEP+30) header2 !old
-
- write(12,rec=(irec-1)*60+(irec-1)*NSTEP*stage_time_scheme+30) header2 !xiezhinan
-
+ header2(2)=NSTEP
+ write(12,rec=(irec-1)*60+(irec-1)*NSTEP+29) header2
+ header2(1)=NINT(deltat*1.0d6) ! deltat (unit: 10^{-6} second)
+ header2(2)=0 ! dummy
+ write(12,rec=(irec-1)*60+(irec-1)*NSTEP+30) header2
if ( seismotype /= 4 .and. seismotype /= 6 .and. p_sv) then
! headers
if (seismo_offset==0) then
-! write(14,rec=(irec-1)*60+(irec-1)*NSTEP+1) irec !old
-! write(14,rec=(irec-1)*60+(irec-1)*NSTEP+10) NINT(st_xval(irec)-x_source) !old
-! write(14,rec=(irec-1)*60+(irec-1)*NSTEP+19) NINT(x_source) !old
-! write(14,rec=(irec-1)*60+(irec-1)*NSTEP+20) NINT(z_source) !old
-! write(14,rec=(irec-1)*60+(irec-1)*NSTEP+21) NINT(st_xval(irec)) !old
-! write(14,rec=(irec-1)*60+(irec-1)*NSTEP+22) NINT(st_zval(irec)) !old
-! if(nrec>1) write(14,rec=(irec-1)*60+(irec-1)*NSTEP+48) SNGL(st_xval(2)-st_xval(1)) !old
-! header2(1)=0 ! dummy !old
-! header2(2)=NSTEP !old
-! write(14,rec=(irec-1)*60+(irec-1)*NSTEP+29) header2 !old
-! header2(1)=NINT(deltat*1.0d6) !old
-! header2(2)=0 ! dummy !old
-! write(14,rec=(irec-1)*60+(irec-1)*NSTEP+30) header2 !old
-
- write(14,rec=(irec-1)*60+(irec-1)*NSTEP*stage_time_scheme+1) irec
- write(14,rec=(irec-1)*60+(irec-1)*NSTEP*stage_time_scheme+10) NINT(st_xval(irec)-x_source)
- write(14,rec=(irec-1)*60+(irec-1)*NSTEP*stage_time_scheme+19) NINT(x_source)
- write(14,rec=(irec-1)*60+(irec-1)*NSTEP*stage_time_scheme+20) NINT(z_source)
- write(14,rec=(irec-1)*60+(irec-1)*NSTEP*stage_time_scheme+21) NINT(st_xval(irec))
- write(14,rec=(irec-1)*60+(irec-1)*NSTEP*stage_time_scheme+22) NINT(st_zval(irec))
- if(nrec>1) write(14,rec=(irec-1)*60+(irec-1)*NSTEP*stage_time_scheme+48) SNGL(st_xval(2)-st_xval(1))
-! header2(1)=0 ! dummy !old
-! header2(2)=NSTEP !old
-! write(14,rec=(irec-1)*60+(irec-1)*NSTEP+29) header2 !old
-! header2(1)=NINT(deltat*1.0d6) !old
-! header2(2)=0 ! dummy !old
-! write(14,rec=(irec-1)*60+(irec-1)*NSTEP+30) header2 !old
-
- header2(1)=0 ! dummy !xiezhinan
- header2(2)=NSTEP*stage_time_scheme !xiezhinan
- write(14,rec=(irec-1)*60+(irec-1)*NSTEP*stage_time_scheme+29) header2 !xiezhinan
- header2(1)=NINT(deltat*1.0d6) !xiezhinan
- header2(2)=0 ! dummy !xiezhinan
- write(14,rec=(irec-1)*60+(irec-1)*NSTEP*stage_time_scheme+30) header2 !xiezhinan
-
+ write(14,rec=(irec-1)*60+(irec-1)*NSTEP+1) irec
+ write(14,rec=(irec-1)*60+(irec-1)*NSTEP+10) NINT(st_xval(irec)-x_source)
+ write(14,rec=(irec-1)*60+(irec-1)*NSTEP+19) NINT(x_source)
+ write(14,rec=(irec-1)*60+(irec-1)*NSTEP+20) NINT(z_source)
+ write(14,rec=(irec-1)*60+(irec-1)*NSTEP+21) NINT(st_xval(irec))
+ write(14,rec=(irec-1)*60+(irec-1)*NSTEP+22) NINT(st_zval(irec))
+ if(nrec>1) write(14,rec=(irec-1)*60+(irec-1)*NSTEP+48) SNGL(st_xval(2)-st_xval(1))
+ header2(1)=0 ! dummy
+ header2(2)=NSTEP
+ write(14,rec=(irec-1)*60+(irec-1)*NSTEP+29) header2
+ header2(1)=NINT(deltat*1.0d6)
+ header2(2)=0 ! dummy
+ write(14,rec=(irec-1)*60+(irec-1)*NSTEP+30) header2
end if
endif
endif
! the "60" in the following corresponds to 240 bytes header (note the reclength is 4 bytes)
do isample = 1, seismo_current
- do i_stage = 1, stage_time_scheme !xiezhinan
-! write(12,rec=irec*60+(irec-1)*NSTEP+seismo_offset+isample) sngl(buffer_binary(isample,1)) !old
-! if ( seismotype /= 4 .and. seismotype /= 6 .and. p_sv) then !old
-! write(14,rec=irec*60+(irec-1)*NSTEP+seismo_offset+isample) sngl(buffer_binary(isample,2)) !old
-! end if
-
- write(12,rec=irec*60+(irec-1)*NSTEP*stage_time_scheme+seismo_offset+isample) sngl(buffer_binary(isample,i_stage,1)) !xiezhinan
- if ( seismotype /= 4 .and. seismotype /= 6 .and. p_sv) then !xiezhinan
- write(14,rec=irec*60+(irec-1)*NSTEP*stage_time_scheme+seismo_offset+isample) sngl(buffer_binary(isample,i_stage,2)) !xiezhinan
- end if !xiezhinan
- enddo !xiezhinan
+ write(12,rec=irec*60+(irec-1)*NSTEP+seismo_offset+isample) sngl(buffer_binary(isample,1))
+ if ( seismotype /= 4 .and. seismotype /= 6 .and. p_sv) then
+ write(14,rec=irec*60+(irec-1)*NSTEP+seismo_offset+isample) sngl(buffer_binary(isample,2))
+ end if
enddo
endif
@@ -477,15 +350,12 @@
else
if ( which_proc_receiver(irec) == myrank ) then
irecloc = irecloc + 1
-! call MPI_SEND(sisux(1,irecloc),NTSTEP_BETWEEN_OUTPUT_SEISMO,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ierror) !old
- call MPI_SEND(sisux(1,1,irecloc),NTSTEP_BETWEEN_OUTPUT_SEISMO,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ierror) !xiezhinan
+ call MPI_SEND(sisux(1,irecloc),NTSTEP_BETWEEN_OUTPUT_SEISMO,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ierror)
if ( number_of_components >= 2 ) then
-! call MPI_SEND(sisuz(1,irecloc),NTSTEP_BETWEEN_OUTPUT_SEISMO,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ierror) !old
- call MPI_SEND(sisuz(1,1,irecloc),NTSTEP_BETWEEN_OUTPUT_SEISMO,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ierror)!xiezhinan
+ call MPI_SEND(sisuz(1,irecloc),NTSTEP_BETWEEN_OUTPUT_SEISMO,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ierror)
end if
if ( number_of_components == 3 ) then
-! call MPI_SEND(siscurl(1,irecloc),NTSTEP_BETWEEN_OUTPUT_SEISMO,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ierror) !old
- call MPI_SEND(siscurl(1,1,irecloc),NTSTEP_BETWEEN_OUTPUT_SEISMO,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ierror)!xiezhinan
+ call MPI_SEND(siscurl(1,irecloc),NTSTEP_BETWEEN_OUTPUT_SEISMO,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ierror)
end if
end if
#endif
More information about the CIG-COMMITS
mailing list