[cig-commits] r8438 - in seismo/2D/SPECFEM2D/trunk: MAILLE90
SPECFEM90
walter at geodynamics.org
walter at geodynamics.org
Fri Dec 7 15:46:22 PST 2007
Author: walter
Date: 2007-12-07 15:46:22 -0800 (Fri, 07 Dec 2007)
New Revision: 8438
Modified:
seismo/2D/SPECFEM2D/trunk/MAILLE90/Par
seismo/2D/SPECFEM2D/trunk/MAILLE90/Par_cours_M2_UPPA
seismo/2D/SPECFEM2D/trunk/MAILLE90/interfaces_cours_M2_UPPA_curved.dat
seismo/2D/SPECFEM2D/trunk/MAILLE90/meshfem2D.f90
seismo/2D/SPECFEM2D/trunk/SPECFEM90/Makefile
seismo/2D/SPECFEM2D/trunk/SPECFEM90/checkgrid.f90
seismo/2D/SPECFEM2D/trunk/SPECFEM90/defarrays.f90
seismo/2D/SPECFEM2D/trunk/SPECFEM90/plotpost.f90
seismo/2D/SPECFEM2D/trunk/SPECFEM90/positrec.f90
seismo/2D/SPECFEM2D/trunk/SPECFEM90/positsource.f90
seismo/2D/SPECFEM2D/trunk/SPECFEM90/specfem2D.f90
seismo/2D/SPECFEM2D/trunk/SPECFEM90/write_seismograms.f90
Log:
Dirac and Gaussian time sources and corresponding convolution routine,
more flexible Par file with any number of comment lines,
Xsu scripts generated automatically,
also cleaned obsolete stuff in 2D solver
Modified: seismo/2D/SPECFEM2D/trunk/MAILLE90/Par
===================================================================
--- seismo/2D/SPECFEM2D/trunk/MAILLE90/Par 2004-12-18 20:09:40 UTC (rev 8437)
+++ seismo/2D/SPECFEM2D/trunk/MAILLE90/Par 2007-12-07 23:46:22 UTC (rev 8438)
@@ -4,7 +4,6 @@
# Put variable names first and actual value after 34th column
#
# ----------------------------------------------------------------
-<- ->
#
# title of job, and file that contains interface data
#
@@ -18,7 +17,7 @@
nx = 80 ! number of elements along X
ngnod = 4 ! noeuds de controle pour blocs (4 ou 9)
initialfield = .false. ! use a plane wave as source or not
-ireadmodel = .false. ! read external earth model or not
+readmodel = .false. ! read external earth model or not
TURN_ANISOTROPY_ON = .false. ! turn anisotropy on or off
TURN_ATTENUATION_ON = .false. ! turn attenuation on or off
#
@@ -36,19 +35,19 @@
#
# source parameters
#
-isource_surf = .true. ! sources dans le volume ou a la surface
+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
-f0 = 10.0 ! central source frequency (Hz)
-t0 = 0.11 ! time delay of the source in seconds
-source_type = 1 ! source type : force=1 or explosion=2
+source_type = 1 ! force = 1 or explosion = 2
+time_function_type = 1 ! Ricker = 1, Gaussian = 2, Dirac = 3
+f0 = 10.0 ! dominant source frequency (Hz) if not Dirac
angle = 0. ! angle of the source (for a force only)
factor = 1.d10 ! amplification factor
#
# receiver line parameters
#
-ienreg_surf = .true. ! enregistrement volume ou surface
-isismostype = 1 ! record 1=displ 2=veloc 3=accel
+enreg_surf = .true. ! enregistrement volume ou surface
+sismostype = 1 ! record 1=displ 2=veloc 3=accel
nrec = 11 ! number of receivers
xdeb = 300. ! first receiver x in meters
zdeb = 2200. ! first receiver z in meters
@@ -59,16 +58,16 @@
# display parameters
#
itaff = 100 ! display frequency in time steps
-ivecttype = 1 ! display 1=displ 2=veloc 3=accel
+vecttype = 1 ! display 1=displ 2=veloc 3=accel
cutvect = 1. ! amplitude min en % pour vector plots
-imeshvect = .true. ! display mesh on vector plots or not
-imodelvect = .false. ! display velocity model on vector plots
-iboundvect = .true. ! display boundary conditions on 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
-iptsdisp = 6 ! points for interpolation of display
-isubsamp = 1 ! subsampling of color snapshots
-ignuplot = .false. ! generate a GNUPLOT file for the grid
-ioutputgrid = .false. ! save the grid in a text file or not
+pointsdisp = 6 ! points for interpolation of display
+subsamp = 1 ! subsampling of color snapshots
+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)
#
Modified: seismo/2D/SPECFEM2D/trunk/MAILLE90/Par_cours_M2_UPPA
===================================================================
--- seismo/2D/SPECFEM2D/trunk/MAILLE90/Par_cours_M2_UPPA 2004-12-18 20:09:40 UTC (rev 8437)
+++ seismo/2D/SPECFEM2D/trunk/MAILLE90/Par_cours_M2_UPPA 2007-12-07 23:46:22 UTC (rev 8438)
@@ -4,7 +4,6 @@
# Put variable names first and actual value after 34th column
#
# ----------------------------------------------------------------
-<- ->
#
# title of job, and file that contains interface data
#
@@ -18,7 +17,7 @@
nx = 80 ! number of elements along X
ngnod = 4 ! noeuds de controle pour blocs (4 ou 9)
initialfield = .false. ! use a plane wave as source or not
-ireadmodel = .false. ! read external earth model or not
+readmodel = .false. ! read external earth model or not
TURN_ANISOTROPY_ON = .false. ! turn anisotropy on or off
TURN_ATTENUATION_ON = .false. ! turn attenuation on or off
#
@@ -36,19 +35,19 @@
#
# source parameters
#
-isource_surf = .true. ! sources dans le volume ou a la surface
+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
-f0 = 10.0 ! central source frequency (Hz)
-t0 = 0.11 ! time delay of the source in seconds
-source_type = 1 ! source type : force=1 or explosion=2
+source_type = 1 ! force = 1 or explosion = 2
+time_function_type = 1 ! Ricker = 1, Gaussian = 2, Dirac = 3
+f0 = 10.0 ! dominant source frequency (Hz) if not Dirac
angle = 0. ! angle of the source (for a force only)
factor = 1.d10 ! amplification factor
#
# receiver line parameters
#
-ienreg_surf = .true. ! enregistrement volume ou surface
-isismostype = 1 ! record 1=displ 2=veloc 3=accel
+enreg_surf = .true. ! enregistrement volume ou surface
+sismostype = 1 ! record 1=displ 2=veloc 3=accel
nrec = 11 ! number of receivers
xdeb = 300. ! first receiver x in meters
zdeb = 2200. ! first receiver z in meters
@@ -59,16 +58,16 @@
# display parameters
#
itaff = 100 ! display frequency in time steps
-ivecttype = 1 ! display 1=displ 2=veloc 3=accel
+vecttype = 1 ! display 1=displ 2=veloc 3=accel
cutvect = 1. ! amplitude min en % pour vector plots
-imeshvect = .true. ! display mesh on vector plots or not
-imodelvect = .false. ! display velocity model on vector plots
-iboundvect = .true. ! display boundary conditions on 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
-iptsdisp = 6 ! points for interpolation of display
-isubsamp = 1 ! subsampling of color snapshots
-ignuplot = .false. ! generate a GNUPLOT file for the grid
-ioutputgrid = .false. ! save the grid in a text file or not
+pointsdisp = 6 ! points for interpolation of display
+subsamp = 1 ! subsampling of color snapshots
+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)
#
Modified: seismo/2D/SPECFEM2D/trunk/MAILLE90/interfaces_cours_M2_UPPA_curved.dat
===================================================================
--- seismo/2D/SPECFEM2D/trunk/MAILLE90/interfaces_cours_M2_UPPA_curved.dat 2004-12-18 20:09:40 UTC (rev 8437)
+++ seismo/2D/SPECFEM2D/trunk/MAILLE90/interfaces_cours_M2_UPPA_curved.dat 2007-12-07 23:46:22 UTC (rev 8438)
@@ -1,19 +1,13 @@
-#
# number of interfaces
-#
4
#
# 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
5000 0
-#
# interface number 2
-#
7
0 1000
1500 1100
@@ -22,9 +16,7 @@
3000 1220
3500 1170
5000 1100
-#
# interface number 3
-#
9
0 2000
500 2000
@@ -35,9 +27,7 @@
3000 2090
3500 2020
5000 2000
-#
# interface number 4 (topography, top of the mesh)
-#
8
0 3000
500 3000
@@ -50,15 +40,9 @@
#
# for each layer, we give the number of spectral elements in the vertical direction
#
-#
# layer number 1 (bottom layer)
-#
20
-#
# layer number 2
-#
20
-#
# layer number 3 (top layer)
-#
20
Modified: seismo/2D/SPECFEM2D/trunk/MAILLE90/meshfem2D.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/MAILLE90/meshfem2D.f90 2004-12-18 20:09:40 UTC (rev 8437)
+++ seismo/2D/SPECFEM2D/trunk/MAILLE90/meshfem2D.f90 2007-12-07 23:46:22 UTC (rev 8438)
@@ -21,10 +21,6 @@
implicit none
-! definir les tableaux pour allocation dynamique
-
- integer, parameter :: ANISOTROPIC_MATERIAL = 1
-
! coordinates of the grid points of the mesh
double precision, dimension(:,:), allocatable :: x,z
@@ -47,22 +43,21 @@
xinterface_top,zinterface_top,coefs_interface_top
! for the source
- integer isource_type
+ integer source_type,time_function_type
double precision xs,zs,f0,t0,angle,factor
! arrays for the receivers
double precision, dimension(:), allocatable :: xrec,zrec
character(len=50) interfacesfile,title
- character(len=34) junk
integer imatnum,inumabs,inumelem
integer nelemabs,npgeo,nspec
integer k,icol,ili,istepx,istepz,ix,iz,irec,i,j
integer ixdebzone,ixfinzone,izdebzone,izfinzone,imodnum
integer izone,imodele,nbzone,nbmodeles
- integer itaff,iptsdisp,isubsamp,nrec
- integer isismostype,ivecttype
+ integer itaff,pointsdisp,subsamp,nrec
+ integer sismostype,vecttype
integer ngnod,nt,nx,nz,nxread,nzread
integer icodematread
@@ -73,14 +68,24 @@
double precision anglerec,xfin,zfin,xdeb,zdeb,xmin,xmax,dt
double precision rhoread,cpread,csread,aniso3read,aniso4read
- logical interpol,ignuplot,ireadmodel,ioutputgrid
+ logical interpol,gnuplot,readmodel,outputgrid
logical abshaut,absbas,absgauche,absdroite
- logical isource_surf,ienreg_surf,imeshvect,initialfield,imodelvect,iboundvect
+ logical source_surf,enreg_surf,meshvect,initialfield,modelvect,boundvect
logical TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON
integer, external :: num
double precision, external :: value_spline
+! flag to indicate an anisotropic material
+ integer, parameter :: ANISOTROPIC_MATERIAL = 1
+
+! file number for input Par file and interface file
+ integer, parameter :: IIN_PAR = 10
+ integer, parameter :: IIN_INTERFACES = 15
+
+! ignore variable name field (junk) at the beginning of each input line
+ logical, parameter :: IGNORE_JUNK = .true.,DONT_IGNORE_JUNK = .false.
+
! ***
! *** read the parameter file
! ***
@@ -90,39 +95,23 @@
open(unit=10,file='Par',status='old')
-! formats
- 1 format(a,f12.5)
- 2 format(a,i8)
- 3 format(a,a)
- 4 format(a,l8)
-
-! read the header
- do i=1,10
- read(10,*)
- enddo
-
! read file names and path for output
- read(10,3) junk,title
- read(10,3) junk,interfacesfile
+ call read_value_string(IIN_PAR,IGNORE_JUNK,title)
+ call read_value_string(IIN_PAR,IGNORE_JUNK,interfacesfile)
write(*,*) 'Titre de la simulation'
write(*,*) title
print *
-! skip comment
- read(10,*)
- read(10,*)
- read(10,*)
-
! read grid parameters
- read(10,1) junk,xmin
- read(10,1) junk,xmax
- read(10,2) junk,nx
- read(10,2) junk,ngnod
- read(10,4) junk,initialfield
- read(10,4) junk,ireadmodel
- read(10,4) junk,TURN_ANISOTROPY_ON
- read(10,4) junk,TURN_ATTENUATION_ON
+ call read_value_double_precision(IIN_PAR,IGNORE_JUNK,xmin)
+ call read_value_double_precision(IIN_PAR,IGNORE_JUNK,xmax)
+ call read_value_integer(IIN_PAR,IGNORE_JUNK,nx)
+ call read_value_integer(IIN_PAR,IGNORE_JUNK,ngnod)
+ call read_value_logical(IIN_PAR,IGNORE_JUNK,initialfield)
+ call read_value_logical(IIN_PAR,IGNORE_JUNK,readmodel)
+ call read_value_logical(IIN_PAR,IGNORE_JUNK,TURN_ANISOTROPY_ON)
+ call read_value_logical(IIN_PAR,IGNORE_JUNK,TURN_ATTENUATION_ON)
! get interface data from external file to count the spectral elements along Z
print *,'Reading interface data from file ',interfacesfile(1:len_trim(interfacesfile)),' to count the spectral elements'
@@ -130,36 +119,21 @@
max_npoints_interface = -1
-! skip header
- read(15,*)
- read(15,*)
- read(15,*)
-
! read number of interfaces
- read(15,*) number_of_interfaces
+ call read_value_integer(IIN_INTERFACES,DONT_IGNORE_JUNK,number_of_interfaces)
if(number_of_interfaces < 2) stop 'not enough interfaces (minimum is 2)'
-! skip header
- read(15,*)
- read(15,*)
- read(15,*)
-
! loop on all the interfaces
do interface_current = 1,number_of_interfaces
-! skip header
- read(15,*)
- read(15,*)
- read(15,*)
-
- read(15,*) npoints_interface_bottom
+ call read_value_integer(IIN_INTERFACES,DONT_IGNORE_JUNK,npoints_interface_bottom)
if(npoints_interface_bottom < 2) stop 'not enough interface points (minimum is 2)'
max_npoints_interface = max(npoints_interface_bottom,max_npoints_interface)
print *,'Reading ',npoints_interface_bottom,' points for interface ',interface_current
! loop on all the points describing this interface
do ipoint_current = 1,npoints_interface_bottom
- read(15,*) xinterface_dummy,zinterface_dummy
+ read(IIN_INTERFACES,*) xinterface_dummy,zinterface_dummy
if(ipoint_current > 1 .and. xinterface_dummy <= xinterface_dummy_previous) &
stop 'interface points must be sorted in increasing X'
xinterface_dummy_previous = xinterface_dummy
@@ -172,21 +146,11 @@
allocate(nz_layer(number_of_layers))
-! skip header
- read(15,*)
- read(15,*)
- read(15,*)
-
! loop on all the layers
do ilayer = 1,number_of_layers
-! skip header
- read(15,*)
- read(15,*)
- read(15,*)
-
! read number of spectral elements in vertical direction in this layer
- read(15,*) nz_layer(ilayer)
+ call read_value_integer(IIN_INTERFACES,DONT_IGNORE_JUNK,nz_layer(ilayer))
if(nz_layer(ilayer) < 1) stop 'not enough spectral elements along Z in layer (minimum is 1)'
print *,'There are ',nz_layer(ilayer),' spectral elements along Z in layer ',ilayer
@@ -210,63 +174,50 @@
nz = nz * 2
endif
-! skip comment
- read(10,*)
- read(10,*)
- read(10,*)
-
! read absorbing boundaries parameters
- read(10,4) junk,abshaut
- read(10,4) junk,absbas
- read(10,4) junk,absgauche
- read(10,4) junk,absdroite
+ call read_value_logical(IIN_PAR,IGNORE_JUNK,abshaut)
+ call read_value_logical(IIN_PAR,IGNORE_JUNK,absbas)
+ call read_value_logical(IIN_PAR,IGNORE_JUNK,absgauche)
+ call read_value_logical(IIN_PAR,IGNORE_JUNK,absdroite)
-! skip comment
- read(10,*)
- read(10,*)
- read(10,*)
-
! read time step parameters
- read(10,2) junk,nt
- read(10,1) junk,dt
+ call read_value_integer(IIN_PAR,IGNORE_JUNK,nt)
+ call read_value_double_precision(IIN_PAR,IGNORE_JUNK,dt)
-! skip comment
- read(10,*)
- read(10,*)
- read(10,*)
-
! read source parameters
- read(10,4) junk,isource_surf
- read(10,1) junk,xs
- read(10,1) junk,zs
- read(10,1) junk,f0
- read(10,1) junk,t0
- read(10,2) junk,isource_type
- read(10,1) junk,angle
- read(10,1) junk,factor
+ call read_value_logical(IIN_PAR,IGNORE_JUNK,source_surf)
+ call read_value_double_precision(IIN_PAR,IGNORE_JUNK,xs)
+ call read_value_double_precision(IIN_PAR,IGNORE_JUNK,zs)
+ call read_value_integer(IIN_PAR,IGNORE_JUNK,source_type)
+ call read_value_integer(IIN_PAR,IGNORE_JUNK,time_function_type)
+ call read_value_double_precision(IIN_PAR,IGNORE_JUNK,f0)
+ call read_value_double_precision(IIN_PAR,IGNORE_JUNK,angle)
+ call read_value_double_precision(IIN_PAR,IGNORE_JUNK,factor)
+! if Dirac source time function, use a very thin Gaussian instead
+ if(time_function_type == 3) f0 = 1.d0 / (5.d0 * dt)
+
+! time delay of the source in seconds, use a 20 % security margin
+ t0 = 1.20d0 / f0
+
print *
print *,'Source:'
print *,'Position xs, zs = ',xs,zs
print *,'Frequency, delay = ',f0,t0
- print *,'Source type (1=force 2=explosion) : ',isource_type
+ print *,'Source type (1=force, 2=explosion): ',source_type
+ print *,'Time function type (1=Ricker, 2=Gaussian, 3=Dirac): ',time_function_type
print *,'Angle of the source if force = ',angle
print *,'Multiplying factor = ',factor
-! skip comment
- read(10,*)
- read(10,*)
- read(10,*)
-
! read receivers line parameters
- read(10,4) junk,ienreg_surf
- read(10,2) junk,isismostype
- read(10,2) junk,nrec
- read(10,1) junk,xdeb
- read(10,1) junk,zdeb
- read(10,1) junk,xfin
- read(10,1) junk,zfin
- read(10,1) junk,anglerec
+ call read_value_logical(IIN_PAR,IGNORE_JUNK,enreg_surf)
+ call read_value_integer(IIN_PAR,IGNORE_JUNK,sismostype)
+ call read_value_integer(IIN_PAR,IGNORE_JUNK,nrec)
+ call read_value_double_precision(IIN_PAR,IGNORE_JUNK,xdeb)
+ call read_value_double_precision(IIN_PAR,IGNORE_JUNK,zdeb)
+ call read_value_double_precision(IIN_PAR,IGNORE_JUNK,xfin)
+ call read_value_double_precision(IIN_PAR,IGNORE_JUNK,zfin)
+ call read_value_double_precision(IIN_PAR,IGNORE_JUNK,anglerec)
allocate(xrec(nrec))
allocate(zrec(nrec))
@@ -280,31 +231,21 @@
zrec(i) = zdeb + dble(i-1)*zspacerec
enddo
-! skip comment
- read(10,*)
- read(10,*)
- read(10,*)
-
! read display parameters
- read(10,2) junk,itaff
- read(10,2) junk,ivecttype
- read(10,1) junk,cutvect
- read(10,4) junk,imeshvect
- read(10,4) junk,imodelvect
- read(10,4) junk,iboundvect
- read(10,4) junk,interpol
- read(10,2) junk,iptsdisp
- read(10,2) junk,isubsamp
- read(10,4) junk,ignuplot
- read(10,4) junk,ioutputgrid
+ call read_value_integer(IIN_PAR,IGNORE_JUNK,itaff)
+ call read_value_integer(IIN_PAR,IGNORE_JUNK,vecttype)
+ call read_value_double_precision(IIN_PAR,IGNORE_JUNK,cutvect)
+ call read_value_logical(IIN_PAR,IGNORE_JUNK,meshvect)
+ call read_value_logical(IIN_PAR,IGNORE_JUNK,modelvect)
+ call read_value_logical(IIN_PAR,IGNORE_JUNK,boundvect)
+ call read_value_logical(IIN_PAR,IGNORE_JUNK,interpol)
+ call read_value_integer(IIN_PAR,IGNORE_JUNK,pointsdisp)
+ call read_value_integer(IIN_PAR,IGNORE_JUNK,subsamp)
+ call read_value_logical(IIN_PAR,IGNORE_JUNK,gnuplot)
+ call read_value_logical(IIN_PAR,IGNORE_JUNK,outputgrid)
-! skip comment
- read(10,*)
- read(10,*)
- read(10,*)
-
! lecture des differents modeles de materiaux
- read(10,2) junk,nbmodeles
+ call read_value_integer(IIN_PAR,IGNORE_JUNK,nbmodeles)
if(nbmodeles <= 0) stop 'Negative number of models not allowed !!'
allocate(icodemat(nbmodeles))
@@ -324,7 +265,7 @@
num_modele(:,:) = 0
do imodele=1,nbmodeles
- read(10,*) i,icodematread,rhoread,cpread,csread,aniso3read,aniso4read
+ read(IIN_PAR,*) i,icodematread,rhoread,cpread,csread,aniso3read,aniso4read
if(i < 1 .or. i > nbmodeles) stop 'Wrong material point number'
icodemat(i) = icodematread
rho(i) = rhoread
@@ -351,7 +292,7 @@
print *
! lecture des numeros de modele des differentes zones
- read(10,2) junk,nbzone
+ call read_value_integer(IIN_PAR,IGNORE_JUNK,nbzone)
if(nbzone <= 0) stop 'Negative number of zones not allowed !!'
@@ -361,7 +302,7 @@
do izone = 1,nbzone
- read(10,*) ixdebzone,ixfinzone,izdebzone,izfinzone,imodnum
+ read(IIN_PAR,*) ixdebzone,ixfinzone,izdebzone,izfinzone,imodnum
if(imodnum < 1) stop 'Negative model number not allowed !!'
if(ixdebzone < 1) stop 'Left coordinate of zone negative !!'
@@ -450,46 +391,26 @@
allocate(zinterface_top(max_npoints_interface))
allocate(coefs_interface_top(max_npoints_interface))
-! skip header
- read(15,*)
- read(15,*)
- read(15,*)
-
! read number of interfaces
- read(15,*) number_of_interfaces
+ call read_value_integer(IIN_INTERFACES,DONT_IGNORE_JUNK,number_of_interfaces)
-! skip header
- read(15,*)
- read(15,*)
- read(15,*)
-
-! skip header
- read(15,*)
- read(15,*)
- read(15,*)
-
! read bottom interface
- read(15,*) npoints_interface_bottom
+ call read_value_integer(IIN_INTERFACES,DONT_IGNORE_JUNK,npoints_interface_bottom)
! loop on all the points describing this interface
do ipoint_current = 1,npoints_interface_bottom
- read(15,*) xinterface_bottom(ipoint_current),zinterface_bottom(ipoint_current)
+ read(IIN_INTERFACES,*) xinterface_bottom(ipoint_current),zinterface_bottom(ipoint_current)
enddo
! boucle sur toutes les couches
do ilayer = 1,number_of_layers
-! skip header
- read(15,*)
- read(15,*)
- read(15,*)
-
! read top interface
- read(15,*) npoints_interface_top
+ call read_value_integer(IIN_INTERFACES,DONT_IGNORE_JUNK,npoints_interface_top)
! loop on all the points describing this interface
do ipoint_current = 1,npoints_interface_top
- read(15,*) xinterface_top(ipoint_current),zinterface_top(ipoint_current)
+ read(IIN_INTERFACES,*) xinterface_top(ipoint_current),zinterface_top(ipoint_current)
enddo
! calculer le spline pour l'interface du bas, imposer la tangente aux deux bords
@@ -508,10 +429,10 @@
if(ilayer == number_of_layers) then
! modifier position de la source si source exactement en surface
- if(isource_surf) zs = value_spline(xs,xinterface_top,zinterface_top,coefs_interface_top,npoints_interface_top)
+ 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(ienreg_surf) then
+ 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
@@ -654,14 +575,14 @@
write(15,*) 'npgeo'
write(15,*) npgeo
- write(15,*) 'ignuplot interpol'
- write(15,*) ignuplot,interpol
+ write(15,*) 'gnuplot interpol'
+ write(15,*) gnuplot,interpol
- write(15,*) 'itaff icolor inumber'
- write(15,*) itaff,1,0
+ write(15,*) 'itaff colors numbers'
+ write(15,*) itaff,' 1 0'
- write(15,*) 'imeshvect imodelvect iboundvect cutvect isubsamp nx_sem_PNM'
- write(15,*) imeshvect,imodelvect,iboundvect,cutvect,isubsamp,nxread
+ write(15,*) 'meshvect modelvect boundvect cutvect subsamp nx_sem_PNM'
+ write(15,*) meshvect,modelvect,boundvect,cutvect,subsamp,nxread
write(15,*) 'nrec anglerec'
write(15,*) nrec,anglerec
@@ -669,17 +590,17 @@
write(15,*) 'initialfield'
write(15,*) initialfield
- write(15,*) 'isismostype ivecttype'
- write(15,*) isismostype,ivecttype
+ write(15,*) 'sismostype vecttype'
+ write(15,*) sismostype,vecttype
- write(15,*) 'ireadmodel ioutputgrid TURN_ANISOTROPY_ON TURN_ATTENUATION_ON'
- write(15,*) ireadmodel,ioutputgrid,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON
+ write(15,*) 'readmodel outputgrid TURN_ANISOTROPY_ON TURN_ATTENUATION_ON'
+ write(15,*) readmodel,outputgrid,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON
write(15,*) 'ncycl dtinc'
write(15,*) nt,dt
write(15,*) 'source'
- write(15,*) isource_type,xs-xmin,zs,f0,t0,factor,angle,' 0'
+ write(15,*) source_type,time_function_type,xs-xmin,zs,f0,t0,factor,angle
write(15,*) 'Receiver positions:'
do irec=1,nrec
@@ -712,8 +633,8 @@
if(abshaut .and. absgauche) nelemabs = nelemabs - 1
if(abshaut .and. absdroite) nelemabs = nelemabs - 1
- write(15,*) 'numat ngnod nspec iptsdisp ielemabs'
- write(15,*) nbmodeles,ngnod,nspec,iptsdisp,nelemabs
+ write(15,*) 'numat ngnod nspec pointsdisp nelemabs'
+ write(15,*) nbmodeles,ngnod,nspec,pointsdisp,nelemabs
write(15,*) 'Material sets (num 0 rho vp vs 0 0) or (num 1 rho c11 c13 c33 c44)'
do i=1,nbmodeles
@@ -909,3 +830,203 @@
end subroutine splint
+
+
+
+
+
+
+
+
+
+
+
+
+!--------------------
+
+ subroutine read_value_integer(iin,ignore_junk,value_to_read)
+
+ implicit none
+
+ integer iin
+ logical ignore_junk
+
+ integer value_to_read
+
+ integer ios
+
+ character(len=100) string_read
+ character(len=34) junk
+
+ do
+ read(unit=iin,fmt=200,iostat=ios) string_read
+ if(ios /= 0) stop 'error while reading input file'
+
+! suppress leading white spaces, if any
+ string_read = adjustl(string_read)
+
+! exit loop when we find the first line that is not a comment
+ if(string_read(1:1) /= '#') exit
+
+ enddo
+
+ if(ignore_junk) then
+ read(string_read,100) junk,value_to_read
+ else
+ read(string_read,*) value_to_read
+ endif
+
+! format
+ 100 format(a,i8)
+ 200 format(a100)
+
+ end subroutine read_value_integer
+
+
+
+
+
+
+
+
+
+
+!--------------------
+
+ subroutine read_value_double_precision(iin,ignore_junk,value_to_read)
+
+ implicit none
+
+ integer iin
+ logical ignore_junk
+
+ double precision value_to_read
+
+ integer ios
+
+ character(len=100) string_read
+ character(len=34) junk
+
+ do
+ read(unit=iin,fmt=200,iostat=ios) string_read
+ if(ios /= 0) stop 'error while reading input file'
+
+! suppress leading white spaces, if any
+ string_read = adjustl(string_read)
+
+! exit loop when we find the first line that is not a comment
+ if(string_read(1:1) /= '#') exit
+
+ enddo
+
+ if(ignore_junk) then
+ read(string_read,100) junk,value_to_read
+ else
+ read(string_read,*) value_to_read
+ endif
+
+! format
+ 100 format(a,f12.5)
+ 200 format(a100)
+
+ end subroutine read_value_double_precision
+
+
+
+
+
+
+
+
+
+!--------------------
+
+ subroutine read_value_logical(iin,ignore_junk,value_to_read)
+
+ implicit none
+
+ integer iin
+ logical ignore_junk
+
+ logical value_to_read
+
+ integer ios
+
+ character(len=100) string_read
+ character(len=34) junk
+
+ do
+ read(unit=iin,fmt=200,iostat=ios) string_read
+ if(ios /= 0) stop 'error while reading input file'
+
+! suppress leading white spaces, if any
+ string_read = adjustl(string_read)
+
+! exit loop when we find the first line that is not a comment
+ if(string_read(1:1) /= '#') exit
+
+ enddo
+
+ if(ignore_junk) then
+ read(string_read,100) junk,value_to_read
+ else
+ read(string_read,*) value_to_read
+ endif
+
+! format
+ 100 format(a,l8)
+ 200 format(a100)
+
+ end subroutine read_value_logical
+
+
+
+
+
+
+
+
+
+!--------------------
+
+ subroutine read_value_string(iin,ignore_junk,value_to_read)
+
+ implicit none
+
+ integer iin
+ logical ignore_junk
+
+ character(len=*) value_to_read
+
+ integer ios
+
+ character(len=100) string_read
+ character(len=34) junk
+
+ do
+ read(unit=iin,fmt=200,iostat=ios) string_read
+ if(ios /= 0) stop 'error while reading input file'
+
+! suppress leading white spaces, if any
+ string_read = adjustl(string_read)
+
+! exit loop when we find the first line that is not a comment
+ if(string_read(1:1) /= '#') exit
+
+ enddo
+
+ if(ignore_junk) then
+ read(string_read,100) junk,value_to_read
+ else
+ read(string_read,*) value_to_read
+ endif
+
+! format
+ 100 format(a34,a)
+ 200 format(a100)
+
+ end subroutine read_value_string
+
+
+
+
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/Makefile
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/Makefile 2004-12-18 20:09:40 UTC (rev 8437)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/Makefile 2007-12-07 23:46:22 UTC (rev 8438)
@@ -37,11 +37,14 @@
default: all
clean:
- /bin/rm -f $(EXEC) $(EXEC).trace $O/*.o *.o $O/*.il *.mod core *.gnu *.ps Ux*.dat Uz*.dat image*.pnm sources;
+ /bin/rm -f $(EXEC) $(EXEC).trace $O/*.o *.o $O/*.il *.mod core *.gnu *.ps Ux*.dat Uz*.dat Ux*.bin Uz*.bin image*.pnm xconvolve_source_timefunction *receiver_line_* sources;
all: $(OBJS)
$(LINK) $(FLAGS) -o $(EXEC) $(OBJS)
+convolve_source_timefunction: $O/convolve_source_timefunction.o
+ ${F90} $(FLAGS) -o xconvolve_source_timefunction $O/convolve_source_timefunction.o
+
$O/checkgrid.o: checkgrid.f90 constants.h
${F90} $(FLAGS) -c -o $O/checkgrid.o checkgrid.f90
@@ -51,6 +54,9 @@
$O/createnum_slow.o: createnum_slow.f90 constants.h
${F90} $(FLAGS) -c -o $O/createnum_slow.o createnum_slow.f90
+$O/convolve_source_timefunction.o: convolve_source_timefunction.f90
+ ${F90} $(FLAGS) -c -o $O/convolve_source_timefunction.o convolve_source_timefunction.f90
+
$O/datim.o: datim.f90 constants.h
${F90} $(FLAGS) -c -o $O/datim.o datim.f90
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/checkgrid.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/checkgrid.f90 2004-12-18 20:09:40 UTC (rev 8437)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/checkgrid.f90 2007-12-07 23:46:22 UTC (rev 8438)
@@ -11,7 +11,7 @@
!
!========================================================================
- subroutine checkgrid(deltat,gltfu,initialfield,rsizemin,rsizemax, &
+ subroutine checkgrid(deltat,f0,t0,initialfield,rsizemin,rsizemax, &
cpoverdxmin,cpoverdxmax,rlamdaSmin,rlamdaSmax,rlamdaPmin,rlamdaPmax)
!
@@ -22,14 +22,12 @@
include "constants.h"
- double precision gltfu(20)
+ double precision f0,t0
double precision deltat,rsizemin,rsizemax,cpoverdxmin,cpoverdxmax, &
rlamdaSmin,rlamdaSmax,rlamdaPmin,rlamdaPmax
logical initialfield
- double precision f0,t0
-
!
!---- verification taille de grille min et max
!
@@ -49,9 +47,6 @@
if(.not. initialfield) then
- f0 = gltfu(5)
- t0 = gltfu(6)
-
print *,' Onset time = ',t0
print *,' Fundamental period = ',1.d0/f0
print *,' Fundamental frequency = ',f0
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/defarrays.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/defarrays.f90 2004-12-18 20:09:40 UTC (rev 8437)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/defarrays.f90 2007-12-07 23:46:22 UTC (rev 8438)
@@ -13,9 +13,9 @@
subroutine defarrays(vpext,vsext,rhoext,density,elastcoef, &
xigll,zigll,xix,xiz,gammax,gammaz,a11,a12, &
- ibool,kmato,coord,gltfu,npoin,rsizemin,rsizemax, &
- cpoverdxmin,cpoverdxmax,rlambdaSmin,rlambdaSmax, &
- rlambdaPmin,rlambdaPmax,vpmin,vpmax,ireadmodel,nspec,numat)
+ ibool,kmato,coord,npoin,rsizemin,rsizemax, &
+ cpoverdxmin,cpoverdxmax,rlambdaSmin,rlambdaSmax,rlambdaPmin,rlambdaPmax, &
+ vpmin,vpmax,readmodel,nspec,numat,source_type,ix_source,iz_source,ispec_source)
! define all the arrays for the variational formulation
@@ -24,7 +24,7 @@
include "constants.h"
integer i,j,ispec,material,ipointnum,npoin,nspec,numat
- integer isourx,isourz,ielems,ir,is
+ integer ix_source,iz_source,ispec_source,ir,is,source_type
integer kmato(nspec),ibool(NGLLX,NGLLX,nspec)
@@ -41,7 +41,6 @@
double precision xigll(NGLLX),zigll(NGLLZ)
- double precision gltfu(20)
double precision vpext(npoin)
double precision vsext(npoin)
double precision rhoext(npoin)
@@ -56,7 +55,7 @@
double precision rsizemin,rsizemax,cpoverdxmin,cpoverdxmax, &
rlambdaSmin,rlambdaSmax,rlambdaPmin,rlambdaPmax,vpmin,vpmax
- logical ireadmodel
+ logical readmodel
double precision, external :: lagrange_deriv_GLL
@@ -104,7 +103,7 @@
do i=1,NGLLX
!--- si formulation heterogene pour un modele de vitesse externe
- if(ireadmodel) then
+ if(readmodel) then
ipointnum = ibool(i,j,ispec)
cploc = vpext(ipointnum)
csloc = vsext(ipointnum)
@@ -182,13 +181,9 @@
print *
! seulement si source explosive
- if(nint(gltfu(2)) == 2) then
+ if(source_type == 2) then
- isourx = nint(gltfu(10))
- isourz = nint(gltfu(11))
- ielems = nint(gltfu(12))
-
- if(isourx == 1 .or. isourx == NGLLX .or. isourz == 1 .or. isourz == NGLLX) &
+ if(ix_source == 1 .or. ix_source == NGLLX .or. iz_source == 1 .or. iz_source == NGLLX) &
stop 'Explosive source on element edge'
!---- definir a11 et a12 - dirac (schema en croix)
@@ -196,15 +191,15 @@
sig0 = one
do ir=1,NGLLX
- flagxprime = lagrange_deriv_GLL(ir-1,isourx-1,xigll,NGLLX)
- a11(ir,isourz) = a11(ir,isourz) + sig0*xix(isourx,isourz,ielems)*flagxprime
- a12(ir,isourz) = a12(ir,isourz) + sig0*xiz(isourx,isourz,ielems)*flagxprime
+ flagxprime = lagrange_deriv_GLL(ir-1,ix_source-1,xigll,NGLLX)
+ a11(ir,iz_source) = a11(ir,iz_source) + sig0*xix(ix_source,iz_source,ispec_source)*flagxprime
+ a12(ir,iz_source) = a12(ir,iz_source) + sig0*xiz(ix_source,iz_source,ispec_source)*flagxprime
enddo
do is=1,NGLLZ
- flagzprime = lagrange_deriv_GLL(is-1,isourz-1,zigll,NGLLZ)
- a11(isourx,is) = a11(isourx,is) + sig0*gammax(isourx,isourz,ielems)*flagzprime
- a12(isourx,is) = a12(isourx,is) + sig0*gammaz(isourx,isourz,ielems)*flagzprime
+ flagzprime = lagrange_deriv_GLL(is-1,iz_source-1,zigll,NGLLZ)
+ a11(ix_source,is) = a11(ix_source,is) + sig0*gammax(ix_source,iz_source,ispec_source)*flagzprime
+ a12(ix_source,is) = a12(ix_source,is) + sig0*gammaz(ix_source,iz_source,ispec_source)*flagzprime
enddo
endif
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/plotpost.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/plotpost.f90 2004-12-18 20:09:40 UTC (rev 8437)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/plotpost.f90 2007-12-07 23:46:22 UTC (rev 8438)
@@ -11,12 +11,12 @@
!
!========================================================================
- subroutine plotpost(displ,coord,vpext,gltfl,posrec,it,dt,coorg, &
+ subroutine plotpost(displ,coord,vpext,iglob_source,iglob_rec,it,dt,coorg, &
xinterp,zinterp,shapeint, &
Uxinterp,Uzinterp,flagrange,density,elastcoef,knods,kmato,ibool, &
numabs,codeabs,anyabs,stitle,npoin,npgeo,vpmin,vpmax,nrec, &
- icolor,inumber,isubsamp,ivecttype,interpol,imeshvect,imodelvect, &
- iboundvect,ireadmodel,cutvect,nelemabs,numat,iptsdisp,nspec,ngnod)
+ colors,numbers,subsamp,vecttype,interpol,meshvect,modelvect, &
+ boundvect,readmodel,cutvect,nelemabs,numat,iptsdisp,nspec,ngnod)
!
! routine affichage postscript
@@ -31,7 +31,7 @@
double precision, dimension(MAXCOLORS) :: red,green,blue
integer it,nrec,nelemabs,numat,iptsdisp,nspec
- integer i,iglobrec,iglobsource,npoin,npgeo,ngnod
+ integer i,iglob_source,npoin,npgeo,ngnod
integer kmato(nspec),knods(ngnod,nspec)
integer ibool(NGLLX,NGLLZ,nspec)
@@ -48,8 +48,7 @@
double precision vpext(npoin)
double precision coorg(NDIME,npgeo)
- double precision gltfl(20)
- double precision posrec(NDIME,nrec)
+ integer iglob_rec(nrec)
integer numabs(nelemabs),codeabs(4,nelemabs)
logical anyabs
@@ -67,8 +66,8 @@
integer k,j,ispec,material,is,ir,nbcols,imat,icol,l,longueur
integer indice,ii,ipoin,in,nnum,ispecabs,ideb,ifin,ibord
- integer icolor,inumber,isubsamp,ivecttype
- logical interpol,imeshvect,imodelvect,iboundvect,ireadmodel
+ integer colors,numbers,subsamp,vecttype
+ logical interpol,meshvect,modelvect,boundvect,readmodel
double precision cutvect
double precision rapp_page,dispmax,xmin,zmin
@@ -305,7 +304,7 @@
write(24,*) '%'
write(24,*) '/Times-Roman findfont'
write(24,*) '.6 CM SCSF'
- if(icolor == 1) write(24,*) '.4 .9 .9 setrgbcolor'
+ if(colors == 1) write(24,*) '.4 .9 .9 setrgbcolor'
write(24,*) '11 CM 1.1 CM MV'
write(24,*) '(X axis) show'
write(24,*) '%'
@@ -316,15 +315,15 @@
write(24,*) '%'
write(24,*) '/Times-Roman findfont'
write(24,*) '.7 CM SCSF'
- if(icolor == 1) write(24,*) '.8 0 .8 setrgbcolor'
+ if(colors == 1) write(24,*) '.8 0 .8 setrgbcolor'
write(24,*) '24.35 CM 18.9 CM MV'
write(24,*) usoffset,' CM 2 div neg 0 MR'
write(24,*) 'currentpoint gsave translate -90 rotate 0 0 moveto'
- if(ivecttype == 1) then
+ if(vecttype == 1) then
write(24,*) '(Displacement vector field) show'
- else if(ivecttype == 2) then
+ else if(vecttype == 2) then
write(24,*) '(Velocity vector field) show'
- else if(ivecttype == 3) then
+ else if(vecttype == 3) then
write(24,*) '(Acceleration vector field) show'
else
stop 'Bad field code in PostScript display'
@@ -358,14 +357,14 @@
!
!---- draw the velocity model in background
!
- if(imodelvect) then
+ if(modelvect) then
do ispec=1,nspec
- do i=1,NGLLX-isubsamp,isubsamp
- do j=1,NGLLX-isubsamp,isubsamp
+ do i=1,NGLLX-subsamp,subsamp
+ do j=1,NGLLX-subsamp,subsamp
if((vpmax-vpmin)/vpmin > 0.02d0) then
- if(ireadmodel) then
+ if(readmodel) then
x1 = (vpext(ibool(i,j,ispec))-vpmin)/ (vpmax-vpmin)
else
material = kmato(ispec)
@@ -395,24 +394,24 @@
zw = zw * centim
write(24,500) xw,zw
- xw = coord(1,ibool(i+isubsamp,j,ispec))
- zw = coord(2,ibool(i+isubsamp,j,ispec))
+ xw = coord(1,ibool(i+subsamp,j,ispec))
+ zw = coord(2,ibool(i+subsamp,j,ispec))
xw = (xw-xmin)*rapp_page + orig_x
zw = (zw-zmin)*rapp_page + orig_z
xw = xw * centim
zw = zw * centim
write(24,499) xw,zw
- xw = coord(1,ibool(i+isubsamp,j+isubsamp,ispec))
- zw = coord(2,ibool(i+isubsamp,j+isubsamp,ispec))
+ xw = coord(1,ibool(i+subsamp,j+subsamp,ispec))
+ zw = coord(2,ibool(i+subsamp,j+subsamp,ispec))
xw = (xw-xmin)*rapp_page + orig_x
zw = (zw-zmin)*rapp_page + orig_z
xw = xw * centim
zw = zw * centim
write(24,499) xw,zw
- xw = coord(1,ibool(i,j+isubsamp,ispec))
- zw = coord(2,ibool(i,j+isubsamp,ispec))
+ xw = coord(1,ibool(i,j+subsamp,ispec))
+ zw = coord(2,ibool(i,j+subsamp,ispec))
xw = (xw-xmin)*rapp_page + orig_x
zw = (zw-zmin)*rapp_page + orig_z
xw = xw * centim
@@ -430,7 +429,7 @@
!---- draw spectral element mesh
!
- if(imeshvect) then
+ if(meshvect) then
write(24,*) '%'
write(24,*) '% spectral element mesh'
@@ -538,7 +537,7 @@
write(24,*) 'CO'
- if(icolor == 1) then
+ if(colors == 1) then
! For the moment 20 different colors max
nbcols = 20
@@ -551,15 +550,14 @@
endif
- if(imodelvect) then
+ if(modelvect) then
write(24,*) 'GC'
else
write(24,*) 'GG'
endif
-! write the element number, the group number and the
-! material number inside the element
- if(inumber == 1) then
+! write the element number, the group number and the material number inside the element
+ if(numbers == 1) then
xw = (coorg(1,knods(1,ispec)) + coorg(1,knods(2,ispec)) + &
coorg(1,knods(3,ispec)) + coorg(1,knods(4,ispec))) / 4.d0
@@ -569,7 +567,7 @@
zw = (zw-zmin)*rapp_page + orig_z
xw = xw * centim
zw = zw * centim
- if(icolor == 1) write(24,*) '1 setgray'
+ if(colors == 1) write(24,*) '1 setgray'
write(24,500) xw,zw
@@ -586,7 +584,7 @@
!---- draw the boundary conditions
!
- if(anyabs .and. iboundvect) then
+ if(anyabs .and. boundvect) then
write(24,*) '%'
write(24,*) '% boundary conditions on the mesh'
@@ -668,7 +666,7 @@
write(24,*) '%'
! fleches en couleur si modele de vitesse en background
- if(imodelvect) then
+ if(modelvect) then
write(24,*) 'Colvects'
else
write(24,*) '0 setgray'
@@ -836,7 +834,7 @@
write(24,*) '0 setgray'
! sources et recepteurs en couleur si modele de vitesse
- if(imodelvect) then
+ if(modelvect) then
write(24,*) 'Colreceiv'
else
write(24,*) '0 setgray'
@@ -845,10 +843,8 @@
!
!---- write position of the source
!
- iglobsource = nint(gltfl(9))
-
- xw = coord(1,iglobsource)
- zw = coord(2,iglobsource)
+ xw = coord(1,iglob_source)
+ zw = coord(2,iglob_source)
xw = (xw-xmin)*rapp_page + orig_x
zw = (zw-zmin)*rapp_page + orig_z
xw = xw * centim
@@ -867,9 +863,8 @@
if(i == 1) write(24,*) '% debut ligne recepteurs'
if(i == nrec) write(24,*) '% fin ligne recepteurs'
- iglobrec = nint(posrec(1,i))
- xw = coord(1,iglobrec)
- zw = coord(2,iglobrec)
+ xw = coord(1,iglob_rec(i))
+ zw = coord(2,iglob_rec(i))
xw = (xw-xmin)*rapp_page + orig_x
zw = (zw-zmin)*rapp_page + orig_z
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/positrec.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/positrec.f90 2004-12-18 20:09:40 UTC (rev 8437)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/positrec.f90 2007-12-07 23:46:22 UTC (rev 8438)
@@ -11,7 +11,7 @@
!
!========================================================================
- subroutine positrec(coord,posrec,npoin,nrec)
+ subroutine positrec(coord,posrec,iglob_rec,npoin,nrec)
!
!---- calculer la position reelle des recepteurs
@@ -24,21 +24,22 @@
integer npoin,nrec
double precision coord(NDIME,npoin)
double precision posrec(NDIME,nrec)
+ integer iglob_rec(nrec)
double precision distminmax,distmin,xs,zs,xp,zp,dist
- integer n,ip,ipoint
+ integer irec,ip
write(iout,200)
distminmax = -HUGEVAL
- do n=1,nrec
+ do irec=1,nrec
distmin = +HUGEVAL
! coordonnees demandees
- xs = posrec(1,n)
- zs = posrec(2,n)
+ xs = posrec(1,irec)
+ zs = posrec(2,irec)
do ip=1,npoin
@@ -51,18 +52,15 @@
! retenir le point pour lequel l'ecart est minimal
if(dist < distmin) then
distmin = dist
- ipoint = ip
+ iglob_rec(irec) = ip
endif
enddo
distminmax = max(distmin,distminmax)
- write(iout,150) n,xs,zs,coord(1,ipoint),coord(2,ipoint),distmin
+ write(iout,150) irec,xs,zs,coord(1,iglob_rec(irec)),coord(2,iglob_rec(irec)),distmin
-! stocker numero global dans premiere coordonnee
- posrec(1,n) = dble(ipoint)
-
enddo
write(iout,160) distminmax
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/positsource.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/positsource.f90 2004-12-18 20:09:40 UTC (rev 8437)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/positsource.f90 2007-12-07 23:46:22 UTC (rev 8438)
@@ -11,7 +11,7 @@
!
!========================================================================
- subroutine positsource(coord,ibool,gltfu,npoin,nspec)
+ subroutine positsource(coord,ibool,npoin,nspec,xs,zs,source_type,ix_source,iz_source,ispec_source,iglob_source)
!
!----- calculer la position reelle de la source
@@ -21,44 +21,42 @@
include "constants.h"
- integer npoin,nspec
- double precision coord(NDIME,npoin)
- double precision gltfu(20)
+ integer npoin,nspec,source_type
integer ibool(NGLLX,NGLLZ,nspec)
- double precision distminmax,distmin,xs,zs,xp,zp,dist
- integer ip,ipoint,ix,iy,numelem,ilowx,ilowy,ihighx,ihighy
+ double precision xs,zs
+ double precision coord(NDIME,npoin)
+ integer ip,ix,iz,numelem,ilowx,ilowz,ihighx,ihighz,ix_source,iz_source,ispec_source,iglob_source
+
+ double precision distminmax,distmin,xp,zp,dist
+
write(iout,200)
distminmax = -HUGEVAL
distmin = +HUGEVAL
-! coordonnees demandees pour la source
- xs = gltfu(3)
- zs = gltfu(4)
-
ilowx = 1
- ilowy = 1
+ ilowz = 1
ihighx = NGLLX
- ihighy = NGLLZ
+ ihighz = NGLLZ
! on ne fait la recherche que sur l'interieur de l'element si source explosive
- if(nint(gltfu(2)) == 2) then
+ if(source_type == 2) then
ilowx = 2
- ilowy = 2
+ ilowz = 2
ihighx = NGLLX-1
- ihighy = NGLLZ-1
+ ihighz = NGLLZ-1
endif
! recherche du point de grille le plus proche
do numelem=1,nspec
do ix=ilowx,ihighx
- do iy=ilowy,ihighy
+ do iz=ilowz,ihighz
! numero global du point
- ip=ibool(ix,iy,numelem)
+ ip=ibool(ix,iz,numelem)
! coordonnees du point de grille
xp = coord(1,ip)
@@ -69,21 +67,19 @@
! retenir le point pour lequel l'ecart est minimal
if(dist < distmin) then
distmin = dist
- gltfu(9) = ip
- gltfu(10) = ix
- gltfu(11) = iy
- gltfu(12) = numelem
+ iglob_source = ip
+ ix_source = ix
+ iz_source = iz
+ ispec_source = numelem
endif
enddo
enddo
enddo
- ipoint = nint(gltfu(9))
-
distminmax = max(distmin,distminmax)
- write(iout,150) xs,zs,coord(1,ipoint),coord(2,ipoint),distmin
+ write(iout,150) xs,zs,coord(1,iglob_source),coord(2,iglob_source),distmin
write(iout,160) distminmax
150 format(1x,f12.3,1x,f12.3,1x,f12.3,1x,f12.3,f12.3)
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/specfem2D.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/specfem2D.f90 2004-12-18 20:09:40 UTC (rev 8437)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/specfem2D.f90 2007-12-07 23:46:22 UTC (rev 8438)
@@ -19,8 +19,10 @@
!
! version 5.1, December 2004 :
+! - Dirac and Gaussian time sources and corresponding convolution routine
! - more general mesher with any number of curved layers
! - color PNM snapshots
+! - more flexible Par file with any number of comment lines
!
! version 5.0, May 2004 :
! - got rid of useless routines, suppressed commons etc.
@@ -45,17 +47,20 @@
character(len=80) datlin
- double precision gltfu(20)
+ integer source_type,time_function_type
+ double precision xs,zs,f0,t0,factor,angleforce
double precision, dimension(:,:), allocatable :: coorg,posrec
double precision, dimension(:), allocatable :: coorgread
double precision, dimension(:), allocatable :: posrecread
+ integer, dimension(:), allocatable :: iglob_rec
+
double precision, dimension(:,:), allocatable :: sisux,sisuz
logical anyabs
- integer i,j,it,irec,iglobrec,ipoin,ip,id,ip1,ip2,in,nnum
+ integer i,j,it,irec,ipoin,ip,id,ip1,ip2,in,nnum
integer nbpoin,inump,n,npoinext,ispec,npoin,npgeo,iglob
double precision valux,valuz,rhoextread,vpextread,vsextread
@@ -113,19 +118,19 @@
integer, dimension(:,:), allocatable :: knods
integer, dimension(:), allocatable :: kmato,numabs
- integer ie,k,iexplo
+ integer ie,k
- integer ielems,iglobsource
- double precision f0,t0,factor,a,angleforce,ricker,displnorm_all
+ integer ispec_source,iglob_source,ix_source,iz_source
+ double precision a,source_time_function,displnorm_all
double precision rsizemin,rsizemax,cpoverdxmin,cpoverdxmax, &
lambdal_Smin,lambdal_Smax,lambdal_Pmin,lambdal_Pmax,vpmin,vpmax
- integer icolor,inumber,isubsamp,ivecttype,itaff,nrec,isismostype
+ integer colors,numbers,subsamp,vecttype,itaff,nrec,sismostype
integer numat,ngnod,nspec,iptsdisp,nelemabs
- logical interpol,imeshvect,imodelvect,iboundvect,ireadmodel,initialfield, &
- ioutputgrid,ignuplot,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON
+ logical interpol,meshvect,modelvect,boundvect,readmodel,initialfield, &
+ outputgrid,gnuplot,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON
double precision cutvect,anglerec
@@ -208,13 +213,13 @@
read(IIN,*) npgeo
read(IIN,40) datlin
- read(IIN,*) ignuplot,interpol
+ read(IIN,*) gnuplot,interpol
read(IIN,40) datlin
- read(IIN,*) itaff,icolor,inumber
+ read(IIN,*) itaff,colors,numbers
read(IIN,40) datlin
- read(IIN,*) imeshvect,imodelvect,iboundvect,cutvect,isubsamp,nx_sem_PNM
+ read(IIN,*) meshvect,modelvect,boundvect,cutvect,subsamp,nx_sem_PNM
cutvect = cutvect / 100.d0
read(IIN,40) datlin
@@ -224,17 +229,17 @@
read(IIN,*) initialfield
read(IIN,40) datlin
- read(IIN,*) isismostype,ivecttype
+ read(IIN,*) sismostype,vecttype
read(IIN,40) datlin
- read(IIN,*) ireadmodel,ioutputgrid,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON
+ read(IIN,*) readmodel,outputgrid,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON
!---- check parameters read
write(IOUT,200) npgeo,NDIME
- write(IOUT,600) itaff,icolor,inumber
- write(IOUT,700) nrec,isismostype,anglerec
- write(IOUT,750) initialfield,ireadmodel,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON,ioutputgrid
- write(IOUT,800) ivecttype,100.d0*cutvect,isubsamp
+ write(IOUT,600) itaff,colors,numbers
+ write(IOUT,700) nrec,sismostype,anglerec
+ write(IOUT,750) initialfield,readmodel,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON,outputgrid
+ write(IOUT,800) vecttype,100.d0*cutvect,subsamp
!---- read time step
read(IIN,40) datlin
@@ -249,29 +254,32 @@
allocate(sisuz(NSTEP,nrec))
allocate(posrec(NDIME,nrec))
allocate(coorg(NDIME,npgeo))
+ allocate(iglob_rec(nrec))
!
!---- read source information
!
read(IIN,40) datlin
- read(IIN,*) (gltfu(k), k=2,9)
+ read(IIN,*) source_type,time_function_type,xs,zs,f0,t0,factor,angleforce
!
!----- check the input
!
if(.not. initialfield) then
- iexplo = nint(gltfu(2))
- if (iexplo == 1) then
- write(IOUT,212) (gltfu(k), k=3,8)
- else if(iexplo == 2) then
- write(IOUT,222) (gltfu(k), k=3,7)
+ if (source_type == 1) then
+ write(IOUT,212) xs,zs,f0,t0,factor,angleforce
+ else if(source_type == 2) then
+ write(IOUT,222) xs,zs,f0,t0,factor
else
stop 'Unknown source type number !'
endif
endif
+! for the source time function
+ a = pi*pi*f0*f0
+
!----- convert angle from degrees to radians
- gltfu(8) = gltfu(8) * pi / 180.d0
+ angleforce = angleforce * pi / 180.d0
!
!---- read receiver locations
@@ -449,7 +457,7 @@
allocate(fglobx(npoin))
allocate(fglobz(npoin))
- if(ireadmodel) then
+ if(readmodel) then
npoinext = npoin
else
npoinext = 1
@@ -486,7 +494,7 @@
!
!--- save the grid of points in a file
!
- if(ioutputgrid) then
+ if(outputgrid) then
print *
print *,'Saving the grid in a text file...'
print *
@@ -501,7 +509,7 @@
!
!----- plot the GLL mesh in a Gnuplot file
!
- if(ignuplot) call plotgll(knods,ibool,coorg,coord,npoin,npgeo,ngnod,nspec)
+ if(gnuplot) call plotgll(knods,ibool,coorg,coord,npoin,npgeo,ngnod,nspec)
!
!---- define coefficients of the Newmark time scheme
@@ -512,13 +520,13 @@
!
!---- definir la position reelle des points source et recepteurs
!
- call positsource(coord,ibool,gltfu,npoin,nspec)
- call positrec(coord,posrec,npoin,nrec)
+ call positsource(coord,ibool,npoin,nspec,xs,zs,source_type,ix_source,iz_source,ispec_source,iglob_source)
+ call positrec(coord,posrec,iglob_rec,npoin,nrec)
!
!---- eventuellement lecture d'un modele externe de vitesse et de densite
!
- if(ireadmodel) then
+ if(readmodel) then
print *
print *,'Reading velocity and density model from external file...'
print *
@@ -540,9 +548,9 @@
!
call defarrays(vpext,vsext,rhoext,density,elastcoef, &
xigll,yigll,xix,xiz,gammax,gammaz,a11,a12, &
- ibool,kmato,coord,gltfu,npoin,rsizemin,rsizemax, &
- cpoverdxmin,cpoverdxmax,lambdal_Smin,lambdal_Smax, &
- lambdal_Pmin,lambdal_Pmax,vpmin,vpmax,ireadmodel,nspec,numat)
+ ibool,kmato,coord,npoin,rsizemin,rsizemax, &
+ cpoverdxmin,cpoverdxmax,lambdal_Smin,lambdal_Smax,lambdal_Pmin,lambdal_Pmax, &
+ vpmin,vpmax,readmodel,nspec,numat,source_type,ix_source,iz_source,ispec_source)
! build the global mass matrix once and for all
rmass(:) = ZERO
@@ -551,7 +559,7 @@
do i = 1,NGLLX
iglob = ibool(i,j,ispec)
!--- if external density model
- if(ireadmodel) then
+ if(readmodel) then
rhol = rhoext(iglob)
else
rhol = density(kmato(ispec))
@@ -566,9 +574,10 @@
!
!---- 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)
!
- call checkgrid(deltat,gltfu,initialfield,rsizemin,rsizemax, &
- cpoverdxmin,cpoverdxmax,lambdal_Smin,lambdal_Smax,lambdal_Pmin,lambdal_Pmax)
+ if(time_function_type /= 3) call checkgrid(deltat,f0,t0,initialfield, &
+ rsizemin,rsizemax,cpoverdxmin,cpoverdxmax,lambdal_Smin,lambdal_Smax,lambdal_Pmin,lambdal_Pmax)
!
!---- for color PNM images
@@ -803,7 +812,7 @@
do i = 1,NGLLX
!--- if external medium, get elastic parameters of current grid point
- if(ireadmodel) then
+ if(readmodel) then
iglob = ibool(i,j,ispec)
cpl = vpext(iglob)
csl = vsext(iglob)
@@ -995,7 +1004,7 @@
zeta = xix(i,j,ispec) * jacobian(i,j,ispec)
! external velocity model
- if(ireadmodel) then
+ if(readmodel) then
cpl = vpext(iglob)
csl = vsext(iglob)
rhol = rhoext(iglob)
@@ -1036,7 +1045,7 @@
zeta = xix(i,j,ispec) * jacobian(i,j,ispec)
! external velocity model
- if(ireadmodel) then
+ if(readmodel) then
cpl = vpext(iglob)
csl = vsext(iglob)
rhol = rhoext(iglob)
@@ -1083,7 +1092,7 @@
xxi = gammaz(i,j,ispec) * jacobian(i,j,ispec)
! external velocity model
- if(ireadmodel) then
+ if(readmodel) then
cpl = vpext(iglob)
csl = vsext(iglob)
rhol = rhoext(iglob)
@@ -1130,7 +1139,7 @@
xxi = gammaz(i,j,ispec) * jacobian(i,j,ispec)
! external velocity model
- if(ireadmodel) then
+ if(readmodel) then
cpl = vpext(iglob)
csl = vsext(iglob)
rhol = rhoext(iglob)
@@ -1167,30 +1176,30 @@
! --- add the source
if(.not. initialfield) then
- f0 = gltfu(5)
- t0 = gltfu(6)
- factor = gltfu(7)
- angleforce = gltfu(8)
+! Ricker source time function
+ if(time_function_type == 1) then
+ source_time_function = - factor * (ONE-TWO*a*(time-t0)**2) * exp(-a*(time-t0)**2)
-! Ricker wavelet for the source time function
- a = pi*pi*f0*f0
- ricker = - factor * (ONE-TWO*a*(time-t0)**2)*exp(-a*(time-t0)**2)
+! Gaussian or Dirac (we use a very thin Gaussian instead) source time function
+ else if(time_function_type == 2 .or. time_function_type == 3) then
+ source_time_function = factor * exp(-a*(time-t0)**2)
+ else
+ stop 'unknown source time function'
+ endif
+
! --- collocated force
- if(nint(gltfu(2)) == 1) then
- iglobsource = nint(gltfu(9))
- accel(1,iglobsource) = accel(1,iglobsource) - dsin(angleforce)*ricker
- accel(2,iglobsource) = accel(2,iglobsource) + dcos(angleforce)*ricker
+ if(source_type == 1) then
+ accel(1,iglob_source) = accel(1,iglob_source) - dsin(angleforce)*source_time_function
+ accel(2,iglob_source) = accel(2,iglob_source) + dcos(angleforce)*source_time_function
!---- explosion
- else if(nint(gltfu(2)) == 2) then
-! recuperer numero d'element de la source
- ielems = nint(gltfu(12))
+ else if(source_type == 2) then
do i=1,NGLLX
do j=1,NGLLX
- iglob = ibool(i,j,ielems)
- accel(1,iglob) = accel(1,iglob) + a11(i,j)*ricker
- accel(2,iglob) = accel(2,iglob) + a12(i,j)*ricker
+ iglob = ibool(i,j,ispec_source)
+ accel(1,iglob) = accel(1,iglob) + a11(i,j)*source_time_function
+ accel(2,iglob) = accel(2,iglob) + a12(i,j)*source_time_function
enddo
enddo
endif
@@ -1321,24 +1330,23 @@
! store the seismograms
do irec=1,nrec
- iglobrec = nint(posrec(1,irec))
- if(isismostype == 1) then
- valux = displ(1,iglobrec)
- valuz = displ(2,iglobrec)
- else if(isismostype == 2) then
- valux = veloc(1,iglobrec)
- valuz = veloc(2,iglobrec)
- else if(isismostype == 3) then
- valux = accel(1,iglobrec)
- valuz = accel(2,iglobrec)
- else
- stop 'Wrong field code for seismogram output'
- endif
+ if(sismostype == 1) then
+ valux = displ(1,iglob_rec(irec))
+ valuz = displ(2,iglob_rec(irec))
+ else if(sismostype == 2) then
+ valux = veloc(1,iglob_rec(irec))
+ valuz = veloc(2,iglob_rec(irec))
+ else if(sismostype == 3) then
+ valux = accel(1,iglob_rec(irec))
+ valuz = accel(2,iglob_rec(irec))
+ else
+ stop 'Wrong field code for seismogram output'
+ endif
! rotation eventuelle des composantes
- sisux(it,irec) = dcosrot*valux + dsinrot*valuz
- sisuz(it,irec) = - dsinrot*valux + dcosrot*valuz
+ sisux(it,irec) = dcosrot*valux + dsinrot*valuz
+ sisuz(it,irec) = - dsinrot*valux + dcosrot*valuz
enddo
@@ -1359,30 +1367,30 @@
!---- affichage postscript
!
write(IOUT,*) 'Dump PostScript'
- if(ivecttype == 1) then
+ if(vecttype == 1) then
write(IOUT,*) 'drawing displacement field...'
- call plotpost(displ,coord,vpext,gltfu,posrec, &
+ call plotpost(displ,coord,vpext,iglob_source,iglob_rec, &
it,deltat,coorg,xinterp,zinterp,shapeint, &
Uxinterp,Uzinterp,flagrange,density,elastcoef,knods,kmato,ibool, &
numabs,codeabs,anyabs,stitle,npoin,npgeo,vpmin,vpmax,nrec, &
- icolor,inumber,isubsamp,ivecttype,interpol,imeshvect,imodelvect, &
- iboundvect,ireadmodel,cutvect,nelemabs,numat,iptsdisp,nspec,ngnod)
- else if(ivecttype == 2) then
+ colors,numbers,subsamp,vecttype,interpol,meshvect,modelvect, &
+ boundvect,readmodel,cutvect,nelemabs,numat,iptsdisp,nspec,ngnod)
+ else if(vecttype == 2) then
write(IOUT,*) 'drawing velocity field...'
- call plotpost(veloc,coord,vpext,gltfu,posrec, &
+ call plotpost(veloc,coord,vpext,iglob_source,iglob_rec, &
it,deltat,coorg,xinterp,zinterp,shapeint, &
Uxinterp,Uzinterp,flagrange,density,elastcoef,knods,kmato,ibool, &
numabs,codeabs,anyabs,stitle,npoin,npgeo,vpmin,vpmax,nrec, &
- icolor,inumber,isubsamp,ivecttype,interpol,imeshvect,imodelvect, &
- iboundvect,ireadmodel,cutvect,nelemabs,numat,iptsdisp,nspec,ngnod)
- else if(ivecttype == 3) then
+ colors,numbers,subsamp,vecttype,interpol,meshvect,modelvect, &
+ boundvect,readmodel,cutvect,nelemabs,numat,iptsdisp,nspec,ngnod)
+ else if(vecttype == 3) then
write(IOUT,*) 'drawing acceleration field...'
- call plotpost(accel,coord,vpext,gltfu,posrec, &
+ call plotpost(accel,coord,vpext,iglob_source,iglob_rec, &
it,deltat,coorg,xinterp,zinterp,shapeint, &
Uxinterp,Uzinterp,flagrange,density,elastcoef,knods,kmato,ibool, &
numabs,codeabs,anyabs,stitle,npoin,npgeo,vpmin,vpmax,nrec, &
- icolor,inumber,isubsamp,ivecttype,interpol,imeshvect,imodelvect, &
- iboundvect,ireadmodel,cutvect,nelemabs,numat,iptsdisp,nspec,ngnod)
+ colors,numbers,subsamp,vecttype,interpol,meshvect,modelvect, &
+ boundvect,readmodel,cutvect,nelemabs,numat,iptsdisp,nspec,ngnod)
else
stop 'wrong field code for PostScript display'
endif
@@ -1407,14 +1415,14 @@
write(IOUT,*) 'Fin creation image PNM'
!---- save temporary seismograms
- call write_seismograms(sisux,sisuz,NSTEP,nrec,deltat)
+ call write_seismograms(sisux,sisuz,NSTEP,nrec,deltat,sismostype,iglob_rec,coord,npoin)
endif
enddo ! end of the main time loop
!---- save final seismograms
- call write_seismograms(sisux,sisuz,NSTEP,nrec,deltat)
+ call write_seismograms(sisux,sisuz,NSTEP,nrec,deltat,sismostype,iglob_rec,coord,npoin)
! print exit banner
call datim(stitle)
@@ -1440,26 +1448,26 @@
'Number of space dimensions . . . . . . . . . (NDIME) =',i8)
600 format(//1x,'C o n t r o l',/1x,34('='),//5x, &
'Display frequency . . . . . . . . . . . . . (itaff) = ',i5/ 5x, &
- 'Color display . . . . . . . . . . . . . . . (icolor) = ',i5/ 5x, &
+ 'Color display . . . . . . . . . . . . . . . (colors) = ',i5/ 5x, &
' == 0 black and white display ', / 5x, &
' == 1 color display ', /5x, &
- 'Numbered mesh . . . . . . . . . . . . . . .(inumber) = ',i5/ 5x, &
+ 'Numbered mesh . . . . . . . . . . . . . . .(numbers) = ',i5/ 5x, &
' == 0 do not number the mesh ', /5x, &
' == 1 number the mesh ')
700 format(//1x,'C o n t r o l',/1x,34('='),//5x, &
'Total number of receivers. . . . . . . . . . .(nrec) = ',i6/5x, &
- 'Seismograms recording type. . . . . . .(isismostype) = ',i6/5x, &
+ 'Seismograms recording type. . . . . . .(sismostype) = ',i6/5x, &
'Angle for first line of receivers. . . . .(anglerec) = ',f6.2)
750 format(//1x,'C o n t r o l',/1x,34('='),//5x, &
'Read external initial field or not . .(initialfield) = ',l6/5x, &
- 'Read external velocity model or not. . .(ireadmodel) = ',l6/5x, &
+ 'Read external velocity model or not. . .(readmodel) = ',l6/5x, &
'Turn anisotropy on or off. . . .(TURN_ANISOTROPY_ON) = ',l6/5x, &
'Turn attenuation on or off. . .(TURN_ATTENUATION_ON) = ',l6/5x, &
- 'Save grid in external file or not . . .(ioutputgrid) = ',l6)
+ 'Save grid in external file or not . . .(outputgrid) = ',l6)
800 format(//1x,'C o n t r o l',/1x,34('='),//5x, &
- 'Vector display type . . . . . . . . . . .(ivecttype) = ',i6/5x, &
+ 'Vector display type . . . . . . . . . . .(vecttype) = ',i6/5x, &
'Percentage of cut for vector plots. . . . .(cutvect) = ',f6.2/5x, &
- 'Subsampling for velocity model display . .(isubsamp) = ',i6)
+ 'Subsampling for velocity model display . .(subsamp) = ',i6)
703 format(//' I t e r a t i o n s '/1x,29('='),//5x, &
'Number of time iterations . . . . .(NSTEP) =',i8,/5x, &
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/write_seismograms.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/write_seismograms.f90 2004-12-18 20:09:40 UTC (rev 8437)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/write_seismograms.f90 2007-12-07 23:46:22 UTC (rev 8438)
@@ -11,22 +11,32 @@
!
!========================================================================
- subroutine write_seismograms(sisux,sisuz,nt,nrec,deltat)
+ subroutine write_seismograms(sisux,sisuz,nt,nrec,deltat,sismostype,iglob_rec,coord,npoin)
! save the seismograms in ASCII format
implicit none
- integer nt,nrec
+ include "constants.h"
+
+ integer nt,nrec,sismostype,npoin
double precision deltat
+ integer iglob_rec(nrec)
+
double precision sisux(nt,nrec)
double precision sisuz(nt,nrec)
+ double precision coord(NDIME,npoin)
integer irec,it
character(len=100) name
+! scaling factor for Seismic Unix xsu dislay
+ double precision, parameter :: FACTORXSU = 1.d0
+
+! write seismograms in ASCII format
+
! X component
do irec=1,nrec
write(name,221) irec
@@ -47,8 +57,81 @@
close(11)
enddo
- 221 format('Ux_file_',i3.3,'.dat')
- 222 format('Uz_file_',i3.3,'.dat')
+!----
+! write seismograms in single precision binary format
+
+! X component
+ open(unit=11,file='Ux_file.bin',status='unknown',access='direct',recl=nt*nrec*4)
+ write(11,rec=1) (sngl(sisux(irec,1)),irec=1,nt*nrec)
+ close(11)
+
+! Z component
+ open(unit=11,file='Uz_file.bin',status='unknown',access='direct',recl=nt*nrec*4)
+ write(11,rec=1) (sngl(sisuz(irec,1)),irec=1,nt*nrec)
+ close(11)
+
+!----
+
+! ligne de recepteurs pour Xsu
+ open(unit=11,file='receiver_line_Xsu_XWindow',status='unknown')
+
+ write(11,110) FACTORXSU,nt,deltat,nrec
+
+ do irec=1,nrec
+ write(11,140) coord(1,iglob_rec(irec))
+ if(irec < nrec) write(11,*) ','
+ enddo
+
+ if(sismostype == 1) then
+ write(11,*) '@title="Ux at displacement@component"@<@Ux_file.bin'
+ else if(sismostype == 2) then
+ write(11,*) '@title="Ux at velocity@component"@<@Ux_file.bin'
+ else
+ write(11,*) '@title="Ux at acceleration@component"@<@Ux_file.bin'
+ endif
+
+ close(11)
+
+! script de visualisation
+ open(unit=11,file='show_receiver_line_Xsu',status='unknown')
+ write(11,100)
+ write(11,*)
+ write(11,*) '/bin/rm -f tempfile receiver_line_Xsu_postscript'
+ write(11,*) '# concatener toutes les lignes'
+ write(11,*) 'tr -d ''\012'' <receiver_line_Xsu_XWindow >tempfile'
+ write(11,*) '# remettre fin de ligne'
+ write(11,*) 'echo " " >> tempfile'
+ write(11,*) '# supprimer espaces, changer arobas, dupliquer'
+ write(11,120)
+ write(11,*) '/bin/rm -f tempfile'
+ write(11,*) '# copier fichier pour sortie postscript'
+ write(11,130)
+ write(11,*) '/bin/rm -f tempfile'
+ write(11,*) 'echo ''rm -f uxpoly.ps uzpoly.ps'' > tempfile'
+ write(11,*) 'cat tempfile receiver_line_Xsu_postscript > tempfile2'
+ write(11,*) '/bin/mv -f tempfile2 receiver_line_Xsu_postscript'
+ write(11,*) '/bin/rm -f tempfile'
+ write(11,*) '# executer commande xsu'
+ write(11,*) 'sh receiver_line_Xsu_XWindow'
+ write(11,*) '/bin/rm -f tempfile tempfile2'
+ close(11)
+
+! formats
+ 100 format('#!/bin/csh -f')
+
+ 110 format('xwigb at xcur=',f8.2,'@n1=',i5,'@d1=',f15.8,'@label1="Time@(s)"@label2="x@(m)"@n2=',i5,'@x2=')
+
+ 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')
+
+ 140 format(f12.5)
+
+ 221 format('Ux_file_',i3.3,'.dat')
+ 222 format('Uz_file_',i3.3,'.dat')
+
end subroutine write_seismograms
More information about the cig-commits
mailing list