[cig-commits] r22572 - seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D

xie.zhinan at geodynamics.org xie.zhinan at geodynamics.org
Fri Jul 12 00:42:46 PDT 2013


Author: xie.zhinan
Date: 2013-07-12 00:42:46 -0700 (Fri, 12 Jul 2013)
New Revision: 22572

Modified:
   seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_element.F90
   seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_crust_mantle_Dev.F90
   seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_crust_mantle_noDev.f90
   seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_inner_core_Dev.F90
   seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_inner_core_noDev.f90
   seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_kernels.f90
   seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part1_classical.f90
   seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part1_undo_att.f90
   seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part2_classical.f90
   seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part2_undo_att.f90
   seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part3_kernel_computation.f90
   seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.F90
Log:
programing the code for flag RECOMPUTE_STRAIN_DO_NOT_STORE   = .false.


Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_element.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_element.F90	2013-07-12 00:33:16 UTC (rev 22571)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_element.F90	2013-07-12 07:42:46 UTC (rev 22572)
@@ -35,7 +35,7 @@
                     R_memory, &
                     one_minus_sum_beta,vnspec, &
                     tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
-                    dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,rho_s_H,PARTIAL_PHYS_DISPERSION_ONLY)
+                    dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,eps_trace_over_3_loc,rho_s_H,PARTIAL_PHYS_DISPERSION_ONLY)
 
 ! this routine is optimized for NGLLX = NGLLY = NGLLZ = 5 using the Deville et al. (2002) inlined matrix-matrix products
 
@@ -87,6 +87,7 @@
 
   real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NGLLZ) :: rho_s_H
   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
   real(kind=CUSTOM_REAL) one_minus_sum_beta_use
@@ -165,6 +166,7 @@
 !ZN Both expressions are fine, but we need to keep in mind that if we put the 1/2 factor here there we need to remove it
 !ZN from the expression in which we use the strain later in the code.
           templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+          eps_trace_over_3_loc(i,j,k) = 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
@@ -355,7 +357,7 @@
                     R_memory, &
                     one_minus_sum_beta,vnspec, &
                     tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
-                    dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,rho_s_H,PARTIAL_PHYS_DISPERSION_ONLY)
+                    dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,eps_trace_over_3_loc,rho_s_H,PARTIAL_PHYS_DISPERSION_ONLY)
 
 ! this routine is optimized for NGLLX = NGLLY = NGLLZ = 5 using the Deville et al. (2002) inlined matrix-matrix products
 
@@ -409,6 +411,7 @@
   real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: dummyx_loc,dummyy_loc,dummyz_loc
   real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NGLLZ) :: rho_s_H
   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
   real(kind=CUSTOM_REAL) one_minus_sum_beta_use
@@ -500,6 +503,7 @@
 !ZN Both expressions are fine, but we need to keep in mind that if we put the 1/2 factor here there we need to remove it
 !ZN from the expression in which we use the strain later in the code.
           templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+          eps_trace_over_3_loc(i,j,k) = 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
@@ -887,7 +891,7 @@
                     R_memory, &
                     one_minus_sum_beta,vnspec, &
                     tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
-                    dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,rho_s_H,PARTIAL_PHYS_DISPERSION_ONLY)
+                    dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,eps_trace_over_3_loc,rho_s_H,PARTIAL_PHYS_DISPERSION_ONLY) 
 
 ! this routine is optimized for NGLLX = NGLLY = NGLLZ = 5 using the Deville et al. (2002) inlined matrix-matrix products
 
@@ -941,6 +945,7 @@
   real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: dummyx_loc,dummyy_loc,dummyz_loc
   real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NGLLZ) :: rho_s_H
   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
   real(kind=CUSTOM_REAL) one_minus_sum_beta_use
@@ -1021,6 +1026,7 @@
 !ZN Both expressions are fine, but we need to keep in mind that if we put the 1/2 factor here there we need to remove it
 !ZN from the expression in which we use the strain later in the code.
           templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+          eps_trace_over_3_loc(i,j,k) = 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
@@ -1688,7 +1694,8 @@
 !--------------------------------------------------------------------------------------------
 !
  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)
+                                       xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,epsilondev_loc_nplus1,&
+                                       eps_trace_over_3_loc_nplus1) 
 
   implicit none
   include "constants.h"
@@ -1702,6 +1709,7 @@
         xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
 
   real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc_nplus1
+  real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: eps_trace_over_3_loc_nplus1
 
 !  local variable
   integer :: i,j,k,iglob
@@ -1854,6 +1862,7 @@
         duzdyl_plus_duydzl = duzdyl + duydzl
 
         templ = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+        eps_trace_over_3_loc_nplus1 = templ
         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

Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_crust_mantle_Dev.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_crust_mantle_Dev.F90	2013-07-12 00:33:16 UTC (rev 22571)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_crust_mantle_Dev.F90	2013-07-12 07:42:46 UTC (rev 22572)
@@ -54,7 +54,8 @@
           ibool,ispec_is_tiso, &
           R_memory,one_minus_sum_beta,deltat,veloc_crust_mantle, &
           alphaval,betaval,gammaval,factor_common,vnspec,PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev,eps_trace_over_3) 
 
 ! this routine is optimized for NGLLX = NGLLY = NGLLZ = 5 using the Deville et al. (2002) inlined matrix-matrix products
 
@@ -102,7 +103,9 @@
   ! 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
-  logical :: PARTIAL_PHYS_DISPERSION_ONLY
+  logical :: PARTIAL_PHYS_DISPERSION_ONLY,RECOMPUTE_STRAIN_DO_NOT_STORE   
+  real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: epsilondev 
+  real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: eps_trace_over_3 
 
   integer :: vnspec
 
@@ -158,7 +161,10 @@
   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(NGLLX,NGLLY,NGLLZ) :: eps_trace_over_3_loc 
   real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc_nplus1
+  real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: eps_trace_over_3_loc_nplus1
+  real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc_nsub1 
   real(kind=CUSTOM_REAL) fac1,fac2,fac3
 
   ! for gravity
