[cig-commits] r8506 - in seismo/2D/SPECFEM2D/trunk: . UTILS
walter at geodynamics.org
walter at geodynamics.org
Fri Dec 7 15:51:50 PST 2007
Author: walter
Date: 2007-12-07 15:51:50 -0800 (Fri, 07 Dec 2007)
New Revision: 8506
Added:
seismo/2D/SPECFEM2D/trunk/UTILS/meshfem2D_circular_canyon.f90
seismo/2D/SPECFEM2D/trunk/UTILS/meshfem2D_non_struct_2.f90
seismo/2D/SPECFEM2D/trunk/UTILS/meshfem2D_non_struct_3.f90
seismo/2D/SPECFEM2D/trunk/UTILS/su_create_Paul.python
Removed:
seismo/2D/SPECFEM2D/trunk/create_earth_model.f90
seismo/2D/SPECFEM2D/trunk/defarrays.f90
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
seismo/2D/SPECFEM2D/trunk/su_create_Paul.python
Log:
deleted two files and moved four old files (meshers for particular geometries or models) to directory UTILS
Copied: seismo/2D/SPECFEM2D/trunk/UTILS/meshfem2D_circular_canyon.f90 (from rev 8505, seismo/2D/SPECFEM2D/trunk/meshfem2D_circular_canyon.f90)
Copied: seismo/2D/SPECFEM2D/trunk/UTILS/meshfem2D_non_struct_2.f90 (from rev 8505, seismo/2D/SPECFEM2D/trunk/meshfem2D_non_struct_2.f90)
Copied: seismo/2D/SPECFEM2D/trunk/UTILS/meshfem2D_non_struct_3.f90 (from rev 8505, seismo/2D/SPECFEM2D/trunk/meshfem2D_non_struct_3.f90)
Copied: seismo/2D/SPECFEM2D/trunk/UTILS/su_create_Paul.python (from rev 8505, seismo/2D/SPECFEM2D/trunk/su_create_Paul.python)
Deleted: seismo/2D/SPECFEM2D/trunk/create_earth_model.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/create_earth_model.f90 2007-03-01 16:00:05 UTC (rev 8505)
+++ seismo/2D/SPECFEM2D/trunk/create_earth_model.f90 2007-12-07 23:51:50 UTC (rev 8506)
@@ -1,72 +0,0 @@
-
-!========================================================================
-!
-! S P E C F E M 2 D Version 5.1
-! ------------------------------
-!
-! Dimitri Komatitsch
-! Universite de Pau et des Pays de l'Adour, France
-!
-! (c) January 2005
-!
-!========================================================================
-
-! modify an external grid file (list of points and coordinates) to include the
-! velocity model (rho, vp, vs, in this order)
-
- program create_earth_model
-
- implicit none
-
- integer ipoin,npoin
-
- double precision rho,vp,vs
-
- double precision, dimension(:), allocatable :: xgrid,zgrid
-
- include "constants.h"
-
-! read the grid from an existing text file
- print *
- print *,'Reading the grid from an existing text file...'
- print *
-
- open(unit=55,file='OUTPUT_FILES/grid_points_and_model.txt',status='old')
-
- read(55,*) npoin
-
- print *,'There are ',npoin,' grid points'
-
- allocate(xgrid(npoin))
- allocate(zgrid(npoin))
-
- do ipoin = 1,npoin
- read(55,*) xgrid(ipoin),zgrid(ipoin)
- enddo
-
- close(55)
-
-! write the velocity model to the same text file
- print *
- print *,'Saving the grid and the velocity model in the same text file...'
- print *
-
- open(unit=55,file='OUTPUT_FILES/grid_points_and_model.txt',status='unknown')
-
- write(55,*) npoin
-
- do ipoin = 1,npoin
-
-! user should change this to assign these values depending on the position of the grid point
- rho = 2200.d0
- vp = 3000.d0
- vs = 1732.d0
-
- write(55,*) xgrid(ipoin),zgrid(ipoin),rho,vp,vs
-
- enddo
-
- close(55)
-
- end program create_earth_model
-
Deleted: seismo/2D/SPECFEM2D/trunk/defarrays.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/defarrays.f90 2007-03-01 16:00:05 UTC (rev 8505)
+++ seismo/2D/SPECFEM2D/trunk/defarrays.f90 2007-12-07 23:51:50 UTC (rev 8506)
@@ -1,166 +0,0 @@
-
-!========================================================================
-!
-! S P E C F E M 2 D Version 5.1
-! ------------------------------
-!
-! Dimitri Komatitsch
-! Universite de Pau et des Pays de l'Adour, France
-!
-! (c) January 2005
-!
-!========================================================================
-
- subroutine defarrays(vpext,vsext,rhoext,density,elastcoef, &
- ibool,kmato,coord,npoin,rsizemin,rsizemax, &
- cpoverdxmax,lambdaSmin,lambdaSmax,lambdaPmin,lambdaPmax, &
- vpmin,vpmax,read_external_model,nspec,numat)
-
-! define all the arrays for the variational formulation
-
- implicit none
-
- include "constants.h"
-
- integer i,j,ispec,material,ipointnum,npoin,nspec,numat
-
- integer kmato(nspec),ibool(NGLLX,NGLLX,nspec)
-
- double precision density(numat),elastcoef(4,numat)
-
- double precision coord(NDIM,npoin)
-
- double precision vpext(npoin)
- double precision vsext(npoin)
- double precision rhoext(npoin)
-
- double precision vsmin,vsmax,densmin,densmax
- double precision lambdaplus2mu,lambda,mu,denst
- double precision kappa,cploc,csloc,x0,z0
- double precision x1,z1,x2,z2,rdist1,rdist2,rapportmin,rapportmax
- double precision lambdamin,lambdamax
-
- double precision rsizemin,rsizemax,cpoverdxmax, &
- lambdaSmin,lambdaSmax,lambdaPmin,lambdaPmax,vpmin,vpmax
-
- logical read_external_model
-
-!
-!-----------------------------------------------------------------------
-!
-
-!---- compute parameters for the spectral elements
-
- vpmin = HUGEVAL
- vsmin = HUGEVAL
- vpmax = -HUGEVAL
- vsmax = -HUGEVAL
- densmin = HUGEVAL
- densmax = -HUGEVAL
-
- rsizemin = HUGEVAL
- rsizemax = -HUGEVAL
-
- cpoverdxmax = -HUGEVAL
-
- lambdaPmin = HUGEVAL
- lambdaSmin = HUGEVAL
- lambdaPmax = -HUGEVAL
- lambdaSmax = -HUGEVAL
-
- do ispec=1,nspec
-
- material = kmato(ispec)
-
- lambda = elastcoef(1,material)
- mu = elastcoef(2,material)
- lambdaplus2mu = elastcoef(3,material)
- denst = density(material)
-
- kappa = lambda + 2.d0*mu/3.d0
-
- cploc = sqrt((kappa + 4.d0*mu/3.d0)/denst)
- csloc = sqrt(mu/denst)
-
- do j=1,NGLLZ
- do i=1,NGLLX
-
-!--- si formulation heterogene pour un modele de vitesse externe
- if(read_external_model) then
- ipointnum = ibool(i,j,ispec)
- cploc = vpext(ipointnum)
- csloc = vsext(ipointnum)
- denst = rhoext(ipointnum)
- mu = denst*csloc*csloc
- lambda = denst*cploc*cploc - 2.d0*mu
- lambdaplus2mu = lambda + 2.d0*mu
- endif
-
-!--- calculer min et max du modele de vitesse et densite
- vpmin = min(vpmin,cploc)
- vpmax = max(vpmax,cploc)
-
- vsmin = min(vsmin,csloc)
- vsmax = max(vsmax,csloc)
-
- densmin = min(densmin,denst)
- densmax = max(densmax,denst)
-
-!--- stocker parametres pour verifs diverses
- if(i < NGLLX .and. j < NGLLZ) then
-
- x0 = coord(1,ibool(i,j,ispec))
- z0 = coord(2,ibool(i,j,ispec))
- x1 = coord(1,ibool(i+1,j,ispec))
- z1 = coord(2,ibool(i+1,j,ispec))
- x2 = coord(1,ibool(i,j+1,ispec))
- z2 = coord(2,ibool(i,j+1,ispec))
-
- rdist1 = sqrt((x1-x0)**2 + (z1-z0)**2)
- rdist2 = sqrt((x2-x0)**2 + (z2-z0)**2)
-
- rsizemin = min(rsizemin,rdist1)
- rsizemin = min(rsizemin,rdist2)
- rsizemax = max(rsizemax,rdist1)
- rsizemax = max(rsizemax,rdist2)
-
- rapportmin = cploc / max(rdist1,rdist2)
- rapportmax = cploc / min(rdist1,rdist2)
- cpoverdxmax = max(cpoverdxmax,rapportmax)
-
- x0 = coord(1,ibool(1,1,ispec))
- z0 = coord(2,ibool(1,1,ispec))
- x1 = coord(1,ibool(NGLLX,1,ispec))
- z1 = coord(2,ibool(NGLLX,1,ispec))
- x2 = coord(1,ibool(1,NGLLZ,ispec))
- z2 = coord(2,ibool(1,NGLLZ,ispec))
-
- rdist1 = sqrt((x1-x0)**2 + (z1-z0)**2)
- rdist2 = sqrt((x2-x0)**2 + (z2-z0)**2)
-
- lambdamin = cploc/max(rdist1,rdist2)
- lambdamax = cploc/min(rdist1,rdist2)
- lambdaPmin = min(lambdaPmin,lambdamin)
- lambdaPmax = max(lambdaPmax,lambdamax)
-
- lambdamin = csloc/max(rdist1,rdist2)
- lambdamax = csloc/min(rdist1,rdist2)
- lambdaSmin = min(lambdaSmin,lambdamin)
- lambdaSmax = max(lambdaSmax,lambdamax)
-
- endif
-
- enddo
- enddo
- enddo
-
- write(IOUT,*)
- write(IOUT,*) '********'
- write(IOUT,*) 'Modele : vitesse P min,max = ',vpmin,vpmax
- write(IOUT,*) 'Modele : vitesse S min,max = ',vsmin,vsmax
- write(IOUT,*) 'Modele : densite min,max = ',densmin,densmax
- write(IOUT,*) '********'
- write(IOUT,*)
-
- end subroutine defarrays
-
Deleted: seismo/2D/SPECFEM2D/trunk/meshfem2D_circular_canyon.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/meshfem2D_circular_canyon.f90 2007-03-01 16:00:05 UTC (rev 8505)
+++ seismo/2D/SPECFEM2D/trunk/meshfem2D_circular_canyon.f90 2007-12-07 23:51:50 UTC (rev 8506)
@@ -1,826 +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
-!
-!=====================================================================
-
-! 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
-
Deleted: seismo/2D/SPECFEM2D/trunk/meshfem2D_non_struct_2.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/meshfem2D_non_struct_2.f90 2007-03-01 16:00:05 UTC (rev 8505)
+++ seismo/2D/SPECFEM2D/trunk/meshfem2D_non_struct_2.f90 2007-12-07 23:51:50 UTC (rev 8506)
@@ -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/meshfem2D_non_struct_3.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/meshfem2D_non_struct_3.f90 2007-03-01 16:00:05 UTC (rev 8505)
+++ seismo/2D/SPECFEM2D/trunk/meshfem2D_non_struct_3.f90 2007-12-07 23:51:50 UTC (rev 8506)
@@ -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
-
Deleted: seismo/2D/SPECFEM2D/trunk/su_create_Paul.python
===================================================================
--- seismo/2D/SPECFEM2D/trunk/su_create_Paul.python 2007-03-01 16:00:05 UTC (rev 8505)
+++ seismo/2D/SPECFEM2D/trunk/su_create_Paul.python 2007-12-07 23:51:50 UTC (rev 8506)
@@ -1,60 +0,0 @@
-#!/usr/bin/env python
-#
-# Python code to generate Seismic Unix files for Ux and Uz and plot them
-#
-from os import *
-import os.path, string
-#
-def createSU():
- SEM=getcwd()
- filename=SEM+'/DATA/Par_file'
- #
- # Variables to be put in SU header
- #
- if path.exists(filename):
- variables=['nt','dt','sismostype']
- else:
- print 'No Par_file found !'
- return
- #
- # Open the file and get the lines
- #
- f = file(filename,'r')
- lignes= f.readlines()
- f.close()
- # Get the title
- for ligne in lignes:
- lsplit=string.split(ligne)
- if lsplit!= []:
- if lsplit[0]=='title':
- title=' '.join(lsplit[2:])
- break
- print '#'*50
- print '# SU file creation for '+title
- print '#'*50
- # Get the variables
- for var in variables:
- for ligne in lignes:
- lsplit=string.split(ligne)
- if lsplit!= []:
- if lsplit[0]==var:
- exec var+'='+string.replace(string.split(''.join(ligne))[2],'d','e')
- break
- #
- print sismostype
- chdir(SEM+'/OUTPUT_FILES')
- labels='label1="Time" label2="Receivers"'
- # Create headers and Su file
- if sismostype==4:
- ordres=['suaddhead < pressure_file_single.bin ns='+str(nt)+' | sushw key=dt a='+str(int(dt*1e6))+' > pressure_file.su']
- ordres.append('suxwigb < pressure_file.su perc=96 '+labels+' title=" Pressure : '+title+'"&')
- else:
- ordres=['suaddhead < Ux_file_single.bin ns='+str(nt)+' | sushw key=dt a='+str(int(dt*1e6))+' > Ux_file.su']
- ordres.append('suaddhead < Uz_file_single.bin ns='+str(nt)+' | sushw key=dt a='+str(int(dt*1e6))+' > Uz_file.su')
- ordres.append('suxwigb < Ux_file.su perc=96 '+labels+' title=" Ux : '+title+'"&')
- ordres.append('suxwigb < Uz_file.su xbox=600 perc=96 '+labels+' title=" Uz : '+title+'"&')
- #
- for i in range(len(ordres)):
- system(ordres[i])
-if __name__=='__main__':
- createSU()
More information about the cig-commits
mailing list