[cig-commits] r8440 - in seismo/2D/SPECFEM2D/trunk: MAILLE90
SPECFEM90
walter at geodynamics.org
walter at geodynamics.org
Fri Dec 7 15:46:33 PST 2007
Author: walter
Date: 2007-12-07 15:46:32 -0800 (Fri, 07 Dec 2007)
New Revision: 8440
Added:
seismo/2D/SPECFEM2D/trunk/MAILLE90/InterfaceSinus.dat
seismo/2D/SPECFEM2D/trunk/MAILLE90/Par_Paul_acoustic
seismo/2D/SPECFEM2D/trunk/MAILLE90/Par_Paul_elastic
seismo/2D/SPECFEM2D/trunk/SPECFEM90/compute_gradient_fluid.f90
Modified:
seismo/2D/SPECFEM2D/trunk/MAILLE90/Makefile
seismo/2D/SPECFEM2D/trunk/MAILLE90/Par
seismo/2D/SPECFEM2D/trunk/MAILLE90/Par_cours_M2_UPPA
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/compute_gradient_attenuation.f90
seismo/2D/SPECFEM2D/trunk/SPECFEM90/createnum_fast.f90
seismo/2D/SPECFEM2D/trunk/SPECFEM90/createnum_slow.f90
seismo/2D/SPECFEM2D/trunk/SPECFEM90/cree_image_PNM.f90
seismo/2D/SPECFEM2D/trunk/SPECFEM90/datim.f90
seismo/2D/SPECFEM2D/trunk/SPECFEM90/defarrays.f90
seismo/2D/SPECFEM2D/trunk/SPECFEM90/define_derivative_matrices.f90
seismo/2D/SPECFEM2D/trunk/SPECFEM90/gll_library.f90
seismo/2D/SPECFEM2D/trunk/SPECFEM90/gmat01.f90
seismo/2D/SPECFEM2D/trunk/SPECFEM90/lagrange_poly.f90
seismo/2D/SPECFEM2D/trunk/SPECFEM90/plotgll.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/q49shape.f90
seismo/2D/SPECFEM2D/trunk/SPECFEM90/q49spec.f90
seismo/2D/SPECFEM2D/trunk/SPECFEM90/specfem2D.f90
seismo/2D/SPECFEM2D/trunk/SPECFEM90/write_seismograms.f90
Log:
modified version 5.1 of SEM_2D_Dimitri to add option for acoustic simulations instead of elastic
Added: seismo/2D/SPECFEM2D/trunk/MAILLE90/InterfaceSinus.dat
===================================================================
--- seismo/2D/SPECFEM2D/trunk/MAILLE90/InterfaceSinus.dat 2005-01-09 13:38:40 UTC (rev 8439)
+++ seismo/2D/SPECFEM2D/trunk/MAILLE90/InterfaceSinus.dat 2007-12-07 23:46:32 UTC (rev 8440)
@@ -0,0 +1,75 @@
+#
+# Number of interfaces
+#
+3
+#
+#
+#
+#
+# Interface 1 (Bottom of the mesh)
+#
+2
+0 0
+4000 0
+#
+# Interface 2
+#
+41
+0 1500
+100 1507
+200 1510
+300 1507
+400 1500
+500 1492
+600 1490
+700 1492
+800 1500
+900 1507
+1000 1510
+1100 1507
+1200 1500
+1300 1492
+1400 1490
+1500 1492
+1600 1500
+1700 1507
+1800 1510
+1900 1507
+2000 1500
+2100 1492
+2200 1490
+2300 1492
+2400 1500
+2500 1507
+2600 1510
+2700 1507
+2800 1500
+2900 1492
+3000 1490
+3100 1492
+3200 1500
+3300 1507
+3400 1510
+3500 1507
+3600 1500
+3700 1492
+3800 1490
+3900 1492
+4000 1500
+#
+# Interface 3
+#
+2
+0 2000
+4000 2000
+#
+#Number of spectral elements
+#
+#
+#
+#
+30
+#
+#
+#
+10
Modified: seismo/2D/SPECFEM2D/trunk/MAILLE90/Makefile
===================================================================
--- seismo/2D/SPECFEM2D/trunk/MAILLE90/Makefile 2005-01-09 13:38:40 UTC (rev 8439)
+++ seismo/2D/SPECFEM2D/trunk/MAILLE90/Makefile 2007-12-07 23:46:32 UTC (rev 8440)
@@ -13,7 +13,7 @@
# Intel Linux
F90 = ifort
-FLAGS=-O0 -implicitnone -warn stderrors -warn truncated_source -warn argument_checking -warn unused -warn declarations -std95 -check bounds -CB
+FLAGS=-O0 -implicitnone -warn stderrors -warn truncated_source -warn argument_checking -warn unused -warn declarations -std95 -check bounds
#FLAGS=-O2 -implicitnone -warn stderrors -warn truncated_source -warn argument_checking -warn unused -warn declarations -std95 -check nobounds
#
Modified: seismo/2D/SPECFEM2D/trunk/MAILLE90/Par
===================================================================
--- seismo/2D/SPECFEM2D/trunk/MAILLE90/Par 2005-01-09 13:38:40 UTC (rev 8439)
+++ seismo/2D/SPECFEM2D/trunk/MAILLE90/Par 2007-12-07 23:46:32 UTC (rev 8440)
@@ -7,8 +7,8 @@
#
# title of job, and file that contains interface data
#
-title = Test pour cours M2 UPPA
-interfacesfile = interfaces_cours_M2_UPPA_curved.dat
+title = Test Sinusoide
+interfacesfile = InterfaceSinus.dat
#
# geometry of the model (origin lower-left corner = 0,0) and mesh description
#
@@ -18,6 +18,7 @@
ngnod = 4 ! noeuds de controle pour blocs (4 ou 9)
initialfield = .false. ! use a plane wave as source or not
readmodel = .false. ! read external earth model or not
+ELASTIC = .false. ! elastic or acoustic simulation
TURN_ANISOTROPY_ON = .false. ! turn anisotropy on or off
TURN_ATTENUATION_ON = .false. ! turn attenuation on or off
#
@@ -30,25 +31,25 @@
#
# time step parameters
#
-nt = 1600 ! nb total de pas de temps
+nt = 1400 ! nb total de pas de temps
dt = 1.d-3 ! valeur du pas de temps
#
# source parameters
#
-source_surf = .true. ! source dans le volume ou a la surface
+source_surf = .false. ! source dans le volume ou a la surface
xs = 2000. ! source location x in meters
-zs = 1500. ! source location z in meters
+zs = 1200. ! source location z in meters
source_type = 1 ! force = 1 or explosion = 2
-time_function_type = 1 ! Ricker = 1, Gaussian = 2, Dirac = 3
+time_function_type = 2 ! Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4
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
+factor = 1.d4 ! amplification factor
#
# receiver line parameters
#
enreg_surf = .true. ! enregistrement volume ou surface
-sismostype = 1 ! record 1=displ 2=veloc 3=accel
-nrec = 11 ! number of receivers
+sismostype = 2 ! record 1=displ 2=veloc 3=accel
+nrec = 100 ! number of receivers
xdeb = 300. ! first receiver x in meters
zdeb = 2200. ! first receiver z in meters
xfin = 3700. ! last receiver x in meters
@@ -58,7 +59,7 @@
# display parameters
#
itaff = 100 ! display frequency in time steps
-vecttype = 1 ! display 1=displ 2=veloc 3=accel
+vecttype = 2 ! display 1=displ 2=veloc 3=accel
cutvect = 1. ! amplitude min en % pour vector plots
meshvect = .true. ! display mesh on vector plots or not
modelvect = .false. ! display velocity model on vector plots
@@ -71,11 +72,11 @@
#
# velocity and density model (nx,nz)
#
-nbmodels = 3 ! nb de modeles differents (rho,vp,vs)
-1 0 2700.d0 3000.d0 1732.051d0 0 0
-2 0 2500.d0 2700.d0 1558.845d0 0 0
-3 0 2200.d0 2500.d0 1443.375d0 0 0
-nbzone = 3 ! nb of zones and model number for each
-1 80 1 20 1
-1 80 21 40 2
-1 80 41 60 3
+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
Added: seismo/2D/SPECFEM2D/trunk/MAILLE90/Par_Paul_acoustic
===================================================================
--- seismo/2D/SPECFEM2D/trunk/MAILLE90/Par_Paul_acoustic 2005-01-09 13:38:40 UTC (rev 8439)
+++ seismo/2D/SPECFEM2D/trunk/MAILLE90/Par_Paul_acoustic 2007-12-07 23:46:32 UTC (rev 8440)
@@ -0,0 +1,82 @@
+# ----------------------------------------------------------------
+#
+# 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
+interfacesfile = InterfaceSinus.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
+ngnod = 4 ! noeuds de controle pour blocs (4 ou 9)
+initialfield = .false. ! use a plane wave as source or not
+readmodel = .false. ! read external earth model or not
+ELASTIC = .false. ! 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 = 1400 ! nb total de pas de temps
+dt = 1.d-3 ! valeur du pas de temps
+#
+# source parameters
+#
+source_surf = .false. ! source dans le volume ou a la surface
+xs = 2000. ! source location x in meters
+zs = 1200. ! source location z in meters
+source_type = 1 ! force = 1 or explosion = 2
+time_function_type = 2 ! Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4
+f0 = 10.0 ! dominant source frequency (Hz) if not Dirac
+angle = 0. ! angle of the source (for a force only)
+factor = 1.d4 ! amplification factor
+#
+# receiver line parameters
+#
+enreg_surf = .true. ! enregistrement volume ou surface
+sismostype = 2 ! record 1=displ 2=veloc 3=accel
+nrec = 100 ! 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
+#
+# display parameters
+#
+itaff = 100 ! display frequency in time steps
+vecttype = 2 ! display 1=displ 2=veloc 3=accel
+cutvect = 1. ! amplitude min en % pour vector plots
+meshvect = .true. ! display mesh on vector plots or not
+modelvect = .false. ! display velocity model on vector plots
+boundvect = .true. ! display boundary conditions on plots
+interpol = .true. ! interpolation of the display or not
+pointsdisp = 6 ! points for interpolation of display
+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)
+#
+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
Added: seismo/2D/SPECFEM2D/trunk/MAILLE90/Par_Paul_elastic
===================================================================
--- seismo/2D/SPECFEM2D/trunk/MAILLE90/Par_Paul_elastic 2005-01-09 13:38:40 UTC (rev 8439)
+++ seismo/2D/SPECFEM2D/trunk/MAILLE90/Par_Paul_elastic 2007-12-07 23:46:32 UTC (rev 8440)
@@ -0,0 +1,82 @@
+# ----------------------------------------------------------------
+#
+# 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
+interfacesfile = InterfaceSinus.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
+ngnod = 4 ! noeuds de controle pour blocs (4 ou 9)
+initialfield = .false. ! use a plane wave as source or not
+readmodel = .false. ! read external earth model or not
+ELASTIC = .true. ! elastic or acoustic simulation
+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 = 1400 ! nb total de pas de temps
+dt = 1.d-3 ! valeur du pas de temps
+#
+# source parameters
+#
+source_surf = .false. ! source dans le volume ou a la surface
+xs = 2000. ! source location x in meters
+zs = 1200. ! source location z in meters
+source_type = 1 ! force = 1 or explosion = 2
+time_function_type = 1 ! Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4
+f0 = 10.0 ! dominant source frequency (Hz) if not Dirac
+angle = 0. ! angle of the source (for a force only)
+factor = 1.d4 ! amplification factor
+#
+# receiver line parameters
+#
+enreg_surf = .true. ! enregistrement volume ou surface
+sismostype = 1 ! record 1=displ 2=veloc 3=accel
+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
+#
+# display parameters
+#
+itaff = 100 ! display frequency in time steps
+vecttype = 1 ! display 1=displ 2=veloc 3=accel
+cutvect = 1. ! amplitude min en % pour vector plots
+meshvect = .true. ! display mesh on vector plots or not
+modelvect = .false. ! display velocity model on vector plots
+boundvect = .true. ! display boundary conditions on plots
+interpol = .true. ! interpolation of the display or not
+pointsdisp = 6 ! points for interpolation of display
+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)
+#
+nbmodels = 2 ! nb de modeles differents (0,rho,vp,vs,0,0)
+1 0 2700.d0 3000.d0 1732.051d0 0 0
+2 0 1800.d0 2000.d0 1100.845d0 0 0
+#3 0 2200.d0 2500.d0 1443.375d0 0 0
+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/MAILLE90/Par_cours_M2_UPPA
===================================================================
--- seismo/2D/SPECFEM2D/trunk/MAILLE90/Par_cours_M2_UPPA 2005-01-09 13:38:40 UTC (rev 8439)
+++ seismo/2D/SPECFEM2D/trunk/MAILLE90/Par_cours_M2_UPPA 2007-12-07 23:46:32 UTC (rev 8440)
@@ -18,6 +18,7 @@
ngnod = 4 ! noeuds de controle pour blocs (4 ou 9)
initialfield = .false. ! use a plane wave as source or not
readmodel = .false. ! read external earth model or not
+ELASTIC = .true. ! elastic or acoustic simulation
TURN_ANISOTROPY_ON = .false. ! turn anisotropy on or off
TURN_ATTENUATION_ON = .false. ! turn attenuation on or off
#
@@ -39,7 +40,7 @@
xs = 2000. ! source location x in meters
zs = 1500. ! source location z in meters
source_type = 1 ! force = 1 or explosion = 2
-time_function_type = 1 ! Ricker = 1, Gaussian = 2, Dirac = 3
+time_function_type = 1 ! Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4
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
@@ -71,7 +72,7 @@
#
# velocity and density model (nx,nz)
#
-nbmodels = 3 ! nb de modeles differents (rho,vp,vs)
+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
3 0 2200.d0 2500.d0 1443.375d0 0 0
Modified: seismo/2D/SPECFEM2D/trunk/MAILLE90/meshfem2D.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/MAILLE90/meshfem2D.f90 2005-01-09 13:38:40 UTC (rev 8439)
+++ seismo/2D/SPECFEM2D/trunk/MAILLE90/meshfem2D.f90 2007-12-07 23:46:32 UTC (rev 8440)
@@ -7,7 +7,7 @@
! Dimitri Komatitsch
! Universite de Pau et des Pays de l'Adour, France
!
-! (c) December 2004
+! (c) January 2005
!
!========================================================================
@@ -51,8 +51,8 @@
character(len=50) interfacesfile,title
- integer imatnum,inumabs,inumelem
- integer nelemabs,npgeo,nspec
+ integer imatnum,inumabs,inumsurface,inumelem
+ integer nelemabs,nelemsurface,npgeo,nspec
integer k,icol,ili,istepx,istepz,ix,iz,irec,i,j
integer ixdebzone,ixfinzone,izdebzone,izfinzone,imodnum
integer izone,imodele,nbzone,nbmodeles
@@ -63,7 +63,7 @@
logical codehaut,codebas,codegauche,codedroite
- double precision tang1,tangN,vpzone,vszone
+ double precision tang1,tangN,vpzone,vszone,poisson_ratio
double precision cutvect,xspacerec,zspacerec
double precision anglerec,xfin,zfin,xdeb,zdeb,xmin,xmax,dt
double precision rhoread,cpread,csread,aniso3read,aniso4read
@@ -71,7 +71,7 @@
logical interpol,gnuplot,readmodel,outputgrid
logical abshaut,absbas,absgauche,absdroite
logical source_surf,enreg_surf,meshvect,initialfield,modelvect,boundvect
- logical TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON
+ logical ELASTIC,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON
integer, external :: num
double precision, external :: value_spline
@@ -110,6 +110,7 @@
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,ELASTIC)
call read_value_logical(IIN_PAR,IGNORE_JUNK,TURN_ANISOTROPY_ON)
call read_value_logical(IIN_PAR,IGNORE_JUNK,TURN_ATTENUATION_ON)
@@ -195,7 +196,7 @@
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)
+ if(time_function_type == 4) f0 = 1.d0 / (5.d0 * dt)
! time delay of the source in seconds, use a 20 % security margin
t0 = 1.20d0 / f0
@@ -205,7 +206,7 @@
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=Gaussian, 3=Dirac): ',time_function_type
+ print *,'Time function type (1=Ricker, 2=First derivative, 3=Gaussian, 4=Dirac): ',time_function_type
print *,'Angle of the source if force = ',angle
print *,'Multiplying factor = ',factor
@@ -266,15 +267,19 @@
do imodele=1,nbmodeles
read(IIN_PAR,*) i,icodematread,rhoread,cpread,csread,aniso3read,aniso4read
- if(i < 1 .or. i > nbmodeles) stop 'Wrong material point number'
+ if(i < 1 .or. i > nbmodeles) stop 'Wrong model number!!'
icodemat(i) = icodematread
rho(i) = rhoread
cp(i) = cpread
cs(i) = csread
+
+ if(rho(i) <= 0.d0 .or. cp(i) <= 0.d0 .or. cs(i) < 0.d0) stop 'negative value of velocity or density'
+
+! check that Cs = 0 if acoustic simulation
+ if(.not. ELASTIC .and. cs(i) > 0.0001) stop 'must have Cs = 0 for acoustic model'
+
aniso3(i) = aniso3read
aniso4(i) = aniso4read
- if(i <= 0) stop 'Negative model number not allowed !!'
- if(rho(i) <= 0.d0 .or. cp(i) <= 0.d0 .or. cs(i) < 0.d0) stop 'negative value of velocity or density'
enddo
print *
@@ -321,8 +326,9 @@
print *,'vp = ',vpzone
print *,'vs = ',vszone
print *,'rho = ',rho(imodnum)
- print *,'Poisson''s ratio = ', &
- 0.5d0*(vpzone*vpzone-2.d0*vszone*vszone) / (vpzone*vpzone-vszone*vszone)
+ poisson_ratio = 0.5d0*(vpzone*vpzone-2.d0*vszone*vszone) / (vpzone*vpzone-vszone*vszone)
+ print *,'Poisson''s ratio = ',poisson_ratio
+ if(poisson_ratio <= -1.00001d0 .or. poisson_ratio >= 0.50001d0) stop 'incorrect value of Poisson''s ratio'
else
print *,'Model # ',imodnum,' anisotrope'
print *,'c11 = ',cp(imodnum)
@@ -372,6 +378,20 @@
!---
+! perform basic checks on parameters read
+
+! for acoustic
+ if(TURN_ANISOTROPY_ON .and. .not. ELASTIC) stop 'currently cannot have anisotropy in acoustic simulation'
+
+ if(TURN_ATTENUATION_ON .and. .not. ELASTIC) stop 'currently cannot have attenuation in acoustic simulation'
+
+ if(source_type == 2 .and. .not. ELASTIC) stop 'currently cannot have moment tensor source in acoustic simulation'
+
+! for attenuation
+ if(TURN_ANISOTROPY_ON .and. TURN_ATTENUATION_ON) stop 'cannot have anisotropy and attenuation both turned on in current version'
+
+!---
+
! allocate arrays for the grid
allocate(x(0:nx,0:nz))
allocate(z(0:nx,0:nz))
@@ -593,8 +613,8 @@
write(15,*) 'sismostype vecttype'
write(15,*) sismostype,vecttype
- write(15,*) 'readmodel outputgrid TURN_ANISOTROPY_ON TURN_ATTENUATION_ON'
- write(15,*) readmodel,outputgrid,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON
+ write(15,*) 'readmodel outputgrid ELASTIC TURN_ANISOTROPY_ON TURN_ATTENUATION_ON'
+ write(15,*) readmodel,outputgrid,ELASTIC,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON
write(15,*) 'ncycl dtinc'
write(15,*) nt,dt
@@ -617,7 +637,6 @@
!
!--- introduction des bords absorbants
!
-
nelemabs = 0
if(absbas) nelemabs = nelemabs + nx
if(abshaut) nelemabs = nelemabs + nx
@@ -633,9 +652,17 @@
if(abshaut .and. absgauche) nelemabs = nelemabs - 1
if(abshaut .and. absdroite) nelemabs = nelemabs - 1
- write(15,*) 'numat ngnod nspec pointsdisp nelemabs'
- write(15,*) nbmodeles,ngnod,nspec,pointsdisp,nelemabs
+!
+!--- introduction de la surface libre si milieu acoustique
+!
+ nelemsurface = nx
+! on a deux fois trop d'elements si elements 9 noeuds
+ if(ngnod == 9) nelemsurface = nelemsurface / 2
+
+ write(15,*) 'numat ngnod nspec pointsdisp nelemabs nelemsurface'
+ write(15,*) nbmodeles,ngnod,nspec,pointsdisp,nelemabs,nelemsurface
+
write(15,*) 'Material sets (num 0 rho vp vs 0 0) or (num 1 rho c11 c13 c33 c44)'
do i=1,nbmodeles
write(15,*) i,icodemat(i),rho(i),cp(i),cs(i),aniso3(i),aniso4(i)
@@ -671,7 +698,6 @@
!
!--- sauvegarde des bords absorbants
!
-
print *
print *,'Au total il y a ',nelemabs,' elements absorbants'
print *
@@ -706,6 +732,29 @@
enddo
endif
+!
+!--- sauvegarde de la surface libre
+!
+ print *
+ print *,'Au total il y a ',nelemsurface,' elements a la surface libre'
+ print *
+
+! generer la liste des elements a la surface libre
+ if(nelemsurface > 0) then
+ write(15,*) 'Liste des elements a la surface libre'
+ write(15,*) abshaut
+ inumsurface = 0
+ do iz = 1,nzread
+ do ix = 1,nxread
+ inumelem = (iz-1)*nxread + ix
+ if(iz == nzread) then
+ inumsurface = inumsurface + 1
+ write(15,*) inumsurface,inumelem
+ endif
+ enddo
+ enddo
+ endif
+
close(15)
10 format('')
@@ -718,7 +767,7 @@
! routines de maillage
! ********************
-! --- numero global du noeud
+!--- numero global du noeud
integer function num(i,j,nx)
@@ -730,7 +779,7 @@
end function num
-! --- representation des interfaces par un spline
+!--- representation des interfaces par un spline
double precision function value_spline(x,xinterface,zinterface,coefs_interface,npoints_interface)
@@ -830,18 +879,6 @@
end subroutine splint
-
-
-
-
-
-
-
-
-
-
-
-
!--------------------
subroutine read_value_integer(iin,ignore_junk,value_to_read)
@@ -882,15 +919,6 @@
end subroutine read_value_integer
-
-
-
-
-
-
-
-
-
!--------------------
subroutine read_value_double_precision(iin,ignore_junk,value_to_read)
@@ -931,14 +959,6 @@
end subroutine read_value_double_precision
-
-
-
-
-
-
-
-
!--------------------
subroutine read_value_logical(iin,ignore_junk,value_to_read)
@@ -979,14 +999,6 @@
end subroutine read_value_logical
-
-
-
-
-
-
-
-
!--------------------
subroutine read_value_string(iin,ignore_junk,value_to_read)
@@ -1027,6 +1039,3 @@
end subroutine read_value_string
-
-
-
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/Makefile
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/Makefile 2005-01-09 13:38:40 UTC (rev 8439)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/Makefile 2007-12-07 23:46:32 UTC (rev 8440)
@@ -13,7 +13,7 @@
# Intel Linux
F90 = ifort
-#FLAGS=-O0 -implicitnone -warn stderrors -warn truncated_source -warn argument_checking -warn unused -warn declarations -std95 -check bounds -CB
+#FLAGS=-O0 -implicitnone -warn stderrors -warn truncated_source -warn argument_checking -warn unused -warn declarations -std95 -check bounds
FLAGS=-O3 -implicitnone -warn stderrors -warn truncated_source -warn argument_checking -warn unused -warn declarations -std95 -check nobounds
#
@@ -32,7 +32,7 @@
OBJS = $O/checkgrid.o $O/datim.o $O/defarrays.o\
$O/lagrange_poly.o $O/gmat01.o $O/gll_library.o $O/plotgll.o $O/define_derivative_matrices.o\
$O/plotpost.o $O/positrec.o $O/positsource.o $O/q49spec.o $O/compute_gradient_attenuation.o\
- $O/specfem2D.o $O/write_seismograms.o $O/createnum_fast.o $O/createnum_slow.o $O/q49shape.o $O/cree_image_PNM.o
+ $O/specfem2D.o $O/write_seismograms.o $O/createnum_fast.o $O/createnum_slow.o $O/q49shape.o $O/cree_image_PNM.o $O/compute_gradient_fluid.o
default: all
@@ -99,6 +99,9 @@
$O/compute_gradient_attenuation.o: compute_gradient_attenuation.f90 constants.h
${F90} $(FLAGS) -c -o $O/compute_gradient_attenuation.o compute_gradient_attenuation.f90
+$O/compute_gradient_fluid.o: compute_gradient_fluid.f90 constants.h
+ ${F90} $(FLAGS) -c -o $O/compute_gradient_fluid.o compute_gradient_fluid.f90
+
$O/cree_image_PNM.o: cree_image_PNM.f90 constants.h
${F90} $(FLAGS) -c -o $O/cree_image_PNM.o cree_image_PNM.f90
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/checkgrid.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/checkgrid.f90 2005-01-09 13:38:40 UTC (rev 8439)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/checkgrid.f90 2007-12-07 23:46:32 UTC (rev 8440)
@@ -7,7 +7,7 @@
! Dimitri Komatitsch
! Universite de Pau et des Pays de l'Adour, France
!
-! (c) December 2004
+! (c) January 2005
!
!========================================================================
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/compute_gradient_attenuation.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/compute_gradient_attenuation.f90 2005-01-09 13:38:40 UTC (rev 8439)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/compute_gradient_attenuation.f90 2007-12-07 23:46:32 UTC (rev 8440)
@@ -7,7 +7,7 @@
! Dimitri Komatitsch
! Universite de Pau et des Pays de l'Adour, France
!
-! (c) December 2004
+! (c) January 2005
!
!========================================================================
Added: seismo/2D/SPECFEM2D/trunk/SPECFEM90/compute_gradient_fluid.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/compute_gradient_fluid.f90 2005-01-09 13:38:40 UTC (rev 8439)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/compute_gradient_fluid.f90 2007-12-07 23:46:32 UTC (rev 8440)
@@ -0,0 +1,83 @@
+
+!========================================================================
+!
+! S P E C F E M 2 D Version 5.1
+! ------------------------------
+!
+! Dimitri Komatitsch
+! Universite de Pau et des Pays de l'Adour, France
+!
+! (c) January 2005
+!
+!========================================================================
+
+ subroutine compute_gradient_fluid(potential,veloc_field_postscript, &
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,NSPEC,npoin)
+
+! compute Grad(potential) in fluid medium
+
+ implicit none
+
+ include "constants.h"
+
+ integer NSPEC,npoin
+
+ integer, dimension(NGLLX,NGLLZ,NSPEC) :: ibool
+
+ double precision, dimension(NGLLX,NGLLZ,NSPEC) :: xix,xiz,gammax,gammaz
+
+ double precision, dimension(NDIME,npoin) :: potential,veloc_field_postscript
+
+! array with derivatives of Lagrange polynomials
+ double precision, dimension(NGLLX,NGLLX) :: hprime_xx
+ double precision, dimension(NGLLZ,NGLLZ) :: hprime_zz
+
+! local variables
+ integer i,j,k,ispec,iglob
+
+! space derivatives
+ double precision tempx1l,tempx2l
+ double precision hp1,hp2
+
+! jacobian
+ double precision xixl,xizl,gammaxl,gammazl
+
+! loop over spectral elements
+ do ispec = 1,NSPEC
+
+! double loop over GLL to compute and store gradients
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+
+! derivative along x
+ tempx1l = ZERO
+ do k = 1,NGLLX
+ hp1 = hprime_xx(k,i)
+ iglob = ibool(k,j,ispec)
+ tempx1l = tempx1l + potential(1,iglob)*hp1
+ enddo
+
+! derivative along z
+ tempx2l = ZERO
+ do k = 1,NGLLZ
+ hp2 = hprime_zz(k,j)
+ iglob = ibool(i,k,ispec)
+ tempx2l = tempx2l + potential(1,iglob)*hp2
+ enddo
+
+ xixl = xix(i,j,ispec)
+ xizl = xiz(i,j,ispec)
+ gammaxl = gammax(i,j,ispec)
+ gammazl = gammaz(i,j,ispec)
+
+! derivatives of velocity potential
+ iglob = ibool(i,j,ispec)
+ veloc_field_postscript(1,iglob) = tempx1l*xixl + tempx2l*gammaxl
+ veloc_field_postscript(2,iglob) = tempx1l*xizl + tempx2l*gammazl
+
+ enddo
+ enddo
+ enddo
+
+ end subroutine compute_gradient_fluid
+
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/createnum_fast.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/createnum_fast.f90 2005-01-09 13:38:40 UTC (rev 8439)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/createnum_fast.f90 2007-12-07 23:46:32 UTC (rev 8440)
@@ -7,7 +7,7 @@
! Dimitri Komatitsch
! Universite de Pau et des Pays de l'Adour, France
!
-! (c) December 2004
+! (c) January 2005
!
!========================================================================
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/createnum_slow.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/createnum_slow.f90 2005-01-09 13:38:40 UTC (rev 8439)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/createnum_slow.f90 2007-12-07 23:46:32 UTC (rev 8440)
@@ -7,7 +7,7 @@
! Dimitri Komatitsch
! Universite de Pau et des Pays de l'Adour, France
!
-! (c) December 2004
+! (c) January 2005
!
!========================================================================
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/cree_image_PNM.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/cree_image_PNM.f90 2005-01-09 13:38:40 UTC (rev 8439)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/cree_image_PNM.f90 2007-12-07 23:46:32 UTC (rev 8440)
@@ -7,7 +7,7 @@
! Dimitri Komatitsch
! Universite de Pau et des Pays de l'Adour, France
!
-! (c) December 2004
+! (c) January 2005
!
!========================================================================
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/datim.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/datim.f90 2005-01-09 13:38:40 UTC (rev 8439)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/datim.f90 2007-12-07 23:46:32 UTC (rev 8440)
@@ -7,7 +7,7 @@
! Dimitri Komatitsch
! Universite de Pau et des Pays de l'Adour, France
!
-! (c) December 2004
+! (c) January 2005
!
!========================================================================
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/defarrays.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/defarrays.f90 2005-01-09 13:38:40 UTC (rev 8439)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/defarrays.f90 2007-12-07 23:46:32 UTC (rev 8440)
@@ -7,14 +7,14 @@
! Dimitri Komatitsch
! Universite de Pau et des Pays de l'Adour, France
!
-! (c) December 2004
+! (c) January 2005
!
!========================================================================
subroutine defarrays(vpext,vsext,rhoext,density,elastcoef, &
xigll,zigll,xix,xiz,gammax,gammaz,a11,a12, &
ibool,kmato,coord,npoin,rsizemin,rsizemax, &
- cpoverdxmin,cpoverdxmax,rlambdaSmin,rlambdaSmax,rlambdaPmin,rlambdaPmax, &
+ cpoverdxmin,cpoverdxmax,lambdaSmin,lambdaSmax,lambdaPmin,lambdaPmax, &
vpmin,vpmax,readmodel,nspec,numat,source_type,ix_source,iz_source,ispec_source)
! define all the arrays for the variational formulation
@@ -46,14 +46,14 @@
double precision rhoext(npoin)
double precision vsmin,vsmax,densmin,densmax
- double precision rKmod,rlambda,rmu,denst
- double precision rKvol,cploc,csloc,x0,z0
+ double precision lambdaplus2mu,lambda,mu,denst
+ double precision kappa,cploc,csloc,x0,z0
double precision x1,z1,x2,z2,rdist1,rdist2,rapportmin,rapportmax
- double precision rlambmin,rlambmax
+ double precision lambdamin,lambdamax
double precision flagxprime,flagzprime,sig0
double precision rsizemin,rsizemax,cpoverdxmin,cpoverdxmax, &
- rlambdaSmin,rlambdaSmax,rlambdaPmin,rlambdaPmax,vpmin,vpmax
+ lambdaSmin,lambdaSmax,lambdaPmin,lambdaPmax,vpmin,vpmax
logical readmodel
@@ -81,24 +81,25 @@
cpoverdxmin = HUGEVAL
cpoverdxmax = -HUGEVAL
- rlambdaPmin = HUGEVAL
- rlambdaSmin = HUGEVAL
- rlambdaPmax = -HUGEVAL
- rlambdaSmax = -HUGEVAL
+ lambdaPmin = HUGEVAL
+ lambdaSmin = HUGEVAL
+ lambdaPmax = -HUGEVAL
+ lambdaSmax = -HUGEVAL
do ispec=1,nspec
- material = kmato(ispec)
+ material = kmato(ispec)
- rlambda = elastcoef(1,material)
- rmu = elastcoef(2,material)
- rKmod = elastcoef(3,material)
- denst = density(material)
+ lambda = elastcoef(1,material)
+ mu = elastcoef(2,material)
+ lambdaplus2mu = elastcoef(3,material)
+ denst = density(material)
- rKvol = rlambda + 2.d0*rmu/3.d0
- cploc = dsqrt((rKvol + 4.d0*rmu/3.d0)/denst)
- csloc = dsqrt(rmu/denst)
+ kappa = lambda + 2.d0*mu/3.d0
+ cploc = sqrt((kappa + 4.d0*mu/3.d0)/denst)
+ csloc = sqrt(mu/denst)
+
do j=1,NGLLZ
do i=1,NGLLX
@@ -108,20 +109,20 @@
cploc = vpext(ipointnum)
csloc = vsext(ipointnum)
denst = rhoext(ipointnum)
- rmu = denst*csloc*csloc
- rlambda = denst*cploc*cploc - 2.d0*rmu
- rKmod = rlambda + 2.d0*rmu
+ mu = denst*csloc*csloc
+ lambda = denst*cploc*cploc - 2.d0*mu
+ lambdaplus2mu = lambda + 2.d0*mu
endif
!--- calculer min et max du modele de vitesse et densite
- vpmin = dmin1(vpmin,cploc)
- vpmax = dmax1(vpmax,cploc)
+ vpmin = min(vpmin,cploc)
+ vpmax = max(vpmax,cploc)
- vsmin = dmin1(vsmin,csloc)
- vsmax = dmax1(vsmax,csloc)
+ vsmin = min(vsmin,csloc)
+ vsmax = max(vsmax,csloc)
- densmin = dmin1(densmin,denst)
- densmax = dmax1(densmax,denst)
+ densmin = min(densmin,denst)
+ densmax = max(densmax,denst)
!--- stocker parametres pour verifs diverses
if(i < NGLLX .and. j < NGLLZ) then
@@ -133,18 +134,18 @@
x2 = coord(1,ibool(i,j+1,ispec))
z2 = coord(2,ibool(i,j+1,ispec))
- rdist1 = dsqrt((x1-x0)**2 + (z1-z0)**2)
- rdist2 = dsqrt((x2-x0)**2 + (z2-z0)**2)
+ rdist1 = sqrt((x1-x0)**2 + (z1-z0)**2)
+ rdist2 = sqrt((x2-x0)**2 + (z2-z0)**2)
- rsizemin = dmin1(rsizemin,rdist1)
- rsizemin = dmin1(rsizemin,rdist2)
- rsizemax = dmax1(rsizemax,rdist1)
- rsizemax = dmax1(rsizemax,rdist2)
+ rsizemin = min(rsizemin,rdist1)
+ rsizemin = min(rsizemin,rdist2)
+ rsizemax = max(rsizemax,rdist1)
+ rsizemax = max(rsizemax,rdist2)
- rapportmin = cploc / dmax1(rdist1,rdist2)
- rapportmax = cploc / dmin1(rdist1,rdist2)
- cpoverdxmin = dmin1(cpoverdxmin,rapportmin)
- cpoverdxmax = dmax1(cpoverdxmax,rapportmax)
+ rapportmin = cploc / max(rdist1,rdist2)
+ rapportmax = cploc / min(rdist1,rdist2)
+ cpoverdxmin = min(cpoverdxmin,rapportmin)
+ cpoverdxmax = max(cpoverdxmax,rapportmax)
x0 = coord(1,ibool(1,1,ispec))
z0 = coord(2,ibool(1,1,ispec))
@@ -153,18 +154,18 @@
x2 = coord(1,ibool(1,NGLLZ,ispec))
z2 = coord(2,ibool(1,NGLLZ,ispec))
- rdist1 = dsqrt((x1-x0)**2 + (z1-z0)**2)
- rdist2 = dsqrt((x2-x0)**2 + (z2-z0)**2)
+ rdist1 = sqrt((x1-x0)**2 + (z1-z0)**2)
+ rdist2 = sqrt((x2-x0)**2 + (z2-z0)**2)
- rlambmin = cploc/dmax1(rdist1,rdist2)
- rlambmax = cploc/dmin1(rdist1,rdist2)
- rlambdaPmin = dmin1(rlambdaPmin,rlambmin)
- rlambdaPmax = dmax1(rlambdaPmax,rlambmax)
+ lambdamin = cploc/max(rdist1,rdist2)
+ lambdamax = cploc/min(rdist1,rdist2)
+ lambdaPmin = min(lambdaPmin,lambdamin)
+ lambdaPmax = max(lambdaPmax,lambdamax)
- rlambmin = csloc/dmax1(rdist1,rdist2)
- rlambmax = csloc/dmin1(rdist1,rdist2)
- rlambdaSmin = dmin1(rlambdaSmin,rlambmin)
- rlambdaSmax = dmax1(rlambdaSmax,rlambmax)
+ lambdamin = csloc/max(rdist1,rdist2)
+ lambdamax = csloc/min(rdist1,rdist2)
+ lambdaSmin = min(lambdaSmin,lambdamin)
+ lambdaSmax = max(lambdaSmax,lambdamax)
endif
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/define_derivative_matrices.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/define_derivative_matrices.f90 2005-01-09 13:38:40 UTC (rev 8439)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/define_derivative_matrices.f90 2007-12-07 23:46:32 UTC (rev 8440)
@@ -7,7 +7,7 @@
! Dimitri Komatitsch
! Universite de Pau et des Pays de l'Adour, France
!
-! (c) December 2004
+! (c) January 2005
!
!========================================================================
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/gll_library.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/gll_library.f90 2005-01-09 13:38:40 UTC (rev 8439)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/gll_library.f90 2007-12-07 23:46:32 UTC (rev 8440)
@@ -123,13 +123,13 @@
gammaf = one
- if (x == -half) gammaf = -two*dsqrt(pi)
- if (x == half) gammaf = dsqrt(pi)
+ if (x == -half) gammaf = -two*sqrt(pi)
+ if (x == half) gammaf = sqrt(pi)
if (x == one ) gammaf = one
if (x == two ) gammaf = one
- if (x == 1.5d0) gammaf = dsqrt(pi)/2.d0
- if (x == 2.5d0) gammaf = 1.5d0*dsqrt(pi)/2.d0
- if (x == 3.5d0) gammaf = 2.5d0*1.5d0*dsqrt(pi)/2.d0
+ if (x == 1.5d0) gammaf = sqrt(pi)/2.d0
+ if (x == 2.5d0) gammaf = 1.5d0*sqrt(pi)/2.d0
+ if (x == 3.5d0) gammaf = 2.5d0*1.5d0*sqrt(pi)/2.d0
if (x == 3.d0 ) gammaf = 2.d0
if (x == 4.d0 ) gammaf = 6.d0
if (x == 5.d0 ) gammaf = 24.d0
@@ -173,15 +173,15 @@
xlast = 0.d0
n = np-1
- dth = 4.d0*datan(1.d0)/(2.d0*dble(n)+2.d0)
+ dth = 4.d0*atan(1.d0)/(2.d0*dble(n)+2.d0)
p = 0.d0
pd = 0.d0
jmin = 0
do j=1,np
if(j == 1) then
- x = dcos((2.d0*(dble(j)-1.d0)+1.d0)*dth)
+ x = cos((2.d0*(dble(j)-1.d0)+1.d0)*dth)
else
- x1 = dcos((2.d0*(dble(j)-1.d0)+1.d0)*dth)
+ x1 = cos((2.d0*(dble(j)-1.d0)+1.d0)*dth)
x2 = xlast
x = (x1+x2)/2.d0
endif
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/gmat01.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/gmat01.f90 2005-01-09 13:38:40 UTC (rev 8439)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/gmat01.f90 2007-12-07 23:46:32 UTC (rev 8440)
@@ -7,11 +7,11 @@
! Dimitri Komatitsch
! Universite de Pau et des Pays de l'Adour, France
!
-! (c) December 2004
+! (c) January 2005
!
!========================================================================
- subroutine gmat01(density,elastcoef,numat)
+ subroutine gmat01(density_array,elastcoef,numat)
! read properties of a 2D isotropic or anisotropic linear elastic element
@@ -20,20 +20,20 @@
include "constants.h"
character(len=80) datlin
- double precision Kmod,Kvol
+ double precision lambdaplus2mu,kappa
integer numat
- double precision density(numat),elastcoef(4,numat)
+ double precision density_array(numat),elastcoef(4,numat)
integer in,n,indic
- double precision young,poisson,denst,cp,cs,amu,a2mu,alam
+ double precision young,poisson,density,cp,cs,mu,two_mu,lambda
double precision val1,val2,val3,val4
double precision c11,c13,c33,c44
!
!---- loop over the different material sets
!
- density(:) = zero
+ density_array(:) = zero
elastcoef(:,:) = zero
write(iout,100) numat
@@ -41,21 +41,32 @@
read(iin ,40) datlin
do in = 1,numat
- read(iin ,*) n,indic,denst,val1,val2,val3,val4
+ read(iin ,*) n,indic,density,val1,val2,val3,val4
if(n<1 .or. n>numat) stop 'Wrong material set number'
!---- materiau isotrope, vitesse P et vitesse S donnees
if(indic == 0) then
+
+! P and S velocity
cp = val1
cs = val2
- amu = denst*cs*cs
- a2mu = 2.d0*amu
- alam = denst*cp*cp - a2mu
- Kmod = alam + a2mu
- Kvol = alam + a2mu/3.d0
- young = 9.d0*Kvol*amu/(3.d0*Kvol + amu)
- poisson = half*(3.d0*Kvol-a2mu)/(3.d0*Kvol+amu)
+
+! Lam'e parameters
+ lambdaplus2mu = density*cp*cp
+ mu = density*cs*cs
+ two_mu = 2.d0*mu
+ lambda = lambdaplus2mu - two_mu
+
+! bulk modulus K
+ kappa = lambda + two_mu/3.d0
+
+! Young modulus
+ young = 9.d0*kappa*mu/(3.d0*kappa + mu)
+
+! Poisson's ratio
+ poisson = half*(3.d0*kappa-two_mu)/(3.d0*kappa+mu)
+
! Poisson's ratio must be between -1 and +1/2
if (poisson < -1.d0 .or. poisson > 0.5d0) stop 'Poisson''s ratio out of range'
@@ -75,9 +86,9 @@
! Transverse anisotropic : c11, c13, c33, c44
!
if(indic == 0 .or. indic == 1) then
- elastcoef(1,n) = alam
- elastcoef(2,n) = amu
- elastcoef(3,n) = Kmod
+ elastcoef(1,n) = lambda
+ elastcoef(2,n) = mu
+ elastcoef(3,n) = lambdaplus2mu
elastcoef(4,n) = zero
else
elastcoef(1,n) = c11
@@ -86,16 +97,16 @@
elastcoef(4,n) = c44
endif
- density(n) = denst
+ density_array(n) = density
!
!---- check the input
!
if(indic == 0 .or. indic == 1) then
- write(iout,200) n,cp,cs,denst,poisson,alam,amu,Kvol,young
+ write(iout,200) n,cp,cs,density,poisson,lambda,mu,kappa,young
else
- write(iout,300) n,c11,c13,c33,c44,denst, &
- sqrt(c33/denst),sqrt(c11/denst),sqrt(c44/denst),sqrt(c44/denst)
+ write(iout,300) n,c11,c13,c33,c44,density, &
+ sqrt(c33/density),sqrt(c11/density),sqrt(c44/density),sqrt(c44/density)
endif
enddo
@@ -113,12 +124,12 @@
'Material set number. . . . . . . . (jmat) =',i5,/5x, &
'P-wave velocity. . . . . . . . . . . (cp) =',1pe15.8,/5x, &
'S-wave velocity. . . . . . . . . . . (cs) =',1pe15.8,/5x, &
- 'Mass density. . . . . . . . . . . (denst) =',1pe15.8,/5x, &
+ 'Mass density. . . . . . . . . . (density) =',1pe15.8,/5x, &
'Poisson''s ratio . . . . . . . .(poisson) =',1pe15.8,/5x, &
- 'First Lame parameter Lambda. . . . (alam) =',1pe15.8,/5x, &
- 'Second Lame parameter Mu. . . . . . (amu) =',1pe15.8,/5x, &
- 'Bulk modulus K . . . . . . . . . . (Kvol) =',1pe15.8,/5x, &
- 'Young''s modulus E . . . . . . . . (young) =',1pe15.8)
+ 'First Lame parameter Lambda. . . (lambda) =',1pe15.8,/5x, &
+ 'Second Lame parameter Mu. . . . . . .(mu) =',1pe15.8,/5x, &
+ 'Bulk modulus K . . . . . . . . . .(kappa) =',1pe15.8,/5x, &
+ 'Young''s modulus E . . . . . . . .(young) =',1pe15.8)
300 format(//5x,'-------------------------------------',/5x, &
'-- Transverse anisotropic material --',/5x, &
'-------------------------------------',/5x, &
@@ -127,7 +138,7 @@
'c13 coefficient (Pascal). . . . . . (c13) =',1pe15.8,/5x, &
'c33 coefficient (Pascal). . . . . . (c33) =',1pe15.8,/5x, &
'c44 coefficient (Pascal). . . . . . (c44) =',1pe15.8,/5x, &
- 'Mass density. . . . . . . . . . . (denst) =',1pe15.8,/5x, &
+ 'Mass density. . . . . . . . . . (density) =',1pe15.8,/5x, &
'Velocity of qP along vertical axis. . . . =',1pe15.8,/5x, &
'Velocity of qP along horizontal axis. . . =',1pe15.8,/5x, &
'Velocity of qSV along vertical axis . . . =',1pe15.8,/5x, &
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/lagrange_poly.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/lagrange_poly.f90 2005-01-09 13:38:40 UTC (rev 8439)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/lagrange_poly.f90 2007-12-07 23:46:32 UTC (rev 8440)
@@ -7,7 +7,7 @@
! Dimitri Komatitsch
! Universite de Pau et des Pays de l'Adour, France
!
-! (c) December 2004
+! (c) January 2005
!
!========================================================================
@@ -32,7 +32,7 @@
EPS = 1.d-5
DZ = Z - ZGLL(I)
- if(dABS(DZ) < EPS) then
+ if(abs(DZ) < EPS) then
HGLL = 1.d0
return
endif
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/plotgll.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/plotgll.f90 2005-01-09 13:38:40 UTC (rev 8439)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/plotgll.f90 2007-12-07 23:46:32 UTC (rev 8440)
@@ -7,7 +7,7 @@
! Dimitri Komatitsch
! Universite de Pau et des Pays de l'Adour, France
!
-! (c) December 2004
+! (c) January 2005
!
!========================================================================
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/plotpost.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/plotpost.f90 2005-01-09 13:38:40 UTC (rev 8439)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/plotpost.f90 2007-12-07 23:46:32 UTC (rev 8440)
@@ -7,7 +7,7 @@
! Dimitri Komatitsch
! Universite de Pau et des Pays de l'Adour, France
!
-! (c) December 2004
+! (c) January 2005
!
!========================================================================
@@ -16,7 +16,7 @@
Uxinterp,Uzinterp,flagrange,density,elastcoef,knods,kmato,ibool, &
numabs,codeabs,anyabs,stitle,npoin,npgeo,vpmin,vpmax,nrec, &
colors,numbers,subsamp,vecttype,interpol,meshvect,modelvect, &
- boundvect,readmodel,cutvect,nelemabs,numat,iptsdisp,nspec,ngnod)
+ boundvect,readmodel,cutvect,nelemabs,numat,iptsdisp,nspec,ngnod,ELASTIC)
!
! routine affichage postscript
@@ -51,7 +51,7 @@
integer iglob_rec(nrec)
integer numabs(nelemabs),codeabs(4,nelemabs)
- logical anyabs
+ logical anyabs,ELASTIC
double precision xmax,zmax,height,xw,zw,usoffset,sizex,sizez,vpmin,vpmax
@@ -337,7 +337,13 @@
write(24,*) '26.45 CM 18.9 CM MV'
write(24,*) usoffset,' CM 2 div neg 0 MR'
write(24,*) 'currentpoint gsave translate -90 rotate 0 0 moveto'
- write(24,*) '(Elastic Wave 2D - Spectral Element Method) show'
+
+ if(ELASTIC) then
+ write(24,*) '(Elastic Wave 2D - Spectral Element Method) show'
+ else
+ write(24,*) '(Acoustic Wave 2D - Spectral Element Method) show'
+ endif
+
write(24,*) 'grestore'
endif
@@ -372,7 +378,7 @@
rmu = elastcoef(2,material)
denst = density(material)
rKvol = rlamda + 2.d0*rmu/3.d0
- cploc = dsqrt((rKvol + 4.d0*rmu/3.d0)/denst)
+ cploc = sqrt((rKvol + 4.d0*rmu/3.d0)/denst)
x1 = (cploc-vpmin)/(vpmax-vpmin)
endif
else
@@ -712,18 +718,18 @@
x2 = Uxinterp(i,j)*sizemax/dispmax
z2 = Uzinterp(i,j)*sizemax/dispmax
- d = dsqrt(x2**2 + z2**2)
+ d = sqrt(x2**2 + z2**2)
! ignorer si vecteur trop petit
if(d > cutvect*sizemax) then
d1 = d * rapport
- d2 = d1 * dcos(angle*convert)
+ d2 = d1 * cos(angle*convert)
dummy = x2/d
if(dummy > 0.9999d0) dummy = 0.9999d0
if(dummy < -0.9999d0) dummy = -0.9999d0
- theta = dacos(dummy)
+ theta = acos(dummy)
if(z2 < 0.d0) theta = 360.d0*convert - theta
thetaup = theta - angle*convert
@@ -734,12 +740,12 @@
z1 = (orig_z+z1) * centim
x2 = x2 * centim
z2 = z2 * centim
- xa = -d2*dcos(thetaup)
- za = -d2*dsin(thetaup)
+ xa = -d2*cos(thetaup)
+ za = -d2*sin(thetaup)
xa = xa * centim
za = za * centim
- xb = -d2*dcos(thetadown)
- zb = -d2*dsin(thetadown)
+ xb = -d2*cos(thetadown)
+ zb = -d2*sin(thetadown)
xb = xb * centim
zb = zb * centim
write(name,700) xb,zb,xa,za,x2,z2,x1,z1
@@ -777,18 +783,18 @@
x2 = displ(1,ipoin)*sizemax/dispmax
z2 = displ(2,ipoin)*sizemax/dispmax
- d = dsqrt(x2**2 + z2**2)
+ d = sqrt(x2**2 + z2**2)
! ignorer si vecteur trop petit
if(d > cutvect*sizemax) then
d1 = d * rapport
- d2 = d1 * dcos(angle*convert)
+ d2 = d1 * cos(angle*convert)
dummy = x2/d
if(dummy > 0.9999d0) dummy = 0.9999d0
if(dummy < -0.9999d0) dummy = -0.9999d0
- theta = dacos(dummy)
+ theta = acos(dummy)
if(z2 < 0.d0) theta = 360.d0*convert - theta
thetaup = theta - angle*convert
@@ -799,12 +805,12 @@
z1 = (orig_z+z1) * centim
x2 = x2 * centim
z2 = z2 * centim
- xa = -d2*dcos(thetaup)
- za = -d2*dsin(thetaup)
+ xa = -d2*cos(thetaup)
+ za = -d2*sin(thetaup)
xa = xa * centim
za = za * centim
- xb = -d2*dcos(thetadown)
- zb = -d2*dsin(thetadown)
+ xb = -d2*cos(thetadown)
+ zb = -d2*sin(thetadown)
xb = xb * centim
zb = zb * centim
write(name,700) xb,zb,xa,za,x2,z2,x1,z1
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/positrec.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/positrec.f90 2005-01-09 13:38:40 UTC (rev 8439)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/positrec.f90 2007-12-07 23:46:32 UTC (rev 8440)
@@ -7,7 +7,7 @@
! Dimitri Komatitsch
! Universite de Pau et des Pays de l'Adour, France
!
-! (c) December 2004
+! (c) January 2005
!
!========================================================================
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/positsource.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/positsource.f90 2005-01-09 13:38:40 UTC (rev 8439)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/positsource.f90 2007-12-07 23:46:32 UTC (rev 8440)
@@ -7,7 +7,7 @@
! Dimitri Komatitsch
! Universite de Pau et des Pays de l'Adour, France
!
-! (c) December 2004
+! (c) January 2005
!
!========================================================================
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/q49shape.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/q49shape.f90 2005-01-09 13:38:40 UTC (rev 8439)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/q49shape.f90 2007-12-07 23:46:32 UTC (rev 8440)
@@ -7,7 +7,7 @@
! Dimitri Komatitsch
! Universite de Pau et des Pays de l'Adour, France
!
-! (c) December 2004
+! (c) January 2005
!
!========================================================================
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/q49spec.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/q49spec.f90 2005-01-09 13:38:40 UTC (rev 8439)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/q49spec.f90 2007-12-07 23:46:32 UTC (rev 8440)
@@ -7,7 +7,7 @@
! Dimitri Komatitsch
! Universite de Pau et des Pays de l'Adour, France
!
-! (c) December 2004
+! (c) January 2005
!
!========================================================================
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/specfem2D.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/specfem2D.f90 2005-01-09 13:38:40 UTC (rev 8439)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/specfem2D.f90 2007-12-07 23:46:32 UTC (rev 8440)
@@ -7,7 +7,7 @@
! Dimitri Komatitsch
! Universite de Pau et des Pays de l'Adour, France
!
-! (c) December 2004
+! (c) January 2005
!
!========================================================================
@@ -18,11 +18,14 @@
!====================================================================================
!
-! version 5.1, December 2004 :
+! version 5.1, January 2005 :
! - Dirac and Gaussian time sources and corresponding convolution routine
! - more general mesher with any number of curved layers
+! - option for acoustic medium instead of elastic
! - color PNM snapshots
! - more flexible Par file with any number of comment lines
+! - Xsu scripts for seismograms
+! - subtract t0 from seismograms
!
! version 5.0, May 2004 :
! - got rid of useless routines, suppressed commons etc.
@@ -65,7 +68,7 @@
double precision valux,valuz,rhoextread,vpextread,vsextread
double precision cpl,csl,rhol
- double precision dcosrot,dsinrot,xcor,zcor
+ double precision cosrot,sinrot,xcor,zcor
! coefficients of the explicit Newmark time scheme
integer NSTEP
@@ -96,13 +99,13 @@
double precision xixl,xizl,gammaxl,gammazl,jacobianl
! material properties of the elastic medium
- double precision mul_relaxed,lambdal_relaxed,lambdalplus2mul_relaxed
+ double precision mul_relaxed,lambdal_relaxed,lambdalplus2mul_relaxed,cpsquare
double precision mul_unrelaxed,lambdal_unrelaxed,lambdalplus2mul_unrelaxed
double precision, dimension(:), allocatable :: xirec,etarec
double precision, dimension(:,:), allocatable :: coord,accel,veloc,displ, &
- flagrange,xinterp,zinterp,Uxinterp,Uzinterp,elastcoef
+ flagrange,xinterp,zinterp,Uxinterp,Uzinterp,elastcoef,vector_field_postscript
double precision, dimension(:), allocatable :: rmass, &
fglobx,fglobz,density,vpext,vsext,rhoext,displread,velocread,accelread
@@ -116,7 +119,7 @@
integer, dimension(:,:,:), allocatable :: ibool
integer, dimension(:,:), allocatable :: knods
- integer, dimension(:), allocatable :: kmato,numabs
+ integer, dimension(:), allocatable :: kmato,numabs,numsurface
integer ie,k
@@ -127,17 +130,17 @@
lambdal_Smin,lambdal_Smax,lambdal_Pmin,lambdal_Pmax,vpmin,vpmax
integer colors,numbers,subsamp,vecttype,itaff,nrec,sismostype
- integer numat,ngnod,nspec,iptsdisp,nelemabs
+ integer numat,ngnod,nspec,iptsdisp,nelemabs,nelemsurface
- logical interpol,meshvect,modelvect,boundvect,readmodel,initialfield, &
- outputgrid,gnuplot,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON
+ logical interpol,meshvect,modelvect,boundvect,readmodel,initialfield,abshaut, &
+ outputgrid,gnuplot,ELASTIC,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON
double precision cutvect,anglerec
-! for absorbing conditions
- integer ispecabs,inum,numabsread,i1abs,i2abs
+! for absorbing and free surface conditions
+ integer ispecabs,ispecsurface,inum,numabsread,numsurfaceread,i1abs,i2abs
logical codeabsread(4)
- double precision nx,nz,vx,vz,vn,rho_vp,rho_vs,tx,tz,weight,xxi,zeta,rKvol
+ double precision nx,nz,vx,vz,vn,rho_vp,rho_vs,tx,tz,weight,xxi,zeta,kappal
logical, dimension(:,:), allocatable :: codeabs
@@ -201,7 +204,7 @@
write(*,*)
write(*,*) '*********************************'
write(*,*) '**** ****'
- write(*,*) '**** SPECFEM2D VERSION 5.0 ****'
+ write(*,*) '**** SPECFEM2D VERSION 5.1 ****'
write(*,*) '**** ****'
write(*,*) '*********************************'
@@ -232,13 +235,13 @@
read(IIN,*) sismostype,vecttype
read(IIN,40) datlin
- read(IIN,*) readmodel,outputgrid,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON
+ read(IIN,*) readmodel,outputgrid,ELASTIC,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON
!---- check parameters read
write(IOUT,200) npgeo,NDIME
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,750) initialfield,readmodel,ELASTIC,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON,outputgrid
write(IOUT,800) vecttype,100.d0*cutvect,subsamp
!---- read time step
@@ -288,7 +291,7 @@
read(IIN,40) datlin
allocate(posrecread(NDIME))
do i=1,nrec
- read(IIN ,*) irec,(posrecread(j),j=1,NDIME)
+ read(IIN,*) irec,(posrecread(j),j=1,NDIME)
if(irec<1 .or. irec>nrec) stop 'Wrong receiver number'
posrec(:,irec) = posrecread
enddo
@@ -310,8 +313,8 @@
!
!---- read the basic properties of the spectral elements
!
- read(IIN ,40) datlin
- read(IIN ,*) numat,ngnod,nspec,iptsdisp,nelemabs
+ read(IIN,40) datlin
+ read(IIN,*) numat,ngnod,nspec,iptsdisp,nelemabs,nelemsurface
!
!---- allocate arrays
@@ -339,6 +342,13 @@
allocate(knods(ngnod,nspec))
allocate(ibool(NGLLX,NGLLZ,nspec))
+! for acoustic
+ if(TURN_ANISOTROPY_ON .and. .not. ELASTIC) stop 'currently cannot have anisotropy in acoustic simulation'
+
+ if(TURN_ATTENUATION_ON .and. .not. ELASTIC) stop 'currently cannot have attenuation in acoustic simulation'
+
+ if(source_type == 2 .and. .not. ELASTIC) stop 'currently cannot have moment tensor source in acoustic simulation'
+
! for attenuation
if(TURN_ANISOTROPY_ON .and. TURN_ATTENUATION_ON) stop 'cannot have anisotropy and attenuation both turned on in current version'
@@ -373,6 +383,10 @@
allocate(numabs(nelemabs))
allocate(codeabs(4,nelemabs))
+! --- allocate array for free surface condition in acoustic medium
+ if(nelemsurface <= 0) nelemsurface = 1
+ allocate(numsurface(nelemsurface))
+
!
!---- print element group main parameters
!
@@ -400,9 +414,9 @@
!---- read absorbing boundary data
!
if(anyabs) then
- read(IIN ,40) datlin
+ read(IIN,40) datlin
do n=1,nelemabs
- read(IIN ,*) inum,numabsread,codeabsread(1),codeabsread(2),codeabsread(3),codeabsread(4)
+ read(IIN,*) inum,numabsread,codeabsread(1),codeabsread(2),codeabsread(3),codeabsread(4)
if(inum < 1 .or. inum > nelemabs) stop 'Wrong absorbing element number'
numabs(inum) = numabsread
codeabs(ITOP,inum) = codeabsread(1)
@@ -415,6 +429,19 @@
endif
!
+!---- read free surface data
+!
+ read(IIN,40) datlin
+ read(IIN,*) abshaut
+ do n=1,nelemsurface
+ read(IIN,*) inum,numsurfaceread
+ if(inum < 1 .or. inum > nelemsurface) stop 'Wrong free surface element number'
+ numsurface(inum) = numsurfaceread
+ enddo
+ write(*,*)
+ write(*,*) 'Number of free surface elements: ',nelemsurface
+
+!
!---- compute the spectral element shape functions and their local derivatives
!
call q49shape(shape,dershape,xigll,yigll,ngnod)
@@ -452,6 +479,13 @@
allocate(displ(NDIME,npoin))
allocate(veloc(NDIME,npoin))
+! for acoustic medium
+ if(ELASTIC) then
+ allocate(vector_field_postscript(NDIME,1))
+ else
+ allocate(vector_field_postscript(NDIME,npoin))
+ endif
+
allocate(rmass(npoin))
allocate(fglobx(npoin))
@@ -558,13 +592,22 @@
do j = 1,NGLLZ
do i = 1,NGLLX
iglob = ibool(i,j,ispec)
-!--- if external density model
+! if external density model
if(readmodel) then
rhol = rhoext(iglob)
+ cpsquare = vpext(iglob)**2
else
rhol = density(kmato(ispec))
+ lambdal_relaxed = elastcoef(1,kmato(ispec))
+ mul_relaxed = elastcoef(2,kmato(ispec))
+ cpsquare = (lambdal_relaxed + 2.d0*mul_relaxed) / rhol
endif
- rmass(iglob) = rmass(iglob) + wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec)
+! for acoustic medium
+ if(ELASTIC) then
+ rmass(iglob) = rmass(iglob) + wxgll(i)*wzgll(j)*rhol*jacobian(i,j,ispec)
+ else
+ rmass(iglob) = rmass(iglob) + wxgll(i)*wzgll(j)*jacobian(i,j,ispec) / cpsquare
+ endif
enddo
enddo
enddo
@@ -576,7 +619,7 @@
!---- verifier le maillage, la stabilite et le nb de points par lambda
!---- seulement si la source en temps n'est pas un Dirac (sinon spectre non defini)
!
- if(time_function_type /= 3) call checkgrid(deltat,f0,t0,initialfield, &
+ if(time_function_type /= 4) call checkgrid(deltat,f0,t0,initialfield, &
rsizemin,rsizemax,cpoverdxmin,cpoverdxmax,lambdal_Smin,lambdal_Smax,lambdal_Pmin,lambdal_Pmax)
!
@@ -700,8 +743,8 @@
sisux = ZERO
sisuz = ZERO
- dcosrot = dcos(anglerec)
- dsinrot = dsin(anglerec)
+ cosrot = cos(anglerec)
+ sinrot = sin(anglerec)
! initialiser les tableaux a zero
accel = ZERO
@@ -799,6 +842,31 @@
veloc(:,:) = veloc(:,:) + deltatover2*accel(:,:)
accel(:,:) = ZERO
+
+!--- free surface for an acoustic medium
+
+! if acoustic, the free surface condition is a Dirichlet condition for the potential,
+! not Neumann, in order to impose zero pressure at the surface. Also check that
+! top absorbing boundary is not set because cannot be both absorbing and free surface
+ if(.not. ELASTIC .and. .not. abshaut) then
+
+ do ispecsurface=1,nelemsurface
+
+ ispec = numsurface(ispecsurface)
+
+ j = NGLLZ
+ do i=1,NGLLX
+ iglob = ibool(i,j,ispec)
+ displ(:,iglob) = ZERO
+ veloc(:,iglob) = ZERO
+ accel(:,iglob) = ZERO
+ enddo
+
+ enddo
+
+ endif ! end of free surface condition for acoustic medium
+
+
! integration over spectral elements
do ispec = 1,NSPEC
@@ -935,6 +1003,12 @@
tempx2(i,j) = jacobianl*(sigma_xx*gammaxl+sigma_zx*gammazl)
tempz2(i,j) = jacobianl*(sigma_xz*gammaxl+sigma_zz*gammazl)
+! for acoustic medium
+ if(.not. ELASTIC) then
+ tempx1(i,j) = jacobianl*(xixl*dUxdxl + xizl*dUxdzl)
+ tempx2(i,j) = jacobianl*(gammaxl*dUxdxl + gammazl*dUxdzl)
+ endif
+
enddo
enddo
@@ -950,7 +1024,7 @@
do k = 1,NGLLX
fac1 = wxgll(k)*hprime_xx(i,k)
tempx1l = tempx1l + tempx1(k,j)*fac1
- tempz1l = tempz1l + tempz1(k,j)*fac1
+ if(ELASTIC) tempz1l = tempz1l + tempz1(k,j)*fac1
enddo
! along z direction
@@ -959,15 +1033,21 @@
do k = 1,NGLLZ
fac2 = wzgll(k)*hprime_zz(j,k)
tempx2l = tempx2l + tempx2(i,k)*fac2
- tempz2l = tempz2l + tempz2(i,k)*fac2
+ if(ELASTIC) tempz2l = tempz2l + tempz2(i,k)*fac2
enddo
+! GLL integration weights
fac1 = wzgll(j)
fac2 = wxgll(i)
+! for acoustic medium
iglob = ibool(i,j,ispec)
accel(1,iglob) = accel(1,iglob) - (fac1*tempx1l + fac2*tempx2l)
- accel(2,iglob) = accel(2,iglob) - (fac1*tempz1l + fac2*tempz2l)
+ if(ELASTIC) then
+ accel(2,iglob) = accel(2,iglob) - (fac1*tempz1l + fac2*tempz2l)
+ else
+ accel(2,iglob) = zero
+ endif
enddo ! second loop over the GLL points
enddo
@@ -987,9 +1067,9 @@
lambdal_relaxed = elastcoef(1,kmato(ispec))
mul_relaxed = elastcoef(2,kmato(ispec))
rhol = density(kmato(ispec))
- rKvol = lambdal_relaxed + TWO*mul_relaxed/3.d0
- cpl = dsqrt((rKvol + 4.d0*mul_relaxed/3.d0)/rhol)
- csl = dsqrt(mul_relaxed/rhol)
+ kappal = lambdal_relaxed + TWO*mul_relaxed/3.d0
+ cpl = sqrt((kappal + 4.d0*mul_relaxed/3.d0)/rhol)
+ csl = sqrt(mul_relaxed/rhol)
!--- left absorbing boundary
@@ -1026,8 +1106,13 @@
weight = zeta*wzgll(j)
- accel(1,iglob) = accel(1,iglob) - tx*weight
- accel(2,iglob) = accel(2,iglob) - tz*weight
+! Clayton-Engquist condition if elastic, Sommerfeld condition if acoustic
+ if(ELASTIC) then
+ accel(1,iglob) = accel(1,iglob) - tx*weight
+ accel(2,iglob) = accel(2,iglob) - tz*weight
+ else
+ accel(1,iglob) = accel(1,iglob) - veloc(1,iglob)*weight/cpl
+ endif
enddo
@@ -1067,8 +1152,13 @@
weight = zeta*wzgll(j)
- accel(1,iglob) = accel(1,iglob) - tx*weight
- accel(2,iglob) = accel(2,iglob) - tz*weight
+! Clayton-Engquist condition if elastic, Sommerfeld condition if acoustic
+ if(ELASTIC) then
+ accel(1,iglob) = accel(1,iglob) - tx*weight
+ accel(2,iglob) = accel(2,iglob) - tz*weight
+ else
+ accel(1,iglob) = accel(1,iglob) - veloc(1,iglob)*weight/cpl
+ endif
enddo
@@ -1114,8 +1204,13 @@
weight = xxi*wxgll(i)
- accel(1,iglob) = accel(1,iglob) - tx*weight
- accel(2,iglob) = accel(2,iglob) - tz*weight
+! Clayton-Engquist condition if elastic, Sommerfeld condition if acoustic
+ if(ELASTIC) then
+ accel(1,iglob) = accel(1,iglob) - tx*weight
+ accel(2,iglob) = accel(2,iglob) - tz*weight
+ else
+ accel(1,iglob) = accel(1,iglob) - veloc(1,iglob)*weight/cpl
+ endif
enddo
@@ -1161,8 +1256,13 @@
weight = xxi*wxgll(i)
- accel(1,iglob) = accel(1,iglob) - tx*weight
- accel(2,iglob) = accel(2,iglob) - tz*weight
+! Clayton-Engquist condition if elastic, Sommerfeld condition if acoustic
+ if(ELASTIC) then
+ accel(1,iglob) = accel(1,iglob) - tx*weight
+ accel(2,iglob) = accel(2,iglob) - tz*weight
+ else
+ accel(1,iglob) = accel(1,iglob) - veloc(1,iglob)*weight/cpl
+ endif
enddo
@@ -1176,24 +1276,34 @@
! --- add the source
if(.not. initialfield) then
-! Ricker source time function
+! Ricker (second derivative of a Gaussian) source time function
if(time_function_type == 1) then
source_time_function = - factor * (ONE-TWO*a*(time-t0)**2) * exp(-a*(time-t0)**2)
+! first derivative of a Gaussian source time function
+ else if(time_function_type == 2) then
+ source_time_function = - factor * TWO*a*(time-t0) * 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
+ else if(time_function_type == 3 .or. time_function_type == 4) then
source_time_function = factor * exp(-a*(time-t0)**2)
else
stop 'unknown source time function'
endif
-! --- collocated force
+! collocated force
+! beware, for acoustic medium, source is a potential, therefore source time function
+! gives shape of velocity, not displacement
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
+ if(ELASTIC) then
+ accel(1,iglob_source) = accel(1,iglob_source) - sin(angleforce)*source_time_function
+ accel(2,iglob_source) = accel(2,iglob_source) + cos(angleforce)*source_time_function
+ else
+ accel(1,iglob_source) = accel(1,iglob_source) + source_time_function
+ endif
-!---- explosion
+! explosion
else if(source_type == 2) then
do i=1,NGLLX
do j=1,NGLLX
@@ -1215,6 +1325,31 @@
! update velocity
veloc(:,:) = veloc(:,:) + deltatover2*accel(:,:)
+
+!--- free surface for an acoustic medium
+
+! if acoustic, the free surface condition is a Dirichlet condition for the potential,
+! not Neumann, in order to impose zero pressure at the surface. Also check that
+! top absorbing boundary is not set because cannot be both absorbing and free surface
+ if(.not. ELASTIC .and. .not. abshaut) then
+
+ do ispecsurface=1,nelemsurface
+
+ ispec = numsurface(ispecsurface)
+
+ j = NGLLZ
+ do i=1,NGLLX
+ iglob = ibool(i,j,ispec)
+ displ(:,iglob) = ZERO
+ veloc(:,iglob) = ZERO
+ accel(:,iglob) = ZERO
+ enddo
+
+ enddo
+
+ endif ! end of free surface condition for acoustic medium
+
+
! implement attenuation
if(TURN_ATTENUATION_ON) then
@@ -1323,30 +1458,54 @@
!---- display max of norm of displacement
if(mod(it,itaff) == 0) then
displnorm_all = maxval(sqrt(displ(1,:)**2 + displ(2,:)**2))
- print *,'Max norm of displacement = ',displnorm_all
+ print *,'Max norm of field = ',displnorm_all
! check stability of the code, exit if unstable
if(displnorm_all > STABILITY_THRESHOLD) stop 'code became unstable and blew up'
endif
! store the seismograms
- do irec=1,nrec
+ if(sismostype < 1 .or. sismostype > 3) stop 'Wrong field code for seismogram output'
+ if(.not. ELASTIC) then
if(sismostype == 1) then
- valux = displ(1,iglob_rec(irec))
- valuz = displ(2,iglob_rec(irec))
+ stop 'cannot store displacement field in acoustic medium because of potential formulation'
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))
+! for acoustic medium, compute gradient for display, displ represents the potential
+ call compute_gradient_fluid(displ,vector_field_postscript, &
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,NSPEC,npoin)
else
- stop 'Wrong field code for seismogram output'
+! for acoustic medium, compute gradient for display, veloc represents the first derivative of the potential
+ call compute_gradient_fluid(veloc,vector_field_postscript, &
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,NSPEC,npoin)
endif
+ endif
+ do irec=1,nrec
+
+ if(ELASTIC) then
+
+ 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
+ valux = accel(1,iglob_rec(irec))
+ valuz = accel(2,iglob_rec(irec))
+ endif
+
+ else
+
+! for acoustic medium
+ valux = vector_field_postscript(1,iglob_rec(irec))
+ valuz = vector_field_postscript(2,iglob_rec(irec))
+
+ endif
+
! rotation eventuelle des composantes
- sisux(it,irec) = dcosrot*valux + dsinrot*valuz
- sisuz(it,irec) = - dsinrot*valux + dcosrot*valuz
+ sisux(it,irec) = cosrot*valux + sinrot*valuz
+ sisuz(it,irec) = - sinrot*valux + cosrot*valuz
enddo
@@ -1367,30 +1526,63 @@
!---- affichage postscript
!
write(IOUT,*) 'Dump PostScript'
- if(vecttype == 1) then
+
+! for elastic medium
+ if(ELASTIC .and. vecttype == 1) then
write(IOUT,*) 'drawing displacement field...'
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, &
colors,numbers,subsamp,vecttype,interpol,meshvect,modelvect, &
- boundvect,readmodel,cutvect,nelemabs,numat,iptsdisp,nspec,ngnod)
- else if(vecttype == 2) then
+ boundvect,readmodel,cutvect,nelemabs,numat,iptsdisp,nspec,ngnod,ELASTIC)
+
+ else if(ELASTIC .and. vecttype == 2) then
write(IOUT,*) 'drawing velocity field...'
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, &
colors,numbers,subsamp,vecttype,interpol,meshvect,modelvect, &
- boundvect,readmodel,cutvect,nelemabs,numat,iptsdisp,nspec,ngnod)
- else if(vecttype == 3) then
+ boundvect,readmodel,cutvect,nelemabs,numat,iptsdisp,nspec,ngnod,ELASTIC)
+
+ else if(ELASTIC .and. vecttype == 3) then
write(IOUT,*) 'drawing acceleration field...'
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, &
colors,numbers,subsamp,vecttype,interpol,meshvect,modelvect, &
- boundvect,readmodel,cutvect,nelemabs,numat,iptsdisp,nspec,ngnod)
+ boundvect,readmodel,cutvect,nelemabs,numat,iptsdisp,nspec,ngnod,ELASTIC)
+
+! for acoustic medium
+ else if(.not. ELASTIC .and. vecttype == 1) then
+ stop 'cannot display displacement field in acoustic medium because of potential formulation'
+
+ else if(.not. ELASTIC .and. vecttype == 2) then
+ write(IOUT,*) 'drawing acoustic velocity field from velocity potential...'
+! for acoustic medium, compute gradient for display, displ represents the potential
+ call compute_gradient_fluid(displ,vector_field_postscript, &
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,NSPEC,npoin)
+ call plotpost(vector_field_postscript,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, &
+ colors,numbers,subsamp,vecttype,interpol,meshvect,modelvect, &
+ boundvect,readmodel,cutvect,nelemabs,numat,iptsdisp,nspec,ngnod,ELASTIC)
+
+ else if(.not. ELASTIC .and. vecttype == 3) then
+ write(IOUT,*) 'drawing acoustic acceleration field from velocity potential...'
+! for acoustic medium, compute gradient for display, veloc represents the first derivative of the potential
+ call compute_gradient_fluid(veloc,vector_field_postscript, &
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,NSPEC,npoin)
+ call plotpost(vector_field_postscript,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, &
+ colors,numbers,subsamp,vecttype,interpol,meshvect,modelvect, &
+ boundvect,readmodel,cutvect,nelemabs,numat,iptsdisp,nspec,ngnod,ELASTIC)
+
else
stop 'wrong field code for PostScript display'
endif
@@ -1405,8 +1597,21 @@
do j = 1,NZ_IMAGE_PNM
do i = 1,NX_IMAGE_PNM
-! afficher la composante verticale du deplacement
- if(iglob_image_PNM_2D(i,j) /= -1) donnees_image_PNM_2D(i,j) = displ(2,iglob_image_PNM_2D(i,j))
+ if(iglob_image_PNM_2D(i,j) /= -1) then
+! display vertical component of vector
+ if(ELASTIC) then
+ if(vecttype == 1) then
+ donnees_image_PNM_2D(i,j) = displ(2,iglob_image_PNM_2D(i,j))
+ else if(vecttype == 2) then
+ donnees_image_PNM_2D(i,j) = veloc(2,iglob_image_PNM_2D(i,j))
+ else
+ donnees_image_PNM_2D(i,j) = accel(2,iglob_image_PNM_2D(i,j))
+ endif
+ else
+! for acoustic medium
+ donnees_image_PNM_2D(i,j) = vector_field_postscript(2,iglob_image_PNM_2D(i,j))
+ endif
+ endif
enddo
enddo
@@ -1415,14 +1620,14 @@
write(IOUT,*) 'Fin creation image PNM'
!---- save temporary seismograms
- call write_seismograms(sisux,sisuz,NSTEP,nrec,deltat,sismostype,iglob_rec,coord,npoin)
+ call write_seismograms(sisux,sisuz,NSTEP,nrec,deltat,sismostype,iglob_rec,coord,npoin,t0)
endif
enddo ! end of the main time loop
!---- save final seismograms
- call write_seismograms(sisux,sisuz,NSTEP,nrec,deltat,sismostype,iglob_rec,coord,npoin)
+ call write_seismograms(sisux,sisuz,NSTEP,nrec,deltat,sismostype,iglob_rec,coord,npoin,t0)
! print exit banner
call datim(stitle)
@@ -1460,10 +1665,11 @@
'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. . .(readmodel) = ',l6/5x, &
+ 'Read external velocity model or not . . .(readmodel) = ',l6/5x, &
+ 'Elastic simulation or acoustic. . . . . . .(ELASTIC) = ',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 . . .(outputgrid) = ',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 . . . . . . . . . . .(vecttype) = ',i6/5x, &
'Percentage of cut for vector plots. . . . .(cutvect) = ',f6.2/5x, &
Modified: seismo/2D/SPECFEM2D/trunk/SPECFEM90/write_seismograms.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/SPECFEM90/write_seismograms.f90 2005-01-09 13:38:40 UTC (rev 8439)
+++ seismo/2D/SPECFEM2D/trunk/SPECFEM90/write_seismograms.f90 2007-12-07 23:46:32 UTC (rev 8440)
@@ -7,11 +7,11 @@
! Dimitri Komatitsch
! Universite de Pau et des Pays de l'Adour, France
!
-! (c) December 2004
+! (c) January 2005
!
!========================================================================
- subroutine write_seismograms(sisux,sisuz,nt,nrec,deltat,sismostype,iglob_rec,coord,npoin)
+ subroutine write_seismograms(sisux,sisuz,nt,nrec,deltat,sismostype,iglob_rec,coord,npoin,t0)
! save the seismograms in ASCII format
@@ -20,7 +20,7 @@
include "constants.h"
integer nt,nrec,sismostype,npoin
- double precision deltat
+ double precision deltat,t0
integer iglob_rec(nrec)
@@ -42,7 +42,7 @@
write(name,221) irec
open(unit=11,file=name,status='unknown')
do it=1,nt
- write(11,*) sngl(dble(it-1)*deltat),' ',sngl(sisux(it,irec))
+ write(11,*) sngl(dble(it-1)*deltat - t0),' ',sngl(sisux(it,irec))
enddo
close(11)
enddo
@@ -52,7 +52,7 @@
write(name,222) irec
open(unit=11,file=name,status='unknown')
do it=1,nt
- write(11,*) sngl(dble(it-1)*deltat),' ',sngl(sisuz(it,irec))
+ write(11,*) sngl(dble(it-1)*deltat - t0),' ',sngl(sisuz(it,irec))
enddo
close(11)
enddo
@@ -62,13 +62,17 @@
! 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)
+ open(unit=11,file='Ux_file.bin',status='unknown',access='direct',recl=nt*4)
+ do irec=1,nrec
+ write(11,rec=irec) (sngl(sisux(it,irec)),it=1,nt)
+ enddo
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)
+ open(unit=11,file='Uz_file.bin',status='unknown',access='direct',recl=nt*4)
+ do irec=1,nrec
+ write(11,rec=irec) (sngl(sisuz(it,irec)),it=1,nt)
+ enddo
close(11)
!----
@@ -76,7 +80,8 @@
! ligne de recepteurs pour Xsu
open(unit=11,file='receiver_line_Xsu_XWindow',status='unknown')
- write(11,110) FACTORXSU,nt,deltat,nrec
+! subtract t0 from seismograms to get correct zero time
+ write(11,110) FACTORXSU,nt,deltat,-t0,nrec
do irec=1,nrec
write(11,140) coord(1,iglob_rec(irec))
@@ -120,7 +125,7 @@
! 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=')
+ 110 format('xwigb at xcur=',f8.2,'@n1=',i5,'@d1=',f15.8,'@f1=',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')
More information about the cig-commits
mailing list