@@ -379,7 +385,7 @@
             R_memory, &
             one_minus_sum_beta,vnspec, &
             tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
-            dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,rho_s_H,PARTIAL_PHYS_DISPERSION_ONLY)
+            dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,eps_trace_over_3_loc,rho_s_H,PARTIAL_PHYS_DISPERSION_ONLY) 
     else
        if( .not. ispec_is_tiso(ispec) ) then
           ! isotropic element
@@ -393,7 +399,7 @@
                R_memory, &
                one_minus_sum_beta,vnspec, &
                tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
-               dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,rho_s_H,PARTIAL_PHYS_DISPERSION_ONLY)
+               dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,eps_trace_over_3_loc,rho_s_H,PARTIAL_PHYS_DISPERSION_ONLY) 
        else
           ! transverse isotropic element
 
@@ -407,7 +413,7 @@
                R_memory, &
                one_minus_sum_beta,vnspec, &
                tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
-               dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,rho_s_H,PARTIAL_PHYS_DISPERSION_ONLY)
+               dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,eps_trace_over_3_loc,rho_s_H,PARTIAL_PHYS_DISPERSION_ONLY) 
        endif ! .not. ispec_is_tiso
     endif
 
@@ -527,18 +533,42 @@
     ! 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_CRUST_MANTLE,NSPEC_CRUST_MANTLE,displ_crust_mantle,veloc_crust_mantle,&
+      if(COMPUTE_AND_STORE_STRAIN .and. (.not. RECOMPUTE_STRAIN_DO_NOT_STORE))then  
+          ! updates R_memory
+          epsilondev_loc_nsub1(1,:,:,:) = epsilondev(1,:,:,:,ispec)
+          epsilondev_loc_nsub1(2,:,:,:) = epsilondev(2,:,:,:,ispec)
+          epsilondev_loc_nsub1(3,:,:,:) = epsilondev(3,:,:,:,ispec)
+          epsilondev_loc_nsub1(4,:,:,:) = epsilondev(4,:,:,:,ispec)
+          epsilondev_loc_nsub1(5,:,:,:) = epsilondev(5,:,:,:,ispec)
+          ! updates R_memory
+          call compute_element_att_memory_cr(ispec,R_memory, &
+                                         vnspec,factor_common, &
+                                         alphaval,betaval,gammaval, &
+                                         c44store,muvstore, &
+                                         epsilondev_loc,epsilondev_loc_nsub1,&  
+                                         istage,R_memory_lddrk,tau_sigma_CUSTOM_REAL,deltat,USE_LDDRK)
+      else
+        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, &
+                                          xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,epsilondev_loc_nplus1,&
+                                          eps_trace_over_3_loc_nplus1)
+        ! updates R_memory
+        call compute_element_att_memory_cr(ispec,R_memory, &
                                          vnspec,factor_common, &
                                          alphaval,betaval,gammaval, &
                                          c44store,muvstore, &
                                          epsilondev_loc_nplus1,epsilondev_loc,&
                                          istage,R_memory_lddrk,tau_sigma_CUSTOM_REAL,deltat,USE_LDDRK)
+      endif
 
+      if(COMPUTE_AND_STORE_STRAIN .and. (.not. RECOMPUTE_STRAIN_DO_NOT_STORE)) then 
+        eps_trace_over_3(:,:,:,ispec) = eps_trace_over_3_loc(:,:,:)
+        epsilondev(1,:,:,:,ispec) = epsilondev_loc(1,:,:,:)
+        epsilondev(2,:,:,:,ispec) = epsilondev_loc(2,:,:,:)
+        epsilondev(3,:,:,:,ispec) = epsilondev_loc(3,:,:,:)
+        epsilondev(4,:,:,:,ispec) = epsilondev_loc(4,:,:,:)
+        epsilondev(5,:,:,:,ispec) = epsilondev_loc(5,:,:,:)  
+      endif
     endif
 
 ! end ispec loop

Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_crust_mantle_noDev.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_crust_mantle_noDev.f90	2013-07-12 00:33:16 UTC (rev 22571)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_crust_mantle_noDev.f90	2013-07-12 07:42:46 UTC (rev 22572)
@@ -54,7 +54,8 @@
           ibool,ispec_is_tiso, &
           R_memory,one_minus_sum_beta,deltat,veloc_crust_mantle, &
           alphaval,betaval,gammaval,factor_common,vnspec,PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev,eps_trace_over_3) 
 
   implicit none
 
@@ -83,7 +84,9 @@
 ! 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
-  logical :: PARTIAL_PHYS_DISPERSION_ONLY
+  logical :: PARTIAL_PHYS_DISPERSION_ONLY,RECOMPUTE_STRAIN_DO_NOT_STORE   
+  real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: epsilondev 
+  real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: eps_trace_over_3 
 
 ! [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
@@ -92,6 +95,7 @@
 
   real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc
   real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc_nplus1
+  real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc_nsub1 
 
 ! arrays with mesh parameters per slice
   integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: ibool
@@ -379,6 +383,9 @@
             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
+            if(.not. RECOMPUTE_STRAIN_DO_NOT_STORE)then  
+              eps_trace_over_3(i,j,k,ispec) = templ
+            endif
           endif
 
           ! precompute terms for attenuation if needed
@@ -906,9 +913,17 @@
 
     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,&
+      if(COMPUTE_AND_STORE_STRAIN .and. (.not. RECOMPUTE_STRAIN_DO_NOT_STORE))then  
+        epsilondev_loc_nsub1(1,:,:,:) = epsilondev(1,:,:,:,ispec)
+        epsilondev_loc_nsub1(2,:,:,:) = epsilondev(2,:,:,:,ispec)
+        epsilondev_loc_nsub1(3,:,:,:) = epsilondev(3,:,:,:,ispec)
+        epsilondev_loc_nsub1(4,:,:,:) = epsilondev(4,:,:,:,ispec)
+        epsilondev_loc_nsub1(5,:,:,:) = epsilondev(5,:,:,:,ispec)
+      else
+        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)
+      endif
 
 ! use Runge-Kutta scheme to march in time
       do i_SLS = 1,N_SLS
