[cig-commits] [commit] : Merged r959:1028 from trunk to v1.4 tag. (bd75705)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Nov 14 20:15:41 PST 2013
Repository : ssh://geoshell/specfem3d
On branch :
Link : https://github.com/geodynamics/specfem2d/compare/1e201257d91c794056b990a43329e05d04f77454...0000000000000000000000000000000000000000
>---------------------------------------------------------------
commit bd75705d977d0ed76049cf6ea4f92050487d93cd
Author: Leif Strand <leif at geodynamics.org>
Date: Tue Oct 24 22:31:34 2006 +0000
Merged r959:1028 from trunk to v1.4 tag.
>---------------------------------------------------------------
bd75705d977d0ed76049cf6ea4f92050487d93cd
Makefile | 10 +-
OUTPUT_FILES/values_from_mesher.h | 51 ---
constants.h | 2 +
create_movie_GMT_basin.f90 | 748 ++++++++++++++++++++++++++++++++++++++
locate_receivers.f90 | 1 +
locate_source.f90 | 1 +
meshfem3D.f90 | 18 +-
specfem3D.f90 | 40 +-
8 files changed, 812 insertions(+), 59 deletions(-)
diff --git a/Makefile b/Makefile
index b7f6dd5..564df32 100644
--- a/Makefile
+++ b/Makefile
@@ -168,7 +168,7 @@ TURN_MPI_ON = -DUSE_MPI
#C_PREPROCESSOR = -qsuffix=cpp=f90
#TURN_MPI_ON = -WF,-DUSE_MPI
-default: meshfem3D create_movie_AVS_DX combine_AVS_DX check_mesh_quality_AVS_DX check_buffers_2D convolve_source_timefunction
+default: meshfem3D create_movie_AVS_DX combine_AVS_DX check_mesh_quality_AVS_DX check_buffers_2D convolve_source_timefunction create_movie_GMT_basin
all: clean default
@@ -515,8 +515,11 @@ check_buffers_2D: constants.h $O/check_buffers_2D.o \
combine_paraview_data: constants.h $O/combine_paraview_data.o $O/write_c_binary.o
${F90} $(FLAGS_CHECK) -o xcombine_paraview_data $O/combine_paraview_data.o $O/write_c_binary.o
+create_movie_GMT_basin : $O/create_movie_GMT_basin.o $O/read_parameter_file.o $O/compute_parameters.o $O/utm_geo.o $O/read_value_parameters.o $O/get_value_parameters.o
+ $(F90) $(FLAGS) -o xcreate_movie_GMT_basin $O/create_movie_GMT_basin.o $O/read_parameter_file.o $O/compute_parameters.o $O/utm_geo.o $O/read_value_parameters.o $O/get_value_parameters.o
+
clean:
- rm -f $O/* *.o *.gnu OUTPUT_FILES/timestamp* OUTPUT_FILES/starttime*txt work.pc* xmeshfem3D xspecfem3D xcombine_AVS_DX xcheck_mesh_quality_AVS_DX xcheck_buffers_2D xconvolve_source_timefunction xcreate_header_file xcreate_movie_AVS_DX xcombine_paraview_data
+ rm -f $O/* *.o *.gnu OUTPUT_FILES/timestamp* OUTPUT_FILES/starttime*txt work.pc* xmeshfem3D xspecfem3D xcombine_AVS_DX xcheck_mesh_quality_AVS_DX xcheck_buffers_2D xconvolve_source_timefunction xcreate_header_file xcreate_movie_AVS_DX xcombine_paraview_data xcreate_movie_GMT_basin
####
#### rule to build each .o file below
@@ -751,6 +754,9 @@ $O/get_attenuation_model.o: constants.h get_attenuation_model.f90
$O/combine_paraview_data.o: constants.h combine_paraview_data.f90
${F90} $(FLAGS_CHECK) -c -o $O/combine_paraview_data.o combine_paraview_data.f90
+$O/create_movie_GMT_basin.o: constants.h create_movie_GMT_basin.f90
+ ${F90} $(FLAGS_CHECK) -c -o $O/create_movie_GMT_basin.o create_movie_GMT_basin.f90
+
###
### C files below
###
diff --git a/OUTPUT_FILES/values_from_mesher.h b/OUTPUT_FILES/values_from_mesher.h
deleted file mode 100644
index 2090aea..0000000
--- a/OUTPUT_FILES/values_from_mesher.h
+++ /dev/null
@@ -1,51 +0,0 @@
-
- !
- ! this is the parameter file for static compilation of the solver
- !
- ! mesh statistics:
- ! ---------------
- !
- ! number of processors = 36
- !
- ! number of ES nodes = 4.500000
- ! percentage of total 640 ES nodes = 0.7031250 %
- ! total memory available on these ES nodes (Gb) = 72.00000
- !
- ! max points per processor = max vector length = 703989
- ! min vector length = 25
- ! min critical vector length = 75
- !
- ! on ES and SX-5, make sure "loopcnt=" parameter
- ! in Makefile is greater than 703989
- !
- ! total elements per AB slice = 10512
- ! total points per AB slice = 703989
- !
- ! total for full mesh:
- ! -------------------
- !
- ! exact total number of spectral elements in entire mesh =
- ! 378432
- ! approximate total number of points in entire mesh =
- ! 25343604.0000000
- ! approximate total number of degrees of freedom in entire mesh =
- ! 76030812.0000000
- !
- ! resolution of the mesh at the surface:
- ! -------------------------------------
- !
- ! spectral elements along X = 288
- ! spectral elements along Y = 288
- ! GLL points along X = 1153
- ! GLL points along Y = 1153
- ! average distance between points along X in m = 239.0178
- ! average distance between points along Y in m = 287.6366
- !
-
- integer, parameter :: NSPEC_AB = 10512
- integer, parameter :: NGLOB_AB = 703989
-
- integer, parameter :: NSPEC_ATTENUATION = 1
- logical, parameter :: ATTENUATION_VAL = .false.
- logical, parameter :: ANISOTROPY_VAL = .false.
-
diff --git a/constants.h b/constants.h
index ef86b10..462b35b 100644
--- a/constants.h
+++ b/constants.h
@@ -50,6 +50,8 @@
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
diff --git a/create_movie_GMT_basin.f90 b/create_movie_GMT_basin.f90
new file mode 100644
index 0000000..6bb1dc9
--- /dev/null
+++ b/create_movie_GMT_basin.f90
@@ -0,0 +1,748 @@
+!=====================================================================
+!
+! 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=100) 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.5_CUSTOM_REAL
+
+ integer it,it1,it2,ivalue,ispec
+ integer nframes,iframe,inorm
+ integer ibool_number
+
+ logical 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
+
+ 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_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,HDUR_MOVIE
+ 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_MESH_FILES,PRINT_SOURCE_TIME_FUNCTION,USE_REGULAR_MESH,SUPPRESS_UTM_PROJECTION
+
+ logical MOVIE_SURFACE,MOVIE_VOLUME,CREATE_SHAKEMAP,SAVE_DISPLACEMENT,USE_HIGHRES_FOR_MOVIES
+ logical SAVE_FORWARD
+ integer NTSTEP_BETWEEN_FRAMES,NTSTEP_BETWEEN_OUTPUT_INFO,SIMULATION_TYPE
+
+ character(len=150) LOCAL_PATH,MODEL
+ ! 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,start_frame)
+ call getarg(3,end_frame)
+ call getarg(4,output_file_prefix)
+ if (trim(movie_data_prefix) == '' &
+ .or. trim(start_frame) == '' .or. trim(end_frame) == '' &
+ .or. trim(output_file_prefix) == '') &
+ stop 'Usage: xcreate_movie_GMT movie_data_prefix 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(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 the movie'
+ else
+ print *,'Vertical velocity will be shown in the movie'
+ endif
+ print *
+
+ 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_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)
+
+! print *, store_val_x(1,0), store_val_y(1,0)
+! 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
+
+! low resolution movie/shakemap
+ 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)
+ 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 *
+
+! comment out the following to try to retain the original value of the shakemap
+! 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
+ !! choose ascii or binary output file format
+ !write(11,rec=igmt) long,lat,field_display(ilocnum+ieoff)
+ write (11,'(2f10.3,e17.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
+
+
+ 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
+
diff --git a/locate_receivers.f90 b/locate_receivers.f90
index 5bbafc3..ce317c7 100644
--- a/locate_receivers.f90
+++ b/locate_receivers.f90
@@ -273,6 +273,7 @@
x_target(irec) = stutm_x(irec)
y_target(irec) = stutm_y(irec)
z_target(irec) = elevation(irec) - stbur(irec)
+ if (myrank == 0) write(IOVTK,*) x_target(irec), y_target(irec), z_target(irec)
! examine top of the elements only (receivers always at the surface)
! k = NGLLZ
diff --git a/locate_source.f90 b/locate_source.f90
index d811b2d..23b4bbd 100644
--- a/locate_source.f90
+++ b/locate_source.f90
@@ -234,6 +234,7 @@
x_target_source = utm_x_source(isource)
y_target_source = utm_y_source(isource)
z_target_source = - depth(isource)*1000.0d0 + elevation(isource)
+ if(myrank == 0) write(IOVTK,*) x_target_source,y_target_source,z_target_source
! set distance to huge initial value
distmin = HUGEVAL
diff --git a/meshfem3D.f90 b/meshfem3D.f90
index 8ce98bf..ec65759 100644
--- a/meshfem3D.f90
+++ b/meshfem3D.f90
@@ -65,7 +65,7 @@
! If you use this code for your own research, please send an email
! to Jeroen Tromp <jtromp AT caltech.edu> for information, and cite:
!
-! @article{KoLiTrSuStSh04,
+! @ARTICLE{KoLiTrSuStSh04,
! author={Dimitri Komatitsch and Qinya Liu and Jeroen Tromp and Peter S\"{u}ss
! and Christiane Stidham and John H. Shaw},
! year=2004,
@@ -73,15 +73,27 @@
! based upon the Spectral-Element Method},
! journal={Bull. Seism. Soc. Am.},
! volume=94,
+! number=1,
! pages={187-206}}
!
-! @article{KoTr99,
+! @ARTICLE{KoTr99,
! author={D. Komatitsch and J. Tromp},
! year=1999,
! title={Introduction to the spectral-element method for 3-{D} seismic wave propagation},
! journal={Geophys. J. Int.},
! volume=139,
-! pages={806-822}}
+! number=3,
+! pages={806-822},
+! doi={10.1046/j.1365-246x.1999.00967.x}}
+!
+! @ARTICLE{KoVi98,
+! author={D. Komatitsch and J. P. Vilotte},
+! title={The spectral-element method: an efficient tool to simulate the seismic response of 2{D} and 3{D} geological structures},
+! journal={Bull. Seismol. Soc. Am.},
+! year=1998,
+! volume=88,
+! number=2,
+! pages={368-392}}
!
! Reference frame - convention:
! ----------------------------
diff --git a/specfem3D.f90 b/specfem3D.f90
index 2a71386..d6176ab 100644
--- a/specfem3D.f90
+++ b/specfem3D.f90
@@ -68,7 +68,7 @@
! If you use this code for your own research, please send an email
! to Jeroen Tromp <jtromp AT caltech.edu> for information, and cite:
!
-! @article{KoLiTrSuStSh04,
+! @ARTICLE{KoLiTrSuStSh04,
! author={Dimitri Komatitsch and Qinya Liu and Jeroen Tromp and Peter S\"{u}ss
! and Christiane Stidham and John H. Shaw},
! year=2004,
@@ -76,15 +76,36 @@
! based upon the Spectral-Element Method},
! journal={Bull. Seism. Soc. Am.},
! volume=94,
+! number=1,
! pages={187-206}}
!
-! @article{KoTr99,
+! @ARTICLE{KoTr99,
! author={D. Komatitsch and J. Tromp},
! year=1999,
! title={Introduction to the spectral-element method for 3-{D} seismic wave propagation},
! journal={Geophys. J. Int.},
! volume=139,
-! pages={806-822}}
+! number=3,
+! pages={806-822},
+! doi={10.1046/j.1365-246x.1999.00967.x}}
+!
+! @ARTICLE{KoVi98,
+! author={D. Komatitsch and J. P. Vilotte},
+! title={The spectral-element method: an efficient tool to simulate the seismic response of 2{D} and 3{D} geological structures},
+! journal={Bull. Seismol. Soc. Am.},
+! year=1998,
+! volume=88,
+! number=2,
+! pages={368-392}}
+!
+! If you use the kernel capabilities of the code, please cite
+!
+! @ARTICLE{Liu06a,
+! AUTHOR = {Q. Liu and J. Tromp},
+! JOURNAL =bssa,
+! TITLE = {{Finite-frequency kernels based upon adjoint methods}},
+! NOTE = {accepted},
+! YEAR = {2006}
!
! Reference frame - convention:
! ----------------------------
@@ -945,6 +966,16 @@
if(nrec < 1) call exit_MPI(myrank,'need at least one receiver')
+! write source and receiver VTK files for Paraview
+ if (myrank == 0) then
+ open(IOVTK,file=trim(OUTPUT_FILES)//'/sr.vtk',status='unknown')
+ write(IOVTK,'(a)') '# vtk DataFile Version 2.0'
+ write(IOVTK,'(a)') 'Source and Receiver VTK file'
+ write(IOVTK,'(a)') 'ASCII'
+ write(IOVTK,'(a)') 'DATASET UNSTRUCTURED_GRID'
+ write(IOVTK, '(a,i6,a)') 'POINTS ', NSOURCES+nrec, ' float'
+ endif
+
! allocate memory for receiver arrays
allocate(islice_selected_rec(nrec))
allocate(ispec_selected_rec(nrec))
@@ -1090,6 +1121,9 @@
nrec_tot_found = nrec_local
#endif
if(myrank == 0) then
+
+ close(IOVTK)
+
write(IMAIN,*)
write(IMAIN,*) 'Total number of samples for seismograms = ',NSTEP
write(IMAIN,*)
More information about the CIG-COMMITS
mailing list