[cig-commits] r20590 - in seismo/3D/SPECFEM3D/trunk/src: decompose_mesh_SCOTCH generate_databases meshfem3D shared specfem3D
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Sat Aug 18 10:02:11 PDT 2012
Author: dkomati1
Date: 2012-08-18 10:02:11 -0700 (Sat, 18 Aug 2012)
New Revision: 20590
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/calc_jacobian.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_absorbing_boundary.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_coupling_surfaces.f90
seismo/3D/SPECFEM3D/trunk/src/generate_databases/memory_eval.f90
seismo/3D/SPECFEM3D/trunk/src/meshfem3D/check_mesh_quality.f90
seismo/3D/SPECFEM3D/trunk/src/meshfem3D/create_regions_mesh.f90
seismo/3D/SPECFEM3D/trunk/src/meshfem3D/define_subregions.f90
seismo/3D/SPECFEM3D/trunk/src/meshfem3D/define_superbrick.f90
seismo/3D/SPECFEM3D/trunk/src/meshfem3D/meshfem3D.f90
seismo/3D/SPECFEM3D/trunk/src/meshfem3D/store_boundaries.f90
seismo/3D/SPECFEM3D/trunk/src/shared/check_mesh_resolution.f90
seismo/3D/SPECFEM3D/trunk/src/shared/detect_surface.f90
seismo/3D/SPECFEM3D/trunk/src/shared/get_element_face.f90
seismo/3D/SPECFEM3D/trunk/src/shared/get_jacobian_boundaries.f90
seismo/3D/SPECFEM3D/trunk/src/shared/get_shape3D.f90
seismo/3D/SPECFEM3D/trunk/src/shared/prepare_assemble_MPI.f90
seismo/3D/SPECFEM3D/trunk/src/specfem3D/PML_init.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/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
Log:
removed old code that was commented out; removed useless white spaces
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 2012-08-18 00:54:58 UTC (rev 20589)
+++ seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 2012-08-18 17:02:11 UTC (rev 20590)
@@ -134,7 +134,6 @@
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)
- !if(num_node /= inode) stop "ERROR : Invalid nodes_coords file."
end do
close(98)
print*, 'total number of nodes: '
@@ -166,32 +165,17 @@
! 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
+ ! here our code assumes that element ordering is:
! 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)
- !read(98,*) num_elmnt, elmnts(5,num_elmnt), elmnts(1,num_elmnt),elmnts(4,num_elmnt), elmnts(8,num_elmnt), &
- ! elmnts(6,num_elmnt), elmnts(2,num_elmnt), elmnts(3,num_elmnt), elmnts(7,num_elmnt)
-
read(98,*) num_elmnt, elmnts(1,num_elmnt), elmnts(2,num_elmnt),elmnts(3,num_elmnt), elmnts(4,num_elmnt), &
elmnts(5,num_elmnt), elmnts(6,num_elmnt), elmnts(7,num_elmnt), elmnts(8,num_elmnt)
if((num_elmnt > nspec) .or. (num_elmnt < 1) ) stop "ERROR : Invalid mesh file."
-
- !outputs info for each element to see ordering
- !print*,'ispec: ',ispec
- !print*,' ',num_elmnt, elmnts(5,num_elmnt), elmnts(1,num_elmnt),elmnts(4,num_elmnt), elmnts(8,num_elmnt), &
- ! elmnts(6,num_elmnt), elmnts(2,num_elmnt), elmnts(3,num_elmnt), elmnts(7,num_elmnt)
- !print*,'elem:',num_elmnt
- !do i=1,8
- ! print*,' i ',i,'val :',elmnts(i,num_elmnt),&
- ! nodes_coords(1,elmnts(i,num_elmnt)),nodes_coords(2,elmnts(i,num_elmnt)),nodes_coords(3,elmnts(i,num_elmnt))
- !enddo
- !print*
-
end do
close(98)
print*, 'total number of spectral elements:'
@@ -207,7 +191,7 @@
do ispec = 1, nspec
! format: # id_element #flag
! note: be aware that elements may not be sorted in materials_file
- read(98,*) num_mat,mat(1,num_mat) !mat(1,ispec)!, mat(2,ispec)
+ read(98,*) num_mat,mat(1,num_mat)
if((num_mat > nspec) .or. (num_mat < 1) ) stop "ERROR : Invalid mat file."
end do
close(98)
@@ -309,7 +293,6 @@
! format: note that we save the arguments in a slightly different order in mat_prop(:,:)
! #(6) material_domain_id #(0) material_id #(1) rho #(2) vp #(3) vs #(4) Q_mu #(5) anisotropy_flag
!
- !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. ier == 0)
@@ -327,8 +310,6 @@
if(idomain_id == 1 .or. idomain_id == 2) then
! material is elastic or acoustic
- !read(98,*) num_mat, mat_prop(1,num_mat),mat_prop(2,num_mat),&
- ! mat_prop(3,num_mat),mat_prop(4,num_mat),mat_prop(5,num_mat)
mat_prop(1,num_mat) = rho
mat_prop(2,num_mat) = vp
mat_prop(3,num_mat) = vs
@@ -374,7 +355,7 @@
! - 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
+ ! reads lines until it reaches a defined material
num_mat = 1
do while( num_mat >= 0 .and. ier == 0 )
read(98,'(A256)',iostat=ier) line
@@ -397,17 +378,8 @@
stop "ERROR: invalid line in nummaterial_velocity_file for undefined material"
endif
- !read(98,*) undef_mat_prop(6,imat),undef_mat_prop(1,imat),undef_mat_prop(2,imat),&
- ! undef_mat_prop(3,imat),undef_mat_prop(4,imat),undef_mat_prop(5,imat)
-
- ! debug output
- !print*,'properties:'
- !print*,undef_mat_prop(:,imat)
- !print*
-
! checks material_id
read(undef_mat_prop(1,imat),*) num_mat
- !print *,'material_id: ',num_mat
if(num_mat > 0 .or. -num_mat > count_undef_mat) &
stop "ERROR : Invalid nummaterial_velocity_file for undefined materials."
if(num_mat /= -imat) &
@@ -491,17 +463,8 @@
! doesn't necessarily have to start on top-rear, then bottom-rear, bottom-front, and finally top-front i.e.:
! point 1 = (0,1,1), point 2 = (0,1,0), point 3 = (0,0,0), point 4 = (0,0,1)
read(98,*) ibelm_xmin(ispec2D), nodes_ibelm_xmin(1,ispec2D), nodes_ibelm_xmin(2,ispec2D), &
- nodes_ibelm_xmin(3,ispec2D), nodes_ibelm_xmin(4,ispec2D)
+ nodes_ibelm_xmin(3,ispec2D), nodes_ibelm_xmin(4,ispec2D) !! DK DK see if we need NGNOD here
- !outputs info for each element for check of ordering
- !print*,'ispec2d:',ispec2d
- !print*,' xmin:', ibelm_xmin(ispec2D), nodes_ibelm_xmin(1,ispec2D), nodes_ibelm_xmin(2,ispec2D), &
- ! nodes_ibelm_xmin(3,ispec2D), nodes_ibelm_xmin(4,ispec2D)
- !do i=1,4
- ! print*,'i',i,'val:',ibelm_xmin(ispec2d),nodes_coords(1,nodes_ibelm_xmin(i,ispec2D)), &
- ! nodes_coords(2,nodes_ibelm_xmin(i,ispec2D)),nodes_coords(3,nodes_ibelm_xmin(i,ispec2D))
- !enddo
- !print*
end do
close(98)
print*, 'absorbing boundaries:'
@@ -891,7 +854,6 @@
!print*, ipart,": nspec_local=",nspec_local, " nnodes_local=", nnodes_loc
! writes out node coordinate locations
- !write(IIN_database,*) nnodes_loc
write(IIN_database) nnodes_loc
@@ -903,7 +865,6 @@
mat_prop, undef_mat_prop)
! writes out spectral element indices
- !write(IIN_database,*) nspec_local
write(IIN_database) nspec_local
call write_partition_database(IIN_database, ipart, nspec_local, nspec, elmnts, &
glob2loc_elmnts, glob2loc_nodes_nparts, &
@@ -928,10 +889,8 @@
! writes out MPI interfaces elements
!print*,' my interfaces:',my_ninterface,maxval(my_nb_interfaces)
if( my_ninterface == 0 ) then
- !write(IIN_database,*) my_ninterface, 0
write(IIN_database) my_ninterface, 0 ! avoids problem with maxval for empty array my_nb_interfaces
else
- !write(IIN_database,*) my_ninterface, maxval(my_nb_interfaces)
write(IIN_database) my_ninterface, maxval(my_nb_interfaces)
endif
@@ -961,3 +920,4 @@
!end program pre_meshfem3D
end module decompose_mesh_SCOTCH
+
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 2012-08-18 00:54:58 UTC (rev 20589)
+++ seismo/3D/SPECFEM3D/trunk/src/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90 2012-08-18 17:02:11 UTC (rev 20590)
@@ -121,14 +121,14 @@
xadj(nodes_elmnts(k+j*nsize)) = xadj(nodes_elmnts(k+j*nsize)) + 1
if (xadj(nodes_elmnts(k+j*nsize))>sup_neighbour) &
- stop 'ERROR : too much neighbours per element, modify the mesh.'
+ stop 'ERROR: too many neighbours per element, modify the mesh.'
adjncy(nodes_elmnts(l+j*nsize)*sup_neighbour &
+ xadj(nodes_elmnts(l+j*nsize))) = nodes_elmnts(k+j*nsize)
xadj(nodes_elmnts(l+j*nsize)) = xadj(nodes_elmnts(l+j*nsize)) + 1
if (xadj(nodes_elmnts(l+j*nsize))>sup_neighbour) &
- stop 'ERROR : too much neighbours per element, modify the mesh.'
+ stop 'ERROR: too many neighbours per element, modify the mesh.'
end if
end if
end do
@@ -149,11 +149,9 @@
xadj(nspec) = nb_edges
-
end subroutine mesh2dual_ncommonnodes
-
!--------------------------------------------------
! construct local numbering for the elements in each partition
!--------------------------------------------------
@@ -185,7 +183,6 @@
num_loc(num_part) = num_loc(num_part) + 1
end do
-
end subroutine build_glob2loc_elmnts
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/calc_jacobian.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/calc_jacobian.f90 2012-08-18 00:54:58 UTC (rev 20589)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/calc_jacobian.f90 2012-08-18 17:02:11 UTC (rev 20590)
@@ -145,220 +145,3 @@
end subroutine calc_jacobian
-!
-!-------------------------------------------------------------------------------------------------
-!
-
-! This subroutine recomputes the 3D jacobian for one element
-! based upon all GLL points
-! Hejun Zhu OCT16,2009
-
-! input: myrank,
-! xstore,ystore,zstore ----- input position
-! xigll,yigll,zigll ----- gll points position
-! ispec,nspec ----- element number
-! ACTUALLY_STORE_ARRAYS ------ save array or not
-
-! output: xixstore,xiystore,xizstore,
-! etaxstore,etaystore,etazstore,
-! gammaxstore,gammaystore,gammazstore ------ parameters used for calculating jacobian
-!
-!
-! subroutine recalc_jacobian_gll3D(myrank,xixstore,xiystore,xizstore, &
-! etaxstore,etaystore,etazstore, &
-! gammaxstore,gammaystore,gammazstore,jacobianstore, &
-! xstore,ystore,zstore, &
-! ispec,nspec, &
-! xigll,yigll,zigll, &
-! ACTUALLY_STORE_ARRAYS)
-!
-! implicit none
-!
-! include "constants.h"
-!
-! ! input parameter
-! integer::myrank,ispec,nspec
-! double precision, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: xstore,ystore,zstore
-! double precision, dimension(NGLLX):: xigll
-! double precision, dimension(NGLLY):: yigll
-! double precision, dimension(NGLLZ):: zigll
-! logical::ACTUALLY_STORE_ARRAYS
-!
-!
-! ! output results
-! real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: &
-! xixstore,xiystore,xizstore,&
-! etaxstore,etaystore,etazstore,&
-! gammaxstore,gammaystore,gammazstore,&
-! jacobianstore
-!
-!
-! ! other parameters for this subroutine
-! integer:: i,j,k,i1,j1,k1
-! double precision:: xxi,xeta,xgamma,yxi,yeta,ygamma,zxi,zeta,zgamma
-! double precision:: xi,eta,gamma
-! double precision,dimension(NGLLX):: hxir,hpxir
-! double precision,dimension(NGLLY):: hetar,hpetar
-! double precision,dimension(NGLLZ):: hgammar,hpgammar
-! double precision:: hlagrange,hlagrange_xi,hlagrange_eta,hlagrange_gamma
-! double precision:: jacobian
-! double precision:: xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
-!
-!
-!
-! ! test parameters which can be deleted
-! double precision:: xmesh,ymesh,zmesh
-! double precision:: sumshape,sumdershapexi,sumdershapeeta,sumdershapegamma
-!
-! ! first go over all 125 gll points
-! do k=1,NGLLZ
-! do j=1,NGLLY
-! do i=1,NGLLX
-!
-! xxi = 0.0
-! xeta = 0.0
-! xgamma = 0.0
-! yxi = 0.0
-! yeta = 0.0
-! ygamma = 0.0
-! zxi = 0.0
-! zeta = 0.0
-! zgamma = 0.0
-!
-! xi = xigll(i)
-! eta = yigll(j)
-! gamma = zigll(k)
-!
-! ! calculate lagrange polynomial and its derivative
-! call lagrange_any(xi,NGLLX,xigll,hxir,hpxir)
-! call lagrange_any(eta,NGLLY,yigll,hetar,hpetar)
-! call lagrange_any(gamma,NGLLZ,zigll,hgammar,hpgammar)
-!
-! ! test parameters
-! sumshape = 0.0
-! sumdershapexi = 0.0
-! sumdershapeeta = 0.0
-! sumdershapegamma = 0.0
-! xmesh = 0.0
-! ymesh = 0.0
-! zmesh = 0.0
-!
-!
-! do k1 = 1,NGLLZ
-! do j1 = 1,NGLLY
-! do i1 = 1,NGLLX
-! hlagrange = hxir(i1)*hetar(j1)*hgammar(k1)
-! hlagrange_xi = hpxir(i1)*hetar(j1)*hgammar(k1)
-! hlagrange_eta = hxir(i1)*hpetar(j1)*hgammar(k1)
-! hlagrange_gamma = hxir(i1)*hetar(j1)*hpgammar(k1)
-!
-!
-! xxi = xxi + xstore(i1,j1,k1,ispec)*hlagrange_xi
-! xeta = xeta + xstore(i1,j1,k1,ispec)*hlagrange_eta
-! xgamma = xgamma + xstore(i1,j1,k1,ispec)*hlagrange_gamma
-!
-! yxi = yxi + ystore(i1,j1,k1,ispec)*hlagrange_xi
-! yeta = yeta + ystore(i1,j1,k1,ispec)*hlagrange_eta
-! ygamma = ygamma + ystore(i1,j1,k1,ispec)*hlagrange_gamma
-!
-! zxi = zxi + zstore(i1,j1,k1,ispec)*hlagrange_xi
-! zeta = zeta + zstore(i1,j1,k1,ispec)*hlagrange_eta
-! zgamma = zgamma + zstore(i1,j1,k1,ispec)*hlagrange_gamma
-!
-! ! test the lagrange polynomial and its derivate
-! xmesh = xmesh + xstore(i1,j1,k1,ispec)*hlagrange
-! ymesh = ymesh + ystore(i1,j1,k1,ispec)*hlagrange
-! zmesh = zmesh + zstore(i1,j1,k1,ispec)*hlagrange
-! sumshape = sumshape + hlagrange
-! sumdershapexi = sumdershapexi + hlagrange_xi
-! sumdershapeeta = sumdershapeeta + hlagrange_eta
-! sumdershapegamma = sumdershapegamma + hlagrange_gamma
-!
-! end do
-! end do
-! end do
-!
-! ! Check the lagrange polynomial and its derivative
-! if (xmesh /=xstore(i,j,k,ispec).or.ymesh/=ystore(i,j,k,ispec).or.zmesh/=zstore(i,j,k,ispec)) then
-! call exit_MPI(myrank,'new mesh positions are wrong in recalc_jacobian_gall3D.f90')
-! end if
-! if(abs(sumshape-one) > TINYVAL) then
-! call exit_MPI(myrank,'error shape functions in recalc_jacobian_gll3D.f90')
-! end if
-! if(abs(sumdershapexi) > TINYVAL) then
-! call exit_MPI(myrank,'error derivative xi shape functions in recalc_jacobian_gll3D.f90')
-! end if
-! if(abs(sumdershapeeta) > TINYVAL) then
-! call exit_MPI(myrank,'error derivative eta shape functions in recalc_jacobian_gll3D.f90')
-! end if
-! if(abs(sumdershapegamma) > TINYVAL) then
-! call exit_MPI(myrank,'error derivative gamma shape functions in recalc_jacobian_gll3D.f90')
-! end if
-!
-!
-! jacobian = xxi*(yeta*zgamma-ygamma*zeta) - &
-! xeta*(yxi*zgamma-ygamma*zxi) + &
-! xgamma*(yxi*zeta-yeta*zxi)
-!
-! ! Check the jacobian
-! if(jacobian <= ZERO) then
-! call exit_MPI(myrank,'3D Jacobian undefined in recalc_jacobian_gll3D.f90')
-! end if
-!
-! ! invert the relation (Fletcher p. 50 vol. 2)
-! xix = (yeta*zgamma-ygamma*zeta) / jacobian
-! xiy = (xgamma*zeta-xeta*zgamma) / jacobian
-! xiz = (xeta*ygamma-xgamma*yeta) / jacobian
-! etax = (ygamma*zxi-yxi*zgamma) / jacobian
-! etay = (xxi*zgamma-xgamma*zxi) / jacobian
-! etaz = (xgamma*yxi-xxi*ygamma) / jacobian
-! gammax = (yxi*zeta-yeta*zxi) / jacobian
-! gammay = (xeta*zxi-xxi*zeta) / jacobian
-! gammaz = (xxi*yeta-xeta*yxi) / jacobian
-!
-!
-! ! compute and store the jacobian for the solver
-! jacobian = 1. / (xix*(etay*gammaz-etaz*gammay) &
-! -xiy*(etax*gammaz-etaz*gammax) &
-! +xiz*(etax*gammay-etay*gammax))
-!
-! ! resave the derivatives and the jacobian
-! ! distinguish between single and double precision for reals
-! if (ACTUALLY_STORE_ARRAYS) then
-!
-! if (myrank == 0) then
-! print*,'xix before',xixstore(i,j,k,ispec),'after',xix
-! print*,'etax before',etaxstore(i,j,k,ispec),'after',etax
-! print*,'gammax before',gammaxstore(i,j,k,ispec),'after',gammax
-! end if
-!
-! if(CUSTOM_REAL == SIZE_REAL) then
-! xixstore(i,j,k,ispec) = sngl(xix)
-! xiystore(i,j,k,ispec) = sngl(xiy)
-! xizstore(i,j,k,ispec) = sngl(xiz)
-! etaxstore(i,j,k,ispec) = sngl(etax)
-! etaystore(i,j,k,ispec) = sngl(etay)
-! etazstore(i,j,k,ispec) = sngl(etaz)
-! gammaxstore(i,j,k,ispec) = sngl(gammax)
-! gammaystore(i,j,k,ispec) = sngl(gammay)
-! gammazstore(i,j,k,ispec) = sngl(gammaz)
-! jacobianstore(i,j,k,ispec) = sngl(jacobian)
-! else
-! xixstore(i,j,k,ispec) = xix
-! xiystore(i,j,k,ispec) = xiy
-! xizstore(i,j,k,ispec) = xiz
-! etaxstore(i,j,k,ispec) = etax
-! etaystore(i,j,k,ispec) = etay
-! etazstore(i,j,k,ispec) = etaz
-! gammaxstore(i,j,k,ispec) = gammax
-! gammaystore(i,j,k,ispec) = gammay
-! gammazstore(i,j,k,ispec) = gammaz
-! jacobianstore(i,j,k,ispec) = jacobian
-! endif
-! end if
-! enddo
-! enddo
-! enddo
-!
-! end subroutine recalc_jacobian_gll3D
-!
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90 2012-08-18 00:54:58 UTC (rev 20589)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/create_regions_mesh.f90 2012-08-18 17:02:11 UTC (rev 20590)
@@ -196,11 +196,6 @@
double precision :: static_memory_size
real(kind=CUSTOM_REAL) :: model_speed_max,min_resolved_period
-! for vtk output
-! character(len=256) prname_file
-! integer,dimension(:),allocatable :: itest_flag
-! integer, dimension(:), allocatable :: elem_flag
-
! initializes arrays
call sync_all()
if( myrank == 0) then
@@ -434,21 +429,6 @@
! local parameters
integer :: ier
-! memory test
-! logical,dimension(:),allocatable :: test_mem
-!
-! tests memory availability (including some small buffer of 10*1024 byte)
-! allocate( test_mem(int(max_static_memory_size)+10*1024),stat=ier)
-! if(ier /= 0) then
-! write(IMAIN,*) 'error: try to increase the available process stack size by'
-! write(IMAIN,*) ' ulimit -s **** '
-! call exit_MPI(myrank,'not enough memory to allocate arrays')
-! endif
-! test_mem(:) = .true.
-! deallocate( test_mem, stat=ier)
-! if(ier /= 0) call exit_MPI(myrank,'error to allocate arrays')
-! call sync_all()
-
allocate( xelm(NGNOD),yelm(NGNOD),zelm(NGNOD),stat=ier)
if( ier /= 0 ) stop 'error allocating array xelm etc.'
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases.f90 2012-08-18 00:54:58 UTC (rev 20589)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/generate_databases.f90 2012-08-18 17:02:11 UTC (rev 20590)
@@ -550,20 +550,6 @@
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')
-! open(unit=55,file=BASEMENT_MAP_FILE,status='old',action='read')
-! do ix=1,NX_BASEMENT
-! do iy=1,NY_BASEMENT
-! read(55,*) iz_basement
-! z_basement(ix,iy) = dble(iz_basement)
-! enddo
-! enddo
-! close(55)
-! endif
-
end subroutine gd_read_topography
!
@@ -585,8 +571,6 @@
! read databases about external mesh simulation
! global node coordinates
call create_name_database(prname,myrank,LOCAL_PATH)
- !open(unit=IIN,file=prname(1:len_trim(prname))//'Database', &
- ! status='old',action='read',form='formatted',iostat=ier)
open(unit=IIN,file=prname(1:len_trim(prname))//'Database', &
status='old',action='read',form='unformatted',iostat=ier)
if( ier /= 0 ) then
@@ -594,19 +578,14 @@
print*,'please make sure file exists'
call exit_mpi(myrank,'error opening database file')
endif
- !read(IIN,*) nnodes_ext_mesh
read(IIN) 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)
-
read(IIN) dummy_node, nodes_coords_ext_mesh(1,inode), nodes_coords_ext_mesh(2,inode), &
nodes_coords_ext_mesh(3,inode)
-
enddo
call sum_all_i(nnodes_ext_mesh,num)
@@ -617,7 +596,6 @@
! read materials' physical properties
! added poroelastic properties and filled with 0 the last 10 entries for elastic/acoustic
- !read(IIN,*) nmat_ext_mesh, nundefMat_ext_mesh
read(IIN) nmat_ext_mesh, nundefMat_ext_mesh
allocate(materials_ext_mesh(16,nmat_ext_mesh),stat=ier)
@@ -626,23 +604,12 @@
! ielastic/acoustic format: #(1) rho #(2) vp #(3) vs #(4) Q_flag #(5) anisotropy_flag #(6) material_domain_id
! and remaining entries are filled with zeros
! poroelastic format: rhos,rhof,phi,tort,eta,material_domain_id,kxx,kxy,kxz,kyy,kyz,kzz,kappas,kappaf,kappafr,mufr
- !read(IIN,*) materials_ext_mesh(1,imat), materials_ext_mesh(2,imat), materials_ext_mesh(3,imat), &
- ! materials_ext_mesh(4,imat), materials_ext_mesh(5,imat), materials_ext_mesh(6,imat), &
- ! materials_ext_mesh(7,imat), materials_ext_mesh(8,imat), materials_ext_mesh(9,imat), &
- ! materials_ext_mesh(10,imat), materials_ext_mesh(11,imat), materials_ext_mesh(12,imat), &
- ! materials_ext_mesh(13,imat), materials_ext_mesh(14,imat), materials_ext_mesh(15,imat), &
- ! materials_ext_mesh(16,imat)
-
read(IIN) materials_ext_mesh(1,imat), materials_ext_mesh(2,imat), materials_ext_mesh(3,imat), &
materials_ext_mesh(4,imat), materials_ext_mesh(5,imat), materials_ext_mesh(6,imat), &
materials_ext_mesh(7,imat), materials_ext_mesh(8,imat), materials_ext_mesh(9,imat), &
materials_ext_mesh(10,imat), materials_ext_mesh(11,imat), materials_ext_mesh(12,imat), &
materials_ext_mesh(13,imat), materials_ext_mesh(14,imat), materials_ext_mesh(15,imat), &
materials_ext_mesh(16,imat)
-
- ! output
- !print*,'materials:',materials_ext_mesh(1,imat), materials_ext_mesh(2,imat), materials_ext_mesh(3,imat), &
- ! materials_ext_mesh(4,imat), materials_ext_mesh(5,imat), materials_ext_mesh(6,imat)
end do
if(myrank == 0) then
@@ -657,15 +624,8 @@
! -1 tomography elastic tomography_model.xyz 1 2
! format example interface:
! -1 interface 14 15 1 2
- !read(IIN,*) undef_mat_prop(1,imat),undef_mat_prop(2,imat),undef_mat_prop(3,imat),undef_mat_prop(4,imat), &
- ! undef_mat_prop(5,imat), undef_mat_prop(6,imat)
-
read(IIN) undef_mat_prop(1,imat),undef_mat_prop(2,imat),undef_mat_prop(3,imat),undef_mat_prop(4,imat), &
undef_mat_prop(5,imat), undef_mat_prop(6,imat)
-
- ! output debug
- !print*,'undefined materials:'
- !print*,undef_mat_prop(:,imat)
end do
if(myrank == 0) then
@@ -674,7 +634,6 @@
call sync_all()
! element indexing
- !read(IIN,*) nelmnts_ext_mesh
read(IIN) nelmnts_ext_mesh
allocate(elmnts_ext_mesh(NGNOD,nelmnts_ext_mesh),stat=ier)
if( ier /= 0 ) stop 'error allocating array elmnts_ext_mesh'
@@ -685,9 +644,6 @@
do ispec = 1, nelmnts_ext_mesh
! format:
! # ispec_local # material_index_1 # material_index_2 # corner_id1 # corner_id2 # ... # corner_id8
- !read(IIN,*) dummy_elmnt, mat_ext_mesh(1,ispec),mat_ext_mesh(2,ispec), &
- ! elmnts_ext_mesh(1,ispec), elmnts_ext_mesh(2,ispec), elmnts_ext_mesh(3,ispec), elmnts_ext_mesh(4,ispec), &
- ! elmnts_ext_mesh(5,ispec), elmnts_ext_mesh(6,ispec), elmnts_ext_mesh(7,ispec), elmnts_ext_mesh(8,ispec)
read(IIN) dummy_elmnt, mat_ext_mesh(1,ispec),mat_ext_mesh(2,ispec), &
elmnts_ext_mesh(1,ispec), elmnts_ext_mesh(2,ispec), elmnts_ext_mesh(3,ispec), elmnts_ext_mesh(4,ispec), &
@@ -705,29 +661,22 @@
endif
call sync_all()
-
! read boundaries
- !read(IIN,*) boundary_number ,nspec2D_xmin
read(IIN) boundary_number ,nspec2D_xmin
if(boundary_number /= 1) stop "Error : invalid database file"
- !read(IIN,*) boundary_number ,nspec2D_xmax
read(IIN) boundary_number ,nspec2D_xmax
if(boundary_number /= 2) stop "Error : invalid database file"
- !read(IIN,*) boundary_number ,nspec2D_ymin
read(IIN) boundary_number ,nspec2D_ymin
if(boundary_number /= 3) stop "Error : invalid database file"
- !read(IIN,*) boundary_number ,nspec2D_ymax
read(IIN) boundary_number ,nspec2D_ymax
if(boundary_number /= 4) stop "Error : invalid database file"
- !read(IIN,*) boundary_number ,nspec2D_bottom_ext
read(IIN) boundary_number ,nspec2D_bottom_ext
if(boundary_number /= 5) stop "Error : invalid database file"
- !read(IIN,*) boundary_number ,nspec2D_top_ext
read(IIN) boundary_number ,nspec2D_top_ext
if(boundary_number /= 6) stop "Error : invalid database file"
@@ -737,42 +686,36 @@
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)
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),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)
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),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)
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),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)
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),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)
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),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)
read(IIN) ibelm_top(ispec2D),(nodes_ibelm_top(j,ispec2D),j=1,4)
end do
@@ -794,7 +737,6 @@
! MPI interfaces between different partitions
if( NPROC > 1 ) then
! format: #number_of_MPI_interfaces #maximum_number_of_elements_on_each_interface
- !read(IIN,*) num_interfaces_ext_mesh, max_interface_size_ext_mesh
read(IIN) num_interfaces_ext_mesh, max_interface_size_ext_mesh
else
num_interfaces_ext_mesh = 0
@@ -819,7 +761,6 @@
! where
! process_interface_id = rank of (neighbor) process to share MPI interface with
! number_of_elements_on_interface = number of interface elements
- !read(IIN,*) my_neighbours_ext_mesh(num_interface), my_nelmnts_neighbours_ext_mesh(num_interface)
read(IIN) my_neighbours_ext_mesh(num_interface), my_nelmnts_neighbours_ext_mesh(num_interface)
! loops over interface elements
@@ -830,10 +771,6 @@
! 1 - corner point only
! 2 - element edge
! 4 - element face
- !read(IIN,*) my_interfaces_ext_mesh(1,ie,num_interface), my_interfaces_ext_mesh(2,ie,num_interface), &
- ! my_interfaces_ext_mesh(3,ie,num_interface), my_interfaces_ext_mesh(4,ie,num_interface), &
- ! my_interfaces_ext_mesh(5,ie,num_interface), my_interfaces_ext_mesh(6,ie,num_interface)
-
read(IIN) my_interfaces_ext_mesh(1,ie,num_interface), my_interfaces_ext_mesh(2,ie,num_interface), &
my_interfaces_ext_mesh(3,ie,num_interface), my_interfaces_ext_mesh(4,ie,num_interface), &
my_interfaces_ext_mesh(5,ie,num_interface), my_interfaces_ext_mesh(6,ie,num_interface)
@@ -849,16 +786,11 @@
! optional moho
if( SAVE_MOHO_MESH ) then
! checks if additional line exists
- !read(IIN,'(a128)',iostat=ier) line
read(IIN,iostat=ier) boundary_number,nspec2D_moho_ext
if( ier /= 0 ) then
! no moho informations given
nspec2D_moho_ext = 0
boundary_number = 7
- !else
- ! ! tries to read in number of moho elements
- ! read(line,*,iostat=ier) boundary_number ,nspec2D_moho_ext
- ! if( ier /= 0 ) call exit_mpi(myrank,'error reading moho mesh in database')
endif
if(boundary_number /= 7) stop "Error : invalid database file"
@@ -871,7 +803,6 @@
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)
read(IIN) ibelm_moho(ispec2D),(nodes_ibelm_moho(j,ispec2D),j=1,4)
end do
@@ -906,7 +837,6 @@
npointot = NSPEC_AB * NGLLCUBE
! use dynamic allocation to allocate memory for arrays
-! allocate(idoubling(NSPEC_AB))
allocate(ibool(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
if( ier /= 0 ) stop 'error allocating array ibool'
allocate(xstore(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
@@ -991,9 +921,9 @@
write(IMAIN,*) 'total number of elements in each slice: ',NSPEC_AB
write(IMAIN,*) 'total number of points in each slice: ',NGLOB_AB
write(IMAIN,*)
- write(IMAIN,*) 'total number of elements in entire mesh: ',nspec_total ! NSPEC_AB*NPROC
- write(IMAIN,*) 'total number of points in entire mesh: ',nglob_total !NGLOB_AB*NPROC
- write(IMAIN,*) 'total number of DOFs in entire mesh: ',nglob_total*NDIM !NGLOB_AB*NPROC*NDIM
+ write(IMAIN,*) 'total number of elements in entire mesh: ',nspec_total
+ write(IMAIN,*) 'total number of points in entire mesh: ',nglob_total
+ write(IMAIN,*) 'total number of DOFs in entire mesh: ',nglob_total*NDIM
write(IMAIN,*)
write(IMAIN,*) 'total number of time steps in the solver will be: ',NSTEP
write(IMAIN,*)
@@ -1042,7 +972,6 @@
! number of surface faces for all partitions together
call sum_all_i(nfaces_surface_ext_mesh,nfaces_surface_glob_ext_mesh)
-
! copy number of elements and points in an include file for the solver
if( myrank == 0 ) then
call save_header_file(NSPEC_AB,NGLOB_AB,NPROC, &
@@ -1070,3 +999,4 @@
call sync_all()
end subroutine gd_finalize
+
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_MPI.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_MPI.f90 2012-08-18 00:54:58 UTC (rev 20589)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_MPI.f90 2012-08-18 17:02:11 UTC (rev 20590)
@@ -143,7 +143,6 @@
write(*,*) ' interface:',iinterface,num_points1,num_points2
call exit_mpi(myrank,'error sorting MPI interface')
endif
- !write(*,*) myrank,'intfc',iinterface,num_points2,nibool_interfaces_ext_mesh_true(iinterface)
! cleanup temporary arrays
deallocate(xp)
@@ -205,8 +204,6 @@
ibool_interfaces_dummy(:,iinterface) = &
ibool_interfaces_ext_mesh(1:max_nibool_interfaces_ext_mesh,iinterface)
count = count + nibool_interfaces_ext_mesh(iinterface)
- !write(*,*) myrank,'interfaces ',iinterface, &
- ! nibool_interfaces_ext_mesh(iinterface),max_nibool_interfaces_ext_mesh
enddo
call sync_all()
@@ -243,3 +240,4 @@
endif
end subroutine get_MPI
+
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_absorbing_boundary.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_absorbing_boundary.f90 2012-08-18 00:54:58 UTC (rev 20589)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_absorbing_boundary.f90 2012-08-18 17:02:11 UTC (rev 20590)
@@ -98,15 +98,12 @@
! sets element
ispec = ibelm_xmin(ispec2D)
- !if(myrank == 0 ) print*,'xmin:',ispec2D,ispec
-
! looks for i,j,k indices of GLL points on boundary face
! determines element face by given CUBIT corners
do icorner=1,NGNOD2D
xcoord(icorner) = nodes_coords_ext_mesh(1,nodes_ibelm_xmin(icorner,ispec2D))
ycoord(icorner) = nodes_coords_ext_mesh(2,nodes_ibelm_xmin(icorner,ispec2D))
zcoord(icorner) = nodes_coords_ext_mesh(3,nodes_ibelm_xmin(icorner,ispec2D))
- !print*,'corner look:',icorner,xcoord(icorner),ycoord(icorner),zcoord(icorner)
enddo
! sets face id of reference element associated with this face
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_coupling_surfaces.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_coupling_surfaces.f90 2012-08-18 00:54:58 UTC (rev 20589)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/get_coupling_surfaces.f90 2012-08-18 17:02:11 UTC (rev 20590)
@@ -24,7 +24,6 @@
!
!=====================================================================
-
subroutine get_coupling_surfaces(myrank, &
nspec,ibool,NPROC, &
nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, &
@@ -208,11 +207,6 @@
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
-
- ! test vtk output
- !integer,dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: gll_data
- !character(len=256):: prname_file
-
! allocates temporary arrays
allocate(tmp_normal(NDIM,NGLLSQUARE,nspec*6),stat=ier)
if( ier /= 0 ) stop 'error allocating array tmp_normal'
@@ -387,11 +381,6 @@
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
-
- ! test vtk output
- !integer,dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: gll_data
- !character(len=256):: prname_file
-
! allocates temporary arrays
allocate(tmp_normal(NDIM,NGLLSQUARE,nspec*6),stat=ier)
if( ier /= 0 ) stop 'error allocating array tmp_normal'
@@ -558,11 +547,6 @@
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
-
- ! test vtk output
- !integer,dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: gll_data
- !character(len=256):: prname_file
-
! allocates temporary arrays
allocate(tmp_normal(NDIM,NGLLSQUARE,nspec*6),stat=ier)
if( ier /= 0 ) stop 'error allocating array tmp_normal'
@@ -623,8 +607,8 @@
if ( (iglob_corners_ref(1) == iglob_corners_ref_el(3)) .and. &
(iglob_corners_ref(3) == iglob_corners_ref_el(1)) ) then
- iface_ref_el = iface_el ![CM]: for some reason this shows a wrong orientation
- ! but the calcul is ok.
+ iface_ref_el = iface_el ![Christina Morency]: for some reason this shows a wrong orientation
+ ! but the calculation is ok.
ispec_ref_el = ispec_el
! gets face GLL points i,j,k indices from poroelastic element face
Modified: seismo/3D/SPECFEM3D/trunk/src/generate_databases/memory_eval.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/generate_databases/memory_eval.f90 2012-08-18 00:54:58 UTC (rev 20589)
+++ seismo/3D/SPECFEM3D/trunk/src/generate_databases/memory_eval.f90 2012-08-18 17:02:11 UTC (rev 20590)
@@ -204,3 +204,4 @@
end subroutine memory_eval_mesher
+
Modified: seismo/3D/SPECFEM3D/trunk/src/meshfem3D/check_mesh_quality.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/meshfem3D/check_mesh_quality.f90 2012-08-18 00:54:58 UTC (rev 20589)
+++ seismo/3D/SPECFEM3D/trunk/src/meshfem3D/check_mesh_quality.f90 2012-08-18 17:02:11 UTC (rev 20590)
@@ -66,15 +66,12 @@
double precision :: distance_min,distance_max
double precision :: distmin,distmax
- !double precision :: equiangle_skewness_min_MPI,edge_aspect_ratio_min_MPI,diagonal_aspect_ratio_min_MPI
double precision :: equiangle_skewness_max_MPI,edge_aspect_ratio_max_MPI,diagonal_aspect_ratio_max_MPI
double precision :: distance_min_MPI,distance_max_MPI
- !double precision :: distmin_MPI,distmax_MPI,skewness_AVS_DX_max_MPI,skewness_AVS_DX_max
- ! double precision :: skewness_AVS_DX_min_MPI,skewness_AVS_DX_min
! for stability
double precision :: dt_suggested,dt_suggested_max,dt_suggested_max_MPI
double precision :: stability,stability_min,stability_max,max_CFL_stability_limit
- double precision :: stability_max_MPI !,max_CFL_stability_limit_MPI,stability_MPI,stability_min_MPI
+ double precision :: stability_max_MPI
! For MPI maxloc reduction
double precision, dimension(2) :: buf_maxloc_send,buf_maxloc_recv
@@ -182,17 +179,13 @@
call max_all_dp(edge_aspect_ratio_max,edge_aspect_ratio_max_MPI)
call max_all_dp(diagonal_aspect_ratio_max,diagonal_aspect_ratio_max_MPI)
-
-
if((myrank == skewness_max_rank) .and. (myrank /= 0)) then
tmp_ispec_max_skewness(1) = ispec_max_skewness
call send_i_t(tmp_ispec_max_skewness,1,0)
end if
-
if(myrank == 0) then
-
if(skewness_max_rank /= myrank) then
call recv_i_t(tmp_ispec_max_skewness_MPI,1,skewness_max_rank)
ispec_max_skewness_MPI = tmp_ispec_max_skewness_MPI(1)
@@ -221,7 +214,6 @@
write(IMAIN,*) '*** max equiangle skewness = ',equiangle_skewness_max_MPI,' in element ',ispec_max_skewness_MPI, &
' of slice ',skewness_max_rank
write(IMAIN,*) '***'
- ! write(IMAIN,*) 'min equiangle skewness = ',equiangle_skewness_min
write(IMAIN,*)
write(IMAIN,*) 'max deviation angle from a right angle (90 degrees) is therefore = ',90.*equiangle_skewness_max_MPI
write(IMAIN,*)
@@ -229,19 +221,15 @@
write(IMAIN,*) 'or ',180. - 90.*(1. - equiangle_skewness_max_MPI),' degrees'
write(IMAIN,*)
write(IMAIN,*) 'max edge aspect ratio = ',edge_aspect_ratio_max_MPI
- ! write(IMAIN,*) 'min edge aspect ratio = ',edge_aspect_ratio_min
write(IMAIN,*)
write(IMAIN,*) 'max diagonal aspect ratio = ',diagonal_aspect_ratio_max_MPI
- ! write(IMAIN,*) 'min diagonal aspect ratio = ',diagonal_aspect_ratio_min
write(IMAIN,*)
- !write(IMAIN,*) 'max stability = ',stability_max_MPI
write(IMAIN,*) '***'
write(IMAIN,'(a50,f13.8)') ' *** Maximum suggested time step for simulation = ',dt_suggested_max_MPI
write(IMAIN,*) '***'
write(IMAIN,*) '*** max stability = ',stability_max_MPI
write(IMAIN,*) '*** computed using VP_MAX = ',VP_MAX
write(IMAIN,*) '***'
- ! write(IMAIN,*) 'min stability = ',stability_min
! max stability CFL value is different in 2D and in 3D
if(NGNOD == 8) then
@@ -266,13 +254,10 @@
endif
write(IMAIN,*)
-
-
! create statistics about mesh quality
write(IMAIN,*) 'creating histogram and statistics of mesh quality'
end if
-
! erase histogram of skewness
classes_skewness(:) = 0
@@ -291,7 +276,6 @@
enddo
-
! sum skewness results in all processes
do iclass = 0,NCLASS-1
call sum_all_i(classes_skewness(iclass),classes_skewnessMPI(iclass))
@@ -299,7 +283,6 @@
call sum_all_i(NSPEC,NSPEC_ALL_SLICES)
-
if(myrank == 0) then
! create histogram of skewness and save in Gnuplot file
write(IMAIN,*)
@@ -373,7 +356,7 @@
enddo
write(66,*) ""
- ! type: hexahedrons
+ ! type: hexahedra
write(66,'(a,i12)') "CELL_TYPES ",nspec
write(66,*) (12,ispec=1,nspec)
write(66,*) ""
@@ -560,4 +543,3 @@
end subroutine create_mesh_quality_data_3D
-
Modified: seismo/3D/SPECFEM3D/trunk/src/meshfem3D/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/meshfem3D/create_regions_mesh.f90 2012-08-18 00:54:58 UTC (rev 20589)
+++ seismo/3D/SPECFEM3D/trunk/src/meshfem3D/create_regions_mesh.f90 2012-08-18 17:02:11 UTC (rev 20590)
@@ -24,12 +24,9 @@
!
!=====================================================================
-
-
module createRegMesh
contains
-
subroutine create_regions_mesh(xgrid,ygrid,zgrid,ibool, &
xstore,ystore,zstore,iproc_xi,iproc_eta,addressing,nspec, &
NGLOB_AB,npointot, &
@@ -244,9 +241,7 @@
nmeshregions = 1
else
call define_superbrick(x_superbrick,y_superbrick,z_superbrick,ibool_superbrick,iboun_sb)
- !call define_superbrick_one_layer(x_superbrick,y_superbrick,z_superbrick,ibool_superbrick,iboun_sb)
nspec_sb = NSPEC_DOUBLING_SUPERBRICK
- !nspec_sb = NSPEC_SUPERBRICK_1L
if(NDOUBLINGS == 1) then
nmeshregions = 3
else if(NDOUBLINGS == 2) then
@@ -278,9 +273,6 @@
ioffset_y = iy+iay*iaddy(ia)
ioffset_z = ir+iar*iaddz(ia)
- !xelm(ia) = xgrid(ir+iar*iaddz(ia),ix+iax*iaddx(ia),iy+iay*iaddy(ia))
- !yelm(ia) = ygrid(ir+iar*iaddz(ia),ix+iax*iaddx(ia),iy+iay*iaddy(ia))
- !zelm(ia) = zgrid(ir+iar*iaddz(ia),ix+iax*iaddx(ia),iy+iay*iaddy(ia))
xelm(ia) = xgrid(ioffset_z,ioffset_x,ioffset_y)
yelm(ia) = ygrid(ioffset_z,ioffset_x,ioffset_y)
zelm(ia) = zgrid(ioffset_z,ioffset_x,ioffset_y)
@@ -447,9 +439,7 @@
ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top,&
NMATERIALS,material_properties)
-
end subroutine create_regions_mesh
end module createRegMesh
-
Modified: seismo/3D/SPECFEM3D/trunk/src/meshfem3D/define_subregions.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/meshfem3D/define_subregions.f90 2012-08-18 00:54:58 UTC (rev 20589)
+++ seismo/3D/SPECFEM3D/trunk/src/meshfem3D/define_subregions.f90 2012-08-18 17:02:11 UTC (rev 20590)
@@ -24,7 +24,6 @@
!
!=====================================================================
-
subroutine define_model_regions(NEX_PER_PROC_XI,NEX_PER_PROC_ETA,iproc_xi,iproc_eta,&
isubregion,nbsubregions,subregions,nblayers,ner_layer,&
iaddx,iaddy,iaddz,ix1,ix2,dix,iy1,iy2,diy,ir1,ir2,dir,iax,iay,iar, &
@@ -319,7 +318,5 @@
end if
-
end subroutine define_mesh_regions
-
Modified: seismo/3D/SPECFEM3D/trunk/src/meshfem3D/define_superbrick.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/meshfem3D/define_superbrick.f90 2012-08-18 00:54:58 UTC (rev 20589)
+++ seismo/3D/SPECFEM3D/trunk/src/meshfem3D/define_superbrick.f90 2012-08-18 17:02:11 UTC (rev 20590)
@@ -2038,4 +2038,6 @@
iboun_sb(8,4) = .true.
iboun_sb(8,5) = .true.
END SELECT
+
end subroutine define_basic_doubling_brick
+
Modified: seismo/3D/SPECFEM3D/trunk/src/meshfem3D/meshfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/meshfem3D/meshfem3D.f90 2012-08-18 00:54:58 UTC (rev 20589)
+++ seismo/3D/SPECFEM3D/trunk/src/meshfem3D/meshfem3D.f90 2012-08-18 17:02:11 UTC (rev 20590)
@@ -270,9 +270,6 @@
double precision Z_DEPTH_BLOCK !,Z_BASEMENT_SURFACE,Z_DEPTH_MOHO
double precision LATITUDE_MIN,LATITUDE_MAX,LONGITUDE_MIN,LONGITUDE_MAX
- !logical TOPOGRAPHY,MOHO_MAP_LUPEI
- !logical BASEMENT_MAP,HAUKSSON_REGIONAL_MODEL,HARVARD_3D_GOCAD_MODEL,IMPOSE_MINIMUM_VP_GOCAD
-
logical SUPPRESS_UTM_PROJECTION,USE_REGULAR_MESH
! Mesh files for visualization
@@ -282,7 +279,7 @@
integer NDOUBLINGS
integer, dimension(2) :: ner_doublings
- character(len=256) OUTPUT_FILES,LOCAL_PATH !,MODEL
+ character(len=256) OUTPUT_FILES,LOCAL_PATH
! parameters deduced from parameters read from file
integer NPROC,NEX_PER_PROC_XI,NEX_PER_PROC_ETA
@@ -296,14 +293,7 @@
NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX
double precision min_elevation,max_elevation
- !double precision min_elevation_all,max_elevation_all
-! for tapered basement map
- !integer icorner_x,icorner_y
- !integer iz_basement
- !double precision x_corner,y_corner
- !character(len=256) BASEMENT_MAP_FILE
-
! interfaces parameters
logical SUPPRESS_UTM_PROJECTION_BOTTOM,SUPPRESS_UTM_PROJECTION_TOP
integer ilayer,interface_current ! ipoint_current
Modified: seismo/3D/SPECFEM3D/trunk/src/meshfem3D/store_boundaries.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/meshfem3D/store_boundaries.f90 2012-08-18 00:54:58 UTC (rev 20589)
+++ seismo/3D/SPECFEM3D/trunk/src/meshfem3D/store_boundaries.f90 2012-08-18 17:02:11 UTC (rev 20590)
@@ -24,7 +24,6 @@
!
!=====================================================================
-
subroutine store_boundaries(myrank,iboun,nspec,&
ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, &
nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax, &
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/check_mesh_resolution.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/check_mesh_resolution.f90 2012-08-18 00:54:58 UTC (rev 20589)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/check_mesh_resolution.f90 2012-08-18 17:02:11 UTC (rev 20590)
@@ -65,12 +65,11 @@
! empirical choice for distorted elements to estimate time step and period resolved:
! courant number for time step estimate
- real(kind=CUSTOM_REAL),parameter :: COURANT_SUGGESTED = 0.3
+ real(kind=CUSTOM_REAL),parameter :: COURANT_SUGGESTED = 0.5 !! DK DK now that Stacey has been fixed, 0.3 is too low
! number of points per minimum wavelength for minimum period estimate
real(kind=CUSTOM_REAL),parameter :: NPTS_PER_WAVELENGTH = 5
!********************************************************************************
- !real(kind=CUSTOM_REAL),parameter :: NELEM_PER_WAVELENGTH = 1.5
logical :: has_vs_zero
real(kind=CUSTOM_REAL),dimension(1) :: tmp_val
@@ -179,7 +178,6 @@
pmax = avg_distance / min( vpmin,vsmin ) * NPTS_PER_WAVELENGTH
pmax_glob = max(pmax_glob,pmax)
-
! old: based on GLL distance, i.e. on maximum ratio ( gridspacing / velocity )
!pmax = distance_max / min( vpmin,vsmin ) * NELEM_PER_WAVELENGTH
!pmax_glob = max(pmax_glob,pmax)
@@ -279,7 +277,6 @@
write(IMAIN,*)
write(IMAIN,*) 'NSPEC_AB_global_min = ',NSPEC_AB_global_min
write(IMAIN,*) 'NSPEC_AB_global_max = ',NSPEC_AB_global_max
- !write(IMAIN,*) 'NSPEC_AB_global_mean = ',NSPEC_AB_global_sum / float(sizeprocs)
write(IMAIN,*) 'NSPEC_AB_global_sum = ',NSPEC_AB_global_sum
write(IMAIN,*)
write(IMAIN,*) 'NGLOB_AB_global_min = ',NGLOB_AB_global_min
@@ -399,12 +396,11 @@
! empirical choice for distorted elements to estimate time step and period resolved:
! courant number for time step estimate
- real(kind=CUSTOM_REAL),parameter :: COURANT_SUGGESTED = 0.3
+ real(kind=CUSTOM_REAL),parameter :: COURANT_SUGGESTED = 0.5 !! DK DK now that Stacey has been fixed, 0.3 is too low
! number of points per minimum wavelength for minimum period estimate
real(kind=CUSTOM_REAL),parameter :: NPTS_PER_WAVELENGTH = 5
!********************************************************************************
- !real(kind=CUSTOM_REAL),parameter :: NELEM_PER_WAVELENGTH = 1.5
logical :: has_vs_zero,has_vp2_zero
real(kind=CUSTOM_REAL),dimension(1) :: tmp_val
@@ -415,7 +411,6 @@
integer:: ier
character(len=256) :: filename,prname
-
! initializations
if( DT <= 0.0d0) then
DT_PRESENT = .false.
@@ -455,7 +450,6 @@
tmp2(:) = 0.0
endif
-
! checks courant number & minimum resolved period for each grid cell
do ispec=1,NSPEC_AB
@@ -641,7 +635,6 @@
write(IMAIN,*)
write(IMAIN,*) 'NSPEC_AB_global_min = ',NSPEC_AB_global_min
write(IMAIN,*) 'NSPEC_AB_global_max = ',NSPEC_AB_global_max
- !write(IMAIN,*) 'NSPEC_AB_global_mean = ',NSPEC_AB_global_sum / float(sizeprocs)
write(IMAIN,*) 'NSPEC_AB_global_sum = ',NSPEC_AB_global_sum
write(IMAIN,*)
write(IMAIN,*) 'NGLOB_AB_global_min = ',NGLOB_AB_global_min
@@ -880,7 +873,6 @@
!-------------------------------------------------------------------------------------------------
!
-
subroutine get_GLL_minmaxdistance(distance_min,distance_max,ispec, &
NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore)
@@ -941,12 +933,10 @@
end subroutine get_GLL_minmaxdistance
-
!
!-------------------------------------------------------------------------------------------------
!
-
subroutine get_elem_minmaxsize(elemsize_min,elemsize_max,ispec, &
NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore)
@@ -1012,3 +1002,4 @@
enddo
end subroutine get_elem_minmaxsize
+
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/detect_surface.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/detect_surface.f90 2012-08-18 00:54:58 UTC (rev 20589)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/detect_surface.f90 2012-08-18 17:02:11 UTC (rev 20590)
@@ -385,8 +385,6 @@
! tries to find closest face if points are inside
do ispec = 1,nspec
- ! checks if already assigned
- !if( ispec_is_surface_external_mesh(ispec) ) cycle
! in case element has still unresolved points in interior,
! we take closest element face to cross-section plane
@@ -437,8 +435,6 @@
! x cross-section plane
- !minface = minloc(midpoint_dist_x)
- !iface = minface(1)
i = iface_midpoint_ijk(1,iface)
j = iface_midpoint_ijk(2,iface)
k = iface_midpoint_ijk(3,iface)
@@ -455,8 +451,6 @@
endif
! y cross-section plane
- !minface = minloc(midpoint_dist_y)
- !iface = minface(1)
i = iface_midpoint_ijk(1,iface)
j = iface_midpoint_ijk(2,iface)
k = iface_midpoint_ijk(3,iface)
@@ -473,8 +467,6 @@
endif
! z cross-section plane
- !minface = minloc(midpoint_dist_z)
- !iface = minface(1)
i = iface_midpoint_ijk(1,iface)
j = iface_midpoint_ijk(2,iface)
k = iface_midpoint_ijk(3,iface)
@@ -772,3 +764,4 @@
num_iglob_image_surface = count
end subroutine detect_surface_PNM_GIF_image
+
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/get_element_face.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/get_element_face.f90 2012-08-18 00:54:58 UTC (rev 20589)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/get_element_face.f90 2012-08-18 17:02:11 UTC (rev 20590)
@@ -24,7 +24,6 @@
!
!=====================================================================
-
subroutine get_element_face_id(ispec,xcoord,ycoord,zcoord,&
ibool,nspec,nglob, &
xstore_dummy,ystore_dummy,zstore_dummy, &
@@ -81,7 +80,6 @@
iface5_corner_ijk,iface6_corner_ijk /),all_faces_shape)
! face orientation
- !real(kind=CUSTOM_REAL) :: face_n(3),face_ntmp(3),tmp
integer :: ifa,icorner,i,j,k,iglob,iloc(1)
! initializes
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/get_jacobian_boundaries.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/get_jacobian_boundaries.f90 2012-08-18 00:54:58 UTC (rev 20589)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/get_jacobian_boundaries.f90 2012-08-18 17:02:11 UTC (rev 20590)
@@ -258,9 +258,12 @@
end subroutine compute_jacobian_2D_face
+!
+!-------------------------
+!
-! This subroutine recompute the 3D jacobian for one element
-! based upon 125 GLL points
+! This subroutine recompute the jacobian for one element
+! based upon the GLL points
! Hejun Zhu, Oct 16, 2009
! input: myrank,
@@ -304,8 +307,6 @@
double precision:: jacobian
double precision:: unx,uny,unz
-
-
! test parameters which can be deleted
double precision:: xmesh,ymesh,zmesh
double precision:: sumshape,sumdershapexi,sumdershapeeta
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/get_shape3D.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/get_shape3D.f90 2012-08-18 00:54:58 UTC (rev 20589)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/get_shape3D.f90 2012-08-18 17:02:11 UTC (rev 20590)
@@ -267,4 +267,3 @@
end subroutine get_shape3D_element_corners
-
Modified: seismo/3D/SPECFEM3D/trunk/src/shared/prepare_assemble_MPI.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/shared/prepare_assemble_MPI.f90 2012-08-18 00:54:58 UTC (rev 20589)
+++ seismo/3D/SPECFEM3D/trunk/src/shared/prepare_assemble_MPI.f90 2012-08-18 17:02:11 UTC (rev 20590)
@@ -24,7 +24,6 @@
!
!=====================================================================
-
subroutine prepare_assemble_MPI (nelmnts,knods, &
ibool,npoin, &
ninterface, max_interface_size, &
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/PML_init.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/PML_init.f90 2012-08-18 00:54:58 UTC (rev 20589)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/PML_init.f90 2012-08-18 17:02:11 UTC (rev 20590)
@@ -1248,3 +1248,4 @@
if(myrank == 0) write(IMAIN,*)
end subroutine PML_output_VTKs
+
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_receivers.f90 2012-08-18 00:54:58 UTC (rev 20589)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_receivers.f90 2012-08-18 17:02:11 UTC (rev 20590)
@@ -316,76 +316,6 @@
altitude_rec(1) = loc_ele
distmin_ele(1) = loc_distmin
-! ! set distance to huge initial value
-! distmin = HUGEVAL
-! if(num_free_surface_faces > 0) then
-! iglob_selected = 1
-! ! loop only on points inside the element
-! ! exclude edges to ensure this point is not shared with other elements
-! imin = 2
-! imax = NGLLX - 1
-! jmin = 2
-! jmax = NGLLY - 1
-! iselected = 0
-! jselected = 0
-! iface_selected = 0
-! do iface=1,num_free_surface_faces
-! do j=jmin,jmax
-! do i=imin,imax
-!
-! ispec = free_surface_ispec(iface)
-! igll = free_surface_ijk(1,(j-1)*NGLLY+i,iface)
-! jgll = free_surface_ijk(2,(j-1)*NGLLY+i,iface)
-! kgll = free_surface_ijk(3,(j-1)*NGLLY+i,iface)
-! iglob = ibool(igll,jgll,kgll,ispec)
-!
-! ! keep this point if it is closer to the receiver
-! dist = (stutm_x(irec)-dble(xstore(iglob)))**2 + &
-! (stutm_y(irec)-dble(ystore(iglob)))**2
-! if(dist < distmin) then
-! distmin = dist
-! iglob_selected = iglob
-! iface_selected = iface
-! iselected = i
-! jselected = j
-! altitude_rec(1) = zstore(iglob_selected)
-! endif
-! enddo
-! enddo
-! ! end of loop on all the elements on the free surface
-! end do
-! ! weighted mean at current point of topography elevation of the four closest nodes
-! ! set distance to huge initial value
-! distmin = HUGEVAL
-! do j=jselected,jselected+1
-! do i=iselected,iselected+1
-! inode = 1
-! do jadjust=0,1
-! do iadjust= 0,1
-! ispec = free_surface_ispec(iface_selected)
-! igll = free_surface_ijk(1,(j-jadjust-1)*NGLLY+i-iadjust,iface_selected)
-! jgll = free_surface_ijk(2,(j-jadjust-1)*NGLLY+i-iadjust,iface_selected)
-! kgll = free_surface_ijk(3,(j-jadjust-1)*NGLLY+i-iadjust,iface_selected)
-! iglob = ibool(igll,jgll,kgll,ispec)
-!
-! elevation_node(inode) = zstore(iglob)
-! dist_node(inode) = dsqrt((stutm_x(irec)-dble(xstore(iglob)))**2 + &
-! (stutm_y(irec)-dble(ystore(iglob)))**2)
-! inode = inode + 1
-! end do
-! end do
-! dist = sum(dist_node)
-! if(dist < distmin) then
-! distmin = dist
-! altitude_rec(1) = (dist_node(1)/dist)*elevation_node(1) + &
-! (dist_node(2)/dist)*elevation_node(2) + &
-! (dist_node(3)/dist)*elevation_node(3) + &
-! (dist_node(4)/dist)*elevation_node(4)
-! endif
-! end do
-! end do
-! end if
-
! MPI communications to determine the best slice
call gather_all_dp(distmin_ele,1,distmin_ele_all,1,NPROC)
call gather_all_dp(altitude_rec,1,elevation_all,1,NPROC)
@@ -425,7 +355,6 @@
! burial depth in STATIONS file given in m
z_target(irec) = elevation(irec) - stbur(irec)
endif
- !if (myrank == 0) write(IOVTK,*) x_target(irec), y_target(irec), z_target(irec)
! reset distance to huge initial value
distmin=HUGEVAL
@@ -983,13 +912,6 @@
write(IMAIN,*) '************************************************************'
endif
- !! write the list of stations and associated epicentral distance
- !open(unit=27,file=trim(OUTPUT_FILES)//'/output_list_stations.txt',status='unknown')
- !do irec=1,nrec
- ! write(27,*) station_name(irec),'.',network_name(irec),' : ',horiz_dist(irec),' km horizontal distance'
- !enddo
- !close(27)
-
! write the locations of stations, so that we can load them and write them to SU headers later
open(unit=IOUT_SU,file=trim(OUTPUT_FILES)//'/output_list_stations.txt', &
status='unknown',action='write',iostat=ios)
@@ -1109,8 +1031,6 @@
! reads in station locations
open(unit=IIN, file=trim(filename), status = 'old', iostat = ios)
- !do irec = 1,nrec
- ! read(IIN,*) station_name,network_name,stlat,stlon,stele,stbur
do while(ios == 0)
read(IIN,"(a256)",iostat = ios) dummystring
if( ios /= 0 ) exit
@@ -1152,11 +1072,7 @@
if( stutm_y >= LATITUDE_MIN .and. stutm_y <= LATITUDE_MAX .and. &
stutm_x >= LONGITUDE_MIN .and. stutm_x <= LONGITUDE_MAX) then
- ! w/out formating
- ! write(IOUT,*) trim(station_name),' ',trim(network_name),' ',sngl(stlat), &
- ! ' ',sngl(stlon), ' ',sngl(stele), ' ',sngl(stbur)
-
- ! w/ specific format
+ ! with specific format
write(IOUT,'(a10,1x,a10,4e18.6)') &
trim(station_name),trim(network_name), &
sngl(stlat),sngl(stlon),sngl(stele),sngl(stbur)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_source.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_source.f90 2012-08-18 00:54:58 UTC (rev 20589)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/locate_source.f90 2012-08-18 17:02:11 UTC (rev 20590)
@@ -69,8 +69,6 @@
integer i,j,k,ispec,iglob,isource
integer imin,imax,jmin,jmax,kmin,kmax
-! integer igll,jgll,kgll,inode,iface,iglob_selected,
-! integer iselected,jselected,iface_selected,iadjust,jadjust
integer iproc(1)
double precision, dimension(NSOURCES) :: utm_x_source,utm_y_source
@@ -106,10 +104,8 @@
double precision,dimension(NPROC) :: distmin_ele_all,elevation_all
-! double precision,dimension(4) :: elevation_node,dist_node
real(kind=CUSTOM_REAL) :: xloc,yloc,loc_ele,loc_distmin
-
integer islice_selected_source(NSOURCES)
! timer MPI
@@ -221,81 +217,6 @@
altitude_source(1) = loc_ele
distmin_ele(1) = loc_distmin
-
-
-! ! set distance to huge initial value
-! distmin = HUGEVAL
-! if(num_free_surface_faces > 0) then
-! iglob_selected = 1
-! ! loop only on points inside the element
-! ! exclude edges to ensure this point is not shared with other elements
-! imin = 2
-! imax = NGLLX - 1
-!
-! jmin = 2
-! jmax = NGLLY - 1
-!
-! iselected = 0
-! jselected = 0
-! iface_selected = 0
-! do iface=1,num_free_surface_faces
-! do j=jmin,jmax
-! do i=imin,imax
-!
-! ispec = free_surface_ispec(iface)
-! igll = free_surface_ijk(1,(j-1)*NGLLY+i,iface)
-! jgll = free_surface_ijk(2,(j-1)*NGLLY+i,iface)
-! kgll = free_surface_ijk(3,(j-1)*NGLLY+i,iface)
-! iglob = ibool(igll,jgll,kgll,ispec)
-!
-! ! keep this point if it is closer to the receiver
-! dist = dsqrt((utm_x_source(isource)-dble(xstore(iglob)))**2 + &
-! (utm_y_source(isource)-dble(ystore(iglob)))**2)
-! if(dist < distmin) then
-! distmin = dist
-! iglob_selected = iglob
-! iface_selected = iface
-! iselected = i
-! jselected = j
-! altitude_source(1) = zstore(iglob_selected)
-! endif
-! enddo
-! enddo
-! ! end of loop on all the elements on the free surface
-! end do
-! ! weighted mean at current point of topography elevation of the four closest nodes
-! ! set distance to huge initial value
-! distmin = HUGEVAL
-! do j=jselected,jselected+1
-! do i=iselected,iselected+1
-! inode = 1
-! do jadjust=0,1
-! do iadjust= 0,1
-! ispec = free_surface_ispec(iface_selected)
-! igll = free_surface_ijk(1,(j-jadjust-1)*NGLLY+i-iadjust,iface_selected)
-! jgll = free_surface_ijk(2,(j-jadjust-1)*NGLLY+i-iadjust,iface_selected)
-! kgll = free_surface_ijk(3,(j-jadjust-1)*NGLLY+i-iadjust,iface_selected)
-! iglob = ibool(igll,jgll,kgll,ispec)
-!
-! elevation_node(inode) = zstore(iglob)
-! dist_node(inode) = dsqrt((utm_x_source(isource)-dble(xstore(iglob)))**2 + &
-! (utm_y_source(isource)-dble(ystore(iglob)))**2)
-! inode = inode + 1
-! end do
-! end do
-! dist = sum(dist_node)
-! if(dist < distmin) then
-! distmin = dist
-! altitude_source(1) = (dist_node(1)/dist)*elevation_node(1) + &
-! (dist_node(2)/dist)*elevation_node(2) + &
-! (dist_node(3)/dist)*elevation_node(3) + &
-! (dist_node(4)/dist)*elevation_node(4)
-! endif
-! end do
-! end do
-! end if
-! distmin_ele(1)= distmin
-
! MPI communications to determine the best slice
call gather_all_dp(distmin_ele,1,distmin_ele_all,1,NPROC)
call gather_all_dp(altitude_source,1,elevation_all,1,NPROC)
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_movie_meshes.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_movie_meshes.f90 2012-08-18 00:54:58 UTC (rev 20589)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_movie_meshes.f90 2012-08-18 17:02:11 UTC (rev 20590)
@@ -311,7 +311,3 @@
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 2012-08-18 00:54:58 UTC (rev 20589)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/setup_sources_receivers.f90 2012-08-18 17:02:11 UTC (rev 20590)
@@ -339,7 +339,7 @@
use specfem_par_acoustic
implicit none
- integer :: irec,isource,ier !,ios
+ integer :: irec,isource,ier
! reads in station file
if (SIMULATION_TYPE == 1) then
@@ -350,14 +350,6 @@
call station_filter(SUPPRESS_UTM_PROJECTION,UTM_PROJECTION_ZONE,myrank,rec_filename,filtered_rec_filename,nrec, &
LATITUDE_MIN, LATITUDE_MAX, LONGITUDE_MIN, LONGITUDE_MAX)
- ! get total number of stations
- !open(unit=IIN,file=rec_filename,iostat=ios,status='old',action='read')
- !nrec = 0
- !do while(ios == 0)
- ! read(IIN,"(a)",iostat=ios) dummystring
- ! if(ios == 0) nrec = nrec + 1
- !enddo
- !close(IIN)
if(nrec < 1) call exit_MPI(myrank,'need at least one receiver')
call sync_all()
@@ -528,10 +520,10 @@
real(kind=CUSTOM_REAL) :: factor_source
real(kind=CUSTOM_REAL) :: junk
integer :: isource,ispec
- integer :: irec !,irec_local
+ integer :: irec
integer :: i,j,k,iglob
integer :: icomp,itime,nadj_files_found,nadj_files_found_tot,ier
- character(len=3),dimension(NDIM) :: comp ! = (/ "BHE", "BHN", "BHZ" /)
+ character(len=3),dimension(NDIM) :: comp
character(len=256) :: filename
! forward simulations
@@ -620,7 +612,7 @@
else
! SIMULATION_TYPE == 2
! allocate dummy array (needed for subroutine calls)
- allocate(sourcearrays(0,0,0,0,0),stat=ier)
+ allocate(sourcearrays(1,1,1,1,1),stat=ier)
if( ier /= 0 ) stop 'error allocating dummy sourcearrays'
endif
@@ -806,10 +798,6 @@
real(kind=CUSTOM_REAL),dimension(NGNOD) :: xelm,yelm,zelm
integer :: ia,ispec,isource,irec,ier,totalpoints
- !INTEGER(kind=4) :: system_command_status
- !integer :: ret
- !integer,external :: system
-
character(len=256) :: filename,filename_new,system_command,system_command1,system_command2
! determines number of points for vtk file
@@ -949,40 +937,21 @@
"('awk ',a1,'{if(NR<5) print $0;if(NR==5)print ',a1,'POINTS',i6,' float',a1,';if(NR>5+',i6,')print $0}',a1,' < ',a,' > ',a)")&
"'",'"',nrec,'"',NSOURCES,"'",trim(filename),trim(filename_new)
-!daniel:
-! gfortran
-! call system(trim(system_command),system_command_status)
-! ifort
-! ret = system(trim(system_command))
-
! extracts source locations
- !"('awk ',a1,'{if(NR< 6 + ',i6,') print $0}END{print}',a1,' < ',a,' > ',a)")&
filename_new = trim(OUTPUT_FILES)//'/source.vtk'
write(system_command1, &
"('awk ',a1,'{if(NR<5) print $0;if(NR==5)print ',a1,'POINTS',i6,' float',a1,';')") &
"'",'"',NSOURCES,'"'
- !daniel
- !print*,'command 1:',trim(system_command1)
-
write(system_command2, &
"('if(NR>5 && NR <6+',i6,')print $0}END{print ',a,'}',a1,' < ',a,' > ',a)") &
NSOURCES,'" "',"'",trim(filename),trim(filename_new)
- !print*,'command 2:',trim(system_command2)
-
system_command = trim(system_command1)//trim(system_command2)
- !print*,'command:',trim(system_command)
-!daniel:
-! gfortran
-! call system(trim(system_command),system_command_status)
-! ifort
-! ret = system(trim(system_command))
-
endif
endif
-
end subroutine setup_sources_receivers_VTKfile
+
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90 2012-08-18 00:54:58 UTC (rev 20589)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/specfem3D_par.f90 2012-08-18 17:02:11 UTC (rev 20590)
@@ -65,8 +65,6 @@
! use integer array to store topography values
integer :: NX_TOPO,NY_TOPO
- !double precision :: ORIG_LAT_TOPO,ORIG_LONG_TOPO,DEGREES_PER_CELL_TOPO
- !character(len=100) :: topo_file
integer, dimension(:,:), allocatable :: itopo_bathy
! absorbing boundary arrays (for all boundaries) - keeps all infos, allowing for irregular surfaces
@@ -215,17 +213,10 @@
! maximum speed in velocity model
real(kind=CUSTOM_REAL):: model_speed_max
-!!!! NL NL REGOLITH : regolith layer for asteroid
-!!$ double precision, external :: materials_ext_mesh
-!!$ logical, dimension(:), allocatable :: ispec_is_regolith
-!!$ real(kind=CUSTOM_REAL) :: weight, jacobianl
-!!!! NL NL REGOLITH
-
! gravity
real(kind=CUSTOM_REAL), dimension(:),allocatable :: minus_deriv_gravity,minus_g
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: wgll_cube
-
! ADJOINT parameters
! time scheme
@@ -269,7 +260,6 @@
end module specfem_par
-
!=====================================================================
module specfem_par_elastic
@@ -505,9 +495,6 @@
C_kl, M_kl, rhob_kl, rhofb_kl, phi_kl, Bb_kl, Cb_kl, Mb_kl, mufrb_kl, &
rhobb_kl, rhofbb_kl, phib_kl, cpI_kl, cpII_kl, cs_kl, ratio_kl
- ! approximate hessian
- !real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: hess_kl
-
! absorbing stacey wavefield parts
real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: b_absorb_fields_poroelastic,b_absorb_fieldw_poroelastic
integer :: b_reclen_fields,b_reclen_fieldw
@@ -539,9 +526,6 @@
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
-
! shakemovies
real(kind=CUSTOM_REAL), dimension(:), allocatable :: store_val_x_external_mesh
real(kind=CUSTOM_REAL), dimension(:), allocatable :: store_val_y_external_mesh
@@ -584,3 +568,4 @@
logical :: MOVIE_SIMULATION
end module specfem_par_movie
+
Modified: seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_movie_output.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_movie_output.f90 2012-08-18 00:54:58 UTC (rev 20589)
+++ seismo/3D/SPECFEM3D/trunk/src/specfem3D/write_movie_output.f90 2012-08-18 17:02:11 UTC (rev 20590)
@@ -113,7 +113,6 @@
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
@@ -165,7 +164,6 @@
ibool,rhostore,GRAVITY)
endif
-
! high-resolution
if (USE_HIGHRES_FOR_MOVIES) then
do ipoin = 1, NGLLX*NGLLY
@@ -1012,9 +1010,6 @@
character(len=3) :: channel
character(len=1) :: compx,compy,compz
- !real(kind=CUSTOM_REAL),dimension(:,:,:,:),allocatable:: tmpdata
- !integer :: i,j,k,iglob
-
! 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
@@ -1132,14 +1127,8 @@
write(27) curl_z
close(27)
- !write(outputname,"('veloc_proc',i6.6,'_it',i6.6,'.bin')") myrank,it
- !open(unit=27,file=trim(LOCAL_PATH)//trim(outputname),status='unknown',form='unformatted')
- !write(27) veloc
- !close(27)
-
endif
-
if( ACOUSTIC_SIMULATION .or. ELASTIC_SIMULATION .or. POROELASTIC_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',iostat=ier)
@@ -1159,63 +1148,6 @@
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',iostat=ier)
- !if( ier /= 0 ) stop 'error opening file velocity_movie'
- !write(27) velocity_movie
- !close(27)
-
- ! norms
- !allocate( tmpdata(NGLLX,NGLLY,NGLLZ,NSPEC_AB),stat=ier)
- !if( ier /= 0 ) stop 'error allocating tmpdata arrays for movie output'
- !
- !if( ELASTIC_SIMULATION ) then
- ! ! norm of displacement
- ! do ispec=1,NSPEC_AB
- ! do k=1,NGLLZ
- ! do j=1,NGLLY
- ! do i=1,NGLLX
- ! iglob = ibool(i,j,k,ispec)
- ! tmpdata(i,j,k,ispec) = sqrt( displ(1,iglob)**2 + displ(2,iglob)**2 + displ(3,iglob)**2 )
- ! enddo
- ! enddo
- ! enddo
- ! enddo
- ! write(outputname,"('/proc',i6.6,'_displ_norm_it',i6.6,'.bin')") myrank,it
- ! 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) tmpdata
- ! close(27)
- !
- ! ! norm of acceleration
- ! do ispec=1,NSPEC_AB
- ! do k=1,NGLLZ
- ! do j=1,NGLLY
- ! do i=1,NGLLX
- ! iglob = ibool(i,j,k,ispec)
- ! tmpdata(i,j,k,ispec) = sqrt( accel(1,iglob)**2 + accel(2,iglob)**2 + accel(3,iglob)**2 )
- ! enddo
- ! enddo
- ! enddo
- ! enddo
- ! write(outputname,"('/proc',i6.6,'_accel_norm_it',i6.6,'.bin')") myrank,it
- ! 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) tmpdata
- ! close(27)
- !endif
- !
- !! norm of velocity
- !tmpdata = sqrt( velocity_x**2 + velocity_y**2 + velocity_z**2)
- !
- !write(outputname,"('/proc',i6.6,'_velocity_norm_it',i6.6,'.bin')") myrank,it
- !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) tmpdata
- !close(27)
- !
- !deallocate(tmpdata)
-
endif
end subroutine wmo_movie_volume_output
@@ -1230,7 +1162,6 @@
hprime_xx,hprime_yy,hprime_zz, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz)
-
! calculates div, curl and velocity
implicit none
@@ -1374,3 +1305,4 @@
enddo
end subroutine wmo_movie_div_curl
+
More information about the CIG-COMMITS
mailing list