updates tomo model to be usable with additional Qp and Qs values and …

…so for gll model to read in also Q values; updates attenuation routines and inserts a criteria from Anderson & Hart (1978) to limit Qs/Qp ratios for bulk attenuation and Olson scaling

src/generate_databases/create_regions_mesh.f90

@@ -343,7 +343,8 @@ subroutine create_regions_mesh()
   deallocate(xixstore,xiystore,xizstore,&
              etaxstore,etaystore,etazstore,&
              gammaxstore,gammaystore,gammazstore)
-  deallocate(jacobianstore,qkappa_attenuation_store,qmu_attenuation_store)
+  deallocate(jacobianstore)
+  deallocate(qkappa_attenuation_store,qmu_attenuation_store)
   deallocate(kappastore,mustore,rho_vp,rho_vs)
   deallocate(rho_vpI,rho_vpII,rho_vsI)
   deallocate(rhoarraystore,kappaarraystore,etastore,phistore,tortstore,permstore)
@@ -403,6 +404,7 @@ subroutine crm_ext_allocate_arrays(nspec,LOCAL_PATH,myrank, &
   allocate(xelm(NGNOD),yelm(NGNOD),zelm(NGNOD),stat=ier)
   if (ier /= 0) stop 'error allocating array xelm etc.'
 
+! attenuation
   allocate(qkappa_attenuation_store(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
   if (ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays')
 

src/generate_databases/get_model.f90

@@ -433,6 +433,12 @@ subroutine get_model_values(materials_ext_mesh,nmat_ext_mesh, &
   iflag_aniso = 0
   idomain_id = IDOMAIN_ELASTIC
 
+  ! attenuation
+  ! shear attenuation: arbitrary value, see maximum in constants.h
+  qmu_atten = ATTENUATION_COMP_MAXIMUM
+  ! bulk attenuation: arbitrary (high) default value
+  qkappa_atten = 9999.0_CUSTOM_REAL
+
   ! selects chosen velocity model
   select case (IMODEL)
 
@@ -450,33 +456,27 @@ subroutine get_model_values(materials_ext_mesh,nmat_ext_mesh, &
   case (IMODEL_1D_PREM )
     ! 1D model profile from PREM
     call model_1D_prem_iso(xmesh,ymesh,zmesh,rho,vp,vs,qmu_atten)
-    qkappa_atten = 9999.  ! undefined in this model
 
   case (IMODEL_1D_PREM_PB )
     ! 1D model profile from PREM modified by Piero
     imaterial_PB = abs(imaterial_id)
     call model_1D_PREM_routine_PB(xmesh,ymesh,zmesh,rho,vp,vs,imaterial_PB)
-    ! attenuation: arbitrary value, see maximum in constants.h
-    qmu_atten = ATTENUATION_COMP_MAXIMUM
+
     ! sets acoustic/elastic domain as given in materials properties
     iundef = - imaterial_id    ! iundef must be positive
     read(undef_mat_prop(6,iundef),*) idomain_id
-    qkappa_atten = 9999.  ! undefined in this model
 
   case (IMODEL_1D_CASCADIA )
     ! 1D model profile for Cascadia region
     call model_1D_cascadia(xmesh,ymesh,zmesh,rho,vp,vs,qmu_atten)
-    qkappa_atten = 9999.  ! undefined in this model
 
   case (IMODEL_1D_SOCAL )
     ! 1D model profile for Southern California
     call model_1D_socal(xmesh,ymesh,zmesh,rho,vp,vs,qmu_atten)
-    qkappa_atten = 9999.  ! undefined in this model
 
   case (IMODEL_SALTON_TROUGH )
     ! gets model values from tomography file
     call model_salton_trough(xmesh,ymesh,zmesh,rho,vp,vs,qmu_atten)
-    qkappa_atten = 9999.  ! undefined in this model
 
   case (IMODEL_TOMO )
     ! gets model values from tomography file

src/generate_databases/model_default.f90

@@ -35,14 +35,15 @@ subroutine model_default(materials_ext_mesh,nmat_ext_mesh, &
                           undef_mat_prop,nundefMat_ext_mesh, &
                           imaterial_id,imaterial_def, &
                           xmesh,ymesh,zmesh, &
-                          rho,vp,vs,iflag_aniso,qkappa_atten,qmu_atten,idomain_id, &
+                          rho,vp,vs, &
+                          iflag_aniso,qkappa_atten,qmu_atten,idomain_id, &
                           rho_s,kappa_s,rho_f,kappa_f,eta_f,kappa_fr,mu_fr, &
                           phi,tort,kxx,kxy,kxz,kyy,kyz,kzz)
 
 ! takes model values specified by mesh properties
 
   use generate_databases_par, only: IDOMAIN_ACOUSTIC,IDOMAIN_ELASTIC,IDOMAIN_POROELASTIC
-  use create_regions_mesh_ext_par
+  use create_regions_mesh_ext_par, only: CUSTOM_REAL,MAX_STRING_LEN
 
   implicit none
 
@@ -71,6 +72,8 @@ subroutine model_default(materials_ext_mesh,nmat_ext_mesh, &
   character(len=MAX_STRING_LEN) :: str_domain
   integer :: ier
 
+  !print *,'model defaults: ',imaterial_id,imaterial_def
+
   ! check if the material is known or unknown
   if (imaterial_id > 0) then
     ! gets velocity model as specified by (cubit) mesh files for elastic & acoustic

src/generate_databases/model_gll.f90

@@ -38,8 +38,12 @@
 
   subroutine model_gll(myrank,nspec,LOCAL_PATH)
 
-  use generate_databases_par,only: NGLLX,NGLLY,NGLLZ,FOUR_THIRDS,IMAIN
-  use create_regions_mesh_ext_par
+  use generate_databases_par,only: NGLLX,NGLLY,NGLLZ,FOUR_THIRDS,IMAIN,MAX_STRING_LEN, &
+    ATTENUATION,FULL_ATTENUATION_SOLID
+
+  use create_regions_mesh_ext_par,only: rhostore,kappastore,mustore,rho_vp,rho_vs, &
+    qkappa_attenuation_store,qmu_attenuation_store
+
   implicit none
 
   integer, intent(in) :: myrank,nspec
@@ -63,9 +67,12 @@ subroutine model_gll(myrank,nspec,LOCAL_PATH)
   !!! use lines for vp only
 
   ! density
-  allocate( rho_read(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
+  allocate(rho_read(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
   if (ier /= 0) stop 'error allocating array rho_read'
 
+  ! user output
+  if (myrank==0) write(IMAIN,*) '     reading in: rho.bin'
+
   filename = prname_lp(1:len_trim(prname_lp))//'rho.bin'
   open(unit=28,file=trim(filename),status='old',action='read',form='unformatted',iostat=ier)
   if (ier /= 0) then
@@ -77,9 +84,12 @@ subroutine model_gll(myrank,nspec,LOCAL_PATH)
   close(28)
 
   ! vp
-  allocate( vp_read(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
+  allocate(vp_read(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
   if (ier /= 0) stop 'error allocating array vp_read'
 
+  ! user output
+  if (myrank==0) write(IMAIN,*) '     reading in: vp.bin'
+
   filename = prname_lp(1:len_trim(prname_lp))//'vp.bin'
   open(unit=28,file=trim(filename),status='old',action='read',form='unformatted',iostat=ier)
   if (ier /= 0) then
@@ -91,9 +101,12 @@ subroutine model_gll(myrank,nspec,LOCAL_PATH)
   close(28)
 
   ! vs
-  allocate( vs_read(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
+  allocate(vs_read(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
   if (ier /= 0) stop 'error allocating array vs_read'
 
+  ! user output
+  if (myrank==0) write(IMAIN,*) '     reading in: vs.bin'
+
   filename = prname_lp(1:len_trim(prname_lp))//'vs.bin'
   open(unit=28,file=trim(filename),status='old',action='read',form='unformatted',iostat=ier)
   if (ier /= 0) then
@@ -127,11 +140,46 @@ subroutine model_gll(myrank,nspec,LOCAL_PATH)
   !!! update arrays that will be saved and used in the solver xspecfem3D
   !!! the following part is neccessary if you uncommented something above
 
-  rhostore    = rho_read
-  kappastore  = rhostore * ( vp_read * vp_read - FOUR_THIRDS * vs_read * vs_read )
-  mustore     = rhostore * vs_read * vs_read
-  rho_vp = rhostore * vp_read
-  rho_vs = rhostore * vs_read
+  rhostore(:,:,:,:) = rho_read(:,:,:,:)
+  kappastore(:,:,:,:) = rhostore(:,:,:,:) * ( vp_read(:,:,:,:) * vp_read(:,:,:,:) &
+                                              - FOUR_THIRDS * vs_read(:,:,:,:) * vs_read(:,:,:,:) )
+  mustore(:,:,:,:) = rhostore(:,:,:,:) * vs_read(:,:,:,:) * vs_read(:,:,:,:)
+
+  rho_vp(:,:,:,:) = rhostore(:,:,:,:) * vp_read(:,:,:,:)
+  rho_vs(:,:,:,:) = rhostore(:,:,:,:) * vs_read(:,:,:,:)
+
+  ! gets attenuation arrays from files
+  if (ATTENUATION) then
+    ! shear attenuation
+    ! user output
+    if (myrank==0) write(IMAIN,*) '     reading in: qmu.bin'
+
+    filename = prname_lp(1:len_trim(prname_lp))//'qmu.bin'
+    open(unit=28,file=trim(filename),status='old',action='read',form='unformatted',iostat=ier)
+    if (ier /= 0) then
+      print *,'Error opening file: ',trim(filename)
+      stop 'Error reading qmu.bin file'
+    endif
+
+    read(28) qmu_attenuation_store
+    close(28)
+
+    ! bulk attenuation
+    if (FULL_ATTENUATION_SOLID) then
+      ! user output
+      if (myrank==0) write(IMAIN,*) '     reading in: qkappa.bin'
+
+      filename = prname_lp(1:len_trim(prname_lp))//'qkappa.bin'
+      open(unit=28,file=trim(filename),status='old',action='read',form='unformatted',iostat=ier)
+      if (ier /= 0) then
+        print *,'error opening file: ',trim(filename)
+        stop 'error reading qkappa.bin file'
+      endif
+
+      read(28) qkappa_attenuation_store
+      close(28)
+    endif
+  endif
 
   ! free memory
   deallocate( rho_read,vp_read,vs_read)

src/generate_databases/model_gll_adios.F90

@@ -38,7 +38,10 @@
 subroutine model_gll_adios(myrank,nspec,LOCAL_PATH)
 
   use adios_read_mod
-  use generate_databases_par,only: NGLLX,NGLLY,NGLLZ,FOUR_THIRDS,IMAIN
+
+  use generate_databases_par,only: NGLLX,NGLLY,NGLLZ,FOUR_THIRDS,IMAIN, &
+    ATTENUATION,FULL_ATTENUATION_SOLID
+
   use create_regions_mesh_ext_par
 
   implicit none
@@ -94,18 +97,33 @@ subroutine model_gll_adios(myrank,nspec,LOCAL_PATH)
   start(1) = local_dim_rho * myrank
   count_ad(1) = NGLLX * NGLLY * NGLLZ * nspec
   call adios_selection_boundingbox(sel, 1, start, count_ad)
+
+  ! wave speeds
   call adios_schedule_read(handle, sel, "rho/array", 0, 1, &
                            rho_read, ier)
   call adios_schedule_read(handle, sel, "vp/array", 0, 1, &
                            vp_read, ier)
   call adios_schedule_read(handle, sel, "vp/array", 0, 1, &
                            vp_read, ier)
 
+  ! attenuation
+  if (ATTENUATION) then
+    ! shear attenuation
+    call adios_schedule_read(handle, sel, "qmu/array", 0, 1, &
+                             qmu_attenuation_store, ier)
+    ! bulk attenuation
+    if (FULL_ATTENUATION_SOLID) then
+      call adios_schedule_read(handle, sel, "qkappa/array", 0, 1, &
+                               qkappa_attenuation_store, ier)
+    endif
+  endif
+
   !---------------------------------------.
   ! Perform read and close the adios file |
   !---------------------------------------'
   call adios_perform_reads(handle, ier)
   if (ier /= 0) call stop_all()
+
   call adios_read_close(handle,ier)
   call adios_read_finalize_method(ADIOS_READ_METHOD_BP, ier)
 
@@ -132,11 +150,12 @@ subroutine model_gll_adios(myrank,nspec,LOCAL_PATH)
   !!! update arrays that will be saved and used in the solver xspecfem3D
   !!! the following part is neccessary if you uncommented something above
 
-  rhostore    = rho_read
-  kappastore  = rhostore * ( vp_read * vp_read - FOUR_THIRDS * vs_read * vs_read )
-  mustore     = rhostore * vs_read * vs_read
-  rho_vp = rhostore * vp_read
-  rho_vs = rhostore * vs_read
+  rhostore(:,:,:,:) = rho_read(:,:,:,:)
+  kappastore(:,:,:,:) = rhostore(:,:,:,:) * ( vp_read(:,:,:,:) * vp_read(:,:,:,:) &
+                          - FOUR_THIRDS * vs_read(:,:,:,:) * vs_read(:,:,:,:) )
+  mustore(:,:,:,:) = rhostore(:,:,:,:) * vs_read(:,:,:,:) * vs_read(:,:,:,:)
+  rho_vp(:,:,:,:) = rhostore(:,:,:,:) * vp_read(:,:,:,:)
+  rho_vs(:,:,:,:) = rhostore(:,:,:,:) * vs_read(:,:,:,:)
 
   ! free memory
   deallocate( rho_read,vp_read,vs_read)