Skip to content

Commit

Permalink
added files for 3D aniso models
Browse files Browse the repository at this point in the history
  • Loading branch information
@bvanderbeek
bvanderbeek committed Jun 30, 2024
1 parent 46dd8d6 commit 8cd3e14
Show file tree
Hide file tree
Showing 6 changed files with 1,622 additions and 4 deletions.
1 change: 1 addition & 0 deletions src/auxiliaries/rules.mk
View file
Original file line number Diff line number Diff line change
Expand Up @@ -584,6 +584,7 @@ xwrite_profile_OBJECTS += \
$O/model_sea1d.check.o \
$O/model_aniso_inner_core.check.o \
$O/model_aniso_mantle.check.o \
$O/model_aniso_mantle_drex.check.o \
$O/model_atten3D_QRFSI12.check.o \
$O/model_attenuation_gll.check.o \
$O/model_attenuation.check.o \
Expand Down
76 changes: 73 additions & 3 deletions src/meshfem3D/meshfem3D_models.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -220,6 +220,10 @@ subroutine meshfem3D_mantle_broadcast()
! Montagner anisotropic model
call model_aniso_mantle_broadcast()

case (THREE_D_MODEL_ANISO_MANTLE_DREX)
! Drex anisotropic mantle model
call model_aniso_mantle_drex_broadcast()

case (THREE_D_MODEL_BKMNS_GLAD)
! GLAD model expansion on block-mantle-spherical-harmonics
if (myrank == 0) write(IMAIN,*) 'setting up both s362ani and bkmns models...'
Expand Down Expand Up @@ -626,6 +630,7 @@ subroutine meshfem3D_models_get3Dmntl_val(iregion_code,r_prem,rho, &
ispec,i,j,k)

use meshfem_models_par
use constants, only: EARTH_R,EARTH_RHOAV,GRAV,PI

implicit none

