supports reading (optional) attenuation Q-values (Qp & Qs) in tomogra…

…phy file

EXAMPLES/tomographic_ocean_model/README

@@ -7,6 +7,14 @@ with topography variations.
 
 The tomographic file will be created by script createTomographyFile.py provided in utils/ directory.
 
+The default tomography file will provide Vp, Vs and density on a data line like:
+  #x #z #vp #vs #rho
+Optionally, attenuation Q-values can be added to these data lines, like:
+  #x #z #vp #vs #rho #Qp #Qs
+to specify a 2D velocity and attenuation model.
+You could for example modify the script createTomographyFile.py accordingly.
+
+
 To run this example, type:
    ./run_this_example.sh
 

setup/constants.h.in

@@ -121,6 +121,9 @@
 ! for viscoacousticity
   logical, parameter :: USE_A_STRONG_FORMULATION_FOR_E1 = .false.
 
+! maximum attenuation Q-value
+  double precision, parameter :: ATTENUATION_COMP_MAXIMUM    = 9999.d0
+
 !!-----------------------------------------------------------
 !!
 !! noise simulations

src/shared/read_material_table.f90

@@ -35,7 +35,8 @@ subroutine read_material_table()
 
 ! reads in material definitions in DATA/Par_file
 
-  use constants, only: IMAIN,TINYVAL,ISOTROPIC_MATERIAL,ANISOTROPIC_MATERIAL,POROELASTIC_MATERIAL
+  use constants, only: IMAIN,TINYVAL,ISOTROPIC_MATERIAL,ANISOTROPIC_MATERIAL,POROELASTIC_MATERIAL, &
+    ATTENUATION_COMP_MAXIMUM
 
   use shared_parameters, only: AXISYM,nbmodels,icodemat,cp,cs, &
                               aniso3,aniso4,aniso5,aniso6,aniso7,aniso8,aniso9,aniso10,aniso11,aniso12, &
@@ -116,8 +117,8 @@ subroutine read_material_table()
   aniso11(:) = 0.d0
 
   comp_g(:) = 0.0d0
-  QKappa(:) = 9999.d0
-  Qmu(:) = 9999.d0
+  QKappa(:) = ATTENUATION_COMP_MAXIMUM
+  Qmu(:) = ATTENUATION_COMP_MAXIMUM
 
   rho_s_read(:) = 0.d0
   rho_f_read(:) = 0.d0

src/specfem2D/define_external_model.f90

@@ -35,7 +35,7 @@ subroutine define_external_model_dummy(coord,material_element,ibool, &
                                          rho,vp,vs,QKappa_attenuation,Qmu_attenuation, &
                                          c11,c12,c13,c15,c23,c25,c33,c35,c55)
 
-  use constants, only: CUSTOM_REAL,NGLLX,NGLLZ,NDIM,IMAIN
+  use constants, only: CUSTOM_REAL,NGLLX,NGLLZ,NDIM,IMAIN,ATTENUATION_COMP_MAXIMUM
 
   use specfem_par, only: myrank,nspec,nglob
 
