[cig-commits] r8472 - in seismo/2D/SPECFEM2D/trunk: . DATA
walter at geodynamics.org
walter at geodynamics.org
Fri Dec 7 15:48:55 PST 2007
Author: walter
Date: 2007-12-07 15:48:54 -0800 (Fri, 07 Dec 2007)
New Revision: 8472
Modified:
seismo/2D/SPECFEM2D/trunk/DATA/Par_file
seismo/2D/SPECFEM2D/trunk/DATA/Par_file_CercleMod50_50_10Hz_Paul
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_Thomas_Lacq
seismo/2D/SPECFEM2D/trunk/DATA/Par_file_cours_M2_UPPA
seismo/2D/SPECFEM2D/trunk/meshfem2D.f90
seismo/2D/SPECFEM2D/trunk/read_value_parameters.f90
seismo/2D/SPECFEM2D/trunk/specfem2D.f90
seismo/2D/SPECFEM2D/trunk/su_create_Paul.python
Log:
added pressure images in PNM format in case of acoustic medium.
added option to have several receiver antennas instead of one.
added option to have only one receiver per antenna if needed (minimum used to be 2).
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file 2005-02-23 17:33:20 UTC (rev 8471)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file 2007-12-07 23:48:54 UTC (rev 8472)
@@ -1,21 +1,19 @@
#
# 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
+title = Test sinusoide acoustique
+interfacesfile = interface_sinus.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
+xmax = 4000.d0 # abscissa of right side of the model
+nx = 80 # 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
+readmodel = .false. # read external earth model or not
+ELASTIC = .false. # elastic or acoustic simulation
TURN_ANISOTROPY_ON = .false. # turn anisotropy on or off
TURN_ATTENUATION_ON = .false. # turn attenuation on or off
@@ -26,53 +24,57 @@
absorbdroite = .true.
# time step parameters
-nt = 10000 # nb total de pas de temps
-dt = 1.d-3 # valeur du pas de temps
+nt = 1400 # nb total de pas de temps
+deltat = 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 = 5 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4, Heaviside = 5
+xs = 2000. # source location x in meters
+zs = 1200. # source location z in meters
+source_type = 1 # force = 1 or moment tensor = 2
+time_function_type = 2 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4, Heaviside = 5
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
+factor = 1.d4 # 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
+sismostype = 4 # record 1=displ 2=veloc 3=accel 4=pressure
+nreceiverlines = 1 # number of receiver lines
anglerec = 0.d0 # angle to rotate components at receivers
+# first receiver line
+nrec = 101 # number of receivers
+xdeb = 300. # first receiver x in meters
+zdeb = 1700. # first receiver z in meters
+xfin = 3700. # last receiver x in meters (ignored if onlyone receiver)
+zfin = 1700. # last receiver z in meters (ignored if onlyone receiver)
+enreg_surf = .false. # enregistrement volume ou surface
+
# display parameters
-itaff = 250 # display frequency in time steps
+itaff = 100 # 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
+output_PNM_image = .true. # output PNM image of the results
+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 = 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
+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 = 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
+nbmodels = 2 # nb de modeles differents (0,rho,vp,vs,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
#1 80 41 60 3
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file_CercleMod50_50_10Hz_Paul
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file_CercleMod50_50_10Hz_Paul 2005-02-23 17:33:20 UTC (rev 8471)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file_CercleMod50_50_10Hz_Paul 2007-12-07 23:48:54 UTC (rev 8472)
@@ -1,8 +1,6 @@
#
# 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 = Test InterfaceCercleMod50_50
@@ -27,13 +25,13 @@
# time step parameters
nt = 6600 # nb total de pas de temps
-dt = 2.d-4 # valeur du pas de temps
+deltat = 2.d-4 # valeur du pas de temps
# source parameters
source_surf = .false. # source dans le volume ou a la surface
xs = 2000. # source location x in meters
zs = 2850. # source location z in meters
-source_type = 1 # force = 1 or explosion = 2
+source_type = 1 # force = 1 or moment tensor = 2
time_function_type = 1 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4, Heaviside = 5
f0 = 10.0 # dominant source frequency (Hz) if not Dirac
angleforce = 0. # angle of the source (for a force only)
@@ -43,14 +41,17 @@
factor = 1.d-6 # amplification factor
# receiver line parameters
-enreg_surf = .false. # enregistrement volume ou surface
sismostype = 2 # record 1=displ 2=veloc 3=accel 4=pressure
+nreceiverlines = 1 # number of receiver lines
+anglerec = 0.d0 # angle to rotate components at receivers
+
+# first receiver line
nrec = 101 # number of receivers
xdeb = 500. # first receiver x in meters
zdeb = 2750. # first receiver z in meters
-xfin = 3500. # last receiver x in meters
-zfin = 2750. # last receiver z in meters
-anglerec = 0.d0 # angle to rotate components at receivers
+xfin = 3500. # last receiver x in meters (ignored if onlyone receiver)
+zfin = 2750. # last receiver z in meters (ignored if onlyone receiver)
+enreg_surf = .false. # enregistrement volume ou surface
# display parameters
itaff = 200 # display frequency in time steps
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Paul_acoustic
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Paul_acoustic 2005-02-23 17:33:20 UTC (rev 8471)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Paul_acoustic 2007-12-07 23:48:54 UTC (rev 8472)
@@ -1,8 +1,6 @@
#
# 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 = Test sinusoide acoustique
@@ -27,13 +25,13 @@
# time step parameters
nt = 1400 # nb total de pas de temps
-dt = 1.d-3 # valeur du pas de temps
+deltat = 1.d-3 # valeur du pas de temps
# source parameters
source_surf = .false. # source dans le volume ou a la surface
xs = 2000. # source location x in meters
zs = 1200. # source location z in meters
-source_type = 1 # force = 1 or explosion = 2
+source_type = 1 # force = 1 or moment tensor = 2
time_function_type = 2 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4, Heaviside = 5
f0 = 10.0 # dominant source frequency (Hz) if not Dirac
angleforce = 0. # angle of the source (for a force only)
@@ -43,14 +41,17 @@
factor = 1.d4 # amplification factor
# receiver line parameters
-enreg_surf = .false. # enregistrement volume ou surface
sismostype = 4 # record 1=displ 2=veloc 3=accel 4=pressure
+nreceiverlines = 1 # number of receiver lines
+anglerec = 0.d0 # angle to rotate components at receivers
+
+# first receiver line
nrec = 101 # number of receivers
xdeb = 300. # first receiver x in meters
zdeb = 1700. # first receiver z in meters
-xfin = 3700. # last receiver x in meters
-zfin = 1700. # last receiver z in meters
-anglerec = 0.d0 # angle to rotate components at receivers
+xfin = 3700. # last receiver x in meters (ignored if onlyone receiver)
+zfin = 1700. # last receiver z in meters (ignored if onlyone receiver)
+enreg_surf = .false. # enregistrement volume ou surface
# display parameters
itaff = 100 # display frequency in time steps
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Paul_elastic
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Paul_elastic 2005-02-23 17:33:20 UTC (rev 8471)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Paul_elastic 2007-12-07 23:48:54 UTC (rev 8472)
@@ -1,8 +1,6 @@
#
# 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 = Test sinusoide elastique
@@ -27,13 +25,13 @@
# time step parameters
nt = 1400 # nb total de pas de temps
-dt = 1.d-3 # valeur du pas de temps
+deltat = 1.d-3 # valeur du pas de temps
# source parameters
source_surf = .false. # source dans le volume ou a la surface
xs = 2000. # source location x in meters
zs = 1200. # source location z in meters
-source_type = 1 # force = 1 or explosion = 2
+source_type = 1 # force = 1 or moment tensor = 2
time_function_type = 1 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4, Heaviside = 5
f0 = 10.0 # dominant source frequency (Hz) if not Dirac
angleforce = 0. # angle of the source (for a force only)
@@ -43,14 +41,17 @@
factor = 1.d12 # amplification factor
# receiver line parameters
-enreg_surf = .true. # enregistrement volume ou surface
sismostype = 1 # record 1=displ 2=veloc 3=accel 4=pressure
+nreceiverlines = 1 # number of receiver lines
+anglerec = 0.d0 # angle to rotate components at receivers
+
+# first receiver line
nrec = 101 # number of receivers
xdeb = 300. # first receiver x in meters
zdeb = 2200. # first receiver z in meters
-xfin = 3700. # last receiver x in meters
-zfin = 2200. # last receiver z in meters
-anglerec = 0.d0 # angle to rotate components at receivers
+xfin = 3700. # last receiver x in meters (ignored if onlyone receiver)
+zfin = 2200. # last receiver z in meters (ignored if onlyone receiver)
+enreg_surf = .true. # enregistrement volume ou surface
# display parameters
itaff = 100 # display frequency in time steps
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Thomas_Lacq
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Thomas_Lacq 2005-02-23 17:33:20 UTC (rev 8471)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file_Thomas_Lacq 2007-12-07 23:48:54 UTC (rev 8472)
@@ -1,8 +1,6 @@
#
# 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
@@ -27,13 +25,13 @@
# time step parameters
nt = 10000 # nb total de pas de temps
-dt = 1.d-3 # valeur du pas de temps
+deltat = 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
+source_type = 1 # force = 1 or moment tensor = 2
time_function_type = 5 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4, Heaviside = 5
f0 = 10.0 # dominant source frequency (Hz) if not Dirac
angleforce = 0. # angle of the source (for a force only)
@@ -43,15 +41,26 @@
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
+sismostype = 2 # record 1=displ 2=veloc 3=accel 4=pressure
+nreceiverlines = 2 # number of receiver lines
+anglerec = 0.d0 # angle to rotate components at receivers
+
+# first receiver line
+nrec = 4 # number of receivers (can be 1 or more)
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
+zdeb = 5000. # first receiver z in meters
+xfin = 20000. # last receiver x in meters (ignored if only one receiver)
+zfin = 5000. # last receiver z in meters (ignored if only one receiver)
+enreg_surf = .true. # enregistrement volume ou surface
+# second receiver line
+nrec = 1 # number of receivers (can be 1 or more)
+xdeb = 6000. # first receiver x in meters
+zdeb = 7000. # first receiver z in meters
+xfin = 99999. # last receiver x in meters (ignored if only one receiver)
+zfin = 99999. # last receiver z in meters (ignored if only one receiver)
+enreg_surf = .false. # enregistrement volume ou surface
+
# display parameters
itaff = 250 # display frequency in time steps
output_postscript_snapshot = .true. # output Postscript image of the results
Modified: seismo/2D/SPECFEM2D/trunk/DATA/Par_file_cours_M2_UPPA
===================================================================
--- seismo/2D/SPECFEM2D/trunk/DATA/Par_file_cours_M2_UPPA 2005-02-23 17:33:20 UTC (rev 8471)
+++ seismo/2D/SPECFEM2D/trunk/DATA/Par_file_cours_M2_UPPA 2007-12-07 23:48:54 UTC (rev 8472)
@@ -1,16 +1,12 @@
#
# 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 = Test pour cours M2 UPPA
interfacesfile = interfaces_cours_M2_UPPA_curved.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 = 4000.d0 # abscissa of right side of the model
nx = 80 # number of elements along X
@@ -20,25 +16,22 @@
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 = 1600 # nb total de pas de temps
-dt = 1.d-3 # valeur du pas de temps
-#
+deltat = 1.d-3 # valeur du pas de temps
+
# source parameters
-#
source_surf = .true. # source dans le volume ou a la surface
xs = 2000. # source location x in meters
zs = 1500. # source location z in meters
-source_type = 1 # force = 1 or explosion = 2
+source_type = 1 # force = 1 or moment tensor = 2
time_function_type = 1 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4, Heaviside = 5
f0 = 10.0 # dominant source frequency (Hz) if not Dirac
angleforce = 0. # angle of the source (for a force only)
@@ -46,20 +39,21 @@
Mzz = 1. # Mzz component (for a moment tensor source only)
Mxz = 0. # Mxz component (for a moment tensor source only)
factor = 1.d10 # amplification factor
-#
+
# receiver line parameters
-#
-enreg_surf = .true. # enregistrement volume ou surface
sismostype = 1 # record 1=displ 2=veloc 3=accel 4=pressure
+nreceiverlines = 1 # number of receiver lines
+anglerec = 0.d0 # angle to rotate components at receivers
+
+# first receiver line
nrec = 11 # number of receivers
xdeb = 300. # first receiver x in meters
zdeb = 2200. # first receiver z in meters
-xfin = 3700. # last receiver x in meters
-zfin = 2200. # last receiver z in meters
-anglerec = 0.d0 # angle to rotate components at receivers
-#
+xfin = 3700. # last receiver x in meters (ignored if onlyone receiver)
+zfin = 2200. # last receiver z in meters (ignored if onlyone receiver)
+enreg_surf = .true. # enregistrement volume ou surface
+
# 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
@@ -74,9 +68,8 @@
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 = 3 # nb de modeles differents (0,rho,vp,vs,0,0)
1 0 2700.d0 3000.d0 1732.051d0 0 0
2 0 2500.d0 2700.d0 1558.845d0 0 0
Modified: seismo/2D/SPECFEM2D/trunk/meshfem2D.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/meshfem2D.f90 2005-02-23 17:33:20 UTC (rev 8471)
+++ seismo/2D/SPECFEM2D/trunk/meshfem2D.f90 2007-12-07 23:48:54 UTC (rev 8472)
@@ -44,13 +44,10 @@
xinterface_bottom,zinterface_bottom,coefs_interface_bottom, &
xinterface_top,zinterface_top,coefs_interface_top
-! for the source
- integer source_type,time_function_type
- double precision xs,zs,f0,t0,angleforce,Mxx,Mzz,Mxz,factor
+! for the source and receivers
+ integer source_type,time_function_type,nrec_total,irec_global_number
+ double precision xs,zs,f0,t0,angleforce,Mxx,Mzz,Mxz,factor,xrec,zrec
-! arrays for the receivers
- double precision, dimension(:), allocatable :: xrec,zrec
-
character(len=50) interfacesfile,title
integer imatnum,inumabs,inumsurface,inumelem
@@ -58,23 +55,26 @@
integer k,icol,ili,istepx,istepz,ix,iz,irec,i,j
integer ixdebzone,ixfinzone,izdebzone,izfinzone,imodnum
integer izone,imodele,nbzone,nbmodeles
- integer itaff,pointsdisp,subsamp,nrec
- integer sismostype,vecttype
- integer ngnod,nt,nx,nz,nxread,nzread
- integer icodematread
+ integer itaff,pointsdisp,subsamp,sismostype,vecttype
+ integer ngnod,nt,nx,nz,nxread,nzread,icodematread,ireceiverlines,nreceiverlines
+ integer, dimension(:), allocatable :: nrec
+
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,sizemax_arrows
- double precision anglerec,xfin,zfin,xdeb,zdeb,xmin,xmax,dt
+ double precision cutvect,sizemax_arrows,anglerec,xmin,xmax,deltat
double precision rhoread,cpread,csread,aniso3read,aniso4read
+ double precision, dimension(:), allocatable :: xdeb,zdeb,xfin,zfin
+
logical interpol,gnuplot,read_external_model,outputgrid
logical abshaut,absbas,absgauche,absdroite
- logical source_surf,enreg_surf,meshvect,initialfield,modelvect,boundvect
+ logical source_surf,meshvect,initialfield,modelvect,boundvect
logical ELASTIC,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON
+ logical, dimension(:), allocatable :: enreg_surf
+
integer, external :: num
double precision, external :: value_spline
@@ -185,7 +185,7 @@
! read time step parameters
call read_value_integer(IIN,IGNORE_JUNK,nt)
- call read_value_double_precision(IIN,IGNORE_JUNK,dt)
+ call read_value_double_precision(IIN,IGNORE_JUNK,deltat)
! read source parameters
call read_value_logical(IIN,IGNORE_JUNK,source_surf)
@@ -202,7 +202,7 @@
! if Dirac source time function, use a very thin Gaussian instead
! if Heaviside source time function, use a very thin error function instead
- if(time_function_type == 4 .or. time_function_type == 5) f0 = 1.d0 / (10.d0 * dt)
+ if(time_function_type == 4 .or. time_function_type == 5) f0 = 1.d0 / (10.d0 * deltat)
! time delay of the source in seconds, use a 20 % security margin (use 2 / f0 if error function)
if(time_function_type == 5) then
@@ -216,33 +216,36 @@
print *,'Position xs, zs = ',xs,zs
print *,'Frequency, delay = ',f0,t0
print *,'Source type (1=force, 2=explosion): ',source_type
- print *,'Time function type (1=Ricker, 2=First derivative, 3=Gaussian, 4=Dirac): ',time_function_type
+ print *,'Time function type (1=Ricker, 2=First derivative, 3=Gaussian, 4=Dirac, 5=Heaviside): ',time_function_type
print *,'Angle of the source if force = ',angleforce
print *,'Mxx of the source if moment tensor = ',Mxx
print *,'Mzz of the source if moment tensor = ',Mzz
print *,'Mxz of the source if moment tensor = ',Mxz
print *,'Multiplying factor = ',factor
-! read receivers line parameters
- call read_value_logical(IIN,IGNORE_JUNK,enreg_surf)
+! read receiver line parameters
call read_value_integer(IIN,IGNORE_JUNK,sismostype)
- call read_value_integer(IIN,IGNORE_JUNK,nrec)
- call read_value_double_precision(IIN,IGNORE_JUNK,xdeb)
- call read_value_double_precision(IIN,IGNORE_JUNK,zdeb)
- call read_value_double_precision(IIN,IGNORE_JUNK,xfin)
- call read_value_double_precision(IIN,IGNORE_JUNK,zfin)
+ call read_value_integer(IIN,IGNORE_JUNK,nreceiverlines)
call read_value_double_precision(IIN,IGNORE_JUNK,anglerec)
- allocate(xrec(nrec))
- allocate(zrec(nrec))
+ if(nreceiverlines < 1) stop 'number of receiver lines must be greater than 1'
- print *
- print *,'There are ',nrec,' receivers'
- xspacerec = (xfin-xdeb) / dble(nrec-1)
- zspacerec = (zfin-zdeb) / dble(nrec-1)
- do i=1,nrec
- xrec(i) = xdeb + dble(i-1)*xspacerec
- zrec(i) = zdeb + dble(i-1)*zspacerec
+! allocate receiver line arrays
+ allocate(nrec(nreceiverlines))
+ allocate(xdeb(nreceiverlines))
+ allocate(zdeb(nreceiverlines))
+ allocate(xfin(nreceiverlines))
+ allocate(zfin(nreceiverlines))
+ allocate(enreg_surf(nreceiverlines))
+
+! loop on all the receiver lines
+ do ireceiverlines = 1,nreceiverlines
+ call read_value_integer(IIN,IGNORE_JUNK,nrec(ireceiverlines))
+ call read_value_double_precision(IIN,IGNORE_JUNK,xdeb(ireceiverlines))
+ call read_value_double_precision(IIN,IGNORE_JUNK,zdeb(ireceiverlines))
+ call read_value_double_precision(IIN,IGNORE_JUNK,xfin(ireceiverlines))
+ call read_value_double_precision(IIN,IGNORE_JUNK,zfin(ireceiverlines))
+ call read_value_logical(IIN,IGNORE_JUNK,enreg_surf(ireceiverlines))
enddo
! read display parameters
@@ -470,20 +473,10 @@
call spline(xinterface_top,zinterface_top,npoints_interface_top,tang1,tangN,coefs_interface_top)
! tester si on est sur la derniere couche, qui contient la topographie
- if(ilayer == number_of_layers) then
+! et modifier position de la source si source exactement en surface
+ if(source_surf .and. ilayer == number_of_layers) &
+ zs = value_spline(xs,xinterface_top,zinterface_top,coefs_interface_top,npoints_interface_top)
-! modifier position de la source si source exactement en surface
- if(source_surf) zs = value_spline(xs,xinterface_top,zinterface_top,coefs_interface_top,npoints_interface_top)
-
-! modifier position des recepteurs si enregistrement exactement en surface
- if(enreg_surf) then
- do irec=1,nrec
- zrec(irec) = value_spline(xrec(irec),xinterface_top,zinterface_top,coefs_interface_top,npoints_interface_top)
- enddo
- endif
-
- endif
-
! calcul de l'offset de cette couche en nombre d'elements spectraux suivant Z
if(ilayer > 1) then
ioffset = sum(nz_layer(1:ilayer-1))
@@ -532,18 +525,6 @@
print *,'Valeurs min et max de Z sur le maillage = ',minval(z),maxval(z)
print *
-! afficher position de la source et des recepteurs
- print *
- print *,'Position (x,z) de la source'
- print *
- print *,'Source = ',xs,zs
- print *
- print *,'Position (x,z) des ',nrec,' recepteurs'
- print *
- do irec=1,nrec
- print *,'Receiver ',irec,' = ',xrec(irec),zrec(irec)
- enddo
-
! ***
! *** generer un fichier Gnuplot pour le controle de la grille ***
! ***
@@ -645,16 +626,16 @@
write(15,*) 'read_external_model outputgrid ELASTIC TURN_ANISOTROPY_ON TURN_ATTENUATION_ON'
write(15,*) read_external_model,outputgrid,ELASTIC,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON
- write(15,*) 'ncycl dtinc'
- write(15,*) nt,dt
+ write(15,*) 'nt deltat'
+ write(15,*) nt,deltat
write(15,*) 'source'
- write(15,*) source_type,time_function_type,xs-xmin,zs,f0,t0,factor,angleforce,Mxx,Mzz,Mxz
+ write(15,*) source_type,time_function_type,xs,zs,f0,t0,factor,angleforce,Mxx,Mzz,Mxz
write(15,*) 'Coordinates of macrobloc mesh (coorg):'
do j=0,nz
do i=0,nx
- write(15,*) num(i,j,nx),x(i,j)-xmin,z(i,j)
+ write(15,*) num(i,j,nx),x(i,j),z(i,j)
enddo
enddo
@@ -781,25 +762,70 @@
close(15)
+! print position of the source
+ print *
+ print *,'Position (x,z) of the source = ',xs,zs
+ print *
-!
-!--- write the STATIONS file
-!
+!--- compute position of the receivers and write the STATIONS file
+
print *
print *,'writing the DATA/STATIONS file'
print *
+! total number of receivers in all the receiver lines
+ nrec_total = sum(nrec)
+
+ print *
+ print *,'There are ',nrec_total,' receivers'
+
+ print *
+ print *,'Position (x,z) of the ',nrec_total,' receivers'
+ print *
+
open(unit=15,file='DATA/STATIONS',status='unknown')
- write(15,*) nrec
- do irec=1,nrec
- write(15,100) irec,xrec(irec)-xmin,zrec(irec)
+ write(15,*) nrec_total
+
+ irec_global_number = 0
+
+! loop on all the receiver lines
+ do ireceiverlines = 1,nreceiverlines
+
+! loop on all the receivers of this receiver line
+ do irec = 1,nrec(ireceiverlines)
+
+! compute global receiver number
+ irec_global_number = irec_global_number + 1
+
+! compute coordinates of the receiver
+ if(nrec(ireceiverlines) > 1) then
+ xrec = xdeb(ireceiverlines) + dble(irec-1)*(xfin(ireceiverlines)-xdeb(ireceiverlines))/dble(nrec(ireceiverlines)-1)
+ zrec = zdeb(ireceiverlines) + dble(irec-1)*(zfin(ireceiverlines)-zdeb(ireceiverlines))/dble(nrec(ireceiverlines)-1)
+ else
+ xrec = xdeb(ireceiverlines)
+ zrec = zdeb(ireceiverlines)
+ endif
+
+! modifier position du recepteur si enregistrement exactement en surface
+ if(enreg_surf(ireceiverlines)) &
+ zrec = value_spline(xrec,xinterface_top,zinterface_top,coefs_interface_top,npoints_interface_top)
+
+! display position of the receiver
+ print *,'Receiver ',irec_global_number,' = ',xrec,zrec
+
+ write(15,100) irec_global_number,xrec,zrec
+
+ enddo
enddo
+
close(15)
+ print *
+
10 format('')
40 format(a50)
- 100 format('S',i3.3,' AA ',f20.7,1x,f20.7,' 0.0 0.0')
+ 100 format('S',i4.4,' AA ',f20.7,1x,f20.7,' 0.0 0.0')
end program meshfem2D
Modified: seismo/2D/SPECFEM2D/trunk/read_value_parameters.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/read_value_parameters.f90 2005-02-23 17:33:20 UTC (rev 8471)
+++ seismo/2D/SPECFEM2D/trunk/read_value_parameters.f90 2007-12-07 23:48:54 UTC (rev 8472)
@@ -84,7 +84,7 @@
logical ignore_junk
character(len=100) string_read
- integer ios,iin
+ integer ios,iin,position_equal_sign
do
read(unit=iin,fmt=200,iostat=ios) string_read
@@ -105,9 +105,16 @@
! suppress trailing comments, if any
if(index(string_read,'#') > 0) string_read = string_read(1:index(string_read,'#')-1)
-! suppress leading junk (first 34 characters) if needed
- if(ignore_junk) string_read = string_read(35:len_trim(string_read))
+! suppress leading junk (until the first equal sign) if needed
+ if(ignore_junk) then
+ position_equal_sign = index(string_read,'=')
+ if(position_equal_sign <= 1) stop 'incorrect syntax detected in DATA/Par_file'
+ string_read = string_read(position_equal_sign+1:len_trim(string_read))
+ endif
+! suppress leading white spaces again, if any, after having suppressed the leading junk
+ string_read = adjustl(string_read)
+
! format
200 format(a100)
Modified: seismo/2D/SPECFEM2D/trunk/specfem2D.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/specfem2D.f90 2005-02-23 17:33:20 UTC (rev 8471)
+++ seismo/2D/SPECFEM2D/trunk/specfem2D.f90 2007-12-07 23:48:54 UTC (rev 8472)
@@ -177,6 +177,7 @@
! for color PNM images
integer :: NX_IMAGE_PNM,NZ_IMAGE_PNM,iplus1,jplus1,iminus1,jminus1,nx_sem_PNM
double precision :: xmin_PNM_image,xmax_PNM_image,zmin_PNM_image,zmax_PNM_image,taille_pixel_horizontal,taille_pixel_vertical
+ integer, dimension(:), allocatable :: ispec_for_PNM_image
integer, dimension(:,:), allocatable :: iglob_image_PNM_2D,copy_iglob_image_PNM_2D
double precision, dimension(:,:), allocatable :: donnees_image_PNM_2D
@@ -800,6 +801,17 @@
write(IOUT,*) 'fin localisation de tous les pixels des images PNM'
+! assign ispec number to be able to determine density for acoustic PNM snapshots
+ allocate(ispec_for_PNM_image(npoin))
+ do ispec = 1,nspec
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+ iglob = ibool(i,j,ispec)
+ ispec_for_PNM_image(iglob) = ispec
+ enddo
+ enddo
+ enddo
+
!
!---- initialiser sismogrammes
!
@@ -1736,20 +1748,34 @@
do j = 1,NZ_IMAGE_PNM
do i = 1,NX_IMAGE_PNM
- if(iglob_image_PNM_2D(i,j) /= -1) then
+
+ iglob = iglob_image_PNM_2D(i,j)
+
+ if(iglob /= -1) then
! 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))
+ donnees_image_PNM_2D(i,j) = displ(2,iglob)
else if(vecttype == 2) then
- donnees_image_PNM_2D(i,j) = veloc(2,iglob_image_PNM_2D(i,j))
+ donnees_image_PNM_2D(i,j) = veloc(2,iglob)
else
- donnees_image_PNM_2D(i,j) = accel(2,iglob_image_PNM_2D(i,j))
+ donnees_image_PNM_2D(i,j) = accel(2,iglob)
endif
+
else
! 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'
+
+! pressure = - rho * Chi_dot
+ material = kmato(ispec_for_PNM_image(iglob))
+ denst = density(material)
+ if(read_external_model) denst = rhoext(iglob)
+
+ donnees_image_PNM_2D(i,j) = - denst * veloc(1,iglob)
+
+! uncomment this for vertical component of velocity vector instead in acoustic medium
+! donnees_image_PNM_2D(i,j) = vector_field_postscript(2,iglob)
+
endif
endif
enddo
Modified: seismo/2D/SPECFEM2D/trunk/su_create_Paul.python
===================================================================
--- seismo/2D/SPECFEM2D/trunk/su_create_Paul.python 2005-02-23 17:33:20 UTC (rev 8471)
+++ seismo/2D/SPECFEM2D/trunk/su_create_Paul.python 2007-12-07 23:48:54 UTC (rev 8472)
@@ -44,7 +44,7 @@
print sismostype
chdir(SEM+'/OUTPUT_FILES')
labels='label1="Time" label2="Receivers"'
- # Create headers ans Su file
+ # Create headers and Su file
if sismostype==4:
ordres=['suaddhead < pressure_file_single.bin ns='+str(nt)+' | sushw key=dt a='+str(int(dt*1e6))+' > pressure_file.su']
ordres.append('suxwigb < pressure_file.su perc=96 '+labels+' title=" Pressure : '+title+'"&')
@@ -57,4 +57,4 @@
for i in range(len(ordres)):
system(ordres[i])
if __name__=='__main__':
- createSU()
\ No newline at end of file
+ createSU()
More information about the cig-commits
mailing list