[cig-commits] r21487 - seismo/2D/SPECFEM2D/trunk/src/specfem2D
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Mon Mar 11 10:45:04 PDT 2013
Author: dkomati1
Date: 2013-03-11 10:45:03 -0700 (Mon, 11 Mar 2013)
New Revision: 21487
Modified:
seismo/2D/SPECFEM2D/trunk/src/specfem2D/calendar.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_poro_fluid.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_poro_solid.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/read_databases.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
Log:
added the full reference of Wang and Texeira (2006) for C-PML;
also removed useless white spaces
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/calendar.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/calendar.f90 2013-03-11 16:22:35 UTC (rev 21486)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/calendar.f90 2013-03-11 17:45:03 UTC (rev 21487)
@@ -429,10 +429,10 @@
!
! Explanation of all internal variables.
! jdref Julian Day on which 1 March begins in the reference year.
-! jmonth Month counter which equals month+1 if month .gt. 2
-! or month+13 if month .le. 2.
-! jyear Year index, jyear=iyear if iyear .gt. 0, jyear=iyear+1
-! if iyear .lt. 0. Thus, jyear does not skip year 0
+! jmonth Month counter which equals month+1 if month > 2
+! or month+13 if month <= 2.
+! jyear Year index, jyear=iyear if iyear > 0, jyear=iyear+1
+! if iyear < 0. Thus, jyear does not skip year 0
! like iyear does between BC and AD years.
! leap =1 if the year is a leap year, =0 if not.
! n1yr Number of complete individual years between iyear and
@@ -456,7 +456,7 @@
! this is 100*365 + 25 = 36525.
! nyrs Number of years from the beginning of yr400
! to the beginning of jyear. (Used for option +/-3).
-! yr400 The largest multiple of 400 years that is .le. jyear.
+! yr400 The largest multiple of 400 years that is <= jyear.
!
!
!----------------------------------------------------------------
@@ -590,7 +590,7 @@
jmonth = month + 1
endif
!
-! Find the closest multiple of 400 years that is .le. jyear.
+! Find the closest multiple of 400 years that is <= jyear.
yr400 = (jyear/400)*400
! = multiple of 400 years at or less than jyear.
if (jyear < yr400) then
@@ -717,7 +717,7 @@
endif
endif
!
-! Adjust the year if it is .le. 0, and hence BC (Before Christ).
+! Adjust the year if it is <= 0, and hence BC (Before Christ).
if (iyear <= 0) then
iyear = iyear - 1
endif
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90 2013-03-11 16:22:35 UTC (rev 21486)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_acoustic.f90 2013-03-11 17:45:03 UTC (rev 21487)
@@ -465,7 +465,7 @@
cpl = vpext(i,j,ispec)
!assuming that in fluid(acoustic) part input cpl is defined by sqrt(kappal/rhol), &
!which is not the same as in cpl input in elastic part
- kappal = rhol*cpl*cpl
+ kappal = rhol*cpl*cpl
else
if(CUSTOM_REAL == SIZE_REAL) then
lambdal_relaxed = sngl(poroelastcoef(1,1,kmato(ispec)))
@@ -594,7 +594,7 @@
rmemory_potential_acoustic(1,i,j,ispec_PML) = rmemory_potential_acoustic(1,i,j,ispec_PML) + &
beta_LDDRK(i_stage) * rmemory_potential_acoust_LDDRK(1,i,j,ispec_PML)
rmemory_potential_acoustic(2,i,j,ispec_PML) =0._CUSTOM_REAL
- end if
+ end if
potential_dot_dot_acoustic_PML(i,j,ispec_PML)= wxgll(i)*wzgll(j)/ kappal*jacobian(i,j,ispec) * &
(A0 * potential_acoustic(iglob) + &
@@ -622,14 +622,14 @@
A4 = -alpha_x_store(i,j,ispec_PML) ** 2*d_x_store(i,j,ispec_PML)*d_z_store(i,j,ispec_PML)
- if(stage_time_scheme == 6) then
+ if(stage_time_scheme == 6) then
A3 = alpha_x_store(i,j,ispec_PML) ** 2*(d_x_store(i,j,ispec_PML) * k_z_store(i,j,ispec_PML)+ &
d_z_store(i,j,ispec_PML) * k_x_store(i,j,ispec_PML)) &
-2._CUSTOM_REAL * alpha_x_store(i,j,ispec_PML)*d_x_store(i,j,ispec_PML)*d_z_store(i,j,ispec_PML)
A4 = alpha_x_store(i,j,ispec_PML) ** 2*d_x_store(i,j,ispec_PML)*d_z_store(i,j,ispec_PML)
- end if
+ end if
- if(stage_time_scheme == 6) then
+ if(stage_time_scheme == 6) then
bb = alpha_z_store(i,j,ispec_PML)
if(bb < 0.0)then
@@ -683,7 +683,7 @@
rmemory_potential_acoustic(2,i,j,ispec_PML) = rmemory_potential_acoustic(2,i,j,ispec_PML) + &
beta_LDDRK(i_stage) * rmemory_potential_acoust_LDDRK(2,i,j,ispec_PML)
- end if
+ end if
potential_dot_dot_acoustic_PML(i,j,ispec_PML)= wxgll(i)*wzgll(j)/ kappal*jacobian(i,j,ispec) * &
(A0 * potential_acoustic(iglob) + &
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_poro_fluid.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_poro_fluid.f90 2013-03-11 16:22:35 UTC (rev 21486)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_poro_fluid.f90 2013-03-11 17:45:03 UTC (rev 21487)
@@ -45,12 +45,12 @@
ispec_selected_source,ispec_selected_rec,is_proc_source,which_proc_receiver,&
source_type,it,NSTEP,anyabs, &
initialfield,ATTENUATION_VISCOELASTIC_SOLID,ATTENUATION_PORO_FLUID_PART,deltat, &
- deltatover2,deltatsquareover2,ibool,kmato,numabs,poroelastic,codeabs, &
+ deltatover2,deltatsquareover2,ibool,kmato,numabs,poroelastic,codeabs, &
accelw_poroelastic,velocw_poroelastic,displw_poroelastic,velocs_poroelastic,displs_poroelastic,&
b_accelw_poroelastic,b_displw_poroelastic,b_displs_poroelastic,&
density,porosity,tortuosity,permeability,poroelastcoef,xix,xiz,gammax,gammaz, &
jacobian,source_time_function,sourcearray,adj_sourcearrays,e11, &
- e13,e11_veloc,e13_veloc,e11_accel,e13_accel,dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n, &
+ e13,e11_veloc,e13_veloc,e11_accel,e13_accel,dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n, &
dvx_dxl_n,dvz_dzl_n,dvz_dxl_n,dvx_dzl_n,hprime_xx,hprimewgll_xx, &
hprime_zz,hprimewgll_zz,wxgll,wzgll,inv_tau_sigma_nu2,&
phi_nu2,Mu_nu2,N_SLS, &
@@ -115,8 +115,8 @@
integer :: N_SLS
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e11,e13
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e11_veloc,e13_veloc
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e11_accel,e13_accel
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e11_veloc,e13_veloc
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e11_accel,e13_accel
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e11_LDDRK,e13_LDDRK
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e11_initial_rk,e13_initial_rk
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS,stage_time_scheme) :: e11_force_RK, e13_force_RK
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_poro_solid.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_poro_solid.f90 2013-03-11 16:22:35 UTC (rev 21486)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_poro_solid.f90 2013-03-11 17:45:03 UTC (rev 21487)
@@ -44,14 +44,14 @@
subroutine compute_forces_poro_solid(nglob,nspec,myrank,nelemabs,numat, &
ispec_selected_source,ispec_selected_rec,is_proc_source,which_proc_receiver,&
source_type,it,NSTEP,anyabs, &
- initialfield,ATTENUATION_VISCOELASTIC_SOLID,ATTENUATION_PORO_FLUID_PART,deltat, &
- deltatover2,deltatsquareover2,ibool,kmato,numabs,poroelastic,codeabs, &
+ initialfield,ATTENUATION_VISCOELASTIC_SOLID,ATTENUATION_PORO_FLUID_PART,deltat, &
+ deltatover2,deltatsquareover2,ibool,kmato,numabs,poroelastic,codeabs, &
accels_poroelastic,velocs_poroelastic,velocw_poroelastic,displs_poroelastic,displw_poroelastic,&
b_accels_poroelastic,b_displs_poroelastic,b_displw_poroelastic,&
density,porosity,tortuosity,permeability,poroelastcoef,xix,xiz,gammax,gammaz, &
jacobian,source_time_function,sourcearray,adj_sourcearrays,e11, &
- e13,e11_veloc,e13_veloc,e11_accel,e13_accel,dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n, &
- dvx_dxl_n,dvz_dzl_n,dvz_dxl_n,dvx_dzl_n,hprime_xx,hprimewgll_xx, &
+ e13,e11_veloc,e13_veloc,e11_accel,e13_accel,dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n, &
+ dvx_dxl_n,dvz_dzl_n,dvz_dxl_n,dvx_dzl_n,hprime_xx,hprimewgll_xx, &
hprime_zz,hprimewgll_zz,wxgll,wzgll,inv_tau_sigma_nu2,&
phi_nu2,Mu_nu2,N_SLS, &
rx_viscous,rz_viscous,theta_e,theta_s,&
@@ -115,8 +115,8 @@
integer :: N_SLS
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e11,e13
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e11_veloc,e13_veloc
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e11_accel,e13_accel
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e11_veloc,e13_veloc
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e11_accel,e13_accel
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e11_LDDRK,e13_LDDRK
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e11_initial_rk,e13_initial_rk
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS,stage_time_scheme) :: e11_force_RK, e13_force_RK
@@ -125,8 +125,8 @@
real(kind=CUSTOM_REAL) :: e11_sum,e13_sum
integer :: i_sls
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: &
- dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n,dvx_dxl_n,dvz_dzl_n,dvz_dxl_n,dvx_dzl_n
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: &
+ dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n,dvx_dxl_n,dvz_dzl_n,dvz_dxl_n,dvx_dzl_n
! viscous attenuation (poroelastic media)
double precision, dimension(NGLLX,NGLLZ,nspec) :: rx_viscous
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90 2013-03-11 16:22:35 UTC (rev 21486)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90 2013-03-11 17:45:03 UTC (rev 21487)
@@ -45,14 +45,14 @@
subroutine compute_forces_viscoelastic(p_sv,nglob,nspec,myrank,nelemabs,numat, &
ispec_selected_source,ispec_selected_rec,is_proc_source,which_proc_receiver, &
source_type,it,NSTEP,anyabs,assign_external_model, &
- initialfield,ATTENUATION_VISCOELASTIC_SOLID,anglesource, &
- deltatover2,deltatsquareover2,ibool,kmato,numabs,elastic,codeabs, &
+ initialfield,ATTENUATION_VISCOELASTIC_SOLID,anglesource, &
+ deltatover2,deltatsquareover2,ibool,kmato,numabs,elastic,codeabs, &
accel_elastic,veloc_elastic,displ_elastic, &
density,poroelastcoef,xix,xiz,gammax,gammaz, &
jacobian,vpext,vsext,rhoext,c11ext,c13ext,c15ext,c33ext,c35ext,c55ext,anisotropic,anisotropy, &
source_time_function,sourcearray,adj_sourcearrays,e1,e11, &
- e13,e1_veloc,e11_veloc,e13_veloc,e1_accel,e11_accel,e13_accel,dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n, &
- dvx_dxl_n,dvz_dzl_n,dvz_dxl_n,dvx_dzl_n,hprime_xx,hprimewgll_xx, &
+ e13,e1_veloc,e11_veloc,e13_veloc,e1_accel,e11_accel,e13_accel,dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n, &
+ dvx_dxl_n,dvz_dzl_n,dvz_dxl_n,dvx_dzl_n,hprime_xx,hprimewgll_xx, &
hprime_zz,hprimewgll_zz,wxgll,wzgll,inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2,Mu_nu1,Mu_nu2,N_SLS, &
deltat,coord,add_Bielak_conditions, &
x0_source, z0_source, A_plane, B_plane, C_plane, anglesource_refl, c_inc, c_refl, time_offset,f0, &
@@ -60,16 +60,16 @@
nleft,nright,nbot,over_critical_angle,NSOURCES,nrec,SIMULATION_TYPE,SAVE_FORWARD,b_absorb_elastic_left,&
b_absorb_elastic_right,b_absorb_elastic_bottom,b_absorb_elastic_top,nspec_left,nspec_right,&
nspec_bottom,nspec_top,ib_left,ib_right,ib_bottom,ib_top,&
- e1_LDDRK,e11_LDDRK,e13_LDDRK,alpha_LDDRK,beta_LDDRK, &
+ e1_LDDRK,e11_LDDRK,e13_LDDRK,alpha_LDDRK,beta_LDDRK, &
e1_initial_rk,e11_initial_rk,e13_initial_rk,e1_force_RK, e11_force_RK, e13_force_RK, &
stage_time_scheme,i_stage,ADD_SPRING_TO_STACEY,x_center_spring,z_center_spring,nadj_rec_local, &
is_PML,nspec_PML,spec_to_PML,region_CPML, &
K_x_store,K_z_store,d_x_store,d_z_store,alpha_x_store,alpha_z_store, &
rmemory_displ_elastic,rmemory_dux_dx,rmemory_dux_dz,rmemory_duz_dx,rmemory_duz_dz, &
rmemory_dux_dx_prime,rmemory_dux_dz_prime,rmemory_duz_dx_prime,rmemory_duz_dz_prime, &
- rmemory_displ_elastic_LDDRK,rmemory_dux_dx_LDDRK,rmemory_dux_dz_LDDRK,&
- rmemory_duz_dx_LDDRK,rmemory_duz_dz_LDDRK, &
- PML_BOUNDARY_CONDITIONS,ROTATE_PML_ACTIVATE,ROTATE_PML_ANGLE,backward_simulation)
+ rmemory_displ_elastic_LDDRK,rmemory_dux_dx_LDDRK,rmemory_dux_dz_LDDRK,&
+ rmemory_duz_dx_LDDRK,rmemory_duz_dz_LDDRK, &
+ PML_BOUNDARY_CONDITIONS,ROTATE_PML_ACTIVATE,ROTATE_PML_ANGLE,backward_simulation)
! compute forces for the elastic elements
@@ -127,8 +127,8 @@
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_veloc,e11_veloc,e13_veloc
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e1_accel,e11_accel,e13_accel
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e1_veloc,e11_veloc,e13_veloc
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e1_accel,e11_accel,e13_accel
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e1_LDDRK,e11_LDDRK,e13_LDDRK
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS) :: e1_initial_rk,e11_initial_rk,e13_initial_rk
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec,N_SLS,stage_time_scheme) :: e1_force_RK, e11_force_RK, e13_force_RK
@@ -164,10 +164,10 @@
! spatial derivatives
real(kind=CUSTOM_REAL) :: dux_dxi,dux_dgamma,duy_dxi,duy_dgamma,duz_dxi,duz_dgamma
real(kind=CUSTOM_REAL) :: dux_dxl,duy_dxl,duz_dxl,dux_dzl,duy_dzl,duz_dzl
- real(kind=CUSTOM_REAL) :: dux_dxl_prime,duz_dxl_prime,dux_dzl_prime,duz_dzl_prime
- real(kind=CUSTOM_REAL) :: theta,ct,st
+ real(kind=CUSTOM_REAL) :: dux_dxl_prime,duz_dxl_prime,dux_dzl_prime,duz_dzl_prime
+ real(kind=CUSTOM_REAL) :: theta,ct,st
real(kind=CUSTOM_REAL) :: sigma_xx,sigma_xy,sigma_xz,sigma_zy,sigma_zz,sigma_zx
- real(kind=CUSTOM_REAL) :: sigma_xx_prime,sigma_xz_prime,sigma_zz_prime,sigma_zx_prime
+ real(kind=CUSTOM_REAL) :: sigma_xx_prime,sigma_xz_prime,sigma_zz_prime,sigma_zx_prime
real(kind=CUSTOM_REAL) :: nx,nz,vx,vy,vz,vn,rho_vp,rho_vs,tx,ty,tz,weight,xxi,zxi,xgamma,zgamma,jacobian1D
real(kind=CUSTOM_REAL) :: displx,disply,displz,displn,spring_position,displtx,displty,displtz
@@ -193,7 +193,7 @@
lambdal_relaxed_viscoelastic,mul_relaxed_viscoelastic,lambdalplusmul_relaxed_viscoel
! for attenuation
- real(kind=CUSTOM_REAL) :: phinu1,phinu2,tauinvnu1,tauinvnu2,theta_n_u,theta_n_v
+ real(kind=CUSTOM_REAL) :: phinu1,phinu2,tauinvnu1,tauinvnu2,theta_n_u,theta_n_v
! for anisotropy
double precision :: c11,c15,c13,c33,c35,c55
@@ -218,14 +218,14 @@
integer, dimension(nspec) :: region_CPML
logical, dimension(nspec) :: is_PML
integer, dimension(nspec) :: spec_to_PML
- logical :: PML_BOUNDARY_CONDITIONS,ROTATE_PML_ACTIVATE
- double precision ROTATE_PML_ANGLE
+ logical :: PML_BOUNDARY_CONDITIONS,ROTATE_PML_ACTIVATE
+ double precision ROTATE_PML_ANGLE
real(kind=CUSTOM_REAL), dimension(2,3,NGLLX,NGLLZ,nspec_PML) :: rmemory_displ_elastic
real(kind=CUSTOM_REAL), dimension(2,3,NGLLX,NGLLZ,nspec_PML) :: rmemory_displ_elastic_LDDRK
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec_PML) :: &
rmemory_dux_dx,rmemory_dux_dz,rmemory_duz_dx,rmemory_duz_dz
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec_PML) :: &
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec_PML) :: &
rmemory_dux_dx_prime,rmemory_dux_dz_prime,rmemory_duz_dx_prime,rmemory_duz_dz_prime
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec_PML) :: &
rmemory_dux_dx_LDDRK,rmemory_dux_dz_LDDRK,rmemory_duz_dx_LDDRK,rmemory_duz_dz_LDDRK
@@ -558,7 +558,12 @@
end if
if(ROTATE_PML_ACTIVATE)then
- !! DK DK new from Wang eq (21)
+
+ ! second-order accurate convolution term calculation from equation (21) of
+ ! Shumin Wang, Robert Lee, and Fernando L. Teixeira,
+ ! Anisotropic-Medium PML for Vector FETD With Modified Basis Functions,
+ ! IEEE Transactions on Antennas and Propagation, vol. 54, no. 1, (2006)
+
rmemory_dux_dx(i,j,ispec_PML) = coef0 * rmemory_dux_dx(i,j,ispec_PML) &
+ PML_dux_dxl_new(i,j,ispec_PML) * coef1 + PML_dux_dxl(i,j,ispec_PML) * coef2
@@ -578,9 +583,9 @@
+ PML_duz_dxl_new(i,j,ispec_PML) * coef1 + PML_duz_dxl(i,j,ispec_PML) * coef2
endif
- endif
+ endif
- if(stage_time_scheme == 6) then
+ if(stage_time_scheme == 6) then
rmemory_dux_dx_LDDRK(i,j,ispec_PML) = alpha_LDDRK(i_stage) * rmemory_dux_dx_LDDRK(i,j,ispec_PML) &
+ deltat * (-bb * rmemory_dux_dx(i,j,ispec_PML) + PML_dux_dxl(i,j,ispec_PML))
@@ -592,8 +597,8 @@
rmemory_duz_dx(i,j,ispec_PML) = rmemory_duz_dx(i,j,ispec_PML) + &
beta_LDDRK(i_stage) * rmemory_duz_dx_LDDRK(i,j,ispec_PML)
- end if
-
+ end if
+
if(ROTATE_PML_ACTIVATE)then
dux_dxl = PML_dux_dxl(i,j,ispec_PML) + A8 * rmemory_dux_dx(i,j,ispec_PML)
dux_dzl = PML_dux_dzl(i,j,ispec_PML) + A8 * rmemory_dux_dz(i,j,ispec_PML)
@@ -623,7 +628,6 @@
end if
if(ROTATE_PML_ACTIVATE)then
- !! DK DK new from Wang eq (21)
rmemory_dux_dx_prime(i,j,ispec_PML) = coef0 * rmemory_dux_dx_prime(i,j,ispec_PML) &
+ PML_dux_dxl_new(i,j,ispec_PML) *coef1 + PML_dux_dxl(i,j,ispec_PML) * coef2
@@ -655,7 +659,7 @@
+ deltat * (-bb * rmemory_duz_dz(i,j,ispec_PML) + PML_duz_dzl(i,j,ispec_PML))
rmemory_duz_dz(i,j,ispec_PML) = rmemory_duz_dz(i,j,ispec_PML) + &
beta_LDDRK(i_stage) * rmemory_duz_dz_LDDRK(i,j,ispec_PML)
- end if
+ end if
if(ROTATE_PML_ACTIVATE)then
dux_dxl_prime = PML_dux_dxl(i,j,ispec_PML) + A5 * rmemory_dux_dx_prime(i,j,ispec_PML)
@@ -683,7 +687,7 @@
bb = d_x_store(i,j,ispec_PML) / k_x_store(i,j,ispec_PML) + alpha_x_store(i,j,ispec_PML)
- if(stage_time_scheme == 1) then
+ if(stage_time_scheme == 1) then
coef0 = exp(- bb * deltat)
@@ -696,7 +700,6 @@
end if
if(ROTATE_PML_ACTIVATE)then
- !! DK DK new from Wang eq (21)
rmemory_dux_dx(i,j,ispec_PML) = coef0 * rmemory_dux_dx(i,j,ispec_PML) &
+ PML_dux_dxl_new(i,j,ispec_PML) * coef1 + PML_dux_dxl(i,j,ispec_PML) * coef2
@@ -718,7 +721,7 @@
end if
- if(stage_time_scheme == 6) then
+ if(stage_time_scheme == 6) then
rmemory_dux_dx_LDDRK(i,j,ispec_PML) = alpha_LDDRK(i_stage) * rmemory_dux_dx_LDDRK(i,j,ispec_PML) &
+ deltat * (-bb * rmemory_dux_dx(i,j,ispec_PML) + PML_dux_dxl(i,j,ispec_PML))
rmemory_dux_dx(i,j,ispec_PML) = rmemory_dux_dx(i,j,ispec_PML) + &
@@ -728,7 +731,7 @@
+ deltat * (-bb * rmemory_duz_dx(i,j,ispec_PML) + PML_duz_dxl(i,j,ispec_PML))
rmemory_duz_dx(i,j,ispec_PML) = rmemory_duz_dx(i,j,ispec_PML) + &
beta_LDDRK(i_stage) * rmemory_duz_dx_LDDRK(i,j,ispec_PML)
- end if
+ end if
if(ROTATE_PML_ACTIVATE)then
dux_dxl = PML_dux_dxl(i,j,ispec_PML) + A8 * rmemory_dux_dx(i,j,ispec_PML)
@@ -749,7 +752,7 @@
bb = d_z_store(i,j,ispec_PML) / k_z_store(i,j,ispec_PML) + alpha_z_store(i,j,ispec_PML)
- if(stage_time_scheme == 1) then
+ if(stage_time_scheme == 1) then
coef0 = exp(- bb * deltat)
@@ -762,7 +765,6 @@
end if
if(ROTATE_PML_ACTIVATE)then
- !! DK DK new from Wang eq (21)
rmemory_dux_dx_prime(i,j,ispec_PML) = coef0 * rmemory_dux_dx_prime(i,j,ispec_PML) &
+ PML_dux_dxl_new(i,j,ispec_PML) *coef1 + PML_dux_dxl(i,j,ispec_PML) * coef2
@@ -785,7 +787,7 @@
end if
- if(stage_time_scheme == 6) then
+ if(stage_time_scheme == 6) then
rmemory_dux_dz_LDDRK(i,j,ispec_PML) = alpha_LDDRK(i_stage) * rmemory_dux_dz_LDDRK(i,j,ispec_PML) &
+ deltat * (-bb * rmemory_dux_dz(i,j,ispec_PML) + PML_dux_dzl(i,j,ispec_PML))
rmemory_dux_dz(i,j,ispec_PML) = rmemory_dux_dz(i,j,ispec_PML) + &
@@ -795,7 +797,7 @@
+ deltat * (-bb * rmemory_duz_dz(i,j,ispec_PML) + PML_duz_dzl(i,j,ispec_PML))
rmemory_duz_dz(i,j,ispec_PML) = rmemory_duz_dz(i,j,ispec_PML) + &
beta_LDDRK(i_stage) * rmemory_duz_dz_LDDRK(i,j,ispec_PML)
- end if
+ end if
if(ROTATE_PML_ACTIVATE)then
dux_dxl_prime = PML_dux_dxl(i,j,ispec_PML) + A5 * rmemory_dux_dx_prime(i,j,ispec_PML)
@@ -817,7 +819,7 @@
A7 = d_z_store(i,j,ispec_PML)
bb = alpha_z_store(i,j,ispec_PML)
- if(stage_time_scheme == 1) then
+ if(stage_time_scheme == 1) then
coef0 = exp(- bb * deltat)
if ( abs( bb ) > 0.001_CUSTOM_REAL) then
@@ -829,7 +831,6 @@
end if
if(ROTATE_PML_ACTIVATE)then
- !! DK DK new from Wang eq (21)
rmemory_dux_dx(i,j,ispec_PML) = coef0 * rmemory_dux_dx(i,j,ispec_PML) &
+ PML_dux_dxl_new(i,j,ispec_PML) * coef1 + PML_dux_dxl(i,j,ispec_PML) * coef2
@@ -849,9 +850,9 @@
+ PML_duz_dxl_new(i,j,ispec_PML) * coef1 + PML_duz_dxl(i,j,ispec_PML) * coef2
endif
- end if
+ end if
- if(stage_time_scheme == 6) then
+ if(stage_time_scheme == 6) then
rmemory_dux_dx_LDDRK(i,j,ispec_PML) = alpha_LDDRK(i_stage) * rmemory_dux_dx_LDDRK(i,j,ispec_PML) &
+ deltat * (-bb * rmemory_dux_dx(i,j,ispec_PML) + PML_dux_dxl(i,j,ispec_PML))
rmemory_dux_dx(i,j,ispec_PML) = rmemory_dux_dx(i,j,ispec_PML) + &
@@ -861,7 +862,7 @@
+ deltat * (-bb * rmemory_duz_dx(i,j,ispec_PML) + PML_duz_dxl(i,j,ispec_PML))
rmemory_duz_dx(i,j,ispec_PML) = rmemory_duz_dx(i,j,ispec_PML) + &
beta_LDDRK(i_stage) * rmemory_duz_dx_LDDRK(i,j,ispec_PML)
- end if
+ end if
if(ROTATE_PML_ACTIVATE)then
dux_dxl = PML_dux_dxl(i,j,ispec_PML) + A7 * rmemory_dux_dx(i,j,ispec_PML)
@@ -879,7 +880,7 @@
A6 = - d_z_store(i,j,ispec_PML) / ( k_z_store(i,j,ispec_PML) ** 2 )
bb = d_z_store(i,j,ispec_PML) / k_z_store(i,j,ispec_PML) + alpha_z_store(i,j,ispec_PML)
- if(stage_time_scheme == 1) then
+ if(stage_time_scheme == 1) then
coef0 = exp(-bb * deltat)
if ( abs(bb) > 0.001_CUSTOM_REAL ) then
@@ -891,7 +892,6 @@
end if
if(ROTATE_PML_ACTIVATE)then
- !! DK DK new from Wang eq (21)
rmemory_dux_dx_prime(i,j,ispec_PML) = coef0 * rmemory_dux_dx_prime(i,j,ispec_PML) &
+ PML_dux_dxl_new(i,j,ispec_PML) *coef1 + PML_dux_dxl(i,j,ispec_PML) * coef2
@@ -913,7 +913,7 @@
end if
- if(stage_time_scheme == 6) then
+ if(stage_time_scheme == 6) then
rmemory_dux_dz_LDDRK(i,j,ispec_PML) = alpha_LDDRK(i_stage) * rmemory_dux_dz_LDDRK(i,j,ispec_PML) &
+ deltat * (-bb * rmemory_dux_dz(i,j,ispec_PML) + PML_dux_dzl(i,j,ispec_PML))
rmemory_dux_dz(i,j,ispec_PML) = rmemory_dux_dz(i,j,ispec_PML) + &
@@ -923,7 +923,7 @@
+ deltat * (-bb * rmemory_duz_dz(i,j,ispec_PML) + PML_duz_dzl(i,j,ispec_PML))
rmemory_duz_dz(i,j,ispec_PML) = rmemory_duz_dz(i,j,ispec_PML) + &
beta_LDDRK(i_stage) * rmemory_duz_dz_LDDRK(i,j,ispec_PML)
- end if
+ end if
if(ROTATE_PML_ACTIVATE)then
dux_dxl_prime = PML_dux_dxl(i,j,ispec_PML) + A6 * rmemory_dux_dx_prime(i,j,ispec_PML)
@@ -1132,7 +1132,7 @@
if (region_CPML(ispec) == CPML_LEFT .or. region_CPML(ispec) == CPML_RIGHT) then
bb = alpha_x_store(i,j,ispec_PML)
- if(stage_time_scheme == 1) then
+ if(stage_time_scheme == 1) then
coef0 = exp(- bb * deltat)
if ( abs( bb ) > 0.001_CUSTOM_REAL) then
@@ -1161,7 +1161,7 @@
!------------------------------------------------------------
bb = alpha_x_store(i,j,ispec_PML)
- if(stage_time_scheme == 1) then
+ if(stage_time_scheme == 1) then
coef0 = exp(- bb * deltat)
if ( abs(bb) > 0.001_CUSTOM_REAL ) then
@@ -1183,7 +1183,7 @@
!------------------------------------------------------------
bb = alpha_z_store(i,j,ispec_PML)
- if(stage_time_scheme == 1) then
+ if(stage_time_scheme == 1) then
coef0 = exp(- bb * deltat)
if ( abs(bb) > 0.001_CUSTOM_REAL ) then
@@ -1260,7 +1260,7 @@
beta_LDDRK(i_stage) * rmemory_displ_elastic_LDDRK(1,3,i,j,ispec_PML)
rmemory_displ_elastic(2,3,i,j,ispec_PML) =0._CUSTOM_REAL
- end if
+ end if
accel_elastic_PML(1,i,j,ispec_PML)= wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec) * &
( &
@@ -1295,14 +1295,14 @@
(it+0.5)*deltat*alpha_x_store(i,j,ispec_PML)**2*d_x_store(i,j,ispec_PML)*d_z_store(i,j,ispec_PML)
A4 = -alpha_x_store(i,j,ispec_PML) ** 2*d_x_store(i,j,ispec_PML)*d_z_store(i,j,ispec_PML)
- if(stage_time_scheme == 6) then
+ if(stage_time_scheme == 6) then
A3 = alpha_x_store(i,j,ispec_PML) ** 2*(d_x_store(i,j,ispec_PML) * k_z_store(i,j,ispec_PML)+ &
d_z_store(i,j,ispec_PML) * k_x_store(i,j,ispec_PML)) &
-2._CUSTOM_REAL * alpha_x_store(i,j,ispec_PML)*d_x_store(i,j,ispec_PML)*d_z_store(i,j,ispec_PML)
A4 = alpha_x_store(i,j,ispec_PML) ** 2*d_x_store(i,j,ispec_PML)*d_z_store(i,j,ispec_PML)
- end if
+ end if
- if(stage_time_scheme == 6) then
+ if(stage_time_scheme == 6) then
bb = alpha_z_store(i,j,ispec_PML)
@@ -1358,7 +1358,7 @@
A3 = 0._CUSTOM_REAL
A4 = d_z_store(i,j,ispec_PML) * alpha_z_store(i,j,ispec_PML) ** 2
- if(stage_time_scheme == 6) then
+ if(stage_time_scheme == 6) then
bb = alpha_z_store(i,j,ispec_PML)
@@ -1692,7 +1692,7 @@
endif
endif
endif !end of backward_simulation
- endif !end of elasitic
+ endif !end of elasitic
enddo
@@ -1835,7 +1835,7 @@
endif
endif !end of backward_simulation
- endif !end of elasitic
+ endif !end of elasitic
enddo
@@ -1970,7 +1970,7 @@
endif
endif !end of backward_simulation
- endif !end of elasitic
+ endif !end of elasitic
enddo
@@ -2122,7 +2122,7 @@
!========================================================================
subroutine compute_forces_viscoelastic_pre_kernel(p_sv,nglob,nspec,displ_elastic,b_displ_elastic,&
- mu_k,kappa_k,elastic,ibool,hprime_xx,hprime_zz,xix,xiz,gammax,gammaz,SIMULATION_TYPE)
+ mu_k,kappa_k,elastic,ibool,hprime_xx,hprime_zz,xix,xiz,gammax,gammaz,SIMULATION_TYPE)
! compute forces for the elastic elements
@@ -2150,7 +2150,7 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx,hprime_zz
real(kind=CUSTOM_REAL), dimension(nglob) :: mu_k,kappa_k
-!ZN logical :: backward_simulation
+!ZN logical :: backward_simulation
do ispec = 1,nspec
if(elastic(ispec))then
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90 2013-03-11 16:22:35 UTC (rev 21486)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/pml_init.F90 2013-03-11 17:45:03 UTC (rev 21487)
@@ -166,18 +166,18 @@
do i=1,NGLLX,NGLLX-1
iglob=ibool(i,j,ispec)
do k=1,ncorner
- if (iglob==icorner_iglob(k)) then
- PML_interior_interface(ibound,ispec) = .true.
- end if
- end do
- end do
- end do
- end if
- end do
+ if (iglob==icorner_iglob(k)) then
+ PML_interior_interface(ibound,ispec) = .true.
+ end if
+ end do
+ end do
+ end do
+ end if
+ end do
- end do !end nelem_thickness loop
+ end do !end nelem_thickness loop
- endif !end of SIMULATION_TYPE == 2
+ endif !end of SIMULATION_TYPE == 2
write(IOUT,*) "number of PML spectral elements on side ", ibound,":", nspec_PML
@@ -190,7 +190,7 @@
.and. (.not. PML_interior_interface(IRIGHT,ispec)) &
.and. (.not. PML_interior_interface(ILEFT,ispec)) &
.and. (.not. which_PML_elem(IRIGHT,ispec)) &
- .and. (.not. which_PML_elem(ILEFT,ispec)))then
+ .and. (.not. which_PML_elem(ILEFT,ispec)))then
nglob_interface = nglob_interface + 5
elseif(PML_interior_interface(ITOP,ispec) &
.and. (.not. PML_interior_interface(IRIGHT,ispec)) &
@@ -202,25 +202,25 @@
.and. (.not. PML_interior_interface(IBOTTOM,ispec)) &
.and. (.not. PML_interior_interface(ITOP,ispec)) &
.and. (.not. which_PML_elem(IBOTTOM,ispec)) &
- .and. (.not. which_PML_elem(ITOP,ispec)))then
+ .and. (.not. which_PML_elem(ITOP,ispec)))then
nglob_interface = nglob_interface + 5
elseif(PML_interior_interface(ILEFT,ispec) &
.and. (.not. PML_interior_interface(IBOTTOM,ispec)) &
.and. (.not. PML_interior_interface(ITOP,ispec)) &
.and. (.not. which_PML_elem(IBOTTOM,ispec)) &
.and. (.not. which_PML_elem(ITOP,ispec)))then
- nglob_interface = nglob_interface + 5
+ nglob_interface = nglob_interface + 5
elseif(PML_interior_interface(ILEFT,ispec) &
- .and. PML_interior_interface(IBOTTOM,ispec))then
+ .and. PML_interior_interface(IBOTTOM,ispec))then
nglob_interface = nglob_interface + 10
elseif(PML_interior_interface(IRIGHT,ispec) &
.and. PML_interior_interface(IBOTTOM,ispec))then
nglob_interface = nglob_interface + 10
elseif(PML_interior_interface(ILEFT,ispec) &
- .and. PML_interior_interface(ITOP,ispec))then
+ .and. PML_interior_interface(ITOP,ispec))then
nglob_interface = nglob_interface + 10
elseif(PML_interior_interface(IRIGHT,ispec) &
- .and. PML_interior_interface(ITOP,ispec))then
+ .and. PML_interior_interface(ITOP,ispec))then
nglob_interface = nglob_interface + 10
endif
enddo
@@ -338,7 +338,7 @@
integer :: nglob_interface, nspec
logical, dimension(4,nspec) :: PML_interior_interface
- logical, dimension(4,nspec) :: which_PML_elem
+ logical, dimension(4,nspec) :: which_PML_elem
integer :: ispec
integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
integer, dimension(nglob_interface) :: point_interface
@@ -376,7 +376,7 @@
point_interface(nglob_interface + 2) = ibool(NGLLX,2,ispec)
point_interface(nglob_interface + 3) = ibool(NGLLX,3,ispec)
point_interface(nglob_interface + 4) = ibool(NGLLX,4,ispec)
- point_interface(nglob_interface + 5) = ibool(NGLLX,5,ispec)
+ point_interface(nglob_interface + 5) = ibool(NGLLX,5,ispec)
nglob_interface = nglob_interface + 5
elseif(PML_interior_interface(ILEFT,ispec) &
.and. (.not. PML_interior_interface(IBOTTOM,ispec)) &
@@ -388,7 +388,7 @@
point_interface(nglob_interface + 3) = ibool(1,3,ispec)
point_interface(nglob_interface + 4) = ibool(1,4,ispec)
point_interface(nglob_interface + 5) = ibool(1,5,ispec)
- nglob_interface = nglob_interface + 5
+ nglob_interface = nglob_interface + 5
elseif(PML_interior_interface(ILEFT,ispec) &
.and. PML_interior_interface(IBOTTOM,ispec))then
point_interface(nglob_interface + 1) = ibool(1,1,ispec)
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/read_databases.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/read_databases.f90 2013-03-11 16:22:35 UTC (rev 21486)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/read_databases.f90 2013-03-11 17:45:03 UTC (rev 21487)
@@ -75,9 +75,9 @@
logical :: meshvect,modelvect,boundvect,initialfield,add_Bielak_conditions, &
assign_external_model,READ_EXTERNAL_SEP_FILE, &
output_grid_ASCII,output_energy,output_wavefield_dumps,p_sv,use_binary_for_wavefield_dumps
- logical :: ATTENUATION_VISCOELASTIC_SOLID,ATTENUATION_PORO_FLUID_PART,PML_BOUNDARY_CONDITIONS,ROTATE_PML_ACTIVATE, &
+ logical :: ATTENUATION_VISCOELASTIC_SOLID,ATTENUATION_PORO_FLUID_PART,PML_BOUNDARY_CONDITIONS,ROTATE_PML_ACTIVATE, &
save_ASCII_seismograms,save_binary_seismograms_single,save_binary_seismograms_double,DRAW_SOURCES_AND_RECEIVERS
- double precision :: ROTATE_PML_ANGLE
+ double precision :: ROTATE_PML_ANGLE
double precision :: cutsnaps,sizemax_arrows,anglerec
double precision :: Q0,freq0
double precision :: deltat
@@ -187,11 +187,11 @@
read(IIN,"(a80)") datlin
read(IIN,*) PML_BOUNDARY_CONDITIONS
- read(IIN,"(a80)") datlin
- read(IIN,*) ROTATE_PML_ACTIVATE
+ read(IIN,"(a80)") datlin
+ read(IIN,*) ROTATE_PML_ACTIVATE
- read(IIN,"(a80)") datlin
- read(IIN,*) ROTATE_PML_ANGLE
+ read(IIN,"(a80)") datlin
+ read(IIN,*) ROTATE_PML_ANGLE
read(IIN,"(a80)") datlin
read(IIN,*) read_external_mesh
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2013-03-11 16:22:35 UTC (rev 21486)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2013-03-11 17:45:03 UTC (rev 21487)
@@ -569,9 +569,9 @@
logical interpol,meshvect,modelvect,boundvect,assign_external_model,initialfield, &
output_grid_ASCII,output_grid_Gnuplot,ATTENUATION_VISCOELASTIC_SOLID,output_postscript_snapshot,output_color_image, &
plot_lowerleft_corner_only,add_Bielak_conditions,output_energy,READ_EXTERNAL_SEP_FILE, &
- output_wavefield_dumps,use_binary_for_wavefield_dumps,PML_BOUNDARY_CONDITIONS,ROTATE_PML_ACTIVATE
+ output_wavefield_dumps,use_binary_for_wavefield_dumps,PML_BOUNDARY_CONDITIONS,ROTATE_PML_ACTIVATE
- double precision :: ROTATE_PML_ANGLE
+ double precision :: ROTATE_PML_ANGLE
!! DK DK for CPML_element_file
logical :: read_external_mesh
@@ -594,9 +594,9 @@
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: e1,e11,e13
!DK DK e1_veloc,e11_veloc,e13_veloc denote first derivative of e1,e11,e13 respect to t
- real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: e1_veloc,e11_veloc,e13_veloc
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: e1_veloc,e11_veloc,e13_veloc
!DK DK e1_accel,e11_accel,e13_accel denote second derivative of e1,e11,e13 respect to t
- real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: e1_accel,e11_accel,e13_accel
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: e1_accel,e11_accel,e13_accel
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: e1_LDDRK,e11_LDDRK,e13_LDDRK
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: e1_initial_rk,e11_initial_rk,e13_initial_rk
real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: e1_force_rk,e11_force_rk,e13_force_rk
@@ -1019,8 +1019,8 @@
real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: &
rmemory_dux_dx,rmemory_duz_dx,rmemory_dux_dz,rmemory_duz_dz
- real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: &
- rmemory_dux_dx_prime,rmemory_duz_dx_prime,rmemory_dux_dz_prime,rmemory_duz_dz_prime
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: &
+ rmemory_dux_dx_prime,rmemory_duz_dx_prime,rmemory_dux_dz_prime,rmemory_duz_dz_prime
real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: &
rmemory_dux_dx_LDDRK,rmemory_duz_dx_LDDRK,rmemory_dux_dz_LDDRK,rmemory_duz_dz_LDDRK
@@ -1028,12 +1028,12 @@
integer, dimension(:), allocatable :: spec_to_PML,icorner_iglob
logical, dimension(:,:), allocatable :: which_PML_elem
-!ZN defined for using PML in elastic simulation
+!ZN defined for using PML in elastic simulation
logical, dimension(:,:), allocatable :: PML_interior_interface
integer, dimension(:), allocatable :: point_interface
- real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: pml_interfeace_history_displ
- real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: pml_interfeace_history_veloc
- real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: pml_interfeace_history_accel
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: pml_interfeace_history_displ
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: pml_interfeace_history_veloc
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: pml_interfeace_history_accel
integer :: nglob_interface
real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: rmemory_displ_elastic
@@ -1261,9 +1261,9 @@
allocate(antecedent_list(nspec))
allocate(perm(nspec))
endif
-
-! DK DK add support for using pml in mpi mode with external mesh
- allocate(region_CPML(nspec))
+
+! DK DK add support for using pml in mpi mode with external mesh
+ allocate(region_CPML(nspec))
call read_databases_mato(ipass,nspec,ngnod,kmato,knods, &
perm,antecedent_list,region_CPML)
@@ -1384,25 +1384,25 @@
allocate(e11(NGLLX,NGLLZ,nspec_allocate,N_SLS))
allocate(e13(NGLLX,NGLLZ,nspec_allocate,N_SLS))
- allocate(e1_veloc(NGLLX,NGLLZ,nspec_allocate,N_SLS))
- allocate(e11_veloc(NGLLX,NGLLZ,nspec_allocate,N_SLS))
- allocate(e13_veloc(NGLLX,NGLLZ,nspec_allocate,N_SLS))
+ allocate(e1_veloc(NGLLX,NGLLZ,nspec_allocate,N_SLS))
+ allocate(e11_veloc(NGLLX,NGLLZ,nspec_allocate,N_SLS))
+ allocate(e13_veloc(NGLLX,NGLLZ,nspec_allocate,N_SLS))
- allocate(e1_accel(NGLLX,NGLLZ,nspec_allocate,N_SLS))
- allocate(e11_accel(NGLLX,NGLLZ,nspec_allocate,N_SLS))
- allocate(e13_accel(NGLLX,NGLLZ,nspec_allocate,N_SLS))
+ allocate(e1_accel(NGLLX,NGLLZ,nspec_allocate,N_SLS))
+ allocate(e11_accel(NGLLX,NGLLZ,nspec_allocate,N_SLS))
+ allocate(e13_accel(NGLLX,NGLLZ,nspec_allocate,N_SLS))
e1(:,:,:,:) = 0._CUSTOM_REAL
e11(:,:,:,:) = 0._CUSTOM_REAL
e13(:,:,:,:) = 0._CUSTOM_REAL
- e1_veloc(:,:,:,:) = 0._CUSTOM_REAL
- e11_veloc(:,:,:,:) = 0._CUSTOM_REAL
- e13_veloc(:,:,:,:) = 0._CUSTOM_REAL
+ e1_veloc(:,:,:,:) = 0._CUSTOM_REAL
+ e11_veloc(:,:,:,:) = 0._CUSTOM_REAL
+ e13_veloc(:,:,:,:) = 0._CUSTOM_REAL
- e1_accel(:,:,:,:) = 0._CUSTOM_REAL
- e11_accel(:,:,:,:) = 0._CUSTOM_REAL
- e13_accel(:,:,:,:) = 0._CUSTOM_REAL
+ e1_accel(:,:,:,:) = 0._CUSTOM_REAL
+ e11_accel(:,:,:,:) = 0._CUSTOM_REAL
+ e13_accel(:,:,:,:) = 0._CUSTOM_REAL
if(time_stepping_scheme == 2)then
allocate(e1_LDDRK(NGLLX,NGLLZ,nspec_allocate,N_SLS))
@@ -1443,10 +1443,10 @@
allocate(duz_dzl_n(NGLLX,NGLLZ,nspec_allocate))
allocate(duz_dxl_n(NGLLX,NGLLZ,nspec_allocate))
allocate(dux_dzl_n(NGLLX,NGLLZ,nspec_allocate))
- allocate(dvx_dxl_n(NGLLX,NGLLZ,nspec_allocate))
- allocate(dvz_dzl_n(NGLLX,NGLLZ,nspec_allocate))
- allocate(dvz_dxl_n(NGLLX,NGLLZ,nspec_allocate))
- allocate(dvx_dzl_n(NGLLX,NGLLZ,nspec_allocate))
+ allocate(dvx_dxl_n(NGLLX,NGLLZ,nspec_allocate))
+ allocate(dvz_dzl_n(NGLLX,NGLLZ,nspec_allocate))
+ allocate(dvz_dxl_n(NGLLX,NGLLZ,nspec_allocate))
+ allocate(dvx_dzl_n(NGLLX,NGLLZ,nspec_allocate))
allocate(Mu_nu1(NGLLX,NGLLZ,nspec))
allocate(Mu_nu2(NGLLX,NGLLZ,nspec))
endif
@@ -2926,14 +2926,14 @@
if( PML_BOUNDARY_CONDITIONS .and. anyabs_glob ) then
!PML code
- allocate(is_PML(nspec))
+ allocate(is_PML(nspec))
allocate(icorner_iglob(nglob))
allocate(which_PML_elem(4,nspec))
allocate(spec_to_PML(nspec))
if(SIMULATION_TYPE == 2 .or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))then
allocate(PML_interior_interface(4,nspec))
- PML_interior_interface = .false.
+ PML_interior_interface = .false.
else
allocate(PML_interior_interface(4,1))
endif
@@ -2941,7 +2941,7 @@
is_PML(:) = .false.
which_PML_elem(:,:) = .false.
-! DK DK add support for using pml in mpi mode with external mesh
+! DK DK add support for using pml in mpi mode with external mesh
call pml_init(nspec,nglob,anyabs,ibool,nelemabs,codeabs,numabs,&
nspec_PML,is_PML,which_PML_elem,spec_to_PML, &
icorner_iglob,NELEM_PML_THICKNESS,&
@@ -2958,7 +2958,7 @@
which_PML_elem,point_interface)
deallocate(PML_interior_interface)
write(outputname,'(a,i6.6,a)') 'pml_elastic',myrank,'.bin'
- open(unit=71,file='OUTPUT_FILES/'//outputname,status='unknown',form='unformatted')
+ open(unit=71,file='OUTPUT_FILES/'//outputname,status='unknown',form='unformatted')
else
allocate(point_interface(1))
allocate(pml_interfeace_history_displ(3,1,1))
@@ -2968,7 +2968,7 @@
if(SIMULATION_TYPE == 2 .and. PML_BOUNDARY_CONDITIONS)then
do it = 1,NSTEP
- do i = 1, nglob_interface
+ do i = 1, nglob_interface
read(71)pml_interfeace_history_accel(1,i,it),&
pml_interfeace_history_accel(2,i,it),&
pml_interfeace_history_accel(3,i,it),&
@@ -2977,9 +2977,9 @@
pml_interfeace_history_veloc(3,i,it),&
pml_interfeace_history_displ(1,i,it),&
pml_interfeace_history_displ(2,i,it),&
- pml_interfeace_history_displ(3,i,it)
- enddo
- enddo
+ pml_interfeace_history_displ(3,i,it)
+ enddo
+ enddo
endif
deallocate(which_PML_elem)
@@ -3027,7 +3027,7 @@
allocate(rmemory_duz_dz(NGLLX,NGLLZ,nspec_PML),stat=ier)
if(ier /= 0) stop 'error: not enough memory to allocate array rmemory_duz_dz'
- if(ROTATE_PML_ACTIVATE)then
+ if(ROTATE_PML_ACTIVATE)then
allocate(rmemory_dux_dx_prime(NGLLX,NGLLZ,nspec_PML),stat=ier)
if(ier /= 0) stop 'error: not enough memory to allocate array rmemory_dux_dx_prime'
allocate(rmemory_dux_dz_prime(NGLLX,NGLLZ,nspec_PML),stat=ier)
@@ -3039,8 +3039,8 @@
endif
if(time_stepping_scheme == 2)then
- allocate(rmemory_displ_elastic_LDDRK(2,3,NGLLX,NGLLZ,nspec_PML),stat=ier)
- if(ier /= 0) stop 'error: not enough memory to allocate array rmemory_displ_elastic'
+ allocate(rmemory_displ_elastic_LDDRK(2,3,NGLLX,NGLLZ,nspec_PML),stat=ier)
+ if(ier /= 0) stop 'error: not enough memory to allocate array rmemory_displ_elastic'
allocate(rmemory_dux_dx_LDDRK(NGLLX,NGLLZ,nspec_PML),stat=ier)
if(ier /= 0) stop 'error: not enough memory to allocate array rmemory_dux_dx'
allocate(rmemory_dux_dz_LDDRK(NGLLX,NGLLZ,nspec_PML),stat=ier)
@@ -3063,7 +3063,7 @@
rmemory_duz_dx(:,:,:) = ZERO
rmemory_duz_dz(:,:,:) = ZERO
- if(ROTATE_PML_ACTIVATE)then
+ if(ROTATE_PML_ACTIVATE)then
rmemory_dux_dx_prime(:,:,:) = ZERO
rmemory_dux_dz_prime(:,:,:) = ZERO
rmemory_duz_dx_prime(:,:,:) = ZERO
@@ -3086,12 +3086,12 @@
allocate(rmemory_duz_dx(1,1,1))
allocate(rmemory_duz_dz(1,1,1))
- allocate(rmemory_dux_dx_prime(1,1,1))
+ allocate(rmemory_dux_dx_prime(1,1,1))
allocate(rmemory_dux_dz_prime(1,1,1))
allocate(rmemory_duz_dx_prime(1,1,1))
allocate(rmemory_duz_dz_prime(1,1,1))
-
+
if(time_stepping_scheme == 2)then
allocate(rmemory_displ_elastic_LDDRK(1,1,1,1,1))
allocate(rmemory_dux_dx_LDDRK(1,1,1))
@@ -3147,7 +3147,7 @@
allocate(rmemory_duz_dx(1,1,1))
allocate(rmemory_duz_dz(1,1,1))
- allocate(rmemory_dux_dx_prime(1,1,1))
+ allocate(rmemory_dux_dx_prime(1,1,1))
allocate(rmemory_dux_dz_prime(1,1,1))
allocate(rmemory_duz_dx_prime(1,1,1))
allocate(rmemory_duz_dz_prime(1,1,1))
@@ -3636,7 +3636,7 @@
!
!----- Files where viscous damping are saved during forward wavefield calculation
!
- if(any_poroelastic .and. (SAVE_FORWARD .or. SIMULATION_TYPE .eq. 2)) then
+ if(any_poroelastic .and. (SAVE_FORWARD .or. SIMULATION_TYPE == 2)) then
allocate(b_viscodampx(nglob))
allocate(b_viscodampz(nglob))
write(outputname,'(a,i6.6,a)') 'viscodampingx',myrank,'.bin'
@@ -5067,7 +5067,7 @@
viscox(:,:,ispec) = viscox_loc(:,:)
viscoz(:,:,ispec) = viscoz_loc(:,:)
endif
-
+
endif ! end of poroelastic element loop
enddo ! end of spectral element loop
@@ -5624,8 +5624,8 @@
call compute_forces_viscoelastic(p_sv,nglob,nspec,myrank,nelemabs,numat, &
ispec_selected_source,ispec_selected_rec,is_proc_source,which_proc_receiver, &
source_type,it,NSTEP,anyabs,assign_external_model, &
- initialfield,ATTENUATION_VISCOELASTIC_SOLID,anglesource, &
- deltatover2,deltatsquareover2,ibool,kmato,numabs,elastic,codeabs, &
+ initialfield,ATTENUATION_VISCOELASTIC_SOLID,anglesource, &
+ deltatover2,deltatsquareover2,ibool,kmato,numabs,elastic,codeabs, &
accel_elastic,veloc_elastic,displ_elastic, &
density,poroelastcoef,xix,xiz,gammax,gammaz, &
jacobian,vpext,vsext,rhoext,c11ext,c13ext,c15ext,c33ext,c35ext,c55ext,anisotropic,anisotropy, &
@@ -5642,7 +5642,7 @@
NSOURCES,nrec,SIMULATION_TYPE,SAVE_FORWARD, &
b_absorb_elastic_left,b_absorb_elastic_right,b_absorb_elastic_bottom,b_absorb_elastic_top, &
nspec_left,nspec_right,nspec_bottom,nspec_top,ib_left,ib_right,ib_bottom,ib_top, &
- e1_LDDRK,e11_LDDRK,e13_LDDRK,alpha_LDDRK,beta_LDDRK, &
+ e1_LDDRK,e11_LDDRK,e13_LDDRK,alpha_LDDRK,beta_LDDRK, &
e1_initial_rk,e11_initial_rk,e13_initial_rk,e1_force_rk, e11_force_rk, e13_force_rk, &
stage_time_scheme,i_stage,ADD_SPRING_TO_STACEY,x_center_spring,z_center_spring,max(1,nadj_rec_local), &
is_PML,nspec_PML,spec_to_PML,region_CPML, &
@@ -5651,11 +5651,11 @@
rmemory_dux_dx_prime,rmemory_dux_dz_prime,rmemory_duz_dx_prime,rmemory_duz_dz_prime, &
rmemory_displ_elastic_LDDRK,rmemory_dux_dx_LDDRK,rmemory_dux_dz_LDDRK,&
rmemory_duz_dx_LDDRK,rmemory_duz_dz_LDDRK, &
- PML_BOUNDARY_CONDITIONS,ROTATE_PML_ACTIVATE,ROTATE_PML_ANGLE,.false.)
+ PML_BOUNDARY_CONDITIONS,ROTATE_PML_ACTIVATE,ROTATE_PML_ANGLE,.false.)
if(SIMULATION_TYPE == 2)then
if(PML_BOUNDARY_CONDITIONS)then
- do ispec = 1,nspec
+ do ispec = 1,nspec
do i = 1, NGLLX
do j = 1, NGLLZ
if(elastic(ispec) .and. is_pml(ispec))then
@@ -5665,11 +5665,11 @@
endif
enddo
enddo
- enddo
+ enddo
endif
if(PML_BOUNDARY_CONDITIONS)then
- do i = 1, nglob_interface
+ do i = 1, nglob_interface
b_accel_elastic(1,point_interface(i)) = pml_interfeace_history_accel(1,i,NSTEP-it+1)
b_accel_elastic(2,point_interface(i)) = pml_interfeace_history_accel(2,i,NSTEP-it+1)
b_accel_elastic(3,point_interface(i)) = pml_interfeace_history_accel(3,i,NSTEP-it+1)
@@ -5678,15 +5678,15 @@
b_veloc_elastic(3,point_interface(i)) = pml_interfeace_history_veloc(3,i,NSTEP-it+1)
b_displ_elastic(1,point_interface(i)) = pml_interfeace_history_displ(1,i,NSTEP-it+1)
b_displ_elastic(2,point_interface(i)) = pml_interfeace_history_displ(2,i,NSTEP-it+1)
- b_displ_elastic(3,point_interface(i)) = pml_interfeace_history_displ(3,i,NSTEP-it+1)
+ b_displ_elastic(3,point_interface(i)) = pml_interfeace_history_displ(3,i,NSTEP-it+1)
enddo
endif
call compute_forces_viscoelastic(p_sv,nglob,nspec,myrank,nelemabs,numat, &
ispec_selected_source,ispec_selected_rec,is_proc_source,which_proc_receiver, &
source_type,it,NSTEP,anyabs,assign_external_model, &
- initialfield,ATTENUATION_VISCOELASTIC_SOLID,anglesource, &
- deltatover2,deltatsquareover2,ibool,kmato,numabs,elastic,codeabs, &
+ initialfield,ATTENUATION_VISCOELASTIC_SOLID,anglesource, &
+ deltatover2,deltatsquareover2,ibool,kmato,numabs,elastic,codeabs, &
b_accel_elastic,b_veloc_elastic,b_displ_elastic, &
density,poroelastcoef,xix,xiz,gammax,gammaz, &
jacobian,vpext,vsext,rhoext,c11ext,c13ext,c15ext,c33ext,c35ext,c55ext,anisotropic,anisotropy, &
@@ -5703,7 +5703,7 @@
NSOURCES,nrec,SIMULATION_TYPE,SAVE_FORWARD, &
b_absorb_elastic_left,b_absorb_elastic_right,b_absorb_elastic_bottom,b_absorb_elastic_top, &
nspec_left,nspec_right,nspec_bottom,nspec_top,ib_left,ib_right,ib_bottom,ib_top, &
- e1_LDDRK,e11_LDDRK,e13_LDDRK,alpha_LDDRK,beta_LDDRK, &
+ e1_LDDRK,e11_LDDRK,e13_LDDRK,alpha_LDDRK,beta_LDDRK, &
e1_initial_rk,e11_initial_rk,e13_initial_rk,e1_force_rk, e11_force_rk, e13_force_rk, &
stage_time_scheme,i_stage,ADD_SPRING_TO_STACEY,x_center_spring,z_center_spring,max(1,nadj_rec_local), &
is_PML,nspec_PML,spec_to_PML,region_CPML, &
@@ -5715,7 +5715,7 @@
!ZN PML_BOUNDARY_CONDITIONS,ROTATE_PML_ACTIVATE,ROTATE_PML_ANGLE,.true.)
.false.,ROTATE_PML_ACTIVATE,ROTATE_PML_ANGLE,.true.)
if(PML_BOUNDARY_CONDITIONS)then
- do ispec = 1,nspec
+ do ispec = 1,nspec
do i = 1, NGLLX
do j = 1, NGLLZ
if(elastic(ispec) .and. is_pml(ispec))then
@@ -5725,10 +5725,10 @@
endif
enddo
enddo
- enddo
+ enddo
endif
- if(PML_BOUNDARY_CONDITIONS)then
+ if(PML_BOUNDARY_CONDITIONS)then
do i = 1, nglob_interface
b_accel_elastic(1,point_interface(i)) = pml_interfeace_history_accel(1,i,NSTEP-it+1)
b_accel_elastic(2,point_interface(i)) = pml_interfeace_history_accel(2,i,NSTEP-it+1)
@@ -5738,7 +5738,7 @@
b_veloc_elastic(3,point_interface(i)) = pml_interfeace_history_veloc(3,i,NSTEP-it+1)
b_displ_elastic(1,point_interface(i)) = pml_interfeace_history_displ(1,i,NSTEP-it+1)
b_displ_elastic(2,point_interface(i)) = pml_interfeace_history_displ(2,i,NSTEP-it+1)
- b_displ_elastic(3,point_interface(i)) = pml_interfeace_history_displ(3,i,NSTEP-it+1)
+ b_displ_elastic(3,point_interface(i)) = pml_interfeace_history_displ(3,i,NSTEP-it+1)
enddo
endif
@@ -6484,14 +6484,14 @@
call compute_forces_poro_solid(nglob,nspec,myrank,nelemabs,numat, &
ispec_selected_source,ispec_selected_rec,is_proc_source,which_proc_receiver,&
source_type,it,NSTEP,anyabs, &
- initialfield,ATTENUATION_VISCOELASTIC_SOLID,ATTENUATION_PORO_FLUID_PART,deltat, &
- deltatover2,deltatsquareover2,ibool,kmato,numabs,poroelastic,codeabs, &
+ initialfield,ATTENUATION_VISCOELASTIC_SOLID,ATTENUATION_PORO_FLUID_PART,deltat, &
+ deltatover2,deltatsquareover2,ibool,kmato,numabs,poroelastic,codeabs, &
accels_poroelastic,velocs_poroelastic,velocw_poroelastic,displs_poroelastic,displw_poroelastic,&
b_accels_poroelastic,b_displs_poroelastic,b_displw_poroelastic,&
density,porosity,tortuosity,permeability,poroelastcoef,xix,xiz,gammax,gammaz, &
jacobian,source_time_function,sourcearray,adj_sourcearrays,e11, &
- e13,e11_veloc,e13_veloc,e11_accel,e13_accel,dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n, &
- dvx_dxl_n,dvz_dzl_n,dvz_dxl_n,dvx_dzl_n,hprime_xx,hprimewgll_xx, &
+ e13,e11_veloc,e13_veloc,e11_accel,e13_accel,dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n, &
+ dvx_dxl_n,dvz_dzl_n,dvz_dxl_n,dvx_dzl_n,hprime_xx,hprimewgll_xx, &
hprime_zz,hprimewgll_zz,wxgll,wzgll,inv_tau_sigma_nu2,&
phi_nu2,Mu_nu2,N_SLS, &
rx_viscous,rz_viscous,theta_e,theta_s,&
@@ -6508,14 +6508,14 @@
call compute_forces_poro_fluid(nglob,nspec,myrank,nelemabs,numat, &
ispec_selected_source,ispec_selected_rec,is_proc_source,which_proc_receiver,&
source_type,it,NSTEP,anyabs, &
- initialfield,ATTENUATION_VISCOELASTIC_SOLID,ATTENUATION_PORO_FLUID_PART,deltat, &
- deltatover2,deltatsquareover2,ibool,kmato,numabs,poroelastic,codeabs, &
+ initialfield,ATTENUATION_VISCOELASTIC_SOLID,ATTENUATION_PORO_FLUID_PART,deltat, &
+ deltatover2,deltatsquareover2,ibool,kmato,numabs,poroelastic,codeabs, &
accelw_poroelastic,velocw_poroelastic,displw_poroelastic,velocs_poroelastic,displs_poroelastic,&
b_accelw_poroelastic,b_displw_poroelastic,b_displs_poroelastic,&
density,porosity,tortuosity,permeability,poroelastcoef,xix,xiz,gammax,gammaz, &
jacobian,source_time_function,sourcearray,adj_sourcearrays,e11, &
- e13,e11_veloc,e13_veloc,e11_accel,e13_accel,dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n, &
- dvx_dxl_n,dvz_dzl_n,dvz_dxl_n,dvx_dzl_n,hprime_xx,hprimewgll_xx, &
+ e13,e11_veloc,e13_veloc,e11_accel,e13_accel,dux_dxl_n,duz_dzl_n,duz_dxl_n,dux_dzl_n, &
+ dvx_dxl_n,dvz_dzl_n,dvz_dxl_n,dvx_dzl_n,hprime_xx,hprimewgll_xx, &
hprime_zz,hprimewgll_zz,wxgll,wzgll,inv_tau_sigma_nu2,&
phi_nu2,Mu_nu2,N_SLS, &
rx_viscous,rz_viscous,theta_e,theta_s,&
@@ -7741,7 +7741,7 @@
call MPI_REDUCE(kinetic_energy, kinetic_energy_total, 1, CUSTOM_MPI_TYPE, MPI_SUM, 0, MPI_COMM_WORLD, ier)
call MPI_REDUCE(potential_energy, potential_energy_total, 1, CUSTOM_MPI_TYPE, MPI_SUM, 0, MPI_COMM_WORLD, ier)
#else
- kinetic_energy_total = kinetic_energy
+ kinetic_energy_total = kinetic_energy
potential_energy_total = potential_energy
#endif
More information about the CIG-COMMITS
mailing list