[cig-commits] r21825 - in seismo/2D/SPECFEM2D/trunk/src: meshfem2D specfem2D

[dkomati1 at geodynamics.org]

Author: dkomati1
Date: 2013-04-11 15:39:22 -0700 (Thu, 11 Apr 2013)
New Revision: 21825

Modified:
   seismo/2D/SPECFEM2D/trunk/src/meshfem2D/read_materials.f90
   seismo/2D/SPECFEM2D/trunk/src/specfem2D/gmat01.f90
Log:
added a test on negative or null values of Q attenuation factor


Modified: seismo/2D/SPECFEM2D/trunk/src/meshfem2D/read_materials.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/meshfem2D/read_materials.f90	2013-04-11 22:27:00 UTC (rev 21824)
+++ seismo/2D/SPECFEM2D/trunk/src/meshfem2D/read_materials.f90	2013-04-11 22:39:22 UTC (rev 21825)
@@ -77,8 +77,8 @@
   aniso6(:) = 0.d0
   aniso7(:) = 0.d0
   aniso8(:) = 0.d0
-  QKappa(:) = 0.d0
-  Qmu(:) = 0.d0
+  QKappa(:) = 9999.d0
+  Qmu(:) = 9999.d0
   rho_s(:) = 0.d0
   rho_f(:) = 0.d0
   phi(:) = 0.d0
@@ -116,8 +116,9 @@
         QKappa(i) = val5read
         Qmu(i) = val6read
 
-        if(rho_s(i) <= 0.d0 .or. cp(i) <= 0.d0 .or. cs(i) < 0.d0) stop 'negative value of velocity or density'
-        if(QKappa(i) <= 0.d0 .or. Qmu(i) <= 0.d0) stop 'non-positive value of QKappa or Qmu'
+        if(rho_s(i) <= 0.00000001d0 .or. cp(i) <= 0.00000001d0 .or. cs(i) < 0.00000001d0) &
+            stop 'negative value of velocity or density'
+        if(QKappa(i) <= 0.00000001d0 .or. Qmu(i) <= 0.00000001d0) stop 'non-positive value of QKappa or Qmu'
 
         aniso3(i) = val3read
         aniso4(i) = val4read
@@ -141,6 +142,11 @@
         Qmu(i) = val8read
         cp(i) = sqrt(val4read/val0read)
         cs(i) = sqrt(val6read/val0read)
+
+        if(rho_s(i) <= 0.00000001d0 .or. cp(i) <= 0.00000001d0 .or. cs(i) < 0.00000001d0) &
+            stop 'negative value of velocity or density'
+        if(QKappa(i) <= 0.00000001d0 .or. Qmu(i) <= 0.00000001d0) stop 'non-positive value of QKappa or Qmu'
+
      else
 
         ! poroelastic materials
@@ -163,8 +169,9 @@
         if(phi(i) <= 0.d0 .or. tortuosity(i) <= 0.d0) stop 'non-positive value of porosity or tortuosity'
         if(kappa_s(i) <= 0.d0 .or. kappa_f(i) <= 0.d0 .or. kappa_fr(i) <= 0.d0 .or. mu_fr(i) <= 0.d0) then
            stop 'non-positive value of modulus'
-        end if
-        if(Qmu(i) <= 0.d0) stop 'non-positive value of Qmu'
+        endif
+        if(Qmu(i) <= 0.00000001d0) stop 'non-positive value of Qmu'
+
      endif
   enddo
 

Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/gmat01.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/gmat01.f90	2013-04-11 22:27:00 UTC (rev 21824)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/gmat01.f90	2013-04-11 22:39:22 UTC (rev 21825)
@@ -83,8 +83,8 @@
   aniso_array(:,:) = zero
   permeability(:,:) = zero
   poroelastcoef(:,:,:) = zero
-  QKappa_array(:) = zero
-  Qmu_array(:) = zero
+  QKappa_array(:) = 9999.
+  Qmu_array(:) = 9999.
 
   if(myrank == 0 .and. ipass == 1) write(IOUT,100) numat
 
@@ -109,8 +109,10 @@
         ! QKappa and Qmu values
         QKappa = val5
         Qmu = val6
+        if(QKappa <= 0.0000001 .or. Qmu <= 0.0000001) &
+          stop 'negative or null values of Q attenuation factor not allowed; set them equal to 9999 to indicate no attenuation'
 
-        ! Lam'e parameters
+        ! Lame parameters
         lambdaplus2mu = density(1)*cp*cp
         mu = density(1)*cs*cs
         two_mu = 2.d0*mu
@@ -145,11 +147,7 @@
         cp = sqrt(c33/density(1))
         cs = sqrt(c55/density(1))
 
-        ! QKappa and Qmu values
-        !QKappa = val9
-        !Qmu = val10
-
-        ! Lam'e parameters
+        ! Lame parameters
         lambdaplus2mu = density(1)*cp*cp
         mu = density(1)*cs*cs
         two_mu = 2.d0*mu
@@ -182,7 +180,7 @@
         ! Frame properties
         kappa_fr = val9
         mu_fr = val11
-        ! Lam'e parameters for the solid phase and the frame
+        ! Lame parameters for the solid phase and the frame
         lambdaplus2mu_s = kappa_s + FOUR_THIRDS*mu_s
         lambda_s = lambdaplus2mu_s - 2.d0*mu_s
         lambdaplus2mu_fr = kappa_fr + FOUR_THIRDS*mu_fr