[cig-commits] r13411 - in seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS: . harness harness/bin harness/chino
leif at geodynamics.org
leif at geodynamics.org
Tue Nov 25 17:37:52 PST 2008
Author: leif
Date: 2008-11-25 17:37:52 -0800 (Tue, 25 Nov 2008)
New Revision: 13411
Added:
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/MANIFEST.in
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/bin/
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/bin/chino
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/__init__.py
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/constants.h
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/exit_mpi.f90
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/harness.py
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/lgndr.f90
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/modeltest.f90
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/moho_stretching.f90
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/precision.h
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/prem_common.f90
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/reduce.f90
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/rthetaphi_xyz.f90
seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/setup.py
Log:
Wrote preliminary test harness.
Added: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/MANIFEST.in
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/MANIFEST.in (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/MANIFEST.in 2008-11-26 01:37:52 UTC (rev 13411)
@@ -0,0 +1,3 @@
+include MANIFEST.in
+include chino/*.f90
+include chino/*.h
Added: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/bin/chino
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/bin/chino (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/bin/chino 2008-11-26 01:37:52 UTC (rev 13411)
@@ -0,0 +1,6 @@
+#!/usr/bin/env python
+
+from chino.harness import Harness
+
+h = Harness("g95", "s20rts/")
+h.run()
Property changes on: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/bin/chino
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/__init__.py
===================================================================
Added: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/constants.h
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/constants.h (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/constants.h 2008-11-26 01:37:52 UTC (rev 13411)
@@ -0,0 +1,433 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 4 . 0
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory, California Institute of Technology, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) California Institute of Technology and University of Pau / CNRS / INRIA
+! February 2008
+!
+! 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.
+!
+!=====================================================================
+
+! constants.h. Generated from constants.h.in by configure.
+
+!
+!--- user can modify parameters below
+!
+
+!
+! 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, SIZE_DOUBLE = 8
+
+! usually the size of integer and logical variables is the same as regular single-precision real variable
+ integer, parameter :: SIZE_INTEGER = SIZE_REAL
+ integer, parameter :: SIZE_LOGICAL = SIZE_REAL
+
+! 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
+
+! if files on a local path on each node are also seen as global with same path
+! set to .true. typically on a shared-memory machine with a common file system
+! set to .false. typically on a cluster of nodes with local disks
+! if running on a cluster of nodes with local disks, also customize global path
+! to local files in create_serial_name_database.f90 ("20 format ...")
+! Flag is used only when one checks the mesh with the serial codes
+! ("xcheck_buffers_1D" etc.), ignore it if you do not plan to use them
+ logical, parameter :: LOCAL_PATH_IS_ALSO_GLOBAL = .false.
+
+! input, output and main MPI I/O files
+ integer, parameter :: ISTANDARD_OUTPUT = 6
+ integer, parameter :: IIN = 40,IOUT = 41,IOUT_SAC = 903
+! local file unit for output of buffers
+ integer, parameter :: IOUT_BUFFERS = 35
+! uncomment this to write messages to a text file
+ integer, parameter :: IMAIN = 42
+! uncomment this to write messages to the screen (slows down the code)
+! integer, parameter :: IMAIN = ISTANDARD_OUTPUT
+! I/O unit for source and receiver vtk file
+ integer, parameter :: IOVTK = 98
+
+
+! R_EARTH is the radius of the bottom of the oceans (radius of Earth in m)
+ double precision, parameter :: R_EARTH = 6371000.d0
+! uncomment line below for PREM with oceans
+! double precision, parameter :: R_EARTH = 6368000.d0
+
+! average density in the full Earth to normalize equation
+ double precision, parameter :: RHOAV = 5514.3d0
+
+! for topography/bathymetry model
+
+!!--- ETOPO5 5-minute model, smoothed Harvard version
+!! size of topography and bathymetry file
+! integer, parameter :: NX_BATHY = 4320,NY_BATHY = 2160
+!! resolution of topography file in minutes
+! integer, parameter :: RESOLUTION_TOPO_FILE = 5
+!! pathname of the topography file
+! character (len=*), parameter :: PATHNAME_TOPO_FILE = 'DATA/topo_bathy/topo_bathy_etopo5_smoothed_Harvard.dat'
+
+!--- ETOPO4 4-minute model created by subsampling and smoothing etopo-2
+! size of topography and bathymetry file
+ integer, parameter :: NX_BATHY = 5400,NY_BATHY = 2700
+! resolution of topography file in minutes
+ integer, parameter :: RESOLUTION_TOPO_FILE = 4
+! pathname of the topography file
+ character (len=*), parameter :: PATHNAME_TOPO_FILE = 'DATA/topo_bathy/topo_bathy_etopo4_smoothed_window_7.dat'
+
+!!--- ETOPO2 2-minute model, not implemented yet
+!! size of topography and bathymetry file
+! integer, parameter :: NX_BATHY = 10800,NY_BATHY = 5400
+!! resolution of topography file in minutes
+! integer, parameter :: RESOLUTION_TOPO_FILE = 2
+!! pathname of the topography file
+! character (len=*), parameter :: PATHNAME_TOPO_FILE = 'DATA/topo_bathy/topo_bathy_etopo2_smoothed_window7.dat'
+
+! maximum depth of the oceans in trenches and height of topo in mountains
+! to avoid taking into account spurious oscillations in global model ETOPO
+ logical, parameter :: USE_MAXIMUM_HEIGHT_TOPO = .false.
+ integer, parameter :: MAXIMUM_HEIGHT_TOPO = +20000
+ logical, parameter :: USE_MAXIMUM_DEPTH_OCEANS = .false.
+ integer, parameter :: MAXIMUM_DEPTH_OCEANS = -20000
+
+! minimum thickness in meters to include the effect of the oceans and topo
+ double precision, parameter :: MINIMUM_THICKNESS_3D_OCEANS = 100.d0
+
+! 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
+
+! flag to exclude elements that are too far from target in source detection
+ logical, parameter :: USE_DISTANCE_CRITERION = .true.
+
+! flag to display detailed information about location of stations
+ logical, parameter :: DISPLAY_DETAILS_STATIONS = .false.
+
+! maximum length of station and network name for receivers
+ integer, parameter :: MAX_LENGTH_STATION_NAME = 32
+ integer, parameter :: MAX_LENGTH_NETWORK_NAME = 8
+
+! we mimic a triangle of half duration equal to half_duration_triangle
+! using a Gaussian having a very close shape, as explained in Figure 4.2
+! of the manual. This source decay rate to mimic an equivalent triangle
+! was found by trial and error
+ double precision, parameter :: SOURCE_DECAY_MIMIC_TRIANGLE = 1.628d0
+
+! maximum number of sources to locate simultaneously
+ integer, parameter :: NSOURCES_SUBSET_MAX = 1000
+
+! distance threshold (in km) above which we consider that a receiver
+! is located outside the mesh and therefore excluded from the station list
+ double precision, parameter :: THRESHOLD_EXCLUDE_STATION = 50.d0
+
+! the first doubling is implemented right below the Moho
+! it seems optimal to implement the three other doublings at these depths
+! in the mantle
+ double precision, parameter :: DEPTH_SECOND_DOUBLING_OPTIMAL = 1650000.d0
+! in the outer core
+ double precision, parameter :: DEPTH_THIRD_DOUBLING_OPTIMAL = 3860000.d0
+! in the outer core
+ double precision, parameter :: DEPTH_FOURTH_DOUBLING_OPTIMAL = 5000000.d0
+
+! Boundary Mesh -- save Moho, 400, 670 km discontinuity topology files (in
+! the mesher) and use them for the computation of boundary kernel (in the solver)
+ logical, parameter :: SAVE_BOUNDARY_MESH = .false.
+
+! this parameter must be set to .true. to compute anisotropic kernels
+! in crust and mantle (related to the 21 Cij in geographical coordinates)
+! default is .false. to compute isotropic kernels (related to alpha and beta)
+ logical, parameter :: ANISOTROPIC_KL = .false.
+
+! print date and time estimate of end of run in another country,
+! in addition to local time.
+! For instance: the code runs at Caltech in California but the person
+! running the code is connected remotely from France, which has 9 hours more.
+! The time difference with that remote location can be positive or negative
+ logical, parameter :: ADD_TIME_ESTIMATE_ELSEWHERE = .false.
+ integer, parameter :: HOURS_TIME_DIFFERENCE = +9
+ integer, parameter :: MINUTES_TIME_DIFFERENCE = +0
+
+!
+!--- debugging flags
+!
+
+! flags to actually assemble with MPI or not
+! and to actually match fluid and solid regions of the Earth or not
+! should always be set to true except when debugging code
+ logical, parameter :: ACTUALLY_ASSEMBLE_MPI_SLICES = .true.
+ logical, parameter :: ACTUALLY_ASSEMBLE_MPI_CHUNKS = .true.
+ logical, parameter :: ACTUALLY_COUPLE_FLUID_CMB = .true.
+ logical, parameter :: ACTUALLY_COUPLE_FLUID_ICB = .true.
+
+!! DK DK UGLY added this in case we are running on MareNostrum in Barcelona
+!! DK DK UGLY because we then need some calls to the system to use GPFS
+ logical, parameter :: RUN_ON_MARENOSTRUM_BARCELONA = .false.
+
+!------------------------------------------------------
+!----------- 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
+ double precision, parameter :: PI_OVER_FOUR = PI / 4.d0
+
+! to convert angles from degrees to radians
+ double precision, parameter :: DEGREES_TO_RADIANS = PI / 180.d0
+
+! 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 with 27 nodes
+ integer, parameter :: NGNOD = 27, NGNOD2D = 9
+
+! gravitational constant
+ double precision, parameter :: GRAV = 6.6723d-11
+
+! 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, &
+ TWO_THIRDS = 2._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
+
+! normalized radius of free surface
+ double precision, parameter :: R_UNIT_SPHERE = ONE
+
+! same radius in km
+ double precision, parameter :: R_EARTH_KM = R_EARTH / 1000.d0
+
+! fixed thickness of 3 km for PREM oceans
+ double precision, parameter :: THICKNESS_OCEANS_PREM = 3000.d0 / R_EARTH
+
+! shortest radius at which crust is implemented (80 km depth)
+! to be constistent with the D80 discontinuity, we impose the crust only above it
+ double precision, parameter :: R_DEEPEST_CRUST = (R_EARTH - 80000.d0) / R_EARTH
+
+! maximum number of chunks (full sphere)
+ integer, parameter :: NCHUNKS_MAX = 6
+
+! define block type based upon chunk number (between 1 and 6)
+! do not change this numbering, chunk AB must be number 1 for central cube
+ integer, parameter :: CHUNK_AB = 1
+ integer, parameter :: CHUNK_AC = 2
+ integer, parameter :: CHUNK_BC = 3
+ integer, parameter :: CHUNK_AC_ANTIPODE = 4
+ integer, parameter :: CHUNK_BC_ANTIPODE = 5
+ integer, parameter :: CHUNK_AB_ANTIPODE = 6
+
+! maximum number of regions in the mesh
+ integer, parameter :: MAX_NUM_REGIONS = 3
+
+! define flag for regions of the global Earth mesh
+ integer, parameter :: IREGION_CRUST_MANTLE = 1
+ integer, parameter :: IREGION_OUTER_CORE = 2
+ integer, parameter :: IREGION_INNER_CORE = 3
+
+! define flag for elements
+ integer, parameter :: IFLAG_CRUST = 1
+
+ integer, parameter :: IFLAG_80_MOHO = 2
+ integer, parameter :: IFLAG_220_80 = 3
+ integer, parameter :: IFLAG_670_220 = 4
+ integer, parameter :: IFLAG_MANTLE_NORMAL = 5
+
+ integer, parameter :: IFLAG_OUTER_CORE_NORMAL = 6
+
+ integer, parameter :: IFLAG_INNER_CORE_NORMAL = 7
+ integer, parameter :: IFLAG_MIDDLE_CENTRAL_CUBE = 8
+ integer, parameter :: IFLAG_BOTTOM_CENTRAL_CUBE = 9
+ integer, parameter :: IFLAG_TOP_CENTRAL_CUBE = 10
+ integer, parameter :: IFLAG_IN_FICTITIOUS_CUBE = 11
+
+ integer, parameter :: NSPEC2D_XI_SUPERBRICK = 8
+ integer, parameter :: NSPEC2D_ETA_SUPERBRICK = 8
+ integer, parameter :: NSPEC2D_XI_SUPERBRICK_1L = 6
+ integer, parameter :: NSPEC2D_ETA_SUPERBRICK_1L = 6
+
+! dummy flag used for mesh display purposes only
+ integer, parameter :: IFLAG_DUMMY = 100
+
+! max number of layers that are used in the radial direction to build the full mesh
+ integer, parameter :: MAX_NUMBER_OF_MESH_LAYERS = 15
+
+! define number of spectral elements and points in basic symmetric mesh doubling superbrick
+ integer, parameter :: NSPEC_DOUBLING_SUPERBRICK = 32
+ integer, parameter :: NGLOB_DOUBLING_SUPERBRICK = 67
+ integer, parameter :: NSPEC_SUPERBRICK_1L = 28
+ integer, parameter :: NGLOB_SUPERBRICK_1L = 58
+ integer, parameter :: NGNOD_EIGHT_CORNERS = 8
+
+! define flag for regions of the global Earth for attenuation
+ integer, parameter :: NUM_REGIONS_ATTENUATION = 5
+
+ integer, parameter :: IREGION_ATTENUATION_INNER_CORE = 1
+ integer, parameter :: IREGION_ATTENUATION_CMB_670 = 2
+ integer, parameter :: IREGION_ATTENUATION_670_220 = 3
+ integer, parameter :: IREGION_ATTENUATION_220_80 = 4
+ integer, parameter :: IREGION_ATTENUATION_80_SURFACE = 5
+ integer, parameter :: IREGION_ATTENUATION_UNDEFINED = 6
+
+! number of standard linear solids for attenuation
+ integer, parameter :: N_SLS = 3
+
+! computation of standard linear solids in meshfem3D
+! ATTENUATION_COMP_RESOLUTION: Number of Digits after decimal
+! ATTENUATION_COMP_MAXIMUM: Maximum Q Value
+ integer, parameter :: ATTENUATION_COMP_RESOLUTION = 1
+ integer, parameter :: ATTENUATION_COMP_MAXIMUM = 5000
+
+! for lookup table for attenuation every 100 m in radial direction of Earth model
+ integer, parameter :: NRAD_ATTENUATION = 70000
+ double precision, parameter :: TABLE_ATTENUATION = R_EARTH_KM * 10.0d0
+
+! for determination of the attenuation period range
+! if this is set to .true. then the hardcoded values will be used
+! otherwise they are computed automatically from the Number of elements
+! This *may* be a useful parameter for Benchmarking against older versions
+ logical, parameter :: ATTENUATION_RANGE_PREDEFINED = .false.
+
+! 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
+
+! flags to select the right corner in each slice
+ integer, parameter :: ILOWERLOWER = 1
+ integer, parameter :: ILOWERUPPER = 2
+ integer, parameter :: IUPPERLOWER = 3
+ integer, parameter :: IUPPERUPPER = 4
+
+! number of points in each AVS or OpenDX quadrangular cell for movies
+ integer, parameter :: NGNOD2D_AVS_DX = 4
+
+! number of faces a given slice can share with other slices
+! this is at most 2, except when there is only once slice per chunk
+! in which case it is 4
+ integer, parameter :: NUMFACES_SHARED = 4
+
+! number of corners a given slice can share with other slices
+! this is at most 1, except when there is only once slice per chunk
+! in which case it is 4
+ integer, parameter :: NUMCORNERS_SHARED = 4
+
+! number of slaves per corner
+ integer, parameter :: NUMSLAVES = 2
+
+! number of layers in PREM
+ integer, parameter :: NR = 640
+
+! smallest real number on many machines = 1.1754944E-38
+! largest real number on many machines = 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 that 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
+
+! small tolerance for conversion from x y z to r theta phi
+ double precision, parameter :: SMALL_VAL_ANGLE = 1.d-10
+
+! geometry tolerance parameter to calculate number of independent grid points
+! sensitive to actual size of model, assumes reference sphere of radius 1
+! this is an absolute value for normalized coordinates in the Earth
+ double precision, parameter :: SMALLVALTOL = 1.d-10
+
+! 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 non linear system for xi and eta
+ integer, parameter :: NUM_ITER = 4
+
+! number of hours per day for rotation rate of the Earth
+ double precision, parameter :: HOURS_PER_DAY = 24.d0
+
+! for lookup table for gravity every 100 m in radial direction of Earth model
+ integer, parameter :: NRAD_GRAVITY = 70000
+
+! to inflate the central cube (set to 0.d0 for a non-inflated cube)
+ double precision, parameter :: CENTRAL_CUBE_INFLATE_FACTOR = 0.41d0
+
+! for the stretching of crustal elements in the case of 3D models
+ double precision, parameter :: MAX_RATIO_CRUST_STRETCHING = 0.6d0
+
+! to suppress the crustal layers (replaced by an extension of the mantle: R_EARTH is not modified, but no more crustal doubling)
+ logical, parameter :: SUPPRESS_CRUSTAL_MESH = .false.
+
+! to add a fourth doubling at the bottom of the outer core
+ logical, parameter :: ADD_4TH_DOUBLING = .false.
+
+!!! parameters for the doubling brick cut
+
+! this for cutting superbrick : 3 possibilities, 4 cases max / possibility
+! 3 possibilities : cut in xi & eta || cut in xi || cut in eta
+! case 1 : ximin & etamin || ximin || etamin
+! case 2 : ximin & etamax || ximax || etamax
+! case 3 : ximax & etamin
+! case 4 : ximax & etamax
+ integer, parameter :: NB_CUT_CASE = 4
+
+! corner 1 : ximin & etamin
+! corner 2 : ximax & etamin
+! corner 3 : ximax & etamax
+! corner 4 : ximin & etamax
+ integer, parameter :: NB_SQUARE_CORNERS = 4
+
+! 2 possibilities : xi || eta
+! face 1 : ximin || etamin
+! face 2 : ximax || etamax
+ integer, parameter :: NB_SQUARE_EDGES_ONEDIR = 2
+
+! this for the geometry of the basic doubling brick
+ integer, parameter :: NSPEC_DOUBLING_BASICBRICK = 8
+ integer, parameter :: NGLOB_DOUBLING_BASICBRICK = 27
+
Copied: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/exit_mpi.f90 (from rev 13385, seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/exit_mpi.f90)
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/exit_mpi.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/exit_mpi.f90 2008-11-26 01:37:52 UTC (rev 13411)
@@ -0,0 +1,66 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 4 . 0
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory, California Institute of Technology, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) California Institute of Technology and University of Pau / CNRS / INRIA
+! February 2008
+!
+! 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.
+!
+!=====================================================================
+
+! end the simulation and exit MPI
+
+! version with rank number printed in the error message
+ subroutine exit_MPI(myrank,error_msg)
+
+ implicit none
+
+ include "constants.h"
+
+ integer myrank
+ character(len=*) error_msg
+
+ write(*,*) error_msg(1:len(error_msg))
+ write(*,*) 'Error detected, aborting MPI... proc ',myrank
+
+ 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
+
+ include "constants.h"
+
+ character(len=*) error_msg
+
+ write(*,*) error_msg(1:len(error_msg))
+ write(*,*) 'Error detected, aborting MPI...'
+
+ stop 'error, program ended in exit_MPI'
+
+ end subroutine exit_MPI_without_rank
+
Added: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/harness.py
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/harness.py (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/harness.py 2008-11-26 01:37:52 UTC (rev 13411)
@@ -0,0 +1,199 @@
+
+harnessSourceFiles = [
+ "exit_mpi.f90",
+ "lgndr.f90",
+ "modeltest.f90",
+ "moho_stretching.f90",
+ "prem_common.f90",
+ "reduce.f90",
+ "rthetaphi_xyz.f90",
+]
+
+harnessHeaderFiles = [
+ "constants.h",
+ "precision.h",
+]
+
+
+class Harness(object):
+
+ def __init__(self, fc, modelDir):
+ from os.path import abspath
+ self.fc = fc
+ self.modelDir = abspath(modelDir)
+
+ def tarUp(self):
+ import tarfile
+ from os import chdir
+ from os.path import split
+
+ tar = tarfile.open(self.model, "w:gz")
+ parent, child = split(self.modelDir)
+ chdir(parent)
+ tar.add(child)
+ tar.close()
+
+ def extractModel(self):
+ import tarfile
+ import os
+ from os.path import dirname
+
+ tgz = tarfile.open(self.model, 'r:gz')
+ path = "model"
+
+ directories = []
+ fortranSourceFiles = []
+
+ for tarinfo in tgz:
+ if tarinfo.isdir():
+ # Extract directory with a safe mode, so that
+ # all files below can be extracted as well.
+ try:
+ os.makedirs(os.path.join(path, tarinfo.name), 0777)
+ except EnvironmentError:
+ pass
+ directories.append(tarinfo)
+ elif tarinfo.name.endswith(".f90"):
+ pathname = os.path.join(path, tarinfo.name)
+ fortranSourceFiles.append(pathname)
+ thisDir = dirname(tarinfo.name) # see bcast_model.c
+ s = tgz.extractfile(tarinfo)
+ f = open(pathname, "w")
+ # Preprocess.
+ for line in s.readlines():
+ line = line.replace('@THIS_DIR@', thisDir)
+ f.write(line)
+ else:
+ tgz.extract(tarinfo, path)
+
+ # Reverse sort directories.
+ directories.sort(lambda a, b: cmp(a.name, b.name))
+ directories.reverse()
+
+ # Set correct owner, mtime and filemode on directories.
+ for tarinfo in directories:
+ path = os.path.join(path, tarinfo.name)
+ try:
+ tgz.chown(tarinfo, path)
+ tgz.utime(tarinfo, path)
+ tgz.chmod(tarinfo, path)
+ except tarfile.ExtractError, e:
+ pass
+
+ self.modelSourceFiles = fortranSourceFiles
+
+ return
+
+ def copyHarnessSourceFiles(self):
+ from os.path import dirname, join
+ from shutil import copyfile
+ from itertools import chain
+ pkgDir = dirname(__file__)
+ for sourceFile in chain(harnessSourceFiles, harnessHeaderFiles):
+ copyfile(join(pkgDir, sourceFile), sourceFile)
+ return
+
+ def generateMakefile(self):
+ from os.path import basename, splitext
+ from itertools import chain
+
+ modelSourceFiles = self.modelSourceFiles
+ s = open("Makefile", "w")
+ print >>s
+ print >>s, "FC = %s" % self.fc
+ print >>s, "FCFLAGS_f90 ="
+ print >>s
+ print >>s, "model_OBJECTS = \\"
+ for sourceFile in modelSourceFiles:
+ base = splitext(basename(sourceFile))[0]
+ print >>s, "\t%s.o \\" % base
+ print >>s, "\t$(empty)"
+ print >>s
+ print >>s, "harness_OBJECTS = \\"
+ for sourceFile in harnessSourceFiles:
+ base = splitext(basename(sourceFile))[0]
+ print >>s, "\t%s.o \\" % base
+ print >>s, "\t$(empty)"
+ print >>s
+ print >>s, "modeltest: $(model_OBJECTS) $(harness_OBJECTS)"
+ print >>s, "\t$(FC) -o modeltest $(model_OBJECTS) $(harness_OBJECTS)"
+ print >>s
+ for sourceFile in chain(modelSourceFiles, harnessSourceFiles):
+ base = splitext(basename(sourceFile))[0]
+ print >>s, "%s.o: constants.h %s" % (base, sourceFile)
+ print >>s, "\t$(FC) -c -o %s.o $(FCFLAGS_f90) %s" % (base, sourceFile)
+ print >>s
+ return
+
+ def spawn(self, *argv):
+ import os, sys
+ status = os.spawnvp(os.P_WAIT, argv[0], argv)
+ if status != 0:
+ sys.exit("%s: exit %d" % (argv[0], status))
+ return
+
+ def ospawn(self, *argv):
+ import os, sys
+ from popen2 import Popen4
+
+ child = Popen4(argv)
+ child.tochild.close()
+ output = child.fromchild.readlines()
+ status = child.wait()
+
+ exitStatus = None
+ if (os.WIFSIGNALED(status)):
+ statusStr = "signal %d" % os.WTERMSIG(status)
+ elif (os.WIFEXITED(status)):
+ exitStatus = os.WEXITSTATUS(status)
+ statusStr = "exit %d" % exitStatus
+ else:
+ statusStr = "status %d" % status
+
+ return output, exitStatus
+
+ def build(self):
+ self.spawn("make")
+
+ def execute(self):
+ import sys
+
+ output, exitStatus = self.ospawn("../modeltest")
+ if exitStatus != 0:
+ for line in output:
+ print line
+ sys.exit("chino: model test failed")
+
+ for line in output:
+ if line.find("Luchooz0") != -1:
+ print line
+ return
+
+ sys.exit("chino: model test failed")
+
+ def run(self):
+ from tempfile import mkdtemp
+ from os import chdir, getcwd, mkdir
+ from os.path import abspath, join
+ from shutil import rmtree
+
+ oldWd = getcwd()
+
+ tempDir = mkdtemp()
+ self.model = abspath(join(tempDir, "model.tgz"))
+
+ self.tarUp()
+
+ chdir(tempDir)
+ self.extractModel()
+ self.copyHarnessSourceFiles()
+ self.generateMakefile()
+ self.build()
+
+ chdir("model")
+ self.execute()
+
+ chdir(oldWd)
+ rmtree(tempDir)
+
+ return
Copied: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/lgndr.f90 (from rev 13385, seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/lgndr.f90)
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/lgndr.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/lgndr.f90 2008-11-26 01:37:52 UTC (rev 13411)
@@ -0,0 +1,152 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 4 . 0
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory, California Institute of Technology, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) California Institute of Technology and University of Pau / CNRS / INRIA
+! February 2008
+!
+! 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 lgndr(l,c,s,x,dx)
+
+! computes Legendre function x(l,m,theta)
+! theta=colatitude,c=cos(theta),s=sin(theta),l=angular order,
+! sin(theta) restricted so that sin(theta) > 1.e-7
+! x(1) contains m=0, x(2) contains m=1, x(k+1) contains m=k
+! m=azimuthal (longitudinal) order 0 <= m <= l
+! dx=dx/dtheta
+!
+! subroutine originally came from Physics Dept. Princeton through
+! Peter Davis, modified by Jeffrey Park
+
+ implicit none
+
+! argument variables
+ integer l
+ double precision x(2*l+1),dx(2*l+1)
+ double precision c,s
+
+! local variables
+ integer i,lp1,lpsafe,lsave
+ integer m,maxsin,mmm,mp1
+
+ double precision sqroot2over2,c1,c2,cot
+ double precision ct,d,f1,f2
+ double precision f3,fac,g1,g2
+ double precision g3,rfpi,sqroot3,sos
+ double precision ss,stom,t,tol
+ double precision v,y
+
+ tol = 1.d-05
+ rfpi = 0.282094791773880d0
+ sqroot3 = 1.73205080756890d0
+ sqroot2over2 = 0.707106781186550d0
+
+ if(s >= 1.0d0-tol) s=1.0d0-tol
+ lsave=l
+ if(l<0) l=-1-l
+ if(l>0) goto 1
+ x(1)=rfpi
+ dx(1)=0.0d0
+ l=lsave
+ return
+ 1 if(l /= 1) goto 2
+ c1=sqroot3*rfpi
+ c2=sqroot2over2*c1
+ x(1)=c1*c
+ x(2)=-c2*s
+ dx(1)=-c1*s
+ dx(2)=-c2*c
+ l=lsave
+ return
+ 2 sos=s
+ if(s<tol) s=tol
+ cot=c/s
+ ct=2.0d0*c
+ ss=s*s
+ lp1=l+1
+ g3=0.0d0
+ g2=1.0d0
+ f3=0.0d0
+
+! evaluate m=l value, sans (sin(theta))**l
+ do i=1,l
+ g2=g2*(1.0d0-1.0d0/(2.0d0*i))
+ enddo
+ g2=rfpi*dsqrt((2*l+1)*g2)
+ f2=l*cot*g2
+ x(lp1)=g2
+ dx(lp1)=f2
+ v=1.0d0
+ y=2.0d0*l
+ d=dsqrt(v*y)
+ t=0.0d0
+ mp1=l
+ m=l-1
+
+! these recursions are similar to ordinary m-recursions, but since we
+! have taken the s**m factor out of the xlm's, the recursion has the powers
+! of sin(theta) instead
+ 3 g1=-(ct*mp1*g2+ss*t*g3)/d
+ f1=(mp1*(2.0d0*s*g2-ct*f2)-t*ss*(f3+cot*g3))/d-cot*g1
+ x(mp1)=g1
+ dx(mp1)=f1
+ if(m == 0) goto 4
+ mp1=m
+ m=m-1
+ v=v+1.0d0
+ y=y-1.0d0
+ t=d
+ d=dsqrt(v*y)
+ g3=g2
+ g2=g1
+ f3=f2
+ f2=f1
+ goto 3
+! explicit conversion to integer added
+ 4 maxsin=int(-72.0d0/log10(s))
+
+! maxsin is the max exponent of sin(theta) without underflow
+ lpsafe=min0(lp1,maxsin)
+ stom=1.0d0
+ fac=sign(1.0d0,dble((l/2)*2-l) + 0.50d0)
+
+! multiply xlm by sin**m
+ do m=1,lpsafe
+ x(m)=fac*x(m)*stom
+ dx(m)=fac*dx(m)*stom
+ stom=stom*s
+ enddo
+
+! set any remaining xlm to zero
+ if(maxsin <= l) then
+ mmm=maxsin+1
+ do m=mmm,lp1
+ x(m)=0.0d0
+ dx(m)=0.0d0
+ enddo
+ endif
+
+ s=sos
+ l=lsave
+
+ end subroutine lgndr
+
Added: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/modeltest.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/modeltest.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/modeltest.f90 2008-11-26 01:37:52 UTC (rev 13411)
@@ -0,0 +1,50 @@
+
+program ModelTest
+
+ logical ANISOTROPIC_3D_MANTLE, &
+ ANISOTROPIC_INNER_CORE, &
+ ATTENUATION_3D, &
+ CASE_3D, &
+ CRUSTAL, &
+ HONOR_1D_SPHERICAL_MOHO, &
+ ISOTROPIC_3D_MANTLE, &
+ ONE_CRUST, &
+ TRANSVERSE_ISOTROPY
+
+ double precision ROCEAN,RMIDDLE_CRUST,RMOHO, &
+ R80,R120,R220,R400,R600,R670,R771, &
+ RTOPDDOUBLEPRIME,RCMB,RICB,RHO_TOP_OC,RHO_BOTTOM_OC
+
+ ANISOTROPIC_3D_MANTLE = .false.
+ ANISOTROPIC_INNER_CORE = .false.
+ ATTENUATION_3D = .false.
+ CASE_3D = .false.
+ CRUSTAL = .false.
+ HONOR_1D_SPHERICAL_MOHO = .false.
+ ISOTROPIC_3D_MANTLE = .false.
+ ONE_CRUST = .false.
+ TRANSVERSE_ISOTROPY = .false.
+
+ call get_model_properties(HONOR_1D_SPHERICAL_MOHO,ONE_CRUST, &
+ TRANSVERSE_ISOTROPY, &
+ ISOTROPIC_3D_MANTLE,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE, &
+ CRUSTAL,CASE_3D, &
+ ATTENUATION_3D)
+
+ call get_reference_1d_model_radii(ROCEAN,RMIDDLE_CRUST,RMOHO, &
+ R80,R120,R220,R400,R600,R670,R771, &
+ RTOPDDOUBLEPRIME,RCMB,RICB,RHO_TOP_OC,RHO_BOTTOM_OC)
+
+ call read_3d_mantle_model()
+
+ if (CRUSTAL) then
+ call read_crust()
+ endif
+
+ call define_reference_1d_model(CRUSTAL)
+
+ ! XXX: Write code to test model here.
+
+ write (*,*) 'Luchooz0 model OK'
+
+end program ModelTest
Copied: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/moho_stretching.f90 (from rev 13385, seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/moho_stretching.f90)
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/moho_stretching.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/moho_stretching.f90 2008-11-26 01:37:52 UTC (rev 13411)
@@ -0,0 +1,301 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 4 . 0
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory, California Institute of Technology, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) California Institute of Technology and University of Pau / CNRS / INRIA
+! February 2008
+!
+! 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 moho_stretching(myrank,xelm,yelm,zelm,RMOHO,R220)
+
+ implicit none
+
+ include "constants.h"
+
+! ocean-continent function maximum spherical harmonic degree
+ integer, parameter :: NL_OCEAN_CONTINENT = 12
+
+! spherical harmonic coefficients of the ocean-continent function (km)
+ double precision A_lm(0:NL_OCEAN_CONTINENT,0:NL_OCEAN_CONTINENT),B_lm(0:NL_OCEAN_CONTINENT,0:NL_OCEAN_CONTINENT)
+
+ common /smooth_moho/ A_lm,B_lm
+
+ integer myrank
+
+ double precision xelm(NGNOD)
+ double precision yelm(NGNOD)
+ double precision zelm(NGNOD)
+
+ double precision RMOHO,R220
+
+ integer ia
+
+ integer l,m
+ double precision r,theta,phi
+ double precision sint,cost,x(2*NL_OCEAN_CONTINENT+1),dx(2*NL_OCEAN_CONTINENT+1)
+ double precision elevation
+ double precision gamma
+
+! we loop on all the points of the element
+ do ia = 1,NGNOD
+
+! convert to r theta phi
+ call xyz_2_rthetaphi_dble(xelm(ia),yelm(ia),zelm(ia),r,theta,phi)
+ call reduce(theta,phi)
+
+ elevation = 0.0d0
+ do l = 0,NL_OCEAN_CONTINENT
+ sint = dsin(theta)
+ cost = dcos(theta)
+ call lgndr(l,cost,sint,x,dx)
+ m = 0
+ elevation = elevation + A_lm(l,m)*x(m+1)
+ do m = 1,l
+ elevation = elevation + (A_lm(l,m)*dcos(dble(m)*phi)+B_lm(l,m)*dsin(dble(m)*phi))*x(m+1)
+ enddo
+ enddo
+ elevation = -0.25d0*elevation/R_EARTH_KM
+
+ gamma = 0.0d0
+ if(r >= RMOHO/R_EARTH) then
+! stretching above the Moho
+ gamma = (1.0d0 - r) / (1.0d0 - RMOHO/R_EARTH)
+ elseif(r>= R220/R_EARTH .and. r< RMOHO/R_EARTH) then
+! stretching between R220 and RMOHO
+ gamma = (r - R220/R_EARTH) / (RMOHO/R_EARTH - R220/R_EARTH)
+ endif
+ if(gamma < -0.0001 .or. gamma > 1.0001) call exit_MPI(myrank,'incorrect value of gamma for Moho topography')
+
+ xelm(ia) = xelm(ia)*(ONE + gamma * elevation / r)
+ yelm(ia) = yelm(ia)*(ONE + gamma * elevation / r)
+ zelm(ia) = zelm(ia)*(ONE + gamma * elevation / r)
+
+ enddo
+
+ end subroutine moho_stretching
+
+ subroutine read_smooth_moho
+
+ implicit none
+
+! ocean-continent function maximum spherical harmonic degree
+ integer, parameter :: NL_OCEAN_CONTINENT = 12
+
+! spherical harmonic coefficients of the ocean-continent function (km)
+ double precision A_lm(0:NL_OCEAN_CONTINENT,0:NL_OCEAN_CONTINENT),B_lm(0:NL_OCEAN_CONTINENT,0:NL_OCEAN_CONTINENT)
+
+ common /smooth_moho/ A_lm,B_lm
+
+! integer l,m
+!
+! ocean-continent function (km)
+! open(unit=10,file='DATA/ocean_continent_function/ocean_continent_function.txt',status='old',action='read')
+! do l=0,NL_OCEAN_CONTINENT
+! read(10,*) A_lm(l,0),(A_lm(l,m),B_lm(l,m),m=1,l)
+! enddo
+! close(10)
+
+ A_lm(0,0) = -3.8201999E-04
+ B_lm(0,0) = 0.
+ A_lm(1,0) = 13.88800
+ B_lm(1,0) = 0.
+ A_lm(1,1) = -15.24000
+ B_lm(1,1) = -9.187200
+ A_lm(2,0) = 11.21500
+ B_lm(2,0) = 0.
+ A_lm(2,1) = -6.754500
+ B_lm(2,1) = -8.516700
+ A_lm(2,2) = -8.327800
+ B_lm(2,2) = -5.029200
+ A_lm(3,0) = -3.614500
+ B_lm(3,0) = 0.
+ A_lm(3,1) = 5.394800
+ B_lm(3,1) = -0.9220800
+ A_lm(3,2) = -10.05100
+ B_lm(3,2) = 13.98100
+ A_lm(3,3) = -2.711200
+ B_lm(3,3) = -13.57100
+ A_lm(4,0) = 7.523300
+ B_lm(4,0) = 0.
+ A_lm(4,1) = 5.156100
+ B_lm(4,1) = 2.184400
+ A_lm(4,2) = -10.67300
+ B_lm(4,2) = 2.640600
+ A_lm(4,3) = -7.786300
+ B_lm(4,3) = 0.3674500
+ A_lm(4,4) = -3.076400
+ B_lm(4,4) = 16.83000
+ A_lm(5,0) = -9.681000
+ B_lm(5,0) = 0.
+ A_lm(5,1) = 0.5026800
+ B_lm(5,1) = 2.111300
+ A_lm(5,2) = -2.931000
+ B_lm(5,2) = -4.329000
+ A_lm(5,3) = -1.766800
+ B_lm(5,3) = -3.621200
+ A_lm(5,4) = 16.08200
+ B_lm(5,4) = -4.493900
+ A_lm(5,5) = -0.3705800
+ B_lm(5,5) = -5.574500
+ A_lm(6,0) = 4.407900
+ B_lm(6,0) = 0.
+ A_lm(6,1) = 0.3799000
+ B_lm(6,1) = 1.589400
+ A_lm(6,2) = -1.886400
+ B_lm(6,2) = -0.5686300
+ A_lm(6,3) = -0.9816800
+ B_lm(6,3) = -5.827800
+ A_lm(6,4) = 3.620600
+ B_lm(6,4) = -2.713100
+ A_lm(6,5) = 1.445600
+ B_lm(6,5) = 3.964100
+ A_lm(6,6) = 1.167400
+ B_lm(6,6) = 2.134100
+ A_lm(7,0) = -4.086100
+ B_lm(7,0) = 0.
+ A_lm(7,1) = 0.5462000
+ B_lm(7,1) = -4.488100
+ A_lm(7,2) = 3.116400
+ B_lm(7,2) = 1.793600
+ A_lm(7,3) = 2.594600
+ B_lm(7,3) = -2.129100
+ A_lm(7,4) = -5.445000
+ B_lm(7,4) = 0.5381500
+ A_lm(7,5) = -2.178100
+ B_lm(7,5) = 1.766700
+ A_lm(7,6) = -1.040000
+ B_lm(7,6) = -5.541000
+ A_lm(7,7) = 1.536500
+ B_lm(7,7) = 3.700600
+ A_lm(8,0) = -2.562200
+ B_lm(8,0) = 0.
+ A_lm(8,1) = 0.3736200
+ B_lm(8,1) = 1.488000
+ A_lm(8,2) = 1.347500
+ B_lm(8,2) = 0.5288200
+ A_lm(8,3) = -0.8493700
+ B_lm(8,3) = -1.626500
+ A_lm(8,4) = 0.2423400
+ B_lm(8,4) = 4.202800
+ A_lm(8,5) = 2.052200
+ B_lm(8,5) = 0.6880400
+ A_lm(8,6) = 2.838500
+ B_lm(8,6) = 2.835700
+ A_lm(8,7) = -4.981400
+ B_lm(8,7) = -1.883100
+ A_lm(8,8) = -1.102800
+ B_lm(8,8) = -1.951700
+ A_lm(9,0) = -1.202100
+ B_lm(9,0) = 0.
+ A_lm(9,1) = 1.020300
+ B_lm(9,1) = 1.371000
+ A_lm(9,2) = -0.3430100
+ B_lm(9,2) = 0.8782800
+ A_lm(9,3) = -0.4462500
+ B_lm(9,3) = -0.3046100
+ A_lm(9,4) = 0.7750700
+ B_lm(9,4) = 2.351600
+ A_lm(9,5) = -2.092600
+ B_lm(9,5) = -2.377100
+ A_lm(9,6) = 0.3126900
+ B_lm(9,6) = 4.996000
+ A_lm(9,7) = -2.284000
+ B_lm(9,7) = 1.183700
+ A_lm(9,8) = 1.445900
+ B_lm(9,8) = 1.080000
+ A_lm(9,9) = 1.146700
+ B_lm(9,9) = 1.457800
+ A_lm(10,0) = -2.516900
+ B_lm(10,0) = 0.
+ A_lm(10,1) = -0.9739500
+ B_lm(10,1) = -0.7195500
+ A_lm(10,2) = -2.846000
+ B_lm(10,2) = -1.464700
+ A_lm(10,3) = 2.720100
+ B_lm(10,3) = 0.8241400
+ A_lm(10,4) = -1.247800
+ B_lm(10,4) = 1.220300
+ A_lm(10,5) = -1.638500
+ B_lm(10,5) = -1.099500
+ A_lm(10,6) = 3.043000
+ B_lm(10,6) = -1.976400
+ A_lm(10,7) = -1.007300
+ B_lm(10,7) = -1.604900
+ A_lm(10,8) = 0.6620500
+ B_lm(10,8) = -1.135000
+ A_lm(10,9) = -3.576800
+ B_lm(10,9) = 0.5554900
+ A_lm(10,10) = 2.418700
+ B_lm(10,10) = -1.482200
+ A_lm(11,0) = 0.7158800
+ B_lm(11,0) = 0.
+ A_lm(11,1) = -3.694800
+ B_lm(11,1) = 0.8491400
+ A_lm(11,2) = 9.3208998E-02
+ B_lm(11,2) = -1.276000
+ A_lm(11,3) = 1.575600
+ B_lm(11,3) = 0.1972100
+ A_lm(11,4) = 0.8989600
+ B_lm(11,4) = -1.063000
+ A_lm(11,5) = -0.6301000
+ B_lm(11,5) = -1.329400
+ A_lm(11,6) = 1.389000
+ B_lm(11,6) = 1.184100
+ A_lm(11,7) = 0.5640700
+ B_lm(11,7) = 2.286200
+ A_lm(11,8) = 1.530300
+ B_lm(11,8) = 0.7677500
+ A_lm(11,9) = 0.8495500
+ B_lm(11,9) = 0.7247500
+ A_lm(11,10) = 2.106800
+ B_lm(11,10) = 0.6588000
+ A_lm(11,11) = 0.6067800
+ B_lm(11,11) = 0.1366800
+ A_lm(12,0) = -2.598700
+ B_lm(12,0) = 0.
+ A_lm(12,1) = -1.150500
+ B_lm(12,1) = -0.8425700
+ A_lm(12,2) = -0.1593300
+ B_lm(12,2) = -1.241400
+ A_lm(12,3) = 1.508600
+ B_lm(12,3) = 0.3385500
+ A_lm(12,4) = -1.941200
+ B_lm(12,4) = 1.120000
+ A_lm(12,5) = -0.4630500
+ B_lm(12,5) = -6.4753003E-02
+ A_lm(12,6) = 0.8967000
+ B_lm(12,6) = 4.7417998E-02
+ A_lm(12,7) = 4.5407999E-02
+ B_lm(12,7) = 0.8876400
+ A_lm(12,8) = -2.444400
+ B_lm(12,8) = 1.172500
+ A_lm(12,9) = -2.593400
+ B_lm(12,9) = 0.1703700
+ A_lm(12,10) = 0.5662700
+ B_lm(12,10) = 0.7050800
+ A_lm(12,11) = -0.1930000
+ B_lm(12,11) = -2.008100
+ A_lm(12,12) = -3.187900
+ B_lm(12,12) = -1.672000
+
+ end subroutine read_smooth_moho
+
Added: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/precision.h
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/precision.h (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/precision.h 2008-11-26 01:37:52 UTC (rev 13411)
@@ -0,0 +1,38 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 4 . 0
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory, California Institute of Technology, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) California Institute of Technology and University of Pau / CNRS / INRIA
+! February 2008
+!
+! 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.
+!
+!=====================================================================
+
+! 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
Copied: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/prem_common.f90 (from rev 13385, seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/prem_common.f90)
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/prem_common.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/prem_common.f90 2008-11-26 01:37:52 UTC (rev 13411)
@@ -0,0 +1,363 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 4 . 0
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory, California Institute of Technology, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) California Institute of Technology and University of Pau / CNRS / INRIA
+! February 2008
+!
+! 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 prem_iso(myrank,x,rho,drhodr,vp,vs,Qkappa,Qmu,idoubling,CRUSTAL, &
+ ONE_CRUST,check_doubling_flag)
+
+ implicit none
+
+ include "constants.h"
+
+! given a normalized radius x, gives the non-dimesionalized density rho,
+! speeds vp and vs, and the quality factors Qkappa and Qmu
+
+ logical CRUSTAL,ONE_CRUST,check_doubling_flag
+
+ integer idoubling,myrank
+
+ double precision x,rho,drhodr,vp,vs,Qkappa,Qmu,RICB,RCMB,RTOPDDOUBLEPRIME, &
+ R600,R670,R220,R771,R400,R80,RMOHO,RMIDDLE_CRUST,ROCEAN
+ double precision R120,RHO_TOP_OC,RHO_BOTTOM_OC
+
+ double precision r,scaleval
+
+ call get_reference_1d_model_radii(ROCEAN,RMIDDLE_CRUST,RMOHO, &
+ R80,R120,R220,R400,R600,R670,R771, &
+ RTOPDDOUBLEPRIME,RCMB,RICB,RHO_TOP_OC,RHO_BOTTOM_OC)
+
+! compute real physical radius in meters
+ r = x * R_EARTH
+
+! check flags to make sure we correctly honor the discontinuities
+! we use strict inequalities since r has been slighly changed in mesher
+
+ if(check_doubling_flag) then
+
+!
+!--- inner core
+!
+
+ if(r >= 0.d0 .and. r < RICB) then
+ if(idoubling /= IFLAG_INNER_CORE_NORMAL .and. &
+ idoubling /= IFLAG_MIDDLE_CENTRAL_CUBE .and. &
+ idoubling /= IFLAG_BOTTOM_CENTRAL_CUBE .and. &
+ idoubling /= IFLAG_TOP_CENTRAL_CUBE .and. &
+ idoubling /= IFLAG_IN_FICTITIOUS_CUBE) &
+ call exit_MPI(myrank,'wrong doubling flag for inner core point')
+!
+!--- outer core
+!
+ else if(r > RICB .and. r < RCMB) then
+ if(idoubling /= IFLAG_OUTER_CORE_NORMAL) &
+ call exit_MPI(myrank,'wrong doubling flag for outer core point')
+!
+!--- D" at the base of the mantle
+!
+ else if(r > RCMB .and. r < RTOPDDOUBLEPRIME) then
+ if(idoubling /= IFLAG_MANTLE_NORMAL) &
+ call exit_MPI(myrank,'wrong doubling flag for D" point')
+!
+!--- mantle: from top of D" to d670
+!
+ else if(r > RTOPDDOUBLEPRIME .and. r < R670) then
+ if(idoubling /= IFLAG_MANTLE_NORMAL) &
+ call exit_MPI(myrank,'wrong doubling flag for top D" -> d670 point')
+
+!
+!--- mantle: from d670 to d220
+!
+ else if(r > R670 .and. r < R220) then
+ if(idoubling /= IFLAG_670_220) &
+ call exit_MPI(myrank,'wrong doubling flag for d670 -> d220 point')
+
+!
+!--- mantle and crust: from d220 to MOHO and then to surface
+!
+ else if(r > R220) then
+ if(idoubling /= IFLAG_220_80 .and. idoubling /= IFLAG_80_MOHO .and. idoubling /= IFLAG_CRUST) &
+ call exit_MPI(myrank,'wrong doubling flag for d220 -> Moho -> surface point')
+
+ endif
+
+ endif
+
+!
+!--- inner core
+!
+ if(r >= 0.d0 .and. r <= RICB) then
+ drhodr=-2.0d0*8.8381d0*x
+ rho=13.0885d0-8.8381d0*x*x
+ vp=11.2622d0-6.3640d0*x*x
+ vs=3.6678d0-4.4475d0*x*x
+ Qmu=84.6d0
+ Qkappa=1327.7d0
+!
+!--- outer core
+!
+ else if(r > RICB .and. r <= RCMB) then
+ drhodr=-1.2638d0-2.0d0*3.6426d0*x-3.0d0*5.5281d0*x*x
+ rho=12.5815d0-1.2638d0*x-3.6426d0*x*x-5.5281d0*x*x*x
+ vp=11.0487d0-4.0362d0*x+4.8023d0*x*x-13.5732d0*x*x*x
+ vs=0.0d0
+ Qmu=0.0d0
+ Qkappa=57827.0d0
+!
+!--- D" at the base of the mantle
+!
+ else if(r > RCMB .and. r <= RTOPDDOUBLEPRIME) then
+ drhodr=-6.4761d0+2.0d0*5.5283d0*x-3.0d0*3.0807d0*x*x
+ rho=7.9565d0-6.4761d0*x+5.5283d0*x*x-3.0807d0*x*x*x
+ vp=15.3891d0-5.3181d0*x+5.5242d0*x*x-2.5514d0*x*x*x
+ vs=6.9254d0+1.4672d0*x-2.0834d0*x*x+0.9783d0*x*x*x
+ Qmu=312.0d0
+ Qkappa=57827.0d0
+!
+!--- mantle: from top of D" to d670
+!
+ else if(r > RTOPDDOUBLEPRIME .and. r <= R771) then
+ drhodr=-6.4761d0+2.0d0*5.5283d0*x-3.0d0*3.0807d0*x*x
+ rho=7.9565d0-6.4761d0*x+5.5283d0*x*x-3.0807d0*x*x*x
+ vp=24.9520d0-40.4673d0*x+51.4832d0*x*x-26.6419d0*x*x*x
+ vs=11.1671d0-13.7818d0*x+17.4575d0*x*x-9.2777d0*x*x*x
+ Qmu=312.0d0
+ Qkappa=57827.0d0
+ else if(r > R771 .and. r <= R670) then
+ drhodr=-6.4761d0+2.0d0*5.5283d0*x-3.0d0*3.0807d0*x*x
+ rho=7.9565d0-6.4761d0*x+5.5283d0*x*x-3.0807d0*x*x*x
+ vp=29.2766d0-23.6027d0*x+5.5242d0*x*x-2.5514d0*x*x*x
+ vs=22.3459d0-17.2473d0*x-2.0834d0*x*x+0.9783d0*x*x*x
+ Qmu=312.0d0
+ Qkappa=57827.0d0
+!
+!--- mantle: above d670
+!
+ else if(r > R670 .and. r <= R600) then
+ drhodr=-1.4836d0
+ rho=5.3197d0-1.4836d0*x
+ vp=19.0957d0-9.8672d0*x
+ vs=9.9839d0-4.9324d0*x
+ Qmu=143.0d0
+ Qkappa=57827.0d0
+ else if(r > R600 .and. r <= R400) then
+ drhodr=-8.0298d0
+ rho=11.2494d0-8.0298d0*x
+ vp=39.7027d0-32.6166d0*x
+ vs=22.3512d0-18.5856d0*x
+ Qmu=143.0d0
+ Qkappa=57827.0d0
+ else if(r > R400 .and. r <= R220) then
+ drhodr=-3.8045d0
+ rho=7.1089d0-3.8045d0*x
+ vp=20.3926d0-12.2569d0*x
+ vs=8.9496d0-4.4597d0*x
+ Qmu=143.0d0
+ Qkappa=57827.0d0
+ else if(r > R220 .and. r <= R80) then
+ drhodr=0.6924d0
+ rho=2.6910d0+0.6924d0*x
+ vp=4.1875d0+3.9382d0*x
+ vs=2.1519d0+2.3481d0*x
+ Qmu=80.0d0
+ Qkappa=57827.0d0
+ else
+ if(CRUSTAL .and. .not. SUPPRESS_CRUSTAL_MESH) then
+! fill with PREM mantle and later add CRUST2.0
+ if(r > R80) then
+ drhodr=0.6924d0
+ rho=2.6910d0+0.6924d0*x
+ vp=4.1875d0+3.9382d0*x
+ vs=2.1519d0+2.3481d0*x
+ Qmu=600.0d0
+ Qkappa=57827.0d0
+ endif
+ else
+! use PREM crust
+ if(r > R80 .and. r <= RMOHO) then
+ drhodr=0.6924d0
+ rho=2.6910d0+0.6924d0*x
+ vp=4.1875d0+3.9382d0*x
+ vs=2.1519d0+2.3481d0*x
+ Qmu=600.0d0
+ Qkappa=57827.0d0
+
+
+ else if (SUPPRESS_CRUSTAL_MESH) then
+!! DK DK extend the Moho up to the surface instead of the crust
+ drhodr=0.6924d0
+ rho = 2.6910d0+0.6924d0*(RMOHO / R_EARTH)
+ vp = 4.1875d0+3.9382d0*(RMOHO / R_EARTH)
+ vs = 2.1519d0+2.3481d0*(RMOHO / R_EARTH)
+ Qmu=600.0d0
+ Qkappa=57827.0d0
+
+ else if(r > RMOHO .and. r <= RMIDDLE_CRUST) then
+ drhodr=0.0d0
+ rho=2.9d0
+ vp=6.8d0
+ vs=3.9d0
+ Qmu=600.0d0
+ Qkappa=57827.0d0
+
+! same properties everywhere in PREM crust if we decide to define only one layer in the crust
+ if(ONE_CRUST) then
+ drhodr=0.0d0
+ rho=2.6d0
+ vp=5.8d0
+ vs=3.2d0
+ Qmu=600.0d0
+ Qkappa=57827.0d0
+ endif
+
+ else if(r > RMIDDLE_CRUST .and. r <= ROCEAN) then
+ drhodr=0.0d0
+ rho=2.6d0
+ vp=5.8d0
+ vs=3.2d0
+ Qmu=600.0d0
+ Qkappa=57827.0d0
+! for density profile for gravity, we do not check that r <= R_EARTH
+ else if(r > ROCEAN) then
+ drhodr=0.0d0
+ rho=2.6d0
+ vp=5.8d0
+ vs=3.2d0
+ Qmu=600.0d0
+ Qkappa=57827.0d0
+
+ endif
+ endif
+ endif
+
+! non-dimensionalize
+! time scaling (s^{-1}) is done with scaleval
+ scaleval=dsqrt(PI*GRAV*RHOAV)
+ drhodr=drhodr*1000.0d0/RHOAV
+ rho=rho*1000.0d0/RHOAV
+ vp=vp*1000.0d0/(R_EARTH*scaleval)
+ vs=vs*1000.0d0/(R_EARTH*scaleval)
+
+ end subroutine prem_iso
+
+!
+!=====================================================================
+!
+
+ subroutine prem_display_outer_core(myrank,x,rho,vp,vs,Qkappa,Qmu,idoubling)
+
+! routine used for AVS or DX display of stability condition
+! and number of points per wavelength only in the fluid outer core
+
+ implicit none
+
+ include "constants.h"
+
+! given a normalized radius x, gives the non-dimesionalized density rho,
+! speeds vp and vs, and the quality factors Qkappa and Qmu
+
+ integer idoubling,myrank
+ double precision x,rho,vp,vs,Qkappa,Qmu
+
+ double precision scaleval
+
+ if(idoubling /= IFLAG_OUTER_CORE_NORMAL) call exit_MPI(myrank,'wrong doubling flag for outer core point')
+
+!
+!--- outer core
+!
+ rho=12.5815d0-1.2638d0*x-3.6426d0*x*x-5.5281d0*x*x*x
+ vp=11.0487d0-4.0362d0*x+4.8023d0*x*x-13.5732d0*x*x*x
+ vs=0.0d0
+ Qmu=0.0d0
+ Qkappa=57827.0d0
+
+! non-dimensionalize
+! time scaling (s^{-1}) is done with scaleval
+ scaleval = dsqrt(PI*GRAV*RHOAV)
+ rho = rho*1000.0d0/RHOAV
+ vp = vp*1000.0d0/(R_EARTH*scaleval)
+ vs = vs*1000.0d0/(R_EARTH*scaleval)
+
+ end subroutine prem_display_outer_core
+
+!
+!=====================================================================
+!
+
+ subroutine prem_density(x,rho,ONE_CRUST,RICB,RCMB,RTOPDDOUBLEPRIME, &
+ R600,R670,R220,R771,R400,R80,RMOHO,RMIDDLE_CRUST,ROCEAN)
+
+ implicit none
+
+ include "constants.h"
+
+ double precision x,rho,RICB,RCMB,RTOPDDOUBLEPRIME, &
+ R600,R670,R220,R771,R400,R80,RMOHO,RMIDDLE_CRUST,ROCEAN
+
+ logical ONE_CRUST
+
+ double precision r
+
+ r = x * R_EARTH
+
+ if(r <= RICB) then
+ rho=13.0885d0-8.8381d0*x*x
+ else if(r > RICB .and. r <= RCMB) then
+ rho=12.5815d0-1.2638d0*x-3.6426d0*x*x-5.5281d0*x*x*x
+ else if(r > RCMB .and. r <= RTOPDDOUBLEPRIME) then
+ rho=7.9565d0-6.4761d0*x+5.5283d0*x*x-3.0807d0*x*x*x
+ else if(r > RTOPDDOUBLEPRIME .and. r <= R771) then
+ rho=7.9565d0-6.4761d0*x+5.5283d0*x*x-3.0807d0*x*x*x
+ else if(r > R771 .and. r <= R670) then
+ rho=7.9565d0-6.4761d0*x+5.5283d0*x*x-3.0807d0*x*x*x
+ else if(r > R670 .and. r <= R600) then
+ rho=5.3197d0-1.4836d0*x
+ else if(r > R600 .and. r <= R400) then
+ rho=11.2494d0-8.0298d0*x
+ else if(r > R400 .and. r <= R220) then
+ rho=7.1089d0-3.8045d0*x
+ else if(r > R220 .and. r <= R80) then
+ rho=2.6910d0+0.6924d0*x
+ else
+ if(r > R80 .and. r <= RMOHO) then
+ rho=2.6910d0+0.6924d0*x
+ else if(r > RMOHO .and. r <= RMIDDLE_CRUST) then
+ if(ONE_CRUST) then
+ rho=2.6d0
+ else
+ rho=2.9d0
+ endif
+ else if(r > RMIDDLE_CRUST .and. r <= ROCEAN) then
+ rho=2.6d0
+ else if(r > ROCEAN) then
+ rho=2.6d0
+ endif
+ endif
+
+ rho=rho*1000.0d0/RHOAV
+
+ end subroutine prem_density
+
Copied: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/reduce.f90 (from rev 13385, seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/reduce.f90)
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/reduce.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/reduce.f90 2008-11-26 01:37:52 UTC (rev 13411)
@@ -0,0 +1,84 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 4 . 0
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory, California Institute of Technology, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) California Institute of Technology and University of Pau / CNRS / INRIA
+! February 2008
+!
+! 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 reduce(theta,phi)
+
+! bring theta between 0 and PI, and phi between 0 and 2*PI
+
+ implicit none
+
+ include "constants.h"
+
+ double precision theta,phi
+
+ integer i
+ double precision th,ph
+
+ th=theta
+ ph=phi
+ i=abs(int(ph/TWO_PI))
+ if(ph<ZERO) then
+ ph=ph+(i+1)*TWO_PI
+ else
+ if(ph>TWO_PI) ph=ph-i*TWO_PI
+ endif
+ phi=ph
+ if(th<ZERO .or. th>PI) then
+ i=int(th/PI)
+ if(th>ZERO) then
+ if(mod(i,2) /= 0) then
+ th=(i+1)*PI-th
+ if(ph<PI) then
+ ph=ph+PI
+ else
+ ph=ph-PI
+ endif
+ else
+ th=th-i*PI
+ endif
+ else
+ if(mod(i,2) == 0) then
+ th=-th+i*PI
+ if(ph<PI) then
+ ph=ph+PI
+ else
+ ph=ph-PI
+ endif
+ else
+ th=th-i*PI
+ endif
+ endif
+ theta=th
+ phi=ph
+ endif
+
+ if(theta<ZERO .or. theta>PI) stop 'theta out of range in reduce'
+
+ if(phi<ZERO .or. phi>TWO_PI) stop 'phi out of range in reduce'
+
+ end subroutine reduce
+
Copied: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/rthetaphi_xyz.f90 (from rev 13385, seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/rthetaphi_xyz.f90)
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/rthetaphi_xyz.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/chino/rthetaphi_xyz.f90 2008-11-26 01:37:52 UTC (rev 13411)
@@ -0,0 +1,119 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 4 . 0
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory, California Institute of Technology, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) California Institute of Technology and University of Pau / CNRS / INRIA
+! February 2008
+!
+! 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 xyz_2_rthetaphi(x,y,z,r,theta,phi)
+
+! convert x y z to r theta phi, single precision call
+
+ implicit none
+
+ include "constants.h"
+
+ real(kind=CUSTOM_REAL) x,y,z,r,theta,phi
+ double precision xmesh,ymesh,zmesh
+
+! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+
+ xmesh = dble(x)
+ ymesh = dble(y)
+ zmesh = dble(z)
+
+ if(zmesh > -SMALL_VAL_ANGLE .and. zmesh <= ZERO) zmesh = -SMALL_VAL_ANGLE
+ if(zmesh < SMALL_VAL_ANGLE .and. zmesh >= ZERO) zmesh = SMALL_VAL_ANGLE
+ theta = sngl(datan2(dsqrt(xmesh*xmesh+ymesh*ymesh),zmesh))
+ if(xmesh > -SMALL_VAL_ANGLE .and. xmesh <= ZERO) xmesh = -SMALL_VAL_ANGLE
+ if(xmesh < SMALL_VAL_ANGLE .and. xmesh >= ZERO) xmesh = SMALL_VAL_ANGLE
+ phi = sngl(datan2(ymesh,xmesh))
+
+ r = sngl(dsqrt(xmesh**2 + ymesh**2 + zmesh**2))
+
+ else
+
+ xmesh = x
+ ymesh = y
+ zmesh = z
+
+ if(zmesh > -SMALL_VAL_ANGLE .and. zmesh <= ZERO) zmesh = -SMALL_VAL_ANGLE
+ if(zmesh < SMALL_VAL_ANGLE .and. zmesh >= ZERO) zmesh = SMALL_VAL_ANGLE
+ theta = datan2(dsqrt(xmesh*xmesh+ymesh*ymesh),zmesh)
+ if(xmesh > -SMALL_VAL_ANGLE .and. xmesh <= ZERO) xmesh = -SMALL_VAL_ANGLE
+ if(xmesh < SMALL_VAL_ANGLE .and. xmesh >= ZERO) xmesh = SMALL_VAL_ANGLE
+ phi = datan2(ymesh,xmesh)
+
+ r = dsqrt(xmesh**2 + ymesh**2 + zmesh**2)
+
+ endif
+
+ end subroutine xyz_2_rthetaphi
+
+!-------------------------------------------------------------
+
+ subroutine xyz_2_rthetaphi_dble(x,y,z,r,theta,phi)
+
+! convert x y z to r theta phi, double precision call
+
+ implicit none
+
+ include "constants.h"
+
+ double precision x,y,z,r,theta,phi
+ double precision xmesh,ymesh,zmesh
+
+ xmesh = x
+ ymesh = y
+ zmesh = z
+
+ if(zmesh > -SMALL_VAL_ANGLE .and. zmesh <= ZERO) zmesh = -SMALL_VAL_ANGLE
+ if(zmesh < SMALL_VAL_ANGLE .and. zmesh >= ZERO) zmesh = SMALL_VAL_ANGLE
+ theta = datan2(dsqrt(xmesh*xmesh+ymesh*ymesh),zmesh)
+ if(xmesh > -SMALL_VAL_ANGLE .and. xmesh <= ZERO) xmesh = -SMALL_VAL_ANGLE
+ if(xmesh < SMALL_VAL_ANGLE .and. xmesh >= ZERO) xmesh = SMALL_VAL_ANGLE
+ phi = datan2(ymesh,xmesh)
+
+ r = dsqrt(xmesh**2 + ymesh**2 + zmesh**2)
+
+ end subroutine xyz_2_rthetaphi_dble
+
+!-------------------------------------------------------------
+
+ subroutine rthetaphi_2_xyz(x,y,z,r,theta,phi)
+
+! convert r theta phi to x y z
+
+ implicit none
+
+ include "constants.h"
+
+ real(kind=CUSTOM_REAL) x,y,z,r,theta,phi
+
+ x = r * sin(theta) * cos(phi)
+ y = r * sin(theta) * sin(phi)
+ z = r * cos(theta)
+
+ end subroutine rthetaphi_2_xyz
+
Added: seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/setup.py
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/setup.py (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/pluggable/MODELS/harness/setup.py 2008-11-26 01:37:52 UTC (rev 13411)
@@ -0,0 +1,18 @@
+
+from distutils.core import setup
+
+from chino.harness import harnessSourceFiles, harnessHeaderFiles
+data = []
+data.extend(harnessSourceFiles)
+data.extend(harnessHeaderFiles)
+
+setup(
+ name = 'chino',
+ version = '1.0.0',
+ url = 'http://www.geodynamics.org/',
+ author = 'Leif Strand',
+ author_email = 'leif at geodynamics.org',
+ packages = ['chino'],
+ package_data = {'chino': data},
+ scripts = ['bin/chino'],
+ )
More information about the CIG-COMMITS
mailing list