[cig-commits] r22310 - in seismo/3D/SPECFEM3D_GLOBE/branches/undo_att: . src/specfem3D
xie.zhinan at geodynamics.org
xie.zhinan at geodynamics.org
Sat Jun 15 08:41:35 PDT 2013
Author: xie.zhinan
Date: 2013-06-15 08:41:35 -0700 (Sat, 15 Jun 2013)
New Revision: 22310
Modified:
seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/part1_classical.f90
seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/part1_undo_att.f90
seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/part2_classical.f90
seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/part2_undo_att.f90
seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/part3_kernel_computation.f90
seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/check_simulation_stability.f90
seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/compute_element.F90
seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/compute_forces_crust_mantle_Dev.F90
seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/compute_forces_crust_mantle_noDev.f90
seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/compute_forces_inner_core_Dev.F90
seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/compute_forces_inner_core_noDev.f90
seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/compute_kernels.f90
seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/compute_seismograms.f90
seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/read_forward_arrays.f90
seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/save_forward_arrays.f90
seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/specfem3D.F90
seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/write_movie_volume.f90
Log:
remove the epsilonedev and eps_trace
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/part1_classical.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/part1_classical.f90 2013-06-15 11:05:55 UTC (rev 22309)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/part1_classical.f90 2013-06-15 15:41:35 UTC (rev 22310)
@@ -64,10 +64,22 @@
! do *NOT* use array syntax for that loop, otherwise you will get a compiler error when MOVIE_VOLUME is off
! because the shape of the arrays will not match (due to some arrays purposely declared with a dummy size of 1)
do ispec = 1,NSPEC_CRUST_MANTLE
+ call compute_element_strain_undo_att_noDev(ispec,nglob_crust_mantle,nspec_crust_mantle,&
+ displ_crust_mantle,hprime_xx,hprime_yy,hprime_zz,ibool_crust_mantle,&
+ xix_crust_mantle,xiy_crust_mantle,xiz_crust_mantle,&
+ etax_crust_mantle,etay_crust_mantle,etaz_crust_mantle,&
+ gammax_crust_mantle,gammay_crust_mantle,gammaz_crust_mantle,&
+ epsilondev_loc_crust_mantle,eps_trace_over_3_loc_crust_mantle)
Iepsilondev_crust_mantle(:,:,:,:,ispec) = Iepsilondev_crust_mantle(:,:,:,:,ispec) &
- + deltat*epsilondev_crust_mantle(:,:,:,:,ispec)
+ + deltat*epsilondev_loc_crust_mantle(:,:,:,:)
Ieps_trace_over_3_crust_mantle(:,:,:,ispec) = Ieps_trace_over_3_crust_mantle(:,:,:,ispec) &
- + deltat*eps_trace_over_3_crust_mantle(:,:,:,ispec)
+ + deltat*eps_trace_over_3_loc_crust_mantle(:,:,:)
+
+!ZN
+!ZN Iepsilondev_crust_mantle(:,:,:,:,ispec) = Iepsilondev_crust_mantle(:,:,:,:,ispec) &
+!ZN + deltat*epsilondev_crust_mantle(:,:,:,:,ispec)
+!ZN Ieps_trace_over_3_crust_mantle(:,:,:,ispec) = Ieps_trace_over_3_crust_mantle(:,:,:,ispec) &
+!ZN + deltat*eps_trace_over_3_crust_mantle(:,:,:,ispec)
enddo
endif
@@ -76,7 +88,7 @@
! and output timestamp file to check that simulation is running fine
if(mod(it,NTSTEP_BETWEEN_OUTPUT_INFO) == 0 .or. it == it_begin+4 .or. it == it_end) then
call check_simulation_stability(it,displ_crust_mantle,displ_inner_core,displ_outer_core, &
- eps_trace_over_3_crust_mantle,epsilondev_crust_mantle, &
+!ZN eps_trace_over_3_crust_mantle,epsilondev_crust_mantle, &
!!!!! DK DK UNDO_ATT SIMULATION_TYPE,OUTPUT_FILES,time_start,DT,t0,NSTEP, &
1,OUTPUT_FILES,time_start,DT,t0,NSTEP, & !!!!! DK DK UNDO_ATT
it_begin,it_end,NUMBER_OF_THIS_RUN,NUMBER_OF_RUNS,myrank)
@@ -558,8 +570,9 @@
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,ispec_is_tiso_crust_mantle, &
- R_memory_crust_mantle,epsilondev_crust_mantle, &
- eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
+!ZN R_memory_crust_mantle,epsilondev_crust_mantle, &
+!ZN eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
+ R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,deltat,veloc_crust_mantle, & !ZN
alphaval,betaval,gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5) )
@@ -599,8 +612,9 @@
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,ispec_is_tiso_crust_mantle, &
- R_memory_crust_mantle,epsilondev_crust_mantle, &
- eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
+!ZN R_memory_crust_mantle,epsilondev_crust_mantle, &
+!ZN eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
+ R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,deltat,veloc_crust_mantle, & !ZN
alphaval,betaval,gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5) )
@@ -727,8 +741,9 @@
kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
c11store_inner_core,c33store_inner_core,c12store_inner_core, &
c13store_inner_core,c44store_inner_core, &
- R_memory_inner_core,epsilondev_inner_core, eps_trace_over_3_inner_core,&
- one_minus_sum_beta_inner_core, &
+!ZN R_memory_inner_core,epsilondev_inner_core, eps_trace_over_3_inner_core,&
+!ZN one_minus_sum_beta_inner_core, &
+ R_memory_inner_core,one_minus_sum_beta_inner_core,deltat,veloc_inner_core,& !ZN
alphaval,betaval,gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
@@ -761,8 +776,9 @@
kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
c11store_inner_core,c33store_inner_core,c12store_inner_core, &
c13store_inner_core,c44store_inner_core, &
- R_memory_inner_core,epsilondev_inner_core, eps_trace_over_3_inner_core,&
- one_minus_sum_beta_inner_core, &
+!ZN R_memory_inner_core,epsilondev_inner_core, eps_trace_over_3_inner_core,&
+!ZN one_minus_sum_beta_inner_core, &
+ R_memory_inner_core,one_minus_sum_beta_inner_core,deltat,veloc_inner_core, &
alphaval,betaval,gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
@@ -1095,8 +1111,9 @@
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,ispec_is_tiso_crust_mantle, &
- R_memory_crust_mantle,epsilondev_crust_mantle, &
- eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
+!ZN R_memory_crust_mantle,epsilondev_crust_mantle, &
+!ZN eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
+ R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,deltat,veloc_crust_mantle, &
alphaval,betaval,gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5) )
@@ -1136,8 +1153,9 @@
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,ispec_is_tiso_crust_mantle, &
- R_memory_crust_mantle,epsilondev_crust_mantle, &
- eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
+!ZN R_memory_crust_mantle,epsilondev_crust_mantle, &
+!ZN eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
+ R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,deltat,veloc_crust_mantle, & !ZN
alphaval,betaval,gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5) )
@@ -1172,8 +1190,9 @@
kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
c11store_inner_core,c33store_inner_core,c12store_inner_core, &
c13store_inner_core,c44store_inner_core, &
- R_memory_inner_core,epsilondev_inner_core, eps_trace_over_3_inner_core,&
- one_minus_sum_beta_inner_core, &
+!ZN R_memory_inner_core,epsilondev_inner_core, eps_trace_over_3_inner_core,&
+!ZN one_minus_sum_beta_inner_core, &
+ R_memory_inner_core,one_minus_sum_beta_inner_core,deltat,veloc_inner_core,&
alphaval,betaval,gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
@@ -1206,8 +1225,9 @@
kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
c11store_inner_core,c33store_inner_core,c12store_inner_core, &
c13store_inner_core,c44store_inner_core, &
- R_memory_inner_core,epsilondev_inner_core, eps_trace_over_3_inner_core,&
- one_minus_sum_beta_inner_core, &
+!ZN R_memory_inner_core,epsilondev_inner_core, eps_trace_over_3_inner_core,&
+!ZN one_minus_sum_beta_inner_core, &
+ R_memory_inner_core,one_minus_sum_beta_inner_core,deltat,veloc_inner_core,&
alphaval,betaval,gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
@@ -1609,7 +1629,7 @@
else if (SIMULATION_TYPE == 2) then
call compute_seismograms_adjoint(NSOURCES,nrec_local,displ_crust_mantle, &
- eps_trace_over_3_crust_mantle,epsilondev_crust_mantle, &
+!ZN eps_trace_over_3_crust_mantle,epsilondev_crust_mantle, &
nu_source,Mxx,Myy,Mzz,Mxy,Mxz,Myz, &
hxir_store,hetar_store,hgammar_store, &
hpxir_store,hpetar_store,hpgammar_store, &
@@ -1959,9 +1979,14 @@
if (MOVIE_VOLUME_TYPE == 1) then ! output strains
call write_movie_volume_strains(myrank,npoints_3dmovie, &
LOCAL_PATH,MOVIE_VOLUME_TYPE,MOVIE_COARSE, &
- it,eps_trace_over_3_crust_mantle,epsilondev_crust_mantle, &
- muvstore_crust_mantle_3dmovie, &
- mask_3dmovie,nu_3dmovie)
+!ZN it,eps_trace_over_3_crust_mantle,epsilondev_crust_mantle, &
+!ZN muvstore_crust_mantle_3dmovie, &
+ it,muvstore_crust_mantle_3dmovie, &
+ mask_3dmovie,nu_3dmovie,& !ZN
+ NSPEC_CRUST_MANTLE,hprime_xx,hprime_yy,hprime_zz,ibool_crust_mantle,&
+ xix_crust_mantle,xiy_crust_mantle,xiz_crust_mantle,&
+ etax_crust_mantle,etay_crust_mantle,etaz_crust_mantle,&
+ gammax_crust_mantle,gammay_crust_mantle,gammaz_crust_mantle) !ZN
else if (MOVIE_VOLUME_TYPE == 2 .or. MOVIE_VOLUME_TYPE == 3) then
! output the Time Integral of Strain, or \mu*TIS
@@ -1972,16 +1997,17 @@
mask_3dmovie,nu_3dmovie)
else if (MOVIE_VOLUME_TYPE == 4) then ! output divergence and curl in whole volume
- call write_movie_volume_divcurl(myrank,it,eps_trace_over_3_crust_mantle,&
- div_displ_outer_core, &
- accel_outer_core,kappavstore_outer_core,rhostore_outer_core,ibool_outer_core, &
- eps_trace_over_3_inner_core, &
- epsilondev_crust_mantle,epsilondev_inner_core, &
- LOCAL_PATH, &
- displ_crust_mantle,displ_inner_core,displ_outer_core, &
- veloc_crust_mantle,veloc_inner_core,veloc_outer_core, &
- accel_crust_mantle,accel_inner_core, &
- ibool_crust_mantle,ibool_inner_core)
+!ZN for undo_att this type of MOVIE is not supported
+!ZN call write_movie_volume_divcurl(myrank,it,eps_trace_over_3_crust_mantle,&
+!ZN div_displ_outer_core, &
+!ZN accel_outer_core,kappavstore_outer_core,rhostore_outer_core,ibool_outer_core, &
+!ZN eps_trace_over_3_inner_core, &
+!ZN epsilondev_crust_mantle,epsilondev_inner_core, &
+!ZN LOCAL_PATH, &
+!ZN displ_crust_mantle,displ_inner_core,displ_outer_core, &
+!ZN veloc_crust_mantle,veloc_inner_core,veloc_outer_core, &
+!ZN accel_crust_mantle,accel_inner_core, &
+!ZN ibool_crust_mantle,ibool_inner_core)
else if (MOVIE_VOLUME_TYPE == 5) then ! output displacement
scalingval = scale_displ
@@ -2005,3 +2031,4 @@
endif
endif ! MOVIE_VOLUME
+
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/part1_undo_att.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/part1_undo_att.f90 2013-06-15 11:05:55 UTC (rev 22309)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/part1_undo_att.f90 2013-06-15 15:41:35 UTC (rev 22310)
@@ -64,10 +64,16 @@
! do *NOT* use array syntax for that loop, otherwise you will get a compiler error when MOVIE_VOLUME is off
! because the shape of the arrays will not match (due to some arrays purposely declared with a dummy size of 1)
do ispec = 1,NSPEC_CRUST_MANTLE
+ call compute_element_strain_undo_att_noDev(ispec,nglob_crust_mantle,nspec_crust_mantle,&
+ displ_crust_mantle,hprime_xx,hprime_yy,hprime_zz,ibool_crust_mantle,&
+ xix_crust_mantle,xiy_crust_mantle,xiz_crust_mantle,&
+ etax_crust_mantle,etay_crust_mantle,etaz_crust_mantle,&
+ gammax_crust_mantle,gammay_crust_mantle,gammaz_crust_mantle,&
+ epsilondev_loc_crust_mantle,eps_trace_over_3_loc_crust_mantle)
Iepsilondev_crust_mantle(:,:,:,:,ispec) = Iepsilondev_crust_mantle(:,:,:,:,ispec) &
- + deltat*epsilondev_crust_mantle(:,:,:,:,ispec)
+ + deltat*epsilondev_loc_crust_mantle(:,:,:,:)
Ieps_trace_over_3_crust_mantle(:,:,:,ispec) = Ieps_trace_over_3_crust_mantle(:,:,:,ispec) &
- + deltat*eps_trace_over_3_crust_mantle(:,:,:,ispec)
+ + deltat*eps_trace_over_3_loc_crust_mantle(:,:,:)
enddo
endif
@@ -76,7 +82,7 @@
! and output timestamp file to check that simulation is running fine
if(mod(it,NTSTEP_BETWEEN_OUTPUT_INFO) == 0 .or. it == it_begin+4 .or. it == it_end) then
call check_simulation_stability(it,displ_crust_mantle,displ_inner_core,displ_outer_core, &
- eps_trace_over_3_crust_mantle,epsilondev_crust_mantle, &
+!ZN eps_trace_over_3_crust_mantle,epsilondev_crust_mantle, &
!!!!! DK DK UNDO_ATT SIMULATION_TYPE,OUTPUT_FILES,time_start,DT,t0,NSTEP, &
1,OUTPUT_FILES,time_start,DT,t0,NSTEP, & !!!!! DK DK UNDO_ATT
it_begin,it_end,NUMBER_OF_THIS_RUN,NUMBER_OF_RUNS,myrank)
@@ -558,8 +564,9 @@
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,ispec_is_tiso_crust_mantle, &
- R_memory_crust_mantle,epsilondev_crust_mantle, &
- eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
+!ZN R_memory_crust_mantle,epsilondev_crust_mantle, &
+!ZN eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
+ R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,deltat,veloc_crust_mantle, & !ZN
alphaval,betaval,gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5) )
@@ -599,8 +606,9 @@
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,ispec_is_tiso_crust_mantle, &
- R_memory_crust_mantle,epsilondev_crust_mantle, &
- eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
+!ZN R_memory_crust_mantle,epsilondev_crust_mantle, &
+!ZN eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
+ R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,deltat,veloc_crust_mantle, & !ZN
alphaval,betaval,gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5) )
@@ -727,8 +735,9 @@
kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
c11store_inner_core,c33store_inner_core,c12store_inner_core, &
c13store_inner_core,c44store_inner_core, &
- R_memory_inner_core,epsilondev_inner_core, eps_trace_over_3_inner_core,&
- one_minus_sum_beta_inner_core, &
+!ZN R_memory_inner_core,epsilondev_inner_core, eps_trace_over_3_inner_core,&
+!ZN one_minus_sum_beta_inner_core, &
+ R_memory_inner_core,one_minus_sum_beta_inner_core,deltat,veloc_inner_core,& !ZN
alphaval,betaval,gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
@@ -761,8 +770,9 @@
kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
c11store_inner_core,c33store_inner_core,c12store_inner_core, &
c13store_inner_core,c44store_inner_core, &
- R_memory_inner_core,epsilondev_inner_core, eps_trace_over_3_inner_core,&
- one_minus_sum_beta_inner_core, &
+!ZN R_memory_inner_core,epsilondev_inner_core, eps_trace_over_3_inner_core,&
+!ZN one_minus_sum_beta_inner_core, &
+ R_memory_inner_core,one_minus_sum_beta_inner_core,deltat,veloc_inner_core, &
alphaval,betaval,gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
@@ -1095,8 +1105,9 @@
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,ispec_is_tiso_crust_mantle, &
- R_memory_crust_mantle,epsilondev_crust_mantle, &
- eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
+!ZN R_memory_crust_mantle,epsilondev_crust_mantle, &
+!ZN eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
+ R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,deltat,veloc_crust_mantle, &
alphaval,betaval,gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5) )
@@ -1136,8 +1147,9 @@
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,ispec_is_tiso_crust_mantle, &
- R_memory_crust_mantle,epsilondev_crust_mantle, &
- eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
+!ZN R_memory_crust_mantle,epsilondev_crust_mantle, &
+!ZN eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
+ R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,deltat,veloc_crust_mantle, & !ZN
alphaval,betaval,gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5) )
@@ -1172,8 +1184,9 @@
kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
c11store_inner_core,c33store_inner_core,c12store_inner_core, &
c13store_inner_core,c44store_inner_core, &
- R_memory_inner_core,epsilondev_inner_core, eps_trace_over_3_inner_core,&
- one_minus_sum_beta_inner_core, &
+!ZN R_memory_inner_core,epsilondev_inner_core, eps_trace_over_3_inner_core,&
+!ZN one_minus_sum_beta_inner_core, &
+ R_memory_inner_core,one_minus_sum_beta_inner_core,deltat,veloc_inner_core,&
alphaval,betaval,gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
@@ -1206,8 +1219,9 @@
kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
c11store_inner_core,c33store_inner_core,c12store_inner_core, &
c13store_inner_core,c44store_inner_core, &
- R_memory_inner_core,epsilondev_inner_core, eps_trace_over_3_inner_core,&
- one_minus_sum_beta_inner_core, &
+!ZN R_memory_inner_core,epsilondev_inner_core, eps_trace_over_3_inner_core,&
+!ZN one_minus_sum_beta_inner_core, &
+ R_memory_inner_core,one_minus_sum_beta_inner_core,deltat,veloc_inner_core,&
alphaval,betaval,gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
@@ -1609,7 +1623,7 @@
else if (SIMULATION_TYPE == 2) then
call compute_seismograms_adjoint(NSOURCES,nrec_local,displ_crust_mantle, &
- eps_trace_over_3_crust_mantle,epsilondev_crust_mantle, &
+!ZN eps_trace_over_3_crust_mantle,epsilondev_crust_mantle, &
nu_source,Mxx,Myy,Mzz,Mxy,Mxz,Myz, &
hxir_store,hetar_store,hgammar_store, &
hpxir_store,hpetar_store,hpgammar_store, &
@@ -1989,9 +2003,14 @@
if (MOVIE_VOLUME_TYPE == 1) then ! output strains
call write_movie_volume_strains(myrank,npoints_3dmovie, &
LOCAL_PATH,MOVIE_VOLUME_TYPE,MOVIE_COARSE, &
- it,eps_trace_over_3_crust_mantle,epsilondev_crust_mantle, &
- muvstore_crust_mantle_3dmovie, &
- mask_3dmovie,nu_3dmovie)
+!ZN it,eps_trace_over_3_crust_mantle,epsilondev_crust_mantle, &
+!ZN muvstore_crust_mantle_3dmovie, &
+ it,muvstore_crust_mantle_3dmovie, &
+ mask_3dmovie,nu_3dmovie,& !ZN
+ hprime_xx,hprime_yy,hprime_zz,ibool_crust_mantle,&
+ xix_crust_mantle,xiy_crust_mantle,xiz_crust_mantle,&
+ etax_crust_mantle,etay_crust_mantle,etaz_crust_mantle,&
+ gammax_crust_mantle,gammay_crust_mantle,gammaz_crust_mantle) !ZN
else if (MOVIE_VOLUME_TYPE == 2 .or. MOVIE_VOLUME_TYPE == 3) then
! output the Time Integral of Strain, or \mu*TIS
@@ -2002,16 +2021,17 @@
mask_3dmovie,nu_3dmovie)
else if (MOVIE_VOLUME_TYPE == 4) then ! output divergence and curl in whole volume
- call write_movie_volume_divcurl(myrank,it,eps_trace_over_3_crust_mantle,&
- div_displ_outer_core, &
- accel_outer_core,kappavstore_outer_core,rhostore_outer_core,ibool_outer_core, &
- eps_trace_over_3_inner_core, &
- epsilondev_crust_mantle,epsilondev_inner_core, &
- LOCAL_PATH, &
- displ_crust_mantle,displ_inner_core,displ_outer_core, &
- veloc_crust_mantle,veloc_inner_core,veloc_outer_core, &
- accel_crust_mantle,accel_inner_core, &
- ibool_crust_mantle,ibool_inner_core)
+!ZN for undo_att this type of MOVIE is not supported
+!ZN call write_movie_volume_divcurl(myrank,it,eps_trace_over_3_crust_mantle,&
+!ZN div_displ_outer_core, &
+!ZN accel_outer_core,kappavstore_outer_core,rhostore_outer_core,ibool_outer_core, &
+!ZN eps_trace_over_3_inner_core, &
+!ZN epsilondev_crust_mantle,epsilondev_inner_core, &
+!ZN LOCAL_PATH, &
+!ZN displ_crust_mantle,displ_inner_core,displ_outer_core, &
+!ZN veloc_crust_mantle,veloc_inner_core,veloc_outer_core, &
+!ZN accel_crust_mantle,accel_inner_core, &
+!ZN ibool_crust_mantle,ibool_inner_core)
else if (MOVIE_VOLUME_TYPE == 5) then ! output displacement
scalingval = scale_displ
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/part2_classical.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/part2_classical.f90 2013-06-15 11:05:55 UTC (rev 22309)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/part2_classical.f90 2013-06-15 15:41:35 UTC (rev 22310)
@@ -39,8 +39,11 @@
! and output timestamp file to check that simulation is running fine
if(mod(it,NTSTEP_BETWEEN_OUTPUT_INFO) == 0 .or. it == it_begin+4 .or. it == it_end) then
if (SIMULATION_TYPE == 3) then
+!ZN call check_simulation_stability(it,b_displ_crust_mantle,b_displ_inner_core,b_displ_outer_core, &
+!ZN eps_trace_over_3_crust_mantle,epsilondev_crust_mantle, &
+!ZN SIMULATION_TYPE,OUTPUT_FILES,time_start,DT,t0,NSTEP, &
+!ZN it_begin,it_end,NUMBER_OF_THIS_RUN,NUMBER_OF_RUNS,myrank)
call check_simulation_stability(it,b_displ_crust_mantle,b_displ_inner_core,b_displ_outer_core, &
- eps_trace_over_3_crust_mantle,epsilondev_crust_mantle, &
SIMULATION_TYPE,OUTPUT_FILES,time_start,DT,t0,NSTEP, &
it_begin,it_end,NUMBER_OF_THIS_RUN,NUMBER_OF_RUNS,myrank)
endif
@@ -60,8 +63,10 @@
! for reconstructing the rotational contributions
if(CUSTOM_REAL == SIZE_REAL) then
time = sngl((dble(NSTEP-it)*DT-t0)*scale_t_inv)
+
else
time = (dble(NSTEP-it)*DT-t0)*scale_t_inv
+
endif
b_iphase = 0 ! do not start any non blocking communications at this stage
@@ -345,8 +350,9 @@
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,ispec_is_tiso_crust_mantle, &
- b_R_memory_crust_mantle,b_epsilondev_crust_mantle, &
- b_eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
+!ZN b_R_memory_crust_mantle,b_epsilondev_crust_mantle, &
+!ZN b_eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
+ b_R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,b_deltat,b_veloc_crust_mantle, &
b_alphaval,b_betaval,b_gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5) )
@@ -386,8 +392,9 @@
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,ispec_is_tiso_crust_mantle, &
- b_R_memory_crust_mantle,b_epsilondev_crust_mantle, &
- b_eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
+!ZN b_R_memory_crust_mantle,b_epsilondev_crust_mantle, &
+!ZN b_eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
+ b_R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,b_deltat,b_displ_crust_mantle, &
b_alphaval,b_betaval,b_gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5) )
@@ -423,8 +430,8 @@
kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
c11store_inner_core,c33store_inner_core,c12store_inner_core, &
c13store_inner_core,c44store_inner_core, &
- b_R_memory_inner_core,b_epsilondev_inner_core, b_eps_trace_over_3_inner_core,&
- one_minus_sum_beta_inner_core, &
+!ZN b_R_memory_inner_core,b_epsilondev_inner_core, b_eps_trace_over_3_inner_core,&
+ b_R_memory_inner_core,one_minus_sum_beta_inner_core,b_deltat,b_veloc_inner_core, &
b_alphaval,b_betaval,b_gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
@@ -457,8 +464,8 @@
kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
c11store_inner_core,c33store_inner_core,c12store_inner_core, &
c13store_inner_core,c44store_inner_core, &
- b_R_memory_inner_core,b_epsilondev_inner_core, b_eps_trace_over_3_inner_core,&
- one_minus_sum_beta_inner_core, &
+!ZN b_R_memory_inner_core,b_epsilondev_inner_core, b_eps_trace_over_3_inner_core,&
+ b_R_memory_inner_core,one_minus_sum_beta_inner_core,b_deltat,b_veloc_inner_core, &
b_alphaval,b_betaval,b_gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
@@ -469,6 +476,7 @@
! Stacey
if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) then
if(SIMULATION_TYPE == 3) then
+
call compute_stacey_crust_mantle_backward(ichunk, &
NSTEP,it,ibool_crust_mantle, &
b_accel_crust_mantle, &
@@ -509,12 +517,13 @@
endif ! Stacey conditions
! add adjoint sources and add sources for backward/reconstructed wavefield
- if (SIMULATION_TYPE == 3) &
+ if (SIMULATION_TYPE == 3) then
call compute_add_sources_backward(myrank,NSOURCES,NSTEP, &
b_accel_crust_mantle,sourcearrays, &
DT,t0,tshift_cmt,hdur_gaussian,ibool_crust_mantle, &
islice_selected_source,ispec_selected_source,it, &
hdur,xi_source,eta_source,gamma_source,nu_source)
+ endif
! NOISE_TOMOGRAPHY
! if ( NOISE_TOMOGRAPHY == 1 ) then
@@ -668,8 +677,9 @@
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,ispec_is_tiso_crust_mantle, &
- b_R_memory_crust_mantle,b_epsilondev_crust_mantle, &
- b_eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
+!ZN b_R_memory_crust_mantle,b_epsilondev_crust_mantle, &
+!ZN b_eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
+ b_R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,b_deltat,b_veloc_crust_mantle, &
b_alphaval,b_betaval,b_gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5) )
@@ -709,8 +719,9 @@
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,ispec_is_tiso_crust_mantle, &
- b_R_memory_crust_mantle,b_epsilondev_crust_mantle, &
- b_eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
+!ZN b_R_memory_crust_mantle,b_epsilondev_crust_mantle, &
+!ZN b_eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
+ b_R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,b_deltat,b_displ_crust_mantle, &
b_alphaval,b_betaval,b_gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5) )
@@ -745,8 +756,8 @@
kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
c11store_inner_core,c33store_inner_core,c12store_inner_core, &
c13store_inner_core,c44store_inner_core, &
- b_R_memory_inner_core,b_epsilondev_inner_core, b_eps_trace_over_3_inner_core,&
- one_minus_sum_beta_inner_core, &
+!ZN b_R_memory_inner_core,b_epsilondev_inner_core, b_eps_trace_over_3_inner_core,&
+ b_R_memory_inner_core,one_minus_sum_beta_inner_core,b_deltat,b_veloc_inner_core, &
b_alphaval,b_betaval,b_gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
@@ -779,8 +790,8 @@
kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
c11store_inner_core,c33store_inner_core,c12store_inner_core, &
c13store_inner_core,c44store_inner_core, &
- b_R_memory_inner_core,b_epsilondev_inner_core, b_eps_trace_over_3_inner_core,&
- one_minus_sum_beta_inner_core, &
+!ZN b_R_memory_inner_core,b_epsilondev_inner_core, b_eps_trace_over_3_inner_core,&
+ b_R_memory_inner_core,one_minus_sum_beta_inner_core,b_deltat,b_veloc_inner_core, &
b_alphaval,b_betaval,b_gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
@@ -895,7 +906,7 @@
b_displ_inner_core,b_veloc_inner_core,b_accel_inner_core, &
b_displ_outer_core,b_veloc_outer_core,b_accel_outer_core, &
b_R_memory_crust_mantle,b_R_memory_inner_core, &
- b_epsilondev_crust_mantle,b_epsilondev_inner_core, &
+!ZN b_epsilondev_crust_mantle,b_epsilondev_inner_core, &
b_A_array_rotation,b_B_array_rotation,LOCAL_PATH)
endif
@@ -977,3 +988,4 @@
!-------------------------------------------------------------------------------------------------
!
+
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/part2_undo_att.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/part2_undo_att.f90 2013-06-15 11:05:55 UTC (rev 22309)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/part2_undo_att.f90 2013-06-15 15:41:35 UTC (rev 22310)
@@ -39,8 +39,11 @@
! and output timestamp file to check that simulation is running fine
if(mod(it,NTSTEP_BETWEEN_OUTPUT_INFO) == 0 .or. it == it_begin+4 .or. it == it_end) then
if (SIMULATION_TYPE == 3) then
+!ZN call check_simulation_stability(it,b_displ_crust_mantle,b_displ_inner_core,b_displ_outer_core, &
+!ZN eps_trace_over_3_crust_mantle,epsilondev_crust_mantle, &
+!ZN SIMULATION_TYPE,OUTPUT_FILES,time_start,DT,t0,NSTEP, &
+!ZN it_begin,it_end,NUMBER_OF_THIS_RUN,NUMBER_OF_RUNS,myrank)
call check_simulation_stability(it,b_displ_crust_mantle,b_displ_inner_core,b_displ_outer_core, &
- eps_trace_over_3_crust_mantle,epsilondev_crust_mantle, &
SIMULATION_TYPE,OUTPUT_FILES,time_start,DT,t0,NSTEP, &
it_begin,it_end,NUMBER_OF_THIS_RUN,NUMBER_OF_RUNS,myrank)
endif
@@ -378,8 +381,9 @@
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,ispec_is_tiso_crust_mantle, &
- b_R_memory_crust_mantle,b_epsilondev_crust_mantle, &
- b_eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
+!ZN b_R_memory_crust_mantle,b_epsilondev_crust_mantle, &
+!ZN b_eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
+ b_R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,b_deltat,b_veloc_crust_mantle, &
b_alphaval,b_betaval,b_gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5) )
@@ -419,8 +423,9 @@
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,ispec_is_tiso_crust_mantle, &
- b_R_memory_crust_mantle,b_epsilondev_crust_mantle, &
- b_eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
+!ZN b_R_memory_crust_mantle,b_epsilondev_crust_mantle, &
+!ZN b_eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
+ b_R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,b_deltat,b_displ_crust_mantle, &
b_alphaval,b_betaval,b_gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5) )
@@ -456,8 +461,8 @@
kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
c11store_inner_core,c33store_inner_core,c12store_inner_core, &
c13store_inner_core,c44store_inner_core, &
- b_R_memory_inner_core,b_epsilondev_inner_core, b_eps_trace_over_3_inner_core,&
- one_minus_sum_beta_inner_core, &
+!ZN b_R_memory_inner_core,b_epsilondev_inner_core, b_eps_trace_over_3_inner_core,&
+ b_R_memory_inner_core,one_minus_sum_beta_inner_core,b_deltat,b_veloc_inner_core, &
b_alphaval,b_betaval,b_gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
@@ -490,8 +495,8 @@
kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
c11store_inner_core,c33store_inner_core,c12store_inner_core, &
c13store_inner_core,c44store_inner_core, &
- b_R_memory_inner_core,b_epsilondev_inner_core, b_eps_trace_over_3_inner_core,&
- one_minus_sum_beta_inner_core, &
+!ZN b_R_memory_inner_core,b_epsilondev_inner_core, b_eps_trace_over_3_inner_core,&
+ b_R_memory_inner_core,one_minus_sum_beta_inner_core,b_deltat,b_veloc_inner_core, &
b_alphaval,b_betaval,b_gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
@@ -736,8 +741,9 @@
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,ispec_is_tiso_crust_mantle, &
- b_R_memory_crust_mantle,b_epsilondev_crust_mantle, &
- b_eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
+!ZN b_R_memory_crust_mantle,b_epsilondev_crust_mantle, &
+!ZN b_eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
+ b_R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,b_deltat,b_veloc_crust_mantle, &
b_alphaval,b_betaval,b_gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5) )
@@ -777,8 +783,9 @@
c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
ibool_crust_mantle,ispec_is_tiso_crust_mantle, &
- b_R_memory_crust_mantle,b_epsilondev_crust_mantle, &
- b_eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
+!ZN b_R_memory_crust_mantle,b_epsilondev_crust_mantle, &
+!ZN b_eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
+ b_R_memory_crust_mantle,one_minus_sum_beta_crust_mantle,b_deltat,b_displ_crust_mantle, &
b_alphaval,b_betaval,b_gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5) )
@@ -813,8 +820,8 @@
kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
c11store_inner_core,c33store_inner_core,c12store_inner_core, &
c13store_inner_core,c44store_inner_core, &
- b_R_memory_inner_core,b_epsilondev_inner_core, b_eps_trace_over_3_inner_core,&
- one_minus_sum_beta_inner_core, &
+!ZN b_R_memory_inner_core,b_epsilondev_inner_core, b_eps_trace_over_3_inner_core,&
+ b_R_memory_inner_core,one_minus_sum_beta_inner_core,b_deltat,b_veloc_inner_core, &
b_alphaval,b_betaval,b_gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
@@ -847,8 +854,8 @@
kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
c11store_inner_core,c33store_inner_core,c12store_inner_core, &
c13store_inner_core,c44store_inner_core, &
- b_R_memory_inner_core,b_epsilondev_inner_core, b_eps_trace_over_3_inner_core,&
- one_minus_sum_beta_inner_core, &
+!ZN b_R_memory_inner_core,b_epsilondev_inner_core, b_eps_trace_over_3_inner_core,&
+ b_R_memory_inner_core,one_minus_sum_beta_inner_core,b_deltat,b_veloc_inner_core, &
b_alphaval,b_betaval,b_gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
@@ -964,7 +971,7 @@
b_displ_inner_core,b_veloc_inner_core,b_accel_inner_core, &
b_displ_outer_core,b_veloc_outer_core,b_accel_outer_core, &
b_R_memory_crust_mantle,b_R_memory_inner_core, &
- b_epsilondev_crust_mantle,b_epsilondev_inner_core, &
+!ZN b_epsilondev_crust_mantle,b_epsilondev_inner_core, &
b_A_array_rotation,b_B_array_rotation,LOCAL_PATH)
endif
endif
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/part3_kernel_computation.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/part3_kernel_computation.f90 2013-06-15 11:05:55 UTC (rev 22309)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/part3_kernel_computation.f90 2013-06-15 15:41:35 UTC (rev 22310)
@@ -9,9 +9,12 @@
rho_kl_crust_mantle,beta_kl_crust_mantle, &
alpha_kl_crust_mantle,cijkl_kl_crust_mantle, &
accel_crust_mantle,b_displ_crust_mantle, &
- epsilondev_crust_mantle,b_epsilondev_crust_mantle, &
- eps_trace_over_3_crust_mantle,b_eps_trace_over_3_crust_mantle, &
- deltat)
+!ZN epsilondev_crust_mantle,b_epsilondev_crust_mantle, &
+!ZN eps_trace_over_3_crust_mantle,b_eps_trace_over_3_crust_mantle, &
+ deltat,displ_crust_mantle,hprime_xx,hprime_xxT,&
+ xix_crust_mantle,xiy_crust_mantle,xiz_crust_mantle,&
+ etax_crust_mantle,etay_crust_mantle,etaz_crust_mantle,&
+ gammax_crust_mantle,gammay_crust_mantle,gammaz_crust_mantle)
! outer core
call compute_kernels_outer_core(ibool_outer_core, &
@@ -34,9 +37,12 @@
rho_kl_inner_core,beta_kl_inner_core, &
alpha_kl_inner_core, &
accel_inner_core,b_displ_inner_core, &
- epsilondev_inner_core,b_epsilondev_inner_core, &
- eps_trace_over_3_inner_core,b_eps_trace_over_3_inner_core, &
- deltat)
+!ZN epsilondev_inner_core,b_epsilondev_inner_core, &
+!ZN eps_trace_over_3_inner_core,b_eps_trace_over_3_inner_core, &
+ deltat,displ_inner_core,hprime_xx,hprime_xxT,&
+ xix_inner_core,xiy_inner_core,xiz_inner_core,&
+ etax_inner_core,etay_inner_core,etaz_inner_core,&
+ gammax_inner_core,gammay_inner_core,gammaz_inner_core)
! NOISE TOMOGRAPHY --- source strength kernel
if (NOISE_TOMOGRAPHY == 3) &
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/check_simulation_stability.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/check_simulation_stability.f90 2013-06-15 11:05:55 UTC (rev 22309)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/check_simulation_stability.f90 2013-06-15 15:41:35 UTC (rev 22310)
@@ -26,7 +26,7 @@
!=====================================================================
subroutine check_simulation_stability(it,displ_crust_mantle,displ_inner_core,displ_outer_core, &
- eps_trace_over_3_crust_mantle,epsilondev_crust_mantle, &
+!ZN eps_trace_over_3_crust_mantle,epsilondev_crust_mantle, &
SIMULATION_TYPE,OUTPUT_FILES,time_start,DT,t0,NSTEP, &
it_begin,it_end,NUMBER_OF_THIS_RUN,NUMBER_OF_RUNS,myrank)
@@ -45,10 +45,10 @@
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_INNER_CORE) :: displ_inner_core
real(kind=CUSTOM_REAL), dimension(NGLOB_OUTER_CORE) :: displ_outer_core
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STRAIN_ONLY) :: &
- eps_trace_over_3_crust_mantle
- real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STR_OR_ATT) :: &
- epsilondev_crust_mantle
+!ZN real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STRAIN_ONLY) :: &
+!ZN eps_trace_over_3_crust_mantle
+!ZN real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STR_OR_ATT) :: &
+!ZN epsilondev_crust_mantle
integer SIMULATION_TYPE
character(len=150) OUTPUT_FILES
@@ -58,7 +58,7 @@
! local parameters
! maximum of the norm of the displacement and of the potential in the fluid
real(kind=CUSTOM_REAL) Usolidnorm,Usolidnorm_all,Ufluidnorm,Ufluidnorm_all
- real(kind=CUSTOM_REAL) Strain_norm,Strain_norm_all,strain2_norm,strain2_norm_all
+!ZN real(kind=CUSTOM_REAL) Strain_norm,Strain_norm_all,strain2_norm,strain2_norm_all
! names of the data files for all the processors in MPI
character(len=150) outputname
! timer MPI
@@ -101,14 +101,14 @@
call MPI_REDUCE(Ufluidnorm,Ufluidnorm_all,1,CUSTOM_MPI_TYPE,MPI_MAX,0, &
MPI_COMM_WORLD,ier)
- if (COMPUTE_AND_STORE_STRAIN) then
- Strain_norm = maxval(abs(eps_trace_over_3_crust_mantle))
- strain2_norm= maxval(abs(epsilondev_crust_mantle))
- call MPI_REDUCE(Strain_norm,Strain_norm_all,1,CUSTOM_MPI_TYPE,MPI_MAX,0, &
- MPI_COMM_WORLD,ier)
- call MPI_REDUCE(Strain2_norm,Strain2_norm_all,1,CUSTOM_MPI_TYPE,MPI_MAX,0, &
- MPI_COMM_WORLD,ier)
- endif
+!ZN if (COMPUTE_AND_STORE_STRAIN) then
+!ZN Strain_norm = maxval(abs(eps_trace_over_3_crust_mantle))
+!ZN strain2_norm= maxval(abs(epsilondev_crust_mantle))
+!ZN call MPI_REDUCE(Strain_norm,Strain_norm_all,1,CUSTOM_MPI_TYPE,MPI_MAX,0, &
+!ZN MPI_COMM_WORLD,ier)
+!ZN call MPI_REDUCE(Strain2_norm,Strain2_norm_all,1,CUSTOM_MPI_TYPE,MPI_MAX,0, &
+!ZN MPI_COMM_WORLD,ier)
+!ZN endif
if(myrank == 0) then
@@ -127,10 +127,10 @@
!! DK DK UNDO_ATT
! if(COMPUTE_AND_STORE_STRAIN) then
- if(SIMULATION_TYPE == 1 .and. COMPUTE_AND_STORE_STRAIN) then
- write(IMAIN,*) 'Max of strain, eps_trace_over_3_crust_mantle =',Strain_norm_all
- write(IMAIN,*) 'Max of strain, epsilondev_crust_mantle =',Strain2_norm_all
- endif
+!ZN if(SIMULATION_TYPE == 1 .and. COMPUTE_AND_STORE_STRAIN) then
+!ZN write(IMAIN,*) 'Max of strain, eps_trace_over_3_crust_mantle =',Strain_norm_all
+!ZN write(IMAIN,*) 'Max of strain, epsilondev_crust_mantle =',Strain2_norm_all
+!ZN endif
! information about the current run only
SHOW_SEPARATE_RUN_INFORMATION = NUMBER_OF_RUNS > 1 .and. NUMBER_OF_THIS_RUN < NUMBER_OF_RUNS
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/compute_element.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/compute_element.F90 2013-06-15 11:05:55 UTC (rev 22309)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/compute_element.F90 2013-06-15 15:41:35 UTC (rev 22310)
@@ -32,7 +32,8 @@
wgll_cube, &
kappavstore,muvstore, &
ibool, &
- R_memory,epsilon_trace_over_3, &
+!ZN R_memory,epsilon_trace_over_3, &
+ R_memory, & !ZN
one_minus_sum_beta,vx,vy,vz,vnspec, &
tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,rho_s_H)
@@ -72,7 +73,7 @@
! memory variables R_ij are stored at the local rather than global level
! to allow for optimization of cache access by compiler
real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: R_memory
- real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: epsilon_trace_over_3
+!ZN real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: epsilon_trace_over_3
integer :: vx,vy,vz,vnspec
real(kind=CUSTOM_REAL), dimension(vx, vy, vz, vnspec) :: one_minus_sum_beta
@@ -110,7 +111,7 @@
double precision factor,sx_l,sy_l,sz_l,gxl,gyl,gzl
double precision Hxxl,Hyyl,Hzzl,Hxyl,Hxzl,Hyzl
- integer :: ispec_strain
+!ZN integer :: ispec_strain
integer :: i,j,k
integer :: int_radius
integer :: iglob1
@@ -159,15 +160,15 @@
! compute deviatoric strain
if (COMPUTE_AND_STORE_STRAIN) then
templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
- if(NSPEC_CRUST_MANTLE_STRAIN_ONLY == 1) then
- ispec_strain = 1
-!$OMP CRITICAL
- epsilon_trace_over_3(i,j,k,ispec_strain) = templ
-!$OMP END CRITICAL
- else
- ispec_strain = ispec
- epsilon_trace_over_3(i,j,k,ispec_strain) = templ
- endif
+!ZN if(NSPEC_CRUST_MANTLE_STRAIN_ONLY == 1) then
+!ZN ispec_strain = 1
+!ZN!$OMP CRITICAL
+!ZN epsilon_trace_over_3(i,j,k,ispec_strain) = templ
+!ZN!$OMP END CRITICAL
+!ZN else
+!ZN ispec_strain = ispec
+!ZN epsilon_trace_over_3(i,j,k,ispec_strain) = templ
+!ZN endif
epsilondev_loc(1,i,j,k) = duxdxl - templ
epsilondev_loc(2,i,j,k) = duydyl - templ
epsilondev_loc(3,i,j,k) = 0.5 * duxdyl_plus_duydxl
@@ -358,7 +359,8 @@
wgll_cube, &
kappavstore,kappahstore,muvstore,muhstore,eta_anisostore, &
ibool, &
- R_memory,epsilon_trace_over_3, &
+!ZN R_memory,epsilon_trace_over_3, &
+ R_memory, & !ZN
one_minus_sum_beta,vx,vy,vz,vnspec, &
tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,rho_s_H)
@@ -400,7 +402,7 @@
! memory variables R_ij are stored at the local rather than global level
! to allow for optimization of cache access by compiler
real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: R_memory
- real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: epsilon_trace_over_3
+!ZN real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: epsilon_trace_over_3
integer vx,vy,vz,vnspec
@@ -454,7 +456,7 @@
double precision Hxxl,Hyyl,Hzzl,Hxyl,Hxzl,Hyzl
real(kind=CUSTOM_REAL) sigma_yx,sigma_zx,sigma_zy
- integer :: ispec_strain
+!ZN integer :: ispec_strain
integer :: i,j,k
integer :: int_radius
integer :: iglob1
@@ -503,15 +505,15 @@
! compute deviatoric strain
if (COMPUTE_AND_STORE_STRAIN) then
templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
- if(NSPEC_CRUST_MANTLE_STRAIN_ONLY == 1) then
- ispec_strain = 1
-!$OMP CRITICAL
- epsilon_trace_over_3(i,j,k,ispec_strain) = templ
-!$OMP END CRITICAL
- else
- ispec_strain = ispec
- epsilon_trace_over_3(i,j,k,ispec_strain) = templ
- endif
+!ZN if(NSPEC_CRUST_MANTLE_STRAIN_ONLY == 1) then
+!ZN ispec_strain = 1
+!ZN!$OMP CRITICAL
+!ZN epsilon_trace_over_3(i,j,k,ispec_strain) = templ
+!ZN!$OMP END CRITICAL
+!ZN else
+!ZN ispec_strain = ispec
+!ZN epsilon_trace_over_3(i,j,k,ispec_strain) = templ
+!ZN endif
epsilondev_loc(1,i,j,k) = duxdxl - templ
epsilondev_loc(2,i,j,k) = duydyl - templ
epsilondev_loc(3,i,j,k) = 0.5 * duxdyl_plus_duydxl
@@ -895,7 +897,8 @@
c23store,c24store,c25store,c26store,c33store,c34store,c35store, &
c36store,c44store,c45store,c46store,c55store,c56store,c66store, &
ibool, &
- R_memory,epsilon_trace_over_3, &
+!ZN R_memory,epsilon_trace_over_3, &
+ R_memory, &
one_minus_sum_beta,vx,vy,vz,vnspec, &
tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,rho_s_H)
@@ -937,7 +940,7 @@
! memory variables R_ij are stored at the local rather than global level
! to allow for optimization of cache access by compiler
real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: R_memory
- real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: epsilon_trace_over_3
+!ZN real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: epsilon_trace_over_3
integer vx,vy,vz,vnspec
@@ -978,7 +981,7 @@
double precision Hxxl,Hyyl,Hzzl,Hxyl,Hxzl,Hyzl
real(kind=CUSTOM_REAL) sigma_yx,sigma_zx,sigma_zy
- integer :: ispec_strain
+!ZN integer :: ispec_strain
integer :: i,j,k
integer :: int_radius
integer :: iglob1
@@ -1027,15 +1030,15 @@
! compute deviatoric strain
if (COMPUTE_AND_STORE_STRAIN) then
templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
- if(NSPEC_CRUST_MANTLE_STRAIN_ONLY == 1) then
- ispec_strain = 1
-!$OMP CRITICAL
- epsilon_trace_over_3(i,j,k,ispec_strain) = templ
-!$OMP END CRITICAL
- else
- ispec_strain = ispec
- epsilon_trace_over_3(i,j,k,ispec_strain) = templ
- endif
+!ZN if(NSPEC_CRUST_MANTLE_STRAIN_ONLY == 1) then
+!ZN ispec_strain = 1
+!ZN!$OMP CRITICAL
+!ZN epsilon_trace_over_3(i,j,k,ispec_strain) = templ
+!ZN!$OMP END CRITICAL
+!ZN else
+!ZN ispec_strain = ispec
+!ZN epsilon_trace_over_3(i,j,k,ispec_strain) = templ
+!ZN endif
epsilondev_loc(1,i,j,k) = duxdxl - templ
epsilondev_loc(2,i,j,k) = duydyl - templ
epsilondev_loc(3,i,j,k) = 0.5 * duxdyl_plus_duydxl
@@ -1293,7 +1296,8 @@
vx,vy,vz,vnspec,factor_common, &
alphaval,betaval,gammaval, &
c44store,muvstore, &
- epsilondev,epsilondev_loc)
+!ZN epsilondev,epsilondev_loc)
+ epsilondev_loc_nplus1,epsilondev_loc)
! crust mantle
! update memory variables based upon the Runge-Kutta scheme
@@ -1331,7 +1335,8 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPECMAX_ANISO_MANTLE) :: c44store
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPECMAX_ISO_MANTLE) :: muvstore
- real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: epsilondev
+!ZN real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: epsilondev
+ real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc_nplus1
real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc
! local parameters
@@ -1357,9 +1362,12 @@
endif
do i_memory = 1,5
+!ZN R_memory(i_memory,i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_memory(i_memory,i_SLS,:,:,:,ispec) &
+!ZN + factor_common_c44_muv(:,:,:) &
+!ZN * (betaval(i_SLS) * epsilondev(i_memory,:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(i_memory,:,:,:))
R_memory(i_memory,i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_memory(i_memory,i_SLS,:,:,:,ispec) &
+ factor_common_c44_muv(:,:,:) &
- * (betaval(i_SLS) * epsilondev(i_memory,:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(i_memory,:,:,:))
+ * (betaval(i_SLS) * epsilondev_loc(i_memory,:,:,:) + gammaval(i_SLS) * epsilondev_loc_nplus1(i_memory,:,:,:))
enddo
enddo ! i_SLS
@@ -1373,7 +1381,8 @@
vx,vy,vz,vnspec,factor_common, &
alphaval,betaval,gammaval, &
muvstore, &
- epsilondev,epsilondev_loc)
+!ZN epsilondev,epsilondev_loc)
+ epsilondev_loc_nplus1,epsilondev_loc) !ZN
! inner core
! update memory variables based upon the Runge-Kutta scheme
@@ -1410,8 +1419,9 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: muvstore
- real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: epsilondev
+!ZN real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: epsilondev
real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc
+ real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc_nplus1
! local parameters
real(kind=CUSTOM_REAL), dimension(NGLLX, NGLLY, NGLLZ) :: factor_common_use
@@ -1431,9 +1441,658 @@
do i_memory = 1,5
R_memory(i_memory,i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_memory(i_memory,i_SLS,:,:,:,ispec) &
+ muvstore(:,:,:,ispec) * factor_common_use(:,:,:) * &
- (betaval(i_SLS) * epsilondev(i_memory,:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(i_memory,:,:,:))
+!ZN (betaval(i_SLS) * epsilondev(i_memory,:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(i_memory,:,:,:))
+ (betaval(i_SLS) * epsilondev_loc_nplus1(i_memory,:,:,:) + gammaval(i_SLS) * epsilondev_loc(i_memory,:,:,:)) !ZN
enddo
enddo
end subroutine compute_element_att_memory_ic
+
+!
+!--------------------------------------------------------------------------------------------
+!
+ subroutine compute_element_strain_undo_att_Dev(ispec,nglob,nspec,displ,ibool,hprime_xx,hprime_xxT,&
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,epsilondev_loc,eps_trace_over_3_loc)
+
+ implicit none
+ include "constants.h"
+
+ integer :: ispec,nglob,nspec
+ real(kind=CUSTOM_REAL), dimension(NDIM,nglob) :: displ
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx,hprime_xxT
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
+
+ real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: eps_trace_over_3_loc
+
+! local variable
+ integer :: i,j,k,iglob
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
+ tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: dummyx_loc,dummyy_loc,dummyz_loc
+ real(kind=CUSTOM_REAL), dimension(NGLLX,m2) :: B1_m1_m2_5points,B2_m1_m2_5points,B3_m1_m2_5points
+ equivalence(dummyx_loc,B1_m1_m2_5points)
+ equivalence(dummyy_loc,B2_m1_m2_5points)
+ equivalence(dummyz_loc,B3_m1_m2_5points)
+ real(kind=CUSTOM_REAL), dimension(m1,m2) :: C1_m1_m2_5points,C2_m1_m2_5points,C3_m1_m2_5points
+ equivalence(tempx1,C1_m1_m2_5points)
+ equivalence(tempy1,C2_m1_m2_5points)
+ equivalence(tempz1,C3_m1_m2_5points)
+ real(kind=CUSTOM_REAL), dimension(m2,m1) :: &
+ C1_mxm_m2_m1_5points,C2_mxm_m2_m1_5points,C3_mxm_m2_m1_5points
+ equivalence(tempx3,C1_mxm_m2_m1_5points)
+ equivalence(tempy3,C2_mxm_m2_m1_5points)
+ equivalence(tempz3,C3_mxm_m2_m1_5points)
+ real(kind=CUSTOM_REAL), dimension(m2,NGLLX) :: &
+ A1_mxm_m2_m1_5points,A2_mxm_m2_m1_5points,A3_mxm_m2_m1_5points
+ equivalence(dummyx_loc,A1_mxm_m2_m1_5points)
+ equivalence(dummyy_loc,A2_mxm_m2_m1_5points)
+ equivalence(dummyz_loc,A3_mxm_m2_m1_5points)
+
+ real(kind=CUSTOM_REAL) xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
+ real(kind=CUSTOM_REAL) duxdxl,duydyl,duzdzl,duxdyl,duydxl,duzdxl,duxdzl,duzdyl,duydzl,&
+ duxdxl_plus_duydyl,duxdxl_plus_duzdzl,duydyl_plus_duzdzl,&
+ duxdyl_plus_duydxl,duzdxl_plus_duxdzl,duzdyl_plus_duydzl
+
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+ dummyx_loc(i,j,k) = displ(1,iglob)
+ dummyy_loc(i,j,k) = displ(2,iglob)
+ dummyz_loc(i,j,k) = displ(3,iglob)
+ enddo
+ enddo
+ enddo
+
+ do j=1,m2
+ do i=1,m1
+ C1_m1_m2_5points(i,j) = hprime_xx(i,1)*B1_m1_m2_5points(1,j) + &
+ hprime_xx(i,2)*B1_m1_m2_5points(2,j) + &
+ hprime_xx(i,3)*B1_m1_m2_5points(3,j) + &
+ hprime_xx(i,4)*B1_m1_m2_5points(4,j) + &
+ hprime_xx(i,5)*B1_m1_m2_5points(5,j)
+
+ C2_m1_m2_5points(i,j) = hprime_xx(i,1)*B2_m1_m2_5points(1,j) + &
+ hprime_xx(i,2)*B2_m1_m2_5points(2,j) + &
+ hprime_xx(i,3)*B2_m1_m2_5points(3,j) + &
+ hprime_xx(i,4)*B2_m1_m2_5points(4,j) + &
+ hprime_xx(i,5)*B2_m1_m2_5points(5,j)
+
+ C3_m1_m2_5points(i,j) = hprime_xx(i,1)*B3_m1_m2_5points(1,j) + &
+ hprime_xx(i,2)*B3_m1_m2_5points(2,j) + &
+ hprime_xx(i,3)*B3_m1_m2_5points(3,j) + &
+ hprime_xx(i,4)*B3_m1_m2_5points(4,j) + &
+ hprime_xx(i,5)*B3_m1_m2_5points(5,j)
+ enddo
+ enddo
+
+ do j=1,m1
+ do i=1,m1
+ ! for efficiency it is better to leave this loop on k inside, it leads to slightly faster code
+ do k = 1,NGLLX
+ tempx2(i,j,k) = dummyx_loc(i,1,k)*hprime_xxT(1,j) + &
+ dummyx_loc(i,2,k)*hprime_xxT(2,j) + &
+ dummyx_loc(i,3,k)*hprime_xxT(3,j) + &
+ dummyx_loc(i,4,k)*hprime_xxT(4,j) + &
+ dummyx_loc(i,5,k)*hprime_xxT(5,j)
+
+ tempy2(i,j,k) = dummyy_loc(i,1,k)*hprime_xxT(1,j) + &
+ dummyy_loc(i,2,k)*hprime_xxT(2,j) + &
+ dummyy_loc(i,3,k)*hprime_xxT(3,j) + &
+ dummyy_loc(i,4,k)*hprime_xxT(4,j) + &
+ dummyy_loc(i,5,k)*hprime_xxT(5,j)
+
+ tempz2(i,j,k) = dummyz_loc(i,1,k)*hprime_xxT(1,j) + &
+ dummyz_loc(i,2,k)*hprime_xxT(2,j) + &
+ dummyz_loc(i,3,k)*hprime_xxT(3,j) + &
+ dummyz_loc(i,4,k)*hprime_xxT(4,j) + &
+ dummyz_loc(i,5,k)*hprime_xxT(5,j)
+ enddo
+ enddo
+ enddo
+
+ do j=1,m1
+ do i=1,m2
+ C1_mxm_m2_m1_5points(i,j) = A1_mxm_m2_m1_5points(i,1)*hprime_xxT(1,j) + &
+ A1_mxm_m2_m1_5points(i,2)*hprime_xxT(2,j) + &
+ A1_mxm_m2_m1_5points(i,3)*hprime_xxT(3,j) + &
+ A1_mxm_m2_m1_5points(i,4)*hprime_xxT(4,j) + &
+ A1_mxm_m2_m1_5points(i,5)*hprime_xxT(5,j)
+
+ C2_mxm_m2_m1_5points(i,j) = A2_mxm_m2_m1_5points(i,1)*hprime_xxT(1,j) + &
+ A2_mxm_m2_m1_5points(i,2)*hprime_xxT(2,j) + &
+ A2_mxm_m2_m1_5points(i,3)*hprime_xxT(3,j) + &
+ A2_mxm_m2_m1_5points(i,4)*hprime_xxT(4,j) + &
+ A2_mxm_m2_m1_5points(i,5)*hprime_xxT(5,j)
+
+ C3_mxm_m2_m1_5points(i,j) = A3_mxm_m2_m1_5points(i,1)*hprime_xxT(1,j) + &
+ A3_mxm_m2_m1_5points(i,2)*hprime_xxT(2,j) + &
+ A3_mxm_m2_m1_5points(i,3)*hprime_xxT(3,j) + &
+ A3_mxm_m2_m1_5points(i,4)*hprime_xxT(4,j) + &
+ A3_mxm_m2_m1_5points(i,5)*hprime_xxT(5,j)
+ enddo
+ enddo
+
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ ! get derivatives of ux, uy and uz with respect to x, y and z
+ xixl = xix(i,j,k,ispec)
+ xiyl = xiy(i,j,k,ispec)
+ xizl = xiz(i,j,k,ispec)
+ etaxl = etax(i,j,k,ispec)
+ etayl = etay(i,j,k,ispec)
+ etazl = etaz(i,j,k,ispec)
+ gammaxl = gammax(i,j,k,ispec)
+ gammayl = gammay(i,j,k,ispec)
+ gammazl = gammaz(i,j,k,ispec)
+
+ ! compute the jacobian
+ jacobianl = 1._CUSTOM_REAL / (xixl*(etayl*gammazl-etazl*gammayl) &
+ - xiyl*(etaxl*gammazl-etazl*gammaxl) &
+ + xizl*(etaxl*gammayl-etayl*gammaxl))
+
+ duxdxl = xixl*tempx1(i,j,k) + etaxl*tempx2(i,j,k) + gammaxl*tempx3(i,j,k)
+ duxdyl = xiyl*tempx1(i,j,k) + etayl*tempx2(i,j,k) + gammayl*tempx3(i,j,k)
+ duxdzl = xizl*tempx1(i,j,k) + etazl*tempx2(i,j,k) + gammazl*tempx3(i,j,k)
+
+ duydxl = xixl*tempy1(i,j,k) + etaxl*tempy2(i,j,k) + gammaxl*tempy3(i,j,k)
+ duydyl = xiyl*tempy1(i,j,k) + etayl*tempy2(i,j,k) + gammayl*tempy3(i,j,k)
+ duydzl = xizl*tempy1(i,j,k) + etazl*tempy2(i,j,k) + gammazl*tempy3(i,j,k)
+
+ duzdxl = xixl*tempz1(i,j,k) + etaxl*tempz2(i,j,k) + gammaxl*tempz3(i,j,k)
+ duzdyl = xiyl*tempz1(i,j,k) + etayl*tempz2(i,j,k) + gammayl*tempz3(i,j,k)
+ duzdzl = xizl*tempz1(i,j,k) + etazl*tempz2(i,j,k) + gammazl*tempz3(i,j,k)
+
+ ! precompute some sums to save CPU time
+ duxdxl_plus_duydyl = duxdxl + duydyl
+ duxdxl_plus_duzdzl = duxdxl + duzdzl
+ duydyl_plus_duzdzl = duydyl + duzdzl
+ duxdyl_plus_duydxl = duxdyl + duydxl
+ duzdxl_plus_duxdzl = duzdxl + duxdzl
+ duzdyl_plus_duydzl = duzdyl + duydzl
+
+ enddo
+ enddo
+ enddo
+
+ eps_trace_over_3_loc(i,j,k) = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ epsilondev_loc(1,i,j,k) = duxdxl - eps_trace_over_3_loc(i,j,k)
+ epsilondev_loc(2,i,j,k) = duydyl - eps_trace_over_3_loc(i,j,k)
+ epsilondev_loc(3,i,j,k) = 0.5 * duxdyl_plus_duydxl
+ epsilondev_loc(4,i,j,k) = 0.5 * duzdxl_plus_duxdzl
+ epsilondev_loc(5,i,j,k) = 0.5 * duzdyl_plus_duydzl
+
+ end subroutine compute_element_strain_undo_att_Dev
+
+!
+!--------------------------------------------------------------------------------------------
+!
+ subroutine compute_element_strain_att_Dev(ispec,nglob,nspec,displ,veloc,deltat,ibool,hprime_xx,hprime_xxT,&
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,epsilondev_loc_nplus1)
+
+ implicit none
+ include "constants.h"
+
+ integer :: ispec,nglob,nspec
+ real(kind=CUSTOM_REAL) :: deltat
+ real(kind=CUSTOM_REAL), dimension(NDIM,nglob) :: displ,veloc
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx,hprime_xxT
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
+
+ real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc_nplus1
+
+! local variable
+ integer :: i,j,k,iglob
+ real(kind=CUSTOM_REAL) :: templ
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
+ tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: dummyx_loc,dummyy_loc,dummyz_loc
+ real(kind=CUSTOM_REAL), dimension(NGLLX,m2) :: B1_m1_m2_5points,B2_m1_m2_5points,B3_m1_m2_5points
+ equivalence(dummyx_loc,B1_m1_m2_5points)
+ equivalence(dummyy_loc,B2_m1_m2_5points)
+ equivalence(dummyz_loc,B3_m1_m2_5points)
+ real(kind=CUSTOM_REAL), dimension(m1,m2) :: C1_m1_m2_5points,C2_m1_m2_5points,C3_m1_m2_5points
+ equivalence(tempx1,C1_m1_m2_5points)
+ equivalence(tempy1,C2_m1_m2_5points)
+ equivalence(tempz1,C3_m1_m2_5points)
+ real(kind=CUSTOM_REAL), dimension(m2,m1) :: &
+ C1_mxm_m2_m1_5points,C2_mxm_m2_m1_5points,C3_mxm_m2_m1_5points
+ equivalence(tempx3,C1_mxm_m2_m1_5points)
+ equivalence(tempy3,C2_mxm_m2_m1_5points)
+ equivalence(tempz3,C3_mxm_m2_m1_5points)
+ real(kind=CUSTOM_REAL), dimension(m2,NGLLX) :: &
+ A1_mxm_m2_m1_5points,A2_mxm_m2_m1_5points,A3_mxm_m2_m1_5points
+ equivalence(dummyx_loc,A1_mxm_m2_m1_5points)
+ equivalence(dummyy_loc,A2_mxm_m2_m1_5points)
+ equivalence(dummyz_loc,A3_mxm_m2_m1_5points)
+
+ real(kind=CUSTOM_REAL) xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
+ real(kind=CUSTOM_REAL) duxdxl,duydyl,duzdzl,duxdyl,duydxl,duzdxl,duxdzl,duzdyl,duydzl,&
+ duxdxl_plus_duydyl,duxdxl_plus_duzdzl,duydyl_plus_duzdzl,&
+ duxdyl_plus_duydxl,duzdxl_plus_duxdzl,duzdyl_plus_duydzl
+
+
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+ dummyx_loc(i,j,k) = displ(1,iglob) + deltat * veloc(1,iglob)
+ dummyy_loc(i,j,k) = displ(2,iglob) + deltat * veloc(2,iglob)
+ dummyz_loc(i,j,k) = displ(3,iglob) + deltat * veloc(3,iglob)
+ enddo
+ enddo
+ enddo
+
+ do j=1,m2
+ do i=1,m1
+ C1_m1_m2_5points(i,j) = hprime_xx(i,1)*B1_m1_m2_5points(1,j) + &
+ hprime_xx(i,2)*B1_m1_m2_5points(2,j) + &
+ hprime_xx(i,3)*B1_m1_m2_5points(3,j) + &
+ hprime_xx(i,4)*B1_m1_m2_5points(4,j) + &
+ hprime_xx(i,5)*B1_m1_m2_5points(5,j)
+
+ C2_m1_m2_5points(i,j) = hprime_xx(i,1)*B2_m1_m2_5points(1,j) + &
+ hprime_xx(i,2)*B2_m1_m2_5points(2,j) + &
+ hprime_xx(i,3)*B2_m1_m2_5points(3,j) + &
+ hprime_xx(i,4)*B2_m1_m2_5points(4,j) + &
+ hprime_xx(i,5)*B2_m1_m2_5points(5,j)
+
+ C3_m1_m2_5points(i,j) = hprime_xx(i,1)*B3_m1_m2_5points(1,j) + &
+ hprime_xx(i,2)*B3_m1_m2_5points(2,j) + &
+ hprime_xx(i,3)*B3_m1_m2_5points(3,j) + &
+ hprime_xx(i,4)*B3_m1_m2_5points(4,j) + &
+ hprime_xx(i,5)*B3_m1_m2_5points(5,j)
+ enddo
+ enddo
+
+ do j=1,m1
+ do i=1,m1
+ ! for efficiency it is better to leave this loop on k inside, it leads to slightly faster code
+ do k = 1,NGLLX
+ tempx2(i,j,k) = dummyx_loc(i,1,k)*hprime_xxT(1,j) + &
+ dummyx_loc(i,2,k)*hprime_xxT(2,j) + &
+ dummyx_loc(i,3,k)*hprime_xxT(3,j) + &
+ dummyx_loc(i,4,k)*hprime_xxT(4,j) + &
+ dummyx_loc(i,5,k)*hprime_xxT(5,j)
+
+ tempy2(i,j,k) = dummyy_loc(i,1,k)*hprime_xxT(1,j) + &
+ dummyy_loc(i,2,k)*hprime_xxT(2,j) + &
+ dummyy_loc(i,3,k)*hprime_xxT(3,j) + &
+ dummyy_loc(i,4,k)*hprime_xxT(4,j) + &
+ dummyy_loc(i,5,k)*hprime_xxT(5,j)
+
+ tempz2(i,j,k) = dummyz_loc(i,1,k)*hprime_xxT(1,j) + &
+ dummyz_loc(i,2,k)*hprime_xxT(2,j) + &
+ dummyz_loc(i,3,k)*hprime_xxT(3,j) + &
+ dummyz_loc(i,4,k)*hprime_xxT(4,j) + &
+ dummyz_loc(i,5,k)*hprime_xxT(5,j)
+ enddo
+ enddo
+ enddo
+
+ do j=1,m1
+ do i=1,m2
+ C1_mxm_m2_m1_5points(i,j) = A1_mxm_m2_m1_5points(i,1)*hprime_xxT(1,j) + &
+ A1_mxm_m2_m1_5points(i,2)*hprime_xxT(2,j) + &
+ A1_mxm_m2_m1_5points(i,3)*hprime_xxT(3,j) + &
+ A1_mxm_m2_m1_5points(i,4)*hprime_xxT(4,j) + &
+ A1_mxm_m2_m1_5points(i,5)*hprime_xxT(5,j)
+
+ C2_mxm_m2_m1_5points(i,j) = A2_mxm_m2_m1_5points(i,1)*hprime_xxT(1,j) + &
+ A2_mxm_m2_m1_5points(i,2)*hprime_xxT(2,j) + &
+ A2_mxm_m2_m1_5points(i,3)*hprime_xxT(3,j) + &
+ A2_mxm_m2_m1_5points(i,4)*hprime_xxT(4,j) + &
+ A2_mxm_m2_m1_5points(i,5)*hprime_xxT(5,j)
+
+ C3_mxm_m2_m1_5points(i,j) = A3_mxm_m2_m1_5points(i,1)*hprime_xxT(1,j) + &
+ A3_mxm_m2_m1_5points(i,2)*hprime_xxT(2,j) + &
+ A3_mxm_m2_m1_5points(i,3)*hprime_xxT(3,j) + &
+ A3_mxm_m2_m1_5points(i,4)*hprime_xxT(4,j) + &
+ A3_mxm_m2_m1_5points(i,5)*hprime_xxT(5,j)
+ enddo
+ enddo
+
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ ! get derivatives of ux, uy and uz with respect to x, y and z
+ xixl = xix(i,j,k,ispec)
+ xiyl = xiy(i,j,k,ispec)
+ xizl = xiz(i,j,k,ispec)
+ etaxl = etax(i,j,k,ispec)
+ etayl = etay(i,j,k,ispec)
+ etazl = etaz(i,j,k,ispec)
+ gammaxl = gammax(i,j,k,ispec)
+ gammayl = gammay(i,j,k,ispec)
+ gammazl = gammaz(i,j,k,ispec)
+
+ ! compute the jacobian
+ jacobianl = 1._CUSTOM_REAL / (xixl*(etayl*gammazl-etazl*gammayl) &
+ - xiyl*(etaxl*gammazl-etazl*gammaxl) &
+ + xizl*(etaxl*gammayl-etayl*gammaxl))
+
+ duxdxl = xixl*tempx1(i,j,k) + etaxl*tempx2(i,j,k) + gammaxl*tempx3(i,j,k)
+ duxdyl = xiyl*tempx1(i,j,k) + etayl*tempx2(i,j,k) + gammayl*tempx3(i,j,k)
+ duxdzl = xizl*tempx1(i,j,k) + etazl*tempx2(i,j,k) + gammazl*tempx3(i,j,k)
+
+ duydxl = xixl*tempy1(i,j,k) + etaxl*tempy2(i,j,k) + gammaxl*tempy3(i,j,k)
+ duydyl = xiyl*tempy1(i,j,k) + etayl*tempy2(i,j,k) + gammayl*tempy3(i,j,k)
+ duydzl = xizl*tempy1(i,j,k) + etazl*tempy2(i,j,k) + gammazl*tempy3(i,j,k)
+
+ duzdxl = xixl*tempz1(i,j,k) + etaxl*tempz2(i,j,k) + gammaxl*tempz3(i,j,k)
+ duzdyl = xiyl*tempz1(i,j,k) + etayl*tempz2(i,j,k) + gammayl*tempz3(i,j,k)
+ duzdzl = xizl*tempz1(i,j,k) + etazl*tempz2(i,j,k) + gammazl*tempz3(i,j,k)
+
+ ! precompute some sums to save CPU time
+ duxdxl_plus_duydyl = duxdxl + duydyl
+ duxdxl_plus_duzdzl = duxdxl + duzdzl
+ duydyl_plus_duzdzl = duydyl + duzdzl
+ duxdyl_plus_duydxl = duxdyl + duydxl
+ duzdxl_plus_duxdzl = duzdxl + duxdzl
+ duzdyl_plus_duydzl = duzdyl + duydzl
+
+ enddo
+ enddo
+ enddo
+
+ templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ epsilondev_loc_nplus1(1,i,j,k) = duxdxl - templ
+ epsilondev_loc_nplus1(2,i,j,k) = duydyl - templ
+ epsilondev_loc_nplus1(3,i,j,k) = 0.5 * duxdyl_plus_duydxl
+ epsilondev_loc_nplus1(4,i,j,k) = 0.5 * duzdxl_plus_duxdzl
+ epsilondev_loc_nplus1(5,i,j,k) = 0.5 * duzdyl_plus_duydzl
+
+ end subroutine compute_element_strain_att_Dev
+!=====================================================================
+
+ subroutine compute_element_strain_undo_att_noDev(ispec,nglob,nspec,displ,hprime_xx,hprime_yy,hprime_zz,ibool,&
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,epsilondev_loc,eps_trace_over_3_loc)
+
+ implicit none
+
+ include "constants.h"
+
+! include values created by the mesher
+! done for performance only using static allocation to allow for loop unrolling
+ include "OUTPUT_FILES/values_from_mesher.h"
+
+ integer ispec,NSPEC,NGLOB
+
+! array with derivatives of Lagrange polynomials and precalculated products
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx
+ real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLY) :: hprime_yy
+ real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ) :: hprime_zz
+
+
+! arrays with mesh parameters per slice
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: ibool
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
+
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB) :: displ
+ real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: eps_trace_over_3_loc
+
+!local parameters
+ integer iglob
+ integer i,j,k,l
+
+ real(kind=CUSTOM_REAL) tempx1l,tempx2l,tempx3l
+ real(kind=CUSTOM_REAL) tempy1l,tempy2l,tempy3l
+ real(kind=CUSTOM_REAL) tempz1l,tempz2l,tempz3l
+
+ real(kind=CUSTOM_REAL) hp1,hp2,hp3
+
+ real(kind=CUSTOM_REAL) xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
+
+ real(kind=CUSTOM_REAL) duxdxl,duxdyl,duxdzl,duydxl,duydyl,duydzl,duzdxl,duzdyl,duzdzl
+
+ real(kind=CUSTOM_REAL) duxdxl_plus_duydyl,duxdxl_plus_duzdzl,duydyl_plus_duzdzl
+ real(kind=CUSTOM_REAL) duxdyl_plus_duydxl,duzdxl_plus_duxdzl,duzdyl_plus_duydzl
+
+
+
+
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+
+ tempx1l = 0._CUSTOM_REAL
+ tempx2l = 0._CUSTOM_REAL
+ tempx3l = 0._CUSTOM_REAL
+
+ tempy1l = 0._CUSTOM_REAL
+ tempy2l = 0._CUSTOM_REAL
+ tempy3l = 0._CUSTOM_REAL
+
+ tempz1l = 0._CUSTOM_REAL
+ tempz2l = 0._CUSTOM_REAL
+ tempz3l = 0._CUSTOM_REAL
+
+ do l=1,NGLLX
+ hp1 = hprime_xx(i,l)
+ iglob = ibool(l,j,k,ispec)
+ tempx1l = tempx1l + displ(1,iglob)*hp1
+ tempy1l = tempy1l + displ(2,iglob)*hp1
+ tempz1l = tempz1l + displ(3,iglob)*hp1
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLY
+ hp2 = hprime_yy(j,l)
+ iglob = ibool(i,l,k,ispec)
+ tempx2l = tempx2l + displ(1,iglob)*hp2
+ tempy2l = tempy2l + displ(2,iglob)*hp2
+ tempz2l = tempz2l + displ(3,iglob)*hp2
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLZ
+ hp3 = hprime_zz(k,l)
+ iglob = ibool(i,j,l,ispec)
+ tempx3l = tempx3l + displ(1,iglob)*hp3
+ tempy3l = tempy3l + displ(2,iglob)*hp3
+ tempz3l = tempz3l + displ(3,iglob)*hp3
+ enddo
+
+! get derivatives of ux, uy and uz with respect to x, y and z
+
+ xixl = xix(i,j,k,ispec)
+ xiyl = xiy(i,j,k,ispec)
+ xizl = xiz(i,j,k,ispec)
+ etaxl = etax(i,j,k,ispec)
+ etayl = etay(i,j,k,ispec)
+ etazl = etaz(i,j,k,ispec)
+ gammaxl = gammax(i,j,k,ispec)
+ gammayl = gammay(i,j,k,ispec)
+ gammazl = gammaz(i,j,k,ispec)
+
+! compute the jacobian
+ jacobianl = 1._CUSTOM_REAL / (xixl*(etayl*gammazl-etazl*gammayl) &
+ - xiyl*(etaxl*gammazl-etazl*gammaxl) &
+ + xizl*(etaxl*gammayl-etayl*gammaxl))
+
+ duxdxl = xixl*tempx1l + etaxl*tempx2l + gammaxl*tempx3l
+ duxdyl = xiyl*tempx1l + etayl*tempx2l + gammayl*tempx3l
+ duxdzl = xizl*tempx1l + etazl*tempx2l + gammazl*tempx3l
+
+ duydxl = xixl*tempy1l + etaxl*tempy2l + gammaxl*tempy3l
+ duydyl = xiyl*tempy1l + etayl*tempy2l + gammayl*tempy3l
+ duydzl = xizl*tempy1l + etazl*tempy2l + gammazl*tempy3l
+
+ duzdxl = xixl*tempz1l + etaxl*tempz2l + gammaxl*tempz3l
+ duzdyl = xiyl*tempz1l + etayl*tempz2l + gammayl*tempz3l
+ duzdzl = xizl*tempz1l + etazl*tempz2l + gammazl*tempz3l
+
+! precompute some sums to save CPU time
+ duxdxl_plus_duydyl = duxdxl + duydyl
+ duxdxl_plus_duzdzl = duxdxl + duzdzl
+ duydyl_plus_duzdzl = duydyl + duzdzl
+ duxdyl_plus_duydxl = duxdyl + duydxl
+ duzdxl_plus_duxdzl = duzdxl + duxdzl
+ duzdyl_plus_duydzl = duzdyl + duydzl
+
+! compute deviatoric strain
+
+ eps_trace_over_3_loc(i,j,k) = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ epsilondev_loc(1,i,j,k) = duxdxl - eps_trace_over_3_loc(i,j,k)
+ epsilondev_loc(2,i,j,k) = duydyl - eps_trace_over_3_loc(i,j,k)
+ epsilondev_loc(3,i,j,k) = 0.5 * duxdyl_plus_duydxl
+ epsilondev_loc(4,i,j,k) = 0.5 * duzdxl_plus_duxdzl
+ epsilondev_loc(5,i,j,k) = 0.5 * duzdyl_plus_duydzl
+
+ enddo ! NGLLX
+ enddo ! NGLLY
+ enddo ! NGLLZ
+
+ end subroutine compute_element_strain_undo_att_noDev
+!=====================================================================
+
+ subroutine compute_element_strain_att_noDev(ispec,nglob,nspec,displ,veloc,deltat,hprime_xx,hprime_yy,hprime_zz,ibool,&
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,epsilondev_loc_nplus1)
+
+ implicit none
+
+ include "constants.h"
+
+! include values created by the mesher
+! done for performance only using static allocation to allow for loop unrolling
+ include "OUTPUT_FILES/values_from_mesher.h"
+
+ integer ispec,NSPEC,NGLOB
+ real(kind=CUSTOM_REAL) deltat
+
+! array with derivatives of Lagrange polynomials and precalculated products
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx
+ real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLY) :: hprime_yy
+ real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ) :: hprime_zz
+
+
+! arrays with mesh parameters per slice
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: ibool
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
+
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB) :: displ,veloc
+
+ real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc_nplus1
+
+!local parameters
+ integer iglob
+ integer i,j,k,l
+ real(kind=CUSTOM_REAL) templ
+
+ real(kind=CUSTOM_REAL) tempx1l,tempx2l,tempx3l
+ real(kind=CUSTOM_REAL) tempy1l,tempy2l,tempy3l
+ real(kind=CUSTOM_REAL) tempz1l,tempz2l,tempz3l
+
+ real(kind=CUSTOM_REAL) hp1,hp2,hp3
+
+ real(kind=CUSTOM_REAL) xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
+
+ real(kind=CUSTOM_REAL) duxdxl,duxdyl,duxdzl,duydxl,duydyl,duydzl,duzdxl,duzdyl,duzdzl
+
+ real(kind=CUSTOM_REAL) duxdxl_plus_duydyl,duxdxl_plus_duzdzl,duydyl_plus_duzdzl
+ real(kind=CUSTOM_REAL) duxdyl_plus_duydxl,duzdxl_plus_duxdzl,duzdyl_plus_duydzl
+
+
+
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+
+ tempx1l = 0._CUSTOM_REAL
+ tempx2l = 0._CUSTOM_REAL
+ tempx3l = 0._CUSTOM_REAL
+
+ tempy1l = 0._CUSTOM_REAL
+ tempy2l = 0._CUSTOM_REAL
+ tempy3l = 0._CUSTOM_REAL
+
+ tempz1l = 0._CUSTOM_REAL
+ tempz2l = 0._CUSTOM_REAL
+ tempz3l = 0._CUSTOM_REAL
+
+ do l=1,NGLLX
+ hp1 = hprime_xx(i,l)
+ iglob = ibool(l,j,k,ispec)
+ tempx1l = tempx1l + (displ(1,iglob) + deltat * veloc(1,iglob))*hp1
+ tempy1l = tempy1l + (displ(2,iglob) + deltat * veloc(2,iglob))*hp1
+ tempz1l = tempz1l + (displ(3,iglob) + deltat * veloc(3,iglob))*hp1
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLY
+ hp2 = hprime_yy(j,l)
+ iglob = ibool(i,l,k,ispec)
+ tempx2l = tempx2l + (displ(1,iglob) + deltat * veloc(1,iglob))*hp2
+ tempy2l = tempy2l + (displ(2,iglob) + deltat * veloc(2,iglob))*hp2
+ tempz2l = tempz2l + (displ(3,iglob) + deltat * veloc(3,iglob))*hp2
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLZ
+ hp3 = hprime_zz(k,l)
+ iglob = ibool(i,j,l,ispec)
+ tempx3l = tempx3l + (displ(1,iglob) + deltat * veloc(1,iglob))*hp3
+ tempy3l = tempy3l + (displ(2,iglob) + deltat * veloc(2,iglob))*hp3
+ tempz3l = tempz3l + (displ(3,iglob) + deltat * veloc(3,iglob))*hp3
+ enddo
+
+! get derivatives of ux, uy and uz with respect to x, y and z
+
+ xixl = xix(i,j,k,ispec)
+ xiyl = xiy(i,j,k,ispec)
+ xizl = xiz(i,j,k,ispec)
+ etaxl = etax(i,j,k,ispec)
+ etayl = etay(i,j,k,ispec)
+ etazl = etaz(i,j,k,ispec)
+ gammaxl = gammax(i,j,k,ispec)
+ gammayl = gammay(i,j,k,ispec)
+ gammazl = gammaz(i,j,k,ispec)
+
+! compute the jacobian
+ jacobianl = 1._CUSTOM_REAL / (xixl*(etayl*gammazl-etazl*gammayl) &
+ - xiyl*(etaxl*gammazl-etazl*gammaxl) &
+ + xizl*(etaxl*gammayl-etayl*gammaxl))
+
+ duxdxl = xixl*tempx1l + etaxl*tempx2l + gammaxl*tempx3l
+ duxdyl = xiyl*tempx1l + etayl*tempx2l + gammayl*tempx3l
+ duxdzl = xizl*tempx1l + etazl*tempx2l + gammazl*tempx3l
+
+ duydxl = xixl*tempy1l + etaxl*tempy2l + gammaxl*tempy3l
+ duydyl = xiyl*tempy1l + etayl*tempy2l + gammayl*tempy3l
+ duydzl = xizl*tempy1l + etazl*tempy2l + gammazl*tempy3l
+
+ duzdxl = xixl*tempz1l + etaxl*tempz2l + gammaxl*tempz3l
+ duzdyl = xiyl*tempz1l + etayl*tempz2l + gammayl*tempz3l
+ duzdzl = xizl*tempz1l + etazl*tempz2l + gammazl*tempz3l
+
+! precompute some sums to save CPU time
+ duxdxl_plus_duydyl = duxdxl + duydyl
+ duxdxl_plus_duzdzl = duxdxl + duzdzl
+ duydyl_plus_duzdzl = duydyl + duzdzl
+ duxdyl_plus_duydxl = duxdyl + duydxl
+ duzdxl_plus_duxdzl = duzdxl + duxdzl
+ duzdyl_plus_duydzl = duzdyl + duydzl
+
+! compute deviatoric strain
+
+ templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ epsilondev_loc_nplus1(1,i,j,k) = duxdxl - templ
+ epsilondev_loc_nplus1(2,i,j,k) = duydyl - templ
+ epsilondev_loc_nplus1(3,i,j,k) = 0.5 * duxdyl_plus_duydxl
+ epsilondev_loc_nplus1(4,i,j,k) = 0.5 * duzdxl_plus_duxdzl
+ epsilondev_loc_nplus1(5,i,j,k) = 0.5 * duzdyl_plus_duydzl
+
+ enddo ! NGLLX
+ enddo ! NGLLY
+ enddo ! NGLLZ
+
+ end subroutine compute_element_strain_att_noDev
+
+
+
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/compute_forces_crust_mantle_Dev.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/compute_forces_crust_mantle_Dev.F90 2013-06-15 11:05:55 UTC (rev 22309)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/compute_forces_crust_mantle_Dev.F90 2013-06-15 15:41:35 UTC (rev 22310)
@@ -52,7 +52,8 @@
c23store,c24store,c25store,c26store,c33store,c34store,c35store, &
c36store,c44store,c45store,c46store,c55store,c56store,c66store, &
ibool,ispec_is_tiso, &
- R_memory,epsilondev,epsilon_trace_over_3,one_minus_sum_beta, &
+!ZN R_memory,epsilondev,epsilon_trace_over_3,one_minus_sum_beta, &
+ R_memory,one_minus_sum_beta,deltat,veloc_crust_mantle, &
alphaval,betaval,gammaval,factor_common,vx,vy,vz,vnspec)
! this routine is optimized for NGLLX = NGLLY = NGLLZ = 5 using the Deville et al. (2002) inlined matrix-matrix products
@@ -66,7 +67,8 @@
include "OUTPUT_FILES/values_from_mesher.h"
! displacement and acceleration
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CRUST_MANTLE) :: displ_crust_mantle,accel_crust_mantle
+ real(kind=CUSTOM_REAL) deltat
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CRUST_MANTLE) :: displ_crust_mantle,accel_crust_mantle,veloc_crust_mantle
! arrays with mesh parameters per slice
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: ibool
@@ -101,8 +103,8 @@
! memory variables R_ij are stored at the local rather than global level
! to allow for optimization of cache access by compiler
real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: R_memory
- real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: epsilondev
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: epsilon_trace_over_3
+!ZN real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: epsilondev
+!ZN real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: epsilon_trace_over_3
integer :: vx,vy,vz,vnspec
@@ -158,6 +160,7 @@
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NGLLZ) :: sum_terms
real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc
+ real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc_nplus1 !ZN
real(kind=CUSTOM_REAL) fac1,fac2,fac3
! for gravity
@@ -411,7 +414,8 @@
c23store,c24store,c25store,c26store,c33store,c34store,c35store, &
c36store,c44store,c45store,c46store,c55store,c56store,c66store, &
ibool, &
- R_memory,epsilon_trace_over_3, &
+!ZN R_memory,epsilon_trace_over_3, &
+ R_memory, &
one_minus_sum_beta,vx,vy,vz,vnspec, &
tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,rho_s_H)
@@ -425,7 +429,8 @@
wgll_cube, &
kappavstore,muvstore, &
ibool, &
- R_memory,epsilon_trace_over_3, &
+!ZN R_memory,epsilon_trace_over_3, &
+ R_memory, & !ZN
one_minus_sum_beta,vx,vy,vz,vnspec, &
tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,rho_s_H)
@@ -438,7 +443,8 @@
wgll_cube, &
kappavstore,kappahstore,muvstore,muhstore,eta_anisostore, &
ibool, &
- R_memory,epsilon_trace_over_3, &
+!ZN R_memory,epsilon_trace_over_3, &
+ R_memory, & !ZN
one_minus_sum_beta,vx,vy,vz,vnspec, &
tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,rho_s_H)
@@ -562,19 +568,23 @@
if(ATTENUATION_VAL .and. ( PARTIAL_PHYS_DISPERSION_ONLY .eqv. .false. ) ) then
+ call compute_element_strain_att_Dev(ispec,NGLOB_CRUST_MANTLE,NSPEC_CRUST_MANTLE,displ_crust_mantle,veloc_crust_mantle,&
+ deltat,ibool,hprime_xx,hprime_xxT,&
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,epsilondev_loc_nplus1)
! updates R_memory
call compute_element_att_memory_cr(ispec,R_memory, &
vx,vy,vz,vnspec,factor_common, &
alphaval,betaval,gammaval, &
c44store,muvstore, &
- epsilondev,epsilondev_loc)
+!ZN epsilondev,epsilondev_loc)
+ epsilondev_loc_nplus1,epsilondev_loc)
endif
! save deviatoric strain for Runge-Kutta scheme
- if(COMPUTE_AND_STORE_STRAIN) then
- epsilondev(:,:,:,:,ispec) = epsilondev_loc(:,:,:,:)
- endif
+!ZN if(COMPUTE_AND_STORE_STRAIN) then
+!ZN epsilondev(:,:,:,:,ispec) = epsilondev_loc(:,:,:,:)
+!ZN endif
! end ispec loop
enddo
!$OMP enddo
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/compute_forces_crust_mantle_noDev.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/compute_forces_crust_mantle_noDev.f90 2013-06-15 11:05:55 UTC (rev 22309)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/compute_forces_crust_mantle_noDev.f90 2013-06-15 15:41:35 UTC (rev 22310)
@@ -52,7 +52,8 @@
c23store,c24store,c25store,c26store,c33store,c34store,c35store, &
c36store,c44store,c45store,c46store,c55store,c56store,c66store, &
ibool,ispec_is_tiso, &
- R_memory,epsilondev,epsilon_trace_over_3,one_minus_sum_beta, &
+!ZN R_memory,epsilondev,epsilon_trace_over_3,one_minus_sum_beta, &
+ R_memory,one_minus_sum_beta,deltat,veloc_crust_mantle, & !ZN
alphaval,betaval,gammaval,factor_common,vx,vy,vz,vnspec)
implicit none
@@ -88,7 +89,7 @@
logical, dimension(NSPEC_CRUST_MANTLE) :: ispec_is_tiso
! displacement and acceleration
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CRUST_MANTLE) :: displ_crust_mantle,accel_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CRUST_MANTLE) :: displ_crust_mantle,accel_crust_mantle,veloc_crust_mantle !ZN
! memory variables for attenuation
! memory variables R_ij are stored at the local rather than global level
@@ -97,13 +98,13 @@
! variable sized array variables for one_minus_sum_beta and factor_common
integer vx, vy, vz, vnspec
- real(kind=CUSTOM_REAL) one_minus_sum_beta_use,minus_sum_beta
+ real(kind=CUSTOM_REAL) one_minus_sum_beta_use,minus_sum_beta,deltat !ZN
real(kind=CUSTOM_REAL), dimension(vx, vy, vz, vnspec) :: one_minus_sum_beta
! for attenuation
real(kind=CUSTOM_REAL) R_xx_val,R_yy_val
real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: R_memory
- real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: epsilondev
+!ZN real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: epsilondev
! [alpha,beta,gamma]val reduced to N_SLS and factor_common to N_SLS*NUM_NODES
real(kind=CUSTOM_REAL), dimension(N_SLS) :: alphaval,betaval,gammaval
@@ -111,7 +112,8 @@
real(kind=CUSTOM_REAL), dimension(NGLLX, NGLLY, NGLLZ) :: factor_common_c44_muv
real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc
- real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: epsilon_trace_over_3
+ real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc_nplus1
+!ZN real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: epsilon_trace_over_3
! arrays with mesh parameters per slice
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: ibool
@@ -147,6 +149,7 @@
integer ispec,iglob,ispec_strain
integer i,j,k,l
+ real(kind=CUSTOM_REAL) templ
! the 21 coefficients for an anisotropic medium in reduced notation
real(kind=CUSTOM_REAL) c11,c22,c33,c44,c55,c66,c12,c13,c23,c14,c24,c34,c15,c25,c35,c45,c16,c26,c36,c46,c56
@@ -380,9 +383,12 @@
else
ispec_strain = ispec
endif
- epsilon_trace_over_3(i,j,k,ispec_strain) = ONE_THIRD * (duxdxl + duydyl + duzdzl)
- epsilondev_loc(1,i,j,k) = duxdxl - epsilon_trace_over_3(i,j,k,ispec_strain)
- epsilondev_loc(2,i,j,k) = duydyl - epsilon_trace_over_3(i,j,k,ispec_strain)
+!ZN epsilon_trace_over_3(i,j,k,ispec_strain) = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+!ZN epsilondev_loc(1,i,j,k) = duxdxl - epsilon_trace_over_3(i,j,k,ispec_strain)
+!ZN epsilondev_loc(2,i,j,k) = duydyl - epsilon_trace_over_3(i,j,k,ispec_strain)
+ templ = ONE_THIRD * (duxdxl + duydyl + duzdzl) !ZN
+ epsilondev_loc(1,i,j,k) = duxdxl - templ !ZN
+ epsilondev_loc(2,i,j,k) = duydyl - templ !ZN
epsilondev_loc(3,i,j,k) = 0.5 * duxdyl_plus_duydxl
epsilondev_loc(4,i,j,k) = 0.5 * duzdxl_plus_duxdzl
epsilondev_loc(5,i,j,k) = 0.5 * duzdyl_plus_duydzl
@@ -911,6 +917,10 @@
if(ATTENUATION_VAL .and. ( PARTIAL_PHYS_DISPERSION_ONLY .eqv. .false. )) then
+ call compute_element_strain_att_noDev(ispec,NGLOB_CRUST_MANTLE,NSPEC_CRUST_MANTLE,displ_crust_mantle,veloc_crust_mantle,&
+ deltat,hprime_xx,hprime_yy,hprime_zz,ibool,&
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,epsilondev_loc_nplus1)
+
! use Runge-Kutta scheme to march in time
do i_SLS = 1,N_SLS
do i_memory = 1,5
@@ -930,24 +940,26 @@
R_memory(i_memory,i_SLS,:,:,:,ispec) = alphaval(i_SLS) * &
R_memory(i_memory,i_SLS,:,:,:,ispec) + &
factor_common_c44_muv * &
- (betaval(i_SLS) * epsilondev(i_memory,:,:,:,ispec) + &
- gammaval(i_SLS) * epsilondev_loc(i_memory,:,:,:))
+!ZN (betaval(i_SLS) * epsilondev(i_memory,:,:,:,ispec) + &
+!ZN gammaval(i_SLS) * epsilondev_loc(i_memory,:,:,:))
+ (betaval(i_SLS) * epsilondev_loc(i_memory,:,:,:) + &
+ gammaval(i_SLS) * epsilondev_loc_nplus1(i_memory,:,:,:))
enddo
enddo
endif
! save deviatoric strain for Runge-Kutta scheme
- if(COMPUTE_AND_STORE_STRAIN) then
- !epsilondev(:,:,:,:,ispec) = epsilondev_loc(:,:,:,:)
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- epsilondev(:,i,j,k,ispec) = epsilondev_loc(:,i,j,k)
- enddo
- enddo
- enddo
- endif
+!ZN if(COMPUTE_AND_STORE_STRAIN) then
+!ZN !epsilondev(:,:,:,:,ispec) = epsilondev_loc(:,:,:,:)
+!ZN do k=1,NGLLZ
+!ZN do j=1,NGLLY
+!ZN do i=1,NGLLX
+!ZN epsilondev(:,i,j,k,ispec) = epsilondev_loc(:,i,j,k)
+!ZN enddo
+!ZN enddo
+!ZN enddo
+!ZN endif
enddo ! spectral element loop NSPEC_CRUST_MANTLE
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/compute_forces_inner_core_Dev.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/compute_forces_inner_core_Dev.F90 2013-06-15 11:05:55 UTC (rev 22309)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/compute_forces_inner_core_Dev.F90 2013-06-15 15:41:35 UTC (rev 22310)
@@ -47,8 +47,9 @@
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
kappavstore,muvstore,ibool,idoubling, &
- c11store,c33store,c12store,c13store,c44store,R_memory,epsilondev,epsilon_trace_over_3,&
- one_minus_sum_beta,alphaval,betaval,gammaval,factor_common, &
+!ZN c11store,c33store,c12store,c13store,c44store,R_memory,epsilondev,epsilon_trace_over_3,&
+ c11store,c33store,c12store,c13store,c44store,R_memory,one_minus_sum_beta,deltat,veloc_inner_core,&
+ alphaval,betaval,gammaval,factor_common, &
vx,vy,vz,vnspec)
! this routine is optimized for NGLLX = NGLLY = NGLLZ = 5 using the Deville et al. (2002) inlined matrix-matrix products
@@ -57,12 +58,14 @@
include "constants.h"
+ real(kind=CUSTOM_REAL) deltat
+
! include values created by the mesher
! done for performance only using static allocation to allow for loop unrolling
include "OUTPUT_FILES/values_from_mesher.h"
! displacement and acceleration
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_INNER_CORE) :: displ_inner_core,accel_inner_core
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_INNER_CORE) :: displ_inner_core,accel_inner_core,veloc_inner_core
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: ibool
real(kind=CUSTOM_REAL), dimension(NGLOB_INNER_CORE) :: xstore,ystore,zstore
@@ -82,8 +85,8 @@
real(kind=CUSTOM_REAL), dimension(N_SLS) :: alphaval,betaval,gammaval
real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ATTENUATION) :: R_memory
- real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: epsilondev
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: epsilon_trace_over_3
+!ZN real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: epsilondev
+!ZN real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: epsilon_trace_over_3
! array with derivatives of Lagrange polynomials and precalculated products
double precision, dimension(NGLLX,NGLLY,NGLLZ) :: wgll_cube
@@ -143,6 +146,7 @@
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NGLLZ) :: sum_terms
real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc
+ real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc_nplus1 !ZN
real(kind=CUSTOM_REAL) xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
real(kind=CUSTOM_REAL) duxdxl,duxdyl,duxdzl,duydxl,duydyl,duydzl,duzdxl,duzdyl,duzdzl
@@ -152,7 +156,7 @@
real(kind=CUSTOM_REAL) sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz
- real(kind=CUSTOM_REAL) fac1,fac2,fac3,templ
+ real(kind=CUSTOM_REAL) fac1,fac2,fac3
real(kind=CUSTOM_REAL) lambdal,mul,lambdalplus2mul
real(kind=CUSTOM_REAL) kappal
@@ -170,7 +174,7 @@
real(kind=CUSTOM_REAL) sigma_yx,sigma_zx,sigma_zy
integer :: int_radius
- integer :: ispec,ispec_strain
+ integer :: ispec
integer :: i,j,k
integer :: iglob1
@@ -404,21 +408,21 @@
duzdxl_plus_duxdzl = duzdxl + duxdzl
duzdyl_plus_duydzl = duzdyl + duydzl
- ! compute deviatoric strain
- if (COMPUTE_AND_STORE_STRAIN) then
- if(NSPEC_INNER_CORE_STRAIN_ONLY == 1) then
- ispec_strain = 1
- else
- ispec_strain = ispec
- endif
- templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
- epsilon_trace_over_3(i,j,k,ispec_strain) = templ
- epsilondev_loc(1,i,j,k) = duxdxl - templ
- epsilondev_loc(2,i,j,k) = duydyl - templ
- epsilondev_loc(3,i,j,k) = 0.5 * duxdyl_plus_duydxl
- epsilondev_loc(4,i,j,k) = 0.5 * duzdxl_plus_duxdzl
- epsilondev_loc(5,i,j,k) = 0.5 * duzdyl_plus_duydzl
- endif
+!ZN ! compute deviatoric strain
+!ZN if (COMPUTE_AND_STORE_STRAIN) then
+!ZN if(NSPEC_INNER_CORE_STRAIN_ONLY == 1) then
+!ZN ispec_strain = 1
+!ZN else
+!ZN ispec_strain = ispec
+!ZN endif
+!ZN templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+!ZN epsilon_trace_over_3(i,j,k,ispec_strain) = templ
+!ZN epsilondev_loc(1,i,j,k) = duxdxl - templ
+!ZN epsilondev_loc(2,i,j,k) = duydyl - templ
+!ZN epsilondev_loc(3,i,j,k) = 0.5 * duxdyl_plus_duydxl
+!ZN epsilondev_loc(4,i,j,k) = 0.5 * duzdxl_plus_duxdzl
+!ZN epsilondev_loc(5,i,j,k) = 0.5 * duzdyl_plus_duydzl
+!ZN endif
if(ATTENUATION_VAL) then
minus_sum_beta = one_minus_sum_beta(i,j,k,ispec) - 1.0
@@ -747,20 +751,25 @@
! we get Q_\alpha = (9 / 4) * Q_\mu = 2.25 * Q_\mu
if(ATTENUATION_VAL .and. ( PARTIAL_PHYS_DISPERSION_ONLY .eqv. .false. ) ) then
+ call compute_element_strain_att_Dev(ispec,NGLOB_INNER_CORE,NSPEC_INNER_CORE,displ_inner_core,&
+ veloc_inner_core,deltat,ibool,hprime_xx,hprime_xxT,&
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,epsilondev_loc_nplus1)
! updates R_memory
call compute_element_att_memory_ic(ispec,R_memory, &
vx,vy,vz,vnspec,factor_common, &
alphaval,betaval,gammaval, &
muvstore, &
- epsilondev,epsilondev_loc)
+!ZN epsilondev,epsilondev_loc)
+ epsilondev_loc_nplus1,epsilondev_loc) !ZN
- endif
- ! save deviatoric strain for Runge-Kutta scheme
- if(COMPUTE_AND_STORE_STRAIN) then
- epsilondev(:,:,:,:,ispec) = epsilondev_loc(:,:,:,:)
endif
+!ZN ! save deviatoric strain for Runge-Kutta scheme
+!ZN if(COMPUTE_AND_STORE_STRAIN) then
+!ZN epsilondev(:,:,:,:,ispec) = epsilondev_loc(:,:,:,:)
+!ZN endif
+
endif ! end test to exclude fictitious elements in central cube
enddo ! spectral element loop
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/compute_forces_inner_core_noDev.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/compute_forces_inner_core_noDev.f90 2013-06-15 11:05:55 UTC (rev 22309)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/compute_forces_inner_core_noDev.f90 2013-06-15 15:41:35 UTC (rev 22310)
@@ -48,8 +48,9 @@
hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
kappavstore,muvstore,ibool,idoubling, &
- c11store,c33store,c12store,c13store,c44store,R_memory,epsilondev,epsilon_trace_over_3,&
- one_minus_sum_beta,alphaval,betaval,gammaval,factor_common, &
+!ZN c11store,c33store,c12store,c13store,c44store,R_memory,epsilondev,epsilon_trace_over_3,&
+ c11store,c33store,c12store,c13store,c44store,R_memory,one_minus_sum_beta,deltat,veloc_inner_core,&
+ alphaval,betaval,gammaval,factor_common, &
vx,vy,vz,vnspec)
implicit none
@@ -60,8 +61,10 @@
! done for performance only using static allocation to allow for loop unrolling
include "OUTPUT_FILES/values_from_mesher.h"
+ real(kind=CUSTOM_REAL) deltat
+
! displacement and acceleration
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_INNER_CORE) :: displ_inner_core,accel_inner_core
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_INNER_CORE) :: displ_inner_core,accel_inner_core,veloc_inner_core
! for attenuation
! memory variables R_ij are stored at the local rather than global level
@@ -79,9 +82,9 @@
real(kind=CUSTOM_REAL), dimension(NGLLX, NGLLY, NGLLZ) :: factor_common_use
real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ATTENUATION) :: R_memory
- real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: epsilondev
- real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: epsilon_trace_over_3
+!ZN real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: epsilondev
+ real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc,epsilondev_loc_nplus1 !ZN
+!ZN real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: epsilon_trace_over_3
! array with the local to global mapping per slice
integer, dimension(NSPEC_INNER_CORE) :: idoubling
@@ -207,6 +210,7 @@
! local to global mapping
integer NSPEC2D_BOTTOM_INNER_CORE
integer, dimension(NSPEC2D_BOTTOM_INNER_CORE) :: ibelm_bottom_inner_core
+ real(kind=CUSTOM_REAL) templ
! ****************************************************
! big loop over all spectral elements in the solid
@@ -337,9 +341,12 @@
else
ispec_strain = ispec
endif
- epsilon_trace_over_3(i,j,k,ispec_strain) = ONE_THIRD * (duxdxl + duydyl + duzdzl)
- epsilondev_loc(1,i,j,k) = duxdxl - epsilon_trace_over_3(i,j,k,ispec_strain)
- epsilondev_loc(2,i,j,k) = duydyl - epsilon_trace_over_3(i,j,k,ispec_strain)
+!ZN epsilon_trace_over_3(i,j,k,ispec_strain) = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+!ZN epsilondev_loc(1,i,j,k) = duxdxl - epsilon_trace_over_3(i,j,k,ispec_strain)
+!ZN epsilondev_loc(2,i,j,k) = duydyl - epsilon_trace_over_3(i,j,k,ispec_strain)
+ templ = ONE_THIRD * (duxdxl + duydyl + duzdzl) !ZN
+ epsilondev_loc(1,i,j,k) = duxdxl - templ !ZN
+ epsilondev_loc(2,i,j,k) = duydyl - templ !ZN
epsilondev_loc(3,i,j,k) = 0.5 * duxdyl_plus_duydxl
epsilondev_loc(4,i,j,k) = 0.5 * duzdxl_plus_duxdzl
epsilondev_loc(5,i,j,k) = 0.5 * duzdyl_plus_duydzl
@@ -644,6 +651,10 @@
if(ATTENUATION_VAL .and. ( PARTIAL_PHYS_DISPERSION_ONLY .eqv. .false. )) then
+ call compute_element_strain_att_noDev(ispec,NGLOB_INNER_CORE,NSPEC_INNER_CORE,displ_inner_core,&
+ veloc_inner_core,deltat,hprime_xx,hprime_yy,hprime_zz,ibool,&
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,epsilondev_loc_nplus1)
+
do i_SLS = 1,N_SLS
factor_common_use = factor_common(i_SLS,:,:,:,ispec)
do i_memory = 1,5
@@ -652,25 +663,25 @@
R_memory(i_memory,i_SLS,:,:,:,ispec) + muvstore(:,:,:,ispec) * &
factor_common_use * &
(betaval(i_SLS) * &
- epsilondev(i_memory,:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(i_memory,:,:,:))
+!ZN epsilondev(i_memory,:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(i_memory,:,:,:))
+ epsilondev_loc_nplus1(i_memory,:,:,:) + gammaval(i_SLS) * epsilondev_loc(i_memory,:,:,:)) !ZN
enddo
enddo
endif
- if (COMPUTE_AND_STORE_STRAIN) then
-! save deviatoric strain for Runge-Kutta scheme
- !epsilondev(:,:,:,:,ispec) = epsilondev_loc(:,:,:,:)
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- epsilondev(:,i,j,k,ispec) = epsilondev_loc(:,i,j,k)
- enddo
- enddo
- enddo
+!ZN if (COMPUTE_AND_STORE_STRAIN) then
+!ZN! save deviatoric strain for Runge-Kutta scheme
+!ZN !epsilondev(:,:,:,:,ispec) = epsilondev_loc(:,:,:,:)
+!ZN do k=1,NGLLZ
+!ZN do j=1,NGLLY
+!ZN do i=1,NGLLX
+!ZN epsilondev(:,i,j,k,ispec) = epsilondev_loc(:,i,j,k)
+!ZN enddo
+!ZN enddo
+!ZN enddo
+!ZN endif
- endif
-
endif ! end test to exclude fictitious elements in central cube
enddo ! spectral element loop
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/compute_kernels.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/compute_kernels.f90 2013-06-15 11:05:55 UTC (rev 22309)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/compute_kernels.f90 2013-06-15 15:41:35 UTC (rev 22310)
@@ -30,9 +30,10 @@
rho_kl_crust_mantle,beta_kl_crust_mantle, &
alpha_kl_crust_mantle,cijkl_kl_crust_mantle, &
accel_crust_mantle,b_displ_crust_mantle, &
- epsilondev_crust_mantle,b_epsilondev_crust_mantle, &
- eps_trace_over_3_crust_mantle,b_eps_trace_over_3_crust_mantle, &
- deltat)
+!ZN epsilondev_crust_mantle,b_epsilondev_crust_mantle, &
+!ZN eps_trace_over_3_crust_mantle,b_eps_trace_over_3_crust_mantle, &
+ deltat,displ_crust_mantle,hprime_xx,hprime_xxT,&
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz)
implicit none
@@ -49,30 +50,48 @@
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CRUST_MANTLE) :: &
- accel_crust_mantle
+ accel_crust_mantle,displ_crust_mantle
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CRUST_MANTLE_ADJOINT) :: &
b_displ_crust_mantle
- real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STR_OR_ATT) :: &
- epsilondev_crust_mantle
- real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT) :: &
- b_epsilondev_crust_mantle
+!ZN real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STR_OR_ATT) :: &
+!ZN epsilondev_crust_mantle
+!ZN real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT) :: &
+!ZN b_epsilondev_crust_mantle
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STRAIN_ONLY) :: &
- eps_trace_over_3_crust_mantle
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT) :: &
- b_eps_trace_over_3_crust_mantle
+!ZN real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STRAIN_ONLY) :: &
+!ZN eps_trace_over_3_crust_mantle
+!ZN real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT) :: &
+!ZN b_eps_trace_over_3_crust_mantle
real(kind=CUSTOM_REAL) deltat
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx,hprime_xxT
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
! local parameters
real(kind=CUSTOM_REAL),dimension(21) :: prod !, cijkl_kl_local
real(kind=CUSTOM_REAL), dimension(5) :: epsilondev_loc
real(kind=CUSTOM_REAL), dimension(5) :: b_epsilondev_loc
+ real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc_matrix,b_epsilondev_loc_matrix
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: eps_trace_over_3_loc_matrix,&
+ b_eps_trace_over_3_loc_matrix
integer :: i,j,k,ispec,iglob
! crust_mantle
do ispec = 1, NSPEC_CRUST_MANTLE
+
+ call compute_element_strain_undo_att_Dev(ispec,NGLOB_CRUST_MANTLE,NSPEC_CRUST_MANTLE,&
+ displ_CRUST_MANTLE,ibool_crust_mantle,hprime_xx,hprime_xxT,&
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,&
+ epsilondev_loc_matrix,eps_trace_over_3_loc_matrix)
+
+ call compute_element_strain_undo_att_Dev(ispec,NGLOB_CRUST_MANTLE,NSPEC_CRUST_MANTLE,&
+ b_displ_CRUST_MANTLE,ibool_crust_mantle,hprime_xx,hprime_xxT,&
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,&
+ b_epsilondev_loc_matrix,b_eps_trace_over_3_loc_matrix)
+
+
do k = 1, NGLLZ
do j = 1, NGLLY
do i = 1, NGLLX
@@ -95,14 +114,19 @@
+ accel_crust_mantle(2,iglob) * b_displ_crust_mantle(2,iglob) &
+ accel_crust_mantle(3,iglob) * b_displ_crust_mantle(3,iglob) )
- epsilondev_loc(:) = epsilondev_crust_mantle(:,i,j,k,ispec)
- b_epsilondev_loc(:) = b_epsilondev_crust_mantle(:,i,j,k,ispec)
+!ZN epsilondev_loc(:) = epsilondev_crust_mantle(:,i,j,k,ispec)
+!ZN b_epsilondev_loc(:) = b_epsilondev_crust_mantle(:,i,j,k,ispec)
+ epsilondev_loc(:) = epsilondev_loc_matrix(:,i,j,k)
+ b_epsilondev_loc(:) = b_epsilondev_loc_matrix(:,i,j,k)
! For anisotropic kernels
if (ANISOTROPIC_KL) then
- call compute_strain_product(prod,eps_trace_over_3_crust_mantle(i,j,k,ispec),epsilondev_loc, &
- b_eps_trace_over_3_crust_mantle(i,j,k,ispec),b_epsilondev_loc)
+!ZN call compute_strain_product(prod,eps_trace_over_3_crust_mantle(i,j,k,ispec),epsilondev_loc, &
+!ZN b_eps_trace_over_3_crust_mantle(i,j,k,ispec),b_epsilondev_loc)
+
+ call compute_strain_product(prod,eps_trace_over_3_loc_matrix(i,j,k),epsilondev_loc, &
+ b_eps_trace_over_3_loc_matrix(i,j,k),b_epsilondev_loc)
cijkl_kl_crust_mantle(:,i,j,k,ispec) = cijkl_kl_crust_mantle(:,i,j,k,ispec) + deltat * prod(:)
else
@@ -118,9 +142,12 @@
! kernel for bulk modulus, see e.g. Tromp et al. (2005), equation (18)
! note: multiplication with kappa(x) will be done after the time loop
+!ZN alpha_kl_crust_mantle(i,j,k,ispec) = alpha_kl_crust_mantle(i,j,k,ispec) &
+!ZN + deltat * (9 * eps_trace_over_3_crust_mantle(i,j,k,ispec) &
+!ZN * b_eps_trace_over_3_crust_mantle(i,j,k,ispec))
alpha_kl_crust_mantle(i,j,k,ispec) = alpha_kl_crust_mantle(i,j,k,ispec) &
- + deltat * (9 * eps_trace_over_3_crust_mantle(i,j,k,ispec) &
- * b_eps_trace_over_3_crust_mantle(i,j,k,ispec))
+ + deltat * (9 * eps_trace_over_3_loc_matrix(i,j,k) &
+ * b_eps_trace_over_3_loc_matrix(i,j,k))
endif
@@ -446,9 +473,10 @@
rho_kl_inner_core,beta_kl_inner_core, &
alpha_kl_inner_core, &
accel_inner_core,b_displ_inner_core, &
- epsilondev_inner_core,b_epsilondev_inner_core, &
- eps_trace_over_3_inner_core,b_eps_trace_over_3_inner_core, &
- deltat)
+!ZN epsilondev_inner_core,b_epsilondev_inner_core, &
+!ZN eps_trace_over_3_inner_core,b_eps_trace_over_3_inner_core, &
+ deltat,displ_inner_core,hprime_xx,hprime_xxT,&
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz)
implicit none
@@ -462,31 +490,48 @@
rho_kl_inner_core, beta_kl_inner_core, alpha_kl_inner_core
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_INNER_CORE) :: &
- accel_inner_core
+ accel_inner_core,displ_inner_core
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_INNER_CORE_ADJOINT) :: &
b_displ_inner_core
- real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_STR_OR_ATT) :: &
- epsilondev_inner_core
- real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ADJOINT) :: &
- b_epsilondev_inner_core
+!ZN real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_STR_OR_ATT) :: &
+!ZN epsilondev_inner_core
+!ZN real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ADJOINT) :: &
+!ZN b_epsilondev_inner_core
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_STRAIN_ONLY) :: &
- eps_trace_over_3_inner_core
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ADJOINT) :: &
- b_eps_trace_over_3_inner_core
+!ZN real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_STRAIN_ONLY) :: &
+!ZN eps_trace_over_3_inner_core
+!ZN real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ADJOINT) :: &
+!ZN b_eps_trace_over_3_inner_core
real(kind=CUSTOM_REAL) deltat
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx,hprime_xxT
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
! local parameters
real(kind=CUSTOM_REAL), dimension(5) :: b_epsilondev_loc
real(kind=CUSTOM_REAL), dimension(5) :: epsilondev_loc
+ real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc_matrix,b_epsilondev_loc_matrix
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: eps_trace_over_3_loc_matrix,&
+ b_eps_trace_over_3_loc_matrix
integer :: i,j,k,ispec,iglob
! inner_core
do ispec = 1, NSPEC_INNER_CORE
+
+ call compute_element_strain_undo_att_Dev(ispec,NGLOB_inner_core,NSPEC_inner_core,&
+ displ_inner_core,ibool_inner_core,hprime_xx,hprime_xxT,&
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,&
+ epsilondev_loc_matrix,eps_trace_over_3_loc_matrix)
+
+ call compute_element_strain_undo_att_Dev(ispec,NGLOB_inner_core,NSPEC_inner_core,&
+ b_displ_inner_core,ibool_inner_core,hprime_xx,hprime_xxT,&
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,&
+ b_epsilondev_loc_matrix,b_eps_trace_over_3_loc_matrix)
+
do k = 1, NGLLZ
do j = 1, NGLLY
do i = 1, NGLLX
@@ -497,16 +542,22 @@
+ accel_inner_core(2,iglob) * b_displ_inner_core(2,iglob) &
+ accel_inner_core(3,iglob) * b_displ_inner_core(3,iglob) )
- epsilondev_loc(:) = epsilondev_inner_core(:,i,j,k,ispec)
- b_epsilondev_loc(:) = b_epsilondev_inner_core(:,i,j,k,ispec)
+!ZN epsilondev_loc(:) = epsilondev_inner_core(:,i,j,k,ispec)
+!ZN b_epsilondev_loc(:) = b_epsilondev_inner_core(:,i,j,k,ispec)
+
+ epsilondev_loc(:) = epsilondev_loc_matrix(:,i,j,k)
+ b_epsilondev_loc(:) = b_epsilondev_loc_matrix(:,i,j,k)
+
beta_kl_inner_core(i,j,k,ispec) = beta_kl_inner_core(i,j,k,ispec) &
+ deltat * (epsilondev_loc(1)*b_epsilondev_loc(1) + epsilondev_loc(2)*b_epsilondev_loc(2) &
+ (epsilondev_loc(1)+epsilondev_loc(2)) * (b_epsilondev_loc(1)+b_epsilondev_loc(2)) &
+ 2 * (epsilondev_loc(3)*b_epsilondev_loc(3) + epsilondev_loc(4)*b_epsilondev_loc(4) &
+ epsilondev_loc(5)*b_epsilondev_loc(5)) )
+!ZN alpha_kl_inner_core(i,j,k,ispec) = alpha_kl_inner_core(i,j,k,ispec) &
+!ZN + deltat * (9 * eps_trace_over_3_inner_core(i,j,k,ispec) * b_eps_trace_over_3_inner_core(i,j,k,ispec))
alpha_kl_inner_core(i,j,k,ispec) = alpha_kl_inner_core(i,j,k,ispec) &
- + deltat * (9 * eps_trace_over_3_inner_core(i,j,k,ispec) * b_eps_trace_over_3_inner_core(i,j,k,ispec))
+ + deltat * (9 * eps_trace_over_3_loc_matrix(i,j,k) * b_eps_trace_over_3_loc_matrix(i,j,k))
enddo
enddo
enddo
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/compute_seismograms.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/compute_seismograms.f90 2013-06-15 11:05:55 UTC (rev 22309)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/compute_seismograms.f90 2013-06-15 15:41:35 UTC (rev 22310)
@@ -191,7 +191,7 @@
!
subroutine compute_seismograms_adjoint(NSOURCES,nrec_local,displ_crust_mantle, &
- eps_trace_over_3_crust_mantle,epsilondev_crust_mantle, &
+!ZN eps_trace_over_3_crust_mantle,epsilondev_crust_mantle, &
nu_source,Mxx,Myy,Mzz,Mxy,Mxz,Myz, &
hxir_store,hetar_store,hgammar_store, &
hpxir_store,hpetar_store,hpgammar_store, &
@@ -213,10 +213,10 @@
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CRUST_MANTLE) :: &
displ_crust_mantle
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STRAIN_ONLY) :: &
- eps_trace_over_3_crust_mantle
- real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STR_OR_ATT) :: &
- epsilondev_crust_mantle
+!ZN real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STRAIN_ONLY) :: &
+!ZN eps_trace_over_3_crust_mantle
+!ZN real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STR_OR_ATT) :: &
+!ZN epsilondev_crust_mantle
double precision, dimension(NDIM,NDIM,NSOURCES) :: nu_source
double precision, dimension(NSOURCES) :: Mxx,Myy,Mzz,Mxy,Mxz,Myz
@@ -266,6 +266,9 @@
double precision, external :: comp_source_time_function
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: eps_trace_over_3_loc_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc_crust_mantle
+
do irec_local = 1,nrec_local
! get global number of that receiver
@@ -284,6 +287,13 @@
dxz = ZERO
dyz = ZERO
+ call compute_element_strain_undo_att_noDev(ispec_selected_source(irec),NGLOB_CRUST_MANTLE,NSPEC_CRUST_MANTLE,&
+ displ_crust_mantle,hprime_xx,hprime_yy,hprime_zz,ibool_crust_mantle,&
+ xix_crust_mantle,xiy_crust_mantle,xiz_crust_mantle, &
+ etax_crust_mantle,etay_crust_mantle,etaz_crust_mantle, &
+ gammax_crust_mantle,gammay_crust_mantle,gammaz_crust_mantle,&
+ epsilondev_loc_crust_mantle,eps_trace_over_3_loc_crust_mantle)
+
do k = 1,NGLLZ
do j = 1,NGLLY
do i = 1,NGLLX
@@ -296,13 +306,20 @@
uyd = uyd + dble(displ_crust_mantle(2,iglob))*hlagrange
uzd = uzd + dble(displ_crust_mantle(3,iglob))*hlagrange
- eps_trace = eps_trace + dble(eps_trace_over_3_crust_mantle(i,j,k,ispec_selected_source(irec)))*hlagrange
- dxx = dxx + dble(epsilondev_crust_mantle(1,i,j,k,ispec_selected_source(irec)))*hlagrange
- dyy = dyy + dble(epsilondev_crust_mantle(2,i,j,k,ispec_selected_source(irec)))*hlagrange
- dxy = dxy + dble(epsilondev_crust_mantle(3,i,j,k,ispec_selected_source(irec)))*hlagrange
- dxz = dxz + dble(epsilondev_crust_mantle(4,i,j,k,ispec_selected_source(irec)))*hlagrange
- dyz = dyz + dble(epsilondev_crust_mantle(5,i,j,k,ispec_selected_source(irec)))*hlagrange
+!ZN eps_trace = eps_trace + dble(eps_trace_over_3_crust_mantle(i,j,k,ispec_selected_source(irec)))*hlagrange
+!ZN dxx = dxx + dble(epsilondev_crust_mantle(1,i,j,k,ispec_selected_source(irec)))*hlagrange
+!ZN dyy = dyy + dble(epsilondev_crust_mantle(2,i,j,k,ispec_selected_source(irec)))*hlagrange
+!ZN dxy = dxy + dble(epsilondev_crust_mantle(3,i,j,k,ispec_selected_source(irec)))*hlagrange
+!ZN dxz = dxz + dble(epsilondev_crust_mantle(4,i,j,k,ispec_selected_source(irec)))*hlagrange
+!ZN dyz = dyz + dble(epsilondev_crust_mantle(5,i,j,k,ispec_selected_source(irec)))*hlagrange
+ eps_trace = eps_trace + dble(eps_trace_over_3_loc_crust_mantle(i,j,k))*hlagrange
+ dxx = dxx + dble(epsilondev_loc_crust_mantle(1,i,j,k))*hlagrange
+ dyy = dyy + dble(epsilondev_loc_crust_mantle(2,i,j,k))*hlagrange
+ dxy = dxy + dble(epsilondev_loc_crust_mantle(3,i,j,k))*hlagrange
+ dxz = dxz + dble(epsilondev_loc_crust_mantle(4,i,j,k))*hlagrange
+ dyz = dyz + dble(epsilondev_loc_crust_mantle(5,i,j,k))*hlagrange
+
displ_s(:,i,j,k) = displ_crust_mantle(:,iglob)
enddo
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/read_forward_arrays.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/read_forward_arrays.f90 2013-06-15 11:05:55 UTC (rev 22309)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/read_forward_arrays.f90 2013-06-15 15:41:35 UTC (rev 22310)
@@ -32,13 +32,13 @@
displ_inner_core,veloc_inner_core,accel_inner_core, &
displ_outer_core,veloc_outer_core,accel_outer_core, &
R_memory_crust_mantle,R_memory_inner_core, &
- epsilondev_crust_mantle,epsilondev_inner_core, &
+!ZN epsilondev_crust_mantle,epsilondev_inner_core, &
A_array_rotation,B_array_rotation, &
b_displ_crust_mantle,b_veloc_crust_mantle,b_accel_crust_mantle, &
b_displ_inner_core,b_veloc_inner_core,b_accel_inner_core, &
b_displ_outer_core,b_veloc_outer_core,b_accel_outer_core, &
b_R_memory_crust_mantle,b_R_memory_inner_core, &
- b_epsilondev_crust_mantle,b_epsilondev_inner_core, &
+!ZN b_epsilondev_crust_mantle,b_epsilondev_inner_core, &
b_A_array_rotation,b_B_array_rotation,LOCAL_PATH)
! reads in saved wavefields
@@ -63,12 +63,12 @@
real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: &
R_memory_crust_mantle
- real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STR_OR_ATT) :: &
- epsilondev_crust_mantle
+!ZN real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STR_OR_ATT) :: &
+!ZN epsilondev_crust_mantle
real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ATTENUATION) :: &
R_memory_inner_core
- real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_STR_OR_ATT) :: &
- epsilondev_inner_core
+!ZN real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_STR_OR_ATT) :: &
+!ZN epsilondev_inner_core
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_OUTER_CORE_ROTATION) :: &
A_array_rotation,B_array_rotation
@@ -82,13 +82,13 @@
real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STR_AND_ATT) :: &
b_R_memory_crust_mantle
- real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT) :: &
- b_epsilondev_crust_mantle
+!ZN real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT) :: &
+!ZN b_epsilondev_crust_mantle
real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_STR_AND_ATT) :: &
b_R_memory_inner_core
- real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ADJOINT) :: &
- b_epsilondev_inner_core
+!ZN real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ADJOINT) :: &
+!ZN b_epsilondev_inner_core
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_OUTER_CORE_ROT_ADJOINT) :: &
b_A_array_rotation,b_B_array_rotation
@@ -125,8 +125,8 @@
read(55) displ_outer_core
read(55) veloc_outer_core
read(55) accel_outer_core
- read(55) epsilondev_crust_mantle
- read(55) epsilondev_inner_core
+!ZN read(55) epsilondev_crust_mantle
+!ZN read(55) epsilondev_inner_core
read(55) A_array_rotation
read(55) B_array_rotation
read(55) R_memory_crust_mantle
@@ -145,8 +145,8 @@
b_displ_outer_core = 0._CUSTOM_REAL
b_veloc_outer_core = 0._CUSTOM_REAL
b_accel_outer_core = 0._CUSTOM_REAL
- b_epsilondev_crust_mantle = 0._CUSTOM_REAL
- b_epsilondev_inner_core = 0._CUSTOM_REAL
+!ZN b_epsilondev_crust_mantle = 0._CUSTOM_REAL
+!ZN b_epsilondev_inner_core = 0._CUSTOM_REAL
if (ROTATION_VAL) then
b_A_array_rotation = 0._CUSTOM_REAL
b_B_array_rotation = 0._CUSTOM_REAL
@@ -168,7 +168,7 @@
b_displ_inner_core,b_veloc_inner_core,b_accel_inner_core, &
b_displ_outer_core,b_veloc_outer_core,b_accel_outer_core, &
b_R_memory_crust_mantle,b_R_memory_inner_core, &
- b_epsilondev_crust_mantle,b_epsilondev_inner_core, &
+!ZN b_epsilondev_crust_mantle,b_epsilondev_inner_core, &
b_A_array_rotation,b_B_array_rotation,LOCAL_PATH)
! reads in saved wavefields
@@ -190,13 +190,13 @@
real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STR_AND_ATT) :: &
b_R_memory_crust_mantle
- real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT) :: &
- b_epsilondev_crust_mantle
+!ZN real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT) :: &
+!ZN b_epsilondev_crust_mantle
real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_STR_AND_ATT) :: &
b_R_memory_inner_core
- real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ADJOINT) :: &
- b_epsilondev_inner_core
+!ZN real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ADJOINT) :: &
+!ZN b_epsilondev_inner_core
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_OUTER_CORE_ROT_ADJOINT) :: &
b_A_array_rotation,b_B_array_rotation
@@ -217,8 +217,8 @@
read(55) b_displ_outer_core
read(55) b_veloc_outer_core
read(55) b_accel_outer_core
- read(55) b_epsilondev_crust_mantle
- read(55) b_epsilondev_inner_core
+!ZN read(55) b_epsilondev_crust_mantle
+!ZN read(55) b_epsilondev_inner_core
if (ROTATION_VAL) then
read(55) b_A_array_rotation
read(55) b_B_array_rotation
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/save_forward_arrays.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/save_forward_arrays.f90 2013-06-15 11:05:55 UTC (rev 22309)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/save_forward_arrays.f90 2013-06-15 15:41:35 UTC (rev 22310)
@@ -31,7 +31,7 @@
displ_inner_core,veloc_inner_core,accel_inner_core, &
displ_outer_core,veloc_outer_core,accel_outer_core, &
R_memory_crust_mantle,R_memory_inner_core, &
- epsilondev_crust_mantle,epsilondev_inner_core, &
+!ZN epsilondev_crust_mantle,epsilondev_inner_core, &
A_array_rotation,B_array_rotation, &
LOCAL_PATH)
@@ -55,12 +55,12 @@
real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: &
R_memory_crust_mantle
- real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STR_OR_ATT) :: &
- epsilondev_crust_mantle
+!ZN real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STR_OR_ATT) :: &
+!ZN epsilondev_crust_mantle
real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ATTENUATION) :: &
R_memory_inner_core
- real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_STR_OR_ATT) :: &
- epsilondev_inner_core
+!ZN real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_STR_OR_ATT) :: &
+!ZN epsilondev_inner_core
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_OUTER_CORE_ROTATION) :: &
A_array_rotation,B_array_rotation
@@ -84,8 +84,8 @@
write(55) displ_outer_core
write(55) veloc_outer_core
write(55) accel_outer_core
- write(55) epsilondev_crust_mantle
- write(55) epsilondev_inner_core
+!ZN write(55) epsilondev_crust_mantle
+!ZN write(55) epsilondev_inner_core
write(55) A_array_rotation
write(55) B_array_rotation
write(55) R_memory_crust_mantle
@@ -106,8 +106,8 @@
write(55) displ_outer_core
write(55) veloc_outer_core
write(55) accel_outer_core
- write(55) epsilondev_crust_mantle
- write(55) epsilondev_inner_core
+!ZN write(55) epsilondev_crust_mantle
+!ZN write(55) epsilondev_inner_core
if (ROTATION_VAL) then
write(55) A_array_rotation
write(55) B_array_rotation
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/specfem3D.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/specfem3D.F90 2013-06-15 11:05:55 UTC (rev 22309)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/specfem3D.F90 2013-06-15 15:41:35 UTC (rev 22310)
@@ -395,23 +395,25 @@
real(kind=CUSTOM_REAL), dimension(N_SLS,ATT1,ATT2,ATT3,ATT5) :: factor_common_inner_core
real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: R_memory_crust_mantle
- real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STR_OR_ATT) :: epsilondev_crust_mantle
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STRAIN_ONLY) :: eps_trace_over_3_crust_mantle
+!ZN real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STR_OR_ATT) :: epsilondev_crust_mantle
+!ZN real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STRAIN_ONLY) :: eps_trace_over_3_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: eps_trace_over_3_loc_crust_mantle
real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ATTENUATION) :: R_memory_inner_core
- real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_STR_OR_ATT) :: epsilondev_inner_core
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_STRAIN_ONLY) :: eps_trace_over_3_inner_core
+!ZN real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_STR_OR_ATT) :: epsilondev_inner_core
+!ZN real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_STRAIN_ONLY) :: eps_trace_over_3_inner_core
! ADJOINT
real(kind=CUSTOM_REAL), dimension(N_SLS) :: b_alphaval, b_betaval, b_gammaval
real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STR_AND_ATT) :: b_R_memory_crust_mantle
- real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT) :: b_epsilondev_crust_mantle
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT) :: b_eps_trace_over_3_crust_mantle
+!ZN real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT) :: b_epsilondev_crust_mantle
+!ZN real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT) :: b_eps_trace_over_3_crust_mantle
real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_STR_AND_ATT) :: b_R_memory_inner_core
- real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ADJOINT) :: b_epsilondev_inner_core
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ADJOINT) :: b_eps_trace_over_3_inner_core
+!ZN real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ADJOINT) :: b_epsilondev_inner_core
+!ZN real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ADJOINT) :: b_eps_trace_over_3_inner_core
! for matching with central cube in inner core
integer, dimension(:), allocatable :: sender_from_slices_to_cube
@@ -2032,16 +2034,16 @@
endif
! initialize to be on the save side for adjoint runs SIMULATION_TYPE==2
- eps_trace_over_3_crust_mantle(:,:,:,:) = 0._CUSTOM_REAL
- epsilondev_crust_mantle(:,:,:,:,:) = 0._CUSTOM_REAL
- eps_trace_over_3_inner_core(:,:,:,:) = 0._CUSTOM_REAL
- epsilondev_inner_core(:,:,:,:,:) = 0._CUSTOM_REAL
- if(FIX_UNDERFLOW_PROBLEM) then
- eps_trace_over_3_crust_mantle(:,:,:,:) = VERYSMALLVAL
- epsilondev_crust_mantle(:,:,:,:,:) = VERYSMALLVAL
- eps_trace_over_3_inner_core(:,:,:,:) = VERYSMALLVAL
- epsilondev_inner_core(:,:,:,:,:) = VERYSMALLVAL
- endif
+!ZN eps_trace_over_3_crust_mantle(:,:,:,:) = 0._CUSTOM_REAL
+!ZN epsilondev_crust_mantle(:,:,:,:,:) = 0._CUSTOM_REAL
+!ZN eps_trace_over_3_inner_core(:,:,:,:) = 0._CUSTOM_REAL
+!ZN epsilondev_inner_core(:,:,:,:,:) = 0._CUSTOM_REAL
+!ZN if(FIX_UNDERFLOW_PROBLEM) then
+!ZN eps_trace_over_3_crust_mantle(:,:,:,:) = VERYSMALLVAL
+!ZN epsilondev_crust_mantle(:,:,:,:,:) = VERYSMALLVAL
+!ZN eps_trace_over_3_inner_core(:,:,:,:) = VERYSMALLVAL
+!ZN epsilondev_inner_core(:,:,:,:,:) = VERYSMALLVAL
+!ZN endif
if (COMPUTE_AND_STORE_STRAIN) then
if(MOVIE_VOLUME .and. (MOVIE_VOLUME_TYPE == 2 .or. MOVIE_VOLUME_TYPE == 3)) then
@@ -2087,13 +2089,13 @@
displ_inner_core,veloc_inner_core,accel_inner_core, &
displ_outer_core,veloc_outer_core,accel_outer_core, &
R_memory_crust_mantle,R_memory_inner_core, &
- epsilondev_crust_mantle,epsilondev_inner_core, &
+!ZN epsilondev_crust_mantle,epsilondev_inner_core, &
A_array_rotation,B_array_rotation, &
b_displ_crust_mantle,b_veloc_crust_mantle,b_accel_crust_mantle, &
b_displ_inner_core,b_veloc_inner_core,b_accel_inner_core, &
b_displ_outer_core,b_veloc_outer_core,b_accel_outer_core, &
b_R_memory_crust_mantle,b_R_memory_inner_core, &
- b_epsilondev_crust_mantle,b_epsilondev_inner_core, &
+!ZN b_epsilondev_crust_mantle,b_epsilondev_inner_core, &
b_A_array_rotation,b_B_array_rotation,LOCAL_PATH)
endif
@@ -2233,7 +2235,6 @@
if(SIMULATION_TYPE == 2)then
!!add this part
-!ZN !ZN we need to be careful to arrange this part
it = 0
do iteration_on_subset = 1, NSTEP / NT_DUMP
@@ -2246,12 +2247,6 @@
include "part1_undo_att.f90"
- ! save source derivatives for adjoint simulations
- if (SIMULATION_TYPE == 2 .and. nrec_local > 0) then
- call save_kernels_source_derivatives(nrec_local,NSOURCES,scale_displ,scale_t, &
- nu_source,moment_der,sloc_der,stshift_der,shdur_der,number_receiver_global)
- endif
-
enddo
enddo
endif
@@ -2323,29 +2318,29 @@
include "part1_undo_att.f90"
- call compute_strain_crust_mantle(displ_crust_mantle,hprime_xx,hprime_yy,hprime_zz,ibool_crust_mantle,&
- xix_crust_mantle,xiy_crust_mantle,xiz_crust_mantle,&
- etax_crust_mantle,etay_crust_mantle,etaz_crust_mantle,&
- gammax_crust_mantle,gammay_crust_mantle,gammaz_crust_mantle, &
- epsilondev_crust_mantle,eps_trace_over_3_crust_mantle)
+!ZN call compute_strain_crust_mantle(displ_crust_mantle,hprime_xx,hprime_yy,hprime_zz,ibool_crust_mantle,&
+!ZN xix_crust_mantle,xiy_crust_mantle,xiz_crust_mantle,&
+!ZN etax_crust_mantle,etay_crust_mantle,etaz_crust_mantle,&
+!ZN gammax_crust_mantle,gammay_crust_mantle,gammaz_crust_mantle, &
+!ZN epsilondev_crust_mantle,eps_trace_over_3_crust_mantle)
- call compute_strain_crust_mantle(b_displ_crust_mantle,hprime_xx,hprime_yy,hprime_zz,ibool_crust_mantle,&
- xix_crust_mantle,xiy_crust_mantle,xiz_crust_mantle,&
- etax_crust_mantle,etay_crust_mantle,etaz_crust_mantle,&
- gammax_crust_mantle,gammay_crust_mantle,gammaz_crust_mantle, &
- b_epsilondev_crust_mantle,b_eps_trace_over_3_crust_mantle)
+!ZN call compute_strain_crust_mantle(b_displ_crust_mantle,hprime_xx,hprime_yy,hprime_zz,ibool_crust_mantle,&
+!ZN xix_crust_mantle,xiy_crust_mantle,xiz_crust_mantle,&
+!ZN etax_crust_mantle,etay_crust_mantle,etaz_crust_mantle,&
+!ZN gammax_crust_mantle,gammay_crust_mantle,gammaz_crust_mantle, &
+!ZN b_epsilondev_crust_mantle,b_eps_trace_over_3_crust_mantle)
- call compute_strain_inner_core(displ_inner_core,hprime_xx,hprime_yy,hprime_zz,ibool_inner_core,&
- xix_inner_core,xiy_inner_core,xiz_inner_core,&
- etax_inner_core,etay_inner_core,etaz_inner_core,&
- gammax_inner_core,gammay_inner_core,gammaz_inner_core, &
- epsilondev_inner_core,eps_trace_over_3_inner_core)
+!ZN call compute_strain_inner_core(displ_inner_core,hprime_xx,hprime_yy,hprime_zz,ibool_inner_core,&
+!ZN xix_inner_core,xiy_inner_core,xiz_inner_core,&
+!ZN etax_inner_core,etay_inner_core,etaz_inner_core,&
+!ZN gammax_inner_core,gammay_inner_core,gammaz_inner_core, &
+!ZN epsilondev_inner_core,eps_trace_over_3_inner_core)
- call compute_strain_inner_core(b_displ_inner_core,hprime_xx,hprime_yy,hprime_zz,ibool_inner_core,&
- xix_inner_core,xiy_inner_core,xiz_inner_core,&
- etax_inner_core,etay_inner_core,etaz_inner_core,&
- gammax_inner_core,gammay_inner_core,gammaz_inner_core, &
- b_epsilondev_inner_core,b_eps_trace_over_3_inner_core)
+!ZN call compute_strain_inner_core(b_displ_inner_core,hprime_xx,hprime_yy,hprime_zz,ibool_inner_core,&
+!ZN xix_inner_core,xiy_inner_core,xiz_inner_core,&
+!ZN etax_inner_core,etay_inner_core,etaz_inner_core,&
+!ZN gammax_inner_core,gammay_inner_core,gammaz_inner_core, &
+!ZN b_epsilondev_inner_core,b_eps_trace_over_3_inner_core)
include "part3_kernel_computation.f90"
@@ -2447,7 +2442,7 @@
displ_inner_core,veloc_inner_core,accel_inner_core, &
displ_outer_core,veloc_outer_core,accel_outer_core, &
R_memory_crust_mantle,R_memory_inner_core, &
- epsilondev_crust_mantle,epsilondev_inner_core, &
+!ZN epsilondev_crust_mantle,epsilondev_inner_core, &
A_array_rotation,B_array_rotation,LOCAL_PATH)
! synchronize all processes
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/write_movie_volume.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/write_movie_volume.f90 2013-06-15 11:05:55 UTC (rev 22309)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/undo_att/src/specfem3D/write_movie_volume.f90 2013-06-15 15:41:35 UTC (rev 22310)
@@ -255,11 +255,18 @@
! ---------------------------------------------
- subroutine write_movie_volume_strains(myrank,npoints_3dmovie,LOCAL_PATH,MOVIE_VOLUME_TYPE,MOVIE_COARSE, &
- it,eps_trace_over_3_crust_mantle,epsilondev_crust_mantle,muvstore_crust_mantle_3dmovie, &
- mask_3dmovie,nu_3dmovie)
+!ZN subroutine write_movie_volume_strains(myrank,npoints_3dmovie,LOCAL_PATH,MOVIE_VOLUME_TYPE,MOVIE_COARSE, &
+!ZN it,eps_trace_over_3_crust_mantle,epsilondev_crust_mantle,muvstore_crust_mantle_3dmovie, &
+!ZN mask_3dmovie,nu_3dmovie)
+ subroutine write_movie_volume_strains(myrank,npoints_3dmovie,LOCAL_PATH,MOVIE_VOLUME_TYPE,MOVIE_COARSE, & !ZN
+ it,muvstore_crust_mantle_3dmovie,mask_3dmovie,nu_3dmovie,& !ZN
+ hprime_xx,hprime_yy,hprime_zz,ibool_crust_mantle,&
+ xix_crust_mantle,xiy_crust_mantle,xiz_crust_mantle,&
+ etax_crust_mantle,etay_crust_mantle,etaz_crust_mantle,&
+ gammax_crust_mantle,gammay_crust_mantle,gammaz_crust_mantle,displ_crust_mantle) !ZN
+
implicit none
include "constants.h"
@@ -267,14 +274,29 @@
! input
integer :: myrank,npoints_3dmovie,MOVIE_VOLUME_TYPE,it
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: eps_trace_over_3_crust_mantle
- real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: epsilondev_crust_mantle
+!ZN real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: eps_trace_over_3_crust_mantle
+!ZN real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: epsilondev_crust_mantle
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: eps_trace_over_3_loc_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc_crust_mantle
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: &
+ xix_crust_mantle,xiy_crust_mantle,xiz_crust_mantle,&
+ etax_crust_mantle,etay_crust_mantle,etaz_crust_mantle, &
+ gammax_crust_mantle,gammay_crust_mantle,gammaz_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: ibool_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx
+ real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLY) :: hprime_yy
+ real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ) :: hprime_zz
+
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_3DMOVIE) :: muvstore_crust_mantle_3dmovie
logical, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: mask_3dmovie
logical :: MOVIE_COARSE
real(kind=CUSTOM_REAL), dimension(3,3,npoints_3dmovie) :: nu_3dmovie
character(len=150) LOCAL_PATH,outputname
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CRUST_MANTLE) :: displ_crust_mantle
+
! variables
character(len=150) prname
integer :: ipoints_3dmovie,i,j,k,ispec,NIT
@@ -310,19 +332,35 @@
write(prname,"('proc',i6.6)") myrank
ipoints_3dmovie=0
do ispec=1,NSPEC_CRUST_MANTLE
+ call compute_element_strain_undo_att_noDev(ispec,nglob_crust_mantle,NSPEC_CRUST_MANTLE,&
+ displ_crust_mantle,hprime_xx,hprime_yy,hprime_zz,ibool_crust_mantle,&
+ xix_crust_mantle,xiy_crust_mantle,xiz_crust_mantle,&
+ etax_crust_mantle,etay_crust_mantle,etaz_crust_mantle,&
+ gammax_crust_mantle,gammay_crust_mantle,gammaz_crust_mantle,&
+ epsilondev_loc_crust_mantle,eps_trace_over_3_loc_crust_mantle)
+
do k=1,NGLLZ,NIT
do j=1,NGLLY,NIT
do i=1,NGLLX,NIT
if(mask_3dmovie(i,j,k,ispec)) then
ipoints_3dmovie=ipoints_3dmovie+1
muv_3dmovie=muvstore_crust_mantle_3dmovie(i,j,k,ispec)
- eps_loc(1,1)=eps_trace_over_3_crust_mantle(i,j,k,ispec) + epsilondev_crust_mantle(1,i,j,k,ispec)
- eps_loc(2,2)=eps_trace_over_3_crust_mantle(i,j,k,ispec) + epsilondev_crust_mantle(2,i,j,k,ispec)
- eps_loc(3,3)=eps_trace_over_3_crust_mantle(i,j,k,ispec)- &
- epsilondev_crust_mantle(1,i,j,k,ispec) - epsilondev_crust_mantle(2,i,j,k,ispec)
- eps_loc(1,2)=epsilondev_crust_mantle(3,i,j,k,ispec)
- eps_loc(1,3)=epsilondev_crust_mantle(4,i,j,k,ispec)
- eps_loc(2,3)=epsilondev_crust_mantle(5,i,j,k,ispec)
+!ZN eps_loc(1,1)=eps_trace_over_3_crust_mantle(i,j,k,ispec) + epsilondev_crust_mantle(1,i,j,k,ispec)
+!ZN eps_loc(2,2)=eps_trace_over_3_crust_mantle(i,j,k,ispec) + epsilondev_crust_mantle(2,i,j,k,ispec)
+!ZN eps_loc(3,3)=eps_trace_over_3_crust_mantle(i,j,k,ispec)- &
+!ZN epsilondev_crust_mantle(1,i,j,k,ispec) - epsilondev_crust_mantle(2,i,j,k,ispec)
+!ZN eps_loc(1,2)=epsilondev_crust_mantle(3,i,j,k,ispec)
+!ZN eps_loc(1,3)=epsilondev_crust_mantle(4,i,j,k,ispec)
+!ZN eps_loc(2,3)=epsilondev_crust_mantle(5,i,j,k,ispec)
+
+ eps_loc(1,1)=eps_trace_over_3_loc_crust_mantle(i,j,k) + epsilondev_loc_crust_mantle(1,i,j,k)
+ eps_loc(2,2)=eps_trace_over_3_loc_crust_mantle(i,j,k) + epsilondev_loc_crust_mantle(2,i,j,k)
+ eps_loc(3,3)=eps_trace_over_3_loc_crust_mantle(i,j,k)- &
+ epsilondev_loc_crust_mantle(1,i,j,k) - epsilondev_loc_crust_mantle(2,i,j,k)
+ eps_loc(1,2)=epsilondev_loc_crust_mantle(3,i,j,k)
+ eps_loc(1,3)=epsilondev_loc_crust_mantle(4,i,j,k)
+ eps_loc(2,3)=epsilondev_loc_crust_mantle(5,i,j,k)
+
eps_loc(2,1)=eps_loc(1,2)
eps_loc(3,1)=eps_loc(1,3)
eps_loc(3,2)=eps_loc(2,3)
More information about the CIG-COMMITS
mailing list