[cig-commits] r22476 - in seismo/3D/SPECFEM3D_GLOBE: branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D trunk/src/meshfem3D
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Sun Jun 30 21:44:25 PDT 2013
Author: dkomati1
Date: 2013-06-30 21:44:25 -0700 (Sun, 30 Jun 2013)
New Revision: 22476
Modified:
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_mass_matrices.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_absorb.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_model.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_doubling_elements.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_mass_matrices.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_regular_elements.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_MPI_cutplanes_eta.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_MPI_cutplanes_xi.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_absorb.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_jacobian_boundaries.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_model.f90
Log:
more obvious merges
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_mass_matrices.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_mass_matrices.f90 2013-07-01 04:21:23 UTC (rev 22475)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_mass_matrices.f90 2013-07-01 04:44:25 UTC (rev 22476)
@@ -29,7 +29,7 @@
iregion_code,xstore,ystore,zstore, &
NSPEC2D_TOP,NSPEC2D_BOTTOM)
-! creates rmassx, rmassy, rmassz and rmass_ocean_load
+ ! creates rmassx, rmassy, rmassz and rmass_ocean_load
use constants
@@ -74,8 +74,13 @@
! initializes matrices
!
+ ! in the case of stacey boundary conditions, add C*delta/2 contribution to the mass matrix
+ ! on the Stacey edges for the crust_mantle and outer_core regions but not for the inner_core region
+ ! thus the mass matrix must be replaced by three mass matrices including the "C" damping matrix
+ !
! if absorbing_conditions are not set or if NCHUNKS=6, only one mass matrix is needed
! for the sake of performance, only "rmassz" array will be filled and "rmassx" & "rmassy" will be obsolete
+
rmassx(:) = 0._CUSTOM_REAL
rmassy(:) = 0._CUSTOM_REAL
rmassz(:) = 0._CUSTOM_REAL
@@ -109,6 +114,7 @@
! definition depends if region is fluid or solid
select case( iregion_code)
+
case( IREGION_CRUST_MANTLE, IREGION_INNER_CORE )
! distinguish between single and double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
@@ -234,7 +240,7 @@
endif
- ! adds C*deltat/2 contribution to the mass matrices on Stacey edges
+ ! add C*deltat/2 contribution to the mass matrices on the Stacey edges
if(NCHUNKS /= 6 .and. ABSORBING_CONDITIONS) then
call create_mass_matrices_Stacey(myrank,nspec,ibool,iregion_code, &
NSPEC2D_BOTTOM)
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_absorb.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_absorb.f90 2013-07-01 04:21:23 UTC (rev 22475)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_absorb.f90 2013-07-01 04:44:25 UTC (rev 22476)
@@ -25,7 +25,7 @@
!
!=====================================================================
- subroutine get_absorb(myrank,prname,iregion, iboun,nspec, &
+ subroutine get_absorb(myrank,prname,iregion,iboun,nspec, &
nimin,nimax,njmin,njmax,nkmin_xi,nkmin_eta, &
NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM)
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_model.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_model.f90 2013-07-01 04:21:23 UTC (rev 22475)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/get_model.f90 2013-07-01 04:44:25 UTC (rev 22476)
@@ -37,7 +37,6 @@
tau_s,tau_e_store,Qmu_store,T_c_source,vx,vy,vz,vnspec, &
elem_in_crust,elem_in_mantle)
-
use meshfem3D_par,only: &
RCMB,RICB,R670,RMOHO,RTOPDDOUBLEPRIME,R600,R220, &
R771,R400,R120,R80,RMIDDLE_CRUST,ROCEAN, &
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_doubling_elements.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_doubling_elements.f90 2013-07-01 04:21:23 UTC (rev 22475)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_doubling_elements.f90 2013-07-01 04:44:25 UTC (rev 22476)
@@ -346,13 +346,14 @@
c23store,c24store,c25store,c26store,c33store,c34store,c35store, &
c36store,c44store,c45store,c46store,c55store,c56store,c66store, &
nspec_ani,nspec_stacey,nspec_att,Qmu_store,tau_e_store,tau_s,T_c_source,&
- rho_vp,rho_vs,ACTUALLY_STORE_ARRAYS,&
+ rho_vp,rho_vs,ACTUALLY_STORE_ARRAYS, &
xigll,yigll,zigll,ispec_is_tiso)
! boundary mesh
if (ipass == 2 .and. SAVE_BOUNDARY_MESH .and. iregion_code == IREGION_CRUST_MANTLE) then
is_superbrick=.true.
- call get_jacobian_discontinuities(myrank,ispec,ix_elem,iy_elem,rmin,rmax,r1,r2,r3,r4,r5,r6,r7,r8, &
+ call get_jacobian_discontinuities(myrank,ispec,ix_elem,iy_elem,rmin,rmax, &
+ r1,r2,r3,r4,r5,r6,r7,r8, &
xstore(:,:,:,ispec),ystore(:,:,:,ispec),zstore(:,:,:,ispec),dershape2D_bottom, &
ibelm_moho_top,ibelm_moho_bot,ibelm_400_top,ibelm_400_bot,ibelm_670_top,ibelm_670_bot, &
normal_moho,normal_400,normal_670,jacobian2D_moho,jacobian2D_400,jacobian2D_670, &
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_mass_matrices.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_mass_matrices.f90 2013-07-01 04:21:23 UTC (rev 22475)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_mass_matrices.f90 2013-07-01 04:44:25 UTC (rev 22476)
@@ -333,7 +333,7 @@
endif
- ! add C*delta/2 contribution to the mass matrices on the Stacey edges
+ ! add C*deltat/2 contribution to the mass matrices on the Stacey edges
if(NCHUNKS /= 6 .and. ABSORBING_CONDITIONS) then
! read arrays for Stacey conditions
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_regular_elements.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_regular_elements.f90 2013-07-01 04:21:23 UTC (rev 22475)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_regular_elements.f90 2013-07-01 04:44:25 UTC (rev 22476)
@@ -272,7 +272,8 @@
if (ipass == 2 .and. SAVE_BOUNDARY_MESH .and. iregion_code == IREGION_CRUST_MANTLE) then
is_superbrick=.false.
ispec_superbrick=0
- call get_jacobian_discontinuities(myrank,ispec,ix_elem,iy_elem,rmin,rmax,r1,r2,r3,r4,r5,r6,r7,r8, &
+ call get_jacobian_discontinuities(myrank,ispec,ix_elem,iy_elem,rmin,rmax, &
+ r1,r2,r3,r4,r5,r6,r7,r8, &
xstore(:,:,:,ispec),ystore(:,:,:,ispec),zstore(:,:,:,ispec),dershape2D_bottom, &
ibelm_moho_top,ibelm_moho_bot,ibelm_400_top,ibelm_400_bot,ibelm_670_top,ibelm_670_bot, &
normal_moho,normal_400,normal_670,jacobian2D_moho,jacobian2D_400,jacobian2D_670, &
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_MPI_cutplanes_eta.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_MPI_cutplanes_eta.f90 2013-07-01 04:21:23 UTC (rev 22475)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_MPI_cutplanes_eta.f90 2013-07-01 04:44:25 UTC (rev 22476)
@@ -26,8 +26,8 @@
!=====================================================================
subroutine get_MPI_cutplanes_eta(myrank,prname,nspec,iMPIcut_eta,ibool, &
- xstore,ystore,zstore,mask_ibool,npointot, &
- NSPEC2D_XI_FACE,iregion,npoin2D_eta)
+ xstore,ystore,zstore,mask_ibool,npointot, &
+ NSPEC2D_XI_FACE,iregion,npoin2D_eta)
! this routine detects cut planes along eta
! In principle the left cut plane of the first slice
@@ -38,7 +38,7 @@
include "constants.h"
- integer nspec,myrank,iregion
+ integer :: nspec,myrank,iregion
integer, dimension(MAX_NUM_REGIONS,NB_SQUARE_EDGES_ONEDIR) :: NSPEC2D_XI_FACE
@@ -64,9 +64,9 @@
! processor identification
character(len=150) prname
-! theoretical number of surface elements in the buffers
-! cut planes along eta=constant correspond to XI faces
- nspec2Dtheor = NSPEC2D_XI_FACE(iregion,1)
+ ! theoretical number of surface elements in the buffers
+ ! cut planes along eta=constant correspond to XI faces
+ nspec2Dtheor = NSPEC2D_XI_FACE(iregion,1)
! write the MPI buffers for the left and right edges of the slice
! and the position of the points to check that the buffers are fine
@@ -78,10 +78,10 @@
! global point number and coordinates left MPI cut-plane
open(unit=10,file=prname(1:len_trim(prname))//'iboolleft_eta.txt',status='unknown')
-! erase the logical mask used to mark points already found
+ ! erase the logical mask used to mark points already found
mask_ibool(:) = .false.
-! nb of global points shared with the other slice
+ ! nb of global points shared with the other slice
npoin2D_eta = 0
! nb of elements in this cut-plane
@@ -120,15 +120,15 @@
!
! determine if the element falls on the right MPI cut plane
!
- nspec2Dtheor = NSPEC2D_XI_FACE(iregion,2)
+ nspec2Dtheor = NSPEC2D_XI_FACE(iregion,2)
! global point number and coordinates right MPI cut-plane
open(unit=10,file=prname(1:len_trim(prname))//'iboolright_eta.txt',status='unknown')
-! erase the logical mask used to mark points already found
+ ! erase the logical mask used to mark points already found
mask_ibool(:) = .false.
-! nb of global points shared with the other slice
+ ! nb of global points shared with the other slice
npoin2D_eta = 0
! nb of elements in this cut-plane
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_MPI_cutplanes_xi.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_MPI_cutplanes_xi.f90 2013-07-01 04:21:23 UTC (rev 22475)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_MPI_cutplanes_xi.f90 2013-07-01 04:44:25 UTC (rev 22476)
@@ -26,8 +26,8 @@
!=====================================================================
subroutine get_MPI_cutplanes_xi(myrank,prname,nspec,iMPIcut_xi,ibool, &
- xstore,ystore,zstore,mask_ibool,npointot, &
- NSPEC2D_ETA_FACE,iregion,npoin2D_xi)
+ xstore,ystore,zstore,mask_ibool,npointot, &
+ NSPEC2D_ETA_FACE,iregion,npoin2D_xi)
! this routine detects cut planes along xi
! In principle the left cut plane of the first slice
@@ -38,7 +38,7 @@
include "constants.h"
- integer nspec,myrank,iregion
+ integer :: nspec,myrank,iregion
integer, dimension(MAX_NUM_REGIONS,NB_SQUARE_EDGES_ONEDIR) :: NSPEC2D_ETA_FACE
logical iMPIcut_xi(2,nspec)
@@ -54,7 +54,7 @@
logical mask_ibool(npointot)
! global element numbering
- integer ispec
+ integer :: ispec
! MPI cut-plane element numbering
integer ispecc1,ispecc2,npoin2D_xi,ix,iy,iz
@@ -64,9 +64,10 @@
! processor identification
character(len=150) prname,errmsg
-! theoretical number of surface elements in the buffers
-! cut planes along xi=constant correspond to ETA faces
- nspec2Dtheor = NSPEC2D_ETA_FACE(iregion,1)
+ ! theoretical number of surface elements in the buffers
+ ! cut planes along xi=constant correspond to ETA faces
+ nspec2Dtheor = NSPEC2D_ETA_FACE(iregion,1)
+
! write the MPI buffers for the left and right edges of the slice
! and the position of the points to check that the buffers are fine
@@ -88,15 +89,15 @@
endif
call exit_mpi(myrank,'error creating iboolleft_xi.txt, please check your Par_file LOCAL_PATH setting')
endif
-! erase the logical mask used to mark points already found
+
+ ! erase the logical mask used to mark points already found
mask_ibool(:) = .false.
-! nb of global points shared with the other slice
+ ! nb of global points shared with the other slice
npoin2D_xi = 0
-! nb of elements in this cut-plane
+ ! nb of elements in this cut-plane
ispecc1=0
-
do ispec=1,nspec
if(iMPIcut_xi(1,ispec)) then
ispecc1=ispecc1+1
@@ -124,7 +125,7 @@
close(10)
-! compare number of surface elements detected to analytical value
+ ! compare number of surface elements detected to analytical value
if(ispecc1 /= nspec2Dtheor) then
write(errmsg,*) 'error MPI cut-planes detection in xi=left T=',nspec2Dtheor,' C=',ispecc1
call exit_MPI(myrank,errmsg)
@@ -132,22 +133,21 @@
!
! determine if the element falls on the right MPI cut plane
!
- nspec2Dtheor = NSPEC2D_ETA_FACE(iregion,2)
+ nspec2Dtheor = NSPEC2D_ETA_FACE(iregion,2)
! global point number and coordinates right MPI cut-plane
open(unit=10,file=prname(1:len_trim(prname))//'iboolright_xi.txt', &
status='unknown',iostat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error creating iboolright_xi.txt for this process')
-! erase the logical mask used to mark points already found
+ ! erase the logical mask used to mark points already found
mask_ibool(:) = .false.
-! nb of global points shared with the other slice
+ ! nb of global points shared with the other slice
npoin2D_xi = 0
-! nb of elements in this cut-plane
+ ! nb of elements in this cut-plane
ispecc2=0
-
do ispec=1,nspec
if(iMPIcut_xi(2,ispec)) then
ispecc2=ispecc2+1
@@ -175,7 +175,7 @@
close(10)
-! compare number of surface elements detected to analytical value
+ ! compare number of surface elements detected to analytical value
if(ispecc2 /= nspec2Dtheor) then
write(errmsg,*) 'error MPI cut-planes detection in xi=right T=',nspec2Dtheor,' C=',ispecc2
call exit_MPI(myrank,errmsg)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_absorb.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_absorb.f90 2013-07-01 04:21:23 UTC (rev 22475)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_absorb.f90 2013-07-01 04:44:25 UTC (rev 22476)
@@ -26,8 +26,8 @@
!=====================================================================
subroutine get_absorb(myrank,prname,iboun,nspec, &
- nimin,nimax,njmin,njmax,nkmin_xi,nkmin_eta, &
- NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM)
+ nimin,nimax,njmin,njmax,nkmin_xi,nkmin_eta, &
+ NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM)
! Stacey, define flags for absorbing boundaries
@@ -35,7 +35,7 @@
include "constants.h"
- integer nspec,myrank
+ integer :: nspec,myrank
integer :: NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM
@@ -43,17 +43,18 @@
integer njmin(2,NSPEC2DMAX_XMIN_XMAX),njmax(2,NSPEC2DMAX_XMIN_XMAX)
integer nkmin_xi(2,NSPEC2DMAX_XMIN_XMAX),nkmin_eta(2,NSPEC2DMAX_YMIN_YMAX)
- logical iboun(6,nspec)
+ logical :: iboun(6,nspec)
! global element numbering
integer ispecg
-! counters to keep track of the number of elements on each of the
-! five absorbing boundaries
- integer ispecb1,ispecb2,ispecb3,ispecb4,ispecb5
+ ! counters to keep track of the number of elements on each of the
+ ! five absorbing boundaries
+ integer :: ispecb1,ispecb2,ispecb3,ispecb4,ispecb5
+ integer :: ier
-! processor identification
- character(len=150) prname
+ ! processor identification
+ character(len=150) :: prname
! initializes
ispecb1=0
@@ -64,16 +65,16 @@
do ispecg=1,nspec
-! determine if the element falls on an absorbing boundary
+ ! determine if the element falls on an absorbing boundary
if(iboun(1,ispecg)) then
-! on boundary 1: xmin
- ispecb1=ispecb1+1
+ ! on boundary 1: xmin
+ ispecb1=ispecb1+1
-! this is useful even if it is constant because it can be zero inside the slices
- njmin(1,ispecb1)=1
- njmax(1,ispecb1)=NGLLY
+ ! this is useful even if it is constant because it can be zero inside the slices
+ njmin(1,ispecb1)=1
+ njmax(1,ispecb1)=NGLLY
! check for ovelap with other boundaries
nkmin_xi(1,ispecb1)=1
@@ -82,12 +83,12 @@
if(iboun(2,ispecg)) then
-! on boundary 2: xmax
- ispecb2=ispecb2+1
+ ! on boundary 2: xmax
+ ispecb2=ispecb2+1
-! this is useful even if it is constant because it can be zero inside the slices
- njmin(2,ispecb2)=1
- njmax(2,ispecb2)=NGLLY
+ ! this is useful even if it is constant because it can be zero inside the slices
+ njmin(2,ispecb2)=1
+ njmax(2,ispecb2)=NGLLY
! check for ovelap with other boundaries
nkmin_xi(2,ispecb2)=1
@@ -96,8 +97,8 @@
if(iboun(3,ispecg)) then
-! on boundary 3: ymin
- ispecb3=ispecb3+1
+ ! on boundary 3: ymin
+ ispecb3=ispecb3+1
! check for ovelap with other boundaries
nimin(1,ispecb3)=1
@@ -110,8 +111,8 @@
if(iboun(4,ispecg)) then
-! on boundary 4: ymax
- ispecb4=ispecb4+1
+ ! on boundary 4: ymax
+ ispecb4=ispecb4+1
! check for ovelap with other boundaries
nimin(2,ispecb4)=1
@@ -127,7 +128,7 @@
enddo
-! check theoretical value of elements at the bottom
+ ! check theoretical value of elements at the bottom
if(ispecb5 /= NSPEC2D_BOTTOM) &
call exit_MPI(myrank,'ispecb5 should equal NSPEC2D_BOTTOM in absorbing boundary detection')
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_jacobian_boundaries.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_jacobian_boundaries.f90 2013-07-01 04:21:23 UTC (rev 22475)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_jacobian_boundaries.f90 2013-07-01 04:44:25 UTC (rev 22476)
@@ -26,9 +26,9 @@
!=====================================================================
subroutine get_jacobian_boundaries(myrank,iboun,nspec,xstore,ystore,zstore, &
- dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
- ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, &
- nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax, &
+ dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
+ ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, &
+ nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax, &
jacobian2D_xmin,jacobian2D_xmax, &
jacobian2D_ymin,jacobian2D_ymax, &
jacobian2D_bottom,jacobian2D_top, &
@@ -42,46 +42,45 @@
include "constants.h"
- integer nspec,myrank
- integer NSPEC2D_BOTTOM,NSPEC2D_TOP,NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX
+ integer :: nspec,myrank
+ integer :: NSPEC2D_BOTTOM,NSPEC2D_TOP,NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX
- integer nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax
- integer ibelm_xmin(NSPEC2DMAX_XMIN_XMAX),ibelm_xmax(NSPEC2DMAX_XMIN_XMAX)
- integer ibelm_ymin(NSPEC2DMAX_YMIN_YMAX),ibelm_ymax(NSPEC2DMAX_YMIN_YMAX)
- integer ibelm_bottom(NSPEC2D_BOTTOM),ibelm_top(NSPEC2D_TOP)
+ integer :: nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax
+ integer :: ibelm_xmin(NSPEC2DMAX_XMIN_XMAX),ibelm_xmax(NSPEC2DMAX_XMIN_XMAX)
+ integer :: ibelm_ymin(NSPEC2DMAX_YMIN_YMAX),ibelm_ymax(NSPEC2DMAX_YMIN_YMAX)
+ integer :: ibelm_bottom(NSPEC2D_BOTTOM)
+ integer :: ibelm_top(NSPEC2D_TOP)
- logical iboun(6,nspec)
+ logical :: iboun(6,nspec)
- double precision xstore(NGLLX,NGLLY,NGLLZ,nspec)
- double precision ystore(NGLLX,NGLLY,NGLLZ,nspec)
- double precision zstore(NGLLX,NGLLY,NGLLZ,nspec)
+ double precision,dimension(NGLLX,NGLLY,NGLLZ,nspec) :: xstore,ystore,zstore
- real(kind=CUSTOM_REAL) jacobian2D_xmin(NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX)
- real(kind=CUSTOM_REAL) jacobian2D_xmax(NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX)
- real(kind=CUSTOM_REAL) jacobian2D_ymin(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)
- real(kind=CUSTOM_REAL) jacobian2D_ymax(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)
- real(kind=CUSTOM_REAL) jacobian2D_bottom(NGLLX,NGLLY,NSPEC2D_BOTTOM)
- real(kind=CUSTOM_REAL) jacobian2D_top(NGLLX,NGLLY,NSPEC2D_TOP)
+ real(kind=CUSTOM_REAL) :: jacobian2D_xmin(NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX)
+ real(kind=CUSTOM_REAL) :: jacobian2D_xmax(NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX)
+ real(kind=CUSTOM_REAL) :: jacobian2D_ymin(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)
+ real(kind=CUSTOM_REAL) :: jacobian2D_ymax(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)
+ real(kind=CUSTOM_REAL) :: jacobian2D_bottom(NGLLX,NGLLY,NSPEC2D_BOTTOM)
+ real(kind=CUSTOM_REAL) :: jacobian2D_top(NGLLX,NGLLY,NSPEC2D_TOP)
- real(kind=CUSTOM_REAL) normal_xmin(NDIM,NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX)
- real(kind=CUSTOM_REAL) normal_xmax(NDIM,NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX)
- real(kind=CUSTOM_REAL) normal_ymin(NDIM,NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)
- real(kind=CUSTOM_REAL) normal_ymax(NDIM,NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)
- real(kind=CUSTOM_REAL) normal_bottom(NDIM,NGLLX,NGLLY,NSPEC2D_BOTTOM)
- real(kind=CUSTOM_REAL) normal_top(NDIM,NGLLX,NGLLY,NSPEC2D_TOP)
+ real(kind=CUSTOM_REAL) :: normal_xmin(NDIM,NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX)
+ real(kind=CUSTOM_REAL) :: normal_xmax(NDIM,NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX)
+ real(kind=CUSTOM_REAL) :: normal_ymin(NDIM,NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)
+ real(kind=CUSTOM_REAL) :: normal_ymax(NDIM,NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)
+ real(kind=CUSTOM_REAL) :: normal_bottom(NDIM,NGLLX,NGLLY,NSPEC2D_BOTTOM)
+ real(kind=CUSTOM_REAL) :: normal_top(NDIM,NGLLX,NGLLY,NSPEC2D_TOP)
- double precision dershape2D_x(NDIM2D,NGNOD2D,NGLLY,NGLLZ)
- double precision dershape2D_y(NDIM2D,NGNOD2D,NGLLX,NGLLZ)
- double precision dershape2D_bottom(NDIM2D,NGNOD2D,NGLLX,NGLLY)
- double precision dershape2D_top(NDIM2D,NGNOD2D,NGLLX,NGLLY)
+ double precision :: dershape2D_x(NDIM2D,NGNOD2D,NGLLY,NGLLZ)
+ double precision :: dershape2D_y(NDIM2D,NGNOD2D,NGLLX,NGLLZ)
+ double precision :: dershape2D_bottom(NDIM2D,NGNOD2D,NGLLX,NGLLY)
+ double precision :: dershape2D_top(NDIM2D,NGNOD2D,NGLLX,NGLLY)
! global element numbering
- integer ispec
+ integer :: ispec
! counters to keep track of number of elements on each of the boundaries
- integer ispecb1,ispecb2,ispecb3,ispecb4,ispecb5,ispecb6
+ integer :: ispecb1,ispecb2,ispecb3,ispecb4,ispecb5,ispecb6
- double precision xelm(NGNOD2D),yelm(NGNOD2D),zelm(NGNOD2D)
+ double precision :: xelm(NGNOD2D),yelm(NGNOD2D),zelm(NGNOD2D)
! Parameters used to calculate 2D Jacobian based upon 25 GLL points
integer:: i,j,k
@@ -421,7 +420,7 @@
zelm(9)=zstore((NGLLX+1)/2,(NGLLY+1)/2,NGLLZ,ispec)
call compute_jacobian_2D(myrank,ispecb6,xelm,yelm,zelm,dershape2D_top, &
- jacobian2D_top,normal_top,NGLLX,NGLLY,NSPEC2D_TOP)
+ jacobian2D_top,normal_top,NGLLX,NGLLY,NSPEC2D_TOP)
else
! get 25 GLL points for zmax
do j = 1,NGLLY
@@ -433,8 +432,8 @@
enddo
! recalcuate jacobian according to 2D gll points
call recalc_jacobian_gll2D(myrank,ispecb6,xelm2D,yelm2D,zelm2D,&
- xigll,yigll,jacobian2D_top,normal_top,&
- NGLLX,NGLLY,NSPEC2D_TOP)
+ xigll,yigll,jacobian2D_top,normal_top,&
+ NGLLX,NGLLY,NSPEC2D_TOP)
endif
@@ -461,7 +460,8 @@
! -------------------------------------------------------
- subroutine compute_jacobian_2D(myrank,ispecb,xelm,yelm,zelm,dershape2D,jacobian2D,normal,NGLLA,NGLLB,NSPEC2DMAX_AB)
+ subroutine compute_jacobian_2D(myrank,ispecb,xelm,yelm,zelm,dershape2D, &
+ jacobian2D,normal,NGLLA,NGLLB,NSPEC2DMAX_AB)
implicit none
@@ -490,6 +490,7 @@
yeta=ZERO
zxi=ZERO
zeta=ZERO
+
do ia=1,NGNOD2D
xxi=xxi+dershape2D(1,ia,i,j)*xelm(ia)
xeta=xeta+dershape2D(2,ia,i,j)*xelm(ia)
@@ -499,16 +500,17 @@
zeta=zeta+dershape2D(2,ia,i,j)*zelm(ia)
enddo
-! calculate the unnormalized normal to the boundary
+ ! calculate the unnormalized normal to the boundary
unx=yxi*zeta-yeta*zxi
uny=zxi*xeta-zeta*xxi
unz=xxi*yeta-xeta*yxi
jacobian=dsqrt(unx**2+uny**2+unz**2)
- if(jacobian == ZERO) call exit_MPI(myrank,'2D Jacobian undefined')
-! normalize normal vector and store surface jacobian
+ if(jacobian <= ZERO) call exit_MPI(myrank,'2D Jacobian undefined')
-! distinguish between single and double precision for reals
+ ! normalize normal vector and store surface jacobian
+
+ ! distinguish between single and double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
jacobian2D(i,j,ispecb)=sngl(jacobian)
normal(1,i,j,ispecb)=sngl(unx/jacobian)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_model.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_model.f90 2013-07-01 04:21:23 UTC (rev 22475)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_model.f90 2013-07-01 04:44:25 UTC (rev 22476)
@@ -44,13 +44,8 @@
integer :: myrank,iregion_code,ispec,nspec,idoubling
- real(kind=CUSTOM_REAL) kappavstore(NGLLX,NGLLY,NGLLZ,nspec)
- real(kind=CUSTOM_REAL) kappahstore(NGLLX,NGLLY,NGLLZ,nspec)
- real(kind=CUSTOM_REAL) muvstore(NGLLX,NGLLY,NGLLZ,nspec)
- real(kind=CUSTOM_REAL) muhstore(NGLLX,NGLLY,NGLLZ,nspec)
- real(kind=CUSTOM_REAL) eta_anisostore(NGLLX,NGLLY,NGLLZ,nspec)
- real(kind=CUSTOM_REAL) rhostore(NGLLX,NGLLY,NGLLZ,nspec)
- real(kind=CUSTOM_REAL) dvpstore(NGLLX,NGLLY,NGLLZ,nspec)
+ real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ,nspec) :: kappavstore,kappahstore, &
+ muvstore,muhstore,eta_anisostore,rhostore,dvpstore
integer :: nspec_ani
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec_ani) :: &
@@ -64,7 +59,7 @@
double precision, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: xstore,ystore,zstore
- double precision rmin,rmax,RCMB,RICB,R670,RMOHO, &
+ double precision :: rmin,rmax,RCMB,RICB,R670,RMOHO, &
RTOPDDOUBLEPRIME,R600,R220,R771,R400,R120,R80,RMIDDLE_CRUST,ROCEAN
! attenuation values
@@ -72,10 +67,10 @@
double precision, dimension(N_SLS) :: tau_s
real(kind=CUSTOM_REAL), dimension(vx, vy, vz, vnspec) :: Qmu_store
real(kind=CUSTOM_REAL), dimension(N_SLS, vx, vy, vz, vnspec) :: tau_e_store
- double precision T_c_source
+ double precision :: T_c_source
logical ABSORBING_CONDITIONS
- logical elem_in_crust,elem_in_mantle
+ logical :: elem_in_crust,elem_in_mantle
! local parameters
double precision :: xmesh,ymesh,zmesh
More information about the CIG-COMMITS
mailing list