@@ -75,8 +75,8 @@ subroutine define_external_model_dummy(coord,material_element,ibool, &
   endif
 
   ! default no attenuation
-  QKappa_attenuation(:,:,:) = 9999.d0
-  Qmu_attenuation(:,:,:) = 9999.d0
+  QKappa_attenuation(:,:,:) = ATTENUATION_COMP_MAXIMUM
+  Qmu_attenuation(:,:,:) = ATTENUATION_COMP_MAXIMUM
 
 ! loop on all the elements of the mesh, and inside each element loop on all the GLL points
   do ispec = 1,nspec
@@ -92,8 +92,8 @@ subroutine define_external_model_dummy(coord,material_element,ibool, &
           rho(i,j,ispec) = 2000.d0 + 10.d0*x - 32.17d0*z**3
           vp(i,j,ispec) = 3000.d0 - 7.12d0*cos(x)
           vs(i,j,ispec) = vp(i,j,ispec) / sqrt(3.d0)
-          QKappa_attenuation(i,j,ispec) = 9999. ! this means no attenuation
-          Qmu_attenuation(i,j,ispec)    = 9999. ! this means no attenuation
+          QKappa_attenuation(i,j,ispec) = ATTENUATION_COMP_MAXIMUM ! this means no attenuation
+          Qmu_attenuation(i,j,ispec)    = ATTENUATION_COMP_MAXIMUM ! this means no attenuation
           c11(i,j,ispec) = 169.d9
           c13(i,j,ispec) = 122.d9
           c15(i,j,ispec) = 0.d0
@@ -142,7 +142,7 @@ subroutine define_external_model(coord,material_element,ibool, &
                                    rho,vp,vs,QKappa_attenuation,Qmu_attenuation, &
                                    c11,c12,c13,c15,c23,c25,c33,c35,c55)
 
-  use constants, only: CUSTOM_REAL,NGLLX,NGLLZ,NDIM,IMAIN
+  use constants, only: CUSTOM_REAL,NGLLX,NGLLZ,NDIM,IMAIN,ATTENUATION_COMP_MAXIMUM
 
   use specfem_par, only: myrank,nspec,nglob
 
@@ -214,8 +214,8 @@ subroutine define_external_model(coord,material_element,ibool, &
   endif
 
   ! default no attenuation
-  QKappa_attenuation(:,:,:) = 9999.d0
-  Qmu_attenuation(:,:,:) = 9999.d0
+  QKappa_attenuation(:,:,:) = ATTENUATION_COMP_MAXIMUM
+  Qmu_attenuation(:,:,:) = ATTENUATION_COMP_MAXIMUM
 
   ! default no anisotropy
   c11(:,:,:) = 0.d0
@@ -1118,8 +1118,8 @@ subroutine define_external_model(coord,material_element,ibool, &
         ! also set fictitious attenuation to a very high value (attenuation is not used in the fluid)
         if (material_element(ispec) == IREGION_OUTER_CORE) then
           vs(i,j,ispec) = 0.d0
-          Qkappa_attenuation(i,j,ispec) = 9999.d0
-          Qmu_attenuation(i,j,ispec) = 9999.d0
+          Qkappa_attenuation(i,j,ispec) = ATTENUATION_COMP_MAXIMUM
+          Qmu_attenuation(i,j,ispec) = ATTENUATION_COMP_MAXIMUM
         endif
 
       enddo

src/specfem2D/define_external_model_from_marmousi.f90

@@ -37,7 +37,7 @@ subroutine define_external_model_from_marmousi(coord,ibool,rho,vp,vs,QKappa_atte
 ! reads in external model using a marmousi format which defines a compaction gradient
 
 
-  use constants, only: CUSTOM_REAL,NGLLX,NGLLZ,NDIM,IMAIN
+  use constants, only: CUSTOM_REAL,NGLLX,NGLLZ,NDIM,IMAIN,ATTENUATION_COMP_MAXIMUM
 
   use specfem_par, only: poroelastcoef,density,kmato,myrank,nspec,nglob
 
@@ -66,8 +66,8 @@ subroutine define_external_model_from_marmousi(coord,ibool,rho,vp,vs,QKappa_atte
   endif
 
   ! default no attenuation
-  QKappa_attenuation(:,:,:) = 9999.d0
-  Qmu_attenuation(:,:,:) = 9999.d0
+  QKappa_attenuation(:,:,:) = ATTENUATION_COMP_MAXIMUM
+  Qmu_attenuation(:,:,:) = ATTENUATION_COMP_MAXIMUM
 
 ! loop on all the elements of the mesh, and inside each element loop on all the GLL points
   do ispec = 1,nspec
@@ -93,8 +93,8 @@ subroutine define_external_model_from_marmousi(coord,ibool,rho,vp,vs,QKappa_atte
   ! assumes Poisson solids
         vs(i,j,ispec) = vp(i,j,ispec) / sqrt(3.d0)
 
-        QKappa_attenuation(i,j,ispec) = 9999. ! this means no attenuation
-        Qmu_attenuation(i,j,ispec)    = 9999. ! this means no attenuation
+        QKappa_attenuation(i,j,ispec) = ATTENUATION_COMP_MAXIMUM ! this means no attenuation
+        Qmu_attenuation(i,j,ispec)    = ATTENUATION_COMP_MAXIMUM ! this means no attenuation
 
         c11(i,j,ispec) = 0.d0   ! this means no anisotropy
         c13(i,j,ispec) = 0.d0