[cig-commits] r22478 - in seismo/3D/SPECFEM3D_GLOBE: branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D trunk/src/auxiliaries trunk/src/meshfem3D trunk/src/shared trunk/src/specfem3D
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Mon Jul 1 05:46:50 PDT 2013
Author: dkomati1
Date: 2013-07-01 05:46:50 -0700 (Mon, 01 Jul 2013)
New Revision: 22478
Modified:
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/check_buffers_1D.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/check_buffers_corners_chunks.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/combine_vol_data.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/compute_coordinates_grid.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_chunk_buffers.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/model_topo_bathy.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_movie_volume.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/check_buffers_corners_chunks.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/check_buffers_faces_chunks.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_AVS_DX.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_vol_data.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/add_topography.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_chunk_buffers.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/save_arrays_solver.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/model_topo_bathy.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_movie_volume.f90
Log:
more obvious merges
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/check_buffers_1D.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/check_buffers_1D.f90 2013-07-01 11:55:54 UTC (rev 22477)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/check_buffers_1D.f90 2013-07-01 12:46:50 UTC (rev 22478)
@@ -110,6 +110,7 @@
integer, dimension(MAX_NUM_REGIONS,NB_SQUARE_CORNERS) :: NGLOB1D_RADIAL_CORNER
integer, dimension(NB_SQUARE_CORNERS) :: NGLOB1D_RADIAL_SPEC_THIS
integer, dimension(NB_SQUARE_CORNERS) :: NGLOB1D_RADIAL_SPEC_OTHER
+
! ************** PROGRAM STARTS HERE **************
print *
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/check_buffers_corners_chunks.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/check_buffers_corners_chunks.f90 2013-07-01 11:55:54 UTC (rev 22477)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/check_buffers_corners_chunks.f90 2013-07-01 12:46:50 UTC (rev 22478)
@@ -102,6 +102,7 @@
logical :: CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA
integer, dimension(NB_SQUARE_CORNERS,NB_CUT_CASE) :: DIFF_NSPEC1D_RADIAL
integer, dimension(NB_SQUARE_EDGES_ONEDIR,NB_CUT_CASE) :: DIFF_NSPEC2D_XI,DIFF_NSPEC2D_ETA
+
! ************** PROGRAM STARTS HERE **************
print *
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/combine_vol_data.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/combine_vol_data.f90 2013-07-01 11:55:54 UTC (rev 22477)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/combine_vol_data.f90 2013-07-01 12:46:50 UTC (rev 22478)
@@ -266,7 +266,6 @@
print *, trim(topo_file)
-
!average data on global points
ibool_count(:) = 0
ibool_dat(:) = 0.0
@@ -303,7 +302,6 @@
num_ibool(:) = 0
numpoin = 0
-
! write point file
do ispec=1,nspec(it)
! checks if element counts
@@ -621,7 +619,6 @@
r(i) = r_ocean+(r_0-r_ocean)*dble(i-634)/dble(6)
enddo
-
! use PREM to get the density profile for ellipticity (fine for other 1D reference models)
do i=1,NR
call prem_density(r(i),rho(i),ONE_CRUST,RICB,RCMB,RTOPDDOUBLEPRIME, &
@@ -762,8 +759,6 @@
integer i,j,n,kdis(28)
integer ndis,nir1
-
-
data kdis/163,323,336,517,530,540,565,590,609,619,626,633,16*0/
ndis = 12
@@ -792,12 +787,12 @@
implicit none
-! Argument variables
+! argument variables
integer kdis(28),n,ndis
double precision r(n),s1(n),s2(n),s3(n)
double precision y(n),yprime(n)
-! Local variables
+! local variables
integer i,j,j1,j2
integer k,nd,ndp
double precision a0,b0,b1
@@ -1013,7 +1008,6 @@
end subroutine spline_evaluation
-
!
! ------------------------------------------------------------------------------------------------
!
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/compute_coordinates_grid.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/compute_coordinates_grid.f90 2013-07-01 11:55:54 UTC (rev 22477)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/compute_coordinates_grid.f90 2013-07-01 12:46:50 UTC (rev 22478)
@@ -312,13 +312,13 @@
fact_y = 2.d0*ratio_y-1.d0
fact_z = 2.d0*ratio_z-1.d0
- xi = PI_OVER_TWO*fact_x;
- eta = PI_OVER_TWO*fact_y;
- gamma = PI_OVER_TWO*fact_z;
+ xi = PI_OVER_TWO*fact_x
+ eta = PI_OVER_TWO*fact_y
+ gamma = PI_OVER_TWO*fact_z
- xgrid_central_cube = radius_cube * fact_x * (1 + (cos(eta)+cos(gamma))*CENTRAL_CUBE_INFLATE_FACTOR / PI);
- ygrid_central_cube = radius_cube * fact_y * (1 + (cos(xi)+cos(gamma))*CENTRAL_CUBE_INFLATE_FACTOR / PI);
- zgrid_central_cube = radius_cube * fact_z * (1 + (cos(xi)+cos(eta))*CENTRAL_CUBE_INFLATE_FACTOR / PI);
+ xgrid_central_cube = radius_cube * fact_x * (1 + (cos(eta)+cos(gamma))*CENTRAL_CUBE_INFLATE_FACTOR / PI)
+ ygrid_central_cube = radius_cube * fact_y * (1 + (cos(xi)+cos(gamma))*CENTRAL_CUBE_INFLATE_FACTOR / PI)
+ zgrid_central_cube = radius_cube * fact_z * (1 + (cos(xi)+cos(eta))*CENTRAL_CUBE_INFLATE_FACTOR / PI)
end subroutine compute_coord_central_cube
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_chunk_buffers.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_chunk_buffers.f90 2013-07-01 11:55:54 UTC (rev 22477)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/create_chunk_buffers.f90 2013-07-01 12:46:50 UTC (rev 22478)
@@ -78,14 +78,18 @@
character(len=150) :: OUTPUT_FILES,ERR_MSG
! mask for ibool to mark points already found
logical, dimension(:), allocatable :: mask_ibool
+
! array to store points selected for the chunk face buffer
integer :: NGLOB2DMAX_XY
integer, dimension(:), allocatable :: ibool_selected
+
double precision, dimension(:), allocatable :: xstore_selected,ystore_selected,zstore_selected
+
! arrays for sorting routine
integer, dimension(:), allocatable :: ind,ninseg,iglob,locval,iwork
logical, dimension(:), allocatable :: ifseg
double precision, dimension(:), allocatable :: work
+
! pairs generated theoretically
! four sides for each of the three types of messages
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/model_topo_bathy.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/model_topo_bathy.f90 2013-07-01 11:55:54 UTC (rev 22477)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/model_topo_bathy.f90 2013-07-01 12:46:50 UTC (rev 22478)
@@ -33,7 +33,6 @@
! by default (constants.h), it uses a smoothed ETOPO 4 dataset
!--------------------------------------------------------------------------------------------------
-
subroutine model_topo_bathy_broadcast(myrank,ibathy_topo,LOCAL_PATH)
! standard routine to setup model
@@ -156,7 +155,6 @@
end subroutine read_topo_bathy_file
-
!
!-------------------------------------------------------------------------------------------------
!
@@ -268,10 +266,10 @@
integer, dimension(NX_BATHY,NY_BATHY),intent(in) :: ibathy_topo
! local parameters
- integer:: iadd1,iel1
- double precision:: samples_per_degree_topo
- double precision:: xlo
- double precision:: lon_corner,lat_corner,ratio_lon,ratio_lat
+ integer :: iadd1,iel1
+ double precision :: samples_per_degree_topo
+ double precision :: xlo
+ double precision :: lon_corner,lat_corner,ratio_lon,ratio_lat
! initializes elevation
value = ZERO
@@ -327,3 +325,7 @@
endif
end subroutine get_topo_bathy
+
+! -------------------------------------------
+
+
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_movie_volume.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_movie_volume.f90 2013-07-01 11:55:54 UTC (rev 22477)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_movie_volume.f90 2013-07-01 12:46:50 UTC (rev 22478)
@@ -25,12 +25,11 @@
!
!=====================================================================
-
- subroutine movie_volume_count_points()
-
! this subroutine counts the number of points and elements within the movie volume
! in this processor slice, and returns arrays that keep track of them, both in global and local indexing schemes
+ subroutine movie_volume_count_points()
+
use specfem_par
use specfem_par_crustmantle
use specfem_par_movie
@@ -282,7 +281,6 @@
!-------------------------------------------------------------------------------------------------
!
-
subroutine write_movie_volume_strains(myrank,npoints_3dmovie, &
LOCAL_TMP_PATH,MOVIE_VOLUME_TYPE,MOVIE_COARSE, &
it,eps_trace_over_3_crust_mantle, &
@@ -291,7 +289,6 @@
muvstore_crust_mantle_3dmovie, &
mask_3dmovie,nu_3dmovie)
-
! outputs strains: MOVIE_VOLUME_TYPE == 1 / 2 / 3
implicit none
@@ -504,7 +501,8 @@
vector_scaled = vector_crust_mantle*scalingval
endif
- ipoints_3dmovie=0
+ ipoints_3dmovie = 0
+
do ispec=1,NSPEC_CRUST_MANTLE
do k=1,NGLLZ,iNIT
do j=1,NGLLY,iNIT
@@ -542,8 +540,7 @@
write(27) store_val3d_Z(1:npoints_3dmovie)
close(27)
- deallocate(store_val3d_N,store_val3d_E,store_val3d_Z, &
- vector_scaled)
+ deallocate(store_val3d_N,store_val3d_E,store_val3d_Z,vector_scaled)
end subroutine write_movie_volume_vector
@@ -596,8 +593,8 @@
! output parameters
logical,parameter :: MOVIE_OUTPUT_DIV = .true. ! divergence
- logical,parameter :: MOVIE_OUTPUT_CURL = .true. ! curl
- logical,parameter :: MOVIE_OUTPUT_CURLNORM = .true. ! frobenius norm of curl
+ logical,parameter :: MOVIE_OUTPUT_CURL = .true. ! curl
+ logical,parameter :: MOVIE_OUTPUT_CURLNORM = .true. ! Frobenius norm of curl
! outputs divergence
if( MOVIE_OUTPUT_DIV ) then
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/check_buffers_corners_chunks.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/check_buffers_corners_chunks.f90 2013-07-01 11:55:54 UTC (rev 22477)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/check_buffers_corners_chunks.f90 2013-07-01 12:46:50 UTC (rev 22478)
@@ -131,7 +131,9 @@
MOVIE_VOLUME,MOVIE_COARSE,ATTENUATION_3D,RECEIVERS_CAN_BE_BURIED, &
PRINT_SOURCE_TIME_FUNCTION,SAVE_MESH_FILES, &
ATTENUATION,REFERENCE_1D_MODEL,THREE_D_MODEL,ABSORBING_CONDITIONS, &
- INCLUDE_CENTRAL_CUBE,INFLATE_CENTRAL_CUBE,LOCAL_PATH,MODEL,SIMULATION_TYPE,SAVE_FORWARD, &
+ INCLUDE_CENTRAL_CUBE,INFLATE_CENTRAL_CUBE, &
+ LOCAL_PATH,MODEL, &
+ SIMULATION_TYPE,SAVE_FORWARD, &
NPROC,NPROCTOT,NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
NSPEC, &
NSPEC2D_XI, &
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/check_buffers_faces_chunks.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/check_buffers_faces_chunks.f90 2013-07-01 11:55:54 UTC (rev 22477)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/check_buffers_faces_chunks.f90 2013-07-01 12:46:50 UTC (rev 22478)
@@ -135,7 +135,9 @@
MOVIE_VOLUME,MOVIE_COARSE,ATTENUATION_3D,RECEIVERS_CAN_BE_BURIED, &
PRINT_SOURCE_TIME_FUNCTION,SAVE_MESH_FILES, &
ATTENUATION,REFERENCE_1D_MODEL,THREE_D_MODEL,ABSORBING_CONDITIONS, &
- INCLUDE_CENTRAL_CUBE,INFLATE_CENTRAL_CUBE,LOCAL_PATH,MODEL,SIMULATION_TYPE,SAVE_FORWARD, &
+ INCLUDE_CENTRAL_CUBE,INFLATE_CENTRAL_CUBE, &
+ LOCAL_PATH,MODEL, &
+ SIMULATION_TYPE,SAVE_FORWARD, &
NPROC,NPROCTOT,NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
NSPEC,NSPEC2D_XI,NSPEC2D_ETA, &
NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_AVS_DX.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_AVS_DX.f90 2013-07-01 11:55:54 UTC (rev 22477)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_AVS_DX.f90 2013-07-01 12:46:50 UTC (rev 22478)
@@ -213,7 +213,9 @@
MOVIE_VOLUME,MOVIE_COARSE,ATTENUATION_3D,RECEIVERS_CAN_BE_BURIED, &
PRINT_SOURCE_TIME_FUNCTION,SAVE_MESH_FILES, &
ATTENUATION,REFERENCE_1D_MODEL,THREE_D_MODEL,ABSORBING_CONDITIONS, &
- INCLUDE_CENTRAL_CUBE,INFLATE_CENTRAL_CUBE,LOCAL_PATH,MODEL,SIMULATION_TYPE,SAVE_FORWARD, &
+ INCLUDE_CENTRAL_CUBE,INFLATE_CENTRAL_CUBE, &
+ LOCAL_PATH,MODEL, &
+ SIMULATION_TYPE,SAVE_FORWARD, &
NPROC,NPROCTOT,NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
NSPEC_COMP,NSPEC2D_XI,NSPEC2D_ETA,NSPEC2DMAX_XMIN_XMAX,&
NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_vol_data.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_vol_data.f90 2013-07-01 11:55:54 UTC (rev 22477)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_vol_data.f90 2013-07-01 12:46:50 UTC (rev 22478)
@@ -192,7 +192,6 @@
print*,'file:',trim(dimension_file)
stop 'Error opening file'
endif
-
read(27,*) nspec(it)
read(27,*) nglob(it)
close(27)
@@ -266,7 +265,6 @@
print *, trim(topo_file)
-
!average data on global points
ibool_count(:) = 0
ibool_dat(:) = 0.0
@@ -303,7 +301,6 @@
num_ibool(:) = 0
numpoin = 0
-
! write point file
do ispec=1,nspec(it)
! checks if element counts
@@ -621,7 +618,6 @@
r(i) = r_ocean+(r_0-r_ocean)*dble(i-634)/dble(6)
enddo
-
! use PREM to get the density profile for ellipticity (fine for other 1D reference models)
do i=1,NR
call prem_density(r(i),rho(i),ONE_CRUST,RICB,RCMB,RTOPDDOUBLEPRIME, &
@@ -762,8 +758,6 @@
integer i,j,n,kdis(28)
integer ndis,nir1
-
-
data kdis/163,323,336,517,530,540,565,590,609,619,626,633,16*0/
ndis = 12
@@ -792,12 +786,12 @@
implicit none
-! Argument variables
+! argument variables
integer kdis(28),n,ndis
double precision r(n),s1(n),s2(n),s3(n)
double precision y(n),yprime(n)
-! Local variables
+! local variables
integer i,j,j1,j2
integer k,nd,ndp
double precision a0,b0,b1
@@ -1013,7 +1007,6 @@
end subroutine spline_evaluation
-
!
! ------------------------------------------------------------------------------------------------
!
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/add_topography.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/add_topography.f90 2013-07-01 11:55:54 UTC (rev 22477)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/add_topography.f90 2013-07-01 12:46:50 UTC (rev 22478)
@@ -31,50 +31,50 @@
include "constants.h"
- double precision xelm(NGNOD)
- double precision yelm(NGNOD)
- double precision zelm(NGNOD)
+ integer :: myrank
- integer myrank
+ double precision,dimension(NGNOD) :: xelm,yelm,zelm
-! use integer array to store values
+ ! use integer array to store values
integer, dimension(NX_BATHY,NY_BATHY) :: ibathy_topo
- integer ia
+ ! local parameters
+ double precision :: lat,lon,elevation
+ double precision :: r,theta,phi,colat
+ double precision :: gamma
+ integer :: ia
- double precision lat,lon,elevation,R220
- double precision r,theta,phi,colat
- double precision gamma
+ double precision R220
-! we loop on all the points of the element
+ ! we loop on all the points of the element
do ia = 1,NGNOD
-! convert to r theta phi
-! slightly move points to avoid roundoff problem when exactly on the polar axis
+ ! gets elevation of point
+ ! convert to r theta phi
+ ! slightly move points to avoid roundoff problem when exactly on the polar axis
call xyz_2_rthetaphi_dble(xelm(ia),yelm(ia),zelm(ia),r,theta,phi)
theta = theta + 0.0000001d0
phi = phi + 0.0000001d0
call reduce(theta,phi)
-! convert the geocentric colatitude to a geographic colatitude
- colat = PI/2.0d0 - datan(1.006760466d0*dcos(theta)/dmax1(TINYVAL,dsin(theta)))
+ ! convert the geocentric colatitude to a geographic colatitude
+ colat = PI_OVER_TWO - datan(1.006760466d0*dcos(theta)/dmax1(TINYVAL,dsin(theta)))
-! get geographic latitude and longitude in degrees
- lat = 90.0d0 - colat*180.0d0/PI
- lon = phi*180.0d0/PI
- elevation = 0.d0
+ ! get geographic latitude and longitude in degrees
+ lat = (PI_OVER_TWO - colat) * RADIANS_TO_DEGREES
+ lon = phi * RADIANS_TO_DEGREES
-! compute elevation at current point
+ ! compute elevation at current point
call get_topo_bathy(lat,lon,elevation,ibathy_topo)
-! non-dimensionalize the elevation, which is in meters
+ ! non-dimensionalize the elevation, which is in meters
elevation = elevation / R_EARTH
-! stretching topography between d220 and the surface
+ ! stretching topography between d220 and the surface
gamma = (r - R220/R_EARTH) / (R_UNIT_SPHERE - R220/R_EARTH)
-! add elevation to all the points of that element
-! also make sure gamma makes sense
+ ! add elevation to all the points of that element
+ ! also make sure gamma makes sense
if(gamma < -0.02 .or. gamma > 1.02) call exit_MPI(myrank,'incorrect value of gamma for topography')
xelm(ia) = xelm(ia)*(ONE + gamma * elevation / r)
@@ -110,7 +110,9 @@
! input parameters
integer:: myrank
integer:: ispec,nspec
+
double precision,dimension(NGLLX,NGLLY,NGLLZ,nspec):: xstore,ystore,zstore
+
integer, dimension(NX_BATHY,NY_BATHY) :: ibathy_topo
double precision:: R220
@@ -133,12 +135,11 @@
! convert the geocentric colatitude to a geographic colatitude
- colat = PI/2.0d0 - datan(1.006760466d0*dcos(theta)/dmax1(TINYVAL,dsin(theta)))
+ colat = PI_OVER_TWO - datan(1.006760466d0*dcos(theta)/dmax1(TINYVAL,dsin(theta)))
! get geographic latitude and longitude in degrees
- lat = 90.0d0 - colat*180.0d0/PI
- lon = phi*180.0d0/PI
- elevation = 0.d0
+ lat = (PI_OVER_TWO - colat) * RADIANS_TO_DEGREES
+ lon = phi * RADIANS_TO_DEGREES
! compute elevation at current point
call get_topo_bathy(lat,lon,elevation,ibathy_topo)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_chunk_buffers.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_chunk_buffers.f90 2013-07-01 11:55:54 UTC (rev 22477)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_chunk_buffers.f90 2013-07-01 12:46:50 UTC (rev 22478)
@@ -28,8 +28,7 @@
! subroutine to create MPI buffers to assemble between chunks
subroutine create_chunk_buffers(iregion_code,nspec,ibool,idoubling, &
- xstore,ystore,zstore, &
- nglob_ori, &
+ xstore,ystore,zstore,nglob_ori, &
NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX, &
NPROC_XI,NPROC_ETA,NPROC,NPROCTOT, &
NGLOB1D_RADIAL_CORNER,NGLOB1D_RADIAL_MAX, &
@@ -54,10 +53,8 @@
integer nspec
integer myrank,NCHUNKS
-! arrays with the mesh
- double precision xstore(NGLLX,NGLLY,NGLLZ,nspec)
- double precision ystore(NGLLX,NGLLY,NGLLZ,nspec)
- double precision zstore(NGLLX,NGLLY,NGLLZ,nspec)
+ ! arrays with the mesh
+ double precision,dimension(NGLLX,NGLLY,NGLLZ,nspec) :: xstore,ystore,zstore
character(len=150) OUTPUT_FILES,LOCAL_PATH,ERR_MSG
@@ -66,22 +63,23 @@
integer idoubling(nspec)
-! mask for ibool to mark points already found
+ ! mask for ibool to mark points already found
logical, dimension(:), allocatable :: mask_ibool
-! array to store points selected for the chunk face buffer
- integer NGLOB2DMAX_XY
+ ! array to store points selected for the chunk face buffer
+ integer :: NGLOB2DMAX_XY
integer, dimension(:), allocatable :: ibool_selected
double precision, dimension(:), allocatable :: xstore_selected,ystore_selected,zstore_selected
-! arrays for sorting routine
+ ! arrays for sorting routine
integer, dimension(:), allocatable :: ind,ninseg,iglob,locval,iwork
logical, dimension(:), allocatable :: ifseg
double precision, dimension(:), allocatable :: work
-! pairs generated theoretically
-! four sides for each of the three types of messages
+ ! pairs generated theoretically
+
+ ! four sides for each of the three types of messages
integer, dimension(:), allocatable :: iproc_sender,iproc_receiver,npoin2D_send,npoin2D_receive
! 1D buffers to remove points belonging to corners
@@ -96,14 +94,14 @@
double precision xdummy,ydummy,zdummy
integer ipoin1D
-! arrays to assemble the corners (3 processors for each corner)
+ ! arrays to assemble the corners (3 processors for each corner)
integer, dimension(:,:), allocatable :: iprocscorners,itypecorner
- integer ichunk_send,iproc_xi_send,iproc_eta_send
- integer ichunk_receive,iproc_xi_receive,iproc_eta_receive
- integer iproc_loop,iproc_xi_loop,iproc_eta_loop
- integer iproc_xi_loop_inv,iproc_eta_loop_inv
- integer imember_corner
+ integer :: ichunk_send,iproc_xi_send,iproc_eta_send
+ integer :: ichunk_receive,iproc_xi_receive,iproc_eta_receive
+ integer :: iproc_loop,iproc_xi_loop,iproc_eta_loop
+ integer :: iproc_xi_loop_inv,iproc_eta_loop_inv
+ integer :: imember_corner
integer iregion_code
@@ -119,7 +117,7 @@
integer i,j,k,ispec,ispec2D,ipoin2D,ier
-! current message number
+ ! current message number
integer imsg
! names of the data files for all the processors in MPI
@@ -174,7 +172,7 @@
ichunk_receive = 0
ichunk_send = 0
-! number of corners and faces shared between chunks and number of message types
+ ! number of corners and faces shared between chunks and number of message types
if(NCHUNKS == 1 .or. NCHUNKS == 2) then
NCORNERSCHUNKS = 1
NUM_FACES = 1
@@ -191,10 +189,10 @@
call exit_MPI(myrank,'number of chunks must be either 1, 2, 3 or 6')
endif
-! if more than one chunk then same number of processors in each direction
+ ! if more than one chunk then same number of processors in each direction
NPROC_ONE_DIRECTION = NPROC_XI
-! total number of messages corresponding to these common faces
+ ! total number of messages corresponding to these common faces
NUMMSGS_FACES = NPROC_ONE_DIRECTION*NUM_FACES*NUM_MSG_TYPES
! check that there is more than one chunk, otherwise nothing to do
@@ -209,7 +207,7 @@
allocate(npoin2D_send(NUMMSGS_FACES))
allocate(npoin2D_receive(NUMMSGS_FACES))
-! allocate array for corners
+ ! allocate array for corners
allocate(iprocscorners(3,NCORNERSCHUNKS))
allocate(itypecorner(3,NCORNERSCHUNKS))
@@ -263,28 +261,28 @@
do iside = 1,NUM_FACES
do iproc_loop = 0,NPROC_ONE_DIRECTION-1
-! create a new message
-! we know there can be no deadlock with this scheme
-! because the three types of messages are independent
+ ! create a new message
+ ! we know there can be no deadlock with this scheme
+ ! because the three types of messages are independent
imsg = imsg + 1
-! check that current message number is correct
+ ! check that current message number is correct
if(imsg > NUMMSGS_FACES) call exit_MPI(myrank,'incorrect message number')
if(myrank == 0) write(IMAIN,*) 'Generating message ',imsg,' for faces out of ',NUMMSGS_FACES
-! we know there is the same number of slices in both directions
+ ! we know there is the same number of slices in both directions
iproc_xi_loop = iproc_loop
iproc_eta_loop = iproc_loop
-! take care of local frame inversions between chunks
+ ! take care of local frame inversions between chunks
iproc_xi_loop_inv = NPROC_ONE_DIRECTION - iproc_loop - 1
iproc_eta_loop_inv = NPROC_ONE_DIRECTION - iproc_loop - 1
-! define the 12 different messages
+ ! define the 12 different messages
-! message type M1
+ ! message type M1
if(imsg_type == 1) then
if(iside == 1) then
@@ -333,7 +331,7 @@
endif
-! message type M2
+ ! message type M2
if(imsg_type == 2) then
if(iside == 1) then
@@ -382,7 +380,7 @@
endif
-! message type M3
+ ! message type M3
if(imsg_type == 3) then
if(iside == 1) then
@@ -442,7 +440,7 @@
! save message type and pair of processors in list of messages
if(myrank == 0) write(IOUT,*) imsg_type,iproc_sender(imsg),iproc_receiver(imsg)
-! loop on sender/receiver (1=sender 2=receiver)
+ ! loop on sender/receiver (1=sender 2=receiver)
do imode_comm=1,2
if(imode_comm == 1) then
@@ -457,7 +455,7 @@
call exit_MPI(myrank,'incorrect communication mode')
endif
-! only do this if current processor is the right one for MPI version
+ ! only do this if current processor is the right one for MPI version
if(iproc == myrank) then
! create the name of the database for each slice
@@ -466,18 +464,19 @@
! open file for 2D buffer
open(unit=IOUT_BUFFERS,file=prname(1:len_trim(prname))//filename_out,status='unknown')
-! determine chunk number and local slice coordinates using addressing
+ ! determine chunk number and local slice coordinates using addressing
ichunk = ichunk_slice(iproc)
iproc_xi = iproc_xi_slice(iproc)
iproc_eta = iproc_eta_slice(iproc)
-! problem if not on edges
+ ! problem if not on edges
if(iproc_xi /= 0 .and. iproc_xi /= NPROC_XI-1 .and. &
- iproc_eta /= 0 .and. iproc_eta /= NPROC_ETA-1) call exit_MPI(myrank,'slice not on any edge')
+ iproc_eta /= 0 .and. iproc_eta /= NPROC_ETA-1) &
+ call exit_MPI(myrank,'slice not on any edge')
nglob=nglob_ori
-! check that iboolmax=nglob
+ ! check that iboolmax=nglob
if(minval(ibool(:,:,:,1:nspec)) /= 1 .or. maxval(ibool(:,:,:,1:nspec)) /= nglob) &
call exit_MPI(myrank,ERR_MSG)
@@ -524,17 +523,16 @@
enddo
close(IIN)
-! erase logical mask
+ ! erase logical mask
mask_ibool(:) = .false.
npoin2D = 0
-! create all the points on each face (no duplicates, but not sorted)
+ ! create all the points on each face (no duplicates, but not sorted)
-! xmin
+ ! xmin
if(iedge == XI_MIN) then
-
-! mark corner points to remove them if needed
+ ! mark corner points to remove them if needed
if(iproc_eta == 0) then
do ipoin1D = 1,NGLOB1D_RADIAL_CORNER(iregion_code,1)
mask_ibool(ibool1D_leftxi_lefteta(ipoin1D)) = .true.
@@ -550,7 +548,7 @@
do ispec2D=1,nspec2D_xmin
ispec=ibelm_xmin(ispec2D)
-! remove central cube for chunk buffers
+ ! remove central cube for chunk buffers
if(idoubling(ispec) == IFLAG_MIDDLE_CENTRAL_CUBE .or. &
idoubling(ispec) == IFLAG_BOTTOM_CENTRAL_CUBE .or. &
idoubling(ispec) == IFLAG_TOP_CENTRAL_CUBE .or. &
@@ -560,10 +558,12 @@
do k=1,NGLLZ
do j=1,NGLLY
if(.not. mask_ibool(ibool(i,j,k,ispec))) then
-! mask and store points found
+ ! mask and store points found
mask_ibool(ibool(i,j,k,ispec)) = .true.
npoin2D = npoin2D + 1
- if(npoin2D > NGLOB2DMAX_XMIN_XMAX) call exit_MPI(myrank,'incorrect 2D point number in xmin')
+ if(npoin2D > NGLOB2DMAX_XMIN_XMAX) &
+ call exit_MPI(myrank,'incorrect 2D point number in xmin')
+
ibool_selected(npoin2D) = ibool(i,j,k,ispec)
xstore_selected(npoin2D) = xstore(i,j,k,ispec)
@@ -574,11 +574,9 @@
enddo
enddo
-! xmax
+ ! xmax
else if(iedge == XI_MAX) then
-
-! mark corner points to remove them if needed
-
+ ! mark corner points to remove them if needed
if(iproc_eta == 0) then
do ipoin1D = 1,NGLOB1D_RADIAL_CORNER(iregion_code,2)
mask_ibool(ibool1D_rightxi_lefteta(ipoin1D)) = .true.
@@ -594,7 +592,7 @@
do ispec2D=1,nspec2D_xmax
ispec=ibelm_xmax(ispec2D)
-! remove central cube for chunk buffers
+ ! remove central cube for chunk buffers
if(idoubling(ispec) == IFLAG_MIDDLE_CENTRAL_CUBE .or. &
idoubling(ispec) == IFLAG_BOTTOM_CENTRAL_CUBE .or. &
idoubling(ispec) == IFLAG_TOP_CENTRAL_CUBE .or. &
@@ -604,10 +602,12 @@
do k=1,NGLLZ
do j=1,NGLLY
if(.not. mask_ibool(ibool(i,j,k,ispec))) then
-! mask and store points found
+ ! mask and store points found
mask_ibool(ibool(i,j,k,ispec)) = .true.
npoin2D = npoin2D + 1
- if(npoin2D > NGLOB2DMAX_XMIN_XMAX) call exit_MPI(myrank,'incorrect 2D point number in xmax')
+ if(npoin2D > NGLOB2DMAX_XMIN_XMAX) &
+ call exit_MPI(myrank,'incorrect 2D point number in xmax')
+
ibool_selected(npoin2D) = ibool(i,j,k,ispec)
xstore_selected(npoin2D) = xstore(i,j,k,ispec)
@@ -618,11 +618,9 @@
enddo
enddo
-! ymin
+ ! ymin
else if(iedge == ETA_MIN) then
-
-! mark corner points to remove them if needed
-
+ ! mark corner points to remove them if needed
if(iproc_xi == 0) then
do ipoin1D = 1,NGLOB1D_RADIAL_CORNER(iregion_code,1)
mask_ibool(ibool1D_leftxi_lefteta(ipoin1D)) = .true.
@@ -638,7 +636,7 @@
do ispec2D=1,nspec2D_ymin
ispec=ibelm_ymin(ispec2D)
-! remove central cube for chunk buffers
+ ! remove central cube for chunk buffers
if(idoubling(ispec) == IFLAG_MIDDLE_CENTRAL_CUBE .or. &
idoubling(ispec) == IFLAG_BOTTOM_CENTRAL_CUBE .or. &
idoubling(ispec) == IFLAG_TOP_CENTRAL_CUBE .or. &
@@ -648,10 +646,12 @@
do k=1,NGLLZ
do i=1,NGLLX
if(.not. mask_ibool(ibool(i,j,k,ispec))) then
-! mask and store points found
+ ! mask and store points found
mask_ibool(ibool(i,j,k,ispec)) = .true.
npoin2D = npoin2D + 1
- if(npoin2D > NGLOB2DMAX_YMIN_YMAX) call exit_MPI(myrank,'incorrect 2D point number in ymin')
+ if(npoin2D > NGLOB2DMAX_YMIN_YMAX) &
+ call exit_MPI(myrank,'incorrect 2D point number in ymin')
+
ibool_selected(npoin2D) = ibool(i,j,k,ispec)
xstore_selected(npoin2D) = xstore(i,j,k,ispec)
@@ -662,11 +662,9 @@
enddo
enddo
-! ymax
+ ! ymax
else if(iedge == ETA_MAX) then
-
-! mark corner points to remove them if needed
-
+ ! mark corner points to remove them if needed
if(iproc_xi == 0) then
do ipoin1D = 1,NGLOB1D_RADIAL_CORNER(iregion_code,4)
mask_ibool(ibool1D_leftxi_righteta(ipoin1D)) = .true.
@@ -682,7 +680,7 @@
do ispec2D=1,nspec2D_ymax
ispec=ibelm_ymax(ispec2D)
-! remove central cube for chunk buffers
+ ! remove central cube for chunk buffers
if(idoubling(ispec) == IFLAG_MIDDLE_CENTRAL_CUBE .or. &
idoubling(ispec) == IFLAG_BOTTOM_CENTRAL_CUBE .or. &
idoubling(ispec) == IFLAG_TOP_CENTRAL_CUBE .or. &
@@ -692,10 +690,12 @@
do k=1,NGLLZ
do i=1,NGLLX
if(.not. mask_ibool(ibool(i,j,k,ispec))) then
-! mask and store points found
+ ! mask and store points found
mask_ibool(ibool(i,j,k,ispec)) = .true.
npoin2D = npoin2D + 1
- if(npoin2D > NGLOB2DMAX_YMIN_YMAX) call exit_MPI(myrank,'incorrect 2D point number in ymax')
+ if(npoin2D > NGLOB2DMAX_YMIN_YMAX) &
+ call exit_MPI(myrank,'incorrect 2D point number in ymax')
+
ibool_selected(npoin2D) = ibool(i,j,k,ispec)
xstore_selected(npoin2D) = xstore(i,j,k,ispec)
@@ -711,16 +711,17 @@
call exit_MPI(myrank,'incorrect edge code')
endif
-! sort buffer obtained to be conforming with neighbor in other chunk
-! sort on x, y and z, the other arrays will be swapped as well
+ ! sort buffer obtained to be conforming with neighbor in other chunk
+ ! sort on x, y and z, the other arrays will be swapped as well
call sort_array_coordinates(npoin2D,xstore_selected,ystore_selected,zstore_selected, &
ibool_selected,iglob,locval,ifseg,nglob,ind,ninseg,iwork,work)
-! check that no duplicate has been detected
+ ! check that no duplicate has been detected
if(nglob /= npoin2D) call exit_MPI(myrank,'duplicates detected in buffer')
-! write list of selected points to output buffer
+ ! write list of selected points to output buffer
+
write(IOUT_BUFFERS,*) npoin2D
do ipoin2D = 1,npoin2D
write(IOUT_BUFFERS,*) ibool_selected(ipoin2D), &
@@ -729,56 +730,57 @@
close(IOUT_BUFFERS)
-! store result to compare number of points for sender and for receiver
+ ! store result to compare number of points for sender and for receiver
if(imode_comm == 1) then
npoin2D_send(imsg) = npoin2D
else
npoin2D_receive(imsg) = npoin2D
endif
-! end of section done only if right processor for MPI
+ ! end of section done only if right processor for MPI
endif
-! end of loop on sender/receiver
+ ! end of loop on sender/receiver
enddo
-! end of loops on all the messages
+ ! end of loops on all the messages
enddo
enddo
enddo
if(myrank == 0) close(IOUT)
-! check that total number of messages is correct
+ ! check that total number of messages is correct
if(imsg /= NUMMSGS_FACES) call exit_MPI(myrank,'incorrect total number of messages')
!
!---- check that number of points detected is the same for sender and receiver
!
-! synchronize all the processes to make sure all the buffers are ready
+ ! synchronize all the processes to make sure all the buffers are ready
call MPI_BARRIER(MPI_COMM_WORLD,ier)
-! gather information about all the messages on all processes
+ ! gather information about all the messages on all processes
do imsg = 1,NUMMSGS_FACES
-
-! gather number of points for sender
+ ! gather number of points for sender
npoin2D_send_local = npoin2D_send(imsg)
call MPI_BCAST(npoin2D_send_local,1,MPI_INTEGER,iproc_sender(imsg),MPI_COMM_WORLD,ier)
if(myrank /= iproc_sender(imsg)) npoin2D_send(imsg) = npoin2D_send_local
-! gather number of points for receiver
+ ! gather number of points for receiver
npoin2D_receive_local = npoin2D_receive(imsg)
call MPI_BCAST(npoin2D_receive_local,1,MPI_INTEGER,iproc_receiver(imsg),MPI_COMM_WORLD,ier)
if(myrank /= iproc_receiver(imsg)) npoin2D_receive(imsg) = npoin2D_receive_local
enddo
-! check the number of points
+ ! check the number of points
do imsg = 1,NUMMSGS_FACES
if(npoin2D_send(imsg) /= npoin2D_receive(imsg)) &
call exit_MPI(myrank,'incorrect number of points for sender/receiver pair detected')
enddo
+
+ ! user output
if(myrank == 0) then
write(IMAIN,*)
write(IMAIN,*) 'all the messages for chunk faces have the right size'
@@ -789,14 +791,14 @@
!---- generate the 8 message patterns sharing a corner of valence 3
!
-! to avoid problem at compile time, use bigger array with fixed dimension
+ ! to avoid problem at compile time, use bigger array with fixed dimension
addressing_big(:,:,:) = 0
addressing_big(1:NCHUNKS,:,:) = addressing(1:NCHUNKS,:,:)
ichunk = 1
iprocscorners(1,ichunk) = addressing_big(CHUNK_AB,0,NPROC_ETA-1)
iprocscorners(2,ichunk) = addressing_big(CHUNK_AC,NPROC_XI-1,NPROC_ETA-1)
-! this line is ok even for NCHUNKS = 2
+ ! this line is ok even for NCHUNKS = 2
iprocscorners(3,ichunk) = addressing_big(CHUNK_BC,NPROC_XI-1,0)
itypecorner(1,ichunk) = ILOWERUPPER
@@ -808,83 +810,83 @@
!! DK DK UGLY formally this is incorrect and should be changed in the future
!! DK DK UGLY in practice this trick works fine
-! this only if more than 3 chunks
+ ! this only if more than 3 chunks
if(NCHUNKS > 3) then
- ichunk = 2
- iprocscorners(1,ichunk) = addressing_big(CHUNK_AB,NPROC_XI-1,0)
- iprocscorners(2,ichunk) = addressing_big(CHUNK_AC_ANTIPODE,0,0)
- iprocscorners(3,ichunk) = addressing_big(CHUNK_BC_ANTIPODE,NPROC_XI-1,0)
+ ichunk = 2
+ iprocscorners(1,ichunk) = addressing_big(CHUNK_AB,NPROC_XI-1,0)
+ iprocscorners(2,ichunk) = addressing_big(CHUNK_AC_ANTIPODE,0,0)
+ iprocscorners(3,ichunk) = addressing_big(CHUNK_BC_ANTIPODE,NPROC_XI-1,0)
- itypecorner(1,ichunk) = IUPPERLOWER
- itypecorner(2,ichunk) = ILOWERLOWER
- itypecorner(3,ichunk) = IUPPERLOWER
+ itypecorner(1,ichunk) = IUPPERLOWER
+ itypecorner(2,ichunk) = ILOWERLOWER
+ itypecorner(3,ichunk) = IUPPERLOWER
- ichunk = 3
- iprocscorners(1,ichunk) = addressing_big(CHUNK_AB,0,0)
- iprocscorners(2,ichunk) = addressing_big(CHUNK_AC,NPROC_XI-1,0)
- iprocscorners(3,ichunk) = addressing_big(CHUNK_BC_ANTIPODE,NPROC_XI-1,NPROC_ETA-1)
+ ichunk = 3
+ iprocscorners(1,ichunk) = addressing_big(CHUNK_AB,0,0)
+ iprocscorners(2,ichunk) = addressing_big(CHUNK_AC,NPROC_XI-1,0)
+ iprocscorners(3,ichunk) = addressing_big(CHUNK_BC_ANTIPODE,NPROC_XI-1,NPROC_ETA-1)
- itypecorner(1,ichunk) = ILOWERLOWER
- itypecorner(2,ichunk) = IUPPERLOWER
- itypecorner(3,ichunk) = IUPPERUPPER
+ itypecorner(1,ichunk) = ILOWERLOWER
+ itypecorner(2,ichunk) = IUPPERLOWER
+ itypecorner(3,ichunk) = IUPPERUPPER
- ichunk = 4
- iprocscorners(1,ichunk) = addressing_big(CHUNK_AB,NPROC_XI-1,NPROC_ETA-1)
- iprocscorners(2,ichunk) = addressing_big(CHUNK_BC,NPROC_XI-1,NPROC_ETA-1)
- iprocscorners(3,ichunk) = addressing_big(CHUNK_AC_ANTIPODE,0,NPROC_ETA-1)
+ ichunk = 4
+ iprocscorners(1,ichunk) = addressing_big(CHUNK_AB,NPROC_XI-1,NPROC_ETA-1)
+ iprocscorners(2,ichunk) = addressing_big(CHUNK_BC,NPROC_XI-1,NPROC_ETA-1)
+ iprocscorners(3,ichunk) = addressing_big(CHUNK_AC_ANTIPODE,0,NPROC_ETA-1)
- itypecorner(1,ichunk) = IUPPERUPPER
- itypecorner(2,ichunk) = IUPPERUPPER
- itypecorner(3,ichunk) = ILOWERUPPER
+ itypecorner(1,ichunk) = IUPPERUPPER
+ itypecorner(2,ichunk) = IUPPERUPPER
+ itypecorner(3,ichunk) = ILOWERUPPER
- ichunk = 5
- iprocscorners(1,ichunk) = addressing_big(CHUNK_AC,0,0)
- iprocscorners(2,ichunk) = addressing_big(CHUNK_BC_ANTIPODE,0,NPROC_ETA-1)
- iprocscorners(3,ichunk) = addressing_big(CHUNK_AB_ANTIPODE,NPROC_XI-1,0)
+ ichunk = 5
+ iprocscorners(1,ichunk) = addressing_big(CHUNK_AC,0,0)
+ iprocscorners(2,ichunk) = addressing_big(CHUNK_BC_ANTIPODE,0,NPROC_ETA-1)
+ iprocscorners(3,ichunk) = addressing_big(CHUNK_AB_ANTIPODE,NPROC_XI-1,0)
- itypecorner(1,ichunk) = ILOWERLOWER
- itypecorner(2,ichunk) = ILOWERUPPER
- itypecorner(3,ichunk) = IUPPERLOWER
+ itypecorner(1,ichunk) = ILOWERLOWER
+ itypecorner(2,ichunk) = ILOWERUPPER
+ itypecorner(3,ichunk) = IUPPERLOWER
- ichunk = 6
- iprocscorners(1,ichunk) = addressing_big(CHUNK_AC_ANTIPODE,NPROC_XI-1,0)
- iprocscorners(2,ichunk) = addressing_big(CHUNK_BC_ANTIPODE,0,0)
- iprocscorners(3,ichunk) = addressing_big(CHUNK_AB_ANTIPODE,0,0)
+ ichunk = 6
+ iprocscorners(1,ichunk) = addressing_big(CHUNK_AC_ANTIPODE,NPROC_XI-1,0)
+ iprocscorners(2,ichunk) = addressing_big(CHUNK_BC_ANTIPODE,0,0)
+ iprocscorners(3,ichunk) = addressing_big(CHUNK_AB_ANTIPODE,0,0)
- itypecorner(1,ichunk) = IUPPERLOWER
- itypecorner(2,ichunk) = ILOWERLOWER
- itypecorner(3,ichunk) = ILOWERLOWER
+ itypecorner(1,ichunk) = IUPPERLOWER
+ itypecorner(2,ichunk) = ILOWERLOWER
+ itypecorner(3,ichunk) = ILOWERLOWER
- ichunk = 7
- iprocscorners(1,ichunk) = addressing_big(CHUNK_AC,0,NPROC_ETA-1)
- iprocscorners(2,ichunk) = addressing_big(CHUNK_BC,0,0)
- iprocscorners(3,ichunk) = addressing_big(CHUNK_AB_ANTIPODE,NPROC_XI-1,NPROC_ETA-1)
+ ichunk = 7
+ iprocscorners(1,ichunk) = addressing_big(CHUNK_AC,0,NPROC_ETA-1)
+ iprocscorners(2,ichunk) = addressing_big(CHUNK_BC,0,0)
+ iprocscorners(3,ichunk) = addressing_big(CHUNK_AB_ANTIPODE,NPROC_XI-1,NPROC_ETA-1)
- itypecorner(1,ichunk) = ILOWERUPPER
- itypecorner(2,ichunk) = ILOWERLOWER
- itypecorner(3,ichunk) = IUPPERUPPER
+ itypecorner(1,ichunk) = ILOWERUPPER
+ itypecorner(2,ichunk) = ILOWERLOWER
+ itypecorner(3,ichunk) = IUPPERUPPER
- ichunk = 8
- iprocscorners(1,ichunk) = addressing_big(CHUNK_BC,0,NPROC_ETA-1)
- iprocscorners(2,ichunk) = addressing_big(CHUNK_AC_ANTIPODE,NPROC_XI-1,NPROC_ETA-1)
- iprocscorners(3,ichunk) = addressing_big(CHUNK_AB_ANTIPODE,0,NPROC_ETA-1)
+ ichunk = 8
+ iprocscorners(1,ichunk) = addressing_big(CHUNK_BC,0,NPROC_ETA-1)
+ iprocscorners(2,ichunk) = addressing_big(CHUNK_AC_ANTIPODE,NPROC_XI-1,NPROC_ETA-1)
+ iprocscorners(3,ichunk) = addressing_big(CHUNK_AB_ANTIPODE,0,NPROC_ETA-1)
- itypecorner(1,ichunk) = ILOWERUPPER
- itypecorner(2,ichunk) = IUPPERUPPER
- itypecorner(3,ichunk) = ILOWERUPPER
+ itypecorner(1,ichunk) = ILOWERUPPER
+ itypecorner(2,ichunk) = IUPPERUPPER
+ itypecorner(3,ichunk) = ILOWERUPPER
endif
! file to store the list of processors for each message for corners
if(myrank == 0) open(unit=IOUT,file=trim(OUTPUT_FILES)//'/list_messages_corners.txt',status='unknown')
-! loop over all the messages to create the addressing
+ ! loop over all the messages to create the addressing
do imsg = 1,NCORNERSCHUNKS
- if(myrank == 0) write(IMAIN,*) 'Generating message ',imsg,' for corners out of ',NCORNERSCHUNKS
+ if(myrank == 0) write(IMAIN,*) 'Generating message ',imsg,' for corners out of ',NCORNERSCHUNKS
-! save triplet of processors in list of messages
+ ! save triplet of processors in list of messages
if(myrank == 0) write(IOUT,*) iprocscorners(1,imsg),iprocscorners(2,imsg),iprocscorners(3,imsg)
! loop on the three processors of a given corner
@@ -954,7 +956,7 @@
if(myrank == 0) close(IOUT)
-! deallocate arrays
+ ! deallocate arrays
deallocate(iproc_sender)
deallocate(iproc_receiver)
deallocate(npoin2D_send)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/save_arrays_solver.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/save_arrays_solver.f90 2013-07-01 11:55:54 UTC (rev 22477)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/save_arrays_solver.f90 2013-07-01 12:46:50 UTC (rev 22478)
@@ -85,7 +85,7 @@
integer ibool(NGLLX,NGLLY,NGLLZ,nspec)
-! doubling mesh flag
+ ! doubling mesh flag
integer, dimension(nspec) :: idoubling
! this for non blocking MPI
@@ -165,6 +165,7 @@
open(unit=27,file=prname(1:len_trim(prname))//'solver_data_1.bin',status='unknown',form='unformatted',action='write')
+ ! local GLL points
write(27) xixstore
write(27) xiystore
write(27) xizstore
@@ -184,7 +185,7 @@
close(29)
endif
-! other terms needed in the solid regions only
+ ! other terms needed in the solid regions only
if(iregion_code /= IREGION_OUTER_CORE) then
! note: muvstore needed for Q_mu shear attenuation in inner core
@@ -341,8 +342,10 @@
close(27)
-! absorbing boundary parameters
- open(unit=27,file=prname(1:len_trim(prname))//'boundary.bin',status='unknown',form='unformatted',action='write')
+ ! absorbing boundary parameters
+ open(unit=27,file=prname(1:len_trim(prname))//'boundary.bin', &
+ status='unknown',form='unformatted',action='write',iostat=ier)
+ if( ier /= 0 ) call exit_mpi(myrank,'error opening boundary.bin file')
write(27) nspec2D_xmin
write(27) nspec2D_xmax
@@ -377,12 +380,15 @@
!> Hejun
! No matter 1D or 3D Attenuation, we save value for gll points
if(ATTENUATION) then
- open(unit=27, file=prname(1:len_trim(prname))//'attenuation.bin', status='unknown', form='unformatted',action='write')
- write(27) tau_s
- write(27) tau_e_store
- write(27) Qmu_store
- write(27) T_c_source
- close(27)
+ open(unit=27, file=prname(1:len_trim(prname))//'attenuation.bin', &
+ status='unknown', form='unformatted',action='write',iostat=ier)
+ if( ier /= 0 ) call exit_mpi(myrank,'error opening attenuation.bin file')
+
+ write(27) tau_s
+ write(27) tau_e_store
+ write(27) Qmu_store
+ write(27) T_c_source
+ close(27)
endif
! uncomment for vp & vs model storage
@@ -400,7 +406,10 @@
write(27) sqrt( muvstore/rhostore )*scaleval1
close(27)
! rho
- open(unit=27,file=prname(1:len_trim(prname))//'rho.bin',status='unknown',form='unformatted',action='write')
+ open(unit=27,file=prname(1:len_trim(prname))//'rho.bin', &
+ status='unknown',form='unformatted',action='write',iostat=ier)
+ if( ier /= 0 ) call exit_mpi(myrank,'error opening rho.bin file')
+
write(27) rhostore*scaleval2
close(27)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/model_topo_bathy.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/model_topo_bathy.f90 2013-07-01 11:55:54 UTC (rev 22477)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/model_topo_bathy.f90 2013-07-01 12:46:50 UTC (rev 22478)
@@ -33,7 +33,6 @@
! by default (constants.h), it uses a smoothed ETOPO 4 dataset
!--------------------------------------------------------------------------------------------------
-
subroutine model_topo_bathy_broadcast(myrank,ibathy_topo)
! standard routine to setup model
@@ -44,10 +43,12 @@
! standard include of the MPI library
include 'mpif.h'
+ integer :: myrank
+
! bathymetry and topography: use integer array to store values
integer, dimension(NX_BATHY,NY_BATHY) :: ibathy_topo
- integer :: myrank
+ ! local parameters
integer :: ier
if(myrank == 0) call read_topo_bathy_file(ibathy_topo)
@@ -71,10 +72,11 @@
character(len=150) topo_bathy_file
-! use integer array to store values
+ ! use integer array to store values
integer, dimension(NX_BATHY,NY_BATHY) :: ibathy_topo
- integer itopo_x,itopo_y,ier
+ ! local parameters
+ integer :: itopo_x,itopo_y,ier
call get_value_string(topo_bathy_file, 'model.topoBathy.PATHNAME_TOPO_FILE', PATHNAME_TOPO_FILE)
@@ -84,6 +86,7 @@
print*,'error opening:',trim(topo_bathy_file)
call exit_mpi(0,'error opening topography data file')
endif
+
! reads in topography array
do itopo_y=1,NY_BATHY
do itopo_x=1,NX_BATHY
@@ -92,7 +95,6 @@
enddo
close(13)
-
! note: we check the limits after reading in the data. this seems to perform sligthly faster
! however, reading ETOPO1.xyz will take ~ 2m 1.2s for a single process
@@ -114,9 +116,11 @@
endif
+ ! user output
+ write(IMAIN,*) " topography/bathymetry: min/max = ",minval(ibathy_topo),maxval(ibathy_topo)
+
end subroutine read_topo_bathy_file
-
!
!-------------------------------------------------------------------------------------------------
!
@@ -131,23 +135,32 @@
include "constants.h"
-! use integer array to store values
- integer, dimension(NX_BATHY,NY_BATHY) :: ibathy_topo
+ ! location latitude/longitude (in degree)
+ double precision,intent(in):: xlat,xlon
- double precision xlat,xlon,value
+ ! returns elevation (in meters)
+ double precision,intent(out):: value
- integer iadd1,iel1
- double precision samples_per_degree_topo
- double precision xlo
- double precision:: lon_corner,lat_corner,ratio_lon,ratio_lat
+ ! use integer array to store values
+ integer, dimension(NX_BATHY,NY_BATHY),intent(in) :: ibathy_topo
+ ! local parameters
+ integer :: iadd1,iel1
+ double precision :: samples_per_degree_topo
+ double precision :: xlo
+ double precision :: lon_corner,lat_corner,ratio_lon,ratio_lat
+
+ ! initializes elevation
+ value = ZERO
+
+ ! longitude within range [0,360] degrees
xlo = xlon
if(xlon < 0.d0) xlo = xlo + 360.d0
-! compute number of samples per degree
+ ! compute number of samples per degree
samples_per_degree_topo = dble(RESOLUTION_TOPO_FILE) / 60.d0
-! compute offset in data file and avoid edge effects
+ ! compute offset in data file and avoid edge effects
iadd1 = 1 + int((90.d0-xlat)/samples_per_degree_topo)
if(iadd1 < 1) iadd1 = 1
if(iadd1 > NY_BATHY) iadd1 = NY_BATHY
@@ -156,7 +169,8 @@
if(iel1 <= 0 .or. iel1 > NX_BATHY) iel1 = NX_BATHY
! Use bilinear interpolation rather nearest point interpolation
-! convert integer value to double precision
+
+ ! convert integer value to double precision
! value = dble(ibathy_topo(iel1,iadd1))
lon_corner=iel1*samples_per_degree_topo
@@ -170,7 +184,7 @@
if(ratio_lat<0.0) ratio_lat=0.0
if(ratio_lat>1.0) ratio_lat=1.0
-! convert integer value to double precision
+ ! convert integer value to double precision
if( iadd1 <= NY_BATHY-1 .and. iel1 <= NX_BATHY-1 ) then
! interpolates for points within boundaries
value = dble(ibathy_topo(iel1,iadd1))*(1-ratio_lon)*(1.-ratio_lat) &
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_movie_volume.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_movie_volume.f90 2013-07-01 11:55:54 UTC (rev 22477)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_movie_volume.f90 2013-07-01 12:46:50 UTC (rev 22478)
@@ -25,10 +25,6 @@
!
!=====================================================================
-! create file OUTPUT_FILES/values_from_mesher.h based upon DATA/Par_file
-! in order to compile the solver with the right array sizes
-
-!---------------------------------------------------------------------------------
! this subroutine counts the number of points and elements within the movie volume
! in this processor slice, and returns arrays that keep track of them, both in global and local indexing schemes
@@ -110,7 +106,10 @@
end subroutine count_points_movie_volume
-! -----------------------------------------------------------------
+!
+!-------------------------------------------------------------------------------------------------
+!
+
! writes meshfiles to merge with solver snapshots for 3D volume movies. Also computes and outputs
! the rotation matrix nu_3dmovie required to transfer to a geographic coordinate system
@@ -255,7 +254,10 @@
end subroutine write_movie_volume_mesh
-! ---------------------------------------------
+!
+!-------------------------------------------------------------------------------------------------
+!
+
subroutine write_movie_volume_strains(myrank,npoints_3dmovie,LOCAL_PATH,MOVIE_VOLUME_TYPE,MOVIE_COARSE, &
it,muvstore_crust_mantle_3dmovie,mask_3dmovie,nu_3dmovie,&
hprime_xx,hprime_yy,hprime_zz,ibool_crust_mantle,&
@@ -263,6 +265,7 @@
etax_crust_mantle,etay_crust_mantle,etaz_crust_mantle,&
gammax_crust_mantle,gammay_crust_mantle,gammaz_crust_mantle,displ_crust_mantle)
+! outputs strains: MOVIE_VOLUME_TYPE == 1 / 2 / 3
implicit none
@@ -309,14 +312,15 @@
allocate(store_val3d_NZ(npoints_3dmovie))
allocate(store_val3d_EZ(npoints_3dmovie))
+ ! check
if(NDIM /= 3) call exit_MPI(myrank, 'write_movie_volume requires NDIM = 3')
if(MOVIE_VOLUME_TYPE == 1) then
- movie_prefix='E' ! strain
+ movie_prefix='E' ! strain
else if(MOVIE_VOLUME_TYPE == 2) then
- movie_prefix='S' ! time integral of strain
+ movie_prefix='S' ! time integral of strain
else if(MOVIE_VOLUME_TYPE == 3) then
- movie_prefix='P' ! potency, or integral of strain x \mu
+ movie_prefix='P' ! potency, or integral of strain x \mu
endif
if(MOVIE_COARSE) then
NIT = NGLLX-1
@@ -404,7 +408,10 @@
end subroutine write_movie_volume_strains
-! ---------------------------------------------
+!
+!-------------------------------------------------------------------------------------------------
+!
+
subroutine write_movie_volume_vector(myrank,it,npoints_3dmovie,LOCAL_PATH,MOVIE_VOLUME_TYPE, &
MOVIE_COARSE,ibool_crust_mantle,vector_crust_mantle, &
scalingval,mask_3dmovie,nu_3dmovie)
@@ -421,6 +428,7 @@
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: ibool_crust_mantle
real(kind=CUSTOM_REAL), dimension(3,NGLOB_CRUST_MANTLE) :: vector_crust_mantle
real(kind=CUSTOM_REAL), dimension(3,3,npoints_3dmovie) :: nu_3dmovie
+
double precision :: scalingval
real(kind=CUSTOM_REAL) :: scalingval_to_use
real(kind=CUSTOM_REAL), dimension(3) :: vector_local,vector_local_new
@@ -437,14 +445,15 @@
! check
if(NDIM /= 3) call exit_MPI(myrank,'write_movie_volume requires NDIM = 3')
+ ! allocates arrays
allocate(store_val3d_N(npoints_3dmovie))
allocate(store_val3d_E(npoints_3dmovie))
allocate(store_val3d_Z(npoints_3dmovie))
if(MOVIE_VOLUME_TYPE == 5) then
- movie_prefix='DI' ! displacement
+ movie_prefix='DI' ! displacement
else if(MOVIE_VOLUME_TYPE == 6) then
- movie_prefix='VE' ! velocity
+ movie_prefix='VE' ! velocity
endif
if(MOVIE_COARSE) then
@@ -506,16 +515,18 @@
!
subroutine write_movie_volume_divcurl(myrank,it,eps_trace_over_3_crust_mantle,&
- div_displ_outer_core, &
- accel_outer_core,kappavstore_outer_core,rhostore_outer_core,ibool_outer_core, &
- eps_trace_over_3_inner_core, &
- epsilondev_crust_mantle,epsilondev_inner_core, &
- LOCAL_PATH, &
- displ_crust_mantle,displ_inner_core,displ_outer_core, &
- veloc_crust_mantle,veloc_inner_core,veloc_outer_core, &
- accel_crust_mantle,accel_inner_core, &
- ibool_crust_mantle,ibool_inner_core)
+ div_displ_outer_core, &
+ accel_outer_core,kappavstore_outer_core,rhostore_outer_core,ibool_outer_core, &
+ eps_trace_over_3_inner_core, &
+ epsilondev_crust_mantle,epsilondev_inner_core, &
+ LOCAL_PATH, &
+ displ_crust_mantle,displ_inner_core,displ_outer_core, &
+ veloc_crust_mantle,veloc_inner_core,veloc_outer_core, &
+ accel_crust_mantle,accel_inner_core, &
+ ibool_crust_mantle,ibool_inner_core)
+! outputs divergence and curl: MOVIE_VOLUME_TYPE == 4
+
implicit none
include "constants.h"
include "OUTPUT_FILES/values_from_mesher.h"
@@ -558,7 +569,7 @@
! output parameters
logical,parameter :: MOVIE_OUTPUT_DIV = .true. ! divergence
logical,parameter :: MOVIE_OUTPUT_CURL = .false. ! curl
- logical,parameter :: MOVIE_OUTPUT_CURLNORM = .true. ! frobenius norm of curl
+ logical,parameter :: MOVIE_OUTPUT_CURLNORM = .true. ! Frobenius norm of curl
logical,parameter :: MOVIE_OUTPUT_DISPLNORM = .false. ! norm of displacement
logical,parameter :: MOVIE_OUTPUT_VELOCNORM = .false. ! norm of velocity
logical,parameter :: MOVIE_OUTPUT_ACCELNORM = .false. ! norm of acceleration
@@ -911,7 +922,6 @@
!-------------------------------------------------------------------------------------------------
-
! external mesh routine for saving vtk files for custom_real values on global points
subroutine write_VTK_data_cr(idoubling,nspec,nglob, &
More information about the CIG-COMMITS
mailing list