[cig-commits] r19219 - in seismo/3D/SPECFEM3D_GLOBE/trunk: . setup src/auxiliaries src/meshfem3D src/shared src/specfem3D
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Sat Nov 19 16:12:06 PST 2011
Author: dkomati1
Date: 2011-11-19 16:12:05 -0800 (Sat, 19 Nov 2011)
New Revision: 19219
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/flags.guess
seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in
seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_vol_data.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/create_movie_AVS_DX.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/create_movie_GMT_global.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_regions_mesh.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_perm_color.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/meshfem3D.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/meshfem3D_models.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/get_model_parameters.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/param_reader.c
seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/write_c_binary.c
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/assemble_MPI_central_cube.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_element.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_crust_mantle_Dev.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_inner_core_Dev.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/noise_tomography.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/read_mesh_databases.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/setup_sources_receivers.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_movie_volume.f90
Log:
added Manh Ha's modifications to save header files for the C / StarSs versions of the code.
changed max number of colors from 10000 to 1000.
removed useless white spaces.
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/flags.guess
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/flags.guess 2011-11-18 23:19:05 UTC (rev 19218)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/flags.guess 2011-11-20 00:12:05 UTC (rev 19219)
@@ -71,7 +71,7 @@
if test x"$FLAGS_CHECK" = x; then
FLAGS_CHECK="\$(FLAGS_NO_CHECK)"
fi
- # useful for debugging... -ffpe-trap=underflow,overflow,zero -fbacktrace
+ # useful for debugging... -ffpe-trap=underflow,overflow,zero -fbacktrace -fbounds-check
#
# older gfortran syntax
# note: -std=f95 can lead to problem: undefined reference to `system_'
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in 2011-11-18 23:19:05 UTC (rev 19218)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in 2011-11-20 00:12:05 UTC (rev 19219)
@@ -60,7 +60,7 @@
! add mesh coloring for the GPU + MPI implementation
logical, parameter :: USE_MESH_COLORING_GPU = .false.
- integer, parameter :: MAX_NUMBER_OF_COLORS = 10000
+ integer, parameter :: MAX_NUMBER_OF_COLORS = 1000
integer, parameter :: NGNOD_HEXAHEDRA = 8
!*********************************************************************************************************
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_vol_data.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_vol_data.f90 2011-11-18 23:19:05 UTC (rev 19218)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_vol_data.f90 2011-11-20 00:12:05 UTC (rev 19219)
@@ -309,7 +309,7 @@
! nothing to do for fictitious elements in central cube
if( idoubling_inner_core(ispec) == IFLAG_IN_FICTITIOUS_CUBE) cycle
endif
-
+
! writes out global point data
do k = 1, NGLLZ, dk
do j = 1, NGLLY, dj
@@ -381,7 +381,7 @@
do ispec = 1, nspec(it)
! checks if element counts
if (ir==3 ) then
- ! inner core
+ ! inner core
! fictitious elements in central cube
if( idoubling_inner_core(ispec) == IFLAG_IN_FICTITIOUS_CUBE) then
! connectivity must be given, otherwise element count would be wrong
@@ -399,12 +399,12 @@
call write_integer_fd(efd,1)
enddo ! i
enddo ! j
- enddo ! k
+ enddo ! k
! takes next element
cycle
- endif
+ endif
endif
-
+
! writes out element connectivity
numpoin = numpoin + 1 ! counts elements
do k = 1, NGLLZ-1, dk
@@ -436,7 +436,7 @@
call write_integer_fd(efd,n8)
enddo ! i
enddo ! j
- enddo ! k
+ enddo ! k
enddo ! ispec
np = np + npoint(it)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/create_movie_AVS_DX.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/create_movie_AVS_DX.f90 2011-11-18 23:19:05 UTC (rev 19218)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/create_movie_AVS_DX.f90 2011-11-20 00:12:05 UTC (rev 19219)
@@ -43,7 +43,7 @@
integer NPROCTOT,NEX_PER_PROC_XI,NEX_PER_PROC_ETA
logical MOVIE_SURFACE
integer USE_COMPONENT
-
+
! ************** PROGRAM STARTS HERE **************
call read_AVS_DX_parameters(NEX_XI,NEX_ETA, &
@@ -75,7 +75,7 @@
print *,'enter component (e.g. 1=Z, 2=N, 3=E)'
read(5,*) USE_COMPONENT
if( USE_COMPONENT < 1 .or. USE_COMPONENT > 3 ) stop 'component must be 1, 2 or 3'
-
+
! run the main program
call create_movie_AVS_DX(iformat,it1,it2, &
NEX_XI,NEX_ETA, &
@@ -134,12 +134,12 @@
real(kind=CUSTOM_REAL) RRval,rhoval
real(kind=CUSTOM_REAL) thetahat_x,thetahat_y,thetahat_z
real(kind=CUSTOM_REAL) phihat_x,phihat_y
-
+
double precision min_field_current,max_field_current,max_absol
logical USE_OPENDX,UNIQUE_FILE,USE_GMT,USE_AVS
integer USE_COMPONENT
-
+
integer iformat,nframes,iframe
character(len=150) outputname
@@ -340,12 +340,12 @@
read(IOUT) store_val_x
read(IOUT) store_val_y
read(IOUT) store_val_z
-
- ! reads in associated values (displacement or velocity..)
+
+ ! reads in associated values (displacement or velocity..)
read(IOUT) store_val_ux
read(IOUT) store_val_uy
read(IOUT) store_val_uz
-
+
close(IOUT)
! clear number of elements kept
@@ -425,7 +425,7 @@
phihat_x = -ycoord / rhoval
phihat_y = xcoord / rhoval
endif
-
+
displn(i,j) = displx*phihat_x + disply*phihat_y
endif
@@ -557,7 +557,7 @@
elseif( USE_COMPONENT == 2) then
write(outputname,"('/DX_movie_',i6.6,'.N.dx')") it
elseif( USE_COMPONENT == 3) then
- write(outputname,"('/DX_movie_',i6.6,'.E.dx')") it
+ write(outputname,"('/DX_movie_',i6.6,'.E.dx')") it
endif
open(unit=11,file=trim(OUTPUT_FILES)//trim(outputname),status='unknown')
write(11,*) 'object 1 class array type float rank 1 shape 3 items ',nglob,' data follows'
@@ -569,31 +569,31 @@
outputname = '/AVS_movie_all.N.inp'
elseif( USE_COMPONENT == 3) then
outputname = '/AVS_movie_all.E.inp'
- endif
+ endif
open(unit=11,file=trim(OUTPUT_FILES)//trim(outputname),status='unknown')
write(11,*) nframes
write(11,*) 'data'
write(11,"('step',i1,' image',i1)") 1,1
write(11,*) nglob,' ',nspectot_AVS_max
else if(.not. UNIQUE_FILE) then
- if( USE_COMPONENT == 1) then
+ if( USE_COMPONENT == 1) then
write(outputname,"('/AVS_movie_',i6.6,'.Z.inp')") it
elseif( USE_COMPONENT == 2) then
write(outputname,"('/AVS_movie_',i6.6,'.N.inp')") it
elseif( USE_COMPONENT == 3) then
write(outputname,"('/AVS_movie_',i6.6,'.E.inp')") it
- endif
+ endif
open(unit=11,file=trim(OUTPUT_FILES)//trim(outputname),status='unknown')
write(11,*) nglob,' ',nspectot_AVS_max,' 1 0 0'
endif
else if(USE_GMT) then
- if( USE_COMPONENT == 1) then
+ if( USE_COMPONENT == 1) then
write(outputname,"('/gmt_movie_',i6.6,'.Z.xyz')") it
elseif( USE_COMPONENT == 2) then
write(outputname,"('/gmt_movie_',i6.6,'.N.xyz')") it
elseif( USE_COMPONENT == 3) then
write(outputname,"('/gmt_movie_',i6.6,'.E.xyz')") it
- endif
+ endif
open(unit=11,file=trim(OUTPUT_FILES)//trim(outputname),status='unknown')
else
stop 'wrong output format selected'
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/create_movie_GMT_global.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/create_movie_GMT_global.f90 2011-11-18 23:19:05 UTC (rev 19218)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/create_movie_GMT_global.f90 2011-11-20 00:12:05 UTC (rev 19219)
@@ -566,7 +566,7 @@
phihat_x = -ycoord / rhoval
phihat_y = xcoord / rhoval
endif
-
+
displn(i,j) = displx*phihat_x + disply*phihat_y
endif
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_regions_mesh.f90 2011-11-18 23:19:05 UTC (rev 19218)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_regions_mesh.f90 2011-11-20 00:12:05 UTC (rev 19219)
@@ -30,6 +30,7 @@
iproc_xi,iproc_eta,ichunk,nspec,nspec_tiso, &
volume_local,area_local_bottom,area_local_top, &
nglob_theor,npointot, &
+ NSTEP,DT,NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX, &
NEX_XI,NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX, &
NSPEC2D_BOTTOM,NSPEC2D_TOP, &
@@ -239,6 +240,16 @@
!****************************************************************************************************
+!///////////////////////////////////////////////////////////////////////////////
+! Manh Ha - 18-11-2011
+! Adding new variables
+
+ integer :: NSTEP
+ double precision :: DT
+ integer :: NGLOB2DMAX_XMIN_XMAX, NGLOB2DMAX_YMIN_YMAX
+
+!///////////////////////////////////////////////////////////////////////////////
+
! Boundary Mesh
integer NSPEC2D_MOHO,NSPEC2D_400,NSPEC2D_670,nex_eta_moho
integer, dimension(:), allocatable :: ibelm_moho_top,ibelm_moho_bot,ibelm_400_top,ibelm_400_bot, &
@@ -874,6 +885,80 @@
endif
+!! DK DK and Manh Ha, Nov 2011: added this to use the new mesher in the CUDA or C / StarSs test codes
+
+ if (myrank == 0 .and. iregion_code == IREGION_CRUST_MANTLE) then
+
+! write a header file for the Fortran version of the solver
+ open(unit=99,file=prname(1:len_trim(prname))//'values_from_mesher_f90.h',status='unknown')
+ write(99,*) 'integer, parameter :: NSPEC = ',nspec
+ write(99,*) 'integer, parameter :: NGLOB = ',nglob
+!!! DK DK use 1000 time steps only for the scaling tests
+ write(99,*) 'integer, parameter :: NSTEP = 1000 !!!!!!!!!!! ',nstep
+ write(99,*) 'real(kind=4), parameter :: deltat = ',DT
+ write(99,*)
+ write(99,*) 'integer, parameter :: NGLOB2DMAX_XMIN_XMAX = ',NGLOB2DMAX_XMIN_XMAX
+ write(99,*) 'integer, parameter :: NGLOB2DMAX_YMIN_YMAX = ',NGLOB2DMAX_YMIN_YMAX
+ write(99,*) 'integer, parameter :: NGLOB2DMAX_ALL = ',max(NGLOB2DMAX_XMIN_XMAX, NGLOB2DMAX_YMIN_YMAX)
+ write(99,*) 'integer, parameter :: NPROC_XI = ',NPROC_XI
+ write(99,*) 'integer, parameter :: NPROC_ETA = ',NPROC_ETA
+ write(99,*)
+ write(99,*) '! element number of the source and of the station'
+ write(99,*) '! after permutation of the elements by mesh coloring'
+ write(99,*) '! and inner/outer set splitting in the mesher'
+ write(99,*) 'integer, parameter :: NSPEC_SOURCE = ',perm(NSPEC/3)
+ write(99,*) 'integer, parameter :: RANK_SOURCE = 0'
+ write(99,*)
+ write(99,*) 'integer, parameter :: RANK_STATION = (NPROC_XI*NPROC_ETA - 1)'
+ write(99,*) 'integer, parameter :: NSPEC_STATION = ',perm(2*NSPEC/3)
+
+! save coordinates of the seismic source
+! write(99,*) xstore(2,2,2,10);
+! write(99,*) ystore(2,2,2,10);
+! write(99,*) zstore(2,2,2,10);
+
+! save coordinates of the seismic station
+! write(99,*) xstore(2,2,2,nspec-10);
+! write(99,*) ystore(2,2,2,nspec-10);
+! write(99,*) zstore(2,2,2,nspec-10);
+ close(99)
+
+!! write a header file for the C version of the solver
+ open(unit=99,file=prname(1:len_trim(prname))//'values_from_mesher_C.h',status='unknown')
+ write(99,*) '#define NSPEC ',nspec
+ write(99,*) '#define NGLOB ',nglob
+!! write(99,*) '#define NSTEP ',nstep
+!!! DK DK use 1000 time steps only for the scaling tests
+ write(99,*) '// #define NSTEP ',nstep
+ write(99,*) '#define NSTEP 1000'
+! put an "f" at the end to force single precision
+ write(99,"('#define deltat ',e18.10,'f')") DT
+ write(99,*) '#define NGLOB2DMAX_XMIN_XMAX ',NGLOB2DMAX_XMIN_XMAX
+ write(99,*) '#define NGLOB2DMAX_YMIN_YMAX ',NGLOB2DMAX_YMIN_YMAX
+ write(99,*) '#define NGLOB2DMAX_ALL ',max(NGLOB2DMAX_XMIN_XMAX, NGLOB2DMAX_YMIN_YMAX)
+ write(99,*) '#define NPROC_XI ',NPROC_XI
+ write(99,*) '#define NPROC_ETA ',NPROC_ETA
+ write(99,*)
+ write(99,*) '// element and MPI slice number of the source and the station'
+ write(99,*) '// after permutation of the elements by mesh coloring'
+ write(99,*) '// and inner/outer set splitting in the mesher'
+ write(99,*) '#define RANK_SOURCE 0'
+ write(99,*) '#define NSPEC_SOURCE ',perm(NSPEC/3)
+ write(99,*)
+ write(99,*) '#define RANK_STATION (NPROC_XI*NPROC_ETA - 1)'
+ write(99,*) '#define NSPEC_STATION ',perm(2*NSPEC/3)
+ close(99)
+
+ open(unit=99,file=prname(1:len_trim(prname))//'values_from_mesher_nspec_outer.h',status='unknown')
+ write(99,*) '#define NSPEC_OUTER ',nspec_outer_max_global
+ write(99,*) '// NSPEC_OUTER_min = ',nspec_outer_min_global
+ write(99,*) '// NSPEC_OUTER_max = ',nspec_outer_max_global
+ close(99)
+
+ endif
+
+!! DK DK and Manh Ha, Nov 2011: added this to use the new mesher in the CUDA or C / StarSs test codes
+
deallocate(perm)
else
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_perm_color.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_perm_color.f90 2011-11-18 23:19:05 UTC (rev 19218)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_perm_color.f90 2011-11-20 00:12:05 UTC (rev 19219)
@@ -14,6 +14,9 @@
include "constants.h"
+! local variables
+ integer nspec, nglob
+
logical, dimension(nspec) :: is_on_a_slice_edge
integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
@@ -22,9 +25,6 @@
integer, dimension(MAX_NUMBER_OF_COLORS) :: first_elem_number_in_this_color
integer :: nb_colors_outer_elements,nb_colors_inner_elements,nspec_outer,myrank
-! local variables
- integer nspec, nglob
-
call get_color_faster(ibool, is_on_a_slice_edge, myrank, nspec, nglob, &
color, nb_colors_outer_elements, nb_colors_inner_elements, nspec_outer)
@@ -52,15 +52,15 @@
include "constants.h"
+! local variables
+ integer nspec,nglob
+
logical, dimension(nspec) :: is_on_a_slice_edge
integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
integer, dimension(nspec) :: color
integer :: nb_colors_outer_elements,nb_colors_inner_elements,myrank,nspec_outer
-! local variables
- integer nspec,nglob
-
integer :: ispec
!! DK DK for mesh coloring GPU Joseph Fourier
@@ -248,6 +248,9 @@
include "constants.h"
+! local variables
+ integer nspec,nglob_GLL_full
+
logical, dimension(nspec) :: is_on_a_slice_edge
integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
@@ -256,12 +259,11 @@
integer, dimension(MAX_NUMBER_OF_COLORS) :: first_elem_number_in_this_color
integer :: nb_colors_outer_elements,nb_colors_inner_elements,nspec_outer,myrank
-! local variables
- integer nspec,nglob_GLL_full
-
! a neighbor of a hexahedral node is a hexahedron that shares a face with it -> max degree of a node = 6
integer, parameter :: MAX_NUMBER_OF_NEIGHBORS = 100
+ integer nglob_eight_corners_only,nglob
+
! global corner numbers that need to be created
integer, dimension(nglob) :: global_corner_number
@@ -269,12 +271,8 @@
integer, dimension(:), allocatable :: ne,np,adj
integer xadj(nspec+1)
- !logical maskel(nspec)
-
integer i,istart,istop,number_of_neighbors
- integer nglob_eight_corners_only,nglob
-
! only count the total size of the array that will be created, or actually create it
logical count_only
integer total_size_ne,total_size_adj
@@ -672,6 +670,8 @@
integer nglob_eight_corners_only
+ integer nspec,iad,ispec,istart,istop,ino,node,jstart,jstop,nelem,jel
+
integer, intent(in) :: mn(nspec*NGNOD_HEXAHEDRA),mp(nspec+1),ne(total_size_ne),np(nglob_eight_corners_only+1)
integer, intent(out) :: adj(total_size_adj),xadj(nspec+1)
@@ -679,8 +679,6 @@
logical maskel(nspec)
integer countel(nspec)
- integer nspec,iad,ispec,istart,istop,ino,node,jstart,jstop,nelem,jel
-
xadj(1) = 1
iad = 1
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/meshfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/meshfem3D.f90 2011-11-18 23:19:05 UTC (rev 19218)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/meshfem3D.f90 2011-11-20 00:12:05 UTC (rev 19219)
@@ -663,6 +663,7 @@
iproc_xi,iproc_eta,ichunk,NSPEC(iregion_code),nspec_tiso, &
volume_local,area_local_bottom,area_local_top, &
nglob(iregion_code),npointot, &
+ NSTEP,DT,NGLOB2DMAX_XMIN_XMAX(iregion_code),NGLOB2DMAX_YMIN_YMAX(iregion_code), &
NEX_XI,NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
NSPEC2DMAX_XMIN_XMAX(iregion_code),NSPEC2DMAX_YMIN_YMAX(iregion_code), &
NSPEC2D_BOTTOM(iregion_code),NSPEC2D_TOP(iregion_code), &
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/meshfem3D_models.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/meshfem3D_models.f90 2011-11-18 23:19:05 UTC (rev 19218)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/meshfem3D_models.f90 2011-11-20 00:12:05 UTC (rev 19219)
@@ -431,7 +431,7 @@
! to create a reference model based on 1D_REF but with 3D crust and 410/660 topography
logical,parameter :: USE_1D_REFERENCE = .false.
-
+
end module meshfem3D_models_par
@@ -905,9 +905,9 @@
dvpv = 0.d0
dvph = 0.d0
dvsv = 0.d0
- dvsh = 0.d0
+ dvsh = 0.d0
endif
-
+
if(TRANSVERSE_ISOTROPY) then
vpv=vpv*(1.0d0+dble(dvpv))
vph=vph*(1.0d0+dble(dvph))
@@ -1085,7 +1085,7 @@
!
!---
found_crust = .false.
-
+
! crustal model can vary for different 3-D models
select case (THREE_D_MODEL )
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/get_model_parameters.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/get_model_parameters.f90 2011-11-18 23:19:05 UTC (rev 19218)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/get_model_parameters.f90 2011-11-20 00:12:05 UTC (rev 19219)
@@ -353,7 +353,7 @@
TRANSVERSE_ISOTROPY = .true. ! to use transverse isotropic PREM 1D ref model
! CRUSTAL = .true. ! with 3D crust: depends on 3D mantle reference model
! THREE_D_MODEL = 0 ! for default crustal model
- ! REFERENCE_1D_MODEL = REFERENCE_MODEL_AK135
+ ! REFERENCE_1D_MODEL = REFERENCE_MODEL_AK135
! TRANSVERSE_ISOTROPY = .false. ! for AK135 ref model
else if(MODEL_ROOT == 'heterogen') then
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/param_reader.c
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/param_reader.c 2011-11-18 23:19:05 UTC (rev 19218)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/param_reader.c 2011-11-20 00:12:05 UTC (rev 19219)
@@ -143,7 +143,7 @@
/* Regular expression for parsing lines from param file.
** Good luck reading this regular expression. Basically, the lines of
** the parameter file should be of the form 'parameter = value',
- ** optionally followed by a #-delimited comment.
+ ** optionally followed by a #-delimited comment.
** 'value' can be any number of space- or tab-separated words. Blank
** lines, lines containing only white space and lines whose first non-
** whitespace character is '#' are ignored. White space is generally
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/write_c_binary.c
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/write_c_binary.c 2011-11-18 23:19:05 UTC (rev 19218)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/write_c_binary.c 2011-11-20 00:12:05 UTC (rev 19219)
@@ -183,7 +183,7 @@
char * blank;
FILE *ft;
int ret;
-
+
// checks filesize
if( *filesize == 0 ){
perror("Error file size for reading");
@@ -206,9 +206,9 @@
ret = setvbuf( ft, work_buffer[*fid], _IOFBF, (size_t)MAX_B );
if( ret != 0 ){
perror("Error setting working buffer");
- exit(EXIT_FAILURE);
+ exit(EXIT_FAILURE);
}
-
+
// stores file index id fid: from 0 to 8
fp_abs[*fid] = ft;
@@ -227,7 +227,7 @@
char * blank;
FILE *ft;
int ret;
-
+
// checks filesize
if( *filesize == 0 ){
perror("Error file size for writing");
@@ -250,9 +250,9 @@
ret = setvbuf( ft, work_buffer[*fid], _IOFBF, (size_t)MAX_B );
if( ret != 0 ){
perror("Error setting working buffer");
- exit(EXIT_FAILURE);
+ exit(EXIT_FAILURE);
}
-
+
// stores file index id fid: from 0 to 8
fp_abs[*fid] = ft;
@@ -330,7 +330,7 @@
perror("Error fseek");
exit(EXIT_FAILURE);
}
-
+
donelen = 0;
remlen = *length;
buf = buffer;
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/assemble_MPI_central_cube.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/assemble_MPI_central_cube.f90 2011-11-18 23:19:05 UTC (rev 19218)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/assemble_MPI_central_cube.f90 2011-11-20 00:12:05 UTC (rev 19219)
@@ -41,13 +41,14 @@
! for matching with central cube in inner core
integer, intent(in) :: ichunk, nb_msgs_theor_in_cube, npoin2D_cube_from_slices
+ integer, intent(in) :: ndim_assemble
+ integer, intent(in) :: receiver_cube_from_slices
integer, intent(inout) :: iphase_CC
integer, dimension(nb_msgs_theor_in_cube), intent(in) :: sender_from_slices_to_cube
double precision, dimension(npoin2D_cube_from_slices,ndim_assemble), intent(inout) :: buffer_slices
double precision, dimension(npoin2D_cube_from_slices,ndim_assemble,nb_msgs_theor_in_cube), intent(inout) :: &
buffer_all_cube_from_slices
integer, dimension(nb_msgs_theor_in_cube,npoin2D_cube_from_slices), intent(in) :: ibool_central_cube
- integer, intent(in) :: receiver_cube_from_slices
! local to global mapping
integer, intent(in) :: NSPEC2D_BOTTOM_INNER_CORE
@@ -56,7 +57,6 @@
integer, dimension(NSPEC2D_BOTTOM_INNER_CORE), intent(in) :: ibelm_bottom_inner_core
! vector
- integer, intent(in) :: ndim_assemble
real(kind=CUSTOM_REAL), dimension(ndim_assemble,NGLOB_INNER_CORE), intent(inout) :: vector_assemble
integer ipoin,idimension, ispec2D, ispec
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_element.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_element.F90 2011-11-18 23:19:05 UTC (rev 19218)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_element.F90 2011-11-20 00:12:05 UTC (rev 19219)
@@ -29,20 +29,20 @@
! #undef _HANDOPT : turns hand-optimized code off
! or compile with: -D_HANDOPT
!#define _HANDOPT
-
+
! note: these hand optimizations should help compilers to pipeline the code and make better use of the cache;
! depending on compilers, it can further decrease the computation time by ~ 30%.
! the original routines are commented with "! way 1", the hand-optimized routines with "! way 2"
subroutine compute_element_iso(ispec, &
minus_gravity_table,density_table,minus_deriv_gravity_table, &
- xstore,ystore,zstore, &
+ xstore,ystore,zstore, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
wgll_cube, &
kappavstore,muvstore, &
ibool, &
R_memory,epsilon_trace_over_3, &
- one_minus_sum_beta,vx,vy,vz,vnspec, &
+ one_minus_sum_beta,vx,vy,vz,vnspec, &
tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,rho_s_H)
@@ -93,7 +93,7 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: dummyx_loc,dummyy_loc,dummyz_loc
-
+
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NGLLZ) :: rho_s_H
real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc
@@ -121,11 +121,11 @@
integer :: ispec_strain
integer :: i,j,k
- integer :: int_radius
- integer :: iglob1
-
+ integer :: int_radius
+ integer :: iglob1
+
! isotropic element
-
+
do k=1,NGLLZ
do j=1,NGLLY
do i=1,NGLLX
@@ -260,7 +260,7 @@
! Cartesian components of gradient of gravitational acceleration
! obtained from spherical components
- minus_g_over_radius = minus_g / radius
+ minus_g_over_radius = minus_g / radius
minus_dg_plus_g_over_radius = minus_dg - minus_g_over_radius
Hxxl = minus_g_over_radius*(cos_phi_sq*cos_theta_sq + sin_phi_sq) + cos_phi_sq*minus_dg*sin_theta_sq
@@ -330,9 +330,9 @@
endif ! end of section with gravity terms
! form dot product with test vector, non-symmetric form
- tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_yx*xiyl + sigma_zx*xizl)
- tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_zy*xizl)
- tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl)
+ tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_yx*xiyl + sigma_zx*xizl)
+ tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_zy*xizl)
+ tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl)
tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_yx*etayl + sigma_zx*etazl)
tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_zy*etazl)
@@ -353,10 +353,10 @@
!
!--------------------------------------------------------------------------------------------------
!
-
-
-
-
+
+
+
+
subroutine compute_element_tiso(ispec, &
minus_gravity_table,density_table,minus_deriv_gravity_table, &
xstore,ystore,zstore, &
@@ -466,7 +466,7 @@
integer :: iglob1
! transverse isotropic element
-
+
do k=1,NGLLZ
do j=1,NGLLY
do i=1,NGLLX
@@ -566,12 +566,12 @@
! use mesh coordinates to get theta and phi
! ystore and zstore contain theta and phi
iglob1 = ibool(i,j,k,ispec)
-
+
theta = ystore(iglob1)
phi = zstore(iglob1)
! precompute some products to reduce the CPU time
-
+
costheta = cos(theta)
sintheta = sin(theta)
cosphi = cos(phi)
@@ -616,11 +616,11 @@
! way 2: pre-compute temporary values
templ1 = four_rhovsvsq - rhovpvsq + twoetaminone*rhovphsq - four_eta_aniso*rhovsvsq
- templ1_cos = rhovphsq - rhovpvsq + costwotheta*templ1
+ templ1_cos = rhovphsq - rhovpvsq + costwotheta*templ1
templ2 = four_rhovsvsq - rhovpvsq - rhovphsq + two_eta_aniso*rhovphsq - four_eta_aniso*rhovsvsq
- templ2_cos = rhovpvsq - rhovphsq + costwotheta*templ2
+ templ2_cos = rhovpvsq - rhovphsq + costwotheta*templ2
templ3 = rhovphsq + rhovpvsq - two_eta_aniso*rhovphsq + four_eta_aniso*rhovsvsq
- templ3_two = templ3 - two_rhovshsq - two_rhovsvsq
+ templ3_two = templ3 - two_rhovshsq - two_rhovsvsq
templ3_cos = templ3_two + costwotheta*templ2
! way 2: reordering operations to facilitate compilation, avoiding divisions, using locality for temporary values
@@ -652,7 +652,7 @@
( 0.5*cosphisq*(templ2_cos + four_rhovshsq - four_rhovsvsq) &
+ sinphisq*(etaminone*rhovphsq + 2.0*(rhovshsq - eta_aniso*rhovsvsq)) )
- ! uses temporary templ2_cos from c14
+ ! uses temporary templ2_cos from c14
c16 = 0.5*cosphi*sinphi*sinthetasq* &
( cosphisq*templ2_cos &
+ 2.0*etaminone*sinphisq*(rhovphsq - two_rhovsvsq) )
@@ -662,18 +662,18 @@
+ sinphifour* &
(rhovphsq*costhetafour + 2.0*costhetasq*sinthetasq*(eta_aniso*rhovphsq &
+ two_rhovsvsq - two_eta_aniso*rhovsvsq) + rhovpvsq*sinthetafour)
-
+
! uses temporary templ1 from c13
c23 = 0.125*sinphisq*(rhovphsq + six_eta_aniso*rhovphsq &
+ rhovpvsq - four_rhovsvsq - 12.0*eta_aniso*rhovsvsq + cosfourtheta*templ1) &
+ cosphisq*(eta_aniso*costhetasq*(rhovphsq - two_rhovsvsq) + sinthetasq*(rhovphsq - two_rhovshsq))
-
- ! uses temporary templ1 from c13
+
+ ! uses temporary templ1 from c13
c24 = costheta*sinphi*sintheta* &
( etaminone*cosphisq*(rhovphsq - two_rhovsvsq) &
+ 0.5*sinphisq*(rhovpvsq - rhovphsq + costwotheta*templ1) )
- ! uses temporary templ2_cos from c14
+ ! uses temporary templ2_cos from c14
c25 = cosphi*costheta*sintheta* &
( cosphisq*(etaminone*rhovphsq + 2.0*(rhovshsq - eta_aniso*rhovsvsq)) &
+ 0.5*sinphisq*(templ2_cos + four_rhovshsq - four_rhovsvsq) )
@@ -681,39 +681,39 @@
! uses temporary templ2_cos from c14
c26 = 0.5*cosphi*sinphi*sinthetasq* &
( 2.0*etaminone*cosphisq*(rhovphsq - two_rhovsvsq) &
- + sinphisq*templ2_cos )
+ + sinphisq*templ2_cos )
c33 = rhovpvsq*costhetafour &
+ 2.0*costhetasq*sinthetasq*(two_rhovsvsq + eta_aniso*(rhovphsq - two_rhovsvsq)) &
- + rhovphsq*sinthetafour
-
+ + rhovphsq*sinthetafour
+
! uses temporary templ1_cos from c13
- c34 = - 0.25*sinphi*sintwotheta*templ1_cos
+ c34 = - 0.25*sinphi*sintwotheta*templ1_cos
- ! uses temporary templ1_cos from c34
+ ! uses temporary templ1_cos from c34
c35 = - 0.25*cosphi*sintwotheta*templ1_cos
-
+
! uses temporary templ1_cos from c34
- c36 = - 0.25*sintwophi*sinthetasq*(templ1_cos - four_rhovshsq + four_rhovsvsq)
+ c36 = - 0.25*sintwophi*sinthetasq*(templ1_cos - four_rhovshsq + four_rhovsvsq)
c44 = cosphisq*(rhovsvsq*costhetasq + rhovshsq*sinthetasq) &
+ sinphisq*(rhovsvsq*costwothetasq + costhetasq*sinthetasq*templ3)
-
- ! uses temporary templ3 from c44
+
+ ! uses temporary templ3 from c44
c46 = - cosphi*costheta*sintheta* &
( cosphisq*(rhovshsq - rhovsvsq) - 0.5*sinphisq*templ3_cos )
- ! uses templ3 from c46
+ ! uses templ3 from c46
c45 = 0.25*sintwophi*sinthetasq* &
(templ3_two + costwotheta*(rhovphsq + rhovpvsq - two_eta_aniso*rhovphsq + 4.0*etaminone*rhovsvsq))
-
+
c55 = sinphisq*(rhovsvsq*costhetasq + rhovshsq*sinthetasq) &
+ cosphisq*(rhovsvsq*costwothetasq &
+ costhetasq*sinthetasq*(rhovphsq - two_eta_aniso*rhovphsq + rhovpvsq + four_eta_aniso*rhovsvsq) )
-
+
! uses temporary templ3_cos from c46
c56 = costheta*sinphi*sintheta* &
- ( 0.5*cosphisq*templ3_cos + sinphisq*(rhovsvsq - rhovshsq) )
+ ( 0.5*cosphisq*templ3_cos + sinphisq*(rhovsvsq - rhovshsq) )
c66 = rhovshsq*costwophisq*costhetasq &
- 2.0*cosphisq*costhetasq*sinphisq*(rhovphsq - two_rhovshsq) &
@@ -722,8 +722,8 @@
+ cos(4.0*phi - 2.0*theta) + cos(2.0*(2.0*phi + theta)) ) &
+ rhovpvsq*cosphisq*sinphisq*sinthetafour &
- 0.5*eta_aniso*sintwophisq*sinthetafour*(rhovphsq - two_rhovsvsq)
-
+
! general expression of stress tensor for full Cijkl with 21 coefficients
sigma_xx = c11*duxdxl + c16*duxdyl_plus_duydxl + c12*duydyl + &
c15*duzdxl_plus_duxdzl + c14*duzdyl_plus_duydzl + c13*duzdzl
@@ -745,11 +745,11 @@
! subtract memory variables if attenuation
if(ATTENUATION_VAL .and. ( USE_ATTENUATION_MIMIC .eqv. .false. ) ) then
-
+
! note: fortran passes pointers to array location, thus R_memory(1,1,...) should be fine
call compute_element_att_stress( R_memory(1,1,i,j,k,ispec), &
sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz)
-
+
endif ! ATTENUATION_VAL
! define symmetric components of sigma for gravity
@@ -879,10 +879,10 @@
enddo ! NGLLX
enddo ! NGLLY
enddo ! NGLLZ
-
-
+
+
end subroutine compute_element_tiso
-
+
!
!--------------------------------------------------------------------------------------------
!
@@ -901,8 +901,8 @@
one_minus_sum_beta,vx,vy,vz,vnspec, &
tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,rho_s_H)
-
+
! this routine is optimized for NGLLX = NGLLY = NGLLZ = 5 using the Deville et al. (2002) inlined matrix-matrix products
implicit none
@@ -984,7 +984,7 @@
integer :: i,j,k
integer :: int_radius
integer :: iglob1
-
+
! anisotropic elements
do k=1,NGLLZ
@@ -1239,8 +1239,8 @@
tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl)
enddo ! NGLLX
enddo ! NGLLY
- enddo ! NGLLZ
-
+ enddo ! NGLLZ
+
end subroutine compute_element_aniso
!
@@ -1250,7 +1250,7 @@
subroutine compute_element_att_stress( R_memory_loc, &
sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz)
-
+
implicit none
include "constants.h"
@@ -1274,7 +1274,7 @@
integer :: imodulo_N_SLS
integer :: i_SLS1,i_SLS2
#endif
-
+
#ifdef _HANDOPT
! way 2:
! note: this should help compilers to pipeline the code and make better use of the cache;
@@ -1284,50 +1284,50 @@
if(imodulo_N_SLS >= 1) then
do i_SLS = 1,imodulo_N_SLS
- R_xx_val1 = R_memory_loc(1,i_SLS)
- R_yy_val1 = R_memory_loc(2,i_SLS)
+ R_xx_val1 = R_memory_loc(1,i_SLS)
+ R_yy_val1 = R_memory_loc(2,i_SLS)
sigma_xx = sigma_xx - R_xx_val1
sigma_yy = sigma_yy - R_yy_val1
sigma_zz = sigma_zz + R_xx_val1 + R_yy_val1
- sigma_xy = sigma_xy - R_memory_loc(3,i_SLS)
- sigma_xz = sigma_xz - R_memory_loc(4,i_SLS)
- sigma_yz = sigma_yz - R_memory_loc(5,i_SLS)
+ sigma_xy = sigma_xy - R_memory_loc(3,i_SLS)
+ sigma_xz = sigma_xz - R_memory_loc(4,i_SLS)
+ sigma_yz = sigma_yz - R_memory_loc(5,i_SLS)
enddo
endif
if(N_SLS >= imodulo_N_SLS+1) then
- ! note: another possibility would be using a reduction example for this loop; was tested but it does not improve,
- ! probably since N_SLS == 3 is too small for a loop benefit
+ ! note: another possibility would be using a reduction example for this loop; was tested but it does not improve,
+ ! probably since N_SLS == 3 is too small for a loop benefit
do i_SLS = imodulo_N_SLS+1,N_SLS,3
- R_xx_val1 = R_memory_loc(1,i_SLS)
- R_yy_val1 = R_memory_loc(2,i_SLS)
+ R_xx_val1 = R_memory_loc(1,i_SLS)
+ R_yy_val1 = R_memory_loc(2,i_SLS)
sigma_xx = sigma_xx - R_xx_val1
sigma_yy = sigma_yy - R_yy_val1
sigma_zz = sigma_zz + R_xx_val1 + R_yy_val1
- sigma_xy = sigma_xy - R_memory_loc(3,i_SLS)
- sigma_xz = sigma_xz - R_memory_loc(4,i_SLS)
- sigma_yz = sigma_yz - R_memory_loc(5,i_SLS)
+ sigma_xy = sigma_xy - R_memory_loc(3,i_SLS)
+ sigma_xz = sigma_xz - R_memory_loc(4,i_SLS)
+ sigma_yz = sigma_yz - R_memory_loc(5,i_SLS)
i_SLS1=i_SLS+1
- R_xx_val2 = R_memory_loc(1,i_SLS1)
- R_yy_val2 = R_memory_loc(2,i_SLS1)
+ R_xx_val2 = R_memory_loc(1,i_SLS1)
+ R_yy_val2 = R_memory_loc(2,i_SLS1)
sigma_xx = sigma_xx - R_xx_val2
sigma_yy = sigma_yy - R_yy_val2
sigma_zz = sigma_zz + R_xx_val2 + R_yy_val2
- sigma_xy = sigma_xy - R_memory_loc(3,i_SLS1)
- sigma_xz = sigma_xz - R_memory_loc(4,i_SLS1)
- sigma_yz = sigma_yz - R_memory_loc(5,i_SLS1)
+ sigma_xy = sigma_xy - R_memory_loc(3,i_SLS1)
+ sigma_xz = sigma_xz - R_memory_loc(4,i_SLS1)
+ sigma_yz = sigma_yz - R_memory_loc(5,i_SLS1)
i_SLS2 =i_SLS+2
- R_xx_val3 = R_memory_loc(1,i_SLS2)
- R_yy_val3 = R_memory_loc(2,i_SLS2)
+ R_xx_val3 = R_memory_loc(1,i_SLS2)
+ R_yy_val3 = R_memory_loc(2,i_SLS2)
sigma_xx = sigma_xx - R_xx_val3
sigma_yy = sigma_yy - R_yy_val3
sigma_zz = sigma_zz + R_xx_val3 + R_yy_val3
- sigma_xy = sigma_xy - R_memory_loc(3,i_SLS2)
- sigma_xz = sigma_xz - R_memory_loc(4,i_SLS2)
- sigma_yz = sigma_yz - R_memory_loc(5,i_SLS2)
- enddo
- endif
+ sigma_xy = sigma_xy - R_memory_loc(3,i_SLS2)
+ sigma_xz = sigma_xz - R_memory_loc(4,i_SLS2)
+ sigma_yz = sigma_yz - R_memory_loc(5,i_SLS2)
+ enddo
+ endif
#else
! way 1:
do i_SLS = 1,N_SLS
@@ -1354,7 +1354,7 @@
alphaval,betaval,gammaval, &
c44store,muvstore, &
epsilondev,epsilondev_loc)
-! crust mantle
+! crust mantle
! update memory variables based upon the Runge-Kutta scheme
! convention for attenuation
@@ -1381,7 +1381,7 @@
! element id
integer :: ispec
-
+
real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: R_memory
integer :: vx,vy,vz,vnspec
@@ -1390,14 +1390,14 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPECMAX_ANISO_MANTLE) :: c44store
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPECMAX_ISO_MANTLE) :: muvstore
-
+
real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: epsilondev
real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc
-! local parameters
+! local parameters
real(kind=CUSTOM_REAL), dimension(NGLLX, NGLLY, NGLLZ) :: factor_common_c44_muv
integer :: i_SLS
-
+
#ifdef _HANDOPT
real(kind=CUSTOM_REAL) :: alphal,betal,gammal
integer :: i,j,k
@@ -1412,12 +1412,12 @@
! IMPROVE we should probably use an average value instead
#ifdef _HANDOPT
-! way 2:
+! way 2:
do i_SLS = 1,N_SLS
alphal = alphaval(i_SLS)
betal = betaval(i_SLS)
- gammal = gammaval(i_SLS)
+ gammal = gammaval(i_SLS)
! reformatted R_memory to handle large factor_common and reduced [alpha,beta,gamma]val
factor_common_c44_muv(:,:,:) = factor_common(i_SLS,:,:,:,ispec)
@@ -1451,7 +1451,7 @@
factor_common_c44_muv(:,:,:) = factor_common_c44_muv(:,:,:) * muvstore(:,:,:,ispec)
endif
- do i_memory = 1,5
+ do i_memory = 1,5
R_memory(i_memory,i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_memory(i_memory,i_SLS,:,:,:,ispec) &
+ factor_common_c44_muv(:,:,:) &
* (betaval(i_SLS) * epsilondev(i_memory,:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(i_memory,:,:,:))
@@ -1461,7 +1461,7 @@
end subroutine compute_element_att_memory_cr
-
+
!
!--------------------------------------------------------------------------------------------
!
@@ -1471,7 +1471,7 @@
alphaval,betaval,gammaval, &
muvstore, &
epsilondev,epsilondev_loc)
-! inner core
+! inner core
! update memory variables based upon the Runge-Kutta scheme
! convention for attenuation
@@ -1498,7 +1498,7 @@
! element id
integer :: ispec
-
+
real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ATTENUATION) :: R_memory
integer :: vx,vy,vz,vnspec
@@ -1506,15 +1506,15 @@
real(kind=CUSTOM_REAL), dimension(N_SLS) :: alphaval,betaval,gammaval
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: muvstore
-
- real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: epsilondev
+
+ real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: epsilondev
real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc
-! local parameters
+! local parameters
real(kind=CUSTOM_REAL), dimension(NGLLX, NGLLY, NGLLZ) :: factor_common_use
-
+
integer :: i_SLS
-
+
#ifdef _HANDOPT
real(kind=CUSTOM_REAL) :: alphal,betal,gammal
integer :: i,j,k
@@ -1534,11 +1534,11 @@
alphal = alphaval(i_SLS)
betal = betaval(i_SLS)
- gammal = gammaval(i_SLS)
+ gammal = gammaval(i_SLS)
! reformatted R_memory to handle large factor_common and reduced [alpha,beta,gamma]val
- factor_common_use(:,:,:) = factor_common(i_SLS,:,:,:,ispec)
-
+ factor_common_use(:,:,:) = factor_common(i_SLS,:,:,:,ispec)
+
factor_common_use(:,:,:) = factor_common_use(:,:,:) * muvstore(:,:,:,ispec)
! this helps to vectorize the inner most loop
@@ -1551,7 +1551,7 @@
enddo
enddo
enddo
-
+
enddo ! i_SLS
#else
! way 1:
@@ -1566,5 +1566,5 @@
#endif
end subroutine compute_element_att_memory_ic
-
-
\ No newline at end of file
+
+
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_crust_mantle_Dev.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_crust_mantle_Dev.F90 2011-11-18 23:19:05 UTC (rev 19218)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_crust_mantle_Dev.F90 2011-11-20 00:12:05 UTC (rev 19219)
@@ -29,7 +29,7 @@
! #undef _HANDOPT : turns hand-optimized code off
! or compile with: -D_HANDOPT
!#define _HANDOPT
-
+
! note: these hand optimizations should help compilers to pipeline the code and make better use of the cache;
! depending on compilers, it can further decrease the computation time by ~ 30%.
! the original routines are commented with "! way 1", the hand-optimized routines with "! way 2"
@@ -169,8 +169,8 @@
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NGLLZ) :: sum_terms
- real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc
- real(kind=CUSTOM_REAL) fac1,fac2,fac3
+ real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc
+ real(kind=CUSTOM_REAL) fac1,fac2,fac3
! for gravity
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NGLLZ) :: rho_s_H
@@ -297,7 +297,7 @@
! way 2:
! since we know that NGLLX = 5, this should help pipelining
iglobv5(:) = ibool(:,j,k,ispec)
-
+
dummyx_loc(1,j,k) = displ_crust_mantle(1,iglobv5(1))
dummyy_loc(1,j,k) = displ_crust_mantle(2,iglobv5(1))
dummyz_loc(1,j,k) = displ_crust_mantle(3,iglobv5(1))
@@ -317,7 +317,7 @@
dummyx_loc(5,j,k) = displ_crust_mantle(1,iglobv5(5))
dummyy_loc(5,j,k) = displ_crust_mantle(2,iglobv5(5))
dummyz_loc(5,j,k) = displ_crust_mantle(3,iglobv5(5))
-
+
#else
! way 1:
do i=1,NGLLX
@@ -413,22 +413,22 @@
R_memory,epsilon_trace_over_3, &
one_minus_sum_beta,vx,vy,vz,vnspec, &
tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
- dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,rho_s_H)
- else
+ dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,rho_s_H)
+ else
if( .not. ispec_is_tiso(ispec) ) then
! isotropic element
call compute_element_iso(ispec, &
minus_gravity_table,density_table,minus_deriv_gravity_table, &
- xstore,ystore,zstore, &
+ xstore,ystore,zstore, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
wgll_cube, &
kappavstore,muvstore, &
ibool, &
R_memory,epsilon_trace_over_3, &
- one_minus_sum_beta,vx,vy,vz,vnspec, &
+ one_minus_sum_beta,vx,vy,vz,vnspec, &
tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
- dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,rho_s_H)
- else
+ dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,rho_s_H)
+ else
! transverse isotropic element
call compute_element_tiso(ispec, &
minus_gravity_table,density_table,minus_deriv_gravity_table, &
@@ -440,10 +440,10 @@
R_memory,epsilon_trace_over_3, &
one_minus_sum_beta,vx,vy,vz,vnspec, &
tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
- dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,rho_s_H)
- endif ! .not. ispec_is_tiso
- endif
-
+ dummyx_loc,dummyy_loc,dummyz_loc,epsilondev_loc,rho_s_H)
+ endif ! .not. ispec_is_tiso
+ endif
+
! subroutines adapted from Deville, Fischer and Mund, High-order methods
! for incompressible fluid flow, Cambridge University Press (2002),
! pages 386 and 389 and Figure 8.3.1
@@ -526,7 +526,7 @@
sum_terms(2,i,j,k) = - (fac1*newtempy1(i,j,k) + fac2*newtempy2(i,j,k) + fac3*newtempy3(i,j,k))
sum_terms(3,i,j,k) = - (fac1*newtempz1(i,j,k) + fac2*newtempz2(i,j,k) + fac3*newtempz3(i,j,k))
- if(GRAVITY_VAL) sum_terms(:,i,j,k) = sum_terms(:,i,j,k) + rho_s_H(:,i,j,k)
+ if(GRAVITY_VAL) sum_terms(:,i,j,k) = sum_terms(:,i,j,k) + rho_s_H(:,i,j,k)
enddo ! NGLLX
@@ -538,7 +538,7 @@
do j=1,NGLLY
#ifdef _HANDOPT
-! way 2:
+! way 2:
iglobv5(:) = ibool(:,j,k,ispec)
accel_crust_mantle(:,iglobv5(1)) = accel_crust_mantle(:,iglobv5(1)) + sum_terms(:,1,j,k)
@@ -546,8 +546,8 @@
accel_crust_mantle(:,iglobv5(3)) = accel_crust_mantle(:,iglobv5(3)) + sum_terms(:,3,j,k)
accel_crust_mantle(:,iglobv5(4)) = accel_crust_mantle(:,iglobv5(4)) + sum_terms(:,4,j,k)
accel_crust_mantle(:,iglobv5(5)) = accel_crust_mantle(:,iglobv5(5)) + sum_terms(:,5,j,k)
-
-#else
+
+#else
! way 1:
do i=1,NGLLX
iglob1 = ibool(i,j,k,ispec)
@@ -580,7 +580,7 @@
alphaval,betaval,gammaval, &
c44store,muvstore, &
epsilondev,epsilondev_loc)
-
+
endif
! save deviatoric strain for Runge-Kutta scheme
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_inner_core_Dev.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_inner_core_Dev.F90 2011-11-18 23:19:05 UTC (rev 19218)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_inner_core_Dev.F90 2011-11-20 00:12:05 UTC (rev 19219)
@@ -29,7 +29,7 @@
! #undef _HANDOPT : turns hand-optimized code off
! or compile with: -D_HANDOPT
!#define _HANDOPT
-
+
! note: these hand optimizations should help compilers to pipeline the code and make better use of the cache;
! depending on compilers, it can further decrease the computation time by ~ 30%.
! the original routines are commented with "! way 1", the hand-optimized routines with "! way 2"
@@ -300,14 +300,14 @@
! subroutines adapted from Deville, Fischer and Mund, High-order methods
! for incompressible fluid flow, Cambridge University Press (2002),
! pages 386 and 389 and Figure 8.3.1
-
+
do k=1,NGLLZ
do j=1,NGLLY
-#ifdef _HANDOPT
+#ifdef _HANDOPT
! way 2:
- ! since we know that NGLLX = 5, this should help pipelining
+ ! since we know that NGLLX = 5, this should help pipelining
iglobv5(:) = ibool(:,j,k,ispec)
-
+
dummyx_loc(1,j,k) = displ_inner_core(1,iglobv5(1))
dummyy_loc(1,j,k) = displ_inner_core(2,iglobv5(1))
dummyz_loc(1,j,k) = displ_inner_core(3,iglobv5(1))
@@ -327,7 +327,7 @@
dummyx_loc(5,j,k) = displ_inner_core(1,iglobv5(5))
dummyy_loc(5,j,k) = displ_inner_core(2,iglobv5(5))
dummyz_loc(5,j,k) = displ_inner_core(3,iglobv5(5))
-
+
#else
! way 1:
do i=1,NGLLX
@@ -339,7 +339,7 @@
#endif
enddo
enddo
-
+
do j=1,m2
do i=1,m1
C1_m1_m2_5points(i,j) = hprime_xx(i,1)*B1_m1_m2_5points(1,j) + &
@@ -765,7 +765,7 @@
! sum contributions from each element to the global mesh and add gravity terms
do k=1,NGLLZ
do j=1,NGLLY
-#ifdef _HANDOPT
+#ifdef _HANDOPT
! way 2:
iglobv5(:) = ibool(:,j,k,ispec)
@@ -774,8 +774,8 @@
accel_inner_core(:,iglobv5(3)) = accel_inner_core(:,iglobv5(3)) + sum_terms(:,3,j,k)
accel_inner_core(:,iglobv5(4)) = accel_inner_core(:,iglobv5(4)) + sum_terms(:,4,j,k)
accel_inner_core(:,iglobv5(5)) = accel_inner_core(:,iglobv5(5)) + sum_terms(:,5,j,k)
-
-#else
+
+#else
! way 1:
do i=1,NGLLX
iglob1 = ibool(i,j,k,ispec)
@@ -800,14 +800,14 @@
! therefore Q_\alpha is not zero; for instance for V_p / V_s = sqrt(3)
! we get Q_\alpha = (9 / 4) * Q_\mu = 2.25 * Q_\mu
if(ATTENUATION_VAL .and. ( USE_ATTENUATION_MIMIC .eqv. .false. ) ) then
-
+
! updates R_memory
call compute_element_att_memory_ic(ispec,R_memory, &
vx,vy,vz,vnspec,factor_common, &
alphaval,betaval,gammaval, &
muvstore, &
epsilondev,epsilondev_loc)
-
+
endif
! save deviatoric strain for Runge-Kutta scheme
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/noise_tomography.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/noise_tomography.f90 2011-11-18 23:19:05 UTC (rev 19218)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/noise_tomography.f90 2011-11-20 00:12:05 UTC (rev 19219)
@@ -298,7 +298,7 @@
! total file size
filesize = reclen
filesize = filesize*NSTEP
-
+
write(outputname,"('/proc',i6.6,'_surface_movie')") myrank
if (NOISE_TOMOGRAPHY==1) call open_file_abs_w(9,trim(LOCAL_PATH)//trim(outputname), &
len_trim(trim(LOCAL_PATH)//trim(outputname)), &
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/read_mesh_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/read_mesh_databases.f90 2011-11-18 23:19:05 UTC (rev 19218)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/read_mesh_databases.f90 2011-11-20 00:12:05 UTC (rev 19219)
@@ -809,7 +809,7 @@
integer NSTEP
! local parameters
- integer(kind=8) :: filesize
+ integer(kind=8) :: filesize
character(len=150) prname
@@ -830,14 +830,14 @@
if (nspec2D_xmin_crust_mantle > 0 .and. (SIMULATION_TYPE == 3 &
.or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))) then
-
+
! size of single record
reclen_xmin_crust_mantle = CUSTOM_REAL * (NDIM * NGLLY * NGLLZ * nspec2D_xmin_crust_mantle)
! total file size
filesize = reclen_xmin_crust_mantle
filesize = filesize*NSTEP
-
+
if (SIMULATION_TYPE == 3) then
! open(unit=51,file=trim(prname)//'absorb_xmin.bin', &
! status='old',action='read',form='unformatted',access='direct', &
@@ -857,14 +857,14 @@
if (nspec2D_xmax_crust_mantle > 0 .and. (SIMULATION_TYPE == 3 &
.or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))) then
-
+
! size of single record
reclen_xmax_crust_mantle = CUSTOM_REAL * (NDIM * NGLLY * NGLLZ * nspec2D_xmax_crust_mantle)
-
+
! total file size
filesize = reclen_xmax_crust_mantle
filesize = filesize*NSTEP
-
+
if (SIMULATION_TYPE == 3) then
! open(unit=52,file=trim(prname)//'absorb_xmax.bin', &
! status='old',action='read',form='unformatted',access='direct', &
@@ -884,7 +884,7 @@
if (nspec2D_ymin_crust_mantle > 0 .and. (SIMULATION_TYPE == 3 &
.or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))) then
-
+
! size of single record
reclen_ymin_crust_mantle = CUSTOM_REAL * (NDIM * NGLLX * NGLLZ * nspec2D_ymin_crust_mantle)
@@ -892,7 +892,7 @@
filesize = reclen_ymin_crust_mantle
filesize = filesize*NSTEP
-
+
if (SIMULATION_TYPE == 3) then
! open(unit=53,file=trim(prname)//'absorb_ymin.bin', &
! status='old',action='read',form='unformatted',access='direct',&
@@ -912,14 +912,14 @@
if (nspec2D_ymax_crust_mantle > 0 .and. (SIMULATION_TYPE == 3 &
.or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))) then
-
+
! size of single record
reclen_ymax_crust_mantle = CUSTOM_REAL * (NDIM * NGLLX * NGLLZ * nspec2D_ymax_crust_mantle)
! total file size
filesize = reclen_ymax_crust_mantle
filesize = filesize*NSTEP
-
+
if (SIMULATION_TYPE == 3) then
! open(unit=54,file=trim(prname)//'absorb_ymax.bin', &
! status='old',action='read',form='unformatted',access='direct',&
@@ -955,14 +955,14 @@
if (nspec2D_xmin_outer_core > 0 .and. (SIMULATION_TYPE == 3 &
.or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))) then
-
+
! size of single record
reclen_xmin_outer_core = CUSTOM_REAL * (NGLLY * NGLLZ * nspec2D_xmin_outer_core)
-
+
! total file size
filesize = reclen_xmin_outer_core
filesize = filesize*NSTEP
-
+
if (SIMULATION_TYPE == 3) then
! open(unit=61,file=trim(prname)//'absorb_xmin.bin', &
! status='old',action='read',form='unformatted',access='direct', &
@@ -982,14 +982,14 @@
if (nspec2D_xmax_outer_core > 0 .and. (SIMULATION_TYPE == 3 &
.or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))) then
-
+
! size of single record
reclen_xmax_outer_core = CUSTOM_REAL * (NGLLY * NGLLZ * nspec2D_xmax_outer_core)
! total file size
filesize = reclen_xmax_outer_core
filesize = filesize*NSTEP
-
+
if (SIMULATION_TYPE == 3) then
! open(unit=62,file=trim(prname)//'absorb_xmax.bin', &
! status='old',action='read',form='unformatted',access='direct', &
@@ -1010,14 +1010,14 @@
if (nspec2D_ymin_outer_core > 0 .and. (SIMULATION_TYPE == 3 &
.or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))) then
-
+
! size of single record
reclen_ymin_outer_core = CUSTOM_REAL * (NGLLX * NGLLZ * nspec2D_ymin_outer_core)
! total file size
filesize = reclen_ymin_outer_core
filesize = filesize*NSTEP
-
+
if (SIMULATION_TYPE == 3) then
! open(unit=63,file=trim(prname)//'absorb_ymin.bin', &
! status='old',action='read',form='unformatted',access='direct',&
@@ -1038,14 +1038,14 @@
if (nspec2D_ymax_outer_core > 0 .and. (SIMULATION_TYPE == 3 &
.or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))) then
-
+
! size of single record
reclen_ymax_outer_core = CUSTOM_REAL * (NGLLX * NGLLZ * nspec2D_ymax_outer_core)
! total file size
filesize = reclen_ymax_outer_core
filesize = filesize*NSTEP
-
+
if (SIMULATION_TYPE == 3) then
! open(unit=64,file=trim(prname)//'absorb_ymax.bin', &
! status='old',action='read',form='unformatted',access='direct',&
@@ -1066,14 +1066,14 @@
if (NSPEC2D_BOTTOM(IREGION_OUTER_CORE) > 0 .and. &
(SIMULATION_TYPE == 3 .or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD)))then
-
+
! size of single record
reclen_zmin = CUSTOM_REAL * (NGLLX * NGLLY * NSPEC2D_BOTTOM(IREGION_OUTER_CORE))
! total file size
filesize = reclen_zmin
filesize = filesize*NSTEP
-
+
if (SIMULATION_TYPE == 3) then
! open(unit=65,file=trim(prname)//'absorb_zmin.bin', &
! status='old',action='read',form='unformatted',access='direct',&
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/setup_sources_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/setup_sources_receivers.f90 2011-11-18 23:19:05 UTC (rev 19218)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/setup_sources_receivers.f90 2011-11-20 00:12:05 UTC (rev 19219)
@@ -111,7 +111,7 @@
character(len=2) :: bic
! makes smaller hdur for movies
logical,parameter :: USE_SMALLER_HDUR_MOVIE = .false.
-
+
! sources
! BS BS moved open statement and writing of first lines into sr.vtk before the
! call to locate_sources, where further write statements to that file follow
@@ -141,14 +141,14 @@
if(abs(minval(tshift_cmt)) > TINYVAL) call exit_MPI(myrank,'one tshift_cmt must be zero, others must be positive')
! filter source time function by Gaussian with hdur = HDUR_MOVIE when outputing movies or shakemaps
- if (MOVIE_SURFACE .or. MOVIE_VOLUME ) then
+ if (MOVIE_SURFACE .or. MOVIE_VOLUME ) then
! smaller hdur_movie will do
if( USE_SMALLER_HDUR_MOVIE ) then
! hdur_movie gets assigned an automatic value based on the simulation resolution
! this will make that a bit smaller to have a higher-frequency movie output
HDUR_MOVIE = 0.5* HDUR_MOVIE
endif
-
+
! new hdur for simulation
hdur = sqrt(hdur**2 + HDUR_MOVIE**2)
if(myrank == 0) then
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.F90 2011-11-18 23:19:05 UTC (rev 19218)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.F90 2011-11-20 00:12:05 UTC (rev 19219)
@@ -31,7 +31,7 @@
! #undef _HANDOPT : turns hand-optimized code off
! or compile with: -D_HANDOPT
!#define _HANDOPT
-
+
! note: these hand optimizations should help compilers to pipeline the code and make better use of the cache;
! depending on compilers, it can further decrease the computation time by ~ 30%.
! the original routines are commented with "! way 1", the hand-optimized routines with "! way 2"
@@ -888,7 +888,7 @@
#ifdef _HANDOPT
integer :: imodulo_NGLOB_CRUST_MANTLE,imodulo_NGLOB_CRUST_MANTLE4, &
- imodulo_NGLOB_INNER_CORE
+ imodulo_NGLOB_INNER_CORE
#endif
! ************** PROGRAM STARTS HERE **************
@@ -942,7 +942,7 @@
! For most compilers and hardware, will result in a significant speedup (> 30% or more, sometimes twice faster).
!
! note 5: a common technique to help compilers enhance pipelining is loop unrolling. We do this here in a simple
-! and straigthforward way, so don't be confused about the do-loop writing. For this to take effect,
+! and straigthforward way, so don't be confused about the do-loop writing. For this to take effect,
! you have to turn the hand-optimization flag on: compile with additional flag -D_HANDOPT
!
! note 6: whenever adding some new code, please make sure to use
@@ -1946,7 +1946,7 @@
! update position in seismograms
seismo_current = seismo_current + 1
- ! Newark time scheme update
+ ! Newark time scheme update
#ifdef _HANDOPT
! way 2:
! One common technique in computational science to help enhance pipelining is loop unrolling
@@ -1994,7 +1994,7 @@
accel_crust_mantle(:,i+2) = 0._CUSTOM_REAL
enddo
- ! outer core
+ ! outer core
do i=1,NGLOB_OUTER_CORE
displ_outer_core(i) = displ_outer_core(i) &
+ deltat*veloc_outer_core(i) + deltatsqover2*accel_outer_core(i)
@@ -2036,8 +2036,8 @@
accel_inner_core(:,i) = 0._CUSTOM_REAL
accel_inner_core(:,i+1) = 0._CUSTOM_REAL
accel_inner_core(:,i+2) = 0._CUSTOM_REAL
- enddo
-
+ enddo
+
#else
! way 1:
! mantle
@@ -2114,7 +2114,7 @@
b_accel_outer_core(i) = 0._CUSTOM_REAL
enddo
- ! inner core
+ ! inner core
if(imodulo_NGLOB_INNER_CORE >= 1) then
do i=1,imodulo_NGLOB_INNER_CORE
b_displ_inner_core(:,i) = b_displ_inner_core(:,i) &
@@ -2143,7 +2143,7 @@
b_accel_inner_core(:,i+1) = 0._CUSTOM_REAL
b_accel_inner_core(:,i+2) = 0._CUSTOM_REAL
enddo
-#else
+#else
! way 1:
! mantle
do i=1,NGLOB_CRUST_MANTLE
@@ -3753,7 +3753,7 @@
veloc_crust_mantle(:,i+3) = veloc_crust_mantle(:,i+3) + deltatover2*accel_crust_mantle(:,i+3)
enddo
- ! inner core
+ ! inner core
if(imodulo_NGLOB_INNER_CORE >= 1) then
do i=1,imodulo_NGLOB_INNER_CORE
accel_inner_core(1,i) = accel_inner_core(1,i)*rmass_inner_core(i) &
@@ -3857,7 +3857,7 @@
b_veloc_inner_core(:,i+2) = b_veloc_inner_core(:,i+2) + b_deltatover2*b_accel_inner_core(:,i+2)
enddo
-#else
+#else
! way 1:
! mantle
do i=1,NGLOB_CRUST_MANTLE
@@ -4295,7 +4295,7 @@
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, &
+ accel_crust_mantle,accel_inner_core, &
ibool_crust_mantle,ibool_inner_core)
else if (MOVIE_VOLUME_TYPE == 5) then !output displacement
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_movie_volume.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_movie_volume.f90 2011-11-18 23:19:05 UTC (rev 19218)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_movie_volume.f90 2011-11-20 00:12:05 UTC (rev 19219)
@@ -478,11 +478,11 @@
veloc_crust_mantle,veloc_inner_core,veloc_outer_core, &
accel_crust_mantle,accel_inner_core, &
ibool_crust_mantle,ibool_inner_core)
-
+
implicit none
include "constants.h"
include "OUTPUT_FILES/values_from_mesher.h"
-
+
integer :: myrank,it
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: eps_trace_over_3_crust_mantle
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_OUTER_CORE_ADJOINT) :: div_displ_outer_core
@@ -492,7 +492,7 @@
rhostore_outer_core,kappavstore_outer_core
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_OUTER_CORE) :: ibool_outer_core
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_STRAIN_ONLY) :: eps_trace_over_3_inner_core
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_STRAIN_ONLY) :: eps_trace_over_3_inner_core
real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: epsilondev_crust_mantle
real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_STR_OR_ATT) :: epsilondev_inner_core
@@ -517,26 +517,26 @@
character(len=150) outputname
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: tmp_data
real(kind=CUSTOM_REAL) :: scale_displ,scale_veloc,scale_accel
-
+
! output parameters
logical,parameter :: MOVIE_OUTPUT_DIV = .true. ! divergence
- logical,parameter :: MOVIE_OUTPUT_CURL = .false. ! curl
+ logical,parameter :: MOVIE_OUTPUT_CURL = .false. ! 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
! outputs divergence
- if( MOVIE_OUTPUT_DIV ) then
+ if( MOVIE_OUTPUT_DIV ) then
! crust_mantle region
! these binary arrays can be converted into mesh format using the utilitiy ./bin/xcombine_vol_data
! old name format: write(outputname,"('proc',i6.6,'_crust_mantle_div_displ_it',i6.6,'.bin')") myrank,it
write(outputname,"('proc',i6.6,'_reg1_div_displ_it',i6.6,'.bin')") myrank,it
open(unit=27,file=trim(LOCAL_PATH)//'/'//trim(outputname),status='unknown',form='unformatted',iostat=ier)
- if( ier /= 0 ) call exit_MPI(myrank,'error opening file '//trim(outputname))
+ if( ier /= 0 ) call exit_MPI(myrank,'error opening file '//trim(outputname))
write(27) eps_trace_over_3_crust_mantle
close(27)
-
+
! outer core
if (NSPEC_OUTER_CORE_ADJOINT /= 1 ) then
write(outputname,"('proc',i6.6,'_reg2_div_displ_it',i6.6,'.bin')") myrank,it
@@ -559,7 +559,7 @@
enddo
enddo
enddo
-
+
! old name format: write(outputname,"('proc',i6.6,'_outer_core_div_displ_it',i6.6,'.bin')") myrank,it
write(outputname,"('proc',i6.6,'_reg2_div_displ_it',i6.6,'.bin')") myrank,it
open(unit=27,file=trim(LOCAL_PATH)//'/'//trim(outputname),status='unknown',form='unformatted',iostat=ier)
@@ -569,7 +569,7 @@
deallocate(div_s_outer_core)
endif
-
+
! inner core
! old name format: write(outputname,"('proc',i6.6,'_inner_core_div_displ_proc_it',i6.6,'.bin')") myrank,it
write(outputname,"('proc',i6.6,'_reg3_div_displ_it',i6.6,'.bin')") myrank,it
@@ -595,10 +595,10 @@
write(27) epsilondev_inner_core
close(27)
endif
-
+
! outputs norm of epsilondev
if( MOVIE_OUTPUT_CURLNORM ) then
- ! these binary arrays can be converted into mesh format using the utilitiy ./bin/xcombine_vol_data
+ ! these binary arrays can be converted into mesh format using the utilitiy ./bin/xcombine_vol_data
! crust_mantle region
write(outputname,"('proc',i6.6,'_reg1_epsdev_displ_it',i6.6,'.bin')") myrank,it
open(unit=27,file=trim(LOCAL_PATH)//'/'//trim(outputname),status='unknown',form='unformatted',iostat=ier)
@@ -608,12 +608,12 @@
do ispec = 1, NSPEC_CRUST_MANTLE
do k = 1, NGLLZ
do j = 1, NGLLY
- do i = 1, NGLLX
+ do i = 1, NGLLX
tmp_data(i,j,k,ispec) = sqrt( epsilondev_crust_mantle(1,i,j,k,ispec)**2 &
+ epsilondev_crust_mantle(2,i,j,k,ispec)**2 &
+ epsilondev_crust_mantle(3,i,j,k,ispec)**2 &
+ epsilondev_crust_mantle(4,i,j,k,ispec)**2 &
- + epsilondev_crust_mantle(5,i,j,k,ispec)**2)
+ + epsilondev_crust_mantle(5,i,j,k,ispec)**2)
enddo
enddo
enddo
@@ -621,7 +621,7 @@
write(27) tmp_data
close(27)
deallocate(tmp_data)
-
+
! alternative: e.g. first component only
!write(27) epsilondev_crust_mantle(1,:,:,:,:)
!close(27)
@@ -635,12 +635,12 @@
do ispec = 1, NSPEC_INNER_CORE
do k = 1, NGLLZ
do j = 1, NGLLY
- do i = 1, NGLLX
+ do i = 1, NGLLX
tmp_data(i,j,k,ispec) = sqrt( epsilondev_inner_core(1,i,j,k,ispec)**2 &
+ epsilondev_inner_core(2,i,j,k,ispec)**2 &
+ epsilondev_inner_core(3,i,j,k,ispec)**2 &
+ epsilondev_inner_core(4,i,j,k,ispec)**2 &
- + epsilondev_inner_core(5,i,j,k,ispec)**2)
+ + epsilondev_inner_core(5,i,j,k,ispec)**2)
enddo
enddo
enddo
@@ -658,11 +658,11 @@
scale_displ = R_EARTH
scale_veloc = scale_displ * sqrt(PI*GRAV*RHOAV)
scale_accel = scale_veloc * dsqrt(PI*GRAV*RHOAV)
-
- ! outputs norm of displacement
+
+ ! outputs norm of displacement
if( MOVIE_OUTPUT_DISPLNORM ) then
! crust mantle
- ! these binary arrays can be converted into mesh format using the utilitiy ./bin/xcombine_vol_data
+ ! these binary arrays can be converted into mesh format using the utilitiy ./bin/xcombine_vol_data
allocate(tmp_data(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE))
do ispec = 1, NSPEC_CRUST_MANTLE
do k = 1, NGLLZ
@@ -691,8 +691,8 @@
do j = 1, NGLLY
do i = 1, NGLLX
iglob = ibool_outer_core(i,j,k,ispec)
- ! norm
- ! note: disp_outer_core is potential, this just outputs the potential,
+ ! norm
+ ! note: disp_outer_core is potential, this just outputs the potential,
! not the actual displacement u = grad(rho * Chi) / rho
tmp_data(i,j,k,ispec) = abs(displ_outer_core(iglob))
enddo
@@ -713,7 +713,7 @@
do j = 1, NGLLY
do i = 1, NGLLX
iglob = ibool_inner_core(i,j,k,ispec)
- ! norm
+ ! norm
tmp_data(i,j,k,ispec) = scale_displ * sqrt( displ_inner_core(1,iglob)**2 &
+ displ_inner_core(2,iglob)**2 &
+ displ_inner_core(3,iglob)**2 )
@@ -730,10 +730,10 @@
endif
- ! outputs norm of velocity
+ ! outputs norm of velocity
if( MOVIE_OUTPUT_VELOCNORM ) then
! crust mantle
- ! these binary arrays can be converted into mesh format using the utilitiy ./bin/xcombine_vol_data
+ ! these binary arrays can be converted into mesh format using the utilitiy ./bin/xcombine_vol_data
allocate(tmp_data(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE))
do ispec = 1, NSPEC_CRUST_MANTLE
do k = 1, NGLLZ
@@ -799,18 +799,18 @@
close(27)
deallocate(tmp_data)
endif
-
+
! outputs norm of acceleration
- if( MOVIE_OUTPUT_ACCELNORM ) then
+ if( MOVIE_OUTPUT_ACCELNORM ) then
! acceleration
- ! these binary arrays can be converted into mesh format using the utilitiy ./bin/xcombine_vol_data
+ ! these binary arrays can be converted into mesh format using the utilitiy ./bin/xcombine_vol_data
allocate(tmp_data(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE))
do ispec = 1, NSPEC_CRUST_MANTLE
do k = 1, NGLLZ
do j = 1, NGLLY
do i = 1, NGLLX
iglob = ibool_crust_mantle(i,j,k,ispec)
- ! norm
+ ! norm
tmp_data(i,j,k,ispec) = scale_accel * sqrt( accel_crust_mantle(1,iglob)**2 &
+ accel_crust_mantle(2,iglob)**2 &
+ accel_crust_mantle(3,iglob)**2 )
@@ -832,9 +832,9 @@
do j = 1, NGLLY
do i = 1, NGLLX
iglob = ibool_outer_core(i,j,k,ispec)
- ! norm
+ ! norm
! note: this outputs only the second time derivative of the potential,
- ! not the actual acceleration or pressure p = - rho * Chi_dot_dot
+ ! not the actual acceleration or pressure p = - rho * Chi_dot_dot
tmp_data(i,j,k,ispec) = abs(accel_outer_core(iglob))
enddo
enddo
More information about the CIG-COMMITS
mailing list