[cig-commits] r12867 - in seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta: . src
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Thu Sep 11 17:57:15 PDT 2008
Author: dkomati1
Date: 2008-09-11 17:57:15 -0700 (Thu, 11 Sep 2008)
New Revision: 12867
Added:
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/assemble_MPI_central_cube_block.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/assemble_MPI_scalar_block.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_forces_CM_IC.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_forces_OC.f90
Removed:
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_forces_crust_mantle.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_forces_inner_core.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_forces_outer_core.f90
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/Makefile
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/attenuation_model.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_element_properties.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/create_regions_mesh.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/get_model.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/main_program.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/meshfem3D.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/specfem3D.F90
Log:
first step to switch to non blocking MPI: split compute_forces in two (inner and outer elements)
and merged crust_mantle and inner_core into a single routine called CM_IC
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/Makefile
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/Makefile 2008-09-11 21:43:27 UTC (rev 12866)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/Makefile 2008-09-12 00:57:15 UTC (rev 12867)
@@ -37,7 +37,7 @@
FC = ifort
MPIFC = mpif90
MPIFLAGS = -DUSE_MPI # -lmpi
-#FLAGS_NO_CHECK = -O1 -vec-report0 -no-heap-arrays -e03 -std03 -implicitnone -warn truncated_source -warn argument_checking -warn unused -warn declarations -warn alignments -warn ignore_loc -warn usage -check nobounds -align sequence -assume byterecl -fpe0 -ftz -traceback -ftrapuv # -mcmodel=medium -shared-intel
+#FLAGS_NO_CHECK = -O1 -vec-report0 -no-heap-arrays -e03 -std03 -implicitnone -warn truncated_source -warn argument_checking -warn unused -warn declarations -warn alignments -warn ignore_loc -warn usage -check all -align sequence -assume byterecl -fpe0 -ftz -traceback -ftrapuv # -mcmodel=medium -shared-intel
FLAGS_NO_CHECK = -O3 -xP -vec-report0 -no-heap-arrays -e03 -std03 -implicitnone -warn truncated_source -warn argument_checking -warn unused -warn declarations -warn alignments -warn ignore_loc -warn usage -check nobounds -align sequence -assume byterecl -fpe3 -ftz # -mcmodel=medium -shared-intel
# we need the -no-heap-arrays flag to force the compiler to allocate memory on the stack
# instead of on the heap to minimize memory fragmentation
@@ -113,6 +113,7 @@
$O/anisotropic_inner_core_model.o \
$O/anisotropic_mantle_model.o \
$O/assemble_MPI_scalar.o \
+ $O/assemble_MPI_scalar_block.o \
$O/assemble_MPI_vector.o \
$O/attenuation_model.o \
$O/calc_jacobian.o \
@@ -191,9 +192,9 @@
# values_from_mesher.h
SOLVER_ARRAY_OBJECTS = \
$O/assemble_MPI_central_cube.o \
- $O/compute_forces_crust_mantle.o \
- $O/compute_forces_inner_core.o \
- $O/compute_forces_outer_core.o \
+ $O/assemble_MPI_central_cube_block.o \
+ $O/compute_forces_CM_IC.o \
+ $O/compute_forces_OC.o \
$O/specfem3D.o \
$(EMPTY_MACRO)
@@ -268,15 +269,12 @@
$O/specfem3D.o: $(SPECINC)/constants.h $(OUTPUT_FILES_INC)/values_from_mesher.h $S/specfem3D.F90
${MPIFCCOMPILE_NO_CHECK} -c -o $O/specfem3D.o ${FCFLAGS_f90} $S/specfem3D.F90
-$O/compute_forces_crust_mantle.o: $(SPECINC)/constants.h $(OUTPUT_FILES_INC)/values_from_mesher.h $S/compute_forces_crust_mantle.f90
- ${FCCOMPILE_NO_CHECK} -c -o $O/compute_forces_crust_mantle.o ${FCFLAGS_f90} $S/compute_forces_crust_mantle.f90
+$O/compute_forces_CM_IC.o: $(SPECINC)/constants.h $(OUTPUT_FILES_INC)/values_from_mesher.h $S/compute_forces_CM_IC.f90
+ ${FCCOMPILE_NO_CHECK} -c -o $O/compute_forces_CM_IC.o ${FCFLAGS_f90} $S/compute_forces_CM_IC.f90
-$O/compute_forces_outer_core.o: $(SPECINC)/constants.h $(OUTPUT_FILES_INC)/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_OC.o: $(SPECINC)/constants.h $(OUTPUT_FILES_INC)/values_from_mesher.h $S/compute_forces_OC.f90
+ ${FCCOMPILE_NO_CHECK} -c -o $O/compute_forces_OC.o ${FCFLAGS_f90} $S/compute_forces_OC.f90
-$O/compute_forces_inner_core.o: $(SPECINC)/constants.h $(OUTPUT_FILES_INC)/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
-
### use MPI here
$O/assemble_MPI_vector.o: $(SPECINC)/constants.h $S/assemble_MPI_vector.F90
${MPIFCCOMPILE_NO_CHECK} -c -o $O/assemble_MPI_vector.o ${FCFLAGS_f90} $S/assemble_MPI_vector.F90
@@ -285,9 +283,15 @@
$O/assemble_MPI_scalar.o: $(SPECINC)/constants.h $S/assemble_MPI_scalar.F90
${MPIFCCOMPILE_NO_CHECK} -c -o $O/assemble_MPI_scalar.o ${FCFLAGS_f90} $S/assemble_MPI_scalar.F90
+$O/assemble_MPI_scalar_block.o: $(SPECINC)/constants.h $S/assemble_MPI_scalar_block.F90
+ ${MPIFCCOMPILE_NO_CHECK} -c -o $O/assemble_MPI_scalar_block.o ${FCFLAGS_f90} $S/assemble_MPI_scalar_block.F90
+
$O/assemble_MPI_central_cube.o: $(SPECINC)/constants.h $(OUTPUT_FILES_INC)/values_from_mesher.h $S/assemble_MPI_central_cube.F90
${MPIFCCOMPILE_NO_CHECK} -c -o $O/assemble_MPI_central_cube.o ${FCFLAGS_f90} $S/assemble_MPI_central_cube.F90
+$O/assemble_MPI_central_cube_block.o: $(SPECINC)/constants.h $(OUTPUT_FILES_INC)/values_from_mesher.h $S/assemble_MPI_central_cube_block.F90
+ ${MPIFCCOMPILE_NO_CHECK} -c -o $O/assemble_MPI_central_cube_block.o ${FCFLAGS_f90} $S/assemble_MPI_central_cube_block.F90
+
###
### regular compilation options here
###
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/assemble_MPI_central_cube_block.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/assemble_MPI_central_cube_block.F90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/assemble_MPI_central_cube_block.F90 2008-09-12 00:57:15 UTC (rev 12867)
@@ -0,0 +1,276 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 4 . 1
+! --------------------------------------------------
+!
+! 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
+! August 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 assemble_MPI_central_cube_block(ichunk,nb_msgs_theor_in_cube, sender_from_slices_to_cube, &
+ npoin2D_cube_from_slices, buffer_all_cube_from_slices, buffer_slices, buffer_slices2, ibool_central_cube, &
+ receiver_cube_from_slices, ibool_inner_core, idoubling_inner_core, NSPEC_INNER_CORE, &
+ ibelm_bottom_inner_core, NSPEC2D_BOTTOM_INNER_CORE,NGLOB_INNER_CORE,vector_assemble,ndim_assemble)
+
+ implicit none
+
+! standard include of the MPI library
+#ifdef USE_MPI
+ include 'mpif.h'
+#endif
+ include 'constants.h'
+
+! for matching with central cube in inner core
+ integer ichunk, nb_msgs_theor_in_cube, npoin2D_cube_from_slices
+ integer, dimension(nb_msgs_theor_in_cube) :: sender_from_slices_to_cube
+ double precision, dimension(npoin2D_cube_from_slices,NDIM) :: buffer_slices,buffer_slices2
+ double precision, dimension(nb_msgs_theor_in_cube,npoin2D_cube_from_slices,NDIM) :: buffer_all_cube_from_slices
+ integer, dimension(nb_msgs_theor_in_cube,npoin2D_cube_from_slices):: ibool_central_cube
+ integer receiver_cube_from_slices
+
+! local to global mapping
+ integer NSPEC_INNER_CORE,NSPEC2D_BOTTOM_INNER_CORE, NGLOB_INNER_CORE
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: ibool_inner_core
+ integer, dimension(NSPEC_INNER_CORE) :: idoubling_inner_core
+ integer, dimension(NSPEC2D_BOTTOM_INNER_CORE) :: ibelm_bottom_inner_core
+
+! vector
+ integer ndim_assemble
+ real(kind=CUSTOM_REAL), dimension(ndim_assemble,NGLOB_INNER_CORE) :: vector_assemble
+
+ integer ipoin,idimension, ispec2D, ispec
+ integer i,j,k
+ integer sender,receiver,imsg
+
+ real(kind=CUSTOM_REAL), dimension(NGLOB_INNER_CORE) :: array_central_cube
+
+! MPI status of messages to be received
+#ifdef USE_MPI
+ integer, dimension(MPI_STATUS_SIZE) :: msg_status
+ integer :: ier
+#endif
+
+! mask
+ logical, dimension(NGLOB_INNER_CORE) :: mask
+
+!---
+!--- now use buffers to assemble mass matrix with central cube once and for all
+!---
+
+! on chunks AB and AB_ANTIPODE, receive all the messages from slices
+ if(ichunk == CHUNK_AB .or. ichunk == CHUNK_AB_ANTIPODE) then
+
+ do imsg = 1,nb_msgs_theor_in_cube-1
+
+! receive buffers from slices
+ sender = sender_from_slices_to_cube(imsg)
+#ifdef USE_MPI
+ call MPI_RECV(buffer_slices, &
+ ndim_assemble*npoin2D_cube_from_slices,MPI_DOUBLE_PRECISION,sender, &
+ itag,MPI_COMM_WORLD,msg_status,ier)
+#endif
+
+! copy buffer in 2D array for each slice
+ buffer_all_cube_from_slices(imsg,:,1:ndim_assemble) = buffer_slices(:,1:ndim_assemble)
+
+ enddo
+ endif
+
+! send info to central cube from all the slices except those in CHUNK_AB & CHUNK_AB_ANTIPODE
+ if(ichunk /= CHUNK_AB .and. ichunk /= CHUNK_AB_ANTIPODE) then
+
+! for bottom elements in contact with central cube from the slices side
+ ipoin = 0
+ do ispec2D = 1,NSPEC2D_BOTTOM_INNER_CORE
+
+ ispec = ibelm_bottom_inner_core(ispec2D)
+
+! only for DOFs exactly on surface of central cube (bottom of these elements)
+ k = 1
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+ ipoin = ipoin + 1
+ buffer_slices(ipoin,1:ndim_assemble) = dble(vector_assemble(1:ndim_assemble,ibool_inner_core(i,j,k,ispec)))
+ enddo
+ enddo
+ enddo
+
+! send buffer to central cube
+ receiver = receiver_cube_from_slices
+#ifdef USE_MPI
+ call MPI_SEND(buffer_slices,ndim_assemble*npoin2D_cube_from_slices, &
+ MPI_DOUBLE_PRECISION,receiver,itag,MPI_COMM_WORLD,ier)
+#endif
+
+ endif ! end sending info to central cube
+
+
+! exchange of their bottom faces between chunks AB and AB_ANTIPODE
+ if(ichunk == CHUNK_AB .or. ichunk == CHUNK_AB_ANTIPODE) then
+
+ ipoin = 0
+ do ispec = NSPEC_INNER_CORE, 1, -1
+ if (idoubling_inner_core(ispec) == IFLAG_BOTTOM_CENTRAL_CUBE) then
+ k = 1
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+ ipoin = ipoin + 1
+ buffer_slices(ipoin,1:ndim_assemble) = dble(vector_assemble(1:ndim_assemble,ibool_inner_core(i,j,k,ispec)))
+ enddo
+ enddo
+ endif
+ enddo
+
+ sender = sender_from_slices_to_cube(nb_msgs_theor_in_cube)
+
+#ifdef USE_MPI
+ call MPI_SENDRECV(buffer_slices,ndim_assemble*npoin2D_cube_from_slices,MPI_DOUBLE_PRECISION,receiver_cube_from_slices, &
+ itag,buffer_slices2,ndim_assemble*npoin2D_cube_from_slices,&
+ MPI_DOUBLE_PRECISION,sender,itag,MPI_COMM_WORLD,msg_status,ier)
+#endif
+
+ buffer_all_cube_from_slices(nb_msgs_theor_in_cube,:,1:ndim_assemble) = buffer_slices2(:,1:ndim_assemble)
+
+ endif
+
+!--- now we need to assemble the contributions
+
+ if(ichunk == CHUNK_AB .or. ichunk == CHUNK_AB_ANTIPODE) then
+
+ do idimension = 1,ndim_assemble
+! erase contributions to central cube array
+ array_central_cube(:) = 0._CUSTOM_REAL
+
+! use indirect addressing to store contributions only once
+! distinguish between single and double precision for reals
+ do imsg = 1,nb_msgs_theor_in_cube-1
+ do ipoin = 1,npoin2D_cube_from_slices
+ if(CUSTOM_REAL == SIZE_REAL) then
+ array_central_cube(ibool_central_cube(imsg,ipoin)) = sngl(buffer_all_cube_from_slices(imsg,ipoin,idimension))
+ else
+ array_central_cube(ibool_central_cube(imsg,ipoin)) = buffer_all_cube_from_slices(imsg,ipoin,idimension)
+ endif
+ enddo
+ enddo
+! add the constribution of AB or AB_ANTIPODE to sum with the external slices on the edges
+! use a mask to avoid taking the same point into account several times.
+ mask(:) = .false.
+ do ipoin = 1,npoin2D_cube_from_slices
+ if (.not. mask(ibool_central_cube(nb_msgs_theor_in_cube,ipoin))) then
+ if(CUSTOM_REAL == SIZE_REAL) then
+ array_central_cube(ibool_central_cube(nb_msgs_theor_in_cube,ipoin)) = &
+ array_central_cube(ibool_central_cube(nb_msgs_theor_in_cube,ipoin)) + &
+ sngl(buffer_all_cube_from_slices(nb_msgs_theor_in_cube,ipoin,idimension))
+ else
+ array_central_cube(ibool_central_cube(nb_msgs_theor_in_cube,ipoin)) = &
+ array_central_cube(ibool_central_cube(nb_msgs_theor_in_cube,ipoin)) + &
+ buffer_all_cube_from_slices(nb_msgs_theor_in_cube,ipoin,idimension)
+ endif
+ mask(ibool_central_cube(nb_msgs_theor_in_cube,ipoin)) = .true.
+ endif
+ enddo
+
+! suppress degrees of freedom already assembled at top of cube on edges
+ do ispec = 1,NSPEC_INNER_CORE
+ if(idoubling_inner_core(ispec) == IFLAG_TOP_CENTRAL_CUBE) then
+ k = NGLLZ
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+ array_central_cube(ibool_inner_core(i,j,k,ispec)) = 0._CUSTOM_REAL
+ enddo
+ enddo
+ endif
+ enddo
+
+! assemble contributions
+ vector_assemble(idimension,:) = vector_assemble(idimension,:) + array_central_cube(:)
+
+! copy sum back
+ do imsg = 1,nb_msgs_theor_in_cube-1
+ do ipoin = 1,npoin2D_cube_from_slices
+ buffer_all_cube_from_slices(imsg,ipoin,idimension) = vector_assemble(idimension,ibool_central_cube(imsg,ipoin))
+ enddo
+ enddo
+
+ enddo
+
+ endif
+
+!----------
+
+! receive info from central cube on all the slices except those in CHUNK_AB & CHUNK_AB_ANTIPODE
+ if(ichunk /= CHUNK_AB .and. ichunk /= CHUNK_AB_ANTIPODE) then
+
+! receive buffers from slices
+ sender = receiver_cube_from_slices
+#ifdef USE_MPI
+ call MPI_RECV(buffer_slices, &
+ ndim_assemble*npoin2D_cube_from_slices,MPI_DOUBLE_PRECISION,sender, &
+ itag,MPI_COMM_WORLD,msg_status,ier)
+#endif
+
+! for bottom elements in contact with central cube from the slices side
+ ipoin = 0
+ do ispec2D = 1,NSPEC2D_BOTTOM_INNER_CORE
+
+ ispec = ibelm_bottom_inner_core(ispec2D)
+
+! only for DOFs exactly on surface of central cube (bottom of these elements)
+ k = 1
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+ ipoin = ipoin + 1
+
+! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ vector_assemble(1:ndim_assemble,ibool_inner_core(i,j,k,ispec)) = sngl(buffer_slices(ipoin,1:ndim_assemble))
+ else
+ vector_assemble(1:ndim_assemble,ibool_inner_core(i,j,k,ispec)) = buffer_slices(ipoin,1:ndim_assemble)
+ endif
+
+ enddo
+ enddo
+ enddo
+
+ endif ! end receiving info from central cube
+
+!------- send info back from central cube to slices
+
+! on chunk AB & CHUNK_AB_ANTIPODE, send all the messages to slices
+ if(ichunk == CHUNK_AB .or. ichunk == CHUNK_AB_ANTIPODE) then
+
+ do imsg = 1,nb_msgs_theor_in_cube-1
+
+! copy buffer in 2D array for each slice
+ buffer_slices(:,1:ndim_assemble) = buffer_all_cube_from_slices(imsg,:,1:ndim_assemble)
+
+! send buffers to slices
+ receiver = sender_from_slices_to_cube(imsg)
+#ifdef USE_MPI
+ call MPI_SEND(buffer_slices,ndim_assemble*npoin2D_cube_from_slices, &
+ MPI_DOUBLE_PRECISION,receiver,itag,MPI_COMM_WORLD,ier)
+#endif
+
+ enddo
+ endif
+
+end subroutine assemble_MPI_central_cube_block
+
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/assemble_MPI_scalar_block.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/assemble_MPI_scalar_block.F90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/assemble_MPI_scalar_block.F90 2008-09-12 00:57:15 UTC (rev 12867)
@@ -0,0 +1,491 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 4 . 1
+! --------------------------------------------------
+!
+! 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
+! August 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.
+!
+!=====================================================================
+
+!----
+!---- assemble the contributions between slices and chunks using MPI
+!----
+
+ subroutine assemble_MPI_scalar_block(myrank,array_val,nglob, &
+ iproc_xi,iproc_eta,ichunk,addressing, &
+ iboolleft_xi,iboolright_xi,iboolleft_eta,iboolright_eta, &
+ npoin2D_faces,npoin2D_xi,npoin2D_eta, &
+ iboolfaces,iboolcorner, &
+ iprocfrom_faces,iprocto_faces,imsg_type, &
+ iproc_master_corners,iproc_worker1_corners,iproc_worker2_corners, &
+ buffer_send_faces_scalar,buffer_received_faces_scalar,npoin2D_max_all, &
+ buffer_send_chunkcorners_scalar,buffer_recv_chunkcorners_scalar, &
+ NUMMSGS_FACES,NUM_MSG_TYPES,NCORNERSCHUNKS, &
+ NPROC_XI,NPROC_ETA,NGLOB1D_RADIAL, &
+ NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX,NGLOB2DMAX_XY,NCHUNKS)
+
+ implicit none
+
+! standard include of the MPI library
+#ifdef USE_MPI
+ include 'mpif.h'
+#endif
+
+ include "constants.h"
+#ifdef USE_MPI
+ include "precision.h"
+#endif
+
+ integer myrank,nglob,NCHUNKS
+
+! array to assemble
+ real(kind=CUSTOM_REAL), dimension(nglob) :: array_val
+
+ integer iproc_xi,iproc_eta,ichunk
+ integer, dimension(NB_SQUARE_EDGES_ONEDIR) :: npoin2D_xi,npoin2D_eta
+ integer npoin2D_faces(NUMFACES_SHARED)
+
+ integer NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX,NGLOB2DMAX_XY
+ integer NPROC_XI,NPROC_ETA,NGLOB1D_RADIAL
+ integer NUMMSGS_FACES,NUM_MSG_TYPES,NCORNERSCHUNKS
+
+! for addressing of the slices
+ integer, dimension(NCHUNKS,0:NPROC_XI-1,0:NPROC_ETA-1) :: addressing
+
+! 2-D addressing and buffers for summation between slices
+ integer, dimension(NGLOB2DMAX_XMIN_XMAX) :: iboolleft_xi,iboolright_xi
+ integer, dimension(NGLOB2DMAX_YMIN_YMAX) :: iboolleft_eta,iboolright_eta
+
+! indirect addressing for each corner of the chunks
+ integer, dimension(NGLOB1D_RADIAL,NUMCORNERS_SHARED) :: iboolcorner
+ integer icount_corners
+
+ integer :: npoin2D_max_all
+ integer, dimension(NGLOB2DMAX_XY,NUMFACES_SHARED) :: iboolfaces
+ real(kind=CUSTOM_REAL), dimension(npoin2D_max_all) :: buffer_send_faces_scalar,buffer_received_faces_scalar
+
+! buffers for send and receive between corners of the chunks
+ real(kind=CUSTOM_REAL), dimension(NGLOB1D_RADIAL) :: buffer_send_chunkcorners_scalar,buffer_recv_chunkcorners_scalar
+
+! ---- arrays to assemble between chunks
+
+! communication pattern for faces between chunks
+ integer, dimension(NUMMSGS_FACES) :: iprocfrom_faces,iprocto_faces,imsg_type
+
+! communication pattern for corners between chunks
+ integer, dimension(NCORNERSCHUNKS) :: iproc_master_corners,iproc_worker1_corners,iproc_worker2_corners
+
+! MPI status of messages to be received
+#ifdef USE_MPI
+ integer, dimension(MPI_STATUS_SIZE) :: msg_status
+#endif
+
+ integer :: ipoin,ipoin2D,ipoin1D
+ integer :: sender,receiver
+ integer :: imsg,imsg_loop
+ integer :: icount_faces,npoin2D_chunks
+
+#ifdef USE_MPI
+ integer :: ier
+#endif
+
+! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+
+! check flag to see if we need to assemble (might be turned off when debugging)
+ if (.not. ACTUALLY_ASSEMBLE_MPI_SLICES) return
+
+! here we have to assemble all the contributions between slices using MPI
+
+!----
+!---- assemble the contributions between slices using MPI
+!----
+
+!----
+!---- first assemble along xi using the 2-D topology
+!----
+
+! assemble along xi only if more than one slice
+ if(NPROC_XI > 1) then
+
+! slices copy the right face into the buffer
+ do ipoin=1,npoin2D_xi(2)
+ buffer_send_faces_scalar(ipoin) = array_val(iboolright_xi(ipoin))
+ enddo
+
+! send messages forward along each row
+ if(iproc_xi == 0) then
+#ifdef USE_MPI
+ sender = MPI_PROC_NULL
+#endif
+ else
+ sender = addressing(ichunk,iproc_xi - 1,iproc_eta)
+ endif
+ if(iproc_xi == NPROC_XI-1) then
+#ifdef USE_MPI
+ receiver = MPI_PROC_NULL
+#endif
+ else
+ receiver = addressing(ichunk,iproc_xi + 1,iproc_eta)
+ endif
+#ifdef USE_MPI
+ call MPI_SENDRECV(buffer_send_faces_scalar,npoin2D_xi(2),CUSTOM_MPI_TYPE,receiver, &
+ itag2,buffer_received_faces_scalar,npoin2D_xi(1),CUSTOM_MPI_TYPE,sender, &
+ itag,MPI_COMM_WORLD,msg_status,ier)
+#endif
+
+! all slices add the buffer received to the contributions on the left face
+ if(iproc_xi > 0) then
+ do ipoin=1,npoin2D_xi(1)
+ array_val(iboolleft_xi(ipoin)) = array_val(iboolleft_xi(ipoin)) + &
+ buffer_received_faces_scalar(ipoin)
+ enddo
+ endif
+
+! the contributions are correctly assembled on the left side of each slice
+! now we have to send the result back to the sender
+! all slices copy the left face into the buffer
+ do ipoin=1,npoin2D_xi(1)
+ buffer_send_faces_scalar(ipoin) = array_val(iboolleft_xi(ipoin))
+ enddo
+
+! send messages backward along each row
+ if(iproc_xi == NPROC_XI-1) then
+#ifdef USE_MPI
+ sender = MPI_PROC_NULL
+#endif
+ else
+ sender = addressing(ichunk,iproc_xi + 1,iproc_eta)
+ endif
+ if(iproc_xi == 0) then
+#ifdef USE_MPI
+ receiver = MPI_PROC_NULL
+#endif
+ else
+ receiver = addressing(ichunk,iproc_xi - 1,iproc_eta)
+ endif
+#ifdef USE_MPI
+ call MPI_SENDRECV(buffer_send_faces_scalar,npoin2D_xi(1),CUSTOM_MPI_TYPE,receiver, &
+ itag2,buffer_received_faces_scalar,npoin2D_xi(2),CUSTOM_MPI_TYPE,sender, &
+ itag,MPI_COMM_WORLD,msg_status,ier)
+#endif
+
+! all slices copy the buffer received to the contributions on the right face
+ if(iproc_xi < NPROC_XI-1) then
+ do ipoin=1,npoin2D_xi(2)
+ array_val(iboolright_xi(ipoin)) = buffer_received_faces_scalar(ipoin)
+ enddo
+ endif
+
+ endif
+
+!----
+!---- then assemble along eta using the 2-D topology
+!----
+
+! assemble along eta only if more than one slice
+ if(NPROC_ETA > 1) then
+
+! slices copy the right face into the buffer
+ do ipoin=1,npoin2D_eta(2)
+ buffer_send_faces_scalar(ipoin) = array_val(iboolright_eta(ipoin))
+ enddo
+
+! send messages forward along each row
+ if(iproc_eta == 0) then
+#ifdef USE_MPI
+ sender = MPI_PROC_NULL
+#endif
+ else
+ sender = addressing(ichunk,iproc_xi,iproc_eta - 1)
+ endif
+ if(iproc_eta == NPROC_ETA-1) then
+#ifdef USE_MPI
+ receiver = MPI_PROC_NULL
+#endif
+ else
+ receiver = addressing(ichunk,iproc_xi,iproc_eta + 1)
+ endif
+#ifdef USE_MPI
+ call MPI_SENDRECV(buffer_send_faces_scalar,npoin2D_eta(2),CUSTOM_MPI_TYPE,receiver, &
+ itag2,buffer_received_faces_scalar,npoin2D_eta(1),CUSTOM_MPI_TYPE,sender, &
+ itag,MPI_COMM_WORLD,msg_status,ier)
+#endif
+
+! all slices add the buffer received to the contributions on the left face
+ if(iproc_eta > 0) then
+ do ipoin=1,npoin2D_eta(1)
+ array_val(iboolleft_eta(ipoin)) = array_val(iboolleft_eta(ipoin)) + &
+ buffer_received_faces_scalar(ipoin)
+ enddo
+ endif
+
+! the contributions are correctly assembled on the left side of each slice
+! now we have to send the result back to the sender
+! all slices copy the left face into the buffer
+ do ipoin=1,npoin2D_eta(1)
+ buffer_send_faces_scalar(ipoin) = array_val(iboolleft_eta(ipoin))
+ enddo
+
+! send messages backward along each row
+ if(iproc_eta == NPROC_ETA-1) then
+#ifdef USE_MPI
+ sender = MPI_PROC_NULL
+#endif
+ else
+ sender = addressing(ichunk,iproc_xi,iproc_eta + 1)
+ endif
+ if(iproc_eta == 0) then
+#ifdef USE_MPI
+ receiver = MPI_PROC_NULL
+#endif
+ else
+ receiver = addressing(ichunk,iproc_xi,iproc_eta - 1)
+ endif
+#ifdef USE_MPI
+ call MPI_SENDRECV(buffer_send_faces_scalar,npoin2D_eta(1),CUSTOM_MPI_TYPE,receiver, &
+ itag2,buffer_received_faces_scalar,npoin2D_eta(2),CUSTOM_MPI_TYPE,sender, &
+ itag,MPI_COMM_WORLD,msg_status,ier)
+#endif
+
+! all slices copy the buffer received to the contributions on the right face
+ if(iproc_eta < NPROC_ETA-1) then
+ do ipoin=1,npoin2D_eta(2)
+ array_val(iboolright_eta(ipoin)) = buffer_received_faces_scalar(ipoin)
+ enddo
+ endif
+
+ endif
+
+!----
+!---- start MPI assembling phase between chunks
+!----
+
+! check flag to see if we need to assemble (might be turned off when debugging)
+! and do not assemble if only one chunk
+ if (.not. ACTUALLY_ASSEMBLE_MPI_CHUNKS .or. NCHUNKS == 1) return
+
+! ***************************************************************
+! transmit messages in forward direction (iprocfrom -> iprocto)
+! ***************************************************************
+
+!---- put slices in receive mode
+!---- a given slice can belong to at most two faces
+
+! use three step scheme that can never deadlock
+! scheme for faces cannot deadlock even if NPROC_XI = NPROC_ETA = 1
+ do imsg_loop = 1,NUM_MSG_TYPES
+
+ icount_faces = 0
+ do imsg = 1,NUMMSGS_FACES
+ if(myrank==iprocfrom_faces(imsg) .or. &
+ myrank==iprocto_faces(imsg)) icount_faces = icount_faces + 1
+ if(myrank==iprocto_faces(imsg) .and. imsg_type(imsg) == imsg_loop) then
+ sender = iprocfrom_faces(imsg)
+ npoin2D_chunks = npoin2D_faces(icount_faces)
+#ifdef USE_MPI
+ call MPI_RECV(buffer_received_faces_scalar, &
+ npoin2D_chunks,CUSTOM_MPI_TYPE,sender, &
+ itag,MPI_COMM_WORLD,msg_status,ier)
+#endif
+ do ipoin2D=1,npoin2D_chunks
+ array_val(iboolfaces(ipoin2D,icount_faces)) = &
+ array_val(iboolfaces(ipoin2D,icount_faces)) + buffer_received_faces_scalar(ipoin2D)
+ enddo
+ endif
+ enddo
+
+!---- put slices in send mode
+!---- a given slice can belong to at most two faces
+ icount_faces = 0
+ do imsg = 1,NUMMSGS_FACES
+ if(myrank==iprocfrom_faces(imsg) .or. &
+ myrank==iprocto_faces(imsg)) icount_faces = icount_faces + 1
+ if(myrank==iprocfrom_faces(imsg) .and. imsg_type(imsg) == imsg_loop) then
+ receiver = iprocto_faces(imsg)
+ npoin2D_chunks = npoin2D_faces(icount_faces)
+ do ipoin2D=1,npoin2D_chunks
+ buffer_send_faces_scalar(ipoin2D) = array_val(iboolfaces(ipoin2D,icount_faces))
+ enddo
+#ifdef USE_MPI
+ call MPI_SEND(buffer_send_faces_scalar,npoin2D_chunks, &
+ CUSTOM_MPI_TYPE,receiver,itag,MPI_COMM_WORLD,ier)
+#endif
+ endif
+ enddo
+
+! *********************************************************************
+! transmit messages back in opposite direction (iprocto -> iprocfrom)
+! *********************************************************************
+
+!---- put slices in receive mode
+!---- a given slice can belong to at most two faces
+
+ icount_faces = 0
+ do imsg = 1,NUMMSGS_FACES
+ if(myrank==iprocfrom_faces(imsg) .or. &
+ myrank==iprocto_faces(imsg)) icount_faces = icount_faces + 1
+ if(myrank==iprocfrom_faces(imsg) .and. imsg_type(imsg) == imsg_loop) then
+ sender = iprocto_faces(imsg)
+ npoin2D_chunks = npoin2D_faces(icount_faces)
+#ifdef USE_MPI
+ call MPI_RECV(buffer_received_faces_scalar, &
+ npoin2D_chunks,CUSTOM_MPI_TYPE,sender, &
+ itag,MPI_COMM_WORLD,msg_status,ier)
+#endif
+ do ipoin2D=1,npoin2D_chunks
+ array_val(iboolfaces(ipoin2D,icount_faces)) = buffer_received_faces_scalar(ipoin2D)
+ enddo
+ endif
+ enddo
+
+!---- put slices in send mode
+!---- a given slice can belong to at most two faces
+ icount_faces = 0
+ do imsg = 1,NUMMSGS_FACES
+ if(myrank==iprocfrom_faces(imsg) .or. &
+ myrank==iprocto_faces(imsg)) icount_faces = icount_faces + 1
+ if(myrank==iprocto_faces(imsg) .and. imsg_type(imsg) == imsg_loop) then
+ receiver = iprocfrom_faces(imsg)
+ npoin2D_chunks = npoin2D_faces(icount_faces)
+ do ipoin2D=1,npoin2D_chunks
+ buffer_send_faces_scalar(ipoin2D) = array_val(iboolfaces(ipoin2D,icount_faces))
+ enddo
+#ifdef USE_MPI
+ call MPI_SEND(buffer_send_faces_scalar,npoin2D_chunks, &
+ CUSTOM_MPI_TYPE,receiver,itag,MPI_COMM_WORLD,ier)
+#endif
+ endif
+ enddo
+
+! end of anti-deadlocking loop
+ enddo
+
+!----
+!---- start MPI assembling corners
+!----
+
+! scheme for corners cannot deadlock even if NPROC_XI = NPROC_ETA = 1
+
+! ***************************************************************
+! transmit messages in forward direction (two workers -> master)
+! ***************************************************************
+
+ icount_corners = 0
+
+ do imsg = 1,NCORNERSCHUNKS
+
+ if(myrank == iproc_master_corners(imsg) .or. &
+ myrank == iproc_worker1_corners(imsg) .or. &
+ (NCHUNKS /= 2 .and. myrank == iproc_worker2_corners(imsg))) icount_corners = icount_corners + 1
+
+!---- receive messages from the two workers on the master
+ if(myrank==iproc_master_corners(imsg)) then
+
+! receive from worker #1 and add to local array
+ sender = iproc_worker1_corners(imsg)
+#ifdef USE_MPI
+ call MPI_RECV(buffer_recv_chunkcorners_scalar,NGLOB1D_RADIAL, &
+ CUSTOM_MPI_TYPE,sender,itag,MPI_COMM_WORLD,msg_status,ier)
+#endif
+ do ipoin1D=1,NGLOB1D_RADIAL
+ array_val(iboolcorner(ipoin1D,icount_corners)) = array_val(iboolcorner(ipoin1D,icount_corners)) + &
+ buffer_recv_chunkcorners_scalar(ipoin1D)
+ enddo
+
+! receive from worker #2 and add to local array
+ if(NCHUNKS /= 2) then
+ sender = iproc_worker2_corners(imsg)
+#ifdef USE_MPI
+ call MPI_RECV(buffer_recv_chunkcorners_scalar,NGLOB1D_RADIAL, &
+ CUSTOM_MPI_TYPE,sender,itag,MPI_COMM_WORLD,msg_status,ier)
+#endif
+ do ipoin1D=1,NGLOB1D_RADIAL
+ array_val(iboolcorner(ipoin1D,icount_corners)) = array_val(iboolcorner(ipoin1D,icount_corners)) + &
+ buffer_recv_chunkcorners_scalar(ipoin1D)
+ enddo
+ endif
+
+ endif
+
+!---- send messages from the two workers to the master
+ if(myrank==iproc_worker1_corners(imsg) .or. &
+ (NCHUNKS /= 2 .and. myrank==iproc_worker2_corners(imsg))) then
+
+ receiver = iproc_master_corners(imsg)
+ do ipoin1D=1,NGLOB1D_RADIAL
+ buffer_send_chunkcorners_scalar(ipoin1D) = array_val(iboolcorner(ipoin1D,icount_corners))
+ enddo
+#ifdef USE_MPI
+ call MPI_SEND(buffer_send_chunkcorners_scalar,NGLOB1D_RADIAL,CUSTOM_MPI_TYPE, &
+ receiver,itag,MPI_COMM_WORLD,ier)
+#endif
+
+ endif
+
+! *********************************************************************
+! transmit messages back in opposite direction (master -> two workers)
+! *********************************************************************
+
+!---- receive messages from the master on the two workers
+ if(myrank==iproc_worker1_corners(imsg) .or. &
+ (NCHUNKS /= 2 .and. myrank==iproc_worker2_corners(imsg))) then
+
+! receive from master and copy to local array
+ sender = iproc_master_corners(imsg)
+#ifdef USE_MPI
+ call MPI_RECV(buffer_recv_chunkcorners_scalar,NGLOB1D_RADIAL, &
+ CUSTOM_MPI_TYPE,sender,itag,MPI_COMM_WORLD,msg_status,ier)
+#endif
+ do ipoin1D=1,NGLOB1D_RADIAL
+ array_val(iboolcorner(ipoin1D,icount_corners)) = buffer_recv_chunkcorners_scalar(ipoin1D)
+ enddo
+
+ endif
+
+!---- send messages from the master to the two workers
+ if(myrank==iproc_master_corners(imsg)) then
+
+ do ipoin1D=1,NGLOB1D_RADIAL
+ buffer_send_chunkcorners_scalar(ipoin1D) = array_val(iboolcorner(ipoin1D,icount_corners))
+ enddo
+
+! send to worker #1
+ receiver = iproc_worker1_corners(imsg)
+#ifdef USE_MPI
+ call MPI_SEND(buffer_send_chunkcorners_scalar,NGLOB1D_RADIAL,CUSTOM_MPI_TYPE, &
+ receiver,itag,MPI_COMM_WORLD,ier)
+#endif
+
+! send to worker #2
+ if(NCHUNKS /= 2) then
+ receiver = iproc_worker2_corners(imsg)
+#ifdef USE_MPI
+ call MPI_SEND(buffer_send_chunkcorners_scalar,NGLOB1D_RADIAL,CUSTOM_MPI_TYPE, &
+ receiver,itag,MPI_COMM_WORLD,ier)
+#endif
+ endif
+
+ endif
+
+ enddo
+
+ end subroutine assemble_MPI_scalar_block
+
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/attenuation_model.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/attenuation_model.F90 2008-09-11 21:43:27 UTC (rev 12866)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/attenuation_model.F90 2008-09-12 00:57:15 UTC (rev 12867)
@@ -83,7 +83,7 @@
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 :: Qone_minus_sum_beta, Qone_minus_sum_beta2 ! 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
@@ -284,7 +284,7 @@
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 :: Qone_minus_sum_beta, Qone_minus_sum_beta2 ! 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
@@ -437,7 +437,7 @@
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 :: Qone_minus_sum_beta, Qone_minus_sum_beta2 ! 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
@@ -512,7 +512,7 @@
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 :: Qone_minus_sum_beta, Qone_minus_sum_beta2 ! 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
@@ -667,7 +667,7 @@
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 :: Qone_minus_sum_beta, Qone_minus_sum_beta2 ! 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
@@ -750,11 +750,11 @@
AM_V%Qr(:) = AM_V%Qr(:) / R_EARTH
allocate(AM_V%Qsf(AM_V%Qn))
- allocate(AM_V%Qomsb(AM_V%Qn))
+ allocate(AM_V%Qone_minus_sum_beta(AM_V%Qn))
allocate(AM_V%Qfc(N_SLS,AM_V%Qn))
allocate(AM_V%Qsf2(AM_V%Qn))
- allocate(AM_V%Qomsb2(AM_V%Qn))
+ allocate(AM_V%Qone_minus_sum_beta2(AM_V%Qn))
allocate(AM_V%Qfc2(N_SLS,AM_V%Qn))
allocate(AM_V%interval_Q(AM_V%Qn))
@@ -764,18 +764,18 @@
do i = 1,AM_V%Qn
if(AM_V%Qmu(i) == 0.0d0) then
- AM_V%Qomsb(i) = 0.0d0
+ AM_V%Qone_minus_sum_beta(i) = 0.0d0
AM_V%Qfc(:,i) = 0.0d0
AM_V%Qsf(i) = 0.0d0
else
- call attenuation_property_values(tau_s, AM_V%Qtau_e(:,i), AM_V%Qfc(:,i), AM_V%Qomsb(i))
+ call attenuation_property_values(tau_s, AM_V%Qtau_e(:,i), AM_V%Qfc(:,i), AM_V%Qone_minus_sum_beta(i))
call attenuation_scale_factor(myrank, AM_V%QT_c_source, AM_V%Qtau_e(:,i), tau_s, AM_V%Qmu(i), AM_V%Qsf(i))
endif
enddo
! Determine the Spline Coefficients or Second Derivatives
call pspline_construction(AM_V%Qr, AM_V%Qsf, AM_V%Qn, Qp1, Qpn, AM_V%Qsf2, AM_V%interval_Q)
- call pspline_construction(AM_V%Qr, AM_V%Qomsb, AM_V%Qn, Qp1, Qpn, AM_V%Qomsb2, AM_V%interval_Q)
+ call pspline_construction(AM_V%Qr, AM_V%Qone_minus_sum_beta, AM_V%Qn, Qp1, Qpn, AM_V%Qone_minus_sum_beta2, AM_V%interval_Q)
do i = 1,N_SLS
! copy the sub-arrays to temporary arrays to avoid a warning by some compilers
! about temporary arrays being created automatically when using this expression
@@ -794,7 +794,8 @@
do i = 1,rmax
call attenuation_lookup_value(i, radius)
call pspline_evaluation(AM_V%Qr, AM_V%Qsf, AM_V%Qsf2, AM_V%Qn, radius, scale_factor(1,1,1,i), AM_V%interval_Q)
- call pspline_evaluation(AM_V%Qr, AM_V%Qomsb, AM_V%Qomsb2, AM_V%Qn, radius, one_minus_sum_beta(1,1,1,i), AM_V%interval_Q)
+ call pspline_evaluation(AM_V%Qr, AM_V%Qone_minus_sum_beta, AM_V%Qone_minus_sum_beta2, &
+ AM_V%Qn, radius, one_minus_sum_beta(1,1,1,i), AM_V%interval_Q)
do j = 1,N_SLS
Qfctmp = AM_V%Qfc(j,:)
Qfc2tmp = AM_V%Qfc2(j,:)
@@ -812,10 +813,10 @@
deallocate(AM_V%Qfc2)
deallocate(AM_V%Qsf2)
- deallocate(AM_V%Qomsb2)
+ deallocate(AM_V%Qone_minus_sum_beta2)
deallocate(AM_V%Qfc)
deallocate(AM_V%Qsf)
- deallocate(AM_V%Qomsb)
+ deallocate(AM_V%Qone_minus_sum_beta)
deallocate(AM_V%Qtau_e)
deallocate(Qfctmp)
deallocate(Qfc2tmp)
@@ -839,7 +840,7 @@
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 :: Qone_minus_sum_beta, Qone_minus_sum_beta2 ! 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
@@ -913,7 +914,7 @@
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 :: Qone_minus_sum_beta, Qone_minus_sum_beta2 ! 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
@@ -994,14 +995,14 @@
include 'constants.h'
integer myrank, vnspec
- double precision, dimension(NGLLX,NGLLY,NGLLZ,vnspec) :: one_minus_sum_beta, scale_factor
+ double precision, dimension(NGLLX,NGLLY,NGLLZ,vnspec) :: one_minus_sum_beta_array, scale_factor
double precision, dimension(N_SLS,NGLLX,NGLLY,NGLLZ,vnspec) :: factor_common
double precision, dimension(N_SLS) :: tau_s
integer i,j,k,ispec
double precision, dimension(N_SLS) :: tau_e, fc
- double precision omsb, Q_mu, sf, T_c_source, scale_t
+ double precision one_minus_sum_beta, Q_mu, sf, T_c_source, scale_t
! All of the following reads use the output parameters as their temporary arrays
! use the filename to determine the actual contents of the read
@@ -1032,10 +1033,10 @@
Q_mu = scale_factor(i,j,k,ispec)
! Determine the factor_common and one_minus_sum_beta from tau_s and tau_e
- call attenuation_property_values(tau_s, tau_e, fc, omsb)
+ call attenuation_property_values(tau_s, tau_e, fc, one_minus_sum_beta)
factor_common(:,i,j,k,ispec) = fc(:)
- one_minus_sum_beta(i,j,k,ispec) = omsb
+ one_minus_sum_beta_array(i,j,k,ispec) = one_minus_sum_beta
! Determine the "scale_factor" from tau_s, tau_e, central source frequency, and Q
call attenuation_scale_factor(myrank, T_c_source, tau_e, tau_s, Q_mu, sf)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_element_properties.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_element_properties.f90 2008-09-11 21:43:27 UTC (rev 12866)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_element_properties.f90 2008-09-12 00:57:15 UTC (rev 12867)
@@ -73,7 +73,7 @@
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 :: Qone_minus_sum_beta, Qone_minus_sum_beta2 ! 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
Copied: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_forces_CM_IC.f90 (from rev 12862, seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_forces_crust_mantle.f90)
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_forces_CM_IC.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_forces_CM_IC.f90 2008-09-12 00:57:15 UTC (rev 12867)
@@ -0,0 +1,959 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 4 . 1
+! --------------------------------------------------
+!
+! 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
+! August 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.
+!
+!=====================================================================
+
+! compute forces in the crust/mantle and in the inner core (i.e., all the solid regions)
+ subroutine compute_forces_CM_IC(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, &
+ kappavstore_crust_mantle,kappahstore_crust_mantle,muvstore_crust_mantle,muhstore_crust_mantle, &
+ eta_anisostore_crust_mantle, &
+ ibool_crust_mantle,idoubling_crust_mantle,R_memory_crust_mantle,epsilondev_crust_mantle,one_minus_sum_beta_crust_mantle, &
+ factor_common_crust_mantle,vx_crust_mantle,vy_crust_mantle,vz_crust_mantle,vnspec_crust_mantle, &
+ is_on_a_slice_edge_crust_mantle, &
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ displ_inner_core,accel_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, &
+ kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
+ R_memory_inner_core,epsilondev_inner_core,&
+ one_minus_sum_beta_inner_core,factor_common_inner_core, &
+ vx_inner_core,vy_inner_core,vz_inner_core,vnspec_inner_core,is_on_a_slice_edge_inner_core, &
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ hprime_xx,hprime_yy,hprime_zz, &
+ hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
+ alphaval,betaval,gammaval, &
+ COMPUTE_AND_STORE_STRAIN,AM_V,icall)
+
+ implicit none
+
+ include "constants.h"
+
+! include values created by the mesher
+! done for performance only using static allocation to allow for loop unrolling
+ include "values_from_mesher.h"
+
+ integer :: icall
+
+ logical, dimension(NSPEC_CRUST_MANTLE) :: is_on_a_slice_edge_crust_mantle
+
+! attenuation_model_variables
+ type attenuation_model_variables
+ sequence
+ double precision min_period, max_period
+ double precision :: QT_c_source ! Source Frequency
+ double precision, dimension(N_SLS) :: 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 :: Qone_minus_sum_beta, Qone_minus_sum_beta2 ! 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_crust_mantle
+
+! displacement and acceleration
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CRUST_MANTLE) :: displ_crust_mantle,accel_crust_mantle
+
+! 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_crust_mantle,vy_crust_mantle,vz_crust_mantle,vnspec_crust_mantle
+
+ real(kind=CUSTOM_REAL) one_minus_sum_beta_use,minus_sum_beta
+ real(kind=CUSTOM_REAL), dimension(vx_crust_mantle,vy_crust_mantle,vz_crust_mantle,vnspec_crust_mantle) :: &
+ one_minus_sum_beta_crust_mantle
+
+ 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_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: epsilondev_crust_mantle
+
+! [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_crust_mantle,vy_crust_mantle,vz_crust_mantle,vnspec_crust_mantle) :: &
+ factor_common_crust_mantle
+
+ real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc
+
+! arrays with mesh parameters per slice
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: ibool_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_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
+
+! 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(NGLLY,NGLLY) :: hprime_yy,hprimewgll_yy
+ real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ) :: hprime_zz,hprimewgll_zz
+ 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_crust_mantle,ystore_crust_mantle,zstore_crust_mantle
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPECMAX_ISO_MANTLE) :: &
+ kappavstore_crust_mantle,muvstore_crust_mantle
+
+! store anisotropic properties only where needed to save memory
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPECMAX_TISO_MANTLE) :: &
+ kappahstore_crust_mantle,muhstore_crust_mantle,eta_anisostore_crust_mantle
+
+ integer ispec,iglob
+ 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
+
+ real(kind=CUSTOM_REAL) tempx1l,tempx2l,tempx3l
+ real(kind=CUSTOM_REAL) tempy1l,tempy2l,tempy3l
+ real(kind=CUSTOM_REAL) tempz1l,tempz2l,tempz3l
+
+ real(kind=CUSTOM_REAL) radius_cr
+
+!=====================================================================
+!=====================================================================
+!=====================================================================
+
+ logical, dimension(NSPEC_OUTER_CORE) :: is_on_a_slice_edge_inner_core
+
+! same attenuation everywhere in the inner core therefore no need to use Brian's routines
+ integer, parameter :: iregion_selected_inner_core = 1
+
+! displacement and acceleration
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_INNER_CORE) :: displ_inner_core,accel_inner_core
+
+! variable lengths for factor_common and one_minus_sum_beta
+ integer vx_inner_core, vy_inner_core, vz_inner_core, vnspec_inner_core
+
+ real(kind=CUSTOM_REAL), dimension(vx_inner_core,vy_inner_core,vz_inner_core,vnspec_inner_core) :: one_minus_sum_beta_inner_core
+
+ real(kind=CUSTOM_REAL), dimension(N_SLS,vx_inner_core,vy_inner_core,vz_inner_core,vnspec_inner_core) :: factor_common_inner_core
+
+ real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ATTENUATION) :: R_memory_inner_core
+ real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: epsilondev_inner_core
+
+! array with the local to global mapping per slice
+ integer, dimension(NSPEC_INNER_CORE) :: idoubling_inner_core
+
+! arrays with mesh parameters per slice
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: ibool_inner_core
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_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
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: kappavstore_inner_core,muvstore_inner_core
+
+! ****************************************************
+! big loop over all spectral elements in the solid
+! ****************************************************
+
+! set acceleration to zero
+ if(icall == 1) accel_crust_mantle(:,:) = 0._CUSTOM_REAL
+
+ do ispec = 1,NSPEC_CRUST_MANTLE
+
+! hide communications by computing the edges first
+ if((icall == 2 .and. is_on_a_slice_edge_crust_mantle(ispec)) .or. &
+ (icall == 1 .and. .not. is_on_a_slice_edge_crust_mantle(ispec))) cycle
+
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+
+ tempx1l = 0._CUSTOM_REAL
+ tempx2l = 0._CUSTOM_REAL
+ tempx3l = 0._CUSTOM_REAL
+
+ tempy1l = 0._CUSTOM_REAL
+ tempy2l = 0._CUSTOM_REAL
+ tempy3l = 0._CUSTOM_REAL
+
+ tempz1l = 0._CUSTOM_REAL
+ tempz2l = 0._CUSTOM_REAL
+ tempz3l = 0._CUSTOM_REAL
+
+ do l=1,NGLLX
+ hp1 = hprime_xx(i,l)
+ iglob = ibool_crust_mantle(l,j,k,ispec)
+ tempx1l = tempx1l + displ_crust_mantle(1,iglob)*hp1
+ tempy1l = tempy1l + displ_crust_mantle(2,iglob)*hp1
+ tempz1l = tempz1l + displ_crust_mantle(3,iglob)*hp1
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLY
+ hp2 = hprime_yy(j,l)
+ iglob = ibool_crust_mantle(i,l,k,ispec)
+ tempx2l = tempx2l + displ_crust_mantle(1,iglob)*hp2
+ tempy2l = tempy2l + displ_crust_mantle(2,iglob)*hp2
+ tempz2l = tempz2l + displ_crust_mantle(3,iglob)*hp2
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLZ
+ hp3 = hprime_zz(k,l)
+ iglob = ibool_crust_mantle(i,j,l,ispec)
+ tempx3l = tempx3l + displ_crust_mantle(1,iglob)*hp3
+ tempy3l = tempy3l + displ_crust_mantle(2,iglob)*hp3
+ tempz3l = tempz3l + displ_crust_mantle(3,iglob)*hp3
+ enddo
+
+! get derivatives of ux, uy and uz with respect to x, y and z
+
+ xixl = xix_crust_mantle(i,j,k,ispec)
+ xiyl = xiy_crust_mantle(i,j,k,ispec)
+ xizl = xiz_crust_mantle(i,j,k,ispec)
+ etaxl = etax_crust_mantle(i,j,k,ispec)
+ etayl = etay_crust_mantle(i,j,k,ispec)
+ etazl = etaz_crust_mantle(i,j,k,ispec)
+ gammaxl = gammax_crust_mantle(i,j,k,ispec)
+ gammayl = gammay_crust_mantle(i,j,k,ispec)
+ gammazl = gammaz_crust_mantle(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*tempx1l + etaxl*tempx2l + gammaxl*tempx3l
+ duxdyl = xiyl*tempx1l + etayl*tempx2l + gammayl*tempx3l
+ duxdzl = xizl*tempx1l + etazl*tempx2l + gammazl*tempx3l
+
+ duydxl = xixl*tempy1l + etaxl*tempy2l + gammaxl*tempy3l
+ duydyl = xiyl*tempy1l + etayl*tempy2l + gammayl*tempy3l
+ duydzl = xizl*tempy1l + etazl*tempy2l + gammazl*tempy3l
+
+ duzdxl = xixl*tempz1l + etaxl*tempz2l + gammaxl*tempz3l
+ duzdyl = xiyl*tempz1l + etayl*tempz2l + gammayl*tempz3l
+ duzdzl = xizl*tempz1l + etazl*tempz2l + gammazl*tempz3l
+
+! 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
+ epsilondev_loc(1,i,j,k) = duxdxl - ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ epsilondev_loc(2,i,j,k) = duydyl - ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ 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
+ radius_cr = xstore_crust_mantle(ibool_crust_mantle(i,j,k,ispec))
+ call get_attenuation_index(idoubling_crust_mantle(ispec), dble(radius_cr), iregion_selected, .FALSE., AM_V)
+ one_minus_sum_beta_use = one_minus_sum_beta_crust_mantle(1,1,1,iregion_selected)
+ minus_sum_beta = one_minus_sum_beta_use - 1.0
+ endif
+
+!
+! compute either isotropic or anisotropic elements
+!
+
+ if(ANISOTROPIC_3D_MANTLE_VAL) then
+
+ else
+
+! do not use transverse isotropy except if element is between d220 and Moho
+ if(.not. (TRANSVERSE_ISOTROPY_VAL .and. (idoubling_crust_mantle(ispec)==IFLAG_220_80 .or. &
+ idoubling_crust_mantle(ispec)==IFLAG_80_MOHO))) then
+
+! layer with no transverse isotropy, use kappav and muv
+ kappal = kappavstore_crust_mantle(i,j,k,ispec)
+ mul = muvstore_crust_mantle(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_crust_mantle(i,j,k,ispec)
+ muvl = muvstore_crust_mantle(i,j,k,ispec)
+
+ kappahl = kappahstore_crust_mantle(i,j,k,ispec)
+ muhl = muhstore_crust_mantle(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_crust_mantle(i,j,k,ispec) !!! that is F / (A - 2 L)
+
+! use mesh coordinates to get theta and phi
+! ystore_crust_mantle and zstore_crust_mantle contain theta and phi
+
+ iglob = ibool_crust_mantle(i,j,k,ispec)
+ theta = ystore_crust_mantle(iglob)
+ phi = zstore_crust_mantle(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_crust_mantle(1,i_sls,i,j,k,ispec)
+ R_yy_val = R_memory_crust_mantle(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_crust_mantle(3,i_sls,i,j,k,ispec)
+ sigma_xz = sigma_xz - R_memory_crust_mantle(4,i_sls,i,j,k,ispec)
+ sigma_yz = sigma_yz - R_memory_crust_mantle(5,i_sls,i,j,k,ispec)
+ enddo
+ endif
+
+! form dot product with test vector, symmetric form
+ tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_xy*xiyl + sigma_xz*xizl)
+ tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_yz*xizl)
+ tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl)
+
+ tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_xy*etayl + sigma_xz*etazl)
+ tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_yz*etazl)
+ tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl)
+
+ tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_xy*gammayl + sigma_xz*gammazl)
+ tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_yz*gammazl)
+ tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl)
+
+ enddo
+ enddo
+ enddo
+
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+
+ tempx1l = 0._CUSTOM_REAL
+ tempy1l = 0._CUSTOM_REAL
+ tempz1l = 0._CUSTOM_REAL
+
+ tempx2l = 0._CUSTOM_REAL
+ tempy2l = 0._CUSTOM_REAL
+ tempz2l = 0._CUSTOM_REAL
+
+ tempx3l = 0._CUSTOM_REAL
+ tempy3l = 0._CUSTOM_REAL
+ tempz3l = 0._CUSTOM_REAL
+
+ do l=1,NGLLX
+ fac1 = hprimewgll_xx(l,i)
+ tempx1l = tempx1l + tempx1(l,j,k)*fac1
+ tempy1l = tempy1l + tempy1(l,j,k)*fac1
+ tempz1l = tempz1l + tempz1(l,j,k)*fac1
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLY
+ fac2 = hprimewgll_yy(l,j)
+ tempx2l = tempx2l + tempx2(i,l,k)*fac2
+ tempy2l = tempy2l + tempy2(i,l,k)*fac2
+ tempz2l = tempz2l + tempz2(i,l,k)*fac2
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLZ
+ fac3 = hprimewgll_zz(l,k)
+ tempx3l = tempx3l + tempx3(i,j,l)*fac3
+ tempy3l = tempy3l + tempy3(i,j,l)*fac3
+ tempz3l = tempz3l + tempz3(i,j,l)*fac3
+ enddo
+
+ fac1 = wgllwgll_yz(j,k)
+ fac2 = wgllwgll_xz(i,k)
+ fac3 = wgllwgll_xy(i,j)
+
+ sum_terms(1,i,j,k) = - (fac1*tempx1l + fac2*tempx2l + fac3*tempx3l)
+ sum_terms(2,i,j,k) = - (fac1*tempy1l + fac2*tempy2l + fac3*tempy3l)
+ sum_terms(3,i,j,k) = - (fac1*tempz1l + fac2*tempz2l + fac3*tempz3l)
+
+ 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_crust_mantle(i,j,k,ispec)
+ accel_crust_mantle(1,iglob) = accel_crust_mantle(1,iglob) + sum_terms(1,i,j,k)
+ accel_crust_mantle(2,iglob) = accel_crust_mantle(2,iglob) + sum_terms(2,i,j,k)
+ accel_crust_mantle(3,iglob) = accel_crust_mantle(3,iglob) + sum_terms(3,i,j,k)
+ enddo
+ enddo
+ enddo
+
+! update memory variables based upon a 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
+ if(ATTENUATION_VAL) then
+ do k = 1,NGLLZ
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+ 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
+ R_memory_crust_mantle(i_memory,i_sls,i,j,k,ispec) = alphaval(i_sls) * &
+ R_memory_crust_mantle(i_memory,i_sls,i,j,k,ispec) + &
+ factor_common_crust_mantle(i_sls,1,1,1,iregion_selected) * muvstore_crust_mantle(i,j,k,ispec) * &
+ (betaval(i_sls) * epsilondev_crust_mantle(i_memory,i,j,k,ispec) + &
+ gammaval(i_sls) * epsilondev_loc(i_memory,i,j,k))
+ enddo
+ enddo
+ enddo
+ enddo
+ enddo
+ endif
+
+! save deviatoric strain for Runge-Kutta scheme
+ if(COMPUTE_AND_STORE_STRAIN) epsilondev_crust_mantle(:,:,:,:,ispec) = epsilondev_loc(:,:,:,:)
+
+ enddo ! spectral element loop
+
+!=====================================================================
+!=====================================================================
+!=====================================================================
+
+! ****************************************************
+! big loop over all spectral elements in the solid
+! ****************************************************
+
+! set acceleration to zero
+ if(icall == 1) accel_inner_core(:,:) = 0._CUSTOM_REAL
+
+ do ispec = 1,NSPEC_INNER_CORE
+
+! hide communications by computing the edges first
+ if((icall == 2 .and. is_on_a_slice_edge_inner_core(ispec)) .or. &
+ (icall == 1 .and. .not. is_on_a_slice_edge_inner_core(ispec))) cycle
+
+! exclude fictitious elements in central cube
+ if(idoubling_inner_core(ispec) /= IFLAG_IN_FICTITIOUS_CUBE) then
+
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+
+ tempx1l = 0._CUSTOM_REAL
+ tempx2l = 0._CUSTOM_REAL
+ tempx3l = 0._CUSTOM_REAL
+
+ tempy1l = 0._CUSTOM_REAL
+ tempy2l = 0._CUSTOM_REAL
+ tempy3l = 0._CUSTOM_REAL
+
+ tempz1l = 0._CUSTOM_REAL
+ tempz2l = 0._CUSTOM_REAL
+ tempz3l = 0._CUSTOM_REAL
+
+ do l=1,NGLLX
+ hp1 = hprime_xx(i,l)
+ iglob = ibool_inner_core(l,j,k,ispec)
+ tempx1l = tempx1l + displ_inner_core(1,iglob)*hp1
+ tempy1l = tempy1l + displ_inner_core(2,iglob)*hp1
+ tempz1l = tempz1l + displ_inner_core(3,iglob)*hp1
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLY
+ hp2 = hprime_yy(j,l)
+ iglob = ibool_inner_core(i,l,k,ispec)
+ tempx2l = tempx2l + displ_inner_core(1,iglob)*hp2
+ tempy2l = tempy2l + displ_inner_core(2,iglob)*hp2
+ tempz2l = tempz2l + displ_inner_core(3,iglob)*hp2
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLZ
+ hp3 = hprime_zz(k,l)
+ iglob = ibool_inner_core(i,j,l,ispec)
+ tempx3l = tempx3l + displ_inner_core(1,iglob)*hp3
+ tempy3l = tempy3l + displ_inner_core(2,iglob)*hp3
+ tempz3l = tempz3l + displ_inner_core(3,iglob)*hp3
+ enddo
+
+! get derivatives of ux, uy and uz with respect to x, y and z
+
+ xixl = xix_inner_core(i,j,k,ispec)
+ xiyl = xiy_inner_core(i,j,k,ispec)
+ xizl = xiz_inner_core(i,j,k,ispec)
+ etaxl = etax_inner_core(i,j,k,ispec)
+ etayl = etay_inner_core(i,j,k,ispec)
+ etazl = etaz_inner_core(i,j,k,ispec)
+ gammaxl = gammax_inner_core(i,j,k,ispec)
+ gammayl = gammay_inner_core(i,j,k,ispec)
+ gammazl = gammaz_inner_core(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*tempx1l + etaxl*tempx2l + gammaxl*tempx3l
+ duxdyl = xiyl*tempx1l + etayl*tempx2l + gammayl*tempx3l
+ duxdzl = xizl*tempx1l + etazl*tempx2l + gammazl*tempx3l
+
+ duydxl = xixl*tempy1l + etaxl*tempy2l + gammaxl*tempy3l
+ duydyl = xiyl*tempy1l + etayl*tempy2l + gammayl*tempy3l
+ duydzl = xizl*tempy1l + etazl*tempy2l + gammazl*tempy3l
+
+ duzdxl = xixl*tempz1l + etaxl*tempz2l + gammaxl*tempz3l
+ duzdyl = xiyl*tempz1l + etayl*tempz2l + gammayl*tempz3l
+ duzdzl = xizl*tempz1l + etazl*tempz2l + gammazl*tempz3l
+
+! 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
+ epsilondev_loc(1,i,j,k) = duxdxl - ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ epsilondev_loc(2,i,j,k) = duydyl - ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ 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
+! same attenuation everywhere in the inner core therefore no need to use Brian's routines
+!!!!! radius_cr = xstore(ibool_inner_core(i,j,k,ispec))
+!!!!! call get_attenuation_index(idoubling_inner_core(ispec), dble(radius_cr), iregion_selected_inner_core, .TRUE., AM_V)
+ minus_sum_beta = one_minus_sum_beta_inner_core(1,1,1,iregion_selected_inner_core) - 1.0
+ endif ! ATTENUATION_VAL
+
+ if(ANISOTROPIC_INNER_CORE_VAL) then
+
+ else
+
+! inner core with no anisotropy, use kappav and muv for instance
+! layer with no anisotropy, use kappav and muv for instance
+ kappal = kappavstore_inner_core(i,j,k,ispec)
+ mul = muvstore_inner_core(i,j,k,ispec)
+
+! use unrelaxed parameters if attenuation
+ if(ATTENUATION_VAL) then
+ mul = mul * one_minus_sum_beta_inner_core(1,1,1,iregion_selected_inner_core)
+ 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_inner_core(1,i_sls,i,j,k,ispec)
+ R_yy_val = R_memory_inner_core(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_inner_core(3,i_sls,i,j,k,ispec)
+ sigma_xz = sigma_xz - R_memory_inner_core(4,i_sls,i,j,k,ispec)
+ sigma_yz = sigma_yz - R_memory_inner_core(5,i_sls,i,j,k,ispec)
+ enddo
+ endif
+
+! form dot product with test vector, symmetric form
+ tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_xy*xiyl + sigma_xz*xizl)
+ tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_yz*xizl)
+ tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl)
+
+ tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_xy*etayl + sigma_xz*etazl)
+ tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_yz*etazl)
+ tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl)
+
+ tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_xy*gammayl + sigma_xz*gammazl)
+ tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_yz*gammazl)
+ tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl)
+
+ enddo
+ enddo
+ enddo
+
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+
+ tempx1l = 0._CUSTOM_REAL
+ tempy1l = 0._CUSTOM_REAL
+ tempz1l = 0._CUSTOM_REAL
+
+ tempx2l = 0._CUSTOM_REAL
+ tempy2l = 0._CUSTOM_REAL
+ tempz2l = 0._CUSTOM_REAL
+
+ tempx3l = 0._CUSTOM_REAL
+ tempy3l = 0._CUSTOM_REAL
+ tempz3l = 0._CUSTOM_REAL
+
+ do l=1,NGLLX
+ fac1 = hprimewgll_xx(l,i)
+ tempx1l = tempx1l + tempx1(l,j,k)*fac1
+ tempy1l = tempy1l + tempy1(l,j,k)*fac1
+ tempz1l = tempz1l + tempz1(l,j,k)*fac1
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLY
+ fac2 = hprimewgll_yy(l,j)
+ tempx2l = tempx2l + tempx2(i,l,k)*fac2
+ tempy2l = tempy2l + tempy2(i,l,k)*fac2
+ tempz2l = tempz2l + tempz2(i,l,k)*fac2
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLZ
+ fac3 = hprimewgll_zz(l,k)
+ tempx3l = tempx3l + tempx3(i,j,l)*fac3
+ tempy3l = tempy3l + tempy3(i,j,l)*fac3
+ tempz3l = tempz3l + tempz3(i,j,l)*fac3
+ enddo
+
+ fac1 = wgllwgll_yz(j,k)
+ fac2 = wgllwgll_xz(i,k)
+ fac3 = wgllwgll_xy(i,j)
+
+ sum_terms(1,i,j,k) = - (fac1*tempx1l + fac2*tempx2l + fac3*tempx3l)
+ sum_terms(2,i,j,k) = - (fac1*tempy1l + fac2*tempy2l + fac3*tempy3l)
+ sum_terms(3,i,j,k) = - (fac1*tempz1l + fac2*tempz2l + fac3*tempz3l)
+
+ 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_inner_core(i,j,k,ispec)
+ accel_inner_core(:,iglob) = accel_inner_core(:,iglob) + sum_terms(:,i,j,k)
+ enddo
+ enddo
+ enddo
+
+! update memory variables based upon a 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
+ if(ATTENUATION_VAL) then
+ do k = 1,NGLLZ
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+ do i_sls = 1,N_SLS
+ do i_memory = 1,5
+ R_memory_inner_core(i_memory,i_sls,i,j,k,ispec) = &
+ alphaval(i_sls) * &
+ R_memory_inner_core(i_memory,i_sls,i,j,k,ispec) + muvstore_inner_core(i,j,k,ispec) * &
+ factor_common_inner_core(i_sls,1,1,1,iregion_selected_inner_core) * &
+ (betaval(i_sls) * &
+ epsilondev_inner_core(i_memory,i,j,k,ispec) + gammaval(i_sls) * epsilondev_loc(i_memory,i,j,k))
+ enddo
+ enddo
+ enddo
+ enddo
+ enddo
+ endif
+
+! save deviatoric strain for Runge-Kutta scheme
+ if(COMPUTE_AND_STORE_STRAIN) epsilondev_inner_core(:,:,:,:,ispec) = epsilondev_loc(:,:,:,:)
+
+ endif ! end test to exclude fictitious elements in central cube
+
+ enddo ! spectral element loop
+
+ end subroutine compute_forces_CM_IC
+
Copied: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_forces_OC.f90 (from rev 12862, seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_forces_outer_core.f90)
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_forces_OC.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_forces_OC.f90 2008-09-12 00:57:15 UTC (rev 12867)
@@ -0,0 +1,234 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 4 . 1
+! --------------------------------------------------
+!
+! 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
+! August 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.
+!
+!=====================================================================
+
+! compute forces in the outer core (i.e., the fluid region)
+ subroutine compute_forces_OC(d_ln_density_dr_table, &
+ displfluid,accelfluid,xstore,ystore,zstore, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ hprime_xx,hprime_yy,hprime_zz, &
+ hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,ibool,icall,is_on_a_slice_edge_outer_core)
+
+ implicit none
+
+ include "constants.h"
+
+! include values created by the mesher
+! done for performance only using static allocation to allow for loop unrolling
+ include "values_from_mesher.h"
+
+ logical, dimension(NSPEC_OUTER_CORE) :: is_on_a_slice_edge_outer_core
+
+ integer :: icall
+
+! displacement and acceleration
+ real(kind=CUSTOM_REAL), dimension(nglob_outer_core) :: displfluid,accelfluid
+
+! 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(NGLLY,NGLLY) :: hprime_yy,hprimewgll_yy
+ real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ) :: hprime_zz,hprimewgll_zz
+ 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
+
+! for gravity
+ integer int_radius
+ double precision radius,theta,phi
+ double precision cos_theta,sin_theta,cos_phi,sin_phi
+ double precision, dimension(NRAD_GRAVITY) :: d_ln_density_dr_table
+ real(kind=CUSTOM_REAL), dimension(nglob_outer_core) :: xstore,ystore,zstore
+
+ 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
+
+ double precision grad_x_ln_rho,grad_y_ln_rho,grad_z_ln_rho
+
+! ****************************************************
+! big loop over all spectral elements in the fluid
+! ****************************************************
+
+! set acceleration to zero
+ if(icall == 1) accelfluid(:) = 0._CUSTOM_REAL
+
+ do ispec = 1,NSPEC_OUTER_CORE
+
+! hide communications by computing the edges first
+ if((icall == 2 .and. is_on_a_slice_edge_outer_core(ispec)) .or. &
+ (icall == 1 .and. .not. is_on_a_slice_edge_outer_core(ispec))) cycle
+
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+
+ iglob = ibool(i,j,k,ispec)
+
+ tempx1l = 0._CUSTOM_REAL
+ tempx2l = 0._CUSTOM_REAL
+ tempx3l = 0._CUSTOM_REAL
+
+ do l=1,NGLLX
+ tempx1l = tempx1l + displfluid(ibool(l,j,k,ispec)) * hprime_xx(i,l)
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLY
+ tempx2l = tempx2l + displfluid(ibool(i,l,k,ispec)) * hprime_yy(j,l)
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLZ
+ tempx3l = tempx3l + displfluid(ibool(i,j,l,ispec)) * hprime_zz(k,l)
+ enddo
+
+! 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*tempx1l + etaxl*tempx2l + gammaxl*tempx3l
+ dpotentialdyl = xiyl*tempx1l + etayl*tempx2l + gammayl*tempx3l
+ dpotentialdzl = xizl*tempx1l + etazl*tempx2l + gammazl*tempx3l
+
+! add (chi/rho)grad(rho) term in no gravity case
+
+! 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
+
+ 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)
+
+! grad(rho)/rho in Cartesian components
+ grad_x_ln_rho = sin_theta * cos_phi * d_ln_density_dr_table(int_radius)
+ grad_y_ln_rho = sin_theta * sin_phi * d_ln_density_dr_table(int_radius)
+ grad_z_ln_rho = cos_theta * d_ln_density_dr_table(int_radius)
+
+! adding (chi/rho)grad(rho)
+ dpotentialdxl = dpotentialdxl + displfluid(iglob) * grad_x_ln_rho
+ dpotentialdyl = dpotentialdyl + displfluid(iglob) * grad_y_ln_rho
+ dpotentialdzl = dpotentialdzl + displfluid(iglob) * grad_z_ln_rho
+
+ tempx1(i,j,k) = jacobianl*(xixl*dpotentialdxl + xiyl*dpotentialdyl + xizl*dpotentialdzl)
+ tempx2(i,j,k) = jacobianl*(etaxl*dpotentialdxl + etayl*dpotentialdyl + etazl*dpotentialdzl)
+ tempx3(i,j,k) = jacobianl*(gammaxl*dpotentialdxl + gammayl*dpotentialdyl + gammazl*dpotentialdzl)
+
+ enddo
+ enddo
+ enddo
+
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+
+ tempx1l = 0._CUSTOM_REAL
+ tempx2l = 0._CUSTOM_REAL
+ tempx3l = 0._CUSTOM_REAL
+
+ do l=1,NGLLX
+ tempx1l = tempx1l + tempx1(l,j,k) * hprimewgll_xx(l,i)
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLY
+ tempx2l = tempx2l + tempx2(i,l,k) * hprimewgll_yy(l,j)
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
+
+!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLZ
+ tempx3l = tempx3l + tempx3(i,j,l) * hprimewgll_zz(l,k)
+ enddo
+
+! sum contributions from each element to the global mesh and add gravity term
+
+ iglob = ibool(i,j,k,ispec)
+ accelfluid(iglob) = accelfluid(iglob) - (wgllwgll_yz(j,k)*tempx1l + wgllwgll_xz(i,k)*tempx2l + wgllwgll_xy(i,j)*tempx3l)
+
+ enddo
+ enddo
+ enddo
+
+ enddo ! spectral element loop
+
+ end subroutine compute_forces_OC
+
Deleted: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_forces_crust_mantle.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_forces_crust_mantle.f90 2008-09-11 21:43:27 UTC (rev 12866)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_forces_crust_mantle.f90 2008-09-12 00:57:15 UTC (rev 12867)
@@ -1,631 +0,0 @@
-!=====================================================================
-!
-! S p e c f e m 3 D G l o b e V e r s i o n 4 . 1
-! --------------------------------------------------
-!
-! 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
-! August 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_crust_mantle(displ,accel,xstore,ystore,zstore, &
- xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
- hprime_xx,hprime_yy,hprime_zz, &
- hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
- kappavstore,kappahstore,muvstore,muhstore,eta_anisostore, &
- ibool,idoubling,R_memory,epsilondev,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 "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(N_SLS) :: 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(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc
-
-! 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(NGLLY,NGLLY) :: hprime_yy,hprimewgll_yy
- real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ) :: hprime_zz,hprimewgll_zz
- 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
-
- integer ispec,iglob
- 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
-
- real(kind=CUSTOM_REAL) tempx1l,tempx2l,tempx3l
- real(kind=CUSTOM_REAL) tempy1l,tempy2l,tempy3l
- real(kind=CUSTOM_REAL) tempz1l,tempz2l,tempz3l
-
- real(kind=CUSTOM_REAL) radius_cr
-
-! ****************************************************
-! big loop over all spectral elements in the solid
-! ****************************************************
-
-! set acceleration to zero
- accel(:,:) = 0._CUSTOM_REAL
-
- do ispec = 1,NSPEC_CRUST_MANTLE
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
-
- tempx1l = 0._CUSTOM_REAL
- tempx2l = 0._CUSTOM_REAL
- tempx3l = 0._CUSTOM_REAL
-
- tempy1l = 0._CUSTOM_REAL
- tempy2l = 0._CUSTOM_REAL
- tempy3l = 0._CUSTOM_REAL
-
- tempz1l = 0._CUSTOM_REAL
- tempz2l = 0._CUSTOM_REAL
- tempz3l = 0._CUSTOM_REAL
-
- do l=1,NGLLX
- hp1 = hprime_xx(i,l)
- iglob = ibool(l,j,k,ispec)
- tempx1l = tempx1l + displ(1,iglob)*hp1
- tempy1l = tempy1l + displ(2,iglob)*hp1
- tempz1l = tempz1l + displ(3,iglob)*hp1
-!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
-
-!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLY
- hp2 = hprime_yy(j,l)
- iglob = ibool(i,l,k,ispec)
- tempx2l = tempx2l + displ(1,iglob)*hp2
- tempy2l = tempy2l + displ(2,iglob)*hp2
- tempz2l = tempz2l + displ(3,iglob)*hp2
-!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
-
-!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLZ
- hp3 = hprime_zz(k,l)
- iglob = ibool(i,j,l,ispec)
- tempx3l = tempx3l + displ(1,iglob)*hp3
- tempy3l = tempy3l + displ(2,iglob)*hp3
- tempz3l = tempz3l + displ(3,iglob)*hp3
- enddo
-
-! 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*tempx1l + etaxl*tempx2l + gammaxl*tempx3l
- duxdyl = xiyl*tempx1l + etayl*tempx2l + gammayl*tempx3l
- duxdzl = xizl*tempx1l + etazl*tempx2l + gammazl*tempx3l
-
- duydxl = xixl*tempy1l + etaxl*tempy2l + gammaxl*tempy3l
- duydyl = xiyl*tempy1l + etayl*tempy2l + gammayl*tempy3l
- duydzl = xizl*tempy1l + etazl*tempy2l + gammazl*tempy3l
-
- duzdxl = xixl*tempz1l + etaxl*tempz2l + gammaxl*tempz3l
- duzdyl = xiyl*tempz1l + etayl*tempz2l + gammayl*tempz3l
- duzdzl = xizl*tempz1l + etazl*tempz2l + gammazl*tempz3l
-
-! 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
- epsilondev_loc(1,i,j,k) = duxdxl - ONE_THIRD * (duxdxl + duydyl + duzdzl)
- epsilondev_loc(2,i,j,k) = duydyl - ONE_THIRD * (duxdxl + duydyl + duzdzl)
- 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
- 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)
- minus_sum_beta = one_minus_sum_beta_use - 1.0
- endif
-
-!
-! compute either isotropic or anisotropic elements
-!
-
- if(ANISOTROPIC_3D_MANTLE_VAL) then
-
- 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
-
-! form dot product with test vector, symmetric form
- tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_xy*xiyl + sigma_xz*xizl)
- tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_yz*xizl)
- tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl)
-
- tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_xy*etayl + sigma_xz*etazl)
- tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_yz*etazl)
- tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl)
-
- tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_xy*gammayl + sigma_xz*gammazl)
- tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_yz*gammazl)
- tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl)
-
- enddo
- enddo
- enddo
-
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
-
- tempx1l = 0._CUSTOM_REAL
- tempy1l = 0._CUSTOM_REAL
- tempz1l = 0._CUSTOM_REAL
-
- tempx2l = 0._CUSTOM_REAL
- tempy2l = 0._CUSTOM_REAL
- tempz2l = 0._CUSTOM_REAL
-
- tempx3l = 0._CUSTOM_REAL
- tempy3l = 0._CUSTOM_REAL
- tempz3l = 0._CUSTOM_REAL
-
- do l=1,NGLLX
- fac1 = hprimewgll_xx(l,i)
- tempx1l = tempx1l + tempx1(l,j,k)*fac1
- tempy1l = tempy1l + tempy1(l,j,k)*fac1
- tempz1l = tempz1l + tempz1(l,j,k)*fac1
-!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
-
-!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLY
- fac2 = hprimewgll_yy(l,j)
- tempx2l = tempx2l + tempx2(i,l,k)*fac2
- tempy2l = tempy2l + tempy2(i,l,k)*fac2
- tempz2l = tempz2l + tempz2(i,l,k)*fac2
-!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
-
-!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLZ
- fac3 = hprimewgll_zz(l,k)
- tempx3l = tempx3l + tempx3(i,j,l)*fac3
- tempy3l = tempy3l + tempy3(i,j,l)*fac3
- tempz3l = tempz3l + tempz3(i,j,l)*fac3
- enddo
-
- fac1 = wgllwgll_yz(j,k)
- fac2 = wgllwgll_xz(i,k)
- fac3 = wgllwgll_xy(i,j)
-
- sum_terms(1,i,j,k) = - (fac1*tempx1l + fac2*tempx2l + fac3*tempx3l)
- sum_terms(2,i,j,k) = - (fac1*tempy1l + fac2*tempy2l + fac3*tempy3l)
- sum_terms(3,i,j,k) = - (fac1*tempz1l + fac2*tempz2l + fac3*tempz3l)
-
- 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(1,iglob) = accel(1,iglob) + sum_terms(1,i,j,k)
- accel(2,iglob) = accel(2,iglob) + sum_terms(2,i,j,k)
- accel(3,iglob) = accel(3,iglob) + sum_terms(3,i,j,k)
- enddo
- enddo
- enddo
-
-! update memory variables based upon a 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
- if(ATTENUATION_VAL) then
- do k = 1,NGLLZ
- do j = 1,NGLLY
- do i = 1,NGLLX
- 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
- R_memory(i_memory,i_sls,i,j,k,ispec) = alphaval(i_sls) * &
- R_memory(i_memory,i_sls,i,j,k,ispec) + &
- factor_common(i_sls,1,1,1,iregion_selected) * muvstore(i,j,k,ispec) * &
- (betaval(i_sls) * epsilondev(i_memory,i,j,k,ispec) + &
- gammaval(i_sls) * epsilondev_loc(i_memory,i,j,k))
- enddo
- enddo
- enddo
- enddo
- enddo
- endif
-
-! save deviatoric strain for Runge-Kutta scheme
- if(COMPUTE_AND_STORE_STRAIN) epsilondev(:,:,:,:,ispec) = epsilondev_loc(:,:,:,:)
-
- enddo ! spectral element loop
-
- end subroutine compute_forces_crust_mantle
-
Deleted: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_forces_inner_core.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_forces_inner_core.f90 2008-09-11 21:43:27 UTC (rev 12866)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_forces_inner_core.f90 2008-09-12 00:57:15 UTC (rev 12867)
@@ -1,376 +0,0 @@
-!=====================================================================
-!
-! S p e c f e m 3 D G l o b e V e r s i o n 4 . 1
-! --------------------------------------------------
-!
-! 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
-! August 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_inner_core(displ,accel, &
- xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
- hprime_xx,hprime_yy,hprime_zz, &
- hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
- kappavstore,muvstore,ibool,idoubling, &
- R_memory,epsilondev,&
- one_minus_sum_beta,alphaval,betaval,gammaval,factor_common, &
- vx,vy,vz,vnspec,COMPUTE_AND_STORE_STRAIN)
-
- implicit none
-
- include "constants.h"
-
-! include values created by the mesher
-! done for performance only using static allocation to allow for loop unrolling
- include "values_from_mesher.h"
-
-! same attenuation everywhere in the inner core therefore no need to use Brian's routines
- integer, parameter :: iregion_selected = 1
-
-! 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(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
-
-! 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(NGLLY,NGLLY) :: hprime_yy,hprimewgll_yy
- real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ) :: hprime_zz,hprimewgll_zz
- 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
-
- 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) 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) tempx1l,tempx2l,tempx3l
- real(kind=CUSTOM_REAL) tempy1l,tempy2l,tempy3l
- real(kind=CUSTOM_REAL) tempz1l,tempz2l,tempz3l
-
-! ****************************************************
-! 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
-
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
-
- tempx1l = 0._CUSTOM_REAL
- tempx2l = 0._CUSTOM_REAL
- tempx3l = 0._CUSTOM_REAL
-
- tempy1l = 0._CUSTOM_REAL
- tempy2l = 0._CUSTOM_REAL
- tempy3l = 0._CUSTOM_REAL
-
- tempz1l = 0._CUSTOM_REAL
- tempz2l = 0._CUSTOM_REAL
- tempz3l = 0._CUSTOM_REAL
-
- do l=1,NGLLX
- hp1 = hprime_xx(i,l)
- iglob = ibool(l,j,k,ispec)
- tempx1l = tempx1l + displ(1,iglob)*hp1
- tempy1l = tempy1l + displ(2,iglob)*hp1
- tempz1l = tempz1l + displ(3,iglob)*hp1
-!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
-
-!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLY
- hp2 = hprime_yy(j,l)
- iglob = ibool(i,l,k,ispec)
- tempx2l = tempx2l + displ(1,iglob)*hp2
- tempy2l = tempy2l + displ(2,iglob)*hp2
- tempz2l = tempz2l + displ(3,iglob)*hp2
-!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
-
-!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLZ
- hp3 = hprime_zz(k,l)
- iglob = ibool(i,j,l,ispec)
- tempx3l = tempx3l + displ(1,iglob)*hp3
- tempy3l = tempy3l + displ(2,iglob)*hp3
- tempz3l = tempz3l + displ(3,iglob)*hp3
- enddo
-
-! 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*tempx1l + etaxl*tempx2l + gammaxl*tempx3l
- duxdyl = xiyl*tempx1l + etayl*tempx2l + gammayl*tempx3l
- duxdzl = xizl*tempx1l + etazl*tempx2l + gammazl*tempx3l
-
- duydxl = xixl*tempy1l + etaxl*tempy2l + gammaxl*tempy3l
- duydyl = xiyl*tempy1l + etayl*tempy2l + gammayl*tempy3l
- duydzl = xizl*tempy1l + etazl*tempy2l + gammazl*tempy3l
-
- duzdxl = xixl*tempz1l + etaxl*tempz2l + gammaxl*tempz3l
- duzdyl = xiyl*tempz1l + etayl*tempz2l + gammayl*tempz3l
- duzdzl = xizl*tempz1l + etazl*tempz2l + gammazl*tempz3l
-
-! 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
- epsilondev_loc(1,i,j,k) = duxdxl - ONE_THIRD * (duxdxl + duydyl + duzdzl)
- epsilondev_loc(2,i,j,k) = duydyl - ONE_THIRD * (duxdxl + duydyl + duzdzl)
- 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
-! same attenuation everywhere in the inner core therefore no need to use Brian's routines
-!!!!! 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_VAL
-
- if(ANISOTROPIC_INNER_CORE_VAL) then
-
- 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
- mul = mul * one_minus_sum_beta(1,1,1,iregion_selected)
- 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
-
-! form dot product with test vector, symmetric form
- tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_xy*xiyl + sigma_xz*xizl)
- tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_yz*xizl)
- tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl)
-
- tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_xy*etayl + sigma_xz*etazl)
- tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_yz*etazl)
- tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl)
-
- tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_xy*gammayl + sigma_xz*gammazl)
- tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_yz*gammazl)
- tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl)
-
- enddo
- enddo
- enddo
-
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
-
- tempx1l = 0._CUSTOM_REAL
- tempy1l = 0._CUSTOM_REAL
- tempz1l = 0._CUSTOM_REAL
-
- tempx2l = 0._CUSTOM_REAL
- tempy2l = 0._CUSTOM_REAL
- tempz2l = 0._CUSTOM_REAL
-
- tempx3l = 0._CUSTOM_REAL
- tempy3l = 0._CUSTOM_REAL
- tempz3l = 0._CUSTOM_REAL
-
- do l=1,NGLLX
- fac1 = hprimewgll_xx(l,i)
- tempx1l = tempx1l + tempx1(l,j,k)*fac1
- tempy1l = tempy1l + tempy1(l,j,k)*fac1
- tempz1l = tempz1l + tempz1(l,j,k)*fac1
-!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
-
-!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLY
- fac2 = hprimewgll_yy(l,j)
- tempx2l = tempx2l + tempx2(i,l,k)*fac2
- tempy2l = tempy2l + tempy2(i,l,k)*fac2
- tempz2l = tempz2l + tempz2(i,l,k)*fac2
-!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
-
-!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLZ
- fac3 = hprimewgll_zz(l,k)
- tempx3l = tempx3l + tempx3(i,j,l)*fac3
- tempy3l = tempy3l + tempy3(i,j,l)*fac3
- tempz3l = tempz3l + tempz3(i,j,l)*fac3
- enddo
-
- fac1 = wgllwgll_yz(j,k)
- fac2 = wgllwgll_xz(i,k)
- fac3 = wgllwgll_xy(i,j)
-
- sum_terms(1,i,j,k) = - (fac1*tempx1l + fac2*tempx2l + fac3*tempx3l)
- sum_terms(2,i,j,k) = - (fac1*tempy1l + fac2*tempy2l + fac3*tempy3l)
- sum_terms(3,i,j,k) = - (fac1*tempz1l + fac2*tempz2l + fac3*tempz3l)
-
- 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
-
-! update memory variables based upon a 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
- if(ATTENUATION_VAL) then
- do k = 1,NGLLZ
- do j = 1,NGLLY
- do i = 1,NGLLX
- do i_sls = 1,N_SLS
- do i_memory = 1,5
- R_memory(i_memory,i_sls,i,j,k,ispec) = &
- alphaval(i_sls) * &
- R_memory(i_memory,i_sls,i,j,k,ispec) + muvstore(i,j,k,ispec) * &
- factor_common(i_sls,1,1,1,iregion_selected) * &
- (betaval(i_sls) * &
- epsilondev(i_memory,i,j,k,ispec) + gammaval(i_sls) * epsilondev_loc(i_memory,i,j,k))
- enddo
- enddo
- enddo
- 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
-
Deleted: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_forces_outer_core.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_forces_outer_core.f90 2008-09-11 21:43:27 UTC (rev 12866)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/compute_forces_outer_core.f90 2008-09-12 00:57:15 UTC (rev 12867)
@@ -1,224 +0,0 @@
-!=====================================================================
-!
-! S p e c f e m 3 D G l o b e V e r s i o n 4 . 1
-! --------------------------------------------------
-!
-! 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
-! August 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(d_ln_density_dr_table, &
- displfluid,accelfluid,xstore,ystore,zstore, &
- xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
- hprime_xx,hprime_yy,hprime_zz, &
- hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,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 "values_from_mesher.h"
-
-! displacement and acceleration
- real(kind=CUSTOM_REAL), dimension(nglob_outer_core) :: displfluid,accelfluid
-
-! 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(NGLLY,NGLLY) :: hprime_yy,hprimewgll_yy
- real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ) :: hprime_zz,hprimewgll_zz
- 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
-
-! for gravity
- integer int_radius
- double precision radius,theta,phi
- double precision cos_theta,sin_theta,cos_phi,sin_phi
- double precision, dimension(NRAD_GRAVITY) :: d_ln_density_dr_table
- real(kind=CUSTOM_REAL), dimension(nglob_outer_core) :: xstore,ystore,zstore
-
- 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
-
- double precision grad_x_ln_rho,grad_y_ln_rho,grad_z_ln_rho
-! ****************************************************
-! big loop over all spectral elements in the fluid
-! ****************************************************
-
-! set acceleration to zero
- accelfluid(:) = 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)
-
- tempx1l = 0._CUSTOM_REAL
- tempx2l = 0._CUSTOM_REAL
- tempx3l = 0._CUSTOM_REAL
-
- do l=1,NGLLX
- tempx1l = tempx1l + displfluid(ibool(l,j,k,ispec)) * hprime_xx(i,l)
-!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
-
-!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLY
- tempx2l = tempx2l + displfluid(ibool(i,l,k,ispec)) * hprime_yy(j,l)
-!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
-
-!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLZ
- tempx3l = tempx3l + displfluid(ibool(i,j,l,ispec)) * hprime_zz(k,l)
- enddo
-
-! 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*tempx1l + etaxl*tempx2l + gammaxl*tempx3l
- dpotentialdyl = xiyl*tempx1l + etayl*tempx2l + gammayl*tempx3l
- dpotentialdzl = xizl*tempx1l + etazl*tempx2l + gammazl*tempx3l
-
-! add (chi/rho)grad(rho) term in no gravity case
-
-! 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
-
- 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)
-
-! grad(rho)/rho in Cartesian components
- grad_x_ln_rho = sin_theta * cos_phi * d_ln_density_dr_table(int_radius)
- grad_y_ln_rho = sin_theta * sin_phi * d_ln_density_dr_table(int_radius)
- grad_z_ln_rho = cos_theta * d_ln_density_dr_table(int_radius)
-
-! adding (chi/rho)grad(rho)
- dpotentialdxl = dpotentialdxl + displfluid(iglob) * grad_x_ln_rho
- dpotentialdyl = dpotentialdyl + displfluid(iglob) * grad_y_ln_rho
- dpotentialdzl = dpotentialdzl + displfluid(iglob) * grad_z_ln_rho
-
- tempx1(i,j,k) = jacobianl*(xixl*dpotentialdxl + xiyl*dpotentialdyl + xizl*dpotentialdzl)
- tempx2(i,j,k) = jacobianl*(etaxl*dpotentialdxl + etayl*dpotentialdyl + etazl*dpotentialdzl)
- tempx3(i,j,k) = jacobianl*(gammaxl*dpotentialdxl + gammayl*dpotentialdyl + gammazl*dpotentialdzl)
-
- enddo
- enddo
- enddo
-
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
-
- tempx1l = 0._CUSTOM_REAL
- tempx2l = 0._CUSTOM_REAL
- tempx3l = 0._CUSTOM_REAL
-
- do l=1,NGLLX
- tempx1l = tempx1l + tempx1(l,j,k) * hprimewgll_xx(l,i)
-!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
-
-!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLY
- tempx2l = tempx2l + tempx2(i,l,k) * hprimewgll_yy(l,j)
-!!! can merge these loops because NGLLX = NGLLY = NGLLZ enddo
-
-!!! can merge these loops because NGLLX = NGLLY = NGLLZ do l=1,NGLLZ
- tempx3l = tempx3l + tempx3(i,j,l) * hprimewgll_zz(l,k)
- enddo
-
-! sum contributions from each element to the global mesh and add gravity term
-
- iglob = ibool(i,j,k,ispec)
- accelfluid(iglob) = accelfluid(iglob) - (wgllwgll_yz(j,k)*tempx1l + wgllwgll_xz(i,k)*tempx2l + wgllwgll_xy(i,j)*tempx3l)
-
- enddo
- enddo
- enddo
-
- enddo ! spectral element loop
-
- end subroutine compute_forces_outer_core
-
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/create_regions_mesh.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/create_regions_mesh.F90 2008-09-11 21:43:27 UTC (rev 12866)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/create_regions_mesh.F90 2008-09-12 00:57:15 UTC (rev 12867)
@@ -50,7 +50,7 @@
xread1D_leftxi_righteta,xread1D_rightxi_righteta,yread1D_leftxi_lefteta,yread1D_rightxi_lefteta,yread1D_leftxi_righteta, &
yread1D_rightxi_righteta,zread1D_leftxi_lefteta,zread1D_rightxi_lefteta,zread1D_leftxi_righteta,zread1D_rightxi_righteta, &
#endif
- rho_vp,rho_vs,Qmu_store,tau_e_store,ifirst_layer_aniso,ilast_layer_aniso,SAVE_MESH_FILES)
+ rho_vp,rho_vs,Qmu_store,tau_e_store,ifirst_layer_aniso,ilast_layer_aniso,SAVE_MESH_FILES,is_on_a_slice_edge)
! create the different regions of the mesh
@@ -69,6 +69,8 @@
!! DK DK to debug the second sorting routine
logical, parameter :: DEBUG = .true.
+ logical, dimension(nspec) :: is_on_a_slice_edge
+
!! DK DK added this for merged version
#ifdef USE_MPI
integer :: npoin2D_xi,npoin2D_eta
@@ -152,7 +154,7 @@
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 :: Qone_minus_sum_beta, Qone_minus_sum_beta2 ! 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
@@ -837,6 +839,9 @@
iboun(5,ispec)= .true.
endif
+ is_on_a_slice_edge(ispec) = iMPIcut_xi(1,ispec) .or. iMPIcut_xi(2,ispec) .or. &
+ iMPIcut_eta(1,ispec) .or. iMPIcut_eta(2,ispec) .or. iboun(5,ispec) .or. iboun(6,ispec)
+
! define the doubling flag of this element
if(iregion_code /= IREGION_OUTER_CORE) idoubling(ispec) = doubling_index(ilayer)
@@ -1021,6 +1026,9 @@
if (ilayer==ifirst_layer) iboun(6,ispec)= iboun_sb(ispec_superbrick,6)
if (ilayer==ilast_layer .and. iz_elem==1) iboun(5,ispec)= iboun_sb(ispec_superbrick,5)
+ is_on_a_slice_edge(ispec) = iMPIcut_xi(1,ispec) .or. iMPIcut_xi(2,ispec) .or. &
+ iMPIcut_eta(1,ispec) .or. iMPIcut_eta(2,ispec) .or. iboun(5,ispec) .or. iboun(6,ispec)
+
! define the doubling flag of this element
if(iregion_code /= IREGION_OUTER_CORE) idoubling(ispec) = doubling_index(ilayer)
@@ -1162,6 +1170,9 @@
if (iproc_eta == NPROC_ETA-1) iboun(4,ispec)= .true.
endif
+ is_on_a_slice_edge(ispec) = iMPIcut_xi(1,ispec) .or. iMPIcut_xi(2,ispec) .or. &
+ iMPIcut_eta(1,ispec) .or. iMPIcut_eta(2,ispec) .or. iboun(5,ispec) .or. iboun(6,ispec)
+
! define the doubling flag of this element
! only two active central cubes, the four others are fictitious
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/get_model.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/get_model.f90 2008-09-11 21:43:27 UTC (rev 12866)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/get_model.f90 2008-09-12 00:57:15 UTC (rev 12867)
@@ -76,7 +76,7 @@
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 :: Qone_minus_sum_beta, Qone_minus_sum_beta2 ! 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
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/main_program.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/main_program.F90 2008-09-11 21:43:27 UTC (rev 12866)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/main_program.F90 2008-09-12 00:57:15 UTC (rev 12867)
@@ -256,6 +256,11 @@
integer npoin2D_faces_inner_core(NUMFACES_SHARED)
#endif
+! for non blocking communications
+ logical, dimension(NSPEC_CRUST_MANTLE) :: is_on_a_slice_edge_crust_mantle
+ logical, dimension(NSPEC_OUTER_CORE) :: is_on_a_slice_edge_outer_core
+ logical, dimension(NSPEC_INNER_CORE) :: is_on_a_slice_edge_inner_core
+
integer, dimension(NB_SQUARE_EDGES_ONEDIR) :: npoin2D_xi_crust_mantle,npoin2D_eta_crust_mantle, &
npoin2D_xi_outer_core,npoin2D_eta_outer_core, &
npoin2D_xi_inner_core,npoin2D_eta_inner_core
@@ -349,7 +354,7 @@
integer, dimension(:), pointer :: Qs ! 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 :: Qone_minus_sum_beta, Qone_minus_sum_beta2 ! 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
@@ -491,7 +496,8 @@
bcast_integer,bcast_double_precision,bcast_logical,MODEL,ner,ratio_sampling_array,doubling_index, &
r_bottom,r_top,rmins,rmaxs,this_layer_has_a_doubling,NSPEC_computed,NSPEC2D_XI,NSPEC2D_ETA, &
NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP,NSPEC1D_RADIAL,NGLOB1D_RADIAL, &
- NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX,NGLOB_computed,DIFF_NSPEC1D_RADIAL,DIFF_NSPEC2D_ETA,DIFF_NSPEC2D_XI)
+ NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX,NGLOB_computed,DIFF_NSPEC1D_RADIAL,DIFF_NSPEC2D_ETA,DIFF_NSPEC2D_XI, &
+ is_on_a_slice_edge_crust_mantle,is_on_a_slice_edge_outer_core,is_on_a_slice_edge_inner_core)
! synchronize all the processes to make sure everybody has finished creating the mesh
#ifdef USE_MPI
@@ -549,7 +555,8 @@
bcast_integer,bcast_double_precision,bcast_logical,MODEL,ner,ratio_sampling_array,doubling_index, &
r_bottom,r_top,rmins,rmaxs,this_layer_has_a_doubling,NSPEC_computed,NSPEC2D_XI,NSPEC2D_ETA, &
NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP,NSPEC1D_RADIAL,NGLOB1D_RADIAL, &
- NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX,NGLOB_computed,DIFF_NSPEC1D_RADIAL,DIFF_NSPEC2D_ETA,DIFF_NSPEC2D_XI)
+ NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX,NGLOB_computed,DIFF_NSPEC1D_RADIAL,DIFF_NSPEC2D_ETA,DIFF_NSPEC2D_XI, &
+ is_on_a_slice_edge_crust_mantle,is_on_a_slice_edge_outer_core,is_on_a_slice_edge_inner_core)
! synchronize all the processes to make sure everybody has finished
#ifdef USE_MPI
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/meshfem3D.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/meshfem3D.F90 2008-09-11 21:43:27 UTC (rev 12866)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/meshfem3D.F90 2008-09-12 00:57:15 UTC (rev 12867)
@@ -59,7 +59,8 @@
bcast_integer,bcast_double_precision,bcast_logical,MODEL,ner,ratio_sampling_array,doubling_index, &
r_bottom,r_top,rmins,rmaxs,this_layer_has_a_doubling,NSPEC,NSPEC2D_XI,NSPEC2D_ETA, &
NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP,NSPEC1D_RADIAL,NGLOB1D_RADIAL, &
- NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX,NGLOB,DIFF_NSPEC1D_RADIAL,DIFF_NSPEC2D_ETA,DIFF_NSPEC2D_XI)
+ NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX,NGLOB,DIFF_NSPEC1D_RADIAL,DIFF_NSPEC2D_ETA,DIFF_NSPEC2D_XI, &
+ is_on_a_slice_edge_crust_mantle,is_on_a_slice_edge_outer_core,is_on_a_slice_edge_inner_core)
use dyn_array
@@ -77,6 +78,11 @@
! include values created by the mesher
include "values_from_mesher.h"
+! for non blocking communications
+ logical, dimension(NSPEC_CRUST_MANTLE) :: is_on_a_slice_edge_crust_mantle
+ logical, dimension(NSPEC_OUTER_CORE) :: is_on_a_slice_edge_outer_core
+ logical, dimension(NSPEC_INNER_CORE) :: is_on_a_slice_edge_inner_core
+
! aniso_mantle_model_variables
type aniso_mantle_model_variables
sequence
@@ -99,7 +105,7 @@
integer, dimension(:), pointer :: Qs ! 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 :: Qone_minus_sum_beta, Qone_minus_sum_beta2 ! 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
@@ -1535,7 +1541,7 @@
xread1D_leftxi_righteta,xread1D_rightxi_righteta,yread1D_leftxi_lefteta,yread1D_rightxi_lefteta,yread1D_leftxi_righteta, &
yread1D_rightxi_righteta,zread1D_leftxi_lefteta,zread1D_rightxi_lefteta,zread1D_leftxi_righteta,zread1D_rightxi_righteta, &
#endif
- rho_vp,rho_vs,Qmu_store,tau_e_store,ifirst_layer_aniso,ilast_layer_aniso,SAVE_MESH_FILES)
+ rho_vp,rho_vs,Qmu_store,tau_e_store,ifirst_layer_aniso,ilast_layer_aniso,SAVE_MESH_FILES,is_on_a_slice_edge_crust_mantle)
else if(iregion_code == IREGION_OUTER_CORE) then
! outer_core
@@ -1570,7 +1576,7 @@
xread1D_leftxi_righteta,xread1D_rightxi_righteta,yread1D_leftxi_lefteta,yread1D_rightxi_lefteta,yread1D_leftxi_righteta, &
yread1D_rightxi_righteta,zread1D_leftxi_lefteta,zread1D_rightxi_lefteta,zread1D_leftxi_righteta,zread1D_rightxi_righteta, &
#endif
- rho_vp,rho_vs,Qmu_store,tau_e_store,ifirst_layer_aniso,ilast_layer_aniso,SAVE_MESH_FILES)
+ rho_vp,rho_vs,Qmu_store,tau_e_store,ifirst_layer_aniso,ilast_layer_aniso,SAVE_MESH_FILES,is_on_a_slice_edge_outer_core)
else if(iregion_code == IREGION_INNER_CORE) then
! inner_core
@@ -1605,7 +1611,7 @@
xread1D_leftxi_righteta,xread1D_rightxi_righteta,yread1D_leftxi_lefteta,yread1D_rightxi_lefteta,yread1D_leftxi_righteta, &
yread1D_rightxi_righteta,zread1D_leftxi_lefteta,zread1D_rightxi_lefteta,zread1D_leftxi_righteta,zread1D_rightxi_righteta, &
#endif
- rho_vp,rho_vs,Qmu_store,tau_e_store,ifirst_layer_aniso,ilast_layer_aniso,SAVE_MESH_FILES)
+ rho_vp,rho_vs,Qmu_store,tau_e_store,ifirst_layer_aniso,ilast_layer_aniso,SAVE_MESH_FILES,is_on_a_slice_edge_inner_core)
else
stop 'DK DK incorrect region in merged code'
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/specfem3D.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/specfem3D.F90 2008-09-11 21:43:27 UTC (rev 12866)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/specfem3D.F90 2008-09-12 00:57:15 UTC (rev 12867)
@@ -59,7 +59,8 @@
bcast_integer,bcast_double_precision,bcast_logical,MODEL,ner,ratio_sampling_array,doubling_index, &
r_bottom,r_top,rmins,rmaxs,this_layer_has_a_doubling,NSPEC_computed,NSPEC2D_XI,NSPEC2D_ETA, &
NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP,NSPEC1D_RADIAL,NGLOB1D_RADIAL, &
- NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX,NGLOB_computed,DIFF_NSPEC1D_RADIAL,DIFF_NSPEC2D_ETA,DIFF_NSPEC2D_XI)
+ NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX,NGLOB_computed,DIFF_NSPEC1D_RADIAL,DIFF_NSPEC2D_ETA,DIFF_NSPEC2D_XI, &
+ is_on_a_slice_edge_crust_mantle,is_on_a_slice_edge_outer_core,is_on_a_slice_edge_inner_core)
use dyn_array
@@ -76,6 +77,13 @@
! include values created by the mesher
include "values_from_mesher.h"
+! for non blocking communications
+ integer :: icall
+ real :: percentage_edge
+ logical, dimension(NSPEC_CRUST_MANTLE) :: is_on_a_slice_edge_crust_mantle
+ logical, dimension(NSPEC_OUTER_CORE) :: is_on_a_slice_edge_outer_core
+ logical, dimension(NSPEC_INNER_CORE) :: is_on_a_slice_edge_inner_core
+
! attenuation_model_variables
type attenuation_model_variables
sequence
@@ -87,7 +95,7 @@
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 :: Qone_minus_sum_beta, Qone_minus_sum_beta2 ! 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
@@ -100,8 +108,8 @@
! memory variables and standard linear solids for attenuation
double precision, dimension(N_SLS) :: tau_sigma_dble
- double precision, dimension(ATT1,ATT2,ATT3,ATT4) :: omsb_crust_mantle_dble, factor_scale_crust_mantle_dble
- double precision, dimension(ATT1,ATT2,ATT3,ATT5) :: omsb_inner_core_dble, factor_scale_inner_core_dble
+ double precision, dimension(ATT1,ATT2,ATT3,ATT4) :: one_minus_sum_beta_CM_dble, factor_scale_crust_mantle_dble
+ double precision, dimension(ATT1,ATT2,ATT3,ATT5) :: one_minus_sum_beta_IC_dble, factor_scale_inner_core_dble
real(kind=CUSTOM_REAL), dimension(ATT1,ATT2,ATT3,ATT4) :: one_minus_sum_beta_crust_mantle, factor_scale_crust_mantle
real(kind=CUSTOM_REAL), dimension(ATT1,ATT2,ATT3,ATT5) :: one_minus_sum_beta_inner_core, factor_scale_inner_core
@@ -809,6 +817,24 @@
write(IMAIN,*) 'smallest and largest possible floating-point numbers are: ',tiny(1._CUSTOM_REAL),huge(1._CUSTOM_REAL)
write(IMAIN,*)
+ write(IMAIN,*) 'for overlapping of communications with calculations:'
+ write(IMAIN,*)
+
+ percentage_edge = 100.*count(is_on_a_slice_edge_crust_mantle(:))/real(NSPEC_CRUST_MANTLE)
+ write(IMAIN,*) 'percentage of edge elements in crust/mantle ',percentage_edge,'%'
+ write(IMAIN,*) 'percentage of volume elements in crust/mantle ',100. - percentage_edge,'%'
+ write(IMAIN,*)
+
+ percentage_edge = 100.*count(is_on_a_slice_edge_outer_core(:))/real(NSPEC_OUTER_CORE)
+ write(IMAIN,*) 'percentage of edge elements in outer core ',percentage_edge,'%'
+ write(IMAIN,*) 'percentage of volume elements in outer core ',100. - percentage_edge,'%'
+ write(IMAIN,*)
+
+ percentage_edge = 100.*count(is_on_a_slice_edge_inner_core(:))/real(NSPEC_INNER_CORE)
+ write(IMAIN,*) 'percentage of edge elements in inner core ',percentage_edge,'%'
+ write(IMAIN,*) 'percentage of volume elements in inner core ',100. - percentage_edge,'%'
+ write(IMAIN,*)
+
endif
! check that the code is running with the requested nb of processes
@@ -1181,7 +1207,7 @@
! ocean load
if (OCEANS) then
- call assemble_MPI_scalar(myrank,rmass_ocean_load,NGLOB_CRUST_MANTLE, &
+ call assemble_MPI_scalar_block(myrank,rmass_ocean_load,NGLOB_CRUST_MANTLE, &
iproc_xi,iproc_eta,ichunk,addressing, &
iboolleft_xi_crust_mantle,iboolright_xi_crust_mantle,iboolleft_eta_crust_mantle,iboolright_eta_crust_mantle, &
npoin2D_faces_crust_mantle,npoin2D_xi_crust_mantle,npoin2D_eta_crust_mantle, &
@@ -1196,7 +1222,7 @@
endif
! crust and mantle
- call assemble_MPI_scalar(myrank,rmass_crust_mantle,NGLOB_CRUST_MANTLE, &
+ call assemble_MPI_scalar_block(myrank,rmass_crust_mantle,NGLOB_CRUST_MANTLE, &
iproc_xi,iproc_eta,ichunk,addressing, &
iboolleft_xi_crust_mantle,iboolright_xi_crust_mantle,iboolleft_eta_crust_mantle,iboolright_eta_crust_mantle, &
npoin2D_faces_crust_mantle,npoin2D_xi_crust_mantle,npoin2D_eta_crust_mantle, &
@@ -1210,7 +1236,7 @@
NGLOB2DMAX_XMIN_XMAX(IREGION_CRUST_MANTLE),NGLOB2DMAX_YMIN_YMAX(IREGION_CRUST_MANTLE),NGLOB2DMAX_XY_VAL_CM,NCHUNKS)
! outer core
- call assemble_MPI_scalar(myrank,rmass_outer_core,NGLOB_OUTER_CORE, &
+ call assemble_MPI_scalar_block(myrank,rmass_outer_core,NGLOB_OUTER_CORE, &
iproc_xi,iproc_eta,ichunk,addressing, &
iboolleft_xi_outer_core,iboolright_xi_outer_core,iboolleft_eta_outer_core,iboolright_eta_outer_core, &
npoin2D_faces_outer_core,npoin2D_xi_outer_core,npoin2D_eta_outer_core, &
@@ -1224,7 +1250,7 @@
NGLOB2DMAX_XMIN_XMAX(IREGION_OUTER_CORE),NGLOB2DMAX_YMIN_YMAX(IREGION_OUTER_CORE),NGLOB2DMAX_XY_VAL_OC,NCHUNKS)
! inner core
- call assemble_MPI_scalar(myrank,rmass_inner_core,NGLOB_INNER_CORE, &
+ call assemble_MPI_scalar_block(myrank,rmass_inner_core,NGLOB_INNER_CORE, &
iproc_xi,iproc_eta,ichunk,addressing, &
iboolleft_xi_inner_core,iboolright_xi_inner_core,iboolleft_eta_inner_core,iboolright_eta_inner_core, &
npoin2D_faces_inner_core,npoin2D_xi_inner_core,npoin2D_eta_inner_core, &
@@ -1308,7 +1334,7 @@
ndim_assemble = 1
! use these buffers to assemble the inner core mass matrix with the central cube
- call assemble_MPI_central_cube(ichunk,nb_msgs_theor_in_cube, sender_from_slices_to_cube, &
+ call assemble_MPI_central_cube_block(ichunk,nb_msgs_theor_in_cube, sender_from_slices_to_cube, &
npoin2D_cube_from_slices, buffer_all_cube_from_slices, buffer_slices, buffer_slices2, ibool_central_cube, &
receiver_cube_from_slices, ibool_inner_core, idoubling_inner_core, NSPEC_INNER_CORE, &
ibelm_bottom_inner_core, NSPEC2D_BOTTOM(IREGION_INNER_CORE),NGLOB_INNER_CORE,rmass_inner_core,ndim_assemble)
@@ -1370,9 +1396,9 @@
! get and store PREM attenuation model
call get_attenuation_model_1D(myrank, IREGION_CRUST_MANTLE, tau_sigma_dble, &
- omsb_crust_mantle_dble, factor_common_crust_mantle_dble, &
+ one_minus_sum_beta_CM_dble, factor_common_crust_mantle_dble, &
factor_scale_crust_mantle_dble, NRAD_ATTENUATION,1,1,1, AM_V)
- omsb_inner_core_dble(:,:,:,1:min(ATT4,ATT5)) = omsb_crust_mantle_dble(:,:,:,1:min(ATT4,ATT5))
+ one_minus_sum_beta_IC_dble(:,:,:,1:min(ATT4,ATT5)) = one_minus_sum_beta_CM_dble(:,:,:,1:min(ATT4,ATT5))
factor_scale_inner_core_dble(:,:,:,1:min(ATT4,ATT5)) = factor_scale_crust_mantle_dble(:,:,:,1:min(ATT4,ATT5))
factor_common_inner_core_dble(:,:,:,:,1:min(ATT4,ATT5)) = factor_common_crust_mantle_dble(:,:,:,:,1:min(ATT4,ATT5))
! Tell the Attenuation Code about the IDOUBLING regions within the Mesh
@@ -1380,19 +1406,19 @@
if(CUSTOM_REAL == SIZE_REAL) then
factor_scale_crust_mantle = sngl(factor_scale_crust_mantle_dble)
- one_minus_sum_beta_crust_mantle = sngl(omsb_crust_mantle_dble)
+ one_minus_sum_beta_crust_mantle = sngl(one_minus_sum_beta_CM_dble)
factor_common_crust_mantle = sngl(factor_common_crust_mantle_dble)
factor_scale_inner_core = sngl(factor_scale_inner_core_dble)
- one_minus_sum_beta_inner_core = sngl(omsb_inner_core_dble)
+ one_minus_sum_beta_inner_core = sngl(one_minus_sum_beta_IC_dble)
factor_common_inner_core = sngl(factor_common_inner_core_dble)
else
factor_scale_crust_mantle = factor_scale_crust_mantle_dble
- one_minus_sum_beta_crust_mantle = omsb_crust_mantle_dble
+ one_minus_sum_beta_crust_mantle = one_minus_sum_beta_CM_dble
factor_common_crust_mantle = factor_common_crust_mantle_dble
factor_scale_inner_core = factor_scale_inner_core_dble
- one_minus_sum_beta_inner_core = omsb_inner_core_dble
+ one_minus_sum_beta_inner_core = one_minus_sum_beta_IC_dble
factor_common_inner_core = factor_common_inner_core_dble
endif
@@ -1947,15 +1973,30 @@
time = (dble(it-1)*DT-t0)/scale_t
endif
-! accel_outer_core, div_displ_outer_core are initialized to zero in the following subroutine.
- call compute_forces_outer_core(d_ln_density_dr_table, &
+ icall = 1
+ call compute_forces_OC(d_ln_density_dr_table, &
displ_outer_core,accel_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,ibool_outer_core)
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,ibool_outer_core,icall,is_on_a_slice_edge_outer_core)
+#ifndef USE_MPI
+!! DK DK put a fictitious source in each region in the case of a serial test if needed
+ if(PUT_SOURCE_IN_EACH_REGION) then
+ stf = 1.d-6 * comp_source_time_function(dble(it-1)*DT-t0,10.d0)
+! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ stf_used = sngl(stf)
+ else
+ stf_used = stf
+ endif
+ iglob = ibool_outer_core(2,2,2,2)
+ accel_outer_core(iglob) = accel_outer_core(iglob) + stf_used
+ endif
+#endif
+
! ****************************************************
! ********** add matching with solid part **********
! ****************************************************
@@ -2058,6 +2099,15 @@
endif
+ icall = 2
+ call compute_forces_OC(d_ln_density_dr_table, &
+ displ_outer_core,accel_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,ibool_outer_core,icall,is_on_a_slice_edge_outer_core)
+
! assemble all the contributions between slices using MPI
! outer core
@@ -2076,21 +2126,6 @@
NGLOB2DMAX_XMIN_XMAX(IREGION_OUTER_CORE),NGLOB2DMAX_YMIN_YMAX(IREGION_OUTER_CORE),NGLOB2DMAX_XY_VAL_OC,NCHUNKS)
#endif
-#ifndef USE_MPI
-!! DK DK put a fictitious source in each region in the case of a serial test if needed
- if(PUT_SOURCE_IN_EACH_REGION) then
- stf = 1.d-6 * comp_source_time_function(dble(it-1)*DT-t0,10.d0)
-! distinguish between single and double precision for reals
- if(CUSTOM_REAL == SIZE_REAL) then
- stf_used = sngl(stf)
- else
- stf_used = stf
- endif
- iglob = ibool_outer_core(2,2,2,2)
- accel_outer_core(iglob) = accel_outer_core(iglob) + stf_used
- endif
-#endif
-
! multiply by the inverse of the mass matrix and update velocity
do i=1,NGLOB_OUTER_CORE
accel_outer_core(i) = accel_outer_core(i)*rmass_outer_core(i)
@@ -2105,34 +2140,33 @@
! for anisotropy and gravity, x y and z contain r theta and phi
- call compute_forces_crust_mantle(displ_crust_mantle,accel_crust_mantle, &
+ icall = 1
+ call compute_forces_CM_IC(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, &
+ 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, &
+ kappavstore_crust_mantle,kappahstore_crust_mantle,muvstore_crust_mantle,muhstore_crust_mantle, &
+ eta_anisostore_crust_mantle, &
+ ibool_crust_mantle,idoubling_crust_mantle,R_memory_crust_mantle,epsilondev_crust_mantle,one_minus_sum_beta_crust_mantle, &
+ 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), &
+ is_on_a_slice_edge_crust_mantle, &
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ displ_inner_core,accel_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, &
+ kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
+ R_memory_inner_core,epsilondev_inner_core,&
+ one_minus_sum_beta_inner_core,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), &
+ is_on_a_slice_edge_inner_core, &
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
hprime_xx,hprime_yy,hprime_zz, &
hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
- kappavstore_crust_mantle,kappahstore_crust_mantle,muvstore_crust_mantle, &
- muhstore_crust_mantle,eta_anisostore_crust_mantle, &
- ibool_crust_mantle,idoubling_crust_mantle, &
- R_memory_crust_mantle,epsilondev_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)
-
- call compute_forces_inner_core(displ_inner_core,accel_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, &
- kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
- R_memory_inner_core,epsilondev_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)
+ COMPUTE_AND_STORE_STRAIN,AM_V,icall)
#ifdef USE_MPI
@@ -2178,6 +2212,21 @@
accel_crust_mantle(3,iglob) = accel_crust_mantle(3,iglob) + stf_used
#endif
+#ifndef USE_MPI
+!! DK DK put a fictitious source in each region in the case of a serial test if needed
+ if(PUT_SOURCE_IN_EACH_REGION) then
+ stf = 1.d-6 * comp_source_time_function(dble(it-1)*DT-t0,10.d0)
+! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ stf_used = sngl(stf)
+ else
+ stf_used = stf
+ endif
+ iglob = ibool_inner_core(2,2,2,2)
+ accel_inner_core(3,iglob) = accel_inner_core(3,iglob) + stf_used
+ endif
+#endif
+
! ****************************************************
! ********** add matching with fluid part **********
! ****************************************************
@@ -2276,6 +2325,34 @@
endif
+ icall = 2
+ call compute_forces_CM_IC(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, &
+ kappavstore_crust_mantle,kappahstore_crust_mantle,muvstore_crust_mantle,muhstore_crust_mantle, &
+ eta_anisostore_crust_mantle, &
+ ibool_crust_mantle,idoubling_crust_mantle,R_memory_crust_mantle,epsilondev_crust_mantle,one_minus_sum_beta_crust_mantle, &
+ 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), &
+ is_on_a_slice_edge_crust_mantle, &
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ displ_inner_core,accel_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, &
+ kappavstore_inner_core,muvstore_inner_core,ibool_inner_core,idoubling_inner_core, &
+ R_memory_inner_core,epsilondev_inner_core,&
+ one_minus_sum_beta_inner_core,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), &
+ is_on_a_slice_edge_inner_core, &
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ hprime_xx,hprime_yy,hprime_zz, &
+ hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
+ alphaval,betaval,gammaval, &
+ COMPUTE_AND_STORE_STRAIN,AM_V,icall)
+
! assemble all the contributions between slices using MPI
! crust/mantle and inner core handled in the same call
@@ -2376,21 +2453,6 @@
veloc_crust_mantle(:,i) = veloc_crust_mantle(:,i) + deltatover2*accel_crust_mantle(:,i)
enddo
-#ifndef USE_MPI
-!! DK DK put a fictitious source in each region in the case of a serial test if needed
- if(PUT_SOURCE_IN_EACH_REGION) then
- stf = 1.d-6 * comp_source_time_function(dble(it-1)*DT-t0,10.d0)
-! distinguish between single and double precision for reals
- if(CUSTOM_REAL == SIZE_REAL) then
- stf_used = sngl(stf)
- else
- stf_used = stf
- endif
- iglob = ibool_inner_core(2,2,2,2)
- accel_inner_core(3,iglob) = accel_inner_core(3,iglob) + stf_used
- endif
-#endif
-
do i=1,NGLOB_INNER_CORE
accel_inner_core(1,i) = accel_inner_core(1,i)*rmass_inner_core(i)
accel_inner_core(2,i) = accel_inner_core(2,i)*rmass_inner_core(i)
More information about the cig-commits
mailing list