[cig-commits] r13217 - in seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta: . DATA OUTPUT_FILES src
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Sat Nov 1 17:24:51 PDT 2008
Author: dkomati1
Date: 2008-11-01 17:24:51 -0700 (Sat, 01 Nov 2008)
New Revision: 13217
Added:
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/DATA/CMTSOLUTION_MOVIE_SICHUAN
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/DATA/Par_file_movie_Sichuan
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/constants_topo.h
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/Makefile
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/OUTPUT_FILES/values_from_mesher.h
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/read_compute_parameters.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/specfem3D.F90
Log:
added support for movies (successfully tested at CINES in France)
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/DATA/CMTSOLUTION_MOVIE_SICHUAN
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/DATA/CMTSOLUTION_MOVIE_SICHUAN (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/DATA/CMTSOLUTION_MOVIE_SICHUAN 2008-11-02 00:24:51 UTC (rev 13217)
@@ -0,0 +1,13 @@
+PDE 2000 12 18 1 19 21.60 -21.1800 -179.1200 628.2 6.4 0.0 SICHUAN, CHINA
+event name: 121800A
+time shift: 0
+half duration: 0
+latitude: 31.4400
+longitude: 104.1000
+depth: 55.7700
+Mrr: 5.690000e+27
+Mtt: 1.040000e+26
+Mpp: -5.790001e+27
+Mrt: -1.990000e+27
+Mrp: 5.610000e+27
+Mtp: -3.480000e+27
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/DATA/Par_file_movie_Sichuan
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/DATA/Par_file_movie_Sichuan (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/DATA/Par_file_movie_Sichuan 2008-11-02 00:24:51 UTC (rev 13217)
@@ -0,0 +1,109 @@
+
+# forward or adjoint simulation
+SIMULATION_TYPE = 1
+SAVE_FORWARD = .false. # save last frame of forward simulation or not
+
+# number of chunks (1,2,3 or 6)
+NCHUNKS = 6
+
+# angular width of the first chunk (not used if full sphere with six chunks)
+ANGULAR_WIDTH_XI_IN_DEGREES = 90.d0 # angular size of a chunk
+ANGULAR_WIDTH_ETA_IN_DEGREES = 90.d0
+CENTER_LATITUDE_IN_DEGREES = 90.d0
+CENTER_LONGITUDE_IN_DEGREES = 0.d0
+GAMMA_ROTATION_AZIMUTH = 0.d0
+
+# number of elements at the surface along the two sides of the first chunk
+# (must be multiple of 16 and also of 8 * multiple of NPROC below)
+NEX_XI = 512
+NEX_ETA = 512
+
+# number of MPI processors along the two sides of the first chunk
+NPROC_XI = 16
+NPROC_ETA = 16
+
+# 1D models with real structure:
+# 1D_isotropic_prem, 1D_transversely_isotropic_prem, 1D_iasp91, 1D_1066a, 1D_ak135, 1D_ref, 1D_ref_iso, 1D_jp3d,1D_sea99
+#
+# 1D models with only one fictitious averaged crustal layer:
+# 1D_isotropic_prem_onecrust, 1D_transversely_isotropic_prem_onecrust, 1D_iasp91_onecrust, 1D_1066a_onecrust, 1D_ak135_onecrust
+#
+# fully 3D models:
+# transversely_isotropic_prem_plus_3D_crust_2.0, 3D_anisotropic, 3D_attenuation,
+# s20rts, s362ani, s362iso, s362wmani, s362ani_prem, s29ea, s29ea,sea99_jp3d1994,sea99,jp3d1994
+MODEL = s20rts
+
+# parameters describing the Earth model
+OCEANS = .false.
+ELLIPTICITY = .true.
+TOPOGRAPHY = .false.
+GRAVITY = .false.
+ROTATION = .false.
+ATTENUATION = .true.
+
+# absorbing boundary conditions for a regional simulation
+ABSORBING_CONDITIONS = .false.
+
+# record length in minutes
+RECORD_LENGTH_IN_MINUTES = 175.0d0
+
+# save AVS or OpenDX movies
+MOVIE_SURFACE = .true.
+MOVIE_VOLUME = .false.
+MOVIE_COARSE = .false.
+NTSTEP_BETWEEN_FRAMES = 50
+HDUR_MOVIE = 10.125d0
+
+# save movie in volume. Will save element if center of element is in prescribed volume
+# top/bottom: depth in KM, use MOVIE_TOP = -100 to make sure the surface is stored.
+# west/east: longitude, degrees East [-180/180] top/bottom: latitute, degrees North [-90/90]
+# start/stop: frames will be stored at MOVIE_START + i*NSTEP_BETWEEN_FRAMES, where i=(0,1,2..) and iNSTEP_BETWEEN_FRAMES <= MOVIE_STOP
+# movie_volume_type: 1=strain, 2=time integral of strain, 3=\mu*time integral of strain
+# type 4 saves the trace and deviatoric stress in the whole volume, 5=displacement, 6=velocity
+#MOVIE_COARSE saves movie only at corners of elements
+MOVIE_VOLUME_TYPE = 2
+MOVIE_TOP_KM = -100.0
+MOVIE_BOTTOM_KM = 1000.0
+MOVIE_WEST_DEG = -90.0
+MOVIE_EAST_DEG = 90.0
+MOVIE_NORTH_DEG = 90.0
+MOVIE_SOUTH_DEG = -90.0
+MOVIE_START = 0
+MOVIE_STOP = 40000
+
+# save mesh files to check the mesh
+SAVE_MESH_FILES = .false.
+
+# restart files
+IT_LAST_VALUE_DUMPED = 0
+INTERVAL_DUMP_FILES = 2000
+
+# interval at which we output time step info and max of norm of displacement
+NTSTEP_BETWEEN_OUTPUT_INFO = 100
+
+# interval in time steps for temporary writing of seismograms
+NTSTEP_BETWEEN_OUTPUT_SEISMOS = 5000000
+NTSTEP_BETWEEN_READ_ADJSRC = 1000
+
+# output format for the seismograms (one can use either or all of the three formats)
+OUTPUT_SEISMOS_ASCII_TEXT = .true.
+OUTPUT_SEISMOS_SAC_ALPHANUM = .false.
+OUTPUT_SEISMOS_SAC_BINARY = .false.
+
+# rotate seismograms to Radial-Transverse-Z or use default North-East-Z reference frame
+ROTATE_SEISMOGRAMS_RT = .false.
+
+# decide if master process writes all the seismograms or if all processes do it in parallel
+WRITE_SEISMOGRAMS_BY_MASTER = .true.
+
+# save all seismograms in one large combined file instead of one file per seismogram
+# to avoid overloading shared non-local file systems such as GPFS for instance
+SAVE_ALL_SEISMOS_IN_ONE_FILE = .true.
+USE_BINARY_FOR_LARGE_FILE = .true.
+
+# flag to impose receivers at the surface or allow them to be buried
+RECEIVERS_CAN_BE_BURIED = .false.
+
+# print source time function
+PRINT_SOURCE_TIME_FUNCTION = .false.
+
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/Makefile
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/Makefile 2008-11-01 23:15:42 UTC (rev 13216)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/Makefile 2008-11-02 00:24:51 UTC (rev 13217)
@@ -136,6 +136,7 @@
$O/count_number_of_sources.o \
$O/create_central_cube_buffers.o \
$O/create_chunk_buffers.o \
+ $O/create_movie_AVS_DX.o \
$O/create_regions_mesh.o \
$O/crustal_model.o \
$O/define_derivation_matrices.o \
@@ -215,6 +216,7 @@
$(OUTPUT_FILES_INC)/values_from_mesher.h \
xspecfem3D \
xcombine_AVS_DX \
+ xcreate_movie_AVS_DX \
$(EMPTY_MACRO)
default: $(DEFAULT)
@@ -243,6 +245,9 @@
xcombine_AVS_DX: $O/combine_AVS_DX.o $(LIBSPECFEM)
${FCCOMPILE_CHECK} -o $(BIN)/xcombine_AVS_DX $O/combine_AVS_DX.o $(LIBSPECFEM)
+xcreate_movie_AVS_DX: $O/create_movie_AVS_DX.o $(LIBSPECFEM)
+ ${FCCOMPILE_CHECK} -o $(BIN)/xcreate_movie_AVS_DX $O/create_movie_AVS_DX.o $(LIBSPECFEM)
+
xconvolve_source_timefunction: $O/convolve_source_timefunction.o
${FCCOMPILE_CHECK} -o $(BIN)/xconvolve_source_timefunction $O/convolve_source_timefunction.o
@@ -250,7 +255,7 @@
${MPIFCCOMPILE_CHECK} -o $(BIN)/xcreate_header_file $O/create_header_file.o $O/exit_mpi.o $O/get_value_parameters.o $O/read_compute_parameters.o $O/memory_eval.o $O/save_header_file.o $O/count_number_of_sources.o $O/read_value_parameters.o $O/euler_angles.o $O/reduce.o $O/rthetaphi_xyz.o $O/auto_ner.o
clean:
- rm -f $O/* *.o work.pc* *.mod $(BIN)/xspecfem3D $(BIN)/xconvolve_source_timefunction $(BIN)/xcreate_header_file $(BIN)/xcombine_AVS_DX
+ rm -f $O/* *.o work.pc* *.mod $(BIN)/xspecfem3D $(BIN)/xconvolve_source_timefunction $(BIN)/xcreate_header_file $(BIN)/xcombine_AVS_DX $(BIN)/xcreate_movie_AVS_DX
###
@@ -312,6 +317,9 @@
$O/comp_source_time_function.o: $S/comp_source_time_function.f90
${FCCOMPILE_CHECK} -c -o $O/comp_source_time_function.o ${FCFLAGS_f90} $S/comp_source_time_function.f90
+$O/create_movie_AVS_DX.o: $(SPECINC)/constants.h $S/create_movie_AVS_DX.f90
+ ${FCCOMPILE_CHECK} -c -o $O/create_movie_AVS_DX.o ${FCFLAGS_f90} $S/create_movie_AVS_DX.f90
+
## leave MPIFLAGS below because we use the C preprocessor for that file
$O/combine_AVS_DX.o: $(SPECINC)/constants.h $S/combine_AVS_DX.F90
${FCCOMPILE_CHECK} $(MPIFLAGS) -c -o $O/combine_AVS_DX.o ${FCFLAGS_f90} $S/combine_AVS_DX.F90
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/OUTPUT_FILES/values_from_mesher.h
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/OUTPUT_FILES/values_from_mesher.h 2008-11-01 23:15:42 UTC (rev 13216)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/OUTPUT_FILES/values_from_mesher.h 2008-11-02 00:24:51 UTC (rev 13217)
@@ -12,23 +12,23 @@
!
! number of processors = 96
!
- ! maximum number of points per region = 576013
+ ! maximum number of points per region = 509449
!
! on NEC SX, make sure "loopcnt=" parameter
- ! in Makefile is greater than max vector length = 1728039
+ ! in Makefile is greater than max vector length = 1528347
!
- ! total elements per slice = 9616
- ! total points per slice = 645095
+ ! total elements per slice = 8592
+ ! total points per slice = 578531
!
! total for full 6-chunk mesh:
! ---------------------------
!
! exact total number of spectral elements in entire mesh =
- ! 902656.000000000
+ ! 804352.000000000
! approximate total number of points in entire mesh =
- ! 60555995.0000000
+ ! 54165851.0000000
! approximate total number of degrees of freedom in entire mesh =
- ! 172454865.000000
+ ! 153284433.000000
!
! resolution of the mesh at the surface:
! -------------------------------------
@@ -39,8 +39,8 @@
! average distance between points in km = 19.54598
! average size of a spectral element in km = 78.18393
!
- ! number of time steps = 7600
- ! value of a time step DT = 0.200000000000000
+ ! number of time steps = 10100
+ ! value of a time step DT = 0.150000000000000
!
! number of seismic sources = 1
!
@@ -48,7 +48,7 @@
! approximate static memory needed by the solver:
! ----------------------------------------------
!
- ! size of static arrays per slice = 0.178467992693186 GB
+ ! size of static arrays per slice = 0.155213307589293 GB
!
! (should be below and typically equal to 80% of 1.5 GB = 1.2 GB on pangu
! at Caltech, and below and typically equal to 85% of 2 GB = 1.7 GB
@@ -56,31 +56,32 @@
! (if significantly more, the job will not run by lack of memory)
! (if significantly less, you waste a significant amount of memory)
!
- ! size of static arrays for all slices = 17.1329272985458 GB
- ! = 1.673137431498617E-002 TB
+ ! size of static arrays for all slices = 14.9004775285721 GB
+ ! = 1.455124758649617E-002 TB
!
integer, parameter :: NEX_XI_VAL = 128
integer, parameter :: NEX_ETA_VAL = 128
- integer, parameter :: NSPEC_CRUST_MANTLE = 8640
+ integer, parameter :: NSPEC_CRUST_MANTLE = 7616
integer, parameter :: NSPEC_OUTER_CORE = 688
integer, parameter :: NSPEC_INNER_CORE = 288
- integer, parameter :: NGLOB_CRUST_MANTLE = 576013
+ integer, parameter :: NGLOB_CRUST_MANTLE = 509449
integer, parameter :: NGLOB_OUTER_CORE = 47985
integer, parameter :: NGLOB_INNER_CORE = 21097
integer, parameter :: NSPECMAX_ANISO_IC = 1
- integer, parameter :: NSPECMAX_ISO_MANTLE = 8640
+ integer, parameter :: NSPECMAX_ISO_MANTLE = 7616
integer, parameter :: NSPECMAX_TISO_MANTLE = 2304
integer, parameter :: NSPECMAX_ANISO_MANTLE = 1
- integer, parameter :: NSPEC_CRUST_MANTLE_ATTENUAT = 8640
+ integer, parameter :: NSPEC_CRUST_MANTLE_ATTENUAT = 7616
+ integer, parameter :: NSPEC_CRUST_MANTLE_ATTENUAT3D = 1
integer, parameter :: NSPEC_INNER_CORE_ATTENUATION = 288
- integer, parameter :: NSPEC_CRUST_MANTLE_STR_OR_ATT = 8640
+ integer, parameter :: NSPEC_CRUST_MANTLE_STR_OR_ATT = 7616
integer, parameter :: NSPEC_INNER_CORE_STR_OR_ATT = 288
integer, parameter :: NSPEC_CRUST_MANTLE_STR_AND_ATT = 1
@@ -100,7 +101,7 @@
integer, parameter :: NSPEC_CRUST_MANTLE_STACEY = 1
integer, parameter :: NSPEC_OUTER_CORE_STACEY = 1
- integer, parameter :: NGLOB_CRUST_MANTLE_OCEANS = 576013
+ integer, parameter :: NGLOB_CRUST_MANTLE_OCEANS = 1
logical, parameter :: TRANSVERSE_ISOTROPY_VAL = .true.
@@ -119,20 +120,20 @@
logical, parameter :: ROTATION_VAL = .false.
integer, parameter :: NSPEC_OUTER_CORE_ROTATION = 1
- integer, parameter :: NGLOB1D_RADIAL_CM = 109
+ integer, parameter :: NGLOB1D_RADIAL_CM = 105
integer, parameter :: NGLOB1D_RADIAL_OC = 65
integer, parameter :: NGLOB1D_RADIAL_IC = 9
- integer, parameter :: NGLOB2DMAX_XMIN_XMAX_CM = 8574
+ integer, parameter :: NGLOB2DMAX_XMIN_XMAX_CM = 7934
integer, parameter :: NGLOB2DMAX_XMIN_XMAX_OC = 2134
integer, parameter :: NGLOB2DMAX_XMIN_XMAX_IC = 1283
- integer, parameter :: NGLOB2DMAX_YMIN_YMAX_CM = 8574
+ integer, parameter :: NGLOB2DMAX_YMIN_YMAX_CM = 7934
integer, parameter :: NGLOB2DMAX_YMIN_YMAX_OC = 2134
integer, parameter :: NGLOB2DMAX_YMIN_YMAX_IC = 1283
integer, parameter :: NPROC_XI_VAL = 4
integer, parameter :: NPROC_ETA_VAL = 4
integer, parameter :: NCHUNKS_VAL = 6
integer, parameter :: NPROCTOT_VAL = 96
- integer, parameter :: NGLOB2DMAX_XY_VAL_CM = 8574
+ integer, parameter :: NGLOB2DMAX_XY_VAL_CM = 7934
integer, parameter :: NGLOB2DMAX_XY_VAL_OC = 2134
integer, parameter :: NGLOB2DMAX_XY_VAL_IC = 1283
integer, parameter :: NUMMSGS_FACES_VAL = 48
@@ -142,8 +143,8 @@
integer, parameter :: ATT3 = 1
integer, parameter :: ATT4 = 70000
integer, parameter :: ATT5 = 70000
- integer, parameter :: NSPEC2DMAX_XMIN_XMAX_CM = 440
- integer, parameter :: NSPEC2DMAX_YMIN_YMAX_CM = 440
+ integer, parameter :: NSPEC2DMAX_XMIN_XMAX_CM = 408
+ integer, parameter :: NSPEC2DMAX_YMIN_YMAX_CM = 408
integer, parameter :: NSPEC2D_BOTTOM_CM = 64
integer, parameter :: NSPEC2D_TOP_CM = 1024
integer, parameter :: NSPEC2DMAX_XMIN_XMAX_IC = 72
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/constants_topo.h
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/constants_topo.h (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/constants_topo.h 2008-11-02 00:24:51 UTC (rev 13217)
@@ -0,0 +1,515 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 4 . 1
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory, California Institute of Technology, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) California Institute of Technology and University of Pau / CNRS / INRIA
+! August 2008
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+! constants.h. Generated from constants.h.in by configure.
+
+!
+!--- user can modify parameters below
+!
+
+! exaggeration factor to display topography in OpenDX
+ double precision, parameter :: EXAGGERATION_FACTOR_TOPO = 22.d0 ! 35.d0
+
+! this for non blocking assembly
+ logical, parameter :: USE_NONBLOCKING_COMMS = .true.
+ integer, parameter :: ELEMENTS_NONBLOCKING_CM_IC = 1500
+ integer, parameter :: ELEMENTS_NONBLOCKING_OC = 3000
+
+ logical, parameter :: DEBUG_NONBLOCKING_COMMS = .false.
+ logical, parameter :: DEBUG_USING_OPENDX = .false.
+
+!! DK DK temporary patch for the large Gordon Bell runs: set RECEIVERS_CAN_BE_BURIED
+!! DK DK to false in all cases etc
+ logical, parameter :: PATCH_FOR_GORDON_BELL = .true.
+
+! save partial seismograms every 10,000 time steps or not
+ logical, parameter :: SAVE_PARTIAL_SEISMOGRAMS = .false.
+
+! (much) faster detection of receivers at high resolution: use grid points only
+ logical, parameter :: FASTER_RECEIVERS_POINTS_ONLY = .true.
+
+! suppress calculation and storage of seismograms if needed
+ logical, parameter :: COMPUTE_STORE_SEISMOGRAMS = .true.
+
+!! DK DK for Gordon Bell
+! integer, parameter :: SEA99_VS_DIM1 = 100, SEA99_VS_DIM2 = 100, SEA99_VS_DIM3 = 100
+ integer, parameter :: SEA99_VS_DIM1 = 1, SEA99_VS_DIM2 = 1, SEA99_VS_DIM3 = 1
+
+!! DK DK for Gordon Bell
+! integer, parameter :: AMM_V_DIM1 = 14, AMM_V_DIM2 = 34, AMM_V_DIM3 = 37, AMM_V_DIM4 = 73
+! use 6 for first index to avoid a warning when compiling
+ integer, parameter :: AMM_V_DIM1 = 6, AMM_V_DIM2 = 1, AMM_V_DIM3 = 1, AMM_V_DIM4 = 1
+
+!
+! solver in single or double precision depending on the machine (4 or 8 bytes)
+!
+! ALSO CHANGE FILE precision.h ACCORDINGLY
+!
+ integer, parameter :: SIZE_REAL = 4, SIZE_DOUBLE = 8
+
+! usually the size of integer and logical variables is the same as regular single-precision real variable
+ integer, parameter :: SIZE_INTEGER = SIZE_REAL
+ integer, parameter :: SIZE_LOGICAL = SIZE_REAL
+
+! set to SIZE_REAL to run in single precision
+! set to SIZE_DOUBLE to run in double precision (increases memory size by 2)
+ integer, parameter :: CUSTOM_REAL = SIZE_REAL
+
+! input, output and main MPI I/O files
+ integer, parameter :: ISTANDARD_OUTPUT = 6
+ integer, parameter :: IIN = 40,IOUT = 41,IOUT_SAC = 903
+! local file unit for output of buffers
+ integer, parameter :: IOUT_BUFFERS = 35
+! uncomment this to write messages to a text file
+ integer, parameter :: IMAIN = 42
+! uncomment this to write messages to the screen (slows down the code)
+! integer, parameter :: IMAIN = ISTANDARD_OUTPUT
+
+! number of values read in Par_file that we need to broadcast
+ integer, parameter :: NVALUES_bcast_integer = 41
+ integer, parameter :: NVALUES_bcast_double_precision = 30
+ integer, parameter :: NVALUES_bcast_logical = 33
+
+! R_EARTH is the radius of the bottom of the oceans (radius of Earth in m)
+ double precision, parameter :: R_EARTH = 6371000.d0
+! uncomment line below for PREM with oceans
+! double precision, parameter :: R_EARTH = 6368000.d0
+
+! average density in the full Earth to normalize equation
+ double precision, parameter :: RHOAV = 5514.3d0
+
+! for topography/bathymetry model
+
+!!--- ETOPO5 5-minute model, smoothed Harvard version
+!! size of topography and bathymetry file
+! integer, parameter :: NX_BATHY = 4320,NY_BATHY = 2160
+!! resolution of topography file in minutes
+! integer, parameter :: RESOLUTION_TOPO_FILE = 5
+!! pathname of the topography file
+! character(len=*), parameter :: PATHNAME_TOPO_FILE = 'DATA/topo_bathy/topo_bathy_etopo5_smoothed_Harvard.dat'
+
+!--- ETOPO4 4-minute model created by subsampling and smoothing etopo-2
+! size of topography and bathymetry file
+ integer, parameter :: NX_BATHY = 5400,NY_BATHY = 2700
+!!!!!!!!!! integer, parameter :: NX_BATHY = 1,NY_BATHY = 1
+! resolution of topography file in minutes
+ integer, parameter :: RESOLUTION_TOPO_FILE = 4
+! pathname of the topography file
+! character(len=*), parameter :: PATHNAME_TOPO_FILE = 'DATA/topo_bathy/topo_bathy_etopo4_smoothed_window_7.dat'
+ character(len=*), parameter :: PATHNAME_TOPO_FILE = 'DATA/topo_bathy/topo_bathy_etopo4_smoothed_window_11_minmax_3500.dat'
+
+!!--- ETOPO2 2-minute model, not implemented yet
+!! size of topography and bathymetry file
+! integer, parameter :: NX_BATHY = 10800,NY_BATHY = 5400
+!! resolution of topography file in minutes
+! integer, parameter :: RESOLUTION_TOPO_FILE = 2
+!! pathname of the topography file
+! character(len=*), parameter :: PATHNAME_TOPO_FILE = 'DATA/topo_bathy/topo_bathy_etopo2_smoothed_window7.dat'
+
+! minimum thickness in meters to include the effect of the oceans and topo
+ double precision, parameter :: MINIMUM_THICKNESS_3D_OCEANS = 100.d0
+
+! number of GLL points in each direction of an element (degree plus one)
+ integer, parameter :: NGLLX = 5
+ integer, parameter :: NGLLY = NGLLX
+ integer, parameter :: NGLLZ = NGLLX
+
+! flag to exclude elements that are too far from target in source detection
+ logical, parameter :: USE_DISTANCE_CRITERION = .true.
+
+! flag to display detailed information about location of stations
+ logical, parameter :: DISPLAY_DETAILS_STATIONS = .false.
+
+! maximum length of station and network name for receivers
+ integer, parameter :: MAX_LENGTH_STATION_NAME = 32
+ integer, parameter :: MAX_LENGTH_NETWORK_NAME = 8
+
+! we mimic a triangle of half duration equal to half_duration_triangle
+! using a Gaussian having a very close shape, as explained in Figure 4.2
+! of the manual. This source decay rate to mimic an equivalent triangle
+! was found by trial and error
+ double precision, parameter :: SOURCE_DECAY_MIMIC_TRIANGLE = 1.628d0
+
+! maximum number of sources to locate simultaneously (do not set it too high otherwise you will use a lot of memory)
+ integer, parameter :: NSOURCES_SUBSET_MAX = 100
+
+! distance threshold (in km) above which we consider that a receiver
+! is located outside the mesh and therefore excluded from the station list
+ double precision, parameter :: THRESHOLD_EXCLUDE_STATION = 50.d0
+
+! the first doubling is implemented right below the Moho
+! it seems optimal to implement the three other doublings at these depths
+! in the mantle
+ double precision, parameter :: DEPTH_SECOND_DOUBLING_OPTIMAL = 1650000.d0
+! in the outer core
+ double precision, parameter :: DEPTH_THIRD_DOUBLING_OPTIMAL = 3860000.d0
+! in the outer core
+ double precision, parameter :: DEPTH_FOURTH_DOUBLING_OPTIMAL = 5000000.d0
+
+! Boundary Mesh -- save Moho, 400, 670 km discontinuity topology files (in
+! the mesher) and use them for the computation of boundary kernel (in the solver)
+ logical, parameter :: SAVE_BOUNDARY_MESH = .false.
+
+! this parameter must be set to .true. to compute anisotropic kernels
+! in crust and mantle (related to the 21 Cij in geographical coordinates)
+! default is .false. to compute isotropic kernels (related to alpha and beta)
+ logical, parameter :: ANISOTROPIC_KL = .false.
+
+! print date and time estimate of end of run in another country,
+! in addition to local time.
+! For instance: the code runs at Caltech in California but the person
+! running the code is connected remotely from France, which has 9 hours more.
+! The time difference with that remote location can be positive or negative
+ logical, parameter :: ADD_TIME_ESTIMATE_ELSEWHERE = .false.
+ integer, parameter :: HOURS_TIME_DIFFERENCE = +9
+ integer, parameter :: MINUTES_TIME_DIFFERENCE = +0
+
+!
+!--- debugging flags
+!
+
+! flags to actually assemble with MPI or not
+! and to actually match fluid and solid regions of the Earth or not
+! should always be set to true except when debugging code
+ logical, parameter :: ACTUALLY_ASSEMBLE_MPI_SLICES = .true.
+ logical, parameter :: ACTUALLY_ASSEMBLE_MPI_CHUNKS = .true.
+ logical, parameter :: ACTUALLY_COUPLE_FLUID_CMB = .true.
+ logical, parameter :: ACTUALLY_COUPLE_FLUID_ICB = .true.
+
+! flag to turn off the conversion of geographic to geocentric coordinates for
+! the seismic source and the stations; i.e. assume a perfect sphere, which
+! can be useful for benchmarks of a spherical Earth with fictitious sources and stations
+ logical, parameter :: ASSUME_PERFECT_SPHERE = .false.
+
+! flag to only create the mesh but not start the solver (for instance to check the mesh obtained)
+ logical, parameter :: MESHER_ONLY = .false.
+
+! flag to put a fictitious source in each region in the case of a serial test
+ logical, parameter :: PUT_SOURCE_IN_EACH_REGION = .false.
+
+!------------------------------------------------------
+!----------- do not modify anything below -------------
+!------------------------------------------------------
+
+! on some processors (e.g. Pentiums) it is necessary to suppress underflows
+! by using a small initial field instead of zero
+ logical, parameter :: FIX_UNDERFLOW_PROBLEM = .true.
+
+! some useful constants
+ double precision, parameter :: PI = 3.141592653589793d0
+ double precision, parameter :: TWO_PI = 2.d0 * PI
+ double precision, parameter :: PI_OVER_FOUR = PI / 4.d0
+
+! to convert angles from degrees to radians
+ double precision, parameter :: DEGREES_TO_RADIANS = PI / 180.d0
+
+! 3-D simulation
+ integer, parameter :: NDIM = 3
+
+! dimension of the boundaries of the slices
+ integer, parameter :: NDIM2D = 2
+
+! number of nodes for 2D and 3D shape functions for hexahedra with 27 nodes
+ integer, parameter :: NGNOD = 27, NGNOD2D = 9
+
+! gravitational constant
+ double precision, parameter :: GRAV = 6.6723d-11
+
+! a few useful constants
+ double precision, parameter :: ZERO = 0.d0,ONE = 1.d0,TWO = 2.d0,HALF = 0.5d0
+
+ real(kind=CUSTOM_REAL), parameter :: &
+ ONE_THIRD = 1._CUSTOM_REAL/3._CUSTOM_REAL, &
+ TWO_THIRDS = 2._CUSTOM_REAL/3._CUSTOM_REAL, &
+ FOUR_THIRDS = 4._CUSTOM_REAL/3._CUSTOM_REAL
+
+! very large and very small values
+ double precision, parameter :: HUGEVAL = 1.d+30,TINYVAL = 1.d-9
+
+! very large real value declared independently of the machine
+ real(kind=CUSTOM_REAL), parameter :: HUGEVAL_SNGL = 1.e+30_CUSTOM_REAL
+
+! very large integer value
+ integer, parameter :: HUGEINT = 100000000
+
+! normalized radius of free surface
+ double precision, parameter :: R_UNIT_SPHERE = ONE
+
+! same radius in km
+ double precision, parameter :: R_EARTH_KM = R_EARTH / 1000.d0
+
+! fixed thickness of 3 km for PREM oceans
+ double precision, parameter :: THICKNESS_OCEANS_PREM = 3000.d0 / R_EARTH
+
+! shortest radius at which crust is implemented (80 km depth)
+! to be constistent with the D80 discontinuity, we impose the crust only above it
+ double precision, parameter :: R_DEEPEST_CRUST = (R_EARTH - 80000.d0) / R_EARTH
+
+! maximum number of chunks (full sphere)
+ integer, parameter :: NCHUNKS_MAX = 6
+
+! define block type based upon chunk number (between 1 and 6)
+! do not change this numbering, chunk AB must be number 1 for central cube
+ integer, parameter :: CHUNK_AB = 1
+ integer, parameter :: CHUNK_AC = 2
+ integer, parameter :: CHUNK_BC = 3
+ integer, parameter :: CHUNK_AC_ANTIPODE = 4
+ integer, parameter :: CHUNK_BC_ANTIPODE = 5
+ integer, parameter :: CHUNK_AB_ANTIPODE = 6
+
+! maximum number of regions in the mesh
+ integer, parameter :: MAX_NUM_REGIONS = 3
+
+! define flag for regions of the global Earth mesh
+ integer, parameter :: IREGION_CRUST_MANTLE = 1
+ integer, parameter :: IREGION_OUTER_CORE = 2
+ integer, parameter :: IREGION_INNER_CORE = 3
+
+! define flag for elements
+ integer, parameter :: IFLAG_CRUST = 1
+
+ integer, parameter :: IFLAG_80_MOHO = 2
+ integer, parameter :: IFLAG_220_80 = 3
+ integer, parameter :: IFLAG_670_220 = 4
+ integer, parameter :: IFLAG_MANTLE_NORMAL = 5
+
+ integer, parameter :: IFLAG_OUTER_CORE_NORMAL = 6
+
+ integer, parameter :: IFLAG_INNER_CORE_NORMAL = 7
+ integer, parameter :: IFLAG_MIDDLE_CENTRAL_CUBE = 8
+ integer, parameter :: IFLAG_BOTTOM_CENTRAL_CUBE = 9
+ integer, parameter :: IFLAG_TOP_CENTRAL_CUBE = 10
+ integer, parameter :: IFLAG_IN_FICTITIOUS_CUBE = 11
+
+ integer, parameter :: NSPEC2D_XI_SUPERBRICK = 8
+ integer, parameter :: NSPEC2D_ETA_SUPERBRICK = 8
+ integer, parameter :: NSPEC2D_XI_SUPERBRICK_1L = 6
+ integer, parameter :: NSPEC2D_ETA_SUPERBRICK_1L = 6
+
+! dummy flag used for mesh display purposes only
+ integer, parameter :: IFLAG_DUMMY = 100
+
+! max number of layers that are used in the radial direction to build the full mesh
+ integer, parameter :: MAX_NUMBER_OF_MESH_LAYERS = 15
+
+! define number of spectral elements and points in basic symmetric mesh doubling superbrick
+ integer, parameter :: NSPEC_DOUBLING_SUPERBRICK = 32
+ integer, parameter :: NGLOB_DOUBLING_SUPERBRICK = 67
+ integer, parameter :: NSPEC_SUPERBRICK_1L = 28
+ integer, parameter :: NGLOB_SUPERBRICK_1L = 58
+ integer, parameter :: NGNOD_EIGHT_CORNERS = 8
+
+! define flag for reference 1D Earth model
+ integer, parameter :: REFERENCE_MODEL_PREM = 1
+ integer, parameter :: REFERENCE_MODEL_IASP91 = 2
+ integer, parameter :: REFERENCE_MODEL_1066A = 3
+ integer, parameter :: REFERENCE_MODEL_AK135 = 4
+ integer, parameter :: REFERENCE_MODEL_REF = 5
+ integer, parameter :: REFERENCE_MODEL_JP1D = 6
+ integer, parameter :: REFERENCE_MODEL_SEA1D = 7
+
+! define flag for 3D Earth model
+ integer, parameter :: THREE_D_MODEL_S20RTS = 1
+ integer, parameter :: THREE_D_MODEL_S362ANI = 2
+ integer, parameter :: THREE_D_MODEL_S362WMANI = 3
+ integer, parameter :: THREE_D_MODEL_S362ANI_PREM = 4
+ integer, parameter :: THREE_D_MODEL_S29EA = 5
+ integer, parameter :: THREE_D_MODEL_SEA99_JP3D = 6
+ integer, parameter :: THREE_D_MODEL_SEA99 = 7
+ integer, parameter :: THREE_D_MODEL_JP3D = 8
+
+! define flag for regions of the global Earth for attenuation
+ integer, parameter :: NUM_REGIONS_ATTENUATION = 5
+
+ integer, parameter :: IREGION_ATTENUATION_INNER_CORE = 1
+ integer, parameter :: IREGION_ATTENUATION_CMB_670 = 2
+ integer, parameter :: IREGION_ATTENUATION_670_220 = 3
+ integer, parameter :: IREGION_ATTENUATION_220_80 = 4
+ integer, parameter :: IREGION_ATTENUATION_80_SURFACE = 5
+ integer, parameter :: IREGION_ATTENUATION_UNDEFINED = 6
+
+! number of standard linear solids for attenuation
+ integer, parameter :: N_SLS = 3
+
+! computation of standard linear solids in meshfem3D
+! ATTENUATION_COMP_RESOLUTION: Number of Digits after decimal
+! ATTENUATION_COMP_MAXIMUM: Maximum Q Value
+ integer, parameter :: ATTENUATION_COMP_RESOLUTION = 1
+ integer, parameter :: ATTENUATION_COMP_MAXIMUM = 5000
+
+! for lookup table for attenuation every 100 m in radial direction of Earth model
+ integer, parameter :: NRAD_ATTENUATION = 70000
+ double precision, parameter :: TABLE_ATTENUATION = R_EARTH_KM * 10.0d0
+
+!!!!!!!!!!!!!! for determination of the attenuation period range
+!!!!!!!!!!!!!! if this is set to .true. then the hardcoded values will be used
+!!!!!!!!!!!!!! otherwise they are computed automatically from the number of elements
+!!!!!!!!!!!!!! This *may* be a useful parameter for benchmarking against older versions
+!!!!!!!!!!!!! logical, parameter :: ATTENUATION_RANGE_PREDEFINED = .false.
+
+! flag for the four edges of each slice and for the bottom edge
+ integer, parameter :: XI_MIN = 1
+ integer, parameter :: XI_MAX = 2
+ integer, parameter :: ETA_MIN = 3
+ integer, parameter :: ETA_MAX = 4
+ integer, parameter :: BOTTOM = 5
+
+! flags to select the right corner in each slice
+ integer, parameter :: ILOWERLOWER = 1
+ integer, parameter :: ILOWERUPPER = 2
+ integer, parameter :: IUPPERLOWER = 3
+ integer, parameter :: IUPPERUPPER = 4
+
+! number of points in each AVS or OpenDX quadrangular cell for movies
+ integer, parameter :: NGNOD2D_AVS_DX = 4
+
+! number of faces a given slice can share with other slices
+! this is at most 2, except when there is only once slice per chunk
+! in which case it is 4
+ integer, parameter :: NUMFACES_SHARED = 2 !!!!! DK DK removed support for one slice only 4
+
+! number of corners a given slice can share with other slices
+! this is at most 1, except when there is only once slice per chunk
+! in which case it is 4
+ integer, parameter :: NUMCORNERS_SHARED = 1 !!!!!! DK DK removed support for one slice only 4
+
+! number of layers in PREM
+ integer, parameter :: NR = 640
+
+! smallest real number on many machines = 1.1754944E-38
+! largest real number on many machines = 3.4028235E+38
+! small negligible initial value to avoid very slow underflow trapping
+! but not too small to avoid trapping on velocity and acceleration in Newmark
+ real(kind=CUSTOM_REAL), parameter :: VERYSMALLVAL = 1.E-24_CUSTOM_REAL
+
+! displacement threshold above which we consider that the code became unstable
+ real(kind=CUSTOM_REAL), parameter :: STABILITY_THRESHOLD = 1.E+25_CUSTOM_REAL
+
+! geometrical tolerance for boundary detection
+ double precision, parameter :: SMALLVAL = 0.00001d0
+
+! small tolerance for conversion from x y z to r theta phi
+ double precision, parameter :: SMALL_VAL_ANGLE = 1.d-10
+
+! geometry tolerance parameter to calculate number of independent grid points
+! sensitive to actual size of model, assumes reference sphere of radius 1
+! this is an absolute value for normalized coordinates in the Earth
+ double precision, parameter :: SMALLVALTOL = 1.d-10
+
+! do not use tags for MPI messages, use dummy tag instead
+ integer, parameter :: itag = 0,itag2 = 0
+
+! for the Gauss-Lobatto-Legendre points and weights
+ double precision, parameter :: GAUSSALPHA = 0.d0,GAUSSBETA = 0.d0
+
+! number of lines per source in CMTSOLUTION file
+ integer, parameter :: NLINES_PER_CMTSOLUTION_SOURCE = 13
+
+! number of iterations to solve the non linear system for xi and eta
+ integer, parameter :: NUM_ITER = 4
+
+! number of hours per day for rotation rate of the Earth
+ double precision, parameter :: HOURS_PER_DAY = 24.d0
+
+! for lookup table for gravity every 100 m in radial direction of Earth model
+ integer, parameter :: NRAD_GRAVITY = 70000
+
+! number of layers in DATA/1066a/1066a.dat
+ integer, parameter :: NR_1066A = 160
+
+! number of layers in DATA/ak135/ak135.dat
+ integer, parameter :: NR_AK135 = 144
+
+! number of layers in DATA/s362ani/REF
+ integer, parameter :: NR_REF = 750
+
+! number of layers in DATA/Lebedev_sea99 1D model
+ integer, parameter :: NR_SEA1D = 163
+
+! three_d_mantle_model_constants
+ integer, parameter :: NK = 20,NS = 20,ND = 1
+
+! Japan 3D model (Zhao, 1994) constants
+ integer, parameter :: MPA=42,MRA=48,MHA=21,MPB=42,MRB=48,MHB=18
+ integer, parameter :: MKA=2101,MKB=2101
+
+! The meaningful range of Zhao et al.'s model (1994) is as follows:
+! latitude : 32 - 45 N
+! longitude: 130-145 E
+! depth : 0 - 500 km
+! The deepest Moho beneath Japan is 40 km
+ double precision,parameter :: LAT_MAX = 45.d0
+ double precision,parameter :: LAT_MIN = 32.d0
+ double precision,parameter :: LON_MAX = 145.d0
+ double precision,parameter :: LON_MIN = 130.d0
+ double precision,parameter :: DEP_MAX = 500.d0
+
+! use sedimentary layers of crust 2.0
+ logical, parameter :: INCLUDE_SEDIMENTS_CRUST = .true.
+
+! number of points per degree in smoothed crust2.0
+ integer, parameter :: NFACTOR_CRUST = 2
+ integer, parameter :: NLON_CRUST = 180 * NFACTOR_CRUST
+ integer, parameter :: NLAT_CRUST = 90 * NFACTOR_CRUST
+
+! to inflate the central cube (set to 0.d0 for a non-inflated cube)
+ double precision, parameter :: CENTRAL_CUBE_INFLATE_FACTOR = 0.41d0
+
+! for the stretching of crustal elements in the case of 3D models
+ double precision, parameter :: MAX_RATIO_CRUST_STRETCHING = 0.6d0
+
+! to suppress the crustal layers (replaced by an extension of the mantle: R_EARTH is not modified, but no more crustal doubling)
+ logical, parameter :: SUPPRESS_CRUSTAL_MESH = .false.
+
+! to add a fourth doubling at the bottom of the outer core
+ logical, parameter :: ADD_4TH_DOUBLING = .false.
+
+! parameters to cut the doubling brick
+
+! this to cut the superbrick: 3 possibilities, 4 cases max / possibility
+! three possibilities: (cut in xi and eta) or (cut in xi) or (cut in eta)
+! case 1: (ximin and etamin) or ximin or etamin
+! case 2: (ximin and etamax) or ximax or etamax
+! case 3: ximax and etamin
+! case 4: ximax and etamax
+ integer, parameter :: NB_CUT_CASE = 4
+
+! corner 1: ximin and etamin
+! corner 2: ximax and etamin
+! corner 3: ximax and etamax
+! corner 4: ximin and etamax
+ integer, parameter :: NB_SQUARE_CORNERS = 4
+
+! two possibilities: xi or eta
+! face 1: ximin or etamin
+! face 2: ximax or etamax
+ integer, parameter :: NB_SQUARE_EDGES_ONEDIR = 2
+
+! this for the geometry of the basic doubling brick
+ integer, parameter :: NSPEC_DOUBLING_BASICBRICK = 8
+ integer, parameter :: NGLOB_DOUBLING_BASICBRICK = 27
+
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/read_compute_parameters.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/read_compute_parameters.F90 2008-11-01 23:15:42 UTC (rev 13216)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/read_compute_parameters.F90 2008-11-02 00:24:51 UTC (rev 13217)
@@ -1162,6 +1162,8 @@
SAVE_ALL_SEISMOS_IN_ONE_FILE = .true.
USE_BINARY_FOR_LARGE_FILE = .true.
RECEIVERS_CAN_BE_BURIED = .false.
+ MOVIE_SURFACE = .false.
+ MOVIE_VOLUME = .false.
! when using restart files, make sure the two intervals coincide
if(USE_RESTART_FILES) then
@@ -1169,7 +1171,6 @@
SAVE_ALL_SEISMOS_IN_ONE_FILE = .false.
USE_BINARY_FOR_LARGE_FILE = .false.
endif
-
endif
if(err_occurred() /= 0) stop 'an error occurred while reading the parameter file'
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/specfem3D.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/specfem3D.F90 2008-11-01 23:15:42 UTC (rev 13216)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/specfem3D.F90 2008-11-02 00:24:51 UTC (rev 13217)
@@ -77,6 +77,15 @@
! include values created by the mesher
include "values_from_mesher.h"
+! to save movie frames
+ integer nmovie_points
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: &
+ store_val_x,store_val_y,store_val_z, &
+ store_val_ux,store_val_uy,store_val_uz
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: &
+ store_val_x_all,store_val_y_all,store_val_z_all, &
+ store_val_ux_all,store_val_uy_all,store_val_uz_all
+
! for non blocking communications
integer :: icall,iphase
real :: percentage_edge
@@ -976,6 +985,16 @@
if(minval(t_cmt) /= 0.) call exit_MPI(myrank,'one t_cmt must be zero, others must be positive')
+! filter source time function by Gaussian with hdur = HDUR_MOVIE when outputing movies or shakemaps
+ if (MOVIE_SURFACE .or. MOVIE_VOLUME ) then
+ hdur = sqrt(hdur**2 + HDUR_MOVIE**2)
+ if(myrank == 0) then
+ write(IMAIN,*)
+ write(IMAIN,*) 'Each source is being convolved with HDUR_MOVIE = ',HDUR_MOVIE
+ write(IMAIN,*)
+ endif
+ endif
+
! convert the half duration for triangle STF to the one for gaussian STF
hdur_gaussian = hdur/SOURCE_DECAY_MIMIC_TRIANGLE
@@ -1595,6 +1614,25 @@
d_ln_density_dr_table(int_radius) = drhodr/rho
enddo
+! allocate files to save movies
+ if(MOVIE_SURFACE) then
+ nmovie_points = NGLLX * NGLLY * NSPEC2D_TOP(IREGION_CRUST_MANTLE)
+
+ allocate(store_val_x(nmovie_points))
+ allocate(store_val_y(nmovie_points))
+ allocate(store_val_z(nmovie_points))
+ allocate(store_val_ux(nmovie_points))
+ allocate(store_val_uy(nmovie_points))
+ allocate(store_val_uz(nmovie_points))
+
+ allocate(store_val_x_all(nmovie_points,0:NPROCTOT-1))
+ allocate(store_val_y_all(nmovie_points,0:NPROCTOT-1))
+ allocate(store_val_z_all(nmovie_points,0:NPROCTOT-1))
+ allocate(store_val_ux_all(nmovie_points,0:NPROCTOT-1))
+ allocate(store_val_uy_all(nmovie_points,0:NPROCTOT-1))
+ allocate(store_val_uz_all(nmovie_points,0:NPROCTOT-1))
+ endif
+
if(myrank == 0) then
write(IMAIN,*)
write(IMAIN,*) ' time step: ',sngl(DT),' s'
@@ -2843,6 +2881,59 @@
endif
+! save movie on surface
+ if(MOVIE_SURFACE .and. mod(it,NTSTEP_BETWEEN_FRAMES) == 0) then
+
+! save velocity here to avoid static offset on displacement for movies
+
+! get coordinates of surface mesh and surface displacement
+ ipoin = 0
+ do ispec2D = 1,NSPEC2D_TOP(IREGION_CRUST_MANTLE)
+ ispec = ibelm_top_crust_mantle(ispec2D)
+ k = NGLLZ
+
+! loop on all the points inside the element
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+ ipoin = ipoin + 1
+ iglob = ibool_crust_mantle(i,j,k,ispec)
+ store_val_x(ipoin) = xstore_crust_mantle(iglob)
+ store_val_y(ipoin) = ystore_crust_mantle(iglob)
+ store_val_z(ipoin) = zstore_crust_mantle(iglob)
+ store_val_ux(ipoin) = veloc_crust_mantle(1,iglob)*scale_veloc
+ store_val_uy(ipoin) = veloc_crust_mantle(2,iglob)*scale_veloc
+ store_val_uz(ipoin) = veloc_crust_mantle(3,iglob)*scale_veloc
+ enddo
+ enddo
+
+ enddo
+
+! gather info on master proc
+ ispec = nmovie_points
+#ifdef USE_MPI
+ call MPI_GATHER(store_val_x,ispec,CUSTOM_MPI_TYPE,store_val_x_all,ispec,CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
+ call MPI_GATHER(store_val_y,ispec,CUSTOM_MPI_TYPE,store_val_y_all,ispec,CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
+ call MPI_GATHER(store_val_z,ispec,CUSTOM_MPI_TYPE,store_val_z_all,ispec,CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
+ call MPI_GATHER(store_val_ux,ispec,CUSTOM_MPI_TYPE,store_val_ux_all,ispec,CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
+ call MPI_GATHER(store_val_uy,ispec,CUSTOM_MPI_TYPE,store_val_uy_all,ispec,CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
+ call MPI_GATHER(store_val_uz,ispec,CUSTOM_MPI_TYPE,store_val_uz_all,ispec,CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
+#endif
+
+! save movie data to disk in home directory
+ if(myrank == 0) then
+ write(outputname,"('/scratch/komatits/moviedata',i6.6)") it
+ open(unit=IOUT,file=outputname,status='unknown',form='unformatted',action='write')
+ write(IOUT) store_val_x_all
+ write(IOUT) store_val_y_all
+ write(IOUT) store_val_z_all
+ write(IOUT) store_val_ux_all
+ write(IOUT) store_val_uy_all
+ write(IOUT) store_val_uz_all
+ close(IOUT)
+ endif
+
+ endif
+
!---- end of time iteration loop
!
enddo ! end of main time loop
More information about the CIG-COMMITS
mailing list