[cig-commits] r18154 - in seismo/3D/SPECFEM3D/trunk: src/decompose_mesh_SCOTCH src/generate_databases src/meshfem3D src/shared src/specfem3D utils/Visualization/Paraview
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Thu Mar 31 04:08:11 PDT 2011
Author: danielpeter
Date: 2011-03-31 04:08:11 -0700 (Thu, 31 Mar 2011)
New Revision: 18154
Modified:
seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90
seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_mass_matrices.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_MPI.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_coupling_surfaces.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_model.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/model_tomography.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90
seismo/3D/SPECFEM3D/trunk/src/meshfem3D/create_regions_mesh.f90
seismo/3D/SPECFEM3D/trunk/src/meshfem3D/get_global.f90
seismo/3D/SPECFEM3D/trunk/src/meshfem3D/meshfem3D.f90
seismo/3D/SPECFEM3D/trunk/src/meshfem3D/save_databases.f90
seismo/3D/SPECFEM3D/trunk/src/shared/assemble_MPI_scalar.f90
seismo/3D/SPECFEM3D/trunk/src/shared/combine_surf_data.f90
seismo/3D/SPECFEM3D/trunk/src/shared/combine_vol_data.f90
seismo/3D/SPECFEM3D/trunk/src/shared/create_movie_shakemap_AVS_DX_GMT.f90
seismo/3D/SPECFEM3D/trunk/src/shared/create_serial_name_database.f90
seismo/3D/SPECFEM3D/trunk/src/shared/get_attenuation_model.f90
seismo/3D/SPECFEM3D/trunk/src/shared/get_element_face.f90
seismo/3D/SPECFEM3D/trunk/src/shared/read_value_parameters.f90
seismo/3D/SPECFEM3D/trunk/src/shared/smooth_vol_data.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/PML_init.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/assemble_MPI_vector.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_elastic.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_PML.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/create_color_image.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/detect_mesh_surfaces.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/initialize_simulation.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_receivers.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_source.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_topography_bathymetry.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/save_adjoint_kernels.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_GLL_points.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_movie_meshes.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_movie_output.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_output_ASCII.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_seismograms.f90
seismo/3D/SPECFEM3D/trunk/utils/Visualization/Paraview/create_slice_VTK.f90
Log:
adds explicit allocation checks; fixes div/curl computation for movie volume runs; updates gatherv calls for collecting movie data; re-adds missing routine create_regions_mesh_save_moho() to create_regions_mesh.f90
Modified: seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -68,7 +68,7 @@
double precision, dimension(SCOTCH_GRAPHDIM) :: scotchgraph
double precision, dimension(SCOTCH_STRATDIM) :: scotchstrat
character(len=256), parameter :: scotch_strategy='b{job=t,map=t,poli=S,sep=h{pass=30}}'
- integer :: ierr,idummy
+ integer :: ier,idummy
!pll
double precision , dimension(:,:), allocatable :: mat_prop
@@ -96,13 +96,14 @@
! reads node coordinates
open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/nodes_coords_file',&
- status='old', form='formatted', iostat = ierr)
- if( ierr /= 0 ) then
+ status='old', form='formatted', iostat = ier)
+ if( ier /= 0 ) then
print*,'could not open file:',localpath_name(1:len_trim(localpath_name))//'/nodes_coords_file'
stop 'error file open'
endif
read(98,*) nnodes
- allocate(nodes_coords(3,nnodes))
+ allocate(nodes_coords(3,nnodes),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array nodes_coords'
do inode = 1, nnodes
! format: #id_node #x_coordinate #y_coordinate #z_coordinate
read(98,*) num_node, nodes_coords(1,num_node), nodes_coords(2,num_node), nodes_coords(3,num_node)
@@ -116,17 +117,19 @@
!(CUBIT calls this the connectivity, guess in the sense that it connects with the points index in
! the global coordinate file "nodes_coords_file"; it doesn't tell you which point is connected with others)
open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/mesh_file', &
- status='old', form='formatted',iostat=ierr)
- if( ierr /= 0 ) stop 'error opening mesh_file'
+ status='old', form='formatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening mesh_file'
read(98,*) nspec
- allocate(elmnts(esize,nspec))
+ allocate(elmnts(esize,nspec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array elmnts'
do ispec = 1, nspec
! format: # element_id #id_node1 ... #id_node8
! note: be aware that here we can have different node ordering for a cube element;
! the ordering from Cubit files might not be consistent for multiple volumes, or uneven, unstructured grids
!
- ! guess here it assumes that spectral elements ordering is like first at the bottom of the element, anticlock-wise, i.e.
+ ! guess here it assumes that spectral elements ordering is like first
+ ! at the bottom of the element, anticlock-wise, i.e.
! point 1 = (0,0,0), point 2 = (0,1,0), point 3 = (1,1,0), point 4 = (1,0,0)
! then top (positive z-direction) of element
! point 5 = (0,0,1), point 6 = (0,1,1), point 7 = (1,1,1), point 8 = (1,0,1)
@@ -158,9 +161,10 @@
! reads material associations
open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/materials_file', &
- status='old', form='formatted',iostat=ierr)
- if( ierr /= 0 ) stop 'error opening materials_file'
- allocate(mat(2,nspec))
+ status='old', form='formatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening materials_file'
+ allocate(mat(2,nspec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array mat'
mat(:,:) = 0
do ispec = 1, nspec
! format: # id_element #flag
@@ -191,14 +195,14 @@
count_def_mat = 0
count_undef_mat = 0
open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/nummaterial_velocity_file',&
- status='old', form='formatted',iostat=ierr)
- if( ierr /= 0 ) stop 'error opening nummaterial_velocity_file'
+ status='old', form='formatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening nummaterial_velocity_file'
! note: format #material_domain_id #material_id #...
- read(98,*,iostat=ierr) idummy,num_mat
+ read(98,*,iostat=ier) idummy,num_mat
print *,'materials:'
! counts materials (defined/undefined)
- do while (ierr == 0)
+ do while (ier == 0)
print*, ' num_mat = ',num_mat
if(num_mat > 0 ) then
! positive materials_id: velocity values will be defined
@@ -207,7 +211,7 @@
! negative materials_id: undefined material properties yet
count_undef_mat = count_undef_mat + 1
end if
- read(98,*,iostat=ierr) idummy,num_mat
+ read(98,*,iostat=ier) idummy,num_mat
end do
close(98)
print*, ' defined = ',count_def_mat, 'undefined = ',count_undef_mat
@@ -218,15 +222,17 @@
print*,' bigger than defined materials in nummaterial_velocity_file:',count_def_mat
stop 'error materials'
endif
- allocate(mat_prop(6,count_def_mat))
- allocate(undef_mat_prop(6,count_undef_mat))
+ allocate(mat_prop(6,count_def_mat),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array mat_prop'
+ allocate(undef_mat_prop(6,count_undef_mat),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array undef_mat_prop'
mat_prop(:,:) = 0.d0
undef_mat_prop(:,:) = ''
! reads in defined material properties
open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/nummaterial_velocity_file', &
- status='old', form='formatted', iostat=ierr)
- if( ierr /= 0 ) stop 'error opening nummaterial_velocity_file'
+ status='old', form='formatted', iostat=ier)
+ if( ier /= 0 ) stop 'error opening nummaterial_velocity_file'
! note: entries in nummaterial_velocity_file can be an unsorted list of all
! defined materials (material_id > 0) and undefined materials (material_id < 0 )
@@ -239,11 +245,11 @@
!read(98,*) idomain_id,num_mat,rho,vp,vs,qmu,aniso_flag
! reads lines unti it reaches a defined material
num_mat = -1
- do while( num_mat < 0 .and. ierr == 0)
- read(98,'(A256)',iostat=ierr) line
+ do while( num_mat < 0 .and. ier == 0)
+ read(98,'(A256)',iostat=ier) line
read(line,*) idomain_id,num_mat
enddo
- if( ierr /= 0 ) stop 'error reading in defined materials in nummaterial_velocity_file'
+ if( ier /= 0 ) stop 'error reading in defined materials in nummaterial_velocity_file'
! reads in defined material properties
read(line,*) idomain_id,num_mat,rho,vp,vs,qmu,aniso_flag
@@ -263,23 +269,24 @@
end do
! reads in undefined material properties
- rewind(98,iostat=ierr) ! back to the beginning of the file
+ rewind(98,iostat=ier) ! back to the beginning of the file
do imat=1,count_undef_mat
! undefined materials: have to be listed in decreasing order of material_id (start with -1, -2, etc...)
! format:
! - for interfaces
- ! #material_domain_id #material_id(<0) #type_name (="interface") #material_id_for_material_below #material_id_for_material_above
+ ! #material_domain_id #material_id(<0) #type_name (="interface")
+ ! #material_id_for_material_below #material_id_for_material_above
! example: 2 -1 interface 1 2
! - for tomography models
! #material_domain_id #material_id (<0) #type_name (="tomography") #block_name
! example: 2 -1 tomography elastic tomography_model.xyz 1
! reads lines unti it reaches a defined material
num_mat = 1
- do while( num_mat >= 0 .and. ierr == 0 )
- read(98,'(A256)',iostat=ierr) line
+ do while( num_mat >= 0 .and. ier == 0 )
+ read(98,'(A256)',iostat=ier) line
read(line,*) idomain_id,num_mat
enddo
- if( ierr /= 0 ) stop 'error reading in undefined materials in nummaterial_velocity_file'
+ if( ier /= 0 ) stop 'error reading in undefined materials in nummaterial_velocity_file'
! checks if interface or tomography definition
read(line,*) undef_mat_prop(6,imat),undef_mat_prop(1,imat),undef_mat_prop(2,imat)
@@ -370,14 +377,16 @@
! reads in absorbing boundary files
open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/absorbing_surface_file_xmin', &
- status='old', form='formatted',iostat=ierr)
- if( ierr /= 0 ) then
+ status='old', form='formatted',iostat=ier)
+ if( ier /= 0 ) then
nspec2D_xmin = 0
else
read(98,*) nspec2D_xmin
endif
- allocate(ibelm_xmin(nspec2D_xmin))
- allocate(nodes_ibelm_xmin(4,nspec2D_xmin))
+ allocate(ibelm_xmin(nspec2D_xmin),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ibelm_xmin'
+ allocate(nodes_ibelm_xmin(4,nspec2D_xmin),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array nodes_ibelm_xmin'
do ispec2D = 1,nspec2D_xmin
! format: #id_(element containing the face) #id_node1_face .. #id_node4_face
! note: ordering for CUBIT seems such that the normal of the face points outward of the element the face belongs to;
@@ -405,14 +414,16 @@
! reads in absorbing boundary files
open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/absorbing_surface_file_xmax', &
- status='old', form='formatted',iostat=ierr)
- if( ierr /= 0 ) then
+ status='old', form='formatted',iostat=ier)
+ if( ier /= 0 ) then
nspec2D_xmax = 0
else
read(98,*) nspec2D_xmax
endif
- allocate(ibelm_xmax(nspec2D_xmax))
- allocate(nodes_ibelm_xmax(4,nspec2D_xmax))
+ allocate(ibelm_xmax(nspec2D_xmax),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ibelm_xmax'
+ allocate(nodes_ibelm_xmax(4,nspec2D_xmax),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array nodes_ibelm_xmax'
do ispec2D = 1,nspec2D_xmax
! format: #id_(element containing the face) #id_node1_face .. #id_node4_face
read(98,*) ibelm_xmax(ispec2D), nodes_ibelm_xmax(1,ispec2D), nodes_ibelm_xmax(2,ispec2D), &
@@ -423,14 +434,16 @@
! reads in absorbing boundary files
open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/absorbing_surface_file_ymin', &
- status='old', form='formatted',iostat=ierr)
- if( ierr /= 0 ) then
+ status='old', form='formatted',iostat=ier)
+ if( ier /= 0 ) then
nspec2D_ymin = 0
else
read(98,*) nspec2D_ymin
endif
- allocate(ibelm_ymin(nspec2D_ymin))
- allocate(nodes_ibelm_ymin(4,nspec2D_ymin))
+ allocate(ibelm_ymin(nspec2D_ymin),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ibelm_ymin'
+ allocate(nodes_ibelm_ymin(4,nspec2D_ymin),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array nodes_ibelm_ymin'
do ispec2D = 1,nspec2D_ymin
! format: #id_(element containing the face) #id_node1_face .. #id_node4_face
read(98,*) ibelm_ymin(ispec2D), nodes_ibelm_ymin(1,ispec2D), nodes_ibelm_ymin(2,ispec2D), &
@@ -441,14 +454,16 @@
! reads in absorbing boundary files
open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/absorbing_surface_file_ymax', &
- status='old', form='formatted',iostat=ierr)
- if( ierr /= 0 ) then
+ status='old', form='formatted',iostat=ier)
+ if( ier /= 0 ) then
nspec2D_ymax = 0
else
read(98,*) nspec2D_ymax
endif
- allocate(ibelm_ymax(nspec2D_ymax))
- allocate(nodes_ibelm_ymax(4,nspec2D_ymax))
+ allocate(ibelm_ymax(nspec2D_ymax),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ibelm_ymax'
+ allocate(nodes_ibelm_ymax(4,nspec2D_ymax),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array nodes_ibelm_ymax'
do ispec2D = 1,nspec2D_ymax
! format: #id_(element containing the face) #id_node1_face .. #id_node4_face
read(98,*) ibelm_ymax(ispec2D), nodes_ibelm_ymax(1,ispec2D), nodes_ibelm_ymax(2,ispec2D), &
@@ -459,14 +474,16 @@
! reads in absorbing boundary files
open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/absorbing_surface_file_bottom', &
- status='old', form='formatted',iostat=ierr)
- if( ierr /= 0 ) then
+ status='old', form='formatted',iostat=ier)
+ if( ier /= 0 ) then
nspec2D_bottom = 0
else
read(98,*) nspec2D_bottom
endif
- allocate(ibelm_bottom(nspec2D_bottom))
- allocate(nodes_ibelm_bottom(4,nspec2D_bottom))
+ allocate(ibelm_bottom(nspec2D_bottom),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ibelm_bottom'
+ allocate(nodes_ibelm_bottom(4,nspec2D_bottom),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array nodes_ibelm_bottom'
do ispec2D = 1,nspec2D_bottom
! format: #id_(element containing the face) #id_node1_face .. #id_node4_face
read(98,*) ibelm_bottom(ispec2D), nodes_ibelm_bottom(1,ispec2D), nodes_ibelm_bottom(2,ispec2D), &
@@ -477,14 +494,16 @@
! reads in free_surface boundary files
open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/free_surface_file', &
- status='old', form='formatted',iostat=ierr)
- if( ierr /= 0 ) then
+ status='old', form='formatted',iostat=ier)
+ if( ier /= 0 ) then
nspec2D_top = 0
else
read(98,*) nspec2D_top
endif
- allocate(ibelm_top(nspec2D_top))
- allocate(nodes_ibelm_top(4,nspec2D_top))
+ allocate(ibelm_top(nspec2D_top),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ibelm_top'
+ allocate(nodes_ibelm_top(4,nspec2D_top),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array nodes_ibelm_top'
do ispec2D = 1,nspec2D_top
! format: #id_(element containing the face) #id_node1_face .. #id_node4_face
read(98,*) ibelm_top(ispec2D), nodes_ibelm_top(1,ispec2D), nodes_ibelm_top(2,ispec2D), &
@@ -495,14 +514,16 @@
! reads in moho_surface boundary files (optional)
open(unit=98, file=localpath_name(1:len_trim(localpath_name))//'/moho_surface_file', &
- status='old', form='formatted',iostat=ierr)
- if( ierr /= 0 ) then
+ status='old', form='formatted',iostat=ier)
+ if( ier /= 0 ) then
nspec2D_moho = 0
else
read(98,*) nspec2D_moho
endif
- allocate(ibelm_moho(nspec2D_moho))
- allocate(nodes_ibelm_moho(4,nspec2D_moho))
+ allocate(ibelm_moho(nspec2D_moho),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ibelm_moho'
+ allocate(nodes_ibelm_moho(4,nspec2D_moho),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array nodes_ibelm_moho'
do ispec2D = 1,nspec2D_moho
! format: #id_(element containing the face) #id_node1_face .. #id_node4_face
read(98,*) ibelm_moho(ispec2D), nodes_ibelm_moho(1,ispec2D), nodes_ibelm_moho(2,ispec2D), &
@@ -519,8 +540,10 @@
subroutine check_valence
- allocate(mask_nodes_elmnts(nnodes))
- allocate(used_nodes_elmnts(nnodes))
+ allocate(mask_nodes_elmnts(nnodes),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array mask_nodes_elmnts'
+ allocate(used_nodes_elmnts(nnodes),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array used_nodes_elmnts'
mask_nodes_elmnts(:) = .false.
used_nodes_elmnts(:) = 0
do ispec = 1, nspec
@@ -549,14 +572,21 @@
subroutine scotch_partitioning
implicit none
-
+ ! local parameters
+ integer, dimension(:),allocatable :: num_material
+ integer :: ier
+
elmnts(:,:) = elmnts(:,:) - 1
! determines maximum neighbors based on 1 common node
- allocate(xadj(1:nspec+1))
- allocate(adjncy(1:sup_neighbour*nspec))
- allocate(nnodes_elmnts(1:nnodes))
- allocate(nodes_elmnts(1:nsize*nnodes))
+ allocate(xadj(1:nspec+1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array xadj'
+ allocate(adjncy(1:sup_neighbour*nspec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array adjncy'
+ allocate(nnodes_elmnts(1:nnodes),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array nnodes_elmnts'
+ allocate(nodes_elmnts(1:nsize*nnodes),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array nodes_elmnts'
call mesh2dual_ncommonnodes(nspec, nnodes, nsize, sup_neighbour, elmnts, xadj, adjncy, nnodes_elmnts, &
nodes_elmnts, max_neighbour, 1)
print*, 'mesh2dual: '
@@ -565,20 +595,29 @@
nb_edges = xadj(nspec+1)
! allocates & initializes partioning of elements
- allocate(part(1:nspec))
+ allocate(part(1:nspec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array part'
part(:) = -1
! initializes
! elements load array
- allocate(elmnts_load(1:nspec))
+ allocate(elmnts_load(1:nspec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array elmnts_load'
! uniform load
elmnts_load(:) = 1
+ ! gets materials id associations
+ allocate(num_material(1:nspec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array num_material'
+ num_material(:) = mat(1,:)
+
! in case of acoustic/elastic simulation, weights elements accordingly
call acoustic_elastic_load(elmnts_load,nspec,count_def_mat,count_undef_mat, &
- mat(1,:),mat_prop,undef_mat_prop)
+ num_material,mat_prop,undef_mat_prop)
+ deallocate(num_material)
+
! SCOTCH partitioning
! we use default strategy for partitioning, thus omit specifing explicit strategy .
@@ -604,18 +643,18 @@
!!as your hand-made strategy did not make use of the
!!multi-level framework.
- call scotchfstratinit (scotchstrat(1), ierr)
- if (ierr /= 0) then
+ call scotchfstratinit (scotchstrat(1), ier)
+ if (ier /= 0) then
stop 'ERROR : MAIN : Cannot initialize strat'
endif
- !call scotchfstratgraphmap (scotchstrat(1), trim(scotch_strategy), ierr)
- ! if (ierr /= 0) then
+ !call scotchfstratgraphmap (scotchstrat(1), trim(scotch_strategy), ier)
+ ! if (ier /= 0) then
! stop 'ERROR : MAIN : Cannot build strat'
!endif
- call scotchfgraphinit (scotchgraph (1), ierr)
- if (ierr /= 0) then
+ call scotchfgraphinit (scotchgraph (1), ier)
+ if (ier /= 0) then
stop 'ERROR : MAIN : Cannot initialize graph'
endif
@@ -629,37 +668,37 @@
xadj (1), xadj (1), &
elmnts_load (1), xadj (1), &
nb_edges, adjncy (1), &
- adjncy (1), ierr)
+ adjncy (1), ier)
! w/out element load, but adjacency array
!call scotchfgraphbuild (scotchgraph (1), 0, nspec, &
! xadj (1), xadj (1), &
! xadj (1), xadj (1), &
! nb_edges, adjncy (1), &
- ! adjncy (1), ierr)
+ ! adjncy (1), ier)
- if (ierr /= 0) then
+ if (ier /= 0) then
stop 'ERROR : MAIN : Cannot build graph'
endif
- call scotchfgraphcheck (scotchgraph (1), ierr)
- if (ierr /= 0) then
+ call scotchfgraphcheck (scotchgraph (1), ier)
+ if (ier /= 0) then
stop 'ERROR : MAIN : Invalid check'
endif
- call scotchfgraphpart (scotchgraph (1), nparts, scotchstrat(1),part(1),ierr)
- if (ierr /= 0) then
+ call scotchfgraphpart (scotchgraph (1), nparts, scotchstrat(1),part(1),ier)
+ if (ier /= 0) then
stop 'ERROR : MAIN : Cannot part graph'
endif
- call scotchfgraphexit (scotchgraph (1), ierr)
- if (ierr /= 0) then
+ call scotchfgraphexit (scotchgraph (1), ier)
+ if (ier /= 0) then
stop 'ERROR : MAIN : Cannot destroy graph'
endif
- call scotchfstratexit (scotchstrat(1), ierr)
- if (ierr /= 0) then
+ call scotchfstratexit (scotchstrat(1), ier)
+ if (ier /= 0) then
stop 'ERROR : MAIN : Cannot destroy strat'
endif
@@ -706,8 +745,13 @@
subroutine write_mesh_databases
- allocate(my_interfaces(0:ninterfaces-1))
- allocate(my_nb_interfaces(0:ninterfaces-1))
+ implicit none
+ !local parameters
+
+ allocate(my_interfaces(0:ninterfaces-1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array my_interfaces'
+ allocate(my_nb_interfaces(0:ninterfaces-1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array my_nb_interfaces'
! writes out Database file for each partition
do ipart = 0, nparts-1
@@ -715,8 +759,8 @@
! opens output file
write(prname, "(i6.6,'_Database')") ipart
open(unit=15,file=outputpath_name(1:len_trim(outputpath_name))//'/proc'//prname,&
- status='unknown', action='write', form='formatted', iostat = ierr)
- if( ierr /= 0 ) then
+ status='unknown', action='write', form='formatted', iostat = ier)
+ if( ier /= 0 ) then
print*,'error file open:',outputpath_name(1:len_trim(outputpath_name))//'/proc'//prname
print*
print*,'check if path exists:',outputpath_name(1:len_trim(outputpath_name))
Modified: seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -143,9 +143,11 @@
integer :: num_glob, num_part, nparts
integer, dimension(0:nparts-1) :: num_loc
-
+ integer :: ier
+
! allocates local numbering array
- allocate(glob2loc_elmnts(0:nelmnts-1))
+ allocate(glob2loc_elmnts(0:nelmnts-1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array glob2loc_elmnts'
! initializes number of local elements per partition
do num_part = 0, nparts-1
@@ -189,9 +191,11 @@
integer :: size_glob2loc_nodes,nparts
integer, dimension(0:nparts-1) :: parts_node
integer, dimension(0:nparts-1) :: num_parts
+ integer :: ier
+
+ allocate(glob2loc_nodes_nparts(0:nnodes),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array glob2loc_nodes_nparts'
- allocate(glob2loc_nodes_nparts(0:nnodes))
-
size_glob2loc_nodes = 0
parts_node(:) = 0
@@ -212,8 +216,10 @@
glob2loc_nodes_nparts(nnodes) = size_glob2loc_nodes
- allocate(glob2loc_nodes_parts(0:glob2loc_nodes_nparts(nnodes)-1))
- allocate(glob2loc_nodes(0:glob2loc_nodes_nparts(nnodes)-1))
+ allocate(glob2loc_nodes_parts(0:glob2loc_nodes_nparts(nnodes)-1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array glob2loc_nodes_parts'
+ allocate(glob2loc_nodes(0:glob2loc_nodes_nparts(nnodes)-1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array glob2loc_nodes'
glob2loc_nodes(0) = 0
@@ -269,9 +275,10 @@
integer, intent(in) :: nparts
! local parameters
- integer :: num_part, num_part_bis, el, el_adj, num_interface, num_edge, ncommon_nodes, &
+ integer :: num_part, num_part_bis, el, el_adj, num_interface, num_edge, ncommon_nodes, &
num_node, num_node_bis
- integer :: i, j
+ integer :: i, j
+ integer :: ier
! counts number of interfaces between partitions
ninterfaces = 0
@@ -281,7 +288,8 @@
end do
end do
- allocate(tab_size_interfaces(0:ninterfaces))
+ allocate(tab_size_interfaces(0:ninterfaces),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array tab_size_interfaces'
tab_size_interfaces(:) = 0
num_interface = 0
@@ -316,7 +324,8 @@
num_interface = 0
num_edge = 0
- allocate(tab_interfaces(0:(tab_size_interfaces(ninterfaces)*7-1)))
+ allocate(tab_interfaces(0:(tab_size_interfaces(ninterfaces)*7-1)),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array tab_interfaces'
tab_interfaces(:) = 0
do num_part = 0, nparts-1
@@ -390,6 +399,7 @@
num_node, num_node_bis
integer :: i, j
logical :: is_acoustic_el, is_acoustic_el_adj
+ integer :: ier
! counts number of interfaces between partitions
ninterfaces = 0
@@ -399,7 +409,8 @@
end do
end do
- allocate(tab_size_interfaces(0:ninterfaces))
+ allocate(tab_size_interfaces(0:ninterfaces),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array tab_size_interfaces'
tab_size_interfaces(:) = 0
num_interface = 0
@@ -455,7 +466,8 @@
num_interface = 0
num_edge = 0
- allocate(tab_interfaces(0:(tab_size_interfaces(ninterfaces)*7-1)))
+ allocate(tab_interfaces(0:(tab_size_interfaces(ninterfaces)*7-1)),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array tab_interfaces'
tab_interfaces(:) = 0
do num_part = 0, nparts-1
@@ -1301,10 +1313,10 @@
integer, dimension(:), allocatable :: nodes_elmnts
integer :: max_neighbour
- integer :: i, iface
+ integer :: i, iface, ier
integer :: el, el_adj
logical :: is_repartitioned
-
+
! sets acoustic/elastic flags for materials
is_acoustic(:) = .false.
is_elastic(:) = .false.
@@ -1318,10 +1330,14 @@
enddo
! gets neighbors by 4 common nodes (face)
- allocate(xadj(0:nelmnts))
- allocate(adjncy(0:sup_neighbour*nelmnts-1))
- allocate(nnodes_elmnts(0:nnodes-1))
- allocate(nodes_elmnts(0:nsize*nnodes-1))
+ allocate(xadj(0:nelmnts),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array xadj'
+ allocate(adjncy(0:sup_neighbour*nelmnts-1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array adjncy'
+ allocate(nnodes_elmnts(0:nnodes-1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array nnodes_elmnts'
+ allocate(nodes_elmnts(0:nsize*nnodes-1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array nodes_elmnts'
!call mesh2dual_ncommonnodes(nelmnts, nnodes, elmnts, xadj, adjncy, nnodes_elmnts, nodes_elmnts,4)
call mesh2dual_ncommonnodes(nelmnts, nnodes, nsize, sup_neighbour, &
elmnts, xadj, adjncy, nnodes_elmnts, &
@@ -1340,7 +1356,8 @@
enddo
! coupled elements
- allocate(faces_coupled(2,nfaces_coupled))
+ allocate(faces_coupled(2,nfaces_coupled),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array faces_coupled'
! stores elements indices
nfaces_coupled = 0
@@ -1423,10 +1440,16 @@
integer :: i, j, iface, inode, ispec2D, counter
integer :: el, el_adj
logical :: is_repartitioned
+ integer :: ier
! temporary flag arrays
- allocate( is_moho(0:nelmnts-1)) ! element ids start from 0
- allocate( node_is_moho(0:nnodes-1) ) ! node ids start from 0
+ ! element ids start from 0
+ allocate( is_moho(0:nelmnts-1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array is_moho'
+ ! node ids start from 0
+ allocate( node_is_moho(0:nnodes-1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array node_is_moho'
+
is_moho(:) = .false.
node_is_moho(:) = .false.
@@ -1468,10 +1491,16 @@
print*,' moho elements = ',counter
! gets neighbors by 4 common nodes (face)
- allocate(xadj(0:nelmnts)) ! contains number of adjacent elements (neighbours)
- allocate(adjncy(0:sup_neighbour*nelmnts-1)) ! contains all element id indices of adjacent elements
- allocate(nnodes_elmnts(0:nnodes-1))
- allocate(nodes_elmnts(0:nsize*nnodes-1))
+ ! contains number of adjacent elements (neighbours)
+ allocate(xadj(0:nelmnts),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array xadj'
+ ! contains all element id indices of adjacent elements
+ allocate(adjncy(0:sup_neighbour*nelmnts-1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array adjncy'
+ allocate(nnodes_elmnts(0:nnodes-1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array nnodes_elmnts'
+ allocate(nodes_elmnts(0:nsize*nnodes-1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array nodes_elmnts'
call mesh2dual_ncommonnodes(nelmnts, nnodes, nsize, sup_neighbour, &
elmnts, xadj, adjncy, nnodes_elmnts, &
@@ -1489,7 +1518,8 @@
enddo
! coupled elements
- allocate(faces_coupled(2,nfaces_coupled))
+ allocate(faces_coupled(2,nfaces_coupled),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array faces_coupled'
! stores elements indices
nfaces_coupled = 0
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_mass_matrices.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_mass_matrices.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_mass_matrices.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -44,8 +44,8 @@
integer :: ispec,i,j,k,iglob,ier
! allocates memory
- allocate(rmass(nglob),stat=ier); if(ier /= 0) stop 'error in allocate'
- allocate(rmass_acoustic(nglob),stat=ier); if(ier /= 0) stop 'error in allocate'
+ allocate(rmass(nglob),stat=ier); if(ier /= 0) stop 'error in allocate rmass'
+ allocate(rmass_acoustic(nglob),stat=ier); if(ier /= 0) stop 'error in allocate rmass_acoustic'
allocate(rmass_solid_poroelastic(nglob),stat=ier); if(ier /= 0) stop 'error in allocate'
allocate(rmass_fluid_poroelastic(nglob),stat=ier); if(ier /= 0) stop 'error in allocate'
@@ -160,13 +160,14 @@
double precision :: long_corner,lat_corner,ratio_xi,ratio_eta
integer :: ix_oceans,iy_oceans,iz_oceans,ispec_oceans,ispec2D,igll,iglobnum
integer :: icornerlong,icornerlat
-
+ integer :: ier
! creates ocean load mass matrix
if(OCEANS) then
! adding ocean load mass matrix at ocean bottom
NGLOB_OCEAN = nglob
- allocate(rmass_ocean_load(NGLOB_OCEAN))
+ allocate(rmass_ocean_load(NGLOB_OCEAN),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array rmass_ocean_load'
! create ocean load mass matrix for degrees of freedom at ocean bottom
rmass_ocean_load(:) = 0._CUSTOM_REAL
@@ -265,7 +266,8 @@
! allocate dummy array if no oceans
NGLOB_OCEAN = 1
- allocate(rmass_ocean_load(NGLOB_OCEAN))
+ allocate(rmass_ocean_load(NGLOB_OCEAN),stat=ier)
+ if( ier /= 0 ) stop 'error allocating dummy array rmass_ocean_load'
endif
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -455,6 +455,7 @@
! call sync_all()
allocate( xelm(NGNOD),yelm(NGNOD),zelm(NGNOD),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array xelm etc.'
allocate( qmu_attenuation_store(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
if(ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays')
@@ -463,25 +464,30 @@
call create_name_database(prname,myrank,LOCAL_PATH)
! Gauss-Lobatto-Legendre points of integration
- allocate(xigll(NGLLX),yigll(NGLLY),zigll(NGLLZ))
+ allocate(xigll(NGLLX),yigll(NGLLY),zigll(NGLLZ),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array xigll etc.'
! Gauss-Lobatto-Legendre weights of integration
- allocate(wxgll(NGLLX),wygll(NGLLY),wzgll(NGLLZ))
+ allocate(wxgll(NGLLX),wygll(NGLLY),wzgll(NGLLZ),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array wxgll etc.'
! 3D shape functions and their derivatives
allocate(shape3D(NGNOD,NGLLX,NGLLY,NGLLZ), &
dershape3D(NDIM,NGNOD,NGLLX,NGLLY,NGLLZ),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array shape3D etc.'
! 2D shape functions and their derivatives
allocate(shape2D_x(NGNOD2D,NGLLY,NGLLZ), &
shape2D_y(NGNOD2D,NGLLX,NGLLZ), &
shape2D_bottom(NGNOD2D,NGLLX,NGLLY), &
- shape2D_top(NGNOD2D,NGLLX,NGLLY), stat=ier)
+ shape2D_top(NGNOD2D,NGLLX,NGLLY),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array shape2D_x etc.'
allocate(dershape2D_x(NDIM2D,NGNOD2D,NGLLY,NGLLZ), &
dershape2D_y(NDIM2D,NGNOD2D,NGLLX,NGLLZ), &
dershape2D_bottom(NDIM2D,NGNOD2D,NGLLX,NGLLY), &
dershape2D_top(NDIM2D,NGNOD2D,NGLLX,NGLLY),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array dershape2D_x etc.'
allocate(wgllwgll_xy(NGLLX,NGLLY), &
wgllwgll_xz(NGLLX,NGLLZ), &
@@ -760,4 +766,320 @@
enddo
enddo
-end subroutine crm_ext_setup_indexing
+ end subroutine crm_ext_setup_indexing
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+
+ subroutine create_regions_mesh_save_moho( myrank,nglob,nspec, &
+ nspec2D_moho_ext,ibelm_moho,nodes_ibelm_moho, &
+ nodes_coords_ext_mesh,nnodes_ext_mesh,ibool )
+
+ use create_regions_mesh_ext_par
+ implicit none
+
+ integer :: nspec2D_moho_ext
+ integer, dimension(nspec2D_moho_ext) :: ibelm_moho
+ integer, dimension(4,nspec2D_moho_ext) :: nodes_ibelm_moho
+
+ integer :: myrank,nglob,nspec
+
+ ! data from the external mesh
+ integer :: nnodes_ext_mesh
+ double precision, dimension(NDIM,nnodes_ext_mesh) :: nodes_coords_ext_mesh
+
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
+
+! local parameters
+ ! Moho mesh
+ real(CUSTOM_REAL), dimension(:,:,:),allocatable :: normal_moho_top
+ real(CUSTOM_REAL), dimension(:,:,:),allocatable :: normal_moho_bot
+ integer,dimension(:,:,:),allocatable :: ijk_moho_top, ijk_moho_bot
+ integer,dimension(:),allocatable :: ibelm_moho_top, ibelm_moho_bot
+ integer :: NSPEC2D_MOHO
+ logical, dimension(:),allocatable :: is_moho_top, is_moho_bot
+
+ real(kind=CUSTOM_REAL),dimension(NGNOD2D) :: xcoord,ycoord,zcoord
+ real(kind=CUSTOM_REAL) :: jacobian2Dw_face(NGLLX,NGLLY)
+ real(kind=CUSTOM_REAL) :: normal_face(NDIM,NGLLX,NGLLY)
+ real(kind=CUSTOM_REAL),dimension(NDIM):: normal
+ integer :: ijk_face(3,NGLLX,NGLLY)
+
+ real(kind=CUSTOM_REAL),dimension(:,:),allocatable :: iglob_normals
+ integer,dimension(:),allocatable:: iglob_is_surface
+
+ integer :: imoho_bot,imoho_top
+ integer :: ispec2D,ispec,icorner,iface,i,j,k,igll,iglob,ier
+ integer :: iglob_midpoint,idirect,counter
+ integer :: imoho_top_all,imoho_bot_all,imoho_all
+
+ ! corners indices of reference cube faces
+ integer,dimension(3,4),parameter :: iface1_corner_ijk = &
+ reshape( (/ 1,1,1, 1,NGLLY,1, 1,NGLLY,NGLLZ, 1,1,NGLLZ /),(/3,4/)) ! xmin
+ integer,dimension(3,4),parameter :: iface2_corner_ijk = &
+ reshape( (/ NGLLX,1,1, NGLLX,NGLLY,1, NGLLX,NGLLY,NGLLZ, NGLLX,1,NGLLZ /),(/3,4/)) ! xmax
+ integer,dimension(3,4),parameter :: iface3_corner_ijk = &
+ reshape( (/ 1,1,1, 1,1,NGLLZ, NGLLX,1,NGLLZ, NGLLX,1,1 /),(/3,4/)) ! ymin
+ integer,dimension(3,4),parameter :: iface4_corner_ijk = &
+ reshape( (/ 1,NGLLY,1, NGLLX,NGLLY,1, NGLLX,NGLLY,NGLLZ, 1,NGLLY,NGLLZ /),(/3,4/)) ! ymax
+ integer,dimension(3,4),parameter :: iface5_corner_ijk = &
+ reshape( (/ 1,1,1, 1,NGLLY,1, NGLLX,NGLLY,1, NGLLX,1,1 /),(/3,4/)) ! bottom
+ integer,dimension(3,4),parameter :: iface6_corner_ijk = &
+ reshape( (/ 1,1,NGLLZ, NGLLX,1,NGLLZ, NGLLX,NGLLY,NGLLZ, 1,NGLLY,NGLLZ /),(/3,4/)) ! top
+ integer,dimension(3,4,6),parameter :: iface_all_corner_ijk = &
+ reshape( (/ iface1_corner_ijk,iface2_corner_ijk, &
+ iface3_corner_ijk,iface4_corner_ijk, &
+ iface5_corner_ijk,iface6_corner_ijk /),(/3,4,6/)) ! all faces
+ ! midpoint indices for each face (xmin,xmax,ymin,ymax,zmin,zmax)
+ integer,dimension(3,6),parameter :: iface_all_midpointijk = &
+ reshape( (/ 1,2,2, NGLLX,2,2, 2,1,2, 2,NGLLY,2, 2,2,1, 2,2,NGLLZ /),(/3,6/)) ! top
+
+ ! temporary arrays for passing information
+ allocate(iglob_is_surface(nglob), &
+ iglob_normals(NDIM,nglob),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array iglob_is_surface'
+
+ iglob_is_surface = 0
+ iglob_normals = 0._CUSTOM_REAL
+
+ ! loops over given moho surface elements
+ do ispec2D=1, nspec2D_moho_ext
+
+ ! gets element id
+ ispec = ibelm_moho(ispec2D)
+
+ ! looks for i,j,k indices of GLL points on boundary face
+ ! determines element face by given CUBIT corners
+ ! (note: uses point locations rather than point indices to find the element face,
+ ! because the indices refer no more to the newly indexed ibool array )
+ do icorner=1,NGNOD2D
+ xcoord(icorner) = nodes_coords_ext_mesh(1,nodes_ibelm_moho(icorner,ispec2D))
+ ycoord(icorner) = nodes_coords_ext_mesh(2,nodes_ibelm_moho(icorner,ispec2D))
+ zcoord(icorner) = nodes_coords_ext_mesh(3,nodes_ibelm_moho(icorner,ispec2D))
+ enddo
+
+ ! sets face id of reference element associated with this face
+ call get_element_face_id(ispec,xcoord,ycoord,zcoord, &
+ ibool,nspec,nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy, &
+ iface)
+
+ ! ijk indices of GLL points for face id
+ call get_element_face_gll_indices(iface,ijk_face,NGLLX,NGLLZ)
+
+ ! weighted jacobian and normal
+ call get_jacobian_boundary_face(myrank,nspec, &
+ xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob,&
+ dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,&
+ ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLZ)
+
+ ! normal convention: points away from element
+ ! switch normal direction if necessary
+ do j=1,NGLLY
+ do i=1,NGLLX
+ call get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, &
+ ibool,nspec,nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy, &
+ normal_face(:,i,j) )
+ enddo
+ enddo
+
+ ! stores information on global points on moho surface
+ igll = 0
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(ijk_face(1,i,j),ijk_face(2,i,j),ijk_face(3,i,j),ispec)
+ ! sets flag
+ iglob_is_surface(iglob) = ispec2D
+ ! sets normals
+ iglob_normals(:,iglob) = normal_face(:,i,j)
+ enddo
+ enddo
+ enddo
+
+ ! stores moho elements
+ NSPEC2D_MOHO = nspec2D_moho_ext
+
+ allocate(ibelm_moho_bot(NSPEC2D_MOHO), &
+ ibelm_moho_top(NSPEC2D_MOHO), &
+ normal_moho_top(NDIM,NGLLSQUARE,NSPEC2D_MOHO), &
+ normal_moho_bot(NDIM,NGLLSQUARE,NSPEC2D_MOHO), &
+ ijk_moho_bot(3,NGLLSQUARE,NSPEC2D_MOHO), &
+ ijk_moho_top(3,NGLLSQUARE,NSPEC2D_MOHO),stat=ier)
+ if( ier /= 0 ) stop 'error allocating ibelm_moho_bot'
+
+ ibelm_moho_bot = 0
+ ibelm_moho_top = 0
+
+ ! element flags
+ allocate(is_moho_top(nspec), &
+ is_moho_bot(nspec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating is_moho_top'
+ is_moho_top = .false.
+ is_moho_bot = .false.
+
+ ! finds spectral elements with moho surface
+ imoho_top = 0
+ imoho_bot = 0
+ do ispec=1,nspec
+
+ ! loops over each face
+ do iface = 1,6
+ ! checks if corners of face on surface
+ counter = 0
+ do icorner = 1,NGNOD2D
+ i = iface_all_corner_ijk(1,icorner,iface)
+ j = iface_all_corner_ijk(2,icorner,iface)
+ k = iface_all_corner_ijk(3,icorner,iface)
+ iglob = ibool(i,j,k,ispec)
+
+ ! checks if point on surface
+ if( iglob_is_surface(iglob) > 0 ) then
+ counter = counter+1
+
+ ! reference corner coordinates
+ xcoord(icorner) = xstore_dummy(iglob)
+ ycoord(icorner) = ystore_dummy(iglob)
+ zcoord(icorner) = zstore_dummy(iglob)
+ endif
+ enddo
+
+ ! stores moho informations
+ if( counter == NGNOD2D ) then
+
+ ! gets face GLL points i,j,k indices from element face
+ call get_element_face_gll_indices(iface,ijk_face,NGLLX,NGLLY)
+
+ ! re-computes face infos
+ ! weighted jacobian and normal
+ call get_jacobian_boundary_face(myrank,nspec, &
+ xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob,&
+ dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,&
+ ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLZ)
+
+ ! normal convention: points away from element
+ ! switch normal direction if necessary
+ do j=1,NGLLZ
+ do i=1,NGLLX
+ call get_element_face_normal(ispec,iface,xcoord,ycoord,zcoord, &
+ ibool,nspec,nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy, &
+ normal_face(:,i,j) )
+ enddo
+ enddo
+
+ ! takes normal stored temporary on a face midpoint
+ i = iface_all_midpointijk(1,iface)
+ j = iface_all_midpointijk(2,iface)
+ k = iface_all_midpointijk(3,iface)
+ iglob_midpoint = ibool(i,j,k,ispec)
+ normal(:) = iglob_normals(:,iglob_midpoint)
+
+ ! determines whether normal points into element or not (top/bottom distinction)
+ call get_element_face_normal_idirect(ispec,iface,xcoord,ycoord,zcoord, &
+ ibool,nspec,nglob, &
+ xstore_dummy,ystore_dummy,zstore_dummy, &
+ normal,idirect )
+
+ ! takes moho surface element id given by id on midpoint
+ ispec2D = iglob_is_surface(iglob_midpoint)
+
+ ! sets face infos for bottom (normal points away from element)
+ if( idirect == 1 ) then
+
+ ! checks validity
+ if( is_moho_bot( ispec) .eqv. .true. ) then
+ print*,'error: moho surface geometry bottom'
+ print*,' does not allow for mulitple element faces in kernel computation'
+ call exit_mpi(myrank,'error moho bottom elements')
+ endif
+
+ imoho_bot = imoho_bot + 1
+ is_moho_bot(ispec) = .true.
+ ibelm_moho_bot(ispec2D) = ispec
+
+ ! stores on surface gll points (assuming NGLLX = NGLLY = NGLLZ)
+ igll = 0
+ do j=1,NGLLZ
+ do i=1,NGLLX
+ igll = igll+1
+ ijk_moho_bot(:,igll,ispec2D) = ijk_face(:,i,j)
+ normal_moho_bot(:,igll,ispec2D) = normal_face(:,i,j)
+ enddo
+ enddo
+
+ ! sets face infos for top element
+ else if( idirect == 2 ) then
+
+ ! checks validity
+ if( is_moho_top( ispec) .eqv. .true. ) then
+ print*,'error: moho surface geometry top'
+ print*,' does not allow for mulitple element faces kernel computation'
+ call exit_mpi(myrank,'error moho top elements')
+ endif
+
+ imoho_top = imoho_top + 1
+ is_moho_top(ispec) = .true.
+ ibelm_moho_top(ispec2D) = ispec
+
+ ! gll points
+ igll = 0
+ do j=1,NGLLZ
+ do i=1,NGLLX
+ igll = igll+1
+ ijk_moho_top(:,igll,ispec) = ijk_face(:,i,j)
+ ! note: top elements have normal pointing into element
+ normal_moho_top(:,igll,ispec) = - normal_face(:,i,j)
+ enddo
+ enddo
+ endif
+
+ endif ! counter
+
+ enddo ! iface
+
+ ! checks validity of top/bottom distinction
+ if( is_moho_top(ispec) .and. is_moho_bot(ispec) ) then
+ print*,'error: moho surface elements confusing'
+ print*,' element:',ispec,'has top and bottom surface'
+ call exit_mpi(myrank,'error moho surface element')
+ endif
+
+ enddo ! ispec2D
+
+ ! note: surface e.g. could be at the free-surface and have no top elements etc...
+ ! user output
+ call sum_all_i( imoho_top, imoho_top_all )
+ call sum_all_i( imoho_bot, imoho_bot_all )
+ call sum_all_i( NSPEC2D_MOHO, imoho_all )
+ if( myrank == 0 ) then
+ write(IMAIN,*) '********'
+ write(IMAIN,*) 'Moho surface:'
+ write(IMAIN,*) ' total surface elements: ',imoho_all
+ write(IMAIN,*) ' top elements :',imoho_top_all
+ write(IMAIN,*) ' bottom elements:',imoho_bot_all
+ write(IMAIN,*) '********'
+ endif
+
+ ! saves moho files: total number of elements, corner points, all points
+ open(unit=27,file=prname(1:len_trim(prname))//'ibelm_moho.bin',status='unknown',form='unformatted')
+ write(27) NSPEC2D_MOHO
+ write(27) ibelm_moho_top
+ write(27) ibelm_moho_bot
+ write(27) ijk_moho_top
+ write(27) ijk_moho_bot
+ close(27)
+ open(unit=27,file=prname(1:len_trim(prname))//'normal_moho.bin',status='unknown',form='unformatted')
+ write(27) normal_moho_top
+ write(27) normal_moho_bot
+ close(27)
+ open(unit=27,file=prname(1:len_trim(prname))//'is_moho.bin',status='unknown',form='unformatted')
+ write(27) is_moho_top
+ write(27) is_moho_bot
+ close(27)
+
+ end subroutine create_regions_mesh_save_moho
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -491,8 +491,7 @@
use generate_databases_par
implicit none
-
-
+
if( OCEANS .and. TOPOGRAPHY ) then
! for Southern California
@@ -503,7 +502,8 @@
DEGREES_PER_CELL_TOPO = DEGREES_PER_CELL_TOPO_SOCAL
topo_file = TOPO_FILE_SOCAL
- allocate(itopo_bathy(NX_TOPO,NY_TOPO))
+ allocate(itopo_bathy(NX_TOPO,NY_TOPO),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array itopo_bathy'
call read_topo_bathy_file(itopo_bathy,NX_TOPO,NY_TOPO,topo_file)
@@ -515,13 +515,15 @@
else
NX_TOPO = 1
NY_TOPO = 1
- allocate(itopo_bathy(NX_TOPO,NY_TOPO))
+ allocate(itopo_bathy(NX_TOPO,NY_TOPO),stat=ier)
+ if( ier /= 0 ) stop 'error allocating dummy array itopo_bathy'
endif
!! read basement map
! if(BASEMENT_MAP) then
-! call get_value_string(BASEMENT_MAP_FILE,'model.BASEMENT_MAP_FILE','../in_data_files/la_basement/reggridbase2_filtered_ascii.dat')
+! call get_value_string(BASEMENT_MAP_FILE,'model.BASEMENT_MAP_FILE', &
+! '../in_data_files/la_basement/reggridbase2_filtered_ascii.dat')
! open(unit=55,file=BASEMENT_MAP_FILE,status='old',action='read')
! do ix=1,NX_BASEMENT
! do iy=1,NY_BASEMENT
@@ -560,7 +562,8 @@
call exit_mpi(myrank,'error opening database file')
endif
read(IIN,*) nnodes_ext_mesh
- allocate(nodes_coords_ext_mesh(NDIM,nnodes_ext_mesh))
+ allocate(nodes_coords_ext_mesh(NDIM,nnodes_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array nodes_coords_ext_mesh'
do inode = 1, nnodes_ext_mesh
read(IIN,*) dummy_node, nodes_coords_ext_mesh(1,inode), nodes_coords_ext_mesh(2,inode), &
nodes_coords_ext_mesh(3,inode)
@@ -574,7 +577,8 @@
! read materials' physical properties
read(IIN,*) nmat_ext_mesh, nundefMat_ext_mesh
- allocate(materials_ext_mesh(6,nmat_ext_mesh))
+ allocate(materials_ext_mesh(6,nmat_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array materials_ext_mesh'
do imat = 1, nmat_ext_mesh
! format: #(1) rho #(2) vp #(3) vs #(4) Q_flag #(5) anisotropy_flag #(6) material_domain_id
read(IIN,*) materials_ext_mesh(1,imat), materials_ext_mesh(2,imat), materials_ext_mesh(3,imat), &
@@ -590,7 +594,8 @@
endif
call sync_all()
- allocate(undef_mat_prop(6,nundefMat_ext_mesh))
+ allocate(undef_mat_prop(6,nundefMat_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array undef_mat_prop'
do imat = 1, nundefMat_ext_mesh
! format example tomography:
! -1 tomography elastic tomography_model.xyz 1 2
@@ -611,8 +616,10 @@
! element indexing
read(IIN,*) nelmnts_ext_mesh
- allocate(elmnts_ext_mesh(esize,nelmnts_ext_mesh))
- allocate(mat_ext_mesh(2,nelmnts_ext_mesh))
+ allocate(elmnts_ext_mesh(esize,nelmnts_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array elmnts_ext_mesh'
+ allocate(mat_ext_mesh(2,nelmnts_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array mat_ext_mesh'
! reads in material association for each spectral element and corner node indices
do ispec = 1, nelmnts_ext_mesh
@@ -652,32 +659,38 @@
NSPEC2D_BOTTOM = nspec2D_bottom_ext
NSPEC2D_TOP = nspec2D_top_ext
- allocate(ibelm_xmin(nspec2D_xmin),nodes_ibelm_xmin(4,nspec2D_xmin))
+ allocate(ibelm_xmin(nspec2D_xmin),nodes_ibelm_xmin(4,nspec2D_xmin),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ibelm_xmin etc.'
do ispec2D = 1,nspec2D_xmin
read(IIN,*) ibelm_xmin(ispec2D),(nodes_ibelm_xmin(j,ispec2D),j=1,4)
end do
- allocate(ibelm_xmax(nspec2D_xmax),nodes_ibelm_xmax(4,nspec2D_xmax))
+ allocate(ibelm_xmax(nspec2D_xmax),nodes_ibelm_xmax(4,nspec2D_xmax),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ibelm_xmax etc.'
do ispec2D = 1,nspec2D_xmax
read(IIN,*) ibelm_xmax(ispec2D),(nodes_ibelm_xmax(j,ispec2D),j=1,4)
end do
- allocate(ibelm_ymin(nspec2D_ymin),nodes_ibelm_ymin(4,nspec2D_ymin))
+ allocate(ibelm_ymin(nspec2D_ymin),nodes_ibelm_ymin(4,nspec2D_ymin),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ibelm_ymin'
do ispec2D = 1,nspec2D_ymin
read(IIN,*) ibelm_ymin(ispec2D),(nodes_ibelm_ymin(j,ispec2D),j=1,4)
end do
- allocate(ibelm_ymax(nspec2D_ymax),nodes_ibelm_ymax(4,nspec2D_ymax))
+ allocate(ibelm_ymax(nspec2D_ymax),nodes_ibelm_ymax(4,nspec2D_ymax),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ibelm_ymax etc.'
do ispec2D = 1,nspec2D_ymax
read(IIN,*) ibelm_ymax(ispec2D),(nodes_ibelm_ymax(j,ispec2D),j=1,4)
end do
- allocate(ibelm_bottom(nspec2D_bottom_ext),nodes_ibelm_bottom(4,nspec2D_bottom_ext))
+ allocate(ibelm_bottom(nspec2D_bottom_ext),nodes_ibelm_bottom(4,nspec2D_bottom_ext),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ibelm_bottom etc.'
do ispec2D = 1,nspec2D_bottom_ext
read(IIN,*) ibelm_bottom(ispec2D),(nodes_ibelm_bottom(j,ispec2D),j=1,4)
end do
- allocate(ibelm_top(nspec2D_top_ext),nodes_ibelm_top(4,nspec2D_top_ext))
+ allocate(ibelm_top(nspec2D_top_ext),nodes_ibelm_top(4,nspec2D_top_ext),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ibelm_top etc.'
do ispec2D = 1,nspec2D_top_ext
read(IIN,*) ibelm_top(ispec2D),(nodes_ibelm_top(j,ispec2D),j=1,4)
end do
@@ -702,11 +715,16 @@
read(IIN,*) num_interfaces_ext_mesh, max_interface_size_ext_mesh
! allocates interfaces
- allocate(my_neighbours_ext_mesh(num_interfaces_ext_mesh))
- allocate(my_nelmnts_neighbours_ext_mesh(num_interfaces_ext_mesh))
- allocate(my_interfaces_ext_mesh(6,max_interface_size_ext_mesh,num_interfaces_ext_mesh))
- allocate(ibool_interfaces_ext_mesh(NGLLX*NGLLX*max_interface_size_ext_mesh,num_interfaces_ext_mesh))
- allocate(nibool_interfaces_ext_mesh(num_interfaces_ext_mesh))
+ allocate(my_neighbours_ext_mesh(num_interfaces_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array my_neighbours_ext_mesh'
+ allocate(my_nelmnts_neighbours_ext_mesh(num_interfaces_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array my_nelmnts_neighbours_ext_mesh'
+ allocate(my_interfaces_ext_mesh(6,max_interface_size_ext_mesh,num_interfaces_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array my_interfaces_ext_mesh'
+ allocate(ibool_interfaces_ext_mesh(NGLLX*NGLLX*max_interface_size_ext_mesh,num_interfaces_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ibool_interfaces_ext_mesh'
+ allocate(nibool_interfaces_ext_mesh(num_interfaces_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array nibool_interfaces_ext_mesh'
! loops over MPI interfaces with other partitions
do num_interface = 1, num_interfaces_ext_mesh
@@ -756,7 +774,8 @@
if( num_moho == 0 ) call exit_mpi(myrank,'error no moho mesh in database')
! reads in element informations
- allocate(ibelm_moho(nspec2D_moho_ext),nodes_ibelm_moho(4,nspec2D_moho_ext))
+ allocate(ibelm_moho(nspec2D_moho_ext),nodes_ibelm_moho(4,nspec2D_moho_ext),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ibelm_moho etc.'
do ispec2D = 1,nspec2D_moho_ext
! format: #element_id #node_id1 #node_id2 #node_id3 #node_id4
read(IIN,*) ibelm_moho(ispec2D),(nodes_ibelm_moho(j,ispec2D),j=1,4)
@@ -792,9 +811,12 @@
! use dynamic allocation to allocate memory for arrays
! allocate(idoubling(nspec))
- allocate(ibool(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(xstore(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(ystore(NGLLX,NGLLY,NGLLZ,nspec))
+ allocate(ibool(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ibool'
+ allocate(xstore(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array xstore'
+ allocate(ystore(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ystore'
allocate(zstore(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
if(ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays')
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_MPI.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_MPI.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_MPI.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -74,7 +74,7 @@
integer :: num_points1, num_points2
! assembly test
- integer :: i,j,k,ispec,iglob,count,inum
+ integer :: i,j,k,ispec,iglob,count,inum,ier
integer :: max_nibool_interfaces_ext_mesh
integer,dimension(:),allocatable :: test_flag
real(kind=CUSTOM_REAL), dimension(:),allocatable :: test_flag_cr
@@ -91,23 +91,34 @@
ibool_interfaces_ext_mesh, &
nibool_interfaces_ext_mesh )
- allocate(nibool_interfaces_ext_mesh_true(num_interfaces_ext_mesh))
+ allocate(nibool_interfaces_ext_mesh_true(num_interfaces_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array nibool_interfaces_ext_mesh_true'
! sorts ibool comm buffers lexicographically for all MPI interfaces
num_points1 = 0
num_points2 = 0
do iinterface = 1, num_interfaces_ext_mesh
- allocate(xp(nibool_interfaces_ext_mesh(iinterface)))
- allocate(yp(nibool_interfaces_ext_mesh(iinterface)))
- allocate(zp(nibool_interfaces_ext_mesh(iinterface)))
- allocate(locval(nibool_interfaces_ext_mesh(iinterface)))
- allocate(ifseg(nibool_interfaces_ext_mesh(iinterface)))
- allocate(reorder_interface_ext_mesh(nibool_interfaces_ext_mesh(iinterface)))
- allocate(ind_ext_mesh(nibool_interfaces_ext_mesh(iinterface)))
- allocate(ninseg_ext_mesh(nibool_interfaces_ext_mesh(iinterface)))
- allocate(iwork_ext_mesh(nibool_interfaces_ext_mesh(iinterface)))
- allocate(work_ext_mesh(nibool_interfaces_ext_mesh(iinterface)))
+ allocate(xp(nibool_interfaces_ext_mesh(iinterface)),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array xp'
+ allocate(yp(nibool_interfaces_ext_mesh(iinterface)),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array yp'
+ allocate(zp(nibool_interfaces_ext_mesh(iinterface)),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array zp'
+ allocate(locval(nibool_interfaces_ext_mesh(iinterface)),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array locval'
+ allocate(ifseg(nibool_interfaces_ext_mesh(iinterface)),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ifseg'
+ allocate(reorder_interface_ext_mesh(nibool_interfaces_ext_mesh(iinterface)),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array reorder_interface_ext_mesh'
+ allocate(ind_ext_mesh(nibool_interfaces_ext_mesh(iinterface)),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ind_ext_mesh'
+ allocate(ninseg_ext_mesh(nibool_interfaces_ext_mesh(iinterface)),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ninseg_ext_mesh'
+ allocate(iwork_ext_mesh(nibool_interfaces_ext_mesh(iinterface)),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array iwork_ext_mesh'
+ allocate(work_ext_mesh(nibool_interfaces_ext_mesh(iinterface)),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array work_ext_mesh'
! gets x,y,z coordinates of global points on MPI interface
do ilocnum = 1, nibool_interfaces_ext_mesh(iinterface)
@@ -158,7 +169,8 @@
endif
! checks with assembly of test fields
- allocate(test_flag(nglob),test_flag_cr(nglob))
+ allocate(test_flag(nglob),test_flag_cr(nglob),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array test_flag etc.'
test_flag(:) = 0
test_flag_cr(:) = 0._CUSTOM_REAL
count = 0
@@ -185,7 +197,8 @@
! collects contributions from different MPI partitions
! sets up MPI communications
max_nibool_interfaces_ext_mesh = maxval( nibool_interfaces_ext_mesh(:) )
- allocate(ibool_interfaces_dummy(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh))
+ allocate(ibool_interfaces_dummy(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ibool_interfaces_dummy'
count = 0
do iinterface = 1, num_interfaces_ext_mesh
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_coupling_surfaces.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_coupling_surfaces.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_coupling_surfaces.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -97,10 +97,14 @@
!character(len=256):: prname_file
! allocates temporary arrays
- allocate(tmp_normal(NDIM,NGLLSQUARE,nspec*6))
- allocate(tmp_jacobian2Dw(NGLLSQUARE,nspec*6))
- allocate(tmp_ijk(3,NGLLSQUARE,nspec*6))
- allocate(tmp_ispec(nspec*6))
+ allocate(tmp_normal(NDIM,NGLLSQUARE,nspec*6),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array tmp_normal'
+ allocate(tmp_jacobian2Dw(NGLLSQUARE,nspec*6),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array tmp_jacobian2Dw'
+ allocate(tmp_ijk(3,NGLLSQUARE,nspec*6),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array tmp_ijk'
+ allocate(tmp_ispec(nspec*6),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array tmp_ispec'
tmp_ispec(:) = 0
tmp_ijk(:,:,:) = 0
tmp_normal(:,:,:) = 0.0
@@ -108,10 +112,13 @@
! sets flags for acoustic / elastic on global points
allocate(elastic_flag(nglob),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array elastic_flag'
allocate(acoustic_flag(nglob),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array acoustic_flag'
allocate(test_flag(nglob),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array test_flag'
allocate(mask_ibool(nglob),stat=ier)
- if( ier /= 0 ) stop 'error allocate flag array'
+ if( ier /= 0 ) stop 'error allocating array mask_ibool'
elastic_flag(:) = 0
acoustic_flag(:) = 0
test_flag(:) = 0
@@ -151,7 +158,7 @@
! sets up MPI communications
max_nibool_interfaces_ext_mesh = maxval( nibool_interfaces_ext_mesh(:) )
allocate(ibool_interfaces_ext_mesh_dummy(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh),stat=ier)
- if( ier /= 0 ) stop 'error allocating array'
+ if( ier /= 0 ) stop 'error allocating array ibool_interfaces_ext_mesh_dummy'
do i = 1, num_interfaces_ext_mesh
ibool_interfaces_ext_mesh_dummy(:,i) = ibool_interfaces_ext_mesh(1:max_nibool_interfaces_ext_mesh,i)
enddo
@@ -288,10 +295,14 @@
! note: no need to store material parameters on these coupling points
! for acoustic-elastic interface
num_coupling_ac_el_faces = inum
- allocate(coupling_ac_el_normal(NDIM,NGLLSQUARE,num_coupling_ac_el_faces))
- allocate(coupling_ac_el_jacobian2Dw(NGLLSQUARE,num_coupling_ac_el_faces))
- allocate(coupling_ac_el_ijk(3,NGLLSQUARE,num_coupling_ac_el_faces))
- allocate(coupling_ac_el_ispec(num_coupling_ac_el_faces))
+ allocate(coupling_ac_el_normal(NDIM,NGLLSQUARE,num_coupling_ac_el_faces),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array coupling_ac_el_normal'
+ allocate(coupling_ac_el_jacobian2Dw(NGLLSQUARE,num_coupling_ac_el_faces),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array coupling_ac_el_jacobian2Dw'
+ allocate(coupling_ac_el_ijk(3,NGLLSQUARE,num_coupling_ac_el_faces),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array coupling_ac_el_ijk'
+ allocate(coupling_ac_el_ispec(num_coupling_ac_el_faces),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array coupling_ac_el_ispec'
do inum = 1,num_coupling_ac_el_faces
coupling_ac_el_normal(:,:,inum) = tmp_normal(:,:,inum)
coupling_ac_el_jacobian2Dw(:,inum) = tmp_jacobian2Dw(:,inum)
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_model.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_model.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_model.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -52,7 +52,7 @@
real(kind=CUSTOM_REAL) :: vp,vs,rho,qmu_atten
real(kind=CUSTOM_REAL) :: c11,c12,c13,c14,c15,c16,c22,c23,c24,c25, &
c26,c33,c34,c35,c36,c44,c45,c46,c55,c56,c66
- integer :: ispec,i,j,k,iundef
+ integer :: ispec,i,j,k,iundef,ier
integer :: iflag,flag_below,flag_above
integer :: iflag_aniso,idomain_id,imaterial_id
@@ -304,21 +304,24 @@
! processors name
write(prname_lp,'(a,i6.6,a)') trim(LOCAL_PATH)//'proc',myrank,'_'
- allocate( rho_read(NGLLX,NGLLY,NGLLZ,nspec) )
+ allocate( rho_read(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array rho_read'
write(prname_lp,'(a,i6.6,a)') trim(LOCAL_PATH)//'proc',myrank,'_'
open(unit=28,file=prname_lp(1:len_trim(prname_lp))//'rho.bin',&
status='unknown',action='read',form='unformatted')
read(28) rho_read
close(28)
- allocate( vp_read(NGLLX,NGLLY,NGLLZ,nspec) )
+ allocate( vp_read(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array vp_read'
write(prname_lp,'(a,i6.6,a)') trim(LOCAL_PATH)//'proc',myrank,'_'
open(unit=28,file=prname_lp(1:len_trim(prname_lp))//'vp.bin',&
status='unknown',action='read',form='unformatted')
read(28) vp_read
close(28)
- allocate( vs_read(NGLLX,NGLLY,NGLLZ,nspec) )
+ allocate( vs_read(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array vs_read'
write(prname_lp,'(a,i6.6,a)') trim(LOCAL_PATH)//'proc',myrank,'_'
open(unit=28,file=prname_lp(1:len_trim(prname_lp))//'vs.bin',&
status='unknown',action='read',form='unformatted')
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/model_tomography.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/model_tomography.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/model_tomography.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -81,7 +81,10 @@
!if(myrank == 0) call read_external_model()
! broadcast the information read on the master to the nodes, e.g.
!call MPI_BCAST(nrecord,1,MPI_INTEGER,0,MPI_COMM_WORLD,ier)
- !if( myrank /= 0 ) allocate( vp_tomography(1:nrecord) )
+ !if( myrank /= 0 ) then
+ ! allocate( vp_tomography(1:nrecord) ,stat=ier)
+ ! if( ier /= 0 ) stop 'error allocating array vp_tomography'
+ !endif
!call MPI_BCAST(vp_tomography,size(vp_tomography),CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
end subroutine model_tomography_broadcast
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/save_arrays_solver.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -267,22 +267,26 @@
! mesh arrays used for example in combine_vol_data.f90
!--- x coordinate
- open(unit=27,file=prname(1:len_trim(prname))//'x.bin',status='unknown',form='unformatted')
+ open(unit=27,file=prname(1:len_trim(prname))//'x.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file x.bin'
write(27) xstore_dummy
close(27)
!--- y coordinate
- open(unit=27,file=prname(1:len_trim(prname))//'y.bin',status='unknown',form='unformatted')
+ open(unit=27,file=prname(1:len_trim(prname))//'y.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file y.bin'
write(27) ystore_dummy
close(27)
!--- z coordinate
- open(unit=27,file=prname(1:len_trim(prname))//'z.bin',status='unknown',form='unformatted')
+ open(unit=27,file=prname(1:len_trim(prname))//'z.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file z.bin'
write(27) zstore_dummy
close(27)
! ibool
- open(unit=27,file=prname(1:len_trim(prname))//'ibool.bin',status='unknown',form='unformatted')
+ open(unit=27,file=prname(1:len_trim(prname))//'ibool.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file ibool.bin'
write(27) ibool
close(27)
@@ -297,7 +301,8 @@
!endif
v_tmp = 0.0
where( rho_vp /= 0._CUSTOM_REAL ) v_tmp = (FOUR_THIRDS * mustore + kappastore) / rho_vp
- open(unit=27,file=prname(1:len_trim(prname))//'vp.bin',status='unknown',form='unformatted')
+ open(unit=27,file=prname(1:len_trim(prname))//'vp.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file vp.bin'
write(27) v_tmp
close(27)
@@ -318,7 +323,8 @@
!endif
v_tmp = 0.0
where( rho_vs /= 0._CUSTOM_REAL ) v_tmp = mustore / rho_vs
- open(unit=27,file=prname(1:len_trim(prname))//'vs.bin',status='unknown',form='unformatted')
+ open(unit=27,file=prname(1:len_trim(prname))//'vs.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file vs.bin'
write(27) v_tmp
close(27)
@@ -332,7 +338,8 @@
! outputs density model for check
v_tmp = 0.0
where( rho_vp /= 0._CUSTOM_REAL ) v_tmp = rho_vp**2 / (FOUR_THIRDS * mustore + kappastore)
- open(unit=27,file=prname(1:len_trim(prname))//'rho.bin',status='unknown',form='unformatted')
+ open(unit=27,file=prname(1:len_trim(prname))//'rho.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file rho.bin'
write(27) v_tmp
close(27)
@@ -347,7 +354,8 @@
! acoustic-elastic domains
if( ACOUSTIC_SIMULATION .and. ELASTIC_SIMULATION ) then
! saves points on acoustic-elastic coupling interface
- allocate( iglob_tmp(NGLLSQUARE*num_coupling_ac_el_faces))
+ allocate( iglob_tmp(NGLLSQUARE*num_coupling_ac_el_faces),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array iglob_tmp'
inum = 0
iglob_tmp(:) = 0
do i=1,num_coupling_ac_el_faces
@@ -366,7 +374,8 @@
filename)
! saves acoustic/elastic flag
- allocate(v_tmp_i(nspec))
+ allocate(v_tmp_i(nspec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array v_tmp_i'
do i=1,nspec
if( ispec_is_acoustic(i) ) then
v_tmp_i(i) = 1
Modified: seismo/3D/SPECFEM3D/trunk/src/meshfem3D/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/meshfem3D/create_regions_mesh.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/meshfem3D/create_regions_mesh.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -166,33 +166,46 @@
call create_name_database(prname,myrank,LOCAL_PATH)
! flag indicating whether point is in the sediments
- allocate(flag_sediments(NGLLX,NGLLY,NGLLZ,nspec))
+ allocate(flag_sediments(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array flag_sediments'
allocate(not_fully_in_bedrock(nspec),stat=ier)
if( ier /= 0 ) call exit_MPI(myrank,'not enough memory to allocate arrays')
! boundary locator
- allocate(iboun(6,nspec))
+ allocate(iboun(6,nspec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array iboun'
! boundary parameters locator
- allocate(ibelm_xmin(NSPEC2DMAX_XMIN_XMAX))
- allocate(ibelm_xmax(NSPEC2DMAX_XMIN_XMAX))
- allocate(ibelm_ymin(NSPEC2DMAX_YMIN_YMAX))
- allocate(ibelm_ymax(NSPEC2DMAX_YMIN_YMAX))
- allocate(ibelm_bottom(NSPEC2D_BOTTOM))
+ allocate(ibelm_xmin(NSPEC2DMAX_XMIN_XMAX),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ibelm_xmin'
+ allocate(ibelm_xmax(NSPEC2DMAX_XMIN_XMAX),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ibelm_xmax'
+ allocate(ibelm_ymin(NSPEC2DMAX_YMIN_YMAX),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ibelm_ymin'
+ allocate(ibelm_ymax(NSPEC2DMAX_YMIN_YMAX),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ibelm_ymax'
+ allocate(ibelm_bottom(NSPEC2D_BOTTOM),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ibelm_bottom'
allocate(ibelm_top(NSPEC2D_TOP),stat=ier)
if( ier /= 0 ) call exit_MPI(myrank,'not enough memory to allocate arrays')
! MPI cut-planes parameters along xi and along eta
- allocate(iMPIcut_xi(2,nspec))
+ allocate(iMPIcut_xi(2,nspec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array iMPIcut_xi'
allocate(iMPIcut_eta(2,nspec),stat=ier)
if( ier /= 0 ) call exit_MPI(myrank,'not enough memory to allocate arrays')
! allocate memory for arrays
- allocate(iglob(npointot))
- allocate(locval(npointot))
- allocate(ifseg(npointot))
- allocate(xp(npointot))
- allocate(yp(npointot))
+ allocate(iglob(npointot),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array iglob'
+ allocate(locval(npointot),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array locval'
+ allocate(ifseg(npointot),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ifseg'
+ allocate(xp(npointot),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array xp'
+ allocate(yp(npointot),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array yp'
allocate(zp(npointot),stat=ier)
if( ier /= 0 ) call exit_MPI(myrank,'not enough memory to allocate arrays')
@@ -370,7 +383,8 @@
call get_global(nspec,xp,yp,zp,iglob,locval,ifseg,nglob,npointot,UTM_X_MIN,UTM_X_MAX)
! put in classical format
- allocate(nodes_coords(nglob,3))
+ allocate(nodes_coords(nglob,3),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array nodes_coords'
do ispec=1,nspec
ieoff = NGLLCUBE*(ispec-1)
Modified: seismo/3D/SPECFEM3D/trunk/src/meshfem3D/get_global.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/meshfem3D/get_global.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/meshfem3D/get_global.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -44,7 +44,7 @@
double precision xp(npointot),yp(npointot),zp(npointot)
double precision UTM_X_MIN,UTM_X_MAX
- integer ispec,i,j
+ integer ispec,i,j,ier
integer ieoff,ilocnum,nseg,ioff,iseg,ig
integer, dimension(:), allocatable :: ind,ninseg,iwork
@@ -58,10 +58,14 @@
SMALLVALTOL = 1.d-10 * dabs(UTM_X_MAX - UTM_X_MIN)
! dynamically allocate arrays
- allocate(ind(npointot))
- allocate(ninseg(npointot))
- allocate(iwork(npointot))
- allocate(work(npointot))
+ allocate(ind(npointot),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ind'
+ allocate(ninseg(npointot),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ninseg'
+ allocate(iwork(npointot),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array iwork'
+ allocate(work(npointot),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array work'
! establish initial pointers
do ispec=1,nspec
Modified: seismo/3D/SPECFEM3D/trunk/src/meshfem3D/meshfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/meshfem3D/meshfem3D.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/meshfem3D/meshfem3D.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -372,7 +372,8 @@
! define number of layers
number_of_layers = number_of_interfaces! - 1
- allocate(ner_layer(number_of_layers))
+ allocate(ner_layer(number_of_layers),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ner_layer'
! loop on all the layers
do ilayer = 1,number_of_layers
@@ -404,16 +405,22 @@
endif
! dynamic allocation of mesh arrays
- allocate(rns(0:2*NER))
+ allocate(rns(0:2*NER),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array rns'
- allocate(xgrid(0:2*NER,0:2*NEX_PER_PROC_XI,0:2*NEX_PER_PROC_ETA))
- allocate(ygrid(0:2*NER,0:2*NEX_PER_PROC_XI,0:2*NEX_PER_PROC_ETA))
+ allocate(xgrid(0:2*NER,0:2*NEX_PER_PROC_XI,0:2*NEX_PER_PROC_ETA),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array xgrid'
+ allocate(ygrid(0:2*NER,0:2*NEX_PER_PROC_XI,0:2*NEX_PER_PROC_ETA),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ygrid'
allocate(zgrid(0:2*NER,0:2*NEX_PER_PROC_XI,0:2*NEX_PER_PROC_ETA),stat=ier)
if( ier /= 0 ) call exit_MPI(myrank,'not enough memory to allocate arrays')
- allocate(addressing(0:NPROC_XI-1,0:NPROC_ETA-1))
- allocate(iproc_xi_slice(0:NPROC-1))
- allocate(iproc_eta_slice(0:NPROC-1))
+ allocate(addressing(0:NPROC_XI-1,0:NPROC_ETA-1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array addressing'
+ allocate(iproc_xi_slice(0:NPROC-1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array iproc_xi_slice'
+ allocate(iproc_eta_slice(0:NPROC-1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array iproc_eta_slice'
! clear arrays
xgrid(:,:,:) = 0.d0
@@ -552,8 +559,10 @@
open(unit=IIN,file=MF_IN_DATA_FILES_PATH(1:len_trim(MF_IN_DATA_FILES_PATH)) &
//INTERFACES_FILE,status='old')
- allocate(interface_bottom(max_npx_interface,max_npy_interface))
- allocate(interface_top(max_npx_interface,max_npy_interface))
+ allocate(interface_bottom(max_npx_interface,max_npy_interface),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array interface_bottom'
+ allocate(interface_top(max_npx_interface,max_npy_interface),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array interface_top'
! read number of interfaces
call read_value_integer(IIN,DONT_IGNORE_JUNK,number_of_interfaces,'NINTERFACES')
@@ -721,12 +730,14 @@
call sync_all()
! use dynamic allocation to allocate memory for arrays
- allocate(ibool(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(xstore(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(ystore(NGLLX,NGLLY,NGLLZ,nspec))
+ allocate(ibool(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ibool'
+ allocate(xstore(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array xstore'
+ allocate(ystore(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ystore'
allocate(zstore(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
-
-! exit if there is not enough memory to allocate all the arrays
+ ! exit if there is not enough memory to allocate all the arrays
if(ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays')
call create_regions_mesh(xgrid,ygrid,zgrid,ibool, &
Modified: seismo/3D/SPECFEM3D/trunk/src/meshfem3D/save_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/meshfem3D/save_databases.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/meshfem3D/save_databases.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -71,7 +71,7 @@
! first dimension : material_id
! second dimension : #rho #vp #vs #Q_flag #anisotropy_flag #domain_id
double precision , dimension(NMATERIALS,6) :: material_properties
-
+ double precision , dimension(6) :: matpropl
integer i,ispec,iglob
! name of the database files
@@ -94,7 +94,9 @@
! Materials properties
write(15,*) NMATERIALS, 0
do idoubl = 1,NMATERIALS
- write(15,*) material_properties(idoubl,:)
+ !write(15,*) material_properties(idoubl,:)
+ matpropl(:) = material_properties(idoubl,:)
+ write(15,*) matpropl
end do
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/assemble_MPI_scalar.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/assemble_MPI_scalar.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/assemble_MPI_scalar.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -64,17 +64,21 @@
integer, dimension(:), allocatable :: request_recv_scalar_ext_mesh
- integer ipoin,iinterface
+ integer ipoin,iinterface,ier
! here we have to assemble all the contributions between partitions using MPI
! assemble only if more than one partition
if(NPROC > 1) then
- allocate(buffer_send_scalar_ext_mesh(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh))
- allocate(buffer_recv_scalar_ext_mesh(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh))
- allocate(request_send_scalar_ext_mesh(num_interfaces_ext_mesh))
- allocate(request_recv_scalar_ext_mesh(num_interfaces_ext_mesh))
+ allocate(buffer_send_scalar_ext_mesh(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array buffer_send_scalar_ext_mesh'
+ allocate(buffer_recv_scalar_ext_mesh(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array buffer_recv_scalar_ext_mesh'
+ allocate(request_send_scalar_ext_mesh(num_interfaces_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array request_send_scalar_ext_mesh'
+ allocate(request_recv_scalar_ext_mesh(num_interfaces_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array request_recv_scalar_ext_mesh'
! partition border copy into the buffer
do iinterface = 1, num_interfaces_ext_mesh
@@ -156,17 +160,21 @@
integer, dimension(:), allocatable :: request_send_scalar_ext_mesh
integer, dimension(:), allocatable :: request_recv_scalar_ext_mesh
- integer :: ipoin,iinterface
+ integer :: ipoin,iinterface,ier
! here we have to assemble all the contributions between partitions using MPI
! assemble only if more than one partition
if(NPROC > 1) then
- allocate(buffer_send_scalar_ext_mesh(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh))
- allocate(buffer_recv_scalar_ext_mesh(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh))
- allocate(request_send_scalar_ext_mesh(num_interfaces_ext_mesh))
- allocate(request_recv_scalar_ext_mesh(num_interfaces_ext_mesh))
+ allocate(buffer_send_scalar_ext_mesh(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array buffer_send_scalar_ext_mesh'
+ allocate(buffer_recv_scalar_ext_mesh(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array buffer_recv_scalar_ext_mesh'
+ allocate(request_send_scalar_ext_mesh(num_interfaces_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array request_send_scalar_ext_mesh'
+ allocate(request_recv_scalar_ext_mesh(num_interfaces_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array request_recv_scalar_ext_mesh'
! partition border copy into the buffer
do iinterface = 1, num_interfaces_ext_mesh
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/combine_surf_data.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/combine_surf_data.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/combine_surf_data.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -59,7 +59,7 @@
integer,dimension(:),allocatable :: num_ibool
logical :: HIGH_RESOLUTION_MESH, FILE_ARRAY_IS_3D
- integer :: ires, nspec_surf, npoint1, npoint2, ispec_surf, inx, iny, idim
+ integer :: ires, nspec_surf, npoint1, npoint2, ispec_surf, inx, iny, idim, ier
integer,dimension(:), allocatable :: ibelm_surf
@@ -163,10 +163,11 @@
nglob = NGLOB_AB
! allocates arrays
- allocate(ibool(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
- allocate(mask_ibool(NGLOB_AB))
- allocate(num_ibool(NGLOB_AB))
- allocate(xstore(NGLOB_AB),ystore(NGLOB_AB),zstore(NGLOB_AB))
+ allocate(ibool(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
+ mask_ibool(NGLOB_AB), &
+ num_ibool(NGLOB_AB), &
+ xstore(NGLOB_AB),ystore(NGLOB_AB),zstore(NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ibool etc.'
! surface file
local_ibool_surf_file = trim(prname) // 'ibelm_' //trim(surfname)// '.bin'
@@ -179,16 +180,21 @@
read(28) npoint1
read(28) npoint2
- if (it == 1) allocate(ibelm_surf(nspec_surf))
+ if (it == 1) then
+ allocate(ibelm_surf(nspec_surf),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ibelm_surf'
+ endif
read(28) ibelm_surf
close(28)
print *, trim(local_ibool_surf_file)
if (it == 1) then
if (FILE_ARRAY_IS_3D) then
- allocate(data_3D(NGLLX,NGLLY,NGLLZ,NSPEC_AB),dat3D(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(data_3D(NGLLX,NGLLY,NGLLZ,NSPEC_AB),dat3D(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array data_3D'
else
- allocate(data_2D(NGLLX,NGLLY,nspec_surf),dat2D(NGLLX,NGLLY,nspec_surf))
+ allocate(data_2D(NGLLX,NGLLY,nspec_surf),dat2D(NGLLX,NGLLY,nspec_surf),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array data_2D'
endif
endif
@@ -316,9 +322,10 @@
nglob = NGLOB_AB
! allocates arrays
- allocate(ibool(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
- allocate(mask_ibool(NGLOB_AB))
- allocate(num_ibool(NGLOB_AB))
+ allocate(ibool(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
+ mask_ibool(NGLOB_AB), &
+ num_ibool(NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ibool etc.'
np = npoint * (it-1)
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/combine_vol_data.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/combine_vol_data.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/combine_vol_data.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -51,7 +51,7 @@
integer :: NSPEC_AB, NGLOB_AB
integer :: numpoin
- integer :: i, ios, it
+ integer :: i, ios, it, ier
integer :: iproc, proc1, proc2, num_node, node_list(300)
integer :: np, ne, npp, nee, nelement, njunk
@@ -185,11 +185,13 @@
read(27) NGLOB_AB
! ibool file
- allocate(ibool(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(ibool(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ibool'
read(27) ibool
! global point arrays
- allocate(xstore(NGLOB_AB),ystore(NGLOB_AB),zstore(NGLOB_AB))
+ allocate(xstore(NGLOB_AB),ystore(NGLOB_AB),zstore(NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array xstore etc.'
read(27) xstore
read(27) ystore
read(27) zstore
@@ -206,13 +208,13 @@
stop
endif
- allocate(dat(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ios)
- if( ios /= 0 ) stop 'error allocating dat array'
+ allocate(dat(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating dat array'
! uses conversion to real values
if( CUSTOM_REAL == SIZE_DOUBLE ) then
- allocate(data(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ios)
- if( ios /= 0 ) stop 'error allocating data array'
+ allocate(data(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating data array'
read(28) data
dat = sngl(data)
deallocate(data)
@@ -265,7 +267,8 @@
read(27) NGLOB_AB
! ibool file
- allocate(ibool(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(ibool(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ibool'
read(27) ibool
close(27)
@@ -325,7 +328,7 @@
integer, dimension(:,:,:,:),allocatable :: ibool
logical, dimension(:),allocatable :: mask_ibool
integer :: NSPEC_AB, NGLOB_AB
- integer :: it,iproc,npoint,nelement,ios,ispec
+ integer :: it,iproc,npoint,nelement,ios,ispec,ier
integer :: iglob1, iglob2, iglob3, iglob4, iglob5, iglob6, iglob7, iglob8
character(len=256) :: prname_lp
@@ -349,7 +352,8 @@
read(27) NGLOB_AB
! gets ibool
if( .not. HIGH_RESOLUTION_MESH ) then
- allocate(ibool(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(ibool(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ibool'
read(27) ibool
endif
close(27)
@@ -368,7 +372,8 @@
else
! mark element corners (global AVS or DX points)
- allocate(mask_ibool(NGLOB_AB))
+ allocate(mask_ibool(NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array mask_ibool'
mask_ibool = .false.
do ispec=1,NSPEC_AB
iglob1=ibool(1,1,1,ispec)
@@ -427,7 +432,7 @@
! local parameters
logical,dimension(:),allocatable :: mask_ibool
real :: x, y, z
- integer :: ispec
+ integer :: ispec,ier
integer :: iglob1, iglob2, iglob3, iglob4, iglob5, iglob6, iglob7, iglob8
! writes out total number of points
@@ -436,7 +441,8 @@
endif
! writes our corner point locations
- allocate(mask_ibool(NGLOB_AB))
+ allocate(mask_ibool(NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array mask_ibool'
mask_ibool(:) = .false.
numpoin = 0
do ispec=1,NSPEC_AB
@@ -562,7 +568,7 @@
! local parameters
logical,dimension(:),allocatable :: mask_ibool
real :: x, y, z
- integer :: ispec,i,j,k,iglob
+ integer :: ispec,i,j,k,iglob,ier
! writes out total number of points
if (it == 1) then
@@ -570,7 +576,9 @@
endif
! writes out point locations and values
- allocate(mask_ibool(NGLOB_AB))
+ allocate(mask_ibool(NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array mask_ibool'
+
mask_ibool(:) = .false.
numpoin = 0
do ispec=1,NSPEC_AB
@@ -613,7 +621,7 @@
! local parameters
logical,dimension(:),allocatable :: mask_ibool
integer,dimension(:),allocatable :: num_ibool
- integer :: ispec
+ integer :: ispec,ier
integer :: iglob1, iglob2, iglob3, iglob4, iglob5, iglob6, iglob7, iglob8
integer :: n1, n2, n3, n4, n5, n6, n7, n8
@@ -623,8 +631,10 @@
end if
! writes out element indices
- allocate(mask_ibool(NGLOB_AB))
- allocate(num_ibool(NGLOB_AB))
+ allocate(mask_ibool(NGLOB_AB), &
+ num_ibool(NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array mask_ibool'
+
mask_ibool(:) = .false.
num_ibool(:) = 0
numpoin = 0
@@ -729,7 +739,7 @@
! local parameters
logical,dimension(:),allocatable :: mask_ibool
integer,dimension(:),allocatable :: num_ibool
- integer :: ispec,i,j,k
+ integer :: ispec,i,j,k,ier
integer :: iglob,iglob1, iglob2, iglob3, iglob4, iglob5, iglob6, iglob7, iglob8
integer :: n1, n2, n3, n4, n5, n6, n7, n8
@@ -740,8 +750,10 @@
endif
! sets numbering num_ibool respecting mask
- allocate(mask_ibool(NGLOB_AB))
- allocate(num_ibool(NGLOB_AB))
+ allocate(mask_ibool(NGLOB_AB), &
+ num_ibool(NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array mask_ibool'
+
mask_ibool(:) = .false.
num_ibool(:) = 0
numpoin = 0
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/create_movie_shakemap_AVS_DX_GMT.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/create_movie_shakemap_AVS_DX_GMT.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/create_movie_shakemap_AVS_DX_GMT.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -249,27 +249,30 @@
npointot = NGNOD2D_AVS_DX * nspectot_AVS_max
! allocate arrays for sorting routine
- allocate(iglob(npointot),loc(npointot))
- allocate(ifseg(npointot))
- allocate(xp(npointot),yp(npointot),zp(npointot))
- allocate(xp_save(npointot),yp_save(npointot),zp_save(npointot))
- allocate(field_display(npointot))
- allocate(mask_point(npointot))
- allocate(ireorder(npointot))
+ allocate(iglob(npointot),loc(npointot), &
+ ifseg(npointot), &
+ xp(npointot),yp(npointot),zp(npointot), &
+ xp_save(npointot),yp_save(npointot),zp_save(npointot), &
+ field_display(npointot), &
+ mask_point(npointot), &
+ ireorder(npointot),stat=ier)
+ if( ier /= 0 ) stop 'error allocating arrays for sorting routine'
! allocates data arrays
- allocate(store_val_x(ilocnum))
- allocate(store_val_y(ilocnum))
- allocate(store_val_z(ilocnum))
- allocate(store_val_ux(ilocnum))
- allocate(store_val_uy(ilocnum))
- allocate(store_val_uz(ilocnum))
+ allocate(store_val_x(ilocnum), &
+ store_val_y(ilocnum), &
+ store_val_z(ilocnum), &
+ store_val_ux(ilocnum), &
+ store_val_uy(ilocnum), &
+ store_val_uz(ilocnum),stat=ier)
+ if( ier /= 0 ) stop 'error allocating arrays for data arrays'
if(USE_HIGHRES_FOR_MOVIES) then
- allocate(x(NGLLX,NGLLY))
- allocate(y(NGLLX,NGLLY))
- allocate(z(NGLLX,NGLLY))
- allocate(display(NGLLX,NGLLY))
+ allocate(x(NGLLX,NGLLY), &
+ y(NGLLX,NGLLY), &
+ z(NGLLX,NGLLY), &
+ display(NGLLX,NGLLY),stat=ier)
+ if( ier /= 0 ) stop 'error allocating arrays for highres'
endif
! user output
@@ -840,7 +843,7 @@
double precision xp(npointot),yp(npointot),zp(npointot)
integer nspec,nglob
- integer ispec,i,j
+ integer ispec,i,j,ier
integer ieoff,ilocnum,nseg,ioff,iseg,ig
integer, dimension(:), allocatable :: ind,ninseg,iwork
@@ -855,10 +858,11 @@
print *, 'SMALLVALTOL', SMALLVALTOL
! dynamically allocate arrays
- allocate(ind(npointot))
- allocate(ninseg(npointot))
- allocate(iwork(npointot))
- allocate(work(npointot))
+ allocate(ind(npointot), &
+ ninseg(npointot), &
+ iwork(npointot), &
+ work(npointot),stat=ier)
+ if( ier /= 0 ) stop 'error allocating arrays ind etc.'
! establish initial pointers
do ispec=1,nspec
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/create_serial_name_database.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/create_serial_name_database.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/create_serial_name_database.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -39,7 +39,7 @@
integer iprocloop,nproc_max_loop
integer, dimension(:), allocatable :: num_active_proc
-
+ integer ier
nproc_max_loop = NPROC-1
! create the name for the database of the current slide and region
@@ -49,7 +49,8 @@
if(.not. LOCAL_PATH_IS_ALSO_GLOBAL) then
! allocate array for active processors
- allocate(num_active_proc(0:nproc_max_loop))
+ allocate(num_active_proc(0:nproc_max_loop),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array num_active_proc'
! read filtered file with name of active machines
open(unit=48,file=trim(OUTPUT_FILES)//'/filtered_machines.txt',status='old',action='read')
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/get_attenuation_model.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/get_attenuation_model.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/get_attenuation_model.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -744,14 +744,17 @@
integer Qtmp
integer, save :: first_time_called = 1
double precision, parameter :: ZERO_TOL = 1.e-5
-
+ integer ier
+
if(first_time_called == 1) then
first_time_called = 0
AM_S%Q_resolution = 10**ATTENUATION_COMP_RESOLUTION
AM_S%Q_max = ATTENUATION_COMP_MAXIMUM
Qtmp = AM_S%Q_resolution * AM_S%Q_max
- allocate(AM_S%tau_eps_storage(N_SLS, Qtmp))
- allocate(AM_S%Qmu_storage(Qtmp))
+
+ allocate(AM_S%tau_eps_storage(N_SLS, Qtmp), &
+ AM_S%Qmu_storage(Qtmp),stat=ier)
+ if( ier /= 0 ) stop 'error allocating arrays for attenuation storage'
AM_S%Qmu_storage(:) = -1
endif
@@ -938,10 +941,12 @@
double precision Q_in
double precision, dimension(nf_in) :: f_in
double precision, dimension(nsls_in) :: tau_s_in
+ integer ier
+
+ allocate(AS_V%f(nf_in), &
+ AS_V%tau_s(nsls_in),stat=ier)
+ if( ier /= 0 ) stop 'error allocating arrays for attenuation simplex'
- allocate(AS_V%f(nf_in))
- allocate(AS_V%tau_s(nsls_in))
-
AS_V%nf = nf_in
AS_V%nsls = nsls_in
AS_V%f = f_in
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/get_element_face.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/get_element_face.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/get_element_face.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -122,8 +122,8 @@
( (xcoord(1)-xcoord(2))**2 &
+ (ycoord(1)-ycoord(2))**2 &
+ (zcoord(1)-zcoord(2))**2 ) ) then
- print*,'error element face midpoint distance:',midpoint_distances(iloc(1)),&
- ( (xcoord(1)-xcoord(2))**2+(ycoord(1)-ycoord(2))**2+(zcoord(1)-zcoord(2))**2 )
+ print*,'error element face midpoint distance:',midpoint_distances(iloc(1)), &
+ (xcoord(1)-xcoord(2))**2+(ycoord(1)-ycoord(2))**2+(zcoord(1)-zcoord(2))**2
! corner locations
do icorner=1,NGNOD2D
i = iface_all_corner_ijk(1,icorner,iloc(1))
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/read_value_parameters.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/read_value_parameters.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/read_value_parameters.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -36,7 +36,7 @@
integer ierr
common /param_err_common/ ierr
- call param_read(string_read, len(string_read), name, len(name), ierr);
+ call param_read(string_read, len(string_read), name, len(name), ierr)
if (ierr .ne. 0) return
read(string_read,*) value_to_read
@@ -54,7 +54,7 @@
integer ierr
common /param_err_common/ ierr
- call param_read(string_read, len(string_read), name, len(name), ierr);
+ call param_read(string_read, len(string_read), name, len(name), ierr)
if (ierr .ne. 0) return
read(string_read,*) value_to_read
@@ -72,7 +72,7 @@
integer ierr
common /param_err_common/ ierr
- call param_read(string_read, len(string_read), name, len(name), ierr);
+ call param_read(string_read, len(string_read), name, len(name), ierr)
if (ierr .ne. 0) return
read(string_read,*) value_to_read
@@ -90,7 +90,7 @@
integer ierr
common /param_err_common/ ierr
- call param_read(string_read, len(string_read), name, len(name), ierr);
+ call param_read(string_read, len(string_read), name, len(name), ierr)
if (ierr .ne. 0) return
value_to_read = string_read
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/smooth_vol_data.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/smooth_vol_data.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/smooth_vol_data.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -256,18 +256,21 @@
read(27) NGLOB_AB
! ibool file
- allocate(ibool(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(ibool(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ibool'
read(27) ibool
! global point arrays
- allocate(xstore(NGLOB_AB),ystore(NGLOB_AB),zstore(NGLOB_AB))
+ allocate(xstore(NGLOB_AB),ystore(NGLOB_AB),zstore(NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array xstore etc.'
read(27) xstore
read(27) ystore
read(27) zstore
! reads in jacobian
allocate(dummy(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
- jacobian(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ jacobian(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array dummy and jacobian'
read(27) dummy ! xix
read(27) dummy ! xiy
read(27) dummy ! xiz
@@ -283,7 +286,8 @@
read(27) dummy ! kappastore
read(27) dummy ! mustore
- allocate(ldummy(NSPEC_AB))
+ allocate(ldummy(NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ldummy'
read(27) ldummy ! ispec_is_acoustic
call any_all_l( ANY(ldummy), ACOUSTIC_SIMULATION )
@@ -295,7 +299,8 @@
deallocate(ldummy)
- allocate(dummy_1(NGLOB_AB))
+ allocate(dummy_1(NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array dummy_1'
if( ACOUSTIC_SIMULATION ) then
read(27) dummy_1 ! rmass_acoustic
read(27) dummy ! rhostore
@@ -312,10 +317,11 @@
deallocate(dummy)
read(27) idummy_a ! num_abs_boundary_faces
- allocate(idummy(idummy_a))
- allocate(idummy_3(3,NGLLSQUARE,idummy_a))
- allocate(dummy_2(NGLLSQUARE,idummy_a))
- allocate(dummy_3(NDIM,NGLLSQUARE,idummy_a))
+ allocate(idummy(idummy_a), &
+ idummy_3(3,NGLLSQUARE,idummy_a), &
+ dummy_2(NGLLSQUARE,idummy_a), &
+ dummy_3(NDIM,NGLLSQUARE,idummy_a),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array idummy etc.'
read(27) idummy ! abs_boundary_ispec
read(27) idummy_3 ! abs_boundary_ijk
read(27) dummy_2 ! abs_boundary_jacobian2Dw
@@ -323,10 +329,11 @@
deallocate( idummy,idummy_3,dummy_2,dummy_3)
read(27) idummy_a ! num_free_surface_faces
- allocate(idummy(idummy_a))
- allocate(idummy_3(3,NGLLSQUARE,idummy_a))
- allocate(dummy_2(NGLLSQUARE,idummy_a))
- allocate(dummy_3(NDIM,NGLLSQUARE,idummy_a))
+ allocate(idummy(idummy_a), &
+ idummy_3(3,NGLLSQUARE,idummy_a), &
+ dummy_2(NGLLSQUARE,idummy_a), &
+ dummy_3(NDIM,NGLLSQUARE,idummy_a),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array idummy etc.'
read(27) idummy ! free_surface_ispec
read(27) idummy_3 ! free_surface_ijk
read(27) dummy_2 ! free_surface_jacobian2Dw
@@ -334,10 +341,11 @@
deallocate( idummy,idummy_3,dummy_2,dummy_3)
read(27) idummy_a ! num_coupling_ac_el_faces
- allocate(idummy(idummy_a))
- allocate(idummy_3(3,NGLLSQUARE,idummy_a))
- allocate(dummy_2(NGLLSQUARE,idummy_a))
- allocate(dummy_3(NDIM,NGLLSQUARE,idummy_a))
+ allocate(idummy(idummy_a), &
+ idummy_3(3,NGLLSQUARE,idummy_a), &
+ dummy_2(NGLLSQUARE,idummy_a), &
+ dummy_3(NDIM,NGLLSQUARE,idummy_a),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array idummy etc.'
read(27) idummy ! coupling_ac_el_ispec
read(27) idummy_3 ! coupling_ac_el_ijk
read(27) dummy_2 ! coupling_ac_el_jacobian2Dw
@@ -346,18 +354,21 @@
read(27) num_interfaces_ext_mesh ! num_interfaces_ext_mesh
read(27) idummy_a ! max_nibool_interfaces_ext_mesh
- allocate(my_neighbours_ext_mesh(num_interfaces_ext_mesh))
+ allocate(my_neighbours_ext_mesh(num_interfaces_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array my_neighbours_ext_mesh'
read(27) my_neighbours_ext_mesh
close(27)
! get the location of the center of the elements and local points
- allocate(xl(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
- allocate(yl(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
- allocate(zl(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
- allocate(cx0(NSPEC_AB))
- allocate(cy0(NSPEC_AB))
- allocate(cz0(NSPEC_AB))
+ allocate(xl(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
+ yl(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
+ zl(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
+ cx0(NSPEC_AB), &
+ cy0(NSPEC_AB), &
+ cz0(NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array xl etc.'
+
do ispec = 1, nspec_AB
do k = 1, NGLLZ
do j = 1, NGLLY
@@ -408,8 +419,10 @@
! loops over slices
! each process reads in his own neighbor slices and gaussian filters the values
- allocate(tk(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
- allocate(bk(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(tk(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
+ bk(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array tk and bk'
+
tk = 0.0_CUSTOM_REAL
bk = 0.0_CUSTOM_REAL
do it=1,num_interfaces_ext_mesh+1
@@ -430,18 +443,21 @@
read(27) NGLOB_N
! ibool file
- allocate(ibool(NGLLX,NGLLY,NGLLZ,NSPEC_N))
+ allocate(ibool(NGLLX,NGLLY,NGLLZ,NSPEC_N),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ibool'
read(27) ibool
! global point arrays
- allocate(xstore(NGLOB_N),ystore(NGLOB_N),zstore(NGLOB_N))
+ allocate(xstore(NGLOB_N),ystore(NGLOB_N),zstore(NGLOB_N),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array xstore etc.'
read(27) xstore
read(27) ystore
read(27) zstore
! reads in jacobian
allocate(dummy(NGLLX,NGLLY,NGLLZ,NSPEC_N), &
- jacobian(NGLLX,NGLLY,NGLLZ,NSPEC_N))
+ jacobian(NGLLX,NGLLY,NGLLZ,NSPEC_N),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array dummy and jacobian'
read(27) dummy ! xix
read(27) dummy ! xiy
read(27) dummy ! xiz
@@ -456,12 +472,14 @@
deallocate(dummy)
! get the location of the center of the elements and local points
- allocate(xx(NGLLX,NGLLY,NGLLZ,NSPEC_N))
- allocate(yy(NGLLX,NGLLY,NGLLZ,NSPEC_N))
- allocate(zz(NGLLX,NGLLY,NGLLZ,NSPEC_N))
- allocate(cx(NSPEC_N))
- allocate(cy(NSPEC_N))
- allocate(cz(NSPEC_N))
+ allocate(xx(NGLLX,NGLLY,NGLLZ,NSPEC_N), &
+ yy(NGLLX,NGLLY,NGLLZ,NSPEC_N), &
+ zz(NGLLX,NGLLY,NGLLZ,NSPEC_N), &
+ cx(NSPEC_N), &
+ cy(NSPEC_N), &
+ cz(NSPEC_N),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array xx etc.'
+
do ispec = 1, nspec_N
do k = 1, NGLLZ
do j = 1, NGLLY
@@ -492,8 +510,8 @@
print *,'Error opening ',trim(local_data_file)
stop
endif
- allocate(dat(NGLLX,NGLLY,NGLLZ,NSPEC_N),stat=ios)
- if( ios /= 0 ) stop 'error allocating dat array'
+ allocate(dat(NGLLX,NGLLY,NGLLZ,NSPEC_N),stat=ier)
+ if( ier /= 0 ) stop 'error allocating dat array'
read(28) dat
close(28)
@@ -559,7 +577,9 @@
! normalizes
!if(myrank==0) print*, 'normalizes values ...'
- allocate(dat_smooth(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(dat_smooth(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array dat_smooth'
+
dat_smooth = 0.0_CUSTOM_REAL
do ispec = 1, nspec_AB
do k = 1, NGLLZ
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/PML_init.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/PML_init.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/PML_init.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -142,7 +142,7 @@
! local parameters
real(kind=CUSTOM_REAL):: d,dprime,d_glob,dprime_glob
real(kind=CUSTOM_REAL) :: dominant_wavelength,hdur_max
- integer :: count,ilayer,sign
+ integer :: count,ilayer,sign,ier
! sets flag
PML = .true.
@@ -155,19 +155,23 @@
endif
! PML element type array: 1 = face, 2 = edge, 3 = corner
- allocate(ispec_is_PML_inum(NSPEC_AB))
+ allocate(ispec_is_PML_inum(NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ispec_is_PML_inum'
num_PML_ispec = 0
! PML interface points between PML and "regular" region
- allocate(iglob_is_PML_interface(NGLOB_AB))
+ allocate(iglob_is_PML_interface(NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array iglob_is_PML_interface'
iglob_is_PML_interface(:) = 0
! PML global points
- allocate(iglob_is_PML(NGLOB_AB))
+ allocate(iglob_is_PML(NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array iglob_is_PML'
iglob_is_PML(:) = 0
! PML ibool mask
- allocate(PML_mask_ibool(NGLOB_AB))
+ allocate(PML_mask_ibool(NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array PML_mask_ibool'
PML_mask_ibool(:) = .false.
! determines dominant wavelength based on maximum model speed
@@ -229,15 +233,18 @@
chi2(NGLLX,NGLLY,NGLLZ,num_PML_ispec),&
chi2_t(NGLLX,NGLLY,NGLLZ,num_PML_ispec),&
chi3(NGLLX,NGLLY,NGLLZ,num_PML_ispec),&
- chi4(NGLLX,NGLLY,NGLLZ,num_PML_ispec))
+ chi4(NGLLX,NGLLY,NGLLZ,num_PML_ispec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array chi1 etc.'
allocate(chi1_dot(NGLLX,NGLLY,NGLLZ,num_PML_ispec),&
chi2_t_dot(NGLLX,NGLLY,NGLLZ,num_PML_ispec),&
chi3_dot(NGLLX,NGLLY,NGLLZ,num_PML_ispec),&
- chi4_dot(NGLLX,NGLLY,NGLLZ,num_PML_ispec))
+ chi4_dot(NGLLX,NGLLY,NGLLZ,num_PML_ispec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array chi1_dot etc.'
allocate(chi1_dot_dot(NGLLX,NGLLY,NGLLZ,num_PML_ispec),&
chi2_t_dot_dot(NGLLX,NGLLY,NGLLZ,num_PML_ispec),&
chi3_dot_dot(NGLLX,NGLLY,NGLLZ,num_PML_ispec),&
- chi4_dot_dot(NGLLX,NGLLY,NGLLZ,num_PML_ispec))
+ chi4_dot_dot(NGLLX,NGLLY,NGLLZ,num_PML_ispec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array chi1_dot_dot etc.'
! potentials
chi1 = 0._CUSTOM_REAL
@@ -310,14 +317,15 @@
! local parameters
real(kind=CUSTOM_REAL),dimension(:,:),allocatable:: temp_ispec_pml_normal
integer,dimension(:),allocatable:: temp_is_pml_elem
- integer:: iface,count,new_elemts,ispec,icorner,igll,iglobf
+ integer:: iface,count,new_elemts,ispec,icorner,igll,iglobf,ier
integer:: i,j,k,iglobcount,iglobcorners(NGNOD)
integer,dimension(3,NGNOD),parameter :: ielem_corner_ijk = &
reshape((/ 1,1,1, 1,NGLLY,1, 1,NGLLY,NGLLZ, 1,1,NGLLZ, &
NGLLX,1,1, NGLLX,NGLLY,1, NGLLX,NGLLY,NGLLZ, NGLLX,1,NGLLZ /),(/3,NGNOD/))
! temporary arrays
- allocate(temp_is_pml_elem(NSPEC_AB))
- allocate(temp_ispec_pml_normal(NDIM,NSPEC_AB))
+ allocate(temp_is_pml_elem(NSPEC_AB), &
+ temp_ispec_pml_normal(NDIM,NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array temp_is_pml_elem'
temp_is_pml_elem(:) = 0
temp_ispec_pml_normal(:,:) = 0._CUSTOM_REAL
@@ -417,7 +425,7 @@
! local parameters
real(kind=CUSTOM_REAL) :: length
- integer :: ispec,ispecPML
+ integer :: ispec,ispecPML,ier
! sets new element type flags
ispec_is_PML_inum(:) = temp_is_pml_elem(:)
@@ -428,8 +436,9 @@
! re-allocates arrays
if( allocated(PML_normal) ) deallocate(PML_normal)
if( allocated(PML_ispec) ) deallocate(PML_ispec)
- allocate(PML_ispec(num_PML_ispec))
- allocate(PML_normal(NDIM,num_PML_ispec))
+ allocate(PML_ispec(num_PML_ispec), &
+ PML_normal(NDIM,num_PML_ispec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array PML_ispec'
! stores PML elements flags and normals
ispecPML = 0
@@ -481,13 +490,14 @@
! local parameters
integer,dimension(:),allocatable:: temp_regulardomain
- integer:: i,j,k,ispec,iglob
+ integer:: i,j,k,ispec,iglob,ier
! PML interface points array
iglob_is_PML_interface(:) = 0
! temporary arrays
- allocate(temp_regulardomain(NGLOB_AB))
+ allocate(temp_regulardomain(NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array temp_regulardomain'
temp_regulardomain(:) = 0
! global PML points
@@ -646,11 +656,12 @@
real(kind=CUSTOM_REAL) :: d
real(kind=CUSTOM_REAL) :: width
- integer:: i,j,k,ispec,iglob,ispecPML,iglobf
+ integer:: i,j,k,ispec,iglob,ispecPML,iglobf,ier
integer:: ispecB,igll,iface
! stores damping coefficient
- allocate( PML_damping_d(NGLLX,NGLLY,NGLLZ,num_PML_ispec))
+ allocate( PML_damping_d(NGLLX,NGLLY,NGLLZ,num_PML_ispec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array PML_damping_d'
PML_damping_d(:,:,:,:) = 0._CUSTOM_REAL
! loops over all PML elements
@@ -765,10 +776,11 @@
real(kind=CUSTOM_REAL) :: xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl
real(kind=CUSTOM_REAL) :: nx,ny,nz
real(kind=CUSTOM_REAL) :: d_dx,d_dy,d_dz,tempd_dx,tempd_dy,tempd_dz
- integer :: ispec,i,j,k,l,ispecPML
+ integer :: ispec,i,j,k,l,ispecPML,ier
! dprime derivatives
- allocate( PML_damping_dprime(NGLLX,NGLLY,NGLLZ,num_PML_ispec))
+ allocate( PML_damping_dprime(NGLLX,NGLLY,NGLLZ,num_PML_ispec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array PML_damping_dprime'
PML_damping_dprime(:,:,:,:) = 0._CUSTOM_REAL
! loops over all PML elements
@@ -847,7 +859,7 @@
real(kind=CUSTOM_REAL),dimension(:,:),allocatable:: iglob_pml_normal
real(kind=CUSTOM_REAL),dimension(:,:),allocatable:: ispec_pml_normal
integer,dimension(:),allocatable:: is_pml_elem
- integer:: i,j,k,iglob,count,ispecPML,ispec,new_elemts
+ integer:: i,j,k,iglob,count,ispecPML,ispec,new_elemts,ier
integer :: iface,icorner,ipmlcorners
integer,dimension(3,4),parameter :: iface1_corner_ijk = &
@@ -872,9 +884,10 @@
logical :: is_done
! temporary arrays
- allocate(is_pml_elem(NSPEC_AB))
- allocate(iglob_pml_normal(NDIM,NGLOB_AB))
- allocate(ispec_pml_normal(NDIM,NSPEC_AB))
+ allocate(is_pml_elem(NSPEC_AB), &
+ iglob_pml_normal(NDIM,NGLOB_AB), &
+ ispec_pml_normal(NDIM,NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array is_pml_elem etc.'
iglob_pml_normal(:,:) = 0._CUSTOM_REAL
ispec_pml_normal(:,:) = 0._CUSTOM_REAL
@@ -1022,7 +1035,7 @@
reshape((/ iface1_corner_ijk,iface2_corner_ijk, &
iface3_corner_ijk,iface4_corner_ijk, &
iface5_corner_ijk,iface6_corner_ijk /),(/3,4,6/)) ! all faces
- integer:: ispecngb,iadj,ipmlinterface
+ integer:: ispecngb,iadj,ipmlinterface,ier
integer :: ispecPMLngb_corner,ispecPMLngb_edge,ispecPMLngb_sngl
integer,dimension(:),allocatable :: iglob_nadj,ispec_is_PML_inum_org
integer,dimension(:,:,:),allocatable :: iglob_adj
@@ -1031,7 +1044,8 @@
! checks normals for elements adjacent to edge/corner elements
! assigns element information to each global point
! (note: mpi partitioning/interface between elements not considered yet)
- allocate(iglob_nadj(NGLOB_AB),iglob_adj(2,32,NGLOB_AB))
+ allocate(iglob_nadj(NGLOB_AB),iglob_adj(2,32,NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array iglob_nadj'
iglob_nadj(:) = 0
iglob_adj(:,:,:) = 0
do ispecPML=1,num_PML_ispec
@@ -1059,7 +1073,8 @@
! finds neighbors based on common nodes and changes type and normal accordingly
! for edges and corners
- allocate(ispec_is_PML_inum_org(NSPEC_AB))
+ allocate(ispec_is_PML_inum_org(NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ispec_is_PML_inum_org'
ispec_is_PML_inum_org(:) = ispec_is_PML_inum(:)
do ispecPML=1,num_PML_ispec
ispec = PML_ispec(ispecPML)
@@ -1169,7 +1184,7 @@
real(kind=CUSTOM_REAL),dimension(:,:),allocatable:: ispec_normal
real(kind=CUSTOM_REAL),dimension(:,:,:,:),allocatable:: temp_gllvalues
integer,dimension(:),allocatable :: temp_iglob
- integer :: count,iglob,ispecPML,ispec
+ integer :: count,iglob,ispecPML,ispec,ier
character(len=256) :: vtkfilename
! element type flags
@@ -1188,7 +1203,8 @@
count = count+1
endif
enddo
- allocate(temp_iglob(count))
+ allocate(temp_iglob(count),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array temp_iglob'
count = 0
do iglob=1,NGLOB_AB
if( iglob_is_PML_interface(iglob) > 0 ) then
@@ -1204,7 +1220,8 @@
! pml normals
if( .false. ) then
- allocate(ispec_normal(3,NSPEC_AB) )
+ allocate(ispec_normal(3,NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ispec_normal'
ispec_normal(:,:) = 0._CUSTOM_REAL
do ispecPML=1,num_PML_ispec
ispec = PML_ispec(ispecPML)
@@ -1218,7 +1235,8 @@
! pml damping coefficients
if( .false. ) then
- allocate(temp_gllvalues(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(temp_gllvalues(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array temp_gllvalues'
temp_gllvalues = 0._CUSTOM_REAL
do ispecPML=1,num_PML_ispec
ispec = PML_ispec(ispecPML)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/assemble_MPI_vector.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/assemble_MPI_vector.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/assemble_MPI_vector.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -60,7 +60,7 @@
integer, dimension(:), allocatable :: request_send_vector_ext_mesh
integer, dimension(:), allocatable :: request_recv_vector_ext_mesh
- integer ipoin,iinterface
+ integer ipoin,iinterface,ier
! here we have to assemble all the contributions between partitions using MPI
@@ -68,10 +68,14 @@
! assemble only if more than one partition
if(NPROC > 1) then
- allocate(buffer_send_vector_ext_mesh(NDIM,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh))
- allocate(buffer_recv_vector_ext_mesh(NDIM,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh))
- allocate(request_send_vector_ext_mesh(num_interfaces_ext_mesh))
- allocate(request_recv_vector_ext_mesh(num_interfaces_ext_mesh))
+ allocate(buffer_send_vector_ext_mesh(NDIM,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array buffer_send_vector_ext_mesh'
+ allocate(buffer_recv_vector_ext_mesh(NDIM,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array buffer_recv_vector_ext_mesh'
+ allocate(request_send_vector_ext_mesh(num_interfaces_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array request_send_vector_ext_mesh'
+ allocate(request_recv_vector_ext_mesh(num_interfaces_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array request_recv_vector_ext_mesh'
! partition border copy into the buffer
do iinterface = 1, num_interfaces_ext_mesh
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_acoustic.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -215,8 +215,10 @@
if (SIMULATION_TYPE == 2 .or. SIMULATION_TYPE == 3) then
! read in adjoint sources block by block (for memory consideration)
- ! e.g., in exploration experiments, both the number of receivers (nrec) and the number of time steps (NSTEP) are huge,
- ! which may cause problems since we have a large array: adj_sourcearrays(nadj_rec_local,NSTEP,NDIM,NGLLX,NGLLY,NGLLZ)
+ ! e.g., in exploration experiments, both the number of receivers (nrec) and
+ ! the number of time steps (NSTEP) are huge,
+ ! which may cause problems since we have a large array:
+ ! adj_sourcearrays(nadj_rec_local,NSTEP,NDIM,NGLLX,NGLLY,NGLLZ)
! figure out if we need to read in a chunk of the adjoint source at this timestep
it_sub_adj = ceiling( dble(it)/dble(NTSTEP_BETWEEN_READ_ADJSRC) ) !chunk_number
@@ -224,7 +226,8 @@
! needs to read in a new chunk/block of the adjoint source
! note that for each partition, we divide it into two parts --- boundaries and interior --- indicated by 'phase_is_inner'
- ! we first do calculations for the boudaries, and then start communication with other partitions while we calculate for the inner part
+ ! we first do calculations for the boudaries, and then start communication
+ ! with other partitions while we calculate for the inner part
! this must be done carefully, otherwise the adjoint sources may be added twice
if (ibool_read_adj_arrays .and. (.not. phase_is_inner)) then
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_elastic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_elastic.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_add_sources_elastic.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -198,8 +198,10 @@
if (SIMULATION_TYPE == 2 .or. SIMULATION_TYPE == 3) then
! read in adjoint sources block by block (for memory consideration)
- ! e.g., in exploration experiments, both the number of receivers (nrec) and the number of time steps (NSTEP) are huge,
- ! which may cause problems since we have a large array: adj_sourcearrays(nadj_rec_local,NSTEP,NDIM,NGLLX,NGLLY,NGLLZ)
+ ! e.g., in exploration experiments, both the number of receivers (nrec) and
+ ! the number of time steps (NSTEP) are huge,
+ ! which may cause problems since we have a large array:
+ ! adj_sourcearrays(nadj_rec_local,NSTEP,NDIM,NGLLX,NGLLY,NGLLZ)
! figure out if we need to read in a chunk of the adjoint source at this timestep
it_sub_adj = ceiling( dble(it)/dble(NTSTEP_BETWEEN_READ_ADJSRC) ) !chunk_number
@@ -207,7 +209,8 @@
! needs to read in a new chunk/block of the adjoint source
! note that for each partition, we divide it into two parts --- boundaries and interior --- indicated by 'phase_is_inner'
- ! we first do calculations for the boudaries, and then start communication with other partitions while calculate for the inner part
+ ! we first do calculations for the boudaries, and then start communication
+ ! with other partitions while calculate for the inner part
! this must be done carefully, otherwise the adjoint sources may be added twice
if (ibool_read_adj_arrays .and. (.not. phase_is_inner)) then
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_PML.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_PML.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/compute_forces_acoustic_PML.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -672,11 +672,13 @@
!local parameters
real(kind=CUSTOM_REAL),dimension(:),allocatable:: contributions,contributions_dot
real(kind=CUSTOM_REAL):: d
- integer :: ispec,ispecPML,i,j,k,iglob
+ integer :: ispec,ispecPML,i,j,k,iglob,ier
! updates local points in PML
- allocate(contributions_dot(NGLOB_AB))
- allocate(contributions(NGLOB_AB))
+ allocate(contributions_dot(NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array contributions_dot'
+ allocate(contributions(NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array contributions'
contributions_dot(:) = 0._CUSTOM_REAL
contributions(:) = 0._CUSTOM_REAL
@@ -1089,9 +1091,10 @@
!local parameters
real(kind=CUSTOM_REAL),dimension(:),allocatable :: contributions_dot_dot,contributions_dot
real(kind=CUSTOM_REAL):: d
- integer :: ispec,ispecPML,i,j,k,iglob
+ integer :: ispec,ispecPML,i,j,k,iglob,ier
- allocate(contributions_dot_dot(NGLOB_AB),contributions_dot(NGLOB_AB))
+ allocate(contributions_dot_dot(NGLOB_AB),contributions_dot(NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array contributions_dot_dot and contributions_dot'
contributions_dot_dot = 0._CUSTOM_REAL
contributions_dot = 0._CUSTOM_REAL
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/create_color_image.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/create_color_image.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/create_color_image.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -361,7 +361,8 @@
! creating and filling array num_pixel_loc with the positions of each colored
! pixel owned by the local process (useful for parallel jobs)
- allocate(num_pixel_loc(nb_pixel_loc))
+ allocate(num_pixel_loc(nb_pixel_loc),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array num_pixel_loc'
nb_pixel_loc = 0
do i = 1, NX_IMAGE_color
do j = 1, NZ_IMAGE_color
@@ -373,12 +374,14 @@
enddo
! filling array iglob_image_color, containing info on which process owns which pixels.
- allocate(nb_pixel_per_proc(0:NPROC-1))
+ allocate(nb_pixel_per_proc(0:NPROC-1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array nb_pixel_per_proc'
call gather_all_i(nb_pixel_loc,1,nb_pixel_per_proc,1,NPROC)
! allocates receiving array
if ( myrank == 0 ) then
- allocate( num_pixel_recv(maxval(nb_pixel_per_proc(:)),0:NPROC-1) )
+ allocate( num_pixel_recv(maxval(nb_pixel_per_proc(:)),0:NPROC-1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array num_pixel_recv'
endif
! fills iglob_image_color index array
if( NPROC > 1 ) then
@@ -407,7 +410,8 @@
image_color_vp_display(:,:) = 0._CUSTOM_REAL
if ( myrank == 0 ) then
- allocate( data_pixel_recv(maxval(nb_pixel_per_proc(:))) )
+ allocate( data_pixel_recv(maxval(nb_pixel_per_proc(:))),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array data_pixel_recv'
data_pixel_recv(:) = 0._CUSTOM_REAL
endif
allocate(data_pixel_send(nb_pixel_loc),stat=ier)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/detect_mesh_surfaces.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/detect_mesh_surfaces.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/detect_mesh_surfaces.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -33,10 +33,12 @@
use specfem_par_acoustic
use specfem_par_elastic
implicit none
-
+ integer :: ier
+
! for mesh surface
- allocate(ispec_is_surface_external_mesh(NSPEC_AB))
- allocate(iglob_is_surface_external_mesh(NGLOB_AB))
+ allocate(ispec_is_surface_external_mesh(NSPEC_AB), &
+ iglob_is_surface_external_mesh(NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array for mesh surface'
! determines model surface
if (.not. RECVS_CAN_BE_BURIED_EXT_MESH .or. &
@@ -98,16 +100,18 @@
if (MOVIE_VOLUME) then
! acoustic
if( ACOUSTIC_SIMULATION .or. ELASTIC_SIMULATION ) then
- allocate(velocity_x(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
- allocate(velocity_y(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
- allocate(velocity_z(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(velocity_x(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
+ velocity_y(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
+ velocity_z(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array movie velocity_x etc.'
endif
! elastic only
if( ELASTIC_SIMULATION ) then
- allocate(div(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
- allocate(curl_x(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
- allocate(curl_y(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
- allocate(curl_z(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(div(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
+ curl_x(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
+ curl_y(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
+ curl_z(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array movie div and curl'
div(:,:,:,:) = 0._CUSTOM_REAL
curl_x(:,:,:,:) = 0._CUSTOM_REAL
curl_y(:,:,:,:) = 0._CUSTOM_REAL
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/initialize_simulation.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/initialize_simulation.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/initialize_simulation.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -134,28 +134,32 @@
endif
! allocate arrays for storing the databases
- allocate(ibool(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
- allocate(xix(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
- allocate(xiy(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
- allocate(xiz(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
- allocate(etax(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
- allocate(etay(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
- allocate(etaz(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
- allocate(gammax(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
- allocate(gammay(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
- allocate(gammaz(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
- allocate(jacobian(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(ibool(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
+ xix(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
+ xiy(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
+ xiz(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
+ etax(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
+ etay(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
+ etaz(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
+ gammax(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
+ gammay(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
+ gammaz(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
+ jacobian(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating arrays for databases'
! mesh node locations
- allocate(xstore(NGLOB_AB))
- allocate(ystore(NGLOB_AB))
- allocate(zstore(NGLOB_AB))
+ allocate(xstore(NGLOB_AB), &
+ ystore(NGLOB_AB), &
+ zstore(NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating arrays for mesh nodes'
! material properties
- allocate(kappastore(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
- allocate(mustore(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(kappastore(NGLLX,NGLLY,NGLLZ,NSPEC_AB), &
+ mustore(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating arrays for material properties'
! material flags
- allocate(ispec_is_acoustic(NSPEC_AB))
- allocate(ispec_is_elastic(NSPEC_AB))
- allocate(ispec_is_poroelastic(NSPEC_AB))
+ allocate(ispec_is_acoustic(NSPEC_AB), &
+ ispec_is_elastic(NSPEC_AB), &
+ ispec_is_poroelastic(NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating arrays for material flags'
! initializes adjoint simulations
call initialize_simulation_adjoint()
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_receivers.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_receivers.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -132,7 +132,7 @@
double precision, allocatable, dimension(:,:) :: xi_receiver_all,eta_receiver_all,gamma_receiver_all
double precision, allocatable, dimension(:,:,:,:) :: nu_all
-
+ integer :: ier
character(len=256) OUTPUT_FILES
! **************
@@ -164,36 +164,35 @@
if (ios /= 0) call exit_mpi(myrank,'error opening file '//trim(rec_filename))
! allocate memory for arrays using number of stations
- allocate(stlat(nrec))
- allocate(stlon(nrec))
- allocate(stele(nrec))
- allocate(stbur(nrec))
- allocate(stutm_x(nrec))
- allocate(stutm_y(nrec))
- allocate(horiz_dist(nrec))
- allocate(elevation(nrec))
+ allocate(stlat(nrec), &
+ stlon(nrec), &
+ stele(nrec), &
+ stbur(nrec), &
+ stutm_x(nrec), &
+ stutm_y(nrec), &
+ horiz_dist(nrec), &
+ elevation(nrec), &
+ ix_initial_guess(nrec), &
+ iy_initial_guess(nrec), &
+ iz_initial_guess(nrec), &
+ x_target(nrec), &
+ y_target(nrec), &
+ z_target(nrec), &
+ x_found(nrec), &
+ y_found(nrec), &
+ z_found(nrec), &
+ final_distance(nrec), &
+ ispec_selected_rec_all(nrec,0:NPROC-1), &
+ xi_receiver_all(nrec,0:NPROC-1), &
+ eta_receiver_all(nrec,0:NPROC-1), &
+ gamma_receiver_all(nrec,0:NPROC-1), &
+ x_found_all(nrec,0:NPROC-1), &
+ y_found_all(nrec,0:NPROC-1), &
+ z_found_all(nrec,0:NPROC-1), &
+ final_distance_all(nrec,0:NPROC-1), &
+ nu_all(3,3,nrec,0:NPROC-1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating arrays for locating receivers'
- allocate(ix_initial_guess(nrec))
- allocate(iy_initial_guess(nrec))
- allocate(iz_initial_guess(nrec))
- allocate(x_target(nrec))
- allocate(y_target(nrec))
- allocate(z_target(nrec))
- allocate(x_found(nrec))
- allocate(y_found(nrec))
- allocate(z_found(nrec))
- allocate(final_distance(nrec))
-
- allocate(ispec_selected_rec_all(nrec,0:NPROC-1))
- allocate(xi_receiver_all(nrec,0:NPROC-1))
- allocate(eta_receiver_all(nrec,0:NPROC-1))
- allocate(gamma_receiver_all(nrec,0:NPROC-1))
- allocate(x_found_all(nrec,0:NPROC-1))
- allocate(y_found_all(nrec,0:NPROC-1))
- allocate(z_found_all(nrec,0:NPROC-1))
- allocate(final_distance_all(nrec,0:NPROC-1))
- allocate(nu_all(3,3,nrec,0:NPROC-1))
-
! loop on all the stations
do irec=1,nrec
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_source.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_source.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_source.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -152,7 +152,7 @@
integer, dimension(3,NGLLSQUARE,num_free_surface_faces) :: free_surface_ijk
integer ix_initial_guess_source,iy_initial_guess_source,iz_initial_guess_source
-
+ integer ier
integer, dimension(NSOURCES) :: idomain
integer, dimension(NGATHER_SOURCES,0:NPROC-1) :: idomain_all
@@ -672,7 +672,8 @@
ispec_selected_source_all(:,:) = -1
! avoids warnings about temporary creations of arrays for function call by compiler
- allocate(tmp_i_local(ng),tmp_i_all_local(ng,0:NPROC-1))
+ allocate(tmp_i_local(ng),tmp_i_all_local(ng,0:NPROC-1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array tmp_i_local'
tmp_i_local(:) = ispec_selected_source(ns:ne)
call gather_all_i(tmp_i_local,ng,tmp_i_all_local,ng,NPROC)
ispec_selected_source_all(1:ng,:) = tmp_i_all_local(:,:)
@@ -685,7 +686,8 @@
deallocate(tmp_i_local,tmp_i_all_local)
! avoids warnings about temporary creations of arrays for function call by compiler
- allocate(tmp_local(ng),tmp_all_local(ng,0:NPROC-1))
+ allocate(tmp_local(ng),tmp_all_local(ng,0:NPROC-1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array tmp_local'
tmp_local(:) = xi_source(ns:ne)
call gather_all_dp(tmp_local,ng,tmp_all_local,ng,NPROC)
xi_source_all(1:ng,:) = tmp_all_local(:,:)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/prepare_timerun.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -36,7 +36,8 @@
implicit none
character(len=256) :: plot_file
-
+ integer :: ier
+
! flag for any movie simulation
if( EXTERNAL_MESH_MOVIE_SURFACE .or. EXTERNAL_MESH_CREATE_SHAKEMAP .or. &
MOVIE_SURFACE .or. CREATE_SHAKEMAP .or. MOVIE_VOLUME .or. PNM_GIF_IMAGE ) then
@@ -145,9 +146,12 @@
! seismograms
if (nrec_local > 0) then
! allocate seismogram array
- allocate(seismograms_d(NDIM,nrec_local,NSTEP))
- allocate(seismograms_v(NDIM,nrec_local,NSTEP))
- allocate(seismograms_a(NDIM,nrec_local,NSTEP))
+ allocate(seismograms_d(NDIM,nrec_local,NSTEP),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array seismograms_d'
+ allocate(seismograms_v(NDIM,nrec_local,NSTEP),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array seismograms_v'
+ allocate(seismograms_a(NDIM,nrec_local,NSTEP),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array seismograms_a'
! initialize seismograms
seismograms_d(:,:,:) = 0._CUSTOM_REAL
@@ -410,7 +414,7 @@
use specfem_par_elastic
use specfem_par_poroelastic
implicit none
-
+ ! local parameters
integer :: ier
! seismograms
@@ -419,7 +423,8 @@
allocate(Mxx_der(nrec_local),Myy_der(nrec_local), &
Mzz_der(nrec_local),Mxy_der(nrec_local), &
Mxz_der(nrec_local),Myz_der(nrec_local), &
- sloc_der(NDIM,nrec_local))
+ sloc_der(NDIM,nrec_local),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array Mxx_der and following arrays'
Mxx_der = 0._CUSTOM_REAL
Myy_der = 0._CUSTOM_REAL
Mzz_der = 0._CUSTOM_REAL
@@ -428,7 +433,8 @@
Myz_der = 0._CUSTOM_REAL
sloc_der = 0._CUSTOM_REAL
- allocate(seismograms_eps(NDIM,NDIM,nrec_local,NSTEP))
+ allocate(seismograms_eps(NDIM,NDIM,nrec_local,NSTEP),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array seismograms_eps'
seismograms_eps(:,:,:,:) = 0._CUSTOM_REAL
endif
@@ -518,7 +524,8 @@
! elastic domains
if( ELASTIC_SIMULATION) then
! allocates wavefield
- allocate(b_absorb_field(NDIM,NGLLSQUARE,b_num_abs_boundary_faces))
+ allocate(b_absorb_field(NDIM,NGLLSQUARE,b_num_abs_boundary_faces),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array b_absorb_field'
b_reclen_field = CUSTOM_REAL * (NDIM * NGLLSQUARE * num_abs_boundary_faces)
@@ -540,7 +547,8 @@
! acoustic domains
if( ACOUSTIC_SIMULATION) then
! allocates wavefield
- allocate(b_absorb_potential(NGLLSQUARE,b_num_abs_boundary_faces))
+ allocate(b_absorb_potential(NGLLSQUARE,b_num_abs_boundary_faces),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array b_absorb_potential'
b_reclen_potential = CUSTOM_REAL * (NGLLSQUARE * num_abs_boundary_faces)
@@ -562,12 +570,16 @@
else
! dummy array
b_num_abs_boundary_faces = 1
- if( ELASTIC_SIMULATION ) &
- allocate(b_absorb_field(NDIM,NGLLSQUARE,b_num_abs_boundary_faces))
-
- if( ACOUSTIC_SIMULATION ) &
- allocate(b_absorb_potential(NGLLSQUARE,b_num_abs_boundary_faces))
-
+ if( ELASTIC_SIMULATION ) then
+ allocate(b_absorb_field(NDIM,NGLLSQUARE,b_num_abs_boundary_faces),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array b_absorb_field'
+ endif
+
+ if( ACOUSTIC_SIMULATION ) then
+ allocate(b_absorb_potential(NGLLSQUARE,b_num_abs_boundary_faces),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array b_absorb_potential'
+ endif
+
endif
endif
@@ -587,7 +599,7 @@
use specfem_par_poroelastic
use specfem_par_movie
implicit none
-
+ ! local parameters
integer :: ier
! for noise simulations
@@ -597,10 +609,14 @@
allocate(noise_sourcearray(NDIM,NGLLX,NGLLY,NGLLZ,NSTEP),stat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error allocating noise source array')
- allocate(normal_x_noise(NGLLX*NGLLY*nfaces_surface_ext_mesh))
- allocate(normal_y_noise(NGLLX*NGLLY*nfaces_surface_ext_mesh))
- allocate(normal_z_noise(NGLLX*NGLLY*nfaces_surface_ext_mesh))
- allocate(mask_noise(NGLLX*NGLLY*nfaces_surface_ext_mesh))
+ allocate(normal_x_noise(NGLLX*NGLLY*nfaces_surface_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array normal_x_noise'
+ allocate(normal_y_noise(NGLLX*NGLLY*nfaces_surface_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array normal_y_noise'
+ allocate(normal_z_noise(NGLLX*NGLLY*nfaces_surface_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array normal_z_noise'
+ allocate(mask_noise(NGLLX*NGLLY*nfaces_surface_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array mask_noise'
! initializes
noise_sourcearray(:,:,:,:,:) = 0._CUSTOM_REAL
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_mesh_databases.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -83,13 +83,18 @@
call any_all_l( ANY(ispec_is_acoustic), ACOUSTIC_SIMULATION )
if( ACOUSTIC_SIMULATION ) then
! potentials
- allocate(potential_acoustic(NGLOB_AB))
- allocate(potential_dot_acoustic(NGLOB_AB))
- allocate(potential_dot_dot_acoustic(NGLOB_AB))
+ allocate(potential_acoustic(NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array potential_acoustic'
+ allocate(potential_dot_acoustic(NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array potential_dot_acoustic'
+ allocate(potential_dot_dot_acoustic(NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array potential_dot_dot_acoustic'
! mass matrix, density
- allocate(rmass_acoustic(NGLOB_AB))
- allocate(rhostore(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(rmass_acoustic(NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array rmass_acoustic'
+ allocate(rhostore(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array rhostore'
read(27) rmass_acoustic
read(27) rhostore
@@ -99,61 +104,73 @@
call any_all_l( ANY(ispec_is_elastic), ELASTIC_SIMULATION )
if( ELASTIC_SIMULATION ) then
! displacement,velocity,acceleration
- allocate(displ(NDIM,NGLOB_AB))
- allocate(veloc(NDIM,NGLOB_AB))
- allocate(accel(NDIM,NGLOB_AB))
+ allocate(displ(NDIM,NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array displ'
+ allocate(veloc(NDIM,NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array veloc'
+ allocate(accel(NDIM,NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array accel'
- allocate(rmass(NGLOB_AB))
- allocate(rho_vp(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
- allocate(rho_vs(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
- allocate(c11store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
- allocate(c12store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
- allocate(c13store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
- allocate(c14store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
- allocate(c15store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
- allocate(c16store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
- allocate(c22store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
- allocate(c23store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
- allocate(c24store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
- allocate(c25store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
- allocate(c26store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
- allocate(c33store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
- allocate(c34store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
- allocate(c35store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
- allocate(c36store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
- allocate(c44store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
- allocate(c45store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
- allocate(c46store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
- allocate(c55store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
- allocate(c56store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
- allocate(c66store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
+ allocate(rmass(NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array rmass'
+ allocate(rho_vp(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array rho_vp'
+ allocate(rho_vs(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array rho_vs'
+ allocate(c11store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c12store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c13store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c14store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c15store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c16store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c22store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c23store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c24store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c25store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c26store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c33store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c34store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c35store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c36store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c44store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c45store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c46store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c55store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c56store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c66store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array c11store etc.'
! note: currently, they need to be defined, as they are used in the routine arguments
! for compute_forces_elastic_Deville()
- allocate(R_xx(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS))
- allocate(R_yy(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS))
- allocate(R_xy(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS))
- allocate(R_xz(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS))
- allocate(R_yz(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS))
+ allocate(R_xx(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS), &
+ R_yy(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS), &
+ R_xy(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS), &
+ R_xz(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS), &
+ R_yz(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array R_xx etc.'
! needed for attenuation and/or kernel computations
- allocate(epsilondev_xx(NGLLX,NGLLY,NGLLZ,NSPEC_STRAIN_ONLY))
- allocate(epsilondev_yy(NGLLX,NGLLY,NGLLZ,NSPEC_STRAIN_ONLY))
- allocate(epsilondev_xy(NGLLX,NGLLY,NGLLZ,NSPEC_STRAIN_ONLY))
- allocate(epsilondev_xz(NGLLX,NGLLY,NGLLZ,NSPEC_STRAIN_ONLY))
- allocate(epsilondev_yz(NGLLX,NGLLY,NGLLZ,NSPEC_STRAIN_ONLY))
+ allocate(epsilondev_xx(NGLLX,NGLLY,NGLLZ,NSPEC_STRAIN_ONLY), &
+ epsilondev_yy(NGLLX,NGLLY,NGLLZ,NSPEC_STRAIN_ONLY), &
+ epsilondev_xy(NGLLX,NGLLY,NGLLZ,NSPEC_STRAIN_ONLY), &
+ epsilondev_xz(NGLLX,NGLLY,NGLLZ,NSPEC_STRAIN_ONLY), &
+ epsilondev_yz(NGLLX,NGLLY,NGLLZ,NSPEC_STRAIN_ONLY),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array epsilondev_xx etc.'
! note: needed for argument of deville routine
- allocate(epsilon_trace_over_3(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT))
+ allocate(epsilon_trace_over_3(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array epsilon_trace_over_3'
! needed for attenuation
- allocate(one_minus_sum_beta(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB))
- allocate(factor_common(N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB))
+ allocate(one_minus_sum_beta(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB), &
+ factor_common(N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array one_minus_sum_beta etc.'
read(27) rmass
if( OCEANS ) then
! ocean mass matrix
- allocate(rmass_ocean_load(NGLOB_AB))
+ allocate(rmass_ocean_load(NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array rmass_ocean_load'
read(27) rmass_ocean_load
endif
!pll
@@ -172,8 +189,10 @@
stop 'not implemented yet: read rmass_solid_poroelastic .. '
- allocate(rmass_solid_poroelastic(NGLOB_AB))
- allocate(rmass_fluid_poroelastic(NGLOB_AB))
+ allocate(rmass_solid_poroelastic(NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array rmass_solid_poroelastic'
+ allocate(rmass_fluid_poroelastic(NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array rmass_fluid_poroelastic'
read(27) rmass_solid_poroelastic
read(27) rmass_fluid_poroelastic
@@ -189,10 +208,11 @@
! absorbing boundary surface
read(27) num_abs_boundary_faces
- allocate(abs_boundary_ispec(num_abs_boundary_faces))
- allocate(abs_boundary_ijk(3,NGLLSQUARE,num_abs_boundary_faces))
- allocate(abs_boundary_jacobian2Dw(NGLLSQUARE,num_abs_boundary_faces))
- allocate(abs_boundary_normal(NDIM,NGLLSQUARE,num_abs_boundary_faces))
+ allocate(abs_boundary_ispec(num_abs_boundary_faces), &
+ abs_boundary_ijk(3,NGLLSQUARE,num_abs_boundary_faces), &
+ abs_boundary_jacobian2Dw(NGLLSQUARE,num_abs_boundary_faces), &
+ abs_boundary_normal(NDIM,NGLLSQUARE,num_abs_boundary_faces),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array abs_boundary_ispec etc.'
read(27) abs_boundary_ispec
read(27) abs_boundary_ijk
read(27) abs_boundary_jacobian2Dw
@@ -200,10 +220,11 @@
! free surface
read(27) num_free_surface_faces
- allocate(free_surface_ispec(num_free_surface_faces))
- allocate(free_surface_ijk(3,NGLLSQUARE,num_free_surface_faces))
- allocate(free_surface_jacobian2Dw(NGLLSQUARE,num_free_surface_faces))
- allocate(free_surface_normal(NDIM,NGLLSQUARE,num_free_surface_faces))
+ allocate(free_surface_ispec(num_free_surface_faces), &
+ free_surface_ijk(3,NGLLSQUARE,num_free_surface_faces), &
+ free_surface_jacobian2Dw(NGLLSQUARE,num_free_surface_faces), &
+ free_surface_normal(NDIM,NGLLSQUARE,num_free_surface_faces),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array free_surface_ispec etc.'
read(27) free_surface_ispec
read(27) free_surface_ijk
read(27) free_surface_jacobian2Dw
@@ -211,10 +232,11 @@
! acoustic-elastic coupling surface
read(27) num_coupling_ac_el_faces
- allocate(coupling_ac_el_normal(NDIM,NGLLSQUARE,num_coupling_ac_el_faces))
- allocate(coupling_ac_el_jacobian2Dw(NGLLSQUARE,num_coupling_ac_el_faces))
- allocate(coupling_ac_el_ijk(3,NGLLSQUARE,num_coupling_ac_el_faces))
- allocate(coupling_ac_el_ispec(num_coupling_ac_el_faces))
+ allocate(coupling_ac_el_normal(NDIM,NGLLSQUARE,num_coupling_ac_el_faces), &
+ coupling_ac_el_jacobian2Dw(NGLLSQUARE,num_coupling_ac_el_faces), &
+ coupling_ac_el_ijk(3,NGLLSQUARE,num_coupling_ac_el_faces), &
+ coupling_ac_el_ispec(num_coupling_ac_el_faces),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array coupling_ac_el_normal etc.'
read(27) coupling_ac_el_ispec
read(27) coupling_ac_el_ijk
read(27) coupling_ac_el_jacobian2Dw
@@ -223,9 +245,10 @@
! MPI interfaces
read(27) num_interfaces_ext_mesh
read(27) max_nibool_interfaces_ext_mesh
- allocate(my_neighbours_ext_mesh(num_interfaces_ext_mesh))
- allocate(nibool_interfaces_ext_mesh(num_interfaces_ext_mesh))
- allocate(ibool_interfaces_ext_mesh(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh))
+ allocate(my_neighbours_ext_mesh(num_interfaces_ext_mesh), &
+ nibool_interfaces_ext_mesh(num_interfaces_ext_mesh), &
+ ibool_interfaces_ext_mesh(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array my_neighbours_ext_mesh etc.'
read(27) my_neighbours_ext_mesh
read(27) nibool_interfaces_ext_mesh
read(27) ibool_interfaces_ext_mesh
@@ -257,18 +280,21 @@
close(27)
! MPI communications
- allocate(buffer_send_vector_ext_mesh(NDIM,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh))
- allocate(buffer_recv_vector_ext_mesh(NDIM,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh))
- allocate(buffer_send_scalar_ext_mesh(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh))
- allocate(buffer_recv_scalar_ext_mesh(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh))
- allocate(request_send_vector_ext_mesh(num_interfaces_ext_mesh))
- allocate(request_recv_vector_ext_mesh(num_interfaces_ext_mesh))
- allocate(request_send_scalar_ext_mesh(num_interfaces_ext_mesh))
- allocate(request_recv_scalar_ext_mesh(num_interfaces_ext_mesh))
+ allocate(buffer_send_vector_ext_mesh(NDIM,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh), &
+ buffer_recv_vector_ext_mesh(NDIM,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh), &
+ buffer_send_scalar_ext_mesh(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh), &
+ buffer_recv_scalar_ext_mesh(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh), &
+ request_send_vector_ext_mesh(num_interfaces_ext_mesh), &
+ request_recv_vector_ext_mesh(num_interfaces_ext_mesh), &
+ request_send_scalar_ext_mesh(num_interfaces_ext_mesh), &
+ request_recv_scalar_ext_mesh(num_interfaces_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array buffer_send_vector_ext_mesh etc.'
! locate inner and outer elements
- allocate(ispec_is_inner(NSPEC_AB))
- allocate(iglob_is_inner(NGLOB_AB))
+ allocate(ispec_is_inner(NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ispec_is_inner'
+ allocate(iglob_is_inner(NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array iglob_is_inner'
ispec_is_inner(:) = .true.
iglob_is_inner(:) = .true.
do iinterface = 1, num_interfaces_ext_mesh
@@ -306,7 +332,8 @@
! stores indices of inner and outer elements for faster(?) computation
num_phase_ispec_acoustic = max(nspec_inner_acoustic,nspec_outer_acoustic)
- allocate( phase_ispec_inner_acoustic(num_phase_ispec_acoustic,2))
+ allocate( phase_ispec_inner_acoustic(num_phase_ispec_acoustic,2),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array phase_ispec_inner_acoustic'
nspec_inner_acoustic = 0
nspec_outer_acoustic = 0
do ispec = 1, NSPEC_AB
@@ -341,7 +368,8 @@
! stores indices of inner and outer elements for faster(?) computation
num_phase_ispec_elastic = max(nspec_inner_elastic,nspec_outer_elastic)
- allocate( phase_ispec_inner_elastic(num_phase_ispec_elastic,2))
+ allocate( phase_ispec_inner_elastic(num_phase_ispec_elastic,2),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array phase_ispec_inner_elastic'
nspec_inner_elastic = 0
nspec_outer_elastic = 0
do ispec = 1, NSPEC_AB
@@ -381,8 +409,10 @@
call check_mesh_resolution(myrank,NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore, &
kappastore,mustore,rho_vp,rho_vs,DT,model_speed_max,min_resolved_period )
else if( ACOUSTIC_SIMULATION ) then
- allocate(rho_vp(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
- allocate(rho_vs(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(rho_vp(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array rho_vp'
+ allocate(rho_vs(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array rho_vs'
rho_vp = sqrt( kappastore / rhostore ) * rhostore
rho_vs = 0.0_CUSTOM_REAL
call check_mesh_resolution(myrank,NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore, &
@@ -414,36 +444,49 @@
! allocates adjoint arrays for elastic simulations
if( ELASTIC_SIMULATION .and. SIMULATION_TYPE == 3 ) then
! backward displacement,velocity,acceleration fields
- allocate(b_displ(NDIM,NGLOB_ADJOINT))
- allocate(b_veloc(NDIM,NGLOB_ADJOINT))
- allocate(b_accel(NDIM,NGLOB_ADJOINT))
+ allocate(b_displ(NDIM,NGLOB_ADJOINT),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array b_displ'
+ allocate(b_veloc(NDIM,NGLOB_ADJOINT),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array b_veloc'
+ allocate(b_accel(NDIM,NGLOB_ADJOINT),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array b_accel'
! adjoint kernels
! primary, isotropic kernels
! density kernel
- allocate(rho_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT))
+ allocate(rho_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array rho_kl'
! shear modulus kernel
- allocate(mu_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT))
+ allocate(mu_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array mu_kl'
! compressional modulus kernel
- allocate(kappa_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT))
+ allocate(kappa_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array kappa_kl'
! derived kernels
! density prime kernel
- allocate(rhop_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT))
+ allocate(rhop_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array rhop_kl'
! vp kernel
- allocate(alpha_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT))
+ allocate(alpha_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array alpha_kl'
! vs kernel
- allocate(beta_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT))
+ allocate(beta_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array beta_kl'
! noise source strength kernel
- if (NOISE_TOMOGRAPHY == 3) allocate(sigma_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT))
+ if (NOISE_TOMOGRAPHY == 3) then
+ allocate(sigma_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array sigma_kl'
+ endif
! MPI handling
- allocate(b_request_send_vector_ext_mesh(num_interfaces_ext_mesh))
- allocate(b_request_recv_vector_ext_mesh(num_interfaces_ext_mesh))
- allocate(b_buffer_send_vector_ext_mesh(NDIM,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh))
- allocate(b_buffer_recv_vector_ext_mesh(NDIM,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh))
+ allocate(b_request_send_vector_ext_mesh(num_interfaces_ext_mesh), &
+ b_request_recv_vector_ext_mesh(num_interfaces_ext_mesh), &
+ b_buffer_send_vector_ext_mesh(NDIM,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh), &
+ b_buffer_recv_vector_ext_mesh(NDIM,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array b_request_send_vector_ext_mesh etc.'
! allocates attenuation solids
if( ATTENUATION ) then
@@ -451,7 +494,8 @@
b_R_yy(NGLLX,NGLLY,NGLLZ,NSPEC_ATT_AND_KERNEL,N_SLS), &
b_R_xy(NGLLX,NGLLY,NGLLZ,NSPEC_ATT_AND_KERNEL,N_SLS), &
b_R_xz(NGLLX,NGLLY,NGLLZ,NSPEC_ATT_AND_KERNEL,N_SLS), &
- b_R_yz(NGLLX,NGLLY,NGLLZ,NSPEC_ATT_AND_KERNEL,N_SLS) )
+ b_R_yz(NGLLX,NGLLY,NGLLZ,NSPEC_ATT_AND_KERNEL,N_SLS),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array b_R_xx etc.'
endif
! note: these arrays are needed for attenuation and/or kernel computations
@@ -459,57 +503,76 @@
b_epsilondev_yy(NGLLX,NGLLY,NGLLZ,NSPEC_STRAIN_ONLY), &
b_epsilondev_xy(NGLLX,NGLLY,NGLLZ,NSPEC_STRAIN_ONLY), &
b_epsilondev_xz(NGLLX,NGLLY,NGLLZ,NSPEC_STRAIN_ONLY), &
- b_epsilondev_yz(NGLLX,NGLLY,NGLLZ,NSPEC_STRAIN_ONLY) )
+ b_epsilondev_yz(NGLLX,NGLLY,NGLLZ,NSPEC_STRAIN_ONLY),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array b_epsilon_dev_xx etc.'
! needed for kernel computations
- allocate(b_epsilon_trace_over_3(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT))
+ allocate(b_epsilon_trace_over_3(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array b_epsilon_trace_over_3'
+ else
+ ! modification: Camille Mazoyer
+ ! dummy allocation
+ allocate(b_displ(1,1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating dummy array b_displ'
+ allocate(b_veloc(1,1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating dummy array b_veloc'
+ allocate(b_accel(1,1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating dummy array b_accel'
+
endif
! allocates adjoint arrays for acoustic simulations
if( ACOUSTIC_SIMULATION .and. SIMULATION_TYPE == 3 ) then
! backward potentials
- allocate(b_potential_acoustic(NGLOB_ADJOINT))
- allocate(b_potential_dot_acoustic(NGLOB_ADJOINT))
- allocate(b_potential_dot_dot_acoustic(NGLOB_ADJOINT))
+ allocate(b_potential_acoustic(NGLOB_ADJOINT), &
+ b_potential_dot_acoustic(NGLOB_ADJOINT), &
+ b_potential_dot_dot_acoustic(NGLOB_ADJOINT),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array b_potential_acoustic etc.'
! kernels
- allocate(rho_ac_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT))
- allocate(rhop_ac_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT))
- allocate(kappa_ac_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT))
- allocate(alpha_ac_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT))
+ allocate(rho_ac_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ rhop_ac_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ kappa_ac_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT), &
+ alpha_ac_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array rho_ac_kl etc.'
! MPI handling
- allocate(b_request_send_scalar_ext_mesh(num_interfaces_ext_mesh))
- allocate(b_request_recv_scalar_ext_mesh(num_interfaces_ext_mesh))
- allocate(b_buffer_send_scalar_ext_mesh(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh))
- allocate(b_buffer_recv_scalar_ext_mesh(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh))
+ allocate(b_request_send_scalar_ext_mesh(num_interfaces_ext_mesh), &
+ b_request_recv_scalar_ext_mesh(num_interfaces_ext_mesh), &
+ b_buffer_send_scalar_ext_mesh(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh), &
+ b_buffer_recv_scalar_ext_mesh(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array b_request_send_scalar_ext_mesh'
else
! backward potentials
- allocate(b_potential_acoustic(1))
- allocate(b_potential_dot_acoustic(1))
- allocate(b_potential_dot_dot_acoustic(1))
+ allocate(b_potential_acoustic(1), &
+ b_potential_dot_acoustic(1), &
+ b_potential_dot_dot_acoustic(1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating dummy array b_potential_acoustic etc.'
! kernels
- allocate(rho_ac_kl(1,1,1,1))
- allocate(rhop_ac_kl(1,1,1,1))
- allocate(kappa_ac_kl(1,1,1,1))
- allocate(alpha_ac_kl(1,1,1,1))
+ allocate(rho_ac_kl(1,1,1,1), &
+ rhop_ac_kl(1,1,1,1), &
+ kappa_ac_kl(1,1,1,1), &
+ alpha_ac_kl(1,1,1,1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating dummy array rho_ac_kl etc.'
! MPI handling
- allocate(b_request_send_scalar_ext_mesh(1))
- allocate(b_request_recv_scalar_ext_mesh(1))
- allocate(b_buffer_send_scalar_ext_mesh(1,1))
- allocate(b_buffer_recv_scalar_ext_mesh(1,1))
+ allocate(b_request_send_scalar_ext_mesh(1), &
+ b_request_recv_scalar_ext_mesh(1), &
+ b_buffer_send_scalar_ext_mesh(1,1), &
+ b_buffer_recv_scalar_ext_mesh(1,1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating dummy array b_request_send_scalar_ext_mesh etc.'
endif
! ADJOINT moho
! moho boundary
if( ELASTIC_SIMULATION ) then
- allocate( is_moho_top(NSPEC_BOUN),is_moho_bot(NSPEC_BOUN) )
+ allocate( is_moho_top(NSPEC_BOUN),is_moho_bot(NSPEC_BOUN),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array is_moho_top etc.'
if( SAVE_MOHO_MESH .and. SIMULATION_TYPE == 3 ) then
@@ -526,12 +589,13 @@
read(27) NSPEC2D_MOHO
! allocates arrays for moho mesh
- allocate(ibelm_moho_bot(NSPEC2D_MOHO))
- allocate(ibelm_moho_top(NSPEC2D_MOHO))
- allocate(normal_moho_top(NDIM,NGLLSQUARE,NSPEC2D_MOHO))
- allocate(normal_moho_bot(NDIM,NGLLSQUARE,NSPEC2D_MOHO))
- allocate(ijk_moho_bot(3,NGLLSQUARE,NSPEC2D_MOHO))
- allocate(ijk_moho_top(3,NGLLSQUARE,NSPEC2D_MOHO))
+ allocate(ibelm_moho_bot(NSPEC2D_MOHO), &
+ ibelm_moho_top(NSPEC2D_MOHO), &
+ normal_moho_top(NDIM,NGLLSQUARE,NSPEC2D_MOHO), &
+ normal_moho_bot(NDIM,NGLLSQUARE,NSPEC2D_MOHO), &
+ ijk_moho_bot(3,NGLLSQUARE,NSPEC2D_MOHO), &
+ ijk_moho_top(3,NGLLSQUARE,NSPEC2D_MOHO),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ibelm_moho_bot etc.'
read(27) ibelm_moho_top
read(27) ibelm_moho_bot
@@ -568,7 +632,8 @@
close(27)
! moho kernel
- allocate( moho_kl(NGLLSQUARE,NSPEC2D_MOHO) )
+ allocate( moho_kl(NGLLSQUARE,NSPEC2D_MOHO),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array moho_kl'
moho_kl = 0._CUSTOM_REAL
else
@@ -576,9 +641,10 @@
endif
allocate( dsdx_top(NDIM,NDIM,NGLLX,NGLLY,NGLLZ,NSPEC2D_MOHO), &
- dsdx_bot(NDIM,NDIM,NGLLX,NGLLY,NGLLZ,NSPEC2D_MOHO), &
- b_dsdx_top(NDIM,NDIM,NGLLX,NGLLY,NGLLZ,NSPEC2D_MOHO), &
- b_dsdx_bot(NDIM,NDIM,NGLLX,NGLLY,NGLLZ,NSPEC2D_MOHO) )
+ dsdx_bot(NDIM,NDIM,NGLLX,NGLLY,NGLLZ,NSPEC2D_MOHO), &
+ b_dsdx_top(NDIM,NDIM,NGLLX,NGLLY,NGLLZ,NSPEC2D_MOHO), &
+ b_dsdx_bot(NDIM,NDIM,NGLLX,NGLLY,NGLLZ,NSPEC2D_MOHO),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array dsdx_top etc.'
endif
end subroutine read_mesh_databases_adjoint
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_topography_bathymetry.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_topography_bathymetry.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/read_topography_bathymetry.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -30,7 +30,8 @@
use specfem_par
implicit none
-
+ integer :: ier
+
! read topography and bathymetry file
if( OCEANS .and. TOPOGRAPHY ) then
@@ -42,7 +43,8 @@
DEGREES_PER_CELL_TOPO = DEGREES_PER_CELL_TOPO_SOCAL
topo_file = TOPO_FILE_SOCAL
- allocate(itopo_bathy(NX_TOPO,NY_TOPO))
+ allocate(itopo_bathy(NX_TOPO,NY_TOPO),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array itopo_bathy'
call read_topo_bathy_file(itopo_bathy,NX_TOPO,NY_TOPO,topo_file)
@@ -56,7 +58,8 @@
else
NX_TOPO = 1
NY_TOPO = 1
- allocate(itopo_bathy(NX_TOPO,NY_TOPO))
+ allocate(itopo_bathy(NX_TOPO,NY_TOPO),stat=ier)
+ if( ier /= 0 ) stop 'error allocating dummy array itopo_bathy'
endif
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/save_adjoint_kernels.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/save_adjoint_kernels.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/save_adjoint_kernels.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -34,7 +34,7 @@
use specfem_par_movie,only: nfaces_surface_ext_mesh
implicit none
- integer:: ispec,i,j,k,iglob
+ integer:: ispec,i,j,k,iglob,ier
real(kind=CUSTOM_REAL) :: rhol,mul,kappal
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: weights_kernel
! flag to save GLL weights
@@ -113,27 +113,34 @@
! save kernels to binary files
if( ELASTIC_SIMULATION ) then
- open(unit=27,file=prname(1:len_trim(prname))//'rho_kernel.bin',status='unknown',form='unformatted')
+ open(unit=27,file=prname(1:len_trim(prname))//'rho_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file rho_kernel.bin'
write(27) rho_kl
close(27)
- open(unit=27,file=prname(1:len_trim(prname))//'mu_kernel.bin',status='unknown',form='unformatted')
+ open(unit=27,file=prname(1:len_trim(prname))//'mu_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file mu_kernel.bin'
write(27) mu_kl
close(27)
- open(unit=27,file=prname(1:len_trim(prname))//'kappa_kernel.bin',status='unknown',form='unformatted')
+ open(unit=27,file=prname(1:len_trim(prname))//'kappa_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file kappa_kernel.bin'
write(27) kappa_kl
close(27)
- open(unit=27,file=prname(1:len_trim(prname))//'rhop_kernel.bin',status='unknown',form='unformatted')
+ open(unit=27,file=prname(1:len_trim(prname))//'rhop_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file rhop_kernel.bin'
write(27) rhop_kl
close(27)
- open(unit=27,file=prname(1:len_trim(prname))//'beta_kernel.bin',status='unknown',form='unformatted')
+ open(unit=27,file=prname(1:len_trim(prname))//'beta_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file beta_kernel.bin'
write(27) beta_kl
close(27)
- open(unit=27,file=prname(1:len_trim(prname))//'alpha_kernel.bin',status='unknown',form='unformatted')
+ open(unit=27,file=prname(1:len_trim(prname))//'alpha_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file alpha_kernel.bin'
write(27) alpha_kl
close(27)
if (SAVE_MOHO_MESH) then
- open(unit=27,file=prname(1:len_trim(prname))//'moho_kernel.bin',status='unknown',form='unformatted')
+ open(unit=27,file=prname(1:len_trim(prname))//'moho_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file moho_kernel.bin'
write(27) moho_kl
close(27)
endif
@@ -143,16 +150,20 @@
! save kernels to binary files
if( ACOUSTIC_SIMULATION ) then
- open(unit=27,file=prname(1:len_trim(prname))//'rho_acoustic_kernel.bin',status='unknown',form='unformatted')
+ open(unit=27,file=prname(1:len_trim(prname))//'rho_acoustic_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file rho_acoustic_kernel.bin'
write(27) rho_ac_kl
close(27)
- open(unit=27,file=prname(1:len_trim(prname))//'kappa_acoustic_kernel.bin',status='unknown',form='unformatted')
+ open(unit=27,file=prname(1:len_trim(prname))//'kappa_acoustic_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file kappa_acoustic_kernel.bin'
write(27) kappa_ac_kl
close(27)
- open(unit=27,file=prname(1:len_trim(prname))//'rhop_acoustic_kernel.bin',status='unknown',form='unformatted')
+ open(unit=27,file=prname(1:len_trim(prname))//'rhop_acoustic_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file rhop_acoustic_kernel.bin'
write(27) rhop_ac_kl
close(27)
- open(unit=27,file=prname(1:len_trim(prname))//'alpha_acoustic_kernel.bin',status='unknown',form='unformatted')
+ open(unit=27,file=prname(1:len_trim(prname))//'alpha_acoustic_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file alpha_acoustic_kernel.bin'
write(27) alpha_ac_kl
close(27)
@@ -160,7 +171,8 @@
! save weights for volume integration, in order to benchmark the kernels with analytical expressions
if( SAVE_WEIGHTS ) then
- allocate(weights_kernel(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(weights_kernel(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array weights_kernel'
do ispec = 1, NSPEC_AB
do k = 1, NGLLZ
do j = 1, NGLLY
@@ -170,7 +182,8 @@
enddo ! j
enddo ! k
enddo ! ispec
- open(unit=27,file=prname(1:len_trim(prname))//'weights_kernel.bin',status='unknown',form='unformatted')
+ open(unit=27,file=prname(1:len_trim(prname))//'weights_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file weights_kernel.bin'
write(27) weights_kernel
close(27)
endif
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_GLL_points.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_GLL_points.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_GLL_points.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -30,7 +30,7 @@
use specfem_par
implicit none
- integer :: i,j
+ integer :: i,j,ier
if(myrank == 0) then
write(IMAIN,*) '******************************************'
@@ -54,12 +54,13 @@
enddo
! allocate 1-D Lagrange interpolators and derivatives
- allocate(hxir(NGLLX))
- allocate(hpxir(NGLLX))
- allocate(hetar(NGLLY))
- allocate(hpetar(NGLLY))
- allocate(hgammar(NGLLZ))
- allocate(hpgammar(NGLLZ))
+ allocate(hxir(NGLLX), &
+ hpxir(NGLLX), &
+ hetar(NGLLY), &
+ hpetar(NGLLY), &
+ hgammar(NGLLZ), &
+ hpgammar(NGLLZ),stat=ier)
+ if( ier /= 0 ) stop 'error allocating arrays for interpolators'
! create name of database
call create_name_database(prname,myrank,LOCAL_PATH)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_movie_meshes.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_movie_meshes.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_movie_meshes.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -34,50 +34,56 @@
use specfem_par_movie
implicit none
- integer :: i,j,k,ispec,iglob
+ integer :: i,j,k,ispec,iglob,ier
integer :: ipoin,nfaces_org
character(len=256):: filename
! initializes mesh arrays for movies and shakemaps
- allocate(nfaces_perproc_surface_ext_mesh(NPROC))
- allocate(faces_surface_offset_ext_mesh(NPROC))
+ allocate(nfaces_perproc_surface_ext_mesh(NPROC), &
+ faces_surface_offset_ext_mesh(NPROC),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array for movie faces'
+
nfaces_org = nfaces_surface_ext_mesh
if (nfaces_surface_ext_mesh == 0) then
! dummy arrays
if (USE_HIGHRES_FOR_MOVIES) then
- allocate(faces_surface_ext_mesh(NGLLX*NGLLY,1))
- allocate(store_val_x_external_mesh(NGLLX*NGLLY*1))
- allocate(store_val_y_external_mesh(NGLLX*NGLLY*1))
- allocate(store_val_z_external_mesh(NGLLX*NGLLY*1))
- allocate(store_val_ux_external_mesh(NGLLX*NGLLY*1))
- allocate(store_val_uy_external_mesh(NGLLX*NGLLY*1))
- allocate(store_val_uz_external_mesh(NGLLX*NGLLY*1))
+ allocate(faces_surface_ext_mesh(NGLLX*NGLLY,1), &
+ store_val_x_external_mesh(NGLLX*NGLLY*1), &
+ store_val_y_external_mesh(NGLLX*NGLLY*1), &
+ store_val_z_external_mesh(NGLLX*NGLLY*1), &
+ store_val_ux_external_mesh(NGLLX*NGLLY*1), &
+ store_val_uy_external_mesh(NGLLX*NGLLY*1), &
+ store_val_uz_external_mesh(NGLLX*NGLLY*1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating dummy arrays for highres movie'
else
- allocate(faces_surface_ext_mesh(NGNOD2D,1))
- allocate(store_val_x_external_mesh(NGNOD2D*1))
- allocate(store_val_y_external_mesh(NGNOD2D*1))
- allocate(store_val_z_external_mesh(NGNOD2D*1))
- allocate(store_val_ux_external_mesh(NGNOD2D*1))
- allocate(store_val_uy_external_mesh(NGNOD2D*1))
- allocate(store_val_uz_external_mesh(NGNOD2D*1))
+ allocate(faces_surface_ext_mesh(NGNOD2D,1), &
+ store_val_x_external_mesh(NGNOD2D*1), &
+ store_val_y_external_mesh(NGNOD2D*1), &
+ store_val_z_external_mesh(NGNOD2D*1), &
+ store_val_ux_external_mesh(NGNOD2D*1), &
+ store_val_uy_external_mesh(NGNOD2D*1), &
+ store_val_uz_external_mesh(NGNOD2D*1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating dummy arrays for lowres movie'
endif
else
if (USE_HIGHRES_FOR_MOVIES) then
- allocate(faces_surface_ext_mesh(NGLLX*NGLLY,nfaces_surface_ext_mesh))
- allocate(store_val_x_external_mesh(NGLLX*NGLLY*nfaces_surface_ext_mesh))
- allocate(store_val_y_external_mesh(NGLLX*NGLLY*nfaces_surface_ext_mesh))
- allocate(store_val_z_external_mesh(NGLLX*NGLLY*nfaces_surface_ext_mesh))
- allocate(store_val_ux_external_mesh(NGLLX*NGLLY*nfaces_surface_ext_mesh))
- allocate(store_val_uy_external_mesh(NGLLX*NGLLY*nfaces_surface_ext_mesh))
- allocate(store_val_uz_external_mesh(NGLLX*NGLLY*nfaces_surface_ext_mesh))
+ allocate(faces_surface_ext_mesh(NGLLX*NGLLY,nfaces_surface_ext_mesh), &
+ store_val_x_external_mesh(NGLLX*NGLLY*nfaces_surface_ext_mesh), &
+ store_val_y_external_mesh(NGLLX*NGLLY*nfaces_surface_ext_mesh), &
+ store_val_z_external_mesh(NGLLX*NGLLY*nfaces_surface_ext_mesh), &
+ store_val_ux_external_mesh(NGLLX*NGLLY*nfaces_surface_ext_mesh), &
+ store_val_uy_external_mesh(NGLLX*NGLLY*nfaces_surface_ext_mesh), &
+ store_val_uz_external_mesh(NGLLX*NGLLY*nfaces_surface_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating arrays for highres movie'
else
- allocate(faces_surface_ext_mesh(NGNOD2D,nfaces_surface_ext_mesh))
- allocate(store_val_x_external_mesh(NGNOD2D*nfaces_surface_ext_mesh))
- allocate(store_val_y_external_mesh(NGNOD2D*nfaces_surface_ext_mesh))
- allocate(store_val_z_external_mesh(NGNOD2D*nfaces_surface_ext_mesh))
- allocate(store_val_ux_external_mesh(NGNOD2D*nfaces_surface_ext_mesh))
- allocate(store_val_uy_external_mesh(NGNOD2D*nfaces_surface_ext_mesh))
- allocate(store_val_uz_external_mesh(NGNOD2D*nfaces_surface_ext_mesh))
+ allocate(faces_surface_ext_mesh(NGNOD2D,nfaces_surface_ext_mesh), &
+ store_val_x_external_mesh(NGNOD2D*nfaces_surface_ext_mesh), &
+ store_val_y_external_mesh(NGNOD2D*nfaces_surface_ext_mesh), &
+ store_val_z_external_mesh(NGNOD2D*nfaces_surface_ext_mesh), &
+ store_val_ux_external_mesh(NGNOD2D*nfaces_surface_ext_mesh), &
+ store_val_uy_external_mesh(NGNOD2D*nfaces_surface_ext_mesh), &
+ store_val_uz_external_mesh(NGNOD2D*nfaces_surface_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating arrays for lowres movie'
endif
endif
store_val_ux_external_mesh(:) = 0._CUSTOM_REAL
@@ -90,19 +96,21 @@
! arrays used for collected/gathered fields
if (myrank == 0) then
if (USE_HIGHRES_FOR_MOVIES) then
- allocate(store_val_x_all_external_mesh(NGLLX*NGLLY*nfaces_surface_glob_ext_mesh))
- allocate(store_val_y_all_external_mesh(NGLLX*NGLLY*nfaces_surface_glob_ext_mesh))
- allocate(store_val_z_all_external_mesh(NGLLX*NGLLY*nfaces_surface_glob_ext_mesh))
- allocate(store_val_ux_all_external_mesh(NGLLX*NGLLY*nfaces_surface_glob_ext_mesh))
- allocate(store_val_uy_all_external_mesh(NGLLX*NGLLY*nfaces_surface_glob_ext_mesh))
- allocate(store_val_uz_all_external_mesh(NGLLX*NGLLY*nfaces_surface_glob_ext_mesh))
+ allocate(store_val_x_all_external_mesh(NGLLX*NGLLY*nfaces_surface_glob_ext_mesh), &
+ store_val_y_all_external_mesh(NGLLX*NGLLY*nfaces_surface_glob_ext_mesh), &
+ store_val_z_all_external_mesh(NGLLX*NGLLY*nfaces_surface_glob_ext_mesh), &
+ store_val_ux_all_external_mesh(NGLLX*NGLLY*nfaces_surface_glob_ext_mesh), &
+ store_val_uy_all_external_mesh(NGLLX*NGLLY*nfaces_surface_glob_ext_mesh), &
+ store_val_uz_all_external_mesh(NGLLX*NGLLY*nfaces_surface_glob_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating arrays for highres movie'
else
- allocate(store_val_x_all_external_mesh(NGNOD2D*nfaces_surface_glob_ext_mesh))
- allocate(store_val_y_all_external_mesh(NGNOD2D*nfaces_surface_glob_ext_mesh))
- allocate(store_val_z_all_external_mesh(NGNOD2D*nfaces_surface_glob_ext_mesh))
- allocate(store_val_ux_all_external_mesh(NGNOD2D*nfaces_surface_glob_ext_mesh))
- allocate(store_val_uy_all_external_mesh(NGNOD2D*nfaces_surface_glob_ext_mesh))
- allocate(store_val_uz_all_external_mesh(NGNOD2D*nfaces_surface_glob_ext_mesh))
+ allocate(store_val_x_all_external_mesh(NGNOD2D*nfaces_surface_glob_ext_mesh), &
+ store_val_y_all_external_mesh(NGNOD2D*nfaces_surface_glob_ext_mesh), &
+ store_val_z_all_external_mesh(NGNOD2D*nfaces_surface_glob_ext_mesh), &
+ store_val_ux_all_external_mesh(NGNOD2D*nfaces_surface_glob_ext_mesh), &
+ store_val_uy_all_external_mesh(NGNOD2D*nfaces_surface_glob_ext_mesh), &
+ store_val_uz_all_external_mesh(NGNOD2D*nfaces_surface_glob_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating arrays for lowres movie'
endif
endif
call gather_all_i(nfaces_surface_ext_mesh,1,nfaces_perproc_surface_ext_mesh,1,NPROC)
@@ -121,7 +129,8 @@
! stores global indices of GLL points on the surface to array faces_surface_ext_mesh
if( EXTERNAL_MESH_MOVIE_SURFACE .or. EXTERNAL_MESH_CREATE_SHAKEMAP ) then
- allocate( faces_surface_ext_mesh_ispec(nfaces_surface_ext_mesh))
+ allocate( faces_surface_ext_mesh_ispec(nfaces_surface_ext_mesh),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array faces_surface_ext_mesh_ispec'
! stores global indices
nfaces_surface_ext_mesh = 0
@@ -287,6 +296,18 @@
endif
+ ! for gathering movie data
+ if (USE_HIGHRES_FOR_MOVIES) then
+ ! hi-res movies output all gll points on surface
+ nfaces_perproc_surface_ext_mesh(:) = nfaces_perproc_surface_ext_mesh(:)*NGLLX*NGLLY
+ nfaces_surface_ext_mesh_points = nfaces_surface_ext_mesh*NGLLX*NGLLY
+ nfaces_surface_glob_em_points = nfaces_surface_glob_ext_mesh*NGLLX*NGLLY
+ else
+ ! low-res movies only output at element corners
+ nfaces_perproc_surface_ext_mesh(:) = nfaces_perproc_surface_ext_mesh(:)*NGNOD2D
+ nfaces_surface_ext_mesh_points = nfaces_surface_ext_mesh*NGNOD2D
+ nfaces_surface_glob_em_points = nfaces_surface_glob_ext_mesh*NGNOD2D
+ endif
end subroutine setup_movie_meshes
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -72,26 +72,27 @@
double precision :: t0_acoustic,min_tshift_cmt_original
integer :: yr,jda,ho,mi
- integer :: isource,ispec
+ integer :: isource,ispec,ier
! allocate arrays for source
- allocate(islice_selected_source(NSOURCES))
- allocate(ispec_selected_source(NSOURCES))
- allocate(Mxx(NSOURCES))
- allocate(Myy(NSOURCES))
- allocate(Mzz(NSOURCES))
- allocate(Mxy(NSOURCES))
- allocate(Mxz(NSOURCES))
- allocate(Myz(NSOURCES))
- allocate(xi_source(NSOURCES))
- allocate(eta_source(NSOURCES))
- allocate(gamma_source(NSOURCES))
- allocate(tshift_cmt(NSOURCES))
- allocate(hdur(NSOURCES))
- allocate(hdur_gaussian(NSOURCES))
- allocate(utm_x_source(NSOURCES))
- allocate(utm_y_source(NSOURCES))
- allocate(nu_source(3,3,NSOURCES))
+ allocate(islice_selected_source(NSOURCES), &
+ ispec_selected_source(NSOURCES), &
+ Mxx(NSOURCES), &
+ Myy(NSOURCES), &
+ Mzz(NSOURCES), &
+ Mxy(NSOURCES), &
+ Mxz(NSOURCES), &
+ Myz(NSOURCES), &
+ xi_source(NSOURCES), &
+ eta_source(NSOURCES), &
+ gamma_source(NSOURCES), &
+ tshift_cmt(NSOURCES), &
+ hdur(NSOURCES), &
+ hdur_gaussian(NSOURCES), &
+ utm_x_source(NSOURCES), &
+ utm_y_source(NSOURCES), &
+ nu_source(3,3,NSOURCES),stat=ier)
+ if( ier /= 0 ) stop 'error allocating arrays for sources'
! locate sources in the mesh
!
@@ -327,7 +328,7 @@
use specfem_par_acoustic
implicit none
- integer :: irec,isource !,ios
+ integer :: irec,isource,ier !,ios
! reads in station file
if (SIMULATION_TYPE == 1) then
@@ -373,14 +374,15 @@
if(nrec < 1) call exit_MPI(myrank,'need at least one receiver')
! allocate memory for receiver arrays
- allocate(islice_selected_rec(nrec))
- allocate(ispec_selected_rec(nrec))
- allocate(xi_receiver(nrec))
- allocate(eta_receiver(nrec))
- allocate(gamma_receiver(nrec))
- allocate(station_name(nrec))
- allocate(network_name(nrec))
- allocate(nu(NDIM,NDIM,nrec))
+ allocate(islice_selected_rec(nrec), &
+ ispec_selected_rec(nrec), &
+ xi_receiver(nrec), &
+ eta_receiver(nrec), &
+ gamma_receiver(nrec), &
+ station_name(nrec), &
+ network_name(nrec), &
+ nu(NDIM,NDIM,nrec),stat=ier)
+ if( ier /= 0 ) stop 'error allocating arrays for receivers'
! locate receivers in the mesh
call locate_receivers(ibool,myrank,NSPEC_AB,NGLOB_AB, &
@@ -518,8 +520,9 @@
! forward simulations
if (SIMULATION_TYPE == 1 .or. SIMULATION_TYPE == 3) then
- allocate(sourcearray(NDIM,NGLLX,NGLLY,NGLLZ))
- allocate(sourcearrays(NSOURCES,NDIM,NGLLX,NGLLY,NGLLZ))
+ allocate(sourcearray(NDIM,NGLLX,NGLLY,NGLLZ), &
+ sourcearrays(NSOURCES,NDIM,NGLLX,NGLLY,NGLLZ),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array sourcearray'
! compute source arrays
do isource = 1,NSOURCES
@@ -625,13 +628,15 @@
! because we may need to read in adjoint sources block by block
! initializes adjoint sources
- allocate(adj_sourcearrays(nadj_rec_local,NTSTEP_BETWEEN_READ_ADJSRC,NDIM,NGLLX,NGLLY,NGLLZ))
+ allocate(adj_sourcearrays(nadj_rec_local,NTSTEP_BETWEEN_READ_ADJSRC,NDIM,NGLLX,NGLLY,NGLLZ),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array adj_sourcearrays'
adj_sourcearrays = 0._CUSTOM_REAL
else
! allocate dummy array in order to be able to use it as a subroutine argument, even if unused
- allocate(adj_sourcearrays(1,1,1,1,1,1))
+ allocate(adj_sourcearrays(1,1,1,1,1,1),stat=ier)
+ if( ier /= 0 ) stop 'error allocating dummy array adj_sourcearrays'
endif
@@ -647,17 +652,19 @@
use specfem_par
implicit none
- integer :: irec,irec_local,isource
+ integer :: irec,irec_local,isource,ier
! stores local receivers interpolation factors
if (nrec_local > 0) then
! allocate Lagrange interpolators for receivers
- allocate(hxir_store(nrec_local,NGLLX))
- allocate(hetar_store(nrec_local,NGLLY))
- allocate(hgammar_store(nrec_local,NGLLZ))
+ allocate(hxir_store(nrec_local,NGLLX), &
+ hetar_store(nrec_local,NGLLY), &
+ hgammar_store(nrec_local,NGLLZ),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array hxir_store etc.'
! define local to global receiver numbering mapping
- allocate(number_receiver_global(nrec_local))
+ allocate(number_receiver_global(nrec_local),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array number_reciever_global'
irec_local = 0
if (SIMULATION_TYPE == 1 .or. SIMULATION_TYPE == 3) then
do irec = 1,nrec
@@ -687,9 +694,10 @@
hgammar_store(irec_local,:) = hgammar(:)
enddo
else
- allocate(hpxir_store(nrec_local,NGLLX))
- allocate(hpetar_store(nrec_local,NGLLY))
- allocate(hpgammar_store(nrec_local,NGLLZ))
+ allocate(hpxir_store(nrec_local,NGLLX), &
+ hpetar_store(nrec_local,NGLLY), &
+ hpgammar_store(nrec_local,NGLLZ),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array hpxir_store'
do irec_local = 1,nrec_local
irec = number_receiver_global(irec_local)
call lagrange_any(xi_source(irec),NGLLX,xigll,hxir,hpxir)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -426,8 +426,8 @@
real(kind=CUSTOM_REAL), dimension(:,:,:,:),allocatable:: div, curl_x, curl_y, curl_z
real(kind=CUSTOM_REAL), dimension(:,:,:,:),allocatable:: velocity_x,velocity_y,velocity_z
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: dvxdxl,dvxdyl,&
- dvxdzl,dvydxl,dvydyl,dvydzl,dvzdxl,dvzdyl,dvzdzl
+! real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: dvxdxl,dvxdyl,&
+! dvxdzl,dvydxl,dvydyl,dvydzl,dvzdxl,dvzdyl,dvzdzl
! shakemovies
real(kind=CUSTOM_REAL), dimension(:), allocatable :: store_val_x_external_mesh
@@ -457,8 +457,8 @@
integer,dimension(:),allocatable :: faces_surface_offset_ext_mesh
integer,dimension(:,:),allocatable :: faces_surface_ext_mesh
integer,dimension(:),allocatable :: faces_surface_ext_mesh_ispec
- integer :: nfaces_surface_ext_mesh
- integer :: nfaces_surface_glob_ext_mesh
+ integer :: nfaces_surface_ext_mesh,nfaces_surface_ext_mesh_points
+ integer :: nfaces_surface_glob_ext_mesh,nfaces_surface_glob_em_points
! face corner indices
integer :: iorderi(NGNOD2D),iorderj(NGNOD2D)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_movie_output.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_movie_output.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_movie_output.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -78,12 +78,20 @@
use specfem_par_movie
implicit none
- real(kind=CUSTOM_REAL),dimension(NDIM,NGLLX,NGLLY,NGLLZ):: &
+ real(kind=CUSTOM_REAL),dimension(:,:,:,:),allocatable:: &
displ_element,veloc_element,accel_element
+ real(kind=CUSTOM_REAL),dimension(1):: dummy
integer :: ipoin,ispec,iglob,ispec2D
- integer :: i,j,k
+ integer :: i,j,k,ier
logical :: is_done
+ ! allocate array for single elements
+ allocate( displ_element(NDIM,NGLLX,NGLLY,NGLLZ), &
+ veloc_element(NDIM,NGLLX,NGLLY,NGLLZ), &
+ accel_element(NDIM,NGLLX,NGLLY,NGLLZ),stat=ier)
+ if( ier /= 0 ) stop 'error allocating arrays for movie elements'
+
+
! initializes arrays for point coordinates
if (it == 1) then
store_val_ux_external_mesh(:) = -HUGEVAL
@@ -254,49 +262,52 @@
! finalizes shakemap: master process collects all info
if (it == NSTEP) then
- if (USE_HIGHRES_FOR_MOVIES) then
- call gatherv_all_cr(store_val_x_external_mesh,nfaces_surface_ext_mesh*NGLLX*NGLLY,&
- store_val_x_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
- call gatherv_all_cr(store_val_y_external_mesh,nfaces_surface_ext_mesh*NGLLX*NGLLY,&
- store_val_y_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
- call gatherv_all_cr(store_val_z_external_mesh,nfaces_surface_ext_mesh*NGLLX*NGLLY,&
- store_val_z_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
- call gatherv_all_cr(store_val_ux_external_mesh,nfaces_surface_ext_mesh*NGLLX*NGLLY,&
- store_val_ux_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
- call gatherv_all_cr(store_val_uy_external_mesh,nfaces_surface_ext_mesh*NGLLX*NGLLY,&
- store_val_uy_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
- call gatherv_all_cr(store_val_uz_external_mesh,nfaces_surface_ext_mesh*NGLLX*NGLLY,&
- store_val_uz_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
+ ! master collects data
+ if( myrank == 0 ) then
+ call gatherv_all_cr(store_val_x_external_mesh,nfaces_surface_ext_mesh_points,&
+ store_val_x_all_external_mesh,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_em_points,NPROC)
+ call gatherv_all_cr(store_val_y_external_mesh,nfaces_surface_ext_mesh_points,&
+ store_val_y_all_external_mesh,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_em_points,NPROC)
+ call gatherv_all_cr(store_val_z_external_mesh,nfaces_surface_ext_mesh_points,&
+ store_val_z_all_external_mesh,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_em_points,NPROC)
+ call gatherv_all_cr(store_val_ux_external_mesh,nfaces_surface_ext_mesh_points,&
+ store_val_ux_all_external_mesh,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_em_points,NPROC)
+ call gatherv_all_cr(store_val_uy_external_mesh,nfaces_surface_ext_mesh_points,&
+ store_val_uy_all_external_mesh,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_em_points,NPROC)
+ call gatherv_all_cr(store_val_uz_external_mesh,nfaces_surface_ext_mesh_points,&
+ store_val_uz_all_external_mesh,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_em_points,NPROC)
else
- call gatherv_all_cr(store_val_x_external_mesh,nfaces_surface_ext_mesh*NGNOD2D,&
- store_val_x_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
- call gatherv_all_cr(store_val_y_external_mesh,nfaces_surface_ext_mesh*NGNOD2D,&
- store_val_y_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
- call gatherv_all_cr(store_val_z_external_mesh,nfaces_surface_ext_mesh*NGNOD2D,&
- store_val_z_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
- call gatherv_all_cr(store_val_ux_external_mesh,nfaces_surface_ext_mesh*NGNOD2D,&
- store_val_ux_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
- call gatherv_all_cr(store_val_uy_external_mesh,nfaces_surface_ext_mesh*NGNOD2D,&
- store_val_uy_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
- call gatherv_all_cr(store_val_uz_external_mesh,nfaces_surface_ext_mesh*NGNOD2D,&
- store_val_uz_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
+ ! all other process just send
+ call gatherv_all_cr(store_val_x_external_mesh,nfaces_surface_ext_mesh_points,&
+ dummy,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ 1,NPROC)
+ call gatherv_all_cr(store_val_y_external_mesh,nfaces_surface_ext_mesh_points,&
+ dummy,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ 1,NPROC)
+ call gatherv_all_cr(store_val_z_external_mesh,nfaces_surface_ext_mesh_points,&
+ dummy,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ 1,NPROC)
+ call gatherv_all_cr(store_val_ux_external_mesh,nfaces_surface_ext_mesh_points,&
+ dummy,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ 1,NPROC)
+ call gatherv_all_cr(store_val_uy_external_mesh,nfaces_surface_ext_mesh_points,&
+ dummy,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ 1,NPROC)
+ call gatherv_all_cr(store_val_uz_external_mesh,nfaces_surface_ext_mesh_points,&
+ dummy,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ 1,NPROC)
endif
! creates shakemap file
if(myrank == 0) then
- open(unit=IOUT,file=trim(OUTPUT_FILES)//'/shakingdata',status='unknown',form='unformatted')
+ open(unit=IOUT,file=trim(OUTPUT_FILES)//'/shakingdata',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file shakingdata'
write(IOUT) store_val_x_all_external_mesh ! x coordinates
write(IOUT) store_val_y_all_external_mesh ! y coordinates
write(IOUT) store_val_z_all_external_mesh ! z coordinates
@@ -307,6 +318,8 @@
endif
endif
+ deallocate(displ_element,veloc_element,accel_element)
+
end subroutine wmo_create_shakemap_em
@@ -322,10 +335,15 @@
use specfem_par_movie
implicit none
- real(kind=CUSTOM_REAL),dimension(NDIM,NGLLX,NGLLY,NGLLZ):: veloc_element
- integer :: ispec2D,ispec,ipoin,iglob,i,j,k
+ real(kind=CUSTOM_REAL),dimension(:,:,:,:),allocatable:: veloc_element
+ real(kind=CUSTOM_REAL),dimension(1):: dummy
+ integer :: ispec2D,ispec,ipoin,iglob,i,j,k,ier
logical :: is_done
+ ! allocate array for single elements
+ allocate( veloc_element(NDIM,NGLLX,NGLLY,NGLLZ),stat=ier)
+ if( ier /= 0 ) stop 'error allocating arrays for movie elements'
+
! initializes arrays for point coordinates
if (it == NTSTEP_BETWEEN_FRAMES ) then
do ispec2D = 1,nfaces_surface_ext_mesh
@@ -430,58 +448,62 @@
enddo
! master process collects all info
- if (USE_HIGHRES_FOR_MOVIES) then
- ! collects locations only once
- if (it == NTSTEP_BETWEEN_FRAMES ) then
- call gatherv_all_cr(store_val_x_external_mesh,nfaces_surface_ext_mesh*NGLLX*NGLLY,&
- store_val_x_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
- call gatherv_all_cr(store_val_y_external_mesh,nfaces_surface_ext_mesh*NGLLX*NGLLY,&
- store_val_y_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
- call gatherv_all_cr(store_val_z_external_mesh,nfaces_surface_ext_mesh*NGLLX*NGLLY,&
- store_val_z_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
+ ! collects locations only once
+ if (it == NTSTEP_BETWEEN_FRAMES ) then
+ ! master collects all
+ if( myrank == 0 ) then
+ call gatherv_all_cr(store_val_x_external_mesh,nfaces_surface_ext_mesh_points,&
+ store_val_x_all_external_mesh,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_em_points,NPROC)
+ call gatherv_all_cr(store_val_y_external_mesh,nfaces_surface_ext_mesh_points,&
+ store_val_y_all_external_mesh,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_em_points,NPROC)
+ call gatherv_all_cr(store_val_z_external_mesh,nfaces_surface_ext_mesh_points,&
+ store_val_z_all_external_mesh,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_em_points,NPROC)
+ else
+ ! slaves just send
+ call gatherv_all_cr(store_val_x_external_mesh,nfaces_surface_ext_mesh_points,&
+ dummy,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ 1,NPROC)
+ call gatherv_all_cr(store_val_y_external_mesh,nfaces_surface_ext_mesh_points,&
+ dummy,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ 1,NPROC)
+ call gatherv_all_cr(store_val_z_external_mesh,nfaces_surface_ext_mesh_points,&
+ dummy,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ 1,NPROC)
endif
- ! updates/gathers velocity field (high-res)
- call gatherv_all_cr(store_val_ux_external_mesh,nfaces_surface_ext_mesh*NGLLX*NGLLY,&
- store_val_ux_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
- call gatherv_all_cr(store_val_uy_external_mesh,nfaces_surface_ext_mesh*NGLLX*NGLLY,&
- store_val_uy_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
- call gatherv_all_cr(store_val_uz_external_mesh,nfaces_surface_ext_mesh*NGLLX*NGLLY,&
- store_val_uz_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
+ endif
+
+ ! updates/gathers velocity field (high-res or low-res)
+ if( myrank == 0 ) then
+ call gatherv_all_cr(store_val_ux_external_mesh,nfaces_surface_ext_mesh_points,&
+ store_val_ux_all_external_mesh,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_em_points,NPROC)
+ call gatherv_all_cr(store_val_uy_external_mesh,nfaces_surface_ext_mesh_points,&
+ store_val_uy_all_external_mesh,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_em_points,NPROC)
+ call gatherv_all_cr(store_val_uz_external_mesh,nfaces_surface_ext_mesh_points,&
+ store_val_uz_all_external_mesh,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_em_points,NPROC)
else
- ! collects locations only once
- if (it == NTSTEP_BETWEEN_FRAMES ) then
- call gatherv_all_cr(store_val_x_external_mesh,nfaces_surface_ext_mesh*NGNOD2D,&
- store_val_x_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
- call gatherv_all_cr(store_val_y_external_mesh,nfaces_surface_ext_mesh*NGNOD2D,&
- store_val_y_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
- call gatherv_all_cr(store_val_z_external_mesh,nfaces_surface_ext_mesh*NGNOD2D,&
- store_val_z_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
- endif
- ! updates/gathers velocity field (low-res)
- call gatherv_all_cr(store_val_ux_external_mesh,nfaces_surface_ext_mesh*NGNOD2D,&
- store_val_ux_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
- call gatherv_all_cr(store_val_uy_external_mesh,nfaces_surface_ext_mesh*NGNOD2D,&
- store_val_uy_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
- call gatherv_all_cr(store_val_uz_external_mesh,nfaces_surface_ext_mesh*NGNOD2D,&
- store_val_uz_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
+ !slaves
+ call gatherv_all_cr(store_val_ux_external_mesh,nfaces_surface_ext_mesh_points,&
+ dummy,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ 1,NPROC)
+ call gatherv_all_cr(store_val_uy_external_mesh,nfaces_surface_ext_mesh_points,&
+ dummy,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ 1,NPROC)
+ call gatherv_all_cr(store_val_uz_external_mesh,nfaces_surface_ext_mesh_points,&
+ dummy,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ 1,NPROC)
endif
! file output
if(myrank == 0) then
write(outputname,"('/moviedata',i6.6)") it
- open(unit=IOUT,file=trim(OUTPUT_FILES)//outputname,status='unknown',form='unformatted')
+ open(unit=IOUT,file=trim(OUTPUT_FILES)//outputname,status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file moviedata'
write(IOUT) store_val_x_all_external_mesh ! x coordinate
write(IOUT) store_val_y_all_external_mesh ! y coordinate
write(IOUT) store_val_z_all_external_mesh ! z coordinate
@@ -491,6 +513,8 @@
close(IOUT)
endif
+ deallocate(veloc_element)
+
end subroutine wmo_create_movie_surface_em
@@ -506,11 +530,16 @@
use specfem_par_movie
implicit none
- real(kind=CUSTOM_REAL),dimension(NDIM,NGLLX,NGLLY,NGLLZ):: val_element
- integer :: ispec,ipoin,iglob,i,j,k
+ real(kind=CUSTOM_REAL),dimension(:,:,:,:),allocatable:: val_element
+ real(kind=CUSTOM_REAL),dimension(1) :: dummy
+ integer :: ispec,ipoin,iglob,i,j,k,ier
integer :: imin,imax,jmin,jmax,kmin,kmax,iface,igll,iloc
logical :: is_done
+ ! allocate array for single elements
+ allocate( val_element(NDIM,NGLLX,NGLLY,NGLLZ),stat=ier)
+ if( ier /= 0 ) stop 'error allocating arrays for movie elements'
+
! initializes arrays for point coordinates
if (it == NTSTEP_BETWEEN_FRAMES ) then
ipoin = 0
@@ -680,54 +709,59 @@
enddo ! iface
! master process collects all info
- if (USE_HIGHRES_FOR_MOVIES) then
- if (it == NTSTEP_BETWEEN_FRAMES ) then
- call gatherv_all_cr(store_val_x_external_mesh,nfaces_surface_ext_mesh*NGLLX*NGLLY,&
- store_val_x_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
- call gatherv_all_cr(store_val_y_external_mesh,nfaces_surface_ext_mesh*NGLLX*NGLLY,&
- store_val_y_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
- call gatherv_all_cr(store_val_z_external_mesh,nfaces_surface_ext_mesh*NGLLX*NGLLY,&
- store_val_z_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
+ if (it == NTSTEP_BETWEEN_FRAMES ) then
+ ! master collects
+ if( myrank == 0 ) then
+ call gatherv_all_cr(store_val_x_external_mesh,nfaces_surface_ext_mesh_points,&
+ store_val_x_all_external_mesh,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_em_points,NPROC)
+ call gatherv_all_cr(store_val_y_external_mesh,nfaces_surface_ext_mesh_points,&
+ store_val_y_all_external_mesh,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_em_points,NPROC)
+ call gatherv_all_cr(store_val_z_external_mesh,nfaces_surface_ext_mesh_points,&
+ store_val_z_all_external_mesh,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_em_points,NPROC)
+ else
+ call gatherv_all_cr(store_val_x_external_mesh,nfaces_surface_ext_mesh_points,&
+ dummy,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ 1,NPROC)
+ call gatherv_all_cr(store_val_y_external_mesh,nfaces_surface_ext_mesh_points,&
+ dummy,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ 1,NPROC)
+ call gatherv_all_cr(store_val_z_external_mesh,nfaces_surface_ext_mesh_points,&
+ dummy,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ 1,NPROC)
endif
- call gatherv_all_cr(store_val_ux_external_mesh,nfaces_surface_ext_mesh*NGLLX*NGLLY,&
- store_val_ux_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
- call gatherv_all_cr(store_val_uy_external_mesh,nfaces_surface_ext_mesh*NGLLX*NGLLY,&
- store_val_uy_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
- call gatherv_all_cr(store_val_uz_external_mesh,nfaces_surface_ext_mesh*NGLLX*NGLLY,&
- store_val_uz_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
+ endif
+
+ ! master collects wavefield
+ if( myrank == 0 ) then
+ call gatherv_all_cr(store_val_ux_external_mesh,nfaces_surface_ext_mesh_points,&
+ store_val_ux_all_external_mesh,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_em_points,NPROC)
+ call gatherv_all_cr(store_val_uy_external_mesh,nfaces_surface_ext_mesh_points,&
+ store_val_uy_all_external_mesh,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_em_points,NPROC)
+ call gatherv_all_cr(store_val_uz_external_mesh,nfaces_surface_ext_mesh_points,&
+ store_val_uz_all_external_mesh,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_em_points,NPROC)
else
- if (it == NTSTEP_BETWEEN_FRAMES ) then
- call gatherv_all_cr(store_val_x_external_mesh,nfaces_surface_ext_mesh*NGNOD2D,&
- store_val_x_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
- call gatherv_all_cr(store_val_y_external_mesh,nfaces_surface_ext_mesh*NGNOD2D,&
- store_val_y_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
- call gatherv_all_cr(store_val_z_external_mesh,nfaces_surface_ext_mesh*NGNOD2D,&
- store_val_z_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
- endif
- call gatherv_all_cr(store_val_ux_external_mesh,nfaces_surface_ext_mesh*NGNOD2D,&
- store_val_ux_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
- call gatherv_all_cr(store_val_uy_external_mesh,nfaces_surface_ext_mesh*NGNOD2D,&
- store_val_uy_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
- call gatherv_all_cr(store_val_uz_external_mesh,nfaces_surface_ext_mesh*NGNOD2D,&
- store_val_uz_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
+ call gatherv_all_cr(store_val_ux_external_mesh,nfaces_surface_ext_mesh_points,&
+ dummy,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ 1,NPROC)
+ call gatherv_all_cr(store_val_uy_external_mesh,nfaces_surface_ext_mesh_points,&
+ dummy,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ 1,NPROC)
+ call gatherv_all_cr(store_val_uz_external_mesh,nfaces_surface_ext_mesh_points,&
+ dummy,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ 1,NPROC)
endif
! file output: note that values are only stored on free surface
if(myrank == 0) then
write(outputname,"('/moviedata',i6.6)") it
- open(unit=IOUT,file=trim(OUTPUT_FILES)//outputname,status='unknown',form='unformatted')
+ open(unit=IOUT,file=trim(OUTPUT_FILES)//outputname,status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file moviedata'
write(IOUT) store_val_x_all_external_mesh ! x coordinate
write(IOUT) store_val_y_all_external_mesh ! y coordinate
write(IOUT) store_val_z_all_external_mesh ! z coordinate
@@ -737,6 +771,8 @@
close(IOUT)
endif
+ deallocate(val_element)
+
end subroutine wmo_movie_surface_output_o
@@ -752,13 +788,20 @@
use specfem_par_movie
implicit none
- real(kind=CUSTOM_REAL),dimension(NDIM,NGLLX,NGLLY,NGLLZ):: &
+ real(kind=CUSTOM_REAL),dimension(:,:,:,:),allocatable:: &
displ_element,veloc_element,accel_element
+ real(kind=CUSTOM_REAL),dimension(1):: dummy
integer :: ipoin,ispec,iglob
integer :: imin,imax,jmin,jmax,kmin,kmax,iface,igll,iloc
- integer :: i,j,k
+ integer :: i,j,k,ier
logical :: is_done
+ ! allocate array for single elements
+ allocate( displ_element(NDIM,NGLLX,NGLLY,NGLLZ), &
+ veloc_element(NDIM,NGLLX,NGLLY,NGLLZ), &
+ accel_element(NDIM,NGLLX,NGLLY,NGLLZ),stat=ier)
+ if( ier /= 0 ) stop 'error allocating arrays for movie elements'
+
! outputs values on free surface
ipoin = 0
do iface=1,num_free_surface_faces
@@ -904,49 +947,51 @@
! saves shakemap only at the end of the simulation
if(it == NSTEP) then
- if (USE_HIGHRES_FOR_MOVIES) then
- call gatherv_all_cr(store_val_x_external_mesh,nfaces_surface_ext_mesh*NGLLX*NGLLY,&
- store_val_x_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
- call gatherv_all_cr(store_val_y_external_mesh,nfaces_surface_ext_mesh*NGLLX*NGLLY,&
- store_val_y_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
- call gatherv_all_cr(store_val_z_external_mesh,nfaces_surface_ext_mesh*NGLLX*NGLLY,&
- store_val_z_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
- call gatherv_all_cr(store_val_ux_external_mesh,nfaces_surface_ext_mesh*NGLLX*NGLLY,&
- store_val_ux_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
- call gatherv_all_cr(store_val_uy_external_mesh,nfaces_surface_ext_mesh*NGLLX*NGLLY,&
- store_val_uy_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
- call gatherv_all_cr(store_val_uz_external_mesh,nfaces_surface_ext_mesh*NGLLX*NGLLY,&
- store_val_uz_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGLLX*NGLLY,NPROC)
+ ! master collects
+ if( myrank == 0 ) then
+ call gatherv_all_cr(store_val_x_external_mesh,nfaces_surface_ext_mesh_points,&
+ store_val_x_all_external_mesh,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_em_points,NPROC)
+ call gatherv_all_cr(store_val_y_external_mesh,nfaces_surface_ext_mesh_points,&
+ store_val_y_all_external_mesh,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_em_points,NPROC)
+ call gatherv_all_cr(store_val_z_external_mesh,nfaces_surface_ext_mesh_points,&
+ store_val_z_all_external_mesh,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_em_points,NPROC)
+ call gatherv_all_cr(store_val_ux_external_mesh,nfaces_surface_ext_mesh_points,&
+ store_val_ux_all_external_mesh,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_em_points,NPROC)
+ call gatherv_all_cr(store_val_uy_external_mesh,nfaces_surface_ext_mesh_points,&
+ store_val_uy_all_external_mesh,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_em_points,NPROC)
+ call gatherv_all_cr(store_val_uz_external_mesh,nfaces_surface_ext_mesh_points,&
+ store_val_uz_all_external_mesh,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ nfaces_surface_glob_em_points,NPROC)
else
- call gatherv_all_cr(store_val_x_external_mesh,nfaces_surface_ext_mesh*NGNOD2D,&
- store_val_x_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
- call gatherv_all_cr(store_val_y_external_mesh,nfaces_surface_ext_mesh*NGNOD2D,&
- store_val_y_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
- call gatherv_all_cr(store_val_z_external_mesh,nfaces_surface_ext_mesh*NGNOD2D,&
- store_val_z_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
- call gatherv_all_cr(store_val_ux_external_mesh,nfaces_surface_ext_mesh*NGNOD2D,&
- store_val_ux_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
- call gatherv_all_cr(store_val_uy_external_mesh,nfaces_surface_ext_mesh*NGNOD2D,&
- store_val_uy_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
- call gatherv_all_cr(store_val_uz_external_mesh,nfaces_surface_ext_mesh*NGNOD2D,&
- store_val_uz_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGNOD2D,faces_surface_offset_ext_mesh,&
- nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
+ call gatherv_all_cr(store_val_x_external_mesh,nfaces_surface_ext_mesh_points,&
+ dummy,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ 1,NPROC)
+ call gatherv_all_cr(store_val_y_external_mesh,nfaces_surface_ext_mesh_points,&
+ dummy,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ 1,NPROC)
+ call gatherv_all_cr(store_val_z_external_mesh,nfaces_surface_ext_mesh_points,&
+ dummy,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ 1,NPROC)
+ call gatherv_all_cr(store_val_ux_external_mesh,nfaces_surface_ext_mesh_points,&
+ dummy,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ 1,NPROC)
+ call gatherv_all_cr(store_val_uy_external_mesh,nfaces_surface_ext_mesh_points,&
+ dummy,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ 1,NPROC)
+ call gatherv_all_cr(store_val_uz_external_mesh,nfaces_surface_ext_mesh_points,&
+ dummy,nfaces_perproc_surface_ext_mesh,faces_surface_offset_ext_mesh,&
+ 1,NPROC)
endif
! creates shakemap file: note that values are only stored on free surface
if(myrank == 0) then
- open(unit=IOUT,file=trim(OUTPUT_FILES)//'/shakingdata',status='unknown',form='unformatted')
+ open(unit=IOUT,file=trim(OUTPUT_FILES)//'/shakingdata',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file shakingdata'
write(IOUT) store_val_x_all_external_mesh ! x coordinates
write(IOUT) store_val_y_all_external_mesh ! y coordinates
write(IOUT) store_val_z_all_external_mesh ! z coordinates
@@ -958,6 +1003,8 @@
endif ! NTSTEP
+ deallocate(displ_element,veloc_element,accel_element)
+
end subroutine wmo_create_shakemap_o
@@ -973,13 +1020,14 @@
use specfem_par_movie
implicit none
- real(kind=CUSTOM_REAL),dimension(NDIM,NGLLX,NGLLY,NGLLZ):: veloc_element
- real(kind=CUSTOM_REAL),dimension(NGLOB_AB):: div_glob,curl_glob ! divergence and curl only in the global nodes
- integer :: ispec,i,j,k,l,iglob
- integer,dimension(NGLOB_AB) :: valency
+ real(kind=CUSTOM_REAL),dimension(:,:,:,:),allocatable:: veloc_element
+ ! divergence and curl only in the global nodes
+ real(kind=CUSTOM_REAL),dimension(:),allocatable:: div_glob,curl_glob
+ integer,dimension(:),allocatable :: valency
+ integer :: ispec,ier
character(len=3) :: channel
character(len=1) :: compx,compy,compz
-
+
! gets component characters: X/Y/Z or E/N/Z
call write_channel_name(1,channel)
compx(1:1) = channel(3:3) ! either X or E
@@ -994,6 +1042,11 @@
velocity_z(:,:,:,:) = 0._CUSTOM_REAL
if( ACOUSTIC_SIMULATION ) then
+
+ ! allocate array for single elements
+ allocate( veloc_element(NDIM,NGLLX,NGLLY,NGLLZ),stat=ier)
+ if( ier /= 0 ) stop 'error allocating arrays for movie elements'
+
! uses div as temporary array to store velocity on all gll points
do ispec=1,NSPEC_AB
if( .not. ispec_is_acoustic(ispec) ) cycle
@@ -1008,130 +1061,64 @@
velocity_y(:,:,:,ispec) = veloc_element(2,:,:,:)
velocity_z(:,:,:,ispec) = veloc_element(3,:,:,:)
enddo
+
+ deallocate(veloc_element)
+
endif ! acoustic
! saves full snapshot data to local disk
if( ELASTIC_SIMULATION ) then
- div_glob=0.0_CUSTOM_REAL
- curl_glob=0.0_CUSTOM_REAL
- do ispec=1,NSPEC_AB
- if( .not. ispec_is_elastic(ispec) ) cycle
-
- ! calculates divergence and curl of velocity field
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- tempx1l = 0._CUSTOM_REAL
- tempx2l = 0._CUSTOM_REAL
- tempx3l = 0._CUSTOM_REAL
- tempy1l = 0._CUSTOM_REAL
- tempy2l = 0._CUSTOM_REAL
- tempy3l = 0._CUSTOM_REAL
- tempz1l = 0._CUSTOM_REAL
- tempz2l = 0._CUSTOM_REAL
- tempz3l = 0._CUSTOM_REAL
-
- do l=1,NGLLX
- hp1 = hprime_xx(i,l)
- iglob = ibool(l,j,k,ispec)
- tempx1l = tempx1l + veloc(1,iglob)*hp1
- tempy1l = tempy1l + veloc(2,iglob)*hp1
- tempz1l = tempz1l + veloc(3,iglob)*hp1
- hp2 = hprime_yy(j,l)
- iglob = ibool(i,l,k,ispec)
- tempx2l = tempx2l + veloc(1,iglob)*hp2
- tempy2l = tempy2l + veloc(2,iglob)*hp2
- tempz2l = tempz2l + veloc(3,iglob)*hp2
- hp3 = hprime_zz(k,l)
- iglob = ibool(i,j,l,ispec)
- tempx3l = tempx3l + veloc(1,iglob)*hp3
- tempy3l = tempy3l + veloc(2,iglob)*hp3
- tempz3l = tempz3l + veloc(3,iglob)*hp3
- enddo
-
- ! get derivatives of ux, uy and uz with respect to x, y and z
- xixl = xix(i,j,k,ispec)
- xiyl = xiy(i,j,k,ispec)
- xizl = xiz(i,j,k,ispec)
- etaxl = etax(i,j,k,ispec)
- etayl = etay(i,j,k,ispec)
- etazl = etaz(i,j,k,ispec)
- gammaxl = gammax(i,j,k,ispec)
- gammayl = gammay(i,j,k,ispec)
- gammazl = gammaz(i,j,k,ispec)
-
- dvxdxl(i,j,k) = xixl*tempx1l + etaxl*tempx2l + gammaxl*tempx3l
- dvxdyl(i,j,k) = xiyl*tempx1l + etayl*tempx2l + gammayl*tempx3l
- dvxdzl(i,j,k) = xizl*tempx1l + etazl*tempx2l + gammazl*tempx3l
-
- dvydxl(i,j,k) = xixl*tempy1l + etaxl*tempy2l + gammaxl*tempy3l
- dvydyl(i,j,k) = xiyl*tempy1l + etayl*tempy2l + gammayl*tempy3l
- dvydzl(i,j,k) = xizl*tempy1l + etazl*tempy2l + gammazl*tempy3l
-
- dvzdxl(i,j,k) = xixl*tempz1l + etaxl*tempz2l + gammaxl*tempz3l
- dvzdyl(i,j,k) = xiyl*tempz1l + etayl*tempz2l + gammayl*tempz3l
- dvzdzl(i,j,k) = xizl*tempz1l + etazl*tempz2l + gammazl*tempz3l
-
- enddo
- enddo
- enddo
-
- do k = 1,NGLLZ
- do j = 1,NGLLY
- do i = 1,NGLLX
- ! divergence \nabla \cdot \bf{v}
- div(i,j,k,ispec) = dvxdxl(i,j,k) + dvydyl(i,j,k) + dvzdzl(i,j,k)
- ! curl
- curl_x(i,j,k,ispec) = dvzdyl(i,j,k) - dvydzl(i,j,k)
- curl_y(i,j,k,ispec) = dvxdzl(i,j,k) - dvzdxl(i,j,k)
- curl_z(i,j,k,ispec) = dvydxl(i,j,k) - dvxdyl(i,j,k)
- ! velocity field
- iglob = ibool(i,j,k,ispec)
- velocity_x(i,j,k,ispec) = veloc(1,iglob)
- velocity_y(i,j,k,ispec) = veloc(2,iglob)
- velocity_z(i,j,k,ispec) = veloc(3,iglob)
-
- valency(iglob)=valency(iglob)+1
-
- div_glob(iglob) = div_glob(iglob) + div(i,j,k,ispec)
- curl_glob(iglob)=curl_glob(iglob)+0.5_CUSTOM_REAL*(curl_x(i,j,k,ispec)+curl_x(i,j,k,ispec)+curl_x(i,j,k,ispec))
- enddo
- enddo
- enddo
- enddo !NSPEC_AB
-
- do i=1,NGLOB_AB
- div_glob(i)=div_glob(i)/valency(i)
- curl_glob(i)=curl_glob(i)/valency(i)
- enddo
-
+ ! allocate array for single elements
+ allocate(div_glob(NGLOB_AB), &
+ curl_glob(NGLOB_AB), &
+ valency(NGLOB_AB), stat=ier)
+ if( ier /= 0 ) stop 'error allocating arrays for movie div and curl'
+
+ ! calculates divergence and curl of velocity field
+ call wmo_movie_div_curl(NSPEC_AB,NGLOB_AB,veloc, &
+ div_glob,curl_glob,valency, &
+ div,curl_x,curl_y,curl_z, &
+ velocity_x,velocity_y,velocity_z, &
+ ibool,ispec_is_elastic, &
+ hprime_xx,hprime_yy,hprime_zz, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz)
+
+ ! writes out div and curl on global points
write(outputname,"('/proc',i6.6,'_div_glob_it',i6.6,'.bin')") myrank,it
- open(unit=27,file=trim(LOCAL_PATH)//trim(outputname),status='unknown',form='unformatted')
+ open(unit=27,file=trim(LOCAL_PATH)//trim(outputname),status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file div_glob'
write(27) div_glob
close(27)
write(outputname,"('/proc',i6.6,'_curl_glob_it',i6.6,'.bin')") myrank,it
- open(unit=27,file=trim(LOCAL_PATH)//trim(outputname),status='unknown',form='unformatted')
+ open(unit=27,file=trim(LOCAL_PATH)//trim(outputname),status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file curl_glob'
write(27) curl_glob
close(27)
+ deallocate(div_glob,curl_glob,valency)
+ ! writes our divergence
write(outputname,"('/proc',i6.6,'_div_it',i6.6,'.bin')") myrank,it
- open(unit=27,file=trim(LOCAL_PATH)//trim(outputname),status='unknown',form='unformatted')
+ open(unit=27,file=trim(LOCAL_PATH)//trim(outputname),status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file div_it'
write(27) div
close(27)
! writes out curl
write(outputname,"('/proc',i6.6,'_curl_',a1,'_it',i6.6,'.bin')") myrank,compx,it
- open(unit=27,file=trim(LOCAL_PATH)//trim(outputname),status='unknown',form='unformatted')
+ open(unit=27,file=trim(LOCAL_PATH)//trim(outputname),status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file curl_x_it '
write(27) curl_x
close(27)
write(outputname,"('/proc',i6.6,'_curl_',a1,'_it',i6.6,'.bin')") myrank,compy,it
- open(unit=27,file=trim(LOCAL_PATH)//trim(outputname),status='unknown',form='unformatted')
+ open(unit=27,file=trim(LOCAL_PATH)//trim(outputname),status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file curl_y_it'
write(27) curl_y
close(27)
write(outputname,"('/proc',i6.6,'_curl_',a1,'_it',i6.6,'.bin')") myrank,compz,it
- open(unit=27,file=trim(LOCAL_PATH)//trim(outputname),status='unknown',form='unformatted')
+ open(unit=27,file=trim(LOCAL_PATH)//trim(outputname),status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file curl_z_it'
write(27) curl_z
close(27)
@@ -1144,26 +1131,184 @@
if( ACOUSTIC_SIMULATION .or. ELASTIC_SIMULATION ) then
write(outputname,"('/proc',i6.6,'_velocity_',a1,'_it',i6.6,'.bin')") myrank,compx,it
- open(unit=27,file=trim(LOCAL_PATH)//trim(outputname),status='unknown',form='unformatted')
+ open(unit=27,file=trim(LOCAL_PATH)//trim(outputname),status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file movie output velocity x'
write(27) velocity_x
close(27)
write(outputname,"('/proc',i6.6,'_velocity_',a1,'_it',i6.6,'.bin')") myrank,compy,it
- open(unit=27,file=trim(LOCAL_PATH)//trim(outputname),status='unknown',form='unformatted')
+ open(unit=27,file=trim(LOCAL_PATH)//trim(outputname),status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file movie output velocity y'
write(27) velocity_y
close(27)
write(outputname,"('/proc',i6.6,'_velocity_',a1,'_it',i6.6,'.bin')") myrank,compz,it
- open(unit=27,file=trim(LOCAL_PATH)//trim(outputname),status='unknown',form='unformatted')
+ open(unit=27,file=trim(LOCAL_PATH)//trim(outputname),status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file movie output velocity z'
write(27) velocity_z
close(27)
!write(outputname,"('/proc',i6.6,'_veloc_it',i6.6,'.bin')") myrank,it
- !open(unit=27,file=trim(LOCAL_PATH)//trim(outputname),status='unknown',form='unformatted')
+ !open(unit=27,file=trim(LOCAL_PATH)//trim(outputname),status='unknown',form='unformatted',iostat=ier)
+ !if( ier /= 0 ) stop 'error opening file velocity_movie'
!write(27) velocity_movie
!close(27)
endif
-
+
end subroutine wmo_movie_volume_output
+!=====================================================================
+
+ subroutine wmo_movie_div_curl(NSPEC_AB,NGLOB_AB,veloc, &
+ div_glob,curl_glob,valency, &
+ div,curl_x,curl_y,curl_z, &
+ velocity_x,velocity_y,velocity_z, &
+ ibool,ispec_is_elastic, &
+ hprime_xx,hprime_yy,hprime_zz, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz)
+
+
+! calculates div, curl and velocity
+
+ implicit none
+ include 'constants.h'
+
+ integer :: NSPEC_AB,NGLOB_AB
+
+ ! velocity field
+ real(kind=CUSTOM_REAL),dimension(NDIM,NGLOB_AB),intent(in) :: veloc
+
+ ! divergence and curl only in the global nodes
+ real(kind=CUSTOM_REAL),dimension(NGLOB_AB) :: div_glob,curl_glob
+ integer,dimension(NGLOB_AB) :: valency
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: div, curl_x, curl_y, curl_z
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: velocity_x,velocity_y,velocity_z
+ integer,dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB):: ibool
+ logical,dimension(NSPEC_AB) :: ispec_is_elastic
+
+ ! array with derivatives of Lagrange polynomials
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx
+ real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLY) :: hprime_yy
+ real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ) :: hprime_zz
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
+
+ ! local parameters
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: dvxdxl,dvxdyl,&
+ dvxdzl,dvydxl,dvydyl,dvydzl,dvzdxl,dvzdyl,dvzdzl
+ real(kind=CUSTOM_REAL) xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl
+ real(kind=CUSTOM_REAL) hp1,hp2,hp3
+ real(kind=CUSTOM_REAL) tempx1l,tempx2l,tempx3l
+ real(kind=CUSTOM_REAL) tempy1l,tempy2l,tempy3l
+ real(kind=CUSTOM_REAL) tempz1l,tempz2l,tempz3l
+ integer :: ispec,i,j,k,l,iglob
+
+ ! initializes
+ div_glob(:) = 0.0_CUSTOM_REAL
+ curl_glob(:) = 0.0_CUSTOM_REAL
+ valency(:) = 0
+
+ ! loops over elements
+ do ispec=1,NSPEC_AB
+ if( .not. ispec_is_elastic(ispec) ) cycle
+
+ ! calculates divergence and curl of velocity field
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ tempx1l = 0._CUSTOM_REAL
+ tempx2l = 0._CUSTOM_REAL
+ tempx3l = 0._CUSTOM_REAL
+ tempy1l = 0._CUSTOM_REAL
+ tempy2l = 0._CUSTOM_REAL
+ tempy3l = 0._CUSTOM_REAL
+ tempz1l = 0._CUSTOM_REAL
+ tempz2l = 0._CUSTOM_REAL
+ tempz3l = 0._CUSTOM_REAL
+
+ do l=1,NGLLX
+ hp1 = hprime_xx(i,l)
+ iglob = ibool(l,j,k,ispec)
+ tempx1l = tempx1l + veloc(1,iglob)*hp1
+ tempy1l = tempy1l + veloc(2,iglob)*hp1
+ tempz1l = tempz1l + veloc(3,iglob)*hp1
+ hp2 = hprime_yy(j,l)
+ iglob = ibool(i,l,k,ispec)
+ tempx2l = tempx2l + veloc(1,iglob)*hp2
+ tempy2l = tempy2l + veloc(2,iglob)*hp2
+ tempz2l = tempz2l + veloc(3,iglob)*hp2
+ hp3 = hprime_zz(k,l)
+ iglob = ibool(i,j,l,ispec)
+ tempx3l = tempx3l + veloc(1,iglob)*hp3
+ tempy3l = tempy3l + veloc(2,iglob)*hp3
+ tempz3l = tempz3l + veloc(3,iglob)*hp3
+ enddo
+
+ ! get derivatives of ux, uy and uz with respect to x, y and z
+ xixl = xix(i,j,k,ispec)
+ xiyl = xiy(i,j,k,ispec)
+ xizl = xiz(i,j,k,ispec)
+ etaxl = etax(i,j,k,ispec)
+ etayl = etay(i,j,k,ispec)
+ etazl = etaz(i,j,k,ispec)
+ gammaxl = gammax(i,j,k,ispec)
+ gammayl = gammay(i,j,k,ispec)
+ gammazl = gammaz(i,j,k,ispec)
+
+ dvxdxl(i,j,k) = xixl*tempx1l + etaxl*tempx2l + gammaxl*tempx3l
+ dvxdyl(i,j,k) = xiyl*tempx1l + etayl*tempx2l + gammayl*tempx3l
+ dvxdzl(i,j,k) = xizl*tempx1l + etazl*tempx2l + gammazl*tempx3l
+
+ dvydxl(i,j,k) = xixl*tempy1l + etaxl*tempy2l + gammaxl*tempy3l
+ dvydyl(i,j,k) = xiyl*tempy1l + etayl*tempy2l + gammayl*tempy3l
+ dvydzl(i,j,k) = xizl*tempy1l + etazl*tempy2l + gammazl*tempy3l
+
+ dvzdxl(i,j,k) = xixl*tempz1l + etaxl*tempz2l + gammaxl*tempz3l
+ dvzdyl(i,j,k) = xiyl*tempz1l + etayl*tempz2l + gammayl*tempz3l
+ dvzdzl(i,j,k) = xizl*tempz1l + etazl*tempz2l + gammazl*tempz3l
+
+ enddo
+ enddo
+ enddo
+
+ do k = 1,NGLLZ
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+ ! divergence \nabla \cdot \bf{v}
+ div(i,j,k,ispec) = dvxdxl(i,j,k) + dvydyl(i,j,k) + dvzdzl(i,j,k)
+ ! curl
+ curl_x(i,j,k,ispec) = dvzdyl(i,j,k) - dvydzl(i,j,k)
+ curl_y(i,j,k,ispec) = dvxdzl(i,j,k) - dvzdxl(i,j,k)
+ curl_z(i,j,k,ispec) = dvydxl(i,j,k) - dvxdyl(i,j,k)
+ ! velocity field
+ iglob = ibool(i,j,k,ispec)
+ velocity_x(i,j,k,ispec) = veloc(1,iglob)
+ velocity_y(i,j,k,ispec) = veloc(2,iglob)
+ velocity_z(i,j,k,ispec) = veloc(3,iglob)
+
+ valency(iglob)=valency(iglob)+1
+
+ div_glob(iglob) = div_glob(iglob) + div(i,j,k,ispec)
+ curl_glob(iglob)= curl_glob(iglob) &
+ + 0.5_CUSTOM_REAL*(curl_x(i,j,k,ispec) &
+ + curl_x(i,j,k,ispec)+curl_x(i,j,k,ispec))
+ enddo
+ enddo
+ enddo
+ enddo !NSPEC_AB
+
+ do i=1,NGLOB_AB
+ ! checks if point has a contribution
+ ! note: might not be the case for points in acoustic elements
+ if( valency(i) /= 0 ) then
+ ! averages by number of contributions
+ div_glob(i) = div_glob(i)/valency(i)
+ curl_glob(i) = curl_glob(i)/valency(i)
+ endif
+
+ enddo
+
+ end subroutine wmo_movie_div_curl
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_output_ASCII.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_output_ASCII.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_output_ASCII.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -50,7 +50,7 @@
character(len=256) sisname,final_LOCAL_PATH
! local parameter
- integer isample,nt_s
+ integer isample,nt_s,ier
real, dimension(:), allocatable :: tr
real(kind=CUSTOM_REAL) :: time_t
@@ -58,7 +58,8 @@
if(SEISMOGRAMS_BINARY)then
! allocate trace
nt_s = min(it,NSTEP)
- allocate(tr(nt_s))
+ allocate(tr(nt_s),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array tr'
! binary format case
open(unit=IOUT, file=final_LOCAL_PATH(1:len_trim(final_LOCAL_PATH))//&
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_seismograms.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_seismograms.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_seismograms.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -47,9 +47,15 @@
do irec_local = 1,nrec_local
- ! get global number of that receiver
+ ! gets global number of that receiver
irec = number_receiver_global(irec_local)
+ ! gets local receiver interpolators
+ ! (1-D Lagrange interpolators)
+ hxir(:) = hxir_store(irec_local,:)
+ hetar(:) = hetar_store(irec_local,:)
+ hgammar(:) = hgammar_store(irec_local,:)
+
! forward simulations
if (SIMULATION_TYPE == 1) then
@@ -62,8 +68,7 @@
call compute_interpolated_dva(displ,veloc,accel,NGLOB_AB, &
ispec,NSPEC_AB,ibool, &
xi_receiver(irec),eta_receiver(irec),gamma_receiver(irec), &
- hxir_store(irec_local,:),hetar_store(irec_local,:), &
- hgammar_store(irec_local,:), &
+ hxir,hetar,hgammar, &
dxd,dyd,dzd,vxd,vyd,vzd,axd,ayd,azd)
endif !elastic
@@ -82,14 +87,13 @@
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
ibool,rhostore)
- ! interpolates displ/veloc/pressure at receiver locations
+ ! interpolates displ/veloc/pressure at receiver locations
call compute_interpolated_dva_ac(displ_element,veloc_element,&
potential_dot_dot_acoustic,potential_dot_acoustic,&
potential_acoustic,NGLOB_AB, &
ispec,NSPEC_AB,ibool, &
xi_receiver(irec),eta_receiver(irec),gamma_receiver(irec), &
- hxir_store(irec_local,:),hetar_store(irec_local,:), &
- hgammar_store(irec_local,:), &
+ hxir,hetar,hgammar, &
dxd,dyd,dzd,vxd,vyd,vzd,axd,ayd,azd)
endif ! acoustic
@@ -105,8 +109,7 @@
call compute_interpolated_dva(displ,veloc,accel,NGLOB_AB, &
ispec,NSPEC_AB,ibool, &
xi_receiver(irec),eta_receiver(irec),gamma_receiver(irec), &
- hxir_store(irec_local,:),hetar_store(irec_local,:), &
- hgammar_store(irec_local,:), &
+ hxir,hetar,hgammar, &
dxd,dyd,dzd,vxd,vyd,vzd,axd,ayd,azd)
! stores elements displacement field
@@ -119,11 +122,15 @@
enddo
enddo
+ ! gets derivatives of local receiver interpolators
+ hpxir(:) = hpxir_store(irec_local,:)
+ hpetar(:) = hpetar_store(irec_local,:)
+ hpgammar(:) = hpgammar_store(irec_local,:)
+
! computes the integrated derivatives of source parameters (M_jk and X_s)
call compute_adj_source_frechet(displ_element,Mxx(irec),Myy(irec),Mzz(irec),&
Mxy(irec),Mxz(irec),Myz(irec),eps_s,eps_m_s, &
- hxir_store(irec_local,:),hetar_store(irec_local,:),hgammar_store(irec_local,:), &
- hpxir_store(irec_local,:),hpetar_store(irec_local,:),hpgammar_store(irec_local,:), &
+ hxir,hetar,hgammar,hpxir,hpetar,hpgammar, &
hprime_xx,hprime_yy,hprime_zz, &
xix(:,:,:,ispec),xiy(:,:,:,ispec),xiz(:,:,:,ispec), &
etax(:,:,:,ispec),etay(:,:,:,ispec),etaz(:,:,:,ispec), &
@@ -162,8 +169,7 @@
potential_acoustic,NGLOB_AB, &
ispec,NSPEC_AB,ibool, &
xi_receiver(irec),eta_receiver(irec),gamma_receiver(irec), &
- hxir_store(irec_local,:),hetar_store(irec_local,:), &
- hgammar_store(irec_local,:), &
+ hxir,hetar,hgammar, &
dxd,dyd,dzd,vxd,vyd,vzd,axd,ayd,azd)
endif ! acoustic
@@ -178,8 +184,7 @@
call compute_interpolated_dva(b_displ,b_veloc,b_accel,NGLOB_ADJOINT,&
ispec,NSPEC_AB,ibool, &
xi_receiver(irec),eta_receiver(irec),gamma_receiver(irec), &
- hxir_store(irec_local,:),hetar_store(irec_local,:), &
- hgammar_store(irec_local,:), &
+ hxir,hetar,hgammar, &
dxd,dyd,dzd,vxd,vyd,vzd,axd,ayd,azd)
endif ! elastic
@@ -204,8 +209,7 @@
b_potential_acoustic,NGLOB_ADJOINT, &
ispec,NSPEC_AB,ibool, &
xi_receiver(irec),eta_receiver(irec),gamma_receiver(irec), &
- hxir_store(irec_local,:),hetar_store(irec_local,:), &
- hgammar_store(irec_local,:), &
+ hxir,hetar,hgammar, &
dxd,dyd,dzd,vxd,vyd,vzd,axd,ayd,azd)
endif ! acoustic
Modified: seismo/3D/SPECFEM3D/trunk/utils/Visualization/Paraview/create_slice_VTK.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/utils/Visualization/Paraview/create_slice_VTK.f90 2011-03-30 23:15:27 UTC (rev 18153)
+++ seismo/3D/SPECFEM3D/trunk/utils/Visualization/Paraview/create_slice_VTK.f90 2011-03-31 11:08:11 UTC (rev 18154)
@@ -57,7 +57,7 @@
character(len=256) :: mesh_file,local_data_file
integer, dimension(300) :: node_list
- integer :: iproc, num_node, i,ios, it
+ integer :: iproc, num_node, i,ios, it, ier
integer :: njunk
! data must be of dimension: (NGLLX,NGLLY,NGLLZ,NSPEC_AB)
@@ -131,11 +131,13 @@
read(27) NGLOB_AB
! ibool file
- allocate(ibool(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(ibool(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array ibool'
read(27) ibool
! global point arrays
- allocate(xstore(NGLOB_AB),ystore(NGLOB_AB),zstore(NGLOB_AB))
+ allocate(xstore(NGLOB_AB),ystore(NGLOB_AB),zstore(NGLOB_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array xstore etc.'
read(27) xstore
read(27) ystore
read(27) zstore
@@ -152,7 +154,8 @@
print *,'Error opening ',trim(local_data_file)
stop
endif
- allocate(data(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(data(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array data'
read(28) data
close(28)
More information about the CIG-COMMITS
mailing list