@@ -951,14 +966,28 @@
            R_memory(i_memory,i_SLS,:,:,:,ispec) = R_memory(i_memory,i_SLS,:,:,:,ispec) + &
                                                BETA_LDDRK(istage) * R_memory_lddrk(i_memory,i_SLS,:,:,:,ispec)
          else
-           R_memory(i_memory,i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_memory(i_memory,i_SLS,:,:,:,ispec) &
+           if(COMPUTE_AND_STORE_STRAIN .and. (.not. RECOMPUTE_STRAIN_DO_NOT_STORE))then  
+             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_loc_nsub1(i_memory,:,:,:) + gammaval(i_SLS) * epsilondev_loc(i_memory,:,:,:)) 
+           else
+             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_loc(i_memory,:,:,:) + gammaval(i_SLS) * epsilondev_loc_nplus1(i_memory,:,:,:))
+           endif
          endif
 
-        enddo
-      enddo
+        enddo ! i_memory = 1,5
+      enddo ! i_SLS = 1,N_SLS
 
+      if(COMPUTE_AND_STORE_STRAIN .and. (.not. RECOMPUTE_STRAIN_DO_NOT_STORE))then  
+        epsilondev(1,:,:,:,ispec) = epsilondev_loc(1,:,:,:) 
+        epsilondev(2,:,:,:,ispec) = epsilondev_loc(2,:,:,:)
+        epsilondev(3,:,:,:,ispec) = epsilondev_loc(3,:,:,:)
+        epsilondev(4,:,:,:,ispec) = epsilondev_loc(4,:,:,:)
+        epsilondev(5,:,:,:,ispec) = epsilondev_loc(5,:,:,:) 
+      endif
+
     endif
 
   enddo   ! spectral element loop NSPEC_CRUST_MANTLE

Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_inner_core_Dev.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_inner_core_Dev.F90	2013-07-12 00:33:16 UTC (rev 22571)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_inner_core_Dev.F90	2013-07-12 07:42:46 UTC (rev 22572)
@@ -50,7 +50,8 @@
           c11store,c33store,c12store,c13store,c44store,R_memory,one_minus_sum_beta,deltat,veloc_inner_core,&
           alphaval,betaval,gammaval,factor_common, &
           vnspec,PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev,eps_trace_over_3) 
 
 ! this routine is optimized for NGLLX = NGLLY = NGLLZ = 5 using the Deville et al. (2002) inlined matrix-matrix products
 
@@ -85,7 +86,9 @@
   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
-  logical :: PARTIAL_PHYS_DISPERSION_ONLY
+  logical :: PARTIAL_PHYS_DISPERSION_ONLY,RECOMPUTE_STRAIN_DO_NOT_STORE   
+  real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ATTENUATION) :: epsilondev 
+  real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ATTENUATION) :: eps_trace_over_3 
 
   ! array with derivatives of Lagrange polynomials and precalculated products
   double precision, dimension(NGLLX,NGLLY,NGLLZ) :: wgll_cube
@@ -147,6 +150,8 @@
   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
+  real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: eps_trace_over_3_loc_nplus1 
+  real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc_nsub1 
 
   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
@@ -423,6 +428,9 @@
               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
+              if(.not. RECOMPUTE_STRAIN_DO_NOT_STORE)then  
+                eps_trace_over_3(i,j,k,ispec) = templ
+              endif
             endif
 
             if(ATTENUATION_3D_VAL) then
@@ -756,17 +764,42 @@
       ! 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, &
+        if(COMPUTE_AND_STORE_STRAIN .and. (.not. RECOMPUTE_STRAIN_DO_NOT_STORE))then  
+          ! updates R_memory
+          epsilondev_loc_nsub1(1,:,:,:) = epsilondev(1,:,:,:,ispec)
+          epsilondev_loc_nsub1(2,:,:,:) = epsilondev(2,:,:,:,ispec)
+          epsilondev_loc_nsub1(3,:,:,:) = epsilondev(3,:,:,:,ispec)
+          epsilondev_loc_nsub1(4,:,:,:) = epsilondev(4,:,:,:,ispec)
+          epsilondev_loc_nsub1(5,:,:,:) = epsilondev(5,:,:,:,ispec)
+
+          call compute_element_att_memory_ic(ispec,R_memory, &
                                       vnspec,factor_common, &
                                       alphaval,betaval,gammaval, &
                                       muvstore, &
+                                      epsilondev_loc,epsilondev_loc_nsub1,& 
+                                      istage,R_memory_lddrk,tau_sigma_CUSTOM_REAL,deltat,USE_LDDRK)
+        else
+          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,&
+                                              eps_trace_over_3_loc_nplus1)
+          ! updates R_memory
+          call compute_element_att_memory_ic(ispec,R_memory, &
+                                      vnspec,factor_common, &
+                                      alphaval,betaval,gammaval, &
+                                      muvstore, &
                                       epsilondev_loc_nplus1,epsilondev_loc,&
                                       istage,R_memory_lddrk,tau_sigma_CUSTOM_REAL,deltat,USE_LDDRK)
+        endif
 
+        if(COMPUTE_AND_STORE_STRAIN .and. (.not. RECOMPUTE_STRAIN_DO_NOT_STORE)) then 
+          epsilondev(1,:,:,:,ispec) = epsilondev_loc(1,:,:,:)
+          epsilondev(2,:,:,:,ispec) = epsilondev_loc(2,:,:,:)
+          epsilondev(3,:,:,:,ispec) = epsilondev_loc(3,:,:,:)
+          epsilondev(4,:,:,:,ispec) = epsilondev_loc(4,:,:,:)
+          epsilondev(5,:,:,:,ispec) = epsilondev_loc(5,:,:,:)  
+        endif
+
       endif
 
     endif   ! end test to exclude fictitious elements in central cube

Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_inner_core_noDev.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_inner_core_noDev.f90	2013-07-12 00:33:16 UTC (rev 22571)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_inner_core_noDev.f90	2013-07-12 07:42:46 UTC (rev 22572)
@@ -51,7 +51,8 @@
           c11store,c33store,c12store,c13store,c44store,R_memory,one_minus_sum_beta,deltat,veloc_inner_core,&
           alphaval,betaval,gammaval,factor_common, &
           vnspec,PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev,eps_trace_over_3) 
 
   implicit none
 
