[cig-commits] r17025 - seismo/3D/SPECFEM3D/trunk
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Mon Jun 28 10:49:55 PDT 2010
Author: danielpeter
Date: 2010-06-28 10:49:54 -0700 (Mon, 28 Jun 2010)
New Revision: 17025
Added:
seismo/3D/SPECFEM3D/trunk/model_interface_bedrock.f90
seismo/3D/SPECFEM3D/trunk/model_tomography.f90
Modified:
seismo/3D/SPECFEM3D/trunk/Makefile.in
seismo/3D/SPECFEM3D/trunk/create_regions_mesh.f90
seismo/3D/SPECFEM3D/trunk/get_model.f90
seismo/3D/SPECFEM3D/trunk/model_external_values.f90
Log:
adds Federica's tomography model routine in new file model_tomography.f90; puts some (old?) outcommented bedrock routines in new file model_interface_bedrock.f90
Modified: seismo/3D/SPECFEM3D/trunk/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D/trunk/Makefile.in 2010-06-26 01:11:50 UTC (rev 17024)
+++ seismo/3D/SPECFEM3D/trunk/Makefile.in 2010-06-28 17:49:54 UTC (rev 17025)
@@ -93,6 +93,7 @@
$O/generate_databases.o \
$O/model_external_values.o \
$O/model_aniso.o \
+ $O/model_tomography.o \
$O/netlib_specfun_erf.o \
$O/param_reader.o \
$O/prepare_assemble_MPI.o \
@@ -391,6 +392,8 @@
$O/model_external_values.o: constants.h model_external_values.f90
${MPIFCCOMPILE_CHECK} -c -o $O/model_external_values.o model_external_values.f90
+$O/model_tomography.o: constants.h model_tomography.f90
+ ${MPIFCCOMPILE_CHECK} -c -o $O/model_tomography.o model_tomography.f90
###
### serial compilation without optimization
Modified: seismo/3D/SPECFEM3D/trunk/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/create_regions_mesh.f90 2010-06-26 01:11:50 UTC (rev 17024)
+++ seismo/3D/SPECFEM3D/trunk/create_regions_mesh.f90 2010-06-28 17:49:54 UTC (rev 17025)
@@ -26,7 +26,7 @@
module create_regions_mesh_ext_par
include 'constants.h'
-
+
! global point coordinates
real(kind=CUSTOM_REAL), dimension(:), allocatable :: xstore_dummy
real(kind=CUSTOM_REAL), dimension(:), allocatable :: ystore_dummy
@@ -46,8 +46,8 @@
etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore,jacobianstore
! for model density, kappa, mu
- real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: rhostore,kappastore,mustore
-
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: rhostore,kappastore,mustore
+
! mass matrix
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass,rmass_acoustic,&
rmass_solid_poroelastic,rmass_fluid_poroelastic
@@ -56,7 +56,7 @@
integer :: NGLOB_OCEAN
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_ocean_load
-! attenuation
+! attenuation
integer, dimension(:,:,:,:), allocatable :: iflag_attenuation_store
! 2D shape functions and their derivatives, weights
@@ -101,7 +101,7 @@
! name of the database file
character(len=256) prname
-
+
end module create_regions_mesh_ext_par
!
@@ -133,7 +133,7 @@
! create the different regions of the mesh
- use create_regions_mesh_ext_par
+ use create_regions_mesh_ext_par
implicit none
!include "constants.h"
@@ -158,15 +158,15 @@
integer, dimension(ESIZE,nelmnts_ext_mesh) :: elmnts_ext_mesh
! static memory size needed by the solver
- double precision :: max_static_memory_size
-
+ double precision :: max_static_memory_size
+
integer, dimension(2,nelmnts_ext_mesh) :: mat_ext_mesh
! material properties
- integer :: nmat_ext_mesh,nundefMat_ext_mesh
- double precision, dimension(6,nmat_ext_mesh) :: materials_ext_mesh
+ integer :: nmat_ext_mesh,nundefMat_ext_mesh
+ double precision, dimension(6,nmat_ext_mesh) :: materials_ext_mesh
character (len=30), dimension(6,nundefMat_ext_mesh):: undef_mat_prop
-
+
! double precision, external :: materials_ext_mesh
! MPI communication
@@ -179,14 +179,14 @@
! absorbing boundaries
integer :: nspec2D_xmin, nspec2D_xmax, nspec2D_ymin, nspec2D_ymax, NSPEC2D_BOTTOM, NSPEC2D_TOP
- integer, dimension(nspec2D_xmin) :: ibelm_xmin
+ integer, dimension(nspec2D_xmin) :: ibelm_xmin
integer, dimension(nspec2D_xmax) :: ibelm_xmax
integer, dimension(nspec2D_ymin) :: ibelm_ymin
integer, dimension(nspec2D_ymax) :: ibelm_ymax
integer, dimension(NSPEC2D_BOTTOM) :: ibelm_bottom
integer, dimension(NSPEC2D_TOP) :: ibelm_top
! node indices of boundary faces
- integer, dimension(4,nspec2D_xmin) :: nodes_ibelm_xmin
+ integer, dimension(4,nspec2D_xmin) :: nodes_ibelm_xmin
integer, dimension(4,nspec2D_xmax) :: nodes_ibelm_xmax
integer, dimension(4,nspec2D_ymin) :: nodes_ibelm_ymin
integer, dimension(4,nspec2D_ymax) :: nodes_ibelm_ymax
@@ -205,12 +205,12 @@
integer :: NX_TOPO,NY_TOPO
double precision :: ORIG_LAT_TOPO,ORIG_LONG_TOPO,DEGREES_PER_CELL_TOPO
integer, dimension(NX_TOPO,NY_TOPO) :: itopo_bathy
-
+
! local parameters
! static memory size needed by the solver
double precision :: static_memory_size
real(kind=CUSTOM_REAL) :: model_speed_max
-
+
! for vtk output
! character(len=256) prname_file
! integer,dimension(:),allocatable :: itest_flag
@@ -233,32 +233,32 @@
! ! store bedrock values
! integer :: icornerlat,icornerlong
! double precision :: lat,long,elevation_bedrock
-! double precision :: lat_corner,long_corner,ratio_xi,ratio_eta
+! double precision :: lat_corner,long_corner,ratio_xi,ratio_eta
!real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: ibedrock
! initializes arrays
call sync_all()
if( myrank == 0) then
- write(IMAIN,*)
+ write(IMAIN,*)
write(IMAIN,*) ' ...allocating arrays '
endif
call crm_ext_allocate_arrays(nspec,LOCAL_PATH,myrank, &
nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax, &
nspec2D_bottom,nspec2D_top,ANISOTROPY)
-
+
! fills location and weights for Gauss-Lobatto-Legendre points, shape and derivations,
! returns jacobianstore,xixstore,...gammazstore
! and GLL-point locations in xstore,ystore,zstore
call sync_all()
if( myrank == 0) then
write(IMAIN,*) ' ...setting up jacobian '
- endif
- call crm_ext_setup_jacobian(myrank, &
+ endif
+ call crm_ext_setup_jacobian(myrank, &
xstore,ystore,zstore,nspec, &
nodes_coords_ext_mesh,nnodes_ext_mesh,&
elmnts_ext_mesh,nelmnts_ext_mesh)
-
+
! creates ibool index array for projection from local to global points
call sync_all()
if( myrank == 0) then
@@ -286,12 +286,12 @@
if( myrank == 0) then
write(IMAIN,*) ' ...determining velocity model'
endif
- call get_model(nspec,ibool,mat_ext_mesh,nelmnts_ext_mesh, &
+ call get_model(myrank,nspec,ibool,mat_ext_mesh,nelmnts_ext_mesh, &
materials_ext_mesh,nmat_ext_mesh, &
undef_mat_prop,nundefMat_ext_mesh, &
ANISOTROPY)
-
-! sets up absorbing/free surface boundaries
+
+! sets up absorbing/free surface boundaries
call sync_all()
if( myrank == 0) then
write(IMAIN,*) ' ...setting up absorbing boundaries '
@@ -303,7 +303,7 @@
nodes_ibelm_bottom,nodes_ibelm_top, &
nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax, &
nspec2D_bottom,nspec2D_top)
-
+
! sets up acoustic-elastic coupling surfaces
call sync_all()
if( myrank == 0) then
@@ -315,14 +315,14 @@
num_interfaces_ext_mesh,max_interface_size_ext_mesh, &
my_neighbours_ext_mesh)
-! creates mass matrix
+! creates mass matrix
call sync_all()
if( myrank == 0) then
write(IMAIN,*) ' ...creating mass matrix '
endif
call create_mass_matrices(nglob,nspec,ibool)
-! creates ocean load mass matrix
+! creates ocean load mass matrix
call sync_all()
if( myrank == 0) then
write(IMAIN,*) ' ...creating ocean load mass matrix '
@@ -350,7 +350,7 @@
free_surface_normal,free_surface_jacobian2Dw, &
free_surface_ijk,free_surface_ispec,num_free_surface_faces, &
coupling_ac_el_normal,coupling_ac_el_jacobian2Dw, &
- coupling_ac_el_ijk,coupling_ac_el_ispec,num_coupling_ac_el_faces, &
+ coupling_ac_el_ijk,coupling_ac_el_ispec,num_coupling_ac_el_faces, &
num_interfaces_ext_mesh,my_neighbours_ext_mesh,nibool_interfaces_ext_mesh, &
max_interface_size_ext_mesh,ibool_interfaces_ext_mesh, &
prname,SAVE_MESH_FILES,ANISOTROPY,NSPEC_ANISO, &
@@ -364,7 +364,7 @@
call memory_eval(nspec,nglob,maxval(nibool_interfaces_ext_mesh),num_interfaces_ext_mesh,static_memory_size)
call max_all_dp(static_memory_size, max_static_memory_size)
-! checks the mesh, stability and resolved period
+! checks the mesh, stability and resolved period
call sync_all()
call check_mesh_resolution(myrank,nspec,nglob,ibool,&
xstore_dummy,ystore_dummy,zstore_dummy, &
@@ -373,7 +373,7 @@
! VTK file output
! if( SAVE_MESH_FILES ) then
-! ! saves material flag assigned for each spectral element into a vtk file
+! ! saves material flag assigned for each spectral element into a vtk file
! prname_file = prname(1:len_trim(prname))//'material_flag'
! allocate(elem_flag(nspec))
! elem_flag(:) = mat_ext_mesh(1,:)
@@ -381,7 +381,7 @@
! xstore_dummy,ystore_dummy,zstore_dummy,ibool, &
! elem_flag,prname_file)
! deallocate(elem_flag)
-!
+!
! !plotting abs boundaries
! ! allocate(itest_flag(nspec))
! ! itest_flag(:) = 0
@@ -393,7 +393,7 @@
! ! xstore_dummy,ystore_dummy,zstore_dummy,ibool, &
! ! itest_flag,prname_file)
! ! deallocate(itest_flag)
-! endif
+! endif
! AVS/DX file output
! create AVS or DX mesh data for the slice, edges and faces
@@ -415,7 +415,7 @@
deallocate(xixstore,xiystore,xizstore,&
etaxstore,etaystore,etazstore,&
gammaxstore,gammaystore,gammazstore)
- deallocate(jacobianstore,iflag_attenuation_store)
+ deallocate(jacobianstore,iflag_attenuation_store)
deallocate(kappastore,mustore,rho_vp,rho_vs)
end subroutine create_regions_mesh_ext
@@ -428,7 +428,7 @@
nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax, &
nspec2D_bottom,nspec2D_top,ANISOTROPY)
- use create_regions_mesh_ext_par
+ use create_regions_mesh_ext_par
implicit none
integer :: nspec,myrank
@@ -436,15 +436,15 @@
nspec2D_bottom,nspec2D_top
character(len=256) :: LOCAL_PATH
-
+
logical :: ANISOTROPY
-! local parameters
+! local parameters
integer :: ier
! memory test
-! logical,dimension(:),allocatable :: test_mem
-!
+! 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
@@ -453,10 +453,10 @@
! call exit_MPI(myrank,'not enough memory to allocate arrays')
! endif
! test_mem(:) = .true.
-! deallocate( test_mem, stat=ier)
+! 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)
allocate( iflag_attenuation_store(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
@@ -480,15 +480,15 @@
shape2D_y(NGNOD2D,NGLLX,NGLLZ), &
shape2D_bottom(NGNOD2D,NGLLX,NGLLY), &
shape2D_top(NGNOD2D,NGLLX,NGLLY), stat=ier)
-
+
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)
-
+
allocate(wgllwgll_xy(NGLLX,NGLLY), &
wgllwgll_xz(NGLLX,NGLLZ), &
- wgllwgll_yz(NGLLY,NGLLZ),stat=ier)
+ wgllwgll_yz(NGLLY,NGLLZ),stat=ier)
if(ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays')
! Stacey
@@ -499,9 +499,9 @@
! array with model density
allocate(rhostore(NGLLX,NGLLY,NGLLZ,nspec), &
kappastore(NGLLX,NGLLY,NGLLZ,nspec), &
- mustore(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
+ mustore(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
!vpstore(NGLLX,NGLLY,NGLLZ,nspec), &
- !vsstore(NGLLX,NGLLY,NGLLZ,nspec),
+ !vsstore(NGLLX,NGLLY,NGLLZ,nspec),
if(ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays')
! arrays with mesh parameters
@@ -517,7 +517,7 @@
jacobianstore(NGLLX,NGLLY,NGLLZ,nspec),stat=ier)
if(ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays')
-! absorbing boundary
+! absorbing boundary
! absorbing faces
num_abs_boundary_faces = nspec2D_xmin + nspec2D_xmax + nspec2D_ymin + nspec2D_ymax + nspec2D_bottom
! adds faces of free surface if it also absorbs
@@ -536,7 +536,7 @@
! no free surface - uses a dummy size
num_free_surface_faces = 1
else
- num_free_surface_faces = nspec2D_top
+ num_free_surface_faces = nspec2D_top
endif
! allocates arrays to store info for each face (assumes NGLLX=NGLLY=NGLLZ)
@@ -580,7 +580,7 @@
ispec_is_elastic(nspec), &
ispec_is_poroelastic(nspec), stat=ier)
if(ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays')
-
+
end subroutine crm_ext_allocate_arrays
@@ -588,7 +588,7 @@
!-------------------------------------------------------------------------------------------------
!
-subroutine crm_ext_setup_jacobian(myrank, &
+subroutine crm_ext_setup_jacobian(myrank, &
xstore,ystore,zstore,nspec, &
nodes_coords_ext_mesh,nnodes_ext_mesh,&
elmnts_ext_mesh,nelmnts_ext_mesh)
@@ -684,7 +684,7 @@
! creates global indexing array ibool
- use create_regions_mesh_ext_par
+ use create_regions_mesh_ext_par
implicit none
! number of spectral elements in each block
@@ -699,7 +699,7 @@
double precision, dimension(NDIM,nnodes_ext_mesh) :: nodes_coords_ext_mesh
! local parameters
-! variables for creating array ibool
+! variables for creating array ibool
double precision, dimension(:), allocatable :: xp,yp,zp
integer, dimension(:), allocatable :: locval
logical, dimension(:), allocatable :: ifseg
@@ -754,8 +754,8 @@
! unique global point locations
allocate(xstore_dummy(nglob), &
ystore_dummy(nglob), &
- zstore_dummy(nglob),stat=ier)
- if(ier /= 0) stop 'error in allocate'
+ zstore_dummy(nglob),stat=ier)
+ if(ier /= 0) stop 'error in allocate'
do ispec = 1, nspec
do k = 1, NGLLZ
do j = 1, NGLLY
@@ -767,153 +767,6 @@
enddo
enddo
enddo
- enddo
+ enddo
end subroutine crm_ext_setup_indexing
-
-!
-!-------------------------------------------------------------------------------------------------
-!
-
-!
-!
-!subroutine bubble_sort( arr, ndim )
-!
-!! sorts values in array arr[ndim] in increasing order
-!
-! implicit none
-!
-! integer :: ndim
-! integer :: arr(ndim)
-!
-! logical :: swapped
-! integer :: j,tmp
-!
-! swapped = .true.
-! do while( swapped )
-! swapped = .false.
-! do j = 1, ndim-1
-! if( arr(j+1) < arr(j) ) then
-! tmp = arr(j)
-! arr(j) = arr(j+1)
-! arr(j+1) = tmp
-! swapped = .true.
-! endif
-! enddo
-! enddo
-!
-!end subroutine bubble_sort
-
-!
-!-------------------------------------------------------------------------------------------------
-!
-
-
-!
-! subroutine interface(iflag,flag_below,flag_above,ispec,nspec,i,j,k,xstore,ystore,zstore,ibedrock)
-
-! implicit none
-
-! include "constants.h"
-
-! integer :: iflag,flag_below,flag_above
-! integer :: ispec,nspec
-! integer :: i,j,k
-! double precision, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: xstore,ystore,zstore
-! real(kind=CUSTOM_REAL), dimension(NX_TOPO_ANT,NY_TOPO_ANT) :: ibedrock
-! integer, parameter :: NUMBER_OF_STATIONS = 1
-! double precision, parameter :: RADIUS_TO_EXCLUDE = 250.d0
-! double precision, dimension(NUMBER_OF_STATIONS) :: utm_x_station,utm_y_station
-
-! !-------------------
-
-! !for Piero
-! logical :: is_around_a_station
-! integer :: istation
-
-! ! store bedrock values
-! integer :: icornerlat,icornerlong
-! double precision :: lat,long,elevation_bedrock
-! double precision :: lat_corner,long_corner,ratio_xi,ratio_eta
-
-
-! !! DK DK store the position of the six stations to be able to
-! !! DK DK exclude circles around each station to make sure they are on the bedrock
-! !! DK DK and not in the ice
-! utm_x_station(1) = 783500.6250000d0
-! utm_y_station(1) = -11828.7519531d0
-
-! utm_x_station(2) = 853644.5000000d0
-! utm_y_station(2) = -114.0138092d0
-
-! utm_x_station(3) = 863406.0000000d0
-! utm_y_station(3) = -53736.1640625d0
-
-! utm_x_station(4) = 823398.8125000d0
-! utm_y_station(4) = 29847.4511719d0
-
-! utm_x_station(5) = 863545.3750000d0
-! utm_y_station(5) = 19669.6621094d0
-
-! utm_x_station(6) = 817099.3750000d0
-! utm_y_station(6) = -24430.2871094d0
-
-! ! since we have suppressed UTM projection for Piero Basini, UTMx is the same as long
-! ! and UTMy is the same as lat
-! long = xstore(i,j,k,ispec)
-! lat = ystore(i,j,k,ispec)
-
-! ! get coordinate of corner in model
-! icornerlong = int((long - ORIG_LONG_TOPO_ANT) / DEGREES_PER_CELL_TOPO_ANT) + 1
-! icornerlat = int((lat - ORIG_LAT_TOPO_ANT) / DEGREES_PER_CELL_TOPO_ANT) + 1
-
-! ! avoid edge effects and extend with identical point if outside model
-! if(icornerlong < 1) icornerlong = 1
-! if(icornerlong > NX_TOPO_ANT-1) icornerlong = NX_TOPO_ANT-1
-! if(icornerlat < 1) icornerlat = 1
-! if(icornerlat > NY_TOPO_ANT-1) icornerlat = NY_TOPO_ANT-1
-
-! ! compute coordinates of corner
-! long_corner = ORIG_LONG_TOPO_ANT + (icornerlong-1)*DEGREES_PER_CELL_TOPO_ANT
-! lat_corner = ORIG_LAT_TOPO_ANT + (icornerlat-1)*DEGREES_PER_CELL_TOPO_ANT
-
-! ! compute ratio for interpolation
-! ratio_xi = (long - long_corner) / DEGREES_PER_CELL_TOPO_ANT
-! ratio_eta = (lat - lat_corner) / DEGREES_PER_CELL_TOPO_ANT
-
-! ! avoid edge effects
-! if(ratio_xi < 0.) ratio_xi = 0.
-! if(ratio_xi > 1.) ratio_xi = 1.
-! if(ratio_eta < 0.) ratio_eta = 0.
-! if(ratio_eta > 1.) ratio_eta = 1.
-
-! ! interpolate elevation at current point
-! elevation_bedrock = &
-! ibedrock(icornerlong,icornerlat)*(1.-ratio_xi)*(1.-ratio_eta) + &
-! ibedrock(icornerlong+1,icornerlat)*ratio_xi*(1.-ratio_eta) + &
-! ibedrock(icornerlong+1,icornerlat+1)*ratio_xi*ratio_eta + &
-! ibedrock(icornerlong,icornerlat+1)*(1.-ratio_xi)*ratio_eta
-
-! !! DK DK exclude circles around each station to make sure they are on the bedrock
-! !! DK DK and not in the ice
-! is_around_a_station = .false.
-! do istation = 1,NUMBER_OF_STATIONS
-! if(sqrt((xstore(i,j,k,ispec) - utm_x_station(istation))**2 + (ystore(i,j,k,ispec) - &
-! utm_y_station(istation))**2) < RADIUS_TO_EXCLUDE) then
-! is_around_a_station = .true.
-! exit
-! endif
-! enddo
-
-! ! we are above the bedrock interface i.e. in the ice, and not too close to a station
-! if(zstore(i,j,k,ispec) >= elevation_bedrock .and. .not. is_around_a_station) then
-! iflag = flag_above
-! !iflag_attenuation_store(i,j,k,ispec) = IATTENUATION_ICE
-! ! we are below the bedrock interface i.e. in the bedrock, or close to a station
-! else
-! iflag = flag_below
-! !iflag_attenuation_store(i,j,k,ispec) = IATTENUATION_BEDROCK
-! endif
-
-
-! end subroutine interface
Modified: seismo/3D/SPECFEM3D/trunk/get_model.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/get_model.f90 2010-06-26 01:11:50 UTC (rev 17024)
+++ seismo/3D/SPECFEM3D/trunk/get_model.f90 2010-06-28 17:49:54 UTC (rev 17025)
@@ -23,103 +23,112 @@
!
!=====================================================================
- subroutine get_model(nspec,ibool,mat_ext_mesh,nelmnts_ext_mesh, &
+ subroutine get_model(myrank,nspec,ibool,mat_ext_mesh,nelmnts_ext_mesh, &
materials_ext_mesh,nmat_ext_mesh, &
undef_mat_prop,nundefMat_ext_mesh, &
ANISOTROPY)
- use create_regions_mesh_ext_par
+ use create_regions_mesh_ext_par
implicit none
! number of spectral elements in each block
- integer :: nspec
+ integer :: myrank,nspec
integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool
! external mesh
integer :: nelmnts_ext_mesh
- integer :: nmat_ext_mesh,nundefMat_ext_mesh
+ integer :: nmat_ext_mesh,nundefMat_ext_mesh
integer, dimension(2,nelmnts_ext_mesh) :: mat_ext_mesh
- double precision, dimension(6,nmat_ext_mesh) :: materials_ext_mesh
+ double precision, dimension(6,nmat_ext_mesh) :: materials_ext_mesh
character (len=30), dimension(6,nundefMat_ext_mesh):: undef_mat_prop
! anisotropy
logical :: ANISOTROPY
! local parameters
- real(kind=CUSTOM_REAL) :: vp,vs,rho
+ real(kind=CUSTOM_REAL) :: vp,vs,rho
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,iflag_atten
integer :: iflag,flag_below,flag_above
integer :: iflag_aniso,idomain_id,imaterial_id
-
-! ! Piero, read bedrock file
-! allocate(ibedrock(NX_TOPO_ANT,NY_TOPO_ANT))
-! if(myrank == 0) then
-! call read_bedrock_file(ibedrock)
-! ! write(IMAIN,*)
-! ! write(IMAIN,*) 'regional bedrock file read ranges in m from ',minval(ibedrock),' to ',maxval(ibedrock)
-! ! write(IMAIN,*)
-! endif
-! ! broadcast the information read on the master to the nodes
-! ! call MPI_BCAST(ibedrock,NX_TOPO_ANT*NY_TOPO_ANT,MPI_REAL,0,MPI_COMM_WORLD,ier)
-! call bcast_all_cr(ibedrock,NX_TOPO_ANT*NY_TOPO_ANT)
+ ! gll point location
+ double precision :: xloc,yloc,zloc
+ integer :: iglob
+
+ ! initializes element domain flags
ispec_is_acoustic(:) = .false.
ispec_is_elastic(:) = .false.
ispec_is_poroelastic(:) = .false.
- ! material properties on all GLL points: taken from material values defined for
+ ! prepares tomography model if needed for elements with undefined material definitions
+ if( nundefMat_ext_mesh > 0 ) then
+ call model_tomography_broadcast(myrank)
+ endif
+
+ ! prepares external model values if needed
+ if( USE_MODEL_EXTERNAL_VALUES ) then
+ call model_external_broadcast(myrank)
+ endif
+
+! ! Piero, read bedrock file
+! in case, see file model_interface_bedrock.f90:
+! call model_bedrock_broadcast(myrank)
+
+
+ ! material properties on all GLL points: taken from material values defined for
! each spectral element in input mesh
do ispec = 1, nspec
do k = 1, NGLLZ
do j = 1, NGLLY
do i = 1, NGLLX
-
+
! material index 1: associated material number
imaterial_id = mat_ext_mesh(1,ispec)
-
+
! check if the material is known or unknown
if( imaterial_id > 0) then
! gets velocity model as specified by (cubit) mesh files
-
- ! density
- ! materials_ext_mesh format:
- ! #index1 = rho #index2 = vp #index3 = vs #index4 = Q_flag #index5 = 0
+
+ ! density
+ ! materials_ext_mesh format:
+ ! #index1 = rho #index2 = vp #index3 = vs #index4 = Q_flag #index5 = 0
rho = materials_ext_mesh(1,imaterial_id)
-
- ! isotropic values: vp, vs
+
+ ! isotropic values: vp, vs
vp = materials_ext_mesh(2,imaterial_id)
vs = materials_ext_mesh(3,imaterial_id)
! attenuation
- iflag_atten = materials_ext_mesh(4,imaterial_id)
+ iflag_atten = materials_ext_mesh(4,imaterial_id)
!change for piero :
!if(mat_ext_mesh(1,ispec) == 1) then
! iflag_attenuation_store(i,j,k,ispec) = 1
!else
! iflag_attenuation_store(i,j,k,ispec) = 2
!endif
-
+
! anisotropy
iflag_aniso = materials_ext_mesh(5,imaterial_id)
-
+
! material domain_id
idomain_id = materials_ext_mesh(6,imaterial_id)
-
+
else if (mat_ext_mesh(2,ispec) == 1) then
-
+
stop 'material: interface not implemented yet'
-
- do iundef = 1,nundefMat_ext_mesh
+
+ do iundef = 1,nundefMat_ext_mesh
if(trim(undef_mat_prop(2,iundef)) == 'interface') then
read(undef_mat_prop(4,iundef),'(1i3)') flag_below
read(undef_mat_prop(5,iundef),'(1i3)') flag_above
endif
enddo
+ ! see file model_interface_bedrock.f90: routine interface()
!call interface(iflag,flag_below,flag_above,ispec,nspec,i,j,k,xstore,ystore,zstore,ibedrock)
iflag = 1
@@ -135,10 +144,27 @@
! endif
iflag_aniso = materials_ext_mesh(5,iflag)
idomain_id = materials_ext_mesh(6,iflag)
+
+ else if ( imaterial_id < 0 ) then
+
+ ! material definition undefined, uses definition from tomography model
+ ! GLL point location
+ iglob = ibool(i,j,k,ispec)
+ xloc = xstore_dummy(iglob)
+ yloc = ystore_dummy(iglob)
+ zloc = zstore_dummy(iglob)
+
+ ! gets model values from tomography file
+ call model_tomography(xloc,yloc,zloc, &
+ rho,vp,vs)
+
+ iflag_atten = 1 ! attenuation: would use IATTENUATION_SEDIMENTS_40
+ iflag_aniso = 0 ! no anisotropy
+ idomain_id = 2 ! elastic domain
+
else
- stop 'material: tomography not implemented yet'
- ! call tomography()
+ stop 'material: not implemented yet'
end if
@@ -153,28 +179,28 @@
c34,c35,c36,c44,c45,c46, &
c55,c56,c66,ANISOTROPY)
endif
-
+
! adds anisotropic default model
if( ANISOTROPY .and. .not. USE_MODEL_EXTERNAL_VALUES ) then
call model_aniso(iflag_aniso,rho,vp,vs,c11,c12,c13,c14,c15,c16, &
c22,c23,c24,c25,c26,c33,c34,c35,c36,c44,c45, &
- c46,c55,c56,c66)
-
+ c46,c55,c56,c66)
+
endif
-
+
! stores velocity model
-
+
! density
rhostore(i,j,k,ispec) = rho
-
+
! kappa, mu
- kappastore(i,j,k,ispec) = rho*( vp*vp - FOUR_THIRDS*vs*vs )
+ kappastore(i,j,k,ispec) = rho*( vp*vp - FOUR_THIRDS*vs*vs )
mustore(i,j,k,ispec) = rho*vs*vs
! attenuation
iflag_attenuation_store(i,j,k,ispec) = iflag_atten
-
- ! Stacey, a completer par la suite
+
+ ! Stacey, a completer par la suite
rho_vp(i,j,k,ispec) = rho*vp
rho_vs(i,j,k,ispec) = rho*vs
!end pll
@@ -182,7 +208,7 @@
! adds anisotropic perturbation to vp, vs
if( ANISOTROPY ) then
!call model_aniso(iflag_aniso,rho,vp,vs,c11,c12,c13,c14,c15,c16, &
- ! c22,c23,c24,c25,c26,c33,c34,c35,c36,c44,c45,c46,c55,c56,c66)
+ ! c22,c23,c24,c25,c26,c33,c34,c35,c36,c44,c45,c46,c55,c56,c66)
c11store(i,j,k,ispec) = c11
c12store(i,j,k,ispec) = c12
c13store(i,j,k,ispec) = c13
@@ -203,13 +229,13 @@
c46store(i,j,k,ispec) = c46
c55store(i,j,k,ispec) = c55
c56store(i,j,k,ispec) = c56
- c66store(i,j,k,ispec) = c66
+ c66store(i,j,k,ispec) = c66
endif
! material domain
- !print*,'velocity model:',ispec,idomain_id
+ !print*,'velocity model:',ispec,idomain_id
if( idomain_id == IDOMAIN_ACOUSTIC ) then
- ispec_is_acoustic(ispec) = .true.
+ ispec_is_acoustic(ispec) = .true.
else if( idomain_id == IDOMAIN_ELASTIC ) then
ispec_is_elastic(ispec) = .true.
else if( idomain_id == IDOMAIN_POROELASTIC ) then
@@ -218,11 +244,11 @@
else
stop 'error material domain index'
endif
-
+
enddo
enddo
enddo
- !print*,myrank,'ispec:',ispec,'rho:',rhostore(1,1,1,ispec),'vp:',vpstore(1,1,1,ispec),'vs:',vsstore(1,1,1,ispec)
+ !print*,myrank,'ispec:',ispec,'rho:',rhostore(1,1,1,ispec),'vp:',vpstore(1,1,1,ispec),'vs:',vsstore(1,1,1,ispec)
enddo
! checks material domains
@@ -233,7 +259,7 @@
print*,'error material domain not assigned to element:',ispec
print*,'acoustic: ',ispec_is_acoustic(ispec)
print*,'elastic: ',ispec_is_elastic(ispec)
- print*,'poroelastic: ',ispec_is_poroelastic(ispec)
+ print*,'poroelastic: ',ispec_is_poroelastic(ispec)
stop 'error material domain index element'
endif
enddo
@@ -242,170 +268,7 @@
! !! DK DK store the position of the six stations to be able to
! !! DK DK exclude circles around each station to make sure they are on the bedrock
! !! DK DK and not in the ice
-! utm_x_station(1) = 783500.6250000d0
-! utm_y_station(1) = -11828.7519531d0
+! in case, see file model_interface_bedrock.f90: routine model_bedrock_store()
-! utm_x_station(2) = 853644.5000000d0
-! utm_y_station(2) = -114.0138092d0
-
-! utm_x_station(3) = 863406.0000000d0
-! utm_y_station(3) = -53736.1640625d0
-
-! utm_x_station(4) = 823398.8125000d0
-! utm_y_station(4) = 29847.4511719d0
-
-! utm_x_station(5) = 863545.3750000d0
-! utm_y_station(5) = 19669.6621094d0
-
-! utm_x_station(6) = 817099.3750000d0
-! utm_y_station(6) = -24430.2871094d0
-
-! print*,myrank,'après store the position of the six stations'
-! call flush(6)
-
-! print*, myrank,minval(nodes_coords_ext_mesh(1,:))
-! call flush(6)
-
-
-! print*, myrank,maxval(nodes_coords_ext_mesh(1,:))
-! call flush(6)
-
-
-! do ispec = 1, nspec
-
-! zmesh = zstore(2,2,2,ispec)
-
-! ! if(doubling_index == IFLAG_ONE_LAYER_TOPOGRAPHY) then
-! if(any(ibelm_top == ispec)) then
-! doubling_value_found_for_Piero = IFLAG_ONE_LAYER_TOPOGRAPHY
-
-! else if(zmesh < Z_23p4km) then
-! doubling_value_found_for_Piero = IFLAG_MANTLE_BELOW_23p4km
-
-! else if(zmesh < Z_14km) then
-! doubling_value_found_for_Piero = IFLAG_14km_to_23p4km
-
-! else
-! doubling_value_found_for_Piero = IFLAG_BEDROCK_down_to_14km
-! endif
-! idoubling(ispec) = doubling_value_found_for_Piero
-
-! do k = 1, NGLLZ
-! do j = 1, NGLLY
-! do i = 1, NGLLX
-
-
-! if(idoubling(ispec) == IFLAG_ONE_LAYER_TOPOGRAPHY .or. &
-! idoubling(ispec) == IFLAG_BEDROCK_down_to_14km) then
-
-! ! since we have suppressed UTM projection for Piero Basini, UTMx is the same as long
-! ! and UTMy is the same as lat
-! long = xstore(i,j,k,ispec)
-! lat = ystore(i,j,k,ispec)
-
-! ! get coordinate of corner in model
-! icornerlong = int((long - ORIG_LONG_TOPO) / DEGREES_PER_CELL_TOPO) + 1
-! icornerlat = int((lat - ORIG_LAT_TOPO) / DEGREES_PER_CELL_TOPO) + 1
-
-! ! avoid edge effects and extend with identical point if outside model
-! if(icornerlong < 1) icornerlong = 1
-! if(icornerlong > NX_TOPO-1) icornerlong = NX_TOPO-1
-! if(icornerlat < 1) icornerlat = 1
-! if(icornerlat > NY_TOPO-1) icornerlat = NY_TOPO-1
-
-! ! compute coordinates of corner
-! long_corner = ORIG_LONG_TOPO + (icornerlong-1)*DEGREES_PER_CELL_TOPO
-! lat_corner = ORIG_LAT_TOPO + (icornerlat-1)*DEGREES_PER_CELL_TOPO
-
-! ! compute ratio for interpolation
-! ratio_xi = (long - long_corner) / DEGREES_PER_CELL_TOPO
-! ratio_eta = (lat - lat_corner) / DEGREES_PER_CELL_TOPO
-
-! ! avoid edge effects
-! if(ratio_xi < 0.) ratio_xi = 0.
-! if(ratio_xi > 1.) ratio_xi = 1.
-! if(ratio_eta < 0.) ratio_eta = 0.
-! if(ratio_eta > 1.) ratio_eta = 1.
-
-! ! interpolate elevation at current point
-! elevation_bedrock = &
-! ibedrock(icornerlong,icornerlat)*(1.-ratio_xi)*(1.-ratio_eta) + &
-! ibedrock(icornerlong+1,icornerlat)*ratio_xi*(1.-ratio_eta) + &
-! ibedrock(icornerlong+1,icornerlat+1)*ratio_xi*ratio_eta + &
-! ibedrock(icornerlong,icornerlat+1)*(1.-ratio_xi)*ratio_eta
-
-! !! DK DK exclude circles around each station to make sure they are on the bedrock
-! !! DK DK and not in the ice
-! is_around_a_station = .false.
-! do istation = 1,NUMBER_OF_STATIONS
-! if(sqrt((long - utm_x_station(istation))**2 + (lat - utm_y_station(istation))**2) < RADIUS_TO_EXCLUDE) then
-! is_around_a_station = .true.
-! exit
-! endif
-! enddo
-
-! ! define elastic parameters in the model
-
-! ! we are above the bedrock interface i.e. in the ice, and not too close to a station
-! if(zmesh >= elevation_bedrock .and. .not. is_around_a_station) then
-! vp = 3800.d0
-! vs = 1900.d0
-! rho = 900.d0
-! iflag_attenuation_store(i,j,k,ispec) = IATTENUATION_ICE
-
-! ! we are below the bedrock interface i.e. in the bedrock, or close to a station
-! else
-! vp = 5800.d0
-! vs = 3200.d0
-! rho = 2600.d0
-! iflag_attenuation_store(i,j,k,ispec) = IATTENUATION_BEDROCK
-! endif
-
-! else if(idoubling(ispec) == IFLAG_14km_to_23p4km) then
-! vp = 6800.d0
-! vs = 3900.d0
-! rho = 2900.d0
-! iflag_attenuation_store(i,j,k,ispec) = IATTENUATION_BEDROCK
-
-! else if(idoubling(ispec) == IFLAG_MANTLE_BELOW_23p4km) then
-! vp = 8100.d0
-! vs = 4480.d0
-! rho = 3380.d0
-! iflag_attenuation_store(i,j,k,ispec) = IATTENUATION_BEDROCK
-
-! endif
-
-! !pll 8/06
-! if(CUSTOM_REAL == SIZE_REAL) then
-! rhostore(i,j,k,ispec) = sngl(rho)
-! vpstore(i,j,k,ispec) = sngl(vp)
-! vsstore(i,j,k,ispec) = sngl(vs)
-! else
-! rhostore(i,j,k,ispec) = rho
-! vpstore(i,j,k,ispec) = vp
-! vsstore(i,j,k,ispec) = vs
-! end if
-
-! kappastore(i,j,k,ispec) = rhostore(i,j,k,ispec)*(vpstore(i,j,k,ispec)*vpstore(i,j,k,ispec) - &
-! 4.d0*vsstore(i,j,k,ispec)*vsstore(i,j,k,ispec)/3.d0)
-! mustore(i,j,k,ispec) = rhostore(i,j,k,ispec)*vsstore(i,j,k,ispec)*&
-! vsstore(i,j,k,ispec)
-
-! ! Stacey, a completer par la suite
-! rho_vp(i,j,k,ispec) = rhostore(i,j,k,ispec)*vpstore(i,j,k,ispec)
-! rho_vs(i,j,k,ispec) = rhostore(i,j,k,ispec)*vsstore(i,j,k,ispec)
-! !end pll
-
-! ! kappastore(i,j,k,ispec) = materials_ext_mesh(1,mat_ext_mesh(ispec))* &
-! ! (materials_ext_mesh(2,mat_ext_mesh(ispec))*materials_ext_mesh(2,mat_ext_mesh(ispec)) - &
-! ! 4.d0*materials_ext_mesh(3,mat_ext_mesh(ispec))*materials_ext_mesh(3,mat_ext_mesh(ispec))/3.d0)
-! ! mustore(i,j,k,ispec) = materials_ext_mesh(1,mat_ext_mesh(ispec))* &
-! materials_ext_mesh(3,mat_ext_mesh(ispec))*&
-! ! x materials_ext_mesh(3,mat_ext_mesh(ispec))
-! enddo
-! enddo
-! enddo
-! enddo
-
end subroutine get_model
Modified: seismo/3D/SPECFEM3D/trunk/model_external_values.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/model_external_values.f90 2010-06-26 01:11:50 UTC (rev 17024)
+++ seismo/3D/SPECFEM3D/trunk/model_external_values.f90 2010-06-28 17:49:54 UTC (rev 17025)
@@ -52,7 +52,66 @@
!-------------------------------------------------------------------------------------------------
!
+ subroutine model_external_broadcast(myrank)
+! standard routine to setup model
+
+ use external_model
+
+ implicit none
+
+ include "constants.h"
+ ! standard include of the MPI library
+ include 'mpif.h'
+
+ integer :: myrank
+
+ ! local parameters
+ integer :: idummy
+
+ ! dummy to ignore compiler warnings
+ idummy = myrank
+
+!---
+!
+! ADD YOUR MODEL HERE
+!
+!---
+
+ ! the variables read are declared and stored in structure MEXT_V
+ !if(myrank == 0) call read_external_model()
+
+ ! broadcast the information read on the master to the nodes
+ !call MPI_BCAST(MEXT_V%dvs,size(MEXT_V%dvs),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+
+ end subroutine model_external_broadcast
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+!
+! subroutine read_external_model()
+!
+! use external_model
+!
+! implicit none
+!
+! include "constants.h"
+!---
+!
+! ADD YOUR MODEL HERE
+!
+!---
+!
+! end subroutine read_external_model
+
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+
subroutine model_external_values(i,j,k,ispec,idomain_id,imaterial_id,&
nspec,ibool, &
iflag_aniso,iflag_atten, &
@@ -156,56 +215,4 @@
end subroutine model_external_values
-!
-!-------------------------------------------------------------------------------------------------
-!
-
- subroutine model_external_broadcast(myrank)
-
-! standard routine to setup model
-
- use external_model
-
- implicit none
-
- include "constants.h"
- ! standard include of the MPI library
- include 'mpif.h'
-
- integer :: myrank
-
- ! local parameters
- integer :: ier
-
- ! dummy to ignore compiler warnings
- ier = myrank
-
-!---
-!
-! ADD YOUR MODEL HERE
-!
-!---
-
- ! the variables read are declared and stored in structure D3MM_V
- !if(myrank == 0) call read_external_model()
-
- ! broadcast the information read on the master to the nodes
- !call MPI_BCAST(MEXT_V%dvs,size(MEXT_V%dvs),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
-
- end subroutine model_external_broadcast
-
-!
-!-------------------------------------------------------------------------------------------------
-!
-
-!
-! subroutine read_external_model()
-!
-! use external_model
-!
-! implicit none
-!
-! include "constants.h"
-!
-! end subroutine read_external_model
\ No newline at end of file
Added: seismo/3D/SPECFEM3D/trunk/model_interface_bedrock.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/model_interface_bedrock.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/model_interface_bedrock.f90 2010-06-28 17:49:54 UTC (rev 17025)
@@ -0,0 +1,390 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology September 2006
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+! interface model file
+! example file only, unused so far
+
+! ! Piero
+! module bedrock
+!
+! real,dimension(:,:),allocatable :: ibedrock
+!
+! end module bedrock
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+! subroutine model_bedrock_broadcast(myrank)
+!
+!! standard routine to setup model
+!
+! use bedrock
+!
+! implicit none
+!
+! include "constants.h"
+! ! standard include of the MPI library
+! include 'mpif.h'
+!
+! integer :: myrank
+!
+! ! local parameters
+! integer :: idummy
+!
+! ! dummy to ignore compiler warnings
+! idummy = myrank
+!
+! allocate(ibedrock(NX_TOPO_ANT,NY_TOPO_ANT))
+
+! if(myrank == 0) then
+! call read_bedrock_file(ibedrock)
+! ! write(IMAIN,*)
+! ! write(IMAIN,*) 'regional bedrock file read ranges in m from ',minval(ibedrock),' to ',maxval(ibedrock)
+! ! write(IMAIN,*)
+! endif
+
+! ! broadcast the information read on the master to the nodes
+! ! call MPI_BCAST(ibedrock,NX_TOPO_ANT*NY_TOPO_ANT,MPI_REAL,0,MPI_COMM_WORLD,ier)
+! call bcast_all_cr(ibedrock,NX_TOPO_ANT*NY_TOPO_ANT)
+
+! end subroutine model_bedrock_broadcast
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+!
+! subroutine read_bedrock_file()
+!
+! use bedrock
+!
+! implicit none
+!
+! include "constants.h"
+!---
+!
+! ADD YOUR MODEL HERE
+!
+!---
+!
+! end subroutine read_external_model
+
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+
+! subroutine model_bedrock_store()
+!
+! use bedrock
+!
+! implicit none
+!
+! !! DK DK store the position of the six stations to be able to
+! !! DK DK exclude circles around each station to make sure they are on the bedrock
+! !! DK DK and not in the ice
+! utm_x_station(1) = 783500.6250000d0
+! utm_y_station(1) = -11828.7519531d0
+
+! utm_x_station(2) = 853644.5000000d0
+! utm_y_station(2) = -114.0138092d0
+
+! utm_x_station(3) = 863406.0000000d0
+! utm_y_station(3) = -53736.1640625d0
+
+! utm_x_station(4) = 823398.8125000d0
+! utm_y_station(4) = 29847.4511719d0
+
+! utm_x_station(5) = 863545.3750000d0
+! utm_y_station(5) = 19669.6621094d0
+
+! utm_x_station(6) = 817099.3750000d0
+! utm_y_station(6) = -24430.2871094d0
+
+! print*,myrank,'après store the position of the six stations'
+! call flush(6)
+
+! print*, myrank,minval(nodes_coords_ext_mesh(1,:))
+! call flush(6)
+
+
+! print*, myrank,maxval(nodes_coords_ext_mesh(1,:))
+! call flush(6)
+
+
+! do ispec = 1, nspec
+
+! zmesh = zstore(2,2,2,ispec)
+
+! ! if(doubling_index == IFLAG_ONE_LAYER_TOPOGRAPHY) then
+! if(any(ibelm_top == ispec)) then
+! doubling_value_found_for_Piero = IFLAG_ONE_LAYER_TOPOGRAPHY
+
+! else if(zmesh < Z_23p4km) then
+! doubling_value_found_for_Piero = IFLAG_MANTLE_BELOW_23p4km
+
+! else if(zmesh < Z_14km) then
+! doubling_value_found_for_Piero = IFLAG_14km_to_23p4km
+
+! else
+! doubling_value_found_for_Piero = IFLAG_BEDROCK_down_to_14km
+! endif
+! idoubling(ispec) = doubling_value_found_for_Piero
+
+! do k = 1, NGLLZ
+! do j = 1, NGLLY
+! do i = 1, NGLLX
+
+
+! if(idoubling(ispec) == IFLAG_ONE_LAYER_TOPOGRAPHY .or. &
+! idoubling(ispec) == IFLAG_BEDROCK_down_to_14km) then
+
+! ! since we have suppressed UTM projection for Piero Basini, UTMx is the same as long
+! ! and UTMy is the same as lat
+! long = xstore(i,j,k,ispec)
+! lat = ystore(i,j,k,ispec)
+
+! ! get coordinate of corner in model
+! icornerlong = int((long - ORIG_LONG_TOPO) / DEGREES_PER_CELL_TOPO) + 1
+! icornerlat = int((lat - ORIG_LAT_TOPO) / DEGREES_PER_CELL_TOPO) + 1
+
+! ! avoid edge effects and extend with identical point if outside model
+! if(icornerlong < 1) icornerlong = 1
+! if(icornerlong > NX_TOPO-1) icornerlong = NX_TOPO-1
+! if(icornerlat < 1) icornerlat = 1
+! if(icornerlat > NY_TOPO-1) icornerlat = NY_TOPO-1
+
+! ! compute coordinates of corner
+! long_corner = ORIG_LONG_TOPO + (icornerlong-1)*DEGREES_PER_CELL_TOPO
+! lat_corner = ORIG_LAT_TOPO + (icornerlat-1)*DEGREES_PER_CELL_TOPO
+
+! ! compute ratio for interpolation
+! ratio_xi = (long - long_corner) / DEGREES_PER_CELL_TOPO
+! ratio_eta = (lat - lat_corner) / DEGREES_PER_CELL_TOPO
+
+! ! avoid edge effects
+! if(ratio_xi < 0.) ratio_xi = 0.
+! if(ratio_xi > 1.) ratio_xi = 1.
+! if(ratio_eta < 0.) ratio_eta = 0.
+! if(ratio_eta > 1.) ratio_eta = 1.
+
+! ! interpolate elevation at current point
+! elevation_bedrock = &
+! ibedrock(icornerlong,icornerlat)*(1.-ratio_xi)*(1.-ratio_eta) + &
+! ibedrock(icornerlong+1,icornerlat)*ratio_xi*(1.-ratio_eta) + &
+! ibedrock(icornerlong+1,icornerlat+1)*ratio_xi*ratio_eta + &
+! ibedrock(icornerlong,icornerlat+1)*(1.-ratio_xi)*ratio_eta
+
+! !! DK DK exclude circles around each station to make sure they are on the bedrock
+! !! DK DK and not in the ice
+! is_around_a_station = .false.
+! do istation = 1,NUMBER_OF_STATIONS
+! if(sqrt((long - utm_x_station(istation))**2 + (lat - utm_y_station(istation))**2) < RADIUS_TO_EXCLUDE) then
+! is_around_a_station = .true.
+! exit
+! endif
+! enddo
+
+! ! define elastic parameters in the model
+
+! ! we are above the bedrock interface i.e. in the ice, and not too close to a station
+! if(zmesh >= elevation_bedrock .and. .not. is_around_a_station) then
+! vp = 3800.d0
+! vs = 1900.d0
+! rho = 900.d0
+! iflag_attenuation_store(i,j,k,ispec) = IATTENUATION_ICE
+
+! ! we are below the bedrock interface i.e. in the bedrock, or close to a station
+! else
+! vp = 5800.d0
+! vs = 3200.d0
+! rho = 2600.d0
+! iflag_attenuation_store(i,j,k,ispec) = IATTENUATION_BEDROCK
+! endif
+
+! else if(idoubling(ispec) == IFLAG_14km_to_23p4km) then
+! vp = 6800.d0
+! vs = 3900.d0
+! rho = 2900.d0
+! iflag_attenuation_store(i,j,k,ispec) = IATTENUATION_BEDROCK
+
+! else if(idoubling(ispec) == IFLAG_MANTLE_BELOW_23p4km) then
+! vp = 8100.d0
+! vs = 4480.d0
+! rho = 3380.d0
+! iflag_attenuation_store(i,j,k,ispec) = IATTENUATION_BEDROCK
+
+! endif
+
+! !pll 8/06
+! if(CUSTOM_REAL == SIZE_REAL) then
+! rhostore(i,j,k,ispec) = sngl(rho)
+! vpstore(i,j,k,ispec) = sngl(vp)
+! vsstore(i,j,k,ispec) = sngl(vs)
+! else
+! rhostore(i,j,k,ispec) = rho
+! vpstore(i,j,k,ispec) = vp
+! vsstore(i,j,k,ispec) = vs
+! end if
+
+! kappastore(i,j,k,ispec) = rhostore(i,j,k,ispec)*(vpstore(i,j,k,ispec)*vpstore(i,j,k,ispec) - &
+! 4.d0*vsstore(i,j,k,ispec)*vsstore(i,j,k,ispec)/3.d0)
+! mustore(i,j,k,ispec) = rhostore(i,j,k,ispec)*vsstore(i,j,k,ispec)*&
+! vsstore(i,j,k,ispec)
+
+! ! Stacey, a completer par la suite
+! rho_vp(i,j,k,ispec) = rhostore(i,j,k,ispec)*vpstore(i,j,k,ispec)
+! rho_vs(i,j,k,ispec) = rhostore(i,j,k,ispec)*vsstore(i,j,k,ispec)
+! !end pll
+
+! ! kappastore(i,j,k,ispec) = materials_ext_mesh(1,mat_ext_mesh(ispec))* &
+! ! (materials_ext_mesh(2,mat_ext_mesh(ispec))*materials_ext_mesh(2,mat_ext_mesh(ispec)) - &
+! ! 4.d0*materials_ext_mesh(3,mat_ext_mesh(ispec))*materials_ext_mesh(3,mat_ext_mesh(ispec))/3.d0)
+! ! mustore(i,j,k,ispec) = materials_ext_mesh(1,mat_ext_mesh(ispec))* &
+! materials_ext_mesh(3,mat_ext_mesh(ispec))*&
+! ! x materials_ext_mesh(3,mat_ext_mesh(ispec))
+! enddo
+! enddo
+! enddo
+! enddo
+!
+! end subroutine
+
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+
+!pll
+! subroutine interface(iflag,flag_below,flag_above,ispec,nspec,i,j,k,xstore,ystore,zstore,ibedrock)
+
+! implicit none
+
+! include "constants.h"
+
+! integer :: iflag,flag_below,flag_above
+! integer :: ispec,nspec
+! integer :: i,j,k
+! double precision, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: xstore,ystore,zstore
+! real(kind=CUSTOM_REAL), dimension(NX_TOPO_ANT,NY_TOPO_ANT) :: ibedrock
+! integer, parameter :: NUMBER_OF_STATIONS = 1
+! double precision, parameter :: RADIUS_TO_EXCLUDE = 250.d0
+! double precision, dimension(NUMBER_OF_STATIONS) :: utm_x_station,utm_y_station
+
+! !-------------------
+
+! !for Piero
+! logical :: is_around_a_station
+! integer :: istation
+
+! ! store bedrock values
+! integer :: icornerlat,icornerlong
+! double precision :: lat,long,elevation_bedrock
+! double precision :: lat_corner,long_corner,ratio_xi,ratio_eta
+
+
+! !! DK DK store the position of the six stations to be able to
+! !! DK DK exclude circles around each station to make sure they are on the bedrock
+! !! DK DK and not in the ice
+! utm_x_station(1) = 783500.6250000d0
+! utm_y_station(1) = -11828.7519531d0
+
+! utm_x_station(2) = 853644.5000000d0
+! utm_y_station(2) = -114.0138092d0
+
+! utm_x_station(3) = 863406.0000000d0
+! utm_y_station(3) = -53736.1640625d0
+
+! utm_x_station(4) = 823398.8125000d0
+! utm_y_station(4) = 29847.4511719d0
+
+! utm_x_station(5) = 863545.3750000d0
+! utm_y_station(5) = 19669.6621094d0
+
+! utm_x_station(6) = 817099.3750000d0
+! utm_y_station(6) = -24430.2871094d0
+
+! ! since we have suppressed UTM projection for Piero Basini, UTMx is the same as long
+! ! and UTMy is the same as lat
+! long = xstore(i,j,k,ispec)
+! lat = ystore(i,j,k,ispec)
+
+! ! get coordinate of corner in model
+! icornerlong = int((long - ORIG_LONG_TOPO_ANT) / DEGREES_PER_CELL_TOPO_ANT) + 1
+! icornerlat = int((lat - ORIG_LAT_TOPO_ANT) / DEGREES_PER_CELL_TOPO_ANT) + 1
+
+! ! avoid edge effects and extend with identical point if outside model
+! if(icornerlong < 1) icornerlong = 1
+! if(icornerlong > NX_TOPO_ANT-1) icornerlong = NX_TOPO_ANT-1
+! if(icornerlat < 1) icornerlat = 1
+! if(icornerlat > NY_TOPO_ANT-1) icornerlat = NY_TOPO_ANT-1
+
+! ! compute coordinates of corner
+! long_corner = ORIG_LONG_TOPO_ANT + (icornerlong-1)*DEGREES_PER_CELL_TOPO_ANT
+! lat_corner = ORIG_LAT_TOPO_ANT + (icornerlat-1)*DEGREES_PER_CELL_TOPO_ANT
+
+! ! compute ratio for interpolation
+! ratio_xi = (long - long_corner) / DEGREES_PER_CELL_TOPO_ANT
+! ratio_eta = (lat - lat_corner) / DEGREES_PER_CELL_TOPO_ANT
+
+! ! avoid edge effects
+! if(ratio_xi < 0.) ratio_xi = 0.
+! if(ratio_xi > 1.) ratio_xi = 1.
+! if(ratio_eta < 0.) ratio_eta = 0.
+! if(ratio_eta > 1.) ratio_eta = 1.
+
+! ! interpolate elevation at current point
+! elevation_bedrock = &
+! ibedrock(icornerlong,icornerlat)*(1.-ratio_xi)*(1.-ratio_eta) + &
+! ibedrock(icornerlong+1,icornerlat)*ratio_xi*(1.-ratio_eta) + &
+! ibedrock(icornerlong+1,icornerlat+1)*ratio_xi*ratio_eta + &
+! ibedrock(icornerlong,icornerlat+1)*(1.-ratio_xi)*ratio_eta
+
+! !! DK DK exclude circles around each station to make sure they are on the bedrock
+! !! DK DK and not in the ice
+! is_around_a_station = .false.
+! do istation = 1,NUMBER_OF_STATIONS
+! if(sqrt((xstore(i,j,k,ispec) - utm_x_station(istation))**2 + (ystore(i,j,k,ispec) - &
+! utm_y_station(istation))**2) < RADIUS_TO_EXCLUDE) then
+! is_around_a_station = .true.
+! exit
+! endif
+! enddo
+
+! ! we are above the bedrock interface i.e. in the ice, and not too close to a station
+! if(zstore(i,j,k,ispec) >= elevation_bedrock .and. .not. is_around_a_station) then
+! iflag = flag_above
+! !iflag_attenuation_store(i,j,k,ispec) = IATTENUATION_ICE
+! ! we are below the bedrock interface i.e. in the bedrock, or close to a station
+! else
+! iflag = flag_below
+! !iflag_attenuation_store(i,j,k,ispec) = IATTENUATION_BEDROCK
+! endif
+
+
+! end subroutine interface
Added: seismo/3D/SPECFEM3D/trunk/model_tomography.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/model_tomography.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/model_tomography.f90 2010-06-28 17:49:54 UTC (rev 17025)
@@ -0,0 +1,355 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology September 2006
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+! generic tomography file
+!
+! note: the idea is to use an external, tomography velocity model
+!
+! most of the routines here are place-holders, please add/implement your own routines
+!
+
+ module tomography
+
+ include "constants.h"
+
+ ! for external tomography....
+ ! (regular spaced, xyz-block file in ascii)
+ character (len=80) :: TOMO_FILENAME = 'DATA/veryfast_tomography_abruzzo_complete.xyz'
+
+ ! model dimensions
+ double precision :: ORIG_X,ORIG_Y,ORIG_Z
+ double precision :: END_X,END_Y,END_Z
+ double precision :: SPACING_X,SPACING_Y,SPACING_Z
+
+ ! model parameter records
+ real(kind=CUSTOM_REAL), dimension (:), allocatable :: vp_tomography,vs_tomography,rho_tomography,z_tomography
+
+ ! model entries
+ integer :: NX,NY,NZ
+ integer :: nrecord
+
+ ! min/max statistics
+ double precision :: VP_MIN,VS_MIN,RHO_MIN,VP_MAX,VS_MAX,RHO_MAX
+
+ end module tomography
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ subroutine model_tomography_broadcast(myrank)
+
+ implicit none
+
+ ! include "constants.h"
+ ! include "precision.h"
+ ! include 'mpif.h'
+ integer :: myrank
+
+ ! all processes read in same file
+ ! note: for a high number of processes this might lead to a bottleneck
+ call read_model_tomography(myrank)
+
+ ! otherwise:
+
+ ! only master reads in model file
+ !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) )
+ !call MPI_BCAST(vp_tomography,size(vp_tomography),CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
+
+ end subroutine model_tomography_broadcast
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ subroutine read_model_tomography(myrank)
+
+! start magnoni 29/11/09
+! read Vp Vs and rho from extracted text file
+
+! assuming that only tomography undefined material is allowed....
+! and all the tomographic regions are collect inside one file called TOMO_FILENAME with homogenous resolution
+! this could be problematic for example if the tomographic regions have different resolution
+! leading to a waste of memory and cpu time in the partitioning process
+
+ use tomography
+
+ implicit none
+
+ integer :: myrank
+
+ ! local parameters
+ real(kind=CUSTOM_REAL) :: x_tomo,y_tomo,z_tomo,vp_tomo,vs_tomo,rho_tomo
+ integer :: irecord,ier
+
+ !TOMO_FILENAME='DATA/veryfast_tomography_abruzzo_complete.xyz'
+ ! probably the simple position for the filename is the constat.h
+ ! but it is also possible to include the name of the file in the material file (therefore in the undef_mat_prop)
+ ! if we want more than one tomofile (Examples: 2 file with a differente resolution
+ ! as in los angeles case we need to loop over mat_ext_mesh(1,ispec)...
+ ! it is a possible solution )
+ ! magnoni 1/12/09
+ open(unit=27,file=TOMO_FILENAME,status='old',iostat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error reading tomography file')
+
+ ! reads in model dimensions
+ read(27,*) ORIG_X, ORIG_Y, ORIG_Z, END_X, END_Y, END_Z
+ read(27,*) SPACING_X, SPACING_Y, SPACING_Z
+ read(27,*) NX, NY, NZ
+ read(27,*) VP_MIN, VP_MAX, VS_MIN, VS_MAX, RHO_MIN, RHO_MAX
+
+ nrecord = NX*NY*NZ
+
+ ! allocates model records
+ allocate(vp_tomography(1:nrecord), &
+ vs_tomography(1:nrecord), &
+ rho_tomography(1:nrecord), &
+ z_tomography(1:nrecord),stat=ier)
+ if(ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays')
+
+ ! reads in record sections
+ do irecord = 1,nrecord
+ read(27,*) x_tomo,y_tomo,z_tomo,vp_tomo,vs_tomo,rho_tomo
+
+ ! stores record values
+ vp_tomography(irecord) = vp_tomo
+ vs_tomography(irecord) = vs_tomo
+ rho_tomography(irecord) = rho_tomo
+ z_tomography(irecord) = z_tomo
+ enddo
+
+ close(27)
+
+ end subroutine read_model_tomography
+
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+
+ subroutine model_tomography(x_eval,y_eval,z_eval, &
+ rho_final,vp_final,vs_final)
+
+ use tomography
+
+ implicit none
+
+ !integer, intent(in) :: NX,NY,NZ
+ !real(kind=CUSTOM_REAL), dimension(1:NX*NY*NZ), intent(in) :: vp_tomography,vs_tomography,rho_tomography,z_tomography
+ !double precision, intent(in) :: ORIG_X,ORIG_Y,ORIG_Z,SPACING_X,SPACING_Y,SPACING_Z
+ !double precision, intent(in) :: VP_MIN,VS_MIN,RHO_MIN,VP_MAX,VS_MAX,RHO_MAX
+
+ double precision, intent(in) :: x_eval,y_eval,z_eval
+ real(kind=CUSTOM_REAL), intent(out) :: vp_final,vs_final,rho_final
+
+ ! local parameters
+ integer :: ix,iy,iz
+ integer :: p0,p1,p2,p3,p4,p5,p6,p7
+
+ double precision :: spac_x,spac_y,spac_z
+ double precision :: gamma_interp_x,gamma_interp_y
+ double precision :: gamma_interp_z1,gamma_interp_z2,gamma_interp_z3, &
+ gamma_interp_z4,gamma_interp_z5,gamma_interp_z6,gamma_interp_z7,gamma_interp_z8
+ real(kind=CUSTOM_REAL) :: vp1,vp2,vp3,vp4,vp5,vp6,vp7,vp8, &
+ vs1,vs2,vs3,vs4,vs5,vs6,vs7,vs8,rho1,rho2,rho3,rho4,rho5,rho6,rho7,rho8
+
+ ! determine spacing and cell for linear interpolation
+ spac_x = (x_eval - ORIG_X) / SPACING_X
+ spac_y = (y_eval - ORIG_Y) / SPACING_Y
+ spac_z = (z_eval - ORIG_Z) / SPACING_Z
+
+ ix = int(spac_x)
+ iy = int(spac_y)
+ iz = int(spac_z)
+
+ gamma_interp_x = spac_x - dble(ix)
+ gamma_interp_y = spac_y - dble(iy)
+
+ ! suppress edge effects for points outside of the model SPOSTARE DOPO
+ if(ix < 0) then
+ ix = 0
+ gamma_interp_x = 0.d0
+ endif
+ if(ix > NX-2) then
+ ix = NX-2
+ gamma_interp_x = 1.d0
+ endif
+
+ if(iy < 0) then
+ iy = 0
+ gamma_interp_y = 0.d0
+ endif
+ if(iy > NY-2) then
+ iy = NY-2
+ gamma_interp_y = 1.d0
+ endif
+
+ if(iz < 0) then
+ iz = 0
+ ! gamma_interp_z = 0.d0
+ endif
+ if(iz > NZ-2) then
+ iz = NZ-2
+ ! gamma_interp_z = 1.d0
+ endif
+
+
+ ! define 8 corners of interpolation element
+ p0 = ix+iy*NX+iz*(NX*NY)
+ p1 = (ix+1)+iy*NX+iz*(NX*NY)
+ p2 = (ix+1)+(iy+1)*NX+iz*(NX*NY)
+ p3 = ix+(iy+1)*NX+iz*(NX*NY)
+ p4 = ix+iy*NX+(iz+1)*(NX*NY)
+ p5 = (ix+1)+iy*NX+(iz+1)*(NX*NY)
+ p6 = (ix+1)+(iy+1)*NX+(iz+1)*(NX*NY)
+ p7 = ix+(iy+1)*NX+(iz+1)*(NX*NY)
+
+ if(z_tomography(p4+1) == z_tomography(p0+1)) then
+ gamma_interp_z1 = 1.d0
+ else
+ gamma_interp_z1 = (z_eval-z_tomography(p0+1))/(z_tomography(p4+1)-z_tomography(p0+1))
+ endif
+ if(gamma_interp_z1 > 1.d0) then
+ gamma_interp_z1 = 1.d0
+ endif
+ if(gamma_interp_z1 < 0.d0) then
+ gamma_interp_z1 = 0.d0
+ endif
+
+
+ if(z_tomography(p5+1) == z_tomography(p1+1)) then
+ gamma_interp_z2 = 1.d0
+ else
+ gamma_interp_z2 = (z_eval-z_tomography(p1+1))/(z_tomography(p5+1)-z_tomography(p1+1))
+ endif
+ if(gamma_interp_z2 > 1.d0) then
+ gamma_interp_z2 = 1.d0
+ endif
+ if(gamma_interp_z2 < 0.d0) then
+ gamma_interp_z2 = 0.d0
+ endif
+
+
+ if(z_tomography(p6+1) == z_tomography(p2+1)) then
+ gamma_interp_z3 = 1.d0
+ else
+ gamma_interp_z3 = (z_eval-z_tomography(p2+1))/(z_tomography(p6+1)-z_tomography(p2+1))
+ endif
+ if(gamma_interp_z3 > 1.d0) then
+ gamma_interp_z3 = 1.d0
+ endif
+ if(gamma_interp_z3 < 0.d0) then
+ gamma_interp_z3 = 0.d0
+ endif
+
+
+ if(z_tomography(p7+1) == z_tomography(p3+1)) then
+ gamma_interp_z4 = 1.d0
+ else
+ gamma_interp_z4 = (z_eval-z_tomography(p3+1))/(z_tomography(p7+1)-z_tomography(p3+1))
+ endif
+ if(gamma_interp_z4 > 1.d0) then
+ gamma_interp_z4 = 1.d0
+ endif
+ if(gamma_interp_z4 < 0.d0) then
+ gamma_interp_z4 = 0.d0
+ endif
+
+ gamma_interp_z5 = 1. - gamma_interp_z1
+ gamma_interp_z6 = 1. - gamma_interp_z2
+ gamma_interp_z7 = 1. - gamma_interp_z3
+ gamma_interp_z8 = 1. - gamma_interp_z4
+
+ vp1 = vp_tomography(p0+1)
+ vp2 = vp_tomography(p1+1)
+ vp3 = vp_tomography(p2+1)
+ vp4 = vp_tomography(p3+1)
+ vp5 = vp_tomography(p4+1)
+ vp6 = vp_tomography(p5+1)
+ vp7 = vp_tomography(p6+1)
+ vp8 = vp_tomography(p7+1)
+
+ vs1 = vs_tomography(p0+1)
+ vs2 = vs_tomography(p1+1)
+ vs3 = vs_tomography(p2+1)
+ vs4 = vs_tomography(p3+1)
+ vs5 = vs_tomography(p4+1)
+ vs6 = vs_tomography(p5+1)
+ vs7 = vs_tomography(p6+1)
+ vs8 = vs_tomography(p7+1)
+
+ rho1 = rho_tomography(p0+1)
+ rho2 = rho_tomography(p1+1)
+ rho3 = rho_tomography(p2+1)
+ rho4 = rho_tomography(p3+1)
+ rho5 = rho_tomography(p4+1)
+ rho6 = rho_tomography(p5+1)
+ rho7 = rho_tomography(p6+1)
+ rho8 = rho_tomography(p7+1)
+
+ ! use trilinear interpolation in cell to define Vp Vs and rho
+ vp_final = &
+ vp1*(1.-gamma_interp_x)*(1.-gamma_interp_y)*(1.-gamma_interp_z1) + &
+ vp2*gamma_interp_x*(1.-gamma_interp_y)*(1.-gamma_interp_z2) + &
+ vp3*gamma_interp_x*gamma_interp_y*(1.-gamma_interp_z3) + &
+ vp4*(1.-gamma_interp_x)*gamma_interp_y*(1.-gamma_interp_z4) + &
+ vp5*(1.-gamma_interp_x)*(1.-gamma_interp_y)*gamma_interp_z1 + &
+ vp6*gamma_interp_x*(1.-gamma_interp_y)*gamma_interp_z2 + &
+ vp7*gamma_interp_x*gamma_interp_y*gamma_interp_z3 + &
+ vp8*(1.-gamma_interp_x)*gamma_interp_y*gamma_interp_z4
+
+ vs_final = &
+ vs1*(1.-gamma_interp_x)*(1.-gamma_interp_y)*(1.-gamma_interp_z1) + &
+ vs2*gamma_interp_x*(1.-gamma_interp_y)*(1.-gamma_interp_z2) + &
+ vs3*gamma_interp_x*gamma_interp_y*(1.-gamma_interp_z3) + &
+ vs4*(1.-gamma_interp_x)*gamma_interp_y*(1.-gamma_interp_z4) + &
+ vs5*(1.-gamma_interp_x)*(1.-gamma_interp_y)*gamma_interp_z1 + &
+ vs6*gamma_interp_x*(1.-gamma_interp_y)*gamma_interp_z2 + &
+ vs7*gamma_interp_x*gamma_interp_y*gamma_interp_z3 + &
+ vs8*(1.-gamma_interp_x)*gamma_interp_y*gamma_interp_z4
+
+ rho_final = &
+ rho1*(1.-gamma_interp_x)*(1.-gamma_interp_y)*(1.-gamma_interp_z1) + &
+ rho2*gamma_interp_x*(1.-gamma_interp_y)*(1.-gamma_interp_z2) + &
+ rho3*gamma_interp_x*gamma_interp_y*(1.-gamma_interp_z3) + &
+ rho4*(1.-gamma_interp_x)*gamma_interp_y*(1.-gamma_interp_z4) + &
+ rho5*(1.-gamma_interp_x)*(1.-gamma_interp_y)*gamma_interp_z1 + &
+ rho6*gamma_interp_x*(1.-gamma_interp_y)*gamma_interp_z2 + &
+ rho7*gamma_interp_x*gamma_interp_y*gamma_interp_z3 + &
+ rho8*(1.-gamma_interp_x)*gamma_interp_y*gamma_interp_z4
+
+ ! impose minimum and maximum velocity and density if needed
+ if(vp_final < VP_MIN) vp_final = VP_MIN
+ if(vs_final < VS_MIN) vs_final = VS_MIN
+ if(rho_final < RHO_MIN) rho_final = RHO_MIN
+ if(vp_final > VP_MAX) vp_final = VP_MAX
+ if(vs_final > VS_MAX) vs_final = VS_MAX
+ if(rho_final > RHO_MAX) rho_final = RHO_MAX
+
+ end subroutine model_tomography
More information about the CIG-COMMITS
mailing list