[cig-commits] r8468 - in seismo/2D/SPECFEM2D/trunk: . DATA
walter at geodynamics.org
walter at geodynamics.org
Fri Dec 7 15:48:35 PST 2007
Author: walter
Date: 2007-12-07 15:48:34 -0800 (Fri, 07 Dec 2007)
New Revision: 8468
Added:
seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Thomas_Lacq
seismo/2D/SPECFEM2D/trunk/DATA/interfaces_thomas_Lacq.dat
Modified:
seismo/2D/SPECFEM2D/trunk/DATA/Par_file
seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Paul_acoustic
seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Paul_elastic
seismo/2D/SPECFEM2D/trunk/DATA/Par_file_cours_M2_UPPA
seismo/2D/SPECFEM2D/trunk/Makefile
seismo/2D/SPECFEM2D/trunk/checkgrid.f90
seismo/2D/SPECFEM2D/trunk/circ.f90
seismo/2D/SPECFEM2D/trunk/compute_gradient_fluid.f90
seismo/2D/SPECFEM2D/trunk/constants.h
seismo/2D/SPECFEM2D/trunk/createnum_fast.f90
seismo/2D/SPECFEM2D/trunk/createnum_slow.f90
seismo/2D/SPECFEM2D/trunk/defarrays.f90
seismo/2D/SPECFEM2D/trunk/maille_non_struct_2.f90
seismo/2D/SPECFEM2D/trunk/maille_non_struct_3.f90
seismo/2D/SPECFEM2D/trunk/meshfem2D.f90
seismo/2D/SPECFEM2D/trunk/plotpost.f90
seismo/2D/SPECFEM2D/trunk/specfem2D.f90
seismo/2D/SPECFEM2D/trunk/write_seismograms.f90
Log:
added Lacq Par_files to 2D code.
added pressure seismograms to 2D code.
added option to display only lowerleft corner of vector plots.
moved sizemax_arrows from constants.h to Par_file.
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file 2005-02-20 00:09:28 UTC (rev 8467)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file 2007-12-07 23:48:34 UTC (rev 8468)
@@ -5,7 +5,7 @@
#-------------------------------->
# title of job, and file that contains interface data
-title = Test sinusoide elastique
+title = Test sinusoide acoustique
interfacesfile = interface_sinus.dat
# geometry of the model (origin lower-left corner = 0,0) and mesh description
@@ -15,7 +15,7 @@
ngnod = 9 # noeuds de controle pour blocs (4 ou 9)
initialfield = .false. # use a plane wave as source or not
readmodel = .false. # read external earth model or not
-ELASTIC = .true. # elastic or acoustic simulation
+ELASTIC = .false. # elastic or acoustic simulation
TURN_ANISOTROPY_ON = .false. # turn anisotropy on or off
TURN_ATTENUATION_ON = .false. # turn attenuation on or off
@@ -34,44 +34,45 @@
xs = 2000. # source location x in meters
zs = 1200. # source location z in meters
source_type = 1 # force = 1 or explosion = 2
-time_function_type = 1 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4
+time_function_type = 2 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4
f0 = 10.0 # dominant source frequency (Hz) if not Dirac
angleforce = 0. # angle of the source (for a force only)
Mxx = 1. # Mxx component (for a moment tensor source only)
Mzz = 1. # Mzz component (for a moment tensor source only)
Mxz = 0. # Mxz component (for a moment tensor source only)
-factor = 1.d12 # amplification factor
+factor = 1.d4 # amplification factor
# receiver line parameters
-enreg_surf = .true. # enregistrement volume ou surface
-sismostype = 1 # record 1=displ 2=veloc 3=accel
+enreg_surf = .false. # enregistrement volume ou surface
+sismostype = 4 # record 1=displ 2=veloc 3=accel 4=pressure
nrec = 101 # number of receivers
xdeb = 300. # first receiver x in meters
-zdeb = 2200. # first receiver z in meters
+zdeb = 1700. # first receiver z in meters
xfin = 3700. # last receiver x in meters
-zfin = 2200. # last receiver z in meters
+zfin = 1700. # last receiver z in meters
anglerec = 0.d0 # angle to rotate components at receivers
# display parameters
itaff = 100 # display frequency in time steps
output_postscript_snapshot = .true. # output Postscript image of the results
output_PNM_image = .true. # output PNM image of the results
-vecttype = 1 # display 1=displ 2=veloc 3=accel
+vecttype = 2 # display 1=displ 2=veloc 3=accel
cutvect = 1. # amplitude min en % pour vector plots
meshvect = .true. # display mesh on vector plots or not
modelvect = .false. # display velocity model on vector plots
boundvect = .true. # display boundary conditions on plots
interpol = .true. # interpolation of the display or not
-pointsdisp = 6 # points for interpolation of display
+pointsdisp = 6 # points for interpolation of display (set to 1 for lower-left corner only)
subsamp = 1 # subsampling of color snapshots
+sizemax_arrows = 1.d0 # maximum size of arrows on vector plots in cm
gnuplot = .false. # generate a GNUPLOT file for the grid
outputgrid = .false. # save the grid in a text file or not
# velocity and density model (nx,nz)
nbmodels = 2 # nb de modeles differents (0,rho,vp,vs,0,0)
-1 0 2700.d0 3000.d0 1732.051d0 0 0
-2 0 1800.d0 2000.d0 1100.845d0 0 0
-#3 0 2200.d0 2500.d0 1443.375d0 0 0
+1 0 2700.d0 3000.d0 0.d0 0 0
+2 0 1800.d0 2000.d0 0.d0 0 0
+#3 0 2200.d0 2500.d0 0.d0 0 0
nbzone = 2 # nb of zones and model number for each
1 80 1 40 1
1 80 31 40 2
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Paul_acoustic
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Paul_acoustic 2005-02-20 00:09:28 UTC (rev 8467)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Paul_acoustic 2007-12-07 23:48:34 UTC (rev 8468)
@@ -50,7 +50,7 @@
# receiver line parameters
#
enreg_surf = .false. # enregistrement volume ou surface
-sismostype = 2 # record 1=displ 2=veloc 3=accel
+sismostype = 2 # record 1=displ 2=veloc 3=accel 4=pressure
nrec = 101 # number of receivers
xdeb = 300. # first receiver x in meters
zdeb = 1700. # first receiver z in meters
@@ -69,8 +69,9 @@
modelvect = .false. # display velocity model on vector plots
boundvect = .true. # display boundary conditions on plots
interpol = .true. # interpolation of the display or not
-pointsdisp = 6 # points for interpolation of display
+pointsdisp = 6 # points for interpolation of display (set to 1 for lower-left corner only)
subsamp = 1 # subsampling of color snapshots
+sizemax_arrows = 1.d0 # maximum size of arrows on vector plots in cm
gnuplot = .false. # generate a GNUPLOT file for the grid
outputgrid = .false. # save the grid in a text file or not
#
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Paul_elastic
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Paul_elastic 2005-02-20 00:09:28 UTC (rev 8467)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Paul_elastic 2007-12-07 23:48:34 UTC (rev 8468)
@@ -44,7 +44,7 @@
# receiver line parameters
enreg_surf = .true. # enregistrement volume ou surface
-sismostype = 1 # record 1=displ 2=veloc 3=accel
+sismostype = 1 # record 1=displ 2=veloc 3=accel 4=pressure
nrec = 101 # number of receivers
xdeb = 300. # first receiver x in meters
zdeb = 2200. # first receiver z in meters
@@ -62,8 +62,9 @@
modelvect = .false. # display velocity model on vector plots
boundvect = .true. # display boundary conditions on plots
interpol = .true. # interpolation of the display or not
-pointsdisp = 6 # points for interpolation of display
+pointsdisp = 6 # points for interpolation of display (set to 1 for lower-left corner only)
subsamp = 1 # subsampling of color snapshots
+sizemax_arrows = 1.d0 # maximum size of arrows on vector plots in cm
gnuplot = .false. # generate a GNUPLOT file for the grid
outputgrid = .false. # save the grid in a text file or not
Added: seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Thomas_Lacq
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Thomas_Lacq 2005-02-20 00:09:28 UTC (rev 8467)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Thomas_Lacq 2007-12-07 23:48:34 UTC (rev 8468)
@@ -0,0 +1,78 @@
+#
+# This is the parameter file
+# Put variable names first and actual value after 34th column
+#
+#-------------------------------->
+
+# title of job, and file that contains interface data
+title = Lacq Thomas
+interfacesfile = interfaces_thomas_Lacq.dat
+
+# geometry of the model (origin lower-left corner = 0,0) and mesh description
+xmin = 0.d0 # abscissa of left side of the model
+xmax = 25000.d0 # abscissa of right side of the model
+nx = 250 # number of elements along X
+ngnod = 9 # noeuds de controle pour blocs (4 ou 9)
+initialfield = .false. # use a plane wave as source or not
+read_external_model = .false. # read external earth model or not
+ELASTIC = .true. # elastic or acoustic simulation
+TURN_ANISOTROPY_ON = .false. # turn anisotropy on or off
+TURN_ATTENUATION_ON = .false. # turn attenuation on or off
+
+# absorbing boundaries parameters
+absorbhaut = .false. # absorbing boundary active or not
+absorbbas = .true.
+absorbgauche = .true.
+absorbdroite = .true.
+
+# time step parameters
+nt = 10000 # nb total de pas de temps
+dt = 1.d-3 # valeur du pas de temps
+
+# source parameters
+source_surf = .false. # source dans le volume ou a la surface
+xs = 12500. # source location x in meters
+zs = 10000. # source location z in meters
+source_type = 1 # force = 1 or explosion = 2
+time_function_type = 3 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4
+f0 = 10.0 # dominant source frequency (Hz) if not Dirac
+angleforce = 0. # angle of the source (for a force only)
+Mxx = 1. # Mxx component (for a moment tensor source only)
+Mzz = 1. # Mzz component (for a moment tensor source only)
+Mxz = 0. # Mxz component (for a moment tensor source only)
+factor = 0.1 # amplification factor
+
+# receiver line parameters
+enreg_surf = .false. # enregistrement volume ou surface
+sismostype = 2 # record 1=displ 2=veloc 3=accel 4=pressure
+nrec = 4 # number of receivers
+xdeb = 5000. # first receiver x in meters
+zdeb = 15000. # first receiver z in meters
+xfin = 20000. # last receiver x in meters
+zfin = 15000. # last receiver z in meters
+anglerec = 0.d0 # angle to rotate components at receivers
+
+# display parameters
+itaff = 250 # display frequency in time steps
+output_postscript_snapshot = .true. # output Postscript image of the results
+output_PNM_image = .false. # output PNM image of the results
+vecttype = 1 # display 1=displ 2=veloc 3=accel
+cutvect = 1.d0 # amplitude min en % pour vector plots
+meshvect = .false. # display mesh on vector plots or not
+modelvect = .true. # display velocity model on vector plots
+boundvect = .true. # display boundary conditions on plots
+interpol = .true. # interpolation of the display or not
+pointsdisp = 1 # points for interpolation of display (set to 1 for lower-left corner only)
+subsamp = 2 # subsampling of color snapshots
+sizemax_arrows = 0.4d0 # maximum size of arrows on vector plots in cm
+gnuplot = .false. # generate a GNUPLOT file for the grid
+outputgrid = .false. # save the grid in a text file or not
+
+# velocity and density model (nx,nz)
+nbmodels = 1 # nb de modeles differents (0,rho,vp,vs,0,0)
+1 0 2700.d0 3000.d0 1732.051d0 0 0
+#2 0 1800.d0 2000.d0 1100.845d0 0 0
+#3 0 2200.d0 2500.d0 1443.375d0 0 0
+nbzone = 1 # nb of zones and model number for each
+1 250 1 150 1
+#1 80 41 60 3
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file_cours_M2_UPPA
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file_cours_M2_UPPA 2005-02-20 00:09:28 UTC (rev 8467)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file_cours_M2_UPPA 2007-12-07 23:48:34 UTC (rev 8468)
@@ -50,7 +50,7 @@
# receiver line parameters
#
enreg_surf = .true. # enregistrement volume ou surface
-sismostype = 1 # record 1=displ 2=veloc 3=accel
+sismostype = 1 # record 1=displ 2=veloc 3=accel 4=pressure
nrec = 11 # number of receivers
xdeb = 300. # first receiver x in meters
zdeb = 2200. # first receiver z in meters
@@ -69,8 +69,9 @@
modelvect = .false. # display velocity model on vector plots
boundvect = .true. # display boundary conditions on plots
interpol = .true. # interpolation of the display or not
-pointsdisp = 6 # points for interpolation of display
+pointsdisp = 6 # points for interpolation of display (set to 1 for lower-left corner only)
subsamp = 1 # subsampling of color snapshots
+sizemax_arrows = 1.d0 # maximum size of arrows on vector plots in cm
gnuplot = .false. # generate a GNUPLOT file for the grid
outputgrid = .false. # save the grid in a text file or not
#
Added: seismo/2D/SPECFEM2D/trunk/DATA/interfaces_thomas_Lacq.dat
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/interfaces_thomas_Lacq.dat 2005-02-20 00:09:28 UTC (rev 8467)
+++ seismo/2D/SPECFEM2D/trunk/DATA/interfaces_thomas_Lacq.dat 2007-12-07 23:48:34 UTC (rev 8468)
@@ -0,0 +1,28 @@
+#
+# number of interfaces
+#
+2
+#
+# for each interface below, we give the number of points and then x,y for each point
+#
+#
+# interface number 1 (bottom of the mesh)
+#
+ 2
+ 0 0
+ 25000 0
+#
+# interface number 4 (topography, top of the mesh)
+#
+ 2
+ 0 15000
+ 25000 15000
+#
+# for each layer, we give the number of spectral elements in the vertical direction
+#
+#
+# layer number 1 (bottom layer)
+#
+150
+#
+#
Modified: seismo/2D/SPECFEM2D/trunk/Makefile
===================================================================
--- seismo/2D/SPECFEM2D/trunk/Makefile 2005-02-20 00:09:28 UTC (rev 8467)
+++ seismo/2D/SPECFEM2D/trunk/Makefile 2007-12-07 23:48:34 UTC (rev 8468)
@@ -14,8 +14,8 @@
# Intel Linux
F90 = ifort
-FLAGS_CHECK=-O0 -implicitnone -warn stderrors -warn truncated_source -warn argument_checking -warn unused -warn declarations -std95 -check bounds
-FLAGS_NOCHECK=-O3 -implicitnone -warn stderrors -warn truncated_source -warn argument_checking -warn unused -warn declarations -std95 -check nobounds
+FLAGS_CHECK=-O0 -implicitnone -warn stderrors -warn truncated_source -warn argument_checking -warn unused -warn declarations -std95 -assume byterecl -check bounds
+FLAGS_NOCHECK=-O3 -implicitnone -warn stderrors -warn truncated_source -warn argument_checking -warn unused -warn declarations -std95 -assume byterecl -check nobounds
#FLAGS_NOCHECK = $(FLAGS_CHECK)
#
Modified: seismo/2D/SPECFEM2D/trunk/checkgrid.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/checkgrid.f90 2005-02-20 00:09:28 UTC (rev 8467)
+++ seismo/2D/SPECFEM2D/trunk/checkgrid.f90 2007-12-07 23:48:34 UTC (rev 8468)
@@ -12,7 +12,7 @@
!========================================================================
subroutine checkgrid(deltat,f0,t0,initialfield,rsizemin,rsizemax, &
- cpoverdxmin,cpoverdxmax,rlamdaSmin,rlamdaSmax,rlamdaPmin,rlamdaPmax)
+ cpoverdxmax,rlamdaSmin,rlamdaSmax,rlamdaPmin,rlamdaPmax)
!
!---- verification taille des mailles, stabilite et nb de points par lambda
@@ -23,7 +23,7 @@
include "constants.h"
double precision f0,t0
- double precision deltat,rsizemin,rsizemax,cpoverdxmin,cpoverdxmax, &
+ double precision deltat,rsizemin,rsizemax,cpoverdxmax, &
rlamdaSmin,rlamdaSmax,rlamdaPmin,rlamdaPmax
logical initialfield
@@ -32,40 +32,39 @@
!---- verification taille de grille min et max
!
- print *
- print *,'******************************************'
- print *,'*** Verification parametres simulation ***'
- print *,'******************************************'
- print *
- print *,'*** Taille max grille = ',rsizemax
- print *,'*** Taille min grille = ',rsizemin
- print *,'*** Rapport max/min = ',rsizemax/rsizemin
- print *
- print *,'*** Stabilite max vitesse P = ',cpoverdxmax*deltat
- print *,'*** Stabilite min vitesse P = ',cpoverdxmin*deltat
- print *
+ write(IOUT,*)
+ write(IOUT,*) '******************************************'
+ write(IOUT,*) '*** Verification parametres simulation ***'
+ write(IOUT,*) '******************************************'
+ write(IOUT,*)
+ write(IOUT,*) '*** Taille max grille = ',rsizemax
+ write(IOUT,*) '*** Taille min grille = ',rsizemin
+ write(IOUT,*) '*** Rapport max/min = ',rsizemax/rsizemin
+ write(IOUT,*)
+ write(IOUT,*) '*** Stabilite max vitesse P = ',cpoverdxmax*deltat
+ write(IOUT,*)
if(.not. initialfield) then
- print *,' Onset time = ',t0
- print *,' Fundamental period = ',1.d0/f0
- print *,' Fundamental frequency = ',f0
+ write(IOUT,*) ' Onset time = ',t0
+ write(IOUT,*) ' Fundamental period = ',1.d0/f0
+ write(IOUT,*) ' Fundamental frequency = ',f0
if(t0 <= 1.d0/f0) then
stop 'Onset time too small'
else
- print *,' --> onset time ok'
+ write(IOUT,*) ' --> onset time ok'
endif
- print *,'----'
- print *,' Nb pts / lambda P max f0 = ',NGLLX*rlamdaPmax/f0
- print *,' Nb pts / lambda P min f0 = ',NGLLX*rlamdaPmin/f0
- print *,' Nb pts / lambda P max fmax = ',NGLLX*rlamdaPmax/(2.5d0*f0)
- print *,' Nb pts / lambda P min fmax = ',NGLLX*rlamdaPmin/(2.5d0*f0)
- print *,'----'
- print *,' Nb pts / lambda S max f0 = ',NGLLX*rlamdaSmax/f0
- print *,' Nb pts / lambda S min f0 = ',NGLLX*rlamdaSmin/f0
- print *,' Nb pts / lambda S max fmax = ',NGLLX*rlamdaSmax/(2.5d0*f0)
- print *,' Nb pts / lambda S min fmax = ',NGLLX*rlamdaSmin/(2.5d0*f0)
- print *,'----'
+ write(IOUT,*) '----'
+ write(IOUT,*) ' Nb pts / lambda P max f0 = ',NGLLX*rlamdaPmax/f0
+ write(IOUT,*) ' Nb pts / lambda P min f0 = ',NGLLX*rlamdaPmin/f0
+ write(IOUT,*) ' Nb pts / lambda P max fmax = ',NGLLX*rlamdaPmax/(2.5d0*f0)
+ write(IOUT,*) ' Nb pts / lambda P min fmax = ',NGLLX*rlamdaPmin/(2.5d0*f0)
+ write(IOUT,*) '----'
+ write(IOUT,*) ' Nb pts / lambda S max f0 = ',NGLLX*rlamdaSmax/f0
+ write(IOUT,*) ' Nb pts / lambda S min f0 = ',NGLLX*rlamdaSmin/f0
+ write(IOUT,*) ' Nb pts / lambda S max fmax = ',NGLLX*rlamdaSmax/(2.5d0*f0)
+ write(IOUT,*) ' Nb pts / lambda S min fmax = ',NGLLX*rlamdaSmin/(2.5d0*f0)
+ write(IOUT,*) '----'
endif
Modified: seismo/2D/SPECFEM2D/trunk/circ.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/circ.f90 2005-02-20 00:09:28 UTC (rev 8467)
+++ seismo/2D/SPECFEM2D/trunk/circ.f90 2007-12-07 23:48:34 UTC (rev 8468)
@@ -36,7 +36,7 @@
integer ndofn,ndime,ngnod,nnode,nbcnd,n1ana
integer nofst,npgeo,nspel,nbmodeles,nbsources,nrec,lquad,isamp,nrec1,nrec2
integer irec,imatnum,netyp,nxgll,nelemperio,nelemabs,nx,nz,i,j
- integer irepr,nrecsur3,nt,niter,itaff,itfirstaff,numerocourant,iptsdisp,isubsamp
+ integer irepr,nrecsur3,nt,niter,itaff,itfirstaff,numerocourant,pointsdisp,isubsamp
double precision R,theta_i,theta_init,delta_theta,eta_j,valseuil,freqmaxrep
double precision f0,t0,xs,zs,angle,factor,dist,xoffs,zoffs
@@ -706,7 +706,7 @@
scalex = 1.
scalez = 1.
sizemax = 1.
- iptsdisp = 7
+ pointsdisp = 7
isubsamp = 2
orig_x = 2.3
orig_z = 3.4
@@ -812,7 +812,7 @@
nelemabs = 0
write(15,*) 'params spectraux'
- write(15,*) netyp,nbmodeles,ngnod,nxgll,nxgll,nspel,iptsdisp, &
+ write(15,*) netyp,nbmodeles,ngnod,nxgll,nxgll,nspel,pointsdisp, &
nelemabs,nelemperio
write(15,*) 'Material sets (num 0 rho vp vs 0 0)'
Modified: seismo/2D/SPECFEM2D/trunk/compute_gradient_fluid.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_gradient_fluid.f90 2005-02-20 00:09:28 UTC (rev 8467)
+++ seismo/2D/SPECFEM2D/trunk/compute_gradient_fluid.f90 2007-12-07 23:48:34 UTC (rev 8468)
@@ -26,6 +26,7 @@
double precision, dimension(NGLLX,NGLLZ,NSPEC) :: xix,xiz,gammax,gammaz
+! for compatibility with elastic arrays, potential is declared as a vector but correctly used below as a scalar
double precision, dimension(NDIM,npoin) :: potential,veloc_field_postscript
! array with derivatives of Lagrange polynomials
Modified: seismo/2D/SPECFEM2D/trunk/constants.h
===================================================================
--- seismo/2D/SPECFEM2D/trunk/constants.h 2005-02-20 00:09:28 UTC (rev 8467)
+++ seismo/2D/SPECFEM2D/trunk/constants.h 2007-12-07 23:48:34 UTC (rev 8468)
@@ -62,9 +62,6 @@
double precision, parameter :: SCALEX = 1.d0
double precision, parameter :: SCALEZ = 1.d0
-! taille de la plus grande fleche en centimetres
- double precision, parameter :: SIZEMAX = 1.d0
-
! US letter paper or European A4
logical, parameter :: US_LETTER = .false.
Modified: seismo/2D/SPECFEM2D/trunk/createnum_fast.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/createnum_fast.f90 2005-02-20 00:09:28 UTC (rev 8467)
+++ seismo/2D/SPECFEM2D/trunk/createnum_fast.f90 2007-12-07 23:48:34 UTC (rev 8468)
@@ -38,10 +38,10 @@
double precision xcor,ycor
!---- create global mesh numbering
- print *
- print *
- print *,'Generating global mesh numbering (fast version)...'
- print *
+ write(IOUT,*)
+ write(IOUT,*)
+ write(IOUT,*) 'Generating global mesh numbering (fast version)...'
+ write(IOUT,*)
nxyz = NGLLX*NGLLZ
ntot = nxyz*nspec
@@ -197,9 +197,9 @@
! verification de la coherence de la numerotation generee
if(minval(ibool) /= 1 .or. maxval(ibool) /= npoin) stop 'Error while generating global numbering'
- print *
- print *,'Total number of points of the global mesh: ',npoin
- print *
+ write(IOUT,*)
+ write(IOUT,*) 'Total number of points of the global mesh: ',npoin
+ write(IOUT,*)
end subroutine createnum_fast
Modified: seismo/2D/SPECFEM2D/trunk/createnum_slow.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/createnum_slow.f90 2005-02-20 00:09:28 UTC (rev 8467)
+++ seismo/2D/SPECFEM2D/trunk/createnum_slow.f90 2007-12-07 23:48:34 UTC (rev 8468)
@@ -31,9 +31,9 @@
integer ngnoddeb(4),ngnodfin(4)
!---- create global mesh numbering
- print *
- print *,'Generating global mesh numbering (slow version)...'
- print *
+ write(IOUT,*)
+ write(IOUT,*) 'Generating global mesh numbering (slow version)...'
+ write(IOUT,*)
npoin = 0
npedge = 0
@@ -277,13 +277,13 @@
! verification de la coherence de la numerotation generee
if(minval(ibool) /= 1 .or. maxval(ibool) /= npoin) stop 'Error while generating global numbering'
- print *,'Total number of points of the global mesh: ',npoin
- print *,'distributed as follows:'
- print *
- print *,'Number of interior points: ',npoin-npedge-npcorn
- print *,'Number of edge points (without corners): ',npedge
- print *,'Number of corner points: ',npcorn
- print *
+ write(IOUT,*) 'Total number of points of the global mesh: ',npoin
+ write(IOUT,*) 'distributed as follows:'
+ write(IOUT,*)
+ write(IOUT,*) 'Number of interior points: ',npoin-npedge-npcorn
+ write(IOUT,*) 'Number of edge points (without corners): ',npedge
+ write(IOUT,*) 'Number of corner points: ',npcorn
+ write(IOUT,*)
end subroutine createnum_slow
Modified: seismo/2D/SPECFEM2D/trunk/defarrays.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/defarrays.f90 2005-02-20 00:09:28 UTC (rev 8467)
+++ seismo/2D/SPECFEM2D/trunk/defarrays.f90 2007-12-07 23:48:34 UTC (rev 8468)
@@ -13,7 +13,7 @@
subroutine defarrays(vpext,vsext,rhoext,density,elastcoef, &
ibool,kmato,coord,npoin,rsizemin,rsizemax, &
- cpoverdxmin,cpoverdxmax,lambdaSmin,lambdaSmax,lambdaPmin,lambdaPmax, &
+ cpoverdxmax,lambdaSmin,lambdaSmax,lambdaPmin,lambdaPmax, &
vpmin,vpmax,read_external_model,nspec,numat)
! define all the arrays for the variational formulation
@@ -40,7 +40,7 @@
double precision x1,z1,x2,z2,rdist1,rdist2,rapportmin,rapportmax
double precision lambdamin,lambdamax
- double precision rsizemin,rsizemax,cpoverdxmin,cpoverdxmax, &
+ double precision rsizemin,rsizemax,cpoverdxmax, &
lambdaSmin,lambdaSmax,lambdaPmin,lambdaPmax,vpmin,vpmax
logical read_external_model
@@ -61,7 +61,6 @@
rsizemin = HUGEVAL
rsizemax = -HUGEVAL
- cpoverdxmin = HUGEVAL
cpoverdxmax = -HUGEVAL
lambdaPmin = HUGEVAL
@@ -92,9 +91,9 @@
cploc = vpext(ipointnum)
csloc = vsext(ipointnum)
denst = rhoext(ipointnum)
- mu = denst*csloc*csloc
- lambda = denst*cploc*cploc - 2.d0*mu
- lambdaplus2mu = lambda + 2.d0*mu
+ mu = denst*csloc*csloc
+ lambda = denst*cploc*cploc - 2.d0*mu
+ lambdaplus2mu = lambda + 2.d0*mu
endif
!--- calculer min et max du modele de vitesse et densite
@@ -127,7 +126,6 @@
rapportmin = cploc / max(rdist1,rdist2)
rapportmax = cploc / min(rdist1,rdist2)
- cpoverdxmin = min(cpoverdxmin,rapportmin)
cpoverdxmax = max(cpoverdxmax,rapportmax)
x0 = coord(1,ibool(1,1,ispec))
@@ -156,13 +154,13 @@
enddo
enddo
- print *
- print *,'********'
- print *,'Modele : vitesse P min,max = ',vpmin,vpmax
- print *,'Modele : vitesse S min,max = ',vsmin,vsmax
- print *,'Modele : densite min,max = ',densmin,densmax
- print *,'********'
- print *
+ write(IOUT,*)
+ write(IOUT,*) '********'
+ write(IOUT,*) 'Modele : vitesse P min,max = ',vpmin,vpmax
+ write(IOUT,*) 'Modele : vitesse S min,max = ',vsmin,vsmax
+ write(IOUT,*) 'Modele : densite min,max = ',densmin,densmax
+ write(IOUT,*) '********'
+ write(IOUT,*)
end subroutine defarrays
Modified: seismo/2D/SPECFEM2D/trunk/maille_non_struct_2.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/maille_non_struct_2.f90 2005-02-20 00:09:28 UTC (rev 8467)
+++ seismo/2D/SPECFEM2D/trunk/maille_non_struct_2.f90 2007-12-07 23:48:34 UTC (rev 8468)
@@ -52,7 +52,7 @@
integer nelemabs,npgeo,nspec,ntopo,nspecvolume,nspecWz
integer k,ix,iz,irec,i,j,iadd
integer imodele,nbmodeles,iaffinfo
- integer itaff,itfirstaff,iptsdisp,isubsamp,nrec,n1ana,n2ana
+ integer itaff,itfirstaff,pointsdisp,isubsamp,nrec,n1ana,n2ana
integer irepr,nrec1,nrec2,isamp,nbsources,isismostype,ivecttype
integer ngnod,nt,niter,idegpoly,nx,nz
integer icodematread
@@ -329,7 +329,7 @@
read(10,4)junk,imodelvect
read(10,4)junk,iboundvect
read(10,4)junk,interpol
- read(10,2)junk,iptsdisp
+ read(10,2)junk,pointsdisp
read(10,2)junk,isubsamp
read(10,1)junk,scalex
read(10,1)junk,scalez
@@ -645,8 +645,8 @@
netyp = 2
nxgll = idegpoly + 1
- write(15,*) 'netyp numat ngnod nxgll nygll nspec iptsdisp ielemabs ielemperio'
- write(15,*) netyp,nbmodeles,ngnod,nxgll,nxgll,nspec,iptsdisp, &
+ write(15,*) 'netyp numat ngnod nxgll nygll nspec pointsdisp ielemabs ielemperio'
+ write(15,*) netyp,nbmodeles,ngnod,nxgll,nxgll,nspec,pointsdisp, &
nelemabs,nelemperio
write(15,*) 'Material sets (num 0 rho vp vs 0 0)'
Modified: seismo/2D/SPECFEM2D/trunk/maille_non_struct_3.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/maille_non_struct_3.f90 2005-02-20 00:09:28 UTC (rev 8467)
+++ seismo/2D/SPECFEM2D/trunk/maille_non_struct_3.f90 2007-12-07 23:48:34 UTC (rev 8468)
@@ -52,7 +52,7 @@
integer nelemabs,npgeo,nspec,ntopo,nspecvolume,nspecWz
integer k,ix,iz,irec,i,j,iadd
integer imodele,nbmodeles,iaffinfo
- integer itaff,itfirstaff,iptsdisp,isubsamp,nrec,n1ana,n2ana
+ integer itaff,itfirstaff,pointsdisp,isubsamp,nrec,n1ana,n2ana
integer irepr,nrec1,nrec2,isamp,nbsources,isismostype,ivecttype
integer ngnod,nt,niter,idegpoly,nx,nz
integer icodematread
@@ -329,7 +329,7 @@
read(10,4)junk,imodelvect
read(10,4)junk,iboundvect
read(10,4)junk,interpol
- read(10,2)junk,iptsdisp
+ read(10,2)junk,pointsdisp
read(10,2)junk,isubsamp
read(10,1)junk,scalex
read(10,1)junk,scalez
@@ -645,8 +645,8 @@
netyp = 2
nxgll = idegpoly + 1
- write(15,*) 'netyp numat ngnod nxgll nygll nspec iptsdisp ielemabs ielemperio'
- write(15,*) netyp,nbmodeles,ngnod,nxgll,nxgll,nspec,iptsdisp, &
+ write(15,*) 'netyp numat ngnod nxgll nygll nspec pointsdisp ielemabs ielemperio'
+ write(15,*) netyp,nbmodeles,ngnod,nxgll,nxgll,nspec,pointsdisp, &
nelemabs,nelemperio
write(15,*) 'Material sets (num 0 rho vp vs 0 0)'
Modified: seismo/2D/SPECFEM2D/trunk/meshfem2D.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/meshfem2D.f90 2005-02-20 00:09:28 UTC (rev 8467)
+++ seismo/2D/SPECFEM2D/trunk/meshfem2D.f90 2007-12-07 23:48:34 UTC (rev 8468)
@@ -63,10 +63,10 @@
integer ngnod,nt,nx,nz,nxread,nzread
integer icodematread
- logical codehaut,codebas,codegauche,codedroite,output_postscript_snapshot,output_PNM_image
+ logical codehaut,codebas,codegauche,codedroite,output_postscript_snapshot,output_PNM_image,plot_lowerleft_corner_only
double precision tang1,tangN,vpzone,vszone,poisson_ratio
- double precision cutvect,xspacerec,zspacerec
+ double precision cutvect,xspacerec,zspacerec,sizemax_arrows
double precision anglerec,xfin,zfin,xdeb,zdeb,xmin,xmax,dt
double precision rhoread,cpread,csread,aniso3read,aniso4read
@@ -252,9 +252,18 @@
call read_value_logical(IIN,IGNORE_JUNK,interpol)
call read_value_integer(IIN,IGNORE_JUNK,pointsdisp)
call read_value_integer(IIN,IGNORE_JUNK,subsamp)
+ call read_value_double_precision(IIN,IGNORE_JUNK,sizemax_arrows)
call read_value_logical(IIN,IGNORE_JUNK,gnuplot)
call read_value_logical(IIN,IGNORE_JUNK,outputgrid)
+! can use only one point to display lower-left corner only for interpolated snapshot
+ if(pointsdisp < 3) then
+ pointsdisp = 3
+ plot_lowerleft_corner_only = .true.
+ else
+ plot_lowerleft_corner_only = .false.
+ endif
+
! lecture des differents modeles de materiaux
call read_value_integer(IIN,IGNORE_JUNK,nbmodeles)
if(nbmodeles <= 0) stop 'Negative number of models not allowed !!'
@@ -611,8 +620,8 @@
write(15,*) 'itaff output_postscript_snapshot output_PNM_image colors numbers'
write(15,*) itaff,output_postscript_snapshot,output_PNM_image,' 1 0'
- write(15,*) 'meshvect modelvect boundvect cutvect subsamp nx_sem_PNM'
- write(15,*) meshvect,modelvect,boundvect,cutvect,subsamp,nxread
+ write(15,*) 'meshvect modelvect boundvect cutvect subsamp sizemax_arrows nx_sem_PNM'
+ write(15,*) meshvect,modelvect,boundvect,cutvect,subsamp,sizemax_arrows,nxread
write(15,*) 'anglerec'
write(15,*) anglerec
@@ -665,8 +674,8 @@
! on a deux fois trop d'elements si elements 9 noeuds
if(ngnod == 9) nelemsurface = nelemsurface / 2
- write(15,*) 'numat ngnod nspec pointsdisp nelemabs nelemsurface'
- write(15,*) nbmodeles,ngnod,nspec,pointsdisp,nelemabs,nelemsurface
+ write(15,*) 'numat ngnod nspec pointsdisp plot_lowerleft_corner_only nelemabs nelemsurface'
+ write(15,*) nbmodeles,ngnod,nspec,pointsdisp,plot_lowerleft_corner_only,nelemabs,nelemsurface
write(15,*) 'Material sets (num 0 rho vp vs 0 0) or (num 1 rho c11 c13 c33 c44)'
do i=1,nbmodeles
Modified: seismo/2D/SPECFEM2D/trunk/plotpost.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/plotpost.f90 2005-02-20 00:09:28 UTC (rev 8467)
+++ seismo/2D/SPECFEM2D/trunk/plotpost.f90 2007-12-07 23:48:34 UTC (rev 8468)
@@ -15,7 +15,8 @@
xinterp,zinterp,shapeint,Uxinterp,Uzinterp,flagrange,density,elastcoef,knods,kmato,ibool, &
numabs,codeabs,anyabs,stitle,npoin,npgeo,vpmin,vpmax,nrec, &
colors,numbers,subsamp,vecttype,interpol,meshvect,modelvect, &
- boundvect,read_external_model,cutvect,nelemabs,numat,iptsdisp,nspec,ngnod,ELASTIC)
+ boundvect,read_external_model,cutvect,sizemax_arrows,nelemabs,numat,pointsdisp,nspec,ngnod,ELASTIC, &
+ plot_lowerleft_corner_only)
!
! routine affichage postscript
@@ -29,17 +30,17 @@
integer, parameter :: MAXCOLORS = 100
double precision, dimension(MAXCOLORS) :: red,green,blue
- integer it,nrec,nelemabs,numat,iptsdisp,nspec
+ integer it,nrec,nelemabs,numat,pointsdisp,pointsdisp_loop,nspec
integer i,npoin,npgeo,ngnod
integer kmato(nspec),knods(ngnod,nspec)
integer ibool(NGLLX,NGLLZ,nspec)
- double precision xinterp(iptsdisp,iptsdisp),zinterp(iptsdisp,iptsdisp)
- double precision shapeint(ngnod,iptsdisp,iptsdisp)
- double precision Uxinterp(iptsdisp,iptsdisp)
- double precision Uzinterp(iptsdisp,iptsdisp)
- double precision flagrange(NGLLX,iptsdisp)
+ double precision xinterp(pointsdisp,pointsdisp),zinterp(pointsdisp,pointsdisp)
+ double precision shapeint(ngnod,pointsdisp,pointsdisp)
+ double precision Uxinterp(pointsdisp,pointsdisp)
+ double precision Uzinterp(pointsdisp,pointsdisp)
+ double precision flagrange(NGLLX,pointsdisp)
double precision density(numat),elastcoef(4,numat)
double precision dt,timeval,x_source,z_source
@@ -50,7 +51,7 @@
double precision, dimension(nrec) :: st_xval,st_zval
integer numabs(nelemabs),codeabs(4,nelemabs)
- logical anyabs,ELASTIC
+ logical anyabs,ELASTIC,plot_lowerleft_corner_only
double precision xmax,zmax,height,xw,zw,usoffset,sizex,sizez,vpmin,vpmax
@@ -67,7 +68,7 @@
integer colors,numbers,subsamp,vecttype
logical interpol,meshvect,modelvect,boundvect,read_external_model
- double precision cutvect
+ double precision cutvect,sizemax_arrows
double precision rapp_page,dispmax,xmin,zmin
@@ -171,8 +172,8 @@
! recherche des positions maximales des points de la grille
xmax=maxval(coord(1,:))
zmax=maxval(coord(2,:))
- write(*,*) 'Max X = ',xmax
- write(*,*) 'Max Z = ',zmax
+ write(IOUT,*) 'Max X = ',xmax
+ write(IOUT,*) 'Max Z = ',zmax
! limite du repere physique
xmin=0.d0
@@ -183,7 +184,7 @@
! recherche de la valeur maximum de la norme du deplacement
dispmax = maxval(sqrt(displ(1,:)**2 + displ(2,:)**2))
- write(*,*) 'Max norme = ',dispmax
+ write(IOUT,*) 'Max norme = ',dispmax
! hauteur des numeros de domaine en CM
height = 0.25d0
@@ -355,7 +356,7 @@
!---- print the spectral elements mesh in PostScript
!
- print *,'Shape functions based on ',ngnod,' control nodes'
+ write(IOUT,*) 'Shape functions based on ',ngnod,' control nodes'
convert = pi/180.d0
@@ -444,8 +445,8 @@
write(24,*) '% elem ',ispec
- do i=1,iptsdisp
- do j=1,iptsdisp
+ do i=1,pointsdisp
+ do j=1,pointsdisp
xinterp(i,j) = 0.d0
zinterp(i,j) = 0.d0
do in = 1,ngnod
@@ -469,22 +470,22 @@
! tracer des droites si elements 4 noeuds
- ir=iptsdisp
+ ir=pointsdisp
x2 = (xinterp(ir,is)-xmin)*rapp_page + orig_x
z2 = (zinterp(ir,is)-zmin)*rapp_page + orig_z
x2 = x2 * centim
z2 = z2 * centim
write(24,681) x2,z2
- ir=iptsdisp
- is=iptsdisp
+ ir=pointsdisp
+ is=pointsdisp
x2 = (xinterp(ir,is)-xmin)*rapp_page + orig_x
z2 = (zinterp(ir,is)-zmin)*rapp_page + orig_z
x2 = x2 * centim
z2 = z2 * centim
write(24,681) x2,z2
- is=iptsdisp
+ is=pointsdisp
ir=1
x2 = (xinterp(ir,is)-xmin)*rapp_page + orig_x
z2 = (zinterp(ir,is)-zmin)*rapp_page + orig_z
@@ -503,7 +504,7 @@
else
! tracer des courbes si elements 9 noeuds
- do ir=2,iptsdisp
+ do ir=2,pointsdisp
x2 = (xinterp(ir,is)-xmin)*rapp_page + orig_x
z2 = (zinterp(ir,is)-zmin)*rapp_page + orig_z
x2 = x2 * centim
@@ -511,8 +512,8 @@
write(24,681) x2,z2
enddo
- ir=iptsdisp
- do is=2,iptsdisp
+ ir=pointsdisp
+ do is=2,pointsdisp
x2 = (xinterp(ir,is)-xmin)*rapp_page + orig_x
z2 = (zinterp(ir,is)-zmin)*rapp_page + orig_z
x2 = x2 * centim
@@ -520,8 +521,8 @@
write(24,681) x2,z2
enddo
- is=iptsdisp
- do ir=iptsdisp-1,1,-1
+ is=pointsdisp
+ do ir=pointsdisp-1,1,-1
x2 = (xinterp(ir,is)-xmin)*rapp_page + orig_x
z2 = (zinterp(ir,is)-zmin)*rapp_page + orig_z
x2 = x2 * centim
@@ -530,7 +531,7 @@
enddo
ir=1
- do is=iptsdisp-1,2,-1
+ do is=pointsdisp-1,2,-1
x2 = (xinterp(ir,is)-xmin)*rapp_page + orig_x
z2 = (zinterp(ir,is)-zmin)*rapp_page + orig_z
x2 = x2 * centim
@@ -662,7 +663,7 @@
! return if the maximum displacement equals zero (no source)
if(dispmax == 0.d0) then
- print *,' null displacement : returning !'
+ write(IOUT,*) ' null displacement : returning !'
return
endif
@@ -679,16 +680,23 @@
if(interpol) then
- print *,'Interpolating the vector field...'
+ write(IOUT,*) 'Interpolating the vector field...'
do ispec=1,nspec
! interpolation sur grille reguliere
- if(mod(ispec,1000) == 0) write(*,*) 'Interpolation uniform grid element ',ispec
+ if(mod(ispec,1000) == 0) write(IOUT,*) 'Interpolation uniform grid element ',ispec
- do i=1,iptsdisp
- do j=1,iptsdisp
+! option to plot only lowerleft corner value to avoid very large files if dense meshes
+ if(plot_lowerleft_corner_only) then
+ pointsdisp_loop = 1
+ else
+ pointsdisp_loop = pointsdisp
+ endif
+ do i=1,pointsdisp_loop
+ do j=1,pointsdisp_loop
+
xinterp(i,j) = 0.d0
zinterp(i,j) = 0.d0
do in = 1,ngnod
@@ -714,13 +722,13 @@
x1 =(xinterp(i,j)-xmin)*rapp_page
z1 =(zinterp(i,j)-zmin)*rapp_page
- x2 = Uxinterp(i,j)*sizemax/dispmax
- z2 = Uzinterp(i,j)*sizemax/dispmax
+ x2 = Uxinterp(i,j)*sizemax_arrows/dispmax
+ z2 = Uzinterp(i,j)*sizemax_arrows/dispmax
d = sqrt(x2**2 + z2**2)
! ignorer si vecteur trop petit
- if(d > cutvect*sizemax) then
+ if(d > cutvect*sizemax_arrows) then
d1 = d * rapport
d2 = d1 * cos(angle*convert)
@@ -779,13 +787,13 @@
x1 =(coord(1,ipoin)-xmin)*rapp_page
z1 =(coord(2,ipoin)-zmin)*rapp_page
- x2 = displ(1,ipoin)*sizemax/dispmax
- z2 = displ(2,ipoin)*sizemax/dispmax
+ x2 = displ(1,ipoin)*sizemax_arrows/dispmax
+ z2 = displ(2,ipoin)*sizemax_arrows/dispmax
d = sqrt(x2**2 + z2**2)
! ignorer si vecteur trop petit
- if(d > cutvect*sizemax) then
+ if(d > cutvect*sizemax_arrows) then
d1 = d * rapport
d2 = d1 * cos(angle*convert)
Modified: seismo/2D/SPECFEM2D/trunk/specfem2D.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/specfem2D.f90 2005-02-20 00:09:28 UTC (rev 8467)
+++ seismo/2D/SPECFEM2D/trunk/specfem2D.f90 2007-12-07 23:48:34 UTC (rev 8468)
@@ -28,6 +28,7 @@
! - more flexible DATA/Par_file with any number of comment lines
! - Xsu scripts for seismograms
! - subtract t0 from seismograms
+! - seismograms and snapshots in pressure in addition to velocity if acoustic medium
!
! version 5.0, May 2004 :
! - got rid of useless routines, suppressed commons etc.
@@ -104,7 +105,7 @@
double precision xixl,xizl,gammaxl,gammazl,jacobianl
! material properties of the elastic medium
- double precision mul_relaxed,lambdal_relaxed,lambdalplus2mul_relaxed,cpsquare
+ double precision mul_relaxed,lambdal_relaxed,lambdalplus2mul_relaxed,cpsquare,denst
double precision mul_unrelaxed,lambdal_unrelaxed,lambdalplus2mul_unrelaxed
double precision, dimension(:,:), allocatable :: coord,accel,veloc,displ, &
@@ -120,22 +121,23 @@
integer, dimension(:,:), allocatable :: knods
integer, dimension(:), allocatable :: kmato,numabs,numsurface
- integer ie,k
+ integer ie,k,material
integer ispec_selected_source,iglob_source,ix_source,iz_source
double precision a,displnorm_all
double precision, dimension(:), allocatable :: source_time_function
- double precision rsizemin,rsizemax,cpoverdxmin,cpoverdxmax, &
+ double precision rsizemin,rsizemax,cpoverdxmax, &
lambdaSmin,lambdaSmax,lambdaPmin,lambdaPmax,vpmin,vpmax
- integer colors,numbers,subsamp,vecttype,itaff,nrec,sismostype
- integer numat,ngnod,nspec,iptsdisp,nelemabs,nelemsurface
+ integer colors,numbers,subsamp,vecttype,IT_AFFICHE,nrec,sismostype
+ integer numat,ngnod,nspec,pointsdisp,nelemabs,nelemsurface
logical interpol,meshvect,modelvect,boundvect,read_external_model,initialfield,abshaut, &
- outputgrid,gnuplot,ELASTIC,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON,output_postscript_snapshot,output_PNM_image
+ outputgrid,gnuplot,ELASTIC,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON,output_postscript_snapshot,output_PNM_image, &
+ plot_lowerleft_corner_only
- double precision cutvect,anglerec,xirec,gammarec
+ double precision cutvect,sizemax_arrows,anglerec,xirec,gammarec
! for absorbing and free surface conditions
integer ispecabs,ispecsurface,inum,numabsread,numsurfaceread,i1abs,i2abs
@@ -211,13 +213,13 @@
!
call datim(stitle)
- write(*,*)
- write(*,*)
- write(*,*) '*********************'
- write(*,*) '**** ****'
- write(*,*) '**** SPECFEM2D ****'
- write(*,*) '**** ****'
- write(*,*) '*********************'
+ write(IOUT,*)
+ write(IOUT,*)
+ write(IOUT,*) '*********************'
+ write(IOUT,*) '**** ****'
+ write(IOUT,*) '**** SPECFEM2D ****'
+ write(IOUT,*) '**** ****'
+ write(IOUT,*) '*********************'
!
!---- read parameters from input file
@@ -230,10 +232,10 @@
read(IIN,*) gnuplot,interpol
read(IIN,40) datlin
- read(IIN,*) itaff,output_postscript_snapshot,output_PNM_image,colors,numbers
+ read(IIN,*) IT_AFFICHE,output_postscript_snapshot,output_PNM_image,colors,numbers
read(IIN,40) datlin
- read(IIN,*) meshvect,modelvect,boundvect,cutvect,subsamp,nx_sem_PNM
+ read(IIN,*) meshvect,modelvect,boundvect,cutvect,subsamp,sizemax_arrows,nx_sem_PNM
cutvect = cutvect / 100.d0
read(IIN,40) datlin
@@ -250,7 +252,7 @@
!---- check parameters read
write(IOUT,200) npgeo,NDIM
- write(IOUT,600) itaff,colors,numbers
+ write(IOUT,600) IT_AFFICHE,colors,numbers
write(IOUT,700) sismostype,anglerec
write(IOUT,750) initialfield,read_external_model,ELASTIC,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON,outputgrid
write(IOUT,800) vecttype,100.d0*cutvect,subsamp
@@ -304,25 +306,25 @@
!---- read the basic properties of the spectral elements
!
read(IIN,40) datlin
- read(IIN,*) numat,ngnod,nspec,iptsdisp,nelemabs,nelemsurface
+ read(IIN,*) numat,ngnod,nspec,pointsdisp,plot_lowerleft_corner_only,nelemabs,nelemsurface
!
!---- allocate arrays
!
allocate(shape2D(ngnod,NGLLX,NGLLZ))
allocate(dershape2D(NDIM,ngnod,NGLLX,NGLLZ))
- allocate(shape2D_display(ngnod,iptsdisp,iptsdisp))
- allocate(dershape2D_display(NDIM,ngnod,iptsdisp,iptsdisp))
+ allocate(shape2D_display(ngnod,pointsdisp,pointsdisp))
+ allocate(dershape2D_display(NDIM,ngnod,pointsdisp,pointsdisp))
allocate(xix(NGLLX,NGLLZ,nspec))
allocate(xiz(NGLLX,NGLLZ,nspec))
allocate(gammax(NGLLX,NGLLZ,nspec))
allocate(gammaz(NGLLX,NGLLZ,nspec))
allocate(jacobian(NGLLX,NGLLZ,nspec))
- allocate(flagrange(NGLLX,iptsdisp))
- allocate(xinterp(iptsdisp,iptsdisp))
- allocate(zinterp(iptsdisp,iptsdisp))
- allocate(Uxinterp(iptsdisp,iptsdisp))
- allocate(Uzinterp(iptsdisp,iptsdisp))
+ allocate(flagrange(NGLLX,pointsdisp))
+ allocate(xinterp(pointsdisp,pointsdisp))
+ allocate(zinterp(pointsdisp,pointsdisp))
+ allocate(Uxinterp(pointsdisp,pointsdisp))
+ allocate(Uzinterp(pointsdisp,pointsdisp))
allocate(density(numat))
allocate(elastcoef(4,numat))
allocate(kmato(nspec))
@@ -380,7 +382,7 @@
!---- print element group main parameters
!
write(IOUT,107)
- write(IOUT,207) nspec,ngnod,NGLLX,NGLLZ,NGLLX*NGLLZ,iptsdisp,numat,nelemabs
+ write(IOUT,207) nspec,ngnod,NGLLX,NGLLZ,NGLLX*NGLLZ,pointsdisp,numat,nelemabs
! set up Gauss-Lobatto-Legendre derivation matrices
call define_derivative_matrices(xigll,zigll,wxgll,wzgll,hprime_xx,hprime_zz)
@@ -413,8 +415,8 @@
codeabs(ILEFT,inum) = codeabsread(3)
codeabs(IRIGHT,inum) = codeabsread(4)
enddo
- write(*,*)
- write(*,*) 'Number of absorbing elements: ',nelemabs
+ write(IOUT,*)
+ write(IOUT,*) 'Number of absorbing elements: ',nelemabs
endif
!
@@ -427,8 +429,8 @@
if(inum < 1 .or. inum > nelemsurface) stop 'Wrong free surface element number'
numsurface(inum) = numsurfaceread
enddo
- write(*,*)
- write(*,*) 'Number of free surface elements: ',nelemsurface
+ write(IOUT,*)
+ write(IOUT,*) 'Number of free surface elements: ',nelemsurface
!
!---- close input file
@@ -456,10 +458,10 @@
endif
!---- compute shape functions and their derivatives for regular !interpolated display grid
- do j = 1,iptsdisp
- do i = 1,iptsdisp
- xirec = 2.d0*dble(i-1)/dble(iptsdisp-1) - 1.d0
- gammarec = 2.d0*dble(j-1)/dble(iptsdisp-1) - 1.d0
+ do j = 1,pointsdisp
+ do i = 1,pointsdisp
+ xirec = 2.d0*dble(i-1)/dble(pointsdisp-1) - 1.d0
+ gammarec = 2.d0*dble(j-1)/dble(pointsdisp-1) - 1.d0
call define_shape_functions(shape2D_display(:,i,j),dershape2D_display(:,:,i,j),xirec,gammarec,ngnod)
enddo
enddo
@@ -467,13 +469,12 @@
!---- compute Lagrange interpolants on a regular interpolated grid in (xi,gamma)
!---- for display (assumes NGLLX = NGLLZ)
do j=1,NGLLX
- do i=1,iptsdisp
- xirec = 2.d0*dble(i-1)/dble(iptsdisp-1) - 1.d0
+ do i=1,pointsdisp
+ xirec = 2.d0*dble(i-1)/dble(pointsdisp-1) - 1.d0
flagrange(j,i) = hgll(j-1,xirec,xigll,NGLLX)
enddo
enddo
-
! read total number of receivers
open(unit=IIN,file='DATA/STATIONS',status='old')
read(IIN,*) nrec
@@ -565,9 +566,9 @@
!--- save the grid of points in a file
!
if(outputgrid) then
- print *
- print *,'Saving the grid in a text file...'
- print *
+ write(IOUT,*)
+ write(IOUT,*) 'Saving the grid in a text file...'
+ write(IOUT,*)
open(unit=55,file='OUTPUT_FILES/grid_points_and_model.txt',status='unknown')
write(55,*) npoin
do n = 1,npoin
@@ -623,9 +624,9 @@
!---- eventuellement lecture d'un modele externe de vitesse et de densite
!
if(read_external_model) then
- print *
- print *,'Reading velocity and density model from external file...'
- print *
+ write(IOUT,*)
+ write(IOUT,*) 'Reading velocity and density model from external file...'
+ write(IOUT,*)
open(unit=55,file='OUTPUT_FILES/grid_points_and_model.txt',status='old')
read(55,*) nbpoin
if(nbpoin /= npoin) stop 'Wrong number of points in input file'
@@ -640,7 +641,7 @@
!
call defarrays(vpext,vsext,rhoext,density,elastcoef, &
ibool,kmato,coord,npoin,rsizemin,rsizemax, &
- cpoverdxmin,cpoverdxmax,lambdaSmin,lambdaSmax,lambdaPmin,lambdaPmax, &
+ cpoverdxmax,lambdaSmin,lambdaSmax,lambdaPmin,lambdaPmax, &
vpmin,vpmax,read_external_model,nspec,numat)
! build the global mass matrix once and for all
@@ -675,7 +676,7 @@
!---- verifier le maillage, la stabilite et le nb de points par lambda
!---- seulement si la source en temps n'est pas un Dirac (sinon spectre non defini)
if(time_function_type /= 4) call checkgrid(deltat,f0,t0,initialfield, &
- rsizemin,rsizemax,cpoverdxmin,cpoverdxmax,lambdaSmin,lambdaSmax,lambdaPmin,lambdaPmax)
+ rsizemin,rsizemax,cpoverdxmax,lambdaSmin,lambdaSmax,lambdaPmin,lambdaPmax)
!
!---- for color PNM images
@@ -704,8 +705,8 @@
allocate(copy_iglob_image_PNM_2D(NX_IMAGE_PNM,NZ_IMAGE_PNM))
! creer tous les pixels
- print *
- print *,'localisation de tous les pixels des images PNM'
+ write(IOUT,*)
+ write(IOUT,*) 'localisation de tous les pixels des images PNM'
taille_pixel_horizontal = (xmax_PNM_image - xmin_PNM_image) / dble(NX_IMAGE_PNM-1)
taille_pixel_vertical = (zmax_PNM_image - zmin_PNM_image) / dble(NZ_IMAGE_PNM-1)
@@ -790,7 +791,7 @@
deallocate(copy_iglob_image_PNM_2D)
- print *,'fin localisation de tous les pixels des images PNM'
+ write(IOUT,*) 'fin localisation de tous les pixels des images PNM'
!
!---- initialiser sismogrammes
@@ -810,9 +811,9 @@
!---- eventuellement lecture des champs initiaux dans un fichier
!
if(initialfield) then
- print *
- print *,'Reading initial fields from external file...'
- print *
+ write(IOUT,*)
+ write(IOUT,*) 'Reading initial fields from external file...'
+ write(IOUT,*)
open(unit=55,file='OUTPUT_FILES/wavefields.txt',status='unknown')
read(55,*) nbpoin
if(nbpoin /= npoin) stop 'Wrong number of points in input file'
@@ -831,7 +832,7 @@
deallocate(velocread)
deallocate(accelread)
close(55)
- print *,'Max norm of initial displacement = ',maxval(sqrt(displ(1,:)**2 + displ(2,:)**2))
+ write(IOUT,*) 'Max norm of initial displacement = ',maxval(sqrt(displ(1,:)**2 + displ(2,:)**2))
endif
! attenuation constants from Carcione 1993 Geophysics volume 58 pages 111 and 112
@@ -872,9 +873,9 @@
allocate(source_time_function(NSTEP))
- print *
- print *,'Saving the source time function in a text file...'
- print *
+ write(IOUT,*)
+ write(IOUT,*) 'Saving the source time function in a text file...'
+ write(IOUT,*)
open(unit=55,file='OUTPUT_FILES/source.txt',status='unknown')
! boucle principale d'evolution en temps
@@ -923,15 +924,6 @@
! compute current time
time = (it-1)*deltat
- if(mod(it,itaff) == 0) then
- write(IOUT,*)
- if(time >= 1.d-3) then
- write(IOUT,100) it,time
- else
- write(IOUT,101) it,time
- endif
- endif
-
! compute Grad(displ) at time step n for attenuation
if(TURN_ATTENUATION_ON) call compute_gradient_attenuation(displ,duxdxl_n,duzdxl_n, &
duxdzl_n,duzdzl_n,xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,NSPEC,npoin)
@@ -1540,17 +1532,29 @@
endif ! end of test on attenuation
-!---- display max of norm of displacement
- if(mod(it,itaff) == 0) then
+!---- display time step max of norm of displacement
+ if(mod(it,IT_AFFICHE) == 0 .or. it == 5) then
+
+ write(IOUT,*)
+ if(time >= 1.d-3) then
+ write(IOUT,100) it,time
+ else
+ write(IOUT,101) it,time
+ endif
+
displnorm_all = maxval(sqrt(displ(1,:)**2 + displ(2,:)**2))
- print *,'Max norm of field = ',displnorm_all
+ write(IOUT,*) 'Max norm of field = ',displnorm_all
! check stability of the code, exit if unstable
if(displnorm_all > STABILITY_THRESHOLD) stop 'code became unstable and blew up'
+
+ write(IOUT,*)
endif
! store the seismograms
- if(sismostype < 1 .or. sismostype > 3) stop 'Wrong field code for seismogram output'
+ if(sismostype < 1 .or. sismostype > 4) stop 'Wrong field code for seismogram output'
+ if(ELASTIC .and. sismostype == 4) stop 'pressure seismograms implemented for an acoustic medium only'
+
if(.not. ELASTIC) then
if(sismostype == 1) then
stop 'cannot store displacement field in acoustic medium because of potential formulation'
@@ -1558,7 +1562,7 @@
! for acoustic medium, compute gradient for display, displ represents the potential
call compute_gradient_fluid(displ,vector_field_postscript, &
xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,NSPEC,npoin)
- else
+ else if(sismostype == 3) then
! for acoustic medium, compute gradient for display, veloc represents the first derivative of the potential
call compute_gradient_fluid(veloc,vector_field_postscript, &
xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,NSPEC,npoin)
@@ -1594,9 +1598,24 @@
else
! for acoustic medium
- dxd = vector_field_postscript(1,iglob)
- dzd = vector_field_postscript(2,iglob)
+! pressure = - rho * Chi_dot
+ if(sismostype == 4) then
+
+ material = kmato(ispec_selected_rec(irec))
+ denst = density(material)
+ if(read_external_model) denst = rhoext(ibool(i,j,ispec_selected_rec(irec)))
+
+ dxd = - denst * veloc(1,iglob)
+ dzd = ZERO
+
+! velocity
+ else
+ dxd = vector_field_postscript(1,iglob)
+ dzd = vector_field_postscript(2,iglob)
+
+ endif
+
endif
! compute interpolated field
@@ -1615,16 +1634,8 @@
!
!---- affichage des resultats a certains pas de temps
!
- if(mod(it,itaff) == 0 .or. it == 5 .or. it == NSTEP) then
+ if(mod(it,IT_AFFICHE) == 0 .or. it == 5 .or. it == NSTEP) then
- write(IOUT,*)
- if(time >= 1.d-3) then
- write(IOUT,110) time
- else
- write(IOUT,111) time
- endif
- write(IOUT,*)
-
!
!---- affichage postscript
!
@@ -1640,7 +1651,8 @@
Uxinterp,Uzinterp,flagrange,density,elastcoef,knods,kmato,ibool, &
numabs,codeabs,anyabs,stitle,npoin,npgeo,vpmin,vpmax,nrec, &
colors,numbers,subsamp,vecttype,interpol,meshvect,modelvect, &
- boundvect,read_external_model,cutvect,nelemabs,numat,iptsdisp,nspec,ngnod,ELASTIC)
+ boundvect,read_external_model,cutvect,sizemax_arrows,nelemabs,numat,pointsdisp,nspec,ngnod,ELASTIC, &
+ plot_lowerleft_corner_only)
else if(ELASTIC .and. vecttype == 2) then
write(IOUT,*) 'drawing velocity field...'
@@ -1649,7 +1661,8 @@
Uxinterp,Uzinterp,flagrange,density,elastcoef,knods,kmato,ibool, &
numabs,codeabs,anyabs,stitle,npoin,npgeo,vpmin,vpmax,nrec, &
colors,numbers,subsamp,vecttype,interpol,meshvect,modelvect, &
- boundvect,read_external_model,cutvect,nelemabs,numat,iptsdisp,nspec,ngnod,ELASTIC)
+ boundvect,read_external_model,cutvect,sizemax_arrows,nelemabs,numat,pointsdisp,nspec,ngnod,ELASTIC, &
+ plot_lowerleft_corner_only)
else if(ELASTIC .and. vecttype == 3) then
write(IOUT,*) 'drawing acceleration field...'
@@ -1658,7 +1671,8 @@
Uxinterp,Uzinterp,flagrange,density,elastcoef,knods,kmato,ibool, &
numabs,codeabs,anyabs,stitle,npoin,npgeo,vpmin,vpmax,nrec, &
colors,numbers,subsamp,vecttype,interpol,meshvect,modelvect, &
- boundvect,read_external_model,cutvect,nelemabs,numat,iptsdisp,nspec,ngnod,ELASTIC)
+ boundvect,read_external_model,cutvect,sizemax_arrows,nelemabs,numat,pointsdisp,nspec,ngnod,ELASTIC, &
+ plot_lowerleft_corner_only)
! for acoustic medium
else if(.not. ELASTIC .and. vecttype == 1) then
@@ -1674,7 +1688,8 @@
Uxinterp,Uzinterp,flagrange,density,elastcoef,knods,kmato,ibool, &
numabs,codeabs,anyabs,stitle,npoin,npgeo,vpmin,vpmax,nrec, &
colors,numbers,subsamp,vecttype,interpol,meshvect,modelvect, &
- boundvect,read_external_model,cutvect,nelemabs,numat,iptsdisp,nspec,ngnod,ELASTIC)
+ boundvect,read_external_model,cutvect,sizemax_arrows,nelemabs,numat,pointsdisp,nspec,ngnod,ELASTIC, &
+ plot_lowerleft_corner_only)
else if(.not. ELASTIC .and. vecttype == 3) then
write(IOUT,*) 'drawing acoustic acceleration field from velocity potential...'
@@ -1686,7 +1701,8 @@
Uxinterp,Uzinterp,flagrange,density,elastcoef,knods,kmato,ibool, &
numabs,codeabs,anyabs,stitle,npoin,npgeo,vpmin,vpmax,nrec, &
colors,numbers,subsamp,vecttype,interpol,meshvect,modelvect, &
- boundvect,read_external_model,cutvect,nelemabs,numat,iptsdisp,nspec,ngnod,ELASTIC)
+ boundvect,read_external_model,cutvect,sizemax_arrows,nelemabs,numat,pointsdisp,nspec,ngnod,ELASTIC, &
+ plot_lowerleft_corner_only)
else
stop 'wrong field code for PostScript display'
@@ -1707,7 +1723,7 @@
do j = 1,NZ_IMAGE_PNM
do i = 1,NX_IMAGE_PNM
if(iglob_image_PNM_2D(i,j) /= -1) then
-! display vertical component of vector
+! display vertical component of vector if elastic medium
if(ELASTIC) then
if(vecttype == 1) then
donnees_image_PNM_2D(i,j) = displ(2,iglob_image_PNM_2D(i,j))
@@ -1717,8 +1733,9 @@
donnees_image_PNM_2D(i,j) = accel(2,iglob_image_PNM_2D(i,j))
endif
else
-! for acoustic medium
+! display pressure if acoustic medium
donnees_image_PNM_2D(i,j) = vector_field_postscript(2,iglob_image_PNM_2D(i,j))
+stop 'DK DK ceci a debugger par DK, mettre pressure au lieu de veloc'
endif
endif
enddo
@@ -1764,7 +1781,7 @@
'Number of spectral element control nodes. . .(npgeo) =',i8/5x, &
'Number of space dimensions. . . . . . . . . . (NDIM) =',i8)
600 format(//1x,'C o n t r o l',/1x,13('='),//5x, &
- 'Display frequency . . . . . . . . . . . . . (itaff) = ',i5/ 5x, &
+ 'Display frequency . . . . . . . . . . . (IT_AFFICHE) = ',i5/ 5x, &
'Color display . . . . . . . . . . . . . . . (colors) = ',i5/ 5x, &
' == 0 black and white display ', / 5x, &
' == 1 color display ', /5x, &
@@ -1798,7 +1815,7 @@
'Number of points in X-direction . . . (NGLLX) =',i7,/5x, &
'Number of points in Y-direction . . . (NGLLZ) =',i7,/5x, &
'Number of points per element. . .(NGLLX*NGLLZ) =',i7,/5x, &
- 'Number of points for display . . . .(iptsdisp) =',i7,/5x, &
+ 'Number of points for display . . .(pointsdisp) =',i7,/5x, &
'Number of element material sets . . . (numat) =',i7,/5x, &
'Number of absorbing elements . . . .(nelemabs) =',i7)
Modified: seismo/2D/SPECFEM2D/trunk/write_seismograms.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/write_seismograms.f90 2005-02-20 00:09:28 UTC (rev 8467)
+++ seismo/2D/SPECFEM2D/trunk/write_seismograms.f90 2007-12-07 23:48:34 UTC (rev 8468)
@@ -29,14 +29,15 @@
character(len=MAX_LENGTH_STATION_NAME), dimension(nrec) :: station_name
character(len=MAX_LENGTH_NETWORK_NAME), dimension(nrec) :: network_name
- integer irec,irecord,length_station_name,length_network_name,iorientation,isample
+ integer irec,irecord,length_station_name,length_network_name,iorientation,isample,number_of_components
character(len=4) chn
character(len=1) component
character(len=150) sisname
-! to write seismograms in single precision SEP binary format
- real(kind=4), dimension(NSTEP*nrec) :: buffer_SEP_binary
+! to write seismograms in single precision SEP and double precision binary format
+ real(kind=4), dimension(NSTEP*nrec) :: buffer_binary_single
+ double precision, dimension(NSTEP*nrec) :: buffer_binary_double
! scaling factor for Seismic Unix xsu dislay
double precision, parameter :: FACTORXSU = 1.d0
@@ -45,21 +46,30 @@
! write seismograms in ASCII format
-! save displacement, velocity or acceleration
+! save displacement, velocity, acceleration or pressure
if(sismostype == 1) then
component = 'd'
else if(sismostype == 2) then
component = 'v'
else if(sismostype == 3) then
component = 'a'
+ else if(sismostype == 4) then
+ component = 'p'
else
stop 'wrong component to save for seismograms'
endif
do irec = 1,nrec
- do iorientation = 1,NDIM
+! only one seismogram if pressurs
+ if(sismostype == 4) then
+ number_of_components = 1
+ else
+ number_of_components = NDIM
+ endif
+ do iorientation = 1,number_of_components
+
if(iorientation == 1) then
chn = 'BHX'
else if(iorientation == 2) then
@@ -68,6 +78,9 @@
stop 'incorrect channel value'
endif
+! in case of pressure, use different abbreviation
+ if(sismostype == 4) chn = 'PRE'
+
! create the name of the seismogram file for each slice
! file name includes the name of the station, the network and the component
length_station_name = len_trim(station_name(irec))
@@ -109,41 +122,78 @@
! write seismograms in single precision SEP binary format
! X component
+
+! delete the old files
+ open(unit=11,file='OUTPUT_FILES/Ux_file_single.bin',status='unknown')
+ close(11,status='delete')
+
+ open(unit=11,file='OUTPUT_FILES/Ux_file_double.bin',status='unknown')
+ close(11,status='delete')
+
+ open(unit=11,file='OUTPUT_FILES/pressure_file_single.bin',status='unknown')
+ close(11,status='delete')
+
+ open(unit=11,file='OUTPUT_FILES/pressure_file_double.bin',status='unknown')
+ close(11,status='delete')
+
irecord = 0
do irec=1,nrec
do isample=1,NSTEP
irecord = irecord + 1
- buffer_SEP_binary(irecord) = sngl(sisux(isample,irec))
+ buffer_binary_single(irecord) = sngl(sisux(isample,irec))
+ buffer_binary_double(irecord) = sisux(isample,irec)
enddo
enddo
-! delete the old file
- open(unit=11,file='OUTPUT_FILES/Ux_file.bin',status='unknown')
- close(11,status='delete')
+! write the new files
+ if(sismostype == 4) then
+ open(unit=11,file='OUTPUT_FILES/pressure_file_single.bin',status='unknown',access='direct',recl=4*NSTEP*nrec)
+ else
+ open(unit=11,file='OUTPUT_FILES/Ux_file_single.bin',status='unknown',access='direct',recl=4*NSTEP*nrec)
+ endif
+ write(11,rec=1) buffer_binary_single
+ close(11)
-! write the new file
- open(unit=11,file='OUTPUT_FILES/Ux_file.bin',status='unknown',access='direct',recl=NSTEP*nrec)
- write(11,rec=1) buffer_SEP_binary
+ if(sismostype == 4) then
+ open(unit=11,file='OUTPUT_FILES/pressure_file_double.bin',status='unknown',access='direct',recl=8*NSTEP*nrec)
+ else
+ open(unit=11,file='OUTPUT_FILES/Ux_file_double.bin',status='unknown',access='direct',recl=8*NSTEP*nrec)
+ endif
+ write(11,rec=1) buffer_binary_double
close(11)
! Z component
+
+! delete the old files
+ open(unit=11,file='OUTPUT_FILES/Uz_file_single.bin',status='unknown')
+ close(11,status='delete')
+
+ open(unit=11,file='OUTPUT_FILES/Uz_file_double.bin',status='unknown')
+ close(11,status='delete')
+
+! no Z component seismogram if pressurs
+ if(sismostype /= 4) then
+
irecord = 0
do irec=1,nrec
do isample=1,NSTEP
irecord = irecord + 1
- buffer_SEP_binary(irecord) = sngl(sisuz(isample,irec))
+ buffer_binary_single(irecord) = sngl(sisuz(isample,irec))
+ buffer_binary_double(irecord) = sisuz(isample,irec)
enddo
enddo
-! delete the old file
- open(unit=11,file='OUTPUT_FILES/Uz_file.bin',status='unknown')
- close(11,status='delete')
+! write the new files
+ open(unit=11,file='OUTPUT_FILES/Uz_file_single.bin',status='unknown',access='direct',recl=4*NSTEP*nrec)
+ write(11,rec=1) buffer_binary_single
+ close(11)
-! write the new file
- open(unit=11,file='OUTPUT_FILES/Uz_file.bin',status='unknown',access='direct',recl=NSTEP*nrec)
- write(11,rec=1) buffer_SEP_binary
+ open(unit=11,file='OUTPUT_FILES/Uz_file_double.bin',status='unknown',access='direct',recl=8*NSTEP*nrec)
+ write(11,rec=1) buffer_binary_double
close(11)
+ endif
+
!----
! ligne de recepteurs pour Xsu
@@ -158,11 +208,11 @@
enddo
if(sismostype == 1) then
- write(11,*) '@title="Ux at displacement@component"@<@Ux_file.bin'
+ write(11,*) '@title="Ux at displacement@component"@<@Ux_file_single.bin'
else if(sismostype == 2) then
- write(11,*) '@title="Ux at velocity@component"@<@Ux_file.bin'
+ write(11,*) '@title="Ux at velocity@component"@<@Ux_file_single.bin'
else
- write(11,*) '@title="Ux at acceleration@component"@<@Ux_file.bin'
+ write(11,*) '@title="Ux at acceleration@component"@<@Ux_file_single.bin'
endif
close(11)
@@ -199,8 +249,8 @@
120 format('sed -e ''1,$s/ //g'' -e ''1,$s/@/ /g'' -e ''1,1p'' -e ''$,$s/Ux/Uz/g'' <tempfile > receiver_line_Xsu_XWindow')
130 format('sed -e ''1,$s/xwigb/pswigp/g'' ', &
- '-e ''1,$s/Ux_file.bin/Ux_file.bin > uxpoly.ps/g'' ', &
- '-e ''1,$s/Uz_file.bin/Uz_file.bin > uzpoly.ps/g'' receiver_line_Xsu_XWindow > receiver_line_Xsu_postscript')
+ '-e ''1,$s/Ux_file_single.bin/Ux_file_single.bin > uxpoly.ps/g'' ', &
+ '-e ''1,$s/Uz_file_single.bin/Uz_file_single.bin > uzpoly.ps/g'' receiver_line_Xsu_XWindow > receiver_line_Xsu_postscript')
140 format(f12.5)
More information about the cig-commits
mailing list