[cig-commits] r16127 - seismo/3D/SPECFEM3D_GLOBE/trunk
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Thu Jan 7 06:04:23 PST 2010
Author: danielpeter
Date: 2010-01-07 06:04:22 -0800 (Thu, 07 Jan 2010)
New Revision: 16127
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_crust_mantle.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_inner_core.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90
Log:
added Deville et al. (2002) routines for elastic solver in specfem3D.f90, compute_forces_crust_mantle.f90 and compute_forces_inner_core.f90 (only valid for NGLLX=NGLLY=NGLLZ=5)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_crust_mantle.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_crust_mantle.f90 2010-01-06 23:41:13 UTC (rev 16126)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_crust_mantle.f90 2010-01-07 14:04:22 UTC (rev 16127)
@@ -849,3 +849,937 @@
end subroutine compute_forces_crust_mantle
+
+
+
+!
+!---------------------------------------------------------------------------------------------------------------------
+!
+
+ subroutine compute_forces_crust_mantle_Dev(minus_gravity_table,density_table,minus_deriv_gravity_table, &
+ displ,accel,xstore,ystore,zstore, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ hprime_xx,hprime_xxT, &
+ hprimewgll_xx,hprimewgll_xxT, &
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
+ kappavstore,kappahstore,muvstore,muhstore,eta_anisostore, &
+ c11store,c12store,c13store,c14store,c15store,c16store,c22store, &
+ c23store,c24store,c25store,c26store,c33store,c34store,c35store, &
+ c36store,c44store,c45store,c46store,c55store,c56store,c66store, &
+ ibool,idoubling,R_memory,epsilondev,epsilon_trace_over_3,one_minus_sum_beta, &
+ alphaval,betaval,gammaval,factor_common,vx,vy,vz,vnspec,COMPUTE_AND_STORE_STRAIN, AM_V)
+
+ implicit none
+
+ include "constants.h"
+
+! include values created by the mesher
+! done for performance only using static allocation to allow for loop unrolling
+ include "OUTPUT_FILES/values_from_mesher.h"
+
+! attenuation_model_variables
+ type attenuation_model_variables
+ sequence
+ double precision min_period, max_period
+ double precision :: QT_c_source ! Source Frequency
+ double precision, dimension(:), pointer :: Qtau_s ! tau_sigma
+ double precision, dimension(:), pointer :: QrDisc ! Discontinutitues Defined
+ double precision, dimension(:), pointer :: Qr ! Radius
+ integer, dimension(:), pointer :: interval_Q ! Steps
+ double precision, dimension(:), pointer :: Qmu ! Shear Attenuation
+ double precision, dimension(:,:), pointer :: Qtau_e ! tau_epsilon
+ double precision, dimension(:), pointer :: Qomsb, Qomsb2 ! one_minus_sum_beta
+ double precision, dimension(:,:), pointer :: Qfc, Qfc2 ! factor_common
+ double precision, dimension(:), pointer :: Qsf, Qsf2 ! scale_factor
+ integer, dimension(:), pointer :: Qrmin ! Max and Mins of idoubling
+ integer, dimension(:), pointer :: Qrmax ! Max and Mins of idoubling
+ integer :: Qn ! Number of points
+ end type attenuation_model_variables
+
+ type (attenuation_model_variables) AM_V
+! attenuation_model_variables
+
+! for forward or backward simulations
+ logical COMPUTE_AND_STORE_STRAIN
+
+! array with the local to global mapping per slice
+ integer, dimension(NSPEC_CRUST_MANTLE) :: idoubling
+
+! displacement and acceleration
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CRUST_MANTLE) :: displ,accel
+
+! memory variables for attenuation
+! memory variables R_ij are stored at the local rather than global level
+! to allow for optimization of cache access by compiler
+ integer i_sls,i_memory
+! variable sized array variables for one_minus_sum_beta and factor_common
+ integer vx, vy, vz, vnspec
+
+ real(kind=CUSTOM_REAL) one_minus_sum_beta_use,minus_sum_beta
+ real(kind=CUSTOM_REAL), dimension(vx, vy, vz, vnspec) :: one_minus_sum_beta
+
+ integer iregion_selected
+
+! for attenuation
+ real(kind=CUSTOM_REAL) R_xx_val,R_yy_val
+ real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: R_memory
+ real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: epsilondev
+
+! [alpha,beta,gamma]val reduced to N_SLS and factor_common to N_SLS*NUM_NODES
+ real(kind=CUSTOM_REAL), dimension(N_SLS) :: alphaval,betaval,gammaval
+ real(kind=CUSTOM_REAL), dimension(N_SLS, vx, vy, vz, vnspec) :: factor_common
+ real(kind=CUSTOM_REAL), dimension(NGLLX, NGLLY, NGLLZ) :: factor_common_c44_muv
+
+ real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc
+ real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: epsilon_trace_over_3
+
+! arrays with mesh parameters per slice
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: ibool
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
+
+! array with derivatives of Lagrange polynomials and precalculated products
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx,hprimewgll_xx
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY) :: wgllwgll_xy
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ) :: wgllwgll_xz
+ real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLZ) :: wgllwgll_yz
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
+ tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3
+
+! x y and z contain r theta and phi
+ real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE) :: xstore,ystore,zstore
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPECMAX_ISO_MANTLE) :: &
+ kappavstore,muvstore
+
+! store anisotropic properties only where needed to save memory
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPECMAX_TISO_MANTLE) :: &
+ kappahstore,muhstore,eta_anisostore
+
+! arrays for full anisotropy only when needed
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPECMAX_ANISO_MANTLE) :: &
+ c11store,c12store,c13store,c14store,c15store,c16store,c22store, &
+ c23store,c24store,c25store,c26store,c33store,c34store,c35store, &
+ c36store,c44store,c45store,c46store,c55store,c56store,c66store
+
+ integer ispec,iglob,ispec_strain
+ integer i,j,k !,l
+
+! the 21 coefficients for an anisotropic medium in reduced notation
+ real(kind=CUSTOM_REAL) c11,c22,c33,c44,c55,c66,c12,c13,c23,c14,c24,c34,c15,c25,c35,c45,c16,c26,c36,c46,c56
+
+ real(kind=CUSTOM_REAL) rhovphsq,sinphifour,cosphisq,sinphisq,costhetasq,rhovsvsq,sinthetasq, &
+ cosphifour,costhetafour,rhovpvsq,sinthetafour,rhovshsq,cosfourphi, &
+ costwotheta,cosfourtheta,sintwophisq,costheta,sinphi,sintheta,cosphi, &
+ sintwotheta,costwophi,sintwophi,costwothetasq,costwophisq,phi,theta
+
+ real(kind=CUSTOM_REAL) two_rhovpvsq,two_rhovphsq,two_rhovsvsq,two_rhovshsq
+ real(kind=CUSTOM_REAL) four_rhovpvsq,four_rhovphsq,four_rhovsvsq,four_rhovshsq
+
+ real(kind=CUSTOM_REAL) twoetaminone,etaminone,eta_aniso
+ real(kind=CUSTOM_REAL) two_eta_aniso,four_eta_aniso,six_eta_aniso
+
+ real(kind=CUSTOM_REAL) xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
+ real(kind=CUSTOM_REAL) duxdxl,duxdyl,duxdzl,duydxl,duydyl,duydzl,duzdxl,duzdyl,duzdzl
+
+ real(kind=CUSTOM_REAL) duxdxl_plus_duydyl,duxdxl_plus_duzdzl,duydyl_plus_duzdzl
+ real(kind=CUSTOM_REAL) duxdyl_plus_duydxl,duzdxl_plus_duxdzl,duzdyl_plus_duydzl
+
+ real(kind=CUSTOM_REAL) sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz
+
+ !real(kind=CUSTOM_REAL) hp1,hp2,hp3
+ real(kind=CUSTOM_REAL) fac1,fac2,fac3
+ real(kind=CUSTOM_REAL) lambdal,mul,lambdalplus2mul
+ real(kind=CUSTOM_REAL) kappal,kappavl,kappahl,muvl,muhl
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NGLLZ) :: sum_terms
+
+! for gravity
+ integer int_radius
+ real(kind=CUSTOM_REAL) sigma_yx,sigma_zx,sigma_zy
+ real(kind=CUSTOM_REAL) radius_cr
+ double precision radius,rho,minus_g,minus_dg
+ double precision minus_g_over_radius,minus_dg_plus_g_over_radius
+ double precision cos_theta,sin_theta,cos_phi,sin_phi
+ double precision cos_theta_sq,sin_theta_sq,cos_phi_sq,sin_phi_sq
+ double precision factor,sx_l,sy_l,sz_l,gxl,gyl,gzl
+ double precision Hxxl,Hyyl,Hzzl,Hxyl,Hxzl,Hyzl
+ double precision, dimension(NRAD_GRAVITY) :: minus_gravity_table,density_table,minus_deriv_gravity_table
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NGLLZ) :: rho_s_H
+ double precision, dimension(NGLLX,NGLLY,NGLLZ) :: wgll_cube
+
+
+! Deville
+! manually inline the calls to the Deville et al. (2002) routines
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: dummyx_loc,dummyy_loc,dummyz_loc, &
+ newtempx1,newtempx2,newtempx3,newtempy1,newtempy2,newtempy3,newtempz1,newtempz2,newtempz3
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xxT,hprimewgll_xxT
+
+ integer, parameter :: m1 = NGLLX, m2 = NGLLX * NGLLY
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,m2) :: B1_m1_m2_5points,B2_m1_m2_5points,B3_m1_m2_5points
+ real(kind=CUSTOM_REAL), dimension(m1,m2) :: C1_m1_m2_5points,C2_m1_m2_5points,C3_m1_m2_5points
+ real(kind=CUSTOM_REAL), dimension(m1,m2) :: E1_m1_m2_5points,E2_m1_m2_5points,E3_m1_m2_5points
+
+ equivalence(dummyx_loc,B1_m1_m2_5points)
+ equivalence(dummyy_loc,B2_m1_m2_5points)
+ equivalence(dummyz_loc,B3_m1_m2_5points)
+ equivalence(tempx1,C1_m1_m2_5points)
+ equivalence(tempy1,C2_m1_m2_5points)
+ equivalence(tempz1,C3_m1_m2_5points)
+ equivalence(newtempx1,E1_m1_m2_5points)
+ equivalence(newtempy1,E2_m1_m2_5points)
+ equivalence(newtempz1,E3_m1_m2_5points)
+
+ real(kind=CUSTOM_REAL), dimension(m2,NGLLX) :: A1_mxm_m2_m1_5points,A2_mxm_m2_m1_5points,A3_mxm_m2_m1_5points
+ real(kind=CUSTOM_REAL), dimension(m2,m1) :: C1_mxm_m2_m1_5points,C2_mxm_m2_m1_5points,C3_mxm_m2_m1_5points
+ real(kind=CUSTOM_REAL), dimension(m2,m1) :: E1_mxm_m2_m1_5points,E2_mxm_m2_m1_5points,E3_mxm_m2_m1_5points
+
+ equivalence(dummyx_loc,A1_mxm_m2_m1_5points)
+ equivalence(dummyy_loc,A2_mxm_m2_m1_5points)
+ equivalence(dummyz_loc,A3_mxm_m2_m1_5points)
+ equivalence(tempx3,C1_mxm_m2_m1_5points)
+ equivalence(tempy3,C2_mxm_m2_m1_5points)
+ equivalence(tempz3,C3_mxm_m2_m1_5points)
+ equivalence(newtempx3,E1_mxm_m2_m1_5points)
+ equivalence(newtempy3,E2_mxm_m2_m1_5points)
+ equivalence(newtempz3,E3_mxm_m2_m1_5points)
+
+! ****************************************************
+! big loop over all spectral elements in the solid
+! ****************************************************
+
+! set acceleration to zero
+ accel(:,:) = 0._CUSTOM_REAL
+
+ do ispec = 1,NSPEC_CRUST_MANTLE
+
+ ! subroutines adapted from Deville, Fischer and Mund, High-order methods
+ ! for incompressible fluid flow, Cambridge University Press (2002),
+ ! pages 386 and 389 and Figure 8.3.1
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+ dummyx_loc(i,j,k) = displ(1,iglob)
+ dummyy_loc(i,j,k) = displ(2,iglob)
+ dummyz_loc(i,j,k) = displ(3,iglob)
+ enddo
+ enddo
+ enddo
+ do j=1,m2
+ do i=1,m1
+ C1_m1_m2_5points(i,j) = hprime_xx(i,1)*B1_m1_m2_5points(1,j) + &
+ hprime_xx(i,2)*B1_m1_m2_5points(2,j) + &
+ hprime_xx(i,3)*B1_m1_m2_5points(3,j) + &
+ hprime_xx(i,4)*B1_m1_m2_5points(4,j) + &
+ hprime_xx(i,5)*B1_m1_m2_5points(5,j)
+
+ C2_m1_m2_5points(i,j) = hprime_xx(i,1)*B2_m1_m2_5points(1,j) + &
+ hprime_xx(i,2)*B2_m1_m2_5points(2,j) + &
+ hprime_xx(i,3)*B2_m1_m2_5points(3,j) + &
+ hprime_xx(i,4)*B2_m1_m2_5points(4,j) + &
+ hprime_xx(i,5)*B2_m1_m2_5points(5,j)
+
+ C3_m1_m2_5points(i,j) = hprime_xx(i,1)*B3_m1_m2_5points(1,j) + &
+ hprime_xx(i,2)*B3_m1_m2_5points(2,j) + &
+ hprime_xx(i,3)*B3_m1_m2_5points(3,j) + &
+ hprime_xx(i,4)*B3_m1_m2_5points(4,j) + &
+ hprime_xx(i,5)*B3_m1_m2_5points(5,j)
+ enddo
+ enddo
+ do j=1,m1
+ do i=1,m1
+ ! for efficiency it is better to leave this loop on k inside, it leads to slightly faster code
+ do k = 1,NGLLX
+ tempx2(i,j,k) = dummyx_loc(i,1,k)*hprime_xxT(1,j) + &
+ dummyx_loc(i,2,k)*hprime_xxT(2,j) + &
+ dummyx_loc(i,3,k)*hprime_xxT(3,j) + &
+ dummyx_loc(i,4,k)*hprime_xxT(4,j) + &
+ dummyx_loc(i,5,k)*hprime_xxT(5,j)
+
+ tempy2(i,j,k) = dummyy_loc(i,1,k)*hprime_xxT(1,j) + &
+ dummyy_loc(i,2,k)*hprime_xxT(2,j) + &
+ dummyy_loc(i,3,k)*hprime_xxT(3,j) + &
+ dummyy_loc(i,4,k)*hprime_xxT(4,j) + &
+ dummyy_loc(i,5,k)*hprime_xxT(5,j)
+
+ tempz2(i,j,k) = dummyz_loc(i,1,k)*hprime_xxT(1,j) + &
+ dummyz_loc(i,2,k)*hprime_xxT(2,j) + &
+ dummyz_loc(i,3,k)*hprime_xxT(3,j) + &
+ dummyz_loc(i,4,k)*hprime_xxT(4,j) + &
+ dummyz_loc(i,5,k)*hprime_xxT(5,j)
+ enddo
+ enddo
+ enddo
+ do j=1,m1
+ do i=1,m2
+ C1_mxm_m2_m1_5points(i,j) = A1_mxm_m2_m1_5points(i,1)*hprime_xxT(1,j) + &
+ A1_mxm_m2_m1_5points(i,2)*hprime_xxT(2,j) + &
+ A1_mxm_m2_m1_5points(i,3)*hprime_xxT(3,j) + &
+ A1_mxm_m2_m1_5points(i,4)*hprime_xxT(4,j) + &
+ A1_mxm_m2_m1_5points(i,5)*hprime_xxT(5,j)
+
+ C2_mxm_m2_m1_5points(i,j) = A2_mxm_m2_m1_5points(i,1)*hprime_xxT(1,j) + &
+ A2_mxm_m2_m1_5points(i,2)*hprime_xxT(2,j) + &
+ A2_mxm_m2_m1_5points(i,3)*hprime_xxT(3,j) + &
+ A2_mxm_m2_m1_5points(i,4)*hprime_xxT(4,j) + &
+ A2_mxm_m2_m1_5points(i,5)*hprime_xxT(5,j)
+
+ C3_mxm_m2_m1_5points(i,j) = A3_mxm_m2_m1_5points(i,1)*hprime_xxT(1,j) + &
+ A3_mxm_m2_m1_5points(i,2)*hprime_xxT(2,j) + &
+ A3_mxm_m2_m1_5points(i,3)*hprime_xxT(3,j) + &
+ A3_mxm_m2_m1_5points(i,4)*hprime_xxT(4,j) + &
+ A3_mxm_m2_m1_5points(i,5)*hprime_xxT(5,j)
+ enddo
+ enddo
+
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ ! get derivatives of ux, uy and uz with respect to x, y and z
+ xixl = xix(i,j,k,ispec)
+ xiyl = xiy(i,j,k,ispec)
+ xizl = xiz(i,j,k,ispec)
+ etaxl = etax(i,j,k,ispec)
+ etayl = etay(i,j,k,ispec)
+ etazl = etaz(i,j,k,ispec)
+ gammaxl = gammax(i,j,k,ispec)
+ gammayl = gammay(i,j,k,ispec)
+ gammazl = gammaz(i,j,k,ispec)
+
+ ! compute the jacobian
+ jacobianl = 1._CUSTOM_REAL / (xixl*(etayl*gammazl-etazl*gammayl) &
+ - xiyl*(etaxl*gammazl-etazl*gammaxl) &
+ + xizl*(etaxl*gammayl-etayl*gammaxl))
+
+ duxdxl = xixl*tempx1(i,j,k) + etaxl*tempx2(i,j,k) + gammaxl*tempx3(i,j,k)
+ duxdyl = xiyl*tempx1(i,j,k) + etayl*tempx2(i,j,k) + gammayl*tempx3(i,j,k)
+ duxdzl = xizl*tempx1(i,j,k) + etazl*tempx2(i,j,k) + gammazl*tempx3(i,j,k)
+
+ duydxl = xixl*tempy1(i,j,k) + etaxl*tempy2(i,j,k) + gammaxl*tempy3(i,j,k)
+ duydyl = xiyl*tempy1(i,j,k) + etayl*tempy2(i,j,k) + gammayl*tempy3(i,j,k)
+ duydzl = xizl*tempy1(i,j,k) + etazl*tempy2(i,j,k) + gammazl*tempy3(i,j,k)
+
+ duzdxl = xixl*tempz1(i,j,k) + etaxl*tempz2(i,j,k) + gammaxl*tempz3(i,j,k)
+ duzdyl = xiyl*tempz1(i,j,k) + etayl*tempz2(i,j,k) + gammayl*tempz3(i,j,k)
+ duzdzl = xizl*tempz1(i,j,k) + etazl*tempz2(i,j,k) + gammazl*tempz3(i,j,k)
+
+ ! precompute some sums to save CPU time
+ duxdxl_plus_duydyl = duxdxl + duydyl
+ duxdxl_plus_duzdzl = duxdxl + duzdzl
+ duydyl_plus_duzdzl = duydyl + duzdzl
+ duxdyl_plus_duydxl = duxdyl + duydxl
+ duzdxl_plus_duxdzl = duzdxl + duxdzl
+ duzdyl_plus_duydzl = duzdyl + duydzl
+
+ ! compute deviatoric strain
+ if (COMPUTE_AND_STORE_STRAIN) then
+ if(NSPEC_CRUST_MANTLE_STRAIN_ONLY == 1) then
+ ispec_strain = 1
+ else
+ ispec_strain = ispec
+ endif
+ epsilon_trace_over_3(i,j,k,ispec_strain) = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ epsilondev_loc(1,i,j,k) = duxdxl - epsilon_trace_over_3(i,j,k,ispec_strain)
+ epsilondev_loc(2,i,j,k) = duydyl - epsilon_trace_over_3(i,j,k,ispec_strain)
+ epsilondev_loc(3,i,j,k) = 0.5 * duxdyl_plus_duydxl
+ epsilondev_loc(4,i,j,k) = 0.5 * duzdxl_plus_duxdzl
+ epsilondev_loc(5,i,j,k) = 0.5 * duzdyl_plus_duydzl
+ endif
+
+ ! precompute terms for attenuation if needed
+ if(ATTENUATION_VAL) then
+ if(ATTENUATION_3D_VAL) then
+ one_minus_sum_beta_use = one_minus_sum_beta(i,j,k,ispec)
+ else
+ radius_cr = xstore(ibool(i,j,k,ispec))
+ call get_attenuation_index(idoubling(ispec), dble(radius_cr), iregion_selected, .FALSE., AM_V)
+ one_minus_sum_beta_use = one_minus_sum_beta(1,1,1,iregion_selected)
+ endif
+ minus_sum_beta = one_minus_sum_beta_use - 1.0
+ endif
+
+ !
+ ! compute either isotropic or anisotropic elements
+ !
+ if(ANISOTROPIC_3D_MANTLE_VAL) then
+
+ c11 = c11store(i,j,k,ispec)
+ c12 = c12store(i,j,k,ispec)
+ c13 = c13store(i,j,k,ispec)
+ c14 = c14store(i,j,k,ispec)
+ c15 = c15store(i,j,k,ispec)
+ c16 = c16store(i,j,k,ispec)
+ c22 = c22store(i,j,k,ispec)
+ c23 = c23store(i,j,k,ispec)
+ c24 = c24store(i,j,k,ispec)
+ c25 = c25store(i,j,k,ispec)
+ c26 = c26store(i,j,k,ispec)
+ c33 = c33store(i,j,k,ispec)
+ c34 = c34store(i,j,k,ispec)
+ c35 = c35store(i,j,k,ispec)
+ c36 = c36store(i,j,k,ispec)
+ c44 = c44store(i,j,k,ispec)
+ c45 = c45store(i,j,k,ispec)
+ c46 = c46store(i,j,k,ispec)
+ c55 = c55store(i,j,k,ispec)
+ c56 = c56store(i,j,k,ispec)
+ c66 = c66store(i,j,k,ispec)
+
+ if(ATTENUATION_VAL) then
+ mul = c44
+ c11 = c11 + FOUR_THIRDS * minus_sum_beta * mul
+ c12 = c12 - TWO_THIRDS * minus_sum_beta * mul
+ c13 = c13 - TWO_THIRDS * minus_sum_beta * mul
+ c22 = c22 + FOUR_THIRDS * minus_sum_beta * mul
+ c23 = c23 - TWO_THIRDS * minus_sum_beta * mul
+ c33 = c33 + FOUR_THIRDS * minus_sum_beta * mul
+ c44 = c44 + minus_sum_beta * mul
+ c55 = c55 + minus_sum_beta * mul
+ c66 = c66 + minus_sum_beta * mul
+ endif
+
+ sigma_xx = c11*duxdxl + c16*duxdyl_plus_duydxl + c12*duydyl + &
+ c15*duzdxl_plus_duxdzl + c14*duzdyl_plus_duydzl + c13*duzdzl
+
+ sigma_yy = c12*duxdxl + c26*duxdyl_plus_duydxl + c22*duydyl + &
+ c25*duzdxl_plus_duxdzl + c24*duzdyl_plus_duydzl + c23*duzdzl
+
+ sigma_zz = c13*duxdxl + c36*duxdyl_plus_duydxl + c23*duydyl + &
+ c35*duzdxl_plus_duxdzl + c34*duzdyl_plus_duydzl + c33*duzdzl
+
+ sigma_xy = c16*duxdxl + c66*duxdyl_plus_duydxl + c26*duydyl + &
+ c56*duzdxl_plus_duxdzl + c46*duzdyl_plus_duydzl + c36*duzdzl
+
+ sigma_xz = c15*duxdxl + c56*duxdyl_plus_duydxl + c25*duydyl + &
+ c55*duzdxl_plus_duxdzl + c45*duzdyl_plus_duydzl + c35*duzdzl
+
+ sigma_yz = c14*duxdxl + c46*duxdyl_plus_duydxl + c24*duydyl + &
+ c45*duzdxl_plus_duxdzl + c44*duzdyl_plus_duydzl + c34*duzdzl
+
+ else
+
+ ! do not use transverse isotropy except if element is between d220 and Moho
+ if(.not. (TRANSVERSE_ISOTROPY_VAL .and. (idoubling(ispec)==IFLAG_220_80 .or. idoubling(ispec)==IFLAG_80_MOHO))) then
+
+ ! layer with no transverse isotropy, use kappav and muv
+ kappal = kappavstore(i,j,k,ispec)
+ mul = muvstore(i,j,k,ispec)
+
+ ! use unrelaxed parameters if attenuation
+ if(ATTENUATION_VAL) mul = mul * one_minus_sum_beta_use
+
+ lambdalplus2mul = kappal + FOUR_THIRDS * mul
+ lambdal = lambdalplus2mul - 2.*mul
+
+ ! compute stress sigma
+ sigma_xx = lambdalplus2mul*duxdxl + lambdal*duydyl_plus_duzdzl
+ sigma_yy = lambdalplus2mul*duydyl + lambdal*duxdxl_plus_duzdzl
+ sigma_zz = lambdalplus2mul*duzdzl + lambdal*duxdxl_plus_duydyl
+
+ sigma_xy = mul*duxdyl_plus_duydxl
+ sigma_xz = mul*duzdxl_plus_duxdzl
+ sigma_yz = mul*duzdyl_plus_duydzl
+
+ else
+
+ ! use Kappa and mu from transversely isotropic model
+ kappavl = kappavstore(i,j,k,ispec)
+ muvl = muvstore(i,j,k,ispec)
+
+ kappahl = kappahstore(i,j,k,ispec)
+ muhl = muhstore(i,j,k,ispec)
+
+ ! use unrelaxed parameters if attenuation
+ ! eta does not need to be shifted since it is a ratio
+ if(ATTENUATION_VAL) then
+ muvl = muvl * one_minus_sum_beta_use
+ muhl = muhl * one_minus_sum_beta_use
+ endif
+
+ rhovpvsq = kappavl + FOUR_THIRDS * muvl !!! that is C
+ rhovphsq = kappahl + FOUR_THIRDS * muhl !!! that is A
+
+ rhovsvsq = muvl !!! that is L
+ rhovshsq = muhl !!! that is N
+
+ eta_aniso = eta_anisostore(i,j,k,ispec) !!! that is F / (A - 2 L)
+
+ ! use mesh coordinates to get theta and phi
+ ! ystore and zstore contain theta and phi
+
+ iglob = ibool(i,j,k,ispec)
+ theta = ystore(iglob)
+ phi = zstore(iglob)
+
+ costheta = cos(theta)
+ sintheta = sin(theta)
+ cosphi = cos(phi)
+ sinphi = sin(phi)
+
+ costhetasq = costheta * costheta
+ sinthetasq = sintheta * sintheta
+ cosphisq = cosphi * cosphi
+ sinphisq = sinphi * sinphi
+
+ costhetafour = costhetasq * costhetasq
+ sinthetafour = sinthetasq * sinthetasq
+ cosphifour = cosphisq * cosphisq
+ sinphifour = sinphisq * sinphisq
+
+ costwotheta = cos(2.*theta)
+ sintwotheta = sin(2.*theta)
+ costwophi = cos(2.*phi)
+ sintwophi = sin(2.*phi)
+
+ cosfourtheta = cos(4.*theta)
+ cosfourphi = cos(4.*phi)
+
+ costwothetasq = costwotheta * costwotheta
+
+ costwophisq = costwophi * costwophi
+ sintwophisq = sintwophi * sintwophi
+
+ etaminone = eta_aniso - 1.
+ twoetaminone = 2. * eta_aniso - 1.
+
+ ! precompute some products to reduce the CPU time
+ two_eta_aniso = 2.*eta_aniso
+ four_eta_aniso = 4.*eta_aniso
+ six_eta_aniso = 6.*eta_aniso
+
+ two_rhovpvsq = 2.*rhovpvsq
+ two_rhovphsq = 2.*rhovphsq
+ two_rhovsvsq = 2.*rhovsvsq
+ two_rhovshsq = 2.*rhovshsq
+
+ four_rhovpvsq = 4.*rhovpvsq
+ four_rhovphsq = 4.*rhovphsq
+ four_rhovsvsq = 4.*rhovsvsq
+ four_rhovshsq = 4.*rhovshsq
+
+ ! the 21 anisotropic coefficients computed using Mathematica
+
+ c11 = rhovphsq*sinphifour + 2.*cosphisq*sinphisq* &
+ (rhovphsq*costhetasq + (eta_aniso*rhovphsq + two_rhovsvsq - two_eta_aniso*rhovsvsq)* &
+ sinthetasq) + cosphifour* &
+ (rhovphsq*costhetafour + 2.*(eta_aniso*rhovphsq + two_rhovsvsq - two_eta_aniso*rhovsvsq)* &
+ costhetasq*sinthetasq + rhovpvsq*sinthetafour)
+
+ c12 = ((rhovphsq - two_rhovshsq)*(3. + cosfourphi)*costhetasq)/4. - &
+ four_rhovshsq*cosphisq*costhetasq*sinphisq + &
+ (rhovphsq*(11. + 4.*costwotheta + cosfourtheta)*sintwophisq)/32. + &
+ eta_aniso*(rhovphsq - two_rhovsvsq)*(cosphifour + &
+ 2.*cosphisq*costhetasq*sinphisq + sinphifour)*sinthetasq + &
+ rhovpvsq*cosphisq*sinphisq*sinthetafour - &
+ rhovsvsq*sintwophisq*sinthetafour
+
+ c13 = (cosphisq*(rhovphsq + six_eta_aniso*rhovphsq + rhovpvsq - four_rhovsvsq - &
+ 12.*eta_aniso*rhovsvsq + (twoetaminone*rhovphsq - rhovpvsq + four_rhovsvsq - &
+ four_eta_aniso*rhovsvsq)*cosfourtheta))/8. + &
+ sinphisq*(eta_aniso*(rhovphsq - two_rhovsvsq)*costhetasq + &
+ (rhovphsq - two_rhovshsq)*sinthetasq)
+
+ c14 = costheta*sinphi*((cosphisq* &
+ (-rhovphsq + rhovpvsq + four_rhovshsq - four_rhovsvsq + &
+ (-rhovphsq + two_eta_aniso*rhovphsq - rhovpvsq + four_rhovsvsq - &
+ four_eta_aniso*rhovsvsq)*costwotheta))/2. + &
+ (etaminone*rhovphsq + 2.*(rhovshsq - eta_aniso*rhovsvsq))*sinphisq)* sintheta
+
+ c15 = cosphi*costheta*((cosphisq* (-rhovphsq + rhovpvsq + &
+ (twoetaminone*rhovphsq - rhovpvsq + four_rhovsvsq - four_eta_aniso*rhovsvsq)* &
+ costwotheta))/2. + etaminone*(rhovphsq - two_rhovsvsq)*sinphisq)*sintheta
+
+ c16 = (cosphi*sinphi*(cosphisq* (-rhovphsq + rhovpvsq + &
+ (-rhovphsq + two_eta_aniso*rhovphsq - rhovpvsq + four_rhovsvsq - &
+ four_eta_aniso*rhovsvsq)*costwotheta) + &
+ 2.*etaminone*(rhovphsq - two_rhovsvsq)*sinphisq)*sinthetasq)/2.
+
+ c22 = rhovphsq*cosphifour + 2.*cosphisq*sinphisq* &
+ (rhovphsq*costhetasq + (eta_aniso*rhovphsq + two_rhovsvsq - two_eta_aniso*rhovsvsq)* &
+ sinthetasq) + sinphifour* &
+ (rhovphsq*costhetafour + 2.*(eta_aniso*rhovphsq + two_rhovsvsq - two_eta_aniso*rhovsvsq)* &
+ costhetasq*sinthetasq + rhovpvsq*sinthetafour)
+
+ c23 = ((rhovphsq + six_eta_aniso*rhovphsq + rhovpvsq - four_rhovsvsq - 12.*eta_aniso*rhovsvsq + &
+ (twoetaminone*rhovphsq - rhovpvsq + four_rhovsvsq - four_eta_aniso*rhovsvsq)* &
+ cosfourtheta)*sinphisq)/8. + &
+ cosphisq*(eta_aniso*(rhovphsq - two_rhovsvsq)*costhetasq + &
+ (rhovphsq - two_rhovshsq)*sinthetasq)
+
+ c24 = costheta*sinphi*(etaminone*(rhovphsq - two_rhovsvsq)*cosphisq + &
+ ((-rhovphsq + rhovpvsq + (twoetaminone*rhovphsq - rhovpvsq + &
+ four_rhovsvsq - four_eta_aniso*rhovsvsq)*costwotheta)*sinphisq)/2.)*sintheta
+
+ c25 = cosphi*costheta*((etaminone*rhovphsq + 2.*(rhovshsq - eta_aniso*rhovsvsq))* &
+ cosphisq + ((-rhovphsq + rhovpvsq + four_rhovshsq - four_rhovsvsq + &
+ (-rhovphsq + two_eta_aniso*rhovphsq - rhovpvsq + four_rhovsvsq - &
+ four_eta_aniso*rhovsvsq)*costwotheta)*sinphisq)/2.)*sintheta
+
+ c26 = (cosphi*sinphi*(2.*etaminone*(rhovphsq - two_rhovsvsq)*cosphisq + &
+ (-rhovphsq + rhovpvsq + (-rhovphsq + two_eta_aniso*rhovphsq - rhovpvsq + four_rhovsvsq - &
+ four_eta_aniso*rhovsvsq)*costwotheta)*sinphisq)*sinthetasq)/2.
+
+ c33 = rhovpvsq*costhetafour + 2.*(eta_aniso*(rhovphsq - two_rhovsvsq) + two_rhovsvsq)* &
+ costhetasq*sinthetasq + rhovphsq*sinthetafour
+
+ c34 = -((rhovphsq - rhovpvsq + (twoetaminone*rhovphsq - rhovpvsq + four_rhovsvsq &
+ - four_eta_aniso*rhovsvsq)*costwotheta)*sinphi*sintwotheta)/4.
+
+ c35 = -(cosphi*(rhovphsq - rhovpvsq + &
+ (twoetaminone*rhovphsq - rhovpvsq + four_rhovsvsq - four_eta_aniso*rhovsvsq)* &
+ costwotheta)*sintwotheta)/4.
+
+ c36 = -((rhovphsq - rhovpvsq - four_rhovshsq + four_rhovsvsq + &
+ (twoetaminone*rhovphsq - rhovpvsq + four_rhovsvsq - four_eta_aniso*rhovsvsq)* &
+ costwotheta)*sintwophi*sinthetasq)/4.
+
+ c44 = cosphisq*(rhovsvsq*costhetasq + rhovshsq*sinthetasq) + &
+ sinphisq*(rhovsvsq*costwothetasq + &
+ (rhovphsq - two_eta_aniso*rhovphsq + rhovpvsq + four_eta_aniso*rhovsvsq)*costhetasq* sinthetasq)
+
+ c45 = ((rhovphsq - two_eta_aniso*rhovphsq + rhovpvsq - two_rhovshsq - two_rhovsvsq + &
+ four_eta_aniso*rhovsvsq + (rhovphsq - two_eta_aniso*rhovphsq + rhovpvsq + &
+ 4.*etaminone*rhovsvsq)*costwotheta)*sintwophi*sinthetasq)/4.
+
+ c46 = -(cosphi*costheta*((rhovshsq - rhovsvsq)*cosphisq - &
+ ((rhovphsq - two_eta_aniso*rhovphsq + rhovpvsq - two_rhovshsq - two_rhovsvsq + &
+ four_eta_aniso*rhovsvsq + (-rhovphsq + two_eta_aniso*rhovphsq - rhovpvsq + &
+ four_rhovsvsq - four_eta_aniso*rhovsvsq)*costwotheta)*sinphisq)/2.)* sintheta)
+
+ c55 = sinphisq*(rhovsvsq*costhetasq + rhovshsq*sinthetasq) + &
+ cosphisq*(rhovsvsq*costwothetasq + &
+ (rhovphsq - two_eta_aniso*rhovphsq + rhovpvsq + four_eta_aniso*rhovsvsq)*costhetasq* sinthetasq)
+
+ c56 = costheta*sinphi*((cosphisq* &
+ (rhovphsq - two_eta_aniso*rhovphsq + rhovpvsq - two_rhovshsq - two_rhovsvsq + &
+ four_eta_aniso*rhovsvsq + (-rhovphsq + two_eta_aniso*rhovphsq - rhovpvsq + &
+ four_rhovsvsq - four_eta_aniso*rhovsvsq)*costwotheta))/2. + &
+ (-rhovshsq + rhovsvsq)*sinphisq)*sintheta
+
+ c66 = rhovshsq*costwophisq*costhetasq - &
+ 2.*(rhovphsq - two_rhovshsq)*cosphisq*costhetasq*sinphisq + &
+ (rhovphsq*(11. + 4.*costwotheta + cosfourtheta)*sintwophisq)/32. - &
+ (rhovsvsq*(-6. - 2.*cosfourphi + cos(4.*phi - 2.*theta) - 2.*costwotheta + &
+ cos(2.*(2.*phi + theta)))*sinthetasq)/8. + &
+ rhovpvsq*cosphisq*sinphisq*sinthetafour - &
+ (eta_aniso*(rhovphsq - two_rhovsvsq)*sintwophisq*sinthetafour)/2.
+
+ ! general expression of stress tensor for full Cijkl with 21 coefficients
+ sigma_xx = c11*duxdxl + c16*duxdyl_plus_duydxl + c12*duydyl + &
+ c15*duzdxl_plus_duxdzl + c14*duzdyl_plus_duydzl + c13*duzdzl
+
+ sigma_yy = c12*duxdxl + c26*duxdyl_plus_duydxl + c22*duydyl + &
+ c25*duzdxl_plus_duxdzl + c24*duzdyl_plus_duydzl + c23*duzdzl
+
+ sigma_zz = c13*duxdxl + c36*duxdyl_plus_duydxl + c23*duydyl + &
+ c35*duzdxl_plus_duxdzl + c34*duzdyl_plus_duydzl + c33*duzdzl
+
+ sigma_xy = c16*duxdxl + c66*duxdyl_plus_duydxl + c26*duydyl + &
+ c56*duzdxl_plus_duxdzl + c46*duzdyl_plus_duydzl + c36*duzdzl
+
+ sigma_xz = c15*duxdxl + c56*duxdyl_plus_duydxl + c25*duydyl + &
+ c55*duzdxl_plus_duxdzl + c45*duzdyl_plus_duydzl + c35*duzdzl
+
+ sigma_yz = c14*duxdxl + c46*duxdyl_plus_duydxl + c24*duydyl + &
+ c45*duzdxl_plus_duxdzl + c44*duzdyl_plus_duydzl + c34*duzdzl
+
+ endif
+
+ endif ! end of test whether isotropic or anisotropic element
+
+ ! subtract memory variables if attenuation
+ if(ATTENUATION_VAL) then
+ do i_sls = 1,N_SLS
+ R_xx_val = R_memory(1,i_sls,i,j,k,ispec)
+ R_yy_val = R_memory(2,i_sls,i,j,k,ispec)
+ sigma_xx = sigma_xx - R_xx_val
+ sigma_yy = sigma_yy - R_yy_val
+ sigma_zz = sigma_zz + R_xx_val + R_yy_val
+ sigma_xy = sigma_xy - R_memory(3,i_sls,i,j,k,ispec)
+ sigma_xz = sigma_xz - R_memory(4,i_sls,i,j,k,ispec)
+ sigma_yz = sigma_yz - R_memory(5,i_sls,i,j,k,ispec)
+ enddo
+ endif
+
+ ! define symmetric components of sigma for gravity
+ sigma_yx = sigma_xy
+ sigma_zx = sigma_xz
+ sigma_zy = sigma_yz
+
+ ! compute non-symmetric terms for gravity
+ if(GRAVITY_VAL) then
+
+ ! use mesh coordinates to get theta and phi
+ ! x y and z contain r theta and phi
+ iglob = ibool(i,j,k,ispec)
+ radius = dble(xstore(iglob))
+ theta = ystore(iglob)
+ phi = zstore(iglob)
+
+ cos_theta = dcos(dble(theta))
+ sin_theta = dsin(dble(theta))
+ cos_phi = dcos(dble(phi))
+ sin_phi = dsin(dble(phi))
+
+ ! get g, rho and dg/dr=dg
+ ! spherical components of the gravitational acceleration
+ ! for efficiency replace with lookup table every 100 m in radial direction
+ int_radius = nint(radius * R_EARTH_KM * 10.d0)
+ minus_g = minus_gravity_table(int_radius)
+ minus_dg = minus_deriv_gravity_table(int_radius)
+ rho = density_table(int_radius)
+
+ ! Cartesian components of the gravitational acceleration
+ gxl = minus_g*sin_theta*cos_phi
+ gyl = minus_g*sin_theta*sin_phi
+ gzl = minus_g*cos_theta
+
+ ! Cartesian components of gradient of gravitational acceleration
+ ! obtained from spherical components
+ minus_g_over_radius = minus_g / radius
+ minus_dg_plus_g_over_radius = minus_dg - minus_g_over_radius
+
+ cos_theta_sq = cos_theta**2
+ sin_theta_sq = sin_theta**2
+ cos_phi_sq = cos_phi**2
+ sin_phi_sq = sin_phi**2
+
+ Hxxl = minus_g_over_radius*(cos_phi_sq*cos_theta_sq + sin_phi_sq) + cos_phi_sq*minus_dg*sin_theta_sq
+ Hyyl = minus_g_over_radius*(cos_phi_sq + cos_theta_sq*sin_phi_sq) + minus_dg*sin_phi_sq*sin_theta_sq
+ Hzzl = cos_theta_sq*minus_dg + minus_g_over_radius*sin_theta_sq
+ Hxyl = cos_phi*minus_dg_plus_g_over_radius*sin_phi*sin_theta_sq
+ Hxzl = cos_phi*cos_theta*minus_dg_plus_g_over_radius*sin_theta
+ Hyzl = cos_theta*minus_dg_plus_g_over_radius*sin_phi*sin_theta
+
+ iglob = ibool(i,j,k,ispec)
+
+ ! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+
+ ! get displacement and multiply by density to compute G tensor
+ sx_l = rho * dble(displ(1,iglob))
+ sy_l = rho * dble(displ(2,iglob))
+ sz_l = rho * dble(displ(3,iglob))
+
+ ! compute G tensor from s . g and add to sigma (not symmetric)
+ sigma_xx = sigma_xx + sngl(sy_l*gyl + sz_l*gzl)
+ sigma_yy = sigma_yy + sngl(sx_l*gxl + sz_l*gzl)
+ sigma_zz = sigma_zz + sngl(sx_l*gxl + sy_l*gyl)
+
+ sigma_xy = sigma_xy - sngl(sx_l * gyl)
+ sigma_yx = sigma_yx - sngl(sy_l * gxl)
+
+ sigma_xz = sigma_xz - sngl(sx_l * gzl)
+ sigma_zx = sigma_zx - sngl(sz_l * gxl)
+
+ sigma_yz = sigma_yz - sngl(sy_l * gzl)
+ sigma_zy = sigma_zy - sngl(sz_l * gyl)
+
+ ! precompute vector
+ factor = dble(jacobianl) * wgll_cube(i,j,k)
+ rho_s_H(1,i,j,k) = sngl(factor * (sx_l * Hxxl + sy_l * Hxyl + sz_l * Hxzl))
+ rho_s_H(2,i,j,k) = sngl(factor * (sx_l * Hxyl + sy_l * Hyyl + sz_l * Hyzl))
+ rho_s_H(3,i,j,k) = sngl(factor * (sx_l * Hxzl + sy_l * Hyzl + sz_l * Hzzl))
+
+ else
+
+ ! get displacement and multiply by density to compute G tensor
+ sx_l = rho * displ(1,iglob)
+ sy_l = rho * displ(2,iglob)
+ sz_l = rho * displ(3,iglob)
+
+ ! compute G tensor from s . g and add to sigma (not symmetric)
+ sigma_xx = sigma_xx + sy_l*gyl + sz_l*gzl
+ sigma_yy = sigma_yy + sx_l*gxl + sz_l*gzl
+ sigma_zz = sigma_zz + sx_l*gxl + sy_l*gyl
+
+ sigma_xy = sigma_xy - sx_l * gyl
+ sigma_yx = sigma_yx - sy_l * gxl
+
+ sigma_xz = sigma_xz - sx_l * gzl
+ sigma_zx = sigma_zx - sz_l * gxl
+
+ sigma_yz = sigma_yz - sy_l * gzl
+ sigma_zy = sigma_zy - sz_l * gyl
+
+ ! precompute vector
+ factor = jacobianl * wgll_cube(i,j,k)
+ rho_s_H(1,i,j,k) = factor * (sx_l * Hxxl + sy_l * Hxyl + sz_l * Hxzl)
+ rho_s_H(2,i,j,k) = factor * (sx_l * Hxyl + sy_l * Hyyl + sz_l * Hyzl)
+ rho_s_H(3,i,j,k) = factor * (sx_l * Hxzl + sy_l * Hyzl + sz_l * Hzzl)
+
+ endif
+
+ endif ! end of section with gravity terms
+
+ ! form dot product with test vector, non-symmetric form
+ tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_yx*xiyl + sigma_zx*xizl)
+ tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_zy*xizl)
+ tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl)
+
+ tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_yx*etayl + sigma_zx*etazl)
+ tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_zy*etazl)
+ tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl)
+
+ tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_yx*gammayl + sigma_zx*gammazl)
+ tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_zy*gammazl)
+ tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl)
+ enddo ! NGLLX
+ enddo ! NGLLY
+ enddo ! NGLLZ
+
+ ! subroutines adapted from Deville, Fischer and Mund, High-order methods
+ ! for incompressible fluid flow, Cambridge University Press (2002),
+ ! pages 386 and 389 and Figure 8.3.1
+ do j=1,m2
+ do i=1,m1
+ E1_m1_m2_5points(i,j) = hprimewgll_xxT(i,1)*C1_m1_m2_5points(1,j) + &
+ hprimewgll_xxT(i,2)*C1_m1_m2_5points(2,j) + &
+ hprimewgll_xxT(i,3)*C1_m1_m2_5points(3,j) + &
+ hprimewgll_xxT(i,4)*C1_m1_m2_5points(4,j) + &
+ hprimewgll_xxT(i,5)*C1_m1_m2_5points(5,j)
+
+ E2_m1_m2_5points(i,j) = hprimewgll_xxT(i,1)*C2_m1_m2_5points(1,j) + &
+ hprimewgll_xxT(i,2)*C2_m1_m2_5points(2,j) + &
+ hprimewgll_xxT(i,3)*C2_m1_m2_5points(3,j) + &
+ hprimewgll_xxT(i,4)*C2_m1_m2_5points(4,j) + &
+ hprimewgll_xxT(i,5)*C2_m1_m2_5points(5,j)
+
+ E3_m1_m2_5points(i,j) = hprimewgll_xxT(i,1)*C3_m1_m2_5points(1,j) + &
+ hprimewgll_xxT(i,2)*C3_m1_m2_5points(2,j) + &
+ hprimewgll_xxT(i,3)*C3_m1_m2_5points(3,j) + &
+ hprimewgll_xxT(i,4)*C3_m1_m2_5points(4,j) + &
+ hprimewgll_xxT(i,5)*C3_m1_m2_5points(5,j)
+ enddo
+ enddo
+ do i=1,m1
+ do j=1,m1
+ ! for efficiency it is better to leave this loop on k inside, it leads to slightly faster code
+ do k = 1,NGLLX
+ newtempx2(i,j,k) = tempx2(i,1,k)*hprimewgll_xx(1,j) + &
+ tempx2(i,2,k)*hprimewgll_xx(2,j) + &
+ tempx2(i,3,k)*hprimewgll_xx(3,j) + &
+ tempx2(i,4,k)*hprimewgll_xx(4,j) + &
+ tempx2(i,5,k)*hprimewgll_xx(5,j)
+
+ newtempy2(i,j,k) = tempy2(i,1,k)*hprimewgll_xx(1,j) + &
+ tempy2(i,2,k)*hprimewgll_xx(2,j) + &
+ tempy2(i,3,k)*hprimewgll_xx(3,j) + &
+ tempy2(i,4,k)*hprimewgll_xx(4,j) + &
+ tempy2(i,5,k)*hprimewgll_xx(5,j)
+
+ newtempz2(i,j,k) = tempz2(i,1,k)*hprimewgll_xx(1,j) + &
+ tempz2(i,2,k)*hprimewgll_xx(2,j) + &
+ tempz2(i,3,k)*hprimewgll_xx(3,j) + &
+ tempz2(i,4,k)*hprimewgll_xx(4,j) + &
+ tempz2(i,5,k)*hprimewgll_xx(5,j)
+ enddo
+ enddo
+ enddo
+ do j=1,m1
+ do i=1,m2
+ E1_mxm_m2_m1_5points(i,j) = C1_mxm_m2_m1_5points(i,1)*hprimewgll_xx(1,j) + &
+ C1_mxm_m2_m1_5points(i,2)*hprimewgll_xx(2,j) + &
+ C1_mxm_m2_m1_5points(i,3)*hprimewgll_xx(3,j) + &
+ C1_mxm_m2_m1_5points(i,4)*hprimewgll_xx(4,j) + &
+ C1_mxm_m2_m1_5points(i,5)*hprimewgll_xx(5,j)
+
+ E2_mxm_m2_m1_5points(i,j) = C2_mxm_m2_m1_5points(i,1)*hprimewgll_xx(1,j) + &
+ C2_mxm_m2_m1_5points(i,2)*hprimewgll_xx(2,j) + &
+ C2_mxm_m2_m1_5points(i,3)*hprimewgll_xx(3,j) + &
+ C2_mxm_m2_m1_5points(i,4)*hprimewgll_xx(4,j) + &
+ C2_mxm_m2_m1_5points(i,5)*hprimewgll_xx(5,j)
+
+ E3_mxm_m2_m1_5points(i,j) = C3_mxm_m2_m1_5points(i,1)*hprimewgll_xx(1,j) + &
+ C3_mxm_m2_m1_5points(i,2)*hprimewgll_xx(2,j) + &
+ C3_mxm_m2_m1_5points(i,3)*hprimewgll_xx(3,j) + &
+ C3_mxm_m2_m1_5points(i,4)*hprimewgll_xx(4,j) + &
+ C3_mxm_m2_m1_5points(i,5)*hprimewgll_xx(5,j)
+ enddo
+ enddo
+
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ fac1 = wgllwgll_yz(j,k)
+ fac2 = wgllwgll_xz(i,k)
+ fac3 = wgllwgll_xy(i,j)
+
+ ! sum contributions
+ sum_terms(1,i,j,k) = - (fac1*newtempx1(i,j,k) + fac2*newtempx2(i,j,k) + fac3*newtempx3(i,j,k))
+ sum_terms(2,i,j,k) = - (fac1*newtempy1(i,j,k) + fac2*newtempy2(i,j,k) + fac3*newtempy3(i,j,k))
+ sum_terms(3,i,j,k) = - (fac1*newtempz1(i,j,k) + fac2*newtempz2(i,j,k) + fac3*newtempz3(i,j,k))
+
+ if(GRAVITY_VAL) sum_terms(:,i,j,k) = sum_terms(:,i,j,k) + rho_s_H(:,i,j,k)
+
+ enddo ! NGLLX
+ enddo ! NGLLY
+ enddo ! NGLLZ
+
+ ! sum contributions from each element to the global mesh and add gravity terms
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+ accel(:,iglob) = accel(:,iglob) + sum_terms(:,i,j,k)
+ enddo
+ enddo
+ enddo
+
+ ! update memory variables based upon the Runge-Kutta scheme
+ ! convention for attenuation
+ ! term in xx = 1
+ ! term in yy = 2
+ ! term in xy = 3
+ ! term in xz = 4
+ ! term in yz = 5
+ ! term in zz not computed since zero trace
+ ! This is because we only implement Q_\mu attenuation and not Q_\kappa.
+ ! Note that this does *NOT* imply that there is no attenuation for P waves
+ ! because for Q_\kappa = infinity one gets (see for instance Dahlen and Tromp (1998)
+ ! equation (9.59) page 350): Q_\alpha = Q_\mu * 3 * (V_p/V_s)^2 / 4
+ ! therefore Q_\alpha is not zero; for instance for V_p / V_s = sqrt(3)
+ ! we get Q_\alpha = (9 / 4) * Q_\mu = 2.25 * Q_\mu
+
+ if(ATTENUATION_VAL) then
+
+ ! use Runge-Kutta scheme to march in time
+ do i_sls = 1,N_SLS
+ do i_memory = 1,5
+
+ ! get coefficients for that standard linear solid
+ ! IMPROVE we use mu_v here even if there is some anisotropy
+ ! IMPROVE we should probably use an average value instead
+
+ ! reformatted R_memory to handle large factor_common and reduced [alpha,beta,gamma]val
+ if(ATTENUATION_3D_VAL) then
+ factor_common_c44_muv = factor_common(i_sls,:,:,:,ispec)
+ else
+ factor_common_c44_muv(:,:,:) = factor_common(i_sls,1,1,1,iregion_selected)
+ endif
+ if(ANISOTROPIC_3D_MANTLE_VAL) then
+ factor_common_c44_muv = factor_common_c44_muv * c44store(:,:,:,ispec)
+ else
+ factor_common_c44_muv = factor_common_c44_muv * muvstore(:,:,:,ispec)
+ endif
+
+ R_memory(i_memory,i_sls,:,:,:,ispec) = alphaval(i_sls) * &
+ R_memory(i_memory,i_sls,:,:,:,ispec) + &
+ factor_common_c44_muv * &
+ (betaval(i_sls) * epsilondev(i_memory,:,:,:,ispec) + &
+ gammaval(i_sls) * epsilondev_loc(i_memory,:,:,:))
+ enddo
+ enddo
+
+ endif
+
+ ! save deviatoric strain for Runge-Kutta scheme
+ if(COMPUTE_AND_STORE_STRAIN) epsilondev(:,:,:,:,ispec) = epsilondev_loc(:,:,:,:)
+
+ enddo ! spectral element loop NSPEC_CRUST_MANTLE
+
+ end subroutine compute_forces_crust_mantle_Dev
+
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_inner_core.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_inner_core.f90 2010-01-06 23:41:13 UTC (rev 16126)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_inner_core.f90 2010-01-07 14:04:22 UTC (rev 16127)
@@ -600,3 +600,686 @@
end subroutine compute_forces_inner_core
+
+!
+!--------------------------------------------------------------------------------------------------
+!
+
+ subroutine compute_forces_inner_core_Dev(minus_gravity_table,density_table,minus_deriv_gravity_table, &
+ displ,accel,xstore,ystore,zstore, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
+ kappavstore,muvstore,ibool,idoubling, &
+ c11store,c33store,c12store,c13store,c44store,R_memory,epsilondev,epsilon_trace_over_3,&
+ one_minus_sum_beta,alphaval,betaval,gammaval,factor_common, &
+ vx,vy,vz,vnspec,COMPUTE_AND_STORE_STRAIN, AM_V)
+
+ implicit none
+
+ include "constants.h"
+
+! include values created by the mesher
+! done for performance only using static allocation to allow for loop unrolling
+ include "OUTPUT_FILES/values_from_mesher.h"
+
+! attenuation_model_variables
+ type attenuation_model_variables
+ sequence
+ double precision min_period, max_period
+ double precision :: QT_c_source ! Source Frequency
+ double precision, dimension(:), pointer :: Qtau_s ! tau_sigma
+ double precision, dimension(:), pointer :: QrDisc ! Discontinutitues Defined
+ double precision, dimension(:), pointer :: Qr ! Radius
+ integer, dimension(:), pointer :: interval_Q ! Steps
+ double precision, dimension(:), pointer :: Qmu ! Shear Attenuation
+ double precision, dimension(:,:), pointer :: Qtau_e ! tau_epsilon
+ double precision, dimension(:), pointer :: Qomsb, Qomsb2 ! one_minus_sum_beta
+ double precision, dimension(:,:), pointer :: Qfc, Qfc2 ! factor_common
+ double precision, dimension(:), pointer :: Qsf, Qsf2 ! scale_factor
+ integer, dimension(:), pointer :: Qrmin ! Max and Mins of idoubling
+ integer, dimension(:), pointer :: Qrmax ! Max and Mins of idoubling
+ integer :: Qn ! Number of points
+ end type attenuation_model_variables
+
+ type (attenuation_model_variables) AM_V
+! attenuation_model_variables
+
+! for forward or backward simulations
+ logical COMPUTE_AND_STORE_STRAIN
+
+! displacement and acceleration
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_INNER_CORE) :: displ,accel
+
+! for attenuation
+! memory variables R_ij are stored at the local rather than global level
+! to allow for optimization of cache access by compiler
+ integer i_sls,i_memory
+ real(kind=CUSTOM_REAL) R_xx_val,R_yy_val
+
+! variable lengths for factor_common and one_minus_sum_beta
+ integer vx, vy, vz, vnspec
+
+ real(kind=CUSTOM_REAL), dimension(vx, vy, vz, vnspec) :: one_minus_sum_beta
+
+ real(kind=CUSTOM_REAL), dimension(N_SLS, vx, vy, vz, vnspec) :: factor_common
+ real(kind=CUSTOM_REAL), dimension(N_SLS) :: alphaval,betaval,gammaval
+ real(kind=CUSTOM_REAL), dimension(NGLLX, NGLLY, NGLLZ) :: factor_common_use
+
+ real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ATTENUATION) :: R_memory
+ real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: epsilondev
+ real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: epsilon_trace_over_3
+
+! array with the local to global mapping per slice
+ integer, dimension(NSPEC_INNER_CORE) :: idoubling
+
+! arrays with mesh parameters per slice
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: ibool
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: xix,xiy,xiz, &
+ etax,etay,etaz,gammax,gammay,gammaz
+
+! array with derivatives of Lagrange polynomials and precalculated products
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx,hprimewgll_xx
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY) :: wgllwgll_xy
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ) :: wgllwgll_xz
+ real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLZ) :: wgllwgll_yz
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
+ tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: kappavstore,muvstore
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: &
+ c11store,c33store,c12store,c13store,c44store
+
+ integer ispec,iglob,ispec_strain
+ integer i,j,k !,l
+
+ real(kind=CUSTOM_REAL) xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
+ real(kind=CUSTOM_REAL) duxdxl,duxdyl,duxdzl,duydxl,duydyl,duydzl,duzdxl,duzdyl,duzdzl
+
+ real(kind=CUSTOM_REAL) duxdxl_plus_duydyl,duxdxl_plus_duzdzl,duydyl_plus_duzdzl
+ real(kind=CUSTOM_REAL) duxdyl_plus_duydxl,duzdxl_plus_duxdzl,duzdyl_plus_duydzl
+
+ real(kind=CUSTOM_REAL) sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz
+
+ !real(kind=CUSTOM_REAL) hp1,hp2,hp3
+ real(kind=CUSTOM_REAL) fac1,fac2,fac3
+ real(kind=CUSTOM_REAL) lambdal,mul,lambdalplus2mul
+ real(kind=CUSTOM_REAL) kappal
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NGLLZ) :: sum_terms
+
+ real(kind=CUSTOM_REAL) minus_sum_beta
+ real(kind=CUSTOM_REAL) c11l,c33l,c12l,c13l,c44l
+
+! for gravity
+ integer int_radius
+ real(kind=CUSTOM_REAL) sigma_yx,sigma_zx,sigma_zy
+ double precision radius,rho,minus_g,minus_dg
+ double precision minus_g_over_radius,minus_dg_plus_g_over_radius
+ double precision cos_theta,sin_theta,cos_phi,sin_phi
+ double precision cos_theta_sq,sin_theta_sq,cos_phi_sq,sin_phi_sq
+ double precision theta,phi,factor,gxl,gyl,gzl,sx_l,sy_l,sz_l
+ double precision Hxxl,Hyyl,Hzzl,Hxyl,Hxzl,Hyzl
+ double precision, dimension(NRAD_GRAVITY) :: minus_gravity_table,density_table,minus_deriv_gravity_table
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NGLLZ) :: rho_s_H
+ double precision, dimension(NGLLX,NGLLY,NGLLZ) :: wgll_cube
+ real(kind=CUSTOM_REAL), dimension(NGLOB_INNER_CORE) :: xstore,ystore,zstore
+ integer iregion_selected
+
+ real(kind=CUSTOM_REAL) radius_cr
+
+! Deville
+! manually inline the calls to the Deville et al. (2002) routines
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: dummyx_loc,dummyy_loc,dummyz_loc, &
+ newtempx1,newtempx2,newtempx3,newtempy1,newtempy2,newtempy3,newtempz1,newtempz2,newtempz3
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xxT,hprimewgll_xxT
+
+ integer, parameter :: m1 = NGLLX, m2 = NGLLX * NGLLY
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,m2) :: B1_m1_m2_5points,B2_m1_m2_5points,B3_m1_m2_5points
+ real(kind=CUSTOM_REAL), dimension(m1,m2) :: C1_m1_m2_5points,C2_m1_m2_5points,C3_m1_m2_5points
+ real(kind=CUSTOM_REAL), dimension(m1,m2) :: E1_m1_m2_5points,E2_m1_m2_5points,E3_m1_m2_5points
+
+ equivalence(dummyx_loc,B1_m1_m2_5points)
+ equivalence(dummyy_loc,B2_m1_m2_5points)
+ equivalence(dummyz_loc,B3_m1_m2_5points)
+ equivalence(tempx1,C1_m1_m2_5points)
+ equivalence(tempy1,C2_m1_m2_5points)
+ equivalence(tempz1,C3_m1_m2_5points)
+ equivalence(newtempx1,E1_m1_m2_5points)
+ equivalence(newtempy1,E2_m1_m2_5points)
+ equivalence(newtempz1,E3_m1_m2_5points)
+
+ real(kind=CUSTOM_REAL), dimension(m2,NGLLX) :: A1_mxm_m2_m1_5points,A2_mxm_m2_m1_5points,A3_mxm_m2_m1_5points
+ real(kind=CUSTOM_REAL), dimension(m2,m1) :: C1_mxm_m2_m1_5points,C2_mxm_m2_m1_5points,C3_mxm_m2_m1_5points
+ real(kind=CUSTOM_REAL), dimension(m2,m1) :: E1_mxm_m2_m1_5points,E2_mxm_m2_m1_5points,E3_mxm_m2_m1_5points
+
+ equivalence(dummyx_loc,A1_mxm_m2_m1_5points)
+ equivalence(dummyy_loc,A2_mxm_m2_m1_5points)
+ equivalence(dummyz_loc,A3_mxm_m2_m1_5points)
+ equivalence(tempx3,C1_mxm_m2_m1_5points)
+ equivalence(tempy3,C2_mxm_m2_m1_5points)
+ equivalence(tempz3,C3_mxm_m2_m1_5points)
+ equivalence(newtempx3,E1_mxm_m2_m1_5points)
+ equivalence(newtempy3,E2_mxm_m2_m1_5points)
+ equivalence(newtempz3,E3_mxm_m2_m1_5points)
+
+
+
+! ****************************************************
+! big loop over all spectral elements in the solid
+! ****************************************************
+
+
+! set acceleration to zero
+ accel(:,:) = 0._CUSTOM_REAL
+
+ do ispec = 1,NSPEC_INNER_CORE
+
+ ! exclude fictitious elements in central cube
+ if(idoubling(ispec) /= IFLAG_IN_FICTITIOUS_CUBE) then
+
+ ! subroutines adapted from Deville, Fischer and Mund, High-order methods
+ ! for incompressible fluid flow, Cambridge University Press (2002),
+ ! pages 386 and 389 and Figure 8.3.1
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+ dummyx_loc(i,j,k) = displ(1,iglob)
+ dummyy_loc(i,j,k) = displ(2,iglob)
+ dummyz_loc(i,j,k) = displ(3,iglob)
+ enddo
+ enddo
+ enddo
+ do j=1,m2
+ do i=1,m1
+ C1_m1_m2_5points(i,j) = hprime_xx(i,1)*B1_m1_m2_5points(1,j) + &
+ hprime_xx(i,2)*B1_m1_m2_5points(2,j) + &
+ hprime_xx(i,3)*B1_m1_m2_5points(3,j) + &
+ hprime_xx(i,4)*B1_m1_m2_5points(4,j) + &
+ hprime_xx(i,5)*B1_m1_m2_5points(5,j)
+
+ C2_m1_m2_5points(i,j) = hprime_xx(i,1)*B2_m1_m2_5points(1,j) + &
+ hprime_xx(i,2)*B2_m1_m2_5points(2,j) + &
+ hprime_xx(i,3)*B2_m1_m2_5points(3,j) + &
+ hprime_xx(i,4)*B2_m1_m2_5points(4,j) + &
+ hprime_xx(i,5)*B2_m1_m2_5points(5,j)
+
+ C3_m1_m2_5points(i,j) = hprime_xx(i,1)*B3_m1_m2_5points(1,j) + &
+ hprime_xx(i,2)*B3_m1_m2_5points(2,j) + &
+ hprime_xx(i,3)*B3_m1_m2_5points(3,j) + &
+ hprime_xx(i,4)*B3_m1_m2_5points(4,j) + &
+ hprime_xx(i,5)*B3_m1_m2_5points(5,j)
+ enddo
+ enddo
+ do j=1,m1
+ do i=1,m1
+ ! for efficiency it is better to leave this loop on k inside, it leads to slightly faster code
+ do k = 1,NGLLX
+ tempx2(i,j,k) = dummyx_loc(i,1,k)*hprime_xxT(1,j) + &
+ dummyx_loc(i,2,k)*hprime_xxT(2,j) + &
+ dummyx_loc(i,3,k)*hprime_xxT(3,j) + &
+ dummyx_loc(i,4,k)*hprime_xxT(4,j) + &
+ dummyx_loc(i,5,k)*hprime_xxT(5,j)
+
+ tempy2(i,j,k) = dummyy_loc(i,1,k)*hprime_xxT(1,j) + &
+ dummyy_loc(i,2,k)*hprime_xxT(2,j) + &
+ dummyy_loc(i,3,k)*hprime_xxT(3,j) + &
+ dummyy_loc(i,4,k)*hprime_xxT(4,j) + &
+ dummyy_loc(i,5,k)*hprime_xxT(5,j)
+
+ tempz2(i,j,k) = dummyz_loc(i,1,k)*hprime_xxT(1,j) + &
+ dummyz_loc(i,2,k)*hprime_xxT(2,j) + &
+ dummyz_loc(i,3,k)*hprime_xxT(3,j) + &
+ dummyz_loc(i,4,k)*hprime_xxT(4,j) + &
+ dummyz_loc(i,5,k)*hprime_xxT(5,j)
+ enddo
+ enddo
+ enddo
+ do j=1,m1
+ do i=1,m2
+ C1_mxm_m2_m1_5points(i,j) = A1_mxm_m2_m1_5points(i,1)*hprime_xxT(1,j) + &
+ A1_mxm_m2_m1_5points(i,2)*hprime_xxT(2,j) + &
+ A1_mxm_m2_m1_5points(i,3)*hprime_xxT(3,j) + &
+ A1_mxm_m2_m1_5points(i,4)*hprime_xxT(4,j) + &
+ A1_mxm_m2_m1_5points(i,5)*hprime_xxT(5,j)
+
+ C2_mxm_m2_m1_5points(i,j) = A2_mxm_m2_m1_5points(i,1)*hprime_xxT(1,j) + &
+ A2_mxm_m2_m1_5points(i,2)*hprime_xxT(2,j) + &
+ A2_mxm_m2_m1_5points(i,3)*hprime_xxT(3,j) + &
+ A2_mxm_m2_m1_5points(i,4)*hprime_xxT(4,j) + &
+ A2_mxm_m2_m1_5points(i,5)*hprime_xxT(5,j)
+
+ C3_mxm_m2_m1_5points(i,j) = A3_mxm_m2_m1_5points(i,1)*hprime_xxT(1,j) + &
+ A3_mxm_m2_m1_5points(i,2)*hprime_xxT(2,j) + &
+ A3_mxm_m2_m1_5points(i,3)*hprime_xxT(3,j) + &
+ A3_mxm_m2_m1_5points(i,4)*hprime_xxT(4,j) + &
+ A3_mxm_m2_m1_5points(i,5)*hprime_xxT(5,j)
+ enddo
+ enddo
+
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+
+ ! get derivatives of ux, uy and uz with respect to x, y and z
+ xixl = xix(i,j,k,ispec)
+ xiyl = xiy(i,j,k,ispec)
+ xizl = xiz(i,j,k,ispec)
+ etaxl = etax(i,j,k,ispec)
+ etayl = etay(i,j,k,ispec)
+ etazl = etaz(i,j,k,ispec)
+ gammaxl = gammax(i,j,k,ispec)
+ gammayl = gammay(i,j,k,ispec)
+ gammazl = gammaz(i,j,k,ispec)
+
+ ! compute the jacobian
+ jacobianl = 1._CUSTOM_REAL / (xixl*(etayl*gammazl-etazl*gammayl) &
+ - xiyl*(etaxl*gammazl-etazl*gammaxl) &
+ + xizl*(etaxl*gammayl-etayl*gammaxl))
+
+ duxdxl = xixl*tempx1(i,j,k) + etaxl*tempx2(i,j,k) + gammaxl*tempx3(i,j,k)
+ duxdyl = xiyl*tempx1(i,j,k) + etayl*tempx2(i,j,k) + gammayl*tempx3(i,j,k)
+ duxdzl = xizl*tempx1(i,j,k) + etazl*tempx2(i,j,k) + gammazl*tempx3(i,j,k)
+
+ duydxl = xixl*tempy1(i,j,k) + etaxl*tempy2(i,j,k) + gammaxl*tempy3(i,j,k)
+ duydyl = xiyl*tempy1(i,j,k) + etayl*tempy2(i,j,k) + gammayl*tempy3(i,j,k)
+ duydzl = xizl*tempy1(i,j,k) + etazl*tempy2(i,j,k) + gammazl*tempy3(i,j,k)
+
+ duzdxl = xixl*tempz1(i,j,k) + etaxl*tempz2(i,j,k) + gammaxl*tempz3(i,j,k)
+ duzdyl = xiyl*tempz1(i,j,k) + etayl*tempz2(i,j,k) + gammayl*tempz3(i,j,k)
+ duzdzl = xizl*tempz1(i,j,k) + etazl*tempz2(i,j,k) + gammazl*tempz3(i,j,k)
+
+ ! precompute some sums to save CPU time
+ duxdxl_plus_duydyl = duxdxl + duydyl
+ duxdxl_plus_duzdzl = duxdxl + duzdzl
+ duydyl_plus_duzdzl = duydyl + duzdzl
+ duxdyl_plus_duydxl = duxdyl + duydxl
+ duzdxl_plus_duxdzl = duzdxl + duxdzl
+ duzdyl_plus_duydzl = duzdyl + duydzl
+
+ ! compute deviatoric strain
+ if (COMPUTE_AND_STORE_STRAIN) then
+ if(NSPEC_INNER_CORE_STRAIN_ONLY == 1) then
+ ispec_strain = 1
+ else
+ ispec_strain = ispec
+ endif
+ epsilon_trace_over_3(i,j,k,ispec_strain) = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ epsilondev_loc(1,i,j,k) = duxdxl - epsilon_trace_over_3(i,j,k,ispec_strain)
+ epsilondev_loc(2,i,j,k) = duydyl - epsilon_trace_over_3(i,j,k,ispec_strain)
+ epsilondev_loc(3,i,j,k) = 0.5 * duxdyl_plus_duydxl
+ epsilondev_loc(4,i,j,k) = 0.5 * duzdxl_plus_duxdzl
+ epsilondev_loc(5,i,j,k) = 0.5 * duzdyl_plus_duydzl
+ endif
+
+ if(ATTENUATION_VAL) then
+ if(ATTENUATION_3D_VAL) then
+ minus_sum_beta = one_minus_sum_beta(i,j,k,ispec) - 1.0
+ else
+ radius_cr = xstore(ibool(i,j,k,ispec))
+ call get_attenuation_index(idoubling(ispec), dble(radius_cr), iregion_selected, .TRUE., AM_V)
+ minus_sum_beta = one_minus_sum_beta(1,1,1,iregion_selected) - 1.0
+ endif ! ATTENUATION_3D_VAL
+ endif ! ATTENUATION_VAL
+
+ if(ANISOTROPIC_INNER_CORE_VAL) then
+ ! elastic tensor for hexagonal symmetry in reduced notation:
+ !
+ ! c11 c12 c13 0 0 0
+ ! c12 c11 c13 0 0 0
+ ! c13 c13 c33 0 0 0
+ ! 0 0 0 c44 0 0
+ ! 0 0 0 0 c44 0
+ ! 0 0 0 0 0 (c11-c12)/2
+ !
+ ! in terms of the A, C, L, N and F of Love (1927):
+ !
+ ! c11 = A
+ ! c12 = A-2N
+ ! c13 = F
+ ! c33 = C
+ ! c44 = L
+ c11l = c11store(i,j,k,ispec)
+ c12l = c12store(i,j,k,ispec)
+ c13l = c13store(i,j,k,ispec)
+ c33l = c33store(i,j,k,ispec)
+ c44l = c44store(i,j,k,ispec)
+
+ ! use unrelaxed parameters if attenuation
+ if(ATTENUATION_VAL) then
+ mul = muvstore(i,j,k,ispec)
+ c11l = c11l + FOUR_THIRDS * minus_sum_beta * mul
+ c12l = c12l - TWO_THIRDS * minus_sum_beta * mul
+ c13l = c13l - TWO_THIRDS * minus_sum_beta * mul
+ c33l = c33l + FOUR_THIRDS * minus_sum_beta * mul
+ c44l = c44l + minus_sum_beta * mul
+ endif
+
+ sigma_xx = c11l*duxdxl + c12l*duydyl + c13l*duzdzl
+ sigma_yy = c12l*duxdxl + c11l*duydyl + c13l*duzdzl
+ sigma_zz = c13l*duxdxl + c13l*duydyl + c33l*duzdzl
+ sigma_xy = 0.5*(c11l-c12l)*duxdyl_plus_duydxl
+ sigma_xz = c44l*duzdxl_plus_duxdzl
+ sigma_yz = c44l*duzdyl_plus_duydzl
+ else
+
+ ! inner core with no anisotropy, use kappav and muv for instance
+ ! layer with no anisotropy, use kappav and muv for instance
+ kappal = kappavstore(i,j,k,ispec)
+ mul = muvstore(i,j,k,ispec)
+
+ ! use unrelaxed parameters if attenuation
+ if(ATTENUATION_VAL) then
+ if(ATTENUATION_3D_VAL) then
+ mul = mul * one_minus_sum_beta(i,j,k,ispec)
+ else
+ mul = mul * one_minus_sum_beta(1,1,1,iregion_selected)
+ endif
+ endif
+
+ lambdalplus2mul = kappal + FOUR_THIRDS * mul
+ lambdal = lambdalplus2mul - 2.*mul
+
+ ! compute stress sigma
+ sigma_xx = lambdalplus2mul*duxdxl + lambdal*duydyl_plus_duzdzl
+ sigma_yy = lambdalplus2mul*duydyl + lambdal*duxdxl_plus_duzdzl
+ sigma_zz = lambdalplus2mul*duzdzl + lambdal*duxdxl_plus_duydyl
+
+ sigma_xy = mul*duxdyl_plus_duydxl
+ sigma_xz = mul*duzdxl_plus_duxdzl
+ sigma_yz = mul*duzdyl_plus_duydzl
+
+ endif
+
+ ! subtract memory variables if attenuation
+ if(ATTENUATION_VAL) then
+ do i_sls = 1,N_SLS
+ R_xx_val = R_memory(1,i_sls,i,j,k,ispec)
+ R_yy_val = R_memory(2,i_sls,i,j,k,ispec)
+ sigma_xx = sigma_xx - R_xx_val
+ sigma_yy = sigma_yy - R_yy_val
+ sigma_zz = sigma_zz + R_xx_val + R_yy_val
+ sigma_xy = sigma_xy - R_memory(3,i_sls,i,j,k,ispec)
+ sigma_xz = sigma_xz - R_memory(4,i_sls,i,j,k,ispec)
+ sigma_yz = sigma_yz - R_memory(5,i_sls,i,j,k,ispec)
+ enddo
+ endif
+
+ ! define symmetric components of sigma for gravity
+ sigma_yx = sigma_xy
+ sigma_zx = sigma_xz
+ sigma_zy = sigma_yz
+
+ ! compute non-symmetric terms for gravity
+ if(GRAVITY_VAL) then
+
+ ! use mesh coordinates to get theta and phi
+ ! x y and z contain r theta and phi
+ iglob = ibool(i,j,k,ispec)
+ radius = dble(xstore(iglob))
+ theta = dble(ystore(iglob))
+ phi = dble(zstore(iglob))
+
+ ! make sure radius is never zero even for points at center of cube
+ ! because we later divide by radius
+ if(radius < 100.d0 / R_EARTH) radius = 100.d0 / R_EARTH
+
+ cos_theta = dcos(theta)
+ sin_theta = dsin(theta)
+ cos_phi = dcos(phi)
+ sin_phi = dsin(phi)
+
+ ! get g, rho and dg/dr=dg
+ ! spherical components of the gravitational acceleration
+ ! for efficiency replace with lookup table every 100 m in radial direction
+ ! make sure we never use zero for point exactly at the center of the Earth
+ int_radius = max(1,nint(radius * R_EARTH_KM * 10.d0))
+ minus_g = minus_gravity_table(int_radius)
+ minus_dg = minus_deriv_gravity_table(int_radius)
+ rho = density_table(int_radius)
+
+ ! Cartesian components of the gravitational acceleration
+ gxl = minus_g*sin_theta*cos_phi
+ gyl = minus_g*sin_theta*sin_phi
+ gzl = minus_g*cos_theta
+
+ ! Cartesian components of gradient of gravitational acceleration
+ ! obtained from spherical components
+ minus_g_over_radius = minus_g / radius
+ minus_dg_plus_g_over_radius = minus_dg - minus_g_over_radius
+
+ cos_theta_sq = cos_theta**2
+ sin_theta_sq = sin_theta**2
+ cos_phi_sq = cos_phi**2
+ sin_phi_sq = sin_phi**2
+
+ Hxxl = minus_g_over_radius*(cos_phi_sq*cos_theta_sq + sin_phi_sq) + cos_phi_sq*minus_dg*sin_theta_sq
+ Hyyl = minus_g_over_radius*(cos_phi_sq + cos_theta_sq*sin_phi_sq) + minus_dg*sin_phi_sq*sin_theta_sq
+ Hzzl = cos_theta_sq*minus_dg + minus_g_over_radius*sin_theta_sq
+ Hxyl = cos_phi*minus_dg_plus_g_over_radius*sin_phi*sin_theta_sq
+ Hxzl = cos_phi*cos_theta*minus_dg_plus_g_over_radius*sin_theta
+ Hyzl = cos_theta*minus_dg_plus_g_over_radius*sin_phi*sin_theta
+
+ iglob = ibool(i,j,k,ispec)
+
+ ! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ ! get displacement and multiply by density to compute G tensor
+ sx_l = rho * dble(displ(1,iglob))
+ sy_l = rho * dble(displ(2,iglob))
+ sz_l = rho * dble(displ(3,iglob))
+
+ ! compute G tensor from s . g and add to sigma (not symmetric)
+ sigma_xx = sigma_xx + sngl(sy_l*gyl + sz_l*gzl)
+ sigma_yy = sigma_yy + sngl(sx_l*gxl + sz_l*gzl)
+ sigma_zz = sigma_zz + sngl(sx_l*gxl + sy_l*gyl)
+
+ sigma_xy = sigma_xy - sngl(sx_l * gyl)
+ sigma_yx = sigma_yx - sngl(sy_l * gxl)
+
+ sigma_xz = sigma_xz - sngl(sx_l * gzl)
+ sigma_zx = sigma_zx - sngl(sz_l * gxl)
+
+ sigma_yz = sigma_yz - sngl(sy_l * gzl)
+ sigma_zy = sigma_zy - sngl(sz_l * gyl)
+
+ ! precompute vector
+ factor = dble(jacobianl) * wgll_cube(i,j,k)
+ rho_s_H(1,i,j,k) = sngl(factor * (sx_l * Hxxl + sy_l * Hxyl + sz_l * Hxzl))
+ rho_s_H(2,i,j,k) = sngl(factor * (sx_l * Hxyl + sy_l * Hyyl + sz_l * Hyzl))
+ rho_s_H(3,i,j,k) = sngl(factor * (sx_l * Hxzl + sy_l * Hyzl + sz_l * Hzzl))
+
+ else
+
+ ! get displacement and multiply by density to compute G tensor
+ sx_l = rho * displ(1,iglob)
+ sy_l = rho * displ(2,iglob)
+ sz_l = rho * displ(3,iglob)
+
+ ! compute G tensor from s . g and add to sigma (not symmetric)
+ sigma_xx = sigma_xx + sy_l*gyl + sz_l*gzl
+ sigma_yy = sigma_yy + sx_l*gxl + sz_l*gzl
+ sigma_zz = sigma_zz + sx_l*gxl + sy_l*gyl
+
+ sigma_xy = sigma_xy - sx_l * gyl
+ sigma_yx = sigma_yx - sy_l * gxl
+
+ sigma_xz = sigma_xz - sx_l * gzl
+ sigma_zx = sigma_zx - sz_l * gxl
+
+ sigma_yz = sigma_yz - sy_l * gzl
+ sigma_zy = sigma_zy - sz_l * gyl
+
+ ! precompute vector
+ factor = jacobianl * wgll_cube(i,j,k)
+ rho_s_H(1,i,j,k) = factor * (sx_l * Hxxl + sy_l * Hxyl + sz_l * Hxzl)
+ rho_s_H(2,i,j,k) = factor * (sx_l * Hxyl + sy_l * Hyyl + sz_l * Hyzl)
+ rho_s_H(3,i,j,k) = factor * (sx_l * Hxzl + sy_l * Hyzl + sz_l * Hzzl)
+
+ endif
+
+ endif ! end of section with gravity terms
+
+ ! form dot product with test vector, non-symmetric form
+ tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_yx*xiyl + sigma_zx*xizl)
+ tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_zy*xizl)
+ tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl)
+
+ tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_yx*etayl + sigma_zx*etazl)
+ tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_zy*etazl)
+ tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl)
+
+ tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_yx*gammayl + sigma_zx*gammazl)
+ tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_zy*gammazl)
+ tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl)
+
+ enddo
+ enddo
+ enddo
+
+ ! subroutines adapted from Deville, Fischer and Mund, High-order methods
+ ! for incompressible fluid flow, Cambridge University Press (2002),
+ ! pages 386 and 389 and Figure 8.3.1
+ do j=1,m2
+ do i=1,m1
+ E1_m1_m2_5points(i,j) = hprimewgll_xxT(i,1)*C1_m1_m2_5points(1,j) + &
+ hprimewgll_xxT(i,2)*C1_m1_m2_5points(2,j) + &
+ hprimewgll_xxT(i,3)*C1_m1_m2_5points(3,j) + &
+ hprimewgll_xxT(i,4)*C1_m1_m2_5points(4,j) + &
+ hprimewgll_xxT(i,5)*C1_m1_m2_5points(5,j)
+
+ E2_m1_m2_5points(i,j) = hprimewgll_xxT(i,1)*C2_m1_m2_5points(1,j) + &
+ hprimewgll_xxT(i,2)*C2_m1_m2_5points(2,j) + &
+ hprimewgll_xxT(i,3)*C2_m1_m2_5points(3,j) + &
+ hprimewgll_xxT(i,4)*C2_m1_m2_5points(4,j) + &
+ hprimewgll_xxT(i,5)*C2_m1_m2_5points(5,j)
+
+ E3_m1_m2_5points(i,j) = hprimewgll_xxT(i,1)*C3_m1_m2_5points(1,j) + &
+ hprimewgll_xxT(i,2)*C3_m1_m2_5points(2,j) + &
+ hprimewgll_xxT(i,3)*C3_m1_m2_5points(3,j) + &
+ hprimewgll_xxT(i,4)*C3_m1_m2_5points(4,j) + &
+ hprimewgll_xxT(i,5)*C3_m1_m2_5points(5,j)
+ enddo
+ enddo
+ do i=1,m1
+ do j=1,m1
+ ! for efficiency it is better to leave this loop on k inside, it leads to slightly faster code
+ do k = 1,NGLLX
+ newtempx2(i,j,k) = tempx2(i,1,k)*hprimewgll_xx(1,j) + &
+ tempx2(i,2,k)*hprimewgll_xx(2,j) + &
+ tempx2(i,3,k)*hprimewgll_xx(3,j) + &
+ tempx2(i,4,k)*hprimewgll_xx(4,j) + &
+ tempx2(i,5,k)*hprimewgll_xx(5,j)
+
+ newtempy2(i,j,k) = tempy2(i,1,k)*hprimewgll_xx(1,j) + &
+ tempy2(i,2,k)*hprimewgll_xx(2,j) + &
+ tempy2(i,3,k)*hprimewgll_xx(3,j) + &
+ tempy2(i,4,k)*hprimewgll_xx(4,j) + &
+ tempy2(i,5,k)*hprimewgll_xx(5,j)
+
+ newtempz2(i,j,k) = tempz2(i,1,k)*hprimewgll_xx(1,j) + &
+ tempz2(i,2,k)*hprimewgll_xx(2,j) + &
+ tempz2(i,3,k)*hprimewgll_xx(3,j) + &
+ tempz2(i,4,k)*hprimewgll_xx(4,j) + &
+ tempz2(i,5,k)*hprimewgll_xx(5,j)
+ enddo
+ enddo
+ enddo
+ do j=1,m1
+ do i=1,m2
+ E1_mxm_m2_m1_5points(i,j) = C1_mxm_m2_m1_5points(i,1)*hprimewgll_xx(1,j) + &
+ C1_mxm_m2_m1_5points(i,2)*hprimewgll_xx(2,j) + &
+ C1_mxm_m2_m1_5points(i,3)*hprimewgll_xx(3,j) + &
+ C1_mxm_m2_m1_5points(i,4)*hprimewgll_xx(4,j) + &
+ C1_mxm_m2_m1_5points(i,5)*hprimewgll_xx(5,j)
+
+ E2_mxm_m2_m1_5points(i,j) = C2_mxm_m2_m1_5points(i,1)*hprimewgll_xx(1,j) + &
+ C2_mxm_m2_m1_5points(i,2)*hprimewgll_xx(2,j) + &
+ C2_mxm_m2_m1_5points(i,3)*hprimewgll_xx(3,j) + &
+ C2_mxm_m2_m1_5points(i,4)*hprimewgll_xx(4,j) + &
+ C2_mxm_m2_m1_5points(i,5)*hprimewgll_xx(5,j)
+
+ E3_mxm_m2_m1_5points(i,j) = C3_mxm_m2_m1_5points(i,1)*hprimewgll_xx(1,j) + &
+ C3_mxm_m2_m1_5points(i,2)*hprimewgll_xx(2,j) + &
+ C3_mxm_m2_m1_5points(i,3)*hprimewgll_xx(3,j) + &
+ C3_mxm_m2_m1_5points(i,4)*hprimewgll_xx(4,j) + &
+ C3_mxm_m2_m1_5points(i,5)*hprimewgll_xx(5,j)
+ enddo
+ enddo
+
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ fac1 = wgllwgll_yz(j,k)
+ fac2 = wgllwgll_xz(i,k)
+ fac3 = wgllwgll_xy(i,j)
+
+ ! sum contributions
+ sum_terms(1,i,j,k) = - (fac1*newtempx1(i,j,k) + fac2*newtempx2(i,j,k) + fac3*newtempx3(i,j,k))
+ sum_terms(2,i,j,k) = - (fac1*newtempy1(i,j,k) + fac2*newtempy2(i,j,k) + fac3*newtempy3(i,j,k))
+ sum_terms(3,i,j,k) = - (fac1*newtempz1(i,j,k) + fac2*newtempz2(i,j,k) + fac3*newtempz3(i,j,k))
+
+ if(GRAVITY_VAL) sum_terms(:,i,j,k) = sum_terms(:,i,j,k) + rho_s_H(:,i,j,k)
+
+ enddo
+ enddo
+ enddo
+
+ ! sum contributions from each element to the global mesh and add gravity terms
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+ accel(:,iglob) = accel(:,iglob) + sum_terms(:,i,j,k)
+ enddo
+ enddo
+ enddo
+
+ ! use Runge-Kutta scheme to march memory variables in time
+ ! convention for attenuation
+ ! term in xx = 1
+ ! term in yy = 2
+ ! term in xy = 3
+ ! term in xz = 4
+ ! term in yz = 5
+ ! term in zz not computed since zero trace
+ ! This is because we only implement Q_\mu attenuation and not Q_\kappa.
+ ! Note that this does *NOT* imply that there is no attenuation for P waves
+ ! because for Q_\kappa = infinity one gets (see for instance Dahlen and Tromp (1998)
+ ! equation (9.59) page 350): Q_\alpha = Q_\mu * 3 * (V_p/V_s)^2 / 4
+ ! therefore Q_\alpha is not zero; for instance for V_p / V_s = sqrt(3)
+ ! we get Q_\alpha = (9 / 4) * Q_\mu = 2.25 * Q_\mu
+ if(ATTENUATION_VAL) then
+
+ do i_sls = 1,N_SLS
+ if(ATTENUATION_3D_VAL) then
+ factor_common_use = factor_common(i_sls,:,:,:,ispec)
+ else
+ factor_common_use(:,:,:) = factor_common(i_sls,1,1,1,iregion_selected)
+ endif
+ do i_memory = 1,5
+ R_memory(i_memory,i_sls,:,:,:,ispec) = &
+ alphaval(i_sls) * &
+ R_memory(i_memory,i_sls,:,:,:,ispec) + muvstore(:,:,:,ispec) * &
+ factor_common_use * &
+ (betaval(i_sls) * &
+ epsilondev(i_memory,:,:,:,ispec) + gammaval(i_sls) * epsilondev_loc(i_memory,:,:,:))
+ enddo
+ enddo
+
+ endif
+
+ ! save deviatoric strain for Runge-Kutta scheme
+ if(COMPUTE_AND_STORE_STRAIN) epsilondev(:,:,:,:,ispec) = epsilondev_loc(:,:,:,:)
+
+ endif ! end test to exclude fictitious elements in central cube
+
+ enddo ! spectral element loop
+
+ end subroutine compute_forces_inner_core_Dev
+
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90 2010-01-06 23:41:13 UTC (rev 16126)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90 2010-01-07 14:04:22 UTC (rev 16127)
@@ -865,6 +865,10 @@
integer, dimension(NB_SQUARE_CORNERS,NB_CUT_CASE) :: DIFF_NSPEC1D_RADIAL
integer, dimension(NB_SQUARE_EDGES_ONEDIR,NB_CUT_CASE) :: DIFF_NSPEC2D_XI,DIFF_NSPEC2D_ETA
+ ! Deville routine
+ logical :: USE_DEVILLE
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xxT,hprimewgll_xxT
+
! ************** PROGRAM STARTS HERE **************
! initialize the MPI communicator and start the NPROCTOT MPI processes.
@@ -1490,6 +1494,20 @@
hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube)
+ if( NGLLX == 5 .and. NGLLY == 5 .and. NGLLZ == 5 ) then
+ USE_DEVILLE = .true.
+
+ ! define transpose of derivation matrix
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+ hprime_xxT(j,i) = hprime_xx(i,j)
+ hprimewgll_xxT(j,i) = hprimewgll_xx(i,j)
+ enddo
+ enddo
+ else
+ USE_DEVILLE = .false.
+ endif
+
! read topography and bathymetry file
if(myrank == 0 .and. (TOPOGRAPHY .or. OCEANS)) call read_topo_bathy_file(ibathy_topo)
! broadcast the information read on the master to the nodes
@@ -3626,12 +3644,37 @@
! for anisotropy and gravity, x y and z contain r theta and phi
- call compute_forces_crust_mantle(minus_gravity_table,density_table,minus_deriv_gravity_table, &
+ if( USE_DEVILLE ) then
+ call compute_forces_crust_mantle_Dev(minus_gravity_table,density_table,minus_deriv_gravity_table, &
displ_crust_mantle,accel_crust_mantle, &
xstore_crust_mantle,ystore_crust_mantle,zstore_crust_mantle, &
xix_crust_mantle,xiy_crust_mantle,xiz_crust_mantle, &
etax_crust_mantle,etay_crust_mantle,etaz_crust_mantle, &
gammax_crust_mantle,gammay_crust_mantle,gammaz_crust_mantle, &
+ hprime_xx,hprime_xxT, &
+ hprimewgll_xx,hprimewgll_xxT, &
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
+ kappavstore_crust_mantle,kappahstore_crust_mantle,muvstore_crust_mantle, &
+ muhstore_crust_mantle,eta_anisostore_crust_mantle, &
+ c11store_crust_mantle,c12store_crust_mantle,c13store_crust_mantle, &
+ c14store_crust_mantle,c15store_crust_mantle,c16store_crust_mantle, &
+ c22store_crust_mantle,c23store_crust_mantle,c24store_crust_mantle, &
+ c25store_crust_mantle,c26store_crust_mantle,c33store_crust_mantle, &
+ c34store_crust_mantle,c35store_crust_mantle,c36store_crust_mantle, &
+ c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
+ c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
+ ibool_crust_mantle,idoubling_crust_mantle, &
+ R_memory_crust_mantle,epsilondev_crust_mantle,eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
+ alphaval,betaval,gammaval,factor_common_crust_mantle, &
+ size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
+ size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),COMPUTE_AND_STORE_STRAIN,AM_V)
+ else
+ call compute_forces_crust_mantle(minus_gravity_table,density_table,minus_deriv_gravity_table, &
+ displ_crust_mantle,accel_crust_mantle, &
+ xstore_crust_mantle,ystore_crust_mantle,zstore_crust_mantle, &
+ xix_crust_mantle,xiy_crust_mantle,xiz_crust_mantle, &
+ etax_crust_mantle,etay_crust_mantle,etaz_crust_mantle, &
+ gammax_crust_mantle,gammay_crust_mantle,gammaz_crust_mantle, &
hprime_xx,hprime_yy,hprime_zz, &
hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
@@ -3649,15 +3692,41 @@
alphaval,betaval,gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),COMPUTE_AND_STORE_STRAIN,AM_V)
-
+ endif
+
if (SIMULATION_TYPE == 3) then
! for anisotropy and gravity, x y and z contain r theta and phi
- call compute_forces_crust_mantle(minus_gravity_table,density_table,minus_deriv_gravity_table, &
+ if( USE_DEVILLE ) then
+ call compute_forces_crust_mantle_Dev(minus_gravity_table,density_table,minus_deriv_gravity_table, &
b_displ_crust_mantle,b_accel_crust_mantle, &
xstore_crust_mantle,ystore_crust_mantle,zstore_crust_mantle, &
xix_crust_mantle,xiy_crust_mantle,xiz_crust_mantle, &
etax_crust_mantle,etay_crust_mantle,etaz_crust_mantle, &
gammax_crust_mantle,gammay_crust_mantle,gammaz_crust_mantle, &
+ hprime_xx,hprime_xxT, &
+ hprimewgll_xx,hprimewgll_xxT, &
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
+ kappavstore_crust_mantle,kappahstore_crust_mantle,muvstore_crust_mantle, &
+ muhstore_crust_mantle,eta_anisostore_crust_mantle, &
+ c11store_crust_mantle,c12store_crust_mantle,c13store_crust_mantle, &
+ c14store_crust_mantle,c15store_crust_mantle,c16store_crust_mantle, &
+ c22store_crust_mantle,c23store_crust_mantle,c24store_crust_mantle, &
+ c25store_crust_mantle,c26store_crust_mantle,c33store_crust_mantle, &
+ c34store_crust_mantle,c35store_crust_mantle,c36store_crust_mantle, &
+ c44store_crust_mantle,c45store_crust_mantle,c46store_crust_mantle, &
+ c55store_crust_mantle,c56store_crust_mantle,c66store_crust_mantle, &
+ ibool_crust_mantle,idoubling_crust_mantle, &
+ b_R_memory_crust_mantle,b_epsilondev_crust_mantle,b_eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
+ b_alphaval,b_betaval,b_gammaval,factor_common_crust_mantle, &
+ size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
+ size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),COMPUTE_AND_STORE_STRAIN,AM_V)
+ else
+ call compute_forces_crust_mantle(minus_gravity_table,density_table,minus_deriv_gravity_table, &
+ b_displ_crust_mantle,b_accel_crust_mantle, &
+ xstore_crust_mantle,ystore_crust_mantle,zstore_crust_mantle, &
+ xix_crust_mantle,xiy_crust_mantle,xiz_crust_mantle, &
+ etax_crust_mantle,etay_crust_mantle,etaz_crust_mantle, &
+ gammax_crust_mantle,gammay_crust_mantle,gammaz_crust_mantle, &
hprime_xx,hprime_yy,hprime_zz, &
hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
@@ -3675,6 +3744,8 @@
b_alphaval,b_betaval,b_gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),COMPUTE_AND_STORE_STRAIN,AM_V)
+
+ endif
endif
@@ -3908,12 +3979,31 @@
endif ! Stacey conditions
- call compute_forces_inner_core(minus_gravity_table,density_table,minus_deriv_gravity_table, &
+ ! deville routine
+ if( USE_DEVILLE ) then
+ call compute_forces_inner_core_Dev(minus_gravity_table,density_table,minus_deriv_gravity_table, &
displ_inner_core,accel_inner_core, &
xstore_inner_core,ystore_inner_core,zstore_inner_core, &
xix_inner_core,xiy_inner_core,xiz_inner_core, &
etax_inner_core,etay_inner_core,etaz_inner_core, &
gammax_inner_core,gammay_inner_core,gammaz_inner_core, &
+ hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
+ kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
+ c11store_inner_core,c33store_inner_core,c12store_inner_core,c13store_inner_core,c44store_inner_core, &
+ R_memory_inner_core,epsilondev_inner_core, eps_trace_over_3_inner_core,&
+ one_minus_sum_beta_inner_core, &
+ alphaval,betaval,gammaval, &
+ factor_common_inner_core, &
+ size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
+ size(factor_common_inner_core,4), size(factor_common_inner_core,5),COMPUTE_AND_STORE_STRAIN,AM_V)
+ else
+ call compute_forces_inner_core(minus_gravity_table,density_table,minus_deriv_gravity_table, &
+ displ_inner_core,accel_inner_core, &
+ xstore_inner_core,ystore_inner_core,zstore_inner_core, &
+ xix_inner_core,xiy_inner_core,xiz_inner_core, &
+ etax_inner_core,etay_inner_core,etaz_inner_core, &
+ gammax_inner_core,gammay_inner_core,gammaz_inner_core, &
hprime_xx,hprime_yy,hprime_zz,hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
@@ -3923,15 +4013,34 @@
alphaval,betaval,gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),COMPUTE_AND_STORE_STRAIN,AM_V)
+ size(factor_common_inner_core,4), size(factor_common_inner_core,5),COMPUTE_AND_STORE_STRAIN,AM_V)
+ endif
if (SIMULATION_TYPE == 3) then
- call compute_forces_inner_core(minus_gravity_table,density_table,minus_deriv_gravity_table, &
+ if( USE_DEVILLE ) then
+ call compute_forces_inner_core_Dev(minus_gravity_table,density_table,minus_deriv_gravity_table, &
b_displ_inner_core,b_accel_inner_core, &
xstore_inner_core,ystore_inner_core,zstore_inner_core, &
xix_inner_core,xiy_inner_core,xiz_inner_core, &
etax_inner_core,etay_inner_core,etaz_inner_core, &
gammax_inner_core,gammay_inner_core,gammaz_inner_core, &
+ hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
+ kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
+ c11store_inner_core,c33store_inner_core,c12store_inner_core,c13store_inner_core,c44store_inner_core, &
+ b_R_memory_inner_core,b_epsilondev_inner_core, b_eps_trace_over_3_inner_core,&
+ one_minus_sum_beta_inner_core, &
+ b_alphaval,b_betaval,b_gammaval, &
+ factor_common_inner_core, &
+ size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
+ size(factor_common_inner_core,4), size(factor_common_inner_core,5),COMPUTE_AND_STORE_STRAIN,AM_V)
+ else
+ call compute_forces_inner_core(minus_gravity_table,density_table,minus_deriv_gravity_table, &
+ b_displ_inner_core,b_accel_inner_core, &
+ xstore_inner_core,ystore_inner_core,zstore_inner_core, &
+ xix_inner_core,xiy_inner_core,xiz_inner_core, &
+ etax_inner_core,etay_inner_core,etaz_inner_core, &
+ gammax_inner_core,gammay_inner_core,gammaz_inner_core, &
hprime_xx,hprime_yy,hprime_zz,hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
@@ -3941,7 +4050,8 @@
b_alphaval,b_betaval,b_gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),COMPUTE_AND_STORE_STRAIN,AM_V)
+ size(factor_common_inner_core,4), size(factor_common_inner_core,5),COMPUTE_AND_STORE_STRAIN,AM_V)
+ endif
endif
! add the sources
More information about the CIG-COMMITS
mailing list