[cig-commits] r8502 - seismo/2D/SPECFEM2D/trunk
walter at geodynamics.org
walter at geodynamics.org
Fri Dec 7 15:51:26 PST 2007
Author: walter
Date: 2007-12-07 15:51:25 -0800 (Fri, 07 Dec 2007)
New Revision: 8502
Added:
seismo/2D/SPECFEM2D/trunk/meshfem2D_circular_canyon.f90
seismo/2D/SPECFEM2D/trunk/meshfem2D_non_struct_2.f90
seismo/2D/SPECFEM2D/trunk/meshfem2D_non_struct_3.f90
Removed:
seismo/2D/SPECFEM2D/trunk/circ.f90
seismo/2D/SPECFEM2D/trunk/maille_non_struct_2.f90
seismo/2D/SPECFEM2D/trunk/maille_non_struct_3.f90
Modified:
seismo/2D/SPECFEM2D/trunk/meshfem2D.f90
seismo/2D/SPECFEM2D/trunk/specfem2D.f90
Log:
renamed obsolete meshers, cleaned the source code of meshfem2D_circular_canyon.f90, and added BibTeX information to the mesher and the solver.
Deleted: seismo/2D/SPECFEM2D/trunk/circ.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/circ.f90 2007-02-15 00:06:00 UTC (rev 8501)
+++ seismo/2D/SPECFEM2D/trunk/circ.f90 2007-12-07 23:51:25 UTC (rev 8502)
@@ -1,839 +0,0 @@
-!
-!=====================================================================
-!
-! P r e m a i l l e u r - 2 D
-! ---------------------------
-!
-! Version 2.1
-! -----------
-!
-! Dimitri Komatitsch
-!
-! Departement de Sismologie
-! Institut de Physique du Globe de Paris
-!
-! (c) Institut de Physique du Globe de Paris, Octobre 1996
-!
-!=====================================================================
-!
-
-! DK DK Mexico August 1999 : mise a jour format base de donnees
-
- program circ
-
- implicit none
-
-! max size of the model in elements
- integer, parameter :: mnx=7,mnz=7
-
- double precision, parameter :: pi=3.141592653589793d0
-
-! seuil pour considerer deux points comme confondus
- double precision, parameter :: rseuil=1.d-2
-
-! declare variables
- integer imaxabs,n2ana,itimetype,isource_type,nump1,nump2,nump3,nump4
- integer ndofn,ndime,ngnod,nnode,nbcnd,n1ana
- integer nofst,npgeo,nspel,nbmodeles,nbsources,nrec,lquad,isamp,nrec1,nrec2
- integer irec,imatnum,netyp,nxgll,nelemperio,nelemabs,nx,nz,i,j
- integer irepr,nrecsur3,nt,niter,itaff,itfirstaff,numerocourant,pointsdisp,isubsamp
-
- double precision R,theta_i,theta_init,delta_theta,eta_j,valseuil,freqmaxrep
- double precision f0,t0,xs,zs,angle,factor,dist,xoffs,zoffs
- double precision xrec,zrec,rho,cp,cs,anglerec
- double precision anglerec2,dt,alphanewm,betanewm,gammanewm
- double precision cutvect,cutcolor,scalex,scalez,sizemax,orig_x,orig_z
- double precision factorana,factorxsu
-
-! stockage de la grille curvi (x et z)
- integer, parameter :: npoinz1=(4*mnx+1)*(mnz+1), nelemz1=(4*mnx)*mnz
- double precision x1(0:4*mnx,0:mnz)
- double precision z1(0:4*mnx,0:mnz)
-
- integer, parameter :: npoinz3=(2*mnx+1)*(4*mnz+1), nelemz3=(2*mnx)*(4*mnz)
- double precision x3(0:2*mnx,0:4*mnz)
- double precision z3(0:2*mnx,0:4*mnz)
-
- integer, parameter :: npoinz4=(2*mnx+1)*(2*mnz+1), nelemz4=(2*mnx)*(2*mnz)
- double precision x4(0:2*mnx,0:2*mnz)
- double precision z4(0:2*mnx,0:2*mnz)
-
- integer, parameter :: npoinz1b=(2*mnx+1)*(mnz+1), nelemz1b=(2*mnx)*mnz
- double precision x1b(0:2*mnx,0:mnz)
- double precision z1b(0:2*mnx,0:mnz)
-
- integer, parameter :: npoinz2b=(mnx+1)*(2*mnz+1), nelemz2b=mnx*(2*mnz)
- double precision x2b(0:mnx,0:2*mnz)
- double precision z2b(0:mnx,0:2*mnz)
-
- integer, parameter :: npoinz3b=(4*mnx+1)*(4*mnz+1), nelemz3b=(4*mnx)*(4*mnz)
- double precision x3b(0:4*mnx,0:4*mnz)
- double precision z3b(0:4*mnx,0:4*mnz)
-
- integer, parameter :: npoinz4b=(2*mnx+1)*(2*mnz+1), nelemz4b=(2*mnx)*(2*mnz)
- double precision x4b(0:2*mnx,0:2*mnz)
- double precision z4b(0:2*mnx,0:2*mnz)
-
-! nombre max de points de maillage, et nombre exact d'elements
- integer, parameter :: npoin = npoinz1+npoinz3+npoinz4+ &
- npoinz1b+npoinz2b+npoinz3b+npoinz4b
- integer, parameter :: nelem = nelemz1+nelemz3+nelemz4+ &
- nelemz1b+nelemz2b+nelemz3b+nelemz4b
-
-! coordonnees geometriques des points
- double precision xpoint(npoin)
- double precision zpoint(npoin)
-
-! coordonnees des sommets de chaque element
- double precision x1e(nelem)
- double precision z1e(nelem)
- double precision x2e(nelem)
- double precision z2e(nelem)
- double precision x3e(nelem)
- double precision z3e(nelem)
- double precision x4e(nelem)
- double precision z4e(nelem)
-
-! numero des points des elements
- integer numpoin1(nelem)
- integer numpoin2(nelem)
- integer numpoin3(nelem)
- integer numpoin4(nelem)
-
-! nom du fichier GNUPLOT contenant la grille
- character(len=50) file1,title
-
- logical iexternal, aleatoire, topoplane, simulate, absstacey
- logical absorbhaut, absorbbas, absorbgauche, sismos
- logical absorbdroite, absorbstacey, absorbmodar, ifullabs
-
- logical display, ignuplot, ivectplot, icolorplot, imeshvect
- logical imeshcolor, imodelvect, iboundvect, interpol, isymbols, initialfield
- logical usletter,compenergy
-
- print *,'Nombre d''elements = ',nelem
- print *,'Nombre max de points = ',npoin
-
- nx = mnx
- nz = mnz
-
- R = 1.
-
-! ***************************************
-! *** ZONE DE DROITE
-! ***************************************
-
-! generer les points de base de l'interpolation lineaire (zone 1)
- theta_init = 3 * pi / 2.
- delta_theta = pi / 2.
- do i=0,4*nx
-
-! --- point de depart
- if(i < 2*nx) then
- x1(i,0) = 2.*R * real(i) / real(2*nx)
- z1(i,0) = - 2.*R
- else
- x1(i,0) = 2.*R
- z1(i,0) = - 2.*R * (1. - real(i - 2*nx) / real(2*nx))
- endif
-
-! --- point d'arrivee
- theta_i = theta_init + delta_theta * real(i) / real(4*nx)
- x1(i,nz) = dcos(theta_i)
- z1(i,nz) = dsin(theta_i)
-
-! --- points intermediaires par interpolation lineaire
- do j=1,nz-1
- eta_j = real(j) / real(nz)
- x1(i,j) = (1.-eta_j)*x1(i,0) + eta_j*x1(i,nz)
- z1(i,j) = (1.-eta_j)*z1(i,0) + eta_j*z1(i,nz)
- enddo
- enddo
-
-! generer zone de gauche (zone 3)
- do i=0,2*nx
- do j=0,4*nz
- x3(i,j) = 5. * real(i) / real(2*nx) + 2.
- if(j <= 2*nz) then
- z3(i,j) = 7. * real(j) / real(2*nz) - 9.
- else
- z3(i,j) = 2. * real(j-2*nz) / real(2*nz) - 2.
- endif
- enddo
- enddo
-
-! generer zone du bas (zone 4)
- do i=0,2*nx
- do j=0,2*nz
- x4(i,j) = 2. * real(i) / real(2*nx)
- z4(i,j) = 7. * real(j) / real(2*nz) - 9.
- enddo
- enddo
-
-! ***************************************
-! *** ZONE DE GAUCHE
-! ***************************************
-
-! generer les points de base de l'interpolation lineaire (zone 1)
- theta_init = pi / 4.
- delta_theta = pi / 4.
- do i=0,2*nx
-! --- point de depart
- x1b(i,0) = 2.*R * (real(i) / real(2*nx) - 1.)
- z1b(i,0) = - 2.*R
-
-! --- point d'arrivee
- theta_i = theta_init + delta_theta * real(i) / real(2*nx)
- x1b(i,nz) = - dcos(theta_i)
- z1b(i,nz) = - dsin(theta_i)
-
-! --- points intermediaires par interpolation lineaire
- do j=1,nz-1
- eta_j = real(j) / real(nz)
- x1b(i,j) = (1.-eta_j)*x1b(i,0) + eta_j*x1b(i,nz)
- z1b(i,j) = (1.-eta_j)*z1b(i,0) + eta_j*z1b(i,nz)
- enddo
- enddo
-
-! generer les points de base de l'interpolation lineaire (zone 2)
- theta_init = pi / 4.
- do j=0,2*nz
-! --- point de depart
- x2b(0,j) = - 2.*R
- z2b(0,j) = 2.*R * (real(j) / real(2*nz) - 1.)
-
-! --- point d'arrivee
- theta_i = theta_init - &
- delta_theta * real(j) / real(2*nz)
- x2b(nx,j) = - dcos(theta_i)
- z2b(nx,j) = - dsin(theta_i)
-
-! --- points intermediaires par interpolation lineaire
- do i=1,nx-1
- eta_j = real(i) / real(nx)
- x2b(i,j) = (1.-eta_j)*x2b(0,j) + eta_j*x2b(nx,j)
- z2b(i,j) = (1.-eta_j)*z2b(0,j) + eta_j*z2b(nx,j)
- enddo
-
- enddo
-
-! generer zone de gauche (zone 3)
- do i=0,4*nx
- do j=0,4*nz
- x3b(i,j) = 10. * real(i) / real(4*nx) - 12.
- if(j <= 2*nz) then
- z3b(i,j) = 7. * real(j) / real(2*nz) - 9.
- else
- z3b(i,j) = 2. * real(j-2*nz) / real(2*nz) - 2.
- endif
- enddo
- enddo
-
-! generer zone du bas (zone 4)
- do i=0,2*nx
- do j=0,2*nz
- x4b(i,j) = 2. * real(i) / real(2*nx) - 2.
- z4b(i,j) = 7. * real(j) / real(2*nz) - 9.
- enddo
- enddo
-
-! ***
-! *** generer un fichier 'GNUPLOT' pour le controle de la grille ***
-! ***
-
- write(*,*)' '
- write(*,*)' Ecriture de la grille format GNUPLOT...'
-
- file1='grid.GNU'
-
- open(unit=20,file=file1,status='unknown')
-
-! *** dessiner la zone 1
- do j=0,nz
- do i=0,4*nx-1
- write(20,*) real(x1(i,j)),real(z1(i,j))
- write(20,*) real(x1(i+1,j)),real(z1(i+1,j))
- write(20,100)
- enddo
- enddo
-
- do i=0,4*nx
- do j=0,nz-1
- write(20,*) real(x1(i,j)),real(z1(i,j))
- write(20,*) real(x1(i,j+1)),real(z1(i,j+1))
- write(20,100)
- enddo
- enddo
-
-! *** dessiner la zone 3
- do j=0,4*nz
- do i=0,2*nx-1
- write(20,*) real(x3(i,j)),real(z3(i,j))
- write(20,*) real(x3(i+1,j)),real(z3(i+1,j))
- write(20,100)
- enddo
- enddo
-
- do i=0,2*nx
- do j=0,4*nz-1
- write(20,*) real(x3(i,j)),real(z3(i,j))
- write(20,*) real(x3(i,j+1)),real(z3(i,j+1))
- write(20,100)
- enddo
- enddo
-
-! *** dessiner la zone 4
- do j=0,2*nz
- do i=0,2*nx-1
- write(20,*) real(x4(i,j)),real(z4(i,j))
- write(20,*) real(x4(i+1,j)),real(z4(i+1,j))
- write(20,100)
- enddo
- enddo
-
- do i=0,2*nx
- do j=0,2*nz-1
- write(20,*) real(x4(i,j)),real(z4(i,j))
- write(20,*) real(x4(i,j+1)),real(z4(i,j+1))
- write(20,100)
- enddo
- enddo
-
-! *** dessiner la zone 1
- do j=0,nz
- do i=0,2*nx-1
- write(20,*) real(x1b(i,j)),real(z1b(i,j))
- write(20,*) real(x1b(i+1,j)),real(z1b(i+1,j))
- write(20,100)
- enddo
- enddo
-
- do i=0,2*nx
- do j=0,nz-1
- write(20,*) real(x1b(i,j)),real(z1b(i,j))
- write(20,*) real(x1b(i,j+1)),real(z1b(i,j+1))
- write(20,100)
- enddo
- enddo
-
-! *** dessiner la zone 2
- do j=0,2*nz
- do i=0,nx-1
- write(20,*) real(x2b(i,j)),real(z2b(i,j))
- write(20,*) real(x2b(i+1,j)),real(z2b(i+1,j))
- write(20,100)
- enddo
- enddo
-
- do i=0,nx
- do j=0,2*nz-1
- write(20,*) real(x2b(i,j)),real(z2b(i,j))
- write(20,*) real(x2b(i,j+1)),real(z2b(i,j+1))
- write(20,100)
- enddo
- enddo
-
-! *** dessiner la zone 3
- do j=0,4*nz
- do i=0,4*nx-1
- write(20,*) real(x3b(i,j)),real(z3b(i,j))
- write(20,*) real(x3b(i+1,j)),real(z3b(i+1,j))
- write(20,100)
- enddo
- enddo
-
- do i=0,4*nx
- do j=0,4*nz-1
- write(20,*) real(x3b(i,j)),real(z3b(i,j))
- write(20,*) real(x3b(i,j+1)),real(z3b(i,j+1))
- write(20,100)
- enddo
- enddo
-
-! *** dessiner la zone 4
- do j=0,2*nz
- do i=0,2*nx-1
- write(20,*) real(x4b(i,j)),real(z4b(i,j))
- write(20,*) real(x4b(i+1,j)),real(z4b(i+1,j))
- write(20,100)
- enddo
- enddo
-
- do i=0,2*nx
- do j=0,2*nz-1
- write(20,*) real(x4b(i,j)),real(z4b(i,j))
- write(20,*) real(x4b(i,j+1)),real(z4b(i,j+1))
- write(20,100)
- enddo
- enddo
-
- close(20)
-
- write(*,*)' Fin ecriture de la grille format GNUPLOT'
- write(*,*)' '
-
- 100 format('')
-
-! ***
-! *** generer la liste des points geometriques
-! ***
-
- numerocourant = 1
-
-! *** zone 1
- do j=0,nz
- do i=0,4*nx
- xpoint(numerocourant) = x1(i,j)
- zpoint(numerocourant) = z1(i,j)
- numerocourant = numerocourant + 1
- enddo
- enddo
-
-! *** zone 3
- do j=0,4*nz
- do i=0,2*nx
- xpoint(numerocourant) = x3(i,j)
- zpoint(numerocourant) = z3(i,j)
- numerocourant = numerocourant + 1
- enddo
- enddo
-
-! *** zone 4
- do j=0,2*nz
- do i=0,2*nx
- xpoint(numerocourant) = x4(i,j)
- zpoint(numerocourant) = z4(i,j)
- numerocourant = numerocourant + 1
- enddo
- enddo
-
-! *** zone 1
- do j=0,nz
- do i=0,2*nx
- xpoint(numerocourant) = x1b(i,j)
- zpoint(numerocourant) = z1b(i,j)
- numerocourant = numerocourant + 1
- enddo
- enddo
-
-! *** zone 2
- do j=0,2*nz
- do i=0,nx
- xpoint(numerocourant) = x2b(i,j)
- zpoint(numerocourant) = z2b(i,j)
- numerocourant = numerocourant + 1
- enddo
- enddo
-
-! *** zone 3
- do j=0,4*nz
- do i=0,4*nx
- xpoint(numerocourant) = x3b(i,j)
- zpoint(numerocourant) = z3b(i,j)
- numerocourant = numerocourant + 1
- enddo
- enddo
-
-! *** zone 4
- do j=0,2*nz
- do i=0,2*nx
- xpoint(numerocourant) = x4b(i,j)
- zpoint(numerocourant) = z4b(i,j)
- numerocourant = numerocourant + 1
- enddo
- enddo
-
- print *,'nb de points stockes = ',numerocourant - 1
-
-! ***
-! *** generer la liste des elements
-! ***
-
- numerocourant = 1
- imaxabs = 0
-
-! *** zone 1
- do j=0,nz-1
- do i=0,4*nx-1
- x1e(numerocourant) = x1(i,j)
- z1e(numerocourant) = z1(i,j)
- x2e(numerocourant) = x1(i+1,j)
- z2e(numerocourant) = z1(i+1,j)
- x3e(numerocourant) = x1(i+1,j+1)
- z3e(numerocourant) = z1(i+1,j+1)
- x4e(numerocourant) = x1(i,j+1)
- z4e(numerocourant) = z1(i,j+1)
- numerocourant = numerocourant + 1
- enddo
- enddo
-
-! *** zone 3
- do j=0,4*nz-1
- do i=0,2*nx-1
- x1e(numerocourant) = x3(i,j)
- z1e(numerocourant) = z3(i,j)
- x2e(numerocourant) = x3(i+1,j)
- z2e(numerocourant) = z3(i+1,j)
- x3e(numerocourant) = x3(i+1,j+1)
- z3e(numerocourant) = z3(i+1,j+1)
- x4e(numerocourant) = x3(i,j+1)
- z4e(numerocourant) = z3(i,j+1)
- numerocourant = numerocourant + 1
- enddo
- enddo
-
-! *** zone 4
- do j=0,2*nz-1
- do i=0,2*nx-1
- x1e(numerocourant) = x4(i,j)
- z1e(numerocourant) = z4(i,j)
- x2e(numerocourant) = x4(i+1,j)
- z2e(numerocourant) = z4(i+1,j)
- x3e(numerocourant) = x4(i+1,j+1)
- z3e(numerocourant) = z4(i+1,j+1)
- x4e(numerocourant) = x4(i,j+1)
- z4e(numerocourant) = z4(i,j+1)
- numerocourant = numerocourant + 1
- enddo
- enddo
-
-! *** zone 1
- do j=0,nz-1
- do i=0,2*nx-1
- x1e(numerocourant) = x1b(i,j)
- z1e(numerocourant) = z1b(i,j)
- x2e(numerocourant) = x1b(i+1,j)
- z2e(numerocourant) = z1b(i+1,j)
- x3e(numerocourant) = x1b(i+1,j+1)
- z3e(numerocourant) = z1b(i+1,j+1)
- x4e(numerocourant) = x1b(i,j+1)
- z4e(numerocourant) = z1b(i,j+1)
- numerocourant = numerocourant + 1
- enddo
- enddo
-
-! *** zone 2
- do j=0,2*nz-1
- do i=0,nx-1
- x1e(numerocourant) = x2b(i,j)
- z1e(numerocourant) = z2b(i,j)
- x2e(numerocourant) = x2b(i+1,j)
- z2e(numerocourant) = z2b(i+1,j)
- x3e(numerocourant) = x2b(i+1,j+1)
- z3e(numerocourant) = z2b(i+1,j+1)
- x4e(numerocourant) = x2b(i,j+1)
- z4e(numerocourant) = z2b(i,j+1)
- numerocourant = numerocourant + 1
- enddo
- enddo
-
-! *** zone 3
- do j=0,4*nz-1
- do i=0,4*nx-1
- x1e(numerocourant) = x3b(i,j)
- z1e(numerocourant) = z3b(i,j)
- x2e(numerocourant) = x3b(i+1,j)
- z2e(numerocourant) = z3b(i+1,j)
- x3e(numerocourant) = x3b(i+1,j+1)
- z3e(numerocourant) = z3b(i+1,j+1)
- x4e(numerocourant) = x3b(i,j+1)
- z4e(numerocourant) = z3b(i,j+1)
- numerocourant = numerocourant + 1
- enddo
- enddo
-
-! *** zone 4
- do j=0,2*nz-1
- do i=0,2*nx-1
- x1e(numerocourant) = x4b(i,j)
- z1e(numerocourant) = z4b(i,j)
- x2e(numerocourant) = x4b(i+1,j)
- z2e(numerocourant) = z4b(i+1,j)
- x3e(numerocourant) = x4b(i+1,j+1)
- z3e(numerocourant) = z4b(i+1,j+1)
- x4e(numerocourant) = x4b(i,j+1)
- z4e(numerocourant) = z4b(i,j+1)
- numerocourant = numerocourant + 1
- enddo
- enddo
-
- print *,'nb d''elements stockes = ',numerocourant - 1
-
-! ***
-! *** creation des elements sous forme topologique
-! ***
-
- write(*,*)' '
- write(*,*)' Creation de la topologie des elements...'
-
- file1='topoelements.txt'
-
- do i=1,nelem
-
-! recherche point 1
- do j=1,npoin
- dist = dsqrt((x1e(i)-xpoint(j))**2 + &
- (z1e(i)-zpoint(j))**2)
- if(dist <= rseuil) then
- nump1 = j
- goto 401
- endif
- enddo
- stop 'point not found !'
- 401 continue
-
-! recherche point 2
- do j=1,npoin
- dist = dsqrt((x2e(i)-xpoint(j))**2 + &
- (z2e(i)-zpoint(j))**2)
- if(dist <= rseuil) then
- nump2 = j
- goto 402
- endif
- enddo
- stop 'point not found !'
- 402 continue
-
-! recherche point 3
- do j=1,npoin
- dist = dsqrt((x3e(i)-xpoint(j))**2 + &
- (z3e(i)-zpoint(j))**2)
- if(dist <= rseuil) then
- nump3 = j
- goto 403
- endif
- enddo
- stop 'point not found !'
- 403 continue
-
-! recherche point 4
- do j=1,npoin
- dist = dsqrt((x4e(i)-xpoint(j))**2 + &
- (z4e(i)-zpoint(j))**2)
- if(dist <= rseuil) then
- nump4 = j
- goto 404
- endif
- enddo
- stop 'point not found !'
- 404 continue
-
- numpoin1(i) = nump1
- numpoin2(i) = nump2
- numpoin3(i) = nump3
- numpoin4(i) = nump4
-
- enddo
-
-! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-! *** generation de la base de donnees
-
- write(*,*)
- write(*,*) 'Generation de la base de donnees...'
-
- open(unit=15,file='OUTPUT_FILES/Database',status='unknown')
-
- title = 'Modele Canyon Paco'
- write(15,*) '#'
- write(15,*) '# Base de Donnees pour Specfem - Premailleur Fortran 90'
- write(15,*) '# ',title
- write(15,*) '# Dimitri Komatitsch, (c) EPS - Harvard February 1998'
- write(15,*) '#'
-
- write(15,*) 'Titre simulation'
- write(15,40) title
-
- ndofn = 2
- ndime = 2
- ngnod = 4
- nnode = 4
- nbcnd = 0
- nofst = 0
- npgeo = npoin
- nspel = nelem
- nbmodeles = 1
- nbsources = 1
- nrec = 150
- lquad = 1
- iexternal = .false.
- aleatoire = .false.
- topoplane = .false.
- simulate = .false.
-
- absorbhaut = .false.
- absorbbas = .false.
- absorbgauche = .false.
- absorbdroite = .false.
- absorbstacey = .true.
- absorbmodar = .false.
- ifullabs = .false.
-
- sismos = .true.
- isamp = 20
- nrec1 = nrec
- nrec2 = 0
- anglerec = 0.
- anglerec2 = 0.
- irepr = 1
- nrecsur3 = nrec / 3
-
- nt = 20000
- dt = 0.625e-3
- niter = 1
- alphanewm = 0.
- betanewm = 0.
- gammanewm = 0.5
- display = .true.
- ignuplot = .false.
- ivectplot = .true.
- icolorplot = .false.
- imeshvect = .true.
- imeshcolor = .false.
- imodelvect = .false.
- iboundvect = .false.
- interpol = .true.
- isymbols = .true.
-
-!! DK DK Mexico August 1999, temporarily suppress external field
- initialfield = .true.
- initialfield = .false.
-
- itaff = 2000
- itfirstaff = 5
- cutvect = 1.
- cutcolor = 2.2
- scalex = 1.
- scalez = 1.
- sizemax = 1.
- pointsdisp = 7
- isubsamp = 2
- orig_x = 2.3
- orig_z = 3.4
- factorana = 50000.
- factorxsu = 3.5
- n1ana = 1
- n2ana = nrec1
-
- write(15,*) 'ndofn ndime npgeo'
- write(15,*) ndofn,ndime,npgeo
-
- write(15,*) 'display ignuplot interpol'
- write(15,*) display,ignuplot,interpol
-
- write(15,*) 'itaff itfirstaff icolor inumber'
- write(15,*) itaff,itfirstaff,0,0
-
- write(15,*) 'ivectplot imeshvect imodelvect iboundvect cutvect isubsamp'
- write(15,*) ivectplot,imeshvect,imodelvect,iboundvect,cutvect,isubsamp
-
- usletter = .true.
- write(15,*) 'scalex scalez sizemax angle rapport USletter'
- write(15,*) scalex,scalez,sizemax,20.,0.40,usletter
-
- write(15,*) 'orig_x orig_z isymbols'
- write(15,*) orig_x,orig_z,isymbols
-
- valseuil = 5.00
- freqmaxrep = 100.
- write(15,*) 'valseuil freqmaxrep'
- write(15,*) valseuil,freqmaxrep
-
- write(15,*) 'sismos nrec nrec1 nrec2 isamp'
- write(15,*) sismos,nrec,nrec1,nrec2,isamp
-
- write(15,*) 'irepr anglerec anglerec2'
- write(15,*) irepr,anglerec,anglerec2
-
- compenergy = .false.
- absstacey = .true.
- write(15,*) 'topoplane absstacey compenergy'
- write(15,*) topoplane,absstacey,compenergy
-
- write(15,*) 'initialfield factorana factorxsu n1ana n2ana'
- write(15,*) initialfield,factorana,factorxsu,n1ana,n2ana
-
- write(15,*) 'isismostype ivecttype iaffinfo'
- write(15,*) '1, 1, 40'
- write(15,*) 'ireadmodel ioutputgrid iavs ivisual3'
- write(15,*) 'F, F, F, F'
-
- write(15,*) 'iexec iecho'
- write(15,*) '1 1'
-
- write(15,*) 'ncycl dtinc niter'
- write(15,*) nt,dt,niter
-
- write(15,*) 'alpha beta gamma'
- write(15,*) alphanewm,betanewm,gammanewm
-
- nbsources = 1
- write(15,*) 'nltfl (number of force or pressure sources)'
- write(15,*) nbsources
-
- itimetype = 6
- isource_type = 2
- f0 = 2.
- t0 = 0.55
- xs = +1.
- zs = -2.
- angle = 0.
- factor = 1.
- xoffs = 12.
- zoffs = 9.
- write(15,*) 'Collocated forces and/or pressure sources:'
- write(15,*) itimetype,isource_type, &
- xs+xoffs,zs+zoffs,f0,t0,factor,angle,0
-
- write(15,*) 'Receivers (number, angle, position in meters)'
- do irec=1,nrec
- if(irec <= nrecsur3) then
- xrec = 2.*dble(irec-1)/dble(nrecsur3-1) + 9.
- zrec = 9.
- else if(irec >= 2*nrecsur3) then
- xrec = 2.*dble(irec-2*nrecsur3)/dble(nrecsur3) + 13.
- zrec = 9.
- else
- angle = pi + pi*dble(irec-nrecsur3)/dble(nrecsur3)
- xrec = 12. + dcos(angle)
- zrec = 9. + dsin(angle)
- endif
- write(15,*) irec,xrec,zrec
- enddo
-
- write(15,*) 'Coordinates of spectral control points'
- do i=1,npoin
- write(15,*) i,xpoint(i)+xoffs,zpoint(i)+zoffs
- enddo
-
- netyp = 2
- nxgll = 6
- nelemperio = 0
- nelemabs = 0
-
- write(15,*) 'params spectraux'
- write(15,*) netyp,nbmodeles,ngnod,nxgll,nxgll,nspel,pointsdisp, &
- nelemabs,nelemperio
-
- write(15,*) 'Material sets (num 0 rho vp vs 0 0)'
- rho = 1.
- cp = 2.
- cs = 1.
- write(15,*) nbmodeles,0,rho,cp,cs,0,0
-
- write(15,*) 'Spectral elements topology'
-
- imatnum = 1
-
- do i=1,nspel
- write(15,*) i,imatnum,numpoin1(i),numpoin2(i),numpoin3(i), &
- numpoin4(i)
- enddo
-
- close(15)
-
- 40 format(a50)
-
-! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
- end program circ
Deleted: seismo/2D/SPECFEM2D/trunk/maille_non_struct_2.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/maille_non_struct_2.f90 2007-02-15 00:06:00 UTC (rev 8501)
+++ seismo/2D/SPECFEM2D/trunk/maille_non_struct_2.f90 2007-12-07 23:51:25 UTC (rev 8502)
@@ -1,962 +0,0 @@
-!=====================================================================
-!
-! P r e m a i l l e u r F o r t r a n 9 0
-! -------------------------------------------
-!
-! Version 3.0
-! -----------
-!
-! Dimitri Komatitsch
-! Department of Earth and Planetary Sciences - Harvard University
-!
-! (c) August 1998
-!
-!=====================================================================
-
-!
-! *** Version optimisee avec maillage non structure Jacques Muller - Elf ***
-! *** Raffinement d'un facteur 2 en surface ***
-!
-
- program maille_non_struct_2
-
- implicit none
-
-! definir les tableaux pour allocation dynamique
-
-! coordinates of the grid points
- double precision, allocatable :: x(:,:),z(:,:)
-
-! variables needed to compute the transformation
- double precision, allocatable :: psi(:),eta(:),absx(:), &
- a00(:),a01(:),valeta(:),bot0(:),top0(:)
-
-! stockage du modele de vitesse et densite
- double precision, allocatable :: rho(:),cp(:),cs(:)
-
-! the topography data
- double precision, allocatable :: xtopo(:),ztopo(:),coefs_topo(:)
-
-! arrays for the source
- double precision, allocatable :: xs(:),zs(:),f0(:),t0(:),angle(:),factor(:)
- integer, allocatable :: isource_type(:),itimetype(:)
-
-! arrays for the receivers
- double precision, allocatable :: xrec(:),zrec(:)
-
- character(len=50) interffile,topofile,title
- character(len=15) junk
-
- integer imatnum,inumabs,inumelem,nelemperio,nxgll,netyp
- integer icodehaut,icodebas,icodegauche,icodedroite
- integer nelemabs,npgeo,nspec,ntopo,nspecvolume,nspecWz
- integer k,ix,iz,irec,i,j,iadd
- integer imodele,nbmodeles,iaffinfo
- integer itaff,itfirstaff,pointsdisp,isubsamp,nrec,n1ana,n2ana
- integer irepr,nrec1,nrec2,isamp,nbsources,isismostype,ivecttype
- integer ngnod,nt,niter,idegpoly,nx,nz
- integer icodematread
-
- double precision valseuil,freqmaxrep,ratio
- double precision tang1,tangN
- double precision orig_x,orig_z,sizemax,cutvect,scalex,scalez
- double precision factorxsu,factorana,xspacerec,zspacerec
- double precision anglerec,anglerec2,xmin,xmax
- double precision xfin,zfin,xfin2,zfin2,xdeb,zdeb,xdeb2,zdeb2
- double precision alphanewm,betanewm,gammanewm,dt
- double precision rhoread,cpread,csread,aniso3read,aniso4read
-
- logical interpol,ignuplot,ireadmodel,iavs,ivisual3,ioutputgrid
- logical abshaut,absbas,absgauche,absdroite,absstacey
- logical periohaut,periogauche
- logical sismos,isources_surf,ienreg_surf,ienreg_surf2,display
- logical ivectplot,imeshvect
- logical topoplane,iexec,initialfield
- logical imodelvect,iboundvect,usletter,compenergy
-
- integer, external :: num
- double precision, external :: bottom,spl,dens
-
- double precision, parameter :: zero = 0.d0, one = 1.d0
-
-! simulation a 2D
- integer, parameter :: ndime = 2
- integer, parameter :: ndofn = 2
-
-! --- code des numeros d'aretes pour les bords absorbants
- integer, parameter :: iaretebas = 1
- integer, parameter :: iaretedroite = 2
- integer, parameter :: iaretehaut = 3
- integer, parameter :: iaretegauche = 4
-
-! DK DK DK ajout Elf : extraction de la topo du fichier SEP
-!! call system('rm -f topo_from_SEP.dat topo_SEP_maille90.dat ; xextract_topo')
-
- print *
- print *,' *** Version optimisee avec maillage non structure ***'
- print *,' *** Raffinement d''un facteur 2 en surface ***'
- print *
-
-! ***
-! *** read the parameter file
-! ***
-
- print *,' Reading the parameter file ... '
- print *
-
- open(unit=10,file='DATA/Par_file',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,topofile
- read(10,3)junk,interffile
-
- 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,nz
- read(10,2)junk,idegpoly
- read(10,2)junk,ngnod
- read(10,1)junk,ratio
- read(10,4)junk,topoplane
- read(10,4)junk,initialfield
- read(10,4)junk,ireadmodel
- read(10,4)junk,iexec
-
-! DK DK forcer pour Elf
- ngnod = 9
- topoplane = .false.
- initialfield = .false.
-
-! pour le non structure, verifier la coherence du maillage
- if(nx < 2) stop 'nx must be greater or equal to 2'
- if(nz < 2) stop 'nz must be greater or equal to 2'
- if(mod(nx,2) /= 0) stop 'nx must be even'
-
-! multiplier par 2 pour implementer le deraffinement non conforme
- nx = nx * 2
- nz = nz * 2
-
-! multiplier par 2 pour elements 9 noeuds
- nx = nx * 2
- nz = nz * 2
-
-! 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
- read(10,4)junk,absstacey
- read(10,4)junk,periohaut
- read(10,4)junk,periogauche
-
-! DK DK forcer pour Elf
- abshaut = .false.
- periohaut = .false.
- periogauche = .false.
-
-! skip comment
- read(10,*)
- read(10,*)
- read(10,*)
-
-! read time step parameters
- read(10,2)junk,nt
- read(10,1)junk,dt
- read(10,2)junk,niter
- read(10,1)junk,alphanewm
- read(10,1)junk,betanewm
- read(10,1)junk,gammanewm
-
-! skip comment
- read(10,*)
- read(10,*)
- read(10,*)
-
-! read source parameters
- read(10,2)junk,nbsources
- read(10,4)junk,isources_surf
- read(10,1)junk,valseuil
- read(10,1)junk,freqmaxrep
- print *,'Nb de sources a lire : ',nbsources
-
- allocate(xs(nbsources))
- allocate(zs(nbsources))
- allocate(f0(nbsources))
- allocate(t0(nbsources))
- allocate(isource_type(nbsources))
- allocate(itimetype(nbsources))
- allocate(angle(nbsources))
- allocate(factor(nbsources))
-
- do i=1,nbsources
- read(10,*)
- read(10,1)junk,xs(i)
- read(10,1)junk,zs(i)
- read(10,1)junk,f0(i)
- read(10,1)junk,t0(i)
- read(10,2)junk,isource_type(i)
- read(10,2)junk,itimetype(i)
- read(10,1)junk,angle(i)
- read(10,1)junk,factor(i)
-
- print *
- print *,' Source #',i
- print *,'Position xs, zs = ',xs(i),zs(i)
- print *,'Frequency, delay = ',f0(i),t0(i)
- print *,'Source type (1=force 2=explo) : ', &
- isource_type(i)
- print *,'Angle of the source if force = ',angle(i)
- print *,'Multiplying factor = ',factor(i)
- enddo
-
-! skip comment
- read(10,*)
- read(10,*)
- read(10,*)
-
-! read receivers line parameters
- read(10,4)junk,sismos
- read(10,2)junk,isamp
- read(10,2)junk,isismostype
- read(10,2)junk,irepr
- read(10,*)
- read(10,2)junk,nrec1
- read(10,1)junk,xdeb
- read(10,1)junk,zdeb
- read(10,1)junk,xfin
- read(10,1)junk,zfin
- read(10,4)junk,ienreg_surf
- read(10,1)junk,anglerec
- read(10,*)
- read(10,2)junk,nrec2
- read(10,1)junk,xdeb2
- read(10,1)junk,zdeb2
- read(10,1)junk,xfin2
- read(10,1)junk,zfin2
- read(10,4)junk,ienreg_surf2
- read(10,1)junk,anglerec2
- read(10,*)
- read(10,1)junk,factorxsu
- read(10,2)junk,n1ana
- read(10,2)junk,n2ana
- read(10,1)junk,factorana
-
-! determination et affichage position ligne de receivers
- if(nrec2 < 0) stop 'negative value of nrec2 !'
-
- if(nrec2 == 0) then
- nrec = nrec1
- else
- nrec = nrec1 + nrec2
- endif
-
-! DK DK forcer pour Elf
- n1ana = 1
- n2ana = nrec
-
- allocate(xrec(nrec))
- allocate(zrec(nrec))
-
- if(nrec2 == 0) then
- print *
- print *,'There are ',nrec,' receivers on a single line'
- xspacerec=(xfin-xdeb)/dble(nrec-1)
- zspacerec=(zfin-zdeb)/dble(nrec-1)
- do i=1,nrec
- xrec(i) = xdeb + dble(i-1)*xspacerec
- zrec(i) = zdeb + dble(i-1)*zspacerec
- enddo
- else
- print *
- print *,'There are ',nrec,' receivers on two lines'
- print *,'First line contains ',nrec1,' receivers'
- print *,'Second line contains ',nrec2,' receivers'
- xspacerec=(xfin-xdeb)/dble(nrec1-1)
- zspacerec=(zfin-zdeb)/dble(nrec1-1)
- do i=1,nrec1
- xrec(i) = xdeb + dble(i-1)*xspacerec
- zrec(i) = zdeb + dble(i-1)*zspacerec
- enddo
- xspacerec=(xfin2-xdeb2)/dble(nrec2-1)
- zspacerec=(zfin2-zdeb2)/dble(nrec2-1)
- do i=1,nrec2
- xrec(i+nrec1) = xdeb2 + dble(i-1)*xspacerec
- zrec(i+nrec1) = zdeb2 + dble(i-1)*zspacerec
- enddo
- endif
-
-! skip comment
- read(10,*)
- read(10,*)
- read(10,*)
-
-! read display parameters
- read(10,4)junk,display
- read(10,2)junk,itaff
- read(10,2)junk,itfirstaff
- read(10,2)junk,iaffinfo
- read(10,4)junk,ivectplot
- 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,pointsdisp
- read(10,2)junk,isubsamp
- read(10,1)junk,scalex
- read(10,1)junk,scalez
- read(10,1)junk,sizemax
- read(10,4)junk,usletter
- read(10,1)junk,orig_x
- read(10,1)junk,orig_z
- read(10,4)junk,ignuplot
- read(10,4)junk,iavs
- read(10,4)junk,ivisual3
- read(10,4)junk,ioutputgrid
- read(10,4)junk,compenergy
-
-! DK DK forcer pour Elf
- ignuplot = .false.
- iavs = .false.
- ivisual3 = .false.
- compenergy = .false.
-
-! skip comment
- read(10,*)
- read(10,*)
- read(10,*)
-
-! lecture des differents modeles de materiaux
-
- read(10,2)junk,nbmodeles
- if(nbmodeles <= 0) stop 'Negative number of models not allowed !!'
-
- allocate(rho(nbmodeles))
- allocate(cp(nbmodeles))
- allocate(cs(nbmodeles))
-
- rho(:) = 0.d0
- cp(:) = 0.d0
- cs(:) = 0.d0
-
- do imodele=1,nbmodeles
- read(10,*) i,icodematread,rhoread,cpread,csread,aniso3read,aniso4read
- if(i<1 .or. i>nbmodeles) stop 'Wrong material set number'
- 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'
- enddo
-
- print *
- print *, 'Nb de modeles de roche = ',nbmodeles
- print *
- do i=1,nbmodeles
- print *,'Modele #',i,' isotrope'
- print *,'rho,cp,cs = ',rho(i),cp(i),cs(i)
- enddo
- print *
-
- close(10)
-
- print *
- print *,' Parameter file successfully read... '
-
-! --------- fin lecture fichier parametres --------------
-
- allocate(psi(0:nx))
- allocate(eta(0:nz))
- allocate(absx(0:nx))
- allocate(a00(0:nz))
- allocate(a01(0:nz))
- allocate(valeta(0:nz))
- allocate(bot0(0:nx))
- allocate(top0(0:nx))
-
-! calcul des points regulierement espaces
- do i=0,nx
- psi(i) = i/dble(nx)
- enddo
- do j=0,nz
- eta(j) = j/dble(nz)
- enddo
-
-! quelques verifications de base a faire
-
- if(ngnod /= 9) stop 'erreur ngnod different de 9 !!'
-
-! calcul du nombre total d'elements spectraux, absorbants et periodiques
- nspecvolume = (nx/2/2)*((nz-4)/2/2)
- nspecWz = 3*(nx/2/2)
- nspec = nspecvolume + nspecWz
- nelemperio = 0
-
- if(absgauche .or. absdroite .or. absbas) then
- nelemabs = 2 * (nz/4 - 2) + nx/4 + 2 + 2
- else
- nelemabs = 0
- endif
-
- print *
- print *,'Le maillage comporte ',nspec,' elements spectraux (nx = ',nx/4, &
- ' nz = ',nz/4,')'
- print *,'soit ',nspecvolume,' elements spectraux dans le volume'
- print *,'et ',nspecWz,' elements spectraux dans la couche Wz'
- print *,'Chaque element comporte ',idegpoly+1,' points dans chaque direction'
- print *,'Le nombre maximum de points theorique est ',nspec*(idegpoly+1)**ndime
- print *,'Les elements de controle sont des elements ',ngnod,' noeuds'
- print *
-
-!------------------------------------------------------
-
- allocate(x(0:nx,0:nz))
- allocate(z(0:nx,0:nz))
-
- x(:,:)=0.d0
- z(:,:)=0.d0
-
-! get topography data from external file
- print *,'Reading topography from file ',topofile
- open(unit=15,file=topofile,status='old')
- read(15,*) ntopo
- if (ntopo < 2) stop 'Not enough topography points (min 2)'
- print *,'Reading ',ntopo,' points from topography file'
- print *
-
- allocate(xtopo(ntopo))
- allocate(ztopo(ntopo))
- allocate(coefs_topo(ntopo))
-
- do i=1,ntopo
- read(15,*) xtopo(i),ztopo(i)
- enddo
- close(15)
-
-! check the values read
- print *
- print *, 'Topography data points (x,z)'
- print *, '----------------------------'
- print *
- print *, 'Topo 1 = (',xtopo(1),',',ztopo(1),')'
- print *, 'Topo ntopo = (',xtopo(ntopo),',',ztopo(ntopo),')'
-
-!--- calculate the spline function for the topography
-!--- imposer les tangentes aux deux bords
- tang1 = (ztopo(2)-ztopo(1))/(xtopo(2)-xtopo(1))
- tangN = (ztopo(ntopo)-ztopo(ntopo-1))/(xtopo(ntopo)-xtopo(ntopo-1))
- call spline(xtopo,ztopo,ntopo,tang1,tangN,coefs_topo)
-
-! *** afficher limites du modele lu
- print *
- print *, 'Limites absolues modele fichier topo :'
- print *
- print *, 'Xmin = ',minval(xtopo),' Xmax = ',maxval(xtopo)
- print *, 'Zmin = ',minval(ztopo),' Zmax = ',maxval(ztopo)
- print *
-
-! *** modifier sources pour position par rapport a la surface
- print *
- print *, 'Position (x,z) des ',nbsources,' sources'
- print *
- do i=1,nbsources
-
-! DK DK DK Elf : position source donnee en profondeur par rapport a la topo
- zs(i) = spl(xs(i),xtopo,ztopo,coefs_topo,ntopo) - zs(i)
-
- if(isources_surf) zs(i) = spl(xs(i),xtopo,ztopo,coefs_topo,ntopo)
- print *, 'Source ',i,' = ',xs(i),zs(i)
- enddo
-
-! *** modifier recepteurs pour enregistrement en surface
- print *
- print *, 'Position (x,z) des ',nrec,' receivers'
- print *
- do irec=1,nrec
-
-! DK DK DK Elf : distinguer les deux lignes de recepteurs
- if(irec <= nrec1) then
- if(ienreg_surf) zrec(irec) = spl(xrec(irec),xtopo,ztopo,coefs_topo,ntopo)
- else
- if(ienreg_surf2) zrec(irec) = spl(xrec(irec),xtopo,ztopo,coefs_topo,ntopo)
- endif
- print *, 'Receiver ',irec,' = ',xrec(irec),zrec(irec)
-
- enddo
-
-!--- definition du maillage suivant X
- do ix=0,nx
- absx(ix) = dens(ix,psi,xmin,xmax,nx)
- enddo
-
-! *** une seule zone
-
- do iz=0,nz
-
-! DK DK DK densification sinusoidale ici en vertical
- valeta(iz) = eta(iz) + ratio * sin(3.14159265 * eta(iz))
- if(valeta(iz) < zero) valeta(iz) = zero
- if(valeta(iz) > one ) valeta(iz) = one
-! DK DK DK densification sinusoidale ici en vertical
-
- a00(iz) = 1-valeta(iz)
- a01(iz) = valeta(iz)
- enddo
-
- do ix=0,nx
- bot0(ix) = bottom(absx(ix))
- top0(ix) = spl(absx(ix),xtopo,ztopo,coefs_topo,ntopo)
- enddo
-
-! valeurs de x et y pour display domaine physique
- do ix=0,nx
- do iz=0,nz
- x(ix,iz) = absx(ix)
- z(ix,iz) = a00(iz)*bot0(ix) + a01(iz)*top0(ix)
- enddo
- enddo
-
-! calculer min et max de X et Z sur la grille
- print *
- print *, 'Valeurs min et max de X sur le maillage = ',minval(x),maxval(x)
- print *, 'Valeurs min et max de Z sur le maillage = ',minval(z),maxval(z)
- print *
-
-! *** generation de la base de donnees
-
- print *
- print *,' Creation de la base de donnees pour SPECFEM...'
- print *
-
- open(unit=15,file='../SPECFEM90/DataBase',status='unknown')
-
- write(15,*) '#'
- write(15,*) '# Base de Donnees pour Specfem - Premailleur Fortran 90'
- write(15,*) '# ',title
- write(15,*) '# Dimitri Komatitsch, (c) EPS - Harvard August 1998'
- write(15,*) '#'
-
- write(15,*) 'Titre simulation'
- write(15,40) title
-
- npgeo = (nx+1)*(nz+1)
- write(15,*) 'ndofn ndime npgeo'
- write(15,*) ndofn,ndime,npgeo
-
- write(15,*) 'display ignuplot interpol'
- write(15,*) display,ignuplot,interpol
-
- write(15,*) 'itaff itfirstaff icolor inumber'
- write(15,*) itaff,itfirstaff,0,0
-
- write(15,*) 'ivectplot imeshvect imodelvect iboundvect cutvect isubsamp'
- write(15,*) ivectplot,imeshvect,imodelvect,iboundvect,cutvect,isubsamp
-
- write(15,*) 'scalex scalez sizemax angle rapport USletter'
- write(15,*) scalex,scalez,sizemax,20.,0.40,usletter
-
- write(15,*) 'orig_x orig_z isymbols'
- write(15,*) orig_x,orig_z,' T'
-
- write(15,*) 'valseuil freqmaxrep'
- write(15,*) valseuil,freqmaxrep
-
- write(15,*) 'sismos nrec nrec1 nrec2 isamp'
- write(15,*) sismos,nrec,nrec1,nrec2,isamp
-
- write(15,*) 'irepr anglerec anglerec2'
- write(15,*) irepr,anglerec,anglerec2
-
- write(15,*) 'topoplane absstacey compenergy'
- write(15,*) topoplane,absstacey,compenergy
-
- write(15,*) 'initialfield factorana factorxsu n1ana n2ana'
- write(15,*) initialfield,factorana,factorxsu,n1ana,n2ana
-
- write(15,*) 'isismostype ivecttype iaffinfo'
- write(15,*) isismostype,ivecttype,iaffinfo
-
- write(15,*) 'ireadmodel ioutputgrid iavs ivisual3'
- write(15,*) ireadmodel,ioutputgrid,iavs,ivisual3
-
- write(15,*) 'iexec iecho'
- if(iexec) then
- write(15,*) '1 1'
- else
- write(15,*) '0 1'
- endif
-
- write(15,*) 'ncycl dtinc niter'
- write(15,*) nt,dt,niter
-
- write(15,*) 'alpha beta gamma (alpha not used for the moment)'
- write(15,*) alphanewm,betanewm,gammanewm
-
- write(15,*) 'nltfl (number of force or pressure sources)'
- write(15,*) nbsources
-
- write(15,*) 'Collocated forces and/or pressure sources:'
- do i=1,nbsources
- write(15,*) itimetype(i),isource_type(i), &
- xs(i)-xmin ,zs(i), &
- f0(i),t0(i),factor(i),angle(i),0
- enddo
-
- write(15,*) 'Receivers positions:'
- do irec=1,nrec
- write(15,*) irec,xrec(irec)-xmin ,zrec(irec)
- enddo
-
- write(15,*) 'Coordinates of macroblocs mesh (coorg):'
- do j=0,nz
- do i=0,nx
- write(15,*) num(i,j,nx),x(i,j)-xmin,z(i,j)
- enddo
- enddo
-
- netyp = 2
- nxgll = idegpoly + 1
-
- write(15,*) 'netyp numat ngnod nxgll nygll nspec pointsdisp ielemabs ielemperio'
- write(15,*) netyp,nbmodeles,ngnod,nxgll,nxgll,nspec,pointsdisp, &
- nelemabs,nelemperio
-
- write(15,*) 'Material sets (num 0 rho vp vs 0 0)'
- do i=1,nbmodeles
- write(15,*) i,0,rho(i),cp(i),cs(i),0,0
- enddo
-
-
- write(15,*) 'Arrays kmato and knods for each bloc:'
-
- imatnum = 1
- k=0
-
-! zone structuree dans le volume
- do j=0,nz-8,4
- do i=0,nx-4,4
- k = k + 1
- write(15,*) k,imatnum,num(i,j,nx),num(i+4,j,nx),num(i+4,j+4,nx), &
- num(i,j+4,nx),num(i+2,j,nx),num(i+4,j+2,nx), &
- num(i+2,j+4,nx),num(i,j+2,nx),num(i+2,j+2,nx)
- enddo
- enddo
-
- if(k /= nspecvolume) stop 'nombre d''elements incoherent dans le volume'
-
-! zone non structuree dans la couche Wz
- j=nz-4
- do i=0,nx-8,8
-
-! element 1 du raccord
- k = k + 1
- write(15,*) k,imatnum,num(i,j,nx),num(i+4,j,nx),num(i+2,j+2,nx), &
- num(i,j+2,nx),num(i+2,j,nx),num(i+3,j+1,nx), &
- num(i+1,j+2,nx),num(i,j+1,nx),num(i+1,j+1,nx)
-
-! element 2 du raccord
- k = k + 1
- write(15,*) k,imatnum,num(i,j+2,nx),num(i+2,j+2,nx),num(i+2,j+4,nx), &
- num(i,j+4,nx),num(i+1,j+2,nx),num(i+2,j+3,nx), &
- num(i+1,j+4,nx),num(i,j+3,nx),num(i+1,j+3,nx)
-
-! element 3 du raccord
- k = k + 1
- write(15,*) k,imatnum,num(i+2,j+2,nx),num(i+4,j,nx),num(i+4,j+4,nx), &
- num(i+2,j+4,nx),num(i+3,j+1,nx),num(i+4,j+2,nx), &
- num(i+3,j+4,nx),num(i+2,j+3,nx),num(i+3,j+3,nx)
-
-! element 4 du raccord
- k = k + 1
- write(15,*) k,imatnum,num(i+4,j,nx),num(i+6,j+2,nx),num(i+6,j+4,nx), &
- num(i+4,j+4,nx),num(i+5,j+1,nx),num(i+6,j+3,nx), &
- num(i+5,j+4,nx),num(i+4,j+2,nx),num(i+5,j+3,nx)
-
-! element 5 du raccord
- k = k + 1
- write(15,*) k,imatnum,num(i+4,j,nx),num(i+8,j,nx),num(i+8,j+2,nx), &
- num(i+6,j+2,nx),num(i+6,j,nx),num(i+8,j+1,nx), &
- num(i+7,j+2,nx),num(i+5,j+1,nx),num(i+7,j+1,nx)
-
-! element 6 du raccord
- k = k + 1
- write(15,*) k,imatnum,num(i+6,j+2,nx),num(i+8,j+2,nx),num(i+8,j+4,nx), &
- num(i+6,j+4,nx),num(i+7,j+2,nx),num(i+8,j+3,nx), &
- num(i+7,j+4,nx),num(i+6,j+3,nx),num(i+7,j+3,nx)
-
- enddo
-
- if(k /= nspec) stop 'nombre d''elements incoherent dans la couche Wz'
-
-!
-!--- sauvegarde des bords absorbants
-!
-
- print *
- print *,'Au total il y a ',nelemabs,' elements absorbants'
- print *
- print *,'Bords absorbants actifs :'
- print *
- print *,'Haut = ',abshaut
- print *,'Bas = ',absbas
- print *,'Gauche = ',absgauche
- print *,'Droite = ',absdroite
- print *
- print *,'Stacey = ',absstacey
- print *
-
-! generer la liste des elements absorbants
- if(nelemabs > 0) then
- write(15,*) 'Liste des elements absorbants (haut bas gauche droite) :'
-
-! repasser aux vrais valeurs de nx et nz
- nx = nx / 4
- nz = nz / 4
-
- inumabs = 0
-
-! bord absorbant du bas sans les coins
- iz = 1
- do ix = 2,nx-1
- inumabs = inumabs + 1
- inumelem = (iz-1)*nx + ix
- icodehaut = 0
- icodebas = iaretebas
- icodegauche = 0
- icodedroite = 0
- write(15,*) inumabs,inumelem,icodehaut,icodebas,icodegauche,icodedroite
- enddo
-
-! coin en bas a gauche
- inumabs = inumabs + 1
- inumelem = 1
- icodehaut = 0
- icodebas = iaretebas
- icodegauche = iaretegauche
- icodedroite = 0
- write(15,*) inumabs,inumelem,icodehaut,icodebas,icodegauche,icodedroite
-
-! coin en bas a droite
- inumabs = inumabs + 1
- inumelem = nx
- icodehaut = 0
- icodebas = iaretebas
- icodegauche = 0
- icodedroite = iaretedroite
- write(15,*) inumabs,inumelem,icodehaut,icodebas,icodegauche,icodedroite
-
-! partie structuree du bord de gauche
- ix = 1
- do iz = 2,nz-1
- inumabs = inumabs + 1
- inumelem = (iz-1)*nx + ix
- icodehaut = 0
- icodebas = 0
- icodegauche = iaretegauche
- icodedroite = 0
- write(15,*) inumabs,inumelem,icodehaut,icodebas,icodegauche,icodedroite
- enddo
-
-! partie structuree du bord de droite
- ix = nx
- do iz = 2,nz-1
- inumabs = inumabs + 1
- inumelem = (iz-1)*nx + ix
- icodehaut = 0
- icodebas = 0
- icodegauche = 0
- icodedroite = iaretedroite
- write(15,*) inumabs,inumelem,icodehaut,icodebas,icodegauche,icodedroite
- enddo
-
-! partie non structuree du bord de gauche (deux elements)
- do iadd = 1,2
- inumabs = inumabs + 1
- inumelem = nx*(nz-1) + iadd
- icodehaut = 0
- icodebas = 0
- icodegauche = iaretegauche
- icodedroite = 0
- write(15,*) inumabs,inumelem,icodehaut,icodebas,icodegauche,icodedroite
- enddo
-
-! partie non structuree du bord de droite (deux elements)
- do iadd = 1,2
- inumabs = inumabs + 1
- inumelem = nspec - iadd + 1
- icodehaut = 0
- icodebas = 0
- icodegauche = 0
- icodedroite = iaretedroite
- write(15,*) inumabs,inumelem,icodehaut,icodebas,icodegauche,icodedroite
- enddo
-
- if(inumabs /= nelemabs) stop 'nombre d''elements absorbants incoherent'
-
- endif
-
-! fermer la base de donnees
-
- close(15)
-
- 40 format(a50)
-
- end program maille_non_struct_2
-
-! *****************
-! routines maillage
-! *****************
-
-! --- numero global du noeud
-
- integer function num(i,j,nx)
- implicit none
- integer i,j,nx
-
- num = j*(nx+1) + i + 1
-
- end function num
-
-! ------- definition des fonctions representant les interfaces -------
-
-!
-! --- bas du modele
-!
-
- double precision function bottom(x)
- implicit none
- double precision x
-
- bottom = 0.d0
-
- end function bottom
-
-!
-! --- representation interfaces par un spline
-!
-
-!--- spline
- double precision function spl(x,xtopo,ztopo,coefs,ntopo)
- implicit none
- integer ntopo
- double precision x,xp
- double precision xtopo(ntopo),ztopo(ntopo)
- double precision coefs(ntopo)
-
- spl = 0.
- xp = x
- if (xp < xtopo(1)) xp = xtopo(1)
- if (xp > xtopo(ntopo)) xp = xtopo(ntopo)
- call splint(xtopo,ztopo,coefs,ntopo,xp,spl)
-
- end function spl
-
-! --- fonction de densification du maillage horizontal
-
- double precision function dens(ix,psi,xmin,xmax,nx)
- implicit none
- integer ix,nx
- double precision psi(0:nx)
- double precision xmin,xmax
-
- dens = xmin + dble(xmax-xmin)*psi(ix)
-
- end function dens
-
-! --------------------------------------
-
-! routine de calcul des coefs du spline (Numerical Recipes)
- subroutine spline(x,y,n,yp1,ypn,y2)
- implicit none
-
- integer n
- double precision x(n),y(n),y2(n)
- double precision, dimension(:), allocatable :: u
- double precision yp1,ypn
-
- integer i,k
- double precision sig,p,qn,un
-
- allocate(u(n))
-
- y2(1)=-0.5
- u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
- do i=2,n-1
- sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
- p=sig*y2(i-1)+2.
- y2(i)=(sig-1.)/p
- u(i)=(6.*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1)) &
- /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
- enddo
- qn=0.5
- un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
- y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
- do k=n-1,1,-1
- y2(k)=y2(k)*y2(k+1)+u(k)
- enddo
-
- deallocate(u)
-
- end subroutine spline
-
-! --------------
-
-! routine d'evaluation du spline (Numerical Recipes)
- SUBROUTINE SPLINT(XA,YA,Y2A,N,X,Y)
- implicit none
-
- integer n
- double precision XA(N),YA(N),Y2A(N)
- double precision x,y
-
- integer k,klo,khi
- double precision h,a,b
-
- KLO=1
- KHI=N
- do while (KHI-KLO > 1)
- K=(KHI+KLO)/2
- IF(XA(K) > X)THEN
- KHI=K
- ELSE
- KLO=K
- ENDIF
- enddo
- H=XA(KHI)-XA(KLO)
- IF (H == 0.d0) stop 'Bad input in spline evaluation'
- A=(XA(KHI)-X)/H
- B=(X-XA(KLO))/H
-
- Y=A*YA(KLO)+B*YA(KHI)+((A**3-A)*Y2A(KLO)+ &
- (B**3-B)*Y2A(KHI))*(H**2)/6.d0
-
- end subroutine SPLINT
-
Deleted: seismo/2D/SPECFEM2D/trunk/maille_non_struct_3.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/maille_non_struct_3.f90 2007-02-15 00:06:00 UTC (rev 8501)
+++ seismo/2D/SPECFEM2D/trunk/maille_non_struct_3.f90 2007-12-07 23:51:25 UTC (rev 8502)
@@ -1,985 +0,0 @@
-!=====================================================================
-!
-! P r e m a i l l e u r F o r t r a n 9 0
-! -------------------------------------------
-!
-! Version 3.0
-! -----------
-!
-! Dimitri Komatitsch
-! Department of Earth and Planetary Sciences - Harvard University
-!
-! (c) August 1998
-!
-!=====================================================================
-
-!
-! *** Version optimisee avec maillage non structure Jacques Muller - Elf ***
-! *** Raffinement d'un facteur 3 en surface ***
-!
-
- program maille_non_struct_3
-
- implicit none
-
-! definir les tableaux pour allocation dynamique
-
-! coordinates of the grid points
- double precision, allocatable :: x(:,:),z(:,:)
-
-! variables needed to compute the transformation
- double precision, allocatable :: psi(:),eta(:),absx(:), &
- a00(:),a01(:),valeta(:),bot0(:),top0(:)
-
-! stockage du modele de vitesse et densite
- double precision, allocatable :: rho(:),cp(:),cs(:)
-
-! the topography data
- double precision, allocatable :: xtopo(:),ztopo(:),coefs_topo(:)
-
-! arrays for the source
- double precision, allocatable :: xs(:),zs(:),f0(:),t0(:),angle(:),factor(:)
- integer, allocatable :: isource_type(:),itimetype(:)
-
-! arrays for the receivers
- double precision, allocatable :: xrec(:),zrec(:)
-
- character(len=50) interffile,topofile,title
- character(len=15) junk
-
- integer imatnum,inumabs,inumelem,nelemperio,nxgll,netyp
- integer icodehaut,icodebas,icodegauche,icodedroite
- integer nelemabs,npgeo,nspec,ntopo,nspecvolume,nspecWz
- integer k,ix,iz,irec,i,j,iadd
- integer imodele,nbmodeles,iaffinfo
- integer itaff,itfirstaff,pointsdisp,isubsamp,nrec,n1ana,n2ana
- integer irepr,nrec1,nrec2,isamp,nbsources,isismostype,ivecttype
- integer ngnod,nt,niter,idegpoly,nx,nz
- integer icodematread
-
- double precision valseuil,freqmaxrep,ratio
- double precision tang1,tangN
- double precision orig_x,orig_z,sizemax,cutvect,scalex,scalez
- double precision factorxsu,factorana,xspacerec,zspacerec
- double precision anglerec,anglerec2,xmin,xmax
- double precision xfin,zfin,xfin2,zfin2,xdeb,zdeb,xdeb2,zdeb2
- double precision alphanewm,betanewm,gammanewm,dt
- double precision rhoread,cpread,csread,aniso3read,aniso4read
-
- logical interpol,ignuplot,ireadmodel,iavs,ivisual3,ioutputgrid
- logical abshaut,absbas,absgauche,absdroite,absstacey
- logical periohaut,periogauche
- logical sismos,isources_surf,ienreg_surf,ienreg_surf2,display
- logical ivectplot,imeshvect
- logical topoplane,iexec,initialfield
- logical imodelvect,iboundvect,usletter,compenergy
-
- integer, external :: num
- double precision, external :: bottom,spl,dens
-
- double precision, parameter :: zero = 0.d0, one = 1.d0
-
-! simulation a 2D
- integer, parameter :: ndime = 2
- integer, parameter :: ndofn = 2
-
-! --- code des numeros d'aretes pour les bords absorbants
- integer, parameter :: iaretebas = 1
- integer, parameter :: iaretedroite = 2
- integer, parameter :: iaretehaut = 3
- integer, parameter :: iaretegauche = 4
-
-! DK DK DK ajout Elf : extraction de la topo du fichier SEP
-!! call system('rm -f topo_from_SEP.dat topo_SEP_maille90.dat ; xextract_topo')
-
- print *
- print *,' *** Version optimisee avec maillage non structure ***'
- print *,' *** Raffinement d''un facteur 3 en surface ***'
- print *
-
-! ***
-! *** read the parameter file
-! ***
-
- print *,' Reading the parameter file ... '
- print *
-
- open(unit=10,file='DATA/Par_file',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,topofile
- read(10,3)junk,interffile
-
- 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,nz
- read(10,2)junk,idegpoly
- read(10,2)junk,ngnod
- read(10,1)junk,ratio
- read(10,4)junk,topoplane
- read(10,4)junk,initialfield
- read(10,4)junk,ireadmodel
- read(10,4)junk,iexec
-
-! DK DK forcer pour Elf
- ngnod = 9
- topoplane = .false.
- initialfield = .false.
-
-! pour le non structure, verifier la coherence du maillage
- if(nx < 2) stop 'nx must be greater or equal to 2'
- if(nz < 2) stop 'nz must be greater or equal to 2'
- if(mod(nx,2) /= 0) stop 'nx must be even'
-
-! multiplier par 3 pour implementer le deraffinement non conforme
- nx = nx * 3
- nz = nz * 3
-
-! multiplier par 2 pour elements 9 noeuds
- nx = nx * 2
- nz = nz * 2
-
-! 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
- read(10,4)junk,absstacey
- read(10,4)junk,periohaut
- read(10,4)junk,periogauche
-
-! DK DK forcer pour Elf
- abshaut = .false.
- periohaut = .false.
- periogauche = .false.
-
-! skip comment
- read(10,*)
- read(10,*)
- read(10,*)
-
-! read time step parameters
- read(10,2)junk,nt
- read(10,1)junk,dt
- read(10,2)junk,niter
- read(10,1)junk,alphanewm
- read(10,1)junk,betanewm
- read(10,1)junk,gammanewm
-
-! skip comment
- read(10,*)
- read(10,*)
- read(10,*)
-
-! read source parameters
- read(10,2)junk,nbsources
- read(10,4)junk,isources_surf
- read(10,1)junk,valseuil
- read(10,1)junk,freqmaxrep
- print *,'Nb de sources a lire : ',nbsources
-
- allocate(xs(nbsources))
- allocate(zs(nbsources))
- allocate(f0(nbsources))
- allocate(t0(nbsources))
- allocate(isource_type(nbsources))
- allocate(itimetype(nbsources))
- allocate(angle(nbsources))
- allocate(factor(nbsources))
-
- do i=1,nbsources
- read(10,*)
- read(10,1)junk,xs(i)
- read(10,1)junk,zs(i)
- read(10,1)junk,f0(i)
- read(10,1)junk,t0(i)
- read(10,2)junk,isource_type(i)
- read(10,2)junk,itimetype(i)
- read(10,1)junk,angle(i)
- read(10,1)junk,factor(i)
-
- print *
- print *,' Source #',i
- print *,'Position xs, zs = ',xs(i),zs(i)
- print *,'Frequency, delay = ',f0(i),t0(i)
- print *,'Source type (1=force 2=explo) : ', &
- isource_type(i)
- print *,'Angle of the source if force = ',angle(i)
- print *,'Multiplying factor = ',factor(i)
- enddo
-
-! skip comment
- read(10,*)
- read(10,*)
- read(10,*)
-
-! read receivers line parameters
- read(10,4)junk,sismos
- read(10,2)junk,isamp
- read(10,2)junk,isismostype
- read(10,2)junk,irepr
- read(10,*)
- read(10,2)junk,nrec1
- read(10,1)junk,xdeb
- read(10,1)junk,zdeb
- read(10,1)junk,xfin
- read(10,1)junk,zfin
- read(10,4)junk,ienreg_surf
- read(10,1)junk,anglerec
- read(10,*)
- read(10,2)junk,nrec2
- read(10,1)junk,xdeb2
- read(10,1)junk,zdeb2
- read(10,1)junk,xfin2
- read(10,1)junk,zfin2
- read(10,4)junk,ienreg_surf2
- read(10,1)junk,anglerec2
- read(10,*)
- read(10,1)junk,factorxsu
- read(10,2)junk,n1ana
- read(10,2)junk,n2ana
- read(10,1)junk,factorana
-
-! determination et affichage position ligne de receivers
- if(nrec2 < 0) stop 'negative value of nrec2 !'
-
- if(nrec2 == 0) then
- nrec = nrec1
- else
- nrec = nrec1 + nrec2
- endif
-
-! DK DK forcer pour Elf
- n1ana = 1
- n2ana = nrec
-
- allocate(xrec(nrec))
- allocate(zrec(nrec))
-
- if(nrec2 == 0) then
- print *
- print *,'There are ',nrec,' receivers on a single line'
- xspacerec=(xfin-xdeb)/dble(nrec-1)
- zspacerec=(zfin-zdeb)/dble(nrec-1)
- do i=1,nrec
- xrec(i) = xdeb + dble(i-1)*xspacerec
- zrec(i) = zdeb + dble(i-1)*zspacerec
- enddo
- else
- print *
- print *,'There are ',nrec,' receivers on two lines'
- print *,'First line contains ',nrec1,' receivers'
- print *,'Second line contains ',nrec2,' receivers'
- xspacerec=(xfin-xdeb)/dble(nrec1-1)
- zspacerec=(zfin-zdeb)/dble(nrec1-1)
- do i=1,nrec1
- xrec(i) = xdeb + dble(i-1)*xspacerec
- zrec(i) = zdeb + dble(i-1)*zspacerec
- enddo
- xspacerec=(xfin2-xdeb2)/dble(nrec2-1)
- zspacerec=(zfin2-zdeb2)/dble(nrec2-1)
- do i=1,nrec2
- xrec(i+nrec1) = xdeb2 + dble(i-1)*xspacerec
- zrec(i+nrec1) = zdeb2 + dble(i-1)*zspacerec
- enddo
- endif
-
-! skip comment
- read(10,*)
- read(10,*)
- read(10,*)
-
-! read display parameters
- read(10,4)junk,display
- read(10,2)junk,itaff
- read(10,2)junk,itfirstaff
- read(10,2)junk,iaffinfo
- read(10,4)junk,ivectplot
- 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,pointsdisp
- read(10,2)junk,isubsamp
- read(10,1)junk,scalex
- read(10,1)junk,scalez
- read(10,1)junk,sizemax
- read(10,4)junk,usletter
- read(10,1)junk,orig_x
- read(10,1)junk,orig_z
- read(10,4)junk,ignuplot
- read(10,4)junk,iavs
- read(10,4)junk,ivisual3
- read(10,4)junk,ioutputgrid
- read(10,4)junk,compenergy
-
-! DK DK forcer pour Elf
- ignuplot = .false.
- iavs = .false.
- ivisual3 = .false.
- compenergy = .false.
-
-! skip comment
- read(10,*)
- read(10,*)
- read(10,*)
-
-! lecture des differents modeles de materiaux
-
- read(10,2)junk,nbmodeles
- if(nbmodeles <= 0) stop 'Negative number of models not allowed !!'
-
- allocate(rho(nbmodeles))
- allocate(cp(nbmodeles))
- allocate(cs(nbmodeles))
-
- rho(:) = 0.d0
- cp(:) = 0.d0
- cs(:) = 0.d0
-
- do imodele=1,nbmodeles
- read(10,*) i,icodematread,rhoread,cpread,csread,aniso3read,aniso4read
- if(i<1 .or. i>nbmodeles) stop 'Wrong material set number'
- 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'
- enddo
-
- print *
- print *, 'Nb de modeles de roche = ',nbmodeles
- print *
- do i=1,nbmodeles
- print *,'Modele #',i,' isotrope'
- print *,'rho,cp,cs = ',rho(i),cp(i),cs(i)
- enddo
- print *
-
- close(10)
-
- print *
- print *,' Parameter file successfully read... '
-
-! --------- fin lecture fichier parametres --------------
-
- allocate(psi(0:nx))
- allocate(eta(0:nz))
- allocate(absx(0:nx))
- allocate(a00(0:nz))
- allocate(a01(0:nz))
- allocate(valeta(0:nz))
- allocate(bot0(0:nx))
- allocate(top0(0:nx))
-
-! calcul des points regulierement espaces
- do i=0,nx
- psi(i) = i/dble(nx)
- enddo
- do j=0,nz
- eta(j) = j/dble(nz)
- enddo
-
-! quelques verifications de base a faire
-
- if(ngnod /= 9) stop 'erreur ngnod different de 9 !!'
-
-! calcul du nombre total d'elements spectraux, absorbants et periodiques
- nspecvolume = (nx/2/3)*((nz-6)/2/3)
- nspecWz = 5*(nx/2/3)
- nspec = nspecvolume + nspecWz
- nelemperio = 0
-
- if(absgauche .or. absdroite .or. absbas) then
- nelemabs = 2 * (nz/6 - 2) + nx/6 + 3 + 3
- else
- nelemabs = 0
- endif
-
- print *
- print *,'Le maillage comporte ',nspec,' elements spectraux (nx = ',nx/6, &
- ' nz = ',nz/6,')'
- print *,'soit ',nspecvolume,' elements spectraux dans le volume'
- print *,'et ',nspecWz,' elements spectraux dans la couche Wz'
- print *,'Chaque element comporte ',idegpoly+1,' points dans chaque direction'
- print *,'Le nombre maximum de points theorique est ',nspec*(idegpoly+1)**ndime
- print *,'Les elements de controle sont des elements ',ngnod,' noeuds'
- print *
-
-!------------------------------------------------------
-
- allocate(x(0:nx,0:nz))
- allocate(z(0:nx,0:nz))
-
- x(:,:)=0.d0
- z(:,:)=0.d0
-
-! get topography data from external file
- print *,'Reading topography from file ',topofile
- open(unit=15,file=topofile,status='old')
- read(15,*) ntopo
- if (ntopo < 2) stop 'Not enough topography points (min 2)'
- print *,'Reading ',ntopo,' points from topography file'
- print *
-
- allocate(xtopo(ntopo))
- allocate(ztopo(ntopo))
- allocate(coefs_topo(ntopo))
-
- do i=1,ntopo
- read(15,*) xtopo(i),ztopo(i)
- enddo
- close(15)
-
-! check the values read
- print *
- print *, 'Topography data points (x,z)'
- print *, '----------------------------'
- print *
- print *, 'Topo 1 = (',xtopo(1),',',ztopo(1),')'
- print *, 'Topo ntopo = (',xtopo(ntopo),',',ztopo(ntopo),')'
-
-!--- calculate the spline function for the topography
-!--- imposer les tangentes aux deux bords
- tang1 = (ztopo(2)-ztopo(1))/(xtopo(2)-xtopo(1))
- tangN = (ztopo(ntopo)-ztopo(ntopo-1))/(xtopo(ntopo)-xtopo(ntopo-1))
- call spline(xtopo,ztopo,ntopo,tang1,tangN,coefs_topo)
-
-! *** afficher limites du modele lu
- print *
- print *, 'Limites absolues modele fichier topo :'
- print *
- print *, 'Xmin = ',minval(xtopo),' Xmax = ',maxval(xtopo)
- print *, 'Zmin = ',minval(ztopo),' Zmax = ',maxval(ztopo)
- print *
-
-! *** modifier sources pour position par rapport a la surface
- print *
- print *, 'Position (x,z) des ',nbsources,' sources'
- print *
- do i=1,nbsources
-
-! DK DK DK Elf : position source donnee en profondeur par rapport a la topo
- zs(i) = spl(xs(i),xtopo,ztopo,coefs_topo,ntopo) - zs(i)
-
- if(isources_surf) zs(i) = spl(xs(i),xtopo,ztopo,coefs_topo,ntopo)
- print *, 'Source ',i,' = ',xs(i),zs(i)
- enddo
-
-! *** modifier recepteurs pour enregistrement en surface
- print *
- print *, 'Position (x,z) des ',nrec,' receivers'
- print *
- do irec=1,nrec
-
-! DK DK DK Elf : distinguer les deux lignes de recepteurs
- if(irec <= nrec1) then
- if(ienreg_surf) zrec(irec) = spl(xrec(irec),xtopo,ztopo,coefs_topo,ntopo)
- else
- if(ienreg_surf2) zrec(irec) = spl(xrec(irec),xtopo,ztopo,coefs_topo,ntopo)
- endif
- print *, 'Receiver ',irec,' = ',xrec(irec),zrec(irec)
-
- enddo
-
-!--- definition du maillage suivant X
- do ix=0,nx
- absx(ix) = dens(ix,psi,xmin,xmax,nx)
- enddo
-
-! *** une seule zone
-
- do iz=0,nz
-
-! DK DK DK densification sinusoidale ici en vertical
- valeta(iz) = eta(iz) + ratio * sin(3.14159265 * eta(iz))
- if(valeta(iz) < zero) valeta(iz) = zero
- if(valeta(iz) > one ) valeta(iz) = one
-! DK DK DK densification sinusoidale ici en vertical
-
- a00(iz) = 1-valeta(iz)
- a01(iz) = valeta(iz)
- enddo
-
- do ix=0,nx
- bot0(ix) = bottom(absx(ix))
- top0(ix) = spl(absx(ix),xtopo,ztopo,coefs_topo,ntopo)
- enddo
-
-! valeurs de x et y pour display domaine physique
- do ix=0,nx
- do iz=0,nz
- x(ix,iz) = absx(ix)
- z(ix,iz) = a00(iz)*bot0(ix) + a01(iz)*top0(ix)
- enddo
- enddo
-
-! calculer min et max de X et Z sur la grille
- print *
- print *, 'Valeurs min et max de X sur le maillage = ',minval(x),maxval(x)
- print *, 'Valeurs min et max de Z sur le maillage = ',minval(z),maxval(z)
- print *
-
-! *** generation de la base de donnees
-
- print *
- print *,' Creation de la base de donnees pour SPECFEM...'
- print *
-
- open(unit=15,file='../SPECFEM90/DataBase',status='unknown')
-
- write(15,*) '#'
- write(15,*) '# Base de Donnees pour Specfem - Premailleur Fortran 90'
- write(15,*) '# ',title
- write(15,*) '# Dimitri Komatitsch, (c) EPS - Harvard August 1998'
- write(15,*) '#'
-
- write(15,*) 'Titre simulation'
- write(15,40) title
-
- npgeo = (nx+1)*(nz+1)
- write(15,*) 'ndofn ndime npgeo'
- write(15,*) ndofn,ndime,npgeo
-
- write(15,*) 'display ignuplot interpol'
- write(15,*) display,ignuplot,interpol
-
- write(15,*) 'itaff itfirstaff icolor inumber'
- write(15,*) itaff,itfirstaff,0,0
-
- write(15,*) 'ivectplot imeshvect imodelvect iboundvect cutvect isubsamp'
- write(15,*) ivectplot,imeshvect,imodelvect,iboundvect,cutvect,isubsamp
-
- write(15,*) 'scalex scalez sizemax angle rapport USletter'
- write(15,*) scalex,scalez,sizemax,20.,0.40,usletter
-
- write(15,*) 'orig_x orig_z isymbols'
- write(15,*) orig_x,orig_z,' T'
-
- write(15,*) 'valseuil freqmaxrep'
- write(15,*) valseuil,freqmaxrep
-
- write(15,*) 'sismos nrec nrec1 nrec2 isamp'
- write(15,*) sismos,nrec,nrec1,nrec2,isamp
-
- write(15,*) 'irepr anglerec anglerec2'
- write(15,*) irepr,anglerec,anglerec2
-
- write(15,*) 'topoplane absstacey compenergy'
- write(15,*) topoplane,absstacey,compenergy
-
- write(15,*) 'initialfield factorana factorxsu n1ana n2ana'
- write(15,*) initialfield,factorana,factorxsu,n1ana,n2ana
-
- write(15,*) 'isismostype ivecttype iaffinfo'
- write(15,*) isismostype,ivecttype,iaffinfo
-
- write(15,*) 'ireadmodel ioutputgrid iavs ivisual3'
- write(15,*) ireadmodel,ioutputgrid,iavs,ivisual3
-
- write(15,*) 'iexec iecho'
- if(iexec) then
- write(15,*) '1 1'
- else
- write(15,*) '0 1'
- endif
-
- write(15,*) 'ncycl dtinc niter'
- write(15,*) nt,dt,niter
-
- write(15,*) 'alpha beta gamma (alpha not used for the moment)'
- write(15,*) alphanewm,betanewm,gammanewm
-
- write(15,*) 'nltfl (number of force or pressure sources)'
- write(15,*) nbsources
-
- write(15,*) 'Collocated forces and/or pressure sources:'
- do i=1,nbsources
- write(15,*) itimetype(i),isource_type(i), &
- xs(i)-xmin ,zs(i), &
- f0(i),t0(i),factor(i),angle(i),0
- enddo
-
- write(15,*) 'Receivers positions:'
- do irec=1,nrec
- write(15,*) irec,xrec(irec)-xmin ,zrec(irec)
- enddo
-
- write(15,*) 'Coordinates of macroblocs mesh (coorg):'
- do j=0,nz
- do i=0,nx
- write(15,*) num(i,j,nx),x(i,j)-xmin,z(i,j)
- enddo
- enddo
-
- netyp = 2
- nxgll = idegpoly + 1
-
- write(15,*) 'netyp numat ngnod nxgll nygll nspec pointsdisp ielemabs ielemperio'
- write(15,*) netyp,nbmodeles,ngnod,nxgll,nxgll,nspec,pointsdisp, &
- nelemabs,nelemperio
-
- write(15,*) 'Material sets (num 0 rho vp vs 0 0)'
- do i=1,nbmodeles
- write(15,*) i,0,rho(i),cp(i),cs(i),0,0
- enddo
-
-
- write(15,*) 'Arrays kmato and knods for each bloc:'
-
- imatnum = 1
- k=0
-
-! zone structuree dans le volume
- do j=0,nz-12,6
- do i=0,nx-6,6
- k = k + 1
- write(15,*) k,imatnum,num(i,j,nx),num(i+6,j,nx),num(i+6,j+6,nx), &
- num(i,j+6,nx),num(i+3,j,nx),num(i+6,j+3,nx), &
- num(i+3,j+6,nx),num(i,j+3,nx),num(i+3,j+3,nx)
- enddo
- enddo
-
- if(k /= nspecvolume) stop 'nombre d''elements incoherent dans le volume'
-
-! zone non structuree dans la couche Wz
- j=nz-6
- do i=0,nx-12,12
-
-! element 1 du raccord
- k = k + 1
- write(15,*) k,imatnum,num(i,j,nx),num(i+6,j,nx),num(i+4,j+2,nx), &
- num(i,j+2,nx),num(i+3,j,nx),num(i+5,j+1,nx), &
- num(i+2,j+2,nx),num(i,j+1,nx),num(i+3,j+1,nx)
-
-! element 2 du raccord
- k = k + 1
- write(15,*) k,imatnum,num(i,j+2,nx),num(i+4,j+2,nx),num(i+2,j+4,nx), &
- num(i,j+4,nx),num(i+2,j+2,nx),num(i+3,j+3,nx), &
- num(i+1,j+4,nx),num(i,j+3,nx),num(i+1,j+3,nx)
-
-! element 3 du raccord
- k = k + 1
- write(15,*) k,imatnum,num(i,j+4,nx),num(i+2,j+4,nx),num(i+2,j+6,nx), &
- num(i,j+6,nx),num(i+1,j+4,nx),num(i+2,j+5,nx), &
- num(i+1,j+6,nx),num(i,j+5,nx),num(i+1,j+5,nx)
-
-! element 4 du raccord
- k = k + 1
- write(15,*) k,imatnum,num(i+2,j+4,nx),num(i+4,j+2,nx),num(i+4,j+6,nx), &
- num(i+2,j+6,nx),num(i+3,j+3,nx),num(i+4,j+4,nx), &
- num(i+3,j+6,nx),num(i+2,j+5,nx),num(i+3,j+5,nx)
-
-! element 5 du raccord
- k = k + 1
- write(15,*) k,imatnum,num(i+4,j+2,nx),num(i+6,j,nx),num(i+6,j+6,nx), &
- num(i+4,j+6,nx),num(i+5,j+1,nx),num(i+6,j+3,nx), &
- num(i+5,j+6,nx),num(i+4,j+4,nx),num(i+5,j+3,nx)
-
-! element 6 du raccord
- k = k + 1
- write(15,*) k,imatnum,num(i+6,j,nx),num(i+8,j+2,nx),num(i+8,j+6,nx), &
- num(i+6,j+6,nx),num(i+7,j+1,nx),num(i+8,j+4,nx), &
- num(i+7,j+6,nx),num(i+6,j+3,nx),num(i+7,j+3,nx)
-
-! element 7 du raccord
- k = k + 1
- write(15,*) k,imatnum,num(i+8,j+2,nx),num(i+10,j+4,nx),num(i+10,j+6,nx), &
- num(i+8,j+6,nx),num(i+9,j+3,nx),num(i+10,j+5,nx), &
- num(i+9,j+6,nx),num(i+8,j+4,nx),num(i+9,j+4,nx)
-
-! element 8 du raccord
- k = k + 1
- write(15,*) k,imatnum,num(i+6,j,nx),num(i+12,j,nx),num(i+12,j+2,nx), &
- num(i+8,j+2,nx),num(i+9,j,nx),num(i+12,j+1,nx), &
- num(i+10,j+2,nx),num(i+7,j+1,nx),num(i+10,j+1,nx)
-
-! element 9 du raccord
- k = k + 1
- write(15,*) k,imatnum,num(i+8,j+2,nx),num(i+12,j+2,nx),num(i+12,j+4,nx), &
- num(i+10,j+4,nx),num(i+10,j+2,nx),num(i+12,j+3,nx), &
- num(i+11,j+4,nx),num(i+9,j+3,nx),num(i+11,j+3,nx)
-
-! element 10 du raccord
- k = k + 1
- write(15,*) k,imatnum,num(i+10,j+4,nx),num(i+12,j+4,nx),num(i+12,j+6,nx),&
- num(i+10,j+6,nx),num(i+11,j+4,nx),num(i+12,j+5,nx), &
- num(i+11,j+6,nx),num(i+10,j+5,nx),num(i+11,j+5,nx)
-
- enddo
-
- if(k /= nspec) stop 'nombre d''elements incoherent dans la couche Wz'
-
-!
-!--- sauvegarde des bords absorbants
-!
-
- print *
- print *,'Au total il y a ',nelemabs,' elements absorbants'
- print *
- print *,'Bords absorbants actifs :'
- print *
- print *,'Haut = ',abshaut
- print *,'Bas = ',absbas
- print *,'Gauche = ',absgauche
- print *,'Droite = ',absdroite
- print *
- print *,'Stacey = ',absstacey
- print *
-
-! generer la liste des elements absorbants
- if(nelemabs > 0) then
- write(15,*) 'Liste des elements absorbants (haut bas gauche droite) :'
-
-! repasser aux vrais valeurs de nx et nz
- nx = nx / 6
- nz = nz / 6
-
- inumabs = 0
-
-! bord absorbant du bas sans les coins
- iz = 1
- do ix = 2,nx-1
- inumabs = inumabs + 1
- inumelem = (iz-1)*nx + ix
- icodehaut = 0
- icodebas = iaretebas
- icodegauche = 0
- icodedroite = 0
- write(15,*) inumabs,inumelem,icodehaut,icodebas,icodegauche,icodedroite
- enddo
-
-! coin en bas a gauche
- inumabs = inumabs + 1
- inumelem = 1
- icodehaut = 0
- icodebas = iaretebas
- icodegauche = iaretegauche
- icodedroite = 0
- write(15,*) inumabs,inumelem,icodehaut,icodebas,icodegauche,icodedroite
-
-! coin en bas a droite
- inumabs = inumabs + 1
- inumelem = nx
- icodehaut = 0
- icodebas = iaretebas
- icodegauche = 0
- icodedroite = iaretedroite
- write(15,*) inumabs,inumelem,icodehaut,icodebas,icodegauche,icodedroite
-
-! partie structuree du bord de gauche
- ix = 1
- do iz = 2,nz-1
- inumabs = inumabs + 1
- inumelem = (iz-1)*nx + ix
- icodehaut = 0
- icodebas = 0
- icodegauche = iaretegauche
- icodedroite = 0
- write(15,*) inumabs,inumelem,icodehaut,icodebas,icodegauche,icodedroite
- enddo
-
-! partie structuree du bord de droite
- ix = nx
- do iz = 2,nz-1
- inumabs = inumabs + 1
- inumelem = (iz-1)*nx + ix
- icodehaut = 0
- icodebas = 0
- icodegauche = 0
- icodedroite = iaretedroite
- write(15,*) inumabs,inumelem,icodehaut,icodebas,icodegauche,icodedroite
- enddo
-
-! partie non structuree du bord de gauche (trois elements)
- do iadd = 1,3
- inumabs = inumabs + 1
- inumelem = nx*(nz-1) + iadd
- icodehaut = 0
- icodebas = 0
- icodegauche = iaretegauche
- icodedroite = 0
- write(15,*) inumabs,inumelem,icodehaut,icodebas,icodegauche,icodedroite
- enddo
-
-! partie non structuree du bord de droite (trois elements)
- do iadd = 1,3
- inumabs = inumabs + 1
- inumelem = nspec - iadd + 1
- icodehaut = 0
- icodebas = 0
- icodegauche = 0
- icodedroite = iaretedroite
- write(15,*) inumabs,inumelem,icodehaut,icodebas,icodegauche,icodedroite
- enddo
-
- if(inumabs /= nelemabs) stop 'nombre d''elements absorbants incoherent'
-
- endif
-
-! fermer la base de donnees
-
- close(15)
-
- 40 format(a50)
-
- end program maille_non_struct_3
-
-! *****************
-! routines maillage
-! *****************
-
-! --- numero global du noeud
-
- integer function num(i,j,nx)
- implicit none
- integer i,j,nx
-
- num = j*(nx+1) + i + 1
-
- end function num
-
-! ------- definition des fonctions representant les interfaces -------
-
-!
-! --- bas du modele
-!
-
- double precision function bottom(x)
- implicit none
- double precision x
-
- bottom = 0.d0
-
- end function bottom
-
-!
-! --- representation interfaces par un spline
-!
-
-!--- spline
- double precision function spl(x,xtopo,ztopo,coefs,ntopo)
- implicit none
- integer ntopo
- double precision x,xp
- double precision xtopo(ntopo),ztopo(ntopo)
- double precision coefs(ntopo)
-
- spl = 0.
- xp = x
- if (xp < xtopo(1)) xp = xtopo(1)
- if (xp > xtopo(ntopo)) xp = xtopo(ntopo)
- call splint(xtopo,ztopo,coefs,ntopo,xp,spl)
-
- end function spl
-
-! --- fonction de densification du maillage horizontal
-
- double precision function dens(ix,psi,xmin,xmax,nx)
- implicit none
- integer ix,nx
- double precision psi(0:nx)
- double precision xmin,xmax
-
- dens = xmin + dble(xmax-xmin)*psi(ix)
-
- end function dens
-
-! --------------------------------------
-
-! routine de calcul des coefs du spline (Numerical Recipes)
- subroutine spline(x,y,n,yp1,ypn,y2)
- implicit none
-
- integer n
- double precision x(n),y(n),y2(n)
- double precision, dimension(:), allocatable :: u
- double precision yp1,ypn
-
- integer i,k
- double precision sig,p,qn,un
-
- allocate(u(n))
-
- y2(1)=-0.5
- u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
- do i=2,n-1
- sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
- p=sig*y2(i-1)+2.
- y2(i)=(sig-1.)/p
- u(i)=(6.*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1)) &
- /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
- enddo
- qn=0.5
- un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
- y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
- do k=n-1,1,-1
- y2(k)=y2(k)*y2(k+1)+u(k)
- enddo
-
- deallocate(u)
-
- end subroutine spline
-
-! --------------
-
-! routine d'evaluation du spline (Numerical Recipes)
- SUBROUTINE SPLINT(XA,YA,Y2A,N,X,Y)
- implicit none
-
- integer n
- double precision XA(N),YA(N),Y2A(N)
- double precision x,y
-
- integer k,klo,khi
- double precision h,a,b
-
- KLO=1
- KHI=N
- do while (KHI-KLO > 1)
- K=(KHI+KLO)/2
- IF(XA(K) > X)THEN
- KHI=K
- ELSE
- KLO=K
- ENDIF
- enddo
- H=XA(KHI)-XA(KLO)
- IF (H == 0.d0) stop 'Bad input in spline evaluation'
- A=(XA(KHI)-X)/H
- B=(X-XA(KLO))/H
-
- Y=A*YA(KLO)+B*YA(KHI)+((A**3-A)*Y2A(KLO)+ (B**3-B)*Y2A(KHI))*(H**2)/6.d0
-
- end subroutine SPLINT
-
Modified: seismo/2D/SPECFEM2D/trunk/meshfem2D.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/meshfem2D.f90 2007-02-15 00:06:00 UTC (rev 8501)
+++ seismo/2D/SPECFEM2D/trunk/meshfem2D.f90 2007-12-07 23:51:25 UTC (rev 8502)
@@ -17,6 +17,27 @@
!
!========================================================================
+! If you use this code for your own research, please cite:
+!
+! @ARTICLE{KoTr99,
+! author={D. Komatitsch and J. Tromp},
+! year=1999,
+! title={Introduction to the spectral-element method for 3-{D} seismic wave propagation},
+! journal={Geophys. J. Int.},
+! volume=139,
+! number=3,
+! pages={806-822},
+! doi={10.1046/j.1365-246x.1999.00967.x}}
+!
+! @ARTICLE{KoVi98,
+! author={D. Komatitsch and J. P. Vilotte},
+! title={The spectral-element method: an efficient tool to simulate the seismic response of 2{D} and 3{D} geological structures},
+! journal={Bull. Seismol. Soc. Am.},
+! year=1998,
+! volume=88,
+! number=2,
+! pages={368-392}}
+
program meshfem2D
implicit none
Copied: seismo/2D/SPECFEM2D/trunk/meshfem2D_circular_canyon.f90 (from rev 8501, seismo/2D/SPECFEM2D/trunk/circ.f90)
===================================================================
--- seismo/2D/SPECFEM2D/trunk/circ.f90 2007-02-15 00:06:00 UTC (rev 8501)
+++ seismo/2D/SPECFEM2D/trunk/meshfem2D_circular_canyon.f90 2007-12-07 23:51:25 UTC (rev 8502)
@@ -0,0 +1,826 @@
+
+!=====================================================================
+!
+! P r e m a i l l e u r - 2 D
+! ---------------------------
+!
+! Version 2.1
+! -----------
+!
+! Dimitri Komatitsch
+!
+! Departement de Sismologie
+! Institut de Physique du Globe de Paris
+!
+! (c) Institut de Physique du Globe de Paris, Octobre 1996
+!
+!=====================================================================
+
+! program to mesh a 2D canyon for a study with Paco Sanchez-Sesma and comparison with Kawase (1988)
+
+! DK DK Mexico August 1999 : mise a jour format base de donnees
+
+ program circular_canyon
+
+ implicit none
+
+! max size of the model in elements
+ integer, parameter :: mnx=7,mnz=7
+
+ double precision, parameter :: pi=3.141592653589793d0
+
+! seuil pour considerer deux points comme confondus
+ double precision, parameter :: rseuil=1.d-2
+
+! declare variables
+ integer imaxabs,n2ana,itimetype,isource_type,nump1,nump2,nump3,nump4
+ integer ndofn,ndime,ngnod,nnode,nbcnd,n1ana
+ integer nofst,npgeo,nspel,nbmodeles,nbsources,nrec,lquad,isamp,nrec1,nrec2
+ integer irec,imatnum,netyp,nxgll,nelemperio,nelemabs,nx,nz,i,j
+ integer irepr,nrecsur3,nt,niter,itaff,itfirstaff,numerocourant,pointsdisp,isubsamp
+
+ double precision R,theta_i,theta_init,delta_theta,eta_j,valseuil,freqmaxrep
+ double precision f0,t0,xs,zs,angle,factor,dist,xoffs,zoffs
+ double precision xrec,zrec,rho,cp,cs,anglerec
+ double precision anglerec2,dt,alphanewm,betanewm,gammanewm
+ double precision cutvect,cutcolor,scalex,scalez,sizemax,orig_x,orig_z
+ double precision factorana,factorxsu
+
+! stockage de la grille curvi (x et z)
+ integer, parameter :: npoinz1=(4*mnx+1)*(mnz+1), nelemz1=(4*mnx)*mnz
+ double precision x1(0:4*mnx,0:mnz)
+ double precision z1(0:4*mnx,0:mnz)
+
+ integer, parameter :: npoinz3=(2*mnx+1)*(4*mnz+1), nelemz3=(2*mnx)*(4*mnz)
+ double precision x3(0:2*mnx,0:4*mnz)
+ double precision z3(0:2*mnx,0:4*mnz)
+
+ integer, parameter :: npoinz4=(2*mnx+1)*(2*mnz+1), nelemz4=(2*mnx)*(2*mnz)
+ double precision x4(0:2*mnx,0:2*mnz)
+ double precision z4(0:2*mnx,0:2*mnz)
+
+ integer, parameter :: npoinz1b=(2*mnx+1)*(mnz+1), nelemz1b=(2*mnx)*mnz
+ double precision x1b(0:2*mnx,0:mnz)
+ double precision z1b(0:2*mnx,0:mnz)
+
+ integer, parameter :: npoinz2b=(mnx+1)*(2*mnz+1), nelemz2b=mnx*(2*mnz)
+ double precision x2b(0:mnx,0:2*mnz)
+ double precision z2b(0:mnx,0:2*mnz)
+
+ integer, parameter :: npoinz3b=(4*mnx+1)*(4*mnz+1), nelemz3b=(4*mnx)*(4*mnz)
+ double precision x3b(0:4*mnx,0:4*mnz)
+ double precision z3b(0:4*mnx,0:4*mnz)
+
+ integer, parameter :: npoinz4b=(2*mnx+1)*(2*mnz+1), nelemz4b=(2*mnx)*(2*mnz)
+ double precision x4b(0:2*mnx,0:2*mnz)
+ double precision z4b(0:2*mnx,0:2*mnz)
+
+! nombre max de points de maillage, et nombre exact d'elements
+ integer, parameter :: npoin = npoinz1+npoinz3+npoinz4+npoinz1b+npoinz2b+npoinz3b+npoinz4b
+ integer, parameter :: nelem = nelemz1+nelemz3+nelemz4+nelemz1b+nelemz2b+nelemz3b+nelemz4b
+
+! coordonnees geometriques des points
+ double precision xpoint(npoin)
+ double precision zpoint(npoin)
+
+! coordonnees des sommets de chaque element
+ double precision x1e(nelem)
+ double precision z1e(nelem)
+ double precision x2e(nelem)
+ double precision z2e(nelem)
+ double precision x3e(nelem)
+ double precision z3e(nelem)
+ double precision x4e(nelem)
+ double precision z4e(nelem)
+
+! numero des points des elements
+ integer numpoin1(nelem)
+ integer numpoin2(nelem)
+ integer numpoin3(nelem)
+ integer numpoin4(nelem)
+
+ character(len=50) title
+
+ logical iexternal, aleatoire, topoplane, simulate, absstacey
+ logical absorbhaut, absorbbas, absorbgauche, sismos
+ logical absorbdroite, absorbstacey, absorbmodar, ifullabs
+
+ logical display, ignuplot, ivectplot, icolorplot, imeshvect
+ logical imeshcolor, imodelvect, iboundvect, interpol, isymbols, initialfield
+ logical usletter,compenergy
+
+ print *,'Nombre d''elements = ',nelem
+ print *,'Nombre max de points = ',npoin
+
+ nx = mnx
+ nz = mnz
+
+ R = 1.
+
+! ***************************************
+! *** ZONE DE DROITE
+! ***************************************
+
+! generer les points de base de l'interpolation lineaire (zone 1)
+ theta_init = 3 * pi / 2.
+ delta_theta = pi / 2.
+ do i=0,4*nx
+
+! --- point de depart
+ if(i < 2*nx) then
+ x1(i,0) = 2.*R * dble(i) / dble(2*nx)
+ z1(i,0) = - 2.*R
+ else
+ x1(i,0) = 2.*R
+ z1(i,0) = - 2.*R * (1. - dble(i - 2*nx) / dble(2*nx))
+ endif
+
+! --- point d'arrivee
+ theta_i = theta_init + delta_theta * dble(i) / dble(4*nx)
+ x1(i,nz) = cos(theta_i)
+ z1(i,nz) = sin(theta_i)
+
+! --- points intermediaires par interpolation lineaire
+ do j=1,nz-1
+ eta_j = dble(j) / dble(nz)
+ x1(i,j) = (1.-eta_j)*x1(i,0) + eta_j*x1(i,nz)
+ z1(i,j) = (1.-eta_j)*z1(i,0) + eta_j*z1(i,nz)
+ enddo
+ enddo
+
+! generer zone de gauche (zone 3)
+ do i=0,2*nx
+ do j=0,4*nz
+ x3(i,j) = 5. * dble(i) / dble(2*nx) + 2.
+ if(j <= 2*nz) then
+ z3(i,j) = 7. * dble(j) / dble(2*nz) - 9.
+ else
+ z3(i,j) = 2. * dble(j-2*nz) / dble(2*nz) - 2.
+ endif
+ enddo
+ enddo
+
+! generer zone du bas (zone 4)
+ do i=0,2*nx
+ do j=0,2*nz
+ x4(i,j) = 2. * dble(i) / dble(2*nx)
+ z4(i,j) = 7. * dble(j) / dble(2*nz) - 9.
+ enddo
+ enddo
+
+! ***************************************
+! *** ZONE DE GAUCHE
+! ***************************************
+
+! generer les points de base de l'interpolation lineaire (zone 1)
+ theta_init = pi / 4.
+ delta_theta = pi / 4.
+
+ do i=0,2*nx
+! --- point de depart
+ x1b(i,0) = 2.*R * (dble(i) / dble(2*nx) - 1.)
+ z1b(i,0) = - 2.*R
+
+! --- point d'arrivee
+ theta_i = theta_init + delta_theta * dble(i) / dble(2*nx)
+ x1b(i,nz) = - cos(theta_i)
+ z1b(i,nz) = - sin(theta_i)
+
+! --- points intermediaires par interpolation lineaire
+ do j=1,nz-1
+ eta_j = dble(j) / dble(nz)
+ x1b(i,j) = (1.-eta_j)*x1b(i,0) + eta_j*x1b(i,nz)
+ z1b(i,j) = (1.-eta_j)*z1b(i,0) + eta_j*z1b(i,nz)
+ enddo
+ enddo
+
+! generer les points de base de l'interpolation lineaire (zone 2)
+ theta_init = pi / 4.
+
+ do j=0,2*nz
+! --- point de depart
+ x2b(0,j) = - 2.*R
+ z2b(0,j) = 2.*R * (dble(j) / dble(2*nz) - 1.)
+
+! --- point d'arrivee
+ theta_i = theta_init - delta_theta * dble(j) / dble(2*nz)
+ x2b(nx,j) = - cos(theta_i)
+ z2b(nx,j) = - sin(theta_i)
+
+! --- points intermediaires par interpolation lineaire
+ do i=1,nx-1
+ eta_j = dble(i) / dble(nx)
+ x2b(i,j) = (1.-eta_j)*x2b(0,j) + eta_j*x2b(nx,j)
+ z2b(i,j) = (1.-eta_j)*z2b(0,j) + eta_j*z2b(nx,j)
+ enddo
+
+ enddo
+
+! generer zone de gauche (zone 3)
+ do i=0,4*nx
+ do j=0,4*nz
+ x3b(i,j) = 10. * dble(i) / dble(4*nx) - 12.
+ if(j <= 2*nz) then
+ z3b(i,j) = 7. * dble(j) / dble(2*nz) - 9.
+ else
+ z3b(i,j) = 2. * dble(j-2*nz) / dble(2*nz) - 2.
+ endif
+ enddo
+ enddo
+
+! generer zone du bas (zone 4)
+ do i=0,2*nx
+ do j=0,2*nz
+ x4b(i,j) = 2. * dble(i) / dble(2*nx) - 2.
+ z4b(i,j) = 7. * dble(j) / dble(2*nz) - 9.
+ enddo
+ enddo
+
+! ***
+! *** generer un fichier 'GNUPLOT' pour le controle de la grille ***
+! ***
+
+ write(*,*)
+ write(*,*) 'Ecriture de la grille format GNUPLOT...'
+
+ open(unit=20,file='grid.gnu',status='unknown')
+
+! *** dessiner la zone 1
+ do j=0,nz
+ do i=0,4*nx-1
+ write(20,*) sngl(x1(i,j)),sngl(z1(i,j))
+ write(20,*) sngl(x1(i+1,j)),sngl(z1(i+1,j))
+ write(20,100)
+ enddo
+ enddo
+
+ do i=0,4*nx
+ do j=0,nz-1
+ write(20,*) sngl(x1(i,j)),sngl(z1(i,j))
+ write(20,*) sngl(x1(i,j+1)),sngl(z1(i,j+1))
+ write(20,100)
+ enddo
+ enddo
+
+! *** dessiner la zone 3
+ do j=0,4*nz
+ do i=0,2*nx-1
+ write(20,*) sngl(x3(i,j)),sngl(z3(i,j))
+ write(20,*) sngl(x3(i+1,j)),sngl(z3(i+1,j))
+ write(20,100)
+ enddo
+ enddo
+
+ do i=0,2*nx
+ do j=0,4*nz-1
+ write(20,*) sngl(x3(i,j)),sngl(z3(i,j))
+ write(20,*) sngl(x3(i,j+1)),sngl(z3(i,j+1))
+ write(20,100)
+ enddo
+ enddo
+
+! *** dessiner la zone 4
+ do j=0,2*nz
+ do i=0,2*nx-1
+ write(20,*) sngl(x4(i,j)),sngl(z4(i,j))
+ write(20,*) sngl(x4(i+1,j)),sngl(z4(i+1,j))
+ write(20,100)
+ enddo
+ enddo
+
+ do i=0,2*nx
+ do j=0,2*nz-1
+ write(20,*) sngl(x4(i,j)),sngl(z4(i,j))
+ write(20,*) sngl(x4(i,j+1)),sngl(z4(i,j+1))
+ write(20,100)
+ enddo
+ enddo
+
+! *** dessiner la zone 1
+ do j=0,nz
+ do i=0,2*nx-1
+ write(20,*) sngl(x1b(i,j)),sngl(z1b(i,j))
+ write(20,*) sngl(x1b(i+1,j)),sngl(z1b(i+1,j))
+ write(20,100)
+ enddo
+ enddo
+
+ do i=0,2*nx
+ do j=0,nz-1
+ write(20,*) sngl(x1b(i,j)),sngl(z1b(i,j))
+ write(20,*) sngl(x1b(i,j+1)),sngl(z1b(i,j+1))
+ write(20,100)
+ enddo
+ enddo
+
+! *** dessiner la zone 2
+ do j=0,2*nz
+ do i=0,nx-1
+ write(20,*) sngl(x2b(i,j)),sngl(z2b(i,j))
+ write(20,*) sngl(x2b(i+1,j)),sngl(z2b(i+1,j))
+ write(20,100)
+ enddo
+ enddo
+
+ do i=0,nx
+ do j=0,2*nz-1
+ write(20,*) sngl(x2b(i,j)),sngl(z2b(i,j))
+ write(20,*) sngl(x2b(i,j+1)),sngl(z2b(i,j+1))
+ write(20,100)
+ enddo
+ enddo
+
+! *** dessiner la zone 3
+ do j=0,4*nz
+ do i=0,4*nx-1
+ write(20,*) sngl(x3b(i,j)),sngl(z3b(i,j))
+ write(20,*) sngl(x3b(i+1,j)),sngl(z3b(i+1,j))
+ write(20,100)
+ enddo
+ enddo
+
+ do i=0,4*nx
+ do j=0,4*nz-1
+ write(20,*) sngl(x3b(i,j)),sngl(z3b(i,j))
+ write(20,*) sngl(x3b(i,j+1)),sngl(z3b(i,j+1))
+ write(20,100)
+ enddo
+ enddo
+
+! *** dessiner la zone 4
+ do j=0,2*nz
+ do i=0,2*nx-1
+ write(20,*) sngl(x4b(i,j)),sngl(z4b(i,j))
+ write(20,*) sngl(x4b(i+1,j)),sngl(z4b(i+1,j))
+ write(20,100)
+ enddo
+ enddo
+
+ do i=0,2*nx
+ do j=0,2*nz-1
+ write(20,*) sngl(x4b(i,j)),sngl(z4b(i,j))
+ write(20,*) sngl(x4b(i,j+1)),sngl(z4b(i,j+1))
+ write(20,100)
+ enddo
+ enddo
+
+ close(20)
+
+ write(*,*) 'Fin ecriture de la grille format GNUPLOT'
+ write(*,*)
+
+ 100 format('')
+
+! ***
+! *** generer la liste des points geometriques
+! ***
+
+ numerocourant = 1
+
+! *** zone 1
+ do j=0,nz
+ do i=0,4*nx
+ xpoint(numerocourant) = x1(i,j)
+ zpoint(numerocourant) = z1(i,j)
+ numerocourant = numerocourant + 1
+ enddo
+ enddo
+
+! *** zone 3
+ do j=0,4*nz
+ do i=0,2*nx
+ xpoint(numerocourant) = x3(i,j)
+ zpoint(numerocourant) = z3(i,j)
+ numerocourant = numerocourant + 1
+ enddo
+ enddo
+
+! *** zone 4
+ do j=0,2*nz
+ do i=0,2*nx
+ xpoint(numerocourant) = x4(i,j)
+ zpoint(numerocourant) = z4(i,j)
+ numerocourant = numerocourant + 1
+ enddo
+ enddo
+
+! *** zone 1
+ do j=0,nz
+ do i=0,2*nx
+ xpoint(numerocourant) = x1b(i,j)
+ zpoint(numerocourant) = z1b(i,j)
+ numerocourant = numerocourant + 1
+ enddo
+ enddo
+
+! *** zone 2
+ do j=0,2*nz
+ do i=0,nx
+ xpoint(numerocourant) = x2b(i,j)
+ zpoint(numerocourant) = z2b(i,j)
+ numerocourant = numerocourant + 1
+ enddo
+ enddo
+
+! *** zone 3
+ do j=0,4*nz
+ do i=0,4*nx
+ xpoint(numerocourant) = x3b(i,j)
+ zpoint(numerocourant) = z3b(i,j)
+ numerocourant = numerocourant + 1
+ enddo
+ enddo
+
+! *** zone 4
+ do j=0,2*nz
+ do i=0,2*nx
+ xpoint(numerocourant) = x4b(i,j)
+ zpoint(numerocourant) = z4b(i,j)
+ numerocourant = numerocourant + 1
+ enddo
+ enddo
+
+ print *,'nb de points stockes = ',numerocourant - 1
+
+! ***
+! *** generer la liste des elements
+! ***
+
+ numerocourant = 1
+ imaxabs = 0
+
+! *** zone 1
+ do j=0,nz-1
+ do i=0,4*nx-1
+ x1e(numerocourant) = x1(i,j)
+ z1e(numerocourant) = z1(i,j)
+ x2e(numerocourant) = x1(i+1,j)
+ z2e(numerocourant) = z1(i+1,j)
+ x3e(numerocourant) = x1(i+1,j+1)
+ z3e(numerocourant) = z1(i+1,j+1)
+ x4e(numerocourant) = x1(i,j+1)
+ z4e(numerocourant) = z1(i,j+1)
+ numerocourant = numerocourant + 1
+ enddo
+ enddo
+
+! *** zone 3
+ do j=0,4*nz-1
+ do i=0,2*nx-1
+ x1e(numerocourant) = x3(i,j)
+ z1e(numerocourant) = z3(i,j)
+ x2e(numerocourant) = x3(i+1,j)
+ z2e(numerocourant) = z3(i+1,j)
+ x3e(numerocourant) = x3(i+1,j+1)
+ z3e(numerocourant) = z3(i+1,j+1)
+ x4e(numerocourant) = x3(i,j+1)
+ z4e(numerocourant) = z3(i,j+1)
+ numerocourant = numerocourant + 1
+ enddo
+ enddo
+
+! *** zone 4
+ do j=0,2*nz-1
+ do i=0,2*nx-1
+ x1e(numerocourant) = x4(i,j)
+ z1e(numerocourant) = z4(i,j)
+ x2e(numerocourant) = x4(i+1,j)
+ z2e(numerocourant) = z4(i+1,j)
+ x3e(numerocourant) = x4(i+1,j+1)
+ z3e(numerocourant) = z4(i+1,j+1)
+ x4e(numerocourant) = x4(i,j+1)
+ z4e(numerocourant) = z4(i,j+1)
+ numerocourant = numerocourant + 1
+ enddo
+ enddo
+
+! *** zone 1
+ do j=0,nz-1
+ do i=0,2*nx-1
+ x1e(numerocourant) = x1b(i,j)
+ z1e(numerocourant) = z1b(i,j)
+ x2e(numerocourant) = x1b(i+1,j)
+ z2e(numerocourant) = z1b(i+1,j)
+ x3e(numerocourant) = x1b(i+1,j+1)
+ z3e(numerocourant) = z1b(i+1,j+1)
+ x4e(numerocourant) = x1b(i,j+1)
+ z4e(numerocourant) = z1b(i,j+1)
+ numerocourant = numerocourant + 1
+ enddo
+ enddo
+
+! *** zone 2
+ do j=0,2*nz-1
+ do i=0,nx-1
+ x1e(numerocourant) = x2b(i,j)
+ z1e(numerocourant) = z2b(i,j)
+ x2e(numerocourant) = x2b(i+1,j)
+ z2e(numerocourant) = z2b(i+1,j)
+ x3e(numerocourant) = x2b(i+1,j+1)
+ z3e(numerocourant) = z2b(i+1,j+1)
+ x4e(numerocourant) = x2b(i,j+1)
+ z4e(numerocourant) = z2b(i,j+1)
+ numerocourant = numerocourant + 1
+ enddo
+ enddo
+
+! *** zone 3
+ do j=0,4*nz-1
+ do i=0,4*nx-1
+ x1e(numerocourant) = x3b(i,j)
+ z1e(numerocourant) = z3b(i,j)
+ x2e(numerocourant) = x3b(i+1,j)
+ z2e(numerocourant) = z3b(i+1,j)
+ x3e(numerocourant) = x3b(i+1,j+1)
+ z3e(numerocourant) = z3b(i+1,j+1)
+ x4e(numerocourant) = x3b(i,j+1)
+ z4e(numerocourant) = z3b(i,j+1)
+ numerocourant = numerocourant + 1
+ enddo
+ enddo
+
+! *** zone 4
+ do j=0,2*nz-1
+ do i=0,2*nx-1
+ x1e(numerocourant) = x4b(i,j)
+ z1e(numerocourant) = z4b(i,j)
+ x2e(numerocourant) = x4b(i+1,j)
+ z2e(numerocourant) = z4b(i+1,j)
+ x3e(numerocourant) = x4b(i+1,j+1)
+ z3e(numerocourant) = z4b(i+1,j+1)
+ x4e(numerocourant) = x4b(i,j+1)
+ z4e(numerocourant) = z4b(i,j+1)
+ numerocourant = numerocourant + 1
+ enddo
+ enddo
+
+ print *,'nb d''elements stockes = ',numerocourant - 1
+
+! ***
+! *** creation des elements sous forme topologique
+! ***
+
+ write(*,*)
+ write(*,*) 'Creation de la topologie des elements...'
+
+ do i=1,nelem
+
+! recherche point 1
+ do j=1,npoin
+ dist = sqrt((x1e(i)-xpoint(j))**2 + (z1e(i)-zpoint(j))**2)
+ if(dist <= rseuil) then
+ nump1 = j
+ goto 401
+ endif
+ enddo
+ stop 'point not found !'
+ 401 continue
+
+! recherche point 2
+ do j=1,npoin
+ dist = sqrt((x2e(i)-xpoint(j))**2 + (z2e(i)-zpoint(j))**2)
+ if(dist <= rseuil) then
+ nump2 = j
+ goto 402
+ endif
+ enddo
+ stop 'point not found !'
+ 402 continue
+
+! recherche point 3
+ do j=1,npoin
+ dist = sqrt((x3e(i)-xpoint(j))**2 + (z3e(i)-zpoint(j))**2)
+ if(dist <= rseuil) then
+ nump3 = j
+ goto 403
+ endif
+ enddo
+ stop 'point not found !'
+ 403 continue
+
+! recherche point 4
+ do j=1,npoin
+ dist = sqrt((x4e(i)-xpoint(j))**2 + (z4e(i)-zpoint(j))**2)
+ if(dist <= rseuil) then
+ nump4 = j
+ goto 404
+ endif
+ enddo
+ stop 'point not found !'
+ 404 continue
+
+ numpoin1(i) = nump1
+ numpoin2(i) = nump2
+ numpoin3(i) = nump3
+ numpoin4(i) = nump4
+
+ enddo
+
+! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+! *** generation de la base de donnees
+
+ write(*,*)
+ write(*,*) 'Generation de la base de donnees...'
+
+ open(unit=15,file='OUTPUT_FILES/Database',status='unknown')
+
+ title = 'Modele Canyon Paco'
+ write(15,*) '#'
+ write(15,*) '# Base de Donnees pour Specfem - Premailleur Fortran 90'
+ write(15,*) '# ',title
+ write(15,*) '# Dimitri Komatitsch, (c) EPS - Harvard February 1998'
+ write(15,*) '#'
+
+ write(15,*) 'Titre simulation'
+ write(15,40) title
+
+ ndofn = 2
+ ndime = 2
+ ngnod = 4
+ nnode = 4
+ nbcnd = 0
+ nofst = 0
+ npgeo = npoin
+ nspel = nelem
+ nbmodeles = 1
+ nbsources = 1
+ nrec = 150
+ lquad = 1
+ iexternal = .false.
+ aleatoire = .false.
+ topoplane = .false.
+ simulate = .false.
+
+ absorbhaut = .false.
+ absorbbas = .false.
+ absorbgauche = .false.
+ absorbdroite = .false.
+ absorbstacey = .true.
+ absorbmodar = .false.
+ ifullabs = .false.
+
+ sismos = .true.
+ isamp = 20
+ nrec1 = nrec
+ nrec2 = 0
+ anglerec = 0.
+ anglerec2 = 0.
+ irepr = 1
+ nrecsur3 = nrec / 3
+
+ nt = 20000
+ dt = 0.625e-3
+ niter = 1
+ alphanewm = 0.
+ betanewm = 0.
+ gammanewm = 0.5
+ display = .true.
+ ignuplot = .false.
+ ivectplot = .true.
+ icolorplot = .false.
+ imeshvect = .true.
+ imeshcolor = .false.
+ imodelvect = .false.
+ iboundvect = .false.
+ interpol = .true.
+ isymbols = .true.
+
+!! DK DK Mexico August 1999, temporarily suppress external field
+ initialfield = .true.
+ initialfield = .false.
+
+ itaff = 2000
+ itfirstaff = 5
+ cutvect = 1.
+ cutcolor = 2.2
+ scalex = 1.
+ scalez = 1.
+ sizemax = 1.
+ pointsdisp = 7
+ isubsamp = 2
+ orig_x = 2.3
+ orig_z = 3.4
+ factorana = 50000.
+ factorxsu = 3.5
+ n1ana = 1
+ n2ana = nrec1
+
+ write(15,*) 'ndofn ndime npgeo'
+ write(15,*) ndofn,ndime,npgeo
+
+ write(15,*) 'display ignuplot interpol'
+ write(15,*) display,ignuplot,interpol
+
+ write(15,*) 'itaff itfirstaff icolor inumber'
+ write(15,*) itaff,itfirstaff,0,0
+
+ write(15,*) 'ivectplot imeshvect imodelvect iboundvect cutvect isubsamp'
+ write(15,*) ivectplot,imeshvect,imodelvect,iboundvect,cutvect,isubsamp
+
+ usletter = .true.
+ write(15,*) 'scalex scalez sizemax angle rapport USletter'
+ write(15,*) scalex,scalez,sizemax,20.,0.40,usletter
+
+ write(15,*) 'orig_x orig_z isymbols'
+ write(15,*) orig_x,orig_z,isymbols
+
+ valseuil = 5.00
+ freqmaxrep = 100.
+ write(15,*) 'valseuil freqmaxrep'
+ write(15,*) valseuil,freqmaxrep
+
+ write(15,*) 'sismos nrec nrec1 nrec2 isamp'
+ write(15,*) sismos,nrec,nrec1,nrec2,isamp
+
+ write(15,*) 'irepr anglerec anglerec2'
+ write(15,*) irepr,anglerec,anglerec2
+
+ compenergy = .false.
+ absstacey = .true.
+ write(15,*) 'topoplane absstacey compenergy'
+ write(15,*) topoplane,absstacey,compenergy
+
+ write(15,*) 'initialfield factorana factorxsu n1ana n2ana'
+ write(15,*) initialfield,factorana,factorxsu,n1ana,n2ana
+
+ write(15,*) 'isismostype ivecttype iaffinfo'
+ write(15,*) '1, 1, 40'
+ write(15,*) 'ireadmodel ioutputgrid iavs ivisual3'
+ write(15,*) 'F, F, F, F'
+
+ write(15,*) 'iexec iecho'
+ write(15,*) '1 1'
+
+ write(15,*) 'ncycl dtinc niter'
+ write(15,*) nt,dt,niter
+
+ write(15,*) 'alpha beta gamma'
+ write(15,*) alphanewm,betanewm,gammanewm
+
+ nbsources = 1
+ write(15,*) 'nltfl (number of force or pressure sources)'
+ write(15,*) nbsources
+
+ itimetype = 6
+ isource_type = 2
+ f0 = 2.
+ t0 = 0.55
+ xs = +1.
+ zs = -2.
+ angle = 0.
+ factor = 1.
+ xoffs = 12.
+ zoffs = 9.
+ write(15,*) 'Collocated forces and/or pressure sources:'
+ write(15,*) itimetype,isource_type,xs+xoffs,zs+zoffs,f0,t0,factor,angle,0
+
+ write(15,*) 'Receivers (number, angle, position in meters)'
+ do irec=1,nrec
+ if(irec <= nrecsur3) then
+ xrec = 2.*dble(irec-1)/dble(nrecsur3-1) + 9.
+ zrec = 9.
+ else if(irec >= 2*nrecsur3) then
+ xrec = 2.*dble(irec-2*nrecsur3)/dble(nrecsur3) + 13.
+ zrec = 9.
+ else
+ angle = pi + pi*dble(irec-nrecsur3)/dble(nrecsur3)
+ xrec = 12. + cos(angle)
+ zrec = 9. + sin(angle)
+ endif
+ write(15,*) irec,xrec,zrec
+ enddo
+
+ write(15,*) 'Coordinates of spectral control points'
+ do i=1,npoin
+ write(15,*) i,xpoint(i)+xoffs,zpoint(i)+zoffs
+ enddo
+
+ netyp = 2
+ nxgll = 6
+ nelemperio = 0
+ nelemabs = 0
+
+ write(15,*) 'params spectraux'
+ write(15,*) netyp,nbmodeles,ngnod,nxgll,nxgll,nspel,pointsdisp,nelemabs,nelemperio
+
+ write(15,*) 'Material sets (num 0 rho vp vs 0 0)'
+ rho = 1.
+ cp = 2.
+ cs = 1.
+ write(15,*) nbmodeles,0,rho,cp,cs,0,0
+
+ write(15,*) 'Spectral-element topology'
+
+ imatnum = 1
+
+ do i=1,nspel
+ write(15,*) i,imatnum,numpoin1(i),numpoin2(i),numpoin3(i),numpoin4(i)
+ enddo
+
+ close(15)
+
+ 40 format(a50)
+
+ end program circular_canyon
+
Copied: seismo/2D/SPECFEM2D/trunk/meshfem2D_non_struct_2.f90 (from rev 8501, seismo/2D/SPECFEM2D/trunk/maille_non_struct_2.f90)
Copied: seismo/2D/SPECFEM2D/trunk/meshfem2D_non_struct_3.f90 (from rev 8501, seismo/2D/SPECFEM2D/trunk/maille_non_struct_3.f90)
Modified: seismo/2D/SPECFEM2D/trunk/specfem2D.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/specfem2D.f90 2007-02-15 00:06:00 UTC (rev 8501)
+++ seismo/2D/SPECFEM2D/trunk/specfem2D.f90 2007-12-07 23:51:25 UTC (rev 8502)
@@ -52,6 +52,27 @@
! Therefore, in order to have a Ricker source in pressure, use the first-derivative
! of a Gaussian for the time source for Chi in the parameter file.
+! If you use this code for your own research, please cite:
+!
+! @ARTICLE{KoTr99,
+! author={D. Komatitsch and J. Tromp},
+! year=1999,
+! title={Introduction to the spectral-element method for 3-{D} seismic wave propagation},
+! journal={Geophys. J. Int.},
+! volume=139,
+! number=3,
+! pages={806-822},
+! doi={10.1046/j.1365-246x.1999.00967.x}}
+!
+! @ARTICLE{KoVi98,
+! author={D. Komatitsch and J. P. Vilotte},
+! title={The spectral-element method: an efficient tool to simulate the seismic response of 2{D} and 3{D} geological structures},
+! journal={Bull. Seismol. Soc. Am.},
+! year=1998,
+! volume=88,
+! number=2,
+! pages={368-392}}
+
program specfem2D
implicit none
More information about the cig-commits
mailing list