@@ -82,8 +83,10 @@
   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) :: epsilondev_loc,epsilondev_loc_nplus1
-  logical :: PARTIAL_PHYS_DISPERSION_ONLY
+  real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc,epsilondev_loc_nplus1,epsilondev_loc_nsub1
+  logical :: PARTIAL_PHYS_DISPERSION_ONLY,RECOMPUTE_STRAIN_DO_NOT_STORE   
+  real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ATTENUATION) :: epsilondev 
+  real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ATTENUATION) :: eps_trace_over_3 
 
 ! array with the local to global mapping per slice
   integer, dimension(NSPEC_INNER_CORE) :: idoubling
@@ -357,6 +360,9 @@
             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
+            if(.not. RECOMPUTE_STRAIN_DO_NOT_STORE)then  
+              eps_trace_over_3(i,j,k,ispec) = templ
+            endif
           endif
 
           if(ATTENUATION_3D_VAL) then
@@ -662,9 +668,17 @@
 
     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,&
+      if(COMPUTE_AND_STORE_STRAIN .and. (.not. RECOMPUTE_STRAIN_DO_NOT_STORE))then  
+        epsilondev_loc_nsub1(1,i,j,k) = epsilondev(1,i,j,k,ispec)
+        epsilondev_loc_nsub1(2,i,j,k) = epsilondev(2,i,j,k,ispec)
+        epsilondev_loc_nsub1(3,i,j,k) = epsilondev(3,i,j,k,ispec)
+        epsilondev_loc_nsub1(4,i,j,k) = epsilondev(4,i,j,k,ispec)
+        epsilondev_loc_nsub1(5,i,j,k) = epsilondev(5,i,j,k,ispec)
+      else
+        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)
+      endif
 
       do i_SLS = 1,N_SLS
 
@@ -695,18 +709,37 @@
                                                BETA_LDDRK(istage) * R_memory_lddrk(i_memory,i_SLS,:,:,:,ispec)
           enddo
         else
-          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_loc_nplus1(i_memory,:,:,:) + gammaval(i_SLS) * epsilondev_loc(i_memory,:,:,:))
-          enddo
+          if(COMPUTE_AND_STORE_STRAIN .and. (.not. RECOMPUTE_STRAIN_DO_NOT_STORE))then  
+            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_loc(i_memory,:,:,:) + gammaval(i_SLS) * epsilondev_loc_nsub1(i_memory,:,:,:)) 
+            enddo 
+          else
+            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_loc_nplus1(i_memory,:,:,:) + gammaval(i_SLS) * epsilondev_loc(i_memory,:,:,:))
+            enddo
+          endif 
         endif
 
-      enddo
+      enddo !do i_SLS = 1,N_SLS
 
+      if(COMPUTE_AND_STORE_STRAIN .and. (.not. RECOMPUTE_STRAIN_DO_NOT_STORE))then  
+         epsilondev(1,:,:,:,ispec) = epsilondev_loc(1,:,:,:) 
+         epsilondev(2,:,:,:,ispec) = epsilondev_loc(2,:,:,:)
+         epsilondev(3,:,:,:,ispec) = epsilondev_loc(3,:,:,:)
+         epsilondev(4,:,:,:,ispec) = epsilondev_loc(4,:,:,:)
+         epsilondev(5,:,:,:,ispec) = epsilondev_loc(5,:,:,:) 
+      endif
+
     endif
 
   endif   ! end test to exclude fictitious elements in central cube

Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_kernels.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_kernels.f90	2013-07-12 00:33:16 UTC (rev 22571)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_kernels.f90	2013-07-12 07:42:46 UTC (rev 22572)
@@ -30,7 +30,9 @@
                           alpha_kl_crust_mantle,cijkl_kl_crust_mantle, &
                           accel_crust_mantle,b_displ_crust_mantle, &
                           deltat,displ_crust_mantle,hprime_xx,hprime_xxT,&
-                          xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,ANISOTROPIC_KL)
+                          xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,ANISOTROPIC_KL,&
+                          RECOMPUTE_STRAIN_DO_NOT_STORE,& 
+                          epsilondev_crust_mantle,eps_trace_over_3_crust_mantle) 
 
   implicit none
 
@@ -57,6 +59,10 @@
 
   logical :: ANISOTROPIC_KL
 
+  logical :: RECOMPUTE_STRAIN_DO_NOT_STORE 
+  real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: epsilondev_crust_mantle 
+  real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: eps_trace_over_3_crust_mantle 
+
   ! local parameters
   real(kind=CUSTOM_REAL),dimension(21) :: prod
   real(kind=CUSTOM_REAL), dimension(5) :: epsilondev_loc
@@ -69,10 +75,19 @@
   ! crust_mantle
   do ispec = 1, NSPEC_CRUST_MANTLE
 
-    call compute_element_strain_undo_att_Dev(ispec,NGLOB_CRUST_MANTLE,NSPEC_CRUST_MANTLE,&
+    if(COMPUTE_AND_STORE_STRAIN .and. RECOMPUTE_STRAIN_DO_NOT_STORE)then  
+        eps_trace_over_3_loc_matrix(:,:,:) = eps_trace_over_3_crust_mantle(:,:,:,ispec)
+        epsilondev_loc_matrix(1,:,:,:) = epsilondev_crust_mantle(1,:,:,:,ispec)
+        epsilondev_loc_matrix(2,:,:,:) = epsilondev_crust_mantle(2,:,:,:,ispec)
+        epsilondev_loc_matrix(3,:,:,:) = epsilondev_crust_mantle(3,:,:,:,ispec)
+        epsilondev_loc_matrix(4,:,:,:) = epsilondev_crust_mantle(4,:,:,:,ispec)
+        epsilondev_loc_matrix(5,:,:,:) = epsilondev_crust_mantle(5,:,:,:,ispec)
+    else
+      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)
+    endif
 
     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,&
