[cig-commits] r17215 - in seismo/3D/SPECFEM3D/trunk: . obsolete_routines_from_SPECFEM3D_v1.4.4
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Fri Sep 24 15:17:01 PDT 2010
Author: dkomati1
Date: 2010-09-24 15:17:01 -0700 (Fri, 24 Sep 2010)
New Revision: 17215
Added:
seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/compute_rho_estimate.f90
seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/define_subregions.f90
seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/define_subregions_heuristic.f90
seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/get_absorb.f90
seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/get_flags_boundaries.f90
seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/hauksson_model.f90
seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/interpolate_gocad_block_HR.f90
seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/interpolate_gocad_block_MR.f90
seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/mesh_vertical.f90
seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/read_arrays_buffers_solver.f90
seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/read_moho_map.f90
seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/salton_trough_gocad.f90
seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/socal_model.f90
Removed:
seismo/3D/SPECFEM3D/trunk/compute_rho_estimate.f90
seismo/3D/SPECFEM3D/trunk/define_subregions.f90
seismo/3D/SPECFEM3D/trunk/define_subregions_heuristic.f90
seismo/3D/SPECFEM3D/trunk/get_absorb.f90
seismo/3D/SPECFEM3D/trunk/get_flags_boundaries.f90
seismo/3D/SPECFEM3D/trunk/hauksson_model.f90
seismo/3D/SPECFEM3D/trunk/interpolate_gocad_block_HR.f90
seismo/3D/SPECFEM3D/trunk/interpolate_gocad_block_MR.f90
seismo/3D/SPECFEM3D/trunk/mesh_vertical.f90
seismo/3D/SPECFEM3D/trunk/read_arrays_buffers_solver.f90
seismo/3D/SPECFEM3D/trunk/read_moho_map.f90
seismo/3D/SPECFEM3D/trunk/salton_trough_gocad.f90
seismo/3D/SPECFEM3D/trunk/socal_model.f90
Log:
moved obsolete routines to directory "obsolete_routines_from_SPECFEM3D_v1.4.4"
Deleted: seismo/3D/SPECFEM3D/trunk/compute_rho_estimate.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/compute_rho_estimate.f90 2010-09-24 22:15:18 UTC (rev 17214)
+++ seismo/3D/SPECFEM3D/trunk/compute_rho_estimate.f90 2010-09-24 22:17:01 UTC (rev 17215)
@@ -1,46 +0,0 @@
-!=====================================================================
-!
-! S p e c f e m 3 D V e r s i o n 1 . 4
-! ---------------------------------------
-!
-! Dimitri Komatitsch and Jeroen Tromp
-! Seismological Laboratory - California Institute of Technology
-! (c) California Institute of Technology September 2006
-!
-! This program is free software; you can redistribute it and/or modify
-! it under the terms of the GNU General Public License as published by
-! the Free Software Foundation; either version 2 of the License, or
-! (at your option) any later version.
-!
-! This program is distributed in the hope that it will be useful,
-! but WITHOUT ANY WARRANTY; without even the implied warranty of
-! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-! GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along
-! with this program; if not, write to the Free Software Foundation, Inc.,
-! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-!
-!=====================================================================
-
- subroutine compute_rho_estimate(rho,vp)
-
-! compute rho estimate in Gocad block and in Hauksson's model
-! based upon Vp
-
- implicit none
-
-! include "constants.h"
- include "constants_gocad.h"
-
- double precision rho,vp
-
-! scale density - use empirical rule from Christiane
- rho = 0.33d0 * vp + 1280.d0
-
-! make sure density estimate is reasonable
- if(rho > DENSITY_MAX) rho = DENSITY_MAX
- if(rho < DENSITY_MIN) rho = DENSITY_MIN
-
- end subroutine compute_rho_estimate
-
Deleted: seismo/3D/SPECFEM3D/trunk/define_subregions.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/define_subregions.f90 2010-09-24 22:15:18 UTC (rev 17214)
+++ seismo/3D/SPECFEM3D/trunk/define_subregions.f90 2010-09-24 22:17:01 UTC (rev 17215)
@@ -1,863 +0,0 @@
-!=====================================================================
-!
-! S p e c f e m 3 D V e r s i o n 1 . 4
-! ---------------------------------------
-!
-! Dimitri Komatitsch and Jeroen Tromp
-! Seismological Laboratory - California Institute of Technology
-! (c) California Institute of Technology September 2006
-!
-! This program is free software; you can redistribute it and/or modify
-! it under the terms of the GNU General Public License as published by
-! the Free Software Foundation; either version 2 of the License, or
-! (at your option) any later version.
-!
-! This program is distributed in the hope that it will be useful,
-! but WITHOUT ANY WARRANTY; without even the implied warranty of
-! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-! GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along
-! with this program; if not, write to the Free Software Foundation, Inc.,
-! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-!
-!=====================================================================
-
- subroutine 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)
-
-! define shape of elements in current subregion of the mesh
-
- implicit none
-
- include "constants.h"
-
- integer myrank
- integer ix1,ix2,dix,iy1,iy2,diy,ir1,ir2,dir
- integer iax,iay,iar
- integer isubregion,doubling_index
- integer npx,npy
-
- integer NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM,NER_SEDIM,NER
-
- logical USE_REGULAR_MESH
-
-! topology of the elements
- integer iaddx(NGNOD)
- integer iaddy(NGNOD)
- integer iaddz(NGNOD)
-
-! **************
-
-!
-!--- case of a regular mesh
-!
- if(USE_REGULAR_MESH) then
-
-! use two layers even for a regular mesh, because the algorithm detects the top of the mesh
-! (the "topography") based on one layer of elements with flag IFLAG_ONE_LAYER_TOPOGRAPHY
- if(isubregion == 2) then
-
- call usual_hex_nodes(iaddx,iaddy,iaddz)
-
- iy1=0
- iy2=npy-2
- diy=2
-
- ix1=0
- ix2=npx-2
- dix=2
-
- ir1=0
- ir2=2*(NER - 2)
- dir=2
-
- iax=1
- iay=1
- iar=1
-
- doubling_index = IFLAG_BASEMENT_TOPO
-
- else if(isubregion == 1) then
-
- call usual_hex_nodes(iaddx,iaddy,iaddz)
-
- iy1=0
- iy2=npy-2
- diy=2
-
- ix1=0
- ix2=npx-2
- dix=2
-
- ir1=2*(NER - 1)
- ir2=ir1
- dir=2
-
- iax=1
- iay=1
- iar=1
-
- doubling_index = IFLAG_ONE_LAYER_TOPOGRAPHY
-
- else
-
- call exit_MPI(myrank,'incorrect subregion code')
-
- endif
-
-!
-!--- case of a non-regular mesh with mesh doublings
-!
- else
-
-! this last region only defined when NER_SEDIM > 1
- if(isubregion == 30) then
-
- call usual_hex_nodes(iaddx,iaddy,iaddz)
-
- iy1=0
- iy2=npy-2
- diy=2
-
- ix1=0
- ix2=npx-2
- dix=2
-
- ir1=2*(NER - NER_SEDIM)
- ir2=2*(NER - 2)
- dir=2
-
- iax=1
- iay=1
- iar=1
-
- doubling_index = IFLAG_BASEMENT_TOPO
-
- else if(isubregion == 29) then
-
- call usual_hex_nodes(iaddx,iaddy,iaddz)
-
- iy1=0
- iy2=npy-2
- diy=2
-
- ix1=0
- ix2=npx-2
- dix=2
-
- ir1=2*(NER - 1)
- ir2=ir1
- dir=2
-
- iax=1
- iay=1
- iar=1
-
- doubling_index = IFLAG_ONE_LAYER_TOPOGRAPHY
-
- else if(isubregion == 28) then
-
- call usual_hex_nodes(iaddx,iaddy,iaddz)
-
- iy1=0
- iy2=npy-8
- diy=8
-
- ix1=0
- ix2=npx-8
- dix=8
-
- ir1= 0
- ir2= 2*NER_BOTTOM_MOHO-8
- dir=8
-
- iax=4
- iay=4
- iar=4
-
- doubling_index= IFLAG_HALFSPACE_MOHO
-
- else if(isubregion == 27) then
-
- call unusual_hex_nodes1(iaddx,iaddy,iaddz)
-
-! generating stage 2 of the mesh doubling below 670
-
- iy1=0
- iy2=npy-8
- diy=8
-
- ix1=0
- ix2=npx-16
- dix=16
-
- dir=4
-
- iax=2
- iay=2
- iar=1
-
- ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 8
- ir2=ir1
-
- doubling_index=IFLAG_16km_BASEMENT
-
- else if(isubregion == 26) then
-
- call unusual_hex_nodes1p(iaddx,iaddy,iaddz)
-
-! generating stage 3 of the mesh doubling below 670
-
- iy1=0
- iy2=npy-8
- diy=8
-
- ix1=8
- ix2=npx-8
- dix=16
-
- dir=4
-
- iax=2
- iay=2
- iar=1
-
- ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 8
- ir2=ir1
-
- doubling_index=IFLAG_16km_BASEMENT
-
- else if(isubregion == 25) then
-
- call unusual_hex_nodes2(iaddx,iaddy,iaddz)
-
-! generating stage 4 of the mesh doubling below 670
-
- iy1=0
- iy2=npy-8
- diy=8
-
- ix1=0
- ix2=npx-16
- dix=16
-
- dir=4
-
- iax=2
- iay=2
- iar=1
-
- ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 8
- ir2=ir1
-
- doubling_index=IFLAG_16km_BASEMENT
-
- else if(isubregion == 24) then
-
- call unusual_hex_nodes2p(iaddx,iaddy,iaddz)
-
-! generating stage 5 of the mesh doubling below 670
-
- iy1=0
- iy2=npy-8
- diy=8
-
- ix1=12
- ix2=npx-4
- dix=16
-
- dir=4
-
- iax=2
- iay=2
- iar=1
-
- ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 6
- ir2=ir1
-
- doubling_index=IFLAG_16km_BASEMENT
-
- else if(isubregion == 23) then
-
- call unusual_hex_nodes3(iaddx,iaddy,iaddz)
-
-! generating stage 6 of the mesh doubling below 670
-
- iy1=0
- iy2=npy-8
- diy=8
-
- ix1=4
- ix2=npx-12
- dix=16
-
- dir=4
-
- iax=2
- iay=2
- iar=1
-
- ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 6
- ir2=ir1
-
- doubling_index=IFLAG_16km_BASEMENT
-
- else if(isubregion == 22) then
-
- call unusual_hex_nodes3(iaddx,iaddy,iaddz)
-
-! generating stage 7 of the mesh doubling below 670
-
- iy1=0
- iy2=npy-8
- diy=8
-
- ix1=8
- ix2=npx-8
- dix=16
-
- dir=4
-
- iax=2
- iay=2
- iar=1
-
- ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 6
- ir2=ir1
-
- doubling_index=IFLAG_16km_BASEMENT
-
- else if(isubregion == 21) then
-
- call unusual_hex_nodes4(iaddx,iaddy,iaddz)
-
-! generating stage 8 of the mesh doubling below 670
-
- iy1=8
- iy2=npy-8
- diy=16
-
- ix1=0
- ix2=npx-4
- dix=4
-
- dir=4
-
- iax=2
- iay=2
- iar=1
-
- ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 4
- ir2=ir1
-
- doubling_index=IFLAG_16km_BASEMENT
-
- else if(isubregion == 20) then
-
- call unusual_hex_nodes4p(iaddx,iaddy,iaddz)
-
-! generating stage 9 of the mesh doubling below 670
-
- iy1=0
- iy2=npy-16
- diy=16
-
- ix1=0
- ix2=npx-4
- dix=4
-
- dir=4
-
- iax=2
- iay=2
- iar=1
-
- ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 4
- ir2=ir1
-
- doubling_index=IFLAG_16km_BASEMENT
-
- else if(isubregion == 19) then
-
- call usual_hex_nodes(iaddx,iaddy,iaddz)
-
-! generating stage 10 of the mesh doubling below 670
-
- iy1=8
- iy2=npy-8
- diy=16
-
- ix1=0
- ix2=npx-4
- dix=4
-
- dir=4
-
- iax=2
- iay=2
- iar=1
-
- ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 2
- ir2=ir1
-
- doubling_index=IFLAG_16km_BASEMENT
-
- else if(isubregion == 18) then
-
- call usual_hex_nodes(iaddx,iaddy,iaddz)
-
-! generating stage 11 of the mesh doubling below 670
-
- iy1=4
- iy2=npy-12
- diy=16
-
- ix1=0
- ix2=npx-4
- dix=4
-
- dir=4
-
- iax=2
- iay=2
- iar=1
-
- ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 2
- ir2=ir1
-
- doubling_index=IFLAG_16km_BASEMENT
-
- else if(isubregion == 17) then
-
- call unusual_hex_nodes6(iaddx,iaddy,iaddz)
-
-! generating stage 12 of the mesh doubling below 670
-
- iy1=12
- iy2=npy-4
- diy=16
-
- ix1=0
- ix2=npx-4
- dix=4
-
- dir=4
-
- iax=2
- iay=2
- iar=1
-
- ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 2
- ir2=ir1
-
- doubling_index=IFLAG_16km_BASEMENT
-
- else if(isubregion == 16) then
-
- call unusual_hex_nodes6p(iaddx,iaddy,iaddz)
-
-! generating stage 13 of the mesh doubling below 670
-
- iy1=0
- iy2=npy-16
- diy=16
-
- ix1=0
- ix2=npx-4
- dix=4
-
- dir=4
-
- iax=2
- iay=2
- iar=1
-
- ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 4
- ir2=ir1
-
- doubling_index=IFLAG_16km_BASEMENT
-
- else if(isubregion == 15) then
-
- call usual_hex_nodes(iaddx,iaddy,iaddz)
-
- iy1=0
- iy2=npy-8
- diy=8
-
- ix1=0
- ix2=npx-8
- dix=8
-
-! honor So-Cal model discontinuity at 16 km for accuracy
- ir1=2*NER_BOTTOM_MOHO
- ir2=2*(NER_BOTTOM_MOHO+NER_MOHO_16) - 4
- dir=4
-
- iax=4
- iay=4
- iar=2
-
- doubling_index = IFLAG_MOHO_16km
-
- else if(isubregion == 14) then
-
- call usual_hex_nodes(iaddx,iaddy,iaddz)
-
- iy1=0
- iy2=npy-8
- diy=8
-
- ix1=0
- ix2=npx-8
- dix=8
-
-! honor So-Cal model discontinuity at 16 km for accuracy
- ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16)
- ir2=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT)-12
- dir=4
-
- iax=4
- iay=4
- iar=2
-
- doubling_index = IFLAG_16km_BASEMENT
-
-
- else if(isubregion == 13) then
-
- call usual_hex_nodes(iaddx,iaddy,iaddz)
-
-! generating stage 1 of the mesh doubling below the Moho
-
- iy1=0
- iy2=npy-4
- diy=4
-
- ix1=0
- ix2=npx-4
- dix=4
-
- ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT)
- ir2=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-12
- dir=4
-
- iax=2
- iay=2
- iar=2
-
- doubling_index=IFLAG_BASEMENT_TOPO
-
- else if(isubregion == 12) then
-
- call unusual_hex_nodes1(iaddx,iaddy,iaddz)
-
-! generating stage 2 of the mesh doubling below the Moho
-
- iy1=0
- iy2=npy-4
- diy=4
-
- ix1=0
- ix2=npx-8
- dix=8
-
- dir=4
-
- iax=1
- iay=1
- iar=1
-
- ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-8
- ir2=ir1
-
- doubling_index=IFLAG_BASEMENT_TOPO
-
- else if(isubregion == 11) then
-
- call unusual_hex_nodes1p(iaddx,iaddy,iaddz)
-
-! generating stage 3 of the mesh doubling below the Moho
-
- iy1=0
- iy2=npy-4
- diy=4
-
- ix1=4
- ix2=npx-4
- dix=8
-
- dir=4
-
- iax=1
- iay=1
- iar=1
-
- ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-8
- ir2=ir1
-
- doubling_index=IFLAG_BASEMENT_TOPO
-
- else if(isubregion == 10) then
-
- call unusual_hex_nodes2(iaddx,iaddy,iaddz)
-
-! generating stage 4 of the mesh doubling below the Moho
-
- iy1=0
- iy2=npy-4
- diy=4
-
- ix1=0
- ix2=npx-8
- dix=8
-
- dir=4
-
- iax=1
- iay=1
- iar=1
-
- ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-8
- ir2=ir1
-
- doubling_index=IFLAG_BASEMENT_TOPO
-
- else if(isubregion == 9) then
-
- call unusual_hex_nodes2p(iaddx,iaddy,iaddz)
-
-! generating stage 5 of the mesh doubling below the Moho
-
- iy1=0
- iy2=npy-4
- diy=4
-
- ix1=6
- ix2=npx-2
- dix=8
-
- dir=4
-
- iax=1
- iay=1
- iar=1
-
- ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-6
- ir2=ir1
-
- doubling_index=IFLAG_BASEMENT_TOPO
-
- else if(isubregion == 8) then
-
- call unusual_hex_nodes3(iaddx,iaddy,iaddz)
-
-! generating stage 6 of the mesh doubling below the Moho
-
- iy1=0
- iy2=npy-4
- diy=4
-
- ix1=2
- ix2=npx-6
- dix=8
-
- dir=4
-
- iax=1
- iay=1
- iar=1
-
- ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-6
- ir2=ir1
-
- doubling_index=IFLAG_BASEMENT_TOPO
-
- else if(isubregion == 7) then
-
- call unusual_hex_nodes3(iaddx,iaddy,iaddz)
-
-! generating stage 7 of the mesh doubling below the Moho
-
- iy1=0
- iy2=npy-4
- diy=4
-
- ix1=4
- ix2=npx-4
- dix=8
-
- dir=4
-
- iax=1
- iay=1
- iar=1
-
- ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-6
- ir2=ir1
-
- doubling_index=IFLAG_BASEMENT_TOPO
-
- else if(isubregion == 6) then
-
- call unusual_hex_nodes4(iaddx,iaddy,iaddz)
-
-! generating stage 8 of the mesh doubling below the Moho
-
- iy1=4
- iy2=npy-4
- diy=8
-
- ix1=0
- ix2=npx-2
- dix=2
-
- dir=4
-
- iax=1
- iay=1
- iar=1
-
- ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-4
- ir2=ir1
-
- doubling_index=IFLAG_BASEMENT_TOPO
-
- else if(isubregion == 5) then
-
- call unusual_hex_nodes4p(iaddx,iaddy,iaddz)
-
-! generating stage 9 of the mesh doubling below the Moho
-
- iy1=0
- iy2=npy-8
- diy=8
-
- ix1=0
- ix2=npx-2
- dix=2
-
- dir=4
-
- iax=1
- iay=1
- iar=1
-
- ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-4
- ir2=ir1
-
- doubling_index=IFLAG_BASEMENT_TOPO
-
- else if(isubregion == 4) then
-
- call usual_hex_nodes(iaddx,iaddy,iaddz)
-
-! generating stage 10 of the mesh doubling below the Moho
-
- iy1=4
- iy2=npy-4
- diy=8
-
- ix1=0
- ix2=npx-2
- dix=2
-
- dir=4
-
- iax=1
- iay=1
- iar=1
-
- ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-2
- ir2=ir1
-
- doubling_index=IFLAG_BASEMENT_TOPO
-
- else if(isubregion == 3) then
-
- call usual_hex_nodes(iaddx,iaddy,iaddz)
-
-! generating stage 11 of the mesh doubling below the Moho
-
- iy1=2
- iy2=npy-6
- diy=8
-
- ix1=0
- ix2=npx-2
- dix=2
-
- dir=4
-
- iax=1
- iay=1
- iar=1
-
- ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-2
- ir2=ir1
-
- doubling_index=IFLAG_BASEMENT_TOPO
-
- else if(isubregion == 2) then
-
- call unusual_hex_nodes6(iaddx,iaddy,iaddz)
-
-! generating stage 12 of the mesh doubling below the Moho
-
- iy1=6
- iy2=npy-2
- diy=8
-
- ix1=0
- ix2=npx-2
- dix=2
-
- dir=4
-
- iax=1
- iay=1
- iar=1
-
- ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-2
- ir2=ir1
-
- doubling_index=IFLAG_BASEMENT_TOPO
-
- else if(isubregion == 1) then
-
- call unusual_hex_nodes6p(iaddx,iaddy,iaddz)
-
-! generating stage 13 of the mesh doubling below the Moho
-
- iy1=0
- iy2=npy-8
- diy=8
-
- ix1=0
- ix2=npx-2
- dix=2
-
- dir=4
-
- iax=1
- iay=1
- iar=1
-
- ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-4
- ir2=ir1
-
- doubling_index=IFLAG_BASEMENT_TOPO
-
- else
-
- call exit_MPI(myrank,'incorrect subregion code')
-
- endif
-
- endif
-
- end subroutine define_subregions
-
Deleted: seismo/3D/SPECFEM3D/trunk/define_subregions_heuristic.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/define_subregions_heuristic.f90 2010-09-24 22:15:18 UTC (rev 17214)
+++ seismo/3D/SPECFEM3D/trunk/define_subregions_heuristic.f90 2010-09-24 22:17:01 UTC (rev 17215)
@@ -1,267 +0,0 @@
-!=====================================================================
-!
-! S p e c f e m 3 D V e r s i o n 1 . 4
-! ---------------------------------------
-!
-! Dimitri Komatitsch and Jeroen Tromp
-! Seismological Laboratory - California Institute of Technology
-! (c) California Institute of Technology September 2006
-!
-! This program is free software; you can redistribute it and/or modify
-! it under the terms of the GNU General Public License as published by
-! the Free Software Foundation; either version 2 of the License, or
-! (at your option) any later version.
-!
-! This program is distributed in the hope that it will be useful,
-! but WITHOUT ANY WARRANTY; without even the implied warranty of
-! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-! GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along
-! with this program; if not, write to the Free Software Foundation, Inc.,
-! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-!
-!=====================================================================
-
- subroutine 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)
-
-! heuristic rule to deform elements to balance angles
-! to 120 degrees in doubling regions
-
- implicit none
-
- include "constants.h"
-
- integer myrank
- integer ix1,ix2,dix,iy1,iy2,diy,ir1,ir2,dir
- integer iax,iay,iar
- integer isubregion,itype_element
- integer npx,npy
-
- integer NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM
-
-! topology of the elements
- integer iaddx(NGNOD)
- integer iaddy(NGNOD)
- integer iaddz(NGNOD)
-
-! type of elements for heuristic rule
- integer, parameter :: ITYPE_UNUSUAL_1 = 1
- integer, parameter :: ITYPE_UNUSUAL_1p = 2
- integer, parameter :: ITYPE_UNUSUAL_4 = 3
- integer, parameter :: ITYPE_UNUSUAL_4p = 4
-
-
-! **************
-
- if(isubregion == 8) then
-
- call unusual_hex_nodes1(iaddx,iaddy,iaddz)
-
-! generating stage 2 of the mesh doubling below 670
-
- iy1=0
- iy2=npy-8
- diy=8
-
- ix1=0
- ix2=npx-16
- dix=16
-
- dir=4
-
- iax=2
- iay=2
- iar=1
-
- ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 8
- ir2=ir1
-
- itype_element = ITYPE_UNUSUAL_1
-
- else if(isubregion == 7) then
-
- call unusual_hex_nodes1p(iaddx,iaddy,iaddz)
-
-! generating stage 3 of the mesh doubling below 670
-
- iy1=0
- iy2=npy-8
- diy=8
-
- ix1=8
- ix2=npx-8
- dix=16
-
- dir=4
-
- iax=2
- iay=2
- iar=1
-
- ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 8
- ir2=ir1
-
- itype_element = ITYPE_UNUSUAL_1p
-
- else if(isubregion == 6) then
-
- call unusual_hex_nodes4(iaddx,iaddy,iaddz)
-
-! generating stage 8 of the mesh doubling below 670
-
- iy1=8
- iy2=npy-8
- diy=16
-
- ix1=0
- ix2=npx-4
- dix=4
-
- dir=4
-
- iax=2
- iay=2
- iar=1
-
- ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 4
- ir2=ir1
-
- itype_element = ITYPE_UNUSUAL_4
-
- else if(isubregion == 5) then
-
- call unusual_hex_nodes4p(iaddx,iaddy,iaddz)
-
-! generating stage 9 of the mesh doubling below 670
-
- iy1=0
- iy2=npy-16
- diy=16
-
- ix1=0
- ix2=npx-4
- dix=4
-
- dir=4
-
- iax=2
- iay=2
- iar=1
-
- ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 4
- ir2=ir1
-
- itype_element = ITYPE_UNUSUAL_4p
-
- else if(isubregion == 4) then
-
- call unusual_hex_nodes1(iaddx,iaddy,iaddz)
-
-! generating stage 2 of the mesh doubling below the Moho
-
- iy1=0
- iy2=npy-4
- diy=4
-
- ix1=0
- ix2=npx-8
- dix=8
-
- dir=4
-
- iax=1
- iay=1
- iar=1
-
- ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-8
- ir2=ir1
-
- itype_element = ITYPE_UNUSUAL_1
-
- else if(isubregion == 3) then
-
- call unusual_hex_nodes1p(iaddx,iaddy,iaddz)
-
-! generating stage 3 of the mesh doubling below the Moho
-
- iy1=0
- iy2=npy-4
- diy=4
-
- ix1=4
- ix2=npx-4
- dix=8
-
- dir=4
-
- iax=1
- iay=1
- iar=1
-
- ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-8
- ir2=ir1
-
- itype_element = ITYPE_UNUSUAL_1p
-
- else if(isubregion == 2) then
-
- call unusual_hex_nodes4(iaddx,iaddy,iaddz)
-
-! generating stage 8 of the mesh doubling below the Moho
-
- iy1=4
- iy2=npy-4
- diy=8
-
- ix1=0
- ix2=npx-2
- dix=2
-
- dir=4
-
- iax=1
- iay=1
- iar=1
-
- ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-4
- ir2=ir1
-
- itype_element = ITYPE_UNUSUAL_4
-
- else if(isubregion == 1) then
-
- call unusual_hex_nodes4p(iaddx,iaddy,iaddz)
-
-! generating stage 9 of the mesh doubling below the Moho
-
- iy1=0
- iy2=npy-8
- diy=8
-
- ix1=0
- ix2=npx-2
- dix=2
-
- dir=4
-
- iax=1
- iay=1
- iar=1
-
- ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-4
- ir2=ir1
-
- itype_element = ITYPE_UNUSUAL_4p
-
- else
-
- call exit_MPI(myrank,'incorrect subregion code')
-
- endif
-
- end subroutine define_subregions_heuristic
-
Deleted: seismo/3D/SPECFEM3D/trunk/get_absorb.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/get_absorb.f90 2010-09-24 22:15:18 UTC (rev 17214)
+++ seismo/3D/SPECFEM3D/trunk/get_absorb.f90 2010-09-24 22:17:01 UTC (rev 17215)
@@ -1,270 +0,0 @@
-!=====================================================================
-!
-! S p e c f e m 3 D V e r s i o n 1 . 4
-! ---------------------------------------
-!
-! Dimitri Komatitsch and Jeroen Tromp
-! Seismological Laboratory - California Institute of Technology
-! (c) California Institute of Technology September 2006
-!
-! This program is free software; you can redistribute it and/or modify
-! it under the terms of the GNU General Public License as published by
-! the Free Software Foundation; either version 2 of the License, or
-! (at your option) any later version.
-!
-! This program is distributed in the hope that it will be useful,
-! but WITHOUT ANY WARRANTY; without even the implied warranty of
-! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-! GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along
-! with this program; if not, write to the Free Software Foundation, Inc.,
-! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-!
-!=====================================================================
-
- subroutine get_absorb(myrank,prname,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
-
-! processor identification
- character(len=256) prname
-
- 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')
-
-! IMPROVE save these temporary arrays for the solver for Stacey conditions
-
- open(unit=27,file=prname(1:len_trim(prname))//'nimin.bin',status='unknown',form='unformatted')
- write(27) nimin
- close(27)
-
- open(unit=27,file=prname(1:len_trim(prname))//'nimax.bin',status='unknown',form='unformatted')
- write(27) nimax
- close(27)
-
- open(unit=27,file=prname(1:len_trim(prname))//'njmin.bin',status='unknown',form='unformatted')
- write(27) njmin
- close(27)
-
- open(unit=27,file=prname(1:len_trim(prname))//'njmax.bin',status='unknown',form='unformatted')
- write(27) njmax
- close(27)
-
- open(unit=27,file=prname(1:len_trim(prname))//'nkmin_xi.bin',status='unknown',form='unformatted')
- write(27) nkmin_xi
- close(27)
-
- open(unit=27,file=prname(1:len_trim(prname))//'nkmin_eta.bin',status='unknown',form='unformatted')
- write(27) nkmin_eta
- close(27)
-
- 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
-
-
-
-
Deleted: seismo/3D/SPECFEM3D/trunk/get_flags_boundaries.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/get_flags_boundaries.f90 2010-09-24 22:15:18 UTC (rev 17214)
+++ seismo/3D/SPECFEM3D/trunk/get_flags_boundaries.f90 2010-09-24 22:17:01 UTC (rev 17215)
@@ -1,162 +0,0 @@
-!=====================================================================
-!
-! S p e c f e m 3 D V e r s i o n 1 . 4
-! ---------------------------------------
-!
-! Dimitri Komatitsch and Jeroen Tromp
-! Seismological Laboratory - California Institute of Technology
-! (c) California Institute of Technology September 2006
-!
-! This program is free software; you can redistribute it and/or modify
-! it under the terms of the GNU General Public License as published by
-! the Free Software Foundation; either version 2 of the License, or
-! (at your option) any later version.
-!
-! This program is distributed in the hope that it will be useful,
-! but WITHOUT ANY WARRANTY; without even the implied warranty of
-! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-! GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along
-! with this program; if not, write to the Free Software Foundation, Inc.,
-! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-!
-!=====================================================================
-
- subroutine get_flags_boundaries(nspec,iproc_xi,iproc_eta,ispec,idoubling, &
- xstore,ystore,zstore,iboun,iMPIcut_xi,iMPIcut_eta, &
- NPROC_XI,NPROC_ETA, &
- UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK)
-
- implicit none
-
- include "constants.h"
- include "constants_gocad.h"
-
- integer nspec
- integer ispec,idoubling
- integer NPROC_XI,NPROC_ETA
-
- double precision UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK
-
- logical iboun(6,nspec)
- logical iMPIcut_xi(2,nspec),iMPIcut_eta(2,nspec)
-
- double precision xstore(NGLLX,NGLLY,NGLLZ)
- double precision ystore(NGLLX,NGLLY,NGLLZ)
- double precision zstore(NGLLX,NGLLY,NGLLZ)
-
-! use iproc_xi and iproc_eta to determine MPI cut planes along xi and eta
- integer iproc_xi,iproc_eta
-
- double precision target,sizeslice,TOLERANCE_METERS
- double precision xelm(8),yelm(8),zelm(8)
-
-! find the coordinates of the eight corner nodes of the element
- xelm(1)=xstore(1,1,1)
- yelm(1)=ystore(1,1,1)
- zelm(1)=zstore(1,1,1)
- xelm(2)=xstore(NGLLX,1,1)
- yelm(2)=ystore(NGLLX,1,1)
- zelm(2)=zstore(NGLLX,1,1)
- xelm(3)=xstore(NGLLX,NGLLY,1)
- yelm(3)=ystore(NGLLX,NGLLY,1)
- zelm(3)=zstore(NGLLX,NGLLY,1)
- xelm(4)=xstore(1,NGLLY,1)
- yelm(4)=ystore(1,NGLLY,1)
- zelm(4)=zstore(1,NGLLY,1)
- xelm(5)=xstore(1,1,NGLLZ)
- yelm(5)=ystore(1,1,NGLLZ)
- zelm(5)=zstore(1,1,NGLLZ)
- xelm(6)=xstore(NGLLX,1,NGLLZ)
- yelm(6)=ystore(NGLLX,1,NGLLZ)
- zelm(6)=zstore(NGLLX,1,NGLLZ)
- xelm(7)=xstore(NGLLX,NGLLY,NGLLZ)
- yelm(7)=ystore(NGLLX,NGLLY,NGLLZ)
- zelm(7)=zstore(NGLLX,NGLLY,NGLLZ)
- xelm(8)=xstore(1,NGLLY,NGLLZ)
- yelm(8)=ystore(1,NGLLY,NGLLZ)
- zelm(8)=zstore(1,NGLLY,NGLLZ)
-
-! compute geometrical tolerance small compared to size of model to detect edges
- TOLERANCE_METERS = dabs(UTM_X_MAX - UTM_X_MIN) / 100000.
-
-! ****************************************************
-! determine if the element falls on a boundary
-! ****************************************************
-
- iboun(:,ispec)=.false.
-
-! on boundary 1: x=xmin
- target= UTM_X_MIN + TOLERANCE_METERS
- if(xelm(1)<target .and. xelm(4)<target .and. xelm(5)<target .and. xelm(8)<target) iboun(1,ispec)=.true.
-
-! on boundary 2: xmax
- target= UTM_X_MAX - TOLERANCE_METERS
- if(xelm(2)>target .and. xelm(3)>target .and. xelm(6)>target .and. xelm(7)>target) iboun(2,ispec)=.true.
-
-! on boundary 3: ymin
- target= UTM_Y_MIN + TOLERANCE_METERS
- if(yelm(1)<target .and. yelm(2)<target .and. yelm(5)<target .and. yelm(6)<target) iboun(3,ispec)=.true.
-
-! on boundary 4: ymax
- target= UTM_Y_MAX - TOLERANCE_METERS
- if(yelm(3)>target .and. yelm(4)>target .and. yelm(7)>target .and. yelm(8)>target) iboun(4,ispec)=.true.
-
-! on boundary 5: bottom
- target = Z_DEPTH_BLOCK + TOLERANCE_METERS
- if(zelm(1)<target .and. zelm(2)<target .and. zelm(3)<target .and. zelm(4)<target) iboun(5,ispec)=.true.
-
-! on boundary 6: top
- if(idoubling == IFLAG_ONE_LAYER_TOPOGRAPHY) iboun(6,ispec)=.true.
-
-! *******************************************************************
-! determine if the element falls on an MPI cut plane along xi
-! *******************************************************************
-
-! detect the MPI cut planes along xi in the cubed sphere
-
- iMPIcut_xi(:,ispec)=.false.
-
-! angular size of a slice along xi
- sizeslice = (UTM_X_MAX-UTM_X_MIN) / NPROC_XI
-
-! left cut-plane in the current slice along X = constant (Xmin of this slice)
-! and add geometrical tolerance
-
- target = UTM_X_MIN + iproc_xi*sizeslice + TOLERANCE_METERS
- if(xelm(1)<target .and. xelm(4)<target .and. xelm(5)<target .and. xelm(8)<target) &
- iMPIcut_xi(1,ispec)=.true.
-
-! right cut-plane in the current slice along X = constant (Xmax of this slice)
-! and add geometrical tolerance
-
- target = UTM_X_MIN + (iproc_xi+1)*sizeslice - TOLERANCE_METERS
- if(xelm(2)>target .and. xelm(3)>target .and. xelm(6)>target .and. xelm(7)>target) &
- iMPIcut_xi(2,ispec)=.true.
-
-! ********************************************************************
-! determine if the element falls on an MPI cut plane along eta
-! ********************************************************************
-
- iMPIcut_eta(:,ispec)=.false.
-
-! angular size of a slice along eta
- sizeslice = (UTM_Y_MAX-UTM_Y_MIN) / NPROC_ETA
-
-! left cut-plane in the current slice along Y = constant (Ymin of this slice)
-! and add geometrical tolerance
-
- target = UTM_Y_MIN + iproc_eta*sizeslice + TOLERANCE_METERS
- if(yelm(1)<target .and. yelm(2)<target .and. yelm(5)<target .and. yelm(6)<target) &
- iMPIcut_eta(1,ispec)=.true.
-
-! right cut-plane in the current slice along Y = constant (Ymax of this slice)
-! and add geometrical tolerance
-
- target = UTM_Y_MIN + (iproc_eta+1)*sizeslice - TOLERANCE_METERS
- if(yelm(3)>target .and. yelm(4)>target .and. yelm(7)>target .and. yelm(8)>target) &
- iMPIcut_eta(2,ispec)=.true.
-
- end subroutine get_flags_boundaries
-
Deleted: seismo/3D/SPECFEM3D/trunk/hauksson_model.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/hauksson_model.f90 2010-09-24 22:15:18 UTC (rev 17214)
+++ seismo/3D/SPECFEM3D/trunk/hauksson_model.f90 2010-09-24 22:17:01 UTC (rev 17215)
@@ -1,227 +0,0 @@
-!=====================================================================
-!
-! S p e c f e m 3 D V e r s i o n 1 . 4
-! ---------------------------------------
-!
-! Dimitri Komatitsch and Jeroen Tromp
-! Seismological Laboratory - California Institute of Technology
-! (c) California Institute of Technology September 2006
-!
-! This program is free software; you can redistribute it and/or modify
-! it under the terms of the GNU General Public License as published by
-! the Free Software Foundation; either version 2 of the License, or
-! (at your option) any later version.
-!
-! This program is distributed in the hope that it will be useful,
-! but WITHOUT ANY WARRANTY; without even the implied warranty of
-! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-! GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along
-! with this program; if not, write to the Free Software Foundation, Inc.,
-! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-!
-!=====================================================================
-
- subroutine hauksson_model(vp,vs,utm_x_eval,utm_y_eval,z_eval,vp_final,vs_final,MOHO_MAP_LUPEI)
-
- implicit none
-
- include "constants.h"
- include "constants_gocad.h"
-
-!! DK DK UGLY one day, we should clarify the issue of merging Hauksson's Moho
-!! DK DK UGLY with our Lupei Moho. Should not be a big issue because in
-!! DK DK UGLY principle Hauksson used Lupei's map to build his model
-
- double precision utm_x_eval,utm_y_eval,z_eval
- double precision vp_final,vs_final
- logical MOHO_MAP_LUPEI
-
- double precision, dimension(NLAYERS_HAUKSSON,NGRID_NEW_HAUKSSON,NGRID_NEW_HAUKSSON) :: vp,vs
- double precision, dimension(NLAYERS_HAUKSSON) :: vp_interp,vs_interp
-
- integer ilayer
- integer icell_interp_x,icell_interp_y
- double precision spacing_x,spacing_y
- double precision utm_x_eval_copy,utm_y_eval_copy
- double precision gamma_interp_x,gamma_interp_y,gamma_interp_z
- double precision v1,v2,v3,v4
- double precision vp_upper,vs_upper,vp_lower,vs_lower,z_upper,z_lower
-
-! copy input values
- utm_x_eval_copy = utm_x_eval
- utm_y_eval_copy = utm_y_eval
-
-! make sure we stay inside Hauksson's grid
- if(utm_x_eval_copy < UTM_X_ORIG_HAUKSSON) utm_x_eval_copy = UTM_X_ORIG_HAUKSSON
- if(utm_y_eval_copy < UTM_Y_ORIG_HAUKSSON) utm_y_eval_copy = UTM_Y_ORIG_HAUKSSON
-
-! determine spacing and cell for linear interpolation
- spacing_x = (utm_x_eval_copy - UTM_X_ORIG_HAUKSSON) / SPACING_UTM_X_HAUKSSON
- spacing_y = (utm_y_eval_copy - UTM_Y_ORIG_HAUKSSON) / SPACING_UTM_Y_HAUKSSON
-
- icell_interp_x = int(spacing_x) + 1
- icell_interp_y = int(spacing_y) + 1
-
- gamma_interp_x = spacing_x - int(spacing_x)
- gamma_interp_y = spacing_y - int(spacing_y)
-
-! suppress edge effects for points outside of Hauksson's model
- if(icell_interp_x < 1) then
- icell_interp_x = 1
- gamma_interp_x = 0.d0
- endif
- if(icell_interp_x > NGRID_NEW_HAUKSSON-1) then
- icell_interp_x = NGRID_NEW_HAUKSSON-1
- gamma_interp_x = 1.d0
- endif
-
- if(icell_interp_y < 1) then
- icell_interp_y = 1
- gamma_interp_y = 0.d0
- endif
- if(icell_interp_y > NGRID_NEW_HAUKSSON-1) then
- icell_interp_y = NGRID_NEW_HAUKSSON-1
- gamma_interp_y = 1.d0
- endif
-
-! make sure interpolation makes sense
- if(gamma_interp_x < -0.001d0 .or. gamma_interp_x > 1.001d0) &
- stop 'interpolation in x is incorrect in Hauksson'
- if(gamma_interp_y < -0.001d0 .or. gamma_interp_y > 1.001d0) &
- stop 'interpolation in y is incorrect in Hauksson'
-
-! interpolate Hauksson's model at right location using bilinear interpolation
- do ilayer = 1,NLAYERS_HAUKSSON
-
-! for Vp
- v1 = vp(ilayer,icell_interp_x,icell_interp_y)
- v2 = vp(ilayer,icell_interp_x+1,icell_interp_y)
- v3 = vp(ilayer,icell_interp_x+1,icell_interp_y+1)
- v4 = vp(ilayer,icell_interp_x,icell_interp_y+1)
-
- vp_interp(ilayer) = v1*(1.-gamma_interp_x)*(1.-gamma_interp_y) + &
- v2*gamma_interp_x*(1.-gamma_interp_y) + &
- v3*gamma_interp_x*gamma_interp_y + &
- v4*(1.-gamma_interp_x)*gamma_interp_y
-
-! for Vs
- v1 = vs(ilayer,icell_interp_x,icell_interp_y)
- v2 = vs(ilayer,icell_interp_x+1,icell_interp_y)
- v3 = vs(ilayer,icell_interp_x+1,icell_interp_y+1)
- v4 = vs(ilayer,icell_interp_x,icell_interp_y+1)
-
- vs_interp(ilayer) = v1*(1.-gamma_interp_x)*(1.-gamma_interp_y) + &
- v2*gamma_interp_x*(1.-gamma_interp_y) + &
- v3*gamma_interp_x*gamma_interp_y + &
- v4*(1.-gamma_interp_x)*gamma_interp_y
-
- enddo
-
-! choose right values depending on depth of target point
- if(z_eval >= Z_HAUKSSON_LAYER_1) then
- vp_final = vp_interp(1)
- vs_final = vs_interp(1)
- return
-
- else if(z_eval <= Z_HAUKSSON_LAYER_9) then
- vp_final = vp_interp(9)
- vs_final = vs_interp(9)
- return
-
- else if(z_eval >= Z_HAUKSSON_LAYER_2) then
- vp_upper = vp_interp(1)
- vs_upper = vs_interp(1)
- z_upper = Z_HAUKSSON_LAYER_1
-
- vp_lower = vp_interp(2)
- vs_lower = vs_interp(2)
- z_lower = Z_HAUKSSON_LAYER_2
-
- else if(z_eval >= Z_HAUKSSON_LAYER_3) then
- vp_upper = vp_interp(2)
- vs_upper = vs_interp(2)
- z_upper = Z_HAUKSSON_LAYER_2
-
- vp_lower = vp_interp(3)
- vs_lower = vs_interp(3)
- z_lower = Z_HAUKSSON_LAYER_3
-
- else if(z_eval >= Z_HAUKSSON_LAYER_4) then
- vp_upper = vp_interp(3)
- vs_upper = vs_interp(3)
- z_upper = Z_HAUKSSON_LAYER_3
-
- vp_lower = vp_interp(4)
- vs_lower = vs_interp(4)
- z_lower = Z_HAUKSSON_LAYER_4
-
- else if(z_eval >= Z_HAUKSSON_LAYER_5) then
- vp_upper = vp_interp(4)
- vs_upper = vs_interp(4)
- z_upper = Z_HAUKSSON_LAYER_4
-
- vp_lower = vp_interp(5)
- vs_lower = vs_interp(5)
- z_lower = Z_HAUKSSON_LAYER_5
-
- else if(z_eval >= Z_HAUKSSON_LAYER_6) then
- vp_upper = vp_interp(5)
- vs_upper = vs_interp(5)
- z_upper = Z_HAUKSSON_LAYER_5
-
- vp_lower = vp_interp(6)
- vs_lower = vs_interp(6)
- z_lower = Z_HAUKSSON_LAYER_6
-
- else if(z_eval >= Z_HAUKSSON_LAYER_7) then
- vp_upper = vp_interp(6)
- vs_upper = vs_interp(6)
- z_upper = Z_HAUKSSON_LAYER_6
-
- vp_lower = vp_interp(7)
- vs_lower = vs_interp(7)
- z_lower = Z_HAUKSSON_LAYER_7
-
- else if(z_eval >= Z_HAUKSSON_LAYER_8) then
- vp_upper = vp_interp(7)
- vs_upper = vs_interp(7)
- z_upper = Z_HAUKSSON_LAYER_7
-
- vp_lower = vp_interp(8)
- vs_lower = vs_interp(8)
- z_lower = Z_HAUKSSON_LAYER_8
-
- else
- if(.not. MOHO_MAP_LUPEI) then
- vp_upper = vp_interp(8)
- vs_upper = vs_interp(8)
- z_upper = Z_HAUKSSON_LAYER_8
-
- vp_lower = vp_interp(9)
- vs_lower = vs_interp(9)
- z_lower = Z_HAUKSSON_LAYER_9
- !!! waiting for better interpolation of Moho maps.
- else
- vp_upper = vp_interp(8)
- vs_upper = vs_interp(8)
- z_upper = Z_HAUKSSON_LAYER_8
-
- vp_lower = vp_interp(9)
- vs_lower = vs_interp(9)
- z_lower = Z_HAUKSSON_LAYER_9
- endif
-
- endif
-
- gamma_interp_z = (z_eval - z_lower) / (z_upper - z_lower)
-
- if(gamma_interp_z < -0.001d0 .or. gamma_interp_z > 1.001d0) &
- stop 'interpolation in z is incorrect in Hauksson'
-
- vp_final = vp_upper * gamma_interp_z + vp_lower * (1.-gamma_interp_z)
- vs_final = vs_upper * gamma_interp_z + vs_lower * (1.-gamma_interp_z)
-
- end subroutine hauksson_model
-
Deleted: seismo/3D/SPECFEM3D/trunk/interpolate_gocad_block_HR.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/interpolate_gocad_block_HR.f90 2010-09-24 22:15:18 UTC (rev 17214)
+++ seismo/3D/SPECFEM3D/trunk/interpolate_gocad_block_HR.f90 2010-09-24 22:17:01 UTC (rev 17215)
@@ -1,161 +0,0 @@
-
- subroutine interpolate_gocad_block_HR(vp_block_gocad_HR,vp_block_gocad_MR, &
- utm_x_eval,utm_y_eval,z_eval,rho_final,vp_final,vs_final,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)
-
- implicit none
-
- include "constants.h"
- include "constants_gocad.h"
-
- double precision vp_block_gocad_HR(0:NX_GOCAD_HR-1,0:NY_GOCAD_HR-1,0:NZ_GOCAD_HR-1)
- double precision vp_block_gocad_MR(0:NX_GOCAD_MR-1,0:NY_GOCAD_MR-1,0:NZ_GOCAD_MR-1)
-
- integer ix,iy,iz
-
- double precision utm_x_eval,utm_y_eval,z_eval
- double precision spacing_x,spacing_y,spacing_z
- double precision gamma_interp_x,gamma_interp_y,gamma_interp_z
- double precision v1,v2,v3,v4,v5,v6,v7,v8
- double precision vp_final,vs_final,rho_final,vp_vs_ratio
- double precision rho_ref_MR,vp_ref_MR,vs_ref_MR
- double precision THICKNESS_TAPER_BLOCK_HR, &
- VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM
-
- logical point_is_in_sediments,dummy_flag,IMPOSE_MINIMUM_VP_GOCAD
-
-! for Hauksson's model
- integer doubling_index
- logical HAUKSSON_REGIONAL_MODEL,MOHO_MAP_LUPEI
- double precision, dimension(NLAYERS_HAUKSSON,NGRID_NEW_HAUKSSON,NGRID_NEW_HAUKSSON) :: vp_hauksson,vs_hauksson
-
-! determine spacing and cell for linear interpolation
- spacing_x = (utm_x_eval - ORIG_X_GOCAD_HR) / SPACING_X_GOCAD_HR
- spacing_y = (utm_y_eval - ORIG_Y_GOCAD_HR) / SPACING_Y_GOCAD_HR
- spacing_z = (z_eval - ORIG_Z_GOCAD_HR) / SPACING_Z_GOCAD_HR
-
- ix = int(spacing_x)
- iy = int(spacing_y)
- iz = int(spacing_z)
-
- gamma_interp_x = spacing_x - dble(ix)
- gamma_interp_y = spacing_y - dble(iy)
- gamma_interp_z = spacing_z - dble(iz)
-
-! this block is smaller than the grid, therefore just exit
-! if the target point is outside of the block
- if(ix < 0 .or. ix > NX_GOCAD_HR-2 .or. iy < 0 .or. iy > NY_GOCAD_HR-2) return
-
-! suppress edge effects in vertical direction
- if(iz < 0) then
- iz = 0
- gamma_interp_z = 0.d0
- endif
- if(iz > NZ_GOCAD_HR-2) then
- iz = NZ_GOCAD_HR-2
- gamma_interp_z = 1.d0
- endif
-
-! define 8 corners of interpolation element
- v1 = vp_block_gocad_HR(ix,iy,iz)
- v2 = vp_block_gocad_HR(ix+1,iy,iz)
- v3 = vp_block_gocad_HR(ix+1,iy+1,iz)
- v4 = vp_block_gocad_HR(ix,iy+1,iz)
-
- v5 = vp_block_gocad_HR(ix,iy,iz+1)
- v6 = vp_block_gocad_HR(ix+1,iy,iz+1)
- v7 = vp_block_gocad_HR(ix+1,iy+1,iz+1)
- v8 = vp_block_gocad_HR(ix,iy+1,iz+1)
-
-! check if element is defined (i.e. is in the sediments in Voxet)
-! do nothing if element is undefined
-! a P-velocity of 20 km/s is used to indicate fictitious elements
- if(v1 < 19000. .and. v2 < 19000. .and. &
- v3 < 19000. .and. v4 < 19000. .and. &
- v5 < 19000. .and. v6 < 19000. .and. &
- v7 < 19000. .and. v8 < 19000.) then
-
-! set flag indicating whether point is in the sediments
- point_is_in_sediments = .true.
-
-! use trilinear interpolation in cell to define Vp
- vp_final = &
- v1*(1.-gamma_interp_x)*(1.-gamma_interp_y)*(1.-gamma_interp_z) + &
- v2*gamma_interp_x*(1.-gamma_interp_y)*(1.-gamma_interp_z) + &
- v3*gamma_interp_x*gamma_interp_y*(1.-gamma_interp_z) + &
- v4*(1.-gamma_interp_x)*gamma_interp_y*(1.-gamma_interp_z) + &
- v5*(1.-gamma_interp_x)*(1.-gamma_interp_y)*gamma_interp_z + &
- v6*gamma_interp_x*(1.-gamma_interp_y)*gamma_interp_z + &
- v7*gamma_interp_x*gamma_interp_y*gamma_interp_z + &
- v8*(1.-gamma_interp_x)*gamma_interp_y*gamma_interp_z
-
-! impose minimum velocity if needed
- if(IMPOSE_MINIMUM_VP_GOCAD .and. vp_final < VP_MIN_GOCAD) vp_final = VP_MIN_GOCAD
-
-! taper edges to make smooth transition between MR and HR blocks
-! get value from edge of medium-resolution block
-! then use linear interpolation from edge of the model
- if(TAPER_GOCAD_TRANSITIONS) then
-
-! x = xmin
- if(utm_x_eval < ORIG_X_GOCAD_HR + THICKNESS_TAPER_BLOCK_HR) then
- gamma_interp_x = (utm_x_eval - ORIG_X_GOCAD_HR) / THICKNESS_TAPER_BLOCK_HR
- call interpolate_gocad_block_MR(vp_block_gocad_MR, &
- ORIG_X_GOCAD_HR,utm_y_eval,z_eval,rho_ref_MR,vp_ref_MR,vs_ref_MR,dummy_flag, &
- 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)
- vp_final = vp_ref_MR * (1. - gamma_interp_x) + vp_final * gamma_interp_x
-
-! x = xmax
- else if(utm_x_eval > END_X_GOCAD_HR - THICKNESS_TAPER_BLOCK_HR) then
- gamma_interp_x = (utm_x_eval - (END_X_GOCAD_HR - THICKNESS_TAPER_BLOCK_HR)) / THICKNESS_TAPER_BLOCK_HR
- call interpolate_gocad_block_MR(vp_block_gocad_MR, &
- END_X_GOCAD_HR,utm_y_eval,z_eval,rho_ref_MR,vp_ref_MR,vs_ref_MR,dummy_flag, &
- 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)
- vp_final = vp_ref_MR * gamma_interp_x + vp_final * (1. - gamma_interp_x)
-
-! y = ymin
- else if(utm_y_eval < ORIG_Y_GOCAD_HR + THICKNESS_TAPER_BLOCK_HR) then
- gamma_interp_y = (utm_y_eval - ORIG_Y_GOCAD_HR) / THICKNESS_TAPER_BLOCK_HR
- call interpolate_gocad_block_MR(vp_block_gocad_MR, &
- utm_x_eval,ORIG_Y_GOCAD_HR,z_eval,rho_ref_MR,vp_ref_MR,vs_ref_MR,dummy_flag, &
- 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)
- vp_final = vp_ref_MR * (1. - gamma_interp_y) + vp_final * gamma_interp_y
-
-! y = ymax
- else if(utm_y_eval > END_Y_GOCAD_HR - THICKNESS_TAPER_BLOCK_HR) then
- gamma_interp_y = (utm_y_eval - (END_Y_GOCAD_HR - THICKNESS_TAPER_BLOCK_HR)) / THICKNESS_TAPER_BLOCK_HR
- call interpolate_gocad_block_MR(vp_block_gocad_MR, &
- utm_x_eval,END_Y_GOCAD_HR,z_eval,rho_ref_MR,vp_ref_MR,vs_ref_MR,dummy_flag, &
- 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)
- vp_final = vp_ref_MR * gamma_interp_y + vp_final * (1. - gamma_interp_y)
-
- endif
-
- endif
-
-! use linear variation of vp/vs ratio with depth, between 0. and 8.5 km
- vp_vs_ratio = VP_VS_RATIO_GOCAD_BOTTOM + &
- (VP_VS_RATIO_GOCAD_TOP - VP_VS_RATIO_GOCAD_BOTTOM) * &
- (z_eval - (-8500.d0)) / (0.d0 - (-8500.d0))
-
-! make sure ratio remains in interval
- if(vp_vs_ratio < VP_VS_RATIO_GOCAD_BOTTOM) vp_vs_ratio = VP_VS_RATIO_GOCAD_BOTTOM
- if(vp_vs_ratio > VP_VS_RATIO_GOCAD_TOP) vp_vs_ratio = VP_VS_RATIO_GOCAD_TOP
-
- vs_final = vp_final / vp_vs_ratio
- call compute_rho_estimate(rho_final,vp_final)
-
- endif
-
- end subroutine interpolate_gocad_block_HR
-
Deleted: seismo/3D/SPECFEM3D/trunk/interpolate_gocad_block_MR.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/interpolate_gocad_block_MR.f90 2010-09-24 22:15:18 UTC (rev 17214)
+++ seismo/3D/SPECFEM3D/trunk/interpolate_gocad_block_MR.f90 2010-09-24 22:17:01 UTC (rev 17215)
@@ -1,186 +0,0 @@
-
- subroutine interpolate_gocad_block_MR(vp_block_gocad_MR, &
- utm_x_eval,utm_y_eval,z_eval,rho_final,vp_final,vs_final,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)
-
- implicit none
-
- include "constants.h"
- include "constants_gocad.h"
-
- double precision vp_block_gocad_MR(0:NX_GOCAD_MR-1,0:NY_GOCAD_MR-1,0:NZ_GOCAD_MR-1)
- double precision VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM,THICKNESS_TAPER_BLOCK_MR
-
- integer ix,iy,iz
-
- double precision utm_x_eval,utm_y_eval,z_eval
- double precision spacing_x,spacing_y,spacing_z
- double precision gamma_interp_x,gamma_interp_y,gamma_interp_z
- double precision v1,v2,v3,v4,v5,v6,v7,v8
- double precision vp_final,vs_final,rho_final,vp_vs_ratio
- double precision xmesh,ymesh,zmesh,vs_dummy,rho_dummy
-
- logical point_is_in_sediments,IMPOSE_MINIMUM_VP_GOCAD
-
-! for Hauksson's model
- integer doubling_index
- logical HAUKSSON_REGIONAL_MODEL,MOHO_MAP_LUPEI
- double precision vp_ref_hauksson
- double precision, dimension(NLAYERS_HAUKSSON,NGRID_NEW_HAUKSSON,NGRID_NEW_HAUKSSON) :: vp_hauksson,vs_hauksson
-
-! determine spacing and cell for linear interpolation
- spacing_x = (utm_x_eval - ORIG_X_GOCAD_MR) / SPACING_X_GOCAD_MR
- spacing_y = (utm_y_eval - ORIG_Y_GOCAD_MR) / SPACING_Y_GOCAD_MR
- spacing_z = (z_eval - ORIG_Z_GOCAD_MR) / SPACING_Z_GOCAD_MR
-
- ix = int(spacing_x)
- iy = int(spacing_y)
- iz = int(spacing_z)
-
- gamma_interp_x = spacing_x - dble(ix)
- gamma_interp_y = spacing_y - dble(iy)
- gamma_interp_z = spacing_z - dble(iz)
-
-! suppress edge effects for points outside of Gocad model
- if(ix < 0) then
- ix = 0
- gamma_interp_x = 0.d0
- endif
- if(ix > NX_GOCAD_MR-2) then
- ix = NX_GOCAD_MR-2
- gamma_interp_x = 1.d0
- endif
-
- if(iy < 0) then
- iy = 0
- gamma_interp_y = 0.d0
- endif
- if(iy > NY_GOCAD_MR-2) then
- iy = NY_GOCAD_MR-2
- gamma_interp_y = 1.d0
- endif
-
- if(iz < 0) then
- iz = 0
- gamma_interp_z = 0.d0
- endif
- if(iz > NZ_GOCAD_MR-2) then
- iz = NZ_GOCAD_MR-2
- gamma_interp_z = 1.d0
- endif
-
-! define 8 corners of interpolation element
- v1 = vp_block_gocad_MR(ix,iy,iz)
- v2 = vp_block_gocad_MR(ix+1,iy,iz)
- v3 = vp_block_gocad_MR(ix+1,iy+1,iz)
- v4 = vp_block_gocad_MR(ix,iy+1,iz)
-
- v5 = vp_block_gocad_MR(ix,iy,iz+1)
- v6 = vp_block_gocad_MR(ix+1,iy,iz+1)
- v7 = vp_block_gocad_MR(ix+1,iy+1,iz+1)
- v8 = vp_block_gocad_MR(ix,iy+1,iz+1)
-
-! check if element is defined (i.e. is in the sediments in Voxet)
-! do nothing if element is undefined
-! a P-velocity of 20 km/s is used to indicate fictitious elements
- if(v1 < 19000. .and. v2 < 19000. .and. &
- v3 < 19000. .and. v4 < 19000. .and. &
- v5 < 19000. .and. v6 < 19000. .and. &
- v7 < 19000. .and. v8 < 19000.) then
-
-! set flag indicating whether point is in the sediments
- point_is_in_sediments = .true.
-
-! use trilinear interpolation in cell to define Vp
- vp_final = &
- v1*(1.-gamma_interp_x)*(1.-gamma_interp_y)*(1.-gamma_interp_z) + &
- v2*gamma_interp_x*(1.-gamma_interp_y)*(1.-gamma_interp_z) + &
- v3*gamma_interp_x*gamma_interp_y*(1.-gamma_interp_z) + &
- v4*(1.-gamma_interp_x)*gamma_interp_y*(1.-gamma_interp_z) + &
- v5*(1.-gamma_interp_x)*(1.-gamma_interp_y)*gamma_interp_z + &
- v6*gamma_interp_x*(1.-gamma_interp_y)*gamma_interp_z + &
- v7*gamma_interp_x*gamma_interp_y*gamma_interp_z + &
- v8*(1.-gamma_interp_x)*gamma_interp_y*gamma_interp_z
-
-! impose minimum velocity if needed
- if(IMPOSE_MINIMUM_VP_GOCAD .and. vp_final < VP_MIN_GOCAD) vp_final = VP_MIN_GOCAD
-
-! taper edges to make smooth transition between Hauksson and MR blocks
-! get value from edge of medium-resolution block
-! then use linear interpolation from edge of the model
- if(TAPER_GOCAD_TRANSITIONS) then
-
-! x = xmin
- if(utm_x_eval < ORIG_X_GOCAD_MR + THICKNESS_TAPER_BLOCK_MR) then
- xmesh = ORIG_X_GOCAD_MR
- ymesh = utm_y_eval
- zmesh = z_eval
- if(HAUKSSON_REGIONAL_MODEL) then
- call hauksson_model(vp_hauksson,vs_hauksson,xmesh,ymesh,zmesh,vp_ref_hauksson,vs_dummy, MOHO_MAP_LUPEI)
- else
- call socal_model(doubling_index,rho_dummy,vp_ref_hauksson,vs_dummy)
- endif
- gamma_interp_x = (utm_x_eval - ORIG_X_GOCAD_MR) / THICKNESS_TAPER_BLOCK_MR
- vp_final = vp_ref_hauksson * (1. - gamma_interp_x) + vp_final * gamma_interp_x
-
-! x = xmax
- else if(utm_x_eval > END_X_GOCAD_MR - THICKNESS_TAPER_BLOCK_MR) then
- xmesh = END_X_GOCAD_MR
- ymesh = utm_y_eval
- zmesh = z_eval
- if(HAUKSSON_REGIONAL_MODEL) then
- call hauksson_model(vp_hauksson,vs_hauksson,xmesh,ymesh,zmesh,vp_ref_hauksson,vs_dummy, MOHO_MAP_LUPEI)
- else
- call socal_model(doubling_index,rho_dummy,vp_ref_hauksson,vs_dummy)
- endif
- gamma_interp_x = (utm_x_eval - (END_X_GOCAD_MR - THICKNESS_TAPER_BLOCK_MR)) / THICKNESS_TAPER_BLOCK_MR
- vp_final = vp_ref_hauksson * gamma_interp_x + vp_final * (1. - gamma_interp_x)
-
-! y = ymin
- else if(utm_y_eval < ORIG_Y_GOCAD_MR + THICKNESS_TAPER_BLOCK_MR) then
- xmesh = utm_x_eval
- ymesh = ORIG_Y_GOCAD_MR
- zmesh = z_eval
- if(HAUKSSON_REGIONAL_MODEL) then
- call hauksson_model(vp_hauksson,vs_hauksson,xmesh,ymesh,zmesh,vp_ref_hauksson,vs_dummy, MOHO_MAP_LUPEI)
- else
- call socal_model(doubling_index,rho_dummy,vp_ref_hauksson,vs_dummy)
- endif
- gamma_interp_y = (utm_y_eval - ORIG_Y_GOCAD_MR) / THICKNESS_TAPER_BLOCK_MR
- vp_final = vp_ref_hauksson * (1. - gamma_interp_y) + vp_final * gamma_interp_y
-
-! y = ymax
- else if(utm_y_eval > END_Y_GOCAD_MR - THICKNESS_TAPER_BLOCK_MR) then
- xmesh = utm_x_eval
- ymesh = END_Y_GOCAD_MR
- zmesh = z_eval
- if(HAUKSSON_REGIONAL_MODEL) then
- call hauksson_model(vp_hauksson,vs_hauksson,xmesh,ymesh,zmesh,vp_ref_hauksson,vs_dummy, MOHO_MAP_LUPEI)
- else
- call socal_model(doubling_index,rho_dummy,vp_ref_hauksson,vs_dummy)
- endif
- gamma_interp_y = (utm_y_eval - (END_Y_GOCAD_MR - THICKNESS_TAPER_BLOCK_MR)) / THICKNESS_TAPER_BLOCK_MR
- vp_final = vp_ref_hauksson * gamma_interp_y + vp_final * (1. - gamma_interp_y)
-
- endif
-
- endif
-
-! use linear variation of vp/vs ratio with depth, between 0. and 8.5 km
- vp_vs_ratio = VP_VS_RATIO_GOCAD_BOTTOM + &
- (VP_VS_RATIO_GOCAD_TOP - VP_VS_RATIO_GOCAD_BOTTOM) * &
- (z_eval - (-8500.d0)) / (0.d0 - (-8500.d0))
-
-! make sure ratio remains in interval
- if(vp_vs_ratio < VP_VS_RATIO_GOCAD_BOTTOM) vp_vs_ratio = VP_VS_RATIO_GOCAD_BOTTOM
- if(vp_vs_ratio > VP_VS_RATIO_GOCAD_TOP) vp_vs_ratio = VP_VS_RATIO_GOCAD_TOP
-
- vs_final = vp_final / vp_vs_ratio
- call compute_rho_estimate(rho_final,vp_final)
-
- endif
-
- end subroutine interpolate_gocad_block_MR
-
Deleted: seismo/3D/SPECFEM3D/trunk/mesh_vertical.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/mesh_vertical.f90 2010-09-24 22:15:18 UTC (rev 17214)
+++ seismo/3D/SPECFEM3D/trunk/mesh_vertical.f90 2010-09-24 22:17:01 UTC (rev 17215)
@@ -1,120 +0,0 @@
-!=====================================================================
-!
-! S p e c f e m 3 D V e r s i o n 1 . 4
-! ---------------------------------------
-!
-! Dimitri Komatitsch and Jeroen Tromp
-! Seismological Laboratory - California Institute of Technology
-! (c) California Institute of Technology September 2006
-!
-! This program is free software; you can redistribute it and/or modify
-! it under the terms of the GNU General Public License as published by
-! the Free Software Foundation; either version 2 of the License, or
-! (at your option) any later version.
-!
-! This program is distributed in the hope that it will be useful,
-! but WITHOUT ANY WARRANTY; without even the implied warranty of
-! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-! GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along
-! with this program; if not, write to the Free Software Foundation, Inc.,
-! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-!
-!=====================================================================
-
- subroutine mesh_vertical(myrank,rn,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)
-
-! create the vertical mesh, honoring the major discontinuities in the model
-
- implicit none
-
- include "constants.h"
-
- integer myrank
- integer NER,NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM,NER_SEDIM
- logical MOHO_MAP_LUPEI
- double precision Z_DEPTH_BLOCK,Z_BASEMENT_SURFACE,Z_DEPTH_MOHO
- double precision rn(0:2*NER)
-
-!! DK DK UGLY modif z_top by Emmanuel Chaljub here
-!! DK DK UGLY modif Manu removed double precision z_top
-
- integer npr,ir
-
- npr = -1
-
-!
-!--- bottom of the mesh (Z_DEPTH_BLOCK) to Moho
-!
- do ir=0,2*NER_BOTTOM_MOHO-1
- npr=npr+1
- rn(npr)=(Z_DEPTH_MOHO-Z_DEPTH_BLOCK)*dble(ir)/dble(2*NER_BOTTOM_MOHO)
- enddo
-
-! do not use d16km when Moho map is honored
- if(MOHO_MAP_LUPEI) then
-
-!
-!--- Moho to modified basement surface
-!
- do ir=0,2*(NER_MOHO_16+NER_16_BASEMENT)-1
- npr=npr+1
- rn(npr)=(Z_DEPTH_MOHO-Z_DEPTH_BLOCK) + (Z_BASEMENT_SURFACE-Z_DEPTH_MOHO)*dble(ir)/dble(2*(NER_MOHO_16+NER_16_BASEMENT))
- enddo
-
- else
-!
-!--- Moho to d16km
-!
- do ir=0,2*NER_MOHO_16-1
- npr=npr+1
- rn(npr)=(Z_DEPTH_MOHO-Z_DEPTH_BLOCK) + (DEPTH_16km_SOCAL-Z_DEPTH_MOHO)*dble(ir)/dble(2*NER_MOHO_16)
- enddo
-!
-!--- d16km to modified basement surface
-!
- do ir=0,2*NER_16_BASEMENT-1
- npr=npr+1
- rn(npr)=(DEPTH_16km_SOCAL-Z_DEPTH_BLOCK) + (Z_BASEMENT_SURFACE-DEPTH_16km_SOCAL)*dble(ir)/dble(2*NER_16_BASEMENT)
- enddo
-
- endif
-
-!
-!--- modified basement surface to surface of model (topography/bathymetry)
-!
-! also create last point exactly at the surface
-! other regions above stop one point below
- do ir=0,2*(NER_BASEMENT_SEDIM+NER_SEDIM) - 0
- npr=npr+1
- rn(npr)=(Z_BASEMENT_SURFACE-Z_DEPTH_BLOCK) + &
-!! DK DK UGLY modif z_top by Emmanuel Chaljub here
-!! DK DK UGLY suppressed Manu's modif and put old code back because better mesh
-!! DK DK UGLY investigate this in detail one day
- (Z_SURFACE-Z_BASEMENT_SURFACE)*dble(ir)/dble(2*(NER_BASEMENT_SEDIM+NER_SEDIM))
-!! DK DK UGLY modif Manu removed (z_top-Z_BASEMENT_SURFACE)*dble(ir)/dble(2*(NER_BASEMENT_SEDIM+NER_SEDIM))
- enddo
-
-! normalize depths
-!! DK DK UGLY modif z_top by Emmanuel Chaljub here
-!! DK DK UGLY suppressed Manu's modif and put old code back because better mesh
-!! DK DK UGLY investigate this in detail one day
-!! DK DK UGLY modif Manu removed rn(:) = rn(:) / (z_top-Z_DEPTH_BLOCK)
-!! DK DK UGLY modif Manu removed
- rn(:) = rn(:) / (Z_SURFACE-Z_DEPTH_BLOCK)
-
-! check that the mesh that has been generated is correct
- if(npr /= 2*NER) call exit_MPI(myrank,'incorrect intervals for model')
-
-! check that vertical spacing makes sense
- do ir=0,2*NER-1
- if(rn(ir+1) < rn(ir)) call exit_MPI(myrank,'incorrect vertical spacing for model')
- enddo
-
- end subroutine mesh_vertical
-
Copied: seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/compute_rho_estimate.f90 (from rev 17213, seismo/3D/SPECFEM3D/trunk/compute_rho_estimate.f90)
===================================================================
--- seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/compute_rho_estimate.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/compute_rho_estimate.f90 2010-09-24 22:17:01 UTC (rev 17215)
@@ -0,0 +1,46 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology September 2006
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+ subroutine compute_rho_estimate(rho,vp)
+
+! compute rho estimate in Gocad block and in Hauksson's model
+! based upon Vp
+
+ implicit none
+
+! include "constants.h"
+ include "constants_gocad.h"
+
+ double precision rho,vp
+
+! scale density - use empirical rule from Christiane
+ rho = 0.33d0 * vp + 1280.d0
+
+! make sure density estimate is reasonable
+ if(rho > DENSITY_MAX) rho = DENSITY_MAX
+ if(rho < DENSITY_MIN) rho = DENSITY_MIN
+
+ end subroutine compute_rho_estimate
+
Copied: seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/define_subregions.f90 (from rev 17213, seismo/3D/SPECFEM3D/trunk/define_subregions.f90)
===================================================================
--- seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/define_subregions.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/define_subregions.f90 2010-09-24 22:17:01 UTC (rev 17215)
@@ -0,0 +1,863 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology September 2006
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+ subroutine 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)
+
+! define shape of elements in current subregion of the mesh
+
+ implicit none
+
+ include "constants.h"
+
+ integer myrank
+ integer ix1,ix2,dix,iy1,iy2,diy,ir1,ir2,dir
+ integer iax,iay,iar
+ integer isubregion,doubling_index
+ integer npx,npy
+
+ integer NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM,NER_SEDIM,NER
+
+ logical USE_REGULAR_MESH
+
+! topology of the elements
+ integer iaddx(NGNOD)
+ integer iaddy(NGNOD)
+ integer iaddz(NGNOD)
+
+! **************
+
+!
+!--- case of a regular mesh
+!
+ if(USE_REGULAR_MESH) then
+
+! use two layers even for a regular mesh, because the algorithm detects the top of the mesh
+! (the "topography") based on one layer of elements with flag IFLAG_ONE_LAYER_TOPOGRAPHY
+ if(isubregion == 2) then
+
+ call usual_hex_nodes(iaddx,iaddy,iaddz)
+
+ iy1=0
+ iy2=npy-2
+ diy=2
+
+ ix1=0
+ ix2=npx-2
+ dix=2
+
+ ir1=0
+ ir2=2*(NER - 2)
+ dir=2
+
+ iax=1
+ iay=1
+ iar=1
+
+ doubling_index = IFLAG_BASEMENT_TOPO
+
+ else if(isubregion == 1) then
+
+ call usual_hex_nodes(iaddx,iaddy,iaddz)
+
+ iy1=0
+ iy2=npy-2
+ diy=2
+
+ ix1=0
+ ix2=npx-2
+ dix=2
+
+ ir1=2*(NER - 1)
+ ir2=ir1
+ dir=2
+
+ iax=1
+ iay=1
+ iar=1
+
+ doubling_index = IFLAG_ONE_LAYER_TOPOGRAPHY
+
+ else
+
+ call exit_MPI(myrank,'incorrect subregion code')
+
+ endif
+
+!
+!--- case of a non-regular mesh with mesh doublings
+!
+ else
+
+! this last region only defined when NER_SEDIM > 1
+ if(isubregion == 30) then
+
+ call usual_hex_nodes(iaddx,iaddy,iaddz)
+
+ iy1=0
+ iy2=npy-2
+ diy=2
+
+ ix1=0
+ ix2=npx-2
+ dix=2
+
+ ir1=2*(NER - NER_SEDIM)
+ ir2=2*(NER - 2)
+ dir=2
+
+ iax=1
+ iay=1
+ iar=1
+
+ doubling_index = IFLAG_BASEMENT_TOPO
+
+ else if(isubregion == 29) then
+
+ call usual_hex_nodes(iaddx,iaddy,iaddz)
+
+ iy1=0
+ iy2=npy-2
+ diy=2
+
+ ix1=0
+ ix2=npx-2
+ dix=2
+
+ ir1=2*(NER - 1)
+ ir2=ir1
+ dir=2
+
+ iax=1
+ iay=1
+ iar=1
+
+ doubling_index = IFLAG_ONE_LAYER_TOPOGRAPHY
+
+ else if(isubregion == 28) then
+
+ call usual_hex_nodes(iaddx,iaddy,iaddz)
+
+ iy1=0
+ iy2=npy-8
+ diy=8
+
+ ix1=0
+ ix2=npx-8
+ dix=8
+
+ ir1= 0
+ ir2= 2*NER_BOTTOM_MOHO-8
+ dir=8
+
+ iax=4
+ iay=4
+ iar=4
+
+ doubling_index= IFLAG_HALFSPACE_MOHO
+
+ else if(isubregion == 27) then
+
+ call unusual_hex_nodes1(iaddx,iaddy,iaddz)
+
+! generating stage 2 of the mesh doubling below 670
+
+ iy1=0
+ iy2=npy-8
+ diy=8
+
+ ix1=0
+ ix2=npx-16
+ dix=16
+
+ dir=4
+
+ iax=2
+ iay=2
+ iar=1
+
+ ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 8
+ ir2=ir1
+
+ doubling_index=IFLAG_16km_BASEMENT
+
+ else if(isubregion == 26) then
+
+ call unusual_hex_nodes1p(iaddx,iaddy,iaddz)
+
+! generating stage 3 of the mesh doubling below 670
+
+ iy1=0
+ iy2=npy-8
+ diy=8
+
+ ix1=8
+ ix2=npx-8
+ dix=16
+
+ dir=4
+
+ iax=2
+ iay=2
+ iar=1
+
+ ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 8
+ ir2=ir1
+
+ doubling_index=IFLAG_16km_BASEMENT
+
+ else if(isubregion == 25) then
+
+ call unusual_hex_nodes2(iaddx,iaddy,iaddz)
+
+! generating stage 4 of the mesh doubling below 670
+
+ iy1=0
+ iy2=npy-8
+ diy=8
+
+ ix1=0
+ ix2=npx-16
+ dix=16
+
+ dir=4
+
+ iax=2
+ iay=2
+ iar=1
+
+ ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 8
+ ir2=ir1
+
+ doubling_index=IFLAG_16km_BASEMENT
+
+ else if(isubregion == 24) then
+
+ call unusual_hex_nodes2p(iaddx,iaddy,iaddz)
+
+! generating stage 5 of the mesh doubling below 670
+
+ iy1=0
+ iy2=npy-8
+ diy=8
+
+ ix1=12
+ ix2=npx-4
+ dix=16
+
+ dir=4
+
+ iax=2
+ iay=2
+ iar=1
+
+ ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 6
+ ir2=ir1
+
+ doubling_index=IFLAG_16km_BASEMENT
+
+ else if(isubregion == 23) then
+
+ call unusual_hex_nodes3(iaddx,iaddy,iaddz)
+
+! generating stage 6 of the mesh doubling below 670
+
+ iy1=0
+ iy2=npy-8
+ diy=8
+
+ ix1=4
+ ix2=npx-12
+ dix=16
+
+ dir=4
+
+ iax=2
+ iay=2
+ iar=1
+
+ ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 6
+ ir2=ir1
+
+ doubling_index=IFLAG_16km_BASEMENT
+
+ else if(isubregion == 22) then
+
+ call unusual_hex_nodes3(iaddx,iaddy,iaddz)
+
+! generating stage 7 of the mesh doubling below 670
+
+ iy1=0
+ iy2=npy-8
+ diy=8
+
+ ix1=8
+ ix2=npx-8
+ dix=16
+
+ dir=4
+
+ iax=2
+ iay=2
+ iar=1
+
+ ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 6
+ ir2=ir1
+
+ doubling_index=IFLAG_16km_BASEMENT
+
+ else if(isubregion == 21) then
+
+ call unusual_hex_nodes4(iaddx,iaddy,iaddz)
+
+! generating stage 8 of the mesh doubling below 670
+
+ iy1=8
+ iy2=npy-8
+ diy=16
+
+ ix1=0
+ ix2=npx-4
+ dix=4
+
+ dir=4
+
+ iax=2
+ iay=2
+ iar=1
+
+ ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 4
+ ir2=ir1
+
+ doubling_index=IFLAG_16km_BASEMENT
+
+ else if(isubregion == 20) then
+
+ call unusual_hex_nodes4p(iaddx,iaddy,iaddz)
+
+! generating stage 9 of the mesh doubling below 670
+
+ iy1=0
+ iy2=npy-16
+ diy=16
+
+ ix1=0
+ ix2=npx-4
+ dix=4
+
+ dir=4
+
+ iax=2
+ iay=2
+ iar=1
+
+ ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 4
+ ir2=ir1
+
+ doubling_index=IFLAG_16km_BASEMENT
+
+ else if(isubregion == 19) then
+
+ call usual_hex_nodes(iaddx,iaddy,iaddz)
+
+! generating stage 10 of the mesh doubling below 670
+
+ iy1=8
+ iy2=npy-8
+ diy=16
+
+ ix1=0
+ ix2=npx-4
+ dix=4
+
+ dir=4
+
+ iax=2
+ iay=2
+ iar=1
+
+ ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 2
+ ir2=ir1
+
+ doubling_index=IFLAG_16km_BASEMENT
+
+ else if(isubregion == 18) then
+
+ call usual_hex_nodes(iaddx,iaddy,iaddz)
+
+! generating stage 11 of the mesh doubling below 670
+
+ iy1=4
+ iy2=npy-12
+ diy=16
+
+ ix1=0
+ ix2=npx-4
+ dix=4
+
+ dir=4
+
+ iax=2
+ iay=2
+ iar=1
+
+ ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 2
+ ir2=ir1
+
+ doubling_index=IFLAG_16km_BASEMENT
+
+ else if(isubregion == 17) then
+
+ call unusual_hex_nodes6(iaddx,iaddy,iaddz)
+
+! generating stage 12 of the mesh doubling below 670
+
+ iy1=12
+ iy2=npy-4
+ diy=16
+
+ ix1=0
+ ix2=npx-4
+ dix=4
+
+ dir=4
+
+ iax=2
+ iay=2
+ iar=1
+
+ ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 2
+ ir2=ir1
+
+ doubling_index=IFLAG_16km_BASEMENT
+
+ else if(isubregion == 16) then
+
+ call unusual_hex_nodes6p(iaddx,iaddy,iaddz)
+
+! generating stage 13 of the mesh doubling below 670
+
+ iy1=0
+ iy2=npy-16
+ diy=16
+
+ ix1=0
+ ix2=npx-4
+ dix=4
+
+ dir=4
+
+ iax=2
+ iay=2
+ iar=1
+
+ ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 4
+ ir2=ir1
+
+ doubling_index=IFLAG_16km_BASEMENT
+
+ else if(isubregion == 15) then
+
+ call usual_hex_nodes(iaddx,iaddy,iaddz)
+
+ iy1=0
+ iy2=npy-8
+ diy=8
+
+ ix1=0
+ ix2=npx-8
+ dix=8
+
+! honor So-Cal model discontinuity at 16 km for accuracy
+ ir1=2*NER_BOTTOM_MOHO
+ ir2=2*(NER_BOTTOM_MOHO+NER_MOHO_16) - 4
+ dir=4
+
+ iax=4
+ iay=4
+ iar=2
+
+ doubling_index = IFLAG_MOHO_16km
+
+ else if(isubregion == 14) then
+
+ call usual_hex_nodes(iaddx,iaddy,iaddz)
+
+ iy1=0
+ iy2=npy-8
+ diy=8
+
+ ix1=0
+ ix2=npx-8
+ dix=8
+
+! honor So-Cal model discontinuity at 16 km for accuracy
+ ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16)
+ ir2=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT)-12
+ dir=4
+
+ iax=4
+ iay=4
+ iar=2
+
+ doubling_index = IFLAG_16km_BASEMENT
+
+
+ else if(isubregion == 13) then
+
+ call usual_hex_nodes(iaddx,iaddy,iaddz)
+
+! generating stage 1 of the mesh doubling below the Moho
+
+ iy1=0
+ iy2=npy-4
+ diy=4
+
+ ix1=0
+ ix2=npx-4
+ dix=4
+
+ ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT)
+ ir2=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-12
+ dir=4
+
+ iax=2
+ iay=2
+ iar=2
+
+ doubling_index=IFLAG_BASEMENT_TOPO
+
+ else if(isubregion == 12) then
+
+ call unusual_hex_nodes1(iaddx,iaddy,iaddz)
+
+! generating stage 2 of the mesh doubling below the Moho
+
+ iy1=0
+ iy2=npy-4
+ diy=4
+
+ ix1=0
+ ix2=npx-8
+ dix=8
+
+ dir=4
+
+ iax=1
+ iay=1
+ iar=1
+
+ ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-8
+ ir2=ir1
+
+ doubling_index=IFLAG_BASEMENT_TOPO
+
+ else if(isubregion == 11) then
+
+ call unusual_hex_nodes1p(iaddx,iaddy,iaddz)
+
+! generating stage 3 of the mesh doubling below the Moho
+
+ iy1=0
+ iy2=npy-4
+ diy=4
+
+ ix1=4
+ ix2=npx-4
+ dix=8
+
+ dir=4
+
+ iax=1
+ iay=1
+ iar=1
+
+ ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-8
+ ir2=ir1
+
+ doubling_index=IFLAG_BASEMENT_TOPO
+
+ else if(isubregion == 10) then
+
+ call unusual_hex_nodes2(iaddx,iaddy,iaddz)
+
+! generating stage 4 of the mesh doubling below the Moho
+
+ iy1=0
+ iy2=npy-4
+ diy=4
+
+ ix1=0
+ ix2=npx-8
+ dix=8
+
+ dir=4
+
+ iax=1
+ iay=1
+ iar=1
+
+ ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-8
+ ir2=ir1
+
+ doubling_index=IFLAG_BASEMENT_TOPO
+
+ else if(isubregion == 9) then
+
+ call unusual_hex_nodes2p(iaddx,iaddy,iaddz)
+
+! generating stage 5 of the mesh doubling below the Moho
+
+ iy1=0
+ iy2=npy-4
+ diy=4
+
+ ix1=6
+ ix2=npx-2
+ dix=8
+
+ dir=4
+
+ iax=1
+ iay=1
+ iar=1
+
+ ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-6
+ ir2=ir1
+
+ doubling_index=IFLAG_BASEMENT_TOPO
+
+ else if(isubregion == 8) then
+
+ call unusual_hex_nodes3(iaddx,iaddy,iaddz)
+
+! generating stage 6 of the mesh doubling below the Moho
+
+ iy1=0
+ iy2=npy-4
+ diy=4
+
+ ix1=2
+ ix2=npx-6
+ dix=8
+
+ dir=4
+
+ iax=1
+ iay=1
+ iar=1
+
+ ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-6
+ ir2=ir1
+
+ doubling_index=IFLAG_BASEMENT_TOPO
+
+ else if(isubregion == 7) then
+
+ call unusual_hex_nodes3(iaddx,iaddy,iaddz)
+
+! generating stage 7 of the mesh doubling below the Moho
+
+ iy1=0
+ iy2=npy-4
+ diy=4
+
+ ix1=4
+ ix2=npx-4
+ dix=8
+
+ dir=4
+
+ iax=1
+ iay=1
+ iar=1
+
+ ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-6
+ ir2=ir1
+
+ doubling_index=IFLAG_BASEMENT_TOPO
+
+ else if(isubregion == 6) then
+
+ call unusual_hex_nodes4(iaddx,iaddy,iaddz)
+
+! generating stage 8 of the mesh doubling below the Moho
+
+ iy1=4
+ iy2=npy-4
+ diy=8
+
+ ix1=0
+ ix2=npx-2
+ dix=2
+
+ dir=4
+
+ iax=1
+ iay=1
+ iar=1
+
+ ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-4
+ ir2=ir1
+
+ doubling_index=IFLAG_BASEMENT_TOPO
+
+ else if(isubregion == 5) then
+
+ call unusual_hex_nodes4p(iaddx,iaddy,iaddz)
+
+! generating stage 9 of the mesh doubling below the Moho
+
+ iy1=0
+ iy2=npy-8
+ diy=8
+
+ ix1=0
+ ix2=npx-2
+ dix=2
+
+ dir=4
+
+ iax=1
+ iay=1
+ iar=1
+
+ ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-4
+ ir2=ir1
+
+ doubling_index=IFLAG_BASEMENT_TOPO
+
+ else if(isubregion == 4) then
+
+ call usual_hex_nodes(iaddx,iaddy,iaddz)
+
+! generating stage 10 of the mesh doubling below the Moho
+
+ iy1=4
+ iy2=npy-4
+ diy=8
+
+ ix1=0
+ ix2=npx-2
+ dix=2
+
+ dir=4
+
+ iax=1
+ iay=1
+ iar=1
+
+ ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-2
+ ir2=ir1
+
+ doubling_index=IFLAG_BASEMENT_TOPO
+
+ else if(isubregion == 3) then
+
+ call usual_hex_nodes(iaddx,iaddy,iaddz)
+
+! generating stage 11 of the mesh doubling below the Moho
+
+ iy1=2
+ iy2=npy-6
+ diy=8
+
+ ix1=0
+ ix2=npx-2
+ dix=2
+
+ dir=4
+
+ iax=1
+ iay=1
+ iar=1
+
+ ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-2
+ ir2=ir1
+
+ doubling_index=IFLAG_BASEMENT_TOPO
+
+ else if(isubregion == 2) then
+
+ call unusual_hex_nodes6(iaddx,iaddy,iaddz)
+
+! generating stage 12 of the mesh doubling below the Moho
+
+ iy1=6
+ iy2=npy-2
+ diy=8
+
+ ix1=0
+ ix2=npx-2
+ dix=2
+
+ dir=4
+
+ iax=1
+ iay=1
+ iar=1
+
+ ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-2
+ ir2=ir1
+
+ doubling_index=IFLAG_BASEMENT_TOPO
+
+ else if(isubregion == 1) then
+
+ call unusual_hex_nodes6p(iaddx,iaddy,iaddz)
+
+! generating stage 13 of the mesh doubling below the Moho
+
+ iy1=0
+ iy2=npy-8
+ diy=8
+
+ ix1=0
+ ix2=npx-2
+ dix=2
+
+ dir=4
+
+ iax=1
+ iay=1
+ iar=1
+
+ ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-4
+ ir2=ir1
+
+ doubling_index=IFLAG_BASEMENT_TOPO
+
+ else
+
+ call exit_MPI(myrank,'incorrect subregion code')
+
+ endif
+
+ endif
+
+ end subroutine define_subregions
+
Copied: seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/define_subregions_heuristic.f90 (from rev 17213, seismo/3D/SPECFEM3D/trunk/define_subregions_heuristic.f90)
===================================================================
--- seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/define_subregions_heuristic.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/define_subregions_heuristic.f90 2010-09-24 22:17:01 UTC (rev 17215)
@@ -0,0 +1,267 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology September 2006
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+ subroutine 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)
+
+! heuristic rule to deform elements to balance angles
+! to 120 degrees in doubling regions
+
+ implicit none
+
+ include "constants.h"
+
+ integer myrank
+ integer ix1,ix2,dix,iy1,iy2,diy,ir1,ir2,dir
+ integer iax,iay,iar
+ integer isubregion,itype_element
+ integer npx,npy
+
+ integer NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM
+
+! topology of the elements
+ integer iaddx(NGNOD)
+ integer iaddy(NGNOD)
+ integer iaddz(NGNOD)
+
+! type of elements for heuristic rule
+ integer, parameter :: ITYPE_UNUSUAL_1 = 1
+ integer, parameter :: ITYPE_UNUSUAL_1p = 2
+ integer, parameter :: ITYPE_UNUSUAL_4 = 3
+ integer, parameter :: ITYPE_UNUSUAL_4p = 4
+
+
+! **************
+
+ if(isubregion == 8) then
+
+ call unusual_hex_nodes1(iaddx,iaddy,iaddz)
+
+! generating stage 2 of the mesh doubling below 670
+
+ iy1=0
+ iy2=npy-8
+ diy=8
+
+ ix1=0
+ ix2=npx-16
+ dix=16
+
+ dir=4
+
+ iax=2
+ iay=2
+ iar=1
+
+ ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 8
+ ir2=ir1
+
+ itype_element = ITYPE_UNUSUAL_1
+
+ else if(isubregion == 7) then
+
+ call unusual_hex_nodes1p(iaddx,iaddy,iaddz)
+
+! generating stage 3 of the mesh doubling below 670
+
+ iy1=0
+ iy2=npy-8
+ diy=8
+
+ ix1=8
+ ix2=npx-8
+ dix=16
+
+ dir=4
+
+ iax=2
+ iay=2
+ iar=1
+
+ ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 8
+ ir2=ir1
+
+ itype_element = ITYPE_UNUSUAL_1p
+
+ else if(isubregion == 6) then
+
+ call unusual_hex_nodes4(iaddx,iaddy,iaddz)
+
+! generating stage 8 of the mesh doubling below 670
+
+ iy1=8
+ iy2=npy-8
+ diy=16
+
+ ix1=0
+ ix2=npx-4
+ dix=4
+
+ dir=4
+
+ iax=2
+ iay=2
+ iar=1
+
+ ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 4
+ ir2=ir1
+
+ itype_element = ITYPE_UNUSUAL_4
+
+ else if(isubregion == 5) then
+
+ call unusual_hex_nodes4p(iaddx,iaddy,iaddz)
+
+! generating stage 9 of the mesh doubling below 670
+
+ iy1=0
+ iy2=npy-16
+ diy=16
+
+ ix1=0
+ ix2=npx-4
+ dix=4
+
+ dir=4
+
+ iax=2
+ iay=2
+ iar=1
+
+ ir1=2*(NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT) - 4
+ ir2=ir1
+
+ itype_element = ITYPE_UNUSUAL_4p
+
+ else if(isubregion == 4) then
+
+ call unusual_hex_nodes1(iaddx,iaddy,iaddz)
+
+! generating stage 2 of the mesh doubling below the Moho
+
+ iy1=0
+ iy2=npy-4
+ diy=4
+
+ ix1=0
+ ix2=npx-8
+ dix=8
+
+ dir=4
+
+ iax=1
+ iay=1
+ iar=1
+
+ ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-8
+ ir2=ir1
+
+ itype_element = ITYPE_UNUSUAL_1
+
+ else if(isubregion == 3) then
+
+ call unusual_hex_nodes1p(iaddx,iaddy,iaddz)
+
+! generating stage 3 of the mesh doubling below the Moho
+
+ iy1=0
+ iy2=npy-4
+ diy=4
+
+ ix1=4
+ ix2=npx-4
+ dix=8
+
+ dir=4
+
+ iax=1
+ iay=1
+ iar=1
+
+ ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-8
+ ir2=ir1
+
+ itype_element = ITYPE_UNUSUAL_1p
+
+ else if(isubregion == 2) then
+
+ call unusual_hex_nodes4(iaddx,iaddy,iaddz)
+
+! generating stage 8 of the mesh doubling below the Moho
+
+ iy1=4
+ iy2=npy-4
+ diy=8
+
+ ix1=0
+ ix2=npx-2
+ dix=2
+
+ dir=4
+
+ iax=1
+ iay=1
+ iar=1
+
+ ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-4
+ ir2=ir1
+
+ itype_element = ITYPE_UNUSUAL_4
+
+ else if(isubregion == 1) then
+
+ call unusual_hex_nodes4p(iaddx,iaddy,iaddz)
+
+! generating stage 9 of the mesh doubling below the Moho
+
+ iy1=0
+ iy2=npy-8
+ diy=8
+
+ ix1=0
+ ix2=npx-2
+ dix=2
+
+ dir=4
+
+ iax=1
+ iay=1
+ iar=1
+
+ ir1=2*(NER_BOTTOM_MOHO+NER_MOHO_16+NER_16_BASEMENT+NER_BASEMENT_SEDIM)-4
+ ir2=ir1
+
+ itype_element = ITYPE_UNUSUAL_4p
+
+ else
+
+ call exit_MPI(myrank,'incorrect subregion code')
+
+ endif
+
+ end subroutine define_subregions_heuristic
+
Copied: seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/get_absorb.f90 (from rev 17213, seismo/3D/SPECFEM3D/trunk/get_absorb.f90)
===================================================================
--- seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/get_absorb.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/get_absorb.f90 2010-09-24 22:17:01 UTC (rev 17215)
@@ -0,0 +1,270 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology September 2006
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+ subroutine get_absorb(myrank,prname,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
+
+! processor identification
+ character(len=256) prname
+
+ 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')
+
+! IMPROVE save these temporary arrays for the solver for Stacey conditions
+
+ open(unit=27,file=prname(1:len_trim(prname))//'nimin.bin',status='unknown',form='unformatted')
+ write(27) nimin
+ close(27)
+
+ open(unit=27,file=prname(1:len_trim(prname))//'nimax.bin',status='unknown',form='unformatted')
+ write(27) nimax
+ close(27)
+
+ open(unit=27,file=prname(1:len_trim(prname))//'njmin.bin',status='unknown',form='unformatted')
+ write(27) njmin
+ close(27)
+
+ open(unit=27,file=prname(1:len_trim(prname))//'njmax.bin',status='unknown',form='unformatted')
+ write(27) njmax
+ close(27)
+
+ open(unit=27,file=prname(1:len_trim(prname))//'nkmin_xi.bin',status='unknown',form='unformatted')
+ write(27) nkmin_xi
+ close(27)
+
+ open(unit=27,file=prname(1:len_trim(prname))//'nkmin_eta.bin',status='unknown',form='unformatted')
+ write(27) nkmin_eta
+ close(27)
+
+ 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
+
+
+
+
Copied: seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/get_flags_boundaries.f90 (from rev 17213, seismo/3D/SPECFEM3D/trunk/get_flags_boundaries.f90)
===================================================================
--- seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/get_flags_boundaries.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/get_flags_boundaries.f90 2010-09-24 22:17:01 UTC (rev 17215)
@@ -0,0 +1,162 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology September 2006
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+ subroutine get_flags_boundaries(nspec,iproc_xi,iproc_eta,ispec,idoubling, &
+ xstore,ystore,zstore,iboun,iMPIcut_xi,iMPIcut_eta, &
+ NPROC_XI,NPROC_ETA, &
+ UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK)
+
+ implicit none
+
+ include "constants.h"
+ include "constants_gocad.h"
+
+ integer nspec
+ integer ispec,idoubling
+ integer NPROC_XI,NPROC_ETA
+
+ double precision UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK
+
+ logical iboun(6,nspec)
+ logical iMPIcut_xi(2,nspec),iMPIcut_eta(2,nspec)
+
+ double precision xstore(NGLLX,NGLLY,NGLLZ)
+ double precision ystore(NGLLX,NGLLY,NGLLZ)
+ double precision zstore(NGLLX,NGLLY,NGLLZ)
+
+! use iproc_xi and iproc_eta to determine MPI cut planes along xi and eta
+ integer iproc_xi,iproc_eta
+
+ double precision target,sizeslice,TOLERANCE_METERS
+ double precision xelm(8),yelm(8),zelm(8)
+
+! find the coordinates of the eight corner nodes of the element
+ xelm(1)=xstore(1,1,1)
+ yelm(1)=ystore(1,1,1)
+ zelm(1)=zstore(1,1,1)
+ xelm(2)=xstore(NGLLX,1,1)
+ yelm(2)=ystore(NGLLX,1,1)
+ zelm(2)=zstore(NGLLX,1,1)
+ xelm(3)=xstore(NGLLX,NGLLY,1)
+ yelm(3)=ystore(NGLLX,NGLLY,1)
+ zelm(3)=zstore(NGLLX,NGLLY,1)
+ xelm(4)=xstore(1,NGLLY,1)
+ yelm(4)=ystore(1,NGLLY,1)
+ zelm(4)=zstore(1,NGLLY,1)
+ xelm(5)=xstore(1,1,NGLLZ)
+ yelm(5)=ystore(1,1,NGLLZ)
+ zelm(5)=zstore(1,1,NGLLZ)
+ xelm(6)=xstore(NGLLX,1,NGLLZ)
+ yelm(6)=ystore(NGLLX,1,NGLLZ)
+ zelm(6)=zstore(NGLLX,1,NGLLZ)
+ xelm(7)=xstore(NGLLX,NGLLY,NGLLZ)
+ yelm(7)=ystore(NGLLX,NGLLY,NGLLZ)
+ zelm(7)=zstore(NGLLX,NGLLY,NGLLZ)
+ xelm(8)=xstore(1,NGLLY,NGLLZ)
+ yelm(8)=ystore(1,NGLLY,NGLLZ)
+ zelm(8)=zstore(1,NGLLY,NGLLZ)
+
+! compute geometrical tolerance small compared to size of model to detect edges
+ TOLERANCE_METERS = dabs(UTM_X_MAX - UTM_X_MIN) / 100000.
+
+! ****************************************************
+! determine if the element falls on a boundary
+! ****************************************************
+
+ iboun(:,ispec)=.false.
+
+! on boundary 1: x=xmin
+ target= UTM_X_MIN + TOLERANCE_METERS
+ if(xelm(1)<target .and. xelm(4)<target .and. xelm(5)<target .and. xelm(8)<target) iboun(1,ispec)=.true.
+
+! on boundary 2: xmax
+ target= UTM_X_MAX - TOLERANCE_METERS
+ if(xelm(2)>target .and. xelm(3)>target .and. xelm(6)>target .and. xelm(7)>target) iboun(2,ispec)=.true.
+
+! on boundary 3: ymin
+ target= UTM_Y_MIN + TOLERANCE_METERS
+ if(yelm(1)<target .and. yelm(2)<target .and. yelm(5)<target .and. yelm(6)<target) iboun(3,ispec)=.true.
+
+! on boundary 4: ymax
+ target= UTM_Y_MAX - TOLERANCE_METERS
+ if(yelm(3)>target .and. yelm(4)>target .and. yelm(7)>target .and. yelm(8)>target) iboun(4,ispec)=.true.
+
+! on boundary 5: bottom
+ target = Z_DEPTH_BLOCK + TOLERANCE_METERS
+ if(zelm(1)<target .and. zelm(2)<target .and. zelm(3)<target .and. zelm(4)<target) iboun(5,ispec)=.true.
+
+! on boundary 6: top
+ if(idoubling == IFLAG_ONE_LAYER_TOPOGRAPHY) iboun(6,ispec)=.true.
+
+! *******************************************************************
+! determine if the element falls on an MPI cut plane along xi
+! *******************************************************************
+
+! detect the MPI cut planes along xi in the cubed sphere
+
+ iMPIcut_xi(:,ispec)=.false.
+
+! angular size of a slice along xi
+ sizeslice = (UTM_X_MAX-UTM_X_MIN) / NPROC_XI
+
+! left cut-plane in the current slice along X = constant (Xmin of this slice)
+! and add geometrical tolerance
+
+ target = UTM_X_MIN + iproc_xi*sizeslice + TOLERANCE_METERS
+ if(xelm(1)<target .and. xelm(4)<target .and. xelm(5)<target .and. xelm(8)<target) &
+ iMPIcut_xi(1,ispec)=.true.
+
+! right cut-plane in the current slice along X = constant (Xmax of this slice)
+! and add geometrical tolerance
+
+ target = UTM_X_MIN + (iproc_xi+1)*sizeslice - TOLERANCE_METERS
+ if(xelm(2)>target .and. xelm(3)>target .and. xelm(6)>target .and. xelm(7)>target) &
+ iMPIcut_xi(2,ispec)=.true.
+
+! ********************************************************************
+! determine if the element falls on an MPI cut plane along eta
+! ********************************************************************
+
+ iMPIcut_eta(:,ispec)=.false.
+
+! angular size of a slice along eta
+ sizeslice = (UTM_Y_MAX-UTM_Y_MIN) / NPROC_ETA
+
+! left cut-plane in the current slice along Y = constant (Ymin of this slice)
+! and add geometrical tolerance
+
+ target = UTM_Y_MIN + iproc_eta*sizeslice + TOLERANCE_METERS
+ if(yelm(1)<target .and. yelm(2)<target .and. yelm(5)<target .and. yelm(6)<target) &
+ iMPIcut_eta(1,ispec)=.true.
+
+! right cut-plane in the current slice along Y = constant (Ymax of this slice)
+! and add geometrical tolerance
+
+ target = UTM_Y_MIN + (iproc_eta+1)*sizeslice - TOLERANCE_METERS
+ if(yelm(3)>target .and. yelm(4)>target .and. yelm(7)>target .and. yelm(8)>target) &
+ iMPIcut_eta(2,ispec)=.true.
+
+ end subroutine get_flags_boundaries
+
Copied: seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/hauksson_model.f90 (from rev 17213, seismo/3D/SPECFEM3D/trunk/hauksson_model.f90)
===================================================================
--- seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/hauksson_model.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/hauksson_model.f90 2010-09-24 22:17:01 UTC (rev 17215)
@@ -0,0 +1,227 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology September 2006
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+ subroutine hauksson_model(vp,vs,utm_x_eval,utm_y_eval,z_eval,vp_final,vs_final,MOHO_MAP_LUPEI)
+
+ implicit none
+
+ include "constants.h"
+ include "constants_gocad.h"
+
+!! DK DK UGLY one day, we should clarify the issue of merging Hauksson's Moho
+!! DK DK UGLY with our Lupei Moho. Should not be a big issue because in
+!! DK DK UGLY principle Hauksson used Lupei's map to build his model
+
+ double precision utm_x_eval,utm_y_eval,z_eval
+ double precision vp_final,vs_final
+ logical MOHO_MAP_LUPEI
+
+ double precision, dimension(NLAYERS_HAUKSSON,NGRID_NEW_HAUKSSON,NGRID_NEW_HAUKSSON) :: vp,vs
+ double precision, dimension(NLAYERS_HAUKSSON) :: vp_interp,vs_interp
+
+ integer ilayer
+ integer icell_interp_x,icell_interp_y
+ double precision spacing_x,spacing_y
+ double precision utm_x_eval_copy,utm_y_eval_copy
+ double precision gamma_interp_x,gamma_interp_y,gamma_interp_z
+ double precision v1,v2,v3,v4
+ double precision vp_upper,vs_upper,vp_lower,vs_lower,z_upper,z_lower
+
+! copy input values
+ utm_x_eval_copy = utm_x_eval
+ utm_y_eval_copy = utm_y_eval
+
+! make sure we stay inside Hauksson's grid
+ if(utm_x_eval_copy < UTM_X_ORIG_HAUKSSON) utm_x_eval_copy = UTM_X_ORIG_HAUKSSON
+ if(utm_y_eval_copy < UTM_Y_ORIG_HAUKSSON) utm_y_eval_copy = UTM_Y_ORIG_HAUKSSON
+
+! determine spacing and cell for linear interpolation
+ spacing_x = (utm_x_eval_copy - UTM_X_ORIG_HAUKSSON) / SPACING_UTM_X_HAUKSSON
+ spacing_y = (utm_y_eval_copy - UTM_Y_ORIG_HAUKSSON) / SPACING_UTM_Y_HAUKSSON
+
+ icell_interp_x = int(spacing_x) + 1
+ icell_interp_y = int(spacing_y) + 1
+
+ gamma_interp_x = spacing_x - int(spacing_x)
+ gamma_interp_y = spacing_y - int(spacing_y)
+
+! suppress edge effects for points outside of Hauksson's model
+ if(icell_interp_x < 1) then
+ icell_interp_x = 1
+ gamma_interp_x = 0.d0
+ endif
+ if(icell_interp_x > NGRID_NEW_HAUKSSON-1) then
+ icell_interp_x = NGRID_NEW_HAUKSSON-1
+ gamma_interp_x = 1.d0
+ endif
+
+ if(icell_interp_y < 1) then
+ icell_interp_y = 1
+ gamma_interp_y = 0.d0
+ endif
+ if(icell_interp_y > NGRID_NEW_HAUKSSON-1) then
+ icell_interp_y = NGRID_NEW_HAUKSSON-1
+ gamma_interp_y = 1.d0
+ endif
+
+! make sure interpolation makes sense
+ if(gamma_interp_x < -0.001d0 .or. gamma_interp_x > 1.001d0) &
+ stop 'interpolation in x is incorrect in Hauksson'
+ if(gamma_interp_y < -0.001d0 .or. gamma_interp_y > 1.001d0) &
+ stop 'interpolation in y is incorrect in Hauksson'
+
+! interpolate Hauksson's model at right location using bilinear interpolation
+ do ilayer = 1,NLAYERS_HAUKSSON
+
+! for Vp
+ v1 = vp(ilayer,icell_interp_x,icell_interp_y)
+ v2 = vp(ilayer,icell_interp_x+1,icell_interp_y)
+ v3 = vp(ilayer,icell_interp_x+1,icell_interp_y+1)
+ v4 = vp(ilayer,icell_interp_x,icell_interp_y+1)
+
+ vp_interp(ilayer) = v1*(1.-gamma_interp_x)*(1.-gamma_interp_y) + &
+ v2*gamma_interp_x*(1.-gamma_interp_y) + &
+ v3*gamma_interp_x*gamma_interp_y + &
+ v4*(1.-gamma_interp_x)*gamma_interp_y
+
+! for Vs
+ v1 = vs(ilayer,icell_interp_x,icell_interp_y)
+ v2 = vs(ilayer,icell_interp_x+1,icell_interp_y)
+ v3 = vs(ilayer,icell_interp_x+1,icell_interp_y+1)
+ v4 = vs(ilayer,icell_interp_x,icell_interp_y+1)
+
+ vs_interp(ilayer) = v1*(1.-gamma_interp_x)*(1.-gamma_interp_y) + &
+ v2*gamma_interp_x*(1.-gamma_interp_y) + &
+ v3*gamma_interp_x*gamma_interp_y + &
+ v4*(1.-gamma_interp_x)*gamma_interp_y
+
+ enddo
+
+! choose right values depending on depth of target point
+ if(z_eval >= Z_HAUKSSON_LAYER_1) then
+ vp_final = vp_interp(1)
+ vs_final = vs_interp(1)
+ return
+
+ else if(z_eval <= Z_HAUKSSON_LAYER_9) then
+ vp_final = vp_interp(9)
+ vs_final = vs_interp(9)
+ return
+
+ else if(z_eval >= Z_HAUKSSON_LAYER_2) then
+ vp_upper = vp_interp(1)
+ vs_upper = vs_interp(1)
+ z_upper = Z_HAUKSSON_LAYER_1
+
+ vp_lower = vp_interp(2)
+ vs_lower = vs_interp(2)
+ z_lower = Z_HAUKSSON_LAYER_2
+
+ else if(z_eval >= Z_HAUKSSON_LAYER_3) then
+ vp_upper = vp_interp(2)
+ vs_upper = vs_interp(2)
+ z_upper = Z_HAUKSSON_LAYER_2
+
+ vp_lower = vp_interp(3)
+ vs_lower = vs_interp(3)
+ z_lower = Z_HAUKSSON_LAYER_3
+
+ else if(z_eval >= Z_HAUKSSON_LAYER_4) then
+ vp_upper = vp_interp(3)
+ vs_upper = vs_interp(3)
+ z_upper = Z_HAUKSSON_LAYER_3
+
+ vp_lower = vp_interp(4)
+ vs_lower = vs_interp(4)
+ z_lower = Z_HAUKSSON_LAYER_4
+
+ else if(z_eval >= Z_HAUKSSON_LAYER_5) then
+ vp_upper = vp_interp(4)
+ vs_upper = vs_interp(4)
+ z_upper = Z_HAUKSSON_LAYER_4
+
+ vp_lower = vp_interp(5)
+ vs_lower = vs_interp(5)
+ z_lower = Z_HAUKSSON_LAYER_5
+
+ else if(z_eval >= Z_HAUKSSON_LAYER_6) then
+ vp_upper = vp_interp(5)
+ vs_upper = vs_interp(5)
+ z_upper = Z_HAUKSSON_LAYER_5
+
+ vp_lower = vp_interp(6)
+ vs_lower = vs_interp(6)
+ z_lower = Z_HAUKSSON_LAYER_6
+
+ else if(z_eval >= Z_HAUKSSON_LAYER_7) then
+ vp_upper = vp_interp(6)
+ vs_upper = vs_interp(6)
+ z_upper = Z_HAUKSSON_LAYER_6
+
+ vp_lower = vp_interp(7)
+ vs_lower = vs_interp(7)
+ z_lower = Z_HAUKSSON_LAYER_7
+
+ else if(z_eval >= Z_HAUKSSON_LAYER_8) then
+ vp_upper = vp_interp(7)
+ vs_upper = vs_interp(7)
+ z_upper = Z_HAUKSSON_LAYER_7
+
+ vp_lower = vp_interp(8)
+ vs_lower = vs_interp(8)
+ z_lower = Z_HAUKSSON_LAYER_8
+
+ else
+ if(.not. MOHO_MAP_LUPEI) then
+ vp_upper = vp_interp(8)
+ vs_upper = vs_interp(8)
+ z_upper = Z_HAUKSSON_LAYER_8
+
+ vp_lower = vp_interp(9)
+ vs_lower = vs_interp(9)
+ z_lower = Z_HAUKSSON_LAYER_9
+ !!! waiting for better interpolation of Moho maps.
+ else
+ vp_upper = vp_interp(8)
+ vs_upper = vs_interp(8)
+ z_upper = Z_HAUKSSON_LAYER_8
+
+ vp_lower = vp_interp(9)
+ vs_lower = vs_interp(9)
+ z_lower = Z_HAUKSSON_LAYER_9
+ endif
+
+ endif
+
+ gamma_interp_z = (z_eval - z_lower) / (z_upper - z_lower)
+
+ if(gamma_interp_z < -0.001d0 .or. gamma_interp_z > 1.001d0) &
+ stop 'interpolation in z is incorrect in Hauksson'
+
+ vp_final = vp_upper * gamma_interp_z + vp_lower * (1.-gamma_interp_z)
+ vs_final = vs_upper * gamma_interp_z + vs_lower * (1.-gamma_interp_z)
+
+ end subroutine hauksson_model
+
Copied: seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/interpolate_gocad_block_HR.f90 (from rev 17213, seismo/3D/SPECFEM3D/trunk/interpolate_gocad_block_HR.f90)
===================================================================
--- seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/interpolate_gocad_block_HR.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/interpolate_gocad_block_HR.f90 2010-09-24 22:17:01 UTC (rev 17215)
@@ -0,0 +1,161 @@
+
+ subroutine interpolate_gocad_block_HR(vp_block_gocad_HR,vp_block_gocad_MR, &
+ utm_x_eval,utm_y_eval,z_eval,rho_final,vp_final,vs_final,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)
+
+ implicit none
+
+ include "constants.h"
+ include "constants_gocad.h"
+
+ double precision vp_block_gocad_HR(0:NX_GOCAD_HR-1,0:NY_GOCAD_HR-1,0:NZ_GOCAD_HR-1)
+ double precision vp_block_gocad_MR(0:NX_GOCAD_MR-1,0:NY_GOCAD_MR-1,0:NZ_GOCAD_MR-1)
+
+ integer ix,iy,iz
+
+ double precision utm_x_eval,utm_y_eval,z_eval
+ double precision spacing_x,spacing_y,spacing_z
+ double precision gamma_interp_x,gamma_interp_y,gamma_interp_z
+ double precision v1,v2,v3,v4,v5,v6,v7,v8
+ double precision vp_final,vs_final,rho_final,vp_vs_ratio
+ double precision rho_ref_MR,vp_ref_MR,vs_ref_MR
+ double precision THICKNESS_TAPER_BLOCK_HR, &
+ VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM
+
+ logical point_is_in_sediments,dummy_flag,IMPOSE_MINIMUM_VP_GOCAD
+
+! for Hauksson's model
+ integer doubling_index
+ logical HAUKSSON_REGIONAL_MODEL,MOHO_MAP_LUPEI
+ double precision, dimension(NLAYERS_HAUKSSON,NGRID_NEW_HAUKSSON,NGRID_NEW_HAUKSSON) :: vp_hauksson,vs_hauksson
+
+! determine spacing and cell for linear interpolation
+ spacing_x = (utm_x_eval - ORIG_X_GOCAD_HR) / SPACING_X_GOCAD_HR
+ spacing_y = (utm_y_eval - ORIG_Y_GOCAD_HR) / SPACING_Y_GOCAD_HR
+ spacing_z = (z_eval - ORIG_Z_GOCAD_HR) / SPACING_Z_GOCAD_HR
+
+ ix = int(spacing_x)
+ iy = int(spacing_y)
+ iz = int(spacing_z)
+
+ gamma_interp_x = spacing_x - dble(ix)
+ gamma_interp_y = spacing_y - dble(iy)
+ gamma_interp_z = spacing_z - dble(iz)
+
+! this block is smaller than the grid, therefore just exit
+! if the target point is outside of the block
+ if(ix < 0 .or. ix > NX_GOCAD_HR-2 .or. iy < 0 .or. iy > NY_GOCAD_HR-2) return
+
+! suppress edge effects in vertical direction
+ if(iz < 0) then
+ iz = 0
+ gamma_interp_z = 0.d0
+ endif
+ if(iz > NZ_GOCAD_HR-2) then
+ iz = NZ_GOCAD_HR-2
+ gamma_interp_z = 1.d0
+ endif
+
+! define 8 corners of interpolation element
+ v1 = vp_block_gocad_HR(ix,iy,iz)
+ v2 = vp_block_gocad_HR(ix+1,iy,iz)
+ v3 = vp_block_gocad_HR(ix+1,iy+1,iz)
+ v4 = vp_block_gocad_HR(ix,iy+1,iz)
+
+ v5 = vp_block_gocad_HR(ix,iy,iz+1)
+ v6 = vp_block_gocad_HR(ix+1,iy,iz+1)
+ v7 = vp_block_gocad_HR(ix+1,iy+1,iz+1)
+ v8 = vp_block_gocad_HR(ix,iy+1,iz+1)
+
+! check if element is defined (i.e. is in the sediments in Voxet)
+! do nothing if element is undefined
+! a P-velocity of 20 km/s is used to indicate fictitious elements
+ if(v1 < 19000. .and. v2 < 19000. .and. &
+ v3 < 19000. .and. v4 < 19000. .and. &
+ v5 < 19000. .and. v6 < 19000. .and. &
+ v7 < 19000. .and. v8 < 19000.) then
+
+! set flag indicating whether point is in the sediments
+ point_is_in_sediments = .true.
+
+! use trilinear interpolation in cell to define Vp
+ vp_final = &
+ v1*(1.-gamma_interp_x)*(1.-gamma_interp_y)*(1.-gamma_interp_z) + &
+ v2*gamma_interp_x*(1.-gamma_interp_y)*(1.-gamma_interp_z) + &
+ v3*gamma_interp_x*gamma_interp_y*(1.-gamma_interp_z) + &
+ v4*(1.-gamma_interp_x)*gamma_interp_y*(1.-gamma_interp_z) + &
+ v5*(1.-gamma_interp_x)*(1.-gamma_interp_y)*gamma_interp_z + &
+ v6*gamma_interp_x*(1.-gamma_interp_y)*gamma_interp_z + &
+ v7*gamma_interp_x*gamma_interp_y*gamma_interp_z + &
+ v8*(1.-gamma_interp_x)*gamma_interp_y*gamma_interp_z
+
+! impose minimum velocity if needed
+ if(IMPOSE_MINIMUM_VP_GOCAD .and. vp_final < VP_MIN_GOCAD) vp_final = VP_MIN_GOCAD
+
+! taper edges to make smooth transition between MR and HR blocks
+! get value from edge of medium-resolution block
+! then use linear interpolation from edge of the model
+ if(TAPER_GOCAD_TRANSITIONS) then
+
+! x = xmin
+ if(utm_x_eval < ORIG_X_GOCAD_HR + THICKNESS_TAPER_BLOCK_HR) then
+ gamma_interp_x = (utm_x_eval - ORIG_X_GOCAD_HR) / THICKNESS_TAPER_BLOCK_HR
+ call interpolate_gocad_block_MR(vp_block_gocad_MR, &
+ ORIG_X_GOCAD_HR,utm_y_eval,z_eval,rho_ref_MR,vp_ref_MR,vs_ref_MR,dummy_flag, &
+ 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)
+ vp_final = vp_ref_MR * (1. - gamma_interp_x) + vp_final * gamma_interp_x
+
+! x = xmax
+ else if(utm_x_eval > END_X_GOCAD_HR - THICKNESS_TAPER_BLOCK_HR) then
+ gamma_interp_x = (utm_x_eval - (END_X_GOCAD_HR - THICKNESS_TAPER_BLOCK_HR)) / THICKNESS_TAPER_BLOCK_HR
+ call interpolate_gocad_block_MR(vp_block_gocad_MR, &
+ END_X_GOCAD_HR,utm_y_eval,z_eval,rho_ref_MR,vp_ref_MR,vs_ref_MR,dummy_flag, &
+ 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)
+ vp_final = vp_ref_MR * gamma_interp_x + vp_final * (1. - gamma_interp_x)
+
+! y = ymin
+ else if(utm_y_eval < ORIG_Y_GOCAD_HR + THICKNESS_TAPER_BLOCK_HR) then
+ gamma_interp_y = (utm_y_eval - ORIG_Y_GOCAD_HR) / THICKNESS_TAPER_BLOCK_HR
+ call interpolate_gocad_block_MR(vp_block_gocad_MR, &
+ utm_x_eval,ORIG_Y_GOCAD_HR,z_eval,rho_ref_MR,vp_ref_MR,vs_ref_MR,dummy_flag, &
+ 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)
+ vp_final = vp_ref_MR * (1. - gamma_interp_y) + vp_final * gamma_interp_y
+
+! y = ymax
+ else if(utm_y_eval > END_Y_GOCAD_HR - THICKNESS_TAPER_BLOCK_HR) then
+ gamma_interp_y = (utm_y_eval - (END_Y_GOCAD_HR - THICKNESS_TAPER_BLOCK_HR)) / THICKNESS_TAPER_BLOCK_HR
+ call interpolate_gocad_block_MR(vp_block_gocad_MR, &
+ utm_x_eval,END_Y_GOCAD_HR,z_eval,rho_ref_MR,vp_ref_MR,vs_ref_MR,dummy_flag, &
+ 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)
+ vp_final = vp_ref_MR * gamma_interp_y + vp_final * (1. - gamma_interp_y)
+
+ endif
+
+ endif
+
+! use linear variation of vp/vs ratio with depth, between 0. and 8.5 km
+ vp_vs_ratio = VP_VS_RATIO_GOCAD_BOTTOM + &
+ (VP_VS_RATIO_GOCAD_TOP - VP_VS_RATIO_GOCAD_BOTTOM) * &
+ (z_eval - (-8500.d0)) / (0.d0 - (-8500.d0))
+
+! make sure ratio remains in interval
+ if(vp_vs_ratio < VP_VS_RATIO_GOCAD_BOTTOM) vp_vs_ratio = VP_VS_RATIO_GOCAD_BOTTOM
+ if(vp_vs_ratio > VP_VS_RATIO_GOCAD_TOP) vp_vs_ratio = VP_VS_RATIO_GOCAD_TOP
+
+ vs_final = vp_final / vp_vs_ratio
+ call compute_rho_estimate(rho_final,vp_final)
+
+ endif
+
+ end subroutine interpolate_gocad_block_HR
+
Copied: seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/interpolate_gocad_block_MR.f90 (from rev 17213, seismo/3D/SPECFEM3D/trunk/interpolate_gocad_block_MR.f90)
===================================================================
--- seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/interpolate_gocad_block_MR.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/interpolate_gocad_block_MR.f90 2010-09-24 22:17:01 UTC (rev 17215)
@@ -0,0 +1,186 @@
+
+ subroutine interpolate_gocad_block_MR(vp_block_gocad_MR, &
+ utm_x_eval,utm_y_eval,z_eval,rho_final,vp_final,vs_final,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)
+
+ implicit none
+
+ include "constants.h"
+ include "constants_gocad.h"
+
+ double precision vp_block_gocad_MR(0:NX_GOCAD_MR-1,0:NY_GOCAD_MR-1,0:NZ_GOCAD_MR-1)
+ double precision VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM,THICKNESS_TAPER_BLOCK_MR
+
+ integer ix,iy,iz
+
+ double precision utm_x_eval,utm_y_eval,z_eval
+ double precision spacing_x,spacing_y,spacing_z
+ double precision gamma_interp_x,gamma_interp_y,gamma_interp_z
+ double precision v1,v2,v3,v4,v5,v6,v7,v8
+ double precision vp_final,vs_final,rho_final,vp_vs_ratio
+ double precision xmesh,ymesh,zmesh,vs_dummy,rho_dummy
+
+ logical point_is_in_sediments,IMPOSE_MINIMUM_VP_GOCAD
+
+! for Hauksson's model
+ integer doubling_index
+ logical HAUKSSON_REGIONAL_MODEL,MOHO_MAP_LUPEI
+ double precision vp_ref_hauksson
+ double precision, dimension(NLAYERS_HAUKSSON,NGRID_NEW_HAUKSSON,NGRID_NEW_HAUKSSON) :: vp_hauksson,vs_hauksson
+
+! determine spacing and cell for linear interpolation
+ spacing_x = (utm_x_eval - ORIG_X_GOCAD_MR) / SPACING_X_GOCAD_MR
+ spacing_y = (utm_y_eval - ORIG_Y_GOCAD_MR) / SPACING_Y_GOCAD_MR
+ spacing_z = (z_eval - ORIG_Z_GOCAD_MR) / SPACING_Z_GOCAD_MR
+
+ ix = int(spacing_x)
+ iy = int(spacing_y)
+ iz = int(spacing_z)
+
+ gamma_interp_x = spacing_x - dble(ix)
+ gamma_interp_y = spacing_y - dble(iy)
+ gamma_interp_z = spacing_z - dble(iz)
+
+! suppress edge effects for points outside of Gocad model
+ if(ix < 0) then
+ ix = 0
+ gamma_interp_x = 0.d0
+ endif
+ if(ix > NX_GOCAD_MR-2) then
+ ix = NX_GOCAD_MR-2
+ gamma_interp_x = 1.d0
+ endif
+
+ if(iy < 0) then
+ iy = 0
+ gamma_interp_y = 0.d0
+ endif
+ if(iy > NY_GOCAD_MR-2) then
+ iy = NY_GOCAD_MR-2
+ gamma_interp_y = 1.d0
+ endif
+
+ if(iz < 0) then
+ iz = 0
+ gamma_interp_z = 0.d0
+ endif
+ if(iz > NZ_GOCAD_MR-2) then
+ iz = NZ_GOCAD_MR-2
+ gamma_interp_z = 1.d0
+ endif
+
+! define 8 corners of interpolation element
+ v1 = vp_block_gocad_MR(ix,iy,iz)
+ v2 = vp_block_gocad_MR(ix+1,iy,iz)
+ v3 = vp_block_gocad_MR(ix+1,iy+1,iz)
+ v4 = vp_block_gocad_MR(ix,iy+1,iz)
+
+ v5 = vp_block_gocad_MR(ix,iy,iz+1)
+ v6 = vp_block_gocad_MR(ix+1,iy,iz+1)
+ v7 = vp_block_gocad_MR(ix+1,iy+1,iz+1)
+ v8 = vp_block_gocad_MR(ix,iy+1,iz+1)
+
+! check if element is defined (i.e. is in the sediments in Voxet)
+! do nothing if element is undefined
+! a P-velocity of 20 km/s is used to indicate fictitious elements
+ if(v1 < 19000. .and. v2 < 19000. .and. &
+ v3 < 19000. .and. v4 < 19000. .and. &
+ v5 < 19000. .and. v6 < 19000. .and. &
+ v7 < 19000. .and. v8 < 19000.) then
+
+! set flag indicating whether point is in the sediments
+ point_is_in_sediments = .true.
+
+! use trilinear interpolation in cell to define Vp
+ vp_final = &
+ v1*(1.-gamma_interp_x)*(1.-gamma_interp_y)*(1.-gamma_interp_z) + &
+ v2*gamma_interp_x*(1.-gamma_interp_y)*(1.-gamma_interp_z) + &
+ v3*gamma_interp_x*gamma_interp_y*(1.-gamma_interp_z) + &
+ v4*(1.-gamma_interp_x)*gamma_interp_y*(1.-gamma_interp_z) + &
+ v5*(1.-gamma_interp_x)*(1.-gamma_interp_y)*gamma_interp_z + &
+ v6*gamma_interp_x*(1.-gamma_interp_y)*gamma_interp_z + &
+ v7*gamma_interp_x*gamma_interp_y*gamma_interp_z + &
+ v8*(1.-gamma_interp_x)*gamma_interp_y*gamma_interp_z
+
+! impose minimum velocity if needed
+ if(IMPOSE_MINIMUM_VP_GOCAD .and. vp_final < VP_MIN_GOCAD) vp_final = VP_MIN_GOCAD
+
+! taper edges to make smooth transition between Hauksson and MR blocks
+! get value from edge of medium-resolution block
+! then use linear interpolation from edge of the model
+ if(TAPER_GOCAD_TRANSITIONS) then
+
+! x = xmin
+ if(utm_x_eval < ORIG_X_GOCAD_MR + THICKNESS_TAPER_BLOCK_MR) then
+ xmesh = ORIG_X_GOCAD_MR
+ ymesh = utm_y_eval
+ zmesh = z_eval
+ if(HAUKSSON_REGIONAL_MODEL) then
+ call hauksson_model(vp_hauksson,vs_hauksson,xmesh,ymesh,zmesh,vp_ref_hauksson,vs_dummy, MOHO_MAP_LUPEI)
+ else
+ call socal_model(doubling_index,rho_dummy,vp_ref_hauksson,vs_dummy)
+ endif
+ gamma_interp_x = (utm_x_eval - ORIG_X_GOCAD_MR) / THICKNESS_TAPER_BLOCK_MR
+ vp_final = vp_ref_hauksson * (1. - gamma_interp_x) + vp_final * gamma_interp_x
+
+! x = xmax
+ else if(utm_x_eval > END_X_GOCAD_MR - THICKNESS_TAPER_BLOCK_MR) then
+ xmesh = END_X_GOCAD_MR
+ ymesh = utm_y_eval
+ zmesh = z_eval
+ if(HAUKSSON_REGIONAL_MODEL) then
+ call hauksson_model(vp_hauksson,vs_hauksson,xmesh,ymesh,zmesh,vp_ref_hauksson,vs_dummy, MOHO_MAP_LUPEI)
+ else
+ call socal_model(doubling_index,rho_dummy,vp_ref_hauksson,vs_dummy)
+ endif
+ gamma_interp_x = (utm_x_eval - (END_X_GOCAD_MR - THICKNESS_TAPER_BLOCK_MR)) / THICKNESS_TAPER_BLOCK_MR
+ vp_final = vp_ref_hauksson * gamma_interp_x + vp_final * (1. - gamma_interp_x)
+
+! y = ymin
+ else if(utm_y_eval < ORIG_Y_GOCAD_MR + THICKNESS_TAPER_BLOCK_MR) then
+ xmesh = utm_x_eval
+ ymesh = ORIG_Y_GOCAD_MR
+ zmesh = z_eval
+ if(HAUKSSON_REGIONAL_MODEL) then
+ call hauksson_model(vp_hauksson,vs_hauksson,xmesh,ymesh,zmesh,vp_ref_hauksson,vs_dummy, MOHO_MAP_LUPEI)
+ else
+ call socal_model(doubling_index,rho_dummy,vp_ref_hauksson,vs_dummy)
+ endif
+ gamma_interp_y = (utm_y_eval - ORIG_Y_GOCAD_MR) / THICKNESS_TAPER_BLOCK_MR
+ vp_final = vp_ref_hauksson * (1. - gamma_interp_y) + vp_final * gamma_interp_y
+
+! y = ymax
+ else if(utm_y_eval > END_Y_GOCAD_MR - THICKNESS_TAPER_BLOCK_MR) then
+ xmesh = utm_x_eval
+ ymesh = END_Y_GOCAD_MR
+ zmesh = z_eval
+ if(HAUKSSON_REGIONAL_MODEL) then
+ call hauksson_model(vp_hauksson,vs_hauksson,xmesh,ymesh,zmesh,vp_ref_hauksson,vs_dummy, MOHO_MAP_LUPEI)
+ else
+ call socal_model(doubling_index,rho_dummy,vp_ref_hauksson,vs_dummy)
+ endif
+ gamma_interp_y = (utm_y_eval - (END_Y_GOCAD_MR - THICKNESS_TAPER_BLOCK_MR)) / THICKNESS_TAPER_BLOCK_MR
+ vp_final = vp_ref_hauksson * gamma_interp_y + vp_final * (1. - gamma_interp_y)
+
+ endif
+
+ endif
+
+! use linear variation of vp/vs ratio with depth, between 0. and 8.5 km
+ vp_vs_ratio = VP_VS_RATIO_GOCAD_BOTTOM + &
+ (VP_VS_RATIO_GOCAD_TOP - VP_VS_RATIO_GOCAD_BOTTOM) * &
+ (z_eval - (-8500.d0)) / (0.d0 - (-8500.d0))
+
+! make sure ratio remains in interval
+ if(vp_vs_ratio < VP_VS_RATIO_GOCAD_BOTTOM) vp_vs_ratio = VP_VS_RATIO_GOCAD_BOTTOM
+ if(vp_vs_ratio > VP_VS_RATIO_GOCAD_TOP) vp_vs_ratio = VP_VS_RATIO_GOCAD_TOP
+
+ vs_final = vp_final / vp_vs_ratio
+ call compute_rho_estimate(rho_final,vp_final)
+
+ endif
+
+ end subroutine interpolate_gocad_block_MR
+
Copied: seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/mesh_vertical.f90 (from rev 17213, seismo/3D/SPECFEM3D/trunk/mesh_vertical.f90)
===================================================================
--- seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/mesh_vertical.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/mesh_vertical.f90 2010-09-24 22:17:01 UTC (rev 17215)
@@ -0,0 +1,120 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology September 2006
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+ subroutine mesh_vertical(myrank,rn,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)
+
+! create the vertical mesh, honoring the major discontinuities in the model
+
+ implicit none
+
+ include "constants.h"
+
+ integer myrank
+ integer NER,NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM,NER_SEDIM
+ logical MOHO_MAP_LUPEI
+ double precision Z_DEPTH_BLOCK,Z_BASEMENT_SURFACE,Z_DEPTH_MOHO
+ double precision rn(0:2*NER)
+
+!! DK DK UGLY modif z_top by Emmanuel Chaljub here
+!! DK DK UGLY modif Manu removed double precision z_top
+
+ integer npr,ir
+
+ npr = -1
+
+!
+!--- bottom of the mesh (Z_DEPTH_BLOCK) to Moho
+!
+ do ir=0,2*NER_BOTTOM_MOHO-1
+ npr=npr+1
+ rn(npr)=(Z_DEPTH_MOHO-Z_DEPTH_BLOCK)*dble(ir)/dble(2*NER_BOTTOM_MOHO)
+ enddo
+
+! do not use d16km when Moho map is honored
+ if(MOHO_MAP_LUPEI) then
+
+!
+!--- Moho to modified basement surface
+!
+ do ir=0,2*(NER_MOHO_16+NER_16_BASEMENT)-1
+ npr=npr+1
+ rn(npr)=(Z_DEPTH_MOHO-Z_DEPTH_BLOCK) + (Z_BASEMENT_SURFACE-Z_DEPTH_MOHO)*dble(ir)/dble(2*(NER_MOHO_16+NER_16_BASEMENT))
+ enddo
+
+ else
+!
+!--- Moho to d16km
+!
+ do ir=0,2*NER_MOHO_16-1
+ npr=npr+1
+ rn(npr)=(Z_DEPTH_MOHO-Z_DEPTH_BLOCK) + (DEPTH_16km_SOCAL-Z_DEPTH_MOHO)*dble(ir)/dble(2*NER_MOHO_16)
+ enddo
+!
+!--- d16km to modified basement surface
+!
+ do ir=0,2*NER_16_BASEMENT-1
+ npr=npr+1
+ rn(npr)=(DEPTH_16km_SOCAL-Z_DEPTH_BLOCK) + (Z_BASEMENT_SURFACE-DEPTH_16km_SOCAL)*dble(ir)/dble(2*NER_16_BASEMENT)
+ enddo
+
+ endif
+
+!
+!--- modified basement surface to surface of model (topography/bathymetry)
+!
+! also create last point exactly at the surface
+! other regions above stop one point below
+ do ir=0,2*(NER_BASEMENT_SEDIM+NER_SEDIM) - 0
+ npr=npr+1
+ rn(npr)=(Z_BASEMENT_SURFACE-Z_DEPTH_BLOCK) + &
+!! DK DK UGLY modif z_top by Emmanuel Chaljub here
+!! DK DK UGLY suppressed Manu's modif and put old code back because better mesh
+!! DK DK UGLY investigate this in detail one day
+ (Z_SURFACE-Z_BASEMENT_SURFACE)*dble(ir)/dble(2*(NER_BASEMENT_SEDIM+NER_SEDIM))
+!! DK DK UGLY modif Manu removed (z_top-Z_BASEMENT_SURFACE)*dble(ir)/dble(2*(NER_BASEMENT_SEDIM+NER_SEDIM))
+ enddo
+
+! normalize depths
+!! DK DK UGLY modif z_top by Emmanuel Chaljub here
+!! DK DK UGLY suppressed Manu's modif and put old code back because better mesh
+!! DK DK UGLY investigate this in detail one day
+!! DK DK UGLY modif Manu removed rn(:) = rn(:) / (z_top-Z_DEPTH_BLOCK)
+!! DK DK UGLY modif Manu removed
+ rn(:) = rn(:) / (Z_SURFACE-Z_DEPTH_BLOCK)
+
+! check that the mesh that has been generated is correct
+ if(npr /= 2*NER) call exit_MPI(myrank,'incorrect intervals for model')
+
+! check that vertical spacing makes sense
+ do ir=0,2*NER-1
+ if(rn(ir+1) < rn(ir)) call exit_MPI(myrank,'incorrect vertical spacing for model')
+ enddo
+
+ end subroutine mesh_vertical
+
Copied: seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/read_arrays_buffers_solver.f90 (from rev 17213, seismo/3D/SPECFEM3D/trunk/read_arrays_buffers_solver.f90)
===================================================================
--- seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/read_arrays_buffers_solver.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/read_arrays_buffers_solver.f90 2010-09-24 22:17:01 UTC (rev 17215)
@@ -0,0 +1,150 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology September 2006
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+ subroutine read_arrays_buffers_solver(myrank, &
+ iboolleft_xi,iboolright_xi,iboolleft_eta,iboolright_eta, &
+ npoin2D_xi,npoin2D_eta, &
+ NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,LOCAL_PATH)
+
+ implicit none
+
+ include "constants.h"
+
+ integer myrank
+
+ integer npoin2D_xi,npoin2D_eta
+ integer NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX
+
+ character(len=256) LOCAL_PATH
+
+ integer, dimension(NPOIN2DMAX_XMIN_XMAX) :: iboolleft_xi,iboolright_xi
+ integer, dimension(NPOIN2DMAX_YMIN_YMAX) :: iboolleft_eta,iboolright_eta
+
+ integer npoin2D_xi_mesher,npoin2D_eta_mesher
+
+ double precision xdummy,ydummy,zdummy
+
+! processor identification
+ character(len=256) prname
+
+! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+
+! create the name for the database of the current slide and region
+ call create_name_database(prname,myrank,LOCAL_PATH)
+
+! read 2-D addressing for summation between slices along xi with MPI
+
+! read iboolleft_xi of this slice
+ open(unit=IIN,file=prname(1:len_trim(prname))//'iboolleft_xi.txt',status='old',action='read')
+ npoin2D_xi = 1
+ 350 continue
+ read(IIN,*) iboolleft_xi(npoin2D_xi),xdummy,ydummy,zdummy
+ if(iboolleft_xi(npoin2D_xi) > 0) then
+ npoin2D_xi = npoin2D_xi + 1
+ goto 350
+ endif
+! subtract the line that contains the flag after the last point
+ npoin2D_xi = npoin2D_xi - 1
+! read nb of points given by the mesher
+ read(IIN,*) npoin2D_xi_mesher
+ if(npoin2D_xi > NPOIN2DMAX_XMIN_XMAX .or. npoin2D_xi /= npoin2D_xi_mesher) &
+ call exit_MPI(myrank,'incorrect iboolleft_xi read')
+ close(IIN)
+
+! read iboolright_xi of this slice
+ open(unit=IIN,file=prname(1:len_trim(prname))//'iboolright_xi.txt',status='old',action='read')
+ npoin2D_xi = 1
+ 360 continue
+ read(IIN,*) iboolright_xi(npoin2D_xi),xdummy,ydummy,zdummy
+ if(iboolright_xi(npoin2D_xi) > 0) then
+ npoin2D_xi = npoin2D_xi + 1
+ goto 360
+ endif
+! subtract the line that contains the flag after the last point
+ npoin2D_xi = npoin2D_xi - 1
+! read nb of points given by the mesher
+ read(IIN,*) npoin2D_xi_mesher
+ if(npoin2D_xi > NPOIN2DMAX_XMIN_XMAX .or. npoin2D_xi /= npoin2D_xi_mesher) &
+ call exit_MPI(myrank,'incorrect iboolright_xi read')
+ close(IIN)
+
+ if(myrank == 0) then
+ write(IMAIN,*)
+ write(IMAIN,*) '# of points in MPI buffers along xi npoin2D_xi = ', &
+ npoin2D_xi
+ write(IMAIN,*) '# of array elements transferred npoin2D_xi*NDIM = ', &
+ npoin2D_xi*NDIM
+ write(IMAIN,*)
+ endif
+
+! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+
+! read 2-D addressing for summation between slices along eta with MPI
+
+! read iboolleft_eta of this slice
+ open(unit=IIN,file=prname(1:len_trim(prname))//'iboolleft_eta.txt',status='old',action='read')
+ npoin2D_eta = 1
+ 370 continue
+ read(IIN,*) iboolleft_eta(npoin2D_eta),xdummy,ydummy,zdummy
+ if(iboolleft_eta(npoin2D_eta) > 0) then
+ npoin2D_eta = npoin2D_eta + 1
+ goto 370
+ endif
+! subtract the line that contains the flag after the last point
+ npoin2D_eta = npoin2D_eta - 1
+! read nb of points given by the mesher
+ read(IIN,*) npoin2D_eta_mesher
+ if(npoin2D_eta > NPOIN2DMAX_YMIN_YMAX .or. npoin2D_eta /= npoin2D_eta_mesher) &
+ call exit_MPI(myrank,'incorrect iboolleft_eta read')
+ close(IIN)
+
+! read iboolright_eta of this slice
+ open(unit=IIN,file=prname(1:len_trim(prname))//'iboolright_eta.txt',status='old',action='read')
+ npoin2D_eta = 1
+ 380 continue
+ read(IIN,*) iboolright_eta(npoin2D_eta),xdummy,ydummy,zdummy
+ if(iboolright_eta(npoin2D_eta) > 0) then
+ npoin2D_eta = npoin2D_eta + 1
+ goto 380
+ endif
+! subtract the line that contains the flag after the last point
+ npoin2D_eta = npoin2D_eta - 1
+! read nb of points given by the mesher
+ read(IIN,*) npoin2D_eta_mesher
+ if(npoin2D_eta > NPOIN2DMAX_YMIN_YMAX .or. npoin2D_eta /= npoin2D_eta_mesher) &
+ call exit_MPI(myrank,'incorrect iboolright_eta read')
+ close(IIN)
+
+ if(myrank == 0) then
+ write(IMAIN,*)
+ write(IMAIN,*) '# of points in MPI buffers along eta npoin2D_eta = ', &
+ npoin2D_eta
+ write(IMAIN,*) '# of array elements transferred npoin2D_eta*NDIM = ', &
+ npoin2D_eta*NDIM
+ write(IMAIN,*)
+ endif
+
+ end subroutine read_arrays_buffers_solver
+
Copied: seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/read_moho_map.f90 (from rev 17213, seismo/3D/SPECFEM3D/trunk/read_moho_map.f90)
===================================================================
--- seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/read_moho_map.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/read_moho_map.f90 2010-09-24 22:17:01 UTC (rev 17215)
@@ -0,0 +1,60 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology September 2006
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+ subroutine read_moho_map(imoho_depth)
+!
+!---- read Lupei Zhu's Moho map of Southern California once and for all
+!
+ implicit none
+
+ include "constants.h"
+
+! use integer array to store Moho depth
+ integer imoho_depth(NX_MOHO,NY_MOHO)
+
+ integer ix,iy
+
+ double precision long,lat,depth_km
+
+ character(len=256) MOHO_MAP_FILE
+
+ imoho_depth(:,:) = 0
+
+ call get_value_string(MOHO_MAP_FILE, &
+ 'model.MOHO_MAP_FILE', &
+ 'DATA/moho_map/moho_lupei_zhu.dat')
+ open(unit=13,file=MOHO_MAP_FILE,status='old',action='read')
+! file starts from North-West corner
+ do iy=NY_MOHO,1,-1
+ do ix=1,NX_MOHO
+ read(13,*) long,lat,depth_km
+! convert depth to meters
+ imoho_depth(ix,iy) = nint(depth_km * 1000.d0)
+ enddo
+ enddo
+ close(13)
+
+ end subroutine read_moho_map
+
Copied: seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/salton_trough_gocad.f90 (from rev 17213, seismo/3D/SPECFEM3D/trunk/salton_trough_gocad.f90)
===================================================================
--- seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/salton_trough_gocad.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/salton_trough_gocad.f90 2010-09-24 22:17:01 UTC (rev 17215)
@@ -0,0 +1,168 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology September 2006
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=======================================================================
+
+subroutine read_salton_sea_model(vp_array)
+
+ implicit none
+
+ include 'constants.h'
+ include 'constants_gocad.h'
+
+ real :: vp_array(GOCAD_ST_NU,GOCAD_ST_NV,GOCAD_ST_NW)
+ integer :: ios, reclen
+
+ character(len=256) SALTON_SEA_MODEL_FILE
+
+ reclen=(GOCAD_ST_NU * GOCAD_ST_NV * GOCAD_ST_NW) * 4
+ call get_value_string(SALTON_SEA_MODEL_FILE, &
+ 'model.SALTON_SEA_MODEL_FILE', &
+ 'DATA/st_3D_block_harvard/regrid3_vel_p.bin')
+ open(11,file=SALTON_SEA_MODEL_FILE,status='old',action='read',form='unformatted',access='direct',recl=reclen,iostat=ios)
+ if (ios /= 0) then
+ print *, 'iostat = ', ios
+ stop 'Error opening file'
+ endif
+ read(11,rec=1,iostat=ios) vp_array
+ if (ios /= 0) stop 'Error reading vp_array'
+ close(11)
+
+end subroutine read_salton_sea_model
+
+
+subroutine vx_xyz2uvw(xmesh, ymesh, zmesh, uc, vc, wc)
+
+
+ implicit none
+ include 'constants.h'
+
+ double precision :: xmesh, ymesh, zmesh, uc, vc, wc
+
+ uc = (GOCAD_ST_NU-1) * ((xmesh - GOCAD_ST_O_X) * GOCAD_ST_V_Y - (ymesh - GOCAD_ST_O_Y) * GOCAD_ST_V_X) &
+ / (GOCAD_ST_U_X * GOCAD_ST_V_Y - GOCAD_ST_U_Y * GOCAD_ST_V_X)
+ vc = (GOCAD_ST_NV-1) * ((ymesh - GOCAD_ST_O_Y) - uc * GOCAD_ST_U_Y/(GOCAD_ST_NU-1) ) / GOCAD_ST_V_Y
+ wc = (GOCAD_ST_NW-1) * (zmesh - GOCAD_ST_O_Z) / GOCAD_ST_W_Z
+
+end subroutine vx_xyz2uvw
+
+
+subroutine vx_xyz_interp(uc,vc,wc, vp, vs, rho, vp_array)
+
+ implicit none
+ include 'constants.h'
+
+ double precision :: uc,vc,wc, vp, vs, rho
+ real :: vp_array(GOCAD_ST_NU,GOCAD_ST_NV,GOCAD_ST_NW)
+
+ integer :: i,j,k,ixi,ieta,iga
+ real :: v1, v2, v3, v4, v5, v6, v7, v8, xi, eta, ga, vi
+ double precision :: zmesh
+ real,parameter :: eps = 1.0e-3
+
+
+ i = uc + 1
+ j = vc + 1
+ k = wc + 1
+
+ xi = uc + 1 - i
+ eta = vc + 1- j
+ ga = wc + 1 -k
+
+ ixi = nint(xi)
+ ieta = nint(eta)
+ iga = nint(ga)
+
+! print *, 'gc = ', i, j, k
+! print *, 'xi, eta, ga = ', xi, eta, ga
+! print *, 'ixi, ieta, iga = ', ixi, ieta, iga
+
+
+ if (i > 0 .or. i < GOCAD_ST_NU .or. j > 0 .or. j < GOCAD_ST_NV .or. k > 0 .or. k < GOCAD_ST_NW) then
+ v1 = vp_array(i,j,k)
+ v2 = vp_array(i+1,j,k)
+ v3 = vp_array(i+1,j+1,k)
+ v4 = vp_array(i,j+1,k)
+ v5 = vp_array(i,j,k+1)
+ v6 = vp_array(i+1,j,k+1)
+ v7 = vp_array(i+1,j+1,k+1)
+ v8 = vp_array(i,j+1,k+1)
+ vi = vp_array(i+ixi,j+ieta,k+iga)
+! print *, v1, v2, v3, v4, v5, v6, v7, v8
+ if ((v1 - GOCAD_ST_NO_DATA_VALUE) > eps .and. &
+ (v2 - GOCAD_ST_NO_DATA_VALUE) > eps .and. &
+ (v3 - GOCAD_ST_NO_DATA_VALUE) > eps .and. &
+ (v4 - GOCAD_ST_NO_DATA_VALUE) > eps .and. &
+ (v5 - GOCAD_ST_NO_DATA_VALUE) > eps .and. &
+ (v6 - GOCAD_ST_NO_DATA_VALUE) > eps .and. &
+ (v7 - GOCAD_ST_NO_DATA_VALUE) > eps .and. &
+ (v8 - GOCAD_ST_NO_DATA_VALUE) > eps ) then
+ vp = dble(&
+ v1 * (1-xi) * (1-eta) * (1-ga) +&
+ v2 * xi * (1-eta) * (1-ga) +&
+ v3 * xi * eta * (1-ga) +&
+ v4 * (1-xi) * eta * (1-ga) +&
+ v5 * (1-xi) * (1-eta) * ga +&
+ v6 * xi * (1-eta) * ga +&
+ v7 * xi * eta * ga +&
+ v8 * (1-xi) * eta * ga)
+ else if ((vi - GOCAD_ST_NO_DATA_VALUE) > eps) then
+ vp = dble(vi)
+! else if ((v1 - GOCAD_ST_NO_DATA_VALUE) > eps) then
+! vp = dble(v1)
+! else if ((v2 - GOCAD_ST_NO_DATA_VALUE) > eps) then
+! vp = dble(v2)
+! else if ((v3 - GOCAD_ST_NO_DATA_VALUE) > eps) then
+! vp = dble(v3)
+! else if ((v4 - GOCAD_ST_NO_DATA_VALUE) > eps) then
+! vp = dble(v4)
+! else if ((v5 - GOCAD_ST_NO_DATA_VALUE) > eps) then
+! vp = dble(v5)
+! else if ((v6 - GOCAD_ST_NO_DATA_VALUE) > eps) then
+! vp = dble(v6)
+! else if ((v7 - GOCAD_ST_NO_DATA_VALUE) > eps) then
+! vp = dble(v7)
+! else if ((v7 - GOCAD_ST_NO_DATA_VALUE) > eps) then
+! vp = dble(v8)
+ else
+ vp = GOCAD_ST_NO_DATA_VALUE
+ endif
+ zmesh = wc / (GOCAD_ST_NW - 1) * GOCAD_ST_W_Z + GOCAD_ST_O_Z
+ if (zmesh > -8500.) then
+ vs = vp / (2 - (0.27*zmesh/(-8500)))
+ else
+ vs = vp/1.73
+ endif
+ if (vp > 2160.) then
+ rho = vp/3 + 1280.
+ else
+ rho = 2000.
+ endif
+ else
+ rho = GOCAD_ST_NO_DATA_VALUE
+ vp = GOCAD_ST_NO_DATA_VALUE
+ vs = GOCAD_ST_NO_DATA_VALUE
+ endif
+
+end subroutine vx_xyz_interp
+
Copied: seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/socal_model.f90 (from rev 17213, seismo/3D/SPECFEM3D/trunk/socal_model.f90)
===================================================================
--- seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/socal_model.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/obsolete_routines_from_SPECFEM3D_v1.4.4/socal_model.f90 2010-09-24 22:17:01 UTC (rev 17215)
@@ -0,0 +1,63 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology September 2006
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+ subroutine socal_model(idoubling,rho,vp,vs)
+
+ implicit none
+
+ include "constants.h"
+ include "constants_gocad.h"
+
+ integer idoubling
+ double precision rho,vp,vs
+
+ if(idoubling == IFLAG_HALFSPACE_MOHO) then
+ vp=7.8d0
+ vs=4.5d0
+ rho=3.0d0
+
+ else if(idoubling == IFLAG_MOHO_16km) then
+ vp=6.7d0
+ vs=3.87d0
+ rho=2.8d0
+
+ else if(idoubling == IFLAG_ONE_LAYER_TOPOGRAPHY .or. idoubling == IFLAG_BASEMENT_TOPO) then
+ vp=5.5d0
+ vs=3.18d0
+ rho=2.4d0
+
+ else
+ vp=6.3d0
+ vs=3.64d0
+ rho=2.67d0
+ endif
+
+! scale to standard units
+ vp = vp * 1000.d0
+ vs = vs * 1000.d0
+ rho = rho * 1000.d0
+
+ end subroutine socal_model
+
Deleted: seismo/3D/SPECFEM3D/trunk/read_arrays_buffers_solver.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/read_arrays_buffers_solver.f90 2010-09-24 22:15:18 UTC (rev 17214)
+++ seismo/3D/SPECFEM3D/trunk/read_arrays_buffers_solver.f90 2010-09-24 22:17:01 UTC (rev 17215)
@@ -1,150 +0,0 @@
-!=====================================================================
-!
-! S p e c f e m 3 D V e r s i o n 1 . 4
-! ---------------------------------------
-!
-! Dimitri Komatitsch and Jeroen Tromp
-! Seismological Laboratory - California Institute of Technology
-! (c) California Institute of Technology September 2006
-!
-! This program is free software; you can redistribute it and/or modify
-! it under the terms of the GNU General Public License as published by
-! the Free Software Foundation; either version 2 of the License, or
-! (at your option) any later version.
-!
-! This program is distributed in the hope that it will be useful,
-! but WITHOUT ANY WARRANTY; without even the implied warranty of
-! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-! GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along
-! with this program; if not, write to the Free Software Foundation, Inc.,
-! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-!
-!=====================================================================
-
- subroutine read_arrays_buffers_solver(myrank, &
- iboolleft_xi,iboolright_xi,iboolleft_eta,iboolright_eta, &
- npoin2D_xi,npoin2D_eta, &
- NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,LOCAL_PATH)
-
- implicit none
-
- include "constants.h"
-
- integer myrank
-
- integer npoin2D_xi,npoin2D_eta
- integer NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX
-
- character(len=256) LOCAL_PATH
-
- integer, dimension(NPOIN2DMAX_XMIN_XMAX) :: iboolleft_xi,iboolright_xi
- integer, dimension(NPOIN2DMAX_YMIN_YMAX) :: iboolleft_eta,iboolright_eta
-
- integer npoin2D_xi_mesher,npoin2D_eta_mesher
-
- double precision xdummy,ydummy,zdummy
-
-! processor identification
- character(len=256) prname
-
-! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
-
-! create the name for the database of the current slide and region
- call create_name_database(prname,myrank,LOCAL_PATH)
-
-! read 2-D addressing for summation between slices along xi with MPI
-
-! read iboolleft_xi of this slice
- open(unit=IIN,file=prname(1:len_trim(prname))//'iboolleft_xi.txt',status='old',action='read')
- npoin2D_xi = 1
- 350 continue
- read(IIN,*) iboolleft_xi(npoin2D_xi),xdummy,ydummy,zdummy
- if(iboolleft_xi(npoin2D_xi) > 0) then
- npoin2D_xi = npoin2D_xi + 1
- goto 350
- endif
-! subtract the line that contains the flag after the last point
- npoin2D_xi = npoin2D_xi - 1
-! read nb of points given by the mesher
- read(IIN,*) npoin2D_xi_mesher
- if(npoin2D_xi > NPOIN2DMAX_XMIN_XMAX .or. npoin2D_xi /= npoin2D_xi_mesher) &
- call exit_MPI(myrank,'incorrect iboolleft_xi read')
- close(IIN)
-
-! read iboolright_xi of this slice
- open(unit=IIN,file=prname(1:len_trim(prname))//'iboolright_xi.txt',status='old',action='read')
- npoin2D_xi = 1
- 360 continue
- read(IIN,*) iboolright_xi(npoin2D_xi),xdummy,ydummy,zdummy
- if(iboolright_xi(npoin2D_xi) > 0) then
- npoin2D_xi = npoin2D_xi + 1
- goto 360
- endif
-! subtract the line that contains the flag after the last point
- npoin2D_xi = npoin2D_xi - 1
-! read nb of points given by the mesher
- read(IIN,*) npoin2D_xi_mesher
- if(npoin2D_xi > NPOIN2DMAX_XMIN_XMAX .or. npoin2D_xi /= npoin2D_xi_mesher) &
- call exit_MPI(myrank,'incorrect iboolright_xi read')
- close(IIN)
-
- if(myrank == 0) then
- write(IMAIN,*)
- write(IMAIN,*) '# of points in MPI buffers along xi npoin2D_xi = ', &
- npoin2D_xi
- write(IMAIN,*) '# of array elements transferred npoin2D_xi*NDIM = ', &
- npoin2D_xi*NDIM
- write(IMAIN,*)
- endif
-
-! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
-
-! read 2-D addressing for summation between slices along eta with MPI
-
-! read iboolleft_eta of this slice
- open(unit=IIN,file=prname(1:len_trim(prname))//'iboolleft_eta.txt',status='old',action='read')
- npoin2D_eta = 1
- 370 continue
- read(IIN,*) iboolleft_eta(npoin2D_eta),xdummy,ydummy,zdummy
- if(iboolleft_eta(npoin2D_eta) > 0) then
- npoin2D_eta = npoin2D_eta + 1
- goto 370
- endif
-! subtract the line that contains the flag after the last point
- npoin2D_eta = npoin2D_eta - 1
-! read nb of points given by the mesher
- read(IIN,*) npoin2D_eta_mesher
- if(npoin2D_eta > NPOIN2DMAX_YMIN_YMAX .or. npoin2D_eta /= npoin2D_eta_mesher) &
- call exit_MPI(myrank,'incorrect iboolleft_eta read')
- close(IIN)
-
-! read iboolright_eta of this slice
- open(unit=IIN,file=prname(1:len_trim(prname))//'iboolright_eta.txt',status='old',action='read')
- npoin2D_eta = 1
- 380 continue
- read(IIN,*) iboolright_eta(npoin2D_eta),xdummy,ydummy,zdummy
- if(iboolright_eta(npoin2D_eta) > 0) then
- npoin2D_eta = npoin2D_eta + 1
- goto 380
- endif
-! subtract the line that contains the flag after the last point
- npoin2D_eta = npoin2D_eta - 1
-! read nb of points given by the mesher
- read(IIN,*) npoin2D_eta_mesher
- if(npoin2D_eta > NPOIN2DMAX_YMIN_YMAX .or. npoin2D_eta /= npoin2D_eta_mesher) &
- call exit_MPI(myrank,'incorrect iboolright_eta read')
- close(IIN)
-
- if(myrank == 0) then
- write(IMAIN,*)
- write(IMAIN,*) '# of points in MPI buffers along eta npoin2D_eta = ', &
- npoin2D_eta
- write(IMAIN,*) '# of array elements transferred npoin2D_eta*NDIM = ', &
- npoin2D_eta*NDIM
- write(IMAIN,*)
- endif
-
- end subroutine read_arrays_buffers_solver
-
Deleted: seismo/3D/SPECFEM3D/trunk/read_moho_map.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/read_moho_map.f90 2010-09-24 22:15:18 UTC (rev 17214)
+++ seismo/3D/SPECFEM3D/trunk/read_moho_map.f90 2010-09-24 22:17:01 UTC (rev 17215)
@@ -1,60 +0,0 @@
-!=====================================================================
-!
-! S p e c f e m 3 D V e r s i o n 1 . 4
-! ---------------------------------------
-!
-! Dimitri Komatitsch and Jeroen Tromp
-! Seismological Laboratory - California Institute of Technology
-! (c) California Institute of Technology September 2006
-!
-! This program is free software; you can redistribute it and/or modify
-! it under the terms of the GNU General Public License as published by
-! the Free Software Foundation; either version 2 of the License, or
-! (at your option) any later version.
-!
-! This program is distributed in the hope that it will be useful,
-! but WITHOUT ANY WARRANTY; without even the implied warranty of
-! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-! GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along
-! with this program; if not, write to the Free Software Foundation, Inc.,
-! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-!
-!=====================================================================
-
- subroutine read_moho_map(imoho_depth)
-!
-!---- read Lupei Zhu's Moho map of Southern California once and for all
-!
- implicit none
-
- include "constants.h"
-
-! use integer array to store Moho depth
- integer imoho_depth(NX_MOHO,NY_MOHO)
-
- integer ix,iy
-
- double precision long,lat,depth_km
-
- character(len=256) MOHO_MAP_FILE
-
- imoho_depth(:,:) = 0
-
- call get_value_string(MOHO_MAP_FILE, &
- 'model.MOHO_MAP_FILE', &
- 'DATA/moho_map/moho_lupei_zhu.dat')
- open(unit=13,file=MOHO_MAP_FILE,status='old',action='read')
-! file starts from North-West corner
- do iy=NY_MOHO,1,-1
- do ix=1,NX_MOHO
- read(13,*) long,lat,depth_km
-! convert depth to meters
- imoho_depth(ix,iy) = nint(depth_km * 1000.d0)
- enddo
- enddo
- close(13)
-
- end subroutine read_moho_map
-
Deleted: seismo/3D/SPECFEM3D/trunk/salton_trough_gocad.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/salton_trough_gocad.f90 2010-09-24 22:15:18 UTC (rev 17214)
+++ seismo/3D/SPECFEM3D/trunk/salton_trough_gocad.f90 2010-09-24 22:17:01 UTC (rev 17215)
@@ -1,168 +0,0 @@
-!=====================================================================
-!
-! S p e c f e m 3 D V e r s i o n 1 . 4
-! ---------------------------------------
-!
-! Dimitri Komatitsch and Jeroen Tromp
-! Seismological Laboratory - California Institute of Technology
-! (c) California Institute of Technology September 2006
-!
-! This program is free software; you can redistribute it and/or modify
-! it under the terms of the GNU General Public License as published by
-! the Free Software Foundation; either version 2 of the License, or
-! (at your option) any later version.
-!
-! This program is distributed in the hope that it will be useful,
-! but WITHOUT ANY WARRANTY; without even the implied warranty of
-! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-! GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along
-! with this program; if not, write to the Free Software Foundation, Inc.,
-! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-!
-!=======================================================================
-
-subroutine read_salton_sea_model(vp_array)
-
- implicit none
-
- include 'constants.h'
- include 'constants_gocad.h'
-
- real :: vp_array(GOCAD_ST_NU,GOCAD_ST_NV,GOCAD_ST_NW)
- integer :: ios, reclen
-
- character(len=256) SALTON_SEA_MODEL_FILE
-
- reclen=(GOCAD_ST_NU * GOCAD_ST_NV * GOCAD_ST_NW) * 4
- call get_value_string(SALTON_SEA_MODEL_FILE, &
- 'model.SALTON_SEA_MODEL_FILE', &
- 'DATA/st_3D_block_harvard/regrid3_vel_p.bin')
- open(11,file=SALTON_SEA_MODEL_FILE,status='old',action='read',form='unformatted',access='direct',recl=reclen,iostat=ios)
- if (ios /= 0) then
- print *, 'iostat = ', ios
- stop 'Error opening file'
- endif
- read(11,rec=1,iostat=ios) vp_array
- if (ios /= 0) stop 'Error reading vp_array'
- close(11)
-
-end subroutine read_salton_sea_model
-
-
-subroutine vx_xyz2uvw(xmesh, ymesh, zmesh, uc, vc, wc)
-
-
- implicit none
- include 'constants.h'
-
- double precision :: xmesh, ymesh, zmesh, uc, vc, wc
-
- uc = (GOCAD_ST_NU-1) * ((xmesh - GOCAD_ST_O_X) * GOCAD_ST_V_Y - (ymesh - GOCAD_ST_O_Y) * GOCAD_ST_V_X) &
- / (GOCAD_ST_U_X * GOCAD_ST_V_Y - GOCAD_ST_U_Y * GOCAD_ST_V_X)
- vc = (GOCAD_ST_NV-1) * ((ymesh - GOCAD_ST_O_Y) - uc * GOCAD_ST_U_Y/(GOCAD_ST_NU-1) ) / GOCAD_ST_V_Y
- wc = (GOCAD_ST_NW-1) * (zmesh - GOCAD_ST_O_Z) / GOCAD_ST_W_Z
-
-end subroutine vx_xyz2uvw
-
-
-subroutine vx_xyz_interp(uc,vc,wc, vp, vs, rho, vp_array)
-
- implicit none
- include 'constants.h'
-
- double precision :: uc,vc,wc, vp, vs, rho
- real :: vp_array(GOCAD_ST_NU,GOCAD_ST_NV,GOCAD_ST_NW)
-
- integer :: i,j,k,ixi,ieta,iga
- real :: v1, v2, v3, v4, v5, v6, v7, v8, xi, eta, ga, vi
- double precision :: zmesh
- real,parameter :: eps = 1.0e-3
-
-
- i = uc + 1
- j = vc + 1
- k = wc + 1
-
- xi = uc + 1 - i
- eta = vc + 1- j
- ga = wc + 1 -k
-
- ixi = nint(xi)
- ieta = nint(eta)
- iga = nint(ga)
-
-! print *, 'gc = ', i, j, k
-! print *, 'xi, eta, ga = ', xi, eta, ga
-! print *, 'ixi, ieta, iga = ', ixi, ieta, iga
-
-
- if (i > 0 .or. i < GOCAD_ST_NU .or. j > 0 .or. j < GOCAD_ST_NV .or. k > 0 .or. k < GOCAD_ST_NW) then
- v1 = vp_array(i,j,k)
- v2 = vp_array(i+1,j,k)
- v3 = vp_array(i+1,j+1,k)
- v4 = vp_array(i,j+1,k)
- v5 = vp_array(i,j,k+1)
- v6 = vp_array(i+1,j,k+1)
- v7 = vp_array(i+1,j+1,k+1)
- v8 = vp_array(i,j+1,k+1)
- vi = vp_array(i+ixi,j+ieta,k+iga)
-! print *, v1, v2, v3, v4, v5, v6, v7, v8
- if ((v1 - GOCAD_ST_NO_DATA_VALUE) > eps .and. &
- (v2 - GOCAD_ST_NO_DATA_VALUE) > eps .and. &
- (v3 - GOCAD_ST_NO_DATA_VALUE) > eps .and. &
- (v4 - GOCAD_ST_NO_DATA_VALUE) > eps .and. &
- (v5 - GOCAD_ST_NO_DATA_VALUE) > eps .and. &
- (v6 - GOCAD_ST_NO_DATA_VALUE) > eps .and. &
- (v7 - GOCAD_ST_NO_DATA_VALUE) > eps .and. &
- (v8 - GOCAD_ST_NO_DATA_VALUE) > eps ) then
- vp = dble(&
- v1 * (1-xi) * (1-eta) * (1-ga) +&
- v2 * xi * (1-eta) * (1-ga) +&
- v3 * xi * eta * (1-ga) +&
- v4 * (1-xi) * eta * (1-ga) +&
- v5 * (1-xi) * (1-eta) * ga +&
- v6 * xi * (1-eta) * ga +&
- v7 * xi * eta * ga +&
- v8 * (1-xi) * eta * ga)
- else if ((vi - GOCAD_ST_NO_DATA_VALUE) > eps) then
- vp = dble(vi)
-! else if ((v1 - GOCAD_ST_NO_DATA_VALUE) > eps) then
-! vp = dble(v1)
-! else if ((v2 - GOCAD_ST_NO_DATA_VALUE) > eps) then
-! vp = dble(v2)
-! else if ((v3 - GOCAD_ST_NO_DATA_VALUE) > eps) then
-! vp = dble(v3)
-! else if ((v4 - GOCAD_ST_NO_DATA_VALUE) > eps) then
-! vp = dble(v4)
-! else if ((v5 - GOCAD_ST_NO_DATA_VALUE) > eps) then
-! vp = dble(v5)
-! else if ((v6 - GOCAD_ST_NO_DATA_VALUE) > eps) then
-! vp = dble(v6)
-! else if ((v7 - GOCAD_ST_NO_DATA_VALUE) > eps) then
-! vp = dble(v7)
-! else if ((v7 - GOCAD_ST_NO_DATA_VALUE) > eps) then
-! vp = dble(v8)
- else
- vp = GOCAD_ST_NO_DATA_VALUE
- endif
- zmesh = wc / (GOCAD_ST_NW - 1) * GOCAD_ST_W_Z + GOCAD_ST_O_Z
- if (zmesh > -8500.) then
- vs = vp / (2 - (0.27*zmesh/(-8500)))
- else
- vs = vp/1.73
- endif
- if (vp > 2160.) then
- rho = vp/3 + 1280.
- else
- rho = 2000.
- endif
- else
- rho = GOCAD_ST_NO_DATA_VALUE
- vp = GOCAD_ST_NO_DATA_VALUE
- vs = GOCAD_ST_NO_DATA_VALUE
- endif
-
-end subroutine vx_xyz_interp
-
Deleted: seismo/3D/SPECFEM3D/trunk/socal_model.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/socal_model.f90 2010-09-24 22:15:18 UTC (rev 17214)
+++ seismo/3D/SPECFEM3D/trunk/socal_model.f90 2010-09-24 22:17:01 UTC (rev 17215)
@@ -1,63 +0,0 @@
-!=====================================================================
-!
-! S p e c f e m 3 D V e r s i o n 1 . 4
-! ---------------------------------------
-!
-! Dimitri Komatitsch and Jeroen Tromp
-! Seismological Laboratory - California Institute of Technology
-! (c) California Institute of Technology September 2006
-!
-! This program is free software; you can redistribute it and/or modify
-! it under the terms of the GNU General Public License as published by
-! the Free Software Foundation; either version 2 of the License, or
-! (at your option) any later version.
-!
-! This program is distributed in the hope that it will be useful,
-! but WITHOUT ANY WARRANTY; without even the implied warranty of
-! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-! GNU General Public License for more details.
-!
-! You should have received a copy of the GNU General Public License along
-! with this program; if not, write to the Free Software Foundation, Inc.,
-! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-!
-!=====================================================================
-
- subroutine socal_model(idoubling,rho,vp,vs)
-
- implicit none
-
- include "constants.h"
- include "constants_gocad.h"
-
- integer idoubling
- double precision rho,vp,vs
-
- if(idoubling == IFLAG_HALFSPACE_MOHO) then
- vp=7.8d0
- vs=4.5d0
- rho=3.0d0
-
- else if(idoubling == IFLAG_MOHO_16km) then
- vp=6.7d0
- vs=3.87d0
- rho=2.8d0
-
- else if(idoubling == IFLAG_ONE_LAYER_TOPOGRAPHY .or. idoubling == IFLAG_BASEMENT_TOPO) then
- vp=5.5d0
- vs=3.18d0
- rho=2.4d0
-
- else
- vp=6.3d0
- vs=3.64d0
- rho=2.67d0
- endif
-
-! scale to standard units
- vp = vp * 1000.d0
- vs = vs * 1000.d0
- rho = rho * 1000.d0
-
- end subroutine socal_model
-
More information about the CIG-COMMITS
mailing list