[cig-commits] r20682 - in seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER: EXAMPLES/noise_examples/global_long EXAMPLES/noise_examples/global_short EXAMPLES/noise_examples/regional EXAMPLES/noise_examples/test_global EXAMPLES/noise_examples/test_regional setup src/auxiliaries src/cuda src/meshfem3D src/shared src/specfem3D
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Tue Sep 4 13:41:20 PDT 2012
Author: danielpeter
Date: 2012-09-04 13:41:19 -0700 (Tue, 04 Sep 2012)
New Revision: 20682
Modified:
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/noise_examples/global_long/NOISE_collect
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/noise_examples/global_short/NOISE_collect
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/noise_examples/regional/NOISE_collect
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/noise_examples/test_global/NOISE_collect
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/noise_examples/test_regional/NOISE_collect
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/setup/constants.h.in
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/combine_surf_data.f90
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_vtk.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/check_fields_cuda.cu
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/model_s362ani.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/save_arrays_solver.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/setup_MPI_interfaces.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/read_compute_parameters.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_outer_core.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_outer_core_Dev.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_receivers.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_sources.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/prepare_timerun.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_arrays_solver.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_mesh_databases.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/setup_sources_receivers.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_movie_volume.f90
Log:
merges solver_data_1.bin, solver_data_2.bin and array_dims.txt into solver_data.bin; re-calculates NSTEP and adds t0 initial time to reach full record length; bug fix in rotation for outer core Deville routine
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/noise_examples/global_long/NOISE_collect
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/noise_examples/global_long/NOISE_collect 2012-09-04 18:44:13 UTC (rev 20681)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/noise_examples/global_long/NOISE_collect 2012-09-04 20:41:19 UTC (rev 20682)
@@ -6,6 +6,5 @@
##### other files can be collected in the same way
mv $DIR_LOCAL/*$iproc*reg1*kernel*.bin $DIR_GLOBAL/
-mv $DIR_LOCAL/*$iproc*reg1_*array_dims.txt $DIR_GLOBAL/
mv $DIR_LOCAL/*$iproc*reg1_solver_data*.bin $DIR_GLOBAL/
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/noise_examples/global_short/NOISE_collect
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/noise_examples/global_short/NOISE_collect 2012-09-04 18:44:13 UTC (rev 20681)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/noise_examples/global_short/NOISE_collect 2012-09-04 20:41:19 UTC (rev 20682)
@@ -6,6 +6,5 @@
##### other files can be collected in the same way
mv $DIR_LOCAL/*$iproc*reg1*kernel*.bin $DIR_GLOBAL/
-mv $DIR_LOCAL/*$iproc*reg1_*array_dims.txt $DIR_GLOBAL/
mv $DIR_LOCAL/*$iproc*reg1_solver_data*.bin $DIR_GLOBAL/
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/noise_examples/regional/NOISE_collect
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/noise_examples/regional/NOISE_collect 2012-09-04 18:44:13 UTC (rev 20681)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/noise_examples/regional/NOISE_collect 2012-09-04 20:41:19 UTC (rev 20682)
@@ -6,6 +6,5 @@
##### other files can be collected in the same way
mv $DIR_LOCAL/*$iproc*reg1*kernel*.bin $DIR_GLOBAL/
-mv $DIR_LOCAL/*$iproc*reg1_*array_dims.txt $DIR_GLOBAL/
mv $DIR_LOCAL/*$iproc*reg1_solver_data*.bin $DIR_GLOBAL/
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/noise_examples/test_global/NOISE_collect
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/noise_examples/test_global/NOISE_collect 2012-09-04 18:44:13 UTC (rev 20681)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/noise_examples/test_global/NOISE_collect 2012-09-04 20:41:19 UTC (rev 20682)
@@ -6,6 +6,5 @@
##### other files can be collected in the same way
mv $DIR_LOCAL/*$iproc*reg1*kernel*.bin $DIR_GLOBAL/
-mv $DIR_LOCAL/*$iproc*reg1_*array_dims.txt $DIR_GLOBAL/
mv $DIR_LOCAL/*$iproc*reg1_solver_data*.bin $DIR_GLOBAL/
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/noise_examples/test_regional/NOISE_collect
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/noise_examples/test_regional/NOISE_collect 2012-09-04 18:44:13 UTC (rev 20681)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/EXAMPLES/noise_examples/test_regional/NOISE_collect 2012-09-04 20:41:19 UTC (rev 20682)
@@ -6,6 +6,5 @@
##### other files can be collected in the same way
mv $DIR_LOCAL/*$iproc*reg1*kernel*.bin $DIR_GLOBAL/
-mv $DIR_LOCAL/*$iproc*reg1_*array_dims.txt $DIR_GLOBAL/
mv $DIR_LOCAL/*$iproc*reg1_solver_data*.bin $DIR_GLOBAL/
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/setup/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/setup/constants.h.in 2012-09-04 18:44:13 UTC (rev 20681)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/setup/constants.h.in 2012-09-04 20:41:19 UTC (rev 20682)
@@ -244,7 +244,7 @@
! using balancing algorithm
! postprocess the colors to balance them if Droux (1993) algorithm is not used
- logical, parameter :: BALANCE_COLORS_SIMPLE_ALGO = .false.
+ logical, parameter :: BALANCE_COLORS_SIMPLE_ALGO = .true.
!!-----------------------------------------------------------
!!
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/combine_surf_data.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/combine_surf_data.f90 2012-09-04 18:44:13 UTC (rev 20681)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/combine_surf_data.f90 2012-09-04 20:41:19 UTC (rev 20682)
@@ -134,8 +134,6 @@
FILE_ARRAY_IS_3D = .true.
endif
- dimen_name = trim(reg_name)//'array_dims.txt'
-
! figure out the total number of points/elements and allocate arrays
write(prname,'(a,i6.6,a)') trim(indir)//'/proc',node_list(1),'_'
nspec2D_file = trim(prname) // trim(belm_name)
@@ -169,14 +167,15 @@
print *, 'total number of elements = ', nelement_total
! ========= write points and elements files ===================
+ dimen_name = trim(reg_name)//'solver_data.bin'
allocate(ibelm_surf(nspec_surf))
do it = 1, num_node
write(prname,'(a,i6.6,a)') trim(indir)//'/proc',node_list(it),'_'
dimension_file = trim(prname) // trim(dimen_name)
- open(unit=27,file=trim(dimension_file),status='old',action='read', iostat = ios)
+ open(unit=27,file=trim(dimension_file),status='old',action='read', iostat = ios, form='unformatted')
if (ios /= 0) stop 'Error opening file'
- read(27,*) nspec(it)
- read(27,*) nglob(it)
+ read(27) nspec(it)
+ read(27) nglob(it)
close(27)
enddo
@@ -255,12 +254,14 @@
close(27)
! ibool file
- ibool_file = trim(prname2) // 'solver_data_2' // '.bin'
+ ibool_file = trim(prname2) // 'solver_data.bin'
print *, trim(ibool_file)
open(unit = 28,file = trim(ibool_file),status='old', iostat = ios, form='unformatted')
if (ios /= 0) then
print *,'Error opening ',trim(ibool_file); stop
endif
+ read(28) nspec(it)
+ read(28) nglob(it)
read(28) xstore(1:nglob(it))
read(28) ystore(1:nglob(it))
read(28) zstore(1:nglob(it))
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 2012-09-04 18:44:13 UTC (rev 20681)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/combine_vol_data.f90 2012-09-04 20:41:19 UTC (rev 20682)
@@ -62,7 +62,7 @@
! note:
! if one wants to remove the topography and ellipticity distortion, you would run the mesher again
! but turning the flags: TOPOGRAPHY and ELLIPTICITY to .false.
- ! then, use those as topo files: proc***_array_dims.txt and proc***_solver_data_2.bin
+ ! then, use those as topo files: proc***_solver_data.bin
! of course, this would also work by just turning ELLIPTICITY to .false. so that the CORRECT_ELLIPTICITY below
! becomes unneccessary
!
@@ -185,17 +185,17 @@
write(prname_file,'(a,i6.6,a,i1,a)') trim(in_file_dir)//'/proc',iproc,'_reg',ir,'_'
- dimension_file = trim(prname_topo) //'array_dims.txt'
- open(unit = 27,file = trim(dimension_file),status='old',action='read', iostat = ios)
+ dimension_file = trim(prname_topo) //'solver_data.bin'
+ open(unit = 27,file = trim(dimension_file),status='old',action='read', iostat = ios, form='unformatted')
if (ios /= 0) then
print*,'error ',ios
print*,'file:',trim(dimension_file)
stop 'Error opening file'
endif
+ read(27) nspec(it)
+ read(27) nglob(it)
+ close(27)
- read(27,*) nspec(it)
- read(27,*) nglob(it)
- close(27)
if (HIGH_RESOLUTION_MESH) then
npoint(it) = nglob(it)
nelement(it) = nspec(it) * (NGLLX-1) * (NGLLY-1) * (NGLLZ-1)
@@ -244,7 +244,7 @@
print *
! topology file
- topo_file = trim(prname_topo) // 'solver_data_2' // '.bin'
+ topo_file = trim(prname_topo) // 'solver_data.bin'
open(unit = 28,file = trim(topo_file),status='old',action='read', iostat = ios, form='unformatted')
if (ios /= 0) then
print*,'error ',ios
@@ -255,6 +255,8 @@
ystore(:) = 0.0
zstore(:) = 0.0
ibool(:,:,:,:) = -1
+ read(28) nspec(it)
+ read(28) nglob(it)
read(28) xstore(1:nglob(it))
read(28) ystore(1:nglob(it))
read(28) zstore(1:nglob(it))
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/combine_vol_data_vtk.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/combine_vol_data_vtk.f90 2012-09-04 18:44:13 UTC (rev 20681)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/auxiliaries/combine_vol_data_vtk.f90 2012-09-04 20:41:19 UTC (rev 20682)
@@ -67,7 +67,7 @@
! note:
! if one wants to remove the topography and ellipticity distortion, you would run the mesher again
! but turning the flags: TOPOGRAPHY and ELLIPTICITY to .false.
- ! then, use those as topo files: proc***_array_dims.txt and proc***_solver_data_2.bin
+ ! then, use those as topo files: proc***_solver_data.bin
! of course, this would also work by just turning ELLIPTICITY to .false. so that the CORRECT_ELLIPTICITY below
! becomes unneccessary
!
@@ -186,16 +186,15 @@
write(prname_file,'(a,i6.6,a,i1,a)') trim(in_file_dir)//'/proc',iproc,'_reg',ir,'_'
- dimension_file = trim(prname_topo) //'array_dims.txt'
- open(unit = 27,file = trim(dimension_file),status='old',action='read', iostat = ios)
+ dimension_file = trim(prname_topo) //'solver_data.bin'
+ open(unit = 28,file = trim(dimension_file),status='old',action='read', iostat = ios, form='unformatted')
if (ios /= 0) then
print*,'error ',ios
print*,'file:',trim(dimension_file)
stop 'Error opening file'
endif
-
- read(27,*) nspec(it)
- read(27,*) nglob(it)
+ read(27) nspec(it)
+ read(27) nglob(it)
close(27)
! check
@@ -286,7 +285,7 @@
print *
! topology file
- topo_file = trim(prname_topo) // 'solver_data_2' // '.bin'
+ topo_file = trim(prname_topo) // 'solver_data.bin'
open(unit = 28,file = trim(topo_file),status='old',action='read', iostat = ios, form='unformatted')
if (ios /= 0) then
print*,'error ',ios
@@ -297,6 +296,8 @@
ystore(:) = 0.0
zstore(:) = 0.0
ibool(:,:,:,:) = -1
+ read(28) nspec(it)
+ read(28) nglob(it)
read(28) xstore(1:nglob(it))
read(28) ystore(1:nglob(it))
read(28) zstore(1:nglob(it))
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/check_fields_cuda.cu
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/check_fields_cuda.cu 2012-09-04 18:44:13 UTC (rev 20681)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/check_fields_cuda.cu 2012-09-04 20:41:19 UTC (rev 20682)
@@ -119,11 +119,13 @@
}
// releases previous contexts
+/* crashes device
#if CUDA_VERSION < 4000
cudaThreadExit();
#else
cudaDeviceReset();
#endif
+*/
// stops program
//free(kernel_name);
@@ -191,11 +193,13 @@
}
// releases previous contexts
+/* crashes device
#if CUDA_VERSION < 4000
cudaThreadExit();
#else
cudaDeviceReset();
#endif
+*/
// stops program
#ifdef WITH_MPI
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu 2012-09-04 18:44:13 UTC (rev 20681)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/cuda/prepare_mesh_constants_cuda.cu 2012-09-04 20:41:19 UTC (rev 20682)
@@ -58,14 +58,14 @@
// we thus pass the address to the pointer above (as void double-pointer) to have it
// pointing to the correct pointer of the array here
print_CUDA_error_if_any(cudaMalloc((void**)d_array_addr_ptr,size*sizeof(int)),
- 10000);
+ 12001);
// copies values onto GPU
//
// note: cudaMemcpy uses the pointer to the array, we thus re-cast the value of
// the double-pointer above to have the correct pointer to the array
print_CUDA_error_if_any(cudaMemcpy((int*) *d_array_addr_ptr,h_array,size*sizeof(int),cudaMemcpyHostToDevice),
- 10001);
+ 12002);
}
/* ----------------------------------------------------------------------------------------------- */
@@ -76,10 +76,10 @@
// allocates memory on GPU
print_CUDA_error_if_any(cudaMalloc((void**)d_array_addr_ptr,size*sizeof(realw)),
- 20000);
+ 22001);
// copies values onto GPU
print_CUDA_error_if_any(cudaMemcpy((realw*) *d_array_addr_ptr,h_array,size*sizeof(realw),cudaMemcpyHostToDevice),
- 20001);
+ 22002);
}
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/model_s362ani.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/model_s362ani.f90 2012-09-04 18:44:13 UTC (rev 20681)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/model_s362ani.f90 2012-09-04 20:41:19 UTC (rev 20682)
@@ -1757,6 +1757,8 @@
subroutine ylm(XLAT,XLON,LMAX,Y,WK1,WK2,WK3)
+ use model_s362ani_par,only : maxl
+
implicit none
complex TEMP,FAC,DFAC
@@ -1770,7 +1772,7 @@
!
real(kind=4) WK1(LMAX+1),WK2(LMAX+1),WK3(LMAX+1)
real(kind=4) XLAT,XLON
- real(kind=4) Y(1) !! Y should go at least from 1 to fac(LMAX)
+ real(kind=4),dimension((maxl+1)**2) :: Y !! Y should go at least from 1 to fac(LMAX)
real(kind=4), parameter :: RADIAN = 57.2957795
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/save_arrays_solver.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/save_arrays_solver.f90 2012-09-04 18:44:13 UTC (rev 20681)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/save_arrays_solver.f90 2012-09-04 20:41:19 UTC (rev 20682)
@@ -74,21 +74,90 @@
! local parameters
integer :: i,j,k,ispec,iglob,ier
+ real(kind=CUSTOM_REAL),dimension(:),allocatable :: tmp_array
+ ! mesh databases
+ open(unit=27,file=prname(1:len_trim(prname))//'solver_data.bin', &
+ status='unknown',form='unformatted',action='write',iostat=ier)
+ if( ier /= 0 ) call exit_mpi(myrank,'error opening solver_data.bin file')
+
! save nspec and nglob, to be used in combine_paraview_data
- open(unit=27,file=prname(1:len_trim(prname))//'array_dims.txt', &
- status='unknown',action='write',iostat=ier)
- if( ier /= 0 ) call exit_mpi(myrank,'error opening array_dims file')
+ write(27) nspec
+ write(27) nglob
- write(27,*) nspec
- write(27,*) nglob
- close(27)
+ ! mesh topology
- ! mesh parameters
- open(unit=27,file=prname(1:len_trim(prname))//'solver_data_1.bin', &
- status='unknown',form='unformatted',action='write',iostat=ier)
- if( ier /= 0 ) call exit_mpi(myrank,'error opening solver_data_1.bin file')
+ ! mesh arrays used in the solver to locate source and receivers
+ ! and for anisotropy and gravity, save in single precision
+ ! use tmp_array for temporary storage to perform conversion
+ allocate(tmp_array(nglob),stat=ier)
+ if( ier /=0 ) call exit_MPI(myrank,'error allocating temporary array for mesh topology')
+ !--- x coordinate
+ tmp_array(:) = 0._CUSTOM_REAL
+ do ispec = 1,nspec
+ do k = 1,NGLLZ
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+ ! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ tmp_array(iglob) = sngl(xstore(i,j,k,ispec))
+ else
+ tmp_array(iglob) = xstore(i,j,k,ispec)
+ endif
+ enddo
+ enddo
+ enddo
+ enddo
+ write(27) tmp_array
+
+ !--- y coordinate
+ tmp_array(:) = 0._CUSTOM_REAL
+ do ispec = 1,nspec
+ do k = 1,NGLLZ
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+ ! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ tmp_array(iglob) = sngl(ystore(i,j,k,ispec))
+ else
+ tmp_array(iglob) = ystore(i,j,k,ispec)
+ endif
+ enddo
+ enddo
+ enddo
+ enddo
+ write(27) tmp_array
+
+ !--- z coordinate
+ tmp_array(:) = 0._CUSTOM_REAL
+ do ispec = 1,nspec
+ do k = 1,NGLLZ
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+ ! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ tmp_array(iglob) = sngl(zstore(i,j,k,ispec))
+ else
+ tmp_array(iglob) = zstore(i,j,k,ispec)
+ endif
+ enddo
+ enddo
+ enddo
+ enddo
+ write(27) tmp_array
+
+ deallocate(tmp_array)
+
+ ! local to global indexing
+ write(27) ibool
+ write(27) idoubling
+ write(27) ispec_is_tiso
+
+ ! local GLL points
write(27) xixstore
write(27) xiystore
write(27) xizstore
@@ -184,82 +253,9 @@
! additional ocean load mass matrix if oceans and if we are in the crust
if(OCEANS .and. iregion_code == IREGION_CRUST_MANTLE) write(27) rmass_ocean_load
- close(27)
+ close(27) ! solver_data.bin
- ! mesh topology
- open(unit=27,file=prname(1:len_trim(prname))//'solver_data_2.bin', &
- status='unknown',form='unformatted',action='write',iostat=ier)
- if( ier /= 0 ) call exit_mpi(myrank,'error opening solver_data_2.bin file')
-! mesh arrays used in the solver to locate source and receivers
-! and for anisotropy and gravity, save in single precision
-! use rmassz for temporary storage to perform conversion, since already saved
-
- !--- x coordinate
- rmassz(:) = 0._CUSTOM_REAL
- do ispec = 1,nspec
- do k = 1,NGLLZ
- do j = 1,NGLLY
- do i = 1,NGLLX
- iglob = ibool(i,j,k,ispec)
- ! distinguish between single and double precision for reals
- if(CUSTOM_REAL == SIZE_REAL) then
- rmassz(iglob) = sngl(xstore(i,j,k,ispec))
- else
- rmassz(iglob) = xstore(i,j,k,ispec)
- endif
- enddo
- enddo
- enddo
- enddo
- write(27) rmassz
-
- !--- y coordinate
- rmassz(:) = 0._CUSTOM_REAL
- do ispec = 1,nspec
- do k = 1,NGLLZ
- do j = 1,NGLLY
- do i = 1,NGLLX
- iglob = ibool(i,j,k,ispec)
- ! distinguish between single and double precision for reals
- if(CUSTOM_REAL == SIZE_REAL) then
- rmassz(iglob) = sngl(ystore(i,j,k,ispec))
- else
- rmassz(iglob) = ystore(i,j,k,ispec)
- endif
- enddo
- enddo
- enddo
- enddo
- write(27) rmassz
-
- !--- z coordinate
- rmassz(:) = 0._CUSTOM_REAL
- do ispec = 1,nspec
- do k = 1,NGLLZ
- do j = 1,NGLLY
- do i = 1,NGLLX
- iglob = ibool(i,j,k,ispec)
- ! distinguish between single and double precision for reals
- if(CUSTOM_REAL == SIZE_REAL) then
- rmassz(iglob) = sngl(zstore(i,j,k,ispec))
- else
- rmassz(iglob) = zstore(i,j,k,ispec)
- endif
- enddo
- enddo
- enddo
- enddo
- write(27) rmassz
-
- write(27) ibool
-
- write(27) idoubling
-
- write(27) ispec_is_tiso
-
- close(27)
-
! absorbing boundary parameters
open(unit=27,file=prname(1:len_trim(prname))//'boundary.bin', &
status='unknown',form='unformatted',action='write',iostat=ier)
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/setup_MPI_interfaces.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/setup_MPI_interfaces.f90 2012-09-04 18:44:13 UTC (rev 20681)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/setup_MPI_interfaces.f90 2012-09-04 20:41:19 UTC (rev 20682)
@@ -175,7 +175,8 @@
buffer_send_chunkcorn_scalar,buffer_recv_chunkcorn_scalar, &
NUMMSGS_FACES,NUM_MSG_TYPES,NCORNERSCHUNKS, &
NPROC_XI,NPROC_ETA,NGLOB1D_RADIAL(IREGION_CRUST_MANTLE), &
- NGLOB2DMAX_XMIN_XMAX(IREGION_CRUST_MANTLE),NGLOB2DMAX_YMIN_YMAX(IREGION_CRUST_MANTLE),NGLOB2DMAX_XY,NCHUNKS)
+ NGLOB2DMAX_XMIN_XMAX(IREGION_CRUST_MANTLE),NGLOB2DMAX_YMIN_YMAX(IREGION_CRUST_MANTLE), &
+ NGLOB2DMAX_XY,NCHUNKS)
! removes own myrank id (+1)
test_flag(:) = test_flag(:) - ( myrank + 1.0)
@@ -304,7 +305,8 @@
buffer_send_chunkcorn_scalar,buffer_recv_chunkcorn_scalar, &
NUMMSGS_FACES,NUM_MSG_TYPES,NCORNERSCHUNKS, &
NPROC_XI,NPROC_ETA,NGLOB1D_RADIAL(IREGION_OUTER_CORE), &
- NGLOB2DMAX_XMIN_XMAX(IREGION_OUTER_CORE),NGLOB2DMAX_YMIN_YMAX(IREGION_OUTER_CORE),NGLOB2DMAX_XY,NCHUNKS)
+ NGLOB2DMAX_XMIN_XMAX(IREGION_OUTER_CORE),NGLOB2DMAX_YMIN_YMAX(IREGION_OUTER_CORE), &
+ NGLOB2DMAX_XY,NCHUNKS)
! removes own myrank id (+1)
@@ -448,7 +450,8 @@
buffer_send_chunkcorn_scalar,buffer_recv_chunkcorn_scalar, &
NUMMSGS_FACES,NUM_MSG_TYPES,NCORNERSCHUNKS, &
NPROC_XI,NPROC_ETA,NGLOB1D_RADIAL(IREGION_INNER_CORE), &
- NGLOB2DMAX_XMIN_XMAX(IREGION_INNER_CORE),NGLOB2DMAX_YMIN_YMAX(IREGION_INNER_CORE),NGLOB2DMAX_XY,NCHUNKS)
+ NGLOB2DMAX_XMIN_XMAX(IREGION_INNER_CORE),NGLOB2DMAX_YMIN_YMAX(IREGION_INNER_CORE), &
+ NGLOB2DMAX_XY,NCHUNKS)
! debug: idoubling inner core
if( DEBUG ) then
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/read_compute_parameters.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/read_compute_parameters.f90 2012-09-04 18:44:13 UTC (rev 20681)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/read_compute_parameters.f90 2012-09-04 20:41:19 UTC (rev 20682)
@@ -226,20 +226,9 @@
ONE_CRUST,HONOR_1D_SPHERICAL_MOHO,CASE_3D,CRUSTAL, &
ANISOTROPIC_INNER_CORE)
- ! compute total number of time steps, rounded to next multiple of 100
+ ! initial guess : compute total number of time steps, rounded to next multiple of 100
NSTEP = 100 * (int(RECORD_LENGTH_IN_MINUTES * 60.d0 / (100.d0*DT)) + 1)
-! if doing benchmark runs to measure scaling of the code for a limited number of time steps only
- if (DO_BENCHMARK_RUN_ONLY) NSTEP = NSTEP_FOR_BENCHMARK
-
- ! noise tomography:
- ! time steps needs to be doubled, due to +/- branches
- if ( NOISE_TOMOGRAPHY /= 0 ) NSTEP = 2*NSTEP-1
-
- ! subsets used to save seismograms must not be larger than the whole time series,
- ! otherwise we waste memory
- if(NTSTEP_BETWEEN_OUTPUT_SEISMOS > NSTEP) NTSTEP_BETWEEN_OUTPUT_SEISMOS = NSTEP
-
! computes a default hdur_movie that creates nice looking movies.
! Sets HDUR_MOVIE as the minimum period the mesh can resolve
if(HDUR_MOVIE <= TINYVAL) &
@@ -251,7 +240,7 @@
call rcp_check_parameters(NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA, &
NCHUNKS,NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES, &
- ATTENUATION_3D,ATTENUATION,ATTENUATION_NEW,ABSORBING_CONDITIONS, &
+ ATTENUATION,ATTENUATION_NEW,ABSORBING_CONDITIONS, &
INCLUDE_CENTRAL_CUBE,OUTPUT_SEISMOS_SAC_ALPHANUM)
! check that mesh can be coarsened in depth three or four times
@@ -362,7 +351,7 @@
subroutine rcp_check_parameters(NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA, &
NCHUNKS,NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES, &
- ATTENUATION_3D,ATTENUATION,ATTENUATION_NEW,ABSORBING_CONDITIONS, &
+ ATTENUATION,ATTENUATION_NEW,ABSORBING_CONDITIONS, &
INCLUDE_CENTRAL_CUBE,OUTPUT_SEISMOS_SAC_ALPHANUM)
implicit none
@@ -373,7 +362,7 @@
double precision ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES
- logical ATTENUATION_3D,ATTENUATION,ATTENUATION_NEW,ABSORBING_CONDITIONS,&
+ logical ATTENUATION,ATTENUATION_NEW,ABSORBING_CONDITIONS,&
INCLUDE_CENTRAL_CUBE,OUTPUT_SEISMOS_SAC_ALPHANUM
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_outer_core.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_outer_core.f90 2012-09-04 18:44:13 UTC (rev 20681)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_outer_core.f90 2012-09-04 20:41:19 UTC (rev 20682)
@@ -160,8 +160,10 @@
sin_two_omega_t = sin(two_omega_earth*time)
! time step deltat of Euler scheme is included in the source
- source_euler_A(i,j,k) = two_omega_deltat * (cos_two_omega_t * dpotentialdyl + sin_two_omega_t * dpotentialdxl)
- source_euler_B(i,j,k) = two_omega_deltat * (sin_two_omega_t * dpotentialdyl - cos_two_omega_t * dpotentialdxl)
+ source_euler_A(i,j,k) = two_omega_deltat &
+ * (cos_two_omega_t * dpotentialdyl + sin_two_omega_t * dpotentialdxl)
+ source_euler_B(i,j,k) = two_omega_deltat &
+ * (sin_two_omega_t * dpotentialdyl - cos_two_omega_t * dpotentialdxl)
A_rotation = A_array_rotation(i,j,k,ispec)
B_rotation = B_array_rotation(i,j,k,ispec)
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_outer_core_Dev.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_outer_core_Dev.f90 2012-09-04 18:44:13 UTC (rev 20681)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_forces_outer_core_Dev.f90 2012-09-04 20:41:19 UTC (rev 20682)
@@ -254,6 +254,9 @@
source_euler_B(i,j,k) = two_omega_deltat &
* (sin_two_omega_t * dpotentialdyl - cos_two_omega_t * dpotentialdxl)
+ A_rotation = A_array_rotation(i,j,k,ispec)
+ B_rotation = B_array_rotation(i,j,k,ispec)
+
ux_rotation = A_rotation*cos_two_omega_t + B_rotation*sin_two_omega_t
uy_rotation = - A_rotation*sin_two_omega_t + B_rotation*cos_two_omega_t
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_receivers.f90 2012-09-04 18:44:13 UTC (rev 20681)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_receivers.f90 2012-09-04 20:41:19 UTC (rev 20682)
@@ -521,13 +521,17 @@
! gamma does not change since we know the receiver is exactly on the surface
dxi = xix*dx + xiy*dy + xiz*dz
deta = etax*dx + etay*dy + etaz*dz
- if(RECEIVERS_CAN_BE_BURIED) dgamma = gammax*dx + gammay*dy + gammaz*dz
! update values
xi = xi + dxi
eta = eta + deta
- if(RECEIVERS_CAN_BE_BURIED) gamma = gamma + dgamma
+ ! buried receivers vary in z depth
+ if(RECEIVERS_CAN_BE_BURIED) then
+ dgamma = gammax*dx + gammay*dy + gammaz*dz
+ gamma = gamma + dgamma
+ endif
+
! impose that we stay in that element
! (useful if user gives a receiver outside the mesh for instance)
! we can go slightly outside the [1,1] segment since with finite elements
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_sources.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_sources.f90 2012-09-04 18:44:13 UTC (rev 20681)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_sources.f90 2012-09-04 20:41:19 UTC (rev 20682)
@@ -37,9 +37,9 @@
use specfem_par,only: &
NSOURCES,myrank, &
tshift_cmt,theta_source,phi_source, &
- NSTEP,DT,hdur,Mxx,Myy,Mzz,Mxy,Mxz,Myz, &
+ DT,hdur,Mxx,Myy,Mzz,Mxy,Mxz,Myz, &
rspl,espl,espl2,nspl,ibathy_topo, &
- PRINT_SOURCE_TIME_FUNCTION,LOCAL_TMP_PATH,SIMULATION_TYPE,TOPOGRAPHY, &
+ LOCAL_TMP_PATH,SIMULATION_TYPE,TOPOGRAPHY, &
xigll,yigll,zigll, &
xi_source,eta_source,gamma_source,nu_source, &
islice_selected_source,ispec_selected_source
@@ -66,7 +66,6 @@
integer i,j,k,ispec,iglob
integer ier
-
double precision ell
double precision elevation
double precision r0,dcost,p20
@@ -95,9 +94,6 @@
double precision x_target_source,y_target_source,z_target_source
double precision r_target_source
- ! timer MPI
- double precision time_start,tCPU
-
integer isources_already_done,isource_in_this_subset
integer, dimension(:), allocatable :: ispec_selected_source_subset
@@ -111,8 +107,6 @@
double precision moment_tensor(6,NSOURCES)
double precision radius
- character(len=150) OUTPUT_FILES
-
double precision, dimension(:), allocatable :: x_found_source,y_found_source,z_found_source
double precision r_found_source
double precision st,ct,sp,cp
@@ -141,9 +135,12 @@
integer,parameter :: MIDY = (NGLLY+1)/2
integer,parameter :: MIDZ = (NGLLZ+1)/2
- ! timing
+ ! timer MPI
+ double precision time_start,tCPU
double precision, external :: wtime
+ character(len=150) :: OUTPUT_FILES
+
! get MPI starting time for all sources
time_start = wtime()
@@ -151,9 +148,6 @@
ispec_selected_source(:) = 0
final_distance_source(:) = HUGEVAL
- ! get the base pathname for output files
- call get_value_string(OUTPUT_FILES, 'OUTPUT_FILES', 'OUTPUT_FILES')
-
! read all the sources
if(myrank == 0) call get_cmt(yr,jda,ho,mi,sec,tshift_cmt,hdur,lat,long,depth,moment_tensor, &
DT,NSOURCES,min_tshift_cmt_original)
@@ -186,6 +180,9 @@
! appends receiver locations to sr.vtk file
if( myrank == 0 ) then
+ ! get the base pathname for output files
+ call get_value_string(OUTPUT_FILES, 'OUTPUT_FILES', 'OUTPUT_FILES')
+
open(IOVTK,file=trim(OUTPUT_FILES)//'/sr_tmp.vtk', &
position='append',status='old',iostat=ier)
if( ier /= 0 ) call exit_MPI(myrank,'Error opening and appending sources to file sr_tmp.vtk')
@@ -707,10 +704,6 @@
write(IMAIN,*) '*****************************************************'
endif
- ! print source time function and spectrum
- if(PRINT_SOURCE_TIME_FUNCTION) call print_stf(NSOURCES,isource,moment_tensor,tshift_cmt,hdur, &
- min_tshift_cmt_original,NSTEP,DT,OUTPUT_FILES)
-
enddo ! end of loop on all the sources within current source subset
endif ! end of section executed by main process only
@@ -858,8 +851,8 @@
!
- subroutine print_stf(NSOURCES,isource,moment_tensor,tshift_cmt,hdur, &
- min_tshift_cmt_original,NSTEP,DT,OUTPUT_FILES)
+ subroutine print_stf(NSOURCES,isource,Mxx,Myy,Mzz,Mxy,Mxz,Myz, &
+ tshift_cmt,hdur,min_tshift_cmt_original,NSTEP,DT)
! prints source time function
@@ -869,17 +862,16 @@
integer :: NSOURCES,isource
- double precision,dimension(6,NSOURCES) :: moment_tensor
+ double precision,dimension(NSOURCES) :: Mxx,Myy,Mzz,Mxy,Mxz,Myz
double precision,dimension(NSOURCES) :: tshift_cmt,hdur
double precision :: min_tshift_cmt_original
integer :: NSTEP
double precision :: DT
- character(len=150) OUTPUT_FILES
! local parameters
- integer :: i,it,iom,ier
+ integer :: it,iom,ier
double precision :: scalar_moment
double precision :: t0, hdur_gaussian(NSOURCES)
double precision :: t_cmt_used(NSOURCES)
@@ -889,7 +881,8 @@
double precision, external :: comp_source_time_function,comp_source_spectrum
double precision, external :: comp_source_time_function_rickr
- character(len=150) plot_file
+ character(len=150) :: OUTPUT_FILES
+ character(len=150) :: plot_file
! number of points to plot the source time function and spectrum
integer, parameter :: NSAMP_PLOT_SOURCE = 1000
@@ -898,6 +891,9 @@
write(IMAIN,*)
write(IMAIN,*) 'printing the source-time function'
+ ! get the base pathname for output files
+ call get_value_string(OUTPUT_FILES, 'OUTPUT_FILES', 'OUTPUT_FILES')
+
! print the source-time function
if(NSOURCES == 1) then
plot_file = '/plot_source_time_function.txt'
@@ -911,14 +907,13 @@
endif
endif
+ ! output file
open(unit=27,file=trim(OUTPUT_FILES)//plot_file, &
status='unknown',iostat=ier)
if( ier /= 0 ) call exit_mpi(0,'error opening plot_source_time_function file')
- scalar_moment = 0.
- do i = 1,6
- scalar_moment = scalar_moment + moment_tensor(i,isource)**2
- enddo
+ scalar_moment = Mxx(isource)**2 + Myy(isource)**2 + Mzz(isource)**2 &
+ + Mxy(isource)**2 + Mxz(isource)**2 + Myz(isource)**2
scalar_moment = dsqrt(scalar_moment/2.0d0)
! define t0 as the earliest start time
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/prepare_timerun.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/prepare_timerun.f90 2012-09-04 18:44:13 UTC (rev 20681)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/prepare_timerun.f90 2012-09-04 20:41:19 UTC (rev 20682)
@@ -75,6 +75,9 @@
! this makes indexing and timing easier to match with adjoint wavefields indexing.
call read_forward_arrays_startrun()
+ ! prepares stacey boundary arrays for re-construction of wavefields
+ if(ABSORBING_CONDITIONS) call prepare_timerun_stacey()
+
! prepares noise simulations
call prepare_timerun_noise()
@@ -1074,6 +1077,290 @@
!-------------------------------------------------------------------------------------------------
!
+ subroutine prepare_timerun_stacey()
+
+! sets up arrays for Stacey conditions
+
+ use specfem_par
+ use specfem_par_crustmantle
+ use specfem_par_outercore
+
+ implicit none
+ ! local parameters
+ integer :: ier
+ integer :: nabs_xmin_cm,nabs_xmax_cm,nabs_ymin_cm,nabs_ymax_cm
+ integer :: nabs_xmin_oc,nabs_xmax_oc,nabs_ymin_oc,nabs_ymax_oc,nabs_zmin_oc
+ integer(kind=8) :: filesize
+
+ ! sets up absorbing boundary buffer arrays
+
+ ! crust_mantle
+
+ ! create name of database
+ call create_name_database(prname,myrank,IREGION_CRUST_MANTLE,LOCAL_PATH)
+
+ ! allocates buffers
+ if (nspec2D_xmin_crust_mantle > 0 .and. (SIMULATION_TYPE == 3 &
+ .or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))) then
+ nabs_xmin_cm = nspec2D_xmin_crust_mantle
+ else
+ nabs_xmin_cm = 1
+ endif
+ allocate(absorb_xmin_crust_mantle(NDIM,NGLLY,NGLLZ,nabs_xmin_cm),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating absorb xmin')
+
+ if (nspec2D_xmax_crust_mantle > 0 .and. (SIMULATION_TYPE == 3 &
+ .or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))) then
+ nabs_xmax_cm = nspec2D_xmax_crust_mantle
+ else
+ nabs_xmax_cm = 1
+ endif
+ allocate(absorb_xmax_crust_mantle(NDIM,NGLLY,NGLLZ,nabs_xmax_cm),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating absorb xmax')
+
+ if (nspec2D_ymin_crust_mantle > 0 .and. (SIMULATION_TYPE == 3 &
+ .or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))) then
+ nabs_ymin_cm = nspec2D_ymin_crust_mantle
+ else
+ nabs_ymin_cm = 1
+ endif
+ allocate(absorb_ymin_crust_mantle(NDIM,NGLLX,NGLLZ,nabs_ymin_cm),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating absorb ymin')
+
+ if (nspec2D_ymax_crust_mantle > 0 .and. (SIMULATION_TYPE == 3 &
+ .or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))) then
+ nabs_ymax_cm = nspec2D_ymax_crust_mantle
+ else
+ nabs_ymax_cm = 1
+ endif
+ allocate(absorb_ymax_crust_mantle(NDIM,NGLLX,NGLLZ,nabs_ymax_cm),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating absorb ymax')
+
+ ! file I/O for re-construction of wavefields
+ 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
+ call open_file_abs_r(0,trim(prname)//'absorb_xmin.bin',len_trim(trim(prname)//'absorb_xmin.bin'), &
+ filesize)
+ else
+ call open_file_abs_w(0,trim(prname)//'absorb_xmin.bin',len_trim(trim(prname)//'absorb_xmin.bin'), &
+ filesize)
+ endif
+ endif
+ 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
+ call open_file_abs_r(1,trim(prname)//'absorb_xmax.bin',len_trim(trim(prname)//'absorb_xmax.bin'), &
+ filesize)
+ else
+ call open_file_abs_w(1,trim(prname)//'absorb_xmax.bin',len_trim(trim(prname)//'absorb_xmax.bin'), &
+ filesize)
+ endif
+ endif
+ 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)
+
+ ! total file size
+ filesize = reclen_ymin_crust_mantle
+ filesize = filesize*NSTEP
+
+
+ if (SIMULATION_TYPE == 3) then
+ call open_file_abs_r(2,trim(prname)//'absorb_ymin.bin',len_trim(trim(prname)//'absorb_ymin.bin'), &
+ filesize)
+ else
+ call open_file_abs_w(2,trim(prname)//'absorb_ymin.bin',len_trim(trim(prname)//'absorb_ymin.bin'), &
+ filesize)
+ endif
+ endif
+ 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
+ call open_file_abs_r(3,trim(prname)//'absorb_ymax.bin',len_trim(trim(prname)//'absorb_ymax.bin'), &
+ filesize)
+ else
+ call open_file_abs_w(3,trim(prname)//'absorb_ymax.bin',len_trim(trim(prname)//'absorb_ymax.bin'), &
+ filesize)
+ endif
+ endif
+
+
+ ! outer_core
+
+ ! create name of database
+ call create_name_database(prname,myrank,IREGION_OUTER_CORE,LOCAL_PATH)
+
+ ! allocates buffers
+ if (nspec2D_xmin_outer_core > 0 .and. (SIMULATION_TYPE == 3 &
+ .or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))) then
+ nabs_xmin_oc = nspec2D_xmin_outer_core
+ else
+ nabs_xmin_oc = 1
+ endif
+ allocate(absorb_xmin_outer_core(NGLLY,NGLLZ,nabs_xmin_oc),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating absorb xmin')
+
+ if (nspec2D_xmax_outer_core > 0 .and. (SIMULATION_TYPE == 3 &
+ .or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))) then
+ nabs_xmax_oc = nspec2D_xmax_outer_core
+ else
+ nabs_xmax_oc = 1
+ endif
+ allocate(absorb_xmax_outer_core(NGLLY,NGLLZ,nabs_xmax_oc),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating absorb xmax')
+
+ if (nspec2D_ymin_outer_core > 0 .and. (SIMULATION_TYPE == 3 &
+ .or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))) then
+ nabs_ymin_oc = nspec2D_ymin_outer_core
+ else
+ nabs_ymin_oc = 1
+ endif
+ allocate(absorb_ymin_outer_core(NGLLX,NGLLZ,nabs_ymin_oc),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating absorb ymin')
+
+ if (nspec2D_ymax_outer_core > 0 .and. (SIMULATION_TYPE == 3 &
+ .or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))) then
+ nabs_ymax_oc = nspec2D_ymax_outer_core
+ else
+ nabs_ymax_oc = 1
+ endif
+ allocate(absorb_ymax_outer_core(NGLLX,NGLLZ,nabs_ymax_oc),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating absorb ymax')
+
+ if (nspec2D_zmin_outer_core > 0 .and. (SIMULATION_TYPE == 3 &
+ .or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))) then
+ nabs_zmin_oc = nspec2D_zmin_outer_core
+ else
+ nabs_zmin_oc = 1
+ endif
+ allocate(absorb_zmin_outer_core(NGLLX,NGLLY,nabs_zmin_oc),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating absorb zmin')
+
+ ! file I/O for re-construction of wavefields
+ 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
+ call open_file_abs_r(4,trim(prname)//'absorb_xmin.bin',len_trim(trim(prname)//'absorb_ymax.bin'), &
+ filesize)
+ else
+ call open_file_abs_w(4,trim(prname)//'absorb_xmin.bin',len_trim(trim(prname)//'absorb_ymax.bin'), &
+ filesize)
+ endif
+ endif
+ 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
+ call open_file_abs_r(5,trim(prname)//'absorb_xmax.bin',len_trim(trim(prname)//'absorb_xmax.bin'), &
+ filesize)
+ else
+ call open_file_abs_w(5,trim(prname)//'absorb_xmax.bin',len_trim(trim(prname)//'absorb_xmax.bin'), &
+ filesize)
+ endif
+ endif
+ 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
+ call open_file_abs_r(6,trim(prname)//'absorb_ymin.bin',len_trim(trim(prname)//'absorb_ymin.bin'), &
+ filesize)
+ else
+ call open_file_abs_w(6,trim(prname)//'absorb_ymin.bin',len_trim(trim(prname)//'absorb_ymin.bin'), &
+ filesize)
+ endif
+ endif
+ 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
+ call open_file_abs_r(7,trim(prname)//'absorb_ymax.bin',len_trim(trim(prname)//'absorb_ymax.bin'), &
+ filesize)
+ else
+ call open_file_abs_w(7,trim(prname)//'absorb_ymax.bin',len_trim(trim(prname)//'absorb_ymax.bin'), &
+ filesize)
+ endif
+ endif
+ if (nspec2D_zmin_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_zmin_outer_core)
+
+ ! total file size
+ filesize = reclen_zmin
+ filesize = filesize*NSTEP
+
+ if (SIMULATION_TYPE == 3) then
+ call open_file_abs_r(8,trim(prname)//'absorb_zmin.bin',len_trim(trim(prname)//'absorb_zmin.bin'), &
+ filesize)
+ else
+ call open_file_abs_w(8,trim(prname)//'absorb_zmin.bin',len_trim(trim(prname)//'absorb_zmin.bin'), &
+ filesize)
+ endif
+ endif
+
+ end subroutine prepare_timerun_stacey
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
subroutine prepare_timerun_noise()
use specfem_par
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_arrays_solver.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_arrays_solver.f90 2012-09-04 18:44:13 UTC (rev 20681)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_arrays_solver.f90 2012-09-04 20:41:19 UTC (rev 20682)
@@ -88,17 +88,46 @@
character(len=150) :: LOCAL_PATH
! local parameters
- integer :: ier
+ integer :: ier,lnspec,lnglob
! processor identification
character(len=150) :: prname
! create the name for the database of the current slide and region
call create_name_database(prname,myrank,iregion_code,LOCAL_PATH)
- open(unit=IIN,file=prname(1:len_trim(prname))//'solver_data_1.bin', &
+ open(unit=IIN,file=prname(1:len_trim(prname))//'solver_data.bin', &
status='old',action='read',form='unformatted',iostat=ier)
- if( ier /= 0 ) call exit_mpi(myrank,'error opening solver_data_1.bin')
+ if( ier /= 0 ) call exit_mpi(myrank,'error opening solver_data.bin')
+ ! read coordinates of the mesh
+
+ read(IIN) lnspec
+ read(IIN) lnglob
+
+ ! checks dimensions
+ if( lnspec /= nspec ) then
+ close(IIN)
+ print*,'error file dimension: nspec in file = ',lnspec,' but nspec desired:',nspec
+ print*,'please check file ',prname(1:len_trim(prname))//'solver_data.bin'
+ call exit_mpi(myrank,'error dimensions in solver_data.bin')
+ endif
+ if( lnglob /= nglob ) then
+ close(IIN)
+ print*,'error file dimension: nglob in file = ',lnglob,' but nglob desired:',nglob
+ print*,'please check file ',prname(1:len_trim(prname))//'solver_data.bin'
+ call exit_mpi(myrank,'error dimensions in solver_data.bin')
+ endif
+
+ ! mesh coordinates
+ read(IIN) xstore
+ read(IIN) ystore
+ read(IIN) zstore
+
+ read(IIN) ibool
+ read(IIN) idoubling
+ read(IIN) ispec_is_tiso
+
+ ! local GLL points
read(IIN) xix
read(IIN) xiy
read(IIN) xiz
@@ -184,24 +213,8 @@
! read additional ocean load mass matrix
if(OCEANS_VAL .and. iregion_code == IREGION_CRUST_MANTLE) read(IIN) rmass_ocean_load
- close(IIN)
+ close(IIN) ! solver_data.bin
- ! read coordinates of the mesh
- open(unit=IIN,file=prname(1:len_trim(prname))//'solver_data_2.bin', &
- status='old',action='read',form='unformatted',iostat=ier)
- if( ier /= 0 ) call exit_mpi(myrank,'error opening file solver_data_2.bin')
- read(IIN) xstore
- read(IIN) ystore
- read(IIN) zstore
-
- read(IIN) ibool
-
- read(IIN) idoubling
-
- read(IIN) ispec_is_tiso
-
- close(IIN)
-
end subroutine read_arrays_solver
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_mesh_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_mesh_databases.f90 2012-09-04 18:44:13 UTC (rev 20681)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_mesh_databases.f90 2012-09-04 20:41:19 UTC (rev 20682)
@@ -44,13 +44,13 @@
! start reading the databases
! read arrays created by the mesher
- ! reads "solver_data_1.bin" & "solver_data_2.bin" files for crust and mantle
+ ! reads "solver_data.bin" files for crust and mantle
call read_mesh_databases_CM()
- ! reads "solver_data_1.bin" & "solver_data_2.bin" files for outer core
+ ! reads "solver_data.bin" files for outer core
call read_mesh_databases_OC()
- ! reads "solver_data_1.bin" & "solver_data_2.bin" files for inner core
+ ! reads "solver_data.bin" files for inner core
call read_mesh_databases_IC()
! reads "boundary.bin" files to couple mantle with outer core and inner core boundaries
@@ -64,7 +64,7 @@
! absorbing boundaries
if(ABSORBING_CONDITIONS) then
- ! reads "stacey.bin" files and sets up arrays for Stacey conditions
+ ! reads "stacey.bin" files
call read_mesh_databases_stacey()
endif
@@ -986,97 +986,8 @@
implicit none
! local parameters
- integer(kind=8) :: filesize
integer :: ier
- integer :: nabs_xmin_cm,nabs_xmax_cm,nabs_ymin_cm,nabs_ymax_cm
- integer :: nabs_xmin_oc,nabs_xmax_oc,nabs_ymin_oc,nabs_ymax_oc,nabs_zmin_oc
- ! sets up absorbing boundary buffer arrays
-
- ! crust_mantle
- if (nspec2D_xmin_crust_mantle > 0 .and. (SIMULATION_TYPE == 3 &
- .or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))) then
- nabs_xmin_cm = nspec2D_xmin_crust_mantle
- else
- nabs_xmin_cm = 1
- endif
- allocate(absorb_xmin_crust_mantle(NDIM,NGLLY,NGLLZ,nabs_xmin_cm),stat=ier)
- if( ier /= 0 ) call exit_MPI(myrank,'error allocating absorb xmin')
-
- if (nspec2D_xmax_crust_mantle > 0 .and. (SIMULATION_TYPE == 3 &
- .or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))) then
- nabs_xmax_cm = nspec2D_xmax_crust_mantle
- else
- nabs_xmax_cm = 1
- endif
- allocate(absorb_xmax_crust_mantle(NDIM,NGLLY,NGLLZ,nabs_xmax_cm),stat=ier)
- if( ier /= 0 ) call exit_MPI(myrank,'error allocating absorb xmax')
-
- if (nspec2D_ymin_crust_mantle > 0 .and. (SIMULATION_TYPE == 3 &
- .or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))) then
- nabs_ymin_cm = nspec2D_ymin_crust_mantle
- else
- nabs_ymin_cm = 1
- endif
- allocate(absorb_ymin_crust_mantle(NDIM,NGLLX,NGLLZ,nabs_ymin_cm),stat=ier)
- if( ier /= 0 ) call exit_MPI(myrank,'error allocating absorb ymin')
-
- if (nspec2D_ymax_crust_mantle > 0 .and. (SIMULATION_TYPE == 3 &
- .or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))) then
- nabs_ymax_cm = nspec2D_ymax_crust_mantle
- else
- nabs_ymax_cm = 1
- endif
- allocate(absorb_ymax_crust_mantle(NDIM,NGLLX,NGLLZ,nabs_ymax_cm),stat=ier)
- if( ier /= 0 ) call exit_MPI(myrank,'error allocating absorb ymax')
-
- ! outer_core
- if (nspec2D_xmin_outer_core > 0 .and. (SIMULATION_TYPE == 3 &
- .or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))) then
- nabs_xmin_oc = nspec2D_xmin_outer_core
- else
- nabs_xmin_oc = 1
- endif
- allocate(absorb_xmin_outer_core(NGLLY,NGLLZ,nabs_xmin_oc),stat=ier)
- if( ier /= 0 ) call exit_MPI(myrank,'error allocating absorb xmin')
-
- if (nspec2D_xmax_outer_core > 0 .and. (SIMULATION_TYPE == 3 &
- .or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))) then
- nabs_xmax_oc = nspec2D_xmax_outer_core
- else
- nabs_xmax_oc = 1
- endif
- allocate(absorb_xmax_outer_core(NGLLY,NGLLZ,nabs_xmax_oc),stat=ier)
- if( ier /= 0 ) call exit_MPI(myrank,'error allocating absorb xmax')
-
- if (nspec2D_ymin_outer_core > 0 .and. (SIMULATION_TYPE == 3 &
- .or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))) then
- nabs_ymin_oc = nspec2D_ymin_outer_core
- else
- nabs_ymin_oc = 1
- endif
- allocate(absorb_ymin_outer_core(NGLLX,NGLLZ,nabs_ymin_oc),stat=ier)
- if( ier /= 0 ) call exit_MPI(myrank,'error allocating absorb ymin')
-
- if (nspec2D_ymax_outer_core > 0 .and. (SIMULATION_TYPE == 3 &
- .or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))) then
- nabs_ymax_oc = nspec2D_ymax_outer_core
- else
- nabs_ymax_oc = 1
- endif
- allocate(absorb_ymax_outer_core(NGLLX,NGLLZ,nabs_ymax_oc),stat=ier)
- if( ier /= 0 ) call exit_MPI(myrank,'error allocating absorb ymax')
-
- if (nspec2D_zmin_outer_core > 0 .and. (SIMULATION_TYPE == 3 &
- .or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD))) then
- nabs_zmin_oc = nspec2D_zmin_outer_core
- else
- nabs_zmin_oc = 1
- endif
- allocate(absorb_zmin_outer_core(NGLLX,NGLLY,nabs_zmin_oc),stat=ier)
- if( ier /= 0 ) call exit_MPI(myrank,'error allocating absorb zmin')
-
-
! crust and mantle
! create name of database
@@ -1084,7 +995,9 @@
! read arrays for Stacey conditions
open(unit=27,file=prname(1:len_trim(prname))//'stacey.bin', &
- status='old',form='unformatted',action='read')
+ status='old',form='unformatted',action='read',iostat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error opening stacey.bin file for crust mantle')
+
read(27) nimin_crust_mantle
read(27) nimax_crust_mantle
read(27) njmin_crust_mantle
@@ -1093,84 +1006,6 @@
read(27) nkmin_eta_crust_mantle
close(27)
- 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
- call open_file_abs_r(0,trim(prname)//'absorb_xmin.bin',len_trim(trim(prname)//'absorb_xmin.bin'), &
- filesize)
- else
- call open_file_abs_w(0,trim(prname)//'absorb_xmin.bin',len_trim(trim(prname)//'absorb_xmin.bin'), &
- filesize)
- endif
- endif
-
- 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
- call open_file_abs_r(1,trim(prname)//'absorb_xmax.bin',len_trim(trim(prname)//'absorb_xmax.bin'), &
- filesize)
- else
- call open_file_abs_w(1,trim(prname)//'absorb_xmax.bin',len_trim(trim(prname)//'absorb_xmax.bin'), &
- filesize)
- endif
- endif
-
- 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)
-
- ! total file size
- filesize = reclen_ymin_crust_mantle
- filesize = filesize*NSTEP
-
-
- if (SIMULATION_TYPE == 3) then
- call open_file_abs_r(2,trim(prname)//'absorb_ymin.bin',len_trim(trim(prname)//'absorb_ymin.bin'), &
- filesize)
- else
- call open_file_abs_w(2,trim(prname)//'absorb_ymin.bin',len_trim(trim(prname)//'absorb_ymin.bin'), &
- filesize)
- endif
- endif
-
- 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
- call open_file_abs_r(3,trim(prname)//'absorb_ymax.bin',len_trim(trim(prname)//'absorb_ymax.bin'), &
- filesize)
- else
- call open_file_abs_w(3,trim(prname)//'absorb_ymax.bin',len_trim(trim(prname)//'absorb_ymax.bin'), &
- filesize)
- endif
- endif
-
-
! outer core
! create name of database
@@ -1178,7 +1013,9 @@
! read arrays for Stacey conditions
open(unit=27,file=prname(1:len_trim(prname))//'stacey.bin', &
- status='old',form='unformatted',action='read')
+ status='old',form='unformatted',action='read',iostat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error opening stacey.bin file for outer core')
+
read(27) nimin_outer_core
read(27) nimax_outer_core
read(27) njmin_outer_core
@@ -1187,103 +1024,5 @@
read(27) nkmin_eta_outer_core
close(27)
- 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
- call open_file_abs_r(4,trim(prname)//'absorb_xmin.bin',len_trim(trim(prname)//'absorb_ymax.bin'), &
- filesize)
- else
- call open_file_abs_w(4,trim(prname)//'absorb_xmin.bin',len_trim(trim(prname)//'absorb_ymax.bin'), &
- filesize)
- endif
- endif
-
- 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
- call open_file_abs_r(5,trim(prname)//'absorb_xmax.bin',len_trim(trim(prname)//'absorb_xmax.bin'), &
- filesize)
- else
- call open_file_abs_w(5,trim(prname)//'absorb_xmax.bin',len_trim(trim(prname)//'absorb_xmax.bin'), &
- filesize)
- endif
-
- endif
-
- 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
- call open_file_abs_r(6,trim(prname)//'absorb_ymin.bin',len_trim(trim(prname)//'absorb_ymin.bin'), &
- filesize)
- else
- call open_file_abs_w(6,trim(prname)//'absorb_ymin.bin',len_trim(trim(prname)//'absorb_ymin.bin'), &
- filesize)
-
- endif
- endif
-
- 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
- call open_file_abs_r(7,trim(prname)//'absorb_ymax.bin',len_trim(trim(prname)//'absorb_ymax.bin'), &
- filesize)
- else
- call open_file_abs_w(7,trim(prname)//'absorb_ymax.bin',len_trim(trim(prname)//'absorb_ymax.bin'), &
- filesize)
-
- endif
- endif
-
- if (nspec2D_zmin_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_zmin_outer_core)
-
- ! total file size
- filesize = reclen_zmin
- filesize = filesize*NSTEP
-
- if (SIMULATION_TYPE == 3) then
- call open_file_abs_r(8,trim(prname)//'absorb_zmin.bin',len_trim(trim(prname)//'absorb_zmin.bin'), &
- filesize)
- else
- call open_file_abs_w(8,trim(prname)//'absorb_zmin.bin',len_trim(trim(prname)//'absorb_zmin.bin'), &
- filesize)
- endif
- endif
-
end subroutine read_mesh_databases_stacey
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/setup_sources_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/setup_sources_receivers.f90 2012-09-04 18:44:13 UTC (rev 20681)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/setup_sources_receivers.f90 2012-09-04 20:41:19 UTC (rev 20682)
@@ -229,6 +229,18 @@
call exit_mpi(myrank,'error negative USER_T0 parameter in constants.h')
endif
+ ! determines number of times steps for simulation
+ call setup_timesteps()
+
+ ! prints source time functions to output files
+ if(PRINT_SOURCE_TIME_FUNCTION .and. myrank == 0) then
+ do isource = 1,NSOURCES
+ ! print source time function and spectrum
+ call print_stf(NSOURCES,isource,Mxx,Myy,Mzz,Mxy,Mxz,Myz, &
+ tshift_cmt,hdur,min_tshift_cmt_original,NSTEP,DT)
+ enddo
+ endif
+
! get information about event name and location
! (e.g. needed for SAC seismograms)
@@ -247,6 +259,39 @@
!-------------------------------------------------------------------------------------------------
!
+ subroutine setup_timesteps()
+
+ use specfem_par
+ implicit none
+
+ ! from intial guess in read_compute_parameters:
+ ! compute total number of time steps, rounded to next multiple of 100
+ ! NSTEP = 100 * (int(RECORD_LENGTH_IN_MINUTES * 60.d0 / (100.d0*DT)) + 1)
+ !
+ ! adds initial t0 time to update number of time steps and reach full record length
+ NSTEP = NSTEP + 100 * (int( abs(t0) / (100.d0*DT)) + 1)
+
+ ! if doing benchmark runs to measure scaling of the code for a limited number of time steps only
+ if (DO_BENCHMARK_RUN_ONLY) NSTEP = NSTEP_FOR_BENCHMARK
+
+ ! noise tomography:
+ ! time steps needs to be doubled, due to +/- branches
+ if ( NOISE_TOMOGRAPHY /= 0 ) NSTEP = 2*NSTEP-1
+
+ ! subsets used to save seismograms must not be larger than the whole time series,
+ ! otherwise we waste memory
+ if(NTSTEP_BETWEEN_OUTPUT_SEISMOS > NSTEP) NTSTEP_BETWEEN_OUTPUT_SEISMOS = NSTEP
+
+ ! re-checks output steps?
+ !if (OUTPUT_SEISMOS_SAC_ALPHANUM .and. (mod(NTSTEP_BETWEEN_OUTPUT_SEISMOS,5)/=0)) &
+ ! stop 'if OUTPUT_SEISMOS_SAC_ALPHANUM = .true. then NTSTEP_BETWEEN_OUTPUT_SEISMOS must be a multiple of 5, check the Par_file'
+
+ end subroutine setup_timesteps
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
subroutine setup_receivers()
use specfem_par
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 2012-09-04 18:44:13 UTC (rev 20681)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_movie_volume.f90 2012-09-04 20:41:19 UTC (rev 20682)
@@ -292,6 +292,8 @@
mask_3dmovie,nu_3dmovie)
+! outputs strains: MOVIE_VOLUME_TYPE == 1 / 2 / 3
+
implicit none
include "constants.h"
@@ -440,6 +442,9 @@
subroutine write_movie_volume_vector(myrank,it,npoints_3dmovie,LOCAL_TMP_PATH,MOVIE_VOLUME_TYPE, &
MOVIE_COARSE,ibool_crust_mantle,vector_crust_mantle, &
scalingval,mask_3dmovie,nu_3dmovie)
+
+! outputs displacement/velocity: MOVIE_VOLUME_TYPE == 5 / 6
+
implicit none
include "constants.h"
@@ -556,6 +561,8 @@
epsilondev_xz_inner_core,epsilondev_yz_inner_core, &
LOCAL_TMP_PATH)
+! outputs divergence and curl: MOVIE_VOLUME_TYPE == 4
+
implicit none
include "constants.h"
include "OUTPUT_FILES/values_from_mesher.h"
@@ -738,6 +745,8 @@
displ_crust_mantle,displ_inner_core,displ_outer_core, &
ibool_crust_mantle,ibool_inner_core,ibool_outer_core)
+! outputs norm of displacement: MOVIE_VOLUME_TYPE == 7
+
implicit none
include "constants.h"
include "OUTPUT_FILES/values_from_mesher.h"
@@ -852,6 +861,8 @@
veloc_crust_mantle,veloc_inner_core,veloc_outer_core, &
ibool_crust_mantle,ibool_inner_core,ibool_outer_core)
+! outputs norm of velocity: MOVIE_VOLUME_TYPE == 8
+
implicit none
include "constants.h"
include "OUTPUT_FILES/values_from_mesher.h"
@@ -966,6 +977,8 @@
accel_crust_mantle,accel_inner_core,accel_outer_core, &
ibool_crust_mantle,ibool_inner_core,ibool_outer_core)
+! outputs norm of acceleration: MOVIE_VOLUME_TYPE == 1
+
implicit none
include "constants.h"
include "OUTPUT_FILES/values_from_mesher.h"
More information about the CIG-COMMITS
mailing list