[cig-commits] r15453 - in seismo/3D/SPECFEM3D_SESAME/trunk: . DATA decompose_mesh_SCOTCH
pieyre at geodynamics.org
pieyre at geodynamics.org
Fri Jul 10 14:02:54 PDT 2009
Author: pieyre
Date: 2009-07-10 14:02:53 -0700 (Fri, 10 Jul 2009)
New Revision: 15453
Modified:
seismo/3D/SPECFEM3D_SESAME/trunk/DATA/Par_file
seismo/3D/SPECFEM3D_SESAME/trunk/Makefile.in
seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces_with_Deville.f90
seismo/3D/SPECFEM3D_SESAME/trunk/constants.h.in
seismo/3D/SPECFEM3D_SESAME/trunk/create_movie_shakemap_AVS_DX_GMT.f90
seismo/3D/SPECFEM3D_SESAME/trunk/create_regions_mesh.f90
seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/compile_all.csh
seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/constants_decompose_mesh_SCOTCH.h
seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90
seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90
seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/scotchf.h
seismo/3D/SPECFEM3D_SESAME/trunk/generate_databases.f90
seismo/3D/SPECFEM3D_SESAME/trunk/get_absorb.f90
seismo/3D/SPECFEM3D_SESAME/trunk/locate_receivers.f90
seismo/3D/SPECFEM3D_SESAME/trunk/locate_source.f90
seismo/3D/SPECFEM3D_SESAME/trunk/save_header_file.f90
seismo/3D/SPECFEM3D_SESAME/trunk/specfem3D.f90
seismo/3D/SPECFEM3D_SESAME/trunk/todo_list_please_dont_remove.txt
Log:
Main modifications : new format for external meshes, attenuation, moment tensor source ... all validated
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/DATA/Par_file
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/DATA/Par_file 2009-07-10 16:39:24 UTC (rev 15452)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/DATA/Par_file 2009-07-10 21:02:53 UTC (rev 15453)
@@ -3,38 +3,38 @@
SIMULATION_TYPE = 1 # 1 = forward, 2 = adjoint, 3 = both simultaneously
SAVE_FORWARD = .false.
-# coordinates of mesh block in latitude/longitude and depth in km
+# coordinates of mesh block in latitude/longitude and depth in km USELESS
LATITUDE_MIN = 33.8d0
LATITUDE_MAX = 34.1d0
LONGITUDE_MIN = -118.4d0
LONGITUDE_MAX = -118.1d0
DEPTH_BLOCK_KM = 60.d0
UTM_PROJECTION_ZONE = 11
-SUPPRESS_UTM_PROJECTION = .false.
+SUPPRESS_UTM_PROJECTION = .true.
-# number of elements at the surface along edges of the mesh at the surface
+# number of elements at the surface along edges of the mesh at the surface USELESS
# (must be 8 * multiple of NPROC below if mesh is not regular and contains mesh doublings)
# (must be multiple of NPROC below if mesh is regular)
NEX_XI = 16
NEX_ETA = 16
-# number of MPI processors along xi and eta (can be different)
+# number of MPI processors along xi and eta (can be different) USELESS
NPROC_XI = 2
NPROC_ETA = 2
-# model (SoCal, Harvard_LA, Min_Chen_anisotropy)
+# model (SoCal, Harvard_LA, Min_Chen_anisotropy)
MODEL = Harvard_LA
# parameters describing the model
-OCEANS = .true.
-TOPOGRAPHY = .true.
-ATTENUATION = .true.
+OCEANS = .false.
+TOPOGRAPHY = .false.
+ATTENUATION = .false.
USE_OLSEN_ATTENUATION = .false.
# absorbing boundary conditions for a regional simulation
-ABSORBING_CONDITIONS = .true.
+ABSORBING_CONDITIONS = .false.
-# record length in seconds
+# record length in seconds USELESS see constants.h
RECORD_LENGTH_IN_SECONDS = 30.0
# save AVS or OpenDX movies
@@ -50,7 +50,7 @@
SAVE_MESH_FILES = .false.
# path to store the local database file on each node
-LOCAL_PATH = /scratch/komatits/DATABASES_MPI
+LOCAL_PATH = ./OUTPUT_FILES
# machine file for MPI
# this is not needed for new cluster and is in go_mesher/go_solver
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/Makefile.in 2009-07-10 16:39:24 UTC (rev 15452)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/Makefile.in 2009-07-10 21:02:53 UTC (rev 15453)
@@ -109,8 +109,6 @@
$O/write_AVS_DX_surface_data.o \
$O/write_seismograms.o \
$O/compute_boundary_kernel.o \
- $O/compute_forces_no_Deville.o \
- $O/compute_forces_with_Deville.o \
$O/memory_eval.o \
$(EMPTY_MACRO)
@@ -120,6 +118,8 @@
$O/assemble_MPI_scalar.o \
$O/assemble_MPI_vector.o \
$O/read_arrays_solver.o \
+ $O/compute_forces_no_Deville.o \
+ $O/compute_forces_with_Deville.o \
$O/specfem3D.o \
$(EMPTY_MACRO)
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces_with_Deville.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces_with_Deville.f90 2009-07-10 16:39:24 UTC (rev 15452)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces_with_Deville.f90 2009-07-10 21:02:53 UTC (rev 15453)
@@ -23,14 +23,26 @@
!
!=====================================================================
-subroutine compute_forces_with_Deville(NSPEC_AB,NGLOB_AB,displ,accel,xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+subroutine compute_forces_with_Deville(NSPEC_AB,NGLOB_AB,ATTENUATION_VAL,displ,accel, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT,wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
kappastore,mustore,jacobian,ibool,ispec_is_inner,phase_is_inner, &
- NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,xi_source,eta_source,gamma_source,nu_source,hdur,dt)
+ NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,xi_source,eta_source,gamma_source,nu_source, &
+ hdur,hdur_gaussian,t_cmt,dt,stf,t0,sourcearrays, & !pll
+ one_minus_sum_beta,factor_common,alphaval,betaval,gammaval,R_xx,R_yy,R_xy,R_xz,R_yz, &
+ epsilondev_xx,epsilondev_yy,epsilondev_xy,epsilondev_xz,epsilondev_yz,iflag_attenuation_store, &
+ ABSORBING_CONDITIONS,SAVE_FORWARD,NSTEP,SIMULATION_TYPE, &
+ nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax,NSPEC2D_BOTTOM,NSPEC2DMAX_XMIN_XMAX_ext,NSPEC2DMAX_YMIN_YMAX_ext,&
+ ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom, &
+ nimin,nimax,njmin,njmax,nkmin_xi,nkmin_eta, &
+ veloc,rho_vp,rho_vs,jacobian2D_xmin,jacobian2D_xmax,jacobian2D_ymin,jacobian2D_ymax,jacobian2D_bottom, &
+ normal_xmin,normal_xmax,normal_ymin,normal_ymax,normal_bottom)
implicit none
include "constants.h"
+! include values created by the mesher
+! include "OUTPUT_FILES/values_from_mesher.h"
integer :: NSPEC_AB,NGLOB_AB
@@ -62,8 +74,9 @@
integer, dimension(NSOURCES) :: islice_selected_source,ispec_selected_source
double precision, dimension(NSOURCES) :: xi_source,eta_source,gamma_source
double precision, dimension(3,3,NSOURCES) :: nu_source
- double precision, dimension(NSOURCES) :: hdur
+ double precision, dimension(NSOURCES) :: hdur,hdur_gaussian,t_cmt
double precision :: dt
+ real(kind=CUSTOM_REAL), dimension(NSOURCES,NDIM,NGLLX,NGLLY,NGLLZ) :: sourcearrays
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3
@@ -116,6 +129,57 @@
integer :: isource
double precision :: t0,f0
+ double precision :: stf
+ real(kind=CUSTOM_REAL) stf_used
+ double precision, external :: comp_source_time_function
+
+! memory variables and standard linear solids for attenuation
+ integer i_SLS
+ integer iselected
+ real(kind=CUSTOM_REAL) R_xx_val,R_yy_val
+ real(kind=CUSTOM_REAL) factor_loc,alphaval_loc,betaval_loc,gammaval_loc,Sn,Snp1
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: epsilondev_xx_loc, &
+ epsilondev_yy_loc, epsilondev_xy_loc, epsilondev_xz_loc, epsilondev_yz_loc
+ real(kind=CUSTOM_REAL) epsilon_trace_over_3
+
+ logical :: ATTENUATION_VAL
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: iflag_attenuation_store
+ real(kind=CUSTOM_REAL), dimension(NUM_REGIONS_ATTENUATION) :: one_minus_sum_beta
+ real(kind=CUSTOM_REAL), dimension(NUM_REGIONS_ATTENUATION,N_SLS) :: factor_common, alphaval,betaval,gammaval
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB,N_SLS) :: &
+ R_xx,R_yy,R_xy,R_xz,R_yz
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
+ epsilondev_xx,epsilondev_yy,epsilondev_xy,epsilondev_xz,epsilondev_yz
+
+! Stacey conditions
+ logical :: ABSORBING_CONDITIONS,SAVE_FORWARD
+ integer :: NSTEP,SIMULATION_TYPE
+ integer :: nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax,NSPEC2D_BOTTOM
+ integer :: NSPEC2DMAX_XMIN_XMAX_ext,NSPEC2DMAX_YMIN_YMAX_ext
+ 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(2,NSPEC2DMAX_YMIN_YMAX_ext) :: nimin,nimax,nkmin_eta
+ integer, dimension(2,NSPEC2DMAX_XMIN_XMAX_ext) :: njmin,njmax,nkmin_xi
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: veloc
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: rho_vp,rho_vs
+ real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLZ,nspec2D_xmin) :: jacobian2D_xmin
+ real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLZ,nspec2D_xmax) :: jacobian2D_xmax
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec2D_ymin) :: jacobian2D_ymin
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec2D_ymax) :: jacobian2D_ymax
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NSPEC2D_BOTTOM) :: jacobian2D_bottom
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLLY,NGLLZ,nspec2D_xmin) :: normal_xmin
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLLY,NGLLZ,nspec2D_xmax) :: normal_xmax
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLLY,NGLLZ,nspec2D_ymin) :: normal_ymin
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLLY,NGLLZ,nspec2D_ymax) :: normal_ymax
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NSPEC2D_BOTTOM) :: normal_bottom
+
+ integer :: ispec2D
+ real(kind=CUSTOM_REAL) vx,vy,vz,nx,ny,nz,tx,ty,tz,vn,weight
+
+
do ispec = 1,NSPEC_AB
if (ispec_is_inner(ispec) .eqv. phase_is_inner) then
@@ -245,6 +309,19 @@
kappal = kappastore(i,j,k,ispec)
mul = mustore(i,j,k,ispec)
+
+ if(ATTENUATION_VAL) then
+ ! compute deviatoric strain
+ epsilon_trace_over_3 = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ epsilondev_xx_loc(i,j,k) = duxdxl - epsilon_trace_over_3
+ epsilondev_yy_loc(i,j,k) = duydyl - epsilon_trace_over_3
+ epsilondev_xy_loc(i,j,k) = 0.5 * duxdyl_plus_duydxl
+ epsilondev_xz_loc(i,j,k) = 0.5 * duzdxl_plus_duxdzl
+ epsilondev_yz_loc(i,j,k) = 0.5 * duzdyl_plus_duydzl
+
+ ! use unrelaxed parameters if attenuation
+ mul = mul * one_minus_sum_beta(iflag_attenuation_store(i,j,k,ispec))
+ endif
lambdalplus2mul = kappal + FOUR_THIRDS * mul
lambdal = lambdalplus2mul - 2.*mul
@@ -258,6 +335,20 @@
sigma_xz = mul*duzdxl_plus_duxdzl
sigma_yz = mul*duzdyl_plus_duydzl
+ ! subtract memory variables if attenuation
+ if(ATTENUATION_VAL) then
+ do i_sls = 1,N_SLS
+ R_xx_val = R_xx(i,j,k,ispec,i_sls)
+ R_yy_val = R_yy(i,j,k,ispec,i_sls)
+ sigma_xx = sigma_xx - R_xx_val
+ sigma_yy = sigma_yy - R_yy_val
+ sigma_zz = sigma_zz + R_xx_val + R_yy_val
+ sigma_xy = sigma_xy - R_xy(i,j,k,ispec,i_sls)
+ sigma_xz = sigma_xz - R_xz(i,j,k,ispec,i_sls)
+ sigma_yz = sigma_yz - R_yz(i,j,k,ispec,i_sls)
+ enddo
+ endif
+
! form dot product with test vector, symmetric form
tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_xy*xiyl + sigma_xz*xizl)
tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_yz*xizl)
@@ -365,53 +456,336 @@
accel(2,iglob) = accel(2,iglob) - fac1*newtempy1(i,j,k) - fac2*newtempy2(i,j,k) - fac3*newtempy3(i,j,k)
accel(3,iglob) = accel(3,iglob) - fac1*newtempz1(i,j,k) - fac2*newtempz2(i,j,k) - fac3*newtempz3(i,j,k)
+ ! update memory variables based upon the Runge-Kutta scheme
+ if(ATTENUATION_VAL) then
+
+ ! use Runge-Kutta scheme to march in time
+ do i_sls = 1,N_SLS
+
+ ! get coefficients for that standard linear solid
+ iselected = iflag_attenuation_store(i,j,k,ispec)
+ factor_loc = mustore(i,j,k,ispec) * factor_common(iselected,i_sls)
+ alphaval_loc = alphaval(iselected,i_sls)
+ betaval_loc = betaval(iselected,i_sls)
+ gammaval_loc = gammaval(iselected,i_sls)
+
+ ! term in xx
+ Sn = factor_loc * epsilondev_xx(i,j,k,ispec)
+ Snp1 = factor_loc * epsilondev_xx_loc(i,j,k)
+ R_xx(i,j,k,ispec,i_sls) = alphaval_loc * R_xx(i,j,k,ispec,i_sls) + betaval_loc * Sn + gammaval_loc * Snp1
+
+ ! term in yy
+ Sn = factor_loc * epsilondev_yy(i,j,k,ispec)
+ Snp1 = factor_loc * epsilondev_yy_loc(i,j,k)
+ R_yy(i,j,k,ispec,i_sls) = alphaval_loc * R_yy(i,j,k,ispec,i_sls) + betaval_loc * Sn + gammaval_loc * Snp1
+
+ ! term in zz not computed since zero trace
+
+ ! term in xy
+ Sn = factor_loc * epsilondev_xy(i,j,k,ispec)
+ Snp1 = factor_loc * epsilondev_xy_loc(i,j,k)
+ R_xy(i,j,k,ispec,i_sls) = alphaval_loc * R_xy(i,j,k,ispec,i_sls) + betaval_loc * Sn + gammaval_loc * Snp1
+
+ ! term in xz
+ Sn = factor_loc * epsilondev_xz(i,j,k,ispec)
+ Snp1 = factor_loc * epsilondev_xz_loc(i,j,k)
+ R_xz(i,j,k,ispec,i_sls) = alphaval_loc * R_xz(i,j,k,ispec,i_sls) + betaval_loc * Sn + gammaval_loc * Snp1
+
+ ! term in yz
+ Sn = factor_loc * epsilondev_yz(i,j,k,ispec)
+ Snp1 = factor_loc * epsilondev_yz_loc(i,j,k)
+ R_yz(i,j,k,ispec,i_sls) = alphaval_loc * R_yz(i,j,k,ispec,i_sls) + betaval_loc * Sn + gammaval_loc * Snp1
+
+ enddo ! end of loop on memory variables
+
+ endif ! end attenuation
+
enddo
enddo
enddo
+ ! save deviatoric strain for Runge-Kutta scheme
+ if(ATTENUATION_VAL) then
+ epsilondev_xx(:,:,:,ispec) = epsilondev_xx_loc(:,:,:)
+ epsilondev_yy(:,:,:,ispec) = epsilondev_yy_loc(:,:,:)
+ epsilondev_xy(:,:,:,ispec) = epsilondev_xy_loc(:,:,:)
+ epsilondev_xz(:,:,:,ispec) = epsilondev_xz_loc(:,:,:)
+ epsilondev_yz(:,:,:,ispec) = epsilondev_yz_loc(:,:,:)
+ endif
+
endif ! if (ispec_is_inner(ispec) .eqv. phase_is_inner)
enddo ! spectral element loop
+
+ ! add Stacey conditions
+ if(ABSORBING_CONDITIONS) then
+
+! xmin
+ do ispec2D=1,nspec2D_xmin
+
+ ispec=ibelm_xmin(ispec2D)
+
+ if (ispec_is_inner(ispec) .eqv. phase_is_inner) then
+
+ ! exclude elements that are not on absorbing edges
+ if(nkmin_xi(1,ispec2D) == 0 .or. njmin(1,ispec2D) == 0) cycle
+
+ i=1
+ do k=nkmin_xi(1,ispec2D),NGLLZ
+ do j=njmin(1,ispec2D),njmax(1,ispec2D)
+
+ iglob=ibool(i,j,k,ispec)
+
+ vx=veloc(1,iglob)
+ vy=veloc(2,iglob)
+ vz=veloc(3,iglob)
+ nx=normal_xmin(1,j,k,ispec2D)
+ ny=normal_xmin(2,j,k,ispec2D)
+ nz=normal_xmin(3,j,k,ispec2D)
+
+ vn=vx*nx+vy*ny+vz*nz
+
+ tx=rho_vp(i,j,k,ispec)*vn*nx+rho_vs(i,j,k,ispec)*(vx-vn*nx)
+ ty=rho_vp(i,j,k,ispec)*vn*ny+rho_vs(i,j,k,ispec)*(vy-vn*ny)
+ tz=rho_vp(i,j,k,ispec)*vn*nz+rho_vs(i,j,k,ispec)*(vz-vn*nz)
+
+ weight=jacobian2D_xmin(j,k,ispec2D)*wgllwgll_yz(j,k)
+
+ accel(1,iglob)=accel(1,iglob) - tx*weight
+ accel(2,iglob)=accel(2,iglob) - ty*weight
+ accel(3,iglob)=accel(3,iglob) - tz*weight
+
+ enddo
+ enddo
+ end if
+ enddo
+
+! xmax
+ do ispec2D=1,nspec2D_xmax
+
+ ispec=ibelm_xmax(ispec2D)
+
+ if (ispec_is_inner(ispec) .eqv. phase_is_inner) then
+
+ ! exclude elements that are not on absorbing edges
+ if(nkmin_xi(2,ispec2D) == 0 .or. njmin(2,ispec2D) == 0) cycle
+
+ i=NGLLX
+ do k=nkmin_xi(2,ispec2D),NGLLZ
+ do j=njmin(2,ispec2D),njmax(2,ispec2D)
+ iglob=ibool(i,j,k,ispec)
+
+ vx=veloc(1,iglob)
+ vy=veloc(2,iglob)
+ vz=veloc(3,iglob)
+
+ nx=normal_xmax(1,j,k,ispec2D)
+ ny=normal_xmax(2,j,k,ispec2D)
+ nz=normal_xmax(3,j,k,ispec2D)
+
+ vn=vx*nx+vy*ny+vz*nz
+
+ tx=rho_vp(i,j,k,ispec)*vn*nx+rho_vs(i,j,k,ispec)*(vx-vn*nx)
+ ty=rho_vp(i,j,k,ispec)*vn*ny+rho_vs(i,j,k,ispec)*(vy-vn*ny)
+ tz=rho_vp(i,j,k,ispec)*vn*nz+rho_vs(i,j,k,ispec)*(vz-vn*nz)
+
+ weight=jacobian2D_xmax(j,k,ispec2D)*wgllwgll_yz(j,k)
+
+ accel(1,iglob)=accel(1,iglob) - tx*weight
+ accel(2,iglob)=accel(2,iglob) - ty*weight
+ accel(3,iglob)=accel(3,iglob) - tz*weight
+
+ enddo
+ enddo
+ end if
+ enddo
+
+! ymin
+ do ispec2D=1,nspec2D_ymin
+
+ ispec=ibelm_ymin(ispec2D)
+
+ if (ispec_is_inner(ispec) .eqv. phase_is_inner) then
+
+ ! exclude elements that are not on absorbing edges
+ if(nkmin_eta(1,ispec2D) == 0 .or. nimin(1,ispec2D) == 0) cycle
+
+ j=1
+ do k=nkmin_eta(1,ispec2D),NGLLZ
+ do i=nimin(1,ispec2D),nimax(1,ispec2D)
+ iglob=ibool(i,j,k,ispec)
+
+ vx=veloc(1,iglob)
+ vy=veloc(2,iglob)
+ vz=veloc(3,iglob)
+
+ nx=normal_ymin(1,i,k,ispec2D)
+ ny=normal_ymin(2,i,k,ispec2D)
+ nz=normal_ymin(3,i,k,ispec2D)
+
+ vn=vx*nx+vy*ny+vz*nz
+
+ tx=rho_vp(i,j,k,ispec)*vn*nx+rho_vs(i,j,k,ispec)*(vx-vn*nx)
+ ty=rho_vp(i,j,k,ispec)*vn*ny+rho_vs(i,j,k,ispec)*(vy-vn*ny)
+ tz=rho_vp(i,j,k,ispec)*vn*nz+rho_vs(i,j,k,ispec)*(vz-vn*nz)
+
+ weight=jacobian2D_ymin(i,k,ispec2D)*wgllwgll_xz(i,k)
+
+ accel(1,iglob)=accel(1,iglob) - tx*weight
+ accel(2,iglob)=accel(2,iglob) - ty*weight
+ accel(3,iglob)=accel(3,iglob) - tz*weight
+
+ enddo
+ enddo
+ endif
+ enddo
+
+! ymax
+ do ispec2D=1,nspec2D_ymax
+
+ ispec=ibelm_ymax(ispec2D)
+
+ if (ispec_is_inner(ispec) .eqv. phase_is_inner) then
+
+ ! exclude elements that are not on absorbing edges
+ if(nkmin_eta(2,ispec2D) == 0 .or. nimin(2,ispec2D) == 0) cycle
+
+ j=NGLLY
+ do k=nkmin_eta(2,ispec2D),NGLLZ
+ do i=nimin(2,ispec2D),nimax(2,ispec2D)
+ iglob=ibool(i,j,k,ispec)
+
+ vx=veloc(1,iglob)
+ vy=veloc(2,iglob)
+ vz=veloc(3,iglob)
+
+ nx=normal_ymax(1,i,k,ispec2D)
+ ny=normal_ymax(2,i,k,ispec2D)
+ nz=normal_ymax(3,i,k,ispec2D)
+
+ vn=vx*nx+vy*ny+vz*nz
+
+ tx=rho_vp(i,j,k,ispec)*vn*nx+rho_vs(i,j,k,ispec)*(vx-vn*nx)
+ ty=rho_vp(i,j,k,ispec)*vn*ny+rho_vs(i,j,k,ispec)*(vy-vn*ny)
+ tz=rho_vp(i,j,k,ispec)*vn*nz+rho_vs(i,j,k,ispec)*(vz-vn*nz)
+
+ weight=jacobian2D_ymax(i,k,ispec2D)*wgllwgll_xz(i,k)
+
+ accel(1,iglob)=accel(1,iglob) - tx*weight
+ accel(2,iglob)=accel(2,iglob) - ty*weight
+ accel(3,iglob)=accel(3,iglob) - tz*weight
+
+ enddo
+ enddo
+ endif
+ enddo
+
+ ! bottom (zmin)
+ do ispec2D=1,NSPEC2D_BOTTOM
+
+ ispec=ibelm_bottom(ispec2D)
+
+ if (ispec_is_inner(ispec) .eqv. phase_is_inner) then
+
+ k=1
+ do j=1,NGLLY
+ do i=1,NGLLX
+
+ iglob=ibool(i,j,k,ispec)
+
+ vx=veloc(1,iglob)
+ vy=veloc(2,iglob)
+ vz=veloc(3,iglob)
+
+ nx=normal_bottom(1,i,j,ispec2D)
+ ny=normal_bottom(2,i,j,ispec2D)
+ nz=normal_bottom(3,i,j,ispec2D)
+
+ vn=vx*nx+vy*ny+vz*nz
+
+ tx=rho_vp(i,j,k,ispec)*vn*nx+rho_vs(i,j,k,ispec)*(vx-vn*nx)
+ ty=rho_vp(i,j,k,ispec)*vn*ny+rho_vs(i,j,k,ispec)*(vy-vn*ny)
+ tz=rho_vp(i,j,k,ispec)*vn*nz+rho_vs(i,j,k,ispec)*(vz-vn*nz)
+
+ weight=jacobian2D_bottom(i,j,ispec2D)*wgllwgll_xy(i,j)
+
+ accel(1,iglob)=accel(1,iglob) - tx*weight
+ accel(2,iglob)=accel(2,iglob) - ty*weight
+ accel(3,iglob)=accel(3,iglob) - tz*weight
+
+ enddo
+ enddo
+ endif
+ enddo
+
+ endif ! end of Stacey conditions
+
+
! adding source
do isource = 1,NSOURCES
if (ispec_is_inner(ispec_selected_source(isource)) .eqv. phase_is_inner) then
- if(USE_FORCE_POINT_SOURCE) then
+ if(USE_FORCE_POINT_SOURCE) then
-! add the source (only if this proc carries the source)
- if(myrank == islice_selected_source(isource)) then
+ ! add the source (only if this proc carries the source)
+ if(myrank == islice_selected_source(isource)) then
+
+ iglob = ibool(nint(xi_source(isource)), &
+ nint(eta_source(isource)), &
+ nint(gamma_source(isource)), &
+ ispec_selected_source(isource))
+ f0 = hdur(isource) !! using hdur as a FREQUENCY just to avoid changing CMTSOLUTION file format
+ t0 = 1.2d0/f0
+
+ if (it == 1 .and. myrank == 0) then
+ print *,'using a source of dominant frequency ',f0
+ print *,'lambda_S at dominant frequency = ',3000./sqrt(3.)/f0
+ print *,'lambda_S at highest significant frequency = ',3000./sqrt(3.)/(2.5*f0)
+ endif
+
+ ! we use nu_source(:,3) here because we want a source normal to the surface.
+ ! This is the expression of a Ricker; should be changed according maybe to the Par_file.
+ !accel(:,iglob) = accel(:,iglob) + &
+ ! sngl(nu_source(:,3,isource) * 10000000.d0 * (1.d0-2.d0*PI*PI*f0*f0*(dble(it-1)*DT-t0)*(dble(it-1)*DT-t0)) * &
+ ! exp(-PI*PI*f0*f0*(dble(it-1)*DT-t0)*(dble(it-1)*DT-t0)))
+ accel(:,iglob) = accel(:,iglob) + &
+ sngl(nu_source(:,3,isource) * 1.d10 * (1.d0-2.d0*PI*PI*f0*f0*(dble(it-1)*DT-t0)*(dble(it-1)*DT-t0)) * &
+ exp(-PI*PI*f0*f0*(dble(it-1)*DT-t0)*(dble(it-1)*DT-t0)))
+
+ endif
- iglob = ibool(nint(xi_source(isource)), &
- nint(eta_source(isource)), &
- nint(gamma_source(isource)), &
- ispec_selected_source(isource))
- f0 = hdur(isource) !! using hdur as a FREQUENCY just to avoid changing CMTSOLUTION file format
- t0 = 1.2d0/f0
+ else
+ ! add the source (only if this proc carries the source)
+ if(myrank == islice_selected_source(isource)) then
+
+ stf = comp_source_time_function(dble(it-1)*DT-t0-t_cmt(isource),hdur_gaussian(isource))
- if (it == 1 .and. myrank == 0) then
- print *,'using a source of dominant frequency ',f0
- print *,'lambda_S at dominant frequency = ',3000./sqrt(3.)/f0
- print *,'lambda_S at highest significant frequency = ',3000./sqrt(3.)/(2.5*f0)
- endif
+ ! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ stf_used = sngl(stf)
+ else
+ stf_used = stf
+ endif
- ! we use nu_source(:,3) here because we want a source normal to the surface.
- ! This is the expression of a Ricker; should be changed according maybe to the Par_file.
- !accel(:,iglob) = accel(:,iglob) + &
- ! sngl(nu_source(:,3,isource) * 10000000.d0 * (1.d0-2.d0*PI*PI*f0*f0*(dble(it-1)*DT-t0)*(dble(it-1)*DT-t0)) * &
- ! exp(-PI*PI*f0*f0*(dble(it-1)*DT-t0)*(dble(it-1)*DT-t0)))
- accel(:,iglob) = accel(:,iglob) + &
- sngl(nu_source(:,3,isource) * 1.d10 * (1.d0-2.d0*PI*PI*f0*f0*(dble(it-1)*DT-t0)*(dble(it-1)*DT-t0)) * &
- exp(-PI*PI*f0*f0*(dble(it-1)*DT-t0)*(dble(it-1)*DT-t0)))
+ ! add source array
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec_selected_source(isource))
+ accel(:,iglob) = accel(:,iglob) + sourcearrays(isource,:,i,j,k)*stf_used
+ enddo
+ enddo
+ enddo
- endif
+ endif
+ endif
+
endif
+
+ enddo
- endif
-
- enddo
-
end subroutine compute_forces_with_Deville
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/constants.h.in 2009-07-10 16:39:24 UTC (rev 15452)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/constants.h.in 2009-07-10 21:02:53 UTC (rev 15453)
@@ -129,9 +129,6 @@
! nlegoff -- Variables that should be read/computed elsewhere.
! Temporarily declared here.
!------------------------------------------------------
-! whether or not an external mesh is used (provided by CUBIT for example)
- logical, parameter :: USE_EXTERNAL_MESH = .true.
-
! no lagrange interpolation on seismograms (we take the value on one NGLL point)
logical, parameter :: FASTER_RECEIVERS_POINTS_ONLY = .false.
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/create_movie_shakemap_AVS_DX_GMT.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/create_movie_shakemap_AVS_DX_GMT.f90 2009-07-10 16:39:24 UTC (rev 15452)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/create_movie_shakemap_AVS_DX_GMT.f90 2009-07-10 21:02:53 UTC (rev 15453)
@@ -157,20 +157,13 @@
endif
print *
- if (USE_EXTERNAL_MESH) then
- NPROC = 1
- if(USE_HIGHRES_FOR_MOVIES) then
- ilocnum = NSPEC_SURFACE_EXT_MESH*NGLLSQUARE
- else
- ilocnum = NSPEC_SURFACE_EXT_MESH*NGNOD2D_AVS_DX
- endif
+ NPROC = 1
+ if(USE_HIGHRES_FOR_MOVIES) then
+ ilocnum = NSPEC_SURFACE_EXT_MESH*NGLLSQUARE
else
- if(USE_HIGHRES_FOR_MOVIES) then
- ilocnum = NGLLSQUARE*NEX_PER_PROC_XI*NEX_PER_PROC_ETA
- else
- ilocnum = NGNOD2D_AVS_DX*NEX_PER_PROC_XI*NEX_PER_PROC_ETA
- endif
+ ilocnum = NSPEC_SURFACE_EXT_MESH*NGNOD2D_AVS_DX
endif
+
allocate(store_val_x(ilocnum,0:NPROC-1))
allocate(store_val_y(ilocnum,0:NPROC-1))
allocate(store_val_z(ilocnum,0:NPROC-1))
@@ -261,20 +254,12 @@
endif
! define the total number of elements at the surface
- if (USE_EXTERNAL_MESH) then
- if(USE_HIGHRES_FOR_MOVIES) then
- nspectot_AVS_max = NSPEC_SURFACE_EXT_MESH * (NGLLX-1) * (NGLLY-1)
- else
- nspectot_AVS_max = NSPEC_SURFACE_EXT_MESH
- endif
+ if(USE_HIGHRES_FOR_MOVIES) then
+ nspectot_AVS_max = NSPEC_SURFACE_EXT_MESH * (NGLLX-1) * (NGLLY-1)
else
- if(USE_HIGHRES_FOR_MOVIES) then
- nspectot_AVS_max = NEX_XI * NEX_ETA * (NGLLX-1) * (NGLLY-1)
- else
- nspectot_AVS_max = NEX_XI * NEX_ETA
- endif
+ nspectot_AVS_max = NSPEC_SURFACE_EXT_MESH
endif
-
+
print *
print *,'there are a total of ',nspectot_AVS_max,' elements at the surface'
print *
@@ -389,12 +374,8 @@
endif
endif
else
- if (USE_EXTERNAL_MESH) then
display(i,j) = sqrt(vectorz**2+vectory**2+vectorx**2)
- else
- display(i,j) = vectorz
- endif
- endif
+ endif
enddo
enddo
@@ -479,11 +460,7 @@
endif
endif
else
- if (USE_EXTERNAL_MESH) then
field_display(ilocnum+ieoff) =sqrt(vectorz**2+vectory**2+vectorx**2)
- else
- field_display(ilocnum+ieoff) = dble(vectorz)
- endif
endif
enddo
@@ -500,12 +477,8 @@
!--- sort the list based upon coordinates to get rid of multiples
print *,'sorting list of points'
- if (USE_EXTERNAL_MESH) then
- call get_global_AVS(nspectot_AVS_max,xp,yp,zp,iglob,loc,ifseg,nglob,npointot, &
- dble(minval(store_val_x(:,0))),dble(maxval(store_val_x(:,0))))
- else
- call get_global_AVS(nspectot_AVS_max,xp,yp,zp,iglob,loc,ifseg,nglob,npointot,UTM_X_MIN,UTM_X_MAX)
- endif
+ call get_global_AVS(nspectot_AVS_max,xp,yp,zp,iglob,loc,ifseg,nglob,npointot, &
+ dble(minval(store_val_x(:,0))),dble(maxval(store_val_x(:,0))))
!--- print total number of points found
print *
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/create_regions_mesh.f90 2009-07-10 16:39:24 UTC (rev 15452)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/create_regions_mesh.f90 2009-07-10 21:02:53 UTC (rev 15453)
@@ -23,1144 +23,17 @@
!
!=====================================================================
- subroutine create_regions_mesh(xgrid,ygrid,zgrid,ibool,idoubling, &
- xstore,ystore,zstore,npx,npy,iproc_xi,iproc_eta,nspec, &
- volume_local,area_local_bottom,area_local_top, &
- NGLOB_AB,npointot, &
- NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM,NER_SEDIM,NER, &
- NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
- NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
- HARVARD_3D_GOCAD_MODEL,NPROC_XI,NPROC_ETA,NSPEC2D_A_XI,NSPEC2D_B_XI, &
- NSPEC2D_A_ETA,NSPEC2D_B_ETA, &
- myrank,LOCAL_PATH,UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK,UTM_PROJECTION_ZONE, &
- HAUKSSON_REGIONAL_MODEL,OCEANS, &
- VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM, &
- IMPOSE_MINIMUM_VP_GOCAD,THICKNESS_TAPER_BLOCK_HR,THICKNESS_TAPER_BLOCK_MR,MOHO_MAP_LUPEI, &
- ANISOTROPY,SAVE_MESH_FILES,SUPPRESS_UTM_PROJECTION, &
- ORIG_LAT_TOPO,ORIG_LONG_TOPO,DEGREES_PER_CELL_TOPO,NX_TOPO,NY_TOPO,USE_REGULAR_MESH)
-! create the different regions of the mesh
-
- implicit none
-
- include "constants.h"
-
-! number of spectral elements in each block
- integer nspec
-
- integer NEX_PER_PROC_XI,NEX_PER_PROC_ETA,UTM_PROJECTION_ZONE
- integer NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM,NER_SEDIM,NER
-
- integer NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP
-
- integer NPROC_XI,NPROC_ETA,NSPEC2D_A_XI,NSPEC2D_B_XI
- integer NSPEC2D_A_ETA,NSPEC2D_B_ETA
- integer NX_TOPO,NY_TOPO
-
- integer npx,npy
- integer npointot
-
- logical HARVARD_3D_GOCAD_MODEL,HAUKSSON_REGIONAL_MODEL
- logical OCEANS,IMPOSE_MINIMUM_VP_GOCAD,USE_REGULAR_MESH
- logical MOHO_MAP_LUPEI,ANISOTROPY,SAVE_MESH_FILES,SUPPRESS_UTM_PROJECTION
-
- double precision UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK
- double precision VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM
- double precision horiz_size,vert_size,THICKNESS_TAPER_BLOCK_HR,THICKNESS_TAPER_BLOCK_MR
- double precision ORIG_LAT_TOPO,ORIG_LONG_TOPO,DEGREES_PER_CELL_TOPO
-
- character(len=150) LOCAL_PATH
-
-! arrays with the mesh
- double precision, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: xstore,ystore,zstore
-
- double precision xstore_local(NGLLX,NGLLY,NGLLZ)
- double precision ystore_local(NGLLX,NGLLY,NGLLZ)
- double precision zstore_local(NGLLX,NGLLY,NGLLZ)
-
- double precision xgrid(0:2*NER,0:2*NEX_PER_PROC_XI,0:2*NEX_PER_PROC_ETA)
- double precision ygrid(0:2*NER,0:2*NEX_PER_PROC_XI,0:2*NEX_PER_PROC_ETA)
- double precision zgrid(0:2*NER,0:2*NEX_PER_PROC_XI,0:2*NEX_PER_PROC_ETA)
-
- double precision xmesh,ymesh,zmesh
-
- integer ibool(NGLLX,NGLLY,NGLLZ,nspec)
-
-! use integer array to store topography values
- integer icornerlat,icornerlong
- double precision lat,long,elevation
- double precision long_corner,lat_corner,ratio_xi,ratio_eta
- integer itopo_bathy(NX_TOPO,NY_TOPO)
-
-! auxiliary variables to generate the mesh
- integer ix,iy,iz,ir,ir1,ir2,dir
- integer ix1,ix2,dix,iy1,iy2,diy
- integer iax,iay,iar
- integer isubregion,nsubregions,doubling_index
-
-! Gauss-Lobatto-Legendre points and weights of integration
- double precision, dimension(:), allocatable :: xigll,yigll,zigll,wxgll,wygll,wzgll
-
-! 3D shape functions and their derivatives
- double precision, dimension(:,:,:,:), allocatable :: shape3D
- double precision, dimension(:,:,:,:,:), allocatable :: dershape3D
-
-! 2D shape functions and their derivatives
- double precision, dimension(:,:,:), allocatable :: shape2D_x,shape2D_y,shape2D_bottom,shape2D_top
- double precision, dimension(:,:,:,:), allocatable :: dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top
-
-! topology of the elements
- integer iaddx(NGNOD)
- integer iaddy(NGNOD)
- integer iaddz(NGNOD)
-
- double precision xelm(NGNOD)
- double precision yelm(NGNOD)
- double precision zelm(NGNOD)
-
-! parameters needed to store the radii of the grid points
-! in the spherically symmetric Earth
- integer idoubling(nspec)
-
-! for model density
- real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: rhostore,kappastore,mustore
- real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: c11store,c12store,c13store,c14store,c15store,c16store,&
- c22store,c23store,c24store,c25store,c26store,c33store,c34store,c35store,c36store,c44store,c45store,c46store,&
- c55store,c56store,c66store
-
-! the jacobian
- real(kind=CUSTOM_REAL) jacobianl
-
-! boundary locator
- logical, dimension(:,:), allocatable :: iboun
-
-! arrays with mesh parameters
- real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: xixstore,xiystore,xizstore, &
- etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore,jacobianstore
-
-! mass matrix and bathymetry for ocean load
- integer ix_oceans,iy_oceans,iz_oceans,ispec_oceans
- integer ispec2D_ocean_bottom
- integer nglob_oceans
- double precision xval,yval
- double precision height_oceans
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_ocean_load
-
-! proc numbers for MPI
- integer myrank
-
-! check area and volume of the final mesh
- double precision weight
- double precision area_local_bottom,area_local_top
- double precision volume_local
-
-! variables for creating array ibool (some arrays also used for AVS or DX files)
- integer, dimension(:), allocatable :: iglob,locval
- logical, dimension(:), allocatable :: ifseg
- double precision, dimension(:), allocatable :: xp,yp,zp
-
- integer nglob,NGLOB_AB
- integer ieoff,ilocnum
-
-! mass matrix
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass
-
-! boundary parameters locator
- integer, dimension(:), allocatable :: ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top
-
-! ---- Moho Vars here ------
-! Moho boundary locator
- integer, dimension(:), allocatable :: ibelm_moho_top, ibelm_moho_bot
- logical, dimension(:), allocatable :: is_moho_top, is_moho_bot
-
-! 2-D jacobian and normals
- real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: jacobian2D_moho
- real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: normal_moho
-
-! number of elements on the boundaries
- integer nspec_moho_top, nspec_moho_bottom
-! ---------------------------
-
-! 2-D jacobians and normals
- real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: &
- jacobian2D_xmin,jacobian2D_xmax, &
- jacobian2D_ymin,jacobian2D_ymax,jacobian2D_bottom,jacobian2D_top
-
- real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: &
- normal_xmin,normal_xmax,normal_ymin,normal_ymax,normal_bottom,normal_top
-
-! MPI cut-planes parameters along xi and along eta
- logical, dimension(:,:), allocatable :: iMPIcut_xi,iMPIcut_eta
-
-! name of the database file
- character(len=150) prname
-
-! number of elements on the boundaries
- integer nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax
-
- integer i,j,k,ia,ispec,iglobnum,itype_element
- integer iproc_xi,iproc_eta
-
- double precision rho,vp,vs
- double precision c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,c33,c34,c35,c36,c44,c45,c46,c55,c56,c66
-
-! for the Harvard 3-D basin model
- double precision vp_block_gocad_MR(0:NX_GOCAD_MR-1,0:NY_GOCAD_MR-1,0:NZ_GOCAD_MR-1)
- double precision vp_block_gocad_HR(0:NX_GOCAD_HR-1,0:NY_GOCAD_HR-1,0:NZ_GOCAD_HR-1)
- integer irecord,nrecord,i_vp
- character(len=150) BASIN_MODEL_3D_MEDIUM_RES_FILE,BASIN_MODEL_3D_HIGH_RES_FILE
-
-! for the harvard 3D salton sea model
- real :: vp_st_gocad(GOCAD_ST_NU,GOCAD_ST_NV,GOCAD_ST_NW)
- double precision :: umesh, vmesh, wmesh, vp_st, vs_st, rho_st
-
-! for Hauksson's model
- double precision, dimension(NLAYERS_HAUKSSON,NGRID_NEW_HAUKSSON,NGRID_NEW_HAUKSSON) :: vp_hauksson,vs_hauksson
- integer ilayer
- character(len=150 ) HAUKSSON_REGIONAL_MODEL_FILE
-
-! Stacey put back
-! indices for Clayton-Engquist absorbing conditions
- integer, dimension(:,:), allocatable :: nimin,nimax,njmin,njmax,nkmin_xi,nkmin_eta
- real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: rho_vp,rho_vs
-
-! flag indicating whether point is in the sediments
- logical point_is_in_sediments
- logical, dimension(:,:,:,:), allocatable :: flag_sediments
- logical, dimension(:), allocatable :: not_fully_in_bedrock
-
-! mask to sort ibool
- integer, dimension(:), allocatable :: mask_ibool
- integer, dimension(:,:,:,:), allocatable :: copy_ibool_ori
- integer :: inumber
-
-! **************
-
-! create the name for the database of the current slide and region
- call create_name_database(prname,myrank,LOCAL_PATH)
-
-! Gauss-Lobatto-Legendre points of integration
- allocate(xigll(NGLLX))
- allocate(yigll(NGLLY))
- allocate(zigll(NGLLZ))
-
-! Gauss-Lobatto-Legendre weights of integration
- allocate(wxgll(NGLLX))
- allocate(wygll(NGLLY))
- allocate(wzgll(NGLLZ))
-
-! 3D shape functions and their derivatives
- allocate(shape3D(NGNOD,NGLLX,NGLLY,NGLLZ))
- allocate(dershape3D(NDIM,NGNOD,NGLLX,NGLLY,NGLLZ))
-
-! 2D shape functions and their derivatives
- allocate(shape2D_x(NGNOD2D,NGLLY,NGLLZ))
- allocate(shape2D_y(NGNOD2D,NGLLX,NGLLZ))
- allocate(shape2D_bottom(NGNOD2D,NGLLX,NGLLY))
- allocate(shape2D_top(NGNOD2D,NGLLX,NGLLY))
- allocate(dershape2D_x(NDIM2D,NGNOD2D,NGLLY,NGLLZ))
- allocate(dershape2D_y(NDIM2D,NGNOD2D,NGLLX,NGLLZ))
- allocate(dershape2D_bottom(NDIM2D,NGNOD2D,NGLLX,NGLLY))
- allocate(dershape2D_top(NDIM2D,NGNOD2D,NGLLX,NGLLY))
-
-! array with model density
- allocate(rhostore(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(kappastore(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(mustore(NGLLX,NGLLY,NGLLZ,nspec))
-
-! array with anisotropy
- allocate(c11store(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(c12store(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(c13store(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(c14store(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(c15store(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(c16store(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(c22store(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(c23store(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(c24store(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(c25store(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(c26store(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(c33store(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(c34store(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(c35store(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(c36store(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(c44store(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(c45store(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(c46store(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(c55store(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(c56store(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(c66store(NGLLX,NGLLY,NGLLZ,nspec))
-
-! Stacey
- allocate(rho_vp(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(rho_vs(NGLLX,NGLLY,NGLLZ,nspec))
-
-! flag indicating whether point is in the sediments
- allocate(flag_sediments(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(not_fully_in_bedrock(nspec))
-
-! boundary locator
- allocate(iboun(6,nspec))
-
-! arrays with mesh parameters
- allocate(xixstore(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(xiystore(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(xizstore(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(etaxstore(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(etaystore(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(etazstore(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(gammaxstore(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(gammaystore(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(gammazstore(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(jacobianstore(NGLLX,NGLLY,NGLLZ,nspec))
-
-! boundary parameters locator
- allocate(ibelm_xmin(NSPEC2DMAX_XMIN_XMAX))
- allocate(ibelm_xmax(NSPEC2DMAX_XMIN_XMAX))
- allocate(ibelm_ymin(NSPEC2DMAX_YMIN_YMAX))
- allocate(ibelm_ymax(NSPEC2DMAX_YMIN_YMAX))
- allocate(ibelm_bottom(NSPEC2D_BOTTOM))
- allocate(ibelm_top(NSPEC2D_TOP))
-
-! 2-D jacobians and normals
- allocate(jacobian2D_xmin(NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX))
- allocate(jacobian2D_xmax(NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX))
- allocate(jacobian2D_ymin(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX))
- allocate(jacobian2D_ymax(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX))
- allocate(jacobian2D_bottom(NGLLX,NGLLY,NSPEC2D_BOTTOM))
- allocate(jacobian2D_top(NGLLX,NGLLY,NSPEC2D_TOP))
-
- allocate(normal_xmin(NDIM,NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX))
- allocate(normal_xmax(NDIM,NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX))
- allocate(normal_ymin(NDIM,NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX))
- allocate(normal_ymax(NDIM,NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX))
- allocate(normal_bottom(NDIM,NGLLX,NGLLY,NSPEC2D_BOTTOM))
- allocate(normal_top(NDIM,NGLLX,NGLLY,NSPEC2D_TOP))
-
-! Moho boundary parameters, 2-D jacobians and normals
- if (SAVE_MOHO_MESH) then
- allocate(ibelm_moho_top(NSPEC2D_BOTTOM))
- allocate(ibelm_moho_bot(NSPEC2D_BOTTOM))
- allocate(is_moho_top(nspec))
- allocate(is_moho_bot(nspec))
- is_moho_top = .false.
- is_moho_bot = .false.
- nspec_moho_top = 0
- nspec_moho_bottom = 0
- allocate(jacobian2D_moho(NGLLX,NGLLY,NSPEC2D_BOTTOM))
- allocate(normal_moho(NDIM,NGLLX,NGLLY,NSPEC2D_BOTTOM))
- endif
-
-
-! Stacey put back
- allocate(nimin(2,NSPEC2DMAX_YMIN_YMAX))
- allocate(nimax(2,NSPEC2DMAX_YMIN_YMAX))
- allocate(njmin(2,NSPEC2DMAX_XMIN_XMAX))
- allocate(njmax(2,NSPEC2DMAX_XMIN_XMAX))
- allocate(nkmin_xi(2,NSPEC2DMAX_XMIN_XMAX))
- allocate(nkmin_eta(2,NSPEC2DMAX_YMIN_YMAX))
-
-! MPI cut-planes parameters along xi and along eta
- allocate(iMPIcut_xi(2,nspec))
- allocate(iMPIcut_eta(2,nspec))
-
-! set up coordinates of the Gauss-Lobatto-Legendre points
- call zwgljd(xigll,wxgll,NGLLX,GAUSSALPHA,GAUSSBETA)
- call zwgljd(yigll,wygll,NGLLY,GAUSSALPHA,GAUSSBETA)
- call zwgljd(zigll,wzgll,NGLLZ,GAUSSALPHA,GAUSSBETA)
-
-! if number of points is odd, the middle abscissa is exactly zero
- if(mod(NGLLX,2) /= 0) xigll((NGLLX-1)/2+1) = ZERO
- if(mod(NGLLY,2) /= 0) yigll((NGLLY-1)/2+1) = ZERO
- if(mod(NGLLZ,2) /= 0) zigll((NGLLZ-1)/2+1) = ZERO
-
-! get the 3-D shape functions
- call get_shape3D(myrank,shape3D,dershape3D,xigll,yigll,zigll)
-
-! get the 2-D shape functions
- call get_shape2D(myrank,shape2D_x,dershape2D_x,yigll,zigll,NGLLY,NGLLZ)
- call get_shape2D(myrank,shape2D_y,dershape2D_y,xigll,zigll,NGLLX,NGLLZ)
- call get_shape2D(myrank,shape2D_bottom,dershape2D_bottom,xigll,yigll,NGLLX,NGLLY)
- call get_shape2D(myrank,shape2D_top,dershape2D_top,xigll,yigll,NGLLX,NGLLY)
-
-! allocate memory for arrays
- allocate(iglob(npointot))
- allocate(locval(npointot))
- allocate(ifseg(npointot))
- allocate(xp(npointot))
- allocate(yp(npointot))
- allocate(zp(npointot))
-
-!--- read Hauksson's model
- if(HAUKSSON_REGIONAL_MODEL) then
- call get_value_string(HAUKSSON_REGIONAL_MODEL_FILE, &
- 'model.HAUKSSON_REGIONAL_MODEL_FILE', &
- 'DATA/hauksson_model/hauksson_final_grid_smooth.dat')
-! call get_value_string(HAUKSSON_REGIONAL_MODEL_FILE, &
-! 'model.HAUKSSON_REGIONAL_MODEL_FILE', &
-! 'DATA/lin_model/lin_final_grid_smooth.dat')
- open(unit=14,file=HAUKSSON_REGIONAL_MODEL_FILE,status='old',action='read')
- do iy = 1,NGRID_NEW_HAUKSSON
- do ix = 1,NGRID_NEW_HAUKSSON
- read(14,*) (vp_hauksson(ilayer,ix,iy),ilayer=1,NLAYERS_HAUKSSON), &
- (vs_hauksson(ilayer,ix,iy),ilayer=1,NLAYERS_HAUKSSON)
- enddo
- enddo
- close(14)
- vp_hauksson(:,:,:) = vp_hauksson(:,:,:) * 1000.d0
- vs_hauksson(:,:,:) = vs_hauksson(:,:,:) * 1000.d0
- endif
-
-!--- read the Harvard 3-D basin model
- if(HARVARD_3D_GOCAD_MODEL) then
-
-! read medium-resolution model
-
-! initialize array to undefined values everywhere
- vp_block_gocad_MR(:,:,:) = 20000.
-
-! read Vp from extracted text file
- call get_value_string(BASIN_MODEL_3D_MEDIUM_RES_FILE, &
- 'model.BASIN_MODEL_3D_MEDIUM_RES_FILE', &
- 'DATA/la_3D_block_harvard/la_3D_medium_res/LA_MR_voxet_extracted.txt')
- open(unit=27,file=BASIN_MODEL_3D_MEDIUM_RES_FILE,status='old',action='read')
- read(27,*) nrecord
- do irecord = 1,nrecord
- read(27,*) ix,iy,iz,i_vp
- if(ix<0 .or. ix>NX_GOCAD_MR-1 .or. iy<0 .or. iy>NY_GOCAD_MR-1 .or. iz<0 .or. iz>NZ_GOCAD_MR-1) &
- stop 'wrong array index read in Gocad medium-resolution file'
- vp_block_gocad_MR(ix,iy,iz) = dble(i_vp)
- enddo
- close(27)
-
-! read high-resolution model
-
-! initialize array to undefined values everywhere
- vp_block_gocad_HR(:,:,:) = 20000.
-
-! read Vp from extracted text file
- call get_value_string(BASIN_MODEL_3D_HIGH_RES_FILE, &
- 'model.BASIN_MODEL_3D_HIGH_RES_FILE', &
- 'DATA/la_3D_block_harvard/la_3D_high_res/LA_HR_voxet_extracted.txt')
- open(unit=27,file=BASIN_MODEL_3D_HIGH_RES_FILE,status='old',action='read')
- read(27,*) nrecord
- do irecord = 1,nrecord
- read(27,*) ix,iy,iz,i_vp
- if(ix<0 .or. ix>NX_GOCAD_HR-1 .or. iy<0 .or. iy>NY_GOCAD_HR-1 .or. iz<0 .or. iz>NZ_GOCAD_HR-1) &
- stop 'wrong array index read in Gocad high-resolution file'
- vp_block_gocad_HR(ix,iy,iz) = dble(i_vp)
- enddo
- close(27)
-
-! read Salton Trough model
- call read_salton_sea_model(vp_st_gocad)
-
- endif
-
-!--- apply heuristic rule to modify doubling regions to balance angles
-
- if(APPLY_HEURISTIC_RULE .and. .not. USE_REGULAR_MESH) then
-
-! define number of subregions affected by heuristic rule in doubling regions
- nsubregions = 8
-
- do isubregion = 1,nsubregions
-
-! define shape of elements for heuristic
- call define_subregions_heuristic(myrank,isubregion,iaddx,iaddy,iaddz, &
- ix1,ix2,dix,iy1,iy2,diy,ir1,ir2,dir,iax,iay,iar, &
- itype_element,npx,npy, &
- NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM)
-
-! loop on all the mesh points in current subregion
- do ir = ir1,ir2,dir
- do iy = iy1,iy2,diy
- do ix = ix1,ix2,dix
-
-! this heuristic rule is only valid for 8-node elements
-! it would not work in the case of 27 nodes
-
-!----
- if(itype_element == ITYPE_UNUSUAL_1) then
-
-! side 1
- horiz_size = xgrid(ir+iar*iaddz(2),ix+iax*iaddx(2),iy+iay*iaddy(2)) &
- - xgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1))
- xgrid(ir+iar*iaddz(5),ix+iax*iaddx(5),iy+iay*iaddy(5)) = &
- xgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1)) + horiz_size * MAGIC_RATIO
-
- vert_size = zgrid(ir+iar*iaddz(5),ix+iax*iaddx(5),iy+iay*iaddy(5)) &
- - zgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1))
- zgrid(ir+iar*iaddz(5),ix+iax*iaddx(5),iy+iay*iaddy(5)) = &
- zgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1)) + vert_size * MAGIC_RATIO / 0.50
-
-! side 2
- horiz_size = xgrid(ir+iar*iaddz(3),ix+iax*iaddx(3),iy+iay*iaddy(3)) &
- - xgrid(ir+iar*iaddz(4),ix+iax*iaddx(4),iy+iay*iaddy(4))
- xgrid(ir+iar*iaddz(8),ix+iax*iaddx(8),iy+iay*iaddy(8)) = &
- xgrid(ir+iar*iaddz(4),ix+iax*iaddx(4),iy+iay*iaddy(4)) + horiz_size * MAGIC_RATIO
-
- vert_size = zgrid(ir+iar*iaddz(8),ix+iax*iaddx(8),iy+iay*iaddy(8)) &
- - zgrid(ir+iar*iaddz(4),ix+iax*iaddx(4),iy+iay*iaddy(4))
- zgrid(ir+iar*iaddz(8),ix+iax*iaddx(8),iy+iay*iaddy(8)) = &
- zgrid(ir+iar*iaddz(4),ix+iax*iaddx(4),iy+iay*iaddy(4)) + vert_size * MAGIC_RATIO / 0.50
-
-!----
- else if(itype_element == ITYPE_UNUSUAL_1p) then
-
-! side 1
- horiz_size = xgrid(ir+iar*iaddz(2),ix+iax*iaddx(2),iy+iay*iaddy(2)) &
- - xgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1))
- xgrid(ir+iar*iaddz(6),ix+iax*iaddx(6),iy+iay*iaddy(6)) = &
- xgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1)) + horiz_size * (1. - MAGIC_RATIO)
-
- vert_size = zgrid(ir+iar*iaddz(5),ix+iax*iaddx(5),iy+iay*iaddy(5)) &
- - zgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1))
- zgrid(ir+iar*iaddz(6),ix+iax*iaddx(6),iy+iay*iaddy(6)) = &
- zgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1)) + vert_size * MAGIC_RATIO / 0.50
-
-! side 2
- horiz_size = xgrid(ir+iar*iaddz(3),ix+iax*iaddx(3),iy+iay*iaddy(3)) &
- - xgrid(ir+iar*iaddz(4),ix+iax*iaddx(4),iy+iay*iaddy(4))
- xgrid(ir+iar*iaddz(7),ix+iax*iaddx(7),iy+iay*iaddy(7)) = &
- xgrid(ir+iar*iaddz(4),ix+iax*iaddx(4),iy+iay*iaddy(4)) + horiz_size * (1. - MAGIC_RATIO)
-
- vert_size = zgrid(ir+iar*iaddz(8),ix+iax*iaddx(8),iy+iay*iaddy(8)) &
- - zgrid(ir+iar*iaddz(4),ix+iax*iaddx(4),iy+iay*iaddy(4))
- zgrid(ir+iar*iaddz(7),ix+iax*iaddx(7),iy+iay*iaddy(7)) = &
- zgrid(ir+iar*iaddz(4),ix+iax*iaddx(4),iy+iay*iaddy(4)) + vert_size * MAGIC_RATIO / 0.50
-
-!----
- else if(itype_element == ITYPE_UNUSUAL_4) then
-
-! side 1
- horiz_size = ygrid(ir+iar*iaddz(3),ix+iax*iaddx(3),iy+iay*iaddy(3)) &
- - ygrid(ir+iar*iaddz(2),ix+iax*iaddx(2),iy+iay*iaddy(2))
- ygrid(ir+iar*iaddz(7),ix+iax*iaddx(7),iy+iay*iaddy(7)) = &
- ygrid(ir+iar*iaddz(2),ix+iax*iaddx(2),iy+iay*iaddy(2)) + horiz_size * (1. - MAGIC_RATIO)
-
- vert_size = zgrid(ir+iar*iaddz(6),ix+iax*iaddx(6),iy+iay*iaddy(6)) &
- - zgrid(ir+iar*iaddz(2),ix+iax*iaddx(2),iy+iay*iaddy(2))
- zgrid(ir+iar*iaddz(7),ix+iax*iaddx(7),iy+iay*iaddy(7)) = &
- zgrid(ir+iar*iaddz(2),ix+iax*iaddx(2),iy+iay*iaddy(2)) + vert_size * MAGIC_RATIO / 0.50
-
-! side 2
- horiz_size = ygrid(ir+iar*iaddz(4),ix+iax*iaddx(4),iy+iay*iaddy(4)) &
- - ygrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1))
- ygrid(ir+iar*iaddz(8),ix+iax*iaddx(8),iy+iay*iaddy(8)) = &
- ygrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1)) + horiz_size * (1. - MAGIC_RATIO)
-
- vert_size = zgrid(ir+iar*iaddz(5),ix+iax*iaddx(5),iy+iay*iaddy(5)) &
- - zgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1))
- zgrid(ir+iar*iaddz(8),ix+iax*iaddx(8),iy+iay*iaddy(8)) = &
- zgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1)) + vert_size * MAGIC_RATIO / 0.50
-
-!----
- else if(itype_element == ITYPE_UNUSUAL_4p) then
-
-! side 1
- horiz_size = ygrid(ir+iar*iaddz(3),ix+iax*iaddx(3),iy+iay*iaddy(3)) &
- - ygrid(ir+iar*iaddz(2),ix+iax*iaddx(2),iy+iay*iaddy(2))
- ygrid(ir+iar*iaddz(6),ix+iax*iaddx(6),iy+iay*iaddy(6)) = &
- ygrid(ir+iar*iaddz(2),ix+iax*iaddx(2),iy+iay*iaddy(2)) + horiz_size * MAGIC_RATIO
-
- vert_size = zgrid(ir+iar*iaddz(6),ix+iax*iaddx(6),iy+iay*iaddy(6)) &
- - zgrid(ir+iar*iaddz(2),ix+iax*iaddx(2),iy+iay*iaddy(2))
- zgrid(ir+iar*iaddz(6),ix+iax*iaddx(6),iy+iay*iaddy(6)) = &
- zgrid(ir+iar*iaddz(2),ix+iax*iaddx(2),iy+iay*iaddy(2)) + vert_size * MAGIC_RATIO / 0.50
-
-! side 2
- horiz_size = ygrid(ir+iar*iaddz(4),ix+iax*iaddx(4),iy+iay*iaddy(4)) &
- - ygrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1))
- ygrid(ir+iar*iaddz(5),ix+iax*iaddx(5),iy+iay*iaddy(5)) = &
- ygrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1)) + horiz_size * MAGIC_RATIO
-
- vert_size = zgrid(ir+iar*iaddz(5),ix+iax*iaddx(5),iy+iay*iaddy(5)) &
- - zgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1))
- zgrid(ir+iar*iaddz(5),ix+iax*iaddx(5),iy+iay*iaddy(5)) = &
- zgrid(ir+iar*iaddz(1),ix+iax*iaddx(1),iy+iay*iaddy(1)) + vert_size * MAGIC_RATIO / 0.50
-
- endif
-
- enddo
- enddo
- enddo
-
- enddo
-
- endif
-
-!---
-
-! generate the elements in all the regions of the mesh
- ispec = 0
-
-! define number of subregions in the mesh
- if(USE_REGULAR_MESH) then
- nsubregions = 2
- else
- if(NER_SEDIM > 1) then
- nsubregions = 30
- else
- nsubregions = 29
- endif
- endif
-
- do isubregion = 1,nsubregions
-
-! define shape of elements
- call define_subregions(myrank,isubregion,iaddx,iaddy,iaddz, &
- ix1,ix2,dix,iy1,iy2,diy,ir1,ir2,dir,iax,iay,iar, &
- doubling_index,npx,npy, &
- NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM,NER_SEDIM,NER,USE_REGULAR_MESH)
-
-! loop on all the mesh points in current subregion
- do ir = ir1,ir2,dir
- do iy = iy1,iy2,diy
- do ix = ix1,ix2,dix
-
-! loop over the NGNOD nodes
- do ia=1,NGNOD
- 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))
- enddo
-
-! add one spectral element to the list and store its material number
- ispec = ispec + 1
- if(ispec > nspec) call exit_MPI(myrank,'ispec greater than nspec in mesh creation')
- idoubling(ispec) = doubling_index
-
-! assign Moho surface element
- if (SAVE_MOHO_MESH) then
- if (isubregion == 15 .and. ir == ir1) then
- nspec_moho_top = nspec_moho_top + 1
- if (nspec_moho_top > NSPEC2D_BOTTOM) call exit_mpi(myrank,"Error counting moho top elements")
- ibelm_moho_top(nspec_moho_top) = ispec
- call compute_jacobian_2D(myrank,nspec_moho_top,xelm(1:NGNOD2D),yelm(1:NGNOD2D),zelm(1:NGNOD2D), &
- dershape2D_bottom,jacobian2D_moho,normal_moho,NGLLX,NGLLY,NSPEC2D_BOTTOM)
- is_moho_top(ispec) = .true.
- else if (isubregion == 28 .and. ir+dir > ir2) then
- nspec_moho_bottom = nspec_moho_bottom + 1
- if (nspec_moho_bottom > NSPEC2D_BOTTOM) call exit_mpi(myrank,"Error counting moho bottom elements")
- ibelm_moho_bot(nspec_moho_bottom) = ispec
- is_moho_bot(ispec) = .true.
- endif
- endif
-
-! initialize flag indicating whether element is in sediments
- not_fully_in_bedrock(ispec) = .false.
-
-! create mesh element
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
-
-! compute mesh coordinates
- xmesh = ZERO
- ymesh = ZERO
- zmesh = ZERO
- do ia=1,NGNOD
- xmesh = xmesh + shape3D(ia,i,j,k)*xelm(ia)
- ymesh = ymesh + shape3D(ia,i,j,k)*yelm(ia)
- zmesh = zmesh + shape3D(ia,i,j,k)*zelm(ia)
- enddo
-
- xstore_local(i,j,k) = xmesh
- ystore_local(i,j,k) = ymesh
- zstore_local(i,j,k) = zmesh
-
-! initialize flag indicating whether point is in the sediments
- point_is_in_sediments = .false.
-
- if(ANISOTROPY) then
- call aniso_model(doubling_index,zmesh,rho,vp,vs,c11,c12,c13,c14,c15,c16,&
- c22,c23,c24,c25,c26,c33,c34,c35,c36,c44,c45,c46,c55,c56,c66)
- else
-! get the regional model parameters
- if(HAUKSSON_REGIONAL_MODEL) then
-! get density from socal model
- call socal_model(doubling_index,rho,vp,vs)
-! get vp and vs from Hauksson
- call hauksson_model(vp_hauksson,vs_hauksson,xmesh,ymesh,zmesh,vp,vs,MOHO_MAP_LUPEI)
-! if Moho map is used, then assume homogeneous medium below the Moho
-! and use bottom layer of Hauksson's model in the halfspace
- if(MOHO_MAP_LUPEI .and. doubling_index == IFLAG_HALFSPACE_MOHO) &
- call socal_model(IFLAG_HALFSPACE_MOHO,rho,vp,vs)
- else
- call socal_model(doubling_index,rho,vp,vs)
-! include attenuation in first SoCal layer if needed
-! uncomment line below to include attenuation in the 1D case
-! if(zmesh >= DEPTH_5p5km_SOCAL) point_is_in_sediments = .true.
- endif
-
-! get the Harvard 3-D basin model
- if(HARVARD_3D_GOCAD_MODEL .and. &
- (doubling_index == IFLAG_ONE_LAYER_TOPOGRAPHY &
- .or. doubling_index == IFLAG_BASEMENT_TOPO) &
- .and. xmesh >= ORIG_X_GOCAD_MR &
- .and. xmesh <= END_X_GOCAD_MR &
- .and. ymesh >= ORIG_Y_GOCAD_MR &
- .and. ymesh <= END_Y_GOCAD_MR) then
-
-! use medium-resolution model first
- call interpolate_gocad_block_MR(vp_block_gocad_MR, &
- xmesh,ymesh,zmesh,rho,vp,vs,point_is_in_sediments, &
- VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM, &
- IMPOSE_MINIMUM_VP_GOCAD,THICKNESS_TAPER_BLOCK_MR, &
- vp_hauksson,vs_hauksson,doubling_index,HAUKSSON_REGIONAL_MODEL,&
- MOHO_MAP_LUPEI)
-
-! then superimpose high-resolution model
- if(xmesh >= ORIG_X_GOCAD_HR &
- .and. xmesh <= END_X_GOCAD_HR &
- .and. ymesh >= ORIG_Y_GOCAD_HR &
- .and. ymesh <= END_Y_GOCAD_HR) &
- call interpolate_gocad_block_HR(vp_block_gocad_HR,vp_block_gocad_MR,&
- xmesh,ymesh,zmesh,rho,vp,vs,point_is_in_sediments, &
- VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM, &
- IMPOSE_MINIMUM_VP_GOCAD,THICKNESS_TAPER_BLOCK_HR, &
- vp_hauksson,vs_hauksson,doubling_index,HAUKSSON_REGIONAL_MODEL, &
- MOHO_MAP_LUPEI)
-
- endif
-! get the Harvard Salton Trough model
- if (HARVARD_3D_GOCAD_MODEL) then
- call vx_xyz2uvw(xmesh, ymesh, zmesh, umesh, vmesh, wmesh)
- if (umesh >= 0 .and. umesh <= GOCAD_ST_NU-1 .and. &
- vmesh >= 0 .and. vmesh <= GOCAD_ST_NV-1 .and. &
- wmesh >= 0 .and. wmesh <= GOCAD_ST_NW-1) then
- call vx_xyz_interp(umesh,vmesh,wmesh, vp_st, vs_st, rho_st, vp_st_gocad)
- if (abs(vp_st - GOCAD_ST_NO_DATA_VALUE) > 1.0d-3) then
- vp = vp_st
- vs = vs_st
- rho = rho_st
- endif
- endif
- endif
- endif
-! store flag indicating whether point is in the sediments
- flag_sediments(i,j,k,ispec) = point_is_in_sediments
- if(point_is_in_sediments) not_fully_in_bedrock(ispec) = .true.
-
-! define elastic parameters in the model
-! distinguish between single and double precision for reals
- if(ANISOTROPY) then
-
- if(CUSTOM_REAL == SIZE_REAL) then
- rhostore(i,j,k,ispec) = sngl(rho)
- kappastore(i,j,k,ispec) = sngl(rho*(vp*vp - 4.d0*vs*vs/3.d0))
- mustore(i,j,k,ispec) = sngl(rho*vs*vs)
- c11store(i,j,k,ispec) = sngl(c11)
- c12store(i,j,k,ispec) = sngl(c12)
- c13store(i,j,k,ispec) = sngl(c13)
- c14store(i,j,k,ispec) = sngl(c14)
- c15store(i,j,k,ispec) = sngl(c15)
- c16store(i,j,k,ispec) = sngl(c16)
- c22store(i,j,k,ispec) = sngl(c22)
- c23store(i,j,k,ispec) = sngl(c23)
- c24store(i,j,k,ispec) = sngl(c24)
- c25store(i,j,k,ispec) = sngl(c25)
- c26store(i,j,k,ispec) = sngl(c26)
- c33store(i,j,k,ispec) = sngl(c33)
- c34store(i,j,k,ispec) = sngl(c34)
- c35store(i,j,k,ispec) = sngl(c35)
- c36store(i,j,k,ispec) = sngl(c36)
- c44store(i,j,k,ispec) = sngl(c44)
- c45store(i,j,k,ispec) = sngl(c45)
- c46store(i,j,k,ispec) = sngl(c46)
- c55store(i,j,k,ispec) = sngl(c55)
- c56store(i,j,k,ispec) = sngl(c56)
- c66store(i,j,k,ispec) = sngl(c66)
-! Stacey
- rho_vp(i,j,k,ispec) = sngl(rho*vp)
- rho_vs(i,j,k,ispec) = sngl(rho*vs)
- else
- rhostore(i,j,k,ispec) = rho
- kappastore(i,j,k,ispec) = rho*(vp*vp - 4.d0*vs*vs/3.d0)
- mustore(i,j,k,ispec) = rho*vs*vs
- c11store(i,j,k,ispec) = c11
- c12store(i,j,k,ispec) = c12
- c13store(i,j,k,ispec) = c13
- c14store(i,j,k,ispec) = c14
- c15store(i,j,k,ispec) = c15
- c16store(i,j,k,ispec) = c16
- c22store(i,j,k,ispec) = c22
- c23store(i,j,k,ispec) = c23
- c24store(i,j,k,ispec) = c24
- c25store(i,j,k,ispec) = c25
- c26store(i,j,k,ispec) = c26
- c33store(i,j,k,ispec) = c33
- c34store(i,j,k,ispec) = c34
- c35store(i,j,k,ispec) = c35
- c36store(i,j,k,ispec) = c36
- c44store(i,j,k,ispec) = c44
- c45store(i,j,k,ispec) = c45
- c46store(i,j,k,ispec) = c46
- c55store(i,j,k,ispec) = c55
- c56store(i,j,k,ispec) = c56
- c66store(i,j,k,ispec) = c66
-! Stacey
- rho_vp(i,j,k,ispec) = rho*vp
- rho_vs(i,j,k,ispec) = rho*vs
- endif
-
-
- else
- if(CUSTOM_REAL == SIZE_REAL) then
- rhostore(i,j,k,ispec) = sngl(rho)
- kappastore(i,j,k,ispec) = sngl(rho*(vp*vp - 4.d0*vs*vs/3.d0))
- mustore(i,j,k,ispec) = sngl(rho*vs*vs)
-
-! Stacey
- rho_vp(i,j,k,ispec) = sngl(rho*vp)
- rho_vs(i,j,k,ispec) = sngl(rho*vs)
- else
- rhostore(i,j,k,ispec) = rho
- kappastore(i,j,k,ispec) = rho*(vp*vp - 4.d0*vs*vs/3.d0)
- mustore(i,j,k,ispec) = rho*vs*vs
-
-! Stacey
- rho_vp(i,j,k,ispec) = rho*vp
- rho_vs(i,j,k,ispec) = rho*vs
- endif
- endif
-
-enddo
-enddo
-enddo
-
-! detect mesh boundaries
- call get_flags_boundaries(nspec,iproc_xi,iproc_eta,ispec,doubling_index, &
- xstore_local,ystore_local,zstore_local, &
- iboun,iMPIcut_xi,iMPIcut_eta,NPROC_XI,NPROC_ETA, &
- UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK)
-
-! compute coordinates and jacobian
- call calc_jacobian(myrank,xixstore,xiystore,xizstore, &
- etaxstore,etaystore,etazstore, &
- gammaxstore,gammaystore,gammazstore,jacobianstore, &
- xstore,ystore,zstore, &
- xelm,yelm,zelm,shape3D,dershape3D,ispec,nspec)
-
-! end of loop on all the mesh points in current subregion
- enddo
- enddo
- enddo
-
-! end of loop on all the subregions of the current region the mesh
- enddo
-
-! check total number of spectral elements created
- if(ispec /= nspec) call exit_MPI(myrank,'ispec should equal nspec')
- if (SAVE_MOHO_MESH) then
- if (nspec_moho_top /= NSPEC2D_BOTTOM .or. nspec_moho_bottom /= NSPEC2D_BOTTOM) &
- call exit_mpi(myrank, "nspec_moho should equal NSPEC2D_BOTTOM")
- endif
-
-
- do ispec=1,nspec
- ieoff = NGLLCUBE*(ispec-1)
- ilocnum = 0
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- ilocnum = ilocnum + 1
- xp(ilocnum+ieoff) = xstore(i,j,k,ispec)
- yp(ilocnum+ieoff) = ystore(i,j,k,ispec)
- zp(ilocnum+ieoff) = zstore(i,j,k,ispec)
- enddo
- enddo
- enddo
- enddo
-
- call get_global(nspec,xp,yp,zp,iglob,locval,ifseg,nglob,npointot,UTM_X_MIN,UTM_X_MAX)
-
-! put in classical format
- do ispec=1,nspec
- ieoff = NGLLCUBE*(ispec-1)
- ilocnum = 0
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- ilocnum = ilocnum + 1
- ibool(i,j,k,ispec) = iglob(ilocnum+ieoff)
- enddo
- enddo
- enddo
- enddo
-
- if(minval(ibool(:,:,:,:)) /= 1 .or. maxval(ibool(:,:,:,:)) /= NGLOB_AB) &
- call exit_MPI(myrank,'incorrect global numbering')
-
-! create a new indirect addressing array instead, to reduce cache misses
-! in memory access in the solver
- allocate(copy_ibool_ori(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(mask_ibool(nglob))
- mask_ibool(:) = -1
- copy_ibool_ori(:,:,:,:) = ibool(:,:,:,:)
-
- inumber = 0
- do ispec=1,nspec
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- if(mask_ibool(copy_ibool_ori(i,j,k,ispec)) == -1) then
-! create a new point
- inumber = inumber + 1
- ibool(i,j,k,ispec) = inumber
- mask_ibool(copy_ibool_ori(i,j,k,ispec)) = inumber
- else
-! use an existing point created previously
- ibool(i,j,k,ispec) = mask_ibool(copy_ibool_ori(i,j,k,ispec))
- endif
- enddo
- enddo
- enddo
- enddo
- deallocate(copy_ibool_ori)
- deallocate(mask_ibool)
-
-! creating mass matrix (will be fully assembled with MPI in the solver)
- allocate(rmass(nglob))
- rmass(:) = 0._CUSTOM_REAL
-
- do ispec=1,nspec
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- weight=wxgll(i)*wygll(j)*wzgll(k)
- iglobnum=ibool(i,j,k,ispec)
-
- jacobianl=jacobianstore(i,j,k,ispec)
-
-! distinguish between single and double precision for reals
- if(CUSTOM_REAL == SIZE_REAL) then
- rmass(iglobnum) = rmass(iglobnum) + &
- sngl(dble(rhostore(i,j,k,ispec)) * dble(jacobianl) * weight)
- else
- rmass(iglobnum) = rmass(iglobnum) + rhostore(i,j,k,ispec) * jacobianl * weight
- endif
-
- enddo
- enddo
- enddo
- enddo
-
- call get_jacobian_boundaries(myrank,iboun,nspec,xstore,ystore,zstore, &
- dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
- ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, &
- nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax, &
- jacobian2D_xmin,jacobian2D_xmax, &
- jacobian2D_ymin,jacobian2D_ymax, &
- jacobian2D_bottom,jacobian2D_top, &
- normal_xmin,normal_xmax, &
- normal_ymin,normal_ymax, &
- normal_bottom,normal_top, &
- NSPEC2D_BOTTOM,NSPEC2D_TOP, &
- NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX)
-
-! create MPI buffers
-! arrays locval(npointot) and ifseg(npointot) used to save memory
- call get_MPI_cutplanes_xi(myrank,prname,nspec,iMPIcut_xi,ibool, &
- xstore,ystore,zstore,ifseg,npointot, &
- NSPEC2D_A_ETA,NSPEC2D_B_ETA)
- call get_MPI_cutplanes_eta(myrank,prname,nspec,iMPIcut_eta,ibool, &
- xstore,ystore,zstore,ifseg,npointot, &
- NSPEC2D_A_XI,NSPEC2D_B_XI)
-
-! Stacey put back
- call get_absorb(myrank,prname,iboun,nspec, &
- nimin,nimax,njmin,njmax,nkmin_xi,nkmin_eta, &
- NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM)
-
-! create AVS or DX mesh data for the slice, edges and faces
- if(SAVE_MESH_FILES) then
- call write_AVS_DX_global_data(myrank,prname,nspec,ibool,idoubling,xstore,ystore,zstore,locval,ifseg,npointot)
- call write_AVS_DX_global_faces_data(myrank,prname,nspec,iMPIcut_xi,iMPIcut_eta,ibool, &
- idoubling,xstore,ystore,zstore,locval,ifseg,npointot)
- call write_AVS_DX_surface_data(myrank,prname,nspec,iboun,ibool, &
- idoubling,xstore,ystore,zstore,locval,ifseg,npointot)
- endif
-
-! create ocean load mass matrix
- if(OCEANS) then
-
-! adding ocean load mass matrix at ocean bottom
- nglob_oceans = nglob
- allocate(rmass_ocean_load(nglob_oceans))
-
-! create ocean load mass matrix for degrees of freedom at ocean bottom
- rmass_ocean_load(:) = 0._CUSTOM_REAL
-
-! add contribution of the oceans for surface elements exactly at ocean bottom
- do ispec2D_ocean_bottom = 1,NSPEC2D_TOP
-
- ispec_oceans = ibelm_top(ispec2D_ocean_bottom)
-
- iz_oceans = NGLLZ
-
- do ix_oceans = 1,NGLLX
- do iy_oceans = 1,NGLLY
-
- iglobnum=ibool(ix_oceans,iy_oceans,iz_oceans,ispec_oceans)
-
-! compute local height of oceans
-
-! get coordinates of current point
- xval = xstore(ix_oceans,iy_oceans,iz_oceans,ispec_oceans)
- yval = ystore(ix_oceans,iy_oceans,iz_oceans,ispec_oceans)
-
-! project x and y in UTM back to long/lat since topo file is in long/lat
- call utm_geo(long,lat,xval,yval,UTM_PROJECTION_ZONE,IUTM2LONGLAT,SUPPRESS_UTM_PROJECTION)
-
-! get coordinate of corner in bathy/topo 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 = &
- itopo_bathy(icornerlong,icornerlat)*(1.-ratio_xi)*(1.-ratio_eta) + &
- itopo_bathy(icornerlong+1,icornerlat)*ratio_xi*(1.-ratio_eta) + &
- itopo_bathy(icornerlong+1,icornerlat+1)*ratio_xi*ratio_eta + &
- itopo_bathy(icornerlong,icornerlat+1)*(1.-ratio_xi)*ratio_eta
-
-! suppress positive elevation, which means no oceans
- if(elevation >= - MINIMUM_THICKNESS_3D_OCEANS) then
- height_oceans = 0.d0
- else
- height_oceans = dabs(elevation)
- endif
-
-! take into account inertia of water column
- weight = wxgll(ix_oceans)*wygll(iy_oceans)*dble(jacobian2D_top(ix_oceans,iy_oceans,ispec2D_ocean_bottom)) &
- * dble(RHO_OCEANS) * height_oceans
-
-! distinguish between single and double precision for reals
- if(CUSTOM_REAL == SIZE_REAL) then
- rmass_ocean_load(iglobnum) = rmass_ocean_load(iglobnum) + sngl(weight)
- else
- rmass_ocean_load(iglobnum) = rmass_ocean_load(iglobnum) + weight
- endif
-
- enddo
- enddo
-
- enddo
-
-! add regular mass matrix to ocean load contribution
- rmass_ocean_load(:) = rmass_ocean_load(:) + rmass(:)
-
- else
-
-! allocate dummy array if no oceans
- nglob_oceans = 1
- allocate(rmass_ocean_load(nglob_oceans))
-
- endif
-
-! save the binary files
- call save_arrays_solver(flag_sediments,not_fully_in_bedrock,rho_vp,rho_vs,prname,xixstore,xiystore,xizstore, &
- etaxstore,etaystore,etazstore, &
- gammaxstore,gammaystore,gammazstore,jacobianstore, &
- xstore,ystore,zstore,kappastore,mustore, &
- ANISOTROPY, &
- c11store,c12store,c13store,c14store,c15store,c16store, &
- c22store,c23store,c24store,c25store,c26store,c33store,c34store,c35store,c36store, &
- c44store,c45store,c46store,c55store,c56store,c66store, &
- ibool,idoubling,rmass,rmass_ocean_load,nglob_oceans, &
- ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, &
- nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax, &
- normal_xmin,normal_xmax,normal_ymin,normal_ymax,normal_bottom,normal_top, &
- jacobian2D_xmin,jacobian2D_xmax, &
- jacobian2D_ymin,jacobian2D_ymax, &
- jacobian2D_bottom,jacobian2D_top, &
- iMPIcut_xi,iMPIcut_eta,nspec,nglob, &
- NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP,OCEANS)
-
-! save Moho mesh arrays
- if (SAVE_MOHO_MESH) then
- open(unit=27,file=prname(1:len_trim(prname))//'ibelm_moho.bin',status='unknown',form='unformatted')
- ! total number of elements, corner points, all points
- write(27) NSPEC2D_BOTTOM
- write(27) (NEX_PER_PROC_XI/4 + 1) * (NEX_PER_PROC_ETA/4 + 1)
- write(27) (NEX_PER_PROC_XI/4 * (NGLLX - 1) + 1) * (NEX_PER_PROC_ETA/4 * (NGLLY - 1) + 1)
- write(27) ibelm_moho_top
- write(27) ibelm_moho_bot
- close(27)
- open(unit=27,file=prname(1:len_trim(prname))//'normal_moho.bin',status='unknown',form='unformatted')
- write(27) normal_moho
- close(27)
- open(unit=27,file=prname(1:len_trim(prname))//'is_moho.bin',status='unknown',form='unformatted')
- write(27) is_moho_top
- write(27) is_moho_bot
- close(27)
- endif
-
-
- do ispec=1,nspec
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- weight=wxgll(i)*wygll(j)*wzgll(k)
- jacobianl=jacobianstore(i,j,k,ispec)
- volume_local = volume_local + dble(jacobianl)*weight
- enddo
- enddo
- enddo
- enddo
-
- do ispec = 1,NSPEC2D_BOTTOM
- do i=1,NGLLX
- do j=1,NGLLY
- weight=wxgll(i)*wygll(j)
- area_local_bottom = area_local_bottom + dble(jacobian2D_bottom(i,j,ispec))*weight
- enddo
- enddo
- enddo
-
- do ispec = 1,NSPEC2D_TOP
- do i=1,NGLLX
- do j=1,NGLLY
- weight=wxgll(i)*wygll(j)
- area_local_top = area_local_top + dble(jacobian2D_top(i,j,ispec))*weight
- enddo
- enddo
- enddo
-
- end subroutine create_regions_mesh
-
-!
-!----
-!
-
subroutine create_regions_mesh_ext_mesh(ibool, &
xstore,ystore,zstore,nspec,npointot,myrank,LOCAL_PATH, &
nnodes_ext_mesh,nelmnts_ext_mesh, &
- nodes_coords_ext_mesh,elmnts_ext_mesh,mat_ext_mesh,max_static_memory_size, &
- ninterface_ext_mesh,max_interface_size_ext_mesh, &
+ nodes_coords_ext_mesh,elmnts_ext_mesh,max_static_memory_size,mat_ext_mesh,materials_ext_mesh, &
+ nmat_ext_mesh,undef_mat_prop,nundefMat_ext_mesh,ninterface_ext_mesh,max_interface_size_ext_mesh, &
my_neighbours_ext_mesh,my_nelmnts_neighbours_ext_mesh,my_interfaces_ext_mesh, &
- ibool_interfaces_ext_mesh,nibool_interfaces_ext_mesh)
+ ibool_interfaces_ext_mesh,nibool_interfaces_ext_mesh, &
+ nspec2D_xmin, nspec2D_xmax, nspec2D_ymin, nspec2D_ymax, NSPEC2D_BOTTOM, NSPEC2D_TOP,&
+ NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX, &
+ ibelm_xmin,ibelm_xmax, ibelm_ymin, ibelm_ymax, ibelm_bottom, ibelm_top)
! create the different regions of the mesh
@@ -1187,14 +60,18 @@
integer :: nnodes_ext_mesh,nelmnts_ext_mesh
double precision, dimension(NDIM,nnodes_ext_mesh) :: nodes_coords_ext_mesh
integer, dimension(ESIZE,nelmnts_ext_mesh) :: elmnts_ext_mesh
- integer, dimension(nelmnts_ext_mesh) :: mat_ext_mesh
- double precision, external :: materials_ext_mesh
+ integer, dimension(2,nelmnts_ext_mesh) :: mat_ext_mesh
+! double precision, external :: materials_ext_mesh
integer :: ninterface_ext_mesh,max_interface_size_ext_mesh
integer, dimension(ninterface_ext_mesh) :: my_neighbours_ext_mesh
integer, dimension(ninterface_ext_mesh) :: my_nelmnts_neighbours_ext_mesh
integer, dimension(6,max_interface_size_ext_mesh,ninterface_ext_mesh) :: my_interfaces_ext_mesh
integer, dimension(NGLLX*NGLLX*max_interface_size_ext_mesh,ninterface_ext_mesh) :: ibool_interfaces_ext_mesh
integer, dimension(ninterface_ext_mesh) :: nibool_interfaces_ext_mesh
+!pll
+ integer :: nmat_ext_mesh,nundefMat_ext_mesh
+ double precision, dimension(5,nmat_ext_mesh) :: materials_ext_mesh
+ character (len=30), dimension(5,nundefMat_ext_mesh):: undef_mat_prop
! for MPI buffers
integer, dimension(:), allocatable :: reorder_interface_ext_mesh,ind_ext_mesh,ninseg_ext_mesh,iwork_ext_mesh
@@ -1225,8 +102,11 @@
etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore,jacobianstore
! for model density
- real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: rhostore,kappastore,mustore
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: rhostore,kappastore,mustore,vpstore,vsstore
+! attenuation
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: iflag_attenuation_store
+
! proc numbers for MPI
integer myrank
@@ -1258,7 +138,75 @@
integer, dimension(:,:,:,:), allocatable :: copy_ibool_ori
integer :: inumber
+! pll
+ integer :: iundef
+ logical, dimension(6,nspec) :: iboun
+! 2D shape functions and their derivatives
+ double precision, dimension(:,:,:), allocatable :: shape2D_x,shape2D_y,shape2D_bottom,shape2D_top
+ double precision, dimension(:,:,:,:), allocatable :: dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top
+
+ integer :: nspec2D_xmin, nspec2D_xmax, nspec2D_ymin, nspec2D_ymax, NSPEC2D_BOTTOM, NSPEC2D_TOP
+ integer :: NSPEC2DMAX_XMIN_XMAX
+ integer :: NSPEC2DMAX_YMIN_YMAX
+ 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
+ integer, dimension(2,NSPEC2DMAX_YMIN_YMAX) :: nimin,nimax
+ integer, dimension(2,NSPEC2DMAX_XMIN_XMAX) :: njmin,njmax
+ integer, dimension(2,NSPEC2DMAX_XMIN_XMAX) :: nkmin_xi
+ integer, dimension(2,NSPEC2DMAX_YMIN_YMAX) :: nkmin_eta
+
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: ibedrock
+
+ integer :: ispec2D
+ integer :: iflag, flag_below, flag_above
+
+ ! 2-D jacobians and normals
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: &
+ jacobian2D_xmin,jacobian2D_xmax, &
+ jacobian2D_ymin,jacobian2D_ymax,jacobian2D_bottom,jacobian2D_top
+
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: &
+ normal_xmin,normal_xmax,normal_ymin,normal_ymax,normal_bottom,normal_top
+
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: rho_vp,rho_vs
+
+
+
+! For Piero Basini :
+! integer :: doubling_value_found_for_Piero
+! double precision :: xmesh,ymesh,zmesh
+! double precision :: rho,vp,vs
+
+! integer,dimension(nspec) :: idoubling
+! integer :: doubling_value_found_for_Piero
+! integer, parameter :: NUMBER_OF_STATIONS = 6
+! double precision, parameter :: RADIUS_TO_EXCLUDE = 250.d0
+! double precision, dimension(NUMBER_OF_STATIONS) :: utm_x_station,utm_y_station
+
+! 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
+
+! ! size of topography and bathymetry file for Piero Basini's model
+! integer, parameter :: NX_TOPO = 787, NY_TOPO = 793
+! double precision, parameter :: ORIG_LAT_TOPO = -102352.d0
+! double precision, parameter :: ORIG_LONG_TOPO = 729806.d0
+! character(len=150), parameter :: TOPO_FILE = 'DATA/piero_model/dem_EV_UTM_regular_250_reordered.dat'
+! ! for Piero Basini's model this is the resolution in meters of the topo file
+! double precision, parameter :: DEGREES_PER_CELL_TOPO = 250.d0
+
+
+
+
! **************
! create the name for the database of the current slide and region
@@ -1278,10 +226,26 @@
allocate(shape3D(NGNOD,NGLLX,NGLLY,NGLLZ))
allocate(dershape3D(NDIM,NGNOD,NGLLX,NGLLY,NGLLZ))
+! pll 2D shape functions and their derivatives
+ allocate(shape2D_x(NGNOD2D,NGLLY,NGLLZ))
+ allocate(shape2D_y(NGNOD2D,NGLLX,NGLLZ))
+ allocate(shape2D_bottom(NGNOD2D,NGLLX,NGLLY))
+ allocate(shape2D_top(NGNOD2D,NGLLX,NGLLY))
+ allocate(dershape2D_x(NDIM2D,NGNOD2D,NGLLY,NGLLZ))
+ allocate(dershape2D_y(NDIM2D,NGNOD2D,NGLLX,NGLLZ))
+ allocate(dershape2D_bottom(NDIM2D,NGNOD2D,NGLLX,NGLLY))
+ allocate(dershape2D_top(NDIM2D,NGNOD2D,NGLLX,NGLLY))
+
+! pll Stacey
+ allocate(rho_vp(NGLLX,NGLLY,NGLLZ,nspec))
+ allocate(rho_vs(NGLLX,NGLLY,NGLLZ,nspec))
+
! array with model density
allocate(rhostore(NGLLX,NGLLY,NGLLZ,nspec))
allocate(kappastore(NGLLX,NGLLY,NGLLZ,nspec))
allocate(mustore(NGLLX,NGLLY,NGLLZ,nspec))
+ allocate(vpstore(NGLLX,NGLLY,NGLLZ,nspec)) !pll
+ allocate(vsstore(NGLLX,NGLLY,NGLLZ,nspec)) !pll
! arrays with mesh parameters
allocate(xixstore(NGLLX,NGLLY,NGLLZ,nspec))
@@ -1295,6 +259,21 @@
allocate(gammazstore(NGLLX,NGLLY,NGLLZ,nspec))
allocate(jacobianstore(NGLLX,NGLLY,NGLLZ,nspec))
+! pll 2-D jacobians and normals
+ allocate(jacobian2D_xmin(NGLLY,NGLLZ,nspec2D_xmin))
+ allocate(jacobian2D_xmax(NGLLY,NGLLZ,nspec2D_xmax))
+ allocate(jacobian2D_ymin(NGLLX,NGLLZ,nspec2D_ymin))
+ allocate(jacobian2D_ymax(NGLLX,NGLLZ,nspec2D_ymax))
+ allocate(jacobian2D_bottom(NGLLX,NGLLY,NSPEC2D_BOTTOM))
+ allocate(jacobian2D_top(NGLLX,NGLLY,NSPEC2D_TOP))
+
+ allocate(normal_xmin(NDIM,NGLLY,NGLLZ,nspec2D_xmin))
+ allocate(normal_xmax(NDIM,NGLLY,NGLLZ,nspec2D_xmax))
+ allocate(normal_ymin(NDIM,NGLLX,NGLLZ,nspec2D_ymin))
+ allocate(normal_ymax(NDIM,NGLLX,NGLLZ,nspec2D_ymax))
+ allocate(normal_bottom(NDIM,NGLLX,NGLLY,NSPEC2D_BOTTOM))
+ allocate(normal_top(NDIM,NGLLX,NGLLY,NSPEC2D_TOP))
+
! set up coordinates of the Gauss-Lobatto-Legendre points
call zwgljd(xigll,wxgll,NGLLX,GAUSSALPHA,GAUSSBETA)
call zwgljd(yigll,wygll,NGLLY,GAUSSALPHA,GAUSSBETA)
@@ -1308,6 +287,12 @@
! get the 3-D shape functions
call get_shape3D(myrank,shape3D,dershape3D,xigll,yigll,zigll)
+! pll get the 2-D shape functions
+ call get_shape2D(myrank,shape2D_x,dershape2D_x,yigll,zigll,NGLLY,NGLLZ)
+ call get_shape2D(myrank,shape2D_y,dershape2D_y,xigll,zigll,NGLLX,NGLLZ)
+ call get_shape2D(myrank,shape2D_bottom,dershape2D_bottom,xigll,yigll,NGLLX,NGLLY)
+ call get_shape2D(myrank,shape2D_top,dershape2D_top,xigll,yigll,NGLLX,NGLLY)
+
! allocate memory for arrays
allocate(iglob(npointot))
allocate(locval(npointot))
@@ -1338,21 +323,239 @@
enddo
+! ! 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)
+
! kappastore and mustore
do ispec = 1, nspec
do k = 1, NGLLZ
do j = 1, NGLLY
do i = 1, NGLLX
- 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))*&
- materials_ext_mesh(3,mat_ext_mesh(ispec))
+! check if the material is known or unknown
+ if (mat_ext_mesh(1,ispec) > 0) then
+ rhostore(i,j,k,ispec) = materials_ext_mesh(1,mat_ext_mesh(1,ispec))
+ vpstore(i,j,k,ispec) = materials_ext_mesh(2,mat_ext_mesh(1,ispec))
+ vsstore(i,j,k,ispec) = materials_ext_mesh(3,mat_ext_mesh(1,ispec))
+ iflag_attenuation_store(i,j,k,ispec) = materials_ext_mesh(4,mat_ext_mesh(1,ispec))
+ !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
+ else if (mat_ext_mesh(2,ispec) == 1) then
+ 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
+ !call interface(iflag,flag_below,flag_above,ispec,nspec,i,j,k,xstore,ystore,zstore,ibedrock)
+ rhostore(i,j,k,ispec) = materials_ext_mesh(1,iflag)
+ vpstore(i,j,k,ispec) = materials_ext_mesh(2,iflag)
+ vsstore(i,j,k,ispec) = materials_ext_mesh(3,iflag)
+ iflag_attenuation_store(i,j,k,ispec) = materials_ext_mesh(4,iflag)
+ !change for piero :
+ ! if(iflag == 1) then
+ ! iflag_attenuation_store(i,j,k,ispec) = 1
+ ! else
+ ! iflag_attenuation_store(i,j,k,ispec) = 2
+ ! endif
+ else
+ ! call tomography()
+ 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
+
enddo
enddo
enddo
enddo
+
+! !! 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
+
locval = 0
ifseg = .false.
xp = 0.d0
@@ -1466,16 +669,50 @@
! distinguish between single and double precision for reals
if(CUSTOM_REAL == SIZE_REAL) then
rmass(iglobnum) = rmass(iglobnum) + &
- sngl(dble(materials_ext_mesh(1,mat_ext_mesh(ispec))) * dble(jacobianl) * weight)
+ sngl((dble(rhostore(i,j,k,ispec))) * dble(jacobianl) * weight)
else
- rmass(iglobnum) = rmass(iglobnum) + materials_ext_mesh(1,mat_ext_mesh(ispec)) * jacobianl * weight
+ rmass(iglobnum) = rmass(iglobnum) + rhostore(i,j,k,ispec) * jacobianl * weight
endif
enddo
enddo
enddo
- enddo
+ enddo
+
+ iboun(:,:) = .false.
+ do ispec2D = 1, nspec2D_xmin
+ iboun(1,ibelm_xmin(ispec2D)) = .true.
+ end do
+ do ispec2D = 1, nspec2D_xmax
+ iboun(2,ibelm_xmax(ispec2D)) = .true.
+ end do
+ do ispec2D = 1, nspec2D_ymin
+ iboun(3,ibelm_ymin(ispec2D)) = .true.
+ end do
+ do ispec2D = 1, nspec2D_ymax
+ iboun(4,ibelm_ymax(ispec2D)) = .true.
+ end do
+ do ispec2D = 1, NSPEC2D_BOTTOM
+ iboun(5,ibelm_bottom(ispec2D)) = .true.
+ end do
+ do ispec2D = 1, NSPEC2D_TOP
+ iboun(6,ibelm_top(ispec2D)) = .true.
+ end do
+
+ call get_jacobian_boundaries(myrank,iboun,nspec,xstore,ystore,zstore, &
+ dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, &
+ ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, &
+ nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax, &
+ jacobian2D_xmin,jacobian2D_xmax, &
+ jacobian2D_ymin,jacobian2D_ymax, &
+ jacobian2D_bottom,jacobian2D_top, &
+ normal_xmin,normal_xmax, &
+ normal_ymin,normal_ymax, &
+ normal_bottom,normal_top, &
+ NSPEC2D_BOTTOM,NSPEC2D_TOP, &
+ NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX)
+
call prepare_assemble_MPI (nelmnts_ext_mesh,ibool, &
elmnts_ext_mesh, ESIZE, &
nglob, &
@@ -1485,6 +722,11 @@
nibool_interfaces_ext_mesh &
)
+ ! Stacey put back
+ call get_absorb_ext_mesh(myrank,iboun,nspec, &
+ nimin,nimax,njmin,njmax,nkmin_xi,nkmin_eta, &
+ NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM)
+
! sort ibool comm buffers lexicographically
allocate(nibool_interfaces_ext_mesh_true(ninterface_ext_mesh))
@@ -1544,6 +786,21 @@
write(IOUT) gammazstore
write(IOUT) jacobianstore
+
+ !pll Stacey
+ write(IOUT) rho_vp
+ write(IOUT) rho_vs
+ write(IOUT) iflag_attenuation_store
+ write(IOUT) NSPEC2DMAX_XMIN_XMAX
+ write(IOUT) NSPEC2DMAX_YMIN_YMAX
+ write(IOUT) nimin
+ write(IOUT) nimax
+ write(IOUT) njmin
+ write(IOUT) njmax
+ write(IOUT) nkmin_xi
+ write(IOUT) nkmin_eta
+ !end pll
+
write(IOUT) kappastore
write(IOUT) mustore
@@ -1555,6 +812,37 @@
write(IOUT) ystore_dummy
write(IOUT) zstore_dummy
+! boundary parameters
+ write(IOUT) nspec2D_xmin
+ write(IOUT) nspec2D_xmax
+ write(IOUT) nspec2D_ymin
+ write(IOUT) nspec2D_ymax
+ write(IOUT) NSPEC2D_BOTTOM
+ write(IOUT) NSPEC2D_TOP
+
+ write(IOUT) ibelm_xmin
+ write(IOUT) ibelm_xmax
+ write(IOUT) ibelm_ymin
+ write(IOUT) ibelm_ymax
+ write(IOUT) ibelm_bottom
+ write(IOUT) ibelm_top
+
+ write(IOUT) normal_xmin
+ write(IOUT) normal_xmax
+ write(IOUT) normal_ymin
+ write(IOUT) normal_ymax
+ write(IOUT) normal_bottom
+ write(IOUT) normal_top
+
+ write(IOUT) jacobian2D_xmin
+ write(IOUT) jacobian2D_xmax
+ write(IOUT) jacobian2D_ymin
+ write(IOUT) jacobian2D_ymax
+ write(IOUT) jacobian2D_bottom
+ write(IOUT) jacobian2D_top
+
+! end boundary parameters
+
write(IOUT) ninterface_ext_mesh
write(IOUT) maxval(nibool_interfaces_ext_mesh)
write(IOUT) my_neighbours_ext_mesh
@@ -1572,45 +860,7 @@
end subroutine create_regions_mesh_ext_mesh
-!
-!----
-!
- double precision function materials_ext_mesh(i,j)
-
- implicit none
-
- integer :: i,j
-
- select case (j)
- case (1)
- select case (i)
- case (1)
- materials_ext_mesh = 2700.d0
- case (2)
- materials_ext_mesh = 3000.d0
- case (3)
- materials_ext_mesh = 1732.051d0
- case default
- call stop_all()
- end select
- case (2)
- select case (i)
- case (1)
- materials_ext_mesh = 2000.d0
- case (2)
- materials_ext_mesh = 900.d0
- case (3)
- materials_ext_mesh = 500.d0
- case default
- call stop_all()
- end select
- case default
- call stop_all()
- end select
-
- end function materials_ext_mesh
-
!
!----
!
@@ -2107,3 +1357,114 @@
end if
end subroutine get_edge
+
+
+
+!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
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/compile_all.csh
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/compile_all.csh 2009-07-10 16:39:24 UTC (rev 15452)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/compile_all.csh 2009-07-10 21:02:53 UTC (rev 15453)
@@ -5,9 +5,11 @@
rm *.o *.mod ./a.out
-gfortran -c part_decompose_mesh_SCOTCH.f90
-gfortran -c decompose_mesh_SCOTCH.f90
+gfortran -Wall -c part_decompose_mesh_SCOTCH.f90
+#gfortran -c decompose_mesh_SCOTCH.f90
+gfortran -Wall -c decompose_mesh_SCOTCH.f90
#gfortran decompose_mesh_SCOTCH.o part_decompose_mesh_SCOTCH.o ~/utils/metis-4.0/libmetis.a
#gfortran decompose_mesh_SCOTCH.o part_decompose_mesh_SCOTCH.o ~/utils/scotch_5.1/lib/libscotchmetis.a ~/utils/scotch_5.1/lib/libscotch.a ~/utils/scotch_5.1/lib/libscotcherr.a
-gfortran decompose_mesh_SCOTCH.o part_decompose_mesh_SCOTCH.o ~/utils/scotch_5.1/lib/libscotch.a ~/utils/scotch_5.1/lib/libscotcherr.a
+#gfortran decompose_mesh_SCOTCH.o part_decompose_mesh_SCOTCH.o ~/utils/scotch_5.1/lib/libscotch.a ~/utils/scotch_5.1/lib/libscotcherr.a
+gfortran -Wall decompose_mesh_SCOTCH.o part_decompose_mesh_SCOTCH.o ~/Download/scotch_5.1/lib/libscotch.a ~/Download/scotch_5.1/lib/libscotcherr.a
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/constants_decompose_mesh_SCOTCH.h
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/constants_decompose_mesh_SCOTCH.h 2009-07-10 16:39:24 UTC (rev 15452)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/constants_decompose_mesh_SCOTCH.h 2009-07-10 21:02:53 UTC (rev 15453)
@@ -1,3 +1,6 @@
+! Number of slices for mesh partitioning
+integer, parameter :: nparts = 2
+
! Useful kind types
integer ,parameter :: short = SELECTED_INT_KIND(4), long = SELECTED_INT_KIND(18)
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 2009-07-10 16:39:24 UTC (rev 15452)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/decompose_mesh_SCOTCH.f90 2009-07-10 21:02:53 UTC (rev 15453)
@@ -1,30 +1,31 @@
+program pre_meshfem3D
- program decompose_mesh_SCOTCH
-
use part_decompose_mesh_SCOTCH
implicit none
-
+
include './constants_decompose_mesh_SCOTCH.h'
- include "./scotchf.h"
+ include './scotchf.h'
- integer, parameter :: nparts = 8
-
integer(long) :: nspec
integer, dimension(:,:), allocatable :: elmnts
- integer, dimension(:), allocatable :: mat
+ integer, dimension(:,:), allocatable :: mat
integer, dimension(:), allocatable :: part
-
+
integer :: nnodes
double precision, dimension(:,:), allocatable :: nodes_coords
-
+
integer, dimension(:), allocatable :: xadj
integer, dimension(:), allocatable :: adjncy
integer, dimension(:), allocatable :: nnodes_elmnts
integer, dimension(:), allocatable :: nodes_elmnts
- integer, dimension(:), allocatable :: vwgt
- integer, dimension(:), allocatable :: adjwgt
- integer, dimension(5) :: metis_options
+! Old Metis variables
+! integer :: wgtflag
+! integer :: num_start
+! integer :: edgecut
+! integer, dimension(:), allocatable :: vwgt
+! integer, dimension(:), allocatable :: adjwgt
+! integer, dimension(5) :: metis_options
integer, dimension(:), pointer :: glob2loc_elmnts
integer, dimension(:), pointer :: glob2loc_nodes_nparts
@@ -36,63 +37,171 @@
integer, dimension(:), allocatable :: my_nb_interfaces
integer :: ninterfaces
integer :: my_ninterface
-
- integer :: nb_materials
- double precision, dimension(:), allocatable :: cs
- integer, dimension(:), allocatable :: num_material
-
- integer(long) :: nsize ! Max number of elements that contain the same node.
- integer :: edgecut
+
+ integer(long) :: nsize ! Max number of elements that contain the same node.
integer :: nb_edges
integer :: ispec, inode
- integer :: wgtflag
- integer :: num_start
integer :: ngnod
- integer :: max_neighbour ! Real maximum number of neighbours per element
+ integer :: max_neighbour ! Real maximum number of neighbours per element
integer(long) :: sup_neighbour ! Majoration of the maximum number of neighbours per element
integer :: ipart, nnodes_loc, nspec_loc
+ integer :: num_elmnt, num_node, num_mat
+
+ ! boundaries
+ integer :: ispec2D
+ integer :: nspec2D_xmin, nspec2D_xmax, nspec2D_ymin, nspec2D_ymax, nspec2D_bottom, nspec2D_top
+ integer, dimension(:), allocatable :: ibelm_xmin, ibelm_xmax, ibelm_ymin, ibelm_ymax, ibelm_bottom, ibelm_top
+ integer, dimension(:,:), allocatable :: nodes_ibelm_xmin, nodes_ibelm_xmax, nodes_ibelm_ymin
+ integer, dimension(:,:), allocatable :: nodes_ibelm_ymax, nodes_ibelm_bottom, nodes_ibelm_top
+
character(len=256) :: prname
logical, dimension(:), allocatable :: mask_nodes_elmnts
integer, dimension(:), allocatable :: used_nodes_elmnts
-!!!! NL NL for SCOTCH partitioner
- double precision, dimension(SCOTCH_GRAPHDIM) :: scotchgraph
- double precision, dimension(SCOTCH_STRATDIM) :: scotchstrat
- character(len=256), parameter :: scotch_strategy='b{job=t,map=t,poli=S,sep=h{pass=30}}'
- integer :: ierr
-!!!! NL NL
+ double precision, dimension(SCOTCH_GRAPHDIM) :: scotchgraph
+ double precision, dimension(SCOTCH_STRATDIM) :: scotchstrat
+ character(len=256), parameter :: scotch_strategy='b{job=t,map=t,poli=S,sep=h{pass=30}}'
+ integer :: ierr
+
+ !pll
+ double precision , dimension(:,:), allocatable :: mat_prop
+ integer :: count_def_mat,count_undef_mat,imat
+ character (len=30), dimension(:,:), allocatable :: undef_mat_prop
ngnod = esize
-! read the elements, material and nodes files
- open(unit=98, file='../model_asteroid_subdivide/mesh', status='old', form='formatted')
+ open(unit=98, file='./OUTPUT_FILES/mesh_file', status='old', form='formatted')
read(98,*) nspec
allocate(elmnts(esize,nspec))
do ispec = 1, nspec
- read(98,*) elmnts(1,ispec), elmnts(2,ispec), elmnts(3,ispec), elmnts(4,ispec), &
- elmnts(5,ispec), elmnts(6,ispec), elmnts(7,ispec), elmnts(8,ispec)
+ 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."
end do
close(98)
+ print*, 'nspec = ', nspec
- open(unit=98, file='../model_asteroid_subdivide/mat', status='old', form='formatted')
- allocate(mat(nspec))
+ open(unit=98, file='./OUTPUT_FILES/materials_file', status='old', form='formatted')
+ allocate(mat(2,nspec))
do ispec = 1, nspec
- read(98,*) mat(ispec)
+ read(98,*) num_mat, mat(1,ispec)!, mat(2,ispec)
+ if((num_mat > nspec) .or. (num_mat < 1) ) stop "ERROR : Invalid mat file."
end do
close(98)
-
- open(unit=98, file='../model_asteroid_subdivide/nodes_coords', status='old', form='formatted')
+!must be changed, if mat(1,i) < 0 1 == interface , 2 == tomography
+ mat(2,:) = 1
+
+ open(unit=98, file='./OUTPUT_FILES/nodes_coords_file', status='old', form='formatted')
read(98,*) nnodes
allocate(nodes_coords(3,nnodes))
do inode = 1, nnodes
- read(98,*) nodes_coords(1,inode), nodes_coords(2,inode), nodes_coords(3,inode)
+ 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*, 'nnodes = ', nnodes
+ count_def_mat = 0
+ count_undef_mat = 0
+ open(unit=98, file='./OUTPUT_FILES/nummaterial_velocity_file', status='old', form='formatted')
+ read(98,*,iostat=ierr) num_mat
+ do while (ierr == 0)
+ print*, 'num_mat = ',num_mat
+ if(num_mat /= -1) then
+ count_def_mat = count_def_mat + 1
+ else
+ count_undef_mat = count_undef_mat + 1
+ end if
+ read(98,*,iostat=ierr) num_mat
+ end do
+ close(98)
+
+ print*, count_def_mat, count_undef_mat
+ allocate(mat_prop(5,count_def_mat))
+ allocate(undef_mat_prop(5,count_undef_mat))
+
+ open(unit=98, file='./OUTPUT_FILES/nummaterial_velocity_file', status='old', form='formatted')
+ do imat=1,count_def_mat
+ 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)
+ if(num_mat < 0 .or. num_mat > count_def_mat) stop "ERROR : Invalid nummaterial_velocity_file file."
+ end do
+ do imat=1,count_undef_mat
+ read(98,'(5A30)') 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)
+ end do
+ close(98)
+
+ !boundary files
+ open(unit=98, file='./OUTPUT_FILES/absorbing_surface_file_xmin', status='old', form='formatted')
+ read(98,*) nspec2D_xmin
+ allocate(ibelm_xmin(nspec2D_xmin))
+ allocate(nodes_ibelm_xmin(4,nspec2D_xmin))
+ do ispec2D = 1,nspec2D_xmin
+ 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)
+ end do
+ close(98)
+ print*, 'nspec2D_xmin = ', nspec2D_xmin
+
+ open(unit=98, file='./OUTPUT_FILES/absorbing_surface_file_xmax', status='old', form='formatted')
+ read(98,*) nspec2D_xmax
+ allocate(ibelm_xmax(nspec2D_xmax))
+ allocate(nodes_ibelm_xmax(4,nspec2D_xmax))
+ do ispec2D = 1,nspec2D_xmax
+ read(98,*) ibelm_xmax(ispec2D), nodes_ibelm_xmax(1,ispec2D), nodes_ibelm_xmax(2,ispec2D), &
+ nodes_ibelm_xmax(3,ispec2D), nodes_ibelm_xmax(4,ispec2D)
+ end do
+ close(98)
+ print*, 'nspec2D_xmax = ', nspec2D_xmax
+
+ open(unit=98, file='./OUTPUT_FILES/absorbing_surface_file_ymin', status='old', form='formatted')
+ read(98,*) nspec2D_ymin
+ allocate(ibelm_ymin(nspec2D_ymin))
+ allocate(nodes_ibelm_ymin(4,nspec2D_ymin))
+ do ispec2D = 1,nspec2D_ymin
+ read(98,*) ibelm_ymin(ispec2D), nodes_ibelm_ymin(1,ispec2D), nodes_ibelm_ymin(2,ispec2D), &
+ nodes_ibelm_ymin(3,ispec2D), nodes_ibelm_ymin(4,ispec2D)
+ end do
+ close(98)
+ print*, 'nspec2D_ymin = ', nspec2D_ymin
+
+ open(unit=98, file='./OUTPUT_FILES/absorbing_surface_file_ymax', status='old', form='formatted')
+ read(98,*) nspec2D_ymax
+ allocate(ibelm_ymax(nspec2D_ymax))
+ allocate(nodes_ibelm_ymax(4,nspec2D_ymax))
+ do ispec2D = 1,nspec2D_ymax
+ read(98,*) ibelm_ymax(ispec2D), nodes_ibelm_ymax(1,ispec2D), nodes_ibelm_ymax(2,ispec2D), &
+ nodes_ibelm_ymax(3,ispec2D), nodes_ibelm_ymax(4,ispec2D)
+ end do
+ close(98)
+ print*, 'nspec2D_ymax = ', nspec2D_ymax
+
+ open(unit=98, file='./OUTPUT_FILES/absorbing_surface_file_bottom', status='old', form='formatted')
+ read(98,*) nspec2D_bottom
+ allocate(ibelm_bottom(nspec2D_bottom))
+ allocate(nodes_ibelm_bottom(4,nspec2D_bottom))
+ do ispec2D = 1,nspec2D_bottom
+ read(98,*) ibelm_bottom(ispec2D), nodes_ibelm_bottom(1,ispec2D), nodes_ibelm_bottom(2,ispec2D), &
+ nodes_ibelm_bottom(3,ispec2D), nodes_ibelm_bottom(4,ispec2D)
+ end do
+ close(98)
+ print*, 'nspec2D_bottom = ', nspec2D_bottom
+
+ open(unit=98, file='./OUTPUT_FILES/free_surface_file', status='old', form='formatted')
+ read(98,*) nspec2D_top
+ allocate(ibelm_top(nspec2D_top))
+ allocate(nodes_ibelm_top(4,nspec2D_top))
+ do ispec2D = 1,nspec2D_top
+ read(98,*) ibelm_top(ispec2D), nodes_ibelm_top(1,ispec2D), nodes_ibelm_top(2,ispec2D), &
+ nodes_ibelm_top(3,ispec2D), nodes_ibelm_top(4,ispec2D)
+ end do
+ close(98)
+ print*, 'nspec2D_top = ', nspec2D_top
+
allocate(mask_nodes_elmnts(nnodes))
allocate(used_nodes_elmnts(nnodes))
mask_nodes_elmnts(:) = .false.
@@ -103,18 +212,15 @@
used_nodes_elmnts(elmnts(inode,ispec)) = used_nodes_elmnts(elmnts(inode,ispec)) + 1
enddo
enddo
- nsize = maxval(used_nodes_elmnts(:))
- sup_neighbour = ngnod * nsize - (ngnod + (ngnod/2 - 1)*nfaces)
- print*, 'nsize = ',nsize, 'sup_neighbour = ', sup_neighbour
-
+ print *, minval(used_nodes_elmnts(:)), maxval(used_nodes_elmnts(:))
do inode = 1, nnodes
if (.not. mask_nodes_elmnts(inode)) then
stop 'ERROR : nodes not used.'
endif
enddo
-! if (maxval(used_nodes_elmnts(:))>nsize) then
-! stop 'ERROR : increase nsize or modify the mesh.'
-! endif
+ nsize = maxval(used_nodes_elmnts(:))
+ sup_neighbour = ngnod * nsize - (ngnod + (ngnod/2 - 1)*nfaces)
+ print*, 'nsize = ',nsize, 'sup_neighbour = ', sup_neighbour
elmnts(:,:) = elmnts(:,:) - 1
@@ -122,16 +228,15 @@
allocate(adjncy(1:sup_neighbour*nspec))
allocate(nnodes_elmnts(1:nnodes))
allocate(nodes_elmnts(1:nsize*nnodes))
-
+
call mesh2dual_ncommonnodes(nspec, nnodes, nsize, sup_neighbour, elmnts, xadj, adjncy, nnodes_elmnts, &
nodes_elmnts, max_neighbour, 1)
print*, 'max_neighbour = ',max_neighbour
-
! elmnts(:,:) = elmnts(:,:) + 1
! adjncy(:) = adjncy(:) + 1
! xadj(:) = xadj(:) + 1
-! allocate(vwgt(0:nspec-1))
+ ! allocate(vwgt(0:nspec-1))
nb_edges = xadj(nspec+1)
! allocate(adjwgt(0:nb_edges-1))
! vwgt(:) = 1
@@ -142,14 +247,13 @@
! metis_options(3) = 1
! metis_options(4) = 1
! metis_options(5) = 0
-
-
+
+
! num_start = 0
! wgtflag = 0
allocate(part(1:nspec))
-
! Old metis partitioning
! call METIS_PartGraphRecursive(nspec, xadj(1), adjncy(1), vwgt(0), adjwgt(0), wgtflag, num_start, nparts, &
! metis_options, edgecut, part(1));
@@ -195,58 +299,62 @@
if (ierr /= 0) then
stop 'ERROR : MAIN : Cannot destroy strat'
endif
-
-
+
! local number of each element for each partition
- call Construct_glob2loc_elmnts(nspec, part, nparts, glob2loc_elmnts)
+ call Construct_glob2loc_elmnts(nspec, part, glob2loc_elmnts)
! local number of each node for each partition
- call Construct_glob2loc_nodes(nspec, nnodes,nsize, nnodes_elmnts, nodes_elmnts, part, nparts, &
+ call Construct_glob2loc_nodes(nspec, nnodes,nsize, nnodes_elmnts, nodes_elmnts, part, &
glob2loc_nodes_nparts, glob2loc_nodes_parts, glob2loc_nodes)
- nb_materials = 1
- allocate(cs(nb_materials))
- allocate(num_material(nspec))
- cs(:) = 1000.d0
- num_material(:) = 1
+ call Construct_interfaces(nspec, sup_neighbour, part, elmnts, xadj, adjncy, tab_interfaces, &
+ tab_size_interfaces, ninterfaces, count_def_mat, mat_prop(3,:), mat(1,:))
- call Construct_interfaces(nspec, nparts, sup_neighbour, part, elmnts, xadj, adjncy, tab_interfaces, &
- tab_size_interfaces, ninterfaces, nb_materials, cs, num_material)
-
allocate(my_interfaces(0:ninterfaces-1))
allocate(my_nb_interfaces(0:ninterfaces-1))
-
-
do ipart = 0, nparts-1
- !write(prname, "('/Database',i5.5)") ipart
write(prname, "(i6.6,'_Database')") ipart
open(unit=15,file='./OUTPUT_FILES/proc'//prname,status='unknown', action='write', form='formatted')
-
+
call write_glob2loc_nodes_database(15, ipart, nnodes_loc, nodes_coords, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
glob2loc_nodes, nnodes, 1)
call write_partition_database(15, ipart, nspec_loc, nspec, elmnts, glob2loc_elmnts, glob2loc_nodes_nparts, &
- glob2loc_nodes_parts, glob2loc_nodes, part, num_material, ngnod, 1)
+ glob2loc_nodes_parts, glob2loc_nodes, part, mat, ngnod, 1)
write(15,*) nnodes_loc
call write_glob2loc_nodes_database(15, ipart, nnodes_loc, nodes_coords, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
glob2loc_nodes, nnodes, 2)
+
+ call write_material_properties_database(15,count_def_mat,count_undef_mat, mat_prop, undef_mat_prop)
+
write(15,*) nspec_loc
call write_partition_database(15, ipart, nspec_loc, nspec, elmnts, glob2loc_elmnts, glob2loc_nodes_nparts, &
- glob2loc_nodes_parts, glob2loc_nodes, part, num_material, ngnod, 2)
+ glob2loc_nodes_parts, glob2loc_nodes, part, mat, ngnod, 2)
- call Write_interfaces_database(15, tab_interfaces, tab_size_interfaces, nparts, ipart, ninterfaces, &
+ call write_boundaries_database(15, ipart, nspec, nspec2D_xmin, nspec2D_xmax, nspec2D_ymin, &
+ nspec2D_ymax, nspec2D_bottom, nspec2D_top, ibelm_xmin, ibelm_xmax, ibelm_ymin, &
+ ibelm_ymax, ibelm_bottom, ibelm_top, nodes_ibelm_xmin, nodes_ibelm_xmax, nodes_ibelm_ymin, &
+ nodes_ibelm_ymax, nodes_ibelm_bottom, nodes_ibelm_top, &
+ glob2loc_elmnts, glob2loc_nodes_nparts, glob2loc_nodes_parts, glob2loc_nodes, part)
+
+ call Write_interfaces_database(15, tab_interfaces, tab_size_interfaces, ipart, ninterfaces, &
my_ninterface, my_interfaces, my_nb_interfaces, glob2loc_elmnts, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
glob2loc_nodes, 1)
write(15,*) my_ninterface, maxval(my_nb_interfaces)
- call Write_interfaces_database(15, tab_interfaces, tab_size_interfaces, nparts, ipart, ninterfaces, &
+ call Write_interfaces_database(15, tab_interfaces, tab_size_interfaces, ipart, ninterfaces, &
my_ninterface, my_interfaces, my_nb_interfaces, glob2loc_elmnts, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
glob2loc_nodes, 2)
-
+
+
close(15)
+
+ end do
- enddo
+
+end program pre_meshfem3D
- end program decompose_mesh_SCOTCH
+
+
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90 2009-07-10 16:39:24 UTC (rev 15452)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/part_decompose_mesh_SCOTCH.f90 2009-07-10 21:02:53 UTC (rev 15453)
@@ -116,12 +116,11 @@
!--------------------------------------------------
! construct local numbering for the elements in each partition
!--------------------------------------------------
- subroutine Construct_glob2loc_elmnts(nelmnts, part, nparts, glob2loc_elmnts)
+ subroutine Construct_glob2loc_elmnts(nelmnts, part, glob2loc_elmnts)
include './constants_decompose_mesh_SCOTCH.h'
integer(long), intent(in) :: nelmnts
- integer, intent(in) :: nparts
integer, dimension(0:nelmnts-1), intent(in) :: part
integer, dimension(:), pointer :: glob2loc_elmnts
@@ -151,13 +150,13 @@
!--------------------------------------------------
! construct local numbering for the nodes in each partition
!--------------------------------------------------
- subroutine Construct_glob2loc_nodes(nelmnts, nnodes, nsize, nnodes_elmnts, nodes_elmnts, part, nparts, &
+ subroutine Construct_glob2loc_nodes(nelmnts, nnodes, nsize, nnodes_elmnts, nodes_elmnts, part, &
glob2loc_nodes_nparts, glob2loc_nodes_parts, glob2loc_nodes)
include './constants_decompose_mesh_SCOTCH.h'
integer(long), intent(in) :: nelmnts, nsize
- integer, intent(in) :: nnodes, nparts
+ integer, intent(in) :: nnodes
integer, dimension(0:nelmnts-1), intent(in) :: part
integer, dimension(0:nnodes-1), intent(in) :: nnodes_elmnts
integer, dimension(0:nsize*nnodes-1), intent(in) :: nodes_elmnts
@@ -237,13 +236,13 @@
! 1/ first element, 2/ second element, 3/ number of common nodes, 4/ first node,
! 5/ second node, if relevant.
! No interface between acoustic and elastic elements.
+ ! Elements with undefined material are considered as elastic elements.
!--------------------------------------------------
- subroutine Construct_interfaces(nelmnts, nparts, sup_neighbour, part, elmnts, xadj, adjncy, &
+ subroutine Construct_interfaces(nelmnts, sup_neighbour, part, elmnts, xadj, adjncy, &
tab_interfaces, tab_size_interfaces, ninterfaces, nb_materials, cs_material, num_material)
include './constants_decompose_mesh_SCOTCH.h'
- integer, intent(in) :: nparts
integer(long), intent(in) :: nelmnts, sup_neighbour
integer, dimension(0:nelmnts-1), intent(in) :: part
integer, dimension(0:esize*nelmnts-1), intent(in) :: elmnts
@@ -278,14 +277,22 @@
do num_part_bis = num_part+1, nparts-1
do el = 0, nelmnts-1
if ( part(el) == num_part ) then
- if ( cs_material(num_material(el+1)) < TINYVAL) then
- is_acoustic_el = .true.
+ if(num_material(el+1) > 0) then
+ if ( cs_material(num_material(el+1)) < TINYVAL) then
+ is_acoustic_el = .true.
+ else
+ is_acoustic_el = .false.
+ end if
else
is_acoustic_el = .false.
end if
do el_adj = xadj(el), xadj(el+1)-1
- if ( cs_material(num_material(adjncy(el_adj)+1)) < TINYVAL) then
- is_acoustic_el_adj = .true.
+ if(num_material(adjncy(el_adj)+1) > 0) then
+ if ( cs_material(num_material(adjncy(el_adj)+1)) < TINYVAL) then
+ is_acoustic_el_adj = .true.
+ else
+ is_acoustic_el_adj = .false.
+ end if
else
is_acoustic_el_adj = .false.
end if
@@ -313,14 +320,22 @@
do num_part_bis = num_part+1, nparts-1
do el = 0, nelmnts-1
if ( part(el) == num_part ) then
- if ( cs_material(num_material(el+1)) < TINYVAL) then
- is_acoustic_el = .true.
+ if(num_material(el+1) > 0) then
+ if ( cs_material(num_material(el+1)) < TINYVAL) then
+ is_acoustic_el = .true.
+ else
+ is_acoustic_el = .false.
+ end if
else
is_acoustic_el = .false.
end if
do el_adj = xadj(el), xadj(el+1)-1
- if ( cs_material(num_material(adjncy(el_adj)+1)) < TINYVAL) then
- is_acoustic_el_adj = .true.
+ if(num_material(adjncy(el_adj)+1) > 0) then
+ if ( cs_material(num_material(adjncy(el_adj)+1)) < TINYVAL) then
+ is_acoustic_el_adj = .true.
+ else
+ is_acoustic_el_adj = .false.
+ end if
else
is_acoustic_el_adj = .false.
end if
@@ -400,8 +415,272 @@
end subroutine Write_glob2loc_nodes_database
+ !--------------------------------------------------
+ ! Write material properties in the Database
+ !--------------------------------------------------
+ subroutine write_material_properties_database(IIN_database,count_def_mat,count_undef_mat, mat_prop, undef_mat_prop)
+ integer, intent(in) :: IIN_database
+ integer, intent(in) :: count_def_mat,count_undef_mat
+ double precision, dimension(5,count_def_mat) :: mat_prop
+ character (len=30), dimension(5,count_undef_mat) :: undef_mat_prop
+ integer :: i
+
+ write(IIN_database,*) count_def_mat,count_undef_mat
+ do i = 1, count_def_mat
+ write(IIN_database,*) mat_prop(1,i), mat_prop(2,i), mat_prop(3,i), mat_prop(4,i), mat_prop(5,i)
+ end do
+ do i = 1, count_undef_mat
+ write(IIN_database,*) trim(undef_mat_prop(1,i)),trim(undef_mat_prop(2,i)),trim(undef_mat_prop(3,i)), &
+ trim(undef_mat_prop(4,i)),trim(undef_mat_prop(5,i))
+ end do
+
+ end subroutine write_material_properties_database
+
+
!--------------------------------------------------
+ ! Write elements on boundaries (and their four nodes on boundaries) pertaining to iproc partition in the corresponding Database
+ !--------------------------------------------------
+ subroutine write_boundaries_database(IIN_database, iproc, nelmnts, nspec2D_xmin, nspec2D_xmax, &
+ nspec2D_ymin, nspec2D_ymax, nspec2D_bottom, nspec2D_top, &
+ ibelm_xmin, ibelm_xmax, ibelm_ymin, ibelm_ymax, ibelm_bottom, ibelm_top, &
+ nodes_ibelm_xmin, nodes_ibelm_xmax, nodes_ibelm_ymin, nodes_ibelm_ymax, nodes_ibelm_bottom, nodes_ibelm_top, &
+ glob2loc_elmnts, glob2loc_nodes_nparts, glob2loc_nodes_parts, glob2loc_nodes, part)
+
+ include './constants_decompose_mesh_SCOTCH.h'
+
+ integer, intent(in) :: IIN_database
+ integer, intent(in) :: iproc
+ integer(long), intent(in) :: nelmnts
+ integer, intent(in) :: nspec2D_xmin, nspec2D_xmax, nspec2D_ymin, nspec2D_ymax, nspec2D_bottom, nspec2D_top
+ integer, dimension(nspec2D_xmin), intent(in) :: ibelm_xmin
+ integer, dimension(nspec2D_xmax), intent(in) :: ibelm_xmax
+ integer, dimension(nspec2D_ymin), intent(in) :: ibelm_ymin
+ integer, dimension(nspec2D_ymax), intent(in) :: ibelm_ymax
+ integer, dimension(nspec2D_bottom), intent(in) :: ibelm_bottom
+ integer, dimension(nspec2D_top), intent(in) :: ibelm_top
+
+ integer, dimension(4,nspec2D_xmin), intent(in) :: nodes_ibelm_xmin
+ integer, dimension(4,nspec2D_xmax), intent(in) :: nodes_ibelm_xmax
+ integer, dimension(4,nspec2D_ymin), intent(in) :: nodes_ibelm_ymin
+ integer, dimension(4,nspec2D_ymax), intent(in) :: nodes_ibelm_ymax
+ integer, dimension(4,nspec2D_bottom), intent(in) :: nodes_ibelm_bottom
+ integer, dimension(4,nspec2D_top), intent(in) :: nodes_ibelm_top
+ integer, dimension(:), pointer :: glob2loc_elmnts
+ integer, dimension(:), pointer :: glob2loc_nodes_nparts
+ integer, dimension(:), pointer :: glob2loc_nodes_parts
+ integer, dimension(:), pointer :: glob2loc_nodes
+ integer, dimension(1:nelmnts) :: part
+
+ integer :: i,j
+ integer :: loc_node1, loc_node2, loc_node3, loc_node4
+ integer :: loc_nspec2D_xmin,loc_nspec2D_xmax,loc_nspec2D_ymin,loc_nspec2D_ymax,loc_nspec2D_bottom,loc_nspec2D_top
+
+
+ loc_nspec2D_xmin = 0
+ do i=1,nspec2D_xmin
+ if(part(ibelm_xmin(i)) == iproc) then
+ loc_nspec2D_xmin = loc_nspec2D_xmin + 1
+ end if
+ end do
+ write(IIN_database,*) 1, loc_nspec2D_xmin
+ loc_nspec2D_xmax = 0
+ do i=1,nspec2D_xmax
+ if(part(ibelm_xmax(i)) == iproc) then
+ loc_nspec2D_xmax = loc_nspec2D_xmax + 1
+ end if
+ end do
+ write(IIN_database,*) 2, loc_nspec2D_xmax
+ loc_nspec2D_ymin = 0
+ do i=1,nspec2D_ymin
+ if(part(ibelm_ymin(i)) == iproc) then
+ loc_nspec2D_ymin = loc_nspec2D_ymin + 1
+ end if
+ end do
+ write(IIN_database,*) 3, loc_nspec2D_ymin
+ loc_nspec2D_ymax = 0
+ do i=1,nspec2D_ymax
+ if(part(ibelm_ymax(i)) == iproc) then
+ loc_nspec2D_ymax = loc_nspec2D_ymax + 1
+ end if
+ end do
+ write(IIN_database,*) 4, loc_nspec2D_ymax
+ loc_nspec2D_bottom = 0
+ do i=1,nspec2D_bottom
+ if(part(ibelm_bottom(i)) == iproc) then
+ loc_nspec2D_bottom = loc_nspec2D_bottom + 1
+ end if
+ end do
+ write(IIN_database,*) 5, loc_nspec2D_bottom
+ loc_nspec2D_top = 0
+ do i=1,nspec2D_top
+ if(part(ibelm_top(i)) == iproc) then
+ loc_nspec2D_top = loc_nspec2D_top + 1
+ end if
+ end do
+ write(IIN_database,*) 6, loc_nspec2D_top
+
+ do i=1,nspec2D_xmin
+ if(part(ibelm_xmin(i)) == iproc) then
+ do j = glob2loc_nodes_nparts(nodes_ibelm_xmin(1,i)-1), glob2loc_nodes_nparts(nodes_ibelm_xmin(1,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node1 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_xmin(2,i)-1), glob2loc_nodes_nparts(nodes_ibelm_xmin(2,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node2 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_xmin(3,i)-1), glob2loc_nodes_nparts(nodes_ibelm_xmin(3,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node3 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_xmin(4,i)-1), glob2loc_nodes_nparts(nodes_ibelm_xmin(4,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node4 = glob2loc_nodes(j)+1
+ end if
+ end do
+ write(IIN_database,*) glob2loc_elmnts(ibelm_xmin(i)-1)+1, loc_node1, loc_node2, loc_node3, loc_node4
+ end if
+ end do
+
+ do i=1,nspec2D_xmax
+ if(part(ibelm_xmax(i)) == iproc) then
+ do j = glob2loc_nodes_nparts(nodes_ibelm_xmax(1,i)-1), glob2loc_nodes_nparts(nodes_ibelm_xmax(1,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node1 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_xmax(2,i)-1), glob2loc_nodes_nparts(nodes_ibelm_xmax(2,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node2 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_xmax(3,i)-1), glob2loc_nodes_nparts(nodes_ibelm_xmax(3,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node3 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_xmax(4,i)-1), glob2loc_nodes_nparts(nodes_ibelm_xmax(4,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node4 = glob2loc_nodes(j)+1
+ end if
+ end do
+ write(IIN_database,*) glob2loc_elmnts(ibelm_xmax(i))+1, loc_node1, loc_node2, loc_node3, loc_node4
+ end if
+ end do
+
+ do i=1,nspec2D_ymin
+ if(part(ibelm_ymin(i)) == iproc) then
+ do j = glob2loc_nodes_nparts(nodes_ibelm_ymin(1,i)-1), glob2loc_nodes_nparts(nodes_ibelm_ymin(1,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node1 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_ymin(2,i)-1), glob2loc_nodes_nparts(nodes_ibelm_ymin(2,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node2 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_ymin(3,i)-1), glob2loc_nodes_nparts(nodes_ibelm_ymin(3,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node3 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_ymin(4,i)-1), glob2loc_nodes_nparts(nodes_ibelm_ymin(4,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node4 = glob2loc_nodes(j)+1
+ end if
+ end do
+ write(IIN_database,*) glob2loc_elmnts(ibelm_ymin(i))+1, loc_node1, loc_node2, loc_node3, loc_node4
+ end if
+ end do
+
+ do i=1,nspec2D_ymax
+ if(part(ibelm_ymax(i)) == iproc) then
+ do j = glob2loc_nodes_nparts(nodes_ibelm_ymax(1,i)-1), glob2loc_nodes_nparts(nodes_ibelm_ymax(1,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node1 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_ymax(2,i)-1), glob2loc_nodes_nparts(nodes_ibelm_ymax(2,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node2 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_ymax(3,i)-1), glob2loc_nodes_nparts(nodes_ibelm_ymax(3,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node3 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_ymax(4,i)-1), glob2loc_nodes_nparts(nodes_ibelm_ymax(4,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node4 = glob2loc_nodes(j)+1
+ end if
+ end do
+ write(IIN_database,*) glob2loc_elmnts(ibelm_ymax(i))+1, loc_node1, loc_node2, loc_node3, loc_node4
+ end if
+ end do
+
+ do i=1,nspec2D_bottom
+ if(part(ibelm_bottom(i)) == iproc) then
+ do j = glob2loc_nodes_nparts(nodes_ibelm_bottom(1,i)-1), glob2loc_nodes_nparts(nodes_ibelm_bottom(1,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node1 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_bottom(2,i)-1), glob2loc_nodes_nparts(nodes_ibelm_bottom(2,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node2 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_bottom(3,i)-1), glob2loc_nodes_nparts(nodes_ibelm_bottom(3,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node3 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_bottom(4,i)-1), glob2loc_nodes_nparts(nodes_ibelm_bottom(4,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node4 = glob2loc_nodes(j)+1
+ end if
+ end do
+ write(IIN_database,*) glob2loc_elmnts(ibelm_bottom(i))+1, loc_node1, loc_node2, loc_node3, loc_node4
+ end if
+ end do
+
+ do i=1,nspec2D_top
+ if(part(ibelm_top(i)) == iproc) then
+ do j = glob2loc_nodes_nparts(nodes_ibelm_top(1,i)-1), glob2loc_nodes_nparts(nodes_ibelm_top(1,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node1 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_top(2,i)-1), glob2loc_nodes_nparts(nodes_ibelm_top(2,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node2 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_top(3,i)-1), glob2loc_nodes_nparts(nodes_ibelm_top(3,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node3 = glob2loc_nodes(j)+1
+ end if
+ end do
+ do j = glob2loc_nodes_nparts(nodes_ibelm_top(4,i)-1), glob2loc_nodes_nparts(nodes_ibelm_top(4,i))-1
+ if (glob2loc_nodes_parts(j) == iproc ) then
+ loc_node4 = glob2loc_nodes(j)+1
+ end if
+ end do
+ write(IIN_database,*) glob2loc_elmnts(ibelm_top(i))+1, loc_node1, loc_node2, loc_node3, loc_node4
+ end if
+ end do
+
+
+ end subroutine write_boundaries_database
+
+
+ !--------------------------------------------------
! Write elements (their nodes) pertaining to iproc partition in the corresponding Database
!--------------------------------------------------
subroutine write_partition_database(IIN_database, iproc, nspec, nelmnts, elmnts, glob2loc_elmnts, glob2loc_nodes_nparts, &
@@ -416,7 +695,7 @@
integer, dimension(0:nelmnts-1) :: part
integer, dimension(0:esize*nelmnts-1) :: elmnts
integer, dimension(:), pointer :: glob2loc_elmnts
- integer, dimension(:) :: num_modele
+ integer, dimension(2,nspec) :: num_modele
integer, dimension(:), pointer :: glob2loc_nodes_nparts
integer, dimension(:), pointer :: glob2loc_nodes_parts
integer, dimension(:), pointer :: glob2loc_nodes
@@ -449,7 +728,7 @@
end do
end do
- write(IIN_database,*) glob2loc_elmnts(i)+1, num_modele(i+1), (loc_nodes(k)+1, k=0,ngnod-1)
+ write(IIN_database,*) glob2loc_elmnts(i)+1, num_modele(1,i+1), num_modele(2,i+1),(loc_nodes(k)+1, k=0,ngnod-1)
end if
end do
end if
@@ -459,17 +738,17 @@
-
!--------------------------------------------------
! Write interfaces (element and common nodes) pertaining to iproc partition in the corresponding Database
!--------------------------------------------------
- subroutine write_interfaces_database(IIN_database, tab_interfaces, tab_size_interfaces, nparts, iproc, ninterfaces, &
+ subroutine write_interfaces_database(IIN_database, tab_interfaces, tab_size_interfaces, iproc, ninterfaces, &
my_ninterface, my_interfaces, my_nb_interfaces, glob2loc_elmnts, glob2loc_nodes_nparts, glob2loc_nodes_parts, &
glob2loc_nodes, num_phase)
+ include './constants_decompose_mesh_SCOTCH.h'
+
integer, intent(in) :: IIN_database
integer, intent(in) :: iproc
- integer, intent(in) :: nparts
integer, intent(in) :: ninterfaces
integer, intent(inout) :: my_ninterface
integer, dimension(:), pointer :: tab_size_interfaces
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/scotchf.h
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/scotchf.h 2009-07-10 16:39:24 UTC (rev 15452)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/scotchf.h 2009-07-10 21:02:53 UTC (rev 15453)
@@ -1,4 +1,61 @@
+! ** Copyright 2004,2007 ENSEIRB, INRIA & CNRS
+! **
+! ** This file is part of the Scotch software package for static mapping,
+! ** graph partitioning and sparse matrix ordering.
+! **
+! ** This software is governed by the CeCILL-C license under French law
+! ** and abiding by the rules of distribution of free software. You can
+! ** use, modify and/or redistribute the software under the terms of the
+! ** CeCILL-C license as circulated by CEA, CNRS and INRIA at the following
+! ** URL: "http://www.cecill.info".
+! **
+! ** As a counterpart to the access to the source code and rights to copy,
+! ** modify and redistribute granted by the license, users are provided
+! ** only with a limited warranty and the software's author, the holder of
+! ** the economic rights, and the successive licensors have only limited
+! ** liability.
+! **
+! ** In this respect, the user's attention is drawn to the risks associated
+! ** with loading, using, modifying and/or developing or reproducing the
+! ** software by the user in light of its specific status of free software,
+! ** that may mean that it is complicated to manipulate, and that also
+! ** therefore means that it is reserved for developers and experienced
+! ** professionals having in-depth computer knowledge. Users are therefore
+! ** encouraged to load and test the software's suitability as regards
+! ** their requirements in conditions enabling the security of their
+! ** systems and/or data to be ensured and, more generally, to use and
+! ** operate it in the same conditions as regards security.
+! **
+! ** The fact that you are presently reading this means that you have had
+! ** knowledge of the CeCILL-C license and that you accept its terms.
+! **
+! ************************************************************
+! ** **
+! ** NAME : ptscotchf.h **
+! ** **
+! ** AUTHOR : Francois PELLEGRINI **
+! ** **
+! ** FUNCTION : FORTRAN declaration file for the **
+! ** LibScotch static mapping and sparse **
+! ** matrix block ordering sequential **
+! ** library. **
+! ** **
+! ** DATES : # Version 3.4 : from : 04 feb 2000 **
+! ** to 22 oct 2001 **
+! ** # Version 4.0 : from : 16 jan 2004 **
+! ** to 16 jan 2004 **
+! ** # Version 5.0 : from : 26 apr 2006 **
+! ** to 26 apr 2006 **
+! ** **
+! ************************************************************
+! ** Size definitions for the SCOTCH opaque
+! ** structures. These structures must be
+! ** allocated as arrays of DOUBLEPRECISION
+! ** values for proper padding. The dummy
+! ** sizes are computed at compile-time by
+! ** program "dummysizes".
+
INTEGER SCOTCH_ARCHDIM
INTEGER SCOTCH_DGRAPHDIM
INTEGER SCOTCH_DORDERDIM
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/generate_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/generate_databases.f90 2009-07-10 16:39:24 UTC (rev 15452)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/generate_databases.f90 2009-07-10 21:02:53 UTC (rev 15453)
@@ -235,10 +235,11 @@
character(len=150) prname
integer :: dummy_node
integer :: dummy_elmnt
- integer :: ispec, inode, num_interface, ie
+ integer :: ispec, inode, num_interface, ie,imat !pll
integer :: nnodes_ext_mesh, nelmnts_ext_mesh
integer :: ninterface_ext_mesh
integer :: max_interface_size_ext_mesh
+ integer :: nmat_ext_mesh, nundefMat_ext_mesh !pll
integer, dimension(:), allocatable :: my_neighbours_ext_mesh
integer, dimension(:), allocatable :: my_nelmnts_neighbours_ext_mesh
integer, dimension(:,:,:), allocatable :: my_interfaces_ext_mesh
@@ -246,8 +247,15 @@
integer, dimension(:), allocatable :: nibool_interfaces_ext_mesh
double precision, dimension(:,:), allocatable :: nodes_coords_ext_mesh
integer, dimension(:,:), allocatable :: elmnts_ext_mesh
- integer, dimension(:), allocatable :: mat_ext_mesh
+ integer, dimension(:,:), allocatable :: mat_ext_mesh
+ ! pll
+ double precision, dimension(:,:), allocatable :: materials_ext_mesh
+ integer, dimension(:), allocatable :: ibelm_xmin,ibelm_xmax, ibelm_ymin, ibelm_ymax, ibelm_bottom, ibelm_top
+ integer :: ispec2D, boundary_number
+ integer :: nspec2D_xmin, nspec2D_xmax, nspec2D_ymin, nspec2D_ymax, nspec2D_bottom_ext, nspec2D_top_ext
+ character (len=30), dimension(:,:), allocatable :: undef_mat_prop
+
! ************** PROGRAM STARTS HERE **************
! sizeprocs returns number of processes started (should be equal to NPROC).
@@ -304,63 +312,12 @@
! info about external mesh simulation
! nlegoff -- should be put in compute_parameters and read_parameter_file for clarity
- if (USE_EXTERNAL_MESH) then
- NPROC = sizeprocs
- endif
+ NPROC = sizeprocs
! check that the code is running with the requested nb of processes
if(sizeprocs /= NPROC) call exit_MPI(myrank,'wrong number of MPI processes')
- if (.not. USE_EXTERNAL_MESH) then
-! dynamic allocation of mesh arrays
- allocate(rns(0:2*NER))
-
- allocate(xgrid(0:2*NER,0:2*NEX_PER_PROC_XI,0:2*NEX_PER_PROC_ETA))
- allocate(ygrid(0:2*NER,0:2*NEX_PER_PROC_XI,0:2*NEX_PER_PROC_ETA))
- allocate(zgrid(0:2*NER,0:2*NEX_PER_PROC_XI,0:2*NEX_PER_PROC_ETA))
-
- allocate(addressing(0:NPROC_XI-1,0:NPROC_ETA-1))
- allocate(iproc_xi_slice(0:NPROC-1))
- allocate(iproc_eta_slice(0:NPROC-1))
-
-! clear arrays
- xgrid(:,:,:) = 0.
- ygrid(:,:,:) = 0.
- zgrid(:,:,:) = 0.
-
- iproc_xi_slice(:) = 0
- iproc_eta_slice(:) = 0
-
-! create global slice addressing for solver
if(myrank == 0) then
- open(unit=IOUT,file=trim(OUTPUT_FILES)//'/addressing.txt',status='unknown')
- write(IMAIN,*) 'creating global slice addressing'
- write(IMAIN,*)
- endif
- do iproc_eta=0,NPROC_ETA-1
- do iproc_xi=0,NPROC_XI-1
- iprocnum = iproc_eta * NPROC_XI + iproc_xi
- iproc_xi_slice(iprocnum) = iproc_xi
- iproc_eta_slice(iprocnum) = iproc_eta
- addressing(iproc_xi,iproc_eta) = iprocnum
- if(myrank == 0) write(IOUT,*) iprocnum,iproc_xi,iproc_eta
- enddo
- enddo
- if(myrank == 0) close(IOUT)
-
- if (myrank == 0) then
- write(IMAIN,*) 'Spatial distribution of slice numbers:'
- do iproc_eta = NPROC_ETA-1, 0, -1
- do iproc_xi = 0, NPROC_XI-1, 1
- write(IMAIN,'(i5)',advance='no') addressing(iproc_xi,iproc_eta)
- enddo
- write(IMAIN,'(a1)',advance='yes') ' '
- enddo
- endif
-
- endif ! end of (.not. USE_EXTERNAL_MESH)
-
- if(myrank == 0) then
write(IMAIN,*) 'This is process ',myrank
write(IMAIN,*) 'There are ',sizeprocs,' MPI processes'
write(IMAIN,*) 'Processes are numbered from 0 to ',sizeprocs-1
@@ -391,30 +348,6 @@
! for the number of standard linear solids for attenuation
if(N_SLS /= 3) call exit_MPI(myrank,'number of SLS must be 3')
- if (.not. USE_EXTERNAL_MESH) then
-! check that Poisson's ratio in Gocad block is fine
- if(VP_VS_RATIO_GOCAD_TOP < sqrt(2.) .or. VP_VS_RATIO_GOCAD_BOTTOM < sqrt(2.))&
- call exit_MPI(myrank,'vp/vs ratio in Gocad block is too small')
-
-! check that number of slices is at least 1 in each direction
- if(NPROC_XI < 1) call exit_MPI(myrank,'NPROC_XI must be greater than 1')
- if(NPROC_ETA < 1) call exit_MPI(myrank,'NPROC_ETA must be greater than 1')
-
-! check that mesh can be cut into the right number of slices
-! also check that mesh can be coarsened in depth twice (block size multiple of 8)
- if(USE_REGULAR_MESH) then
- if(mod(NEX_XI,NPROC_XI) /= 0) call exit_MPI(myrank,'NEX_XI must be a multiple of NPROC_XI for a regular mesh')
- if(mod(NEX_ETA,NPROC_ETA) /= 0) call exit_MPI(myrank,'NEX_ETA must be a multiple of NPROC_ETA for a regular mesh')
- else
- if(mod(NEX_XI,8) /= 0) call exit_MPI(myrank,'NEX_XI must be a multiple of 8 for a non-regular mesh')
- if(mod(NEX_ETA,8) /= 0) call exit_MPI(myrank,'NEX_ETA must be a multiple of 8 for a non-regular mesh')
-
- if(mod(NEX_XI/8,NPROC_XI) /= 0) call exit_MPI(myrank,'NEX_XI must be a multiple of 8*NPROC_XI for a non-regular mesh')
- if(mod(NEX_ETA/8,NPROC_ETA) /= 0) call exit_MPI(myrank,'NEX_ETA must be a multiple of 8*NPROC_ETA for a non-regular mesh')
- endif
-
- endif ! end of (.not. USE_EXTERNAL_MESH)
-
if(myrank == 0) then
write(IMAIN,*) 'region selected:'
@@ -530,184 +463,6 @@
close(55)
endif
- if (.not. USE_EXTERNAL_MESH) then
-
-! get addressing for this process
- iproc_xi = iproc_xi_slice(myrank)
- iproc_eta = iproc_eta_slice(myrank)
-
-! number of elements in each slice
- npx = 2*NEX_PER_PROC_XI
- npy = 2*NEX_PER_PROC_ETA
-
- min_elevation = +HUGEVAL
- max_elevation = -HUGEVAL
-
-! fill the region between the cutoff depth and the free surface
- do iy=0,npy
- do ix=0,npx
-
-! define the mesh points on the top and the bottom
-
- xin=dble(ix)/dble(npx)
- x_current = UTM_X_MIN + (dble(iproc_xi)+xin)*(UTM_X_MAX-UTM_X_MIN)/dble(NPROC_XI)
-
- etan=dble(iy)/dble(npy)
- y_current = UTM_Y_MIN + (dble(iproc_eta)+etan)*(UTM_Y_MAX-UTM_Y_MIN)/dble(NPROC_ETA)
-
-! define model between topography surface and fictitious bottom
- if(TOPOGRAPHY) then
-
-! project x and y in UTM back to long/lat since topo file is in long/lat
- call utm_geo(long,lat,x_current,y_current,UTM_PROJECTION_ZONE,IUTM2LONGLAT,SUPPRESS_UTM_PROJECTION)
-
-! get coordinate of corner in bathy/topo 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 = &
- itopo_bathy(icornerlong,icornerlat)*(1.-ratio_xi)*(1.-ratio_eta) + &
- itopo_bathy(icornerlong+1,icornerlat)*ratio_xi*(1.-ratio_eta) + &
- itopo_bathy(icornerlong+1,icornerlat+1)*ratio_xi*ratio_eta + &
- itopo_bathy(icornerlong,icornerlat+1)*(1.-ratio_xi)*ratio_eta
-
- else
-
- elevation = 0.d0
-
- endif
-
- z_top = Z_SURFACE + elevation
- z_bot = - dabs(Z_DEPTH_BLOCK)
-
-! compute global min and max of elevation
- min_elevation = dmin1(min_elevation,elevation)
- max_elevation = dmax1(max_elevation,elevation)
-
-! create vertical point distribution at current horizontal point
- if(BASEMENT_MAP) then
-
-! get coordinate of corner in bathy/topo model
- icorner_x = int((x_current - ORIG_X_BASEMENT) / SPACING_X_BASEMENT) + 1
- icorner_y = int((y_current - ORIG_Y_BASEMENT) / SPACING_Y_BASEMENT) + 1
-
-! avoid edge effects and extend with identical point if outside model
- if(icorner_x < 1) icorner_x = 1
- if(icorner_x > NX_BASEMENT-1) icorner_x = NX_BASEMENT-1
- if(icorner_y < 1) icorner_y = 1
- if(icorner_y > NY_BASEMENT-1) icorner_y = NY_BASEMENT-1
-
-! compute coordinates of corner
- x_corner = ORIG_X_BASEMENT + (icorner_x-1)*SPACING_X_BASEMENT
- y_corner = ORIG_Y_BASEMENT + (icorner_y-1)*SPACING_Y_BASEMENT
-
-! compute ratio for interpolation
- ratio_xi = (x_current - x_corner) / SPACING_X_BASEMENT
- ratio_eta = (y_current - y_corner) / SPACING_Y_BASEMENT
-
-! 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 basement surface at current point
- Z_BASEMENT_SURFACE = &
- z_basement(icorner_x,icorner_y)*(1.-ratio_xi)*(1.-ratio_eta) + &
- z_basement(icorner_x+1,icorner_y)*ratio_xi*(1.-ratio_eta) + &
- z_basement(icorner_x+1,icorner_y+1)*ratio_xi*ratio_eta + &
- z_basement(icorner_x,icorner_y+1)*(1.-ratio_xi)*ratio_eta
-
- else
- Z_BASEMENT_SURFACE = DEPTH_5p5km_SOCAL
- endif
-
-! honor Lupei Zhu's Moho map
- if(MOHO_MAP_LUPEI) then
-
-! project x and y in UTM back to long/lat since topo file is in long/lat
- call utm_geo(long,lat,x_current,y_current,UTM_PROJECTION_ZONE,IUTM2LONGLAT,SUPPRESS_UTM_PROJECTION)
-
-! get coordinate of corner in Moho map
- icornerlong = int((long - ORIG_LONG_MOHO) / DEGREES_PER_CELL_MOHO) + 1
- icornerlat = int((lat - ORIG_LAT_MOHO) / DEGREES_PER_CELL_MOHO) + 1
-
-! avoid edge effects and extend with identical point if outside model
- if(icornerlong < 1) icornerlong = 1
- if(icornerlong > NX_MOHO-1) icornerlong = NX_MOHO-1
- if(icornerlat < 1) icornerlat = 1
- if(icornerlat > NY_MOHO-1) icornerlat = NY_MOHO-1
-
-! compute coordinates of corner
- long_corner = ORIG_LONG_MOHO + (icornerlong-1)*DEGREES_PER_CELL_MOHO
- lat_corner = ORIG_LAT_MOHO + (icornerlat-1)*DEGREES_PER_CELL_MOHO
-
-! compute ratio for interpolation
- ratio_xi = (long - long_corner) / DEGREES_PER_CELL_MOHO
- ratio_eta = (lat - lat_corner) / DEGREES_PER_CELL_MOHO
-
-! 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 Moho depth at current point
- Z_DEPTH_MOHO = &
- - (imoho_depth(icornerlong,icornerlat)*(1.-ratio_xi)*(1.-ratio_eta) + &
- imoho_depth(icornerlong+1,icornerlat)*ratio_xi*(1.-ratio_eta) + &
- imoho_depth(icornerlong+1,icornerlat+1)*ratio_xi*ratio_eta + &
- imoho_depth(icornerlong,icornerlat+1)*(1.-ratio_xi)*ratio_eta)
-
- else
- Z_DEPTH_MOHO = DEPTH_MOHO_SOCAL
- endif
-
-! define vertical spacing of the mesh in case of a non-regular mesh with mesh doublings
- if(.not. USE_REGULAR_MESH) call mesh_vertical(myrank,rns,NER,NER_BOTTOM_MOHO,NER_MOHO_16, &
- NER_16_BASEMENT,NER_BASEMENT_SEDIM,NER_SEDIM, &
-!! DK DK UGLY modif z_top by Emmanuel Chaljub here
-!! DK DK UGLY modif Manu removed z_top, &
- Z_DEPTH_BLOCK,Z_BASEMENT_SURFACE,Z_DEPTH_MOHO,MOHO_MAP_LUPEI)
-
-! fill the volume
- do ir = 0,2*NER
- if(USE_REGULAR_MESH) then
- rn = dble(ir) / dble(2*NER)
- else
- rn = rns(ir)
- endif
- xgrid(ir,ix,iy) = x_current
- ygrid(ir,ix,iy) = y_current
- zgrid(ir,ix,iy) = z_bot*(ONE-rn) + z_top*rn
- enddo
-
- enddo
- enddo
-
- endif ! end of (.not. USE_EXTERNAL_MESH)
-
if(myrank == 0) then
write(IMAIN,*)
write(IMAIN,*) '**************************'
@@ -722,45 +477,105 @@
area_local_top = ZERO
! read databases about external mesh simulation
-! nlegoff --
- if (USE_EXTERNAL_MESH) then
- call create_name_database(prname,myrank,LOCAL_PATH)
- open(unit=IIN,file=prname(1:len_trim(prname))//'Database',status='old',action='read',form='formatted')
- read(IIN,*) nnodes_ext_mesh
- allocate(nodes_coords_ext_mesh(NDIM,nnodes_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)
- enddo
+
+ call create_name_database(prname,myrank,LOCAL_PATH)
+ open(unit=IIN,file=prname(1:len_trim(prname))//'Database',status='old',action='read',form='formatted')
+ read(IIN,*) nnodes_ext_mesh
+ allocate(nodes_coords_ext_mesh(NDIM,nnodes_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)
+ enddo
+
+! read materials' physical properties
+ read(IIN,*) nmat_ext_mesh, nundefMat_ext_mesh
+ allocate(materials_ext_mesh(5,nmat_ext_mesh))
+ allocate(undef_mat_prop(5,nundefMat_ext_mesh))
+ do imat = 1, nmat_ext_mesh
+ 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)
+ end do
+
+ do imat = 1, nundefMat_ext_mesh
+ 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)
+ end do
- read(IIN,*) nelmnts_ext_mesh
- allocate(elmnts_ext_mesh(esize,nelmnts_ext_mesh))
- allocate(mat_ext_mesh(nelmnts_ext_mesh))
- do ispec = 1, nelmnts_ext_mesh
- read(IIN,*) dummy_elmnt, mat_ext_mesh(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)
- enddo
- NSPEC_AB = nelmnts_ext_mesh
+ read(IIN,*) nelmnts_ext_mesh
+ allocate(elmnts_ext_mesh(esize,nelmnts_ext_mesh))
+ allocate(mat_ext_mesh(2,nelmnts_ext_mesh))
+ do ispec = 1, nelmnts_ext_mesh
+ 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)
+ enddo
+ NSPEC_AB = nelmnts_ext_mesh
+
+! read boundaries
+ read(IIN,*) boundary_number ,nspec2D_xmin
+ if(boundary_number /= 1) stop "Error : invalid database file"
+ read(IIN,*) boundary_number ,nspec2D_xmax
+ if(boundary_number /= 2) stop "Error : invalid database file"
+ read(IIN,*) boundary_number ,nspec2D_ymin
+ if(boundary_number /= 3) stop "Error : invalid database file"
+ read(IIN,*) boundary_number ,nspec2D_ymax
+ if(boundary_number /= 4) stop "Error : invalid database file"
+ read(IIN,*) boundary_number ,nspec2D_bottom_ext
+ if(boundary_number /= 5) stop "Error : invalid database file"
+ read(IIN,*) boundary_number ,nspec2D_top_ext
+ if(boundary_number /= 6) stop "Error : invalid database file"
+ NSPEC2DMAX_XMIN_XMAX = max(nspec2D_xmin,nspec2D_xmax)
+ NSPEC2DMAX_YMIN_YMAX = max(nspec2D_ymin,nspec2D_ymax)
+ NSPEC2D_BOTTOM = nspec2D_bottom_ext
+ NSPEC2D_TOP = nspec2D_top_ext
+
+ allocate(ibelm_xmin(nspec2D_xmin))
+ do ispec2D = 1,nspec2D_xmin
+ read(IIN,*) ibelm_xmin(ispec2D)
+ end do
+
+ allocate(ibelm_xmax(nspec2D_xmax))
+ do ispec2D = 1,nspec2D_xmax
+ read(IIN,*) ibelm_xmax(ispec2D)
+ end do
+
+ allocate(ibelm_ymin(nspec2D_ymin))
+ do ispec2D = 1,nspec2D_ymin
+ read(IIN,*) ibelm_ymin(ispec2D)
+ end do
+
+ allocate(ibelm_ymax(nspec2D_ymax))
+ do ispec2D = 1,nspec2D_ymax
+ read(IIN,*) ibelm_ymax(ispec2D)
+ end do
+
+ allocate(ibelm_bottom(nspec2D_bottom_ext))
+ do ispec2D = 1,nspec2D_bottom_ext
+ read(IIN,*) ibelm_bottom(ispec2D)
+ end do
- read(IIN,*) ninterface_ext_mesh, max_interface_size_ext_mesh
- allocate(my_neighbours_ext_mesh(ninterface_ext_mesh))
- allocate(my_nelmnts_neighbours_ext_mesh(ninterface_ext_mesh))
- allocate(my_interfaces_ext_mesh(6,max_interface_size_ext_mesh,ninterface_ext_mesh))
- allocate(ibool_interfaces_ext_mesh(NGLLX*NGLLX*max_interface_size_ext_mesh,ninterface_ext_mesh))
- allocate(nibool_interfaces_ext_mesh(ninterface_ext_mesh))
- do num_interface = 1, ninterface_ext_mesh
- read(IIN,*) my_neighbours_ext_mesh(num_interface), my_nelmnts_neighbours_ext_mesh(num_interface)
- do ie = 1, my_nelmnts_neighbours_ext_mesh(num_interface)
+ allocate(ibelm_top(nspec2D_top_ext))
+ do ispec2D = 1,nspec2D_top_ext
+ read(IIN,*) ibelm_top(ispec2D)
+ end do
+
+ read(IIN,*) ninterface_ext_mesh, max_interface_size_ext_mesh
+ allocate(my_neighbours_ext_mesh(ninterface_ext_mesh))
+ allocate(my_nelmnts_neighbours_ext_mesh(ninterface_ext_mesh))
+ allocate(my_interfaces_ext_mesh(6,max_interface_size_ext_mesh,ninterface_ext_mesh))
+ allocate(ibool_interfaces_ext_mesh(NGLLX*NGLLX*max_interface_size_ext_mesh,ninterface_ext_mesh))
+ allocate(nibool_interfaces_ext_mesh(ninterface_ext_mesh))
+ do num_interface = 1, ninterface_ext_mesh
+ read(IIN,*) my_neighbours_ext_mesh(num_interface), my_nelmnts_neighbours_ext_mesh(num_interface)
+ do ie = 1, my_nelmnts_neighbours_ext_mesh(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)
- enddo
- enddo
+ enddo
+ enddo
+
+ close(IIN)
- close(IIN)
- endif
-
! assign theoretical number of elements
nspec = NSPEC_AB
@@ -781,35 +596,18 @@
if(ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays')
! create all the regions of the mesh
- if (USE_EXTERNAL_MESH) then
call create_regions_mesh_ext_mesh(ibool, &
- xstore,ystore,zstore,nspec, &
- npointot,myrank,LOCAL_PATH, &
- nnodes_ext_mesh,nelmnts_ext_mesh, &
- nodes_coords_ext_mesh,elmnts_ext_mesh,mat_ext_mesh,static_memory_size, &
- ninterface_ext_mesh,max_interface_size_ext_mesh, &
- my_neighbours_ext_mesh,my_nelmnts_neighbours_ext_mesh,my_interfaces_ext_mesh, &
- ibool_interfaces_ext_mesh,nibool_interfaces_ext_mesh)
+ xstore,ystore,zstore,nspec,npointot,myrank,LOCAL_PATH, &
+ nnodes_ext_mesh,nelmnts_ext_mesh, &
+ nodes_coords_ext_mesh,elmnts_ext_mesh,static_memory_size,mat_ext_mesh,materials_ext_mesh, &
+ nmat_ext_mesh,undef_mat_prop,nundefMat_ext_mesh,ninterface_ext_mesh,max_interface_size_ext_mesh, &
+ my_neighbours_ext_mesh,my_nelmnts_neighbours_ext_mesh,my_interfaces_ext_mesh, &
+ ibool_interfaces_ext_mesh,nibool_interfaces_ext_mesh, &
+ nspec2D_xmin, nspec2D_xmax, nspec2D_ymin, nspec2D_ymax, NSPEC2D_BOTTOM, NSPEC2D_TOP,&
+ NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX, &
+ ibelm_xmin,ibelm_xmax, ibelm_ymin, ibelm_ymax, ibelm_bottom, ibelm_top)
+
- else
- call create_regions_mesh(xgrid,ygrid,zgrid,ibool,idoubling, &
- xstore,ystore,zstore,npx,npy, &
- iproc_xi,iproc_eta,nspec, &
- volume_local,area_local_bottom,area_local_top, &
- NGLOB_AB,npointot, &
- NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM,NER_SEDIM,NER, &
- NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
- NSPEC2DMAX_XMIN_XMAX, &
- NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
- HARVARD_3D_GOCAD_MODEL,NPROC_XI,NPROC_ETA,NSPEC2D_A_XI,NSPEC2D_B_XI, &
- NSPEC2D_A_ETA,NSPEC2D_B_ETA,myrank,LOCAL_PATH, &
- UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK,UTM_PROJECTION_ZONE, &
- HAUKSSON_REGIONAL_MODEL,OCEANS, &
- VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM, &
- IMPOSE_MINIMUM_VP_GOCAD,THICKNESS_TAPER_BLOCK_HR,THICKNESS_TAPER_BLOCK_MR,MOHO_MAP_LUPEI, &
- ANISOTROPY,SAVE_MESH_FILES,SUPPRESS_UTM_PROJECTION, &
- ORIG_LAT_TOPO,ORIG_LONG_TOPO,DEGREES_PER_CELL_TOPO,NX_TOPO,NY_TOPO,USE_REGULAR_MESH)
- endif
! print min and max of topography included
if(TOPOGRAPHY) then
@@ -900,51 +698,53 @@
UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,ATTENUATION,ANISOTROPY,NSTEP, &
NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,SIMULATION_TYPE,static_memory_size)
- call get_value_string(rec_filename, 'solver.STATIONS', 'DATA/STATIONS')
- call get_value_string(filtered_rec_filename, 'solver.STATIONS_FILTERED', 'DATA/STATIONS_FILTERED')
+! call get_value_string(rec_filename, 'solver.STATIONS', 'DATA/STATIONS')
+! call get_value_string(filtered_rec_filename, 'solver.STATIONS_FILTERED', 'DATA/STATIONS_FILTERED')
! 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)
+! 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)
! filter list of stations, only retain stations that are in the model
- nrec_filtered = 0
- open(unit=IIN,file=rec_filename,status='old',action='read')
- do irec = 1,nrec
- read(IIN,*) station_name,network_name,stlat,stlon,stele,stbur
- if((stlat >= LATITUDE_MIN .and. stlat <= LATITUDE_MAX .and. stlon >= LONGITUDE_MIN .and. stlon <= LONGITUDE_MAX) &
- .or. USE_EXTERNAL_MESH) &
- nrec_filtered = nrec_filtered + 1
- enddo
- close(IIN)
+! nrec_filtered = 0
+! open(unit=IIN,file=rec_filename,status='old',action='read')
+! do irec = 1,nrec
+! read(IIN,*) station_name,network_name,stlat,stlon,stele,stbur
+! if((stlat >= LATITUDE_MIN .and. stlat <= LATITUDE_MAX .and. stlon >= LONGITUDE_MIN .and. stlon <= LONGITUDE_MAX) &
+! .or. USE_EXTERNAL_MESH) &
+! nrec_filtered = nrec_filtered + 1
+! enddo
+! close(IIN)
- write(IMAIN,*)
- write(IMAIN,*) 'there are ',nrec,' stations in file ', trim(rec_filename)
- write(IMAIN,*) 'saving ',nrec_filtered,' stations inside the model in file ', trim(filtered_rec_filename)
- write(IMAIN,*) 'excluding ',nrec - nrec_filtered,' stations located outside the model'
- write(IMAIN,*)
+! write(IMAIN,*)
+! write(IMAIN,*) 'there are ',nrec,' stations in file ', trim(rec_filename)
+! write(IMAIN,*) 'saving ',nrec_filtered,' stations inside the model in file ', trim(filtered_rec_filename)
+! write(IMAIN,*) 'excluding ',nrec - nrec_filtered,' stations located outside the model'
+! write(IMAIN,*)
- if(nrec_filtered < 1) call exit_MPI(myrank,'need at least one station in the model')
+! if(nrec_filtered < 1) call exit_MPI(myrank,'need at least one station in the model')
- open(unit=IIN,file=rec_filename,status='old',action='read')
- open(unit=IOUT,file=filtered_rec_filename,status='unknown')
+! if(nrec < 1) call exit_MPI(myrank,'need at least one station in the model')
- do irec = 1,nrec
- read(IIN,*) station_name,network_name,stlat,stlon,stele,stbur
- if((stlat >= LATITUDE_MIN .and. stlat <= LATITUDE_MAX .and. stlon >= LONGITUDE_MIN .and. stlon <= LONGITUDE_MAX) &
- .or. USE_EXTERNAL_MESH) &
- write(IOUT,*) station_name(1:len_trim(station_name)),' ',network_name(1:len_trim(network_name)),' ', &
- sngl(stlat),' ',sngl(stlon), ' ', sngl(stele), ' ', sngl(stbur)
- enddo
+! open(unit=IIN,file=rec_filename,status='old',action='read')
+! open(unit=IOUT,file=filtered_rec_filename,status='unknown')
- close(IIN)
- close(IOUT)
+! do irec = 1,nrec
+! read(IIN,*) station_name,network_name,stlat,stlon,stele,stbur
+! if((stlat >= LATITUDE_MIN .and. stlat <= LATITUDE_MAX .and. stlon >= LONGITUDE_MIN .and. stlon <= LONGITUDE_MAX) &
+! .or. USE_EXTERNAL_MESH) &
+! write(IOUT,*) station_name(1:len_trim(station_name)),' ',network_name(1:len_trim(network_name)),' ', &
+! sngl(stlat),' ',sngl(stlon), ' ', sngl(stele), ' ', sngl(stbur)
+! enddo
+! close(IIN)
+! close(IOUT)
+
endif ! end of section executed by main process only
! elapsed time since beginning of mesh generation
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/get_absorb.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/get_absorb.f90 2009-07-10 16:39:24 UTC (rev 15452)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/get_absorb.f90 2009-07-10 21:02:53 UTC (rev 15453)
@@ -156,3 +156,113 @@
end subroutine get_absorb
+
+
+
+ subroutine get_absorb_ext_mesh(myrank,iboun,nspec, &
+ nimin,nimax,njmin,njmax,nkmin_xi,nkmin_eta, &
+ NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM)
+
+! put Stacey back, here define overlap flags
+
+ implicit none
+
+ include "constants.h"
+
+ integer nspec,myrank
+
+ integer NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM
+
+ integer nimin(2,NSPEC2DMAX_YMIN_YMAX),nimax(2,NSPEC2DMAX_YMIN_YMAX)
+ integer njmin(2,NSPEC2DMAX_XMIN_XMAX),njmax(2,NSPEC2DMAX_XMIN_XMAX)
+ integer nkmin_xi(2,NSPEC2DMAX_XMIN_XMAX),nkmin_eta(2,NSPEC2DMAX_YMIN_YMAX)
+
+ logical iboun(6,nspec)
+
+! global element numbering
+ integer ispecg
+
+! counters to keep track of the number of elements on each of the
+! five absorbing boundaries
+ integer ispecb1,ispecb2,ispecb3,ispecb4,ispecb5
+
+ ispecb1=0
+ ispecb2=0
+ ispecb3=0
+ ispecb4=0
+ ispecb5=0
+
+ do ispecg=1,nspec
+
+! determine if the element falls on an absorbing boundary
+
+ if(iboun(1,ispecg)) then
+
+! on boundary 1: xmin
+ ispecb1=ispecb1+1
+
+! this is useful even if it is constant because it can be zero inside the slices
+ njmin(1,ispecb1)=1
+ njmax(1,ispecb1)=NGLLY
+
+! check for ovelap with other boundaries
+ nkmin_xi(1,ispecb1)=1
+ if(iboun(5,ispecg)) nkmin_xi(1,ispecb1)=2
+ endif
+
+ if(iboun(2,ispecg)) then
+
+! on boundary 2: xmax
+ ispecb2=ispecb2+1
+
+! this is useful even if it is constant because it can be zero inside the slices
+ njmin(2,ispecb2)=1
+ njmax(2,ispecb2)=NGLLY
+
+! check for ovelap with other boundaries
+ nkmin_xi(2,ispecb2)=1
+ if(iboun(5,ispecg)) nkmin_xi(2,ispecb2)=2
+ endif
+
+ if(iboun(3,ispecg)) then
+
+! on boundary 3: ymin
+ ispecb3=ispecb3+1
+
+! check for ovelap with other boundaries
+ nimin(1,ispecb3)=1
+ if(iboun(1,ispecg)) nimin(1,ispecb3)=2
+ nimax(1,ispecb3)=NGLLX
+ if(iboun(2,ispecg)) nimax(1,ispecb3)=NGLLX-1
+ nkmin_eta(1,ispecb3)=1
+ if(iboun(5,ispecg)) nkmin_eta(1,ispecb3)=2
+ endif
+
+ if(iboun(4,ispecg)) then
+
+! on boundary 4: ymax
+ ispecb4=ispecb4+1
+
+! check for ovelap with other boundaries
+ nimin(2,ispecb4)=1
+ if(iboun(1,ispecg)) nimin(2,ispecb4)=2
+ nimax(2,ispecb4)=NGLLX
+ if(iboun(2,ispecg)) nimax(2,ispecb4)=NGLLX-1
+ nkmin_eta(2,ispecb4)=1
+ if(iboun(5,ispecg)) nkmin_eta(2,ispecb4)=2
+ endif
+
+! on boundary 5: bottom
+ if(iboun(5,ispecg)) ispecb5=ispecb5+1
+
+ enddo
+
+! check theoretical value of elements at the bottom
+ if(ispecb5 /= NSPEC2D_BOTTOM) &
+ call exit_MPI(myrank,'ispecb5 should equal NSPEC2D_BOTTOM in absorbing boundary detection')
+
+ end subroutine get_absorb_ext_mesh
+
+
+
+
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/locate_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/locate_receivers.f90 2009-07-10 16:39:24 UTC (rev 15452)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/locate_receivers.f90 2009-07-10 21:02:53 UTC (rev 15453)
@@ -201,7 +201,7 @@
! convert station location to UTM
call utm_geo(stlon(irec),stlat(irec),stutm_x(irec),stutm_y(irec),UTM_PROJECTION_ZONE,ILONGLAT2UTM, &
- (SUPPRESS_UTM_PROJECTION .or. USE_EXTERNAL_MESH))
+ .true.)
! compute horizontal distance between source and receiver in km
horiz_dist(irec) = dsqrt((stutm_y(irec)-utm_y_source)**2 + (stutm_x(irec)-utm_x_source)**2) / 1000.
@@ -231,61 +231,14 @@
nu(3,2,irec) = 0.d0
nu(3,3,irec) = 1.d0
- if (.not. USE_EXTERNAL_MESH) then
-! compute elevation of topography at the receiver location
-! we assume that receivers are always at the surface i.e. not buried
- if(TOPOGRAPHY) then
-! get coordinate of corner in bathy/topo model
- icornerlong = int((stlon(irec) - ORIG_LONG_TOPO) / DEGREES_PER_CELL_TOPO) + 1
- icornerlat = int((stlat(irec) - 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 = (stlon(irec) - long_corner) / DEGREES_PER_CELL_TOPO
- ratio_eta = (stlat(irec) - 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 = &
- itopo_bathy(icornerlong,icornerlat)*(1.-ratio_xi)*(1.-ratio_eta) + &
- itopo_bathy(icornerlong+1,icornerlat)*ratio_xi*(1.-ratio_eta) + &
- itopo_bathy(icornerlong+1,icornerlat+1)*ratio_xi*ratio_eta + &
- itopo_bathy(icornerlong,icornerlat+1)*(1.-ratio_xi)*ratio_eta
-
- else
- elevation(irec) = 0.d0
- endif
-
-! compute the Cartesian position of the receiver
x_target(irec) = stutm_x(irec)
y_target(irec) = stutm_y(irec)
- z_target(irec) = elevation(irec) - stbur(irec)
+ z_target(irec) = stbur(irec)
if (myrank == 0) write(IOVTK,*) x_target(irec), y_target(irec), z_target(irec)
- else
+
- x_target(irec) = stutm_x(irec)
- y_target(irec) = stutm_y(irec)
- z_target(irec) = stbur(irec)
- if (myrank == 0) write(IOVTK,*) x_target(irec), y_target(irec), z_target(irec)
-
- endif ! of if (.not. USE_EXTERNAL_MESH)
-
! examine top of the elements only (receivers always at the surface)
! k = NGLLZ
@@ -323,7 +276,7 @@
iglob = ibool(i,j,k,ispec)
- if (USE_EXTERNAL_MESH .and. (.not. RECVS_CAN_BE_BURIED_EXT_MESH)) then
+ if (.not. RECVS_CAN_BE_BURIED_EXT_MESH) then
if ((.not. iglob_is_surface_external_mesh(iglob)) .or. (.not. ispec_is_surface_external_mesh(ispec))) then
cycle
endif
@@ -366,7 +319,7 @@
endif
! get normal to the face of the hexaedra if receiver is on the surface
- if (USE_EXTERNAL_MESH .and. (.not. RECVS_CAN_BE_BURIED_EXT_MESH) .and. &
+ if ((.not. RECVS_CAN_BE_BURIED_EXT_MESH) .and. &
.not. (ispec_selected_rec(irec) == 0)) then
pt0_ix = -1
pt0_iy = -1
@@ -517,7 +470,7 @@
nu(3,3,irec) = 1.d0
endif
- endif ! of if (USE_EXTERNAL_MESH .and. (.not. RECVS_CAN_BE_BURIED_EXT_MESH))
+ endif ! of if (.not. RECVS_CAN_BE_BURIED_EXT_MESH)
! end of loop on all the stations
enddo
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/locate_source.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/locate_source.f90 2009-07-10 16:39:24 UTC (rev 15452)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/locate_source.f90 2009-07-10 21:02:53 UTC (rev 15453)
@@ -167,14 +167,6 @@
! loop on all the sources
do isource = 1,NSOURCES
- if (.not. USE_EXTERNAL_MESH) then
-! check that the current source is inside the model
- if(lat(isource) < LATITUDE_MIN .or. lat(isource) > LATITUDE_MAX .or. long(isource) < LONGITUDE_MIN &
- .or. long(isource) > LONGITUDE_MAX) call exit_MPI(myrank,'the current source is outside the model')
-
- if(depth(isource) >= dabs(Z_DEPTH_BLOCK/1000.d0)) &
- call exit_MPI(myrank,'the current source is below the bottom of the model')
- endif
!
! r -> z, theta -> -y, phi -> x
!
@@ -194,7 +186,7 @@
Mxy(isource) = - moment_tensor(6,isource)
call utm_geo(long(isource),lat(isource),utm_x_source(isource),utm_y_source(isource), &
- UTM_PROJECTION_ZONE,ILONGLAT2UTM,(SUPPRESS_UTM_PROJECTION .or. USE_EXTERNAL_MESH))
+ UTM_PROJECTION_ZONE,ILONGLAT2UTM,.true.)
! orientation consistent with the UTM projection
@@ -213,62 +205,12 @@
nu_source(3,2,isource) = 0.d0
nu_source(3,3,isource) = 1.d0
- if (.not. USE_EXTERNAL_MESH) then
+ x_target_source = utm_x_source(isource)
+ y_target_source = utm_y_source(isource)
+ z_target_source = depth(isource)
+ if (myrank == 0) write(IOVTK,*) x_target_source, y_target_source, z_target_source
+
-! compute elevation of topography at the epicenter
- if(TOPOGRAPHY) then
-
-! get coordinate of corner in bathy/topo model
- icornerlong = int((long(isource) - ORIG_LONG_TOPO) / DEGREES_PER_CELL_TOPO) + 1
- icornerlat = int((lat(isource) - 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(isource) - long_corner) / DEGREES_PER_CELL_TOPO
- ratio_eta = (lat(isource) - 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(isource) = &
- itopo_bathy(icornerlong,icornerlat)*(1.-ratio_xi)*(1.-ratio_eta) + &
- itopo_bathy(icornerlong+1,icornerlat)*ratio_xi*(1.-ratio_eta) + &
- itopo_bathy(icornerlong+1,icornerlat+1)*ratio_xi*ratio_eta + &
- itopo_bathy(icornerlong,icornerlat+1)*(1.-ratio_xi)*ratio_eta
-
- else
- elevation = 0.d0
- endif
-
-! compute the Cartesian position of the source
-! take elevation of the surface into account
- x_target_source = utm_x_source(isource)
- y_target_source = utm_y_source(isource)
- z_target_source = - depth(isource)*1000.0d0 + elevation(isource)
- if(myrank == 0) write(IOVTK,*) x_target_source,y_target_source,z_target_source
-
- else
-
- x_target_source = utm_x_source(isource)
- y_target_source = utm_y_source(isource)
- z_target_source = depth(isource)
- if (myrank == 0) write(IOVTK,*) x_target_source, y_target_source, z_target_source
-
- endif ! of if (.not. USE_EXTERNAL_MESH)
-
! set distance to huge initial value
distmin = HUGEVAL
@@ -278,7 +220,7 @@
! define the interval in which we look for points
- if(USE_FORCE_POINT_SOURCE .and. USE_EXTERNAL_MESH) then
+ if(USE_FORCE_POINT_SOURCE) then
imin = 1
imax = NGLLX
@@ -307,7 +249,7 @@
iglob = ibool(i,j,k,ispec)
- if (USE_EXTERNAL_MESH .and. (.not. SOURCES_CAN_BE_BURIED_EXT_MESH)) then
+ if (.not. SOURCES_CAN_BE_BURIED_EXT_MESH) then
if ((.not. iglob_is_surface_external_mesh(iglob)) .or. (.not. ispec_is_surface_external_mesh(ispec))) then
cycle
endif
@@ -350,7 +292,7 @@
endif
! get normal to the face of the hexaedra if receiver is on the surface
- if (USE_EXTERNAL_MESH .and. (.not. SOURCES_CAN_BE_BURIED_EXT_MESH) .and. &
+ if ((.not. SOURCES_CAN_BE_BURIED_EXT_MESH) .and. &
.not. (ispec_selected_source(isource) == 0)) then
pt0_ix = -1
pt0_iy = -1
@@ -484,7 +426,7 @@
nu_source(3,2,isource) = v_vector(3)
nu_source(3,3,isource) = w_vector(3)
- endif ! of if (USE_EXTERNAL_MESH .and. (.not. RECEIVERS_CAN_BE_BURIED_EXT_MESH))
+ endif ! of if (.not. RECEIVERS_CAN_BE_BURIED_EXT_MESH)
! *******************************************
! find the best (xi,eta,gamma) for the source
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/save_header_file.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/save_header_file.f90 2009-07-10 16:39:24 UTC (rev 15452)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/save_header_file.f90 2009-07-10 21:02:53 UTC (rev 15453)
@@ -123,7 +123,6 @@
! if attenuation is off, set dummy size of arrays to one
if(ATTENUATION) then
- stop 'ATTENUATION not supported yet in the CUBIT + SCOTCH version because of arrays of constant size defined'
write(IOUT,*) 'integer, parameter :: NSPEC_ATTENUATION = ', NSPEC_AB
write(IOUT,*) 'logical, parameter :: ATTENUATION_VAL = .true.'
else
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/specfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/specfem3D.f90 2009-07-10 16:39:24 UTC (rev 15452)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/specfem3D.f90 2009-07-10 21:02:53 UTC (rev 15453)
@@ -141,11 +141,14 @@
integer iattenuation
double precision scale_factor
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION,N_SLS) :: &
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: &
R_xx,R_yy,R_xy,R_xz,R_yz
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION) :: &
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: &
epsilondev_xx,epsilondev_yy,epsilondev_xy,epsilondev_xz,epsilondev_yz
+ integer :: NSPEC_ATTENUATION_AB
+ integer, dimension(:,:,:,:),allocatable :: iflag_attenuation_store
+
! ADJOINT
real(kind=CUSTOM_REAL), dimension(NUM_REGIONS_ATTENUATION,N_SLS) :: b_alphaval, b_betaval, b_gammaval
!! DK DK array not created yet for CUBIT
@@ -155,14 +158,26 @@
! b_epsilondev_yy,b_epsilondev_xy,b_epsilondev_xz,b_epsilondev_yz
! ADJOINT
-! integer NPOIN2DMAX_XY
-
! 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
+ integer :: NSPEC2DMAX_XMIN_XMAX_ext,NSPEC2DMAX_YMIN_YMAX_ext
+ integer, dimension(:), allocatable :: ibelm_xmin,ibelm_xmax
+ integer, dimension(:), allocatable :: ibelm_ymin,ibelm_ymax
+ integer, dimension(:), allocatable :: ibelm_bottom
+ integer, dimension(:), allocatable :: ibelm_top
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: jacobian2D_xmin,jacobian2D_xmax
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: jacobian2D_ymin,jacobian2D_ymax
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: jacobian2D_bottom
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: jacobian2D_top
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: normal_xmin,normal_xmax
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: normal_ymin,normal_ymax
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: normal_bottom
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: normal_top
+
!! DK DK array not created yet for CUBIT
! integer, dimension(NSPEC2D_TOP_VAL) :: ibelm_top
! real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NSPEC2D_TOP_VAL) :: normal_top
@@ -228,12 +243,14 @@
! real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_ADJOINT) :: b_displ, b_veloc, b_accel
! real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT) :: rho_kl, mu_kl, kappa_kl, &
! rhop_kl, beta_kl, alpha_kl
+! real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: absorb_xmin, absorb_xmax, &
+! absorb_ymin, absorb_ymax, absorb_zmin ! for absorbing b.c.
+! integer reclen_xmin, reclen_xmax, reclen_ymin, reclen_ymax, reclen_zmin
+
real(kind=CUSTOM_REAL) b_deltat, b_deltatover2, b_deltatsqover2
! ADJOINT
integer l
- real(kind=CUSTOM_REAL) vs_val,Q_mu
- integer iattenuation_sediments,int_Q_mu
! Moho kernel
! integer ispec2D_moho_top, ispec2D_moho_bot, k_top, k_bot, ispec_top, ispec_bot, iglob_top, iglob_bot
@@ -371,6 +388,8 @@
! Stacey conditions put back
integer nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax,ispec2D
real(kind=CUSTOM_REAL) nx,ny,nz
+ integer, dimension(:,:),allocatable :: nimin,nimax,nkmin_eta
+ integer, dimension(:,:),allocatable :: njmin,njmax,nkmin_xi
! to save movie frames
integer ipoin, nmovie_points, iloc, iorderi(NGNOD2D), iorderj(NGNOD2D)
@@ -491,11 +510,6 @@
! get the base pathname for output files
call get_value_string(OUTPUT_FILES, 'OUTPUT_FILES', 'OUTPUT_FILES')
-! check that we use an external mesh, because this new version of the solver
-! (SPECFEM3D_SESAME) is for external meshes only (created for instance using CUBIT
-! and decomposed with METIS or SCOTCH)
- if (.not. USE_EXTERNAL_MESH) stop 'SPECFEM3D_SESAME is for external meshes only'
-
! check that optimized routines from Deville et al. (2002) can be used
if(NGLLX /= 5 .or. NGLLY /= 5 .or. NGLLZ /= 5) &
stop 'optimized routines from Deville et al. (2002) such as mxm_m1_m2_5points can only be used if NGLL = 5'
@@ -509,6 +523,8 @@
open(unit=27,file=prname(1:len_trim(prname))//'external_mesh.bin',status='old',action='read',form='unformatted')
read(27) NSPEC_AB
read(27) NGLOB_AB
+ !pll
+ NSPEC_ATTENUATION_AB = NSPEC_AB
close(27)
! open main output file, only written to by process 0
@@ -594,6 +610,7 @@
allocate(displ(NDIM,NGLOB_AB))
allocate(veloc(NDIM,NGLOB_AB))
allocate(accel(NDIM,NGLOB_AB))
+ allocate(iflag_attenuation_store(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
! info about external mesh simulation
! nlegoff -- should be put in read_arrays_solver and read_arrays_buffer_solver for clarity
@@ -611,6 +628,23 @@
read(27) gammay
read(27) gammaz
read(27) jacobian
+
+ !pll
+ read(27) rho_vp
+ read(27) rho_vs
+ read(27) iflag_attenuation_store
+ read(27) NSPEC2DMAX_XMIN_XMAX_ext
+ read(27) NSPEC2DMAX_YMIN_YMAX_ext
+ allocate(nimin(2,NSPEC2DMAX_YMIN_YMAX_ext),nimax(2,NSPEC2DMAX_YMIN_YMAX_ext),nkmin_eta(2,NSPEC2DMAX_YMIN_YMAX_ext))
+ allocate(njmin(2,NSPEC2DMAX_XMIN_XMAX_ext),njmax(2,NSPEC2DMAX_XMIN_XMAX_ext),nkmin_xi(2,NSPEC2DMAX_XMIN_XMAX_ext))
+ read(27) nimin
+ read(27) nimax
+ read(27) njmin
+ read(27) njmax
+ read(27) nkmin_xi
+ read(27) nkmin_eta
+ !end pll
+
read(27) kappastore
read(27) mustore
read(27) rmass
@@ -619,6 +653,51 @@
read(27) ystore
read(27) zstore
+ !pll
+ read(27) nspec2D_xmin
+ read(27) nspec2D_xmax
+ read(27) nspec2D_ymin
+ read(27) nspec2D_ymax
+ read(27) NSPEC2D_BOTTOM
+ read(27) NSPEC2D_TOP
+ allocate(ibelm_xmin(nspec2D_xmin))
+ allocate(ibelm_xmax(nspec2D_xmax))
+ allocate(ibelm_ymin(nspec2D_ymin))
+ allocate(ibelm_ymax(nspec2D_ymax))
+ allocate(ibelm_bottom(NSPEC2D_BOTTOM))
+ allocate(ibelm_top(NSPEC2D_TOP))
+ allocate(jacobian2D_xmin(NGLLY,NGLLZ,nspec2D_xmin))
+ allocate(jacobian2D_xmax(NGLLY,NGLLZ,nspec2D_xmax))
+ allocate(jacobian2D_ymin(NGLLX,NGLLZ,nspec2D_ymin))
+ allocate(jacobian2D_ymax(NGLLX,NGLLZ,nspec2D_ymax))
+ allocate(jacobian2D_bottom(NGLLX,NGLLY,NSPEC2D_BOTTOM))
+ allocate(jacobian2D_top(NGLLX,NGLLY,NSPEC2D_TOP))
+ allocate(normal_xmin(NDIM,NGLLY,NGLLZ,nspec2D_xmin))
+ allocate(normal_xmax(NDIM,NGLLY,NGLLZ,nspec2D_xmax))
+ allocate(normal_ymin(NDIM,NGLLX,NGLLZ,nspec2D_ymin))
+ allocate(normal_ymax(NDIM,NGLLX,NGLLZ,nspec2D_ymax))
+ allocate(normal_bottom(NDIM,NGLLX,NGLLY,NSPEC2D_BOTTOM))
+ allocate(normal_top(NDIM,NGLLX,NGLLY,NSPEC2D_TOP))
+ read(27) ibelm_xmin
+ read(27) ibelm_xmax
+ read(27) ibelm_ymin
+ read(27) ibelm_ymax
+ read(27) ibelm_bottom
+ read(27) ibelm_top
+ read(27) normal_xmin
+ read(27) normal_xmax
+ read(27) normal_ymin
+ read(27) normal_ymax
+ read(27) normal_bottom
+ read(27) normal_top
+ read(27) jacobian2D_xmin
+ read(27) jacobian2D_xmax
+ read(27) jacobian2D_ymin
+ read(27) jacobian2D_ymax
+ read(27) jacobian2D_bottom
+ read(27) jacobian2D_top
+ !end pll
+
read(27) ninterfaces_ext_mesh
read(27) max_nibool_interfaces_ext_mesh
allocate(my_neighbours_ext_mesh(ninterfaces_ext_mesh))
@@ -1162,17 +1241,16 @@
!$$$$$$$$$$$$$$$$$$ RECEIVERS $$$$$$$$$$$$$$$$$$$$$
if (SIMULATION_TYPE == 1) then
- call get_value_string(filtered_rec_filename, 'solver.STATIONS_FILTERED', 'DATA/STATIONS_FILTERED')
+ call get_value_string(rec_filename, 'solver.STATIONS', 'DATA/STATIONS')
! get total number of stations
- open(unit=IIN,file=filtered_rec_filename,iostat=ios,status='old',action='read')
+ 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')
else
@@ -1208,7 +1286,7 @@
! locate receivers in the mesh
call locate_receivers(ibool,myrank,NSPEC_AB,NGLOB_AB, &
- xstore,ystore,zstore,xigll,yigll,zigll,filtered_rec_filename, &
+ xstore,ystore,zstore,xigll,yigll,zigll,rec_filename, &
nrec,islice_selected_rec,ispec_selected_rec, &
xi_receiver,eta_receiver,gamma_receiver,station_name,network_name,nu, &
NPROC,utm_x_source(1),utm_y_source(1), &
@@ -1457,74 +1535,87 @@
! rescale shear modulus according to attenuation model
+!pll
+! do ispec = 1,NSPEC_AB
+! if(not_fully_in_bedrock(ispec)) then
+! do k=1,NGLLZ
+! do j=1,NGLLY
+! do i=1,NGLLX
+!
+!! distinguish attenuation factors
+! if(flag_sediments(i,j,k,ispec)) then
+!
+!! use constant attenuation of Q = 90
+!! or use scaling rule similar to Olsen et al. (2003)
+!! We might need to fix the attenuation part for the anisotropy case
+!! At this stage, we turn the ATTENUATION flag off always, and still keep mustore
+! if(USE_OLSEN_ATTENUATION) then
+! vs_val = mustore(i,j,k,ispec) / rho_vs(i,j,k,ispec)
+!! use rule Q_mu = constant * v_s
+! Q_mu = OLSEN_ATTENUATION_RATIO * vs_val
+! int_Q_mu = 10 * nint(Q_mu / 10.)
+! if(int_Q_mu < 40) int_Q_mu = 40
+! if(int_Q_mu > 150) int_Q_mu = 150
+!
+! if(int_Q_mu == 40) then
+! iattenuation_sediments = IATTENUATION_SEDIMENTS_40
+! else if(int_Q_mu == 50) then
+! iattenuation_sediments = IATTENUATION_SEDIMENTS_50
+! else if(int_Q_mu == 60) then
+! iattenuation_sediments = IATTENUATION_SEDIMENTS_60
+! else if(int_Q_mu == 70) then
+! iattenuation_sediments = IATTENUATION_SEDIMENTS_70
+! else if(int_Q_mu == 80) then
+! iattenuation_sediments = IATTENUATION_SEDIMENTS_80
+! else if(int_Q_mu == 90) then
+! iattenuation_sediments = IATTENUATION_SEDIMENTS_90
+! else if(int_Q_mu == 100) then
+! iattenuation_sediments = IATTENUATION_SEDIMENTS_100
+! else if(int_Q_mu == 110) then
+! iattenuation_sediments = IATTENUATION_SEDIMENTS_110
+! else if(int_Q_mu == 120) then
+! iattenuation_sediments = IATTENUATION_SEDIMENTS_120
+! else if(int_Q_mu == 130) then
+! iattenuation_sediments = IATTENUATION_SEDIMENTS_130
+! else if(int_Q_mu == 140) then
+! iattenuation_sediments = IATTENUATION_SEDIMENTS_140
+! else if(int_Q_mu == 150) then
+! iattenuation_sediments = IATTENUATION_SEDIMENTS_150
+! else
+! stop 'incorrect attenuation coefficient'
+! endif
+!
+! else
+! iattenuation_sediments = IATTENUATION_SEDIMENTS_90
+! endif
+!
+! scale_factor = factor_scale(iattenuation_sediments)
+! else
+! scale_factor = factor_scale(IATTENUATION_BEDROCK)
+! endif
+!
+! mustore(i,j,k,ispec) = mustore(i,j,k,ispec) * scale_factor
+!
+! enddo
+! enddo
+! enddo
+! endif
+! enddo
+
+ !pll
do ispec = 1,NSPEC_AB
- if(not_fully_in_bedrock(ispec)) then
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
-
-! distinguish attenuation factors
- if(flag_sediments(i,j,k,ispec)) then
-
-! use constant attenuation of Q = 90
-! or use scaling rule similar to Olsen et al. (2003)
-! We might need to fix the attenuation part for the anisotropy case
-! At this stage, we turn the ATTENUATION flag off always, and still keep mustore
- if(USE_OLSEN_ATTENUATION) then
- vs_val = mustore(i,j,k,ispec) / rho_vs(i,j,k,ispec)
-! use rule Q_mu = constant * v_s
- Q_mu = OLSEN_ATTENUATION_RATIO * vs_val
- int_Q_mu = 10 * nint(Q_mu / 10.)
- if(int_Q_mu < 40) int_Q_mu = 40
- if(int_Q_mu > 150) int_Q_mu = 150
-
- if(int_Q_mu == 40) then
- iattenuation_sediments = IATTENUATION_SEDIMENTS_40
- else if(int_Q_mu == 50) then
- iattenuation_sediments = IATTENUATION_SEDIMENTS_50
- else if(int_Q_mu == 60) then
- iattenuation_sediments = IATTENUATION_SEDIMENTS_60
- else if(int_Q_mu == 70) then
- iattenuation_sediments = IATTENUATION_SEDIMENTS_70
- else if(int_Q_mu == 80) then
- iattenuation_sediments = IATTENUATION_SEDIMENTS_80
- else if(int_Q_mu == 90) then
- iattenuation_sediments = IATTENUATION_SEDIMENTS_90
- else if(int_Q_mu == 100) then
- iattenuation_sediments = IATTENUATION_SEDIMENTS_100
- else if(int_Q_mu == 110) then
- iattenuation_sediments = IATTENUATION_SEDIMENTS_110
- else if(int_Q_mu == 120) then
- iattenuation_sediments = IATTENUATION_SEDIMENTS_120
- else if(int_Q_mu == 130) then
- iattenuation_sediments = IATTENUATION_SEDIMENTS_130
- else if(int_Q_mu == 140) then
- iattenuation_sediments = IATTENUATION_SEDIMENTS_140
- else if(int_Q_mu == 150) then
- iattenuation_sediments = IATTENUATION_SEDIMENTS_150
- else
- stop 'incorrect attenuation coefficient'
- endif
-
- else
- iattenuation_sediments = IATTENUATION_SEDIMENTS_90
- endif
-
- scale_factor = factor_scale(iattenuation_sediments)
- else
- scale_factor = factor_scale(IATTENUATION_BEDROCK)
- endif
-
- mustore(i,j,k,ispec) = mustore(i,j,k,ispec) * scale_factor
-
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ scale_factor = factor_scale(iflag_attenuation_store(i,j,k,ispec))
+ mustore(i,j,k,ispec) = mustore(i,j,k,ispec) * scale_factor
+ enddo
enddo
- enddo
- enddo
- endif
+ enddo
enddo
+
+ endif
- endif
-
! allocate seismogram array
if (nrec_local > 0) then
allocate(seismograms_d(NDIM,nrec_local,NSTEP))
@@ -1659,9 +1750,22 @@
endif
endif
+
+ !pll, to put elsewhere
+ allocate(R_xx(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS))
+ allocate(R_yy(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS))
+ allocate(R_xy(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS))
+ allocate(R_xz(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS))
+ allocate(R_yz(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS))
+ allocate(epsilondev_xx(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB))
+ allocate(epsilondev_yy(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB))
+ allocate(epsilondev_xy(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB))
+ allocate(epsilondev_xz(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB))
+ allocate(epsilondev_yz(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB))
+
! clear memory variables if attenuation
if(ATTENUATION) then
-
+
! initialize memory variables for attenuation
epsilondev_xx(:,:,:,:) = 0._CUSTOM_REAL
epsilondev_yy(:,:,:,:) = 0._CUSTOM_REAL
@@ -1839,11 +1943,22 @@
! endif
! assemble all the contributions between slices using MPI
+
+
if(USE_DEVILLE_PRODUCTS) then
- call compute_forces_with_Deville(NSPEC_AB,NGLOB_AB,displ,accel,xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ call compute_forces_with_Deville(NSPEC_AB,NGLOB_AB,ATTENUATION,displ,accel,xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT,wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
kappastore,mustore,jacobian,ibool,ispec_is_inner_ext_mesh,.false., &
- NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,xi_source,eta_source,gamma_source,nu_source,hdur,dt)
+ NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,xi_source,eta_source,gamma_source,nu_source, &
+ hdur,hdur_gaussian,t_cmt,dt,stf,t0,sourcearrays, &
+ one_minus_sum_beta,factor_common,alphaval,betaval,gammaval,R_xx,R_yy,R_xy,R_xz,R_yz, &
+ epsilondev_xx,epsilondev_yy,epsilondev_xy,epsilondev_xz,epsilondev_yz,iflag_attenuation_store, &
+ ABSORBING_CONDITIONS,SAVE_FORWARD,NSTEP,SIMULATION_TYPE, &
+ nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax,NSPEC2D_BOTTOM,NSPEC2DMAX_XMIN_XMAX_ext,NSPEC2DMAX_YMIN_YMAX_ext, &
+ ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom, &
+ nimin,nimax,njmin,njmax,nkmin_xi,nkmin_eta, &
+ veloc,rho_vp,rho_vs,jacobian2D_xmin,jacobian2D_xmax,jacobian2D_ymin,jacobian2D_ymax,jacobian2D_bottom, &
+ normal_xmin,normal_xmax,normal_ymin,normal_ymax,normal_bottom)
else
call compute_forces_no_Deville(NSPEC_AB,NGLOB_AB,displ,accel,xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
hprime_xx,hprime_yy,hprime_zz,hprimewgll_xx,hprimewgll_yy,hprimewgll_zz,wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
@@ -1858,10 +1973,19 @@
request_send_vector_ext_mesh,request_recv_vector_ext_mesh)
if(USE_DEVILLE_PRODUCTS) then
- call compute_forces_with_Deville(NSPEC_AB,NGLOB_AB,displ,accel,xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ call compute_forces_with_Deville(NSPEC_AB,NGLOB_AB,ATTENUATION,displ,accel,xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT,wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
kappastore,mustore,jacobian,ibool,ispec_is_inner_ext_mesh,.true., &
- NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,xi_source,eta_source,gamma_source,nu_source,hdur,dt)
+ NSOURCES,myrank,it,islice_selected_source,ispec_selected_source,xi_source,eta_source,gamma_source,nu_source, &
+ hdur,hdur_gaussian,t_cmt,dt,stf,t0,sourcearrays, &
+ one_minus_sum_beta,factor_common,alphaval,betaval,gammaval,R_xx,R_yy,R_xy,R_xz,R_yz, &
+ epsilondev_xx,epsilondev_yy,epsilondev_xy,epsilondev_xz,epsilondev_yz,iflag_attenuation_store, &
+ ABSORBING_CONDITIONS,SAVE_FORWARD,NSTEP,SIMULATION_TYPE, &
+ nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax,NSPEC2D_BOTTOM,NSPEC2DMAX_XMIN_XMAX_ext,NSPEC2DMAX_YMIN_YMAX_ext, &
+ ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom, &
+ nimin,nimax,njmin,njmax,nkmin_xi,nkmin_eta, &
+ veloc,rho_vp,rho_vs,jacobian2D_xmin,jacobian2D_xmax,jacobian2D_ymin,jacobian2D_ymax,jacobian2D_bottom, &
+ normal_xmin,normal_xmax,normal_ymin,normal_ymax,normal_bottom)
else
call compute_forces_no_Deville(NSPEC_AB,NGLOB_AB,displ,accel,xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
hprime_xx,hprime_yy,hprime_zz,hprimewgll_xx,hprimewgll_yy,hprimewgll_zz,wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/todo_list_please_dont_remove.txt
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/todo_list_please_dont_remove.txt 2009-07-10 16:39:24 UTC (rev 15452)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/todo_list_please_dont_remove.txt 2009-07-10 21:02:53 UTC (rev 15453)
@@ -43,6 +43,40 @@
- etudier comment lire un maillage CUBIT stocke au format natif de CUBIT (Exodus, base sur netCDF) depuis SPECFEM3D. On pourrait utiliser la commande "ncdump" dans CUBIT si necessaire d'après Emanuele Casarotti. Une autre option serait d'utiliser le format de stockage ABAQUS dans CUBIT.
+
+by Nicolas Le Goff :
+--------------------
+
+DESIGN :
+
+ - The subdivision of the mesh should be done inside meshfem3D (there is actually no point in doing it prior to pre_meshfem). The subdivision itself is already implemented, but we also have to subdivide the communications and any other similar data that concerns elements (absorbing boundaries, PML,...).
+
+ - Nearly all constants declared in constants.h concerning external meshes (look for string 'nlegoff') should be read/computed from a kind of Par_file. The same goes for the declaration of materials, done with a function in meshfem3D.f90 (and also right now in specfem3D.f90 in case of associating materials after pre_meshfem but we should get rid of it). There are also some variables declared as parameters used in create_movie_AVS_DX.f90.
+
+ - Materials should be associated to elements before pre_meshfem, but it is not always possible (that was the case for Celine Blitz's asteroid mesh with regolith because of limitations in the CUBIT software with huge meshes). An example of how we can isolate elements in contact with the envelope and associate them with another material is commented in specfem3D.f90 (look for string 'NL NL REGOLITH') tough it should in fact be inside meshfem3D, not specfem3D. Identifying the envelope is based on the NGLL points, so there is no problem with tranfering it to meshfem3D.
+
+ - Calculations on elements in specfem3D for regular meshes should be merged inside compute_forces (the basic calculation is okay, but not attenuation, absorbing boundaries, ...) as there is no reason to differenciate those.
+
+MISC :
+
+ - The following default options in flags.guess -qlanglvl=2003pure (for xlf) and -std=f2003 (for gfortran) causes some problems with intrinsics. There might be a better way than to simply remove these options.
+
+ - Reading the file DATA/STATIONS_FILTERED with pgf90 returned an error, because it considered that the lines were truncated, thus failed to read some data even when DATA/STATIONS_FILTERED was also generated by pgf90. Go figures. There might be a better workaround than simply removing writing of blanks ' ' (see rev).
+
+Add-ONS :
+
+ - Receivers distance to source along the surface computed for an asteroid can be done by computing the shortest path in a graph. No problem in serial, doing so in parallel might be tricky.
+
+ - Partitioning should be done in meshfem3D, using a parallel partitioner such as ParMETIS and/or SCOTCH (PT-SCOTCH). Reading the (distributed or not) mesh, building the graph, distributing, and partitioning is not a problem, but redistributing the graph/mesh according to the partitioning obtained might well be.
+
+ - Adding absorbing boundaries, attenuation, PML, moment sources, ... for external meshes.
+
+ - Differentiate between elements in contact with the envelope, and those that are in contact with a fracture : source, receivers, ... anything based on the envelope detection will have to take into account those fractures. Really tricky without prior information on the mesh.
+
+ - Settle the issue of normal recording for the receivers. The normal is unique, but the plane can be described using two vectors at random (while conserving an orthonormal reference), so we have no reason to choose a pair of vectors compared to the others.i
+
+
+
---------------------------------------------------------------------------------
---------------------------------------------------------------------------------
---------------------------------------------------------------------------------
@@ -53,18 +87,19 @@
- Here is how to run the code:
-1/ creer un maillage CUBIT
+1/ create a CUBIT mesh
+Possibly decimate it with "decimate_mesh"
-2/ controler sa qualite avec mon code Fortran
+2/ check his quality with Dimitri's Fortran code
"check_mesh_quality_CUBIT_Abaqus"
-2/ le decomposer avec "decompose_mesh_SCOTCH"
+3/ export it with Emanuele's scripts, in decompose_mesh_SCOTCH/OUTPUT_FILES
-3/ eventuellement le decimer avec "decimate_mesh"
+2/ decompose it with "decompose_mesh_SCOTCH" and move the output in your "LOCAL_PATH" (see DATA/Par_file)
-4/ creer les bases de donnees et les matrices avec "generate_databases.f90"
+4/ create databases and matrices with "generate_databases"
-5/ tourner le solveur "specfem3D.f90"
+5/ compile and run the solver "specfem3D"
- remove old things from SPECFEM3D_BASIN in the manual and add new things
related to CUBIT, METIS, SCOTCH etc. Add new figures showing the asteroid
@@ -78,3 +113,6 @@
Done:
-----
+ - Added attenuation
+
+ - Added moment sources
More information about the CIG-COMMITS
mailing list