[cig-commits] r16310 - seismo/3D/SPECFEM3D_GLOBE/trunk
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Mon Feb 22 15:17:53 PST 2010
Author: danielpeter
Date: 2010-02-22 15:17:53 -0800 (Mon, 22 Feb 2010)
New Revision: 16310
Added:
seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_outer_core_Dev.f90
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/Makefile.in
seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_outer_core.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/constants.h.in
seismo/3D/SPECFEM3D_GLOBE/trunk/model_topo_bathy.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90
Log:
added Deville et al. (2002) routines for fluid outer core
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/Makefile.in 2010-02-22 21:58:51 UTC (rev 16309)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/Makefile.in 2010-02-22 23:17:53 UTC (rev 16310)
@@ -210,6 +210,7 @@
$O/compute_forces_inner_core.o \
$O/compute_forces_inner_core_Dev.o \
$O/compute_forces_outer_core.o \
+ $O/compute_forces_outer_core_Dev.o \
$O/compute_kernels.o \
$O/compute_seismograms.o \
$O/compute_stacey_crust_mantle.o \
@@ -412,6 +413,9 @@
$O/compute_forces_outer_core.o: constants.h OUTPUT_FILES/values_from_mesher.h $S/compute_forces_outer_core.f90
${FCCOMPILE_NO_CHECK} -c -o $O/compute_forces_outer_core.o ${FCFLAGS_f90} $S/compute_forces_outer_core.f90
+$O/compute_forces_outer_core_Dev.o: constants.h OUTPUT_FILES/values_from_mesher.h $S/compute_forces_outer_core_Dev.f90
+ ${FCCOMPILE_NO_CHECK} -c -o $O/compute_forces_outer_core_Dev.o ${FCFLAGS_f90} $S/compute_forces_outer_core_Dev.f90
+
$O/compute_forces_inner_core.o: constants.h OUTPUT_FILES/values_from_mesher.h $S/compute_forces_inner_core.f90
${FCCOMPILE_NO_CHECK} -c -o $O/compute_forces_inner_core.o ${FCFLAGS_f90} $S/compute_forces_inner_core.f90
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_outer_core.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_outer_core.f90 2010-02-22 21:58:51 UTC (rev 16309)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_outer_core.f90 2010-02-22 23:17:53 UTC (rev 16310)
@@ -106,7 +106,6 @@
do j=1,NGLLY
do i=1,NGLLX
-
tempx1l = 0._CUSTOM_REAL
tempx2l = 0._CUSTOM_REAL
tempx3l = 0._CUSTOM_REAL
@@ -170,39 +169,38 @@
! add (chi/rho)grad(rho) term in no gravity case
if(.not. GRAVITY_VAL) then
+ ! With regards to the non-gravitating case: we cannot set N^2 = 0 *and* let g = 0.
+ ! We can *either* assume N^2 = 0 but keep gravity g, *or* we can assume that gravity
+ ! is negligible to begin with, as in our GJI 2002a, in which case N does not arise.
+ ! We get:
+ !
+ ! \ddot\chi = \rho^{-1}\kappa\bdel\cdot(\bdel\chi+\chi\bdel\ln\rho)
+ !
+ ! Then the displacement is
+ !
+ ! \bu = \bdel\chi+\chi\bdel\ln\rho = \rho^{-1}\bdel(\rho\chi)
+ !
+ ! and the pressure is
+ !
+ ! p = -\rho\ddot{\chi}
+ !
+ ! Thus in our 2002b GJI paper eqn (21) is wrong, and equation (41)
+ ! in our AGU monograph is incorrect; these equations should be replaced by
+ !
+ ! \ddot\chi = \rho^{-1}\kappa\bdel\cdot(\bdel\chi+\chi\bdel\ln\rho)
+ !
+ ! Note that the fluid potential we use in GJI 2002a differs from the one used here:
+ !
+ ! \chi_GJI2002a = \rho\partial\t\chi
+ !
+ ! such that
+ !
+ ! \bv = \partial_t\bu=\rho^{-1}\bdel\chi_GJI2002a (GJI 2002a eqn 20)
+ !
+ ! p = - \partial_t\chi_GJI2002a (GJI 2002a eqn 19)
-! With regards to the non-gravitating case: we cannot set N^2 = 0 *and* let g = 0.
-! We can *either* assume N^2 = 0 but keep gravity g, *or* we can assume that gravity
-! is negligible to begin with, as in our GJI 2002a, in which case N does not arise.
-! We get:
-!
-! \ddot\chi = \rho^{-1}\kappa\bdel\cdot(\bdel\chi+\chi\bdel\ln\rho)
-!
-! Then the displacement is
-!
-! \bu = \bdel\chi+\chi\bdel\ln\rho = \rho^{-1}\bdel(\rho\chi)
-!
-! and the pressure is
-!
-! p = -\rho\ddot{\chi}
-!
-! Thus in our 2002b GJI paper eqn (21) is wrong, and equation (41)
-! in our AGU monograph is incorrect; these equations should be replaced by
-!
-! \ddot\chi = \rho^{-1}\kappa\bdel\cdot(\bdel\chi+\chi\bdel\ln\rho)
-!
-! Note that the fluid potential we use in GJI 2002a differs from the one used here:
-!
-! \chi_GJI2002a = \rho\partial\t\chi
-!
-! such that
-!
-! \bv = \partial_t\bu=\rho^{-1}\bdel\chi_GJI2002a (GJI 2002a eqn 20)
-!
-! p = - \partial_t\chi_GJI2002a (GJI 2002a eqn 19)
-
-! use mesh coordinates to get theta and phi
-! x y z contain r theta phi
+ ! use mesh coordinates to get theta and phi
+ ! x y z contain r theta phi
iglob = ibool(i,j,k,ispec)
radius = dble(xstore(iglob))
@@ -222,8 +220,8 @@
grad_z_ln_rho = cos_theta * d_ln_density_dr_table(int_radius)
! adding (chi/rho)grad(rho)
- dpotentialdx_with_rot = dpotentialdxl + displfluid(iglob) * grad_x_ln_rho
- dpotentialdy_with_rot = dpotentialdyl + displfluid(iglob) * grad_y_ln_rho
+ dpotentialdx_with_rot = dpotentialdx_with_rot + displfluid(iglob) * grad_x_ln_rho
+ dpotentialdy_with_rot = dpotentialdx_with_rot + displfluid(iglob) * grad_y_ln_rho
dpotentialdzl = dpotentialdzl + displfluid(iglob) * grad_z_ln_rho
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_outer_core_Dev.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_outer_core_Dev.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_outer_core_Dev.f90 2010-02-22 23:17:53 UTC (rev 16310)
@@ -0,0 +1,409 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 4 . 0
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory, California Institute of Technology, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) California Institute of Technology and University of Pau / CNRS / INRIA
+! February 2008
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+ subroutine compute_forces_outer_core_Dev(time,deltat,two_omega_earth, &
+ A_array_rotation,B_array_rotation, &
+ d_ln_density_dr_table, &
+ minus_rho_g_over_kappa_fluid,displfluid,accelfluid, &
+ div_displfluid, &
+ 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, &
+ ibool)
+
+ 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"
+
+! displacement and acceleration
+ real(kind=CUSTOM_REAL), dimension(nglob_outer_core) :: displfluid,accelfluid
+
+! divergence of displacement
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec_outer_core) :: div_displfluid
+
+! arrays with mesh parameters per slice
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nspec_outer_core) :: ibool
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec_outer_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,NGLLX) :: hprime_xxT,hprimewgll_xxT
+ 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
+ double precision, dimension(NGLLX,NGLLY,NGLLZ) :: wgll_cube
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: tempx1,tempx2,tempx3
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: newtempx1,newtempx2,newtempx3
+
+! for gravity
+ integer int_radius
+ double precision radius,theta,phi,gxl,gyl,gzl
+ double precision cos_theta,sin_theta,cos_phi,sin_phi
+ double precision, dimension(NRAD_GRAVITY) :: minus_rho_g_over_kappa_fluid
+ double precision, dimension(NRAD_GRAVITY) :: d_ln_density_dr_table
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: gravity_term
+ real(kind=CUSTOM_REAL), dimension(nglob_outer_core) :: xstore,ystore,zstore
+
+! for the Euler scheme for rotation
+ real(kind=CUSTOM_REAL) time,deltat,two_omega_earth
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_OUTER_CORE_ROTATION) :: &
+ A_array_rotation,B_array_rotation
+
+ real(kind=CUSTOM_REAL) two_omega_deltat,cos_two_omega_t,sin_two_omega_t,A_rotation,B_rotation, &
+ ux_rotation,uy_rotation,dpotentialdx_with_rot,dpotentialdy_with_rot
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: source_euler_A,source_euler_B
+
+ integer ispec,iglob
+ integer i,j,k,l
+
+ real(kind=CUSTOM_REAL) xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
+ real(kind=CUSTOM_REAL) dpotentialdxl,dpotentialdyl,dpotentialdzl
+ real(kind=CUSTOM_REAL) tempx1l,tempx2l,tempx3l,sum_terms
+
+ double precision grad_x_ln_rho,grad_y_ln_rho,grad_z_ln_rho
+
+ ! manually inline the calls to the Deville et al. (2002) routines
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: dummyx_loc
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,m2) :: B1_m1_m2_5points
+ real(kind=CUSTOM_REAL), dimension(m1,m2) :: C1_m1_m2_5points
+ real(kind=CUSTOM_REAL), dimension(m1,m2) :: E1_m1_m2_5points
+
+ equivalence(dummyx_loc,B1_m1_m2_5points)
+ equivalence(tempx1,C1_m1_m2_5points)
+ equivalence(newtempx1,E1_m1_m2_5points)
+
+ real(kind=CUSTOM_REAL), dimension(m2,NGLLX) :: A1_mxm_m2_m1_5points
+ real(kind=CUSTOM_REAL), dimension(m2,m1) :: C1_mxm_m2_m1_5points
+ real(kind=CUSTOM_REAL), dimension(m2,m1) :: E1_mxm_m2_m1_5points
+
+ equivalence(dummyx_loc,A1_mxm_m2_m1_5points)
+ equivalence(tempx3,C1_mxm_m2_m1_5points)
+ equivalence(newtempx3,E1_mxm_m2_m1_5points)
+
+ double precision, dimension(NGLLX,NGLLY,NGLLZ) :: temp_gxl,temp_gyl,temp_gzl
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
+ displ_times_grad_x_ln_rho,displ_times_grad_y_ln_rho,displ_times_grad_z_ln_rho
+
+! ****************************************************
+! big loop over all spectral elements in the fluid
+! ****************************************************
+
+ if (NSPEC_OUTER_CORE_ADJOINT /= 1) div_displfluid(:,:,:,:) = 0._CUSTOM_REAL
+
+ do ispec = 1,NSPEC_OUTER_CORE
+
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+
+ ! stores "displacement"
+ dummyx_loc(i,j,k) = displfluid(iglob)
+
+ ! pre-computes factors
+ ! use mesh coordinates to get theta and phi
+ ! x y z contain r theta phi
+ radius = dble(xstore(iglob))
+ theta = dble(ystore(iglob))
+ phi = dble(zstore(iglob))
+
+ cos_theta = dcos(theta)
+ sin_theta = dsin(theta)
+ cos_phi = dcos(phi)
+ sin_phi = dsin(phi)
+
+ int_radius = nint(radius * R_EARTH_KM * 10.d0)
+
+ if( .not. GRAVITY_VAL ) then
+ ! grad(rho)/rho in Cartesian components
+ displ_times_grad_x_ln_rho(i,j,k) = dummyx_loc(i,j,k) &
+ * sngl(sin_theta * cos_phi * d_ln_density_dr_table(int_radius))
+ displ_times_grad_y_ln_rho(i,j,k) = dummyx_loc(i,j,k) &
+ * sngl(sin_theta * sin_phi * d_ln_density_dr_table(int_radius))
+ displ_times_grad_z_ln_rho(i,j,k) = dummyx_loc(i,j,k) &
+ * sngl(cos_theta * d_ln_density_dr_table(int_radius))
+ else
+ ! Cartesian components of the gravitational acceleration
+ ! integrate and multiply by rho / Kappa
+ temp_gxl(i,j,k) = sin_theta*cos_phi
+ temp_gyl(i,j,k) = sin_theta*sin_phi
+ temp_gzl(i,j,k) = cos_theta
+ endif
+
+ 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
+ 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)
+ enddo
+ enddo
+ do k = 1,NGLLX
+ do j=1,m1
+ do i=1,m1
+ 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)
+ 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)
+ enddo
+ enddo
+
+
+
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+
+ ! get derivatives of velocity potential 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))
+
+ dpotentialdxl = xixl*tempx1(i,j,k) + etaxl*tempx2(i,j,k) + gammaxl*tempx3(i,j,k)
+ dpotentialdyl = xiyl*tempx1(i,j,k) + etayl*tempx2(i,j,k) + gammayl*tempx3(i,j,k)
+ dpotentialdzl = xizl*tempx1(i,j,k) + etazl*tempx2(i,j,k) + gammazl*tempx3(i,j,k)
+
+ ! compute contribution of rotation and add to gradient of potential
+ ! this term has no Z component
+ if(ROTATION_VAL) then
+
+ ! store the source for the Euler scheme for A_rotation and B_rotation
+ two_omega_deltat = deltat * two_omega_earth
+
+ cos_two_omega_t = cos(two_omega_earth*time)
+ sin_two_omega_t = sin(two_omega_earth*time)
+
+ ! time step deltat of Euler scheme is included in the source
+ source_euler_A(i,j,k) = two_omega_deltat &
+ * (cos_two_omega_t * dpotentialdyl + sin_two_omega_t * dpotentialdxl)
+ source_euler_B(i,j,k) = two_omega_deltat &
+ * (sin_two_omega_t * dpotentialdyl - cos_two_omega_t * dpotentialdxl)
+
+ A_rotation = A_array_rotation(i,j,k,ispec)
+ B_rotation = B_array_rotation(i,j,k,ispec)
+
+ ux_rotation = A_rotation*cos_two_omega_t + B_rotation*sin_two_omega_t
+ uy_rotation = - A_rotation*sin_two_omega_t + B_rotation*cos_two_omega_t
+
+ dpotentialdx_with_rot = dpotentialdxl + ux_rotation
+ dpotentialdy_with_rot = dpotentialdyl + uy_rotation
+
+ else
+
+ dpotentialdx_with_rot = dpotentialdxl
+ dpotentialdy_with_rot = dpotentialdyl
+
+ endif ! end of section with rotation
+
+ ! add (chi/rho)grad(rho) term in no gravity case
+ if(.not. GRAVITY_VAL) then
+
+ ! With regards to the non-gravitating case: we cannot set N^2 = 0 *and* let g = 0.
+ ! We can *either* assume N^2 = 0 but keep gravity g, *or* we can assume that gravity
+ ! is negligible to begin with, as in our GJI 2002a, in which case N does not arise.
+ ! We get:
+ !
+ ! \ddot\chi = \rho^{-1}\kappa\bdel\cdot(\bdel\chi+\chi\bdel\ln\rho)
+ !
+ ! Then the displacement is
+ !
+ ! \bu = \bdel\chi+\chi\bdel\ln\rho = \rho^{-1}\bdel(\rho\chi)
+ !
+ ! and the pressure is
+ !
+ ! p = -\rho\ddot{\chi}
+ !
+ ! Thus in our 2002b GJI paper eqn (21) is wrong, and equation (41)
+ ! in our AGU monograph is incorrect; these equations should be replaced by
+ !
+ ! \ddot\chi = \rho^{-1}\kappa\bdel\cdot(\bdel\chi+\chi\bdel\ln\rho)
+ !
+ ! Note that the fluid potential we use in GJI 2002a differs from the one used here:
+ !
+ ! \chi_GJI2002a = \rho\partial\t\chi
+ !
+ ! such that
+ !
+ ! \bv = \partial_t\bu=\rho^{-1}\bdel\chi_GJI2002a (GJI 2002a eqn 20)
+ !
+ ! p = - \partial_t\chi_GJI2002a (GJI 2002a eqn 19)
+
+ ! use mesh coordinates to get theta and phi
+ ! x y z contain r theta phi
+ dpotentialdx_with_rot = dpotentialdx_with_rot + displ_times_grad_x_ln_rho(i,j,k)
+ dpotentialdy_with_rot = dpotentialdx_with_rot + displ_times_grad_y_ln_rho(i,j,k)
+ dpotentialdzl = dpotentialdzl + displ_times_grad_z_ln_rho(i,j,k)
+
+ else ! if gravity is turned on
+
+ ! compute divergence of displacment
+ gxl = temp_gxl(i,j,k)
+ gyl = temp_gyl(i,j,k)
+ gzl = temp_gzl(i,j,k)
+
+ ! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ gravity_term(i,j,k) = &
+ sngl( minus_rho_g_over_kappa_fluid(int_radius) &
+ * dble(jacobianl) * wgll_cube(i,j,k) &
+ * (dble(dpotentialdx_with_rot) * gxl &
+ + dble(dpotentialdy_with_rot) * gyl &
+ + dble(dpotentialdzl) * gzl) )
+ else
+ gravity_term(i,j,k) = minus_rho_g_over_kappa_fluid(int_radius) * &
+ jacobianl * wgll_cube(i,j,k) &
+ * (dpotentialdx_with_rot * gxl &
+ + dpotentialdy_with_rot * gyl &
+ + dpotentialdzl * gzl)
+ endif
+
+ ! divergence of displacement field with gravity on
+ ! note: these calculations are only considered for SIMULATION_TYPE == 1 .and. SAVE_FORWARD
+ ! and one has set MOVIE_VOLUME_TYPE == 4 when MOVIE_VOLUME is .true.;
+ ! in case of SIMULATION_TYPE == 3, it gets overwritten by compute_kernels_outer_core()
+ if (NSPEC_OUTER_CORE_ADJOINT /= 1) then
+ div_displfluid(i,j,k,ispec) = &
+ minus_rho_g_over_kappa_fluid(int_radius) &
+ * (dpotentialdx_with_rot * gxl &
+ + dpotentialdy_with_rot * gyl &
+ + dpotentialdzl * gzl)
+ endif
+
+ endif
+
+ tempx1(i,j,k) = jacobianl*(xixl*dpotentialdx_with_rot &
+ + xiyl*dpotentialdy_with_rot + xizl*dpotentialdzl)
+ tempx2(i,j,k) = jacobianl*(etaxl*dpotentialdx_with_rot &
+ + etayl*dpotentialdy_with_rot + etazl*dpotentialdzl)
+ tempx3(i,j,k) = jacobianl*(gammaxl*dpotentialdx_with_rot &
+ + gammayl*dpotentialdy_with_rot + gammazl*dpotentialdzl)
+
+ 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)
+ enddo
+ enddo
+ do k = 1,NGLLX
+ do j=1,m1
+ do i=1,m1
+ 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)
+ 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)
+ enddo
+ enddo
+
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+
+ ! sum contributions from each element to the global mesh and add gravity term
+ sum_terms = - (wgllwgll_yz(j,k)*newtempx1(i,j,k) &
+ + wgllwgll_xz(i,k)*newtempx2(i,j,k) &
+ + wgllwgll_xy(i,j)*newtempx3(i,j,k))
+
+ if(GRAVITY_VAL) sum_terms = sum_terms + gravity_term(i,j,k)
+
+ iglob = ibool(i,j,k,ispec)
+ accelfluid(iglob) = accelfluid(iglob) + sum_terms
+
+ enddo
+ enddo
+ enddo
+
+ ! update rotation term with Euler scheme
+ if(ROTATION_VAL) then
+ ! use the source saved above
+ A_array_rotation(:,:,:,ispec) = A_array_rotation(:,:,:,ispec) + source_euler_A(:,:,:)
+ B_array_rotation(:,:,:,ispec) = B_array_rotation(:,:,:,ispec) + source_euler_B(:,:,:)
+ endif
+
+ enddo ! spectral element loop
+
+ end subroutine compute_forces_outer_core_Dev
+
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/constants.h.in 2010-02-22 21:58:51 UTC (rev 16309)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/constants.h.in 2010-02-22 23:17:53 UTC (rev 16310)
@@ -109,8 +109,6 @@
! integer, parameter :: RESOLUTION_TOPO_FILE = 1
!! pathname of the topography file (un-smoothed)
! character (len=*), parameter :: PATHNAME_TOPO_FILE = 'DATA/topo_bathy/ETOPO1.xyz'
-!! pathname of the topography file (smoothed)
-! character (len=*), parameter :: PATHNAME_TOPO_FILE = 'DATA/topo_bathy/topo_bathy_etopo1_smoothed_window_7.dat'
! Use GLL points to capture TOPOGRAPHY and ELLIPTICITY
logical,parameter :: USE_GLL = .false.
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/model_topo_bathy.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/model_topo_bathy.f90 2010-02-22 21:58:51 UTC (rev 16309)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/model_topo_bathy.f90 2010-02-22 23:17:53 UTC (rev 16310)
@@ -74,12 +74,17 @@
! use integer array to store values
integer, dimension(NX_BATHY,NY_BATHY) :: ibathy_topo
- integer itopo_x,itopo_y
+ integer itopo_x,itopo_y,ier
call get_value_string(topo_bathy_file, 'model.topoBathy.PATHNAME_TOPO_FILE', PATHNAME_TOPO_FILE)
! reads in topography values from file
- open(unit=13,file=topo_bathy_file,status='old',action='read')
+ open(unit=13,file=trim(topo_bathy_file),status='old',action='read',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening:',trim(topo_bathy_file)
+ call exit_mpi(0,'error opening topography data file')
+ endif
+ ! reads in topography array
do itopo_y=1,NY_BATHY
do itopo_x=1,NX_BATHY
read(13,*) ibathy_topo(itopo_x,itopo_y)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90 2010-02-22 21:58:51 UTC (rev 16309)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90 2010-02-22 23:17:53 UTC (rev 16310)
@@ -1898,56 +1898,57 @@
time = (dble(it-1)*DT-t0)*scale_t_inv
endif
- ! div_displ_outer_core is initialized to zero in the following subroutine.
- call compute_forces_outer_core(time,deltat,two_omega_earth, &
+ if( USE_DEVILLE_VAL ) then
+ ! uses deville et al. (2002) routine
+ call compute_forces_outer_core_Dev(time,deltat,two_omega_earth, &
A_array_rotation,B_array_rotation,d_ln_density_dr_table, &
minus_rho_g_over_kappa_fluid,displ_outer_core,accel_outer_core,div_displ_outer_core, &
xstore_outer_core,ystore_outer_core,zstore_outer_core, &
xix_outer_core,xiy_outer_core,xiz_outer_core, &
etax_outer_core,etay_outer_core,etaz_outer_core, &
gammax_outer_core,gammay_outer_core,gammaz_outer_core, &
+ hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
+ ibool_outer_core)
+ else
+ ! div_displ_outer_core is initialized to zero in the following subroutine.
+ call compute_forces_outer_core(time,deltat,two_omega_earth, &
+ A_array_rotation,B_array_rotation,d_ln_density_dr_table, &
+ minus_rho_g_over_kappa_fluid,displ_outer_core,accel_outer_core,div_displ_outer_core, &
+ xstore_outer_core,ystore_outer_core,zstore_outer_core, &
+ xix_outer_core,xiy_outer_core,xiz_outer_core, &
+ etax_outer_core,etay_outer_core,etaz_outer_core, &
+ gammax_outer_core,gammay_outer_core,gammaz_outer_core, &
hprime_xx,hprime_yy,hprime_zz,hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
ibool_outer_core)
+ endif
-! ! uses deville et al. (2002) routine
-! call compute_forces_outer_core_Dev(time,deltat,two_omega_earth, &
-! A_array_rotation,B_array_rotation,d_ln_density_dr_table, &
-! minus_rho_g_over_kappa_fluid,displ_outer_core,accel_outer_core,div_displ_outer_core, &
-! xstore_outer_core,ystore_outer_core,zstore_outer_core, &
-! xix_outer_core,xiy_outer_core,xiz_outer_core, &
-! etax_outer_core,etay_outer_core,etaz_outer_core, &
-! gammax_outer_core,gammay_outer_core,gammaz_outer_core, &
-! hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
-! wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
-! ibool_outer_core)
-
-
if (SIMULATION_TYPE == 3) then
- call compute_forces_outer_core(time,b_deltat,b_two_omega_earth, &
+ if( USE_DEVILLE_VAL ) then
+ ! uses deville et al. (2002) routine
+ call compute_forces_outer_core_Dev(time,b_deltat,b_two_omega_earth, &
b_A_array_rotation,b_B_array_rotation,d_ln_density_dr_table, &
minus_rho_g_over_kappa_fluid,b_displ_outer_core,b_accel_outer_core,b_div_displ_outer_core, &
xstore_outer_core,ystore_outer_core,zstore_outer_core, &
xix_outer_core,xiy_outer_core,xiz_outer_core, &
etax_outer_core,etay_outer_core,etaz_outer_core, &
gammax_outer_core,gammay_outer_core,gammaz_outer_core, &
+ hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
+ ibool_outer_core)
+ else
+ call compute_forces_outer_core(time,b_deltat,b_two_omega_earth, &
+ b_A_array_rotation,b_B_array_rotation,d_ln_density_dr_table, &
+ minus_rho_g_over_kappa_fluid,b_displ_outer_core,b_accel_outer_core,b_div_displ_outer_core, &
+ xstore_outer_core,ystore_outer_core,zstore_outer_core, &
+ xix_outer_core,xiy_outer_core,xiz_outer_core, &
+ etax_outer_core,etay_outer_core,etaz_outer_core, &
+ gammax_outer_core,gammay_outer_core,gammaz_outer_core, &
hprime_xx,hprime_yy,hprime_zz,hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
ibool_outer_core)
-
-! ! uses deville et al. (2002) routine
-! call compute_forces_outer_core_Dev(time,b_deltat,b_two_omega_earth, &
-! b_A_array_rotation,b_B_array_rotation,d_ln_density_dr_table, &
-! minus_rho_g_over_kappa_fluid,b_displ_outer_core,b_accel_outer_core,b_div_displ_outer_core, &
-! xstore_outer_core,ystore_outer_core,zstore_outer_core, &
-! xix_outer_core,xiy_outer_core,xiz_outer_core, &
-! etax_outer_core,etay_outer_core,etaz_outer_core, &
-! gammax_outer_core,gammay_outer_core,gammaz_outer_core, &
-! hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
-! wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
-! ibool_outer_core)
-
-
+ endif
endif
! Stacey absorbing boundaries
More information about the CIG-COMMITS
mailing list