[cig-commits] r14505 - in seismo/3D/ADJOINT_TOMO/iterate_adj/pangu: . model_slice
carltape at geodynamics.org
carltape at geodynamics.org
Fri Mar 27 17:36:33 PDT 2009
Author: carltape
Date: 2009-03-27 17:36:31 -0700 (Fri, 27 Mar 2009)
New Revision: 14505
Added:
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/OUTPUT_FILES/
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/README
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/compile
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/constants.h
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/exit_mpi.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/exit_mpi.o
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/ftags
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/interp.bash
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/locate_receivers.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/mpif.h
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/precision.h
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/read_basin_topo_bathy_file.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/read_basin_topo_bathy_file.o
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/run.lsf
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/sem_model_slice
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/sem_model_slice.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/sem_model_slice.o
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/utm_geo.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/utm_geo.o
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/values_from_mesher.h
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_001_input
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_002_input
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_vb_m00_001.gmt
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_vb_m00_002.gmt
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_vb_m16_001.gmt
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_vb_m16_002.gmt
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_vs_m00_001.gmt
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_vs_m00_002.gmt
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_vs_m16_001.gmt
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_vs_m16_002.gmt
Log:
Adding a set of fortran codes written by Qinya Liu to read in a set of target points, and then extract the nearest GLL point from a pre-stored set of topology files describing a mesh. The motivation is to obtain high-resolution cross sections directly from the pre-stored models. These output values are then used within GMT plotting scripts.
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/README
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/README (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/README 2009-03-28 00:36:31 UTC (rev 14505)
@@ -0,0 +1,23 @@
+Qinya Liu, 30-Nov-2008
+
+This program is intended to compute the model values for given spatial points (as input), and the output file can be used directly in GMT plotting program such as grdimage to plot the model cross-section.
+
+Use an old version of the basin code that Carl uses for socal adjoint (subspace) tomo inversion.
+
+Includes an option (TOPOGRAPHY = .true.) for assigning negative values to 'air points' that are above the topography. This is useful for vertical cross sections.
+
+Max number of input points is set at NMAXPTS = 100000.
+
+--------------------------
+Carl Tape
+27-March-2009
+
+The example is executed by run.lsf and shown in interp.bash.
+
+There are two input files, vert_xc_001_input and vert_xc_002_input, each with two points.
+
+The output files are listed for four different models: vb_m00, vs_m00, vb_m16, vs_m16.
+
+Of course, in the future it would be best to pre-determine the target cross section points, and then compute (once-and-for-all) the nearest GLL point for each target point.
+
+==========================
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/compile
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/compile (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/compile 2009-03-28 00:36:31 UTC (rev 14505)
@@ -0,0 +1 @@
+mpif90 -O3 -o sem_model_slice sem_model_slice.f90 exit_mpi.f90 read_basin_topo_bathy_file.f90 utm_geo.f90
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/compile
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/constants.h
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/constants.h (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/constants.h 2009-03-28 00:36:31 UTC (rev 14505)
@@ -0,0 +1,373 @@
+!=====================================================================
+!
+! S p e c f e m 3 D B a s i n 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 July 2005
+!
+! A signed non-commercial agreement is required to use this program.
+! Please check http://www.gps.caltech.edu/research/jtromp for details.
+! Free for non-commercial academic research ONLY.
+! This program is distributed WITHOUT ANY WARRANTY whatsoever.
+! Do not redistribute this program without written permission.
+!
+!=====================================================================
+
+! constants.h. Generated from constants.h.in by configure.
+
+!
+! solver in single or double precision depending on the machine (4 or 8 bytes)
+!
+! ALSO CHANGE FILE precision.h ACCORDINGLY
+!
+ integer, parameter :: SIZE_REAL = 4
+ integer, parameter :: SIZE_DOUBLE = 8
+
+! set to SIZE_REAL to run in single precision
+! set to SIZE_DOUBLE to run in double precision (increases memory size by 2)
+ integer, parameter :: CUSTOM_REAL = SIZE_REAL
+
+!----------- parameters that can be changed by the user -----------
+
+! set to .false. if running on a Beowulf-type machine with local disks
+! set to .true. if running on a shared-memory machine with common file system
+! if running on a Beowulf, also modify name of nodes in filter_machine_file.f90
+ logical, parameter :: LOCAL_PATH_IS_ALSO_GLOBAL = .false.
+
+! apply heuristic rule to modify doubling regions to balance angles
+ logical, parameter :: APPLY_HEURISTIC_RULE = .true.
+
+! number of GLL points in each direction of an element (degree plus one)
+ integer, parameter :: NGLLX = 5
+ integer, parameter :: NGLLY = NGLLX
+ integer, parameter :: NGLLZ = NGLLX
+
+! input, output and main MPI I/O files
+ integer, parameter :: ISTANDARD_OUTPUT = 6
+ integer, parameter :: IIN = 40,IOUT = 41
+! uncomment this to write messages to a text file
+ integer, parameter :: IMAIN = 42
+! uncomment this to write messages to the screen
+! integer, parameter :: IMAIN = ISTANDARD_OUTPUT
+! I/O unit for source and receiver vtk file
+ integer, parameter :: IOVTK = 98
+
+! minimum thickness in meters to include the effect of the oceans
+! to avoid taking into account spurious oscillations in topography model
+ double precision, parameter :: MINIMUM_THICKNESS_3D_OCEANS = 10.d0
+
+! min and max density in the model (based on Brocher, 2005, eq. 1)
+ double precision, parameter :: DENSITY_MAX = 3476.d0 ! 3000.d0 or 3476.d0
+ double precision, parameter :: DENSITY_MIN = 1635.d0 ! 2000.d0 or 1635.d0
+
+! density of sea water
+ real(kind=CUSTOM_REAL), parameter :: RHO_OCEANS = 1020.0
+
+! extend basin model below threshold and above topography to make sure
+! there is no small gap between interpolated maps and sediments
+ logical, parameter :: EXTEND_VOXET_BELOW_BASEMENT = .true.
+ logical, parameter :: EXTEND_VOXET_ABOVE_TOPO = .true.
+ double precision, parameter :: DISTMAX_ASSUME_SEDIMENTS = 210.d0
+ integer, parameter :: NCELLS_EXTEND = 8
+
+! depth at which we start to honor the basement interface
+ double precision, parameter :: Z_THRESHOLD_HONOR_BASEMENT = -4700.d0
+
+! flag to print the details of source location
+ logical, parameter :: SHOW_DETAILS_LOCATE_SOURCE = .false.
+
+! maximum length of station and network name for receivers
+ integer, parameter :: MAX_LENGTH_STATION_NAME = 32
+ integer, parameter :: MAX_LENGTH_NETWORK_NAME = 8
+
+! number of sources to be gathered by MPI_Gather
+ integer, parameter :: NGATHER_SOURCES = 10000
+
+! source decay rate
+ double precision, parameter :: SOURCE_DECAY_RATE = 1.628d0
+
+! ---------------------------------------------------------------------------------------
+! LQY -- Following 3 variables stays here temporarily,
+! we need to move them to Par_file at a proper time
+! ---------------------------------------------------------------------------------------
+! save moho mesh and compute Moho boundary kernels
+ logical, parameter :: SAVE_MOHO_MESH = .true.
+
+! number of steps to save the state variables in the forward simulation,
+! to be used in the backward reconstruction in the presence of attenuation
+ integer, parameter :: NSTEP_Q_SAVE = 200
+
+! the scratch disk to save the state variables saved in the forward
+! simulation, this can be a global scratch disk in case you run out of
+! space on the local scratch disk
+ character(len=150), parameter :: LOCAL_PATH_Q = '/ibrixfs1/scratch/lqy'
+
+!------------------------------------------------------
+!----------- do not modify anything below -------------
+!------------------------------------------------------
+
+! on some processors (e.g. Pentiums) it is necessary to suppress underflows
+! by using a small initial field instead of zero
+ logical, parameter :: FIX_UNDERFLOW_PROBLEM = .true.
+
+! some useful constants
+ double precision, parameter :: PI = 3.141592653589793d0
+ double precision, parameter :: TWO_PI = 2.d0 * PI
+
+! 3-D simulation
+ integer, parameter :: NDIM = 3
+
+! dimension of the boundaries of the slices
+ integer, parameter :: NDIM2D = 2
+
+! number of nodes for 2D and 3D shape functions for hexahedra
+! we use 8-node mesh bricks, which are more stable than 27-node elements
+ integer, parameter :: NGNOD = 8, NGNOD2D = 4
+
+! a few useful constants
+ double precision, parameter :: ZERO = 0.d0,ONE = 1.d0,TWO = 2.d0,HALF = 0.5d0
+
+ real(kind=CUSTOM_REAL), parameter :: &
+ ONE_THIRD = 1._CUSTOM_REAL/3._CUSTOM_REAL, &
+ FOUR_THIRDS = 4._CUSTOM_REAL/3._CUSTOM_REAL
+
+! very large and very small values
+ double precision, parameter :: HUGEVAL = 1.d+30,TINYVAL = 1.d-9
+
+! very large real value declared independently of the machine
+ real(kind=CUSTOM_REAL), parameter :: HUGEVAL_SNGL = 1.e+30_CUSTOM_REAL
+
+! very large integer value
+ integer, parameter :: HUGEINT = 100000000
+
+! define flag for elements
+ integer, parameter :: IFLAG_ONE_LAYER_TOPOGRAPHY = 1
+ integer, parameter :: IFLAG_BASEMENT_TOPO = 2
+ integer, parameter :: IFLAG_16km_BASEMENT = 3
+ integer, parameter :: IFLAG_MOHO_16km = 4
+ integer, parameter :: IFLAG_HALFSPACE_MOHO = 5
+
+! define flag for regions for attenuation
+ integer, parameter :: NUM_REGIONS_ATTENUATION = 13
+
+ integer, parameter :: IATTENUATION_SEDIMENTS_40 = 1
+ integer, parameter :: IATTENUATION_SEDIMENTS_50 = 2
+ integer, parameter :: IATTENUATION_SEDIMENTS_60 = 3
+ integer, parameter :: IATTENUATION_SEDIMENTS_70 = 4
+ integer, parameter :: IATTENUATION_SEDIMENTS_80 = 5
+ integer, parameter :: IATTENUATION_SEDIMENTS_90 = 6
+ integer, parameter :: IATTENUATION_SEDIMENTS_100 = 7
+ integer, parameter :: IATTENUATION_SEDIMENTS_110 = 8
+ integer, parameter :: IATTENUATION_SEDIMENTS_120 = 9
+ integer, parameter :: IATTENUATION_SEDIMENTS_130 = 10
+ integer, parameter :: IATTENUATION_SEDIMENTS_140 = 11
+ integer, parameter :: IATTENUATION_SEDIMENTS_150 = 12
+ integer, parameter :: IATTENUATION_BEDROCK = 13
+
+! Olsen's constant for Q_mu = constant * v_s attenuation rule
+ real, parameter :: OLSEN_ATTENUATION_RATIO = 0.05
+
+! number of standard linear solids in parallel for attenuation
+ integer, parameter :: N_SLS = 3
+
+! flag for the four edges of each slice and for the bottom edge
+ integer, parameter :: XI_MIN = 1
+ integer, parameter :: XI_MAX = 2
+ integer, parameter :: ETA_MIN = 3
+ integer, parameter :: ETA_MAX = 4
+ integer, parameter :: BOTTOM = 5
+
+! number of points per surface element
+ integer, parameter :: NGLLSQUARE = NGLLX * NGLLY
+
+! number of points per spectral element
+ integer, parameter :: NGLLCUBE = NGLLX * NGLLY * NGLLZ
+
+! for vectorization of loops
+ integer, parameter :: NGLLSQUARE_NDIM = NGLLSQUARE * NDIM
+ integer, parameter :: NGLLCUBE_NDIM = NGLLCUBE * NDIM
+
+! flag for projection from latitude/longitude to UTM, and back
+ integer, parameter :: ILONGLAT2UTM = 0, IUTM2LONGLAT = 1
+
+! smallest real number on the Pentium and the SGI = 1.1754944E-38
+! largest real number on the Pentium and the SGI = 3.4028235E+38
+! small negligible initial value to avoid very slow underflow trapping
+! but not too small to avoid trapping on velocity and acceleration in Newmark
+ real(kind=CUSTOM_REAL), parameter :: VERYSMALLVAL = 1.E-24_CUSTOM_REAL
+
+! displacement threshold above which we consider the code became unstable
+ real(kind=CUSTOM_REAL), parameter :: STABILITY_THRESHOLD = 1.E+25_CUSTOM_REAL
+
+! geometrical tolerance for boundary detection
+ double precision, parameter :: SMALLVAL = 0.00001d0
+
+! do not use tags for MPI messages, use dummy tag instead
+ integer, parameter :: itag = 0,itag2 = 0
+
+! for the Gauss-Lobatto-Legendre points and weights
+ double precision, parameter :: GAUSSALPHA = 0.d0,GAUSSBETA = 0.d0
+
+! number of lines per source in CMTSOLUTION file
+ integer, parameter :: NLINES_PER_CMTSOLUTION_SOURCE = 13
+
+! number of iterations to solve the system for xi and eta
+ integer, parameter :: NUM_ITER = 4
+
+! size of topography and bathymetry file for Southern California
+ integer, parameter :: NX_TOPO_SOCAL = 1401,NY_TOPO_SOCAL = 1001
+ double precision, parameter :: ORIG_LAT_TOPO_SOCAL = 32.d0
+ double precision, parameter :: ORIG_LONG_TOPO_SOCAL = -121.d0
+ double precision, parameter :: DEGREES_PER_CELL_TOPO_SOCAL = 5.d0 / 1000.d0
+ character(len=100), parameter :: TOPO_FILE_SOCAL = 'DATA/la_topography/topo_bathy_final.dat'
+
+! size of topography and bathymetry file for Lacq (France) gas field
+ integer, parameter :: NX_TOPO_LACQ = 499,NY_TOPO_LACQ = 401
+ double precision, parameter :: ORIG_LAT_TOPO_LACQ = 100000.d0
+ double precision, parameter :: ORIG_LONG_TOPO_LACQ = 340400.d0
+ double precision, parameter :: DEGREES_PER_CELL_TOPO_LACQ = 100.d0
+ character(len=100), parameter :: TOPO_FILE_LACQ = 'DATA/lacq_thomas/mnt_Lacq_Lambert_final_dimitri.dat'
+
+! size of Lupei Zhu's Moho map file for Southern California
+ integer, parameter :: NX_MOHO = 71,NY_MOHO = 51
+ double precision, parameter :: ORIG_LAT_MOHO = 32.d0
+ double precision, parameter :: ORIG_LONG_MOHO = -121.d0
+ double precision, parameter :: DEGREES_PER_CELL_MOHO = 0.1d0
+
+! size of basement map file
+ integer, parameter :: NX_BASEMENT = 161,NY_BASEMENT = 144
+ double precision, parameter :: ORIG_X_BASEMENT = 316000.
+ double precision, parameter :: ORIG_Y_BASEMENT = 3655000.
+ double precision, parameter :: SPACING_X_BASEMENT = 1000.
+ double precision, parameter :: SPACING_Y_BASEMENT = 1000.
+
+!
+! new Gocad Voxets Peter July 29, 2002 - high-res and medium-res blocks
+!
+
+! size of the medium-resolution Gocad voxet
+ integer, parameter :: NX_GOCAD_MR = 194, NY_GOCAD_MR = 196, NZ_GOCAD_MR = 100
+
+ double precision, parameter :: ORIG_X_GOCAD_MR = 283000.
+ double precision, parameter :: ORIG_Y_GOCAD_MR = 3655000.
+ double precision, parameter :: ORIG_Z_GOCAD_MR = -15000.
+
+ double precision, parameter :: SPACING_X_GOCAD_MR = 1000.
+ double precision, parameter :: SPACING_Y_GOCAD_MR = 1000.
+ double precision, parameter :: SPACING_Z_GOCAD_MR = 200.
+
+! maximum size of model for tapering of transition between Hauksson and MR
+ double precision, parameter :: END_X_GOCAD_MR = ORIG_X_GOCAD_MR + SPACING_X_GOCAD_MR * (NX_GOCAD_MR - 1)
+ double precision, parameter :: END_Y_GOCAD_MR = ORIG_Y_GOCAD_MR + SPACING_Y_GOCAD_MR * (NY_GOCAD_MR - 1)
+
+! size of the high-resolution Gocad voxet
+ integer, parameter :: NX_GOCAD_HR = 185, NY_GOCAD_HR = 196, NZ_GOCAD_HR = 100
+
+ double precision, parameter :: ORIG_X_GOCAD_HR = 371052.25
+ double precision, parameter :: ORIG_Y_GOCAD_HR = 3725250.
+ double precision, parameter :: ORIG_Z_GOCAD_HR = -9500.
+
+ double precision, parameter :: SPACING_X_GOCAD_HR = 250.
+ double precision, parameter :: SPACING_Y_GOCAD_HR = 250.
+ double precision, parameter :: SPACING_Z_GOCAD_HR = 100.
+
+! maximum size of model for tapering of transition between HR and MR
+ double precision, parameter :: END_X_GOCAD_HR = ORIG_X_GOCAD_HR + SPACING_X_GOCAD_HR * (NX_GOCAD_HR - 1)
+ double precision, parameter :: END_Y_GOCAD_HR = ORIG_Y_GOCAD_HR + SPACING_Y_GOCAD_HR * (NY_GOCAD_HR - 1)
+
+! implement smooth transition between Hauksson, HR and MR Gocad blocks
+ logical, parameter :: TAPER_GOCAD_TRANSITIONS = .true.
+
+! Salton Sea Gocad voxet
+ integer, parameter :: GOCAD_ST_NU = 638, GOCAD_ST_NV = 219, GOCAD_ST_NW = 76
+ double precision, parameter :: GOCAD_ST_O_X = 720844.0, GOCAD_ST_O_Y = 3401799.250, &
+ GOCAD_ST_O_Z = -6354.334
+ double precision, parameter :: GOCAD_ST_U_X = -209197.89, GOCAD_ST_U_Y = 320741.71
+ double precision, parameter :: GOCAD_ST_V_X = 109670.74, GOCAD_ST_V_Y = 71530.72
+ double precision, parameter :: GOCAD_ST_W_Z = 7666.334
+ double precision, parameter :: GOCAD_ST_NO_DATA_VALUE = -99999
+
+!
+!--- larger Hauksson model for entire So-Cal, 15 km resolution
+!
+
+!! Hauksson (2000)
+!! number of non-constant layers
+! integer, parameter :: NLAYERS_HAUKSSON = 9
+!! depth of layers
+! double precision, parameter :: Z_HAUKSSON_LAYER_1 = -1000.d0
+! double precision, parameter :: Z_HAUKSSON_LAYER_2 = -4000.d0
+! double precision, parameter :: Z_HAUKSSON_LAYER_3 = -6000.d0
+! double precision, parameter :: Z_HAUKSSON_LAYER_4 = -10000.d0
+! double precision, parameter :: Z_HAUKSSON_LAYER_5 = -15000.d0
+! double precision, parameter :: Z_HAUKSSON_LAYER_6 = -17000.d0
+! double precision, parameter :: Z_HAUKSSON_LAYER_7 = -22000.d0
+! double precision, parameter :: Z_HAUKSSON_LAYER_8 = -31000.d0
+! double precision, parameter :: Z_HAUKSSON_LAYER_9 = -33000.d0
+!
+! integer, parameter :: NGRID_NEW_HAUKSSON = 201
+!
+!! corners of new Hauksson's interpolated grid
+! double precision, parameter :: UTM_X_ORIG_HAUKSSON = 122035.012d0
+! double precision, parameter :: UTM_X_END_HAUKSSON = 766968.628d0
+! double precision, parameter :: UTM_Y_ORIG_HAUKSSON = 3547232.986d0
+! double precision, parameter :: UTM_Y_END_HAUKSSON = 4098868.501d0
+
+! Lin-Shearer-Hauksson-Thurber (2007)
+! number of non-constant layers
+ integer, parameter :: NLAYERS_HAUKSSON = 8
+! depth of layers
+ double precision, parameter :: Z_HAUKSSON_LAYER_1 = 0.d0
+ double precision, parameter :: Z_HAUKSSON_LAYER_2 = -3000.d0
+ double precision, parameter :: Z_HAUKSSON_LAYER_3 = -6000.d0
+ double precision, parameter :: Z_HAUKSSON_LAYER_4 = -10000.d0
+ double precision, parameter :: Z_HAUKSSON_LAYER_5 = -15000.d0
+ double precision, parameter :: Z_HAUKSSON_LAYER_6 = -17000.d0
+ double precision, parameter :: Z_HAUKSSON_LAYER_7 = -22000.d0
+ double precision, parameter :: Z_HAUKSSON_LAYER_8 = -31000.d0
+
+ integer, parameter :: NGRID_NEW_HAUKSSON = 241
+
+! corners of new Hauksson's interpolated grid
+ double precision, parameter :: UTM_X_ORIG_HAUKSSON = 88021.568d0
+ double precision, parameter :: UTM_X_END_HAUKSSON = 861517.886d0
+ double precision, parameter :: UTM_Y_ORIG_HAUKSSON = 3404059.875d0
+ double precision, parameter :: UTM_Y_END_HAUKSSON = 4180234.582d0
+
+ double precision, parameter :: SPACING_UTM_X_HAUKSSON = (UTM_X_END_HAUKSSON - UTM_X_ORIG_HAUKSSON) / (NGRID_NEW_HAUKSSON-1.d0)
+ double precision, parameter :: SPACING_UTM_Y_HAUKSSON = (UTM_Y_END_HAUKSSON - UTM_Y_ORIG_HAUKSSON) / (NGRID_NEW_HAUKSSON-1.d0)
+
+! layers in the So-Cal regional model
+! DEPTH_MOHO_SOCAL = -35 km was based on Dreger and Helmberger (1990)
+! and is (July 2007) the preferred Moho depth for Dreger.
+! The depth of 32 km is used in the standard processing (Wald et al., 1995)
+! of SoCal events and is the value in the original Kanamori-Hadley (1975) model.
+ double precision, parameter :: DEPTH_5p5km_SOCAL = -5500.d0
+ double precision, parameter :: DEPTH_16km_SOCAL = -16000.d0
+ double precision, parameter :: DEPTH_MOHO_SOCAL = -32000.d0
+! double precision, parameter :: DEPTH_5p5km_SOCAL = -7700.d0
+! double precision, parameter :: DEPTH_16km_SOCAL = -18200.d0
+! double precision, parameter :: DEPTH_MOHO_SOCAL = -34200.d0
+
+! layer for Lacq (France) gas field
+ double precision, parameter :: DEPTH_INTERFACE_LACQ = -7000.d0
+
+! reference surface of the model before adding topography
+ double precision, parameter :: Z_SURFACE = 0.d0
+
+! number of points in each AVS or OpenDX quadrangular cell for movies
+ integer, parameter :: NGNOD2D_AVS_DX = 4
+
+! magic ratio for heuristic rule
+! this gives 120 degree angles in doubling
+! standard value 0.5 gives 135-135-90, which is not optimal
+ double precision, parameter :: MAGIC_RATIO = 0.6056d0
+
+! 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
+
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/exit_mpi.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/exit_mpi.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/exit_mpi.f90 2009-03-28 00:36:31 UTC (rev 14505)
@@ -0,0 +1,91 @@
+!=====================================================================
+!
+! S p e c f e m 3 D B a s i n 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
+!
+! A signed non-commercial agreement is required to use this program.
+! Please check http://www.gps.caltech.edu/research/jtromp for details.
+! Free for non-commercial academic research ONLY.
+! This program is distributed WITHOUT ANY WARRANTY whatsoever.
+! Do not redistribute this program without written permission.
+!
+!=====================================================================
+
+! end the simulation and exit MPI
+
+ subroutine exit_MPI(myrank,error_msg)
+
+ implicit none
+
+! standard include of the MPI library
+ include 'mpif.h'
+
+ include "constants.h"
+ include "precision.h"
+
+! identifier for error message file
+ integer, parameter :: IERROR = 30
+
+ integer myrank
+ character(len=*) error_msg
+
+ integer ier
+ character(len=80) outputname
+ character(len=150) OUTPUT_FILES
+
+! write error message to screen
+ write(*,*) error_msg(1:len(error_msg))
+ write(*,*) 'Error detected, aborting MPI... proc ',myrank
+
+! write error message to file
+ write(outputname,"('/error_message',i6.6,'.txt')") myrank
+ open(unit=IERROR,file='OUTPUT_FILES/'//outputname,status='unknown')
+ write(IERROR,*) error_msg(1:len(error_msg))
+ write(IERROR,*) 'Error detected, aborting MPI... proc ',myrank
+ close(IERROR)
+
+! close output file
+ if(myrank == 0 .and. IMAIN /= ISTANDARD_OUTPUT) close(IMAIN)
+
+! stop all the MPI processes, and exit
+! on some machines, MPI_FINALIZE needs to be called before MPI_ABORT
+ call MPI_FINALIZE(ier)
+ call MPI_ABORT(ier)
+ stop 'error, program ended in exit_MPI'
+
+ end subroutine exit_MPI
+
+!
+!----
+!
+
+! version without rank number printed in the error message
+ subroutine exit_MPI_without_rank(error_msg)
+
+ implicit none
+
+! standard include of the MPI library
+ include 'mpif.h'
+
+ include "constants.h"
+
+ character(len=*) error_msg
+
+ integer ier
+
+! write error message to screen
+ write(*,*) error_msg(1:len(error_msg))
+ write(*,*) 'Error detected, aborting MPI...'
+
+! stop all the MPI processes, and exit
+! on some machines, MPI_FINALIZE needs to be called before MPI_ABORT
+ call MPI_FINALIZE(ier)
+ call MPI_ABORT(ier)
+ stop 'error, program ended in exit_MPI'
+
+ end subroutine exit_MPI_without_rank
+
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/exit_mpi.o
===================================================================
(Binary files differ)
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/exit_mpi.o
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/ftags
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/ftags (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/ftags 2009-03-28 00:36:31 UTC (rev 14505)
@@ -0,0 +1,2 @@
+001
+002
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/interp.bash
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/interp.bash (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/interp.bash 2009-03-28 00:36:31 UTC (rev 14505)
@@ -0,0 +1,27 @@
+#!/bin/bash -v
+#BSUB -o OUTPUT_FILES/%J.o
+#BSUB -a mpich_gm
+#BSUB -J Model_slice_1
+
+echo "$LSB_MCPU_HOSTS" > OUTPUT_FILES/lsf_machines
+echo "$LSB_JOBID" > OUTPUT_FILES/jobid
+remap_lsf_machines.pl OUTPUT_FILES/lsf_machines >OUTPUT_FILES/machines
+
+for ftag in `cat ftags`; do
+
+#--------------------
+# vertical cross sections
+
+mpirun.lsf --gm-no-shmem --gm-copy-env sem_model_slice vert_xc_${ftag}_input /ibrixfs1/home/carltape/ADJOINT_TOMO_EXTRA/topo /ibrixfs1/home/carltape/ADJOINT_TOMO_EXTRA/KERNELS_MODELS/models/m16 vs_m16 vert_xc_vs_m16_${ftag}.gmt
+sleep 5s
+
+mpirun.lsf --gm-no-shmem --gm-copy-env sem_model_slice vert_xc_${ftag}_input /ibrixfs1/home/carltape/ADJOINT_TOMO_EXTRA/topo /ibrixfs1/home/carltape/ADJOINT_TOMO_EXTRA/KERNELS_MODELS/models/m00 vs_m00 vert_xc_vs_m00_${ftag}.gmt
+sleep 5s
+
+mpirun.lsf --gm-no-shmem --gm-copy-env sem_model_slice vert_xc_${ftag}_input /ibrixfs1/home/carltape/ADJOINT_TOMO_EXTRA/topo /ibrixfs1/home/carltape/ADJOINT_TOMO_EXTRA/KERNELS_MODELS/models/m16 vb_m16 vert_xc_vb_m16_${ftag}.gmt
+sleep 5s
+
+mpirun.lsf --gm-no-shmem --gm-copy-env sem_model_slice vert_xc_${ftag}_input /ibrixfs1/home/carltape/ADJOINT_TOMO_EXTRA/topo /ibrixfs1/home/carltape/ADJOINT_TOMO_EXTRA/KERNELS_MODELS/models/m00 vb_m00 vert_xc_vb_m00_${ftag}.gmt
+sleep 5s
+
+done
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/interp.bash
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/locate_receivers.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/locate_receivers.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/locate_receivers.f90 2009-03-28 00:36:31 UTC (rev 14505)
@@ -0,0 +1,664 @@
+!=====================================================================
+!
+! S p e c f e m 3 D B a s i n 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
+!
+! A signed non-commercial agreement is required to use this program.
+! Please check http://www.gps.caltech.edu/research/jtromp for details.
+! Free for non-commercial academic research ONLY.
+! This program is distributed WITHOUT ANY WARRANTY whatsoever.
+! Do not redistribute this program without written permission.
+!
+!=====================================================================
+
+!----
+!---- locate_receivers finds the correct position of the receivers
+!----
+ subroutine locate_receivers(ibool,myrank,NSPEC_AB,NGLOB_AB, &
+ xstore,ystore,zstore,xigll,yigll,zigll,rec_filename, &
+ nrec,islice_selected_rec,ispec_selected_rec, &
+ xi_receiver,eta_receiver,gamma_receiver,station_name,network_name,nu, &
+ NPROC,utm_x_source,utm_y_source, &
+ TOPOGRAPHY,itopo_bathy_basin,UTM_PROJECTION_ZONE,SUPPRESS_UTM_PROJECTION, &
+ NX_TOPO,NY_TOPO,ORIG_LAT_TOPO,ORIG_LONG_TOPO,DEGREES_PER_CELL_TOPO)
+
+ implicit none
+
+#ifdef USE_MPI
+! standard include of the MPI library
+ include 'mpif.h'
+#endif
+
+ include "constants.h"
+#ifdef USE_MPI
+ include "precision.h"
+#endif
+
+ integer NPROC,UTM_PROJECTION_ZONE,NX_TOPO,NY_TOPO
+
+ logical TOPOGRAPHY,SUPPRESS_UTM_PROJECTION
+
+ integer nrec,myrank
+
+ integer NSPEC_AB,NGLOB_AB
+
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
+
+ double precision ORIG_LAT_TOPO,ORIG_LONG_TOPO,DEGREES_PER_CELL_TOPO
+
+! arrays containing coordinates of the points
+ real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: xstore,ystore,zstore
+
+! use integer array to store topography values
+ integer itopo_bathy_basin(NX_TOPO,NY_TOPO)
+ double precision long_corner,lat_corner,ratio_xi,ratio_eta
+
+ integer, allocatable, dimension(:) :: ix_initial_guess,iy_initial_guess,iz_initial_guess
+
+ integer iprocloop
+ integer nrec_dummy
+ integer ios
+
+ double precision, allocatable, dimension(:) :: x_target,y_target,z_target
+ double precision, allocatable, dimension(:) :: horiz_dist,elevation
+ double precision, allocatable, dimension(:) :: x_found,y_found,z_found
+ double precision, allocatable, dimension(:,:) :: x_found_all,y_found_all,z_found_all
+
+ integer irec
+ integer i,j,k,ispec,iglob
+#ifdef USE_MPI
+ integer ier
+#endif
+
+ integer icornerlong,icornerlat
+ double precision utm_x_source,utm_y_source
+ double precision dist
+ double precision xi,eta,gamma,dx,dy,dz,dxi,deta,dgamma
+
+! Gauss-Lobatto-Legendre points of integration
+ double precision xigll(NGLLX)
+ double precision yigll(NGLLY)
+ double precision zigll(NGLLZ)
+
+! input receiver file name
+ character(len=*) rec_filename
+
+! topology of the control points of the surface element
+ integer iax,iay,iaz
+ integer iaddx(NGNOD),iaddy(NGNOD),iaddz(NGNOD)
+
+! coordinates of the control points of the surface element
+ double precision xelm(NGNOD),yelm(NGNOD),zelm(NGNOD)
+
+ integer iter_loop,ispec_iterate
+
+ integer ia
+ double precision x,y,z
+ double precision xix,xiy,xiz
+ double precision etax,etay,etaz
+ double precision gammax,gammay,gammaz
+
+! timer MPI
+ double precision time_start,tCPU
+
+! use dynamic allocation
+ double precision, dimension(:), allocatable :: final_distance
+ double precision, dimension(:,:), allocatable :: final_distance_all
+ double precision distmin,final_distance_max
+
+! receiver information
+! timing information for the stations
+! station information for writing the seismograms
+
+ integer, dimension(nrec) :: islice_selected_rec,ispec_selected_rec
+ double precision, dimension(nrec) :: xi_receiver,eta_receiver,gamma_receiver
+ double precision, dimension(3,3,nrec) :: nu
+ character(len=MAX_LENGTH_STATION_NAME), dimension(nrec) :: station_name
+ character(len=MAX_LENGTH_NETWORK_NAME), dimension(nrec) :: network_name
+
+ integer, allocatable, dimension(:,:) :: ispec_selected_rec_all
+ double precision, allocatable, dimension(:) :: stlat,stlon,stele,stbur,stutm_x,stutm_y
+ double precision, allocatable, dimension(:,:) :: xi_receiver_all,eta_receiver_all,gamma_receiver_all
+
+ character(len=150) OUTPUT_FILES
+
+! **************
+
+! get MPI starting time
+#ifdef USE_MPI
+ time_start = MPI_WTIME()
+#else
+ time_start = 0.d0
+#endif
+
+ if(myrank == 0) then
+ write(IMAIN,*)
+ write(IMAIN,*) '********************'
+ write(IMAIN,*) ' locating receivers'
+ write(IMAIN,*) '********************'
+ write(IMAIN,*)
+ endif
+
+! define topology of the control element
+ call usual_hex_nodes(iaddx,iaddy,iaddz)
+
+ if(myrank == 0) then
+ write(IMAIN,*)
+ write(IMAIN,*) '*****************************************************************'
+ write(IMAIN,'(1x,a,a,a)') 'reading receiver information from ', trim(rec_filename), ' file'
+ write(IMAIN,*) '*****************************************************************'
+ endif
+
+! get number of stations from receiver file
+ open(unit=1,file=trim(rec_filename),status='old',action='read',iostat=ios)
+ if (ios /= 0) call exit_mpi(myrank,'error opening file '//trim(rec_filename))
+ read(1,*) nrec_dummy
+
+ if(nrec_dummy /= nrec) call exit_MPI(myrank,'problem with number of receivers')
+
+! allocate memory for arrays using number of stations
+ allocate(stlat(nrec))
+ allocate(stlon(nrec))
+ allocate(stele(nrec))
+ allocate(stbur(nrec))
+ allocate(stutm_x(nrec))
+ allocate(stutm_y(nrec))
+ allocate(horiz_dist(nrec))
+ allocate(elevation(nrec))
+
+ allocate(ix_initial_guess(nrec))
+ allocate(iy_initial_guess(nrec))
+ allocate(iz_initial_guess(nrec))
+ allocate(x_target(nrec))
+ allocate(y_target(nrec))
+ allocate(z_target(nrec))
+ allocate(x_found(nrec))
+ allocate(y_found(nrec))
+ allocate(z_found(nrec))
+ allocate(final_distance(nrec))
+
+ allocate(ispec_selected_rec_all(nrec,0:NPROC-1))
+ allocate(xi_receiver_all(nrec,0:NPROC-1))
+ allocate(eta_receiver_all(nrec,0:NPROC-1))
+ allocate(gamma_receiver_all(nrec,0:NPROC-1))
+ allocate(x_found_all(nrec,0:NPROC-1))
+ allocate(y_found_all(nrec,0:NPROC-1))
+ allocate(z_found_all(nrec,0:NPROC-1))
+ allocate(final_distance_all(nrec,0:NPROC-1))
+
+! loop on all the stations
+ do irec=1,nrec
+
+! set distance to huge initial value
+ distmin=HUGEVAL
+
+ read(1,*,iostat=ios) station_name(irec),network_name(irec),stlat(irec),stlon(irec),stele(irec),stbur(irec)
+ if (ios /= 0) call exit_mpi(myrank, 'Error reading station file '//trim(rec_filename))
+
+! convert station location to UTM
+ call utm_geo(stlon(irec),stlat(irec),stutm_x(irec),stutm_y(irec),UTM_PROJECTION_ZONE,ILONGLAT2UTM,SUPPRESS_UTM_PROJECTION)
+
+! compute horizontal distance between source and receiver in km
+ horiz_dist(irec) = dsqrt((stutm_y(irec)-utm_y_source)**2 + (stutm_x(irec)-utm_x_source)**2) / 1000.
+
+! print some information about stations
+ if(myrank == 0) &
+ write(IMAIN,*) 'Station #',irec,': ',station_name(irec)(1:len_trim(station_name(irec))), &
+ '.',network_name(irec)(1:len_trim(network_name(irec))), &
+ ' horizontal distance: ',sngl(horiz_dist(irec)),' km'
+
+! get the Cartesian components of n in the model: nu
+
+! orientation consistent with the UTM projection
+
+! East
+ nu(1,1,irec) = 1.d0
+ nu(1,2,irec) = 0.d0
+ nu(1,3,irec) = 0.d0
+
+! North
+ nu(2,1,irec) = 0.d0
+ nu(2,2,irec) = 1.d0
+ nu(2,3,irec) = 0.d0
+
+! Vertical
+ nu(3,1,irec) = 0.d0
+ nu(3,2,irec) = 0.d0
+ nu(3,3,irec) = 1.d0
+
+! compute elevation of topography at the receiver location
+! we assume that receivers are always at the surface i.e. not buried
+ if(TOPOGRAPHY) then
+
+! get coordinate of corner in bathy/topo model
+ icornerlong = int((stlon(irec) - ORIG_LONG_TOPO) / DEGREES_PER_CELL_TOPO) + 1
+ icornerlat = int((stlat(irec) - ORIG_LAT_TOPO) / DEGREES_PER_CELL_TOPO) + 1
+
+! avoid edge effects and extend with identical point if outside model
+ if(icornerlong < 1) icornerlong = 1
+ if(icornerlong > NX_TOPO-1) icornerlong = NX_TOPO-1
+ if(icornerlat < 1) icornerlat = 1
+ if(icornerlat > NY_TOPO-1) icornerlat = NY_TOPO-1
+
+! compute coordinates of corner
+ long_corner = ORIG_LONG_TOPO + (icornerlong-1)*DEGREES_PER_CELL_TOPO
+ lat_corner = ORIG_LAT_TOPO + (icornerlat-1)*DEGREES_PER_CELL_TOPO
+
+! compute ratio for interpolation
+ ratio_xi = (stlon(irec) - long_corner) / DEGREES_PER_CELL_TOPO
+ ratio_eta = (stlat(irec) - lat_corner) / DEGREES_PER_CELL_TOPO
+
+! avoid edge effects
+ if(ratio_xi < 0.) ratio_xi = 0.
+ if(ratio_xi > 1.) ratio_xi = 1.
+ if(ratio_eta < 0.) ratio_eta = 0.
+ if(ratio_eta > 1.) ratio_eta = 1.
+
+! interpolate elevation at current point
+ elevation = &
+ itopo_bathy_basin(icornerlong,icornerlat)*(1.-ratio_xi)*(1.-ratio_eta) + &
+ itopo_bathy_basin(icornerlong+1,icornerlat)*ratio_xi*(1.-ratio_eta) + &
+ itopo_bathy_basin(icornerlong+1,icornerlat+1)*ratio_xi*ratio_eta + &
+ itopo_bathy_basin(icornerlong,icornerlat+1)*(1.-ratio_xi)*ratio_eta
+
+ else
+ elevation(irec) = 0.d0
+ endif
+
+! compute the Cartesian position of the receiver
+ x_target(irec) = stutm_x(irec)
+ y_target(irec) = stutm_y(irec)
+ z_target(irec) = elevation(irec) - stbur(irec)
+ if (myrank == 0) write(IOVTK,*) x_target(irec), y_target(irec), z_target(irec)
+
+! examine top of the elements only (receivers always at the surface)
+! k = NGLLZ
+
+ do ispec=1,NSPEC_AB
+
+! modification by Qinya Liu: idoubling is no longer used because receivers can now be in depth
+! (in previous versions, receivers were always assumed to be at the surface)
+
+! loop only on points inside the element
+! exclude edges to ensure this point is not shared with other elements
+ do k=2,NGLLZ-1
+ do j=2,NGLLY-1
+ do i=2,NGLLX-1
+
+ iglob = ibool(i,j,k,ispec)
+ dist = dsqrt((x_target(irec)-dble(xstore(iglob)))**2 &
+ +(y_target(irec)-dble(ystore(iglob)))**2 &
+ +(z_target(irec)-dble(zstore(iglob)))**2)
+
+! keep this point if it is closer to the receiver
+ if(dist < distmin) then
+ distmin = dist
+ ispec_selected_rec(irec) = ispec
+ ix_initial_guess(irec) = i
+ iy_initial_guess(irec) = j
+ iz_initial_guess(irec) = k
+ endif
+
+ enddo
+ enddo
+ enddo
+! endif
+
+! end of loop on all the spectral elements in current slice
+ enddo
+
+! end of loop on all the stations
+ enddo
+
+! close receiver file
+ close(1)
+
+! ****************************************
+! find the best (xi,eta,gamma) for each receiver
+! ****************************************
+
+! loop on all the receivers to iterate in that slice
+ do irec = 1,nrec
+
+ ispec_iterate = ispec_selected_rec(irec)
+
+! use initial guess in xi and eta
+ xi = xigll(ix_initial_guess(irec))
+ eta = yigll(iy_initial_guess(irec))
+ gamma = zigll(iz_initial_guess(irec))
+
+! define coordinates of the control points of the element
+
+ do ia=1,NGNOD
+
+ if(iaddx(ia) == 0) then
+ iax = 1
+ else if(iaddx(ia) == 1) then
+ iax = (NGLLX+1)/2
+ else if(iaddx(ia) == 2) then
+ iax = NGLLX
+ else
+ call exit_MPI(myrank,'incorrect value of iaddx')
+ endif
+
+ if(iaddy(ia) == 0) then
+ iay = 1
+ else if(iaddy(ia) == 1) then
+ iay = (NGLLY+1)/2
+ else if(iaddy(ia) == 2) then
+ iay = NGLLY
+ else
+ call exit_MPI(myrank,'incorrect value of iaddy')
+ endif
+
+ if(iaddz(ia) == 0) then
+ iaz = 1
+ else if(iaddz(ia) == 1) then
+ iaz = (NGLLZ+1)/2
+ else if(iaddz(ia) == 2) then
+ iaz = NGLLZ
+ else
+ call exit_MPI(myrank,'incorrect value of iaddz')
+ endif
+
+ iglob = ibool(iax,iay,iaz,ispec_iterate)
+ xelm(ia) = dble(xstore(iglob))
+ yelm(ia) = dble(ystore(iglob))
+ zelm(ia) = dble(zstore(iglob))
+
+ enddo
+
+! iterate to solve the non linear system
+ do iter_loop = 1,NUM_ITER
+
+! impose receiver exactly at the surface
+! gamma = 1.d0
+
+! recompute jacobian for the new point
+ call recompute_jacobian(xelm,yelm,zelm,xi,eta,gamma,x,y,z, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz)
+
+! compute distance to target location
+ dx = - (x - x_target(irec))
+ dy = - (y - y_target(irec))
+ dz = - (z - z_target(irec))
+
+! compute increments
+! gamma does not change since we know the receiver is exactly on the surface
+ dxi = xix*dx + xiy*dy + xiz*dz
+ deta = etax*dx + etay*dy + etaz*dz
+ dgamma = gammax*dx + gammay*dy + gammaz*dz
+
+! update values
+ xi = xi + dxi
+ eta = eta + deta
+ gamma = gamma + dgamma
+
+! impose that we stay in that element
+! (useful if user gives a receiver outside the mesh for instance)
+! we can go slightly outside the [1,1] segment since with finite elements
+! the polynomial solution is defined everywhere
+! this can be useful for convergence of itertive scheme with distorted elements
+ if (xi > 1.10d0) xi = 1.10d0
+ if (xi < -1.10d0) xi = -1.10d0
+ if (eta > 1.10d0) eta = 1.10d0
+ if (eta < -1.10d0) eta = -1.10d0
+ if (gamma > 1.10d0) gamma = 1.10d0
+ if (gamma < -1.10d0) gamma = -1.10d0
+
+! end of non linear iterations
+ enddo
+
+! impose receiver exactly at the surface after final iteration
+! gamma = 1.d0
+
+! compute final coordinates of point found
+ call recompute_jacobian(xelm,yelm,zelm,xi,eta,gamma,x,y,z, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz)
+
+! store xi,eta and x,y,z of point found
+ xi_receiver(irec) = xi
+ eta_receiver(irec) = eta
+ gamma_receiver(irec) = gamma
+ x_found(irec) = x
+ y_found(irec) = y
+ z_found(irec) = z
+
+! compute final distance between asked and found (converted to km)
+ final_distance(irec) = dsqrt((x_target(irec)-x_found(irec))**2 + &
+ (y_target(irec)-y_found(irec))**2 + (z_target(irec)-z_found(irec))**2)
+
+ enddo
+
+! synchronize all the processes to make sure all the estimates are available
+#ifdef USE_MPI
+ call MPI_BARRIER(MPI_COMM_WORLD,ier)
+#endif
+
+! for MPI version, gather information from all the nodes
+ ispec_selected_rec_all(:,:) = -1
+#ifdef USE_MPI
+ call MPI_GATHER(ispec_selected_rec,nrec,MPI_INTEGER,ispec_selected_rec_all,nrec,MPI_INTEGER,0,MPI_COMM_WORLD,ier)
+
+ call MPI_GATHER(xi_receiver,nrec,MPI_DOUBLE_PRECISION,xi_receiver_all,nrec,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+ call MPI_GATHER(eta_receiver,nrec,MPI_DOUBLE_PRECISION,eta_receiver_all,nrec,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+ call MPI_GATHER(gamma_receiver,nrec,MPI_DOUBLE_PRECISION,gamma_receiver_all,nrec,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+ call MPI_GATHER(final_distance,nrec,MPI_DOUBLE_PRECISION,final_distance_all,nrec,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+ call MPI_GATHER(x_found,nrec,MPI_DOUBLE_PRECISION,x_found_all,nrec,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+ call MPI_GATHER(y_found,nrec,MPI_DOUBLE_PRECISION,y_found_all,nrec,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+ call MPI_GATHER(z_found,nrec,MPI_DOUBLE_PRECISION,z_found_all,nrec,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+#else
+ ispec_selected_rec_all(:,0) = ispec_selected_rec
+
+ xi_receiver_all(:,0) = xi_receiver(:)
+ eta_receiver_all(:,0) = eta_receiver(:)
+ gamma_receiver_all(:,0) = gamma_receiver(:)
+ final_distance_all(:,0) = final_distance(:)
+ x_found_all(:,0) = x_found(:)
+ y_found_all(:,0) = y_found(:)
+ z_found_all(:,0) = z_found(:)
+#endif
+
+! this is executed by main process only
+ if(myrank == 0) then
+
+! check that the gather operation went well
+ if(any(ispec_selected_rec_all(:,:) == -1)) call exit_MPI(myrank,'gather operation failed for receivers')
+
+! MPI loop on all the results to determine the best slice
+ islice_selected_rec(:) = -1
+ do irec = 1,nrec
+ distmin = HUGEVAL
+ do iprocloop = 0,NPROC-1
+ if(final_distance_all(irec,iprocloop) < distmin) then
+ distmin = final_distance_all(irec,iprocloop)
+ islice_selected_rec(irec) = iprocloop
+ ispec_selected_rec(irec) = ispec_selected_rec_all(irec,iprocloop)
+ xi_receiver(irec) = xi_receiver_all(irec,iprocloop)
+ eta_receiver(irec) = eta_receiver_all(irec,iprocloop)
+ gamma_receiver(irec) = gamma_receiver_all(irec,iprocloop)
+ x_found(irec) = x_found_all(irec,iprocloop)
+ y_found(irec) = y_found_all(irec,iprocloop)
+ z_found(irec) = z_found_all(irec,iprocloop)
+ endif
+ enddo
+ final_distance(irec) = distmin
+ enddo
+
+ do irec=1,nrec
+
+ write(IMAIN,*)
+ write(IMAIN,*) 'station # ',irec,' ',station_name(irec),network_name(irec)
+
+ if(final_distance(irec) == HUGEVAL) call exit_MPI(myrank,'error locating receiver')
+
+ write(IMAIN,*) ' original latitude: ',sngl(stlat(irec))
+ write(IMAIN,*) ' original longitude: ',sngl(stlon(irec))
+ write(IMAIN,*) ' original UTM x: ',sngl(stutm_x(irec))
+ write(IMAIN,*) ' original UTM y: ',sngl(stutm_y(irec))
+ write(IMAIN,*) ' horizontal distance: ',sngl(horiz_dist(irec))
+ if(TOPOGRAPHY) write(IMAIN,*) ' topography elevation: ',sngl(elevation(irec))
+ write(IMAIN,*) ' target x, y, z: ',sngl(x_target(irec)),sngl(y_target(irec)),sngl(z_target(irec))
+
+ write(IMAIN,*) 'closest estimate found: ',sngl(final_distance(irec)),' m away'
+ write(IMAIN,*) ' in slice ',islice_selected_rec(irec),' in element ',ispec_selected_rec(irec)
+ write(IMAIN,*) ' at xi,eta,gamma coordinates = ',xi_receiver(irec),eta_receiver(irec),gamma_receiver(irec)
+
+! add warning if estimate is poor
+! (usually means receiver outside the mesh given by the user)
+ if(final_distance(irec) > 3000.d0) then
+ write(IMAIN,*) '*******************************************************'
+ write(IMAIN,*) '***** WARNING: receiver location estimate is poor *****'
+ write(IMAIN,*) '*******************************************************'
+ endif
+
+ write(IMAIN,*)
+
+ enddo
+
+! compute maximal distance for all the receivers
+ final_distance_max = maxval(final_distance(:))
+
+! display maximum error for all the receivers
+ write(IMAIN,*) 'maximum error in location of all the receivers: ',sngl(final_distance_max),' m'
+
+! add warning if estimate is poor
+! (usually means receiver outside the mesh given by the user)
+ if(final_distance_max > 1000.d0) then
+ write(IMAIN,*)
+ write(IMAIN,*) '************************************************************'
+ write(IMAIN,*) '************************************************************'
+ write(IMAIN,*) '***** WARNING: at least one receiver is poorly located *****'
+ write(IMAIN,*) '************************************************************'
+ write(IMAIN,*) '************************************************************'
+ endif
+
+! get the base pathname for output files
+ call get_value_string(OUTPUT_FILES, 'OUTPUT_FILES', 'OUTPUT_FILES')
+
+! write the list of stations and associated epicentral distance
+ open(unit=27,file=trim(OUTPUT_FILES)//'/output_list_stations.txt',status='unknown')
+ do irec=1,nrec
+ write(27,*) station_name(irec),'.',network_name(irec),' : ',horiz_dist(irec),' km horizontal distance'
+ enddo
+ close(27)
+
+! elapsed time since beginning of mesh generation
+#ifdef USE_MPI
+ tCPU = MPI_WTIME() - time_start
+#else
+ tCPU = 0.d0
+#endif
+ write(IMAIN,*)
+ write(IMAIN,*) 'Elapsed time for receiver detection in seconds = ',tCPU
+ write(IMAIN,*)
+ write(IMAIN,*) 'End of receiver detection - done'
+ write(IMAIN,*)
+
+ endif ! end of section executed by main process only
+
+#ifdef USE_MPI
+! main process broadcasts the results to all the slices
+ call MPI_BCAST(islice_selected_rec,nrec,MPI_INTEGER,0,MPI_COMM_WORLD,ier)
+ call MPI_BCAST(ispec_selected_rec,nrec,MPI_INTEGER,0,MPI_COMM_WORLD,ier)
+ call MPI_BCAST(xi_receiver,nrec,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+ call MPI_BCAST(eta_receiver,nrec,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+ call MPI_BCAST(gamma_receiver,nrec,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+
+! synchronize all the processes to make sure everybody has finished
+ call MPI_BARRIER(MPI_COMM_WORLD,ier)
+#endif
+
+! deallocate arrays
+ deallocate(stlat)
+ deallocate(stlon)
+ deallocate(stele)
+ deallocate(stbur)
+ deallocate(stutm_x)
+ deallocate(stutm_y)
+ deallocate(horiz_dist)
+ deallocate(ix_initial_guess)
+ deallocate(iy_initial_guess)
+ deallocate(iz_initial_guess)
+ deallocate(x_target)
+ deallocate(y_target)
+ deallocate(z_target)
+ deallocate(x_found)
+ deallocate(y_found)
+ deallocate(z_found)
+ deallocate(final_distance)
+ deallocate(ispec_selected_rec_all)
+ deallocate(xi_receiver_all)
+ deallocate(eta_receiver_all)
+ deallocate(gamma_receiver_all)
+ deallocate(x_found_all)
+ deallocate(y_found_all)
+ deallocate(z_found_all)
+ deallocate(final_distance_all)
+
+ end subroutine locate_receivers
+
+!===========================
+
+ subroutine station_filter(myrank,filename,filtered_filename,nfilter, &
+ LATITUDE_MIN, LATITUDE_MAX, LONGITUDE_MIN, LONGITUDE_MAX)
+
+ implicit none
+
+ include 'constants.h'
+
+! input
+ integer :: myrank
+ character(len=*) filename,filtered_filename
+ double precision LATITUDE_MIN,LATITUDE_MAX,LONGITUDE_MIN,LONGITUDE_MAX
+
+! output
+ integer nfilter
+
+ integer irec, nrec, nrec_filtered, ios
+
+ double precision stlat,stlon,stele,stbur
+ character(len=MAX_LENGTH_STATION_NAME) station_name
+ character(len=MAX_LENGTH_NETWORK_NAME) network_name
+
+ nrec_filtered = 0
+
+ open(unit=IIN, file=trim(filename), status = 'old', iostat = ios)
+ if (ios /= 0) call exit_mpi(myrank, 'No file '//trim(filename)//', exit')
+ read(IIN, *) nrec
+ do irec = 1, nrec
+ read(IIN, *) station_name, network_name, stlat, stlon, stele, stbur
+ if(stlat >= LATITUDE_MIN .and. stlat <= LATITUDE_MAX .and. stlon >= LONGITUDE_MIN .and. stlon <= LONGITUDE_MAX) &
+ nrec_filtered = nrec_filtered + 1
+ enddo
+ close(IIN)
+
+ if (myrank == 0) then
+ open(unit=IIN,file=trim(filename),status='old',action='read')
+ open(unit=IOUT,file=trim(filtered_filename),status='unknown')
+ read(IIN,*) nrec
+ write(IOUT,*) nrec_filtered
+ do irec = 1,nrec
+ read(IIN,*) station_name,network_name,stlat,stlon,stele,stbur
+ if(stlat >= LATITUDE_MIN .and. stlat <= LATITUDE_MAX .and. stlon >= LONGITUDE_MIN .and. stlon <= LONGITUDE_MAX) &
+ write(IOUT,*) station_name,' ',network_name,' ',sngl(stlat),' ',sngl(stlon), ' ',sngl(stele), ' ',sngl(stbur)
+ enddo
+ close(IIN)
+ close(IOUT)
+
+ write(IMAIN,*)
+ write(IMAIN,*) 'there are ',nrec,' stations in file ', trim(filename)
+ write(IMAIN,*) 'saving ',nrec_filtered,' stations inside the model in file ', trim(filtered_filename)
+ write(IMAIN,*) 'excluding ',nrec - nrec_filtered,' stations located outside the model'
+ write(IMAIN,*)
+
+ endif
+
+ nfilter = nrec_filtered
+
+ end subroutine station_filter
+
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/mpif.h
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/mpif.h (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/mpif.h 2009-03-28 00:36:31 UTC (rev 14505)
@@ -0,0 +1,265 @@
+! mpif.h. Generated from /opt/mpich/myrinet/intel//include/mpif.h by configure.
+
+!
+!
+! (C) 1993 by Argonne National Laboratory and Mississipi State University.
+! All rights reserved. See COPYRIGHT in top-level directory.
+!
+!
+! user include file for MPI programs, with no dependencies
+!
+! It really is not possible to make a perfect include file that can
+! be used by both F77 and F90 compilers, but this is close. We have removed
+! continuation lines (allows free form input in F90); systems whose
+! Fortran compilers support ! instead of just C or * for comments can
+! globally replace a C in the first column with !; the resulting file
+! should work for both Fortran 77 and Fortran 90.
+!
+! If your Fortran compiler supports ! for comments, you can run this
+! through sed with
+! sed -e 's/^C/\!/g'
+!
+! We have also removed the use of contractions (involving the single quote)
+! character because some users use .F instead of .f files (to invoke the
+! cpp preprocessor) and further, their preprocessor is determined to find
+! matching single quote pairs (and probably double quotes; given the
+! different rules in C and Fortran, this sounds like a disaster). Rather than
+! take the position that the poor users should get a better system, we
+! have removed the text that caused problems. Of course, the users SHOULD
+! get a better system...
+!
+! return codes
+ INTEGER MPI_SUCCESS,MPI_ERR_BUFFER,MPI_ERR_COUNT,MPI_ERR_TYPE
+ INTEGER MPI_ERR_TAG,MPI_ERR_COMM,MPI_ERR_RANK,MPI_ERR_ROOT
+ INTEGER MPI_ERR_GROUP
+ INTEGER MPI_ERR_OP,MPI_ERR_TOPOLOGY,MPI_ERR_DIMS,MPI_ERR_ARG
+ INTEGER MPI_ERR_UNKNOWN,MPI_ERR_TRUNCATE,MPI_ERR_OTHER
+ INTEGER MPI_ERR_INTERN,MPI_ERR_IN_STATUS,MPI_ERR_PENDING
+ INTEGER MPI_ERR_REQUEST, MPI_ERR_LASTCODE
+ PARAMETER (MPI_SUCCESS=0,MPI_ERR_BUFFER=1,MPI_ERR_COUNT=2)
+ PARAMETER (MPI_ERR_TYPE=3,MPI_ERR_TAG=4,MPI_ERR_COMM=5)
+ PARAMETER (MPI_ERR_RANK=6,MPI_ERR_ROOT=7,MPI_ERR_GROUP=8)
+ PARAMETER (MPI_ERR_OP=9,MPI_ERR_TOPOLOGY=10,MPI_ERR_DIMS=11)
+ PARAMETER (MPI_ERR_ARG=12,MPI_ERR_UNKNOWN=13)
+ PARAMETER (MPI_ERR_TRUNCATE=14,MPI_ERR_OTHER=15)
+ PARAMETER (MPI_ERR_INTERN=16,MPI_ERR_IN_STATUS=17)
+ PARAMETER (MPI_ERR_PENDING=18,MPI_ERR_REQUEST=19)
+ PARAMETER (MPI_ERR_LASTCODE=1073741823)
+!
+ INTEGER MPI_UNDEFINED
+ parameter (MPI_UNDEFINED = (-32766))
+!
+ INTEGER MPI_GRAPH, MPI_CART
+ PARAMETER (MPI_GRAPH = 1, MPI_CART = 2)
+ INTEGER MPI_PROC_NULL
+ PARAMETER ( MPI_PROC_NULL = (-1) )
+!
+ INTEGER MPI_BSEND_OVERHEAD
+ PARAMETER ( MPI_BSEND_OVERHEAD = 512 )
+
+ INTEGER MPI_SOURCE, MPI_TAG, MPI_ERROR
+ PARAMETER(MPI_SOURCE=2, MPI_TAG=3, MPI_ERROR=4)
+ INTEGER MPI_STATUS_SIZE
+ PARAMETER (MPI_STATUS_SIZE=4)
+ INTEGER MPI_MAX_PROCESSOR_NAME, MPI_MAX_ERROR_STRING
+ PARAMETER (MPI_MAX_PROCESSOR_NAME=256)
+ PARAMETER (MPI_MAX_ERROR_STRING=512)
+ INTEGER MPI_MAX_NAME_STRING
+ PARAMETER (MPI_MAX_NAME_STRING=63)
+ INTEGER MPI_MAX_PORT_NAME
+ PARAMETER (MPI_MAX_PORT_NAME=256)
+!
+ INTEGER MPI_COMM_NULL
+ PARAMETER (MPI_COMM_NULL=0)
+!
+ INTEGER MPI_DATATYPE_NULL
+ PARAMETER (MPI_DATATYPE_NULL = 0)
+
+ INTEGER MPI_ERRHANDLER_NULL
+ PARAMETER (MPI_ERRHANDLER_NULL = 0)
+
+ INTEGER MPI_GROUP_NULL
+ PARAMETER (MPI_GROUP_NULL = 0)
+
+ INTEGER MPI_KEYVAL_INVALID
+ PARAMETER (MPI_KEYVAL_INVALID = 0)
+
+ INTEGER MPI_REQUEST_NULL
+ PARAMETER (MPI_REQUEST_NULL = 0)
+!
+ INTEGER MPI_IDENT, MPI_CONGRUENT, MPI_SIMILAR, MPI_UNEQUAL
+ PARAMETER (MPI_IDENT=0, MPI_CONGRUENT=1, MPI_SIMILAR=2)
+ PARAMETER (MPI_UNEQUAL=3)
+!
+! MPI_BOTTOM needs to be a known address; here we put it at the
+! beginning of the common block. The point-to-point and collective
+! routines know about MPI_BOTTOM, but MPI_TYPE_STRUCT as yet does not.
+!
+! MPI_STATUS_IGNORE and MPI_STATUSES_IGNORE are similar objects
+! Until the underlying MPI library implements the C version of these
+! (a null pointer), these are declared as arrays of MPI_STATUS_SIZE
+!
+! The types MPI_INTEGER1,2,4 and MPI_REAL4,8 are OPTIONAL.
+! Their values are zero if they are not available. Note that
+! using these reduces the portability of code (though may enhance
+! portability between Crays and other systems)
+!
+ INTEGER MPI_TAG_UB, MPI_HOST, MPI_IO
+ INTEGER MPI_BOTTOM
+ INTEGER MPI_STATUS_IGNORE(MPI_STATUS_SIZE)
+ INTEGER MPI_STATUSES_IGNORE(MPI_STATUS_SIZE)
+ INTEGER MPI_INTEGER, MPI_REAL, MPI_DOUBLE_PRECISION
+ INTEGER MPI_COMPLEX, MPI_DOUBLE_COMPLEX,MPI_LOGICAL
+ INTEGER MPI_CHARACTER, MPI_BYTE, MPI_2INTEGER, MPI_2REAL
+ INTEGER MPI_2DOUBLE_PRECISION, MPI_2COMPLEX, MPI_2DOUBLE_COMPLEX
+ INTEGER MPI_UB, MPI_LB
+ INTEGER MPI_PACKED, MPI_WTIME_IS_GLOBAL
+ INTEGER MPI_COMM_WORLD, MPI_COMM_SELF, MPI_GROUP_EMPTY
+ INTEGER MPI_SUM, MPI_MAX, MPI_MIN, MPI_PROD, MPI_LAND, MPI_BAND
+ INTEGER MPI_LOR, MPI_BOR, MPI_LXOR, MPI_BXOR, MPI_MINLOC
+ INTEGER MPI_MAXLOC
+ INTEGER MPI_OP_NULL
+ INTEGER MPI_ERRORS_ARE_FATAL, MPI_ERRORS_RETURN
+!
+ PARAMETER (MPI_ERRORS_ARE_FATAL=119)
+ PARAMETER (MPI_ERRORS_RETURN=120)
+!
+ PARAMETER (MPI_COMPLEX=23,MPI_DOUBLE_COMPLEX=24,MPI_LOGICAL=25)
+ PARAMETER (MPI_REAL=26,MPI_DOUBLE_PRECISION=27,MPI_INTEGER=28)
+ PARAMETER (MPI_2INTEGER=29,MPI_2COMPLEX=30,MPI_2DOUBLE_COMPLEX=31)
+ PARAMETER (MPI_2REAL=32,MPI_2DOUBLE_PRECISION=33,MPI_CHARACTER=1)
+ PARAMETER (MPI_BYTE=3,MPI_UB=16,MPI_LB=15,MPI_PACKED=14)
+
+ INTEGER MPI_ORDER_C, MPI_ORDER_FORTRAN
+ PARAMETER (MPI_ORDER_C=56, MPI_ORDER_FORTRAN=57)
+ INTEGER MPI_DISTRIBUTE_BLOCK, MPI_DISTRIBUTE_CYCLIC
+ INTEGER MPI_DISTRIBUTE_NONE, MPI_DISTRIBUTE_DFLT_DARG
+ PARAMETER (MPI_DISTRIBUTE_BLOCK=121, MPI_DISTRIBUTE_CYCLIC=122)
+ PARAMETER (MPI_DISTRIBUTE_NONE=123)
+ PARAMETER (MPI_DISTRIBUTE_DFLT_DARG=-49767)
+ INTEGER MPI_MAX_INFO_KEY, MPI_MAX_INFO_VAL
+ PARAMETER (MPI_MAX_INFO_KEY=255, MPI_MAX_INFO_VAL=1024)
+ INTEGER MPI_INFO_NULL
+ PARAMETER (MPI_INFO_NULL=0)
+
+!
+! Optional Fortran Types. Configure attempts to determine these.
+!
+ INTEGER MPI_INTEGER1, MPI_INTEGER2, MPI_INTEGER4, MPI_INTEGER8
+ INTEGER MPI_INTEGER16
+ INTEGER MPI_REAL4, MPI_REAL8, MPI_REAL16
+ INTEGER MPI_COMPLEX8, MPI_COMPLEX16, MPI_COMPLEX32
+ PARAMETER (MPI_INTEGER1=1,MPI_INTEGER2=4)
+ PARAMETER (MPI_INTEGER4=6)
+ PARAMETER (MPI_INTEGER8=8)
+ PARAMETER (MPI_INTEGER16=0)
+ PARAMETER (MPI_REAL4=10)
+ PARAMETER (MPI_REAL8=11)
+ PARAMETER (MPI_REAL16=12)
+ PARAMETER (MPI_COMPLEX8=23)
+ PARAMETER (MPI_COMPLEX16=24)
+ PARAMETER (MPI_COMPLEX32=0)
+!
+! This is now handled with either the "pointer" extension or this same
+! code, appended at the end.
+! COMMON /MPIPRIV/ MPI_BOTTOM,MPI_STATUS_IGNORE,MPI_STATUSES_IGNORE
+!C
+!C Without this save, some Fortran implementations may make the common
+!C dynamic!
+!C
+!C For a Fortran90 module, we might replace /MPIPRIV/ with a simple
+!C SAVE MPI_BOTTOM
+!C
+! SAVE /MPIPRIV/
+!
+ PARAMETER (MPI_MAX=100,MPI_MIN=101,MPI_SUM=102,MPI_PROD=103)
+ PARAMETER (MPI_LAND=104,MPI_BAND=105,MPI_LOR=106,MPI_BOR=107)
+ PARAMETER (MPI_LXOR=108,MPI_BXOR=109,MPI_MINLOC=110)
+ PARAMETER (MPI_MAXLOC=111, MPI_OP_NULL=0)
+!
+ PARAMETER (MPI_GROUP_EMPTY=90,MPI_COMM_WORLD=91,MPI_COMM_SELF=92)
+ PARAMETER (MPI_TAG_UB=80,MPI_HOST=82,MPI_IO=84)
+ PARAMETER (MPI_WTIME_IS_GLOBAL=86)
+!
+ INTEGER MPI_ANY_SOURCE
+ PARAMETER (MPI_ANY_SOURCE = (-2))
+ INTEGER MPI_ANY_TAG
+ PARAMETER (MPI_ANY_TAG = (-1))
+!
+ INTEGER MPI_VERSION, MPI_SUBVERSION
+ PARAMETER (MPI_VERSION = 1, MPI_SUBVERSION = 2)
+!
+! There are additional MPI-2 constants
+ INTEGER MPI_ADDRESS_KIND, MPI_OFFSET_KIND
+ PARAMETER (MPI_ADDRESS_KIND=8)
+ PARAMETER (MPI_OFFSET_KIND=8)
+!
+! All other MPI routines are subroutines
+! This may cause some Fortran compilers to complain about defined and
+! not used. Such compilers should be improved.
+!
+! Some Fortran compilers will not link programs that contain
+! external statements to routines that are not provided, even if
+! the routine is never called. Remove PMPI_WTIME and PMPI_WTICK
+! if you have trouble with them.
+!
+ DOUBLE PRECISION MPI_WTIME, MPI_WTICK,PMPI_WTIME,PMPI_WTICK
+ EXTERNAL MPI_WTIME, MPI_WTICK,PMPI_WTIME,PMPI_WTICK
+!
+! The attribute copy/delete subroutines are symbols that can be passed
+! to MPI routines
+!
+ EXTERNAL MPI_NULL_COPY_FN, MPI_NULL_DELETE_FN, MPI_DUP_FN
+ COMMON /MPIPRIV/ MPI_BOTTOM,MPI_STATUS_IGNORE,MPI_STATUSES_IGNORE
+!
+! Without this save, some Fortran implementations may make the common
+! dynamic!
+!
+! For a Fortran90 module, we might replace /MPIPRIV/ with a simple
+! SAVE MPI_BOTTOM
+!
+ SAVE /MPIPRIV/
+!
+! $Id: mpiof.h.in,v 1.3 1999/08/06 18:33:09 thakur Exp $
+!
+! Copyright (C) 1997 University of Chicago.
+! See COPYRIGHT notice in top-level directory.
+!
+!
+! user include file for Fortran MPI-IO programs
+!
+ INTEGER MPI_MODE_RDONLY, MPI_MODE_RDWR, MPI_MODE_WRONLY
+ INTEGER MPI_MODE_DELETE_ON_CLOSE, MPI_MODE_UNIQUE_OPEN
+ INTEGER MPI_MODE_CREATE, MPI_MODE_EXCL
+ INTEGER MPI_MODE_APPEND, MPI_MODE_SEQUENTIAL
+ PARAMETER (MPI_MODE_RDONLY=2, MPI_MODE_RDWR=8, MPI_MODE_WRONLY=4)
+ PARAMETER (MPI_MODE_CREATE=1, MPI_MODE_DELETE_ON_CLOSE=16)
+ PARAMETER (MPI_MODE_UNIQUE_OPEN=32, MPI_MODE_EXCL=64)
+ PARAMETER (MPI_MODE_APPEND=128, MPI_MODE_SEQUENTIAL=256)
+!
+ INTEGER MPI_FILE_NULL
+ PARAMETER (MPI_FILE_NULL=0)
+!
+ INTEGER MPI_MAX_DATAREP_STRING
+ PARAMETER (MPI_MAX_DATAREP_STRING=128)
+!
+ INTEGER MPI_SEEK_SET, MPI_SEEK_CUR, MPI_SEEK_END
+ PARAMETER (MPI_SEEK_SET=600, MPI_SEEK_CUR=602, MPI_SEEK_END=604)
+!
+ INTEGER MPIO_REQUEST_NULL
+ PARAMETER (MPIO_REQUEST_NULL=0)
+!
+!
+!
+
+
+
+
+
+
+
+!
+!
+!
+!
+!
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/precision.h
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/precision.h (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/precision.h 2009-03-28 00:36:31 UTC (rev 14505)
@@ -0,0 +1,29 @@
+!=====================================================================
+!
+! S p e c f e m 3 D B a s i n 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 July 2005
+!
+! A signed non-commercial agreement is required to use this program.
+! Please check http://www.gps.caltech.edu/research/jtromp for details.
+! Free for non-commercial academic research ONLY.
+! This program is distributed WITHOUT ANY WARRANTY whatsoever.
+! Do not redistribute this program without written permission.
+!
+!=====================================================================
+
+! precision.h. Generated from precision.h.in by configure.
+
+!
+! solver in single or double precision depending on the machine
+!
+! set to MPI_REAL to run in single precision
+! set to MPI_DOUBLE_PRECISION to run in double precision
+!
+! ALSO CHANGE FILE constants.h ACCORDINGLY
+!
+ integer, parameter :: CUSTOM_MPI_TYPE = MPI_REAL
+ integer, parameter :: CUSTOM_MPI_2REAL = MPI_2REAL
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/read_basin_topo_bathy_file.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/read_basin_topo_bathy_file.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/read_basin_topo_bathy_file.f90 2009-03-28 00:36:31 UTC (rev 14505)
@@ -0,0 +1,46 @@
+!=====================================================================
+!
+! S p e c f e m 3 D B a s i n 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
+!
+! A signed non-commercial agreement is required to use this program.
+! Please check http://www.gps.caltech.edu/research/jtromp for details.
+! Free for non-commercial academic research ONLY.
+! This program is distributed WITHOUT ANY WARRANTY whatsoever.
+! Do not redistribute this program without written permission.
+!
+!=====================================================================
+
+ subroutine read_basin_topo_bathy_file(itopo_bathy_basin,NX_TOPO,NY_TOPO,topo_file)
+!
+!---- read basin topography and bathymetry file once and for all
+!
+ implicit none
+
+ include "constants.h"
+
+ integer NX_TOPO,NY_TOPO
+
+! use integer array to store topography values
+ integer itopo_bathy_basin(NX_TOPO,NY_TOPO)
+
+ character(len=100) topo_file
+
+ integer ix,iy
+
+ itopo_bathy_basin(:,:) = 0
+
+ open(unit=13,file=topo_file,status='old',action='read')
+ do iy=1,NY_TOPO
+ do ix=1,NX_TOPO
+ read(13,*) itopo_bathy_basin(ix,iy)
+ enddo
+ enddo
+ close(13)
+
+ end subroutine read_basin_topo_bathy_file
+
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/read_basin_topo_bathy_file.o
===================================================================
(Binary files differ)
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/read_basin_topo_bathy_file.o
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/run.lsf
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/run.lsf (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/run.lsf 2009-03-28 00:36:31 UTC (rev 14505)
@@ -0,0 +1,9 @@
+#!/bin/bash
+compile
+mkdir -p OUTPUT_FILES
+rm -f OUTPUT_FILES/*
+date
+#bsub -q debug -W 30 -n 168 < interp.bash
+#bsub -q normal -W 60 -n 168 < interp.bash
+bsub -q normal -W 1440 -n 168 < interp.bash
+date
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/run.lsf
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/sem_model_slice
===================================================================
(Binary files differ)
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/sem_model_slice
___________________________________________________________________
Name: svn:executable
+ *
Name: svn:mime-type
+ application/octet-stream
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/sem_model_slice.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/sem_model_slice.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/sem_model_slice.f90 2009-03-28 00:36:31 UTC (rev 14505)
@@ -0,0 +1,230 @@
+program sem_model_slice
+
+ implicit none
+
+ include 'mpif.h'
+ include 'constants.h'
+ include "precision.h"
+ include 'values_from_mesher.h'
+
+ integer, parameter :: NMAXPTS = 100000
+ integer :: ier,sizeprocs,myrank,ios, i,j, k, ispec,iglob,ipt, npts
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
+ character(len=150) :: xyz_infile,topo_dir,model_dir,data_name,gmt_outfile,&
+ local_data_file, prname
+ real(kind=CUSTOM_REAL),dimension(NMAXPTS) :: x, y, z, v, vmin, &
+ vall,distmin, dist, distall
+ integer,dimension(NMAXPTS) :: ispec_min,ix_min,iy_min,iz_min
+ real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: vstore
+ real(kind=CUSTOM_REAL),dimension(NGLOB_AB) :: xstore,ystore,zstore
+ real(kind=CUSTOM_REAL),dimension(2,NMAXPTS) :: in, out
+
+ ! true --> replace "air" points with NaN (vertical cross sections)
+ ! false --> take the closest value to the "air" points (horizontal cross section)
+ logical, parameter :: TOPOGRAPHY = .true.
+
+ character(len=100) :: topo_file
+ integer, dimension(NX_TOPO_SOCAL,NY_TOPO_SOCAL) :: itopo_bathy_basin
+ double precision :: elevation
+
+ ! MPI initialization
+ call MPI_INIT(ier)
+ call MPI_COMM_SIZE(MPI_COMM_WORLD,sizeprocs,ier)
+ call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ier)
+
+ ! input arguments
+ call getarg(1,xyz_infile)
+ call getarg(2,topo_dir)
+ call getarg(3,model_dir)
+ call getarg(4,data_name)
+ call getarg(5,gmt_outfile)
+
+ ! read points to be interpolated
+ open(11,file=xyz_infile,iostat=ios)
+ i=0
+ do while (1 == 1)
+ i=i+1
+ read(11,*,iostat=ios) x(i),y(i),z(i)
+ if (ios /= 0) exit
+ enddo
+ close(11)
+ npts=i-1
+ if (myrank == 0) then
+ if (npts > NMAXPTS .or. npts < 0) call exit_mpi(myrank,'Npts error ...')
+ write(*,*) 'Total number of points = ', npts
+ endif
+
+ ! read data and topo files
+ write(prname,'(a,i6.6,a)') trim(model_dir)//'/proc',myrank,'_'
+ local_data_file = trim(prname) // trim(data_name)//'.bin'
+ open(unit = 27,file=local_data_file,status='old', iostat = ios,form ='unformatted')
+ if (ios /= 0) call exit_mpi(myrank, 'Error reading model data file')
+ read(27) vstore
+ close(27)
+
+ write(prname,'(a,i6.6,a)') trim(topo_dir)//'/proc',myrank,'_'
+ open(unit = 27,file = trim(prname)//'x.bin',status='old',action='read', iostat = ios,form ='unformatted')
+ if (ios /= 0) call exit_mpi(myrank,'Error reading local x file')
+ read(27) xstore
+ close(27)
+ open(unit = 27,file = trim(prname)//'y.bin',status='old',action='read', iostat = ios,form ='unformatted')
+ if (ios /= 0) call exit_mpi(myrank,'Error reading local y file')
+ read(27) ystore
+ close(27)
+ open(unit = 27,file = trim(prname)//'z.bin',status='old',action='read', iostat = ios,form ='unformatted')
+ if (ios /= 0) call exit_mpi(myrank,'Error reading local z file')
+ read(27) zstore
+ close(27)
+ open(unit = 27,file = trim(prname)//'ibool.bin',status='old',action='read', iostat = ios,form ='unformatted')
+ if (ios /= 0) call exit_mpi(myrank,'Error reading local z file')
+ read(27) ibool
+ close(27)
+
+
+ ! search for local minimum-distance point
+ distmin(1:npts) = HUGEVAL
+
+ do ispec=1,NSPEC_AB
+
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec)
+ dist(1:npts)=dsqrt((x(1:npts)-dble(xstore(iglob)))**2 &
+ +(y(1:npts)-dble(ystore(iglob)))**2 &
+ +(z(1:npts)-dble(zstore(iglob)))**2)
+ do ipt=1,npts
+ if(dist(ipt) < distmin(ipt)) then
+ distmin(ipt)=dist(ipt)
+ ispec_min(ipt)=ispec
+ ix_min(ipt)=i; iy_min(ipt)=j; iz_min(ipt)=k
+ vmin(ipt)=vstore(i,j,k,ispec)
+ endif
+ enddo
+
+ enddo
+ enddo
+ enddo
+
+ ! end of loop on all the elements in current slice
+ enddo
+
+ call MPI_BARRIER(MPI_COMM_WORLD,ier)
+ if (myrank == 0) print *, 'Done looping over global points ...'
+
+ ! choose the minimum value
+
+ ! write(prname,'(a,i6.6,a)') 'OUTPUT_FILES/in_',myrank,'.txt'
+ ! open(33,file=prname)
+ ! do i = 1, npts
+ ! write(33,*) i, myrank, distmin(i), ispec_min(i), ix_min(i), iy_min(i), iz_min(i),vmin(i)
+ ! enddo
+ ! close(33)
+
+ do i=1, npts
+ in(1,i) = distmin(i)
+ in(2,i) = myrank ! myrank is coerced to a double
+ enddo
+ call MPI_REDUCE(in,out,npts,CUSTOM_MPI_2REAL,MPI_MINLOC,0,MPI_COMM_WORLD,ier)
+
+ ! if (myrank == 0) then
+ ! open(33,file='OUTPUT_FILES/out.txt')
+ ! do i = 1, npts
+ ! write(33,*) i, out(1,i), out(2,i)
+ ! enddo
+ ! close(33)
+ ! endif
+
+ call MPI_BCAST(out,2*npts,CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
+
+ v(1:npts) = 0
+ dist(1:npts) = 0.
+
+ do i = 1, npts
+ if (myrank == nint(out(2,i))) then
+ v(i) = vmin(i)
+! if (GLL_INTERPOLATION) call xeg_search(x(i),y(i),z(i),&
+! ispec_min(i),ix_min(i),iy_min(i),iz_min(i),v(i))
+ dist(i) = distmin(i)
+ endif
+ enddo
+
+ call MPI_REDUCE(v,vall,npts,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call MPI_REDUCE(dist,distall,npts,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+
+ if (myrank == 0) then
+
+ if (TOPOGRAPHY) then
+ topo_file='/ibrixfs1/home/lqy/cmt3d/test_dir/'//trim(TOPO_FILE_SOCAL)
+ call read_basin_topo_bathy_file(itopo_bathy_basin,NX_TOPO_SOCAL,NY_TOPO_SOCAL,topo_file)
+ endif
+
+ print *, 'Writing out gmt file ...'
+ open(12,file=gmt_outfile,status='unknown')
+ do i = 1, npts
+ if (TOPOGRAPHY) then
+ call topo_value(itopo_bathy_basin,dble(x(i)),dble(y(i)),elevation)
+ if (elevation < z(i)) vall(i)= -1000.00
+ endif
+ write(12,*) x(i), y(i), z(i), vall(i), distall(i)
+ enddo
+ close(12)
+ endif
+
+ call MPI_FINALIZE(ier)
+
+end program sem_model_slice
+
+!------------------------------------------------
+
+subroutine topo_value(itopo_bathy_basin,x,y,elevation)
+
+ implicit none
+ include 'constants.h'
+
+ integer,dimension(NX_TOPO_SOCAL,NY_TOPO_SOCAL) :: itopo_bathy_basin
+ double precision :: x, y
+ double precision elevation
+
+ double precision :: long, lat
+ integer :: icornerlong, icornerlat
+ double precision :: long_corner, lat_corner,ratio_xi,ratio_eta
+ integer, parameter :: UTM_PROJECTION_ZONE = 11
+ logical, parameter :: SUPPRESS_UTM_PROJECTION = .false.
+
+
+ call utm_geo(long,lat,x,y,UTM_PROJECTION_ZONE,IUTM2LONGLAT,SUPPRESS_UTM_PROJECTION)
+
+ ! get coordinate of corner in bathy/topo model
+ icornerlong = int((long - ORIG_LONG_TOPO_SOCAL) / DEGREES_PER_CELL_TOPO_SOCAL) + 1
+ icornerlat = int((lat - ORIG_LAT_TOPO_SOCAL) / DEGREES_PER_CELL_TOPO_SOCAL) + 1
+
+ ! avoid edge effects and extend with identical point if outside model
+ if(icornerlong < 1) icornerlong = 1
+ if(icornerlong > NX_TOPO_SOCAL-1) icornerlong = NX_TOPO_SOCAL-1
+ if(icornerlat < 1) icornerlat = 1
+ if(icornerlat > NY_TOPO_SOCAL-1) icornerlat = NY_TOPO_SOCAL-1
+
+ ! compute coordinates of corner
+ long_corner = ORIG_LONG_TOPO_SOCAL + (icornerlong-1)*DEGREES_PER_CELL_TOPO_SOCAL
+ lat_corner = ORIG_LAT_TOPO_SOCAL + (icornerlat-1)*DEGREES_PER_CELL_TOPO_SOCAL
+
+ ! compute ratio for interpolation
+ ratio_xi = (long - long_corner) / DEGREES_PER_CELL_TOPO_SOCAL
+ ratio_eta = (lat - lat_corner) / DEGREES_PER_CELL_TOPO_SOCAL
+
+ ! avoid edge effects
+ if(ratio_xi < 0.) ratio_xi = 0.
+ if(ratio_xi > 1.) ratio_xi = 1.
+ if(ratio_eta < 0.) ratio_eta = 0.
+ if(ratio_eta > 1.) ratio_eta = 1.
+
+ ! interpolate elevation at current point
+ elevation = &
+ itopo_bathy_basin(icornerlong,icornerlat)*(1.-ratio_xi)*(1.-ratio_eta) + &
+ itopo_bathy_basin(icornerlong+1,icornerlat)*ratio_xi*(1.-ratio_eta) + &
+ itopo_bathy_basin(icornerlong+1,icornerlat+1)*ratio_xi*ratio_eta + &
+ itopo_bathy_basin(icornerlong,icornerlat+1)*(1.-ratio_xi)*ratio_eta
+
+
+end subroutine topo_value
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/sem_model_slice.o
===================================================================
(Binary files differ)
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/sem_model_slice.o
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/utm_geo.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/utm_geo.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/utm_geo.f90 2009-03-28 00:36:31 UTC (rev 14505)
@@ -0,0 +1,198 @@
+!=====================================================================
+!
+! UTM (Universal Transverse Mercator) projection from the USGS
+!
+!=====================================================================
+
+ subroutine utm_geo(rlon,rlat,rx,ry,UTM_PROJECTION_ZONE,iway,SUPPRESS_UTM_PROJECTION)
+
+! convert geodetic longitude and latitude to UTM, and back
+! use iway = ILONGLAT2UTM for long/lat to UTM, IUTM2LONGLAT for UTM to lat/long
+! a list of UTM zones of the world is available at www.dmap.co.uk/utmworld.htm
+
+ implicit none
+
+ include "constants.h"
+
+!
+!-----CAMx v2.03
+!
+! UTM_GEO performs UTM to geodetic (long/lat) translation, and back.
+!
+! This is a Fortran version of the BASIC program "Transverse Mercator
+! Conversion", Copyright 1986, Norman J. Berls (Stefan Musarra, 2/94)
+! Based on algorithm taken from "Map Projections Used by the USGS"
+! by John P. Snyder, Geological Survey Bulletin 1532, USDI.
+!
+! Input/Output arguments:
+!
+! rlon Longitude (deg, negative for West)
+! rlat Latitude (deg)
+! rx UTM easting (m)
+! ry UTM northing (m)
+! UTM_PROJECTION_ZONE UTM zone
+! iway Conversion type
+! ILONGLAT2UTM = geodetic to UTM
+! IUTM2LONGLAT = UTM to geodetic
+!
+
+ integer UTM_PROJECTION_ZONE,iway
+ double precision rx,ry,rlon,rlat
+ logical SUPPRESS_UTM_PROJECTION
+
+ double precision, parameter :: degrad=PI/180.d0, raddeg=180.d0/PI
+ double precision, parameter :: semimaj=6378206.4d0, semimin=6356583.8d0
+ double precision, parameter :: scfa=0.9996d0
+ double precision, parameter :: north=0.d0, east=500000.d0
+
+ double precision e2,e4,e6,ep2,xx,yy,dlat,dlon,zone,cm,cmr,delam
+ double precision f1,f2,f3,f4,rm,rn,t,c,a,e1,u,rlat1,dlat1,c1,t1,rn1,r1,d
+ double precision rx_save,ry_save,rlon_save,rlat_save
+
+ if(SUPPRESS_UTM_PROJECTION) then
+ if (iway == ILONGLAT2UTM) then
+ rx = rlon
+ ry = rlat
+ else
+ rlon = rx
+ rlat = ry
+ endif
+ return
+ endif
+
+! save original parameters
+ rlon_save = rlon
+ rlat_save = rlat
+ rx_save = rx
+ ry_save = ry
+
+! define parameters of reference ellipsoid
+ e2=1.0-(semimin/semimaj)**2.0
+ e4=e2*e2
+ e6=e2*e4
+ ep2=e2/(1.-e2)
+
+ if (iway == IUTM2LONGLAT) then
+ xx = rx
+ yy = ry
+ else
+ dlon = rlon
+ dlat = rlat
+ endif
+!
+!----- Set Zone parameters
+!
+ zone = dble(UTM_PROJECTION_ZONE)
+ cm = zone*6.0 - 183.0
+ cmr = cm*degrad
+!
+!---- Lat/Lon to UTM conversion
+!
+ if (iway == ILONGLAT2UTM) then
+
+ rlon = degrad*dlon
+ rlat = degrad*dlat
+
+ delam = dlon - cm
+ if (delam < -180.) delam = delam + 360.
+ if (delam > 180.) delam = delam - 360.
+ delam = delam*degrad
+
+ f1 = (1. - e2/4. - 3.*e4/64. - 5.*e6/256)*rlat
+ f2 = 3.*e2/8. + 3.*e4/32. + 45.*e6/1024.
+ f2 = f2*sin(2.*rlat)
+ f3 = 15.*e4/256.*45.*e6/1024.
+ f3 = f3*sin(4.*rlat)
+ f4 = 35.*e6/3072.
+ f4 = f4*sin(6.*rlat)
+ rm = semimaj*(f1 - f2 + f3 - f4)
+ if (dlat == 90. .or. dlat == -90.) then
+ xx = 0.
+ yy = scfa*rm
+ else
+ rn = semimaj/sqrt(1. - e2*sin(rlat)**2)
+ t = tan(rlat)**2
+ c = ep2*cos(rlat)**2
+ a = cos(rlat)*delam
+
+ f1 = (1. - t + c)*a**3/6.
+ f2 = 5. - 18.*t + t**2 + 72.*c - 58.*ep2
+ f2 = f2*a**5/120.
+ xx = scfa*rn*(a + f1 + f2)
+ f1 = a**2/2.
+ f2 = 5. - t + 9.*c + 4.*c**2
+ f2 = f2*a**4/24.
+ f3 = 61. - 58.*t + t**2 + 600.*c - 330.*ep2
+ f3 = f3*a**6/720.
+ yy = scfa*(rm + rn*tan(rlat)*(f1 + f2 + f3))
+ endif
+ xx = xx + east
+ yy = yy + north
+
+!
+!---- UTM to Lat/Lon conversion
+!
+ else
+
+ xx = xx - east
+ yy = yy - north
+ e1 = sqrt(1. - e2)
+ e1 = (1. - e1)/(1. + e1)
+ rm = yy/scfa
+ u = 1. - e2/4. - 3.*e4/64. - 5.*e6/256.
+ u = rm/(semimaj*u)
+
+ f1 = 3.*e1/2. - 27.*e1**3./32.
+ f1 = f1*sin(2.*u)
+ f2 = 21.*e1**2/16. - 55.*e1**4/32.
+ f2 = f2*sin(4.*u)
+ f3 = 151.*e1**3./96.
+ f3 = f3*sin(6.*u)
+ rlat1 = u + f1 + f2 + f3
+ dlat1 = rlat1*raddeg
+ if (dlat1 >= 90. .or. dlat1 <= -90.) then
+ dlat1 = dmin1(dlat1,dble(90.) )
+ dlat1 = dmax1(dlat1,dble(-90.) )
+ dlon = cm
+ else
+ c1 = ep2*cos(rlat1)**2.
+ t1 = tan(rlat1)**2.
+ f1 = 1. - e2*sin(rlat1)**2.
+ rn1 = semimaj/sqrt(f1)
+ r1 = semimaj*(1. - e2)/sqrt(f1**3)
+ d = xx/(rn1*scfa)
+
+ f1 = rn1*tan(rlat1)/r1
+ f2 = d**2/2.
+ f3 = 5.*3.*t1 + 10.*c1 - 4.*c1**2 - 9.*ep2
+ f3 = f3*d**2*d**2/24.
+ f4 = 61. + 90.*t1 + 298.*c1 + 45.*t1**2. - 252.*ep2 - 3.*c1**2
+ f4 = f4*(d**2)**3./720.
+ rlat = rlat1 - f1*(f2 - f3 + f4)
+ dlat = rlat*raddeg
+
+ f1 = 1. + 2.*t1 + c1
+ f1 = f1*d**2*d/6.
+ f2 = 5. - 2.*c1 + 28.*t1 - 3.*c1**2 + 8.*ep2 + 24.*t1**2.
+ f2 = f2*(d**2)**2*d/120.
+ rlon = cmr + (d - f1 + f2)/cos(rlat1)
+ dlon = rlon*raddeg
+ if (dlon < -180.) dlon = dlon + 360.
+ if (dlon > 180.) dlon = dlon - 360.
+ endif
+ endif
+
+ if (iway == IUTM2LONGLAT) then
+ rlon = dlon
+ rlat = dlat
+ rx = rx_save
+ ry = ry_save
+ else
+ rx = xx
+ ry = yy
+ rlon = rlon_save
+ rlat = rlat_save
+ endif
+
+ end subroutine utm_geo
+
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/utm_geo.o
===================================================================
(Binary files differ)
Property changes on: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/utm_geo.o
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/values_from_mesher.h
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/values_from_mesher.h (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/values_from_mesher.h 2009-03-28 00:36:31 UTC (rev 14505)
@@ -0,0 +1,54 @@
+
+ !
+ ! this is the parameter file for static compilation of the solver
+ !
+ ! mesh statistics:
+ ! ---------------
+ !
+ ! number of processors = 168
+ !
+ ! number of ES nodes = 21.00000
+ ! percentage of total 640 ES nodes = 3.281250 %
+ ! total memory available on these ES nodes (Gb) = 336.0000
+ !
+ ! max points per processor = max vector length = 163941
+ ! min vector length = 25
+ ! min critical vector length = 75
+ !
+ ! on ES and SX-5, make sure "loopcnt=" parameter
+ ! in Makefile is greater than 163941
+ !
+ ! total elements per AB slice = 2412
+ ! total points per AB slice = 163941
+ !
+ ! total for full mesh:
+ ! -------------------
+ !
+ ! exact total number of spectral elements in entire mesh =
+ ! 405216
+ ! approximate total number of points in entire mesh =
+ ! 27542088.0000000
+ ! approximate total number of degrees of freedom in entire mesh =
+ ! 82626264.0000000
+ !
+ ! resolution of the mesh at the surface:
+ ! -------------------------------------
+ !
+ ! spectral elements along X = 336
+ ! spectral elements along Y = 288
+ ! GLL points along X = 1345
+ ! GLL points along Y = 1153
+ ! average distance between points along X in m = 475.4175
+ ! average distance between points along Y in m = 436.8476
+ !
+
+ integer, parameter :: NSPEC_AB = 2412
+ integer, parameter :: NGLOB_AB = 163941
+
+ !
+ ! number of time steps = 18200
+ !
+ integer, parameter :: NSPEC_ATTENUATION = NSPEC_AB
+ logical, parameter :: ATTENUATION_VAL = .true.
+ logical, parameter :: ANISOTROPY_VAL = .false.
+
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_001_input
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_001_input (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_001_input 2009-03-28 00:36:31 UTC (rev 14505)
@@ -0,0 +1,2 @@
+ 2.273952e+05 4.024866e+06 -4.500000e+04 -2.000000e+05
+ 2.273952e+05 4.024866e+06 -4.449495e+04 -2.000000e+05
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_002_input
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_002_input (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_002_input 2009-03-28 00:36:31 UTC (rev 14505)
@@ -0,0 +1,2 @@
+ 2.807261e+05 4.053913e+06 -4.500000e+04 -2.000000e+05
+ 2.807261e+05 4.053913e+06 -4.449495e+04 -2.000000e+05
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_vb_m00_001.gmt
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_vb_m00_001.gmt (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_vb_m00_001.gmt 2009-03-28 00:36:31 UTC (rev 14505)
@@ -0,0 +1,2 @@
+ 227395.2 4024866. -45000.00 5817.216 991.2888
+ 227395.2 4024866. -44494.95 5817.216 493.4351
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_vb_m00_002.gmt
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_vb_m00_002.gmt (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_vb_m00_002.gmt 2009-03-28 00:36:31 UTC (rev 14505)
@@ -0,0 +1,2 @@
+ 280726.1 4053913. -45000.00 5817.216 1008.704
+ 280726.1 4053913. -44494.95 5817.216 522.5519
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_vb_m16_001.gmt
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_vb_m16_001.gmt (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_vb_m16_001.gmt 2009-03-28 00:36:31 UTC (rev 14505)
@@ -0,0 +1,2 @@
+ 227395.2 4024866. -45000.00 5816.855 991.2888
+ 227395.2 4024866. -44494.95 5816.855 493.4351
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_vb_m16_002.gmt
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_vb_m16_002.gmt (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_vb_m16_002.gmt 2009-03-28 00:36:31 UTC (rev 14505)
@@ -0,0 +1,2 @@
+ 280726.1 4053913. -45000.00 5817.097 1008.704
+ 280726.1 4053913. -44494.95 5817.097 522.5519
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_vs_m00_001.gmt
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_vs_m00_001.gmt (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_vs_m00_001.gmt 2009-03-28 00:36:31 UTC (rev 14505)
@@ -0,0 +1,2 @@
+ 227395.2 4024866. -45000.00 4500.000 991.2888
+ 227395.2 4024866. -44494.95 4500.000 493.4351
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_vs_m00_002.gmt
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_vs_m00_002.gmt (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_vs_m00_002.gmt 2009-03-28 00:36:31 UTC (rev 14505)
@@ -0,0 +1,2 @@
+ 280726.1 4053913. -45000.00 4500.000 1008.704
+ 280726.1 4053913. -44494.95 4500.000 522.5519
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_vs_m16_001.gmt
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_vs_m16_001.gmt (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_vs_m16_001.gmt 2009-03-28 00:36:31 UTC (rev 14505)
@@ -0,0 +1,2 @@
+ 227395.2 4024866. -45000.00 4502.412 991.2888
+ 227395.2 4024866. -44494.95 4502.412 493.4351
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_vs_m16_002.gmt
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_vs_m16_002.gmt (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_slice/vert_xc_vs_m16_002.gmt 2009-03-28 00:36:31 UTC (rev 14505)
@@ -0,0 +1,2 @@
+ 280726.1 4053913. -45000.00 4501.273 1008.704
+ 280726.1 4053913. -44494.95 4501.273 522.5519
More information about the CIG-COMMITS
mailing list