Expand All @@ -645,8 +650,8 @@ subroutine meshfem3D_models_get3Dmntl_val(iregion_code,r_prem,rho, &
integer, intent(in) :: ispec, i, j, k

! local parameters
double precision :: r_used
double precision :: dvp,dvs,drho,vp,vs,moho,sediment
double precision :: scaleval,scale_GPa,r_used,r_dummy,min_lon,max_lon,min_colat,max_colat,lon_step,colat_step
double precision :: dvp,dvs,drho,vp,vs,moho,sediment,Nparam,Lparam,E_Chi,PwaveMod,SwaveMod
double precision :: dvpv,dvph,dvsv,dvsh,deta
double precision :: lat,lon
double precision :: A,C,L,N,F
Expand Down Expand Up @@ -894,6 +899,71 @@ subroutine meshfem3D_models_get3Dmntl_val(iregion_code,r_prem,rho, &
c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26, &
c33,c34,c35,c36,c44,c45,c46,c55,c56,c66)

case (THREE_D_MODEL_ANISO_MANTLE_DREX)
! Drex anisotropic model

! Following Elodie's advice 04/09/2020
if (r_prem>RCMB/EARTH_R .and. r_prem<RMOHO/EARTH_R) then !!ELODIE
! 14/07/2020 Rappisi ---> if you don't use r_used=r_prem you get problems at discontinuities
r_used = r_prem
else
r_used=RMOHO/EARTH_R
endif

call model_aniso_mantle_drex(r_used,theta,phi,&
rho,c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26, &
c33,c34,c35,c36,c44,c45,c46,c55,c56,c66)

! calculate velocities for visualization (must be in local
! (radial) coordinates system
vpv = sqrt(c33/rho)
vph = sqrt(c11/rho)
vsv = sqrt(c44/rho)
vsh = sqrt(0.5d0*((c11-c12)/rho))
eta_aniso = c13/(c11 - 2.d0*c44)

! 14/10/2020 just for the plot of isotropic velocities otherwise
! the code plots vsv
vp = sqrt( ((8.d0+4.d0*eta_aniso)*vph*vph + 3.d0*vpv*vpv &
+ (8.d0 - 8.d0*eta_aniso)*vsv*vsv)/15.d0 )
vs = sqrt( ((1.d0-2.d0*eta_aniso)*vph*vph + vpv*vpv &
+ 5.d0*vsh*vsh + (6.d0+4.d0*eta_aniso)*vsv*vsv)/15.d0)

! activate these 5 lines if you want to plot the isotropic velocities
!vpv = vp
!vph = vp
!vsv = vs
!vsh = vs
!eta_aniso = 1.d0

Nparam = (1.d0/8.d0)*(c11+c22)-(1.d0/4.d0)*(c12)+(1.d0/2.d0)*(c66)
Lparam = (1.d0/2.d0)*(c44+c55)
scaleval = dsqrt(PI*GRAV*EARTH_RHOAV)
E_Chi = (Nparam/Lparam)
! put vpv=E_Chi if you want to plot radial anisotropy instead of vpv
! vpv = E_Chi*1000.0d0/(scaleval*R_EARTH)

! 9/03/2021 rotate from radial (local) to global, needed to calculate
! seismograms; NB: this rotation is needed only if the tensor is
! in radial system
call rotate_tensor_radial_to_global(theta,phi,c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26, &
c33,c34,c35,c36,c44,c45,c46,c55,c56,c66, &
c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26, &
c33,c34,c35,c36,c44,c45,c46,c55,c56,c66)


!04/06/2020 - Find isotropic P and S waves moduli
PwaveMod = (3.d0/15.d0)*(c11+c22+c33)+(2.d0/15.d0)*(c23+c13+c12)+(4.d0/15.d0)*(c44+c55+c66)
SwaveMod = (1.d0/15.d0)*(c11+c22+c33)-(1.d0/15.d0)*(c23+c13+c12)+(3.d0/15.d0)*(c44+c55+c66)

! calculate isotropic velocities from anisotropic tensor if necessary
!vph = sqrt(PwaveMod/rho)
!vpv = vph
!vsh = sqrt(SwaveMod/rho)
!vsv = vsh
!eta_aniso = 1.d0


case (THREE_D_MODEL_BKMNS_GLAD)
! GLAD model expansion on block-mantle-spherical-harmonics
! model takes actual position (includes topography between 80km depth and surface topo)
Expand Down Expand Up @@ -1011,7 +1081,7 @@ subroutine meshfem3D_models_get3Dmntl_val(iregion_code,r_prem,rho, &

! special case for fully anisotropic models
! anisotropic models defined between the Moho and 670 km (change to CMB if desired)
if (THREE_D_MODEL == THREE_D_MODEL_ANISO_MANTLE) then
if (THREE_D_MODEL == THREE_D_MODEL_ANISO_MANTLE .or. THREE_D_MODEL == THREE_D_MODEL_ANISO_MANTLE_DREX) then
! special case for model_aniso_mantle()
if (suppress_mantle_extension) then
! point is in crust, and CRUSTAL == .false. (no 3D crustal model); model_aniso_mantle() hasn't been called.
Expand Down
2 changes: 1 addition & 1 deletion src/meshfem3D/meshfem3D_par.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,7 @@ module meshfem_models_par
THREE_D_MODEL_SEA99,THREE_D_MODEL_JP3D,THREE_D_MODEL_JP3D, &
THREE_D_MODEL_S362ANI,THREE_D_MODEL_S362WMANI,THREE_D_MODEL_S362ANI_PREM,THREE_D_MODEL_S29EA, &
THREE_D_MODEL_PPM,THREE_D_MODEL_GAPP2,THREE_D_MODEL_SGLOBE,THREE_D_MODEL_SGLOBE_ISO, &
THREE_D_MODEL_ANISO_MANTLE,THREE_D_MODEL_GLL, &
THREE_D_MODEL_ANISO_MANTLE,THREE_D_MODEL_ANISO_MANTLE_DREX,THREE_D_MODEL_GLL, &
THREE_D_MODEL_BKMNS_GLAD, &
THREE_D_MODEL_SPIRAL,THREE_D_MODEL_HETEROGEN_PREM, &
THREE_D_MODEL_SH_MARS, &
Expand Down
Loading

0 comments on commit 8cd3e14

Please sign in to comment.