[cig-commits] r14820 - in seismo/3D/SPECFEM3D_SESAME: tags/v1.4.4_last_BASIN trunk trunk/USER_MANUAL
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Wed Apr 29 15:50:18 PDT 2009
Author: dkomati1
Date: 2009-04-29 15:50:18 -0700 (Wed, 29 Apr 2009)
New Revision: 14820
Added:
seismo/3D/SPECFEM3D_SESAME/trunk/create_movie_shakemap_AVS_DX_GMT.f90
Removed:
seismo/3D/SPECFEM3D_SESAME/trunk/create_movie_AVS_DX.f90
Modified:
seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/create_movie_shakemap_AVS_DX_GMT.f90
seismo/3D/SPECFEM3D_SESAME/trunk/USER_MANUAL/manual_SPECFEM3D.tex
Log:
renamed create_movie_AVS_DX.f90 to create_movie_shakemap_AVS_DX_GMT.f90.
Updated the acknowledgements section of the manual.
Modified: seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/create_movie_shakemap_AVS_DX_GMT.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/create_movie_shakemap_AVS_DX_GMT.f90 2009-04-29 22:43:27 UTC (rev 14819)
+++ seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/create_movie_shakemap_AVS_DX_GMT.f90 2009-04-29 22:50:18 UTC (rev 14820)
@@ -24,11 +24,12 @@
!=====================================================================
!
-!--- create a movie of vertical component of surface displacement or velocity
-!--- in AVS or OpenDX format
+!--- create a movie of the vertical component of surface displacement or velocity
+!--- or a ShakeMap(R) (i.e. map of the maximum absolute value of the two horizontal components
+!--- of the velocity vector) in AVS, OpenDX or GMT format
!
- program create_movie_AVS_DX
+ program create_movie_shakemap
implicit none
@@ -46,7 +47,7 @@
integer iformat,nframes,iframe,inumber,inorm,iscaling_shake
integer ibool_number,ibool_number1,ibool_number2,ibool_number3,ibool_number4
- logical USE_OPENDX,USE_AVS,USE_GMT,UNIQUE_FILE,plot_shaking_map
+ logical USE_OPENDX,USE_AVS,USE_GMT,plot_shaking_map
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: x,y,z,display
real(kind=CUSTOM_REAL) xcoord,ycoord,zcoord
@@ -164,34 +165,28 @@
allocate(store_val_uz(ilocnum,0:NPROC-1))
print *,'1 = create files in OpenDX format'
- print *,'2 = create files in AVS UCD format with individual files'
- print *,'3 = create files in AVS UCD format with one time-dependent file'
- print *,'4 = create files in GMT xyz Ascii long/lat/Uz format'
+ print *,'2 = create files in AVS UCD format'
+ print *,'3 = create files in GMT xyz Ascii long/lat/Uz format'
print *,'any other value = exit'
print *
+
print *,'enter value:'
read(5,*) iformat
- if(iformat<1 .or. iformat>4) stop 'exiting...'
+
+ if(iformat<1 .or. iformat>3) stop 'exiting...'
+
if(iformat == 1) then
USE_OPENDX = .true.
USE_AVS = .false.
USE_GMT = .false.
- UNIQUE_FILE = .false.
else if(iformat == 2) then
USE_OPENDX = .false.
USE_AVS = .true.
USE_GMT = .false.
- UNIQUE_FILE = .false.
- else if(iformat == 3) then
- USE_OPENDX = .false.
- USE_AVS = .true.
- USE_GMT = .false.
- UNIQUE_FILE = .true.
else
USE_OPENDX = .false.
USE_AVS = .false.
USE_GMT = .true.
- UNIQUE_FILE = .false.
endif
if(.not. USE_GMT) then
@@ -577,7 +572,6 @@
open(unit=11,file=trim(OUTPUT_FILES)//outputname,status='unknown')
write(11,*) 'object 1 class array type float rank 1 shape 3 items ',nglob,' data follows'
else if(USE_AVS) then
- if(UNIQUE_FILE) stop 'cannot use unique file AVS option for shaking map'
write(outputname,"('/AVS_shaking_map.inp')")
open(unit=11,file=trim(OUTPUT_FILES)//outputname,status='unknown')
write(11,*) nglob,' ',nspectot_AVS_max,' 1 0 0'
@@ -595,17 +589,9 @@
open(unit=11,file=trim(OUTPUT_FILES)//outputname,status='unknown')
write(11,*) 'object 1 class array type float rank 1 shape 3 items ',nglob,' data follows'
else if(USE_AVS) then
- if(UNIQUE_FILE .and. iframe == 1) then
- open(unit=11,file=trim(OUTPUT_FILES)//'/AVS_movie_all.inp',status='unknown')
- write(11,*) nframes
- write(11,*) 'data'
- write(11,"('step',i1,' image',i1)") 1,1
- write(11,*) nglob,' ',nspectot_AVS_max
- else if(.not. UNIQUE_FILE) then
- write(outputname,"('/AVS_movie_',i6.6,'.inp')") ivalue
- open(unit=11,file=trim(OUTPUT_FILES)//outputname,status='unknown')
- write(11,*) nglob,' ',nspectot_AVS_max,' 1 0 0'
- endif
+ write(outputname,"('/AVS_movie_',i6.6,'.inp')") ivalue
+ open(unit=11,file=trim(OUTPUT_FILES)//outputname,status='unknown')
+ write(11,*) nglob,' ',nspectot_AVS_max,' 1 0 0'
else if(USE_GMT) then
write(outputname,"('/gmt_movie_',i6.6,'.xyz')") ivalue
open(unit=11,file=trim(OUTPUT_FILES)//outputname,status='unknown')
@@ -635,9 +621,6 @@
else
-! if unique file, output geometry only once
- if(.not. UNIQUE_FILE .or. iframe == 1) then
-
! output list of points
mask_point = .false.
ipoin = 0
@@ -681,28 +664,11 @@
endif
enddo
- endif
-
if(USE_OPENDX) then
write(11,*) 'attribute "element type" string "quads"'
write(11,*) 'attribute "ref" string "positions"'
write(11,*) 'object 3 class array type float rank 0 items ',nglob,' data follows'
else
- if(UNIQUE_FILE) then
-! step number for AVS multistep file
- if(iframe > 1) then
- if(iframe < 10) then
- write(11,"('step',i1,' image',i1)") iframe,iframe
- else if(iframe < 100) then
- write(11,"('step',i2,' image',i2)") iframe,iframe
- else if(iframe < 1000) then
- write(11,"('step',i3,' image',i3)") iframe,iframe
- else
- write(11,"('step',i4,' image',i4)") iframe,iframe
- endif
- endif
- write(11,*) '1 0'
- endif
! dummy text for labels
write(11,*) '1 1'
write(11,*) 'a, b'
@@ -749,14 +715,12 @@
! end of test for GMT format
endif
- if(.not. UNIQUE_FILE) close(11)
+ close(11)
! end of loop and test on all the time steps for all the movie images
endif
enddo
- if(UNIQUE_FILE) close(11)
-
print *
print *,'done creating movie or shaking map'
print *
@@ -789,7 +753,7 @@
deallocate(display)
endif
- end program create_movie_AVS_DX
+ end program create_movie_shakemap
!
!=====================================================================
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/USER_MANUAL/manual_SPECFEM3D.tex
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/USER_MANUAL/manual_SPECFEM3D.tex 2009-04-29 22:43:27 UTC (rev 14819)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/USER_MANUAL/manual_SPECFEM3D.tex 2009-04-29 22:50:18 UTC (rev 14820)
@@ -1553,7 +1553,7 @@
To report bugs or suggest improvements to the code, please send an
e-mail to the CIG Computational Seismology Mailing List \url{cig-seismo at geodynamics.org}
-or Jeroen Tromp \url{jtromp-AT-gps.caltech.edu}, and/or use our online
+or Jeroen Tromp \url{jtromp-AT-princeton.edu}, and/or use our online
bug tracking system Roundup \url{www.geodynamics.org/roundup}.
@@ -1566,10 +1566,10 @@
clarity).
The Gauss-Lobatto-Legendre subroutines in \texttt{gll\_library.f90}
-are based in part on software libraries from Massachusetts Institute
-of Technology, Department of Mechanical Engineering, Cambridge, Massachusetts.
+are based in part on software libraries from the Massachusetts Institute
+of Technology, Department of Mechanical Engineering (Cambridge, Massachusetts, USA).
The non-structured global numbering software was provided by Paul
-F. Fischer.
+F. Fischer (Brown University, Providence, Rhode Island, USA, now at Argonne National Laboratory, USA).
OpenDX \url{http://www.opendx.org} is open-source based on IBM Data
Explorer, AVS \url{http://www.avs.com} is a trademark of Advanced
@@ -1580,10 +1580,9 @@
Komatitsch, Jeroen Tromp, and Qinya Liu. The following individuals
(listed in alphabetical order) have also contributed to the development
or improvement of the source code: Min Chen, Vala Hjörleifsdóttir,
-Jes\'us Labarta, Brian Savage,
-and Leif Strand. The following individuals (listed in alphabetical
+Jes\'us Labarta, and Leif Strand. The following individuals (listed in alphabetical
order) contributed to this manual: Min Chen, Vala Hjörleifsdóttir,
-Sue Kientz, Dimitri Komatitsch, Qinya Liu, Alessia Maggi, Brian Savage,
+Sue Kientz, Dimitri Komatitsch, Qinya Liu, Alessia Maggi,
Carl Tape, and Jeroen Tromp. The manual's cover graphic was created
by Santiago Lombeyda from Caltech's Center for Advanced Computing Research (CACR) \url{http://www.cacr.caltech.edu/}.
Deleted: seismo/3D/SPECFEM3D_SESAME/trunk/create_movie_AVS_DX.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/create_movie_AVS_DX.f90 2009-04-29 22:43:27 UTC (rev 14819)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/create_movie_AVS_DX.f90 2009-04-29 22:50:18 UTC (rev 14820)
@@ -1,1029 +0,0 @@
-!=====================================================================
-!
-! 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.
-!
-!=====================================================================
-
-!
-!--- create a movie of vertical component of surface displacement or velocity
-!--- in AVS or OpenDX format
-!
-
- program create_movie_AVS_DX
-
- implicit none
-
- include "constants.h"
-
-! threshold in percent of the maximum below which we cut the amplitude
- logical, parameter :: APPLY_THRESHOLD = .true.
- 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 = .true.
- real(kind=CUSTOM_REAL), parameter :: POWER_SCALING = 0.13_CUSTOM_REAL
-
- integer it,it1,it2,ivalue,nspectot_AVS_max,ispec
- integer iformat,nframes,iframe,inumber,inorm,iscaling_shake
- integer ibool_number,ibool_number1,ibool_number2,ibool_number3,ibool_number4
-
- logical USE_OPENDX,USE_AVS,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
-
-! 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
-
- logical MOVIE_SURFACE,MOVIE_VOLUME,CREATE_SHAKEMAP,SAVE_DISPLACEMENT, &
- USE_HIGHRES_FOR_MOVIES,SUPPRESS_UTM_PROJECTION,USE_REGULAR_MESH
- integer NTSTEP_BETWEEN_FRAMES,NTSTEP_BETWEEN_OUTPUT_INFO
-
- 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,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 ANISOTROPY,SAVE_MESH_FILES,PRINT_SOURCE_TIME_FUNCTION
-
- character(len=150) OUTPUT_FILES,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
-
-!!!! NL NL for external meshes
- real(kind=CUSTOM_REAL), parameter :: RADIUS_TO_MUTE = 1000._CUSTOM_REAL
- logical, parameter :: MUTE_SOURCE = .true.
- real(kind=CUSTOM_REAL), parameter :: X_SOURCE_EXT_MESH = -9023.021484375
- real(kind=CUSTOM_REAL), parameter :: Y_SOURCE_EXT_MESH = 6123.611328125
- real(kind=CUSTOM_REAL), parameter :: Z_SOURCE_EXT_MESH = 17.96331405639648
- integer, parameter :: NSPEC_SURFACE_EXT_MESH = 15808*4
-
-!!!! NL NL
-
-! ************** PROGRAM STARTS HERE **************
-
- print *
- print *,'Recombining all movie frames to create a movie'
- print *
-
- print *
- print *,'reading parameter file'
- print *
-
-! 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)
-
-! get the base pathname for output files
- call get_value_string(OUTPUT_FILES, 'OUTPUT_FILES', 'OUTPUT_FILES')
-
- print *
- print *,'There are ',NPROC,' slices numbered from 0 to ',NPROC-1
- 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_EXTERNAL_MESH) then
- NPROC = 1
- if(USE_HIGHRES_FOR_MOVIES) then
- ilocnum = NSPEC_SURFACE_EXT_MESH*NGLLSQUARE
- else
- ilocnum = NSPEC_SURFACE_EXT_MESH*NGNOD2D_AVS_DX
- endif
- else
- 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
- 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))
-
-
- print *,'1 = create files in OpenDX format'
- print *,'2 = create files in AVS UCD format with individual files'
- print *,'3 = create files in AVS UCD format with one time-dependent file'
- print *,'4 = create files in GMT xyz Ascii long/lat/Uz format'
- print *,'any other value = exit'
- print *
- print *,'enter value:'
- read(5,*) iformat
- if(iformat<1 .or. iformat>4) stop 'exiting...'
- if(iformat == 1) then
- USE_OPENDX = .true.
- USE_AVS = .false.
- UNIQUE_FILE = .false.
- else if(iformat == 2) then
- USE_OPENDX = .false.
- USE_AVS = .true.
- UNIQUE_FILE = .false.
- else if(iformat == 3) then
- USE_OPENDX = .false.
- USE_AVS = .true.
- UNIQUE_FILE = .true.
- else
- USE_OPENDX = .false.
- USE_AVS = .false.
- UNIQUE_FILE = .false.
- endif
-
- print *,'movie frames have been saved every ',NTSTEP_BETWEEN_FRAMES,' time steps'
- print *
-
- plot_shaking_map = .false.
- print *,'enter first time step of movie (e.g. 1, enter -1 for shaking map)'
- read(5,*) it1
- if(it1 == -1) plot_shaking_map = .true.
-
- if(.not. plot_shaking_map) then
-
- print *,'enter last time step of movie (e.g. ',NSTEP,')'
- read(5,*) it2
-
- print *
- print *,'1 = define file names using frame number'
- print *,'2 = define file names using time step number'
- print *,'any other value = exit'
- print *
- print *,'enter value:'
- read(5,*) inumber
- if(inumber<1 .or. inumber>2) stop 'exiting...'
-
- print *
- print *,'looping from ',it1,' to ',it2,' every ',NTSTEP_BETWEEN_FRAMES,' time steps'
-
-! count number of movie frames
- nframes = 0
- do it = it1,it2
- if(mod(it,NTSTEP_BETWEEN_FRAMES) == 0) nframes = nframes + 1
- enddo
- print *
- print *,'total number of frames will be ',nframes
- if(nframes == 0) stop 'null number of frames'
-
- else
-
-! only one frame if shaking map
- nframes = 1
- it1 = 1
- it2 = 1
- endif
-
- iscaling_shake = 0
- if(plot_shaking_map) then
- print *
- print *,'norm to display in shaking map:'
- print *,'1=displacement 2=velocity 3=acceleration'
- print *
- read(5,*) inorm
- if(inorm < 1 .or. inorm > 3) stop 'incorrect value of inorm'
-
- print *
- print *,'apply non-linear scaling to shaking map:'
- print *,'1=non-linear 2=no scaling'
- print *
- read(5,*) iscaling_shake
- if(iscaling_shake < 1 .or. iscaling_shake > 2) stop 'incorrect value of iscaling_shake'
- endif
-
-! define the total number of elements at the surface
- if (USE_EXTERNAL_MESH) then
- if(USE_HIGHRES_FOR_MOVIES) then
- nspectot_AVS_max = NSPEC_SURFACE_EXT_MESH * (NGLLX-1) * (NGLLY-1)
- else
- nspectot_AVS_max = NSPEC_SURFACE_EXT_MESH
- endif
- else
- 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
- 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))
-
- print *
- if(APPLY_THRESHOLD .and. .not. plot_shaking_map) print *,'Will apply a threshold to amplitude below ',100.*THRESHOLD,' %'
- if(NONLINEAR_SCALING .and. (.not. plot_shaking_map .or. iscaling_shake == 1)) &
- print *,'Will apply a non linear scaling with coef ',POWER_SCALING
-
-! --------------------------------------
-
- 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
- if(mod(it,NTSTEP_BETWEEN_FRAMES) == 0 .or. plot_shaking_map) then
-
- iframe = iframe + 1
-
- print *
- if(plot_shaking_map) then
- print *,'reading shaking map snapshot'
- else
- print *,'reading snapshot time step ',it,' out of ',NSTEP
- endif
- print *
-
-! read all the elements from the same file
- if(plot_shaking_map) then
- write(outputname,"('/shakingdata')")
- else
- write(outputname,"('/moviedata',i6.6)") it
- endif
- open(unit=IOUT,file=trim(OUTPUT_FILES)//outputname,status='old',action='read',form='unformatted')
- 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)
-
-! 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
- do ispecloc = 1, NSPEC_SURFACE_EXT_MESH
-
- 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)
-
- 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
-!!!! NL NL mute value near source
- if ( (sqrt(((x(i,j) - (X_SOURCE_EXT_MESH))**2 + &
- (y(i,j) - (Y_SOURCE_EXT_MESH))**2 + &
- (z(i,j) - (Z_SOURCE_EXT_MESH))**2)) < RADIUS_TO_MUTE) &
- .and. MUTE_SOURCE) then
-
- display(i,j) = 0.
- else
-
- if(inorm == 1) then
- display(i,j) = vectorx
- else if(inorm == 2) then
- display(i,j) = vectory
- else
- display(i,j) = vectorz
- endif
- endif
- else
- if (USE_EXTERNAL_MESH) then
- display(i,j) = sqrt(vectorz**2+vectory**2+vectorx**2)
- else
- display(i,j) = vectorz
- endif
- 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)
-
- 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
-!!!! NL NL mute value near source
- if ( (sqrt(((dble(xcoord) - (X_SOURCE_EXT_MESH))**2 + &
- (dble(ycoord) - (Y_SOURCE_EXT_MESH))**2 + &
- (dble(zcoord) - (Z_SOURCE_EXT_MESH))**2)) < RADIUS_TO_MUTE) &
- .and. MUTE_SOURCE) then
-
- field_display(ilocnum+ieoff) = 0.
- else
-
-
- 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
- endif
- else
- if (USE_EXTERNAL_MESH) then
- field_display(ilocnum+ieoff) =sqrt(vectorz**2+vectory**2+vectorx**2)
- else
- field_display(ilocnum+ieoff) = dble(vectorz)
- endif
- 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'
- if (USE_EXTERNAL_MESH) then
- call get_global_AVS(nspectot_AVS_max,xp,yp,zp,iglob,loc,ifseg,nglob,npointot, &
- dble(minval(store_val_x(:,0))),dble(maxval(store_val_x(:,0))))
- else
- call get_global_AVS(nspectot_AVS_max,xp,yp,zp,iglob,loc,ifseg,nglob,npointot,UTM_X_MIN,UTM_X_MAX)
- endif
-
-!--- 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(:))
-
-! 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
-
-! compute min and max of data value to normalize
- min_field_current = minval(field_display(:))
- max_field_current = maxval(field_display(:))
-
-! print minimum and maximum amplitude in current snapshot
- print *
- print *,'minimum amplitude in current snapshot after removal = ',min_field_current
- print *,'maximum amplitude in current snapshot after removal = ',max_field_current
- print *
-
- endif
-
-!-----------------------------------------
-
-
-! apply scaling in all cases for movies
- if(.not. plot_shaking_map) then
-
-! make sure range is always symmetric and center is in zero
-! this assumption works only for fields that can be negative
-! would not work for norm of vector for instance
-! (we would lose half of the color palette if no negative values)
- max_absol = max(abs(min_field_current),abs(max_field_current))
- min_field_current = - max_absol
- max_field_current = + max_absol
-
-! normalize field to [0:1]
- field_display(:) = (field_display(:) - min_field_current) / (max_field_current - min_field_current)
-
-! rescale to [-1,1]
- field_display(:) = 2.*field_display(:) - 1.
-
-! apply threshold to normalized field
- if(APPLY_THRESHOLD) &
- where(abs(field_display(:)) <= THRESHOLD) field_display = 0.
-
-! apply non linear scaling to normalized field if needed
- if(NONLINEAR_SCALING) then
- where(field_display(:) >= 0.)
- field_display = field_display ** POWER_SCALING
- elsewhere
- field_display = - abs(field_display) ** POWER_SCALING
- endwhere
- endif
-
-! map back to [0,1]
- field_display(:) = (field_display(:) + 1.) / 2.
-
-! map field to [0:255] for AVS color scale
- field_display(:) = 255. * field_display(:)
-
-
-! apply scaling only if selected for shaking map
-
- else if(NONLINEAR_SCALING .and. iscaling_shake == 1) 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 ******
-
- if(inumber == 1) then
- ivalue = iframe
- else
- ivalue = it
- endif
-
-! create file name and open file
- if(plot_shaking_map) then
-
- if(USE_OPENDX) then
- write(outputname,"('/DX_shaking_map.dx')")
- open(unit=11,file=trim(OUTPUT_FILES)//outputname,status='unknown')
- write(11,*) 'object 1 class array type float rank 1 shape 3 items ',nglob,' data follows'
- else if(USE_AVS) then
- if(UNIQUE_FILE) stop 'cannot use unique file AVS option for shaking map'
- write(outputname,"('/AVS_shaking_map.inp')")
- open(unit=11,file=trim(OUTPUT_FILES)//outputname,status='unknown')
- write(11,*) nglob,' ',nspectot_AVS_max,' 1 0 0'
- else
- stop 'wrong output format selected'
- endif
-
- else
-
- if(USE_OPENDX) then
- write(outputname,"('/DX_movie_',i6.6,'.dx')") ivalue
- open(unit=11,file=trim(OUTPUT_FILES)//outputname,status='unknown')
- write(11,*) 'object 1 class array type float rank 1 shape 3 items ',nglob,' data follows'
- else if(USE_AVS) then
- if(UNIQUE_FILE .and. iframe == 1) then
- open(unit=11,file=trim(OUTPUT_FILES)//'/AVS_movie_all.inp',status='unknown')
- write(11,*) nframes
- write(11,*) 'data'
- write(11,"('step',i1,' image',i1)") 1,1
- write(11,*) nglob,' ',nspectot_AVS_max
- else if(.not. UNIQUE_FILE) then
- write(outputname,"('/AVS_movie_',i6.6,'.inp')") ivalue
- open(unit=11,file=trim(OUTPUT_FILES)//outputname,status='unknown')
- write(11,*) nglob,' ',nspectot_AVS_max,' 1 0 0'
- endif
- else
- stop 'wrong output format selected'
- endif
-
- endif
-
- if(.false.) then
-
- else
-
-! if unique file, output geometry only once
- if(.not. UNIQUE_FILE .or. iframe == 1) then
-
-! output list of points
- mask_point = .false.
- ipoin = 0
- 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
- ipoin = ipoin + 1
- ireorder(ibool_number) = ipoin
- if(USE_OPENDX) then
- write(11,*) sngl(xp_save(ilocnum+ieoff)),sngl(yp_save(ilocnum+ieoff)),sngl(zp_save(ilocnum+ieoff))
- else if(USE_AVS) then
- write(11,*) ireorder(ibool_number),sngl(xp_save(ilocnum+ieoff)), &
- sngl(yp_save(ilocnum+ieoff)),sngl(zp_save(ilocnum+ieoff))
- endif
- endif
- mask_point(ibool_number) = .true.
- enddo
- enddo
-
- if(USE_OPENDX) &
- write(11,*) 'object 2 class array type int rank 1 shape 4 items ',nspectot_AVS_max,' data follows'
-
-! output list of elements
- do ispec=1,nspectot_AVS_max
- ieoff = NGNOD2D_AVS_DX*(ispec-1)
-! four points for each element
- ibool_number1 = iglob(ieoff + 1)
- ibool_number2 = iglob(ieoff + 2)
- ibool_number3 = iglob(ieoff + 3)
- ibool_number4 = iglob(ieoff + 4)
- if(USE_OPENDX) then
-! point order in OpenDX is 1,4,2,3 *not* 1,2,3,4 as in AVS
- write(11,"(i10,1x,i10,1x,i10,1x,i10)") ireorder(ibool_number1)-1, &
- ireorder(ibool_number4)-1,ireorder(ibool_number2)-1,ireorder(ibool_number3)-1
- else
- write(11,"(i10,' 1 quad ',i10,1x,i10,1x,i10,1x,i10)") ispec,ireorder(ibool_number1), &
- ireorder(ibool_number4),ireorder(ibool_number2),ireorder(ibool_number3)
- endif
- enddo
-
- endif
-
- if(USE_OPENDX) then
- write(11,*) 'attribute "element type" string "quads"'
- write(11,*) 'attribute "ref" string "positions"'
- write(11,*) 'object 3 class array type float rank 0 items ',nglob,' data follows'
- else
- if(UNIQUE_FILE) then
-! step number for AVS multistep file
- if(iframe > 1) then
- if(iframe < 10) then
- write(11,"('step',i1,' image',i1)") iframe,iframe
- else if(iframe < 100) then
- write(11,"('step',i2,' image',i2)") iframe,iframe
- else if(iframe < 1000) then
- write(11,"('step',i3,' image',i3)") iframe,iframe
- else
- write(11,"('step',i4,' image',i4)") iframe,iframe
- endif
- endif
- write(11,*) '1 0'
- endif
-! dummy text for labels
- write(11,*) '1 1'
- write(11,*) 'a, b'
- endif
-
-! output data values
- mask_point = .false.
-
-! output point data
- 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
- if(USE_OPENDX) then
- if(plot_shaking_map) then
- write(11,*) sngl(field_display(ilocnum+ieoff))
- else
- write(11,"(f7.2)") field_display(ilocnum+ieoff)
- endif
- else
- if(plot_shaking_map) then
- write(11,*) ireorder(ibool_number),field_display(ilocnum+ieoff)
- else
- write(11,"(i10,1x,f7.2)") ireorder(ibool_number),field_display(ilocnum+ieoff)
- endif
- endif
- endif
- mask_point(ibool_number) = .true.
- enddo
- enddo
-
-! define OpenDX field
- if(USE_OPENDX) then
- write(11,*) 'attribute "dep" string "positions"'
- write(11,*) 'object "irregular positions irregular connections" class field'
- write(11,*) 'component "positions" value 1'
- write(11,*) 'component "connections" value 2'
- write(11,*) 'component "data" value 3'
- write(11,*) 'end'
- endif
-
-! end of test for GMT format
- endif
-
- if(.not. UNIQUE_FILE) close(11)
-
-! end of loop and test on all the time steps for all the movie images
- endif
- enddo
-
- if(UNIQUE_FILE) close(11)
-
- print *
- print *,'done creating movie or shaking map'
- print *
- if(USE_OPENDX) print *,'DX files are stored in ', trim(OUTPUT_FILES), '/DX_*.dx'
- if(USE_AVS) print *,'AVS files are stored in ', trim(OUTPUT_FILES), '/AVS_*.inp'
-
- 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_AVS_DX
-
-!
-!=====================================================================
-!
-
- 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)
- print *, 'UTM_X_MAX', UTM_X_MAX
- print *, 'UTM_X_MIN', UTM_X_MIN
- print *, 'SMALLVALTOL', SMALLVALTOL
-
-! 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
-
Copied: seismo/3D/SPECFEM3D_SESAME/trunk/create_movie_shakemap_AVS_DX_GMT.f90 (from rev 14763, seismo/3D/SPECFEM3D_SESAME/trunk/create_movie_AVS_DX.f90)
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/create_movie_shakemap_AVS_DX_GMT.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/create_movie_shakemap_AVS_DX_GMT.f90 2009-04-29 22:50:18 UTC (rev 14820)
@@ -0,0 +1,995 @@
+!=====================================================================
+!
+! 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.
+!
+!=====================================================================
+
+!
+!--- create a movie of vertical component of surface displacement or velocity
+!--- in AVS or OpenDX format
+!
+
+ program create_movie_AVS_DX
+
+ implicit none
+
+ include "constants.h"
+
+! threshold in percent of the maximum below which we cut the amplitude
+ logical, parameter :: APPLY_THRESHOLD = .true.
+ 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 = .true.
+ real(kind=CUSTOM_REAL), parameter :: POWER_SCALING = 0.13_CUSTOM_REAL
+
+ integer it,it1,it2,ivalue,nspectot_AVS_max,ispec
+ integer iformat,nframes,iframe,inumber,inorm,iscaling_shake
+ integer ibool_number,ibool_number1,ibool_number2,ibool_number3,ibool_number4
+
+ logical USE_OPENDX,USE_AVS,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
+
+! 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
+
+ logical MOVIE_SURFACE,MOVIE_VOLUME,CREATE_SHAKEMAP,SAVE_DISPLACEMENT, &
+ USE_HIGHRES_FOR_MOVIES,SUPPRESS_UTM_PROJECTION,USE_REGULAR_MESH
+ integer NTSTEP_BETWEEN_FRAMES,NTSTEP_BETWEEN_OUTPUT_INFO
+
+ 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,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 ANISOTROPY,SAVE_MESH_FILES,PRINT_SOURCE_TIME_FUNCTION
+
+ character(len=150) OUTPUT_FILES,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
+
+!!!! NL NL for external meshes
+ real(kind=CUSTOM_REAL), parameter :: RADIUS_TO_MUTE = 1000._CUSTOM_REAL
+ logical, parameter :: MUTE_SOURCE = .true.
+ real(kind=CUSTOM_REAL), parameter :: X_SOURCE_EXT_MESH = -9023.021484375
+ real(kind=CUSTOM_REAL), parameter :: Y_SOURCE_EXT_MESH = 6123.611328125
+ real(kind=CUSTOM_REAL), parameter :: Z_SOURCE_EXT_MESH = 17.96331405639648
+ integer, parameter :: NSPEC_SURFACE_EXT_MESH = 15808*4
+
+!!!! NL NL
+
+! ************** PROGRAM STARTS HERE **************
+
+ print *
+ print *,'Recombining all movie frames to create a movie'
+ print *
+
+ print *
+ print *,'reading parameter file'
+ print *
+
+! 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)
+
+! get the base pathname for output files
+ call get_value_string(OUTPUT_FILES, 'OUTPUT_FILES', 'OUTPUT_FILES')
+
+ print *
+ print *,'There are ',NPROC,' slices numbered from 0 to ',NPROC-1
+ 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_EXTERNAL_MESH) then
+ NPROC = 1
+ if(USE_HIGHRES_FOR_MOVIES) then
+ ilocnum = NSPEC_SURFACE_EXT_MESH*NGLLSQUARE
+ else
+ ilocnum = NSPEC_SURFACE_EXT_MESH*NGNOD2D_AVS_DX
+ endif
+ else
+ 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
+ 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))
+
+
+ print *,'1 = create files in OpenDX format'
+ print *,'2 = create files in AVS UCD format'
+ print *,'3 = create files in GMT xyz Ascii long/lat/Uz format'
+ print *,'any other value = exit'
+ print *
+
+ print *,'enter value:'
+ read(5,*) iformat
+
+ if(iformat < 1 .or. iformat > 3) stop 'exiting...'
+
+ if(iformat == 1) then
+ USE_OPENDX = .true.
+ USE_AVS = .false.
+ else if(iformat == 2) then
+ USE_OPENDX = .false.
+ USE_AVS = .true.
+ else
+ USE_OPENDX = .false.
+ USE_AVS = .false.
+ endif
+
+ print *,'movie frames have been saved every ',NTSTEP_BETWEEN_FRAMES,' time steps'
+ print *
+
+ plot_shaking_map = .false.
+ print *,'enter first time step of movie (e.g. 1, enter -1 for shaking map)'
+ read(5,*) it1
+ if(it1 == -1) plot_shaking_map = .true.
+
+ if(.not. plot_shaking_map) then
+
+ print *,'enter last time step of movie (e.g. ',NSTEP,')'
+ read(5,*) it2
+
+ print *
+ print *,'1 = define file names using frame number'
+ print *,'2 = define file names using time step number'
+ print *,'any other value = exit'
+ print *
+ print *,'enter value:'
+ read(5,*) inumber
+ if(inumber<1 .or. inumber>2) stop 'exiting...'
+
+ print *
+ print *,'looping from ',it1,' to ',it2,' every ',NTSTEP_BETWEEN_FRAMES,' time steps'
+
+! count number of movie frames
+ nframes = 0
+ do it = it1,it2
+ if(mod(it,NTSTEP_BETWEEN_FRAMES) == 0) nframes = nframes + 1
+ enddo
+ print *
+ print *,'total number of frames will be ',nframes
+ if(nframes == 0) stop 'null number of frames'
+
+ else
+
+! only one frame if shaking map
+ nframes = 1
+ it1 = 1
+ it2 = 1
+ endif
+
+ iscaling_shake = 0
+ if(plot_shaking_map) then
+ print *
+ print *,'norm to display in shaking map:'
+ print *,'1=displacement 2=velocity 3=acceleration'
+ print *
+ read(5,*) inorm
+ if(inorm < 1 .or. inorm > 3) stop 'incorrect value of inorm'
+
+ print *
+ print *,'apply non-linear scaling to shaking map:'
+ print *,'1=non-linear 2=no scaling'
+ print *
+ read(5,*) iscaling_shake
+ if(iscaling_shake < 1 .or. iscaling_shake > 2) stop 'incorrect value of iscaling_shake'
+ endif
+
+! define the total number of elements at the surface
+ if (USE_EXTERNAL_MESH) then
+ if(USE_HIGHRES_FOR_MOVIES) then
+ nspectot_AVS_max = NSPEC_SURFACE_EXT_MESH * (NGLLX-1) * (NGLLY-1)
+ else
+ nspectot_AVS_max = NSPEC_SURFACE_EXT_MESH
+ endif
+ else
+ 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
+ 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))
+
+ print *
+ if(APPLY_THRESHOLD .and. .not. plot_shaking_map) print *,'Will apply a threshold to amplitude below ',100.*THRESHOLD,' %'
+ if(NONLINEAR_SCALING .and. (.not. plot_shaking_map .or. iscaling_shake == 1)) &
+ print *,'Will apply a non linear scaling with coef ',POWER_SCALING
+
+! --------------------------------------
+
+ 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
+ if(mod(it,NTSTEP_BETWEEN_FRAMES) == 0 .or. plot_shaking_map) then
+
+ iframe = iframe + 1
+
+ print *
+ if(plot_shaking_map) then
+ print *,'reading shaking map snapshot'
+ else
+ print *,'reading snapshot time step ',it,' out of ',NSTEP
+ endif
+ print *
+
+! read all the elements from the same file
+ if(plot_shaking_map) then
+ write(outputname,"('/shakingdata')")
+ else
+ write(outputname,"('/moviedata',i6.6)") it
+ endif
+ open(unit=IOUT,file=trim(OUTPUT_FILES)//outputname,status='old',action='read',form='unformatted')
+ 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)
+
+! 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
+ do ispecloc = 1, NSPEC_SURFACE_EXT_MESH
+
+ 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)
+
+ 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
+!!!! NL NL mute value near source
+ if ( (sqrt(((x(i,j) - (X_SOURCE_EXT_MESH))**2 + &
+ (y(i,j) - (Y_SOURCE_EXT_MESH))**2 + &
+ (z(i,j) - (Z_SOURCE_EXT_MESH))**2)) < RADIUS_TO_MUTE) &
+ .and. MUTE_SOURCE) then
+
+ display(i,j) = 0.
+ else
+
+ if(inorm == 1) then
+ display(i,j) = vectorx
+ else if(inorm == 2) then
+ display(i,j) = vectory
+ else
+ display(i,j) = vectorz
+ endif
+ endif
+ else
+ if (USE_EXTERNAL_MESH) then
+ display(i,j) = sqrt(vectorz**2+vectory**2+vectorx**2)
+ else
+ display(i,j) = vectorz
+ endif
+ 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)
+
+ 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
+!!!! NL NL mute value near source
+ if ( (sqrt(((dble(xcoord) - (X_SOURCE_EXT_MESH))**2 + &
+ (dble(ycoord) - (Y_SOURCE_EXT_MESH))**2 + &
+ (dble(zcoord) - (Z_SOURCE_EXT_MESH))**2)) < RADIUS_TO_MUTE) &
+ .and. MUTE_SOURCE) then
+
+ field_display(ilocnum+ieoff) = 0.
+ else
+
+
+ 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
+ endif
+ else
+ if (USE_EXTERNAL_MESH) then
+ field_display(ilocnum+ieoff) =sqrt(vectorz**2+vectory**2+vectorx**2)
+ else
+ field_display(ilocnum+ieoff) = dble(vectorz)
+ endif
+ 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'
+ if (USE_EXTERNAL_MESH) then
+ call get_global_AVS(nspectot_AVS_max,xp,yp,zp,iglob,loc,ifseg,nglob,npointot, &
+ dble(minval(store_val_x(:,0))),dble(maxval(store_val_x(:,0))))
+ else
+ call get_global_AVS(nspectot_AVS_max,xp,yp,zp,iglob,loc,ifseg,nglob,npointot,UTM_X_MIN,UTM_X_MAX)
+ endif
+
+!--- 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(:))
+
+! 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
+
+! compute min and max of data value to normalize
+ min_field_current = minval(field_display(:))
+ max_field_current = maxval(field_display(:))
+
+! print minimum and maximum amplitude in current snapshot
+ print *
+ print *,'minimum amplitude in current snapshot after removal = ',min_field_current
+ print *,'maximum amplitude in current snapshot after removal = ',max_field_current
+ print *
+
+ endif
+
+!-----------------------------------------
+
+
+! apply scaling in all cases for movies
+ if(.not. plot_shaking_map) then
+
+! make sure range is always symmetric and center is in zero
+! this assumption works only for fields that can be negative
+! would not work for norm of vector for instance
+! (we would lose half of the color palette if no negative values)
+ max_absol = max(abs(min_field_current),abs(max_field_current))
+ min_field_current = - max_absol
+ max_field_current = + max_absol
+
+! normalize field to [0:1]
+ field_display(:) = (field_display(:) - min_field_current) / (max_field_current - min_field_current)
+
+! rescale to [-1,1]
+ field_display(:) = 2.*field_display(:) - 1.
+
+! apply threshold to normalized field
+ if(APPLY_THRESHOLD) &
+ where(abs(field_display(:)) <= THRESHOLD) field_display = 0.
+
+! apply non linear scaling to normalized field if needed
+ if(NONLINEAR_SCALING) then
+ where(field_display(:) >= 0.)
+ field_display = field_display ** POWER_SCALING
+ elsewhere
+ field_display = - abs(field_display) ** POWER_SCALING
+ endwhere
+ endif
+
+! map back to [0,1]
+ field_display(:) = (field_display(:) + 1.) / 2.
+
+! map field to [0:255] for AVS color scale
+ field_display(:) = 255. * field_display(:)
+
+
+! apply scaling only if selected for shaking map
+
+ else if(NONLINEAR_SCALING .and. iscaling_shake == 1) 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 ******
+
+ if(.not. plot_shaking_map) then
+ if(inumber == 1) then
+ ivalue = iframe
+ else
+ ivalue = it
+ endif
+ endif
+
+! create file name and open file
+ if(plot_shaking_map) then
+
+ if(USE_OPENDX) then
+ write(outputname,"('/DX_shaking_map.dx')")
+ open(unit=11,file=trim(OUTPUT_FILES)//outputname,status='unknown')
+ write(11,*) 'object 1 class array type float rank 1 shape 3 items ',nglob,' data follows'
+ else if(USE_AVS) then
+ write(outputname,"('/AVS_shaking_map.inp')")
+ open(unit=11,file=trim(OUTPUT_FILES)//outputname,status='unknown')
+ write(11,*) nglob,' ',nspectot_AVS_max,' 1 0 0'
+ else
+ stop 'wrong output format selected'
+ endif
+
+ else
+
+ if(USE_OPENDX) then
+ write(outputname,"('/DX_movie_',i6.6,'.dx')") ivalue
+ open(unit=11,file=trim(OUTPUT_FILES)//outputname,status='unknown')
+ write(11,*) 'object 1 class array type float rank 1 shape 3 items ',nglob,' data follows'
+ else if(USE_AVS) then
+ write(outputname,"('/AVS_movie_',i6.6,'.inp')") ivalue
+ open(unit=11,file=trim(OUTPUT_FILES)//outputname,status='unknown')
+ write(11,*) nglob,' ',nspectot_AVS_max,' 1 0 0'
+ else
+ stop 'wrong output format selected'
+ endif
+
+ endif
+
+ if(.false.) then
+
+ else
+
+! output list of points
+ mask_point = .false.
+ ipoin = 0
+ 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
+ ipoin = ipoin + 1
+ ireorder(ibool_number) = ipoin
+ if(USE_OPENDX) then
+ write(11,*) xp_save(ilocnum+ieoff),yp_save(ilocnum+ieoff),zp_save(ilocnum+ieoff)
+ else if(USE_AVS) then
+ write(11,*) ireorder(ibool_number),xp_save(ilocnum+ieoff), &
+ yp_save(ilocnum+ieoff),zp_save(ilocnum+ieoff)
+ endif
+ endif
+ mask_point(ibool_number) = .true.
+ enddo
+ enddo
+
+ if(USE_OPENDX) &
+ write(11,*) 'object 2 class array type int rank 1 shape 4 items ',nspectot_AVS_max,' data follows'
+
+! output list of elements
+ do ispec=1,nspectot_AVS_max
+ ieoff = NGNOD2D_AVS_DX*(ispec-1)
+! four points for each element
+ ibool_number1 = iglob(ieoff + 1)
+ ibool_number2 = iglob(ieoff + 2)
+ ibool_number3 = iglob(ieoff + 3)
+ ibool_number4 = iglob(ieoff + 4)
+ if(USE_OPENDX) then
+! point order in OpenDX is 1,4,2,3 *not* 1,2,3,4 as in AVS
+ write(11,"(i10,1x,i10,1x,i10,1x,i10)") ireorder(ibool_number1)-1, &
+ ireorder(ibool_number4)-1,ireorder(ibool_number2)-1,ireorder(ibool_number3)-1
+ else
+ write(11,"(i10,' 1 quad ',i10,1x,i10,1x,i10,1x,i10)") ispec,ireorder(ibool_number1), &
+ ireorder(ibool_number4),ireorder(ibool_number2),ireorder(ibool_number3)
+ endif
+ enddo
+
+ if(USE_OPENDX) then
+ write(11,*) 'attribute "element type" string "quads"'
+ write(11,*) 'attribute "ref" string "positions"'
+ write(11,*) 'object 3 class array type float rank 0 items ',nglob,' data follows'
+ else
+! dummy text for labels
+ write(11,*) '1 1'
+ write(11,*) 'a, b'
+ endif
+
+! output data values
+ mask_point = .false.
+
+! output point data
+ 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
+ if(USE_OPENDX) then
+ if(plot_shaking_map) then
+ write(11,*) sngl(field_display(ilocnum+ieoff))
+ else
+ write(11,"(f7.2)") field_display(ilocnum+ieoff)
+ endif
+ else
+ if(plot_shaking_map) then
+ write(11,*) ireorder(ibool_number),field_display(ilocnum+ieoff)
+ else
+ write(11,"(i10,1x,f7.2)") ireorder(ibool_number),field_display(ilocnum+ieoff)
+ endif
+ endif
+ endif
+ mask_point(ibool_number) = .true.
+ enddo
+ enddo
+
+! define OpenDX field
+ if(USE_OPENDX) then
+ write(11,*) 'attribute "dep" string "positions"'
+ write(11,*) 'object "irregular positions irregular connections" class field'
+ write(11,*) 'component "positions" value 1'
+ write(11,*) 'component "connections" value 2'
+ write(11,*) 'component "data" value 3'
+ write(11,*) 'end'
+ endif
+
+! end of test for GMT format
+ endif
+
+ close(11)
+
+! end of loop and test on all the time steps for all the movie images
+ endif
+ enddo
+
+ print *
+ print *,'done creating movie or shaking map'
+ print *
+ if(USE_OPENDX) print *,'DX files are stored in ', trim(OUTPUT_FILES), '/DX_*.dx'
+ if(USE_AVS) print *,'AVS files are stored in ', trim(OUTPUT_FILES), '/AVS_*.inp'
+
+ 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_AVS_DX
+
+!
+!=====================================================================
+!
+
+ 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)
+ print *, 'UTM_X_MAX', UTM_X_MAX
+ print *, 'UTM_X_MIN', UTM_X_MIN
+ print *, 'SMALLVALTOL', SMALLVALTOL
+
+! 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
+
Property changes on: seismo/3D/SPECFEM3D_SESAME/trunk/create_movie_shakemap_AVS_DX_GMT.f90
___________________________________________________________________
Name: svn:keywords
+ Author Date Id Revision
Name: svn:mergeinfo
+
Name: svn:eol-style
+ native
More information about the CIG-COMMITS
mailing list