[cig-commits] r19592 - seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Tue Feb 7 08:23:34 PST 2012
Author: xie.zhinan
Date: 2012-02-07 08:23:33 -0800 (Tue, 07 Feb 2012)
New Revision: 19592
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/prepare_absorb.f90
seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/specfem2D.F90
seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/write_seismograms.F90
Log:
add LDDRK time scheme for pure elastic adjoint 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-07 01:19:09 UTC (rev 19591)
+++ seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/compute_arrays_source.f90 2012-02-07 16:23:33 UTC (rev 19592)
@@ -129,23 +129,26 @@
! ------------------------------------------------------------------------------------------------------
+! 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
+! 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
! Gauss-Lobatto-Legendre points of integration and weights
double precision, dimension(NGLLX) :: xigll
@@ -153,7 +156,8 @@
double precision :: hxir(NGLLX), hpxir(NGLLX), hgammar(NGLLZ), hpgammar(NGLLZ)
- real(kind=CUSTOM_REAL) :: adj_src_s(NSTEP,3)
+! real(kind=CUSTOM_REAL) :: adj_src_s(NSTEP,3) !old
+ real(kind=CUSTOM_REAL) :: adj_src_s(NSTEP,stage_time_scheme,3) !xiezhinan
integer icomp, itime, i, k, ios
double precision :: junk
@@ -163,7 +167,8 @@
call lagrange_any(xi_receiver,NGLLX,xigll,hxir,hpxir)
call lagrange_any(gamma_receiver,NGLLZ,zigll,hgammar,hpgammar)
- adj_sourcearray(:,:,:,:) = 0.
+! adj_sourcearray(:,:,:,:) = 0. !old
+ adj_sourcearray(:,:,:,:,:) = 0. !xiezhinan
comp = (/"BXX","BXY","BXZ"/)
@@ -174,7 +179,10 @@
if (ios /= 0) call exit_MPI(' file '//trim(filename)//'does not exist')
do itime = 1, NSTEP
- read(IIN,*) junk, adj_src_s(itime,icomp)
+ 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
enddo
close(IIN)
@@ -182,7 +190,8 @@
do k = 1, NGLLZ
do i = 1, NGLLX
- adj_sourcearray(:,:,i,k) = hxir(i) * hgammar(k) * adj_src_s(:,:)
+! adj_sourcearray(:,:,i,k) = hxir(i) * hgammar(k) * adj_src_s(:,:) !old
+ adj_sourcearray(:,:,:,i,k) = hxir(i) * hgammar(k) * adj_src_s(:,:,:) !xiezhinan
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-07 01:19:09 UTC (rev 19591)
+++ seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/compute_forces_viscoelastic.f90 2012-02-07 16:23:33 UTC (rev 19592)
@@ -112,13 +112,19 @@
real(kind=CUSTOM_REAL), dimension(NSOURCES,NDIM,NGLLX,NGLLZ) :: sourcearray
real(kind=CUSTOM_REAL), dimension(3,nglob) :: b_accel_elastic,b_displ_elastic
- real(kind=CUSTOM_REAL), dimension(nrec,NSTEP,3,NGLLX,NGLLZ) :: adj_sourcearrays
+! 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(nglob) :: mu_k,kappa_k
- 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) :: 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,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) :: e1_LDDRK,e11_LDDRK,e13_LDDRK
@@ -205,7 +211,7 @@
! compute Grad(displ_elastic) at time step n for attenuation
if(ATTENUATION_VISCOELASTIC_SOLID) then
- temp_displ_elastic = displ_elastic + deltat * veloc_elastic + 0.5 * deltat**2 * accel_elastic !xiezhinan
+ 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)
endif
@@ -579,20 +585,33 @@
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
- b_absorb_elastic_left(3,j,ib_left(ispecabs),it) = (tz + traction_z_t0)*weight
+! 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
else !SH (membrane) waves
- b_absorb_elastic_left(2,j,ib_left(ispecabs),it) = ty*weight
+! 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
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_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)
+! 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
else !SH (membrane) waves
- b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - &
- b_absorb_elastic_left(2,j,ib_left(ispecabs),NSTEP-it+1)
+! 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
+
endif
endif
@@ -670,20 +689,35 @@
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
- b_absorb_elastic_right(3,j,ib_right(ispecabs),it) = (tz - traction_z_t0)*weight
+! 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
+
else! SH (membrane) waves
- b_absorb_elastic_right(2,j,ib_right(ispecabs),it) = ty*weight
+! 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
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_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)
+! 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
+
else! SH (membrane) waves
- b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - &
- b_absorb_elastic_right(2,j,ib_right(ispecabs),NSTEP-it+1)
+! 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
+
endif
endif
@@ -776,20 +810,34 @@
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
- b_absorb_elastic_bottom(3,i,ib_bottom(ispecabs),it) = (tz + traction_z_t0)*weight
+! 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
+
else!SH (membrane) waves
- b_absorb_elastic_bottom(2,i,ib_bottom(ispecabs),it) = ty*weight
+! 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
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_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)
+! 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
+
else!SH (membrane) waves
- b_accel_elastic(2,iglob) = b_accel_elastic(2,iglob) - &
- b_absorb_elastic_bottom(2,i,ib_bottom(ispecabs),NSTEP-it+1)
+! 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
endif
endif
@@ -874,17 +922,30 @@
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
- b_absorb_elastic_top(3,i,ib_top(ispecabs),it) = (tz- traction_z_t0)*weight
+! 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
+
else!SH (membrane) waves
- b_absorb_elastic_top(2,i,ib_top(ispecabs),it) = ty*weight
+! 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
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)
- b_accel_elastic(3,iglob) = b_accel_elastic(3,iglob) - b_absorb_elastic_top(3,i,ib_top(ispecabs),NSTEP-it+1)
+! 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
+
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)
+! 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
endif
endif
@@ -929,6 +990,7 @@
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
@@ -951,11 +1013,18 @@
do j=1,NGLLZ
do i=1,NGLLX
iglob = ibool(i,j,ispec_selected_rec(irec))
- 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)
+ 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
else !SH (membrane) waves
- accel_elastic(2,iglob) = accel_elastic(2,iglob) + adj_sourcearrays(irec_local,NSTEP-it+1,2,i,j)
+! 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
endif
enddo
enddo
@@ -990,18 +1059,18 @@
! evolution e1
if(stage_time_scheme == 1) then
- Un = e1(i,j,ispec,i_sls)
- tauinv = - inv_tau_sigma_nu1(i,j,ispec,i_sls)
- tauinvsquare = tauinv * tauinv
- tauinvcube = tauinvsquare * tauinv
- tauinvUn = tauinv * Un
- Sn = theta_n * phi_nu1(i,j,ispec,i_sls)
- Snp1 = theta_np1 * phi_nu1(i,j,ispec,i_sls)
- Unp1 = Un + (deltatfourth*tauinvcube*(Sn + tauinvUn) + &
- twelvedeltat*(Sn + Snp1 + 2*tauinvUn) + &
- fourdeltatsquare*tauinv*(2*Sn + Snp1 + 3*tauinvUn) + &
- deltatcube*tauinvsquare*(3*Sn + Snp1 + 4*tauinvUn))* ONE_OVER_24
- e1(i,j,ispec,i_sls) = Unp1
+ Un = e1(i,j,ispec,i_sls)
+ tauinv = - inv_tau_sigma_nu1(i,j,ispec,i_sls)
+ tauinvsquare = tauinv * tauinv
+ tauinvcube = tauinvsquare * tauinv
+ tauinvUn = tauinv * Un
+ Sn = theta_n * phi_nu1(i,j,ispec,i_sls)
+ Snp1 = theta_np1 * phi_nu1(i,j,ispec,i_sls)
+ Unp1 = Un + (deltatfourth*tauinvcube*(Sn + tauinvUn) + &
+ twelvedeltat*(Sn + Snp1 + 2*tauinvUn) + &
+ fourdeltatsquare*tauinv*(2*Sn + Snp1 + 3*tauinvUn) + &
+ deltatcube*tauinvsquare*(3*Sn + Snp1 + 4*tauinvUn))* ONE_OVER_24
+ e1(i,j,ispec,i_sls) = Unp1
endif
if(stage_time_scheme == 6) then
@@ -1009,158 +1078,129 @@
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_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * e1_LDDRK(i,j,ispec,i_sls) + &
+ e1_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * e1_LDDRK(i,j,ispec,i_sls) + &
deltat * (theta_n * temp_time_scheme - Un * tauinv)
- e1(i,j,ispec,i_sls) = Un + beta_LDDRK(i_stage) * e1_LDDRK(i,j,ispec,i_sls)
+ e1(i,j,ispec,i_sls) = Un + beta_LDDRK(i_stage) * 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)
+ 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 .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
+ if(i_stage==1)then
+ e1_initial_rk(i,j,ispec,i_sls) = e1(i,j,ispec,i_sls)
+ endif
- e1_initial_rk(i,j,ispec,i_sls) = e1(i,j,ispec,i_sls)
+ 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)
- endif
+ elseif(i_stage==4)then
- 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)
+ 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
- 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))
-
+ ! evolution e11
+ if(stage_time_scheme == 1) then
+ Un = e11(i,j,ispec,i_sls)
+ tauinv = - inv_tau_sigma_nu2(i,j,ispec,i_sls)
+ tauinvsquare = tauinv * tauinv
+ tauinvcube = tauinvsquare * tauinv
+ tauinvUn = tauinv * Un
+ Sn = (dux_dxl_n(i,j,ispec) - theta_n/TWO) * phi_nu2(i,j,ispec,i_sls)
+ Snp1 = (dux_dxl_np1(i,j,ispec) - theta_np1/TWO) * phi_nu2(i,j,ispec,i_sls)
+ Unp1 = Un + (deltatfourth*tauinvcube*(Sn + tauinvUn) + &
+ twelvedeltat*(Sn + Snp1 + 2*tauinvUn) + &
+ fourdeltatsquare*tauinv*(2*Sn + Snp1 + 3*tauinvUn) + &
+ deltatcube*tauinvsquare*(3*Sn + Snp1 + 4*tauinvUn))* ONE_OVER_24
+ e11(i,j,ispec,i_sls) = Unp1
endif
+ if(stage_time_scheme == 6) then
+ temp_time_scheme = dux_dxl_n(i,j,ispec)-theta_n/TWO
+ temper_time_scheme = phi_nu2(i,j,ispec,i_sls)
+ e11_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * e11_LDDRK(i,j,ispec,i_sls) &
+ + deltat * (temp_time_scheme * temper_time_scheme)- &
+ deltat * (e11(i,j,ispec,i_sls) * inv_tau_sigma_nu2(i,j,ispec,i_sls))
+ e11(i,j,ispec,i_sls) = e11(i,j,ispec,i_sls)+beta_LDDRK(i_stage)*e11_LDDRK(i,j,ispec,i_sls)
endif
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- ! evolution e11
- if(stage_time_scheme == 1) then
- Un = e11(i,j,ispec,i_sls)
- tauinv = - inv_tau_sigma_nu2(i,j,ispec,i_sls)
- tauinvsquare = tauinv * tauinv
- tauinvcube = tauinvsquare * tauinv
- tauinvUn = tauinv * Un
- Sn = (dux_dxl_n(i,j,ispec) - theta_n/TWO) * phi_nu2(i,j,ispec,i_sls)
- Snp1 = (dux_dxl_np1(i,j,ispec) - theta_np1/TWO) * phi_nu2(i,j,ispec,i_sls)
- Unp1 = Un + (deltatfourth*tauinvcube*(Sn + tauinvUn) + &
- twelvedeltat*(Sn + Snp1 + 2*tauinvUn) + &
- fourdeltatsquare*tauinv*(2*Sn + Snp1 + 3*tauinvUn) + &
- deltatcube*tauinvsquare*(3*Sn + Snp1 + 4*tauinvUn))* ONE_OVER_24
- e11(i,j,ispec,i_sls) = Unp1
- endif
-
- if(stage_time_scheme == 6) then
- temp_time_scheme = dux_dxl_n(i,j,ispec)-theta_n/TWO
- temper_time_scheme = phi_nu2(i,j,ispec,i_sls)
- e11_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * e11_LDDRK(i,j,ispec,i_sls) &
- + deltat * (temp_time_scheme * temper_time_scheme)- &
- deltat * (e11(i,j,ispec,i_sls) * inv_tau_sigma_nu2(i,j,ispec,i_sls))
- e11(i,j,ispec,i_sls) = e11(i,j,ispec,i_sls)+beta_LDDRK(i_stage)*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))
+ 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
-
+ 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 = e13(i,j,ispec,i_sls)
- tauinv = - inv_tau_sigma_nu2(i,j,ispec,i_sls)
- tauinvsquare = tauinv * tauinv
- tauinvcube = tauinvsquare * tauinv
- tauinvUn = tauinv * Un
- Sn = (dux_dzl_n(i,j,ispec) + duz_dxl_n(i,j,ispec)) * phi_nu2(i,j,ispec,i_sls)
- Snp1 = (dux_dzl_np1(i,j,ispec) + duz_dxl_np1(i,j,ispec)) * phi_nu2(i,j,ispec,i_sls)
- Unp1 = Un + (deltatfourth*tauinvcube*(Sn + tauinvUn) + &
- twelvedeltat*(Sn + Snp1 + 2*tauinvUn) + &
- fourdeltatsquare*tauinv*(2*Sn + Snp1 + 3*tauinvUn) + &
- deltatcube*tauinvsquare*(3*Sn + Snp1 + 4*tauinvUn))* ONE_OVER_24
- e13(i,j,ispec,i_sls) = Unp1
- endif
+ if(stage_time_scheme == 1) then
+ Un = e13(i,j,ispec,i_sls)
+ tauinv = - inv_tau_sigma_nu2(i,j,ispec,i_sls)
+ tauinvsquare = tauinv * tauinv
+ tauinvcube = tauinvsquare * tauinv
+ tauinvUn = tauinv * Un
+ Sn = (dux_dzl_n(i,j,ispec) + duz_dxl_n(i,j,ispec)) * phi_nu2(i,j,ispec,i_sls)
+ Snp1 = (dux_dzl_np1(i,j,ispec) + duz_dxl_np1(i,j,ispec)) * phi_nu2(i,j,ispec,i_sls)
+ Unp1 = Un + (deltatfourth*tauinvcube*(Sn + tauinvUn) + &
+ twelvedeltat*(Sn + Snp1 + 2*tauinvUn) + &
+ fourdeltatsquare*tauinv*(2*Sn + Snp1 + 3*tauinvUn) + &
+ deltatcube*tauinvsquare*(3*Sn + Snp1 + 4*tauinvUn))* ONE_OVER_24
+ e13(i,j,ispec,i_sls) = Unp1
+ endif
if(stage_time_scheme == 6) 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_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * e13_LDDRK(i,j,ispec,i_sls)+&
+ 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_LDDRK(i,j,ispec,i_sls) = alpha_LDDRK(i_stage) * e13_LDDRK(i,j,ispec,i_sls)+&
deltat * (temp_time_scheme*temper_time_scheme)- &
deltat * (e13(i,j,ispec,i_sls) * inv_tau_sigma_nu2(i,j,ispec,i_sls))
- e13(i,j,ispec,i_sls) = e13(i,j,ispec,i_sls)+beta_LDDRK(i_stage) * e13_LDDRK(i,j,ispec,i_sls)
- endif
+ e13(i,j,ispec,i_sls) = e13(i,j,ispec,i_sls)+beta_LDDRK(i_stage) * 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- &
+ 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
-
+ 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
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-07 01:19:09 UTC (rev 19591)
+++ seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/prepare_absorb.f90 2012-02-07 16:23:33 UTC (rev 19592)
@@ -254,45 +254,67 @@
!-------------------------------------------------------------------------------------------------
!
- 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, & !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, & !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
- 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)
+ integer :: stage_time_scheme,i_stage !xiezhinan
+! 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)
+! read(35) b_absorb_elastic_left(1,i,ispec,it) !old
+ read(35) b_absorb_elastic_left(1,i,ispec,it,i_stage) !xiezhinan
enddo
do i=1,NGLLZ
- read(35) b_absorb_elastic_left(3,i,ispec,it)
+! read(35) b_absorb_elastic_left(3,i,ispec,it) !old
+ read(35) b_absorb_elastic_left(3,i,ispec,it,i_stage) !xiezhinan
enddo
- b_absorb_elastic_left(2,:,ispec,it) = ZERO
+! b_absorb_elastic_left(2,:,ispec,it) = ZERO !old
+ b_absorb_elastic_left(2,:,ispec,it,i_stage) = ZERO !xiezhinan
else!SH (membrane) waves
do i=1,NGLLZ
- read(35) b_absorb_elastic_left(2,i,ispec,it)
+! read(35) b_absorb_elastic_left(2,i,ispec,it) !old
+ read(35) b_absorb_elastic_left(2,i,ispec,it,i_stage) !xiezhinan
enddo
- b_absorb_elastic_left(1,:,ispec,it) = ZERO
- b_absorb_elastic_left(3,:,ispec,it) = ZERO
+! 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
+
endif
enddo
@@ -304,18 +326,26 @@
if(p_sv)then!P-SV waves
do i=1,NGLLZ
- read(36) b_absorb_elastic_right(1,i,ispec,it)
+! read(36) b_absorb_elastic_right(1,i,ispec,it) !old
+ read(36) b_absorb_elastic_right(1,i,ispec,it,i_stage) !xiezhinan
enddo
do i=1,NGLLZ
- read(36) b_absorb_elastic_right(3,i,ispec,it)
+! read(36) b_absorb_elastic_right(3,i,ispec,it) !old
+ read(36) b_absorb_elastic_right(3,i,ispec,it,i_stage) !xiezhinan
enddo
- b_absorb_elastic_right(2,:,ispec,it) = ZERO
+! b_absorb_elastic_right(2,:,ispec,it) = ZERO !old
+ b_absorb_elastic_right(2,:,ispec,it,i_stage) = ZERO !xiezhinan
else!SH (membrane) waves
do i=1,NGLLZ
- read(36) b_absorb_elastic_right(2,i,ispec,it)
+! read(36) b_absorb_elastic_right(2,i,ispec,it) !old
+ read(36) b_absorb_elastic_right(2,i,ispec,it,i_stage) !xiezhinan
enddo
- b_absorb_elastic_right(1,:,ispec,it) = ZERO
- b_absorb_elastic_right(3,:,ispec,it) = ZERO
+! 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
+
endif
enddo
@@ -327,18 +357,24 @@
if(p_sv)then!P-SV waves
do i=1,NGLLX
- read(37) b_absorb_elastic_bottom(1,i,ispec,it)
+! read(37) b_absorb_elastic_bottom(1,i,ispec,it) !old
+ read(37) b_absorb_elastic_bottom(1,i,ispec,it,i_stage) !xiezhinan
enddo
do i=1,NGLLX
- read(37) b_absorb_elastic_bottom(3,i,ispec,it)
+! read(37) b_absorb_elastic_bottom(3,i,ispec,it) !old
+ read(37) b_absorb_elastic_bottom(3,i,ispec,it,i_stage) !xiezhinan
enddo
- b_absorb_elastic_bottom(2,:,ispec,it) = ZERO
+! b_absorb_elastic_bottom(2,:,ispec,it) = ZERO !old
+ b_absorb_elastic_bottom(2,:,ispec,it,i_stage) = ZERO !xiezhinan
else!SH (membrane) waves
do i=1,NGLLZ
- read(37) b_absorb_elastic_bottom(2,i,ispec,it)
+! read(37) b_absorb_elastic_bottom(2,i,ispec,it) !old
+ read(37) b_absorb_elastic_bottom(2,i,ispec,it,i_stage) !xiezhinan
enddo
- b_absorb_elastic_bottom(1,:,ispec,it) = ZERO
- b_absorb_elastic_bottom(3,:,ispec,it) = ZERO
+! 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
endif
enddo
@@ -350,23 +386,31 @@
if(p_sv)then!P-SV waves
do i=1,NGLLX
- read(38) b_absorb_elastic_top(1,i,ispec,it)
+! read(38) b_absorb_elastic_top(1,i,ispec,it) !old
+ read(38) b_absorb_elastic_top(1,i,ispec,it,i_stage) !xiezhinan
enddo
do i=1,NGLLX
- read(38) b_absorb_elastic_top(3,i,ispec,it)
+! read(38) b_absorb_elastic_top(3,i,ispec,it) !old
+ read(38) b_absorb_elastic_top(3,i,ispec,it,i_stage) !xiezhinan
enddo
- b_absorb_elastic_top(2,:,ispec,it) = ZERO
+! b_absorb_elastic_top(2,:,ispec,it) = ZERO !old
+ b_absorb_elastic_top(2,:,ispec,it,i_stage) = ZERO !xiezhinan
else!SH (membrane) waves
do i=1,NGLLZ
- read(38) b_absorb_elastic_top(2,i,ispec,it)
+! read(38) b_absorb_elastic_top(2,i,ispec,it) !old
+ read(38) b_absorb_elastic_top(2,i,ispec,it,i_stage) !xiezhinan
enddo
- b_absorb_elastic_top(1,:,ispec,it) = ZERO
- b_absorb_elastic_top(3,:,ispec,it) = ZERO
+! 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
+
endif
enddo
endif
-
+ enddo !xiezhinan
enddo
end subroutine prepare_absorb_elastic
Modified: seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/specfem2D.F90 2012-02-07 01:19:09 UTC (rev 19591)
+++ seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/specfem2D.F90 2012-02-07 16:23:33 UTC (rev 19592)
@@ -387,7 +387,8 @@
character(len=150) dummystring
! for seismograms
- double precision, dimension(:,:), allocatable :: sisux,sisuz,siscurl
+ double precision, dimension(:,:,:), allocatable :: sisux,sisuz,siscurl !xiezhinan
+! double precision, dimension(:,:), allocatable :: sisux,sisuz,siscurl !old
integer :: seismo_offset, seismo_current
! vector field in an element
@@ -633,6 +634,8 @@
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
@@ -655,16 +658,29 @@
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
- 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 :: 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 :: &
b_absorb_acoustic_left,b_absorb_acoustic_right,&
b_absorb_acoustic_bottom, b_absorb_acoustic_top
@@ -863,7 +879,8 @@
!<SU_FORMAT
integer(kind=4) :: r4head(60)
character(len=512) :: filename
- real(kind=4),dimension(:,:),allocatable :: adj_src_s
+! real(kind=4),dimension(:,:),allocatable :: adj_src_s !old
+ real(kind=4),dimension(:,:,:),allocatable :: adj_src_s !xiezhinan
integer(kind=2) :: header2(2)
!>SU_FORMAT
@@ -969,6 +986,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(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
@@ -1423,15 +1442,27 @@
! 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))
- 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))
+! 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
+
else
- 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))
+! 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
+
endif
if(any_poroelastic .and. (SAVE_FORWARD .or. SIMULATION_TYPE == 2)) then
allocate(b_absorb_poro_s_left(NDIM,NGLLZ,nspec_left,NSTEP))
@@ -1468,10 +1499,15 @@
else
if(.not. allocated(b_absorb_elastic_left)) then
- 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))
+! 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
endif
if(.not. allocated(b_absorb_poro_s_left)) then
@@ -2084,11 +2120,14 @@
nadj_rec_local = nadj_rec_local + 1
endif
enddo
- if(ipass == 1) allocate(adj_sourcearray(NSTEP,3,NGLLX,NGLLZ))
+! 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 (nadj_rec_local > 0 .and. ipass == 1) then
- allocate(adj_sourcearrays(nadj_rec_local,NSTEP,3,NGLLX,NGLLZ))
+! 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
else if (ipass == 1) then
- allocate(adj_sourcearrays(1,1,1,1,1))
+! allocate(adj_sourcearrays(1,1,1,1,1)) !old
+ allocate(adj_sourcearrays(1,1,1,1,1,1)) !xiezhinan
endif
if (.not. SU_FORMAT) then
@@ -2098,30 +2137,42 @@
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, &
- xi_receiver(irec), gamma_receiver(irec), &
- adj_sourcearray, xigll,zigll,NSTEP)
- adj_sourcearrays(irec_local,:,:,:,:) = adj_sourcearray(:,:,:,:)
+! 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
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)
+! 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
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)
+! 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
if (ios /= 0) call exit_MPI(' file '//trim(filename)//'does not exist')
- allocate(adj_src_s(NSTEP,3))
+! allocate(adj_src_s(NSTEP,3)) !old
+ allocate(adj_src_s(NSTEP,stage_time_scheme,3)) !xiezhinan
do irec = 1, nrec
if(myrank == which_proc_receiver(irec))then
irec_local = irec_local + 1
- adj_sourcearray(:,:,:,:) = 0.0
- read(111,rec=irec,iostat=ios) r4head, adj_src_s(:,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
if (ios /= 0) call exit_MPI(' file '//trim(filename)//' read error')
- read(113,rec=irec,iostat=ios) r4head, adj_src_s(:,3)
+! read(113,rec=irec,iostat=ios) r4head, adj_src_s(:,3) !old
+ read(113,rec=irec,iostat=ios) r4head, adj_src_s(:,:,3)!xiezhinan
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)
@@ -2129,10 +2180,12 @@
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(:,:)
+! adj_sourcearray(:,:,i,k) = hxir(i) * hgammar(k) * adj_src_s(:,:) !old
+ adj_sourcearray(:,:,:,i,k) = hxir(i) * hgammar(k) * adj_src_s(:,:,:) !xiezhinan
enddo
enddo
- adj_sourcearrays(irec_local,:,:,:,:) = adj_sourcearray(:,:,:,:)
+! adj_sourcearrays(irec_local,:,:,:,:) = adj_sourcearray(:,:,:,:) !old
+ adj_sourcearrays(irec_local,:,:,:,:,:) = adj_sourcearray(:,:,:,:,:) !xiezhinan
endif ! if(myrank == which_proc_receiver(irec))
enddo ! irec
close(111)
@@ -2140,7 +2193,8 @@
deallocate(adj_src_s)
endif
else if (ipass == 1) then
- allocate(adj_sourcearrays(1,1,1,1,1))
+! allocate(adj_sourcearrays(1,1,1,1,1)) !old
+ allocate(adj_sourcearrays(1,1,1,1,1,1)) !xiezhinan
endif
if (ipass == 1) then
@@ -2364,9 +2418,13 @@
! allocate seismogram arrays
if(ipass == 1) then
- allocate(sisux(NTSTEP_BETWEEN_OUTPUT_SEISMO,nrecloc))
- allocate(sisuz(NTSTEP_BETWEEN_OUTPUT_SEISMO,nrecloc))
- allocate(siscurl(NTSTEP_BETWEEN_OUTPUT_SEISMO,nrecloc))
+ 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
endif
! check if acoustic receiver is exactly on the free surface because pressure is zero there
@@ -2458,6 +2516,12 @@
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))
@@ -2478,6 +2542,12 @@
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))
@@ -3159,11 +3229,16 @@
! reads in absorbing boundary data
if(any_elastic) then
- 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, & !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, & !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, &
@@ -3755,11 +3830,19 @@
! solid/porous edge detection
! the two elements forming an edge are already known (computed in meshfem2D),
! the common nodes forming the edge are computed here
- if(coupled_elastic_poro) then
+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) &
stop 'Attenuation not supported for mixed elastic/poroelastic simulations'
+ if(time_stepping_scheme == 2.or. time_stepping_scheme == 3) &
+ stop 'RK and LDDRK time scheme not supported for mixed elastic/poroelastic simulations'
+
if ( myrank == 0 ) then
print *
print *,'Mixed elastic/poroelastic simulation'
@@ -4329,10 +4412,13 @@
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
@@ -4356,12 +4442,19 @@
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
@@ -4926,8 +5019,10 @@
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) &
- + adj_sourcearrays(irec_local,NSTEP-it+1,1,i,j)
+! 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
enddo
enddo
endif ! if element acoustic
@@ -5100,14 +5195,17 @@
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)
+! write(35) b_absorb_elastic_left(1,i,ispec,it) !old
+ write(35) b_absorb_elastic_left(1,i,ispec,it,i_stage) !xiezhinan
enddo
do i=1,NGLLZ
- write(35) b_absorb_elastic_left(3,i,ispec,it)
+! write(35) b_absorb_elastic_left(3,i,ispec,it) !old
+ write(35) b_absorb_elastic_left(3,i,ispec,it,i_stage) !xiezhinan
enddo
else!SH (membrane) waves
do i=1,NGLLZ
- write(35) b_absorb_elastic_left(2,i,ispec,it)
+! write(35) b_absorb_elastic_left(2,i,ispec,it) !old
+ write(35) b_absorb_elastic_left(2,i,ispec,it,i_stage) !xiezhinan
enddo
endif
enddo
@@ -5118,14 +5216,17 @@
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)
+! write(36) b_absorb_elastic_right(1,i,ispec,it) !old
+ write(36) b_absorb_elastic_right(1,i,ispec,it,i_stage) !xiezhinan
enddo
do i=1,NGLLZ
- write(36) b_absorb_elastic_right(3,i,ispec,it)
+! write(36) b_absorb_elastic_right(3,i,ispec,it) !old
+ write(36) b_absorb_elastic_right(3,i,ispec,it,i_stage) !xiezhinan
enddo
else!SH (membrane) waves
do i=1,NGLLZ
- write(36) b_absorb_elastic_right(2,i,ispec,it)
+! write(36) b_absorb_elastic_right(2,i,ispec,it) !old
+ write(36) b_absorb_elastic_right(2,i,ispec,it,i_stage) !xiezhinan
enddo
endif
enddo
@@ -5136,14 +5237,18 @@
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)
+! write(37) b_absorb_elastic_bottom(1,i,ispec,it) !old
+ write(37) b_absorb_elastic_bottom(1,i,ispec,it,i_stage) !xiezhinan
enddo
do i=1,NGLLX
- write(37) b_absorb_elastic_bottom(3,i,ispec,it)
+! write(37) b_absorb_elastic_bottom(3,i,ispec,it) !old
+ write(37) b_absorb_elastic_bottom(3,i,ispec,it,i_stage) !xiezhinan
+
enddo
else!SH (membrane) waves
do i=1,NGLLX
- write(37) b_absorb_elastic_bottom(2,i,ispec,it)
+! write(37) b_absorb_elastic_bottom(2,i,ispec,it) !old
+ write(37) b_absorb_elastic_bottom(2,i,ispec,it,i_stage) !xiezhinan
enddo
endif
enddo
@@ -5154,14 +5259,17 @@
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)
+! write(38) b_absorb_elastic_top(1,i,ispec,it) !old
+ write(38) b_absorb_elastic_top(1,i,ispec,it,i_stage)!xiezhinan
enddo
do i=1,NGLLX
- write(38) b_absorb_elastic_top(3,i,ispec,it)
+! write(38) b_absorb_elastic_top(3,i,ispec,it) !old
+ write(38) b_absorb_elastic_top(3,i,ispec,it,i_stage)!xiezhinan
enddo
else!SH (membrane) waves
do i=1,NGLLX
- write(38) b_absorb_elastic_top(2,i,ispec,it)
+! write(38) b_absorb_elastic_top(2,i,ispec,it) !old
+ write(38) b_absorb_elastic_top(2,i,ispec,it,i_stage) !xiezhinan
enddo
endif
enddo
@@ -5768,6 +5876,8 @@
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(:)
@@ -5775,6 +5885,21 @@
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)
@@ -6607,223 +6732,223 @@
endif
- if(time_stepping_scheme == 2)then
+! if(time_stepping_scheme == 2)then
! recovering original velocities and accelerations on boundaries (elastic side)
- veloc_elastic = veloc_elastic - beta_LDDRK(i_stage) * veloc_elastic_LDDRK
- displ_elastic = displ_elastic - beta_LDDRK(i_stage) * displ_elastic_LDDRK
- veloc_elastic_LDDRK = (veloc_elastic_LDDRK - deltat * accel_elastic) / alpha_LDDRK(i_stage)
- displ_elastic_LDDRK = (displ_elastic_LDDRK - deltat * veloc_elastic) / alpha_LDDRK(i_stage)
- accel_elastic(1,iglob) = accel_elastic(1,iglob) / rmass_inverse_elastic(iglob)
- accel_elastic(3,iglob) = accel_elastic(3,iglob) / rmass_inverse_elastic(iglob)
+! veloc_elastic = veloc_elastic - beta_LDDRK(i_stage) * veloc_elastic_LDDRK
+! displ_elastic = displ_elastic - beta_LDDRK(i_stage) * displ_elastic_LDDRK
+! veloc_elastic_LDDRK = (veloc_elastic_LDDRK - deltat * accel_elastic) / alpha_LDDRK(i_stage)
+! displ_elastic_LDDRK = (displ_elastic_LDDRK - deltat * veloc_elastic) / alpha_LDDRK(i_stage)
+! accel_elastic(1,iglob) = accel_elastic(1,iglob) / rmass_inverse_elastic(iglob)
+! accel_elastic(3,iglob) = accel_elastic(3,iglob) / rmass_inverse_elastic(iglob)
! recovering original velocities and accelerations on boundaries (poro side)
- velocs_poroelastic = velocs_poroelastic - beta_LDDRK(i_stage) * velocs_poroelastic_LDDRK
- displs_poroelastic = displs_poroelastic - beta_LDDRK(i_stage) * displs_poroelastic_LDDRK
- velocs_poroelastic_LDDRK = (velocs_poroelastic_LDDRK - deltat * accels_poroelastic) / alpha_LDDRK(i_stage)
- displs_poroelastic_LDDRK = (velocs_poroelastic_LDDRK - deltat * velocs_poroelastic) / alpha_LDDRK(i_stage)
- accels_poroelastic(1,iglob) = accels_poroelastic(1,iglob) / rmass_s_inverse_poroelastic(iglob)
- accels_poroelastic(2,iglob) = accels_poroelastic(2,iglob) / rmass_s_inverse_poroelastic(iglob)
+! velocs_poroelastic = velocs_poroelastic - beta_LDDRK(i_stage) * velocs_poroelastic_LDDRK
+! displs_poroelastic = displs_poroelastic - beta_LDDRK(i_stage) * displs_poroelastic_LDDRK
+! velocs_poroelastic_LDDRK = (velocs_poroelastic_LDDRK - deltat * accels_poroelastic) / alpha_LDDRK(i_stage)
+! displs_poroelastic_LDDRK = (velocs_poroelastic_LDDRK - deltat * velocs_poroelastic) / alpha_LDDRK(i_stage)
+! accels_poroelastic(1,iglob) = accels_poroelastic(1,iglob) / rmass_s_inverse_poroelastic(iglob)
+! accels_poroelastic(2,iglob) = accels_poroelastic(2,iglob) / rmass_s_inverse_poroelastic(iglob)
! assembling accelerations
- accel_elastic(1,iglob) = ( accel_elastic(1,iglob) + accels_poroelastic(1,iglob) ) / &
- ( 1.0/rmass_inverse_elastic(iglob) +1.0/rmass_s_inverse_poroelastic(iglob) )
- accel_elastic(3,iglob) = ( accel_elastic(3,iglob) + accels_poroelastic(2,iglob) ) / &
- ( 1.0/rmass_inverse_elastic(iglob) +1.0/rmass_s_inverse_poroelastic(iglob) )
- accels_poroelastic(1,iglob) = accel_elastic(1,iglob)
- accels_poroelastic(2,iglob) = accel_elastic(3,iglob)
+! accel_elastic(1,iglob) = ( accel_elastic(1,iglob) + accels_poroelastic(1,iglob) ) / &
+! ( 1.0/rmass_inverse_elastic(iglob) +1.0/rmass_s_inverse_poroelastic(iglob) )
+! accel_elastic(3,iglob) = ( accel_elastic(3,iglob) + accels_poroelastic(2,iglob) ) / &
+! ( 1.0/rmass_inverse_elastic(iglob) +1.0/rmass_s_inverse_poroelastic(iglob) )
+! accels_poroelastic(1,iglob) = accel_elastic(1,iglob)
+! accels_poroelastic(2,iglob) = accel_elastic(3,iglob)
! updating velocities
! updating velocities(elastic side)
- veloc_elastic_LDDRK = alpha_LDDRK(i_stage) * veloc_elastic_LDDRK + deltat * accel_elastic
- displ_elastic_LDDRK = alpha_LDDRK(i_stage) * displ_elastic_LDDRK + deltat * veloc_elastic
- veloc_elastic = veloc_elastic + beta_LDDRK(i_stage) * veloc_elastic_LDDRK
- displ_elastic = displ_elastic + beta_LDDRK(i_stage) * displ_elastic_LDDRK
+! veloc_elastic_LDDRK = alpha_LDDRK(i_stage) * veloc_elastic_LDDRK + deltat * accel_elastic
+! displ_elastic_LDDRK = alpha_LDDRK(i_stage) * displ_elastic_LDDRK + deltat * veloc_elastic
+! veloc_elastic = veloc_elastic + beta_LDDRK(i_stage) * veloc_elastic_LDDRK
+! displ_elastic = displ_elastic + beta_LDDRK(i_stage) * displ_elastic_LDDRK
! updating velocities(poro side)
- velocs_poroelastic_LDDRK = alpha_LDDRK(i_stage) * velocs_poroelastic_LDDRK + deltat * accels_poroelastic
- displs_poroelastic_LDDRK = alpha_LDDRK(i_stage) * displs_poroelastic_LDDRK + deltat * velocs_poroelastic
- velocs_poroelastic = velocs_poroelastic + beta_LDDRK(i_stage) * velocs_poroelastic_LDDRK
- displs_poroelastic = displs_poroelastic + beta_LDDRK(i_stage) * displs_poroelastic_LDDRK
+! velocs_poroelastic_LDDRK = alpha_LDDRK(i_stage) * velocs_poroelastic_LDDRK + deltat * accels_poroelastic
+! displs_poroelastic_LDDRK = alpha_LDDRK(i_stage) * displs_poroelastic_LDDRK + deltat * velocs_poroelastic
+! velocs_poroelastic = velocs_poroelastic + beta_LDDRK(i_stage) * velocs_poroelastic_LDDRK
+! displs_poroelastic = displs_poroelastic + beta_LDDRK(i_stage) * displs_poroelastic_LDDRK
! zeros w
- accelw_poroelastic(1,iglob) = ZERO
- accelw_poroelastic(2,iglob) = ZERO
- velocw_poroelastic(1,iglob) = ZERO
- velocw_poroelastic(2,iglob) = ZERO
- endif
+! accelw_poroelastic(1,iglob) = ZERO
+! accelw_poroelastic(2,iglob) = ZERO
+! velocw_poroelastic(1,iglob) = ZERO
+! velocw_poroelastic(2,iglob) = ZERO
+! endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
- if(time_stepping_scheme == 3)then
+! if(time_stepping_scheme == 3)then
! recovering original velocities and accelerations on boundaries (elastic side)
- if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then
+! 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)weight_rk = 0.5d0
+! if(i_stage == 2)weight_rk = 0.5d0
+! if(i_stage == 3)weight_rk = 1.0d0
- veloc_elastic(1,iglob) = veloc_elastic_initial_rk(1,iglob) - weight_rk * accel_elastic_rk(1,iglob,i_stage)
- veloc_elastic(3,iglob) = veloc_elastic_initial_rk(3,iglob) - weight_rk * accel_elastic_rk(3,iglob,i_stage)
- displ_elastic(1,iglob) = displ_elastic_initial_rk(1,iglob) - weight_rk * veloc_elastic_rk(1,iglob,i_stage)
- displ_elastic(3,iglob) = displ_elastic_initial_rk(3,iglob) - weight_rk * veloc_elastic_rk(3,iglob,i_stage)
+! veloc_elastic(1,iglob) = veloc_elastic_initial_rk(1,iglob) - weight_rk * accel_elastic_rk(1,iglob,i_stage)
+! veloc_elastic(3,iglob) = veloc_elastic_initial_rk(3,iglob) - weight_rk * accel_elastic_rk(3,iglob,i_stage)
+! displ_elastic(1,iglob) = displ_elastic_initial_rk(1,iglob) - weight_rk * veloc_elastic_rk(1,iglob,i_stage)
+! displ_elastic(3,iglob) = displ_elastic_initial_rk(3,iglob) - weight_rk * veloc_elastic_rk(3,iglob,i_stage)
- elseif(i_stage==4)then
+! elseif(i_stage==4)then
- veloc_elastic(1,iglob) = veloc_elastic_initial_rk(1,iglob) - 1.0d0 / 6.0d0 * &
- (accel_elastic_rk(1,iglob,1) + 2.0d0 * accel_elastic_rk(1,iglob,2) + &
- 2.0d0 * accel_elastic_rk(1,iglob,3) + accel_elastic_rk(1,iglob,4))
+! veloc_elastic(1,iglob) = veloc_elastic_initial_rk(1,iglob) - 1.0d0 / 6.0d0 * &
+! (accel_elastic_rk(1,iglob,1) + 2.0d0 * accel_elastic_rk(1,iglob,2) + &
+! 2.0d0 * accel_elastic_rk(1,iglob,3) + accel_elastic_rk(1,iglob,4))
- veloc_elastic(3,iglob) = veloc_elastic_initial_rk(3,iglob) - 1.0d0 / 6.0d0 * &
- (accel_elastic_rk(3,iglob,1) + 2.0d0 * accel_elastic_rk(3,iglob,2) + &
- 2.0d0 * accel_elastic_rk(3,iglob,3) + accel_elastic_rk(3,iglob,4))
+! veloc_elastic(3,iglob) = veloc_elastic_initial_rk(3,iglob) - 1.0d0 / 6.0d0 * &
+! (accel_elastic_rk(3,iglob,1) + 2.0d0 * accel_elastic_rk(3,iglob,2) + &
+! 2.0d0 * accel_elastic_rk(3,iglob,3) + accel_elastic_rk(3,iglob,4))
- displ_elastic(1,iglob) = displ_elastic_initial_rk(1,iglob) - 1.0d0 / 6.0d0 * &
- (veloc_elastic_rk(1,iglob,1) + 2.0d0 * veloc_elastic_rk(1,iglob,2) + &
- 2.0d0 * veloc_elastic_rk(1,iglob,3) + veloc_elastic_rk(1,iglob,4))
+! displ_elastic(1,iglob) = displ_elastic_initial_rk(1,iglob) - 1.0d0 / 6.0d0 * &
+! (veloc_elastic_rk(1,iglob,1) + 2.0d0 * veloc_elastic_rk(1,iglob,2) + &
+! 2.0d0 * veloc_elastic_rk(1,iglob,3) + veloc_elastic_rk(1,iglob,4))
- displ_elastic(3,iglob) = displ_elastic_initial_rk(3,iglob) - 1.0d0 / 6.0d0 * &
- (veloc_elastic_rk(3,iglob,1) + 2.0d0 * veloc_elastic_rk(3,iglob,2) + &
- 2.0d0 * veloc_elastic_rk(3,iglob,3) + veloc_elastic_rk(3,iglob,4))
+! displ_elastic(3,iglob) = displ_elastic_initial_rk(3,iglob) - 1.0d0 / 6.0d0 * &
+! (veloc_elastic_rk(3,iglob,1) + 2.0d0 * veloc_elastic_rk(3,iglob,2) + &
+! 2.0d0 * veloc_elastic_rk(3,iglob,3) + veloc_elastic_rk(3,iglob,4))
- endif
+! endif
- accel_elastic(1,iglob) = accel_elastic(1,iglob) / rmass_inverse_elastic(iglob)
- accel_elastic(3,iglob) = accel_elastic(3,iglob) / rmass_inverse_elastic(iglob)
+! accel_elastic(1,iglob) = accel_elastic(1,iglob) / rmass_inverse_elastic(iglob)
+! accel_elastic(3,iglob) = accel_elastic(3,iglob) / rmass_inverse_elastic(iglob)
- accel_elastic_rk(1,iglob,i_stage) = accel_elastic(1,iglob) / deltat
- accel_elastic_rk(3,iglob,i_stage) = accel_elastic(3,iglob) / deltat
- veloc_elastic_rk(1,iglob,i_stage) = veloc_elastic(1,iglob) / deltat
- veloc_elastic_rk(3,iglob,i_stage) = veloc_elastic(3,iglob) / deltat
+! accel_elastic_rk(1,iglob,i_stage) = accel_elastic(1,iglob) / deltat
+! accel_elastic_rk(3,iglob,i_stage) = accel_elastic(3,iglob) / deltat
+! veloc_elastic_rk(1,iglob,i_stage) = veloc_elastic(1,iglob) / deltat
+! veloc_elastic_rk(3,iglob,i_stage) = veloc_elastic(3,iglob) / deltat
! recovering original velocities and accelerations on boundaries (poro side)
- if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then
+! 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)weight_rk = 0.5d0
+! if(i_stage == 2)weight_rk = 0.5d0
+! if(i_stage == 3)weight_rk = 1.0d0
- velocs_poroelastic(1,iglob) = velocs_poroelastic_initial_rk(1,iglob) - weight_rk * accels_poroelastic_rk(1,iglob,i_stage)
- velocs_poroelastic(2,iglob) = velocs_poroelastic_initial_rk(2,iglob) - weight_rk * accels_poroelastic_rk(2,iglob,i_stage)
- displs_poroelastic(1,iglob) = displs_poroelastic_initial_rk(1,iglob) - weight_rk * velocs_poroelastic_rk(1,iglob,i_stage)
- displs_poroelastic(2,iglob) = displs_poroelastic_initial_rk(2,iglob) - weight_rk * velocs_poroelastic_rk(2,iglob,i_stage)
+! velocs_poroelastic(1,iglob) = velocs_poroelastic_initial_rk(1,iglob) - weight_rk * accels_poroelastic_rk(1,iglob,i_stage)
+! velocs_poroelastic(2,iglob) = velocs_poroelastic_initial_rk(2,iglob) - weight_rk * accels_poroelastic_rk(2,iglob,i_stage)
+! displs_poroelastic(1,iglob) = displs_poroelastic_initial_rk(1,iglob) - weight_rk * velocs_poroelastic_rk(1,iglob,i_stage)
+! displs_poroelastic(2,iglob) = displs_poroelastic_initial_rk(2,iglob) - weight_rk * velocs_poroelastic_rk(2,iglob,i_stage)
- elseif(i_stage==4)then
+! elseif(i_stage==4)then
- velocs_poroelastic(1,iglob) = velocs_poroelastic_initial_rk(1,iglob) - 1.0d0 / 6.0d0 * &
- (accels_poroelastic_rk(1,iglob,1) + 2.0d0 * accels_poroelastic_rk(1,iglob,2) + &
- 2.0d0 * accels_poroelastic_rk(1,iglob,3) + accels_poroelastic_rk(1,iglob,4))
+! velocs_poroelastic(1,iglob) = velocs_poroelastic_initial_rk(1,iglob) - 1.0d0 / 6.0d0 * &
+! (accels_poroelastic_rk(1,iglob,1) + 2.0d0 * accels_poroelastic_rk(1,iglob,2) + &
+! 2.0d0 * accels_poroelastic_rk(1,iglob,3) + accels_poroelastic_rk(1,iglob,4))
- velocs_poroelastic(2,iglob) = velocs_poroelastic_initial_rk(2,iglob) - 1.0d0 / 6.0d0 * &
- (accels_poroelastic_rk(2,iglob,1) + 2.0d0 * accels_poroelastic_rk(2,iglob,2) + &
- 2.0d0 * accels_poroelastic_rk(2,iglob,3) + accels_poroelastic_rk(2,iglob,4))
+! velocs_poroelastic(2,iglob) = velocs_poroelastic_initial_rk(2,iglob) - 1.0d0 / 6.0d0 * &
+! (accels_poroelastic_rk(2,iglob,1) + 2.0d0 * accels_poroelastic_rk(2,iglob,2) + &
+! 2.0d0 * accels_poroelastic_rk(2,iglob,3) + accels_poroelastic_rk(2,iglob,4))
- displs_poroelastic(1,iglob) = displs_poroelastic_initial_rk(1,iglob) - 1.0d0 / 6.0d0 * &
- (velocs_poroelastic_rk(1,iglob,1) + 2.0d0 * velocs_poroelastic_rk(1,iglob,2) + &
- 2.0d0 * velocs_poroelastic_rk(1,iglob,3) + velocs_poroelastic_rk(1,iglob,4))
+! displs_poroelastic(1,iglob) = displs_poroelastic_initial_rk(1,iglob) - 1.0d0 / 6.0d0 * &
+! (velocs_poroelastic_rk(1,iglob,1) + 2.0d0 * velocs_poroelastic_rk(1,iglob,2) + &
+! 2.0d0 * velocs_poroelastic_rk(1,iglob,3) + velocs_poroelastic_rk(1,iglob,4))
- displs_poroelastic(2,iglob) = displs_poroelastic_initial_rk(2,iglob) - 1.0d0 / 6.0d0 * &
- (velocs_poroelastic_rk(2,iglob,1) + 2.0d0 * velocs_poroelastic_rk(2,iglob,2) + &
- 2.0d0 * velocs_poroelastic_rk(2,iglob,3) + velocs_poroelastic_rk(2,iglob,4))
+! displs_poroelastic(2,iglob) = displs_poroelastic_initial_rk(2,iglob) - 1.0d0 / 6.0d0 * &
+! (velocs_poroelastic_rk(2,iglob,1) + 2.0d0 * velocs_poroelastic_rk(2,iglob,2) + &
+! 2.0d0 * velocs_poroelastic_rk(2,iglob,3) + velocs_poroelastic_rk(2,iglob,4))
- endif
+! endif
- accels_poroelastic(1,iglob) = accels_poroelastic(1,iglob) / rmass_s_inverse_poroelastic(iglob)
- accels_poroelastic(2,iglob) = accels_poroelastic(2,iglob) / rmass_s_inverse_poroelastic(iglob)
+! accels_poroelastic(1,iglob) = accels_poroelastic(1,iglob) / rmass_s_inverse_poroelastic(iglob)
+! accels_poroelastic(2,iglob) = accels_poroelastic(2,iglob) / rmass_s_inverse_poroelastic(iglob)
- accels_poroelastic_rk(1,iglob,i_stage) = accels_poroelastic(1,iglob) / deltat
- accels_poroelastic_rk(2,iglob,i_stage) = accels_poroelastic(2,iglob) / deltat
- velocs_poroelastic_rk(1,iglob,i_stage) = velocs_poroelastic(1,iglob) / deltat
- velocs_poroelastic_rk(2,iglob,i_stage) = velocs_poroelastic(2,iglob) / deltat
+! accels_poroelastic_rk(1,iglob,i_stage) = accels_poroelastic(1,iglob) / deltat
+! accels_poroelastic_rk(2,iglob,i_stage) = accels_poroelastic(2,iglob) / deltat
+! velocs_poroelastic_rk(1,iglob,i_stage) = velocs_poroelastic(1,iglob) / deltat
+! velocs_poroelastic_rk(2,iglob,i_stage) = velocs_poroelastic(2,iglob) / deltat
! assembling accelerations
- accel_elastic(1,iglob) = ( accel_elastic(1,iglob) + accels_poroelastic(1,iglob) ) / &
- ( 1.0/rmass_inverse_elastic(iglob) +1.0/rmass_s_inverse_poroelastic(iglob) )
- accel_elastic(3,iglob) = ( accel_elastic(3,iglob) + accels_poroelastic(2,iglob) ) / &
- ( 1.0/rmass_inverse_elastic(iglob) +1.0/rmass_s_inverse_poroelastic(iglob) )
- accels_poroelastic(1,iglob) = accel_elastic(1,iglob)
- accels_poroelastic(2,iglob) = accel_elastic(3,iglob)
+! accel_elastic(1,iglob) = ( accel_elastic(1,iglob) + accels_poroelastic(1,iglob) ) / &
+! ( 1.0/rmass_inverse_elastic(iglob) +1.0/rmass_s_inverse_poroelastic(iglob) )
+! accel_elastic(3,iglob) = ( accel_elastic(3,iglob) + accels_poroelastic(2,iglob) ) / &
+! ( 1.0/rmass_inverse_elastic(iglob) +1.0/rmass_s_inverse_poroelastic(iglob) )
+! accels_poroelastic(1,iglob) = accel_elastic(1,iglob)
+! accels_poroelastic(2,iglob) = accel_elastic(3,iglob)
- ! updating velocities
+ ! updating velocities
! updating velocities(elastic side)
- accel_elastic_rk(1,iglob,i_stage) = accel_elastic(1,iglob) * deltat
- accel_elastic_rk(3,iglob,i_stage) = accel_elastic(3,iglob) * deltat
+ ! accel_elastic_rk(1,iglob,i_stage) = accel_elastic(1,iglob) * deltat
+ ! accel_elastic_rk(3,iglob,i_stage) = accel_elastic(3,iglob) * deltat
- if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then
+ ! 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)weight_rk = 0.5d0
+ ! if(i_stage == 2)weight_rk = 0.5d0
+ ! if(i_stage == 3)weight_rk = 1.0d0
- veloc_elastic(1,iglob) = veloc_elastic_initial_rk(1,iglob) + weight_rk * accel_elastic_rk(1,iglob,i_stage)
- veloc_elastic(3,iglob) = veloc_elastic_initial_rk(3,iglob) + weight_rk * accel_elastic_rk(3,iglob,i_stage)
- displ_elastic(1,iglob) = displ_elastic_initial_rk(1,iglob) + weight_rk * veloc_elastic_rk(1,iglob,i_stage)
- displ_elastic(3,iglob) = displ_elastic_initial_rk(3,iglob) + weight_rk * veloc_elastic_rk(3,iglob,i_stage)
+ ! veloc_elastic(1,iglob) = veloc_elastic_initial_rk(1,iglob) + weight_rk * accel_elastic_rk(1,iglob,i_stage)
+ ! veloc_elastic(3,iglob) = veloc_elastic_initial_rk(3,iglob) + weight_rk * accel_elastic_rk(3,iglob,i_stage)
+ ! displ_elastic(1,iglob) = displ_elastic_initial_rk(1,iglob) + weight_rk * veloc_elastic_rk(1,iglob,i_stage)
+ ! displ_elastic(3,iglob) = displ_elastic_initial_rk(3,iglob) + weight_rk * veloc_elastic_rk(3,iglob,i_stage)
- elseif(i_stage==4)then
+ ! elseif(i_stage==4)then
- veloc_elastic(1,iglob) = veloc_elastic_initial_rk(1,iglob) + 1.0d0 / 6.0d0 * &
- (accel_elastic_rk(1,iglob,1) + 2.0d0 * accel_elastic_rk(1,iglob,2) + &
- 2.0d0 * accel_elastic_rk(1,iglob,3) + accel_elastic_rk(1,iglob,4))
+ ! veloc_elastic(1,iglob) = veloc_elastic_initial_rk(1,iglob) + 1.0d0 / 6.0d0 * &
+ ! (accel_elastic_rk(1,iglob,1) + 2.0d0 * accel_elastic_rk(1,iglob,2) + &
+ ! 2.0d0 * accel_elastic_rk(1,iglob,3) + accel_elastic_rk(1,iglob,4))
+!
+ ! veloc_elastic(3,iglob) = veloc_elastic_initial_rk(3,iglob) + 1.0d0 / 6.0d0 * &
+ ! (accel_elastic_rk(3,iglob,1) + 2.0d0 * accel_elastic_rk(3,iglob,2) + &
+ ! 2.0d0 * accel_elastic_rk(3,iglob,3) + accel_elastic_rk(3,iglob,4))
- veloc_elastic(3,iglob) = veloc_elastic_initial_rk(3,iglob) + 1.0d0 / 6.0d0 * &
- (accel_elastic_rk(3,iglob,1) + 2.0d0 * accel_elastic_rk(3,iglob,2) + &
- 2.0d0 * accel_elastic_rk(3,iglob,3) + accel_elastic_rk(3,iglob,4))
+ ! displ_elastic(1,iglob) = displ_elastic_initial_rk(1,iglob) + 1.0d0 / 6.0d0 * &
+ ! (veloc_elastic_rk(1,iglob,1) + 2.0d0 * veloc_elastic_rk(1,iglob,2) + &
+ ! 2.0d0 * veloc_elastic_rk(1,iglob,3) + veloc_elastic_rk(1,iglob,4))
- displ_elastic(1,iglob) = displ_elastic_initial_rk(1,iglob) + 1.0d0 / 6.0d0 * &
- (veloc_elastic_rk(1,iglob,1) + 2.0d0 * veloc_elastic_rk(1,iglob,2) + &
- 2.0d0 * veloc_elastic_rk(1,iglob,3) + veloc_elastic_rk(1,iglob,4))
+ ! displ_elastic(3,iglob) = displ_elastic_initial_rk(3,iglob) + 1.0d0 / 6.0d0 * &
+ ! (veloc_elastic_rk(3,iglob,1) + 2.0d0 * veloc_elastic_rk(3,iglob,2) + &
+ ! 2.0d0 * veloc_elastic_rk(3,iglob,3) + veloc_elastic_rk(3,iglob,4))
- displ_elastic(3,iglob) = displ_elastic_initial_rk(3,iglob) + 1.0d0 / 6.0d0 * &
- (veloc_elastic_rk(3,iglob,1) + 2.0d0 * veloc_elastic_rk(3,iglob,2) + &
- 2.0d0 * veloc_elastic_rk(3,iglob,3) + veloc_elastic_rk(3,iglob,4))
-
- endif
+ ! endif
! updating velocities(poro side)
- accels_poroelastic_rk(1,iglob,i_stage) = deltat * accels_poroelastic(1,iglob)
- accels_poroelastic_rk(2,iglob,i_stage) = deltat * accels_poroelastic(2,iglob)
- velocs_poroelastic_rk(1,iglob,i_stage) = deltat * velocs_poroelastic(1,iglob)
- velocs_poroelastic_rk(2,iglob,i_stage) = deltat * velocs_poroelastic(2,iglob)
+ ! accels_poroelastic_rk(1,iglob,i_stage) = deltat * accels_poroelastic(1,iglob)
+ ! accels_poroelastic_rk(2,iglob,i_stage) = deltat * accels_poroelastic(2,iglob)
+ ! velocs_poroelastic_rk(1,iglob,i_stage) = deltat * velocs_poroelastic(1,iglob)
+ ! velocs_poroelastic_rk(2,iglob,i_stage) = deltat * velocs_poroelastic(2,iglob)
- if(i_stage==1 .or. i_stage==2 .or. i_stage==3)then
+ ! 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)weight_rk = 0.5d0
+ ! if(i_stage == 2)weight_rk = 0.5d0
+ ! if(i_stage == 3)weight_rk = 1.0d0
- velocs_poroelastic(1,iglob) = velocs_poroelastic_initial_rk(1,iglob) + weight_rk * accels_poroelastic_rk(1,iglob,i_stage)
- velocs_poroelastic(2,iglob) = velocs_poroelastic_initial_rk(2,iglob) + weight_rk * accels_poroelastic_rk(2,iglob,i_stage)
- displs_poroelastic(1,iglob) = displs_poroelastic_initial_rk(1,iglob) + weight_rk * velocs_poroelastic_rk(1,iglob,i_stage)
- displs_poroelastic(2,iglob) = displs_poroelastic_initial_rk(2,iglob) + weight_rk * velocs_poroelastic_rk(2,iglob,i_stage)
+ ! velocs_poroelastic(1,iglob) = velocs_poroelastic_initial_rk(1,iglob) + weight_rk * accels_poroelastic_rk(1,iglob,i_stage)
+ ! velocs_poroelastic(2,iglob) = velocs_poroelastic_initial_rk(2,iglob) + weight_rk * accels_poroelastic_rk(2,iglob,i_stage)
+ ! displs_poroelastic(1,iglob) = displs_poroelastic_initial_rk(1,iglob) + weight_rk * velocs_poroelastic_rk(1,iglob,i_stage)
+ ! displs_poroelastic(2,iglob) = displs_poroelastic_initial_rk(2,iglob) + weight_rk * velocs_poroelastic_rk(2,iglob,i_stage)
- elseif(i_stage==4)then
+ ! elseif(i_stage==4)then
- velocs_poroelastic(1,iglob) = velocs_poroelastic_initial_rk(1,iglob) + 1.0d0 / 6.0d0 * &
- (accels_poroelastic_rk(1,iglob,1) + 2.0d0 * accels_poroelastic_rk(1,iglob,2) + &
- 2.0d0 * accels_poroelastic_rk(1,iglob,3) + accels_poroelastic_rk(1,iglob,4))
+ ! velocs_poroelastic(1,iglob) = velocs_poroelastic_initial_rk(1,iglob) + 1.0d0 / 6.0d0 * &
+ ! (accels_poroelastic_rk(1,iglob,1) + 2.0d0 * accels_poroelastic_rk(1,iglob,2) + &
+ ! 2.0d0 * accels_poroelastic_rk(1,iglob,3) + accels_poroelastic_rk(1,iglob,4))
- velocs_poroelastic(2,iglob) = velocs_poroelastic_initial_rk(2,iglob) + 1.0d0 / 6.0d0 * &
- (accels_poroelastic_rk(2,iglob,1) + 2.0d0 * accels_poroelastic_rk(2,iglob,2) + &
- 2.0d0 * accels_poroelastic_rk(2,iglob,3) + accels_poroelastic_rk(2,iglob,4))
+ ! velocs_poroelastic(2,iglob) = velocs_poroelastic_initial_rk(2,iglob) + 1.0d0 / 6.0d0 * &
+ ! (accels_poroelastic_rk(2,iglob,1) + 2.0d0 * accels_poroelastic_rk(2,iglob,2) + &
+ ! 2.0d0 * accels_poroelastic_rk(2,iglob,3) + accels_poroelastic_rk(2,iglob,4))
+!
+ ! displs_poroelastic(1,iglob) = displs_poroelastic_initial_rk(1,iglob) + 1.0d0 / 6.0d0 * &
+ ! (velocs_poroelastic_rk(1,iglob,1) + 2.0d0 * velocs_poroelastic_rk(1,iglob,2) + &
+ ! 2.0d0 * velocs_poroelastic_rk(1,iglob,3) + velocs_poroelastic_rk(1,iglob,4))
+!
+ ! displs_poroelastic(2,iglob) = displs_poroelastic_initial_rk(2,iglob) + 1.0d0 / 6.0d0 * &
+ ! (velocs_poroelastic_rk(2,iglob,1) + 2.0d0 * velocs_poroelastic_rk(2,iglob,2) + &
+ ! 2.0d0 * velocs_poroelastic_rk(2,iglob,3) + velocs_poroelastic_rk(2,iglob,4))
- displs_poroelastic(1,iglob) = displs_poroelastic_initial_rk(1,iglob) + 1.0d0 / 6.0d0 * &
- (velocs_poroelastic_rk(1,iglob,1) + 2.0d0 * velocs_poroelastic_rk(1,iglob,2) + &
- 2.0d0 * velocs_poroelastic_rk(1,iglob,3) + velocs_poroelastic_rk(1,iglob,4))
+ ! endif
- displs_poroelastic(2,iglob) = displs_poroelastic_initial_rk(2,iglob) + 1.0d0 / 6.0d0 * &
- (velocs_poroelastic_rk(2,iglob,1) + 2.0d0 * velocs_poroelastic_rk(2,iglob,2) + &
- 2.0d0 * velocs_poroelastic_rk(2,iglob,3) + velocs_poroelastic_rk(2,iglob,4))
-
- endif
-
- endif
+ ! endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
if(SIMULATION_TYPE == 2) then
@@ -6862,7 +6987,7 @@
enddo
endif
- enddo !LDDRK or RK
+ ! enddo !LDDRK or RK !initially set the end of LDDRK
! ********************************************************************************************
! reading lastframe for adjoint/kernels calculation
@@ -6961,68 +7086,8 @@
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
@@ -7158,21 +7223,80 @@
! 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,irecloc) = cosrot_irec(irecloc)*valux + sinrot_irec(irecloc)*valuz
- sisuz(seismo_current,irecloc) = - sinrot_irec(irecloc)*valux + cosrot_irec(irecloc)*valuz
+ 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
else
- sisux(seismo_current,irecloc) = valuy
- sisuz(seismo_current,irecloc) = ZERO
+ sisux(seismo_current,i_stage,irecloc) = valuy !xiezhinan
+ sisuz(seismo_current,i_stage,irecloc) = ZERO !xiezhinan
endif
else
- sisux(seismo_current,irecloc) = valux
- sisuz(seismo_current,irecloc) = ZERO
+ sisux(seismo_current,i_stage,irecloc) = valux !xiezhinan
+ sisuz(seismo_current,i_stage,irecloc) = ZERO !xiezhinan
endif
- siscurl(seismo_current,irecloc) = valcurl
+ siscurl(seismo_current,i_stage,irecloc) = valcurl !xiezhinan
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
@@ -7890,8 +8014,13 @@
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)
+ st_zval,x_source(1),z_source(1),SU_FORMAT,stage_time_scheme,SIMULATION_TYPE,c_LDDRK) !xiezhinan
+! 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-07 01:19:09 UTC (rev 19591)
+++ seismo/2D/SPECFEM2D/branches/Xie_Zhinan/src/specfem2D/write_seismograms.F90 2012-02-07 16:23:33 UTC (rev 19592)
@@ -44,10 +44,15 @@
! 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)
+ st_zval,x_source,z_source,SU_FORMAT,stage_time_scheme,SIMULATION_TYPE,c_LDDRK) !xiezhinan
implicit none
@@ -58,6 +63,9 @@
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)
@@ -68,8 +76,10 @@
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
+! double precision, dimension(NTSTEP_BETWEEN_OUTPUT_SEISMO,nrecloc), intent(in) :: sisux,sisuz,siscurl !old
+ 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
@@ -82,7 +92,8 @@
character(len=150) sisname
! to write seismograms in single precision SEP and double precision binary format
- double precision, dimension(:,:), allocatable :: buffer_binary
+! double precision, dimension(:,:), allocatable :: buffer_binary !old
+ double precision, dimension(:,:,:), allocatable :: buffer_binary
! scaling factor for Seismic Unix xsu dislay
double precision, parameter :: FACTORXSU = 1.d0
@@ -129,7 +140,8 @@
number_of_components = NDIM
endif
- allocate(buffer_binary(NTSTEP_BETWEEN_OUTPUT_SEISMO,number_of_components))
+! 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
if ( myrank == 0 .and. seismo_offset == 0 ) then
@@ -204,12 +216,17 @@
if ( which_proc_receiver(irec) == myrank ) then
irecloc = irecloc + 1
- buffer_binary(:,1) = sisux(:,irecloc)
+! buffer_binary(:,1) = sisux(:,irecloc) !old
+ buffer_binary(:,:,1) = sisux(:,:,irecloc) !xiezhinan
if ( number_of_components == 2 ) then
- buffer_binary(:,2) = sisuz(:,irecloc)
+! buffer_binary(:,2) = sisuz(:,irecloc) !old
+ buffer_binary(:,:,2) = sisuz(:,:,irecloc) !xiezhinan
else if ( number_of_components == 3 ) then
- buffer_binary(:,2) = sisuz(:,irecloc)
- buffer_binary(:,3) = siscurl(:,irecloc)
+! buffer_binary(:,2) = sisuz(:,irecloc) !old
+! buffer_binary(:,3) = siscurl(:,irecloc) !old
+
+ buffer_binary(:,:,2) = sisuz(:,:,irecloc) !xiezhinan
+ buffer_binary(:,:,3) = siscurl(:,:,irecloc) !xiezhinan
end if
#ifdef USE_MPI
@@ -226,6 +243,8 @@
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
@@ -278,71 +297,179 @@
! 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),' ', &
- sngl(buffer_binary(isample,iorientation))
+! 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
+
else
- write(11,*) sngl(dble(seismo_offset+isample-1)*deltat - t0),' ', &
- sngl(buffer_binary(isample,iorientation))
+! 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
+
endif
+ enddo
enddo
close(11)
end do
! write binary seismogram
do isample = 1, seismo_current
- 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)
+ 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
+
if ( seismotype /= 4 .and. seismotype /= 6 .and. p_sv) then
- 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)
+! 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
+
end if
+
if ( seismotype == 5 ) then
- 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)
+! 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
+
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
- 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
+! 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
+
header2(1)=0 ! dummy
- 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
+! 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
+
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
- 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
+! 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
+
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
- 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
+ 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
enddo
endif
More information about the CIG-COMMITS
mailing list