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

[dkomati1 at geodynamics.org]

Author: dkomati1
Date: 2012-05-30 14:31:29 -0700 (Wed, 30 May 2012)
New Revision: 20233

Modified:
   seismo/2D/SPECFEM2D/trunk/src/specfem2D/define_external_model.f90
   seismo/2D/SPECFEM2D/trunk/src/specfem2D/read_external_model.f90
Log:
minor improvements made in the context of the definition of the Munk velocity profile


Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/define_external_model.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/define_external_model.f90	2012-05-30 20:47:41 UTC (rev 20232)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/define_external_model.f90	2012-05-30 21:31:29 UTC (rev 20233)
@@ -61,7 +61,9 @@
   double precision, intent(out) :: QKappa_attenuation,Qmu_attenuation
   double precision, intent(out) :: c11,c15,c13,c33,c35,c55
 
-! dummy routine here, just to demonstrate how the model can be assigned
+! completely dummy routine here, just to demonstrate how the model can be assigned
+! and how such a routine can be written
+
    if(iflag_element == 1 .or. x < 1700.d0 .or. y >= 2300.d0) then
      rho = 2000.d0
      vp = 3000.d0
@@ -74,18 +76,22 @@
      c33 = c11
      c35 = 0.d0
      c55 = 75.3d9
-   else
+
+   else if(iflag_element == 2) then
      rho = 2500.d0
      vp = 3600.d0
      vs = vp / 2.d0
-     QKappa_attenuation = 60.
-     Qmu_attenuation = 60.
+     QKappa_attenuation = 120.
+     Qmu_attenuation = 120.
      c11 = 0.d0   ! this means no anisotropy
      c13 = 0.d0
      c15 = 0.d0
      c33 = 0.d0
      c35 = 0.d0
      c55 = 0.d0
+
+   else
+     stop 'wrong flag number in external model'
    endif
 
   end subroutine define_external_model

Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/read_external_model.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/read_external_model.f90	2012-05-30 20:47:41 UTC (rev 20232)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/read_external_model.f90	2012-05-30 21:31:29 UTC (rev 20233)
@@ -123,9 +123,9 @@
                                     c11ext(i,j,ispec),c13ext(i,j,ispec),c15ext(i,j,ispec), &
                                     c33ext(i,j,ispec),c35ext(i,j,ispec),c55ext(i,j,ispec))
 
-          if((c11ext(i,j,ispec) /= 0) .or. (c13ext(i,j,ispec) /= 0) .or. (c15ext(i,j,ispec) /= 0) .or. &
-            (c33ext(i,j,ispec) /= 0) .or. (c35ext(i,j,ispec) /= 0) .or. (c55ext(i,j,ispec) /= 0)) then
-            ! vp, vs : dummy values, trick to avoid floating point errors
+          if(c11ext(i,j,ispec) > TINYVAL .or. c13ext(i,j,ispec) > TINYVAL .or. c15ext(i,j,ispec) > TINYVAL .or. &
+             c33ext(i,j,ispec) > TINYVAL .or. c35ext(i,j,ispec) > TINYVAL .or. c55ext(i,j,ispec) > TINYVAL) then
+            ! vp, vs : assign dummy values, trick to avoid floating point errors in the case of an anisotropic medium
             vpext(i,j,ispec) = 20.d0
             vsext(i,j,ispec) = 10.d0
           end if
@@ -153,15 +153,15 @@
           (vsext(i,j,ispec) < TINYVAL .and. previous_vsext >= TINYVAL)))  &
           call exit_MPI('external velocity model cannot be both fluid and solid inside the same spectral element')
 
-        if((c11ext(i,j,ispec) /= 0) .or. (c13ext(i,j,ispec) /= 0) .or. (c15ext(i,j,ispec) /= 0) .or. &
-          (c33ext(i,j,ispec) /= 0) .or. (c35ext(i,j,ispec) /= 0) .or. (c55ext(i,j,ispec) /= 0)) then
+        if(c11ext(i,j,ispec) > TINYVAL .or. c13ext(i,j,ispec) > TINYVAL .or. c15ext(i,j,ispec) > TINYVAL .or. &
+           c33ext(i,j,ispec) > TINYVAL .or. c35ext(i,j,ispec) > TINYVAL .or. c55ext(i,j,ispec) > TINYVAL) then
           anisotropic(ispec) = .true.
           poroelastic(ispec) = .false.
           elastic(ispec) = .true.
           any_elastic = .true.
           QKappa_attenuationext(i,j,ispec) = 10.d0
           Qmu_attenuationext(i,j,ispec) = 10.d0
-        elseif(vsext(i,j,ispec) < TINYVAL) then
+        else if(vsext(i,j,ispec) < TINYVAL) then
           elastic(ispec) = .false.
           poroelastic(ispec) = .false.
           any_acoustic = .true.