[cig-commits] r12537 - seismo/3D/SPECFEM3D_GLOBE/trunk
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Wed Aug 6 03:28:46 PDT 2008
Author: dkomati1
Date: 2008-08-06 03:28:46 -0700 (Wed, 06 Aug 2008)
New Revision: 12537
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/Makefile.in
seismo/3D/SPECFEM3D_GLOBE/trunk/attenuation_model_QRFSI12.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/compute_element_properties.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/constants.h.in
seismo/3D/SPECFEM3D_GLOBE/trunk/create_regions_mesh.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/get_MPI_1D_buffers.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/get_MPI_cutplanes_eta.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/get_MPI_cutplanes_xi.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/get_jacobian_boundaries.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/heterogen_mantle_model.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/write_movie_volume.f90
Log:
completely suppressed Cuthill-McKee sorting because of all the stability issues we had in the last few months.
Also removed useless white spaces in the process.
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/Makefile.in 2008-08-06 09:31:17 UTC (rev 12536)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/Makefile.in 2008-08-06 10:28:46 UTC (rev 12537)
@@ -96,7 +96,6 @@
$O/get_jacobian_boundaries.o \
$O/get_jacobian_discontinuities.o \
$O/get_model.o \
- $O/get_perm_cuthill_mckee.o \
$O/get_shape2D.o \
$O/get_shape3D.o \
$O/get_value_parameters.o \
@@ -543,9 +542,6 @@
$O/create_serial_name_database.o: constants.h $S/create_serial_name_database.f90
${FCCOMPILE_CHECK} -c -o $O/create_serial_name_database.o ${FCFLAGS_f90} $S/create_serial_name_database.f90
-$O/get_perm_cuthill_mckee.o: constants.h $S/get_perm_cuthill_mckee.f90
- ${FCCOMPILE_CHECK} -c -o $O/get_perm_cuthill_mckee.o ${FCFLAGS_f90} $S/get_perm_cuthill_mckee.f90
-
## use MPI here
$O/read_arrays_buffers_solver.o: constants.h $S/read_arrays_buffers_solver.f90
${MPIFCCOMPILE_CHECK} -c -o $O/read_arrays_buffers_solver.o ${FCFLAGS_f90} $S/read_arrays_buffers_solver.f90
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/attenuation_model_QRFSI12.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/attenuation_model_QRFSI12.f90 2008-08-06 09:31:17 UTC (rev 12536)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/attenuation_model_QRFSI12.f90 2008-08-06 10:28:46 UTC (rev 12537)
@@ -1,6 +1,6 @@
!=====================================================================
!
-! This file contains subroutines to read in and get values for
+! This file contains subroutines to read in and get values for
! 3-D attenuation model QRFSI12 (Dalton, Ekstrom, & Dziewonski, 2008)
!
! Last edit: Colleen Dalton, March 25, 2008
@@ -35,10 +35,10 @@
character(len=150) QRFSI12,QRFSI12_ref
! read in QRFSI12
-! hard-wire for now
+! hard-wire for now
QRFSI12='DATA/QRFSI12/QRFSI12.dat'
QRFSI12_ref='DATA/QRFSI12/ref_QRFSI12'
-
+
! get the dq model coefficients
open(unit=10,file=QRFSI12,status='old',action='read')
do k=1,NKQ
@@ -61,8 +61,8 @@
enddo
enddo
close(10)
-
-! get the depths (km) of the spline knots
+
+! get the depths (km) of the spline knots
QRFSI12_Q%spknt(1) = 24.4
QRFSI12_Q%spknt(2) = 75.0
QRFSI12_Q%spknt(3) = 150.0
@@ -71,25 +71,25 @@
QRFSI12_Q%spknt(6) = 410.0
QRFSI12_Q%spknt(7) = 530.0
QRFSI12_Q%spknt(8) = 650.0
-
+
! get the depths and 1/Q values of the reference model
open(11,file=QRFSI12_ref,status='old',action='read')
do j=1,NDEPTHS_REFQ
read(11,*)QRFSI12_Q%refdepth(j),QRFSI12_Q%refqmu(j)
enddo
close(11)
-
-
+
+
end subroutine read_atten_model_3D_QRFSI12
!----------------------------------
!----------------------------------
-
+
subroutine attenuation_model_3D_QRFSI12(radius,theta,phi,Qmu,QRFSI12_Q,idoubling)
implicit none
- include "constants.h"
+ include "constants.h"
! atten_model_QRFSI12_variables
type atten_model_QRFSI12_variables
@@ -112,16 +112,16 @@
real(kind=4) xlmvec(NSQ),wk1(NSQ),wk2(NSQ),wk3(NSQ)
double precision, parameter :: rmoho = 6371.0-24.4
double precision, parameter :: rcmb = 3480.0
-
- !in Colleen's original code theta refers to the latitude. Here we have redefined theta to be colatitude
- !to agree with the rest of specfem
+
+ !in Colleen's original code theta refers to the latitude. Here we have redefined theta to be colatitude
+ !to agree with the rest of specfem
! print *,'entering QRFSI12 subroutine'
ylat=90.0d0-theta
- xlon=phi
+ xlon=phi
if(idoubling == IFLAG_CRUST .or. radius >= rmoho) then
- Qmu = 600.0d0
+ Qmu = 600.0d0
! print *,'QRFSI12: we are in the crust'
else if(idoubling == IFLAG_INNER_CORE_NORMAL .or. idoubling == IFLAG_MIDDLE_CENTRAL_CUBE .or. &
idoubling == IFLAG_BOTTOM_CENTRAL_CUBE .or. idoubling == IFLAG_TOP_CENTRAL_CUBE .or. &
@@ -143,7 +143,7 @@
if(ifnd == 0)then
write(6,"('problem finding reference Q value at depth: ',f8.3)") depth
stop
- endif
+ endif
smallq_ref=QRFSI12_Q%refqmu(ifnd)
smallq = smallq_ref
@@ -175,12 +175,12 @@
if(Qmu >= ATTENUATION_COMP_MAXIMUM) Qmu = 0.99d0*ATTENUATION_COMP_MAXIMUM
endif
-
+
end subroutine attenuation_model_3D_QRFSI12
!----------------------------------
!----------------------------------
-
+
!!$ subroutine vbspl(x,np,xarr,splcon,splcond)
!!$!
!!$!---- this subroutine returns the spline contributions at a particular value of x
@@ -584,7 +584,7 @@
!!$ end subroutine ylm
!!$ subroutine legndr(THETA,L,M,X,XP,XCOSEC)
-!!$ implicit none
+!!$ implicit none
!!$
!!$ integer :: L,M,i,k,LP1,MP1
!!$ real(kind=4) :: THETA,X,XP,XCOSEC,SFL3
@@ -658,4 +658,4 @@
!!$ enddo
!!$ RETURN
!!$ end subroutine legndr
-
+
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/compute_element_properties.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/compute_element_properties.f90 2008-08-06 09:31:17 UTC (rev 12536)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/compute_element_properties.f90 2008-08-06 10:28:46 UTC (rev 12537)
@@ -29,7 +29,7 @@
subroutine compute_element_properties(ispec,iregion_code,idoubling, &
xstore,ystore,zstore,nspec, &
nspl,rspl,espl,espl2,ELLIPTICITY,TOPOGRAPHY,TRANSVERSE_ISOTROPY, &
- ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE,ISOTROPIC_3D_MANTLE, &
+ ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE,ISOTROPIC_3D_MANTLE, &
HETEROGEN_3D_MANTLE,CRUSTAL,ONE_CRUST, &
myrank,ibathy_topo,ATTENUATION,ATTENUATION_3D, &
ABSORBING_CONDITIONS,REFERENCE_1D_MODEL,THREE_D_MODEL, &
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/constants.h.in 2008-08-06 09:31:17 UTC (rev 12536)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/constants.h.in 2008-08-06 10:28:46 UTC (rev 12537)
@@ -497,11 +497,3 @@
integer, parameter :: NSPEC_DOUBLING_BASICBRICK = 8
integer, parameter :: NGLOB_DOUBLING_BASICBRICK = 27
-! for Cuthill-McKee (1969) permutation
- logical, parameter :: PERFORM_CUTHILL_MCKEE = .false.
- integer, parameter :: NGNOD_HEXAHEDRA = 8
-! perform classical or multi-level Cuthill-McKee ordering
- logical, parameter :: CMcK_MULTI = .false.
-! maximum size if multi-level Cuthill-McKee ordering
- integer, parameter :: LIMIT_MULTI_CUTHILL = 50
-
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/create_regions_mesh.f90 2008-08-06 09:31:17 UTC (rev 12536)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/create_regions_mesh.f90 2008-08-06 10:28:46 UTC (rev 12537)
@@ -36,7 +36,7 @@
ELLIPTICITY,TOPOGRAPHY,TRANSVERSE_ISOTROPY, &
ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE,ISOTROPIC_3D_MANTLE,HETEROGEN_3D_MANTLE,CRUSTAL,ONE_CRUST, &
NPROC_XI,NPROC_ETA,NSPEC2D_XI_FACE, &
- NSPEC2D_ETA_FACE,NSPEC1D_RADIAL_CORNER,NGLOB1D_RADIAL_CORNER,NGLOB2DMAX_XY, &
+ NSPEC2D_ETA_FACE,NSPEC1D_RADIAL_CORNER,NGLOB1D_RADIAL_CORNER, &
myrank,LOCAL_PATH,OCEANS,ibathy_topo, &
rotation_matrix,ANGULAR_WIDTH_XI_RAD,ANGULAR_WIDTH_ETA_RAD,&
ATTENUATION,ATTENUATION_3D,SAVE_MESH_FILES, &
@@ -507,8 +507,6 @@
integer NUMBER_OF_MESH_LAYERS,layer_shift,cpt,first_layer_aniso,last_layer_aniso,FIRST_ELT_NON_ANISO
double precision, dimension(:,:), allocatable :: stretch_tab
- integer :: NGLOB2DMAX_XY
-
integer :: nb_layer_above_aniso,FIRST_ELT_ABOVE_ANISO
integer, parameter :: maxker=200
@@ -560,14 +558,6 @@
double precision r_moho,r_400,r_670
logical :: is_superbrick
-! added for Cuthill McKee permutation
- integer, dimension(:), allocatable :: perm,perm_tmp,temp_array_1D_int
- logical, dimension(:,:), allocatable :: temp_array_2D_log
- integer, dimension(:,:,:,:), allocatable :: temp_array_int
- real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: temp_array_real
- double precision, dimension(:,:,:,:), allocatable :: temp_array_dble
- double precision, dimension(:,:,:,:,:), allocatable :: temp_array_dble_5dim
-
! the height at which the central cube is cut
integer :: nz_inf_limit
@@ -1468,10 +1458,12 @@
! arrays locval(npointot) and ifseg(npointot) used to save memory
call get_MPI_cutplanes_xi(myrank,prname,nspec,iMPIcut_xi,ibool, &
xstore,ystore,zstore,ifseg,npointot, &
- NSPEC2D_ETA_FACE,iregion_code,NGLOB2DMAX_XY)
+ NSPEC2D_ETA_FACE,iregion_code)
+
call get_MPI_cutplanes_eta(myrank,prname,nspec,iMPIcut_eta,ibool, &
xstore,ystore,zstore,ifseg,npointot, &
- NSPEC2D_XI_FACE,iregion_code,NGLOB2DMAX_XY)
+ NSPEC2D_XI_FACE,iregion_code)
+
call get_MPI_1D_buffers(myrank,prname,nspec,iMPIcut_xi,iMPIcut_eta,ibool,idoubling, &
xstore,ystore,zstore,ifseg,npointot, &
NSPEC1D_RADIAL_CORNER,NGLOB1D_RADIAL_CORNER,iregion_code)
@@ -1509,140 +1501,6 @@
! should be zero in all the regions except in the mantle
nspec_tiso = count(idoubling(1:nspec) == IFLAG_220_80) + count(idoubling(1:nspec) == IFLAG_80_MOHO)
-! ***************************************************
-! Cuthill McKee permutation
-! ***************************************************
- if (PERFORM_CUTHILL_MCKEE) then
- allocate(perm(nspec))
- if(iregion_code == IREGION_CRUST_MANTLE) then
- ! do not permute anisotropic elements
- perm(1:FIRST_ELT_NON_ANISO-1) = (/ (i,i=1,FIRST_ELT_NON_ANISO-1) /)
-
- ! no more connectivity between layers below and above the anisotropic layers => 2 permutations
- ! permute the bottom of the region : below the aniso layers
- allocate(perm_tmp(FIRST_ELT_ABOVE_ANISO-FIRST_ELT_NON_ANISO))
- call get_perm(ibool(:,:,:,FIRST_ELT_NON_ANISO:FIRST_ELT_ABOVE_ANISO-1),perm_tmp,LIMIT_MULTI_CUTHILL,&
- (FIRST_ELT_ABOVE_ANISO-FIRST_ELT_NON_ANISO),nglob,.true.,.false.)
- perm(FIRST_ELT_NON_ANISO:FIRST_ELT_ABOVE_ANISO-1) = perm_tmp(:)+(FIRST_ELT_NON_ANISO-1)
- deallocate(perm_tmp)
-
- ! permute the top of the region : above the aniso layers
- allocate(perm_tmp(nspec-FIRST_ELT_ABOVE_ANISO+1))
- call get_perm(ibool(:,:,:,FIRST_ELT_ABOVE_ANISO:nspec),perm_tmp,LIMIT_MULTI_CUTHILL,&
- (nspec-FIRST_ELT_ABOVE_ANISO+1),nglob,.true.,.false.)
- perm(FIRST_ELT_ABOVE_ANISO:nspec) = perm_tmp(:)+(FIRST_ELT_ABOVE_ANISO-1)
- deallocate(perm_tmp)
- else
- ! the 3 last parameters are : PERFORM_CUTHILL_MCKEE,INVERSE,FACE
- call get_perm(ibool,perm,LIMIT_MULTI_CUTHILL,nspec,nglob,.true.,.false.)
- endif
-
- ! permutation of xstore, ystore, zstore, rhostore, kappavstore, kappahstore,
- ! muvstore, muhstore, eta_anisostore, xixstore, xiystore, xizstore, etaxstore,
- ! etaystore, etazstore, gammaxstore, gammaystore, gammazstore, no more jacobianstore
-
- allocate(temp_array_dble(NGLLX,NGLLY,NGLLZ,nspec))
- if(ATTENUATION .and. ATTENUATION_3D) then
- call permute_elements_dble(Qmu_store,temp_array_dble,perm,nspec)
- allocate(temp_array_dble_5dim(N_SLS,NGLLX,NGLLY,NGLLZ,nspec))
- temp_array_dble_5dim(:,:,:,:,:) = tau_e_store(:,:,:,:,:)
- do i = 1,nspec
- tau_e_store(:,:,:,:,perm(i)) = temp_array_dble_5dim(:,:,:,:,i)
- enddo
- deallocate(temp_array_dble_5dim)
- endif
- call permute_elements_dble(xstore,temp_array_dble,perm,nspec)
- call permute_elements_dble(ystore,temp_array_dble,perm,nspec)
- call permute_elements_dble(zstore,temp_array_dble,perm,nspec)
- deallocate(temp_array_dble)
-
-
- allocate(temp_array_real(NGLLX,NGLLY,NGLLZ,nspec))
- if(NCHUNKS /= 6) then
- call permute_elements_real(rho_vp,temp_array_real,perm,nspec)
- call permute_elements_real(rho_vs,temp_array_real,perm,nspec)
- endif
- if((ANISOTROPIC_INNER_CORE .and. iregion_code == IREGION_INNER_CORE) .or. &
- (ANISOTROPIC_3D_MANTLE .and. iregion_code == IREGION_CRUST_MANTLE)) then
- call permute_elements_real(c11store,temp_array_real,perm,nspec)
- call permute_elements_real(c12store,temp_array_real,perm,nspec)
- call permute_elements_real(c13store,temp_array_real,perm,nspec)
- call permute_elements_real(c14store,temp_array_real,perm,nspec)
- call permute_elements_real(c15store,temp_array_real,perm,nspec)
- call permute_elements_real(c16store,temp_array_real,perm,nspec)
- call permute_elements_real(c22store,temp_array_real,perm,nspec)
- call permute_elements_real(c23store,temp_array_real,perm,nspec)
- call permute_elements_real(c24store,temp_array_real,perm,nspec)
- call permute_elements_real(c25store,temp_array_real,perm,nspec)
- call permute_elements_real(c26store,temp_array_real,perm,nspec)
- call permute_elements_real(c33store,temp_array_real,perm,nspec)
- call permute_elements_real(c34store,temp_array_real,perm,nspec)
- call permute_elements_real(c35store,temp_array_real,perm,nspec)
- call permute_elements_real(c36store,temp_array_real,perm,nspec)
- call permute_elements_real(c44store,temp_array_real,perm,nspec)
- call permute_elements_real(c45store,temp_array_real,perm,nspec)
- call permute_elements_real(c46store,temp_array_real,perm,nspec)
- call permute_elements_real(c55store,temp_array_real,perm,nspec)
- call permute_elements_real(c56store,temp_array_real,perm,nspec)
- call permute_elements_real(c66store,temp_array_real,perm,nspec)
- endif
- call permute_elements_real(rhostore,temp_array_real,perm,nspec)
- call permute_elements_real(kappavstore,temp_array_real,perm,nspec)
- call permute_elements_real(kappahstore,temp_array_real,perm,nspec)
- call permute_elements_real(muvstore,temp_array_real,perm,nspec)
- call permute_elements_real(muhstore,temp_array_real,perm,nspec)
- call permute_elements_real(eta_anisostore,temp_array_real,perm,nspec)
- call permute_elements_real(xixstore,temp_array_real,perm,nspec)
- call permute_elements_real(xiystore,temp_array_real,perm,nspec)
- call permute_elements_real(xizstore,temp_array_real,perm,nspec)
- call permute_elements_real(etaxstore,temp_array_real,perm,nspec)
- call permute_elements_real(etaystore,temp_array_real,perm,nspec)
- call permute_elements_real(etazstore,temp_array_real,perm,nspec)
- call permute_elements_real(gammaxstore,temp_array_real,perm,nspec)
- call permute_elements_real(gammaystore,temp_array_real,perm,nspec)
- call permute_elements_real(gammazstore,temp_array_real,perm,nspec)
- deallocate(temp_array_real)
-
- ! permutation of ibool
- allocate(temp_array_int(NGLLX,NGLLY,NGLLZ,nspec))
- call permute_elements_integer(ibool,temp_array_int,perm,nspec)
- deallocate(temp_array_int)
-
- ! permutation of iMPIcut_*
- allocate(temp_array_2D_log(2,nspec))
- temp_array_2D_log(:,:) = iMPIcut_xi(:,:)
- do i = 1,nspec
- iMPIcut_xi(:,perm(i)) = temp_array_2D_log(:,i)
- enddo
- temp_array_2D_log(:,:) = iMPIcut_eta(:,:)
- do i = 1,nspec
- iMPIcut_eta(:,perm(i)) = temp_array_2D_log(:,i)
- enddo
- deallocate(temp_array_2D_log)
-
- ! permutation of iboun
- allocate(temp_array_2D_log(6,nspec))
- temp_array_2D_log(:,:) = iboun(:,:)
- do i = 1,nspec
- iboun(:,perm(i)) = temp_array_2D_log(:,i)
- enddo
- deallocate(temp_array_2D_log)
-
- ! permutation of idoubling
- allocate(temp_array_1D_int(nspec))
- temp_array_1D_int(:) = idoubling(:)
- do i = 1,nspec
- idoubling(perm(i)) = temp_array_1D_int(i)
- enddo
- deallocate(temp_array_1D_int)
-
- deallocate(perm)
- endif
-
-! ***************************************************
-! end of Cuthill McKee permutation
-! ***************************************************
-
call 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, &
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/get_MPI_1D_buffers.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/get_MPI_1D_buffers.f90 2008-08-06 09:31:17 UTC (rev 12536)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/get_MPI_1D_buffers.f90 2008-08-06 10:28:46 UTC (rev 12537)
@@ -35,7 +35,7 @@
include "constants.h"
- integer nspec,myrank,nglob,ipoin1D,iregion
+ integer nspec,myrank,iregion
integer, dimension(MAX_NUM_REGIONS,NB_SQUARE_CORNERS) :: NSPEC1D_RADIAL_CORNER,NGLOB1D_RADIAL_CORNER
logical iMPIcut_xi(2,nspec)
@@ -62,29 +62,6 @@
! processor identification
character(len=150) prname
-! arrays for sorting routine
- integer, dimension(:), allocatable :: ind,ninseg,iglob,locval,iwork
- logical, dimension(:), allocatable :: ifseg
- double precision, dimension(:), allocatable :: work
- integer, dimension(:), allocatable :: ibool_selected
- double precision, dimension(:), allocatable :: xstore_selected,ystore_selected,zstore_selected
-
-! allocate arrays for message buffers with maximum size
-! define maximum size for message buffers
- if (PERFORM_CUTHILL_MCKEE) then
- allocate(ibool_selected(maxval(NGLOB1D_RADIAL_CORNER(iregion,:))))
- allocate(xstore_selected(maxval(NGLOB1D_RADIAL_CORNER(iregion,:))))
- allocate(ystore_selected(maxval(NGLOB1D_RADIAL_CORNER(iregion,:))))
- allocate(zstore_selected(maxval(NGLOB1D_RADIAL_CORNER(iregion,:))))
- allocate(ind(maxval(NGLOB1D_RADIAL_CORNER(iregion,:))))
- allocate(ninseg(maxval(NGLOB1D_RADIAL_CORNER(iregion,:))))
- allocate(iglob(maxval(NGLOB1D_RADIAL_CORNER(iregion,:))))
- allocate(locval(maxval(NGLOB1D_RADIAL_CORNER(iregion,:))))
- allocate(ifseg(maxval(NGLOB1D_RADIAL_CORNER(iregion,:))))
- allocate(iwork(maxval(NGLOB1D_RADIAL_CORNER(iregion,:))))
- allocate(work(maxval(NGLOB1D_RADIAL_CORNER(iregion,:))))
- endif
-
! 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
@@ -123,29 +100,13 @@
if(.not. mask_ibool(ibool(ix,iy,iz,ispec))) then
mask_ibool(ibool(ix,iy,iz,ispec)) = .true.
npoin1D = npoin1D + 1
- if (PERFORM_CUTHILL_MCKEE) then
- ibool_selected(npoin1D) = ibool(ix,iy,iz,ispec)
- xstore_selected(npoin1D) = xstore(ix,iy,iz,ispec)
- ystore_selected(npoin1D) = ystore(ix,iy,iz,ispec)
- zstore_selected(npoin1D) = zstore(ix,iy,iz,ispec)
- else
- write(10,*) ibool(ix,iy,iz,ispec), xstore(ix,iy,iz,ispec), &
- ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
- endif
+ write(10,*) ibool(ix,iy,iz,ispec), xstore(ix,iy,iz,ispec), &
+ ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
endif
enddo
endif
enddo
- if (PERFORM_CUTHILL_MCKEE) then
- call sort_array_coordinates(npoin1D,xstore_selected,ystore_selected,zstore_selected, &
- ibool_selected,iglob,locval,ifseg,nglob,ind,ninseg,iwork,work)
-
- do ipoin1D=1,npoin1D
- write(10,*) ibool_selected(ipoin1D), xstore_selected(ipoin1D), &
- ystore_selected(ipoin1D),zstore_selected(ipoin1D)
- enddo
- endif
! put flag to indicate end of the list of points
write(10,*) '0 0 0. 0. 0.'
@@ -188,30 +149,13 @@
if(.not. mask_ibool(ibool(ix,iy,iz,ispec))) then
mask_ibool(ibool(ix,iy,iz,ispec)) = .true.
npoin1D = npoin1D + 1
- if (PERFORM_CUTHILL_MCKEE) then
- ibool_selected(npoin1D) = ibool(ix,iy,iz,ispec)
- xstore_selected(npoin1D) = xstore(ix,iy,iz,ispec)
- ystore_selected(npoin1D) = ystore(ix,iy,iz,ispec)
- zstore_selected(npoin1D) = zstore(ix,iy,iz,ispec)
- else
- write(10,*) ibool(ix,iy,iz,ispec), xstore(ix,iy,iz,ispec), &
- ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
- endif
+ write(10,*) ibool(ix,iy,iz,ispec), xstore(ix,iy,iz,ispec), &
+ ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
endif
enddo
endif
enddo
- if (PERFORM_CUTHILL_MCKEE) then
- call sort_array_coordinates(npoin1D,xstore_selected,ystore_selected,zstore_selected, &
- ibool_selected,iglob,locval,ifseg,nglob,ind,ninseg,iwork,work)
-
- do ipoin1D=1,npoin1D
- write(10,*) ibool_selected(ipoin1D), xstore_selected(ipoin1D), &
- ystore_selected(ipoin1D),zstore_selected(ipoin1D)
- enddo
- endif
-
! put flag to indicate end of the list of points
write(10,*) '0 0 0. 0. 0.'
@@ -264,30 +208,13 @@
if(.not. mask_ibool(ibool(ix,iy,iz,ispec))) then
mask_ibool(ibool(ix,iy,iz,ispec)) = .true.
npoin1D = npoin1D + 1
- if (PERFORM_CUTHILL_MCKEE) then
- ibool_selected(npoin1D) = ibool(ix,iy,iz,ispec)
- xstore_selected(npoin1D) = xstore(ix,iy,iz,ispec)
- ystore_selected(npoin1D) = ystore(ix,iy,iz,ispec)
- zstore_selected(npoin1D) = zstore(ix,iy,iz,ispec)
- else
- write(10,*) ibool(ix,iy,iz,ispec), xstore(ix,iy,iz,ispec), &
- ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
- endif
+ write(10,*) ibool(ix,iy,iz,ispec), xstore(ix,iy,iz,ispec), &
+ ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
endif
enddo
endif
enddo
- if (PERFORM_CUTHILL_MCKEE) then
- call sort_array_coordinates(npoin1D,xstore_selected,ystore_selected,zstore_selected, &
- ibool_selected,iglob,locval,ifseg,nglob,ind,ninseg,iwork,work)
-
- do ipoin1D=1,npoin1D
- write(10,*) ibool_selected(ipoin1D), xstore_selected(ipoin1D), &
- ystore_selected(ipoin1D),zstore_selected(ipoin1D)
- enddo
- endif
-
! put flag to indicate end of the list of points
write(10,*) '0 0 0. 0. 0.'
@@ -336,30 +263,13 @@
if(.not. mask_ibool(ibool(ix,iy,iz,ispec))) then
mask_ibool(ibool(ix,iy,iz,ispec)) = .true.
npoin1D = npoin1D + 1
- if (PERFORM_CUTHILL_MCKEE) then
- ibool_selected(npoin1D) = ibool(ix,iy,iz,ispec)
- xstore_selected(npoin1D) = xstore(ix,iy,iz,ispec)
- ystore_selected(npoin1D) = ystore(ix,iy,iz,ispec)
- zstore_selected(npoin1D) = zstore(ix,iy,iz,ispec)
- else
- write(10,*) ibool(ix,iy,iz,ispec), xstore(ix,iy,iz,ispec), &
- ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
- endif
+ write(10,*) ibool(ix,iy,iz,ispec), xstore(ix,iy,iz,ispec), &
+ ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
endif
enddo
endif
enddo
- if (PERFORM_CUTHILL_MCKEE) then
- call sort_array_coordinates(npoin1D,xstore_selected,ystore_selected,zstore_selected, &
- ibool_selected,iglob,locval,ifseg,nglob,ind,ninseg,iwork,work)
-
- do ipoin1D=1,npoin1D
- write(10,*) ibool_selected(ipoin1D), xstore_selected(ipoin1D), &
- ystore_selected(ipoin1D),zstore_selected(ipoin1D)
- enddo
- endif
-
! put flag to indicate end of the list of points
write(10,*) '0 0 0. 0. 0.'
@@ -372,19 +282,5 @@
if(ispeccount /= NSPEC1D_RADIAL_CORNER(iregion,3) .or. npoin1D /= NGLOB1D_RADIAL_CORNER(iregion,3)) &
call exit_MPI(myrank,'error MPI 1D buffer detection in xi=right')
- if (PERFORM_CUTHILL_MCKEE) then
- deallocate(ibool_selected)
- deallocate(xstore_selected)
- deallocate(ystore_selected)
- deallocate(zstore_selected)
- deallocate(ind)
- deallocate(ninseg)
- deallocate(iglob)
- deallocate(locval)
- deallocate(ifseg)
- deallocate(iwork)
- deallocate(work)
- endif
-
end subroutine get_MPI_1D_buffers
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/get_MPI_cutplanes_eta.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/get_MPI_cutplanes_eta.f90 2008-08-06 09:31:17 UTC (rev 12536)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/get_MPI_cutplanes_eta.f90 2008-08-06 10:28:46 UTC (rev 12537)
@@ -27,7 +27,7 @@
subroutine get_MPI_cutplanes_eta(myrank,prname,nspec,iMPIcut_eta,ibool, &
xstore,ystore,zstore,mask_ibool,npointot, &
- NSPEC2D_XI_FACE,iregion,NGLOB2DMAX_XY)
+ NSPEC2D_XI_FACE,iregion)
! 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,nglob,ipoin2D,NGLOB2DMAX_XY,iregion
+ integer nspec,myrank,iregion
integer, dimension(MAX_NUM_REGIONS,NB_SQUARE_EDGES_ONEDIR) :: NSPEC2D_XI_FACE
@@ -64,30 +64,6 @@
! processor identification
character(len=150) prname
-! arrays for sorting routine
- integer, dimension(:), allocatable :: ind,ninseg,iglob,locval,iwork
- logical, dimension(:), allocatable :: ifseg
- double precision, dimension(:), allocatable :: work
- integer, dimension(:), allocatable :: ibool_selected
- double precision, dimension(:), allocatable :: xstore_selected,ystore_selected,zstore_selected
-
-
-! allocate arrays for message buffers with maximum size
-! define maximum size for message buffers
- if (PERFORM_CUTHILL_MCKEE) then
- allocate(ibool_selected(NGLOB2DMAX_XY))
- allocate(xstore_selected(NGLOB2DMAX_XY))
- allocate(ystore_selected(NGLOB2DMAX_XY))
- allocate(zstore_selected(NGLOB2DMAX_XY))
- allocate(ind(NGLOB2DMAX_XY))
- allocate(ninseg(NGLOB2DMAX_XY))
- allocate(iglob(NGLOB2DMAX_XY))
- allocate(locval(NGLOB2DMAX_XY))
- allocate(ifseg(NGLOB2DMAX_XY))
- allocate(iwork(NGLOB2DMAX_XY))
- allocate(work(NGLOB2DMAX_XY))
- endif
-
! theoretical number of surface elements in the buffers
! cut planes along eta=constant correspond to XI faces
nspec2Dtheor = NSPEC2D_XI_FACE(iregion,1)
@@ -122,31 +98,14 @@
if(.not. mask_ibool(ibool(ix,iy,iz,ispec))) then
mask_ibool(ibool(ix,iy,iz,ispec)) = .true.
npoin2D_eta = npoin2D_eta + 1
- if (PERFORM_CUTHILL_MCKEE) then
- ibool_selected(npoin2D_eta) = ibool(ix,iy,iz,ispec)
- xstore_selected(npoin2D_eta) = xstore(ix,iy,iz,ispec)
- ystore_selected(npoin2D_eta) = ystore(ix,iy,iz,ispec)
- zstore_selected(npoin2D_eta) = zstore(ix,iy,iz,ispec)
- else
- write(10,*) ibool(ix,iy,iz,ispec), xstore(ix,iy,iz,ispec), &
- ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
- endif
+ write(10,*) ibool(ix,iy,iz,ispec), xstore(ix,iy,iz,ispec), &
+ ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
endif
enddo
enddo
endif
enddo
- if (PERFORM_CUTHILL_MCKEE) then
- call sort_array_coordinates(npoin2D_eta,xstore_selected,ystore_selected,zstore_selected, &
- ibool_selected,iglob,locval,ifseg,nglob,ind,ninseg,iwork,work)
-
- do ipoin2D=1,npoin2D_eta
- write(10,*) ibool_selected(ipoin2D), xstore_selected(ipoin2D), &
- ystore_selected(ipoin2D),zstore_selected(ipoin2D)
- enddo
- endif
-
! put flag to indicate end of the list of points
write(10,*) '0 0 0. 0. 0.'
@@ -186,31 +145,14 @@
if(.not. mask_ibool(ibool(ix,iy,iz,ispec))) then
mask_ibool(ibool(ix,iy,iz,ispec)) = .true.
npoin2D_eta = npoin2D_eta + 1
- if (PERFORM_CUTHILL_MCKEE) then
- ibool_selected(npoin2D_eta) = ibool(ix,iy,iz,ispec)
- xstore_selected(npoin2D_eta) = xstore(ix,iy,iz,ispec)
- ystore_selected(npoin2D_eta) = ystore(ix,iy,iz,ispec)
- zstore_selected(npoin2D_eta) = zstore(ix,iy,iz,ispec)
- else
- write(10,*) ibool(ix,iy,iz,ispec), xstore(ix,iy,iz,ispec), &
- ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
- endif
+ write(10,*) ibool(ix,iy,iz,ispec), xstore(ix,iy,iz,ispec), &
+ ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
endif
enddo
enddo
endif
enddo
- if (PERFORM_CUTHILL_MCKEE) then
- call sort_array_coordinates(npoin2D_eta,xstore_selected,ystore_selected,zstore_selected, &
- ibool_selected,iglob,locval,ifseg,nglob,ind,ninseg,iwork,work)
-
- do ipoin2D=1,npoin2D_eta
- write(10,*) ibool_selected(ipoin2D), xstore_selected(ipoin2D), &
- ystore_selected(ipoin2D),zstore_selected(ipoin2D)
- enddo
- endif
-
! put flag to indicate end of the list of points
write(10,*) '0 0 0. 0. 0.'
@@ -222,19 +164,5 @@
! compare number of surface elements detected to analytical value
if(ispecc2 /= nspec2Dtheor) call exit_MPI(myrank,'error MPI cut-planes detection in eta=right')
- if (PERFORM_CUTHILL_MCKEE) then
- deallocate(ibool_selected)
- deallocate(xstore_selected)
- deallocate(ystore_selected)
- deallocate(zstore_selected)
- deallocate(ind)
- deallocate(ninseg)
- deallocate(iglob)
- deallocate(locval)
- deallocate(ifseg)
- deallocate(iwork)
- deallocate(work)
- endif
-
end subroutine get_MPI_cutplanes_eta
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/get_MPI_cutplanes_xi.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/get_MPI_cutplanes_xi.f90 2008-08-06 09:31:17 UTC (rev 12536)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/get_MPI_cutplanes_xi.f90 2008-08-06 10:28:46 UTC (rev 12537)
@@ -27,7 +27,7 @@
subroutine get_MPI_cutplanes_xi(myrank,prname,nspec,iMPIcut_xi,ibool, &
xstore,ystore,zstore,mask_ibool,npointot, &
- NSPEC2D_ETA_FACE,iregion,NGLOB2DMAX_XY)
+ NSPEC2D_ETA_FACE,iregion)
! 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,nglob,ipoin2D,iregion
+ integer nspec,myrank,iregion
integer, dimension(MAX_NUM_REGIONS,NB_SQUARE_EDGES_ONEDIR) :: NSPEC2D_ETA_FACE
logical iMPIcut_xi(2,nspec)
@@ -63,31 +63,6 @@
! processor identification
character(len=150) prname,errmsg
-! arrays for sorting routine
- integer, dimension(:), allocatable :: ind,ninseg,iglob,locval,iwork
- logical, dimension(:), allocatable :: ifseg
- double precision, dimension(:), allocatable :: work
- integer NGLOB2DMAX_XY
- integer, dimension(:), allocatable :: ibool_selected
- double precision, dimension(:), allocatable :: xstore_selected,ystore_selected,zstore_selected
-
-! allocate arrays for message buffers with maximum size
-! define maximum size for message buffers
- if (PERFORM_CUTHILL_MCKEE) then
- allocate(ibool_selected(NGLOB2DMAX_XY))
- allocate(xstore_selected(NGLOB2DMAX_XY))
- allocate(ystore_selected(NGLOB2DMAX_XY))
- allocate(zstore_selected(NGLOB2DMAX_XY))
- allocate(ind(NGLOB2DMAX_XY))
- allocate(ninseg(NGLOB2DMAX_XY))
- allocate(iglob(NGLOB2DMAX_XY))
- allocate(locval(NGLOB2DMAX_XY))
- allocate(ifseg(NGLOB2DMAX_XY))
- allocate(iwork(NGLOB2DMAX_XY))
- allocate(work(NGLOB2DMAX_XY))
- endif
-
-
! theoretical number of surface elements in the buffers
! cut planes along xi=constant correspond to ETA faces
nspec2Dtheor = NSPEC2D_ETA_FACE(iregion,1)
@@ -121,31 +96,14 @@
if(.not. mask_ibool(ibool(ix,iy,iz,ispec))) then
mask_ibool(ibool(ix,iy,iz,ispec)) = .true.
npoin2D_xi = npoin2D_xi + 1
- if (PERFORM_CUTHILL_MCKEE) then
- ibool_selected(npoin2D_xi) = ibool(ix,iy,iz,ispec)
- xstore_selected(npoin2D_xi) = xstore(ix,iy,iz,ispec)
- ystore_selected(npoin2D_xi) = ystore(ix,iy,iz,ispec)
- zstore_selected(npoin2D_xi) = zstore(ix,iy,iz,ispec)
- else
- write(10,*) ibool(ix,iy,iz,ispec), xstore(ix,iy,iz,ispec), &
- ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
- endif
+ write(10,*) ibool(ix,iy,iz,ispec), xstore(ix,iy,iz,ispec), &
+ ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
endif
enddo
enddo
endif
enddo
- if (PERFORM_CUTHILL_MCKEE) then
- call sort_array_coordinates(npoin2D_xi,xstore_selected,ystore_selected,zstore_selected, &
- ibool_selected,iglob,locval,ifseg,nglob,ind,ninseg,iwork,work)
-
- do ipoin2D=1,npoin2D_xi
- write(10,*) ibool_selected(ipoin2D), xstore_selected(ipoin2D), &
- ystore_selected(ipoin2D),zstore_selected(ipoin2D)
- enddo
- endif
-
! put flag to indicate end of the list of points
write(10,*) '0 0 0. 0. 0.'
@@ -187,32 +145,14 @@
if(.not. mask_ibool(ibool(ix,iy,iz,ispec))) then
mask_ibool(ibool(ix,iy,iz,ispec)) = .true.
npoin2D_xi = npoin2D_xi + 1
- if (PERFORM_CUTHILL_MCKEE) then
- ibool_selected(npoin2D_xi) = ibool(ix,iy,iz,ispec)
- xstore_selected(npoin2D_xi) = xstore(ix,iy,iz,ispec)
- ystore_selected(npoin2D_xi) = ystore(ix,iy,iz,ispec)
- zstore_selected(npoin2D_xi) = zstore(ix,iy,iz,ispec)
- else
- write(10,*) ibool(ix,iy,iz,ispec), xstore(ix,iy,iz,ispec), &
- ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
- endif
+ write(10,*) ibool(ix,iy,iz,ispec), xstore(ix,iy,iz,ispec), &
+ ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec)
endif
enddo
enddo
endif
enddo
- if (PERFORM_CUTHILL_MCKEE) then
- call sort_array_coordinates(npoin2D_xi,xstore_selected,ystore_selected,zstore_selected, &
- ibool_selected,iglob,locval,ifseg,nglob,ind,ninseg,iwork,work)
-
- do ipoin2D=1,npoin2D_xi
- write(10,*) ibool_selected(ipoin2D), xstore_selected(ipoin2D), &
- ystore_selected(ipoin2D),zstore_selected(ipoin2D)
- enddo
- endif
-
-
! put flag to indicate end of the list of points
write(10,*) '0 0 0. 0. 0.'
@@ -227,19 +167,5 @@
call exit_MPI(myrank,errmsg)
endif
- if (PERFORM_CUTHILL_MCKEE) then
- deallocate(ibool_selected)
- deallocate(xstore_selected)
- deallocate(ystore_selected)
- deallocate(zstore_selected)
- deallocate(ind)
- deallocate(ninseg)
- deallocate(iglob)
- deallocate(locval)
- deallocate(ifseg)
- deallocate(iwork)
- deallocate(work)
- endif
-
end subroutine get_MPI_cutplanes_xi
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/get_jacobian_boundaries.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/get_jacobian_boundaries.f90 2008-08-06 09:31:17 UTC (rev 12536)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/get_jacobian_boundaries.f90 2008-08-06 10:28:46 UTC (rev 12537)
@@ -137,10 +137,6 @@
call compute_jacobian_2D(myrank,ispecb1,xelm,yelm,zelm,dershape2D_x, &
jacobian2D_xmin,normal_xmin,NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX)
- if (PERFORM_CUTHILL_MCKEE) then
- call sort_arrays_for_cuthill (ispecb1,xstore,ystore,zstore,ibelm_xmin,normal_xmin,&
- jacobian2D_xmin,NSPEC2DMAX_XMIN_XMAX,NGLLY,NGLLZ,nspec)
- endif
endif
! on boundary: xmax
@@ -182,10 +178,6 @@
call compute_jacobian_2D(myrank,ispecb2,xelm,yelm,zelm,dershape2D_x, &
jacobian2D_xmax,normal_xmax,NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX)
- if (PERFORM_CUTHILL_MCKEE) then
- call sort_arrays_for_cuthill (ispecb2,xstore,ystore,zstore,ibelm_xmax,normal_xmax,&
- jacobian2D_xmax,NSPEC2DMAX_XMIN_XMAX,NGLLY,NGLLZ,nspec)
- endif
endif
! on boundary: ymin
@@ -227,10 +219,6 @@
call compute_jacobian_2D(myrank,ispecb3,xelm,yelm,zelm,dershape2D_y, &
jacobian2D_ymin,normal_ymin,NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)
- if (PERFORM_CUTHILL_MCKEE) then
- call sort_arrays_for_cuthill (ispecb3,xstore,ystore,zstore,ibelm_ymin,normal_ymin,&
- jacobian2D_ymin,NSPEC2DMAX_YMIN_YMAX,NGLLX,NGLLZ,nspec)
- endif
endif
! on boundary: ymax
@@ -272,10 +260,6 @@
call compute_jacobian_2D(myrank,ispecb4,xelm,yelm,zelm,dershape2D_y, &
jacobian2D_ymax,normal_ymax,NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)
- if (PERFORM_CUTHILL_MCKEE) then
- call sort_arrays_for_cuthill (ispecb4,xstore,ystore,zstore,ibelm_ymax,normal_ymax,&
- jacobian2D_ymax,NSPEC2DMAX_YMIN_YMAX,NGLLX,NGLLZ,nspec)
- endif
endif
! on boundary: bottom
@@ -316,10 +300,6 @@
call compute_jacobian_2D(myrank,ispecb5,xelm,yelm,zelm,dershape2D_bottom, &
jacobian2D_bottom,normal_bottom,NGLLX,NGLLY,NSPEC2D_BOTTOM)
- if (PERFORM_CUTHILL_MCKEE) then
- call sort_arrays_for_cuthill (ispecb5,xstore,ystore,zstore,ibelm_bottom,normal_bottom,&
- jacobian2D_bottom,NSPEC2D_BOTTOM,NGLLX,NGLLY,nspec)
- endif
endif
! on boundary: top
@@ -360,15 +340,10 @@
call compute_jacobian_2D(myrank,ispecb6,xelm,yelm,zelm,dershape2D_top, &
jacobian2D_top,normal_top,NGLLX,NGLLY,NSPEC2D_TOP)
- if (PERFORM_CUTHILL_MCKEE) then
- call sort_arrays_for_cuthill (ispecb6,xstore,ystore,zstore,ibelm_top,normal_top,&
- jacobian2D_top,NSPEC2D_TOP,NGLLX,NGLLY,nspec)
- endif
endif
enddo
-
! check theoretical value of elements at the bottom
if(ispecb5 /= NSPEC2D_BOTTOM) then
call exit_MPI(myrank,'ispecb5 should equal NSPEC2D_BOTTOM')
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/heterogen_mantle_model.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/heterogen_mantle_model.f90 2008-08-06 09:31:17 UTC (rev 12536)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/heterogen_mantle_model.f90 2008-08-06 10:28:46 UTC (rev 12537)
@@ -29,7 +29,7 @@
!
subroutine read_heterogen_mantle_model(HMM)
-
+
implicit none
include "constants.h"
@@ -51,7 +51,7 @@
form='formatted',recl=20,status='old',action='read')
j = N_R*N_THETA*N_PHI
-
+
do i = 1,j
read(10,rec=i,fmt='(F20.15)') HMM%rho_in(i)
end do
@@ -93,18 +93,18 @@
radius = radius*R_EARTH
r_inner = 3.500d6 !lower bound for heterogeneity zone
-! NOTE: r_outer NEEDS TO BE (just) SMALLER THAN R_EARTH!!!!!!!!
+! NOTE: r_outer NEEDS TO BE (just) SMALLER THAN R_EARTH!!!!!!!!
r_outer = R_EARTH-1.0d1 !6.300d6 !upper bound for heterogeneity zone (lower mantle: e.g. 4.500d6)
delta = 2.*R_EARTH/(real(N_R-1))
- delta2 = 2.*R_EARTH/(real(N_R-2))
- !delta2 = 2.*R_EARTH/(real(N_R))
+ delta2 = 2.*R_EARTH/(real(N_R-2))
+ !delta2 = 2.*R_EARTH/(real(N_R))
if ((radius >= r_inner) .and. (radius <= r_outer)) then
! convert spherical point to cartesian point, move origin to corner
x = R_EARTH + radius*sin(theta)*cos(phi)
y = R_EARTH + radius*sin(theta)*sin(phi)
- z = R_EARTH + radius*cos(theta)
+ z = R_EARTH + radius*cos(theta)
! determine which points to search for in heterogen.dat
! find x_low,y_low,z_low etc.
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D.f90 2008-08-06 09:31:17 UTC (rev 12536)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D.f90 2008-08-06 10:28:46 UTC (rev 12537)
@@ -1317,7 +1317,7 @@
call MPI_BCAST(AM_V%Qtau_s(1), 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ier)
call MPI_BCAST(AM_V%Qtau_s(2), 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ier)
call MPI_BCAST(AM_V%Qtau_s(3), 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ier)
- if(myrank == 0) call read_atten_model_3D_QRFSI12(QRFSI12_Q)
+ if(myrank == 0) call read_atten_model_3D_QRFSI12(QRFSI12_Q)
call MPI_BCAST(QRFSI12_Q%dqmu, NKQ*NSQ,MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ier)
call MPI_BCAST(QRFSI12_Q%spknt, NKQ,MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ier)
call MPI_BCAST(QRFSI12_Q%refdepth, NDEPTHS_REFQ,MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ier)
@@ -1419,7 +1419,6 @@
ANISOTROPIC_INNER_CORE,ISOTROPIC_3D_MANTLE,HETEROGEN_3D_MANTLE,CRUSTAL,ONE_CRUST, &
NPROC_XI,NPROC_ETA,NSPEC2D_XI_FACE, &
NSPEC2D_ETA_FACE,NSPEC1D_RADIAL_CORNER,NGLOB1D_RADIAL_CORNER, &
- max(NGLOB2DMAX_XMIN_XMAX(iregion_code),NGLOB2DMAX_YMIN_YMAX(iregion_code)), &
myrank,LOCAL_PATH,OCEANS,ibathy_topo, &
rotation_matrix,ANGULAR_WIDTH_XI_RAD,ANGULAR_WIDTH_ETA_RAD, &
ATTENUATION,ATTENUATION_3D,SAVE_MESH_FILES, &
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90 2008-08-06 09:31:17 UTC (rev 12536)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90 2008-08-06 10:28:46 UTC (rev 12537)
@@ -4399,7 +4399,7 @@
dxy = dxy + dble(epsilondev_crust_mantle(3,i,j,k,ispec_selected_source(irec)))*hlagrange
dxz = dxz + dble(epsilondev_crust_mantle(4,i,j,k,ispec_selected_source(irec)))*hlagrange
dyz = dyz + dble(epsilondev_crust_mantle(5,i,j,k,ispec_selected_source(irec)))*hlagrange
-
+
displ_s(:,i,j,k) = displ_crust_mantle(:,iglob)
enddo
@@ -4450,7 +4450,7 @@
xix_crust_mantle(:,:,:,ispec),xiy_crust_mantle(:,:,:,ispec),xiz_crust_mantle(:,:,:,ispec), &
etax_crust_mantle(:,:,:,ispec),etay_crust_mantle(:,:,:,ispec),etaz_crust_mantle(:,:,:,ispec), &
gammax_crust_mantle(:,:,:,ispec),gammay_crust_mantle(:,:,:,ispec),gammaz_crust_mantle(:,:,:,ispec))
-
+
stf = comp_source_time_function(dble(NSTEP-it)*DT-t0-t_cmt(irec),hdur_gaussian(irec))
stf_deltat = stf * deltat
moment_der(:,:,irec_local) = moment_der(:,:,irec_local) + eps_s(:,:) * stf_deltat
@@ -5154,7 +5154,7 @@
! save source derivatives for adjoint simulations
if (SIMULATION_TYPE == 2 .and. nrec_local > 0) then
- scale_mass = RHOAV * (R_EARTH**3)
+ scale_mass = RHOAV * (R_EARTH**3)
do irec_local = 1, nrec_local
! rotate and scale the location derivatives to correspond to dn,de,dz
@@ -5162,7 +5162,7 @@
! rotate scale the moment derivatives to correspond to M[n,e,z][n,e,z]
moment_der(:,:,irec_local) = matmul(matmul(nu_source(:,:,irec),moment_der(:,:,irec_local)),&
transpose(nu_source(:,:,irec))) * scale_t ** 3 / scale_mass
-
+
write(outputname,'(a,i5.5)') 'OUTPUT_FILES/src_frechet.',number_receiver_global(irec_local)
open(unit=27,file=trim(outputname),status='unknown')
!
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/write_movie_volume.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/write_movie_volume.f90 2008-08-06 09:31:17 UTC (rev 12536)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/write_movie_volume.f90 2008-08-06 10:28:46 UTC (rev 12537)
@@ -430,7 +430,7 @@
ipoints_3dmovie=ipoints_3dmovie+1
iglob = ibool_crust_mantle(i,j,k,ispec)
vector_local(:) = vector_scaled(:,iglob)
-
+
! rotate eps_loc to spherical coordinates
vector_local_new(:) = matmul(nu_3dmovie(:,:,ipoints_3dmovie), vector_local(:))
store_val3d_N(ipoints_3dmovie)=vector_local_new(1)
More information about the cig-commits
mailing list