[cig-commits] r21297 - seismo/2D/SPECFEM2D/trunk/src/specfem2D

[xie.zhinan at geodynamics.org]

Author: xie.zhinan
Date: 2013-01-25 08:35:47 -0800 (Fri, 25 Jan 2013)
New Revision: 21297

Modified:
   seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90
Log:
clean the code and change the computation of relaxed modulus according the Qkappa  and Qmu definition in Par_file


Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90	2013-01-25 15:22:23 UTC (rev 21296)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/compute_forces_viscoelastic.F90	2013-01-25 16:35:47 UTC (rev 21297)
@@ -198,7 +198,7 @@
   ! material properties of the elastic medium
   real(kind=CUSTOM_REAL) :: mul_unrelaxed_elastic,lambdal_unrelaxed_elastic, &
     lambdaplus2mu_unrelaxed_elastic,kappal,cpl,csl,rhol, &
-    lambdal_relaxed_viscoelastic,mul_relaxed_viscoelastic,lambdalplus2mul_relaxed_viscoel
+    lambdal_relaxed_viscoelastic,mul_relaxed_viscoelastic,lambdalplusmul_relaxed_viscoel
 
   ! for attenuation
   real(kind=CUSTOM_REAL) :: phinu1,phinu2,tauinvnu1,tauinvnu2,theta_n_u,theta_n_v 
@@ -1004,10 +1004,11 @@
                  ! and porous media, Elsevier, p. 124-125, 2007
 
                  ! compute unrelaxed elastic coefficients from formulas in Carcione 2007 page 125
-                 lambdal_relaxed_viscoelastic = (lambdal_unrelaxed_elastic + mul_unrelaxed_elastic) / Mu_nu1(i,j,ispec) &
-                                                - mul_unrelaxed_elastic / Mu_nu2(i,j,ispec)
+                 lambdal_relaxed_viscoelastic = (lambdal_unrelaxed_elastic + 2._CUSTOM_REAL*mul_unrelaxed_elastic/3._CUSTOM_REAL)&
+                                                / Mu_nu1(i,j,ispec) &
+                                                - (2._CUSTOM_REAL*mul_unrelaxed_elastic/3._CUSTOM_REAL) / Mu_nu2(i,j,ispec)
                  mul_relaxed_viscoelastic = mul_unrelaxed_elastic / Mu_nu2(i,j,ispec)
-                 lambdalplus2mul_relaxed_viscoel = lambdal_relaxed_viscoelastic + TWO*mul_relaxed_viscoelastic
+                 lambdalplusmul_relaxed_viscoel = lambdal_relaxed_viscoelastic + mul_relaxed_viscoelastic
 
                  ! compute the stress using the unrelaxed Lame parameters (Carcione 2007 page 125)
                  sigma_xx = lambdaplus2mu_unrelaxed_elastic*dux_dxl + lambdal_unrelaxed_elastic*duz_dzl
@@ -1026,11 +1027,9 @@
                     e13_sum = e13_sum + e13(i,j,ispec,i_sls)
                  enddo
 
-                 sigma_xx = sigma_xx + (lambdal_relaxed_viscoelastic + mul_relaxed_viscoelastic) * e1_sum &
-                                     + TWO * mul_relaxed_viscoelastic * e11_sum
+                 sigma_xx = sigma_xx + lambdalplusmul_relaxed_viscoel * e1_sum + TWO * mul_relaxed_viscoelastic * e11_sum
                  sigma_xz = sigma_xz + mul_relaxed_viscoelastic * e13_sum
-                 sigma_zz = sigma_zz + (lambdal_relaxed_viscoelastic + mul_relaxed_viscoelastic) * e1_sum &
-                                     - TWO * mul_relaxed_viscoelastic * e11_sum
+                 sigma_zz = sigma_zz + lambdalplusmul_relaxed_viscoel * e1_sum - TWO * mul_relaxed_viscoelastic * e11_sum
                  sigma_zx = sigma_xz
 
                  if(PML_BOUNDARY_CONDITIONS .and. is_PML(ispec)) then