[cig-commits] r17232 - in seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT: . create_movie_GMT create_movie_GMT/DATA
liuqy at geodynamics.org
liuqy at geodynamics.org
Fri Oct 1 08:00:21 PDT 2010
Author: liuqy
Date: 2010-10-01 08:00:21 -0700 (Fri, 01 Oct 2010)
New Revision: 17232
Added:
seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/
seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/DATA/
seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/DATA/CMTSOLUTION
seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/OUTPUT_FILES/
seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/Par_file
seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/compute_parameters.f90
seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/constants.h
seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/convert_movie_real_to_double.f90
seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/convert_nonlinear.f90
seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/convert_tiff_topo.pl
seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/create_movie_GMT.f90
seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/get_value_parameters.f90
seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/makefile
seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/movie2gif.pl
seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/obj/
seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/read_parameter_file.f90
seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/read_value_parameters.f90
seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/readme.txt
seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/utm_geo.f90
Log:
Add old create_movie_GMT package for movie2gif.pl
Added: seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/DATA/CMTSOLUTION
===================================================================
--- seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/DATA/CMTSOLUTION (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/DATA/CMTSOLUTION 2010-10-01 15:00:21 UTC (rev 17232)
@@ -0,0 +1,13 @@
+PDE 2010 6 15 4 26 58.00 32.6978 -115.9235 8.0 5.7 5.7 10704013 from SCSN
+event name: 10704013
+time shift: 0.0000
+half duration: 2.5
+latitude: 32.6978
+longitude: -115.9235
+depth: 8.0000
+Mrr: -9.238200e+23
+Mtt: -4.114500e+24
+Mpp: 5.038300e+24
+Mrt: -8.087700e+23
+Mrp: 2.265500e+23
+Mtp: -9.079000e+23
Added: seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/Par_file
===================================================================
--- seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/Par_file (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/Par_file 2010-10-01 15:00:21 UTC (rev 17232)
@@ -0,0 +1,68 @@
+
+# forward or adjoint simulation
+SIMULATION_TYPE = 1 # 1 = forward, 2 = adjoint, 3 = both simultaneously
+SAVE_FORWARD = .false.
+
+# coordinates of mesh block in latitude/longitude and depth in km
+LATITUDE_MIN = 32.2d0
+LATITUDE_MAX = 36.8d0
+LONGITUDE_MIN = -120.3d0
+LONGITUDE_MAX = -114.7d0
+DEPTH_BLOCK_KM = 60.d0
+UTM_PROJECTION_ZONE = 11
+SUPPRESS_UTM_PROJECTION = .false.
+
+# number of elements at the surface along edges of the mesh at the surface
+# (must be 8 * multiple of NPROC below if mesh is not regular and contains mesh doublings)
+# (must be multiple of NPROC below if mesh is regular)
+NEX_XI = 288
+NEX_ETA = 288
+
+# number of MPI processors along xi and eta (can be different)
+NPROC_XI = 12
+NPROC_ETA = 12
+
+# model (SoCal, Harvard_LA, Min_Chen_anisotropy)
+MODEL = Harvard_LA
+
+# parameters describing the model
+OCEANS = .true.
+TOPOGRAPHY = .true.
+ATTENUATION = .true.
+USE_OLSEN_ATTENUATION = .false.
+
+# absorbing boundary conditions for a regional simulation
+ABSORBING_CONDITIONS = .true.
+
+# record length in seconds
+RECORD_LENGTH_IN_SECONDS = 180
+
+# save AVS or OpenDX movies
+MOVIE_SURFACE = .true.
+MOVIE_VOLUME = .false.
+NTSTEP_BETWEEN_FRAMES = 200
+CREATE_SHAKEMAP = .true.
+SAVE_DISPLACEMENT = .false.
+USE_HIGHRES_FOR_MOVIES = .false.
+HDUR_MOVIE = 0.0
+
+# save AVS or OpenDX mesh files to check the mesh
+SAVE_MESH_FILES = .false.
+
+# path to store the local database file on each node
+LOCAL_PATH = /state/partition1/scratch/friberg/DATABASES_MPI.123319
+
+# machine file for MPI
+# this is not needed for new cluster and is in go_mesher/go_solver
+# MACHINE_FILE = mymachines
+
+# interval at which we output time step info and max of norm of displacement
+NTSTEP_BETWEEN_OUTPUT_INFO = 100
+
+# interval in time steps for writing of seismograms
+NTSTEP_BETWEEN_OUTPUT_SEISMOS = 200000
+
+# print source time function
+PRINT_SOURCE_TIME_FUNCTION = .false.
+
+
Added: seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/compute_parameters.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/compute_parameters.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/compute_parameters.f90 2010-10-01 15:00:21 UTC (rev 17232)
@@ -0,0 +1,264 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology September 2006
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+ subroutine compute_parameters(NER,NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA, &
+ NPROC,NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
+ NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM,NER_SEDIM, &
+ NSPEC_AB,NSPEC2D_A_XI,NSPEC2D_B_XI, &
+ NSPEC2D_A_ETA,NSPEC2D_B_ETA, &
+ NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
+ NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NGLOB_AB,USE_REGULAR_MESH)
+
+ implicit none
+
+ include "constants.h"
+
+! parameters read from parameter file
+ integer NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA
+ integer NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM,NER_SEDIM
+
+! parameters to be computed based upon parameters above read from file
+ integer NPROC,NEX_PER_PROC_XI,NEX_PER_PROC_ETA,NER
+
+ integer NSPEC_AB,NSPEC2D_A_XI,NSPEC2D_B_XI, &
+ NSPEC2D_A_ETA,NSPEC2D_B_ETA, &
+ NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
+ NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NGLOB_AB
+
+ integer NEX_DOUBLING_SEDIM_XI,NEX_DOUBLING_SEDIM_ETA
+ integer NEX_DOUBLING_SEDIM_PER_PROC_XI,NEX_DOUBLING_SEDIM_PER_PROC_ETA
+ integer NSPEC2D_DOUBLING_A_XI,NSPEC2D_DOUBLING_A_ETA
+ integer NSPEC2D_DOUBLING_B_XI,NSPEC2D_DOUBLING_B_ETA
+ integer NSPEC_DOUBLING_AB
+ integer NUM_DOUBLING_BRICKS
+ integer NUM2D_DOUBLING_BRICKS_XI,NUM2D_DOUBLING_BRICKS_ETA
+ integer nglob_no_doubling_volume,nglob_no_doubling_surface
+ integer nblocks_xi,nblocks_eta
+ integer nglob_surface_typeA,nglob_surface_typeB
+ integer NSPEC1D_RADIAL_BEDROCK,NPOIN1D_RADIAL_BEDROCK
+
+ integer NSPEC_NO_DOUBLING,NSPEC2D_NO_DOUBLING_XI,NSPEC2D_NO_DOUBLING_ETA
+
+ logical USE_REGULAR_MESH
+
+!
+!--- case of a regular mesh
+!
+ if(USE_REGULAR_MESH) then
+
+! total number of spectral elements along radius
+ NER = NER_SEDIM
+
+! number of elements horizontally in each slice (i.e. per processor)
+! these two values MUST be equal in all cases
+ NEX_PER_PROC_XI = NEX_XI / NPROC_XI
+ NEX_PER_PROC_ETA = NEX_ETA / NPROC_ETA
+
+! total number of processors in each of the six chunks
+ NPROC = NPROC_XI * NPROC_ETA
+
+! exact number of spectral elements without doubling layers
+ NSPEC_NO_DOUBLING = NEX_XI*NEX_ETA*NER_SEDIM
+
+! %%%%%%%%%%%%%% surface elements %%%%%%%%%%%%%%%%%%%
+
+! exact number of surface elements for a chunk without doubling layers
+
+ NSPEC2D_NO_DOUBLING_XI = NEX_PER_PROC_XI*NER_SEDIM
+
+ NSPEC2D_NO_DOUBLING_ETA = NEX_PER_PROC_ETA*NER_SEDIM
+
+! exact number of spectral elements
+ NSPEC_AB = NSPEC_NO_DOUBLING / NPROC
+
+! exact number of surface elements for faces A and B along XI and ETA
+ NSPEC2D_A_XI = NSPEC2D_NO_DOUBLING_XI
+ NSPEC2D_B_XI = NSPEC2D_NO_DOUBLING_XI
+ NSPEC2D_A_ETA = NSPEC2D_NO_DOUBLING_ETA
+ NSPEC2D_B_ETA = NSPEC2D_NO_DOUBLING_ETA
+
+! exact number of surface elements on the bottom and top boundaries
+! and theoretical number of spectral elements in radial direction
+
+ NSPEC2D_TOP = NEX_XI*NEX_ETA / NPROC
+ NSPEC2D_BOTTOM = NSPEC2D_TOP
+
+ NSPEC1D_RADIAL_BEDROCK = NER
+
+! face with max number of elements is type B here
+! maximum number of surface elements on vertical boundaries of the slices
+ NSPEC2DMAX_XMIN_XMAX = NSPEC2D_B_ETA
+ NSPEC2DMAX_YMIN_YMAX = NSPEC2D_B_XI
+
+! theoretical number of Gauss-Lobatto points in radial direction
+ NPOIN1D_RADIAL_BEDROCK = NSPEC1D_RADIAL_BEDROCK*(NGLLZ-1)+1
+
+! 2-D addressing and buffers for summation between slices
+! we add one to number of points because of the flag after the last point
+ NPOIN2DMAX_XMIN_XMAX = NSPEC2DMAX_XMIN_XMAX*NGLLY*NGLLZ + 1
+ NPOIN2DMAX_YMIN_YMAX = NSPEC2DMAX_YMIN_YMAX*NGLLX*NGLLZ + 1
+
+! exact number of global points
+ NGLOB_AB = (NEX_PER_PROC_XI*(NGLLX-1)+1) * (NEX_PER_PROC_ETA*(NGLLY-1)+1) * (NER*(NGLLZ-1)+1)
+
+!
+!--- case of a non-regular mesh with mesh doublings
+!
+ else
+
+! total number of spectral elements along radius
+ NER = NER_BOTTOM_MOHO + NER_MOHO_16 + NER_16_BASEMENT + NER_BASEMENT_SEDIM + NER_SEDIM
+
+! number of elements horizontally in each slice (i.e. per processor)
+! these two values MUST be equal in all cases
+ NEX_PER_PROC_XI = NEX_XI / NPROC_XI
+ NEX_PER_PROC_ETA = NEX_ETA / NPROC_ETA
+
+! total number of processors in each of the six chunks
+ NPROC = NPROC_XI * NPROC_ETA
+
+! number of spectral elements at the bottom of the doubling below the moho
+ NEX_DOUBLING_SEDIM_XI=NEX_XI/2
+ NEX_DOUBLING_SEDIM_ETA=NEX_ETA/2
+ NEX_DOUBLING_SEDIM_PER_PROC_XI=NEX_PER_PROC_XI/2
+ NEX_DOUBLING_SEDIM_PER_PROC_ETA=NEX_PER_PROC_ETA/2
+
+! exact number of spectral elements without doubling layers
+ NSPEC_NO_DOUBLING = &
+ (NEX_DOUBLING_SEDIM_XI*NEX_DOUBLING_SEDIM_ETA*(NER_BASEMENT_SEDIM/2-3) &
+ +(NEX_XI/4)*(NEX_ETA/4)*(NER_16_BASEMENT/2-3) &
+ +(NEX_XI/4)*(NEX_ETA/4)*(NER_MOHO_16/2) &
+ +(NEX_XI/4)*(NEX_ETA/4)*(NER_BOTTOM_MOHO/4)) + NEX_XI*NEX_ETA*NER_SEDIM
+
+! exact number of spectral elements in the doubling regions
+
+! number of elementary bricks in the two regions with doubling
+ NUM_DOUBLING_BRICKS = ((NEX_XI/4)*(NEX_ETA/4) &
+ +NEX_DOUBLING_SEDIM_XI*NEX_DOUBLING_SEDIM_ETA)/4
+
+! for type AB, each doubling brick contains 40 elements on 3 levels
+ NSPEC_DOUBLING_AB=40*NUM_DOUBLING_BRICKS
+
+! %%%%%%%%%%%%%% surface elements %%%%%%%%%%%%%%%%%%%
+
+! exact number of surface elements for a chunk without doubling layers
+
+ NSPEC2D_NO_DOUBLING_XI = &
+ NEX_DOUBLING_SEDIM_PER_PROC_XI*(NER_BASEMENT_SEDIM/2-3) &
+ +(NEX_PER_PROC_XI/4)*(NER_16_BASEMENT/2-3) &
+ +(NEX_PER_PROC_XI/4)*(NER_MOHO_16/2) &
+ +(NEX_PER_PROC_XI/4)*(NER_BOTTOM_MOHO/4) + NEX_PER_PROC_XI*NER_SEDIM
+
+ NSPEC2D_NO_DOUBLING_ETA = &
+ NEX_DOUBLING_SEDIM_PER_PROC_ETA*(NER_BASEMENT_SEDIM/2-3) &
+ +(NEX_PER_PROC_ETA/4)*(NER_16_BASEMENT/2-3) &
+ +(NEX_PER_PROC_ETA/4)*(NER_MOHO_16/2) &
+ +(NEX_PER_PROC_ETA/4)*(NER_BOTTOM_MOHO/4) + NEX_PER_PROC_ETA*NER_SEDIM
+
+! exact number of surface elements in the doubling regions
+
+! number of elementary bricks in the two regions with doubling
+ NUM2D_DOUBLING_BRICKS_XI = ((NEX_PER_PROC_XI/4) &
+ +NEX_DOUBLING_SEDIM_PER_PROC_XI)/2
+
+ NUM2D_DOUBLING_BRICKS_ETA = ((NEX_PER_PROC_ETA/4) &
+ +NEX_DOUBLING_SEDIM_PER_PROC_ETA)/2
+
+! for type A, each doubling brick contains 10 elements on 3 levels
+ NSPEC2D_DOUBLING_A_XI=10*NUM2D_DOUBLING_BRICKS_XI
+ NSPEC2D_DOUBLING_A_ETA=10*NUM2D_DOUBLING_BRICKS_ETA
+
+! for type B, each doubling brick contains 12 elements on 3 levels
+ NSPEC2D_DOUBLING_B_XI=12*NUM2D_DOUBLING_BRICKS_XI
+ NSPEC2D_DOUBLING_B_ETA=12*NUM2D_DOUBLING_BRICKS_ETA
+
+! exact number of spectral elements
+ NSPEC_AB = (NSPEC_NO_DOUBLING + NSPEC_DOUBLING_AB) / NPROC
+
+! exact number of surface elements for faces A and B
+! along XI and ETA for doubling region
+ NSPEC2D_A_XI = NSPEC2D_NO_DOUBLING_XI + NSPEC2D_DOUBLING_A_XI
+ NSPEC2D_B_XI = NSPEC2D_NO_DOUBLING_XI + NSPEC2D_DOUBLING_B_XI
+ NSPEC2D_A_ETA = NSPEC2D_NO_DOUBLING_ETA + NSPEC2D_DOUBLING_A_ETA
+ NSPEC2D_B_ETA = NSPEC2D_NO_DOUBLING_ETA + NSPEC2D_DOUBLING_B_ETA
+
+! exact number of surface elements on the bottom and top boundaries
+! and theoretical number of spectral elements in radial direction
+
+ NSPEC2D_TOP = NEX_XI*NEX_ETA / NPROC
+ NSPEC2D_BOTTOM = (NEX_XI/4)*(NEX_ETA/4) / NPROC
+
+ NSPEC1D_RADIAL_BEDROCK = (NER_BASEMENT_SEDIM+NER_16_BASEMENT+NER_MOHO_16)/2 + NER_BOTTOM_MOHO/4
+
+! face with max number of elements is type B here
+! maximum number of surface elements on vertical boundaries of the slices
+ NSPEC2DMAX_XMIN_XMAX = NSPEC2D_B_ETA
+ NSPEC2DMAX_YMIN_YMAX = NSPEC2D_B_XI
+
+! theoretical number of Gauss-Lobatto points in radial direction
+ NPOIN1D_RADIAL_BEDROCK = NSPEC1D_RADIAL_BEDROCK*(NGLLZ-1)+1
+
+! 2-D addressing and buffers for summation between slices
+! we add one to number of points because of the flag after the last point
+ NPOIN2DMAX_XMIN_XMAX = NSPEC2DMAX_XMIN_XMAX*NGLLY*NGLLZ + 1
+ NPOIN2DMAX_YMIN_YMAX = NSPEC2DMAX_YMIN_YMAX*NGLLX*NGLLZ + 1
+
+! exact number of global points
+
+! case of the doubling regions
+! formulas computed using Mathematica for the basic 3D doubling brick
+
+! compute number of points in blocks with no doubling
+! exclude the three surfaces in contact with the doubling regions
+ nglob_no_doubling_volume = (4*(NGLLX-1)+1)*(4*(NGLLX-1)+1)*((NER_BASEMENT_SEDIM/2-3 )*(NGLLX-1)-1) &
+ +(2*(NGLLX-1)+1)*(2*(NGLLX-1)+1)*(((NER_16_BASEMENT/2+NER_MOHO_16/2+NER_BOTTOM_MOHO/4)-3)*(NGLLX-1)+0)
+
+! number of basic blocks in each slice
+ nblocks_xi = NEX_PER_PROC_XI / 8
+ nblocks_eta = NEX_PER_PROC_ETA / 8
+
+ NGLOB_AB = nblocks_xi*nblocks_eta*(200*NGLLX**3 - 484*NGLLX**2 + 392*NGLLX - 106 + nglob_no_doubling_volume)
+
+! same thing for 2D surfaces for the three types of faces
+ nglob_no_doubling_surface = (4*(NGLLX-1)+1)*((NER_BASEMENT_SEDIM/2-3)*(NGLLX-1)-1) &
+ +(2*(NGLLX-1)+1)*(((NER_16_BASEMENT/2+NER_MOHO_16/2+NER_BOTTOM_MOHO/4)-3)*(NGLLX-1)+0)
+
+ nglob_surface_typeA = 30*NGLLX**2 - 45 * NGLLX + 17
+ nglob_surface_typeB = 36*NGLLX**2 - 57 * NGLLX + 23
+
+! final number of points in volume obtained by removing planes counted twice
+ NGLOB_AB = NGLOB_AB &
+ - (nblocks_xi-1)*nblocks_eta*(nglob_surface_typeA + nglob_no_doubling_surface) &
+ - (nblocks_eta-1)*nblocks_xi*(nglob_surface_typeB + nglob_no_doubling_surface) &
+ + (nblocks_eta-1)*(nblocks_xi-1)*NPOIN1D_RADIAL_BEDROCK
+
+! add number of points in the sediments
+ NGLOB_AB = NGLOB_AB + (NEX_PER_PROC_XI*(NGLLX-1)+1) &
+ *(NEX_PER_PROC_ETA*(NGLLY-1)+1)*(NER_SEDIM*(NGLLZ-1)+0)
+
+ endif ! end of section for non-regular mesh with doublings
+
+ end subroutine compute_parameters
+
Added: seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/constants.h
===================================================================
--- seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/constants.h (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/constants.h 2010-10-01 15:00:21 UTC (rev 17232)
@@ -0,0 +1,357 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology July 2005
+!
+! 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.
+
+!
+! 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.
+
+! use inlined products of Deville et al. (2002) to speedup the calculations to compute internal forces
+! logical, parameter :: USE_DEVILLE_PRODUCTS = .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
+
+! for optimized routines by Deville et al. (2002)
+ integer, parameter :: m1 = NGLLX, m2 = NGLLX * NGLLY
+
+! input, output and main MPI I/O files
+ integer, parameter :: ISTANDARD_OUTPUT = 6
+ integer, parameter :: IIN = 40,IOUT = 41
+ integer, parameter :: IIN_NOISE = 43,IOUT_NOISE = 44
+ integer, parameter :: NOISE_TOMOGRAPHY = 0
+! 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
+ double precision, parameter :: DENSITY_MAX = 3000.d0
+ double precision, parameter :: DENSITY_MIN = 2000.d0
+
+! density of sea water
+ real(kind=CUSTOM_REAL), parameter :: RHO_OCEANS = 1020.0
+
+! extend 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
+
+! 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
+
+! ---------------------------------------------------------------------------------------
+! 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 = .false.
+
+! 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/DATABASES_MPI_Q/'
+
+!------------------------------------------------------
+!----------- 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 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
+!
+
+! 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
+
+ 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
+
+! 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/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/convert_movie_real_to_double.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/convert_movie_real_to_double.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/convert_movie_real_to_double.f90 2010-10-01 15:00:21 UTC (rev 17232)
@@ -0,0 +1,205 @@
+program convert_movie_real_to_double
+
+ implicit none
+
+ include 'constants.h'
+
+ character(len=100) par_file, movie_data_prefix, start_frame, end_frame, &
+ output_file_prefix
+ integer ios1, ios2, nspectot_AVS_max
+ ! threshold in percent of the maximum below which we cut the amplitude
+ logical, parameter :: APPLY_THRESHOLD = .false.
+ real(kind=CUSTOM_REAL), parameter :: THRESHOLD = 1._CUSTOM_REAL / 100._CUSTOM_REAL
+
+ ! coefficient of power law used for non linear scaling
+ logical, parameter :: NONLINEAR_SCALING = .false.
+ real(kind=CUSTOM_REAL), parameter :: POWER_SCALING = 0.50_CUSTOM_REAL
+
+ integer it,it1,it2,ivalue,ispec
+ integer iformat,nframes,iframe,inumber,inorm,iscaling_shake
+ integer ibool_number,ibool_number1,ibool_number2,ibool_number3,ibool_number4
+
+ logical UNIQUE_FILE,plot_shaking_map
+
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: x,y,z,display
+ real(kind=CUSTOM_REAL) xcoord,ycoord,zcoord
+ real(kind=CUSTOM_REAL) vectorx,vectory,vectorz
+
+ double precision min_field_current,max_field_current,max_absol
+
+ character(len=150) outputname,outputname1
+
+ integer iproc,ipoin
+
+ ! GMT
+ double precision lat,long,zscaling
+ integer igmt
+
+ ! for sorting routine
+ integer npointot,ilocnum,nglob,i,j,ielm,ieoff,ispecloc
+ integer, dimension(:), allocatable :: iglob,loc,ireorder
+ logical, dimension(:), allocatable :: ifseg,mask_point
+ double precision, dimension(:), allocatable :: xp,yp,zp,xp_save,yp_save,zp_save,field_display
+
+ ! movie files stored by solver
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: store_val
+ real*8, dimension(:,:), allocatable :: store_val_double
+
+! parameters read from parameter file
+ integer NER_SEDIM,NER_BASEMENT_SEDIM,NER_16_BASEMENT, &
+ NER_MOHO_16,NER_BOTTOM_MOHO,NEX_ETA,NEX_XI, &
+ NPROC_ETA,NPROC_XI,NTSTEP_BETWEEN_OUTPUT_SEISMOS,NSTEP,UTM_PROJECTION_ZONE
+ integer NSOURCES
+
+ double precision UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK
+ double precision DT,LATITUDE_MIN,LATITUDE_MAX,LONGITUDE_MIN,LONGITUDE_MAX
+ double precision THICKNESS_TAPER_BLOCK_HR,THICKNESS_TAPER_BLOCK_MR,VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM
+
+ logical HARVARD_3D_GOCAD_MODEL,TOPOGRAPHY,ATTENUATION,USE_OLSEN_ATTENUATION, &
+ OCEANS,IMPOSE_MINIMUM_VP_GOCAD,HAUKSSON_REGIONAL_MODEL, &
+ BASEMENT_MAP,MOHO_MAP_LUPEI,ABSORBING_CONDITIONS
+ logical ANISOTROPY,SAVE_AVS_DX_MESH_FILES,PRINT_SOURCE_TIME_FUNCTION
+
+ logical MOVIE_SURFACE,MOVIE_VOLUME,CREATE_SHAKEMAP,SAVE_DISPLACEMENT,USE_HIGHRES_FOR_MOVIES
+ integer NTSTEP_BETWEEN_FRAMES,NTSTEP_BETWEEN_OUTPUT_INFO
+
+ character(len=150) LOCAL_PATH,clean_LOCAL_PATH,final_LOCAL_PATH,prname
+ ! parameters deduced from parameters read from file
+ integer NPROC,NEX_PER_PROC_XI,NEX_PER_PROC_ETA
+ integer NER
+
+ integer NSPEC_AB,NSPEC2D_A_XI,NSPEC2D_B_XI, &
+ NSPEC2D_A_ETA,NSPEC2D_B_ETA, &
+ NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX, &
+ NSPEC2D_BOTTOM,NSPEC2D_TOP, &
+ NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NGLOB_AB
+
+ double precision max_all_frames
+ ! ************** PROGRAM STARTS HERE **************
+
+
+ call getarg(1,movie_data_prefix)
+ call getarg(2,par_file)
+ call getarg(3,start_frame)
+ call getarg(4,end_frame)
+ if (trim(movie_data_prefix) == '' .or. trim(par_file) == '' &
+ .or. trim(start_frame) == '' .or. trim(end_frame) == '' ) &
+ stop 'Usage: xcreate_movie_GMT movie_data_prefix par_file start_frame end_frame output_file_prefix'
+
+ read(start_frame, *,iostat=ios1) it1
+ read(end_frame, *, iostat=ios2) it2
+ if (ios1 /= 0 .or. ios2 /= 0) stop 'Error reading start_frame and end_frame'
+
+! read the parameter file
+ call read_parameter_file(par_file,LATITUDE_MIN,LATITUDE_MAX,LONGITUDE_MIN,LONGITUDE_MAX, &
+ UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK, &
+ NER_SEDIM,NER_BASEMENT_SEDIM,NER_16_BASEMENT,NER_MOHO_16,NER_BOTTOM_MOHO, &
+ NEX_ETA,NEX_XI,NPROC_ETA,NPROC_XI,NTSTEP_BETWEEN_OUTPUT_SEISMOS,NSTEP,UTM_PROJECTION_ZONE,DT, &
+ ATTENUATION,USE_OLSEN_ATTENUATION,HARVARD_3D_GOCAD_MODEL,TOPOGRAPHY,LOCAL_PATH,NSOURCES, &
+ THICKNESS_TAPER_BLOCK_HR,THICKNESS_TAPER_BLOCK_MR,VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM, &
+ OCEANS,IMPOSE_MINIMUM_VP_GOCAD,HAUKSSON_REGIONAL_MODEL,ANISOTROPY, &
+ BASEMENT_MAP,MOHO_MAP_LUPEI,ABSORBING_CONDITIONS, &
+ MOVIE_SURFACE,MOVIE_VOLUME,CREATE_SHAKEMAP,SAVE_DISPLACEMENT, &
+ NTSTEP_BETWEEN_FRAMES,USE_HIGHRES_FOR_MOVIES, &
+ SAVE_AVS_DX_MESH_FILES,PRINT_SOURCE_TIME_FUNCTION,NTSTEP_BETWEEN_OUTPUT_INFO)
+
+! compute other parameters based upon values read
+ call compute_parameters(NER,NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA, &
+ NPROC,NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
+ NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM,NER_SEDIM, &
+ NSPEC_AB,NSPEC2D_A_XI,NSPEC2D_B_XI, &
+ NSPEC2D_A_ETA,NSPEC2D_B_ETA, &
+ NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
+ NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NGLOB_AB)
+
+ if(USE_HIGHRES_FOR_MOVIES) then
+ ilocnum = NGLLSQUARE*NEX_PER_PROC_XI*NEX_PER_PROC_ETA
+ else
+ ilocnum = NGNOD2D_AVS_DX*NEX_PER_PROC_XI*NEX_PER_PROC_ETA
+ endif
+ allocate(store_val(ilocnum,0:NPROC-1))
+ allocate(store_val_double(ilocnum,0:NPROC-1))
+ print *, 'data = ', ilocnum*NPROC*8/1024/1024,' MB'
+
+ if(it1 == -1) then
+ plot_shaking_map = .true.
+ nframes = 1
+ it1 = 1
+ inorm = it2
+ if(inorm < 1 .or. inorm > 3) stop 'incorrect value of inorm'
+ it2 = 1
+ else
+ plot_shaking_map = .false.
+ print *,'movie frames have been saved every ',NTSTEP_BETWEEN_FRAMES,' time steps'
+ print *
+ print *
+ print *,'looping from ',it1,' to ',it2,' frame'
+ ! count number of movie frames
+ nframes = 0
+ do it = it1,it2
+ nframes = nframes + 1
+ enddo
+ print *
+ print *,'total number of frames will be ',nframes
+ if(nframes == 0) stop 'null number of frames'
+ max_all_frames = -100.0
+ endif
+
+ ! define the total number of elements at the surface
+ if(USE_HIGHRES_FOR_MOVIES) then
+ nspectot_AVS_max = NEX_XI * NEX_ETA * (NGLLX-1) * (NGLLY-1)
+ else
+ nspectot_AVS_max = NEX_XI * NEX_ETA
+ endif
+
+
+ ! loop on all the time steps in the range entered
+ do it = it1,it2
+
+ ! check if time step corresponds to a movie frame
+ iframe = iframe + 1
+ ivalue = it * NTSTEP_BETWEEN_FRAMES
+ print *, 'ivalue = ' ,ivalue
+ if(plot_shaking_map) then
+ print *,'reading shaking map snapshot'
+ else
+ print *,'reading snapshot frame ',it,' out of ',NSTEP/NTSTEP_BETWEEN_FRAMES
+ endif
+ print *
+
+ ! read all the elements from the same file
+ if(plot_shaking_map) then
+ write(outputname,"('/shakingdata')")
+ else
+ write(outputname,"('/moviedata',i6.6)") ivalue
+ endif
+
+ outputname = trim(movie_data_prefix) // trim(outputname)
+ outputname1 = trim(outputname)//'.new'
+ open(unit=IIN,file=outputname,status='old',form='unformatted',iostat=ios1)
+ open(unit=IOUT,file=outputname1,form='unformatted',iostat=ios2)
+
+ if (ios1 /= 0 .or. ios2 /= 0) then
+ print *, 'Error opening file ', trim(outputname), ' and ', trim(outputname1)
+ stop
+ else
+ print *, 'reading file ', trim(outputname)
+ print *, 'writing file ', trim(outputname1)
+ endif
+
+ do i = 1, 6
+ read(IIN,iostat=ios1) store_val_double
+ if (ios1 /= 0) stop 'Error reading store_val_double file'
+ store_val(:,:) = sngl(store_val_double(:,:))
+ write(IOUT,iostat=ios1) store_val
+ if (ios1 /= 0) stop 'Error writing store_val file'
+ enddo
+
+ close(IIN)
+ close(IOUT)
+
+ end do
+
+
+
+ end program convert_movie_real_to_double
Added: seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/convert_nonlinear.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/convert_nonlinear.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/convert_nonlinear.f90 2010-10-01 15:00:21 UTC (rev 17232)
@@ -0,0 +1,68 @@
+program convert_nonlinear
+
+ ! reads in the output of shakemaps
+ ! and convert it to non-linear scale
+
+ implicit none
+
+ character(len=150) :: infile,outfile,ch_power,ch_zmax
+ real*8 :: power,zmax,zmin
+ logical :: calculate_zmax
+ integer :: ios,i,j
+ real*8 :: lon,lat,zv
+
+
+
+ call getarg(1,infile)
+ call getarg(2,outfile)
+ call getarg(3,ch_power)
+ call getarg(4,ch_zmax)
+
+ if (trim(infile) == '' .or. trim(outfile) == '' ) then
+ stop 'Usage: convert_nonlinear infile outfile power [zmax]'
+ endif
+ if (trim(ch_power) == '') then
+ power = 0.25
+ else
+ read(ch_power,*) power
+ endif
+ if (trim(ch_zmax) == '') then
+ calculate_zmax = .true.
+ else
+ calculate_zmax = .false.
+ read(ch_zmax,*) zmax
+ endif
+
+
+ open(11,file = infile,status = 'old',iostat = ios)
+ if (ios /= 0) stop 'Error open input xyz file'
+ open(12,file = outfile,status='unknown',iostat = ios)
+ if (ios /= 0) stop 'Error open output xyz file'
+ i = 0
+ if (calculate_zmax) then
+ do while (ios == 0)
+ read(11,*,iostat=ios) lon,lat,zv
+ if (ios /= 0) exit
+ i = i + 1
+ if (i==1) then
+ zmax = zv
+ else if (zmax < zv) then
+ zmax = zv
+ endif
+ enddo
+ print *, 'The maximum value of the dataset is ', sngl(zmax)
+ close(11)
+ open(11,file = infile,status = 'old',iostat = ios)
+ if (ios /= 0) stop 'Error open input xyz file'
+ endif
+
+ do while (ios == 0)
+ read(11,*,iostat=ios) lon,lat,zv
+ if (ios /= 0) exit
+ write(12,'(f18.12,f16.12,e19.12)',iostat=ios) lon,lat,(zv/zmax)**power
+ if (ios /= 0) stop 'Error writing'
+ enddo
+ close(12)
+ close(11)
+
+end program convert_nonlinear
Added: seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/convert_tiff_topo.pl
===================================================================
--- seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/convert_tiff_topo.pl (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/convert_tiff_topo.pl 2010-10-01 15:00:21 UTC (rev 17232)
@@ -0,0 +1,49 @@
+#!/usr/bin/perl -w
+use lib '/opt/seismo-util/lib/perl';
+use CMT_TOOLS;
+
+ at frames = (6..6);
+
+$dx = 0.01458333333333;
+$dy = 0.01197916666667;
+$R = "-R-120.3/-114.7/32.2/36.8";
+$JM = "-JM5i";
+$B = " -B2/2wesn ";
+$null = 127.5000;
+
+$datalib = "/opt/seismo-util/data/datalib_SC";
+$int_file = "$datalib/topo_cal.int";
+$fault_file = "$datalib/jennings.xy";
+$cmt_file = "CMTSOLUTION";
+$cpt_file = "disp.cpt";
+
+($elat,$elon) = get_cmt_location($cmt_file);
+
+open(CSH,">convert_tiff.csh");
+
+print CSH "gmtset BASEMAP_TYPE plain ANOT_FONT_SIZE 9 HEADER_FONT_SIZE 15\n";
+print CSH "makecpt -Cpolar -T0/255/2 -Z -V > $cpt_file\n";
+foreach $frame (@frames) {
+ if ($frame < 10) {$frame = "0$frame";}
+ $file = "gmt_movie_0000$frame";
+ $time = $frame * 200 * 0.009;
+ $xyz_file = "$file.xyz";
+ $grd_file = "$file.grd";
+ $sample_grd_file = "$grd_file.sample";
+ $cut_int_file = "$file.int";
+ $ps_file = "$file.ps";
+ $tiff_file = "$file.tiff";
+ print CSH "\n# frame $frame\n";
+ print CSH "awk '\$3 != $null {print \$0}' $xyz_file | xyz2grd -I$dx/$dy $R -G$grd_file -N128 -V \n";
+ print CSH "grdsample $grd_file -G$sample_grd_file -I0.00833333/0.00833333 -F -V\n";
+ print CSH "grdcut $int_file -G$cut_int_file $R -V \n";
+ print CSH "grdimage $sample_grd_file $JM $R -I$cut_int_file -C$cpt_file $B -K -V -P > $ps_file\n";
+ print CSH "pscoast $JM $R -W4 -Na -Dh -K -O -P -V >> $ps_file \n";
+ print CSH "psxy $JM $R -M -W2 -K -O -P -V $fault_file >> $ps_file\n";
+ print CSH "psxy $JM $R -Sa0.10 -W1 -G0/0/0 -K -O -P -V <<EOF >> $ps_file\n$elon $elat\nEOF\n";
+ print CSH "pstext $JM $R -N -W255/255/255 -O -P -V <<EOF >>$ps_file \n";
+ print CSH "-114.8 36.5 12 0 0 RT time = $time s \nEOF\n";
+ print CSH "convert $ps_file $tiff_file\n";
+}
+close(CSH);
+
Property changes on: seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/convert_tiff_topo.pl
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/create_movie_GMT.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/create_movie_GMT.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/create_movie_GMT.f90 2010-10-01 15:00:21 UTC (rev 17232)
@@ -0,0 +1,746 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 3 . 3
+! --------------------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology September 2002
+!
+! 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.
+!
+!=====================================================================
+
+!
+!--- create a movie of vertical component of surface displacement or velocity
+!--- in AVS or OpenDX format
+!
+
+program create_movie_GMT
+
+ implicit none
+
+ include 'constants.h'
+
+ character(len=200) par_file, movie_data_prefix, start_frame, end_frame, &
+ output_file_prefix
+ integer ios1, ios2, nspectot_AVS_max
+ ! threshold in percent of the maximum below which we cut the amplitude
+ logical, parameter :: APPLY_THRESHOLD = .false.
+ real(kind=CUSTOM_REAL), parameter :: THRESHOLD = 1._CUSTOM_REAL / 100._CUSTOM_REAL
+
+ ! coefficient of power law used for non linear scaling
+ logical, parameter :: NONLINEAR_SCALING = .false.
+ real(kind=CUSTOM_REAL), parameter :: POWER_SCALING = 0.50_CUSTOM_REAL
+
+ integer it,it1,it2,ivalue,ispec
+ integer iformat,nframes,iframe,inumber,inorm,iscaling_shake
+ integer ibool_number,ibool_number1,ibool_number2,ibool_number3,ibool_number4
+
+ logical UNIQUE_FILE,plot_shaking_map
+
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: x,y,z,display
+ real(kind=CUSTOM_REAL) xcoord,ycoord,zcoord
+ real(kind=CUSTOM_REAL) vectorx,vectory,vectorz
+
+ double precision min_field_current,max_field_current,max_absol
+
+ character(len=150) outputname
+
+ integer iproc,ipoin
+
+ ! GMT
+ double precision lat,long,zscaling
+ integer igmt
+
+ ! for sorting routine
+ integer npointot,ilocnum,nglob,i,j,ielm,ieoff,ispecloc
+ integer, dimension(:), allocatable :: iglob,loc,ireorder
+ logical, dimension(:), allocatable :: ifseg,mask_point
+ double precision, dimension(:), allocatable :: xp,yp,zp,xp_save,yp_save,zp_save,field_display
+
+ ! movie files stored by solver
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: &
+ store_val_x,store_val_y,store_val_z, &
+ store_val_ux,store_val_uy,store_val_uz
+
+! parameters read from parameter file
+ integer NER_SEDIM,NER_BASEMENT_SEDIM,NER_16_BASEMENT,NER_MOHO_16,NER_BOTTOM_MOHO, &
+ NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA,NTSTEP_BETWEEN_OUTPUT_SEISMOS,NSTEP,UTM_PROJECTION_ZONE,SIMULATION_TYPE
+ integer NSOURCES,NTSTEP_BETWEEN_FRAMES,NTSTEP_BETWEEN_OUTPUT_INFO
+
+ double precision UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK, UTM_MAX
+ double precision LATITUDE_MIN,LATITUDE_MAX,LONGITUDE_MIN,LONGITUDE_MAX,DT,HDUR_MOVIE
+ double precision THICKNESS_TAPER_BLOCK_HR,THICKNESS_TAPER_BLOCK_MR,VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM
+
+ logical HARVARD_3D_GOCAD_MODEL,TOPOGRAPHY,ATTENUATION,USE_OLSEN_ATTENUATION, &
+ OCEANS,IMPOSE_MINIMUM_VP_GOCAD,HAUKSSON_REGIONAL_MODEL, &
+ BASEMENT_MAP,MOHO_MAP_LUPEI,ABSORBING_CONDITIONS,SAVE_FORWARD
+ logical MOVIE_SURFACE,MOVIE_VOLUME,CREATE_SHAKEMAP,SAVE_DISPLACEMENT,USE_HIGHRES_FOR_MOVIES
+ logical ANISOTROPY,SAVE_MESH_FILES,PRINT_SOURCE_TIME_FUNCTION,SUPPRESS_UTM_PROJECTION,USE_REGULAR_MESH
+
+ character(len=150) LOCAL_PATH,MODEL,CMTSOLUTION
+
+
+ ! parameters deduced from parameters read from file
+ integer NPROC,NEX_PER_PROC_XI,NEX_PER_PROC_ETA
+ integer NER
+
+ integer NSPEC_AB,NSPEC2D_A_XI,NSPEC2D_B_XI, &
+ NSPEC2D_A_ETA,NSPEC2D_B_ETA, &
+ NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX, &
+ NSPEC2D_BOTTOM,NSPEC2D_TOP, &
+ NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NGLOB_AB
+
+ double precision max_all_frames
+ ! ************** PROGRAM STARTS HERE **************
+
+
+ call getarg(1,movie_data_prefix)
+ call getarg(2,par_file)
+ call getarg(3,start_frame)
+ call getarg(4,end_frame)
+ call getarg(5,output_file_prefix)
+ if (trim(movie_data_prefix) == '' .or. trim(par_file) == '' &
+ .or. trim(start_frame) == '' .or. trim(end_frame) == '' &
+ .or. trim(output_file_prefix) == '') &
+ stop 'Usage: xcreate_movie_GMT movie_data_prefix par_file start_frame end_frame output_file_prefix'
+
+ read(start_frame, *,iostat=ios1) it1
+ read(end_frame, *, iostat=ios2) it2
+ if (ios1 /= 0 .or. ios2 /= 0) stop 'Error reading start_frame and end_frame'
+
+! read the parameter file
+ call read_parameter_file(par_file,LATITUDE_MIN,LATITUDE_MAX,LONGITUDE_MIN,LONGITUDE_MAX, &
+ UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK, &
+ NER_SEDIM,NER_BASEMENT_SEDIM,NER_16_BASEMENT,NER_MOHO_16,NER_BOTTOM_MOHO, &
+ NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA,NTSTEP_BETWEEN_OUTPUT_SEISMOS,NSTEP,UTM_PROJECTION_ZONE,DT, &
+ ATTENUATION,USE_OLSEN_ATTENUATION,HARVARD_3D_GOCAD_MODEL,TOPOGRAPHY,LOCAL_PATH,NSOURCES, &
+ THICKNESS_TAPER_BLOCK_HR,THICKNESS_TAPER_BLOCK_MR,VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM, &
+ OCEANS,IMPOSE_MINIMUM_VP_GOCAD,HAUKSSON_REGIONAL_MODEL,ANISOTROPY, &
+ BASEMENT_MAP,MOHO_MAP_LUPEI,ABSORBING_CONDITIONS, &
+ MOVIE_SURFACE,MOVIE_VOLUME,CREATE_SHAKEMAP,SAVE_DISPLACEMENT, &
+ NTSTEP_BETWEEN_FRAMES,USE_HIGHRES_FOR_MOVIES,HDUR_MOVIE, &
+ SAVE_MESH_FILES,PRINT_SOURCE_TIME_FUNCTION,NTSTEP_BETWEEN_OUTPUT_INFO, &
+ SUPPRESS_UTM_PROJECTION,MODEL,USE_REGULAR_MESH,SIMULATION_TYPE,SAVE_FORWARD)
+
+! compute other parameters based upon values read
+ call compute_parameters(NER,NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA, &
+ NPROC,NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
+ NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM,NER_SEDIM, &
+ NSPEC_AB,NSPEC2D_A_XI,NSPEC2D_B_XI, &
+ NSPEC2D_A_ETA,NSPEC2D_B_ETA, &
+ NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
+ NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NGLOB_AB,USE_REGULAR_MESH)
+
+ print *
+ print *,'There are ',NPROC,' slices numbered from 0 to ',NPROC-1
+ print *
+
+ print *, 'DT = ', DT , ' NSTEP = ', NSTEP
+ print *
+
+ if(SAVE_DISPLACEMENT) then
+ print *,'Vertical displacement will be shown in movie'
+ else
+ print *,'Vertical velocity will be shown in movie'
+ endif
+ print *
+
+ if(USE_HIGHRES_FOR_MOVIES) then
+ print *, 'Movie is in high-resolution'
+ ilocnum = NGLLSQUARE*NEX_PER_PROC_XI*NEX_PER_PROC_ETA
+ else
+ print *, 'Movie is in low-resolution'
+ ilocnum = NGNOD2D_AVS_DX*NEX_PER_PROC_XI*NEX_PER_PROC_ETA
+ endif
+ allocate(store_val_x(ilocnum,0:NPROC-1))
+ allocate(store_val_y(ilocnum,0:NPROC-1))
+ allocate(store_val_z(ilocnum,0:NPROC-1))
+ allocate(store_val_ux(ilocnum,0:NPROC-1))
+ allocate(store_val_uy(ilocnum,0:NPROC-1))
+ allocate(store_val_uz(ilocnum,0:NPROC-1))
+
+ zscaling = 0.
+
+ if(it1 == -1) then
+ plot_shaking_map = .true.
+ nframes = 1
+ it1 = 1
+ inorm = it2
+ if(inorm < 1 .or. inorm > 3) stop 'incorrect value of inorm'
+ it2 = 1
+ else
+ plot_shaking_map = .false.
+ print *,'movie frames have been saved every ',NTSTEP_BETWEEN_FRAMES,' time steps'
+ print *
+ print *
+ print *,'looping from ',it1,' to ',it2,' frame'
+ ! count number of movie frames
+ nframes = 0
+ do it = it1,it2
+ nframes = nframes + 1
+ enddo
+ print *
+ print *,'total number of frames will be ',nframes
+ if(nframes == 0) stop 'null number of frames'
+ max_all_frames = -100.0
+ endif
+
+
+ ! define the total number of elements at the surface
+ if(USE_HIGHRES_FOR_MOVIES) then
+ nspectot_AVS_max = NEX_XI * NEX_ETA * (NGLLX-1) * (NGLLY-1)
+ else
+ nspectot_AVS_max = NEX_XI * NEX_ETA
+ endif
+
+ print *
+ print *,'there are a total of ',nspectot_AVS_max,' elements at the surface'
+ print *
+
+ ! maximum theoretical number of points at the surface
+ npointot = NGNOD2D_AVS_DX * nspectot_AVS_max
+
+ ! allocate arrays for sorting routine
+ allocate(iglob(npointot),loc(npointot))
+ allocate(ifseg(npointot))
+ allocate(xp(npointot),yp(npointot),zp(npointot))
+ allocate(xp_save(npointot),yp_save(npointot),zp_save(npointot))
+ allocate(field_display(npointot))
+ allocate(mask_point(npointot))
+ allocate(ireorder(npointot))
+
+ ! --------------------------------------
+
+ if(USE_HIGHRES_FOR_MOVIES) then
+ allocate(x(NGLLX,NGLLY))
+ allocate(y(NGLLX,NGLLY))
+ allocate(z(NGLLX,NGLLY))
+ allocate(display(NGLLX,NGLLY))
+ endif
+
+ iframe = 0
+
+ ! loop on all the time steps in the range entered
+ do it = it1,it2
+
+ ! check if time step corresponds to a movie frame
+ iframe = iframe + 1
+ ivalue = it * NTSTEP_BETWEEN_FRAMES
+! print *, 'ivalue = ' ,ivalue
+ if(plot_shaking_map) then
+ print *,'reading shaking map snapshot'
+ else
+ print *,'reading snapshot frame ',it,' out of ',NSTEP/NTSTEP_BETWEEN_FRAMES
+ endif
+ print *
+
+ ! read all the elements from the same file
+ if(plot_shaking_map) then
+ write(outputname,"('/shakingdata')")
+ else
+ write(outputname,"('/moviedata',i6.6)") ivalue
+ endif
+ outputname = trim(movie_data_prefix) // trim(outputname)
+ open(unit=IOUT,file=outputname,status='old',form='unformatted',iostat=ios1)
+ if (ios1 /= 0) then
+ print *, 'Error opening file ', trim(outputname)
+ stop
+ else
+ print *, 'reading file ', trim(outputname)
+ endif
+ read(IOUT) store_val_x
+ read(IOUT) store_val_y
+ read(IOUT) store_val_z
+ read(IOUT) store_val_ux
+ read(IOUT) store_val_uy
+ read(IOUT) store_val_uz
+ close(IOUT)
+
+ call utm_geo(long,lat,dble(minval(store_val_x)),dble(minval(store_val_y)), &
+ UTM_PROJECTION_ZONE,IUTM2LONGLAT,SUPPRESS_UTM_PROJECTION)
+! print *
+! print *, 'min long = ', long, ' min y = ', lat
+ call utm_geo(long,lat,dble(maxval(store_val_x)),dble(maxval(store_val_y)), &
+ UTM_PROJECTION_ZONE,IUTM2LONGLAT,SUPPRESS_UTM_PROJECTION)
+! print *
+! print *, 'max long = ', long, ' max lat = ', lat
+
+ ! clear number of elements kept
+ ispec = 0
+
+ ! read points for all the slices
+ do iproc = 0,NPROC-1
+
+ ! reset point number
+ ipoin = 0
+
+ do ispecloc = 1,NEX_PER_PROC_XI*NEX_PER_PROC_ETA
+
+ if(USE_HIGHRES_FOR_MOVIES) then
+ ! assign the OpenDX "elements"
+
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+
+ ipoin = ipoin + 1
+
+ xcoord = store_val_x(ipoin,iproc)
+ ycoord = store_val_y(ipoin,iproc)
+ zcoord = store_val_z(ipoin,iproc)
+
+ ! amplify topography, otherwise too flat to see anything
+ zcoord = zcoord * zscaling
+
+ vectorx = store_val_ux(ipoin,iproc)
+ vectory = store_val_uy(ipoin,iproc)
+ vectorz = store_val_uz(ipoin,iproc)
+
+ x(i,j) = xcoord
+ y(i,j) = ycoord
+ z(i,j) = zcoord
+
+ if(plot_shaking_map) then
+ if(inorm == 1) then
+ display(i,j) = vectorx
+ else if(inorm == 2) then
+ display(i,j) = vectory
+ else
+ display(i,j) = vectorz
+ endif
+ else
+ display(i,j) = vectorz
+ endif
+
+ enddo
+ enddo
+
+ ! assign the values of the corners of the OpenDX "elements"
+ ispec = ispec + 1
+ ielm = (NGLLX-1)*(NGLLY-1)*(ispec-1)
+
+ do j = 1,NGLLY-1
+ do i = 1,NGLLX-1
+ ieoff = NGNOD2D_AVS_DX*(ielm+(i-1)+(j-1)*(NGLLX-1))
+ do ilocnum = 1,NGNOD2D_AVS_DX
+
+ if(ilocnum == 1) then
+ xp(ieoff+ilocnum) = dble(x(i,j))
+ yp(ieoff+ilocnum) = dble(y(i,j))
+ zp(ieoff+ilocnum) = dble(z(i,j))
+ field_display(ieoff+ilocnum) = dble(display(i,j))
+ elseif(ilocnum == 2) then
+ xp(ieoff+ilocnum) = dble(x(i+1,j))
+ yp(ieoff+ilocnum) = dble(y(i+1,j))
+ zp(ieoff+ilocnum) = dble(z(i+1,j))
+ field_display(ieoff+ilocnum) = dble(display(i+1,j))
+ elseif(ilocnum == 3) then
+ xp(ieoff+ilocnum) = dble(x(i+1,j+1))
+ yp(ieoff+ilocnum) = dble(y(i+1,j+1))
+ zp(ieoff+ilocnum) = dble(z(i+1,j+1))
+ field_display(ieoff+ilocnum) = dble(display(i+1,j+1))
+ else
+ xp(ieoff+ilocnum) = dble(x(i,j+1))
+ yp(ieoff+ilocnum) = dble(y(i,j+1))
+ zp(ieoff+ilocnum) = dble(z(i,j+1))
+ field_display(ieoff+ilocnum) = dble(display(i,j+1))
+ endif
+
+ enddo
+ enddo
+ enddo
+
+ else
+
+ ispec = ispec + 1
+ ieoff = NGNOD2D_AVS_DX*(ispec-1)
+
+ ! four points for each element
+ do ilocnum = 1,NGNOD2D_AVS_DX
+
+ ipoin = ipoin + 1
+
+ xcoord = store_val_x(ipoin,iproc)
+ ycoord = store_val_y(ipoin,iproc)
+ zcoord = store_val_z(ipoin,iproc)
+
+ ! amplify topography, otherwise too flat to see anything
+ zcoord = zcoord * zscaling
+
+ vectorx = store_val_ux(ipoin,iproc)
+ vectory = store_val_uy(ipoin,iproc)
+ vectorz = store_val_uz(ipoin,iproc)
+
+ xp(ilocnum+ieoff) = dble(xcoord)
+ yp(ilocnum+ieoff) = dble(ycoord)
+ zp(ilocnum+ieoff) = dble(zcoord)
+
+ ! show vertical component of displacement or velocity in the movie
+ ! or show norm of vector if shaking map
+ ! for shaking map, norm of U stored in ux, V in uy and A in uz
+ if(plot_shaking_map) then
+ if(inorm == 1) then
+ field_display(ilocnum+ieoff) = dble(vectorx)
+ else if(inorm == 2) then
+ field_display(ilocnum+ieoff) = dble(vectory)
+ else
+ field_display(ilocnum+ieoff) = dble(vectorz)
+ endif
+ else
+ field_display(ilocnum+ieoff) = dble(vectorz) ! only plotting z component for now.
+ !field_display(ilocnum+ieoff) = dble(vectorx * 0.87 - vectory * 0.5) ! plot hroizontal movie
+ endif
+
+ enddo
+
+ endif
+
+ enddo
+ enddo
+
+ ! copy coordinate arrays since the sorting routine does not preserve them
+ xp_save(:) = xp(:)
+ yp_save(:) = yp(:)
+ zp_save(:) = zp(:)
+
+ !--- sort the list based upon coordinates to get rid of multiples
+ print *,'sorting list of points'
+ call get_global_AVS(nspectot_AVS_max,xp,yp,zp,iglob,loc,ifseg,nglob,npointot,UTM_X_MIN,UTM_X_MAX)
+
+ !--- print total number of points found
+ print *
+ print *,'found a total of ',nglob,' points'
+ print *,'initial number of points (with multiples) was ',npointot
+
+
+ !--- normalize and scale vector field
+
+ ! compute min and max of data value to normalize
+ min_field_current = minval(field_display(:))
+ max_field_current = maxval(field_display(:))
+
+
+ if (max_field_current > max_all_frames) max_all_frames = max_field_current
+ ! print minimum and maximum amplitude in current snapshot
+ print *
+ print *,'minimum amplitude in current snapshot = ',min_field_current
+ print *,'maximum amplitude in current snapshot = ',max_field_current
+ print *
+
+ if(plot_shaking_map) then
+
+ ! normalize field to [0:1]
+ field_display(:) = field_display(:) / max_field_current
+
+ ! apply non linear scaling to normalized field
+ field_display = field_display ** POWER_SCALING
+
+ ! map field to [0:255] for AVS color scale
+ field_display(:) = 255. * field_display(:)
+
+ endif
+
+ !--- ****** create AVS file using sorted list ******
+
+
+ ! create file name and open file
+ if(plot_shaking_map) then
+ write(outputname,"('/gmt_shaking_map.xyz')")
+ open(unit=11,file=trim(output_file_prefix)//outputname,status='unknown',iostat=ios1)
+ else
+ write(outputname,"('/gmt_movie_',i6.6,'.xyz')") it
+ outputname = trim(output_file_prefix) //trim(outputname)
+ call system('\\rm -f ' // trim(outputname))
+ open(unit=11,file=outputname,status='unknown')!,form='unformatted',access='direct',recl = 3 * 8,iostat=ios1)
+ endif
+ if (ios1 /= 0) stop 'Error opening output file'
+
+ igmt = 1
+ ! output list of points
+ mask_point = .false.
+ do ispec=1,nspectot_AVS_max
+ ieoff = NGNOD2D_AVS_DX*(ispec-1)
+ ! four points for each element
+ do ilocnum = 1,NGNOD2D_AVS_DX
+ ibool_number = iglob(ilocnum+ieoff)
+ if(.not. mask_point(ibool_number)) then
+ call utm_geo(long,lat,xp_save(ilocnum+ieoff),yp_save(ilocnum+ieoff), &
+ UTM_PROJECTION_ZONE,IUTM2LONGLAT,SUPPRESS_UTM_PROJECTION)
+ if (plot_shaking_map) then
+ write(11,*) long,lat,field_display(ilocnum+ieoff)
+ else
+ !write(11,rec=igmt) long,lat,field_display(ilocnum+ieoff)
+ write (11,'(3e17.6)') long, lat, field_display(ilocnum+ieoff)
+ igmt = igmt + 1
+ endif
+ endif
+ mask_point(ibool_number) = .true.
+ enddo
+ enddo
+ print *, 'total number of record is ', igmt - 1
+ print *, ' '
+ print *, ' '
+
+ enddo
+
+ ! step number for AVS multistep file
+401 format('step',i1,' image',i1)
+402 format('step',i2,' image',i2)
+403 format('step',i3,' image',i3)
+404 format('step',i4,' image',i4)
+
+ print *,'GMT files are stored in OUTPUT_FILES/gmt_*.xyz'
+ print *
+ print *, 'The maximum absolute value of all frames is ', max_all_frames
+ print *
+ print *
+
+
+ deallocate(store_val_x)
+ deallocate(store_val_y)
+ deallocate(store_val_z)
+ deallocate(store_val_ux)
+ deallocate(store_val_uy)
+ deallocate(store_val_uz)
+
+ ! deallocate arrays for sorting routine
+ deallocate(iglob,loc)
+ deallocate(ifseg)
+ deallocate(xp,yp,zp)
+ deallocate(xp_save,yp_save,zp_save)
+ deallocate(field_display)
+ deallocate(mask_point)
+ deallocate(ireorder)
+
+ if(USE_HIGHRES_FOR_MOVIES) then
+ deallocate(x)
+ deallocate(y)
+ deallocate(z)
+ deallocate(display)
+ endif
+
+
+
+
+end program create_movie_GMT
+
+!
+!=====================================================================
+!
+
+ subroutine get_global_AVS(nspec,xp,yp,zp,iglob,loc,ifseg,nglob,npointot,UTM_X_MIN,UTM_X_MAX)
+
+! this routine MUST be in double precision to avoid sensitivity
+! to roundoff errors in the coordinates of the points
+
+! leave sorting subroutines in same source file to allow for inlining
+
+ implicit none
+
+ include "constants.h"
+
+! geometry tolerance parameter to calculate number of independent grid points
+! small value for double precision and to avoid sensitivity to roundoff
+ double precision SMALLVALTOL
+
+ integer npointot
+ integer iglob(npointot),loc(npointot)
+ logical ifseg(npointot)
+ double precision xp(npointot),yp(npointot),zp(npointot)
+ integer nspec,nglob
+
+ integer ispec,i,j
+ integer ieoff,ilocnum,nseg,ioff,iseg,ig
+
+ integer, dimension(:), allocatable :: ind,ninseg,iwork
+ double precision, dimension(:), allocatable :: work
+
+ double precision UTM_X_MIN,UTM_X_MAX
+
+! define geometrical tolerance based upon typical size of the model
+ SMALLVALTOL = 1.d-10 * dabs(UTM_X_MAX - UTM_X_MIN)
+
+! dynamically allocate arrays
+ allocate(ind(npointot))
+ allocate(ninseg(npointot))
+ allocate(iwork(npointot))
+ allocate(work(npointot))
+
+! establish initial pointers
+ do ispec=1,nspec
+ ieoff=NGNOD2D_AVS_DX*(ispec-1)
+ do ilocnum=1,NGNOD2D_AVS_DX
+ loc(ilocnum+ieoff)=ilocnum+ieoff
+ enddo
+ enddo
+
+ ifseg(:)=.false.
+
+ nseg=1
+ ifseg(1)=.true.
+ ninseg(1)=npointot
+
+ do j=1,NDIM
+
+! sort within each segment
+ ioff=1
+ do iseg=1,nseg
+ if(j == 1) then
+ call rank(xp(ioff),ind,ninseg(iseg))
+ else if(j == 2) then
+ call rank(yp(ioff),ind,ninseg(iseg))
+ else
+ call rank(zp(ioff),ind,ninseg(iseg))
+ endif
+ call swap_all(loc(ioff),xp(ioff),yp(ioff),zp(ioff),iwork,work,ind,ninseg(iseg))
+ ioff=ioff+ninseg(iseg)
+ enddo
+
+! check for jumps in current coordinate
+! compare the coordinates of the points within a small tolerance
+ if(j == 1) then
+ do i=2,npointot
+ if(dabs(xp(i)-xp(i-1)) > SMALLVALTOL) ifseg(i)=.true.
+ enddo
+ else if(j == 2) then
+ do i=2,npointot
+ if(dabs(yp(i)-yp(i-1)) > SMALLVALTOL) ifseg(i)=.true.
+ enddo
+ else
+ do i=2,npointot
+ if(dabs(zp(i)-zp(i-1)) > SMALLVALTOL) ifseg(i)=.true.
+ enddo
+ endif
+
+! count up number of different segments
+ nseg=0
+ do i=1,npointot
+ if(ifseg(i)) then
+ nseg=nseg+1
+ ninseg(nseg)=1
+ else
+ ninseg(nseg)=ninseg(nseg)+1
+ endif
+ enddo
+ enddo
+
+! assign global node numbers (now sorted lexicographically)
+ ig=0
+ do i=1,npointot
+ if(ifseg(i)) ig=ig+1
+ iglob(loc(i))=ig
+ enddo
+
+ nglob=ig
+
+! deallocate arrays
+ deallocate(ind)
+ deallocate(ninseg)
+ deallocate(iwork)
+ deallocate(work)
+
+ end subroutine get_global_AVS
+
+! -----------------------------------
+
+! sorting routines put in same file to allow for inlining
+
+ subroutine rank(A,IND,N)
+!
+! Use Heap Sort (Numerical Recipes)
+!
+ implicit none
+
+ integer n
+ double precision A(n)
+ integer IND(n)
+
+ integer i,j,l,ir,indx
+ double precision q
+
+ do j=1,n
+ IND(j)=j
+ enddo
+
+ if (n == 1) return
+
+ L=n/2+1
+ ir=n
+ 100 CONTINUE
+ IF (l>1) THEN
+ l=l-1
+ indx=ind(l)
+ q=a(indx)
+ ELSE
+ indx=ind(ir)
+ q=a(indx)
+ ind(ir)=ind(1)
+ ir=ir-1
+ if (ir == 1) then
+ ind(1)=indx
+ return
+ endif
+ ENDIF
+ i=l
+ j=l+l
+ 200 CONTINUE
+ IF (J <= IR) THEN
+ IF (J<IR) THEN
+ IF ( A(IND(j))<A(IND(j+1)) ) j=j+1
+ ENDIF
+ IF (q<A(IND(j))) THEN
+ IND(I)=IND(J)
+ I=J
+ J=J+J
+ ELSE
+ J=IR+1
+ ENDIF
+ goto 200
+ ENDIF
+ IND(I)=INDX
+ goto 100
+ end subroutine rank
+
+! ------------------------------------------------------------------
+
+ subroutine swap_all(IA,A,B,C,IW,W,ind,n)
+!
+! swap arrays IA, A, B and C according to addressing in array IND
+!
+ implicit none
+
+ integer n
+
+ integer IND(n)
+ integer IA(n),IW(n)
+ double precision A(n),B(n),C(n),W(n)
+
+ integer i
+
+ IW(:) = IA(:)
+ W(:) = A(:)
+
+ do i=1,n
+ IA(i)=IW(ind(i))
+ A(i)=W(ind(i))
+ enddo
+
+ W(:) = B(:)
+
+ do i=1,n
+ B(i)=W(ind(i))
+ enddo
+
+ W(:) = C(:)
+
+ do i=1,n
+ C(i)=W(ind(i))
+ enddo
+
+ end subroutine swap_all
+
Added: seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/get_value_parameters.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/get_value_parameters.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/get_value_parameters.f90 2010-10-01 15:00:21 UTC (rev 17232)
@@ -0,0 +1,82 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology September 2006
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+ subroutine get_value_integer(value_to_get, name, default_value)
+
+ implicit none
+
+ integer value_to_get, default_value
+ character(len=*) name
+
+ call unused_string(name)
+
+ value_to_get = default_value
+
+ end subroutine get_value_integer
+
+!--------------------
+
+ subroutine get_value_double_precision(value_to_get, name, default_value)
+
+ implicit none
+
+ double precision value_to_get, default_value
+ character(len=*) name
+
+ call unused_string(name)
+
+ value_to_get = default_value
+
+ end subroutine get_value_double_precision
+
+!--------------------
+
+ subroutine get_value_logical(value_to_get, name, default_value)
+
+ implicit none
+
+ logical value_to_get, default_value
+ character(len=*) name
+
+ call unused_string(name)
+
+ value_to_get = default_value
+
+ end subroutine get_value_logical
+
+!--------------------
+
+ subroutine get_value_string(value_to_get, name, default_value)
+
+ implicit none
+
+ character(len=*) value_to_get, default_value
+ character(len=*) name
+
+ call unused_string(name)
+
+ value_to_get = default_value
+
+ end subroutine get_value_string
Added: seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/makefile
===================================================================
--- seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/makefile (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/makefile 2010-10-01 15:00:21 UTC (rev 17232)
@@ -0,0 +1,28 @@
+F90 = gfortran
+FLAGS = #-O3
+LIB =
+OBJ = obj/
+BIN = .
+
+SRC_FILES = create_movie_GMT \
+ read_parameter_file compute_parameters \
+ utm_geo read_value_parameters get_value_parameters
+
+OBJ_FILES = $(patsubst %, $(OBJ)%.o, $(SRC_FILES))
+
+all : movie convert_nonlinear
+
+movie : $(OBJ_FILES)
+ $(F90) $(FLAGS) -o $(BIN)/xcreate_movie_GMT $(OBJ_FILES) $(LIB)
+
+$(OBJ)%.o : %.f90
+ $(F90) $(FLAGS) -o $(OBJ)$*.o -c $*.f90
+
+convert_nonlinear: convert_nonlinear.f90
+ $(F90) $(FLAGS) -o $(BIN)/xconvert_nonlinear convert_nonlinear.f90 $(LIB)
+
+all : movie convert_nonlinear
+
+clean:
+ \rm -f $(OBJ)*.o *~ xcreate_movie_GMT xconvert_nonlinear
+
Added: seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/movie2gif.pl
===================================================================
--- seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/movie2gif.pl (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/movie2gif.pl 2010-10-01 15:00:21 UTC (rev 17232)
@@ -0,0 +1,198 @@
+#!/usr/bin/perl -w
+
+use Getopt::Std;
+use lib '/opt/seismo-util/lib/perl';
+use CMT_TOOLS;
+use GMT_PLOT;
+use vars qw($opt_t $opt_p $opt_n $opt_g $opt_x $opt_2 $opt_s $opt_d $opt_R);
+sub Usage {
+ print STDERR <<EOF;
+
+ This script converts moviedata to xyz files, and then plot them in GMT
+ Usage : movie2gif.pl -R -m cmt_file -g -f start/end -p -2 -n -d dist(km) -s sta_file -x
+ -m cmt file
+ -R min_lon/max_lon/min_lat/max_lat (default -120.3/-114.7/32.2/36.8)
+ -g convert the movie data to gmt xyz files, you can choose to skip
+ this step if you have already done it.
+ -f start and end frame
+ -p add topography to the plot
+ -2 2d plot, otherwise 3d plot
+ -s station name file to plot
+ -n -- run nearneighbor command to interpolate xyz file into grd file,
+ since this is the step that takes most of the time, you can choose
+ to skip this step if you have already run it.
+ -d -- distance to average for nearneighbor command option -S (default 5 km)
+ -x -- executing the c shell script
+
+EOF
+exit(1);
+}
+#$bindir = "/opt/seismo-util/source/basin_inversion/bin";
+
+
+if (@ARGV == 0) {Usage();}
+if (not getopts('m:f:pg2xns:R:d:')) {die("Check options\n");}
+if ($opt_m) {$cmt_file = $opt_m;
+ ($elat,$elon) = get_cmt_location($cmt_file);}
+
+if (not defined $opt_f) {die("give the start and end of frame\n");}
+($start,$end) = split(/\//,$opt_f);
+#$transparent = 1;
+if ($opt_p) {$no_topo = 0;} else {$no_topo = 1;}
+if ($opt_g) {$tran_to_gmt = 1;} else {$tran_to_gmt = 0;}
+if ($opt_2) {$two_d = 1; $E = "";} else {$two_d = 0; $E = "-E190/60";}
+if ($opt_n) {$near = 1;} else {$near = 0;}
+if ($opt_s) {
+ $sta_file = $opt_s;
+ system("filelist.pl -I $sta_file -f $sta_file -n sta.tmp");}
+
+if (not defined $opt_n) {$opt_n = 1;}
+
+if ($opt_R) {$R = "-R$opt_R";} else {$R = "-R-120.3/-114.7/32.2/36.8";} #-120.3
+$R2 = "${R}/0/0.8";
+$R3 = "${R}/-1000/4000";
+$JM = "-JM4.5i";
+$JM3 = "${JM} -JZ1.5i";
+$B = " -B2/2wesn ";
+$B2 = "-B2/2/0.5WSen";
+$B3 = "-B2/2/1000WSen";
+
+#$Itopo = "-I1m/1m";
+$I = "-I0.5m";
+if ($opt_d) {$S = "-S${opt_d}k";} else {$S = "-S5k";}
+
+$datalib = "/opt/seismo-util/data/datalib_SC";
+$fault_file = "$datalib/jennings.xy";
+$fault_file_3d = "$datalib/jennings.xyz";
+$cpt_file = "disp.cpt";
+$topo_cpt = "$datalib/topo_cal_color.cpt";
+$csh_file = "movie2gif.csh";
+$big_grd_file = "$datalib/big_sc.grd";
+
+$factor = 10000;
+
+
+# *******************
+# convert from movie data to gmt xyz files
+if ($tran_to_gmt) {
+ print "Transfer the movie data files to gmt xyz files \n";
+ $cmdline="xcreate_movie_GMT . Par_file $start $end . > movie.log";
+ print "running command:\n $cmdline ... \n";
+ system($cmdline);
+ if ($? != 0) {die("Error create movie xyz files\n"); }
+ (@junk) = split(" ", `grep '^ DT' movie.log`); $dt = $junk[2];
+ (@junk) = split(" ", `grep '^NTSTEP_BETWEEN_FRAMES' Par_file`); $nframe = $junk[2];
+} else {
+ if (not -f "movie.log") {die(" Check if movie log file exist or not\n");}
+ $dt = 0.009; $nframe = 200;
+}
+# figure out the maximum value
+(@junk) = split(" ", `grep 'maximum absolute value' movie.log`);
+$max = $junk[-1] / 1;
+print " The maximum value of all frames are : $max \n";
+
+# *********************
+open(CSH,">$csh_file");
+
+print CSH "gmtset BASEMAP_TYPE plain ANOT_FONT_SIZE 18 HEADER_FONT_SIZE 20\n";
+# make cpt file, manipulate the white color part
+if ($no_topo) {
+ system("makecpt -T-1/1/0.1 -Cpolar > temp.cpt");
+ open(SED,">sed.in");
+ print SED "/^N/s/^.*\$/N\t255\t255\t255/ \n";
+ print SED "/^F/s/^.*\$/F\t255\t0\t0/ \n";
+ print SED "/^B/s/^.*\$/B\t0\t0\t255/ \n";
+ print SED "/^-0.1/s/^.*\$/-0.1\t255\t255\t255\t0\t255\t255\t255/ \n";
+ print SED "/^5.55112e-17/s/^.*\$/0\t255\t255\t255\t0.1\t255\t255\t255/ \n";
+ close(SED);
+ system("sed -f sed.in temp.cpt > temp1.cpt\n");
+ convert_linear("temp1.cpt",$cpt_file,3);
+} else {
+ if (not -f $big_grd_file) {die("No such $big_grd_file \n");}
+ print CSH "grdcut $big_grd_file $R -Gtopo.grd \n";
+}
+
+foreach $frame ($start .. $end) {
+ if ($frame < 10) {$frame = sprintf("00%02d",$frame);}
+ if ($frame < 100) {$frame = sprintf("0%03d",$frame);}
+ if ($frame < 1000) {$frame = sprintf("%04d",$frame);}
+ $file = "gmt_movie_00$frame";
+ $time = $frame * $nframe * $dt;
+ $xyz_file = "$file.xyz"; if (not -f $xyz_file) {die("check $xyz_file\n");}
+ $ps_file = "$file.ps";
+ $gif_file = "$file.gif";
+# $tran_gif_file = "${file}_tran.gif";
+ $grd_file = "$file.grd";
+
+ print CSH "\n# frame $frame\n";
+ if ($near) {
+ print CSH "nearneighbor -F $R $I $S -G$grd_file $xyz_file -E0 \n";
+ print CSH "grdmath $grd_file $max DIV = field1.grd\n";}
+# else {
+## print CSH "awk '{print \$1,\$2,\$3/$max}' $xyz_file | pscontour $R $JM $B -A- -C$cpt_file -I -P> $ps_file\n";} next;
+
+ if ($no_topo) {
+ print CSH "grdgradient field1.grd -V -A225 -Nt -Gfield.int\n";
+ if ($two_d) { # grdimage for no topo
+ print CSH "grdimage $JM $R field1.grd -Ifield.int -C$cpt_file -K -P -V> $ps_file\n";
+ } else { # grdview for no topo
+ print CSH "psbasemap -P $E $R3 $B2 $JM3 -K > $ps_file\n";
+ print CSH "grdview field1.grd -Qs -P $E $R2 $JM3 -C$cpt_file -Ifield.int -K -O >> $ps_file\n"; }
+ } else {
+ print CSH "grdmath field1.grd $factor MUL topo.grd ADD = field2.grd \n";
+ print CSH "grdgradient field2.grd -V -A190/60 -Nt -Gfield.int\n";
+ if ($two_d) { # for grdimage topo
+ print CSH "grdimage field2.grd $JM $R -C$topo_cpt $B -Ifield.int -K -P > $ps_file\n";}
+ else { # grdview for topo
+ print CSH "psbasemap -P $E $R3 $B3 $JM3 -K > $ps_file\n";
+ print CSH "grdview field2.grd -Qs -P -Gtopo.grd $E $R3 $JM3 -C$topo_cpt -Ifield.int -K -O >> $ps_file \n";}
+ }
+# if ($no_topo) {
+ if ($two_d) {
+ print CSH "pscoast $JM $R -W4 -Na -Dh -K -O -P -V -S205/255/255 >> $ps_file \n";
+ print CSH "psxy $JM $R -M -W2 -K -O -P -V $fault_file >> $ps_file\n";
+ if ($opt_m) {print CSH "psxy $JM $R -Sa0.15 -W1 -G255/0/0 -K -O -P -V <<EOF >> $ps_file\n$elon $elat\nEOF\n";}
+ if ($opt_s) {
+ print CSH "awk '{print \$4, \$3}' sta.tmp | psxy $JM $R -St0.12 -W0.8 -G0/255/0 -K -O -P -V >> $ps_file\n";
+ print CSH "awk '{print \$4, \$3+0.1, 12, 0, 4, \"CM\", \$1}' sta.tmp | pstext $JM $R -G0/0/255 -N -P -K -O -V >> $ps_file \n";
+ }
+ print CSH "pstext $JM $R -N -W255/255/255 -O -P -V <<EOF >>$ps_file \n -114.8 36.5 12 0 0 RT time = $time s \nEOF\n";
+ } else {
+ print CSH "pscoast -JM -JZ -R -W5 -Na -Dh -K -O -P -V $E -S205/255/255 >> $ps_file \n";
+ print CSH "psclip -JM -JZ -R -K -O -P -Z0 >> $ps_file <<EOF\n-120.0 32.2\n-120.0 36.4 \n-114.7 36.4\n-114.7 32.2\n-120.0 32.2\nEOF\n";
+ print CSH "psxyz -JM -JZ -R -M -W2 -K -O -P -V $E $fault_file_3d >> $ps_file\n";
+ print CSH "psclip -C -K -O >> $ps_file\n";
+ if ($opt_m) {print CSH "psxyz -JM -JZ -R -Sa0.15 -W1 -G255/0/0 -K -O -P -V $E <<EOF >> $ps_file\n$elon $elat 0 \nEOF\n";}
+ if ($opt_s) {
+ print CSH "awk '{print \$4, \$3, 0}' sta.tmp | psxyz $E -JM -JZ -R -St0.12 -W0.8 -G0/255/0 -K -O -P -V >> $ps_file\n";
+ print CSH "awk '{print \$4, \$3+0.15, 15, 0, 4, \"CM\", \$1}' sta.tmp | pstext -JM -JZ -G0/0/255 $E -R -N -P -K -O -V -Z0>> $ps_file \n";
+ }
+ print CSH "pstext -JM -JZ -R -N -W255/255/255 -O -P -V -Z0<<EOF >>$ps_file \n -114.8 36.5 12 0 0 RT time = $time s \nEOF\n";
+ }
+#}
+ print CSH "convert $ps_file $gif_file\n";
+# if ($no_topo and $transparent) {print CSH "giftrans -t\\\#fefefe $gif_file > $tran_gif_file\n";}
+}
+
+close(CSH);
+if ($opt_x) {system("csh -fv $csh_file");}
+
+
+#---------------------------
+
+sub convert_linear {
+
+ my ($cpt,$new_cpt,$power) = @_;
+
+ open(CPT,"$cpt"); @old_cpt = <CPT>;close(CPT);
+ open(CPT2,">$new_cpt");
+ foreach $line (@old_cpt) {
+ if ($line =~/^\d+/) {
+ ($v1,$r1,$g1,$b1,$v2,$r2,$g2,$b2) = split(" ",$line);
+ $v1 = $v1 ** $power, $v2 = $v2 ** $power;
+ printf CPT2 ("%-15.10f\t%s\t%s\t%s\t%-15.10f\t%s\t%s\t%s\n",
+ $v1,$r1,$g1,$b1,$v2,$r2,$g2,$b2);
+ }else {print CPT2 "$line";}
+ }
+ close(CPT2);
+}
Property changes on: seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/movie2gif.pl
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/read_parameter_file.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/read_parameter_file.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/read_parameter_file.f90 2010-10-01 15:00:21 UTC (rev 17232)
@@ -0,0 +1,304 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology September 2006
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+ subroutine read_parameter_file(par_file,LATITUDE_MIN,LATITUDE_MAX,LONGITUDE_MIN,LONGITUDE_MAX, &
+ UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK, &
+ NER_SEDIM,NER_BASEMENT_SEDIM,NER_16_BASEMENT,NER_MOHO_16,NER_BOTTOM_MOHO, &
+ NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA,NTSTEP_BETWEEN_OUTPUT_SEISMOS,NSTEP,UTM_PROJECTION_ZONE,DT, &
+ ATTENUATION,USE_OLSEN_ATTENUATION,HARVARD_3D_GOCAD_MODEL,TOPOGRAPHY,LOCAL_PATH,NSOURCES, &
+ THICKNESS_TAPER_BLOCK_HR,THICKNESS_TAPER_BLOCK_MR,VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM, &
+ OCEANS,IMPOSE_MINIMUM_VP_GOCAD,HAUKSSON_REGIONAL_MODEL,ANISOTROPY, &
+ BASEMENT_MAP,MOHO_MAP_LUPEI,ABSORBING_CONDITIONS, &
+ MOVIE_SURFACE,MOVIE_VOLUME,CREATE_SHAKEMAP,SAVE_DISPLACEMENT, &
+ NTSTEP_BETWEEN_FRAMES,USE_HIGHRES_FOR_MOVIES,HDUR_MOVIE, &
+ SAVE_MESH_FILES,PRINT_SOURCE_TIME_FUNCTION,NTSTEP_BETWEEN_OUTPUT_INFO, &
+ SUPPRESS_UTM_PROJECTION,MODEL,USE_REGULAR_MESH,SIMULATION_TYPE,SAVE_FORWARD)
+
+ implicit none
+
+ include "constants.h"
+
+ character(len=*) :: par_file
+ integer NER_SEDIM,NER_BASEMENT_SEDIM,NER_16_BASEMENT,NER_MOHO_16,NER_BOTTOM_MOHO, &
+ NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA,NTSTEP_BETWEEN_OUTPUT_SEISMOS,NSTEP,UTM_PROJECTION_ZONE,SIMULATION_TYPE
+ integer NSOURCES,NTSTEP_BETWEEN_FRAMES,NTSTEP_BETWEEN_OUTPUT_INFO
+
+ double precision UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK, UTM_MAX
+ double precision LATITUDE_MIN,LATITUDE_MAX,LONGITUDE_MIN,LONGITUDE_MAX,DT,HDUR_MOVIE
+ double precision THICKNESS_TAPER_BLOCK_HR,THICKNESS_TAPER_BLOCK_MR,VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM
+
+ logical HARVARD_3D_GOCAD_MODEL,TOPOGRAPHY,ATTENUATION,USE_OLSEN_ATTENUATION, &
+ OCEANS,IMPOSE_MINIMUM_VP_GOCAD,HAUKSSON_REGIONAL_MODEL, &
+ BASEMENT_MAP,MOHO_MAP_LUPEI,ABSORBING_CONDITIONS,SAVE_FORWARD
+ logical MOVIE_SURFACE,MOVIE_VOLUME,CREATE_SHAKEMAP,SAVE_DISPLACEMENT,USE_HIGHRES_FOR_MOVIES
+ logical ANISOTROPY,SAVE_MESH_FILES,PRINT_SOURCE_TIME_FUNCTION,SUPPRESS_UTM_PROJECTION,USE_REGULAR_MESH
+
+ character(len=150) LOCAL_PATH,MODEL,CMTSOLUTION
+
+! local variables
+ integer ios,icounter,isource,idummy,NEX_MAX
+
+ double precision DEPTH_BLOCK_KM,RECORD_LENGTH_IN_SECONDS,hdur,minval_hdur
+
+ character(len=150) dummystring
+
+ integer, external :: err_occurred
+
+ open(unit=IIN,file=trim(par_file),status='old',action='read')
+
+ call read_value_integer(SIMULATION_TYPE, 'solver.SIMULATION_TYPE')
+ if(err_occurred() /= 0) return
+ call read_value_logical(SAVE_FORWARD, 'solver.SAVE_FORWARD')
+ if(err_occurred() /= 0) return
+
+ call read_value_double_precision(LATITUDE_MIN, 'mesher.LATITUDE_MIN')
+ if(err_occurred() /= 0) return
+ call read_value_double_precision(LATITUDE_MAX, 'mesher.LATITUDE_MAX')
+ if(err_occurred() /= 0) return
+ call read_value_double_precision(LONGITUDE_MIN, 'mesher.LONGITUDE_MIN')
+ if(err_occurred() /= 0) return
+ call read_value_double_precision(LONGITUDE_MAX, 'mesher.LONGITUDE_MAX')
+ if(err_occurred() /= 0) return
+ call read_value_double_precision(DEPTH_BLOCK_KM, 'mesher.DEPTH_BLOCK_KM')
+ if(err_occurred() /= 0) return
+ call read_value_integer(UTM_PROJECTION_ZONE, 'mesher.UTM_PROJECTION_ZONE')
+ if(err_occurred() /= 0) return
+ call read_value_logical(SUPPRESS_UTM_PROJECTION, 'mesher.SUPPRESS_UTM_PROJECTION')
+ if(err_occurred() /= 0) return
+
+ call read_value_integer(NEX_XI, 'mesher.NEX_XI')
+ if(err_occurred() /= 0) return
+ call read_value_integer(NEX_ETA, 'mesher.NEX_ETA')
+ if(err_occurred() /= 0) return
+ call read_value_integer(NPROC_XI, 'mesher.NPROC_XI')
+ if(err_occurred() /= 0) return
+ call read_value_integer(NPROC_ETA, 'mesher.NPROC_ETA')
+ if(err_occurred() /= 0) return
+
+! convert model size to UTM coordinates and depth of mesh to meters
+ call utm_geo(LONGITUDE_MIN,LATITUDE_MIN,UTM_X_MIN,UTM_Y_MIN,UTM_PROJECTION_ZONE,ILONGLAT2UTM,SUPPRESS_UTM_PROJECTION)
+ call utm_geo(LONGITUDE_MAX,LATITUDE_MAX,UTM_X_MAX,UTM_Y_MAX,UTM_PROJECTION_ZONE,ILONGLAT2UTM,SUPPRESS_UTM_PROJECTION)
+
+ Z_DEPTH_BLOCK = - dabs(DEPTH_BLOCK_KM) * 1000.d0
+
+! check that parameters computed are consistent
+ if(UTM_X_MIN >= UTM_X_MAX) stop 'horizontal dimension of UTM block incorrect'
+ if(UTM_Y_MIN >= UTM_Y_MAX) stop 'vertical dimension of UTM block incorrect'
+
+! set time step and radial distribution of elements
+! right distribution is determined based upon maximum value of NEX
+
+ NEX_MAX = max(NEX_XI,NEX_ETA)
+ UTM_MAX = max(UTM_Y_MAX-UTM_Y_MIN, UTM_X_MAX-UTM_X_MIN)/1000.0 ! in KM
+
+ call read_value_string(MODEL, 'model.name')
+ if(err_occurred() /= 0) return
+ call read_value_logical(OCEANS, 'model.OCEANS')
+ if(err_occurred() /= 0) return
+ call read_value_logical(TOPOGRAPHY, 'model.TOPOGRAPHY')
+ if(err_occurred() /= 0) return
+ call read_value_logical(ATTENUATION, 'model.ATTENUATION')
+ if(err_occurred() /= 0) return
+ call read_value_logical(USE_OLSEN_ATTENUATION, 'model.USE_OLSEN_ATTENUATION')
+ if(err_occurred() /= 0) return
+
+ if(dabs(DEPTH_BLOCK_KM) <= DEPTH_MOHO_SOCAL) &
+ stop 'bottom of mesh must be deeper than deepest regional layer for Southern California'
+
+! standard mesh for Southern California on Caltech cluster
+ if (UTM_MAX/NEX_MAX >= 1.5) then
+
+! time step in seconds
+ DT = 0.011d0
+
+! number of elements in the vertical direction
+ NER_SEDIM = 1
+ NER_BASEMENT_SEDIM = 2
+ NER_16_BASEMENT = 2
+ NER_MOHO_16 = 3
+ NER_BOTTOM_MOHO = 3
+
+ else if(UTM_MAX/NEX_MAX >= 1.0) then
+
+! time step in seconds
+ DT = 0.009d0
+
+! number of elements in the vertical direction
+ NER_SEDIM = 1
+ NER_BASEMENT_SEDIM = 2
+ NER_16_BASEMENT = 3
+ NER_MOHO_16 = 4
+ NER_BOTTOM_MOHO = 7
+
+ else
+ stop 'this value of NEX_MAX is not in the database, edit read_parameter_file.f90 and recompile'
+ endif
+
+! multiply parameters read by mesh scaling factor
+ NER_BASEMENT_SEDIM = NER_BASEMENT_SEDIM * 2
+ NER_16_BASEMENT = NER_16_BASEMENT * 2
+ NER_MOHO_16 = NER_MOHO_16 * 2
+ NER_BOTTOM_MOHO = NER_BOTTOM_MOHO * 4
+
+ if(MODEL == 'SoCal') then
+
+ BASEMENT_MAP = .false.
+ MOHO_MAP_LUPEI = .false.
+ HAUKSSON_REGIONAL_MODEL = .false.
+ HARVARD_3D_GOCAD_MODEL = .false.
+ THICKNESS_TAPER_BLOCK_HR = 1.d0
+ THICKNESS_TAPER_BLOCK_MR = 1.d0
+ IMPOSE_MINIMUM_VP_GOCAD = .false.
+ VP_MIN_GOCAD = 1.d0
+ VP_VS_RATIO_GOCAD_TOP = 2.0d0
+ VP_VS_RATIO_GOCAD_BOTTOM = 1.732d0
+ ANISOTROPY = .false.
+ USE_REGULAR_MESH = .false.
+
+ else if(MODEL == 'Harvard_LA') then
+
+ BASEMENT_MAP = .true.
+ MOHO_MAP_LUPEI = .true.
+ HAUKSSON_REGIONAL_MODEL = .true.
+ HARVARD_3D_GOCAD_MODEL = .true.
+ THICKNESS_TAPER_BLOCK_HR = 3000.d0
+ THICKNESS_TAPER_BLOCK_MR = 15000.d0
+ IMPOSE_MINIMUM_VP_GOCAD = .true.
+ VP_MIN_GOCAD = 750.d0
+ VP_VS_RATIO_GOCAD_TOP = 2.0d0
+ VP_VS_RATIO_GOCAD_BOTTOM = 1.732d0
+ ANISOTROPY = .false.
+ USE_REGULAR_MESH = .false.
+
+ else if(MODEL == 'Min_Chen_anisotropy') then
+
+ BASEMENT_MAP = .false.
+ MOHO_MAP_LUPEI = .false.
+ HAUKSSON_REGIONAL_MODEL = .false.
+ HARVARD_3D_GOCAD_MODEL = .false.
+ THICKNESS_TAPER_BLOCK_HR = 1.d0
+ THICKNESS_TAPER_BLOCK_MR = 1.d0
+ IMPOSE_MINIMUM_VP_GOCAD = .false.
+ VP_MIN_GOCAD = 1.d0
+ VP_VS_RATIO_GOCAD_TOP = 2.0d0
+ VP_VS_RATIO_GOCAD_BOTTOM = 1.732d0
+ ANISOTROPY = .true.
+ USE_REGULAR_MESH = .false.
+
+ else
+ stop 'model not implemented, edit read_parameter_file.f90 and recompile'
+ endif
+
+! check that Poisson's ratio is positive
+ if(VP_VS_RATIO_GOCAD_TOP <= sqrt(2.d0) .or. VP_VS_RATIO_GOCAD_BOTTOM <= sqrt(2.d0)) &
+ stop 'wrong value of Poisson''s ratio for Gocad Vs block'
+
+ call read_value_logical(ABSORBING_CONDITIONS, 'solver.ABSORBING_CONDITIONS')
+ if(err_occurred() /= 0) return
+ call read_value_double_precision(RECORD_LENGTH_IN_SECONDS, 'solver.RECORD_LENGTH_IN_SECONDS')
+ if(err_occurred() /= 0) return
+
+! compute total number of time steps, rounded to next multiple of 100
+ NSTEP = 100 * (int(RECORD_LENGTH_IN_SECONDS / (100.d0*DT)) + 1)
+ if ( NOISE_TOMOGRAPHY /= 0 ) NSTEP = 2*NSTEP-1 ! time steps needs to be doubled, due to +/- branches
+
+! compute the total number of sources in the CMTSOLUTION file
+! there are NLINES_PER_CMTSOLUTION_SOURCE lines per source in that file
+ call get_value_string(CMTSOLUTION, 'solver.CMTSOLUTION', 'DATA/CMTSOLUTION')
+ open(unit=1,file=CMTSOLUTION,iostat=ios,status='old',action='read')
+ if(ios /= 0) stop 'error opening CMTSOLUTION file'
+ icounter = 0
+ do while(ios == 0)
+ read(1,"(a)",iostat=ios) dummystring
+ if(ios == 0) icounter = icounter + 1
+ enddo
+ close(1)
+ if(mod(icounter,NLINES_PER_CMTSOLUTION_SOURCE) /= 0) &
+ stop 'total number of lines in CMTSOLUTION file should be a multiple of NLINES_PER_CMTSOLUTION_SOURCE'
+ NSOURCES = icounter / NLINES_PER_CMTSOLUTION_SOURCE
+ if(NSOURCES < 1) stop 'need at least one source in CMTSOLUTION file'
+
+ call read_value_logical(MOVIE_SURFACE, 'solver.MOVIE_SURFACE')
+ if(err_occurred() /= 0) return
+ call read_value_logical(MOVIE_VOLUME, 'solver.MOVIE_VOLUME')
+ if(err_occurred() /= 0) return
+ call read_value_integer(NTSTEP_BETWEEN_FRAMES, 'solver.NTSTEP_BETWEEN_FRAMES')
+ if(err_occurred() /= 0) return
+ call read_value_logical(CREATE_SHAKEMAP, 'solver.CREATE_SHAKEMAP')
+ if(err_occurred() /= 0) return
+ call read_value_logical(SAVE_DISPLACEMENT, 'solver.SAVE_DISPLACEMENT')
+ if(err_occurred() /= 0) return
+ call read_value_logical(USE_HIGHRES_FOR_MOVIES, 'solver.USE_HIGHRES_FOR_MOVIES')
+ if(err_occurred() /= 0) return
+ call read_value_double_precision(HDUR_MOVIE, 'solver.HDUR_MOVIE')
+ if(err_occurred() /= 0) return
+! computes a default hdur_movie that creates nice looking movies.
+! Sets HDUR_MOVIE as the minimum period the mesh can resolve for Southern California model
+ if(HDUR_MOVIE <=TINYVAL .and. (MODEL == 'Harvard_LA' .or. MODEL == 'SoCal')) &
+ HDUR_MOVIE = max(384/NEX_XI*2.4,384/NEX_ETA*2.4)
+
+! compute the minimum value of hdur in CMTSOLUTION file
+ open(unit=1,file=CMTSOLUTION,status='old',action='read')
+ minval_hdur = HUGEVAL
+ do isource = 1,NSOURCES
+
+! skip other information
+ do idummy = 1,3
+ read(1,"(a)") dummystring
+ enddo
+
+! read half duration and compute minimum
+ read(1,"(a)") dummystring
+ read(dummystring(15:len_trim(dummystring)),*) hdur
+ minval_hdur = min(minval_hdur,hdur)
+
+! skip other information
+ do idummy = 1,9
+ read(1,"(a)") dummystring
+ enddo
+
+ enddo
+ close(1)
+
+! one cannot use a Heaviside source for the movies
+! if((MOVIE_SURFACE .or. MOVIE_VOLUME) .and. sqrt(minval_hdur**2 + HDUR_MOVIE**2) < TINYVAL) &
+! stop 'hdur too small for movie creation, movies do not make sense for Heaviside source'
+
+ call read_value_logical(SAVE_MESH_FILES, 'mesher.SAVE_MESH_FILES')
+ if(err_occurred() /= 0) return
+ call read_value_string(LOCAL_PATH, 'LOCAL_PATH')
+ if(err_occurred() /= 0) return
+ call read_value_integer(NTSTEP_BETWEEN_OUTPUT_INFO, 'solver.NTSTEP_BETWEEN_OUTPUT_INFO')
+ if(err_occurred() /= 0) return
+ call read_value_integer(NTSTEP_BETWEEN_OUTPUT_SEISMOS, 'solver.NTSTEP_BETWEEN_OUTPUT_SEISMOS')
+ if(err_occurred() /= 0) return
+ call read_value_logical(PRINT_SOURCE_TIME_FUNCTION, 'solver.PRINT_SOURCE_TIME_FUNCTION')
+ if(err_occurred() /= 0) return
+
+! close parameter file
+ close(IIN)
+
+ end subroutine read_parameter_file
+
Added: seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/read_value_parameters.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/read_value_parameters.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/read_value_parameters.f90 2010-10-01 15:00:21 UTC (rev 17232)
@@ -0,0 +1,177 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology September 2006
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+! read values from parameter file, ignoring white lines and comments
+
+ subroutine read_value_integer(value_to_read, name)
+
+ implicit none
+
+ integer value_to_read
+ character(len=*) name
+ character(len=100) string_read
+
+ call unused_string(name)
+
+ call read_next_line(string_read)
+ read(string_read,*) value_to_read
+
+ end subroutine read_value_integer
+
+!--------------------
+
+ subroutine read_value_double_precision(value_to_read, name)
+
+ implicit none
+
+ double precision value_to_read
+ character(len=*) name
+ character(len=100) string_read
+
+ call unused_string(name)
+
+ call read_next_line(string_read)
+ read(string_read,*) value_to_read
+
+ end subroutine read_value_double_precision
+
+!--------------------
+
+ subroutine read_value_logical(value_to_read, name)
+
+ implicit none
+
+ logical value_to_read
+ character(len=*) name
+ character(len=100) string_read
+
+ call unused_string(name)
+
+ call read_next_line(string_read)
+ read(string_read,*) value_to_read
+
+ end subroutine read_value_logical
+
+!--------------------
+
+ subroutine read_value_string(value_to_read, name)
+
+ implicit none
+
+ character(len=*) value_to_read
+ character(len=*) name
+ character(len=100) string_read
+
+ call unused_string(name)
+
+ call read_next_line(string_read)
+ value_to_read = string_read
+
+ end subroutine read_value_string
+
+!--------------------
+
+ subroutine read_next_line(string_read)
+
+ implicit none
+
+ include "constants.h"
+
+ character(len=100) string_read
+
+ integer index_equal_sign,ios
+
+ do
+ read(unit=IIN,fmt="(a100)",iostat=ios) string_read
+ if(ios /= 0) stop 'error while reading parameter file'
+
+! suppress leading white spaces, if any
+ string_read = adjustl(string_read)
+
+! suppress trailing carriage return (ASCII code 13) if any (e.g. if input text file coming from Windows/DOS)
+ if(index(string_read,achar(13)) > 0) string_read = string_read(1:index(string_read,achar(13))-1)
+
+! exit loop when we find the first line that is not a comment or a white line
+ if(len_trim(string_read) == 0) cycle
+ if(string_read(1:1) /= '#') exit
+
+ enddo
+
+! suppress trailing white spaces, if any
+ string_read = string_read(1:len_trim(string_read))
+
+! suppress trailing comments, if any
+ if(index(string_read,'#') > 0) string_read = string_read(1:index(string_read,'#')-1)
+
+! suppress leading junk (up to the first equal sign, included)
+ index_equal_sign = index(string_read,'=')
+ if(index_equal_sign <= 1 .or. index_equal_sign == len_trim(string_read)) stop 'incorrect syntax detected in DATA/Par_file'
+ string_read = string_read(index_equal_sign + 1:len_trim(string_read))
+
+! suppress leading and trailing white spaces again, if any, after having suppressed the leading junk
+ string_read = adjustl(string_read)
+ string_read = string_read(1:len_trim(string_read))
+
+ end subroutine read_next_line
+
+!--------------------
+
+ subroutine open_parameter_file
+
+ include "constants.h"
+
+ open(unit=IIN,file='DATA/Par_file',status='old',action='read')
+
+ end subroutine open_parameter_file
+
+!--------------------
+
+ subroutine close_parameter_file
+
+ include "constants.h"
+
+ close(IIN)
+
+ end subroutine close_parameter_file
+
+!--------------------
+
+ integer function err_occurred()
+
+ err_occurred = 0
+
+ end function err_occurred
+
+!--------------------
+
+! dummy subroutine to avoid warnings about variable not used in other subroutines
+ subroutine unused_string(s)
+
+ character(len=*) s
+
+ if (len(s) == 1) continue
+
+ end subroutine unused_string
+
Added: seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/readme.txt
===================================================================
--- seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/readme.txt (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/readme.txt 2010-10-01 15:00:21 UTC (rev 17232)
@@ -0,0 +1 @@
+This is an ancient version that still uses the old read/compute_parameter_file subroutines. To compile properly, you need to copy over the right subroutines you used for your runs and modify them accordingly in the main program.
Added: seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/utm_geo.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/utm_geo.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/UTILS/Visualization/GMT/create_movie_GMT/utm_geo.f90 2010-10-01 15:00:21 UTC (rev 17232)
@@ -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., raddeg=180./PI
+ double precision, parameter :: semimaj=6378206.4d0, semimin=6356583.8d0
+ double precision, parameter :: scfa=.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
+
More information about the CIG-COMMITS
mailing list