adds support for model SEMUCB_A3d_3dQ (SEMUCB w/ 3D Berkeley attenuation model)

Author: danielpeter

src/auxiliaries/rules.mk

@@ -586,6 +586,7 @@ xwrite_profile_OBJECTS += \
 	$O/model_aniso_inner_core.check.o \
 	$O/model_aniso_mantle.check.o \
 	$O/model_atten3D_QRFSI12.check.o \
+	$O/model_atten3D_berkeley.check.o \
 	$O/model_attenuation_gll.check.o \
 	$O/model_attenuation.check.o \
 	$O/model_berkeley.check.o \

src/auxiliaries/write_profile.f90

@@ -687,6 +687,7 @@ subroutine write_profile_setup()
     if (ATTENUATION) then
       write(IMAIN,*) '  incorporating attenuation using ',N_SLS,' standard linear solids'
       if (ATTENUATION_3D) write(IMAIN,*) '  using 3D attenuation model'
+      if (ATTENUATION_3D_BERKELEY) write(IMAIN,*) '  using 3D Berkeley attenuation model'
     else
       write(IMAIN,*) '  no attenuation'
     endif
@@ -999,10 +1000,11 @@ subroutine write_profile_model_values(r,r_prem,theta,phi,iregion_code,idoubling,
   !
   !note:  only Qmu attenuation considered, Qkappa attenuation not used so far...
   if (ATTENUATION) then
-    call meshfem3D_models_getatten_val(idoubling,r_prem,theta,phi, &
+    call meshfem3D_models_getatten_val(iregion_code,idoubling, &
+                                       r_prem,theta,phi, &
                                        ispec, i, j, k, &
                                        tau_e,tau_s, &
-                                       moho,Qmu,Qkappa,elem_in_crust)
+                                       moho,Qmu,Qkappa,elem_in_crust,rho)
   endif
 !> end GET_MODEL
 

src/meshfem3D/get_model.F90

@@ -47,7 +47,7 @@ subroutine get_model(iregion_code,ispec,nspec,idoubling, &
 
   use meshfem_models_par, only: &
     ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE, &
-    ATTENUATION,ATTENUATION_3D,ATTENUATION_1D_WITH_3D_STORAGE, &
+    ATTENUATION,ATTENUATION_3D,ATTENUATION_3D_BERKELEY,ATTENUATION_1D_WITH_3D_STORAGE, &
     CEM_ACCEPT,CRUSTAL
 
   use regions_mesh_par2, only: &
@@ -239,10 +239,11 @@ subroutine get_model(iregion_code,ispec,nspec,idoubling, &
         !
         !note:  only Qmu attenuation considered, Qkappa attenuation not used so far...
         if (ATTENUATION) then
-          call meshfem3D_models_getatten_val(idoubling,r_prem,theta,phi, &
+          call meshfem3D_models_getatten_val(iregion_code,idoubling, &
+                                             r_prem,theta,phi, &
                                              ispec, i, j, k, &
                                              tau_e,tau_s, &
-                                             moho,Qmu,Qkappa,elem_in_crust)
+                                             moho,Qmu,Qkappa,elem_in_crust,rho)
         endif
 
         ! define elastic parameters in the model
@@ -318,13 +319,12 @@ subroutine get_model(iregion_code,ispec,nspec,idoubling, &
 
         ! stores attenuation arrays
         if (ATTENUATION) then
-          if (ATTENUATION_3D .or. ATTENUATION_1D_WITH_3D_STORAGE) then
+          if (ATTENUATION_3D .or. ATTENUATION_3D_BERKELEY .or. ATTENUATION_1D_WITH_3D_STORAGE) then
             ! distinguish between single and double precision for reals
             do i_sls = 1,N_SLS
               tau_e_store(i,j,k,i_sls,ispec) = real(tau_e(i_sls), kind=CUSTOM_REAL)
             enddo
             Qmu_store(i,j,k,ispec) = real(Qmu, kind=CUSTOM_REAL)
-
           else
             ! single node per element
             ! distinguish between single and double precision for reals
@@ -335,7 +335,6 @@ subroutine get_model(iregion_code,ispec,nspec,idoubling, &
               enddo
               Qmu_store(1,1,1,ispec) = real(Qmu, kind=CUSTOM_REAL)
             endif
-
           endif
         endif
 

src/meshfem3D/meshfem3D_models.F90

@@ -81,6 +81,9 @@ subroutine meshfem3D_models_broadcast()
     else if (ATTENUATION_3D) then
       ! Colleen's model defined originally between 24.4km and 650km
       call model_atten3D_QRFSI12_broadcast()
+    else if (ATTENUATION_3D_BERKELEY) then
+      ! Berkeley 3D attenuation model
+      call model_atten3D_berkeley_broadcast()
     else
       ! sets up attenuation coefficients according to the chosen, "pure" 1D model
       ! (including their 1D-crustal profiles)
@@ -1647,10 +1650,12 @@ end subroutine meshfem3D_model_crust
 !-------------------------------------------------------------------------------------------------
 !
 
-  subroutine meshfem3D_models_getatten_val(idoubling,r_prem,theta,phi, &
+  subroutine meshfem3D_models_getatten_val(iregion_code,idoubling, &
+                                           r_prem,theta,phi, &
                                            ispec, i, j, k, &
                                            tau_e,tau_s, &
-                                           moho,Qmu,Qkappa,elem_in_crust)
+                                           moho,Qmu,Qkappa, &
+                                           elem_in_crust,rho)
 
 ! sets attenuation values tau_e and Qmu for a given point
 !
@@ -1667,7 +1672,7 @@ subroutine meshfem3D_models_getatten_val(idoubling,r_prem,theta,phi, &
 
   implicit none
 
-  integer,intent(in) :: idoubling
+  integer,intent(in) :: iregion_code,idoubling
 
   double precision,intent(in) :: r_prem
   double precision,intent(in) :: theta,phi
@@ -1681,11 +1686,11 @@ subroutine meshfem3D_models_getatten_val(idoubling,r_prem,theta,phi, &
   double precision, dimension(N_SLS),intent(inout) :: tau_s, tau_e
 
   logical,intent(in) :: elem_in_crust
+  double precision, intent(in) :: rho
 
   ! local parameters
   double precision :: theta_degrees,phi_degrees
   double precision :: r_used
-
   ! geographical values
   double precision :: dist, theta_c, phi_c, dist_c, edge, sponge
 
@@ -1705,6 +1710,7 @@ subroutine meshfem3D_models_getatten_val(idoubling,r_prem,theta,phi, &
     call model_attenuation_gll(ispec, i, j, k, Qmu)
 
   else if (ATTENUATION_3D) then
+    ! 3D attenuation model
     ! used for models: s362ani_3DQ, s362iso_3DQ, 3D_attenuation, SPiRal
 
     ! gets spherical coordinates in degrees
@@ -1713,7 +1719,7 @@ subroutine meshfem3D_models_getatten_val(idoubling,r_prem,theta,phi, &
 
     ! in case models incorporate a 3D crust, attenuation values for mantle
     ! get expanded up to surface, and for the crustal points Qmu for PREM crust is imposed
-    r_used = r_prem*R_PLANET_KM
+    r_used = r_prem * R_PLANET_KM
     if (CRUSTAL) then
       if (r_prem > (ONE-moho) .or. elem_in_crust) then
         ! points in actual crust: puts point radius into prem crust
@@ -1728,8 +1734,19 @@ subroutine meshfem3D_models_getatten_val(idoubling,r_prem,theta,phi, &
     ! gets value according to radius/theta/phi location and idoubling flag
     call model_atten3D_QRFSI12(r_used,theta_degrees,phi_degrees,Qmu,idoubling)
 
-  else
+  else if (ATTENUATION_3D_BERKELEY) then
+    ! 3D Berkeley attenuation model
+    if (elem_in_crust) then
+      ! fixes Q for crust
+      Qmu = 300.0d0
+      Qkappa = 99900.d0  !  not used so far...
+    else
+      r_used = r_prem
+      call model_atten3D_berkeley(r_used,theta,phi,rho,Qkappa,Qmu,iregion_code,CRUSTAL)
+    endif
 
+  else
+    ! default 1D reference (attenuation) model
     select case (REFERENCE_1D_MODEL)
 
       ! case (REFERENCE_MODEL_PREM,REFERENCE_MODEL_PREM2)

src/meshfem3D/meshfem3D_par.f90

@@ -88,7 +88,7 @@ module meshfem_models_par
     MODEL_GLL, &
     HONOR_1D_SPHERICAL_MOHO,CRUSTAL,ONE_CRUST,CASE_3D,TRANSVERSE_ISOTROPY, &
     MODEL_3D_MANTLE_PERTUBATIONS,HETEROGEN_3D_MANTLE,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE, &
-    ATTENUATION_3D, ATTENUATION_GLL, &
+    ATTENUATION_3D,ATTENUATION_3D_BERKELEY,ATTENUATION_GLL, &
     CEM_REQUEST, CEM_ACCEPT, &
     EMC_MODEL,EMC_MODEL_QMU, &
     NX_BATHY,NY_BATHY, &

src/meshfem3D/model_1dberkeley.f90

@@ -63,7 +63,7 @@ module model_1dberkeley_par
     Mref_V_Qmu_berkeley
 
   ! Berkeley 1D model
-  character (len=100) :: berkeley_file_model1D = trim(A3d_folder) // 'model1D.dat'
+  character(len=*), parameter :: berkeley_file_model1D = trim(A3d_folder) // 'model1D.dat'
 
   ! moho layer index
   integer :: moho1D_layer_index = -1