[cig-commits] r8461 - seismo/2D/SPECFEM2D/trunk
walter at geodynamics.org
walter at geodynamics.org
Fri Dec 7 15:47:57 PST 2007
Author: walter
Date: 2007-12-07 15:47:57 -0800 (Fri, 07 Dec 2007)
New Revision: 8461
Modified:
seismo/2D/SPECFEM2D/trunk/Makefile
seismo/2D/SPECFEM2D/trunk/defarrays.f90
seismo/2D/SPECFEM2D/trunk/meshfem2D.f90
seismo/2D/SPECFEM2D/trunk/plotpost.f90
seismo/2D/SPECFEM2D/trunk/specfem2D.f90
Log:
improved handling of external velocity model in 2D code
Modified: seismo/2D/SPECFEM2D/trunk/Makefile
===================================================================
--- seismo/2D/SPECFEM2D/trunk/Makefile 2005-02-10 17:31:34 UTC (rev 8460)
+++ seismo/2D/SPECFEM2D/trunk/Makefile 2007-12-07 23:47:57 UTC (rev 8461)
@@ -46,7 +46,7 @@
all: default
clean:
- /bin/rm -r -f xmeshfem2D xmeshfem2D.trace xspecfem2D xspecfem2D.trace $O/*.o *.o $O/*.il *.mod core *.gnu *.ps Ux*.bin Uz*.bin image*.pnm xconvolve_source_timefunction *receiver_line_* plotgnu source.txt *.sem* OUTPUT_FILES/*
+ /bin/rm -r -f xmeshfem2D xmeshfem2D.trace xspecfem2D xspecfem2D.trace $O/*.o *.o $O/*.il *.mod core *.gnu *.ps Ux*.bin Uz*.bin image*.pnm xconvolve_source_timefunction *receiver_line_* plotgnu source.txt *.sem* xcreate_earth_model
meshfem2D: $(OBJS_MESHFEM2D)
$(LINK) $(FLAGS_CHECK) -o xmeshfem2D $(OBJS_MESHFEM2D)
@@ -58,6 +58,9 @@
convolve_source_timefunction: $O/convolve_source_timefunction.o
${F90} $(FLAGS_CHECK) -o xconvolve_source_timefunction $O/convolve_source_timefunction.o
+create_earth_model: $O/create_earth_model.o
+ ${F90} $(FLAGS_CHECK) -o xcreate_earth_model $O/create_earth_model.o
+
$O/checkgrid.o: checkgrid.f90 constants.h
${F90} $(FLAGS_CHECK) -c -o $O/checkgrid.o checkgrid.f90
@@ -73,6 +76,9 @@
$O/convolve_source_timefunction.o: convolve_source_timefunction.f90
${F90} $(FLAGS_CHECK) -c -o $O/convolve_source_timefunction.o convolve_source_timefunction.f90
+$O/create_earth_model.o: create_earth_model.f90
+ ${F90} $(FLAGS_CHECK) -c -o $O/create_earth_model.o create_earth_model.f90
+
$O/datim.o: datim.f90 constants.h
${F90} $(FLAGS_CHECK) -c -o $O/datim.o datim.f90
Modified: seismo/2D/SPECFEM2D/trunk/defarrays.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/defarrays.f90 2005-02-10 17:31:34 UTC (rev 8460)
+++ seismo/2D/SPECFEM2D/trunk/defarrays.f90 2007-12-07 23:47:57 UTC (rev 8461)
@@ -14,7 +14,7 @@
subroutine defarrays(vpext,vsext,rhoext,density,elastcoef, &
ibool,kmato,coord,npoin,rsizemin,rsizemax, &
cpoverdxmin,cpoverdxmax,lambdaSmin,lambdaSmax,lambdaPmin,lambdaPmax, &
- vpmin,vpmax,readmodel,nspec,numat)
+ vpmin,vpmax,read_external_model,nspec,numat)
! define all the arrays for the variational formulation
@@ -43,7 +43,7 @@
double precision rsizemin,rsizemax,cpoverdxmin,cpoverdxmax, &
lambdaSmin,lambdaSmax,lambdaPmin,lambdaPmax,vpmin,vpmax
- logical readmodel
+ logical read_external_model
!
!-----------------------------------------------------------------------
@@ -87,7 +87,7 @@
do i=1,NGLLX
!--- si formulation heterogene pour un modele de vitesse externe
- if(readmodel) then
+ if(read_external_model) then
ipointnum = ibool(i,j,ispec)
cploc = vpext(ipointnum)
csloc = vsext(ipointnum)
Modified: seismo/2D/SPECFEM2D/trunk/meshfem2D.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/meshfem2D.f90 2005-02-10 17:31:34 UTC (rev 8460)
+++ seismo/2D/SPECFEM2D/trunk/meshfem2D.f90 2007-12-07 23:47:57 UTC (rev 8461)
@@ -68,7 +68,7 @@
double precision anglerec,xfin,zfin,xdeb,zdeb,xmin,xmax,dt
double precision rhoread,cpread,csread,aniso3read,aniso4read
- logical interpol,gnuplot,readmodel,outputgrid
+ logical interpol,gnuplot,read_external_model,outputgrid
logical abshaut,absbas,absgauche,absdroite
logical source_surf,enreg_surf,meshvect,initialfield,modelvect,boundvect
logical ELASTIC,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON
@@ -109,7 +109,7 @@
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,read_external_model)
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)
@@ -620,8 +620,8 @@
write(15,*) 'sismostype vecttype'
write(15,*) sismostype,vecttype
- write(15,*) 'readmodel outputgrid ELASTIC TURN_ANISOTROPY_ON TURN_ATTENUATION_ON'
- write(15,*) readmodel,outputgrid,ELASTIC,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON
+ write(15,*) 'read_external_model outputgrid ELASTIC TURN_ANISOTROPY_ON TURN_ATTENUATION_ON'
+ write(15,*) read_external_model,outputgrid,ELASTIC,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON
write(15,*) 'ncycl dtinc'
write(15,*) nt,dt
Modified: seismo/2D/SPECFEM2D/trunk/plotpost.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/plotpost.f90 2005-02-10 17:31:34 UTC (rev 8460)
+++ seismo/2D/SPECFEM2D/trunk/plotpost.f90 2007-12-07 23:47:57 UTC (rev 8461)
@@ -15,7 +15,7 @@
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)
+ boundvect,read_external_model,cutvect,nelemabs,numat,iptsdisp,nspec,ngnod,ELASTIC)
!
! routine affichage postscript
@@ -66,7 +66,7 @@
integer indice,ii,ipoin,in,nnum,ispecabs,ideb,ifin,ibord
integer colors,numbers,subsamp,vecttype
- logical interpol,meshvect,modelvect,boundvect,readmodel
+ logical interpol,meshvect,modelvect,boundvect,read_external_model
double precision cutvect
double precision rapp_page,dispmax,xmin,zmin
@@ -369,7 +369,7 @@
do j=1,NGLLX-subsamp,subsamp
if((vpmax-vpmin)/vpmin > 0.02d0) then
- if(readmodel) then
+ if(read_external_model) then
x1 = (vpext(ibool(i,j,ispec))-vpmin)/ (vpmax-vpmin)
else
material = kmato(ispec)
Modified: seismo/2D/SPECFEM2D/trunk/specfem2D.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/specfem2D.f90 2005-02-10 17:31:34 UTC (rev 8460)
+++ seismo/2D/SPECFEM2D/trunk/specfem2D.f90 2007-12-07 23:47:57 UTC (rev 8461)
@@ -70,7 +70,7 @@
integer i,j,it,irec,ipoin,ip,id
integer nbpoin,inump,n,npoinext,ispec,npoin,npgeo,iglob
- double precision dxd,dzd,valux,valuz,hlagrange,rhoextread,vpextread,vsextread
+ double precision dxd,dzd,valux,valuz,hlagrange,xdummy,zdummy
double precision cpl,csl,rhol
double precision cosrot,sinrot
double precision xi,gamma,x,z
@@ -132,7 +132,7 @@
integer colors,numbers,subsamp,vecttype,itaff,nrec,sismostype
integer numat,ngnod,nspec,iptsdisp,nelemabs,nelemsurface
- logical interpol,meshvect,modelvect,boundvect,readmodel,initialfield,abshaut, &
+ logical interpol,meshvect,modelvect,boundvect,read_external_model,initialfield,abshaut, &
outputgrid,gnuplot,ELASTIC,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON
double precision cutvect,anglerec,xirec,gammarec
@@ -246,13 +246,13 @@
read(IIN,*) sismostype,vecttype
read(IIN,40) datlin
- read(IIN,*) readmodel,outputgrid,ELASTIC,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON
+ read(IIN,*) read_external_model,outputgrid,ELASTIC,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON
!---- check parameters read
write(IOUT,200) npgeo,NDIM
write(IOUT,600) itaff,colors,numbers
write(IOUT,700) sismostype,anglerec
- write(IOUT,750) initialfield,readmodel,ELASTIC,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON,outputgrid
+ write(IOUT,750) initialfield,read_external_model,ELASTIC,TURN_ANISOTROPY_ON,TURN_ATTENUATION_ON,outputgrid
write(IOUT,800) vecttype,100.d0*cutvect,subsamp
!---- read time step
@@ -329,6 +329,8 @@
allocate(knods(ngnod,nspec))
allocate(ibool(NGLLX,NGLLZ,nspec))
+ if(read_external_model .and. outputgrid) stop 'cannot output the grid and read external model at the same time'
+
! for acoustic
if(TURN_ANISOTROPY_ON .and. .not. ELASTIC) stop 'currently cannot have anisotropy in acoustic simulation'
@@ -525,7 +527,7 @@
allocate(rmass(npoin))
- if(readmodel) then
+ if(read_external_model) then
npoinext = npoin
else
npoinext = 1
@@ -566,10 +568,10 @@
print *
print *,'Saving the grid in a text file...'
print *
- open(unit=55,file='OUTPUT_FILES/gridpoints.txt',status='unknown')
+ open(unit=55,file='OUTPUT_FILES/grid_points_and_model.txt',status='unknown')
write(55,*) npoin
do n = 1,npoin
- write(55,*) n,(coord(i,n), i=1,NDIM)
+ write(55,*) (coord(i,n), i=1,NDIM)
enddo
close(55)
endif
@@ -620,19 +622,15 @@
!
!---- eventuellement lecture d'un modele externe de vitesse et de densite
!
- if(readmodel) then
+ if(read_external_model) then
print *
print *,'Reading velocity and density model from external file...'
print *
- open(unit=55,file='OUTPUT_FILES/extmodel.txt',status='unknown')
+ open(unit=55,file='OUTPUT_FILES/grid_points_and_model.txt',status='old')
read(55,*) nbpoin
if(nbpoin /= npoin) stop 'Wrong number of points in input file'
do n = 1,npoin
- read(55,*) inump,rhoextread,vpextread,vsextread
- if(inump<1 .or. inump>npoin) stop 'Wrong point number'
- rhoext(inump) = rhoextread
- vpext(inump) = vpextread
- vsext(inump) = vsextread
+ read(55,*) xdummy,zdummy,rhoext(n),vpext(n),vsext(n)
enddo
close(55)
endif
@@ -643,7 +641,7 @@
call defarrays(vpext,vsext,rhoext,density,elastcoef, &
ibool,kmato,coord,npoin,rsizemin,rsizemax, &
cpoverdxmin,cpoverdxmax,lambdaSmin,lambdaSmax,lambdaPmin,lambdaPmax, &
- vpmin,vpmax,readmodel,nspec,numat)
+ vpmin,vpmax,read_external_model,nspec,numat)
! build the global mass matrix once and for all
rmass(:) = ZERO
@@ -652,7 +650,7 @@
do i = 1,NGLLX
iglob = ibool(i,j,ispec)
! if external density model
- if(readmodel) then
+ if(read_external_model) then
rhol = rhoext(iglob)
cpsquare = vpext(iglob)**2
else
@@ -981,7 +979,7 @@
do i = 1,NGLLX
!--- if external medium, get elastic parameters of current grid point
- if(readmodel) then
+ if(read_external_model) then
iglob = ibool(i,j,ispec)
cpl = vpext(iglob)
csl = vsext(iglob)
@@ -1185,7 +1183,7 @@
zgamma = xix(i,j,ispec) * jacobian(i,j,ispec)
! external velocity model
- if(readmodel) then
+ if(read_external_model) then
cpl = vpext(iglob)
csl = vsext(iglob)
rhol = rhoext(iglob)
@@ -1231,7 +1229,7 @@
zgamma = xix(i,j,ispec) * jacobian(i,j,ispec)
! external velocity model
- if(readmodel) then
+ if(read_external_model) then
cpl = vpext(iglob)
csl = vsext(iglob)
rhol = rhoext(iglob)
@@ -1283,7 +1281,7 @@
xxi = gammaz(i,j,ispec) * jacobian(i,j,ispec)
! external velocity model
- if(readmodel) then
+ if(read_external_model) then
cpl = vpext(iglob)
csl = vsext(iglob)
rhol = rhoext(iglob)
@@ -1335,7 +1333,7 @@
xxi = gammaz(i,j,ispec) * jacobian(i,j,ispec)
! external velocity model
- if(readmodel) then
+ if(read_external_model) then
cpl = vpext(iglob)
csl = vsext(iglob)
rhol = rhoext(iglob)
@@ -1640,7 +1638,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,ELASTIC)
+ boundvect,read_external_model,cutvect,nelemabs,numat,iptsdisp,nspec,ngnod,ELASTIC)
else if(ELASTIC .and. vecttype == 2) then
write(IOUT,*) 'drawing velocity field...'
@@ -1649,7 +1647,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,ELASTIC)
+ boundvect,read_external_model,cutvect,nelemabs,numat,iptsdisp,nspec,ngnod,ELASTIC)
else if(ELASTIC .and. vecttype == 3) then
write(IOUT,*) 'drawing acceleration field...'
@@ -1658,7 +1656,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,ELASTIC)
+ boundvect,read_external_model,cutvect,nelemabs,numat,iptsdisp,nspec,ngnod,ELASTIC)
! for acoustic medium
else if(.not. ELASTIC .and. vecttype == 1) then
@@ -1674,7 +1672,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,ELASTIC)
+ boundvect,read_external_model,cutvect,nelemabs,numat,iptsdisp,nspec,ngnod,ELASTIC)
else if(.not. ELASTIC .and. vecttype == 3) then
write(IOUT,*) 'drawing acoustic acceleration field from velocity potential...'
@@ -1686,7 +1684,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,ELASTIC)
+ boundvect,read_external_model,cutvect,nelemabs,numat,iptsdisp,nspec,ngnod,ELASTIC)
else
stop 'wrong field code for PostScript display'
@@ -1769,8 +1767,8 @@
'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,13('='),//5x, &
- 'Read external initial field or not . .(initialfield) = ',l6/5x, &
- 'Read external velocity model or not . . .(readmodel) = ',l6/5x, &
+ 'Read external initial field. . . . . .(initialfield) = ',l6/5x, &
+ 'Read external velocity model. .(read_external_model) = ',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, &
More information about the cig-commits
mailing list