@@ -306,7 +321,9 @@
                           alpha_kl_inner_core, &
                           accel_inner_core,b_displ_inner_core, &
                           deltat,displ_inner_core,hprime_xx,hprime_xxT,&
-                          xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz)
+                          xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,&
+                          RECOMPUTE_STRAIN_DO_NOT_STORE,& 
+                          epsilondev_inner_core,eps_trace_over_3_inner_core) 
 
   implicit none
 
@@ -328,6 +345,10 @@
   real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: &
         xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
 
+  logical :: RECOMPUTE_STRAIN_DO_NOT_STORE 
+  real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: epsilondev_inner_core 
+  real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: eps_trace_over_3_inner_core 
+
   ! local parameters
   real(kind=CUSTOM_REAL), dimension(5) :: b_epsilondev_loc
   real(kind=CUSTOM_REAL), dimension(5) :: epsilondev_loc
@@ -340,10 +361,19 @@
   ! inner_core
   do ispec = 1, NSPEC_INNER_CORE
 
-    call compute_element_strain_undo_att_Dev(ispec,NGLOB_inner_core,NSPEC_inner_core,&
+    if(COMPUTE_AND_STORE_STRAIN .and. RECOMPUTE_STRAIN_DO_NOT_STORE)then  
+        eps_trace_over_3_loc_matrix(:,:,:) = eps_trace_over_3_inner_core(:,:,:,ispec)
+        epsilondev_loc_matrix(1,:,:,:) = epsilondev_inner_core(1,:,:,:,ispec)
+        epsilondev_loc_matrix(2,:,:,:) = epsilondev_inner_core(2,:,:,:,ispec)
+        epsilondev_loc_matrix(3,:,:,:) = epsilondev_inner_core(3,:,:,:,ispec)
+        epsilondev_loc_matrix(4,:,:,:) = epsilondev_inner_core(4,:,:,:,ispec)
+        epsilondev_loc_matrix(5,:,:,:) = epsilondev_inner_core(5,:,:,:,ispec)
+    else
+      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)
+    endif
 
     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,&

Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part1_classical.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part1_classical.f90	2013-07-12 00:33:16 UTC (rev 22571)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part1_classical.f90	2013-07-12 07:42:46 UTC (rev 22572)
@@ -320,7 +320,8 @@
           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,5),PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev_crust_mantle,eps_trace_over_3_crust_mantle) 
     else
       call compute_forces_crust_mantle(minus_gravity_table,density_table,minus_deriv_gravity_table, &
           displ_crust_mantle,accel_crust_mantle, &
@@ -360,7 +361,8 @@
           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,5),PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev_crust_mantle,eps_trace_over_3_crust_mantle) 
     endif
 
     ! Deville routine
@@ -396,7 +398,8 @@
           alphaval,betaval,gammaval, &
           factor_common_inner_core, &
           size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev_inner_core,eps_trace_over_3_inner_core) 
     else
       call compute_forces_inner_core(minus_gravity_table,density_table,minus_deriv_gravity_table, &
           displ_inner_core,accel_inner_core, &
@@ -429,7 +432,8 @@
           alphaval,betaval,gammaval, &
           factor_common_inner_core, &
           size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev_inner_core,eps_trace_over_3_inner_core) 
     endif
 
     ! Stacey
@@ -612,7 +616,8 @@
           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,5),PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev_crust_mantle,eps_trace_over_3_crust_mantle) 
       else
         call compute_forces_crust_mantle(minus_gravity_table,density_table,minus_deriv_gravity_table, &
           displ_crust_mantle,accel_crust_mantle, &
@@ -652,7 +657,8 @@
           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,5),PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev_crust_mantle,eps_trace_over_3_crust_mantle) 
       endif
 
       ! Deville routine
@@ -688,7 +694,8 @@
           alphaval,betaval,gammaval, &
           factor_common_inner_core, &
           size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev_inner_core,eps_trace_over_3_inner_core) 
       else
         call compute_forces_inner_core(minus_gravity_table,density_table,minus_deriv_gravity_table, &
           displ_inner_core,accel_inner_core, &
@@ -721,7 +728,8 @@
           alphaval,betaval,gammaval, &
           factor_common_inner_core, &
           size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev_inner_core,eps_trace_over_3_inner_core) 
       endif
 
 ! assemble all the contributions between slices using MPI

Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part1_undo_att.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part1_undo_att.f90	2013-07-12 00:33:16 UTC (rev 22571)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part1_undo_att.f90	2013-07-12 07:42:46 UTC (rev 22572)
@@ -349,7 +349,8 @@
           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,5),PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev_crust_mantle,eps_trace_over_3_crust_mantle) 
     else
       call compute_forces_crust_mantle(minus_gravity_table,density_table,minus_deriv_gravity_table, &
           displ_crust_mantle,accel_crust_mantle, &
@@ -389,7 +390,8 @@
           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,5),PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev_crust_mantle,eps_trace_over_3_crust_mantle) 
     endif
 
     ! Deville routine
@@ -425,7 +427,8 @@
           alphaval,betaval,gammaval, &
           factor_common_inner_core, &
           size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev_inner_core,eps_trace_over_3_inner_core) 
     else
       call compute_forces_inner_core(minus_gravity_table,density_table,minus_deriv_gravity_table, &
           displ_inner_core,accel_inner_core, &
@@ -458,7 +461,8 @@
           alphaval,betaval,gammaval, &
           factor_common_inner_core, &
           size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev_inner_core,eps_trace_over_3_inner_core) 
     endif
 
     ! Stacey
@@ -637,7 +641,8 @@
           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,5),PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev_crust_mantle,eps_trace_over_3_crust_mantle) 
       else
         call compute_forces_crust_mantle(minus_gravity_table,density_table,minus_deriv_gravity_table, &
           displ_crust_mantle,accel_crust_mantle, &
@@ -677,7 +682,8 @@
           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,5),PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev_crust_mantle,eps_trace_over_3_crust_mantle) 
       endif
 
       ! Deville routine
@@ -713,7 +719,8 @@
           alphaval,betaval,gammaval, &
           factor_common_inner_core, &
           size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev_inner_core,eps_trace_over_3_inner_core) 
       else
         call compute_forces_inner_core(minus_gravity_table,density_table,minus_deriv_gravity_table, &
           displ_inner_core,accel_inner_core, &
@@ -746,7 +753,8 @@
           alphaval,betaval,gammaval, &
           factor_common_inner_core, &
           size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev_inner_core,eps_trace_over_3_inner_core) 
       endif
 
 ! assemble all the contributions between slices using MPI

Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part2_classical.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part2_classical.f90	2013-07-12 00:33:16 UTC (rev 22571)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part2_classical.f90	2013-07-12 07:42:46 UTC (rev 22572)
@@ -329,7 +329,8 @@
           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,5),PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev_crust_mantle,eps_trace_over_3_crust_mantle) 
       else
         call compute_forces_crust_mantle(minus_gravity_table,density_table,minus_deriv_gravity_table, &
           b_displ_crust_mantle,b_accel_crust_mantle, &
@@ -369,7 +370,8 @@
           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,5),PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev_crust_mantle,eps_trace_over_3_crust_mantle) 
       endif
     endif
 
@@ -406,7 +408,8 @@
           b_alphaval,b_betaval,b_gammaval, &
           factor_common_inner_core, &
           size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev_inner_core,eps_trace_over_3_inner_core) 
       else
         call compute_forces_inner_core(minus_gravity_table,density_table,minus_deriv_gravity_table, &
           b_displ_inner_core,b_accel_inner_core, &
@@ -439,7 +442,8 @@
           b_alphaval,b_betaval,b_gammaval, &
           factor_common_inner_core, &
           size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev_inner_core,eps_trace_over_3_inner_core) 
       endif
     endif
 
@@ -605,7 +609,8 @@
           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,5),PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev_crust_mantle,eps_trace_over_3_crust_mantle) 
         else
           call compute_forces_crust_mantle(minus_gravity_table,density_table,minus_deriv_gravity_table, &
           b_displ_crust_mantle,b_accel_crust_mantle, &
@@ -645,7 +650,8 @@
           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,5),PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev_crust_mantle,eps_trace_over_3_crust_mantle) 
         endif
 
         ! Deville routine
@@ -681,7 +687,8 @@
           b_alphaval,b_betaval,b_gammaval, &
           factor_common_inner_core, &
           size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev_inner_core,eps_trace_over_3_inner_core) 
         else
           call compute_forces_inner_core(minus_gravity_table,density_table,minus_deriv_gravity_table, &
           b_displ_inner_core,b_accel_inner_core, &
@@ -714,7 +721,8 @@
           b_alphaval,b_betaval,b_gammaval, &
           factor_common_inner_core, &
           size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev_inner_core,eps_trace_over_3_inner_core) 
         endif
 
 ! assemble all the contributions between slices using MPI

Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part2_undo_att.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part2_undo_att.f90	2013-07-12 00:33:16 UTC (rev 22571)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part2_undo_att.f90	2013-07-12 07:42:46 UTC (rev 22572)
@@ -360,7 +360,8 @@
           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,5),PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev_crust_mantle,eps_trace_over_3_crust_mantle) 
       else
         call compute_forces_crust_mantle(minus_gravity_table,density_table,minus_deriv_gravity_table, &
           b_displ_crust_mantle,b_accel_crust_mantle, &
@@ -400,7 +401,8 @@
           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,5),PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev_crust_mantle,eps_trace_over_3_crust_mantle) 
       endif
     endif
 
@@ -437,7 +439,8 @@
           b_alphaval,b_betaval,b_gammaval, &
           factor_common_inner_core, &
           size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev_inner_core,eps_trace_over_3_inner_core) 
       else
         call compute_forces_inner_core(minus_gravity_table,density_table,minus_deriv_gravity_table, &
           b_displ_inner_core,b_accel_inner_core, &
@@ -470,7 +473,8 @@
           b_alphaval,b_betaval,b_gammaval, &
           factor_common_inner_core, &
           size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev_inner_core,eps_trace_over_3_inner_core) 
       endif
     endif
 
@@ -670,7 +674,8 @@
           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,5),PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev_crust_mantle,eps_trace_over_3_crust_mantle) 
         else
           call compute_forces_crust_mantle(minus_gravity_table,density_table,minus_deriv_gravity_table, &
           b_displ_crust_mantle,b_accel_crust_mantle, &
@@ -710,7 +715,8 @@
           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,5),PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_crust_mantle_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev_crust_mantle,eps_trace_over_3_crust_mantle) 
         endif
 
         ! Deville routine
@@ -746,7 +752,8 @@
           b_alphaval,b_betaval,b_gammaval, &
           factor_common_inner_core, &
           size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev_inner_core,eps_trace_over_3_inner_core) 
         else
           call compute_forces_inner_core(minus_gravity_table,density_table,minus_deriv_gravity_table, &
           b_displ_inner_core,b_accel_inner_core, &
@@ -779,7 +786,8 @@
           b_alphaval,b_betaval,b_gammaval, &
           factor_common_inner_core, &
           size(factor_common_inner_core,5),PARTIAL_PHYS_DISPERSION_ONLY,&
-          istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK)
+          istage,R_memory_inner_core_lddrk,tau_sigma_CUSTOM_REAL,USE_LDDRK,&
+          RECOMPUTE_STRAIN_DO_NOT_STORE,epsilondev_inner_core,eps_trace_over_3_inner_core) 
         endif
 
 ! assemble all the contributions between slices using MPI

Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part3_kernel_computation.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part3_kernel_computation.f90	2013-07-12 00:33:16 UTC (rev 22571)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/part3_kernel_computation.f90	2013-07-12 07:42:46 UTC (rev 22572)
@@ -11,7 +11,9 @@
                           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,ANISOTROPIC_KL)
+                          gammax_crust_mantle,gammay_crust_mantle,gammaz_crust_mantle,ANISOTROPIC_KL,&
+                          RECOMPUTE_STRAIN_DO_NOT_STORE,& 
+                          epsilondev_crust_mantle,eps_trace_over_3_crust_mantle) 
 
     ! outer core
     call compute_kernels_outer_core(ibool_outer_core, &
@@ -36,7 +38,9 @@
                           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)
+                          gammax_inner_core,gammay_inner_core,gammaz_inner_core,&
+                          RECOMPUTE_STRAIN_DO_NOT_STORE,& 
+                          epsilondev_inner_core,eps_trace_over_3_inner_core) 
 
     ! NOISE TOMOGRAPHY --- source strength kernel
     if (NOISE_TOMOGRAPHY == 3)  &

Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.F90	2013-07-12 00:33:16 UTC (rev 22571)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.F90	2013-07-12 07:42:46 UTC (rev 22572)
@@ -393,10 +393,16 @@
   real(kind=CUSTOM_REAL), dimension(N_SLS,ATT1_VAL,ATT2_VAL,ATT3_VAL,ATT5_VAL) :: 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(:,:,:,:,:), allocatable :: epsilondev_crust_mantle 
+  real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: eps_trace_over_3_crust_mantle 
+  real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc 
+  real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: eps_trace_over_3_loc 
   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(:,:,:,:,:), allocatable :: epsilondev_inner_core 
+  real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: eps_trace_over_3_inner_core 
 
 ! to save a significant amount of memory space, we equivalence these two arrays that are never used simultaneously
 ! (ibathy_topo is only used before the beginning of the time loop, while R_memory_crust_mantle is only used inside the time loop)
@@ -1110,6 +1116,47 @@
 ! in exact undoing of attenuation
   undo_att_sim_type_3 = .false.
 
+! ZN if we want to storing the strain to acclerate the code but cost more memory then
+  if(ATTENUATION_VAL .and. COMPUTE_AND_STORE_STRAIN .and. (.not. RECOMPUTE_STRAIN_DO_NOT_STORE))then
+    allocate(epsilondev_crust_mantle(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT),stat=ier)
+    if( ier /= 0 ) call exit_MPI(myrank,'error allocating epsilondev_crust_mantle')
+    allocate(eps_trace_over_3_crust_mantle(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT),stat=ier)
+    if( ier /= 0 ) call exit_MPI(myrank,'error allocating eps_trace_over_3_crust_mantle')
+    allocate(epsilondev_inner_core(5,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ATTENUATION),stat=ier)
+    if( ier /= 0 ) call exit_MPI(myrank,'error allocating epsilondev_inner_core')
+    allocate(eps_trace_over_3_inner_core(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ATTENUATION),stat=ier)
+    if( ier /= 0 ) call exit_MPI(myrank,'error allocating eps_trace_over_3_inner_core')
+    epsilondev_crust_mantle(:,:,:,:,:) = 0._CUSTOM_REAL
+    eps_trace_over_3_crust_mantle(:,:,:,:) = 0._CUSTOM_REAL
+    epsilondev_inner_core(:,:,:,:,:) = 0._CUSTOM_REAL
+    eps_trace_over_3_inner_core(:,:,:,:) = 0._CUSTOM_REAL
+    if(FIX_UNDERFLOW_PROBLEM) then
+      epsilondev_crust_mantle(:,:,:,:,:) = VERYSMALLVAL
+      eps_trace_over_3_crust_mantle(:,:,:,:) = VERYSMALLVAL
+      epsilondev_inner_core(:,:,:,:,:) = VERYSMALLVAL
+      eps_trace_over_3_inner_core(:,:,:,:) = VERYSMALLVAL
+    endif
+  else
+    allocate(epsilondev_crust_mantle(5,NGLLX,NGLLY,NGLLZ,1),stat=ier)
+    if( ier /= 0 ) call exit_MPI(myrank,'error allocating epsilondev_crust_mantle')
+    allocate(eps_trace_over_3_crust_mantle(NGLLX,NGLLY,NGLLZ,1),stat=ier)
+    if( ier /= 0 ) call exit_MPI(myrank,'error allocating eps_trace_over_3_crust_mantle')
+    allocate(epsilondev_inner_core(5,NGLLX,NGLLY,NGLLZ,1),stat=ier)
+    if( ier /= 0 ) call exit_MPI(myrank,'error allocating epsilondev_inner_core')
+    allocate(eps_trace_over_3_inner_core(NGLLX,NGLLY,NGLLZ,1),stat=ier)
+    if( ier /= 0 ) call exit_MPI(myrank,'error allocating eps_trace_over_3_inner_core')
+    epsilondev_crust_mantle(:,:,:,:,:) = 0._CUSTOM_REAL
+    eps_trace_over_3_crust_mantle(:,:,:,:) = 0._CUSTOM_REAL
+    epsilondev_inner_core(:,:,:,:,:) = 0._CUSTOM_REAL
+    eps_trace_over_3_inner_core(:,:,:,:) = 0._CUSTOM_REAL
+    if(FIX_UNDERFLOW_PROBLEM) then
+      epsilondev_crust_mantle(:,:,:,:,:) = VERYSMALLVAL
+      eps_trace_over_3_crust_mantle(:,:,:,:) = VERYSMALLVAL
+      epsilondev_inner_core(:,:,:,:,:) = VERYSMALLVAL
+      eps_trace_over_3_inner_core(:,:,:,:) = VERYSMALLVAL
+    endif
+  endif
+
 !
 !-------------------------------------------------------------------------------------------------
 !-------------------------------------------------------------------------------------------------
@@ -2387,6 +2434,41 @@
       it_temp = it
       seismo_current_temp = seismo_current
 
+        if(COMPUTE_AND_STORE_STRAIN .and. (.not. RECOMPUTE_STRAIN_DO_NOT_STORE))then  
+          if(.not. USE_DEVILLE_PRODUCTS_VAL) &
+             call exit_MPI(myrank,'COMPUTE_AND_STORE_STRAIN .and. (.not. RECOMPUTE_STRAIN_DO_NOT_STORE) is .true.'// &
+                           'is not implemented without USE_DEVILLE_PRODUCTS_VAL') 
+          do ispec = 1, NSPEC_INNER_CORE
+            call compute_element_strain_att_Dev(ispec,NGLOB_INNER_CORE,NSPEC_INNER_CORE,b_displ_inner_core,&
+                                            b_veloc_inner_core,0._CUSTOM_REAL,ibool_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,&
+                                            epsilondev_loc,eps_trace_over_3_loc)
+            eps_trace_over_3_inner_core(:,:,:,ispec) = eps_trace_over_3_loc(:,:,:)
+            epsilondev_inner_core(1,:,:,:,ispec) = epsilondev_loc(1,:,:,:)
+            epsilondev_inner_core(2,:,:,:,ispec) = epsilondev_loc(2,:,:,:)
+            epsilondev_inner_core(3,:,:,:,ispec) = epsilondev_loc(3,:,:,:)
+            epsilondev_inner_core(4,:,:,:,ispec) = epsilondev_loc(4,:,:,:)
+            epsilondev_inner_core(5,:,:,:,ispec) = epsilondev_loc(5,:,:,:)          
+          enddo
+
+          do ispec = 1, NSPEC_crust_mantle
+            call compute_element_strain_att_Dev(ispec,NGLOB_crust_mantle,NSPEC_crust_mantle,b_displ_crust_mantle,&
+                                            b_veloc_crust_mantle,0._CUSTOM_REAL,ibool_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,&
+                                            epsilondev_loc,eps_trace_over_3_loc)
+            eps_trace_over_3_crust_mantle(:,:,:,ispec) = eps_trace_over_3_loc(:,:,:)
+            epsilondev_crust_mantle(1,:,:,:,ispec) = epsilondev_loc(1,:,:,:)
+            epsilondev_crust_mantle(2,:,:,:,ispec) = epsilondev_loc(2,:,:,:)
+            epsilondev_crust_mantle(3,:,:,:,ispec) = epsilondev_loc(3,:,:,:)
+            epsilondev_crust_mantle(4,:,:,:,ispec) = epsilondev_loc(4,:,:,:)
+            epsilondev_crust_mantle(5,:,:,:,ispec) = epsilondev_loc(5,:,:,:)          
+          enddo
+        endif
+
       do it_of_this_subset = 1, NT_DUMP_ATTENUATION
 
         it = it + 1
@@ -2403,6 +2485,41 @@
       it = it_temp
       seismo_current = seismo_current_temp
 
+      if(COMPUTE_AND_STORE_STRAIN .and. (.not. RECOMPUTE_STRAIN_DO_NOT_STORE))then  
+        if(.not. USE_DEVILLE_PRODUCTS_VAL) &
+           call exit_MPI(myrank,'COMPUTE_AND_STORE_STRAIN .and. (.not. RECOMPUTE_STRAIN_DO_NOT_STORE) is .true.'// &
+                         'is not implemented without USE_DEVILLE_PRODUCTS_VAL') 
+        do ispec = 1, NSPEC_INNER_CORE
+          call compute_element_strain_att_Dev(ispec,NGLOB_INNER_CORE,NSPEC_INNER_CORE,displ_inner_core,&
+                                            veloc_inner_core,0._CUSTOM_REAL,ibool_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,&
+                                            epsilondev_loc,eps_trace_over_3_loc)
+          eps_trace_over_3_inner_core(:,:,:,ispec) = eps_trace_over_3_loc(:,:,:)
+          epsilondev_inner_core(1,:,:,:,ispec) = epsilondev_loc(1,:,:,:)
+          epsilondev_inner_core(2,:,:,:,ispec) = epsilondev_loc(2,:,:,:)
+          epsilondev_inner_core(3,:,:,:,ispec) = epsilondev_loc(3,:,:,:)
+          epsilondev_inner_core(4,:,:,:,ispec) = epsilondev_loc(4,:,:,:)
+          epsilondev_inner_core(5,:,:,:,ispec) = epsilondev_loc(5,:,:,:)          
+        enddo
+
+        do ispec = 1, NSPEC_crust_mantle
+          call compute_element_strain_att_Dev(ispec,NGLOB_crust_mantle,NSPEC_crust_mantle,displ_crust_mantle,&
+                                          veloc_crust_mantle,0._CUSTOM_REAL,ibool_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,&
+                                          epsilondev_loc,eps_trace_over_3_loc)
+          eps_trace_over_3_crust_mantle(:,:,:,ispec) = eps_trace_over_3_loc(:,:,:)
+          epsilondev_crust_mantle(1,:,:,:,ispec) = epsilondev_loc(1,:,:,:)
+          epsilondev_crust_mantle(2,:,:,:,ispec) = epsilondev_loc(2,:,:,:)
+          epsilondev_crust_mantle(3,:,:,:,ispec) = epsilondev_loc(3,:,:,:)
+          epsilondev_crust_mantle(4,:,:,:,ispec) = epsilondev_loc(4,:,:,:)
+          epsilondev_crust_mantle(5,:,:,:,ispec) = epsilondev_loc(5,:,:,:)          
+        enddo
+      endif
+
       do it_of_this_subset = 1, NT_DUMP_ATTENUATION
         do i = 1, NDIM
           do j =1,NGLOB_CRUST_MANTLE_ADJOINT



More information about the CIG-COMMITS mailing list