[cig-commits] r21168 - in seismo/2D/SPECFEM2D/trunk/EXAMPLES: . global_Earth_ak135f
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Fri Dec 14 17:15:45 PST 2012
Author: dkomati1
Date: 2012-12-14 17:15:45 -0800 (Fri, 14 Dec 2012)
New Revision: 21168
Added:
seismo/2D/SPECFEM2D/trunk/EXAMPLES/global_Earth_ak135f/
seismo/2D/SPECFEM2D/trunk/EXAMPLES/global_Earth_ak135f/Par_file
seismo/2D/SPECFEM2D/trunk/EXAMPLES/global_Earth_ak135f/SOURCE
seismo/2D/SPECFEM2D/trunk/EXAMPLES/global_Earth_ak135f/create_mesh_AK135F_2D_with_central_cube_no_PML.F90
seismo/2D/SPECFEM2D/trunk/EXAMPLES/global_Earth_ak135f/make_all_check.csh
seismo/2D/SPECFEM2D/trunk/EXAMPLES/global_Earth_ak135f/make_all_fast.csh
Log:
added EXAMPLES/global_Earth_ak135f
Added: seismo/2D/SPECFEM2D/trunk/EXAMPLES/global_Earth_ak135f/Par_file
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/global_Earth_ak135f/Par_file (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/global_Earth_ak135f/Par_file 2012-12-15 01:15:45 UTC (rev 21168)
@@ -0,0 +1,168 @@
+# title of job
+title = AK135F for half a disk for axisymmetric runs
+
+# forward or adjoint simulation
+SIMULATION_TYPE = 1 # 1 = forward, 2 = adjoint + kernels
+NOISE_TOMOGRAPHY = 0 # 0 = earthquake simulation, 1/2/3 = noise simulation
+SAVE_FORWARD = .false. # save the last frame, needed for adjoint simulation
+
+# parameters concerning partitioning
+nproc = 4 # number of processes
+partitioning_method = 3 # ascending order = 1, Scotch = 3
+PERFORM_CUTHILL_MCKEE = .false. # perform inverse Cuthill-McKee (1969) optimization/permutation for mesh numbering (can be very costly and not very useful)
+
+ngnod = 9 # number of control nodes per element (4 or 9)
+initialfield = .false. # use a plane wave as source or not
+add_Bielak_conditions = .false. # add Bielak conditions or not if initial plane wave
+assign_external_model = .true. # define external earth model or not
+READ_EXTERNAL_SEP_FILE = .false. # Read external SEP file from DATA/model_velocity.dat_input, or use routine
+ATTENUATION_VISCOELASTIC_SOLID = .false. # turn attenuation (viscoelasticity) on or off for non-poroelastic solid parts of the model
+ATTENUATION_PORO_FLUID_PART = .false. # turn viscous attenuation on or off for the fluid part of poroelastic parts of the model
+Q0 = 1 # quality factor for viscous attenuation
+freq0 = 10 # frequency for viscous attenuation
+p_sv = .true. # set the type of calculation (P-SV or SH/membrane waves)
+
+# time step parameters
+nt = 8000 # total number of time steps
+deltat = 0.25 # duration of a time step (see section "How to choose the time step" of the manual for how to do this)
+USER_T0 = 0.0d0 # use this t0 as earliest starting time rather than the automatically calculated one
+time_stepping_scheme = 1 # 1 = Newmark (2nd order), 2 = LDDRK4-6 (4th-order 6-stage low storage Runge-Kutta), 3 = classical 4th-order 4-stage Runge-Kutta
+
+# source parameters
+NSOURCES = 1 # number of sources [source info read in CMTSOLUTION file]
+force_normal_to_surface = .false. # angleforce normal to surface (external mesh and curve file needed)
+
+# constants for attenuation
+N_SLS = 2 # number of standard linear solids for attenuation
+f0_attenuation = 5.196152422706633 # (Hz) relevant only if source is a Dirac or a Heaviside, else it is f0
+
+# receiver set parameters for seismograms
+seismotype = 1 # record 1=displ 2=veloc 3=accel 4=pressure 5=curl of displ 6=the fluid potential
+NSTEP_BETWEEN_OUTPUT_SEISMOS = 5000000 # every how many time steps we save the seismograms (costly, do not use a very small value; if you use a very large value that is larger than the total number of time steps of the run, the seismograms will automatically be saved once at the end of the run anyway)
+save_ASCII_seismograms = .true. # save seismograms in ASCII format or not
+save_binary_seismograms_single = .false. # save seismograms in single precision binary format or not (can be used jointly with ASCII above to save both)
+save_binary_seismograms_double = .false. # save seismograms in double precision binary format or not (can be used jointly with both flags above to save all)
+SU_FORMAT = .false. # output single precision binary seismograms in Seismic Unix format (adjoint traces will be read in the same format)
+subsamp_seismos = 1 # subsampling of the seismograms to create smaller files (but less accurately sampled in time)
+generate_STATIONS = .true. # creates a new STATION file in ./DATA from this file or use the existing one
+nreceiversets = 1 # number of receiver sets
+anglerec = 0.d0 # angle to rotate components at receivers
+rec_normal_to_surface = .false. # replace anglerec with the normal to the surface (external mesh and curve file needed, anglerec above ignored)
+
+# first receiver set
+nrec = 2 # number of receivers
+xdeb = 6050000. # first receiver x in meters
+zdeb = 900000. # first receiver z in meters
+xfin = 6050000. # last receiver x in meters (ignored if onlyone receiver)
+zfin = -900000. # last receiver z in meters (ignored if onlyone receiver)
+enreg_surf_same_vertical = .false. # receivers inside the medium or at the surface
+
+# display parameters
+NSTEP_BETWEEN_OUTPUT_INFO = 200 # every how many time steps we display information about the simulation (costly, do not use a very small value)
+NSTEP_BETWEEN_OUTPUT_IMAGES = 200 # every how many time steps we draw JPEG or PostScript pictures of the simulation (costly, do not use a very small value)
+cutsnaps = 1. # minimum amplitude kept in % for the JPEG and PostScript snapshots; amplitudes below that are muted
+
+#### for JPEG color images ####
+output_color_image = .true. # output JPEG color image of the results every NSTEP_BETWEEN_OUTPUT_IMAGES time steps or not
+imagetype_JPEG = 2 # display 1=displ_Ux 2=displ_Uz 3=displ_norm 4=veloc_Vx 5=veloc_Vz 6=veloc_norm 7=accel_Ax 8=accel_Az 9=accel_norm 10=pressure
+factor_subsample_image = 1 # factor to subsample color images output by the code (useful for very large models)
+POWER_DISPLAY_COLOR = 0.30d0 # non linear display to enhance small amplitudes in color images
+DRAW_SOURCES_AND_RECEIVERS = .true. # display sources as orange crosses and receivers as green squares in JPEG images or not
+DRAW_WATER_IN_BLUE = .true. # display acoustic layers as constant blue in JPEG images, because they likely correspond to water in the case of ocean acoustics or in the case of offshore oil industry experiments (if off, display them as greyscale, as for elastic or poroelastic elements, for instance for acoustic-only oil industry models of solid media)
+USE_SNAPSHOT_NUMBER_IN_FILENAME = .false. # use snapshot number in the file name of JPEG color snapshots instead of the time step
+
+#### for PostScript snapshots ####
+output_postscript_snapshot = .false. # output Postscript snapshot of the results
+imagetype_postscript = 1 # display 1=displ vector 2=veloc vector 3=accel vector; small arrows are displayed for the vectors
+meshvect = .true. # display mesh on vector plots or not
+modelvect = .false. # display velocity model on vector plots
+boundvect = .true. # display boundary conditions on plots
+interpol = .true. # interpolation of the display or not
+pointsdisp = 6 # points for interpolation of display (set to 1 for lower-left corner only)
+subsamp_postscript = 1 # subsampling of color snapshots
+sizemax_arrows = 1.d0 # maximum size of arrows on vector plots in cm
+US_LETTER = .false. # US letter paper or European A4
+
+#### for wavefield dumps ####
+NSTEP_BETWEEN_OUTPUT_WAVE_DUMPS = 100 # every how many time steps we dump results of the simulation as ASCII or binary files (costly, do not use a very small value)
+output_wavefield_dumps = .false. # output wave field to a text file every NSTEP_BETWEEN_OUTPUT_TEXT_DUMPS time steps (creates very big files)
+imagetype_wavefield_dumps = 1 # display 1=displ vector 2=veloc vector 3=accel vector 4=pressure
+use_binary_for_wavefield_dumps = .false. # use ASCII or single-precision binary format for the wave field dumps
+####
+output_grid_Gnuplot = .false. # save the grid in a text file or not
+output_grid_ASCII = .false. # dump the grid in an ASCII text file consisting of a set of X,Y,Z points or not
+output_energy = .false. # compute and output acoustic and elastic energy (slows down the code significantly)
+
+# velocity and density models
+nbmodels = 4 # nb of different models
+# define models as I: (model_number,1,rho,Vp,Vs,0,0,QKappa Qmu) or II: (model_number,2,rho,c11,c13,c15,c33,c35,c55) or III: (model_number,3,rhos,rhof,phi,c,kxx,kxz,kzz,Ks,Kf,Kfr,etaf,mufr,Qmu).
+# For istropic elastic/acoustic material use I and set Vs to zero to make a given model acoustic, for anisotropic elastic use II,
+# and for isotropic poroelastic material use III. The mesh can contain acoustic, elastic, & poroelastic models simultaneously
+1 1 3320.d0 8040.d0 4480.0d0 0 0 9999.d0 9999.d0 0 0 0 0 0 0
+2 1 4985.d0 12503.d0 6805.0d0 0 0 9999.d0 9999.d0 0 0 0 0 0 0
+3 1 11539.d0 9751.d0 0.0d0 0 0 9999.d0 9999.d0 0 0 0 0 0 0
+4 1 13004.d0 11256.d0 3663.0d0 0 0 9999.d0 9999.d0 0 0 0 0 0 0
+#3 3 2650.d0 1000.0d0 0.4 1.25 1d-10 0.0 1d-10 3.6d10 2.25d9 2d9 0.0d-4 3.204d9 10.d0
+#2 1 2500.d0 5000.d0 2500.0d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+#3 1 2200.d0 2500.d0 1443.375d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+#4 3 2200.d0 786.3d0 0.4 2.0 1d-11 0.0 1d-11 5.341d9 2d9 3d9 0.0d-4 3.204d9 10.d0
+#5 2 2500.d0 169.d9 122.d9 0.d0 169.d9 0.d0 75.3d9 0 0 0 0 0 0
+#4 1 2200.d0 2200.d0 1343.375d0 0 0 10.d0 10.d0 0 0 0 0 0 0
+
+# external mesh or not
+read_external_mesh = .true.
+
+# absorbing boundary active or not
+PML_BOUNDARY_CONDITIONS = .false.
+NELEM_PML_THICKNESS = 3
+STACEY_ABSORBING_CONDITIONS = .false.
+ADD_SPRING_TO_STACEY = .false.
+
+# for horizontal periodic conditions: detect common points between left and right edges
+ADD_PERIODIC_CONDITIONS = .false.
+
+# horizontal periodicity distance for periodic conditions
+PERIODIC_horiz_dist = 0.3597d0
+
+# grid point detection tolerance for periodic conditions
+PERIODIC_DETECT_TOL = 3.3334d-6
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR EXTERNAL MESHING
+
+# data concerning mesh, when generated using third-party app (more info in README)
+# (see also absorbing_conditions above)
+mesh_file = ./DATA/Mesh_AK135F_NO_MUD # file containing the mesh
+nodes_coords_file = ./DATA/Nodes_AK135F_NO_MUD # file containing the nodes coordinates
+materials_file = ./DATA/Material_AK135F_NO_MUD # file containing the material number for each element
+free_surface_file = ./DATA/Surf_free_AK135F_NO_MUD # file containing the free surface
+absorbing_surface_file = ./DATA/Surf_abs_AK135F_NO_MUD # file containing the absorbing surface
+CPML_element_file = ./DATA/Elements_CPML_list # file containing the CPML element numbers
+tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the curve delimiting the velocity model
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR INTERNAL MESHING
+
+# file containing interfaces for internal mesh
+interfacesfile = Interface_flat_ASM_DGA_119_62kHz.dat
+
+# geometry of the model (origin lower-left corner = 0,0) and mesh description
+xmin = -1.d0 # abscissa of left side of the model
+xmax = 4.d0 # abscissa of right side of the model# file containing interfaces for internal mesh
+nx = 835 # number of elements along X
+
+# absorbing boundaries parameters
+absorbbottom = .true.
+absorbright = .true.
+absorbtop = .true.
+absorbleft = .true.
+
+# define the different regions of the model in the (nx,nz) spectral element mesh
+nbregions = 2 # nb of regions and model number for each
+1 835 1 71 1
+1 835 72 96 2
+#1 80 21 40 4
+#1 80 41 60 3
+#60 80 21 40 4
+#30 40 50 60 2
+#35 40 50 60 5
Added: seismo/2D/SPECFEM2D/trunk/EXAMPLES/global_Earth_ak135f/SOURCE
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/global_Earth_ak135f/SOURCE (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/global_Earth_ak135f/SOURCE 2012-12-15 01:15:45 UTC (rev 21168)
@@ -0,0 +1,13 @@
+#source 1. The components of a moment tensor source must be given in N.m, not in dyne.cm as in the DATA/CMTSOLUTION source file of the 3D version of the code.
+source_surf = .false. # source inside the medium or at the surface
+xs = 3500000. # source location x in meters
+zs = 5000000. # source location z in meters
+source_type = 1 # elastic force or acoustic pressure = 1 or moment tensor = 2
+time_function_type = 1 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4, Heaviside = 5
+f0 = 0.02d0 # dominant source frequency (Hz) if not Dirac or Heaviside
+t0 = 0.0 # time shift when multi sources (if one source, must be zero)
+anglesource = 0.0 # angle of the source (for a force only)
+Mxx = 1. # Mxx component (for a moment tensor source only)
+Mzz = 1. # Mzz component (for a moment tensor source only)
+Mxz = 0. # Mxz component (for a moment tensor source only)
+factor = 1.d10 # amplification factor
Added: seismo/2D/SPECFEM2D/trunk/EXAMPLES/global_Earth_ak135f/create_mesh_AK135F_2D_with_central_cube_no_PML.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/global_Earth_ak135f/create_mesh_AK135F_2D_with_central_cube_no_PML.F90 (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/global_Earth_ak135f/create_mesh_AK135F_2D_with_central_cube_no_PML.F90 2012-12-15 01:15:45 UTC (rev 21168)
@@ -0,0 +1,1449 @@
+
+ program generate_mesh_PREM
+
+ implicit none
+
+! dans le code SPECFEM3D_GLOBE on fait :
+ ! un doubling sous le Moho
+ ! un doubling un peu plus bas que sous la d670
+ ! un doubling un peu au dessus du milieu du outer core, entre la CMB et l'ICB
+
+! ici on fait un doubling de moins, et l'un des doublings n'est pas tout a fait au meme endroit :
+ ! un doubling juste sous la d670
+ ! un doubling un peu au dessus du milieu du outer core, entre la CMB et l'ICB
+
+! in the case of very large meshes, this option can be useful to switch from ASCII to binary for the mesh files
+#define USE_BINARY_FOR_EXTERNAL_MESH_DATABASE
+
+! 1 = generate the whole mesh 2 = generate only half the mesh
+ integer, parameter :: factor_divide_mesh = 2
+
+! multiplication factor for the mesh to be accurate at 1 second of minimum seismic period
+! 1 is accurate down to periods of about 45 seconds
+! 2 is accurate down to periods of about 22.5 seconds
+! 22 is accurate down to periods of about 2 seconds
+! 45 is accurate down to periods of about 1 second
+! 64 is accurate down to periods of about 0.7 second
+! 90 is accurate down to periods of about 0.5 second
+ integer, parameter :: multiplication_factor = 2 ! 1 ! 22 ! 45
+
+! output a Gnuplot grid or not (set to false when creating a very large mesh for very short periods)
+ logical, parameter :: output_gnuplot_grid = .false. ! .true.
+
+! The routine uses 9 control nodes defined as follows:
+!
+! 4 . . . . 7 . . . . 3
+! . .
+! . y .
+! . .
+! 8 9 x 6
+! . .
+! . .
+! . .
+! 1 . . . . 5 . . . . 2
+
+! nombre de dimensions du probleme physique
+ integer, parameter :: NDIM = 2
+
+! nombre de noeuds de controle par element
+ integer, parameter :: ngnod = 9
+
+! double du nombre d'elements a la surface
+! double du nombre d'elements spectraux dans la direction radiale
+ integer, parameter :: nspec_surf_whole_circle = 512 * multiplication_factor
+ integer, parameter :: nspec_rad_670_surf = 10 * multiplication_factor
+ integer, parameter :: nspec_rad_CMB_670 = 56 * multiplication_factor
+ integer, parameter :: nspec_rad_doubling_OC_to_CMB = 24 * multiplication_factor
+ integer, parameter :: nspec_rad_ICB_to_doubling_OC = 40 * multiplication_factor
+ integer, parameter :: nspec_rad_Cube_ICB = 16 * multiplication_factor
+
+! taille du cube a l'interieur de la graine
+ double precision, parameter :: radius_cube = 650000.d0
+
+! to inflate the central cube (set to 0.d0 for a non-inflated cube)
+ double precision, parameter :: CENTRAL_CUBE_INFLATE_FACTOR = 0.41d0
+
+! flags for material properties
+ integer, parameter :: IREGION_MANTLE_CRUST_ABOVE_d670 = 1
+ integer, parameter :: IREGION_MANTLE_BELOW_d670 = 2
+ integer, parameter :: IREGION_OUTER_CORE = 3
+ integer, parameter :: IREGION_INNER_CORE = 4
+
+! flags for symmetry-axis edges
+ integer, parameter :: IBOTTOM = 1
+ integer, parameter :: IRIGHT = 2
+ integer, parameter :: ITOP = 3
+ integer, parameter :: ILEFT = 4
+
+ double precision :: ratio_x, ratio_y, fact_x, fact_y, xi, eta
+
+! prevision exacte du nombre max d'elements et prevision grossiere du nombre max de points
+ integer, parameter :: nspec_max_670_surf = (nspec_surf_whole_circle/2)*(nspec_rad_670_surf/2)
+ integer, parameter :: nspec_max_CMB_670 = (nspec_surf_whole_circle/4)*(nspec_rad_CMB_670/4 - 1) + &
+ (nspec_surf_whole_circle/4)*3
+ integer, parameter :: nspec_max_doubling_OC_CMB = (nspec_surf_whole_circle/4)*(nspec_rad_doubling_OC_to_CMB/4)
+ integer, parameter :: nspec_max_ICB_doubling_OC = (nspec_surf_whole_circle/8)*(nspec_rad_ICB_to_doubling_OC/4 - 1) + &
+ (nspec_surf_whole_circle/8)*3
+ integer, parameter :: nspec_max_Cube_ICB = (nspec_surf_whole_circle/8)*(nspec_rad_Cube_ICB/4)
+ integer, parameter :: nspec_max_inside_cube = (nspec_surf_whole_circle/32)*(nspec_surf_whole_circle/32)
+
+ integer, parameter :: nspec_exact = (nspec_max_670_surf + nspec_max_CMB_670 + nspec_max_doubling_OC_CMB + &
+ nspec_max_ICB_doubling_OC + nspec_max_Cube_ICB + nspec_max_inside_cube) / factor_divide_mesh
+
+! prevision surestimee du nombre max de points pour routine Fischer
+ integer, parameter :: npoin_max = nspec_exact * ngnod
+
+! topologie des elements du maillage et coordonnees des points de controle
+ integer, dimension(ngnod,nspec_exact) :: ibool
+ integer, dimension(nspec_exact) :: idoubling
+ double precision, dimension(ngnod,nspec_exact) :: xcoord,ycoord
+
+ double precision, dimension(npoin_max) :: xp,yp,work
+ integer, dimension(npoin_max) :: loc,ind,ninseg,iglob
+ logical, dimension(npoin_max) :: ifseg
+
+! save memory in the case of very large meshes to generate by using two different shapes for the same array
+ equivalence(ibool,iglob)
+
+ double precision, parameter :: R_EARTH = 6371000.d0
+
+! values for AK135F_NO_MUD (modified version of AK135)
+ double precision, parameter :: R670 = 5711000.d0
+ double precision, parameter :: RCMB = 3479500.d0
+ double precision, parameter :: RICB = 1217500.d0
+ double precision, parameter :: R_DOUBLING_OUTER_CORE = RICB + 0.56*(RCMB - RICB)
+
+ double precision, parameter :: PI = 3.141592653589793d0
+ double precision, parameter :: PI_OVER_TWO = PI / 2.d0
+
+ double precision zero,one
+ parameter(zero=0.d0, one=1.d0)
+
+! to be able to estimate the shortest seismic period that can be resolved by the mesh
+ double precision, parameter :: cs_min_in_crust_in_km_per_s = 3.460000d0
+
+ logical :: generate_only_half_the_mesh
+ double precision :: delta_theta,size_of_a_surface_element_in_km,radius_interf
+ integer :: nspec,ispec,ipoin,ia1,ia2,icentral_cube1,icentral_cube2,ispec_count
+
+! %%%%% pour le maillage elements spectraux %%%%%
+
+!---- zone d670 -> surface
+
+ double precision x1surf(0:nspec_surf_whole_circle)
+ double precision y1surf(0:nspec_surf_whole_circle)
+
+ double precision x1bot(0:nspec_surf_whole_circle)
+ double precision y1bot(0:nspec_surf_whole_circle)
+
+ double precision x1vol(0:nspec_surf_whole_circle,0:nspec_rad_670_surf)
+ double precision y1vol(0:nspec_surf_whole_circle,0:nspec_rad_670_surf)
+
+!---- zone CMB -> d670
+
+ double precision x2surf(0:nspec_surf_whole_circle)
+ double precision y2surf(0:nspec_surf_whole_circle)
+
+ double precision x2bot(0:nspec_surf_whole_circle)
+ double precision y2bot(0:nspec_surf_whole_circle)
+
+ double precision x2vol(0:nspec_surf_whole_circle,0:nspec_rad_CMB_670)
+ double precision y2vol(0:nspec_surf_whole_circle,0:nspec_rad_CMB_670)
+
+!---- zone ICB -> CMB
+
+ double precision x3surf(0:nspec_surf_whole_circle)
+ double precision y3surf(0:nspec_surf_whole_circle)
+
+ double precision x3bot(0:nspec_surf_whole_circle)
+ double precision y3bot(0:nspec_surf_whole_circle)
+
+ double precision x3vol(0:nspec_surf_whole_circle,0:nspec_rad_ICB_to_doubling_OC)
+ double precision y3vol(0:nspec_surf_whole_circle,0:nspec_rad_ICB_to_doubling_OC)
+
+!---- zone Cube -> ICB
+
+ double precision x4surf(0:nspec_surf_whole_circle)
+ double precision y4surf(0:nspec_surf_whole_circle)
+
+ double precision x4bot(0:nspec_surf_whole_circle)
+ double precision y4bot(0:nspec_surf_whole_circle)
+
+ double precision x4vol(0:nspec_surf_whole_circle,0:nspec_rad_Cube_ICB)
+ double precision y4vol(0:nspec_surf_whole_circle,0:nspec_rad_Cube_ICB)
+
+!---- zone interieur du Cube
+
+ double precision x5vol(0:nspec_surf_whole_circle/16,0:nspec_surf_whole_circle/16)
+ double precision y5vol(0:nspec_surf_whole_circle/16,0:nspec_surf_whole_circle/16)
+
+! save a lot of memory in the case of a very dense mesh by using the same memory block for the five arrays
+ equivalence(x1vol,x2vol,x3vol,x4vol,x5vol)
+ equivalence(y1vol,y2vol,y3vol,y4vol,y5vol)
+
+ integer ix,irad,ia
+ double precision xicoord,radcoord,thetacoord,xlincoord
+
+! %%%%% pour le maillage elements spectraux %%%%%
+
+! %%%%% pour routine numbering Fischer %%%%%
+! pour numerotation locale -> globale
+ integer ie,nseg,ioff,iseg,ig,ieoff
+ integer ilocnum,i,j,npoin
+
+! mesh tolerance for fast global numbering
+ double precision, parameter :: SMALLVALTOL = 1.d-10
+
+! very large and very small values
+ double precision, parameter :: HUGEVAL = 1.d+30
+
+ double precision xmaxval,xminval,ymaxval,yminval,xtol,xtypdist
+
+! elements 9 noeuds donc on a seulement la moitie des elements
+ delta_theta = 2. * pi / dble(nspec_surf_whole_circle/2)
+ size_of_a_surface_element_in_km = delta_theta*R_EARTH/1000.
+
+ print *
+ print *,'Nombre d''elements en surface du cercle complet = ',nspec_surf_whole_circle
+ print *
+ print *,'Taille d''un element spectral en surface = ',size_of_a_surface_element_in_km,' km'
+ print *
+ print *,'Since minimum S velocity in the crust is about ',cs_min_in_crust_in_km_per_s,' km/s'
+ print *,'and since one spectral element is needed to resolve the shortest wavelength,'
+ print *,'this mesh is accurate approximately down to a minimum period of ', &
+ size_of_a_surface_element_in_km/cs_min_in_crust_in_km_per_s,' seconds'
+ print *
+
+! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+!
+!---- generation de la grille pour SPECFEM90
+!
+
+ print *,'Generation de la grille pour Specfem...'
+
+! generate only half the mesh or the whole mesh
+ if(factor_divide_mesh < 1 .or. factor_divide_mesh > 2) stop 'incorrect value of factor_divide_mesh'
+
+ if(factor_divide_mesh == 1) then
+ generate_only_half_the_mesh = .false.
+ else
+ generate_only_half_the_mesh = .true.
+ endif
+
+ if(generate_only_half_the_mesh) then
+ print *,'generating only half the mesh'
+ else
+ print *,'generating the whole mesh'
+ endif
+ print *
+
+! verification de la coherence des parametres entres
+
+! --- pour pouvoir generer des elements 9 noeuds
+ if(mod(nspec_rad_670_surf,2) /= 0) stop 'nspec_rad_670_surf should be a multiple of 2'
+ if(mod(nspec_rad_CMB_670,4) /= 0) stop 'nspec_rad_CMB_670 should be a multiple of 4'
+ if(mod(nspec_rad_ICB_to_doubling_OC,4) /= 0) stop 'nspec_rad_ICB_to_doubling_OC should be a multiple of 4'
+ if(mod(nspec_rad_Cube_ICB,4) /= 0) stop 'nspec_rad_Cube_ICB should be a multiple of 4'
+
+! --- pour permettre le deraffinement
+ if(mod(nspec_surf_whole_circle,16) /= 0) stop 'nspec_surf_whole_circle should be a multiple of 16'
+ if(nspec_rad_CMB_670 < 6) stop 'nspec_rad_CMB_670 should be greater than 6'
+ if(nspec_rad_ICB_to_doubling_OC < 6) stop 'nspec_rad_ICB_to_doubling_OC should be greater than 6'
+
+! %%%% zone d670 -> surface %%%%
+
+! generation maillage de la surface
+
+ do ix=0,nspec_surf_whole_circle
+
+ xicoord = dble(ix)/dble(nspec_surf_whole_circle)
+
+! parcours de l'angle dans le sens negatif
+ thetacoord = +pi/2 - 2*pi*xicoord
+
+! coordonnees cartesiennes correspondantes
+
+!---- surface
+
+ radius_interf = R_EARTH
+
+ x1surf(ix) = radius_interf * dcos(thetacoord)
+ y1surf(ix) = radius_interf * dsin(thetacoord)
+
+!---- d670
+
+ radius_interf = R670
+
+ x1bot(ix) = radius_interf * dcos(thetacoord)
+ y1bot(ix) = radius_interf * dsin(thetacoord)
+
+!---- volume
+
+ do irad=0,nspec_rad_670_surf
+ radcoord = dble(irad)/dble(nspec_rad_670_surf)
+ x1vol(ix,irad) = x1surf(ix) * radcoord + x1bot(ix) * (one - radcoord)
+ y1vol(ix,irad) = y1surf(ix) * radcoord + y1bot(ix) * (one - radcoord)
+ enddo
+
+ enddo
+
+! %%% ecrire la grille dans un tableau pour routine de Fischer
+ ispec = 0
+
+! %%% bloc d670 -> surface
+ do irad=0,nspec_rad_670_surf-2,2
+ do ix=0,nspec_surf_whole_circle/factor_divide_mesh-2,2
+
+ ispec = ispec + 1
+
+ idoubling(ispec) = IREGION_MANTLE_CRUST_ABOVE_d670
+
+ xcoord(1,ispec) = x1vol(ix ,irad)
+ xcoord(5,ispec) = x1vol(ix+1,irad)
+ xcoord(2,ispec) = x1vol(ix+2,irad)
+ xcoord(8,ispec) = x1vol(ix ,irad+1)
+ xcoord(9,ispec) = x1vol(ix+1,irad+1)
+ xcoord(6,ispec) = x1vol(ix+2,irad+1)
+ xcoord(4,ispec) = x1vol(ix ,irad+2)
+ xcoord(7,ispec) = x1vol(ix+1,irad+2)
+ xcoord(3,ispec) = x1vol(ix+2,irad+2)
+
+ ycoord(1,ispec) = y1vol(ix ,irad)
+ ycoord(5,ispec) = y1vol(ix+1,irad)
+ ycoord(2,ispec) = y1vol(ix+2,irad)
+ ycoord(8,ispec) = y1vol(ix ,irad+1)
+ ycoord(9,ispec) = y1vol(ix+1,irad+1)
+ ycoord(6,ispec) = y1vol(ix+2,irad+1)
+ ycoord(4,ispec) = y1vol(ix ,irad+2)
+ ycoord(7,ispec) = y1vol(ix+1,irad+2)
+ ycoord(3,ispec) = y1vol(ix+2,irad+2)
+ enddo
+ enddo
+
+! %%%% zone CMB -> d670 %%%%
+
+! generation maillage de la surface
+
+ do ix=0,nspec_surf_whole_circle
+
+ xicoord = dble(ix)/dble(nspec_surf_whole_circle)
+
+! parcours de l'angle dans le sens negatif
+ thetacoord = +pi/2 - 2*pi*xicoord
+
+! coordonnees cartesiennes correspondantes
+
+!---- d670
+
+ radius_interf = R670
+
+ x2surf(ix) = radius_interf * dcos(thetacoord)
+ y2surf(ix) = radius_interf * dsin(thetacoord)
+
+!---- CMB
+
+ radius_interf = RCMB
+
+ x2bot(ix) = radius_interf * dcos(thetacoord)
+ y2bot(ix) = radius_interf * dsin(thetacoord)
+
+!---- volume
+
+ do irad=0,nspec_rad_CMB_670
+ radcoord = dble(irad)/dble(nspec_rad_CMB_670)
+ x2vol(ix,irad) = x2surf(ix) * radcoord + x2bot(ix) * (one - radcoord)
+ y2vol(ix,irad) = y2surf(ix) * radcoord + y2bot(ix) * (one - radcoord)
+ enddo
+
+ enddo
+
+! --- bloc principal
+ do irad=0,nspec_rad_CMB_670-8,4
+ do ix=0,nspec_surf_whole_circle/factor_divide_mesh-4,4
+
+ ispec = ispec + 1
+
+ idoubling(ispec) = IREGION_MANTLE_BELOW_d670
+
+ xcoord(1,ispec) = x2vol(ix ,irad)
+ xcoord(5,ispec) = x2vol(ix+2,irad)
+ xcoord(2,ispec) = x2vol(ix+4,irad)
+ xcoord(8,ispec) = x2vol(ix ,irad+2)
+ xcoord(9,ispec) = x2vol(ix+2,irad+2)
+ xcoord(6,ispec) = x2vol(ix+4,irad+2)
+ xcoord(4,ispec) = x2vol(ix ,irad+4)
+ xcoord(7,ispec) = x2vol(ix+2,irad+4)
+ xcoord(3,ispec) = x2vol(ix+4,irad+4)
+
+ ycoord(1,ispec) = y2vol(ix ,irad)
+ ycoord(5,ispec) = y2vol(ix+2,irad)
+ ycoord(2,ispec) = y2vol(ix+4,irad)
+ ycoord(8,ispec) = y2vol(ix ,irad+2)
+ ycoord(9,ispec) = y2vol(ix+2,irad+2)
+ ycoord(6,ispec) = y2vol(ix+4,irad+2)
+ ycoord(4,ispec) = y2vol(ix ,irad+4)
+ ycoord(7,ispec) = y2vol(ix+2,irad+4)
+ ycoord(3,ispec) = y2vol(ix+4,irad+4)
+ enddo
+ enddo
+
+! --- zone de raccord geometrique conforme
+ irad=nspec_rad_CMB_670-4
+ do ix=0,nspec_surf_whole_circle/factor_divide_mesh-8,8
+
+ ispec = ispec + 1
+
+ idoubling(ispec) = IREGION_MANTLE_BELOW_d670
+
+ xcoord(1,ispec) = x2vol(ix ,irad+2)
+ xcoord(5,ispec) = x2vol(ix+1,irad+2)
+ xcoord(2,ispec) = x2vol(ix+2,irad+2)
+ xcoord(8,ispec) = x2vol(ix ,irad+3)
+ xcoord(9,ispec) = x2vol(ix+1,irad+3)
+ xcoord(6,ispec) = x2vol(ix+2,irad+3)
+ xcoord(4,ispec) = x2vol(ix ,irad+4)
+ xcoord(7,ispec) = x2vol(ix+1,irad+4)
+ xcoord(3,ispec) = x2vol(ix+2,irad+4)
+
+ ycoord(1,ispec) = y2vol(ix ,irad+2)
+ ycoord(5,ispec) = y2vol(ix+1,irad+2)
+ ycoord(2,ispec) = y2vol(ix+2,irad+2)
+ ycoord(8,ispec) = y2vol(ix ,irad+3)
+ ycoord(9,ispec) = y2vol(ix+1,irad+3)
+ ycoord(6,ispec) = y2vol(ix+2,irad+3)
+ ycoord(4,ispec) = y2vol(ix ,irad+4)
+ ycoord(7,ispec) = y2vol(ix+1,irad+4)
+ ycoord(3,ispec) = y2vol(ix+2,irad+4)
+
+ ispec = ispec + 1
+
+ idoubling(ispec) = IREGION_MANTLE_BELOW_d670
+
+ xcoord(1,ispec) = x2vol(ix ,irad)
+ xcoord(5,ispec) = x2vol(ix+2,irad)
+ xcoord(2,ispec) = x2vol(ix+4,irad)
+ xcoord(8,ispec) = x2vol(ix ,irad+1)
+ xcoord(9,ispec) = (x2vol(ix+1,irad+1) + x2vol(ix+2,irad+1))/2.
+ xcoord(6,ispec) = x2vol(ix+3,irad+1)
+ xcoord(4,ispec) = x2vol(ix ,irad+2)
+ xcoord(7,ispec) = x2vol(ix+1,irad+2)
+ xcoord(3,ispec) = x2vol(ix+2,irad+2)
+
+ ycoord(1,ispec) = y2vol(ix ,irad)
+ ycoord(5,ispec) = y2vol(ix+2,irad)
+ ycoord(2,ispec) = y2vol(ix+4,irad)
+ ycoord(8,ispec) = y2vol(ix ,irad+1)
+ ycoord(9,ispec) = (y2vol(ix+1,irad+1) + y2vol(ix+2,irad+1))/2.
+ ycoord(6,ispec) = y2vol(ix+3,irad+1)
+ ycoord(4,ispec) = y2vol(ix ,irad+2)
+ ycoord(7,ispec) = y2vol(ix+1,irad+2)
+ ycoord(3,ispec) = y2vol(ix+2,irad+2)
+
+ ispec = ispec + 1
+
+ idoubling(ispec) = IREGION_MANTLE_BELOW_d670
+
+ xcoord(1,ispec) = x2vol(ix+2,irad+2)
+ xcoord(5,ispec) = x2vol(ix+3,irad+1)
+ xcoord(2,ispec) = x2vol(ix+4,irad)
+ xcoord(8,ispec) = x2vol(ix+2,irad+3)
+ xcoord(9,ispec) = (x2vol(ix+3,irad+2) + x2vol(ix+3,irad+3))/2.
+ xcoord(6,ispec) = x2vol(ix+4,irad+2)
+ xcoord(4,ispec) = x2vol(ix+2,irad+4)
+ xcoord(7,ispec) = x2vol(ix+3,irad+4)
+ xcoord(3,ispec) = x2vol(ix+4,irad+4)
+
+ ycoord(1,ispec) = y2vol(ix+2,irad+2)
+ ycoord(5,ispec) = y2vol(ix+3,irad+1)
+ ycoord(2,ispec) = y2vol(ix+4,irad)
+ ycoord(8,ispec) = y2vol(ix+2,irad+3)
+ ycoord(9,ispec) = (y2vol(ix+3,irad+2) + y2vol(ix+3,irad+3))/2.
+ ycoord(6,ispec) = y2vol(ix+4,irad+2)
+ ycoord(4,ispec) = y2vol(ix+2,irad+4)
+ ycoord(7,ispec) = y2vol(ix+3,irad+4)
+ ycoord(3,ispec) = y2vol(ix+4,irad+4)
+
+ enddo
+
+! --- zone de raccord geometrique conforme inverse
+ irad=nspec_rad_CMB_670-4
+ do ix=4,nspec_surf_whole_circle/factor_divide_mesh-4,8
+
+ ispec = ispec + 1
+
+ idoubling(ispec) = IREGION_MANTLE_BELOW_d670
+
+ xcoord(1,ispec) = x2vol(ix+2,irad+2)
+ xcoord(5,ispec) = x2vol(ix+3,irad+2)
+ xcoord(2,ispec) = x2vol(ix+4,irad+2)
+ xcoord(8,ispec) = x2vol(ix+2,irad+3)
+ xcoord(9,ispec) = x2vol(ix+3,irad+3)
+ xcoord(6,ispec) = x2vol(ix+4,irad+3)
+ xcoord(4,ispec) = x2vol(ix+2,irad+4)
+ xcoord(7,ispec) = x2vol(ix+3,irad+4)
+ xcoord(3,ispec) = x2vol(ix+4,irad+4)
+
+ ycoord(1,ispec) = y2vol(ix+2,irad+2)
+ ycoord(5,ispec) = y2vol(ix+3,irad+2)
+ ycoord(2,ispec) = y2vol(ix+4,irad+2)
+ ycoord(8,ispec) = y2vol(ix+2,irad+3)
+ ycoord(9,ispec) = y2vol(ix+3,irad+3)
+ ycoord(6,ispec) = y2vol(ix+4,irad+3)
+ ycoord(4,ispec) = y2vol(ix+2,irad+4)
+ ycoord(7,ispec) = y2vol(ix+3,irad+4)
+ ycoord(3,ispec) = y2vol(ix+4,irad+4)
+
+ ispec = ispec + 1
+
+ idoubling(ispec) = IREGION_MANTLE_BELOW_d670
+
+ xcoord(1,ispec) = x2vol(ix ,irad)
+ xcoord(5,ispec) = x2vol(ix+2,irad)
+ xcoord(2,ispec) = x2vol(ix+4,irad)
+ xcoord(8,ispec) = x2vol(ix+1,irad+1)
+ xcoord(9,ispec) = (x2vol(ix+2,irad+1) + x2vol(ix+3,irad+1))/2.
+ xcoord(6,ispec) = x2vol(ix+4,irad+1)
+ xcoord(4,ispec) = x2vol(ix+2,irad+2)
+ xcoord(7,ispec) = x2vol(ix+3,irad+2)
+ xcoord(3,ispec) = x2vol(ix+4,irad+2)
+
+ ycoord(1,ispec) = y2vol(ix ,irad)
+ ycoord(5,ispec) = y2vol(ix+2,irad)
+ ycoord(2,ispec) = y2vol(ix+4,irad)
+ ycoord(8,ispec) = y2vol(ix+1,irad+1)
+ ycoord(9,ispec) = (y2vol(ix+2,irad+1) + y2vol(ix+3,irad+1))/2.
+ ycoord(6,ispec) = y2vol(ix+4,irad+1)
+ ycoord(4,ispec) = y2vol(ix+2,irad+2)
+ ycoord(7,ispec) = y2vol(ix+3,irad+2)
+ ycoord(3,ispec) = y2vol(ix+4,irad+2)
+
+ ispec = ispec + 1
+
+ idoubling(ispec) = IREGION_MANTLE_BELOW_d670
+
+ xcoord(1,ispec) = x2vol(ix,irad)
+ xcoord(5,ispec) = x2vol(ix+1,irad+1)
+ xcoord(2,ispec) = x2vol(ix+2,irad+2)
+ xcoord(8,ispec) = x2vol(ix,irad+2)
+ xcoord(9,ispec) = (x2vol(ix+1,irad+2) + x2vol(ix+1,irad+3))/2.
+ xcoord(6,ispec) = x2vol(ix+2,irad+3)
+ xcoord(4,ispec) = x2vol(ix,irad+4)
+ xcoord(7,ispec) = x2vol(ix+1,irad+4)
+ xcoord(3,ispec) = x2vol(ix+2,irad+4)
+
+ ycoord(1,ispec) = y2vol(ix,irad)
+ ycoord(5,ispec) = y2vol(ix+1,irad+1)
+ ycoord(2,ispec) = y2vol(ix+2,irad+2)
+ ycoord(8,ispec) = y2vol(ix,irad+2)
+ ycoord(9,ispec) = (y2vol(ix+1,irad+2) + y2vol(ix+1,irad+3))/2.
+ ycoord(6,ispec) = y2vol(ix+2,irad+3)
+ ycoord(4,ispec) = y2vol(ix,irad+4)
+ ycoord(7,ispec) = y2vol(ix+1,irad+4)
+ ycoord(3,ispec) = y2vol(ix+2,irad+4)
+
+ enddo
+
+! %%%% zone doubling in the CMB -> CMB %%%%
+
+! generation maillage de la surface
+
+ do ix=0,nspec_surf_whole_circle
+
+ xicoord = dble(ix)/dble(nspec_surf_whole_circle)
+
+! parcours de l'angle dans le sens negatif
+ thetacoord = +pi/2 - 2*pi*xicoord
+
+! coordonnees cartesiennes correspondantes
+
+!---- CMB
+
+ radius_interf = RCMB
+
+ x3surf(ix) = radius_interf * dcos(thetacoord)
+ y3surf(ix) = radius_interf * dsin(thetacoord)
+
+!---- doubling in the outer core
+
+ radius_interf = R_DOUBLING_OUTER_CORE
+
+ x3bot(ix) = radius_interf * dcos(thetacoord)
+ y3bot(ix) = radius_interf * dsin(thetacoord)
+
+!---- volume
+
+ do irad=0,nspec_rad_doubling_OC_to_CMB
+ radcoord = dble(irad)/dble(nspec_rad_doubling_OC_to_CMB)
+ x3vol(ix,irad) = x3surf(ix) * radcoord + x3bot(ix) * (one - radcoord)
+ y3vol(ix,irad) = y3surf(ix) * radcoord + y3bot(ix) * (one - radcoord)
+ enddo
+
+ enddo
+
+! --- bloc principal
+ do irad=0,nspec_rad_doubling_OC_to_CMB-4,4
+ do ix=0,nspec_surf_whole_circle/factor_divide_mesh-4,4
+
+ ispec = ispec + 1
+
+ idoubling(ispec) = IREGION_OUTER_CORE
+
+ xcoord(1,ispec) = x3vol(ix ,irad)
+ xcoord(5,ispec) = x3vol(ix+2,irad)
+ xcoord(2,ispec) = x3vol(ix+4,irad)
+ xcoord(8,ispec) = x3vol(ix ,irad+2)
+ xcoord(9,ispec) = x3vol(ix+2,irad+2)
+ xcoord(6,ispec) = x3vol(ix+4,irad+2)
+ xcoord(4,ispec) = x3vol(ix ,irad+4)
+ xcoord(7,ispec) = x3vol(ix+2,irad+4)
+ xcoord(3,ispec) = x3vol(ix+4,irad+4)
+
+ ycoord(1,ispec) = y3vol(ix ,irad)
+ ycoord(5,ispec) = y3vol(ix+2,irad)
+ ycoord(2,ispec) = y3vol(ix+4,irad)
+ ycoord(8,ispec) = y3vol(ix ,irad+2)
+ ycoord(9,ispec) = y3vol(ix+2,irad+2)
+ ycoord(6,ispec) = y3vol(ix+4,irad+2)
+ ycoord(4,ispec) = y3vol(ix ,irad+4)
+ ycoord(7,ispec) = y3vol(ix+2,irad+4)
+ ycoord(3,ispec) = y3vol(ix+4,irad+4)
+ enddo
+ enddo
+
+! %%%% zone ICB -> doubling in the CMB %%%%
+
+! generation maillage de la surface
+
+ do ix=0,nspec_surf_whole_circle
+
+ xicoord = dble(ix)/dble(nspec_surf_whole_circle)
+
+! parcours de l'angle dans le sens negatif
+ thetacoord = +pi/2 - 2*pi*xicoord
+
+! coordonnees cartesiennes correspondantes
+
+!---- doubling in the outer core
+
+ radius_interf = R_DOUBLING_OUTER_CORE
+
+ x3surf(ix) = radius_interf * dcos(thetacoord)
+ y3surf(ix) = radius_interf * dsin(thetacoord)
+
+!---- ICB
+
+ radius_interf = RICB
+
+ x3bot(ix) = radius_interf * dcos(thetacoord)
+ y3bot(ix) = radius_interf * dsin(thetacoord)
+
+!---- volume
+
+ do irad=0,nspec_rad_ICB_to_doubling_OC
+ radcoord = dble(irad)/dble(nspec_rad_ICB_to_doubling_OC)
+ x3vol(ix,irad) = x3surf(ix) * radcoord + x3bot(ix) * (one - radcoord)
+ y3vol(ix,irad) = y3surf(ix) * radcoord + y3bot(ix) * (one - radcoord)
+ enddo
+
+ enddo
+
+! --- bloc principal
+ do irad=0,nspec_rad_ICB_to_doubling_OC-8,4
+ do ix=0,nspec_surf_whole_circle/factor_divide_mesh-8,8
+
+ ispec = ispec + 1
+
+ idoubling(ispec) = IREGION_OUTER_CORE
+
+ xcoord(1,ispec) = x3vol(ix ,irad)
+ xcoord(5,ispec) = x3vol(ix+4,irad)
+ xcoord(2,ispec) = x3vol(ix+8,irad)
+ xcoord(8,ispec) = x3vol(ix ,irad+2)
+ xcoord(9,ispec) = x3vol(ix+4,irad+2)
+ xcoord(6,ispec) = x3vol(ix+8,irad+2)
+ xcoord(4,ispec) = x3vol(ix ,irad+4)
+ xcoord(7,ispec) = x3vol(ix+4,irad+4)
+ xcoord(3,ispec) = x3vol(ix+8,irad+4)
+
+ ycoord(1,ispec) = y3vol(ix ,irad)
+ ycoord(5,ispec) = y3vol(ix+4,irad)
+ ycoord(2,ispec) = y3vol(ix+8,irad)
+ ycoord(8,ispec) = y3vol(ix ,irad+2)
+ ycoord(9,ispec) = y3vol(ix+4,irad+2)
+ ycoord(6,ispec) = y3vol(ix+8,irad+2)
+ ycoord(4,ispec) = y3vol(ix ,irad+4)
+ ycoord(7,ispec) = y3vol(ix+4,irad+4)
+ ycoord(3,ispec) = y3vol(ix+8,irad+4)
+ enddo
+ enddo
+
+! --- zone de raccord geometrique conforme
+ irad=nspec_rad_ICB_to_doubling_OC-4
+ do ix=0,nspec_surf_whole_circle/factor_divide_mesh-16,16
+
+ ispec = ispec + 1
+
+ idoubling(ispec) = IREGION_OUTER_CORE
+
+ xcoord(1,ispec) = x3vol(ix ,irad+2)
+ xcoord(5,ispec) = x3vol(ix+2,irad+2)
+ xcoord(2,ispec) = x3vol(ix+4,irad+2)
+ xcoord(8,ispec) = x3vol(ix ,irad+3)
+ xcoord(9,ispec) = x3vol(ix+2,irad+3)
+ xcoord(6,ispec) = x3vol(ix+4,irad+3)
+ xcoord(4,ispec) = x3vol(ix ,irad+4)
+ xcoord(7,ispec) = x3vol(ix+2,irad+4)
+ xcoord(3,ispec) = x3vol(ix+4,irad+4)
+
+ ycoord(1,ispec) = y3vol(ix ,irad+2)
+ ycoord(5,ispec) = y3vol(ix+2,irad+2)
+ ycoord(2,ispec) = y3vol(ix+4,irad+2)
+ ycoord(8,ispec) = y3vol(ix ,irad+3)
+ ycoord(9,ispec) = y3vol(ix+2,irad+3)
+ ycoord(6,ispec) = y3vol(ix+4,irad+3)
+ ycoord(4,ispec) = y3vol(ix ,irad+4)
+ ycoord(7,ispec) = y3vol(ix+2,irad+4)
+ ycoord(3,ispec) = y3vol(ix+4,irad+4)
+
+ ispec = ispec + 1
+
+ idoubling(ispec) = IREGION_OUTER_CORE
+
+ xcoord(1,ispec) = x3vol(ix ,irad)
+ xcoord(5,ispec) = x3vol(ix+4,irad)
+ xcoord(2,ispec) = x3vol(ix+8,irad)
+ xcoord(8,ispec) = x3vol(ix ,irad+1)
+ xcoord(9,ispec) = (x3vol(ix+2,irad+1) + x3vol(ix+4,irad+1))/2.
+ xcoord(6,ispec) = x3vol(ix+6,irad+1)
+ xcoord(4,ispec) = x3vol(ix ,irad+2)
+ xcoord(7,ispec) = x3vol(ix+2,irad+2)
+ xcoord(3,ispec) = x3vol(ix+4,irad+2)
+
+ ycoord(1,ispec) = y3vol(ix ,irad)
+ ycoord(5,ispec) = y3vol(ix+4,irad)
+ ycoord(2,ispec) = y3vol(ix+8,irad)
+ ycoord(8,ispec) = y3vol(ix ,irad+1)
+ ycoord(9,ispec) = (y3vol(ix+2,irad+1) + y3vol(ix+4,irad+1))/2.
+ ycoord(6,ispec) = y3vol(ix+6,irad+1)
+ ycoord(4,ispec) = y3vol(ix ,irad+2)
+ ycoord(7,ispec) = y3vol(ix+2,irad+2)
+ ycoord(3,ispec) = y3vol(ix+4,irad+2)
+
+ ispec = ispec + 1
+
+ idoubling(ispec) = IREGION_OUTER_CORE
+
+ xcoord(1,ispec) = x3vol(ix+4,irad+2)
+ xcoord(5,ispec) = x3vol(ix+6,irad+1)
+ xcoord(2,ispec) = x3vol(ix+8,irad)
+ xcoord(8,ispec) = x3vol(ix+4,irad+3)
+ xcoord(9,ispec) = (x3vol(ix+6,irad+2) + x3vol(ix+6,irad+3))/2.
+ xcoord(6,ispec) = x3vol(ix+8,irad+2)
+ xcoord(4,ispec) = x3vol(ix+4,irad+4)
+ xcoord(7,ispec) = x3vol(ix+6,irad+4)
+ xcoord(3,ispec) = x3vol(ix+8,irad+4)
+
+ ycoord(1,ispec) = y3vol(ix+4,irad+2)
+ ycoord(5,ispec) = y3vol(ix+6,irad+1)
+ ycoord(2,ispec) = y3vol(ix+8,irad)
+ ycoord(8,ispec) = y3vol(ix+4,irad+3)
+ ycoord(9,ispec) = (y3vol(ix+6,irad+2) + y3vol(ix+6,irad+3))/2.
+ ycoord(6,ispec) = y3vol(ix+8,irad+2)
+ ycoord(4,ispec) = y3vol(ix+4,irad+4)
+ ycoord(7,ispec) = y3vol(ix+6,irad+4)
+ ycoord(3,ispec) = y3vol(ix+8,irad+4)
+
+ enddo
+
+! --- zone de raccord geometrique conforme inverse
+ irad=nspec_rad_ICB_to_doubling_OC-4
+ do ix=8,nspec_surf_whole_circle/factor_divide_mesh-8,16
+
+ ispec = ispec + 1
+
+ idoubling(ispec) = IREGION_OUTER_CORE
+
+ xcoord(1,ispec) = x3vol(ix+4,irad+2)
+ xcoord(5,ispec) = x3vol(ix+6,irad+2)
+ xcoord(2,ispec) = x3vol(ix+8,irad+2)
+ xcoord(8,ispec) = x3vol(ix+4,irad+3)
+ xcoord(9,ispec) = x3vol(ix+6,irad+3)
+ xcoord(6,ispec) = x3vol(ix+8,irad+3)
+ xcoord(4,ispec) = x3vol(ix+4,irad+4)
+ xcoord(7,ispec) = x3vol(ix+6,irad+4)
+ xcoord(3,ispec) = x3vol(ix+8,irad+4)
+
+ ycoord(1,ispec) = y3vol(ix+4,irad+2)
+ ycoord(5,ispec) = y3vol(ix+6,irad+2)
+ ycoord(2,ispec) = y3vol(ix+8,irad+2)
+ ycoord(8,ispec) = y3vol(ix+4,irad+3)
+ ycoord(9,ispec) = y3vol(ix+6,irad+3)
+ ycoord(6,ispec) = y3vol(ix+8,irad+3)
+ ycoord(4,ispec) = y3vol(ix+4,irad+4)
+ ycoord(7,ispec) = y3vol(ix+6,irad+4)
+ ycoord(3,ispec) = y3vol(ix+8,irad+4)
+
+ ispec = ispec + 1
+
+ idoubling(ispec) = IREGION_OUTER_CORE
+
+ xcoord(1,ispec) = x3vol(ix ,irad)
+ xcoord(5,ispec) = x3vol(ix+4,irad)
+ xcoord(2,ispec) = x3vol(ix+8,irad)
+ xcoord(8,ispec) = x3vol(ix+2,irad+1)
+ xcoord(9,ispec) = (x3vol(ix+4,irad+1) + x3vol(ix+6,irad+1))/2.
+ xcoord(6,ispec) = x3vol(ix+8,irad+1)
+ xcoord(4,ispec) = x3vol(ix+4,irad+2)
+ xcoord(7,ispec) = x3vol(ix+6,irad+2)
+ xcoord(3,ispec) = x3vol(ix+8,irad+2)
+
+ ycoord(1,ispec) = y3vol(ix ,irad)
+ ycoord(5,ispec) = y3vol(ix+4,irad)
+ ycoord(2,ispec) = y3vol(ix+8,irad)
+ ycoord(8,ispec) = y3vol(ix+2,irad+1)
+ ycoord(9,ispec) = (y3vol(ix+4,irad+1) + y3vol(ix+6,irad+1))/2.
+ ycoord(6,ispec) = y3vol(ix+8,irad+1)
+ ycoord(4,ispec) = y3vol(ix+4,irad+2)
+ ycoord(7,ispec) = y3vol(ix+6,irad+2)
+ ycoord(3,ispec) = y3vol(ix+8,irad+2)
+
+ ispec = ispec + 1
+
+ idoubling(ispec) = IREGION_OUTER_CORE
+
+ xcoord(1,ispec) = x3vol(ix,irad)
+ xcoord(5,ispec) = x3vol(ix+2,irad+1)
+ xcoord(2,ispec) = x3vol(ix+4,irad+2)
+ xcoord(8,ispec) = x3vol(ix,irad+2)
+ xcoord(9,ispec) = (x3vol(ix+2,irad+2) + x3vol(ix+2,irad+3))/2.
+ xcoord(6,ispec) = x3vol(ix+4,irad+3)
+ xcoord(4,ispec) = x3vol(ix,irad+4)
+ xcoord(7,ispec) = x3vol(ix+2,irad+4)
+ xcoord(3,ispec) = x3vol(ix+4,irad+4)
+
+ ycoord(1,ispec) = y3vol(ix,irad)
+ ycoord(5,ispec) = y3vol(ix+2,irad+1)
+ ycoord(2,ispec) = y3vol(ix+4,irad+2)
+ ycoord(8,ispec) = y3vol(ix,irad+2)
+ ycoord(9,ispec) = (y3vol(ix+2,irad+2) + y3vol(ix+2,irad+3))/2.
+ ycoord(6,ispec) = y3vol(ix+4,irad+3)
+ ycoord(4,ispec) = y3vol(ix,irad+4)
+ ycoord(7,ispec) = y3vol(ix+2,irad+4)
+ ycoord(3,ispec) = y3vol(ix+4,irad+4)
+
+ enddo
+
+! %%%% zone Cube -> ICB %%%%
+
+! generation maillage de la surface
+
+ do ix=0,nspec_surf_whole_circle
+
+ xicoord = dble(ix)/dble(nspec_surf_whole_circle)
+
+! parcours de l'angle dans le sens negatif with a different starting point
+ thetacoord = +pi/2 - 2*pi*xicoord
+
+! coordonnees cartesiennes correspondantes
+
+!---- ICB
+
+ radius_interf = RICB
+
+ x4surf(ix) = radius_interf * dcos(thetacoord)
+ y4surf(ix) = radius_interf * dsin(thetacoord)
+
+!---- Cube
+
+ if(thetacoord > pi/4) thetacoord = thetacoord - 2*pi
+
+! Xmax face (right)
+ if(thetacoord >= - pi/4 .and. thetacoord <= + pi/4) then
+ ratio_x = +1
+ ratio_y = (thetacoord - (-pi/4)) / (pi/2)
+
+! Ymin face (bottom)
+ else if(thetacoord >= - 3*pi/4 .and. thetacoord <= - pi/4) then
+ ratio_x = (thetacoord - (-3*pi/4)) / (pi/2)
+ ratio_y = 0
+
+! Xmin face (left)
+ else if(thetacoord >= - 5*pi/4 .and. thetacoord <= - 3*pi/4) then
+ ratio_x = 0
+ ratio_y = 1 - (thetacoord - (-5*pi/4)) / (pi/2)
+
+! Ymax face (top)
+ else
+ ratio_x = 1 - (thetacoord - (-7*pi/4)) / (pi/2)
+ ratio_y = +1
+
+ endif
+
+!---- volume
+
+! use a "flat" cubed sphere to create the central cube
+
+! map ratio to [-1,1] and then map to real radius
+! then add deformation
+ fact_x = 2.d0*ratio_x-1.d0
+ fact_y = 2.d0*ratio_y-1.d0
+
+ xi = PI_OVER_TWO*fact_x
+ eta = PI_OVER_TWO*fact_y
+
+ x4bot(ix) = radius_cube * fact_x * (1 + cos(eta)*CENTRAL_CUBE_INFLATE_FACTOR / PI)
+ y4bot(ix) = radius_cube * fact_y * (1 + cos(xi)*CENTRAL_CUBE_INFLATE_FACTOR / PI)
+
+!---- volume
+
+ do irad=0,nspec_rad_Cube_ICB
+ radcoord = dble(irad)/dble(nspec_rad_Cube_ICB)
+ x4vol(ix,irad) = x4surf(ix) * radcoord + x4bot(ix) * (one - radcoord)
+ y4vol(ix,irad) = y4surf(ix) * radcoord + y4bot(ix) * (one - radcoord)
+ enddo
+
+ enddo
+
+ do irad=0,nspec_rad_Cube_ICB-4,4
+ do ix=0,nspec_surf_whole_circle/factor_divide_mesh-8,8
+
+ ispec = ispec + 1
+
+ idoubling(ispec) = IREGION_INNER_CORE
+
+ xcoord(1,ispec) = x4vol(ix ,irad)
+ xcoord(5,ispec) = x4vol(ix+4,irad)
+ xcoord(2,ispec) = x4vol(ix+8,irad)
+ xcoord(8,ispec) = x4vol(ix ,irad+2)
+ xcoord(9,ispec) = x4vol(ix+4,irad+2)
+ xcoord(6,ispec) = x4vol(ix+8,irad+2)
+ xcoord(4,ispec) = x4vol(ix ,irad+4)
+ xcoord(7,ispec) = x4vol(ix+4,irad+4)
+ xcoord(3,ispec) = x4vol(ix+8,irad+4)
+
+ ycoord(1,ispec) = y4vol(ix ,irad)
+ ycoord(5,ispec) = y4vol(ix+4,irad)
+ ycoord(2,ispec) = y4vol(ix+8,irad)
+ ycoord(8,ispec) = y4vol(ix ,irad+2)
+ ycoord(9,ispec) = y4vol(ix+4,irad+2)
+ ycoord(6,ispec) = y4vol(ix+8,irad+2)
+ ycoord(4,ispec) = y4vol(ix ,irad+4)
+ ycoord(7,ispec) = y4vol(ix+4,irad+4)
+ ycoord(3,ispec) = y4vol(ix+8,irad+4)
+
+ enddo
+ enddo
+
+!----
+!---- generer l'interieur du cube
+!----
+
+ do ix=0,nspec_surf_whole_circle/16
+
+ xlincoord = dble(ix)/dble(nspec_surf_whole_circle/16)
+
+!---- volume
+
+ do irad=0,nspec_surf_whole_circle/16
+ radcoord = dble(irad)/dble(nspec_surf_whole_circle/16)
+
+! use a "flat" cubed sphere to create the central cube
+ ratio_x = xlincoord
+ ratio_y = radcoord
+
+! map ratio to [-1,1] and then map to real radius
+! then add deformation
+ fact_x = 2.d0*ratio_x-1.d0
+ fact_y = 2.d0*ratio_y-1.d0
+
+ xi = PI_OVER_TWO*fact_x
+ eta = PI_OVER_TWO*fact_y
+
+ x5vol(ix,irad) = radius_cube * fact_x * (1 + cos(eta)*CENTRAL_CUBE_INFLATE_FACTOR / PI)
+ y5vol(ix,irad) = radius_cube * fact_y * (1 + cos(xi)*CENTRAL_CUBE_INFLATE_FACTOR / PI)
+
+ enddo
+
+ enddo
+
+ if(generate_only_half_the_mesh) then
+ icentral_cube1 = nspec_surf_whole_circle/factor_divide_mesh/16
+ icentral_cube2 = nspec_surf_whole_circle/16-2
+ else
+ icentral_cube1 = 0
+ icentral_cube2 = nspec_surf_whole_circle/16-2
+ endif
+
+ do irad=0,nspec_surf_whole_circle/16-2,2
+ do ix=icentral_cube1,icentral_cube2,2
+
+ ispec = ispec + 1
+
+ idoubling(ispec) = IREGION_INNER_CORE
+
+ xcoord(1,ispec) = x5vol(ix ,irad)
+ xcoord(5,ispec) = x5vol(ix+1,irad)
+ xcoord(2,ispec) = x5vol(ix+2,irad)
+ xcoord(8,ispec) = x5vol(ix ,irad+1)
+ xcoord(9,ispec) = x5vol(ix+1,irad+1)
+ xcoord(6,ispec) = x5vol(ix+2,irad+1)
+ xcoord(4,ispec) = x5vol(ix ,irad+2)
+ xcoord(7,ispec) = x5vol(ix+1,irad+2)
+ xcoord(3,ispec) = x5vol(ix+2,irad+2)
+
+ ycoord(1,ispec) = y5vol(ix ,irad)
+ ycoord(5,ispec) = y5vol(ix+1,irad)
+ ycoord(2,ispec) = y5vol(ix+2,irad)
+ ycoord(8,ispec) = y5vol(ix ,irad+1)
+ ycoord(9,ispec) = y5vol(ix+1,irad+1)
+ ycoord(6,ispec) = y5vol(ix+2,irad+1)
+ ycoord(4,ispec) = y5vol(ix ,irad+2)
+ ycoord(7,ispec) = y5vol(ix+1,irad+2)
+ ycoord(3,ispec) = y5vol(ix+2,irad+2)
+
+ enddo
+ enddo
+
+! stocker le nombre total d'elements spectraux generes
+ nspec = ispec
+
+ print *
+ print *,'Total number of spectral elements = ',nspec
+ print *
+
+ if(nspec /= nspec_exact) stop 'incorrect number of spectral elements generated'
+
+!
+!---- generation de la numerotation pour SPECFEM90
+!
+
+ print *,'Generation de la numerotation pour Specfem...'
+
+! get coordinates of the grid points
+ xp(:) = 0.d0
+ yp(:) = 0.d0
+ do ispec=1,nspec
+ ieoff = ngnod*(ispec - 1)
+ ilocnum = 0
+
+ do ia = 1,ngnod
+ ilocnum = ilocnum + 1
+ xp(ilocnum + ieoff) = xcoord(ia,ispec)
+ yp(ilocnum + ieoff) = ycoord(ia,ispec)
+ enddo
+
+ enddo
+
+! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+! Establish initial pointers
+ do ie=1,nspec
+ ieoff = ngnod*(ie -1)
+ do ix=1,ngnod
+ loc (ix+ieoff) = ix+ieoff
+ enddo
+ enddo
+
+! set up local geometric tolerances
+
+ xtypdist=+HUGEVAL
+
+ do ispec=1,nspec
+
+ xminval=+HUGEVAL
+ yminval=+HUGEVAL
+ xmaxval=-HUGEVAL
+ ymaxval=-HUGEVAL
+ ieoff=ngnod*(ispec-1)
+ do ilocnum=1,ngnod
+ xmaxval=max(xp(ieoff+ilocnum),xmaxval)
+ xminval=min(xp(ieoff+ilocnum),xminval)
+ ymaxval=max(yp(ieoff+ilocnum),ymaxval)
+ yminval=min(yp(ieoff+ilocnum),yminval)
+ enddo
+
+! compute the minimum typical "size" of an element in the mesh
+ xtypdist = min(xtypdist,xmaxval-xminval)
+ xtypdist = min(xtypdist,ymaxval-yminval)
+
+ enddo
+
+! define a tolerance, small with respect to the minimum size
+ xtol = SMALLVALTOL * xtypdist
+
+ ifseg(:) = .false.
+ nseg = 1
+ ifseg(1) = .true.
+ ninseg(1) = npoin_max
+
+ do j=1,NDIM
+! Sort within each segment
+ ioff=1
+ do iseg=1,nseg
+ if (j==1) then
+ call rank(xp(ioff),ind,ninseg(iseg))
+ else
+ call rank(yp(ioff),ind,ninseg(iseg))
+ endif
+ call swap(xp(ioff),work,ind,ninseg(iseg))
+ call swap(yp(ioff),work,ind,ninseg(iseg))
+ call iswap(loc(ioff),work,ind,ninseg(iseg))
+ ioff=ioff+ninseg(iseg)
+ enddo
+! Check for jumps in current coordinate
+ if (j==1) then
+ do i=2,npoin_max
+ if (abs(xp(i)-xp(i-1)) > xtol) ifseg(i)=.true.
+ enddo
+ else
+ do i=2,npoin_max
+ if (abs(yp(i)-yp(i-1)) > xtol) ifseg(i)=.true.
+ enddo
+ endif
+! Count up number of different segments
+ nseg = 0
+ do i=1,npoin_max
+ if (ifseg(i)) then
+ nseg = nseg+1
+ ninseg(nseg) = 1
+ else
+ ninseg(nseg) = ninseg(nseg) + 1
+ endif
+ enddo
+ enddo
+!
+! Assign global node numbers (now sorted lexicographically!)
+!
+ ig = 0
+ do i=1,npoin_max
+ if (ifseg(i)) ig=ig+1
+ iglob(loc(i)) = ig
+ enddo
+
+ npoin = ig
+
+! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+! verifier la coherence du nombre de points generes
+ if(npoin > npoin_max) stop 'Too many points generated'
+
+! verification de la coherence de la numerotation generee
+! conversion from iglob() to ibool() is handled automatically by the "equivalence" statement
+ if(minval(ibool) /= 1 .or. maxval(ibool) /= npoin) stop 'Error when generating global numbering'
+
+ print *,'Total number of points of the global mesh: ',npoin
+ print *
+
+! generer les coordonnees des points du maillage global
+! in global numbering
+ print *,'Generating the coordinates of the points of the global mesh...'
+ do ispec=1,nspec
+ do ia = 1,ngnod
+ xp(ibool(ia,ispec)) = xcoord(ia,ispec)
+ yp(ibool(ia,ispec)) = ycoord(ia,ispec)
+ enddo
+ enddo
+ print *
+
+! ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+!
+!---- write the external mesh database for SPECFEM2D
+!
+
+ print *,'Writing the external mesh database for SPECFEM2D...'
+
+#ifdef USE_BINARY_FOR_EXTERNAL_MESH_DATABASE
+ print *,'using binary file format to store the mesh databases'
+#else
+ print *,'using ASCII file format to store the mesh databases'
+#endif
+
+! write the mesh points
+#ifdef USE_BINARY_FOR_EXTERNAL_MESH_DATABASE
+ open(unit=22,file='Nodes_AK135F_NO_MUD',form='unformatted',status='unknown',action='write')
+ write(22) npoin
+#else
+ open(unit=22,file='Nodes_AK135F_NO_MUD',form='formatted',status='unknown',action='write')
+ write(22,*) npoin
+#endif
+ do ipoin = 1,npoin
+#ifdef USE_BINARY_FOR_EXTERNAL_MESH_DATABASE
+ write(22) xp(ipoin),yp(ipoin)
+#else
+ write(22,*) xp(ipoin),yp(ipoin)
+#endif
+ enddo
+ close(22)
+
+! write the mesh elements
+#ifdef USE_BINARY_FOR_EXTERNAL_MESH_DATABASE
+ open(unit=22,file='Mesh_AK135F_NO_MUD',form='unformatted',status='unknown',action='write')
+ write(22) nspec
+#else
+ open(unit=22,file='Mesh_AK135F_NO_MUD',form='formatted',status='unknown',action='write')
+ write(22,*) nspec
+#endif
+ do ispec = 1,nspec
+#ifdef USE_BINARY_FOR_EXTERNAL_MESH_DATABASE
+ write(22) (ibool(ia,ispec), ia=1,ngnod)
+#else
+ write(22,*) (ibool(ia,ispec), ia=1,ngnod)
+#endif
+ enddo
+ close(22)
+
+! write the material properties
+#ifdef USE_BINARY_FOR_EXTERNAL_MESH_DATABASE
+ open(unit=22,file='Material_AK135F_NO_MUD',form='unformatted',status='unknown',action='write')
+#else
+ open(unit=22,file='Material_AK135F_NO_MUD',form='formatted',status='unknown',action='write')
+#endif
+ do ispec = 1,nspec
+#ifdef USE_BINARY_FOR_EXTERNAL_MESH_DATABASE
+ write(22) idoubling(ispec)
+#else
+ write(22,*) idoubling(ispec)
+#endif
+ enddo
+ close(22)
+
+! write empty file for the absorbing boundary
+#ifdef USE_BINARY_FOR_EXTERNAL_MESH_DATABASE
+ open(unit=22,file='Surf_abs_AK135F_NO_MUD',form='unformatted',status='unknown',action='write')
+ write(22) 0
+#else
+ open(unit=22,file='Surf_abs_AK135F_NO_MUD',form='formatted',status='unknown',action='write')
+ write(22,*) 0
+#endif
+ close(22)
+
+! write empty file for the acoustic free surface boundary (since there is no acoustic free surface in the 1D Earth without oceans)
+#ifdef USE_BINARY_FOR_EXTERNAL_MESH_DATABASE
+ open(unit=22,file='Surf_free_AK135F_NO_MUD',form='unformatted',status='unknown',action='write')
+ write(22) 0
+#else
+ open(unit=22,file='Surf_free_AK135F_NO_MUD',form='formatted',status='unknown',action='write')
+ write(22,*) 0
+#endif
+ close(22)
+
+ print *,'Done writing the external mesh database for SPECFEM2D'
+ print *
+
+! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+! list all the elements that are in contact with the symmetry axis, and by what edge
+
+ ispec_count = 0
+
+! count the number of elements that are in contact with the symmetry axis
+ do ispec=1,nspec
+ i = 0
+
+ if(xcoord(1,ispec) < 0.001d0) i = i + 1
+ if(xcoord(2,ispec) < 0.001d0) i = i + 1
+ if(xcoord(3,ispec) < 0.001d0) i = i + 1
+ if(xcoord(4,ispec) < 0.001d0) i = i + 1
+
+ if(i > 0) then
+! we know in advance from the way the mesh is designed that single points can never be detected, only a full edge,
+! and a single edge can of course be detected otherwise that element would have a 180-degree angle
+ if(i /= 2) stop 'error: element in contact with the symmetry axis by incorrect number of points'
+ ispec_count = ispec_count + 1
+ endif
+ enddo
+
+ print *,'number of elements in contact with the symmetry axis = ',ispec_count
+
+! then save them to a file and also save the (only) edge that is in contact
+#ifdef USE_BINARY_FOR_EXTERNAL_MESH_DATABASE
+ open(unit=22,file='Symmetry_axis_elements_and_edges_AK135F_NO_MUD',form='unformatted',status='unknown',action='write')
+ write(22) ispec_count
+#else
+ open(unit=22,file='Symmetry_axis_elements_and_edges_AK135F_NO_MUD',form='formatted',status='unknown',action='write')
+ write(22,*) ispec_count
+#endif
+
+ do ispec=1,nspec
+#ifdef USE_BINARY_FOR_EXTERNAL_MESH_DATABASE
+ if(xcoord(1,ispec) < 0.001d0 .and. xcoord(2,ispec) < 0.001d0) write(22) ispec,IBOTTOM
+ if(xcoord(2,ispec) < 0.001d0 .and. xcoord(3,ispec) < 0.001d0) write(22) ispec,IRIGHT
+ if(xcoord(3,ispec) < 0.001d0 .and. xcoord(4,ispec) < 0.001d0) write(22) ispec,ITOP
+ if(xcoord(4,ispec) < 0.001d0 .and. xcoord(1,ispec) < 0.001d0) write(22) ispec,ILEFT
+#else
+ if(xcoord(1,ispec) < 0.001d0 .and. xcoord(2,ispec) < 0.001d0) write(22,*) ispec,IBOTTOM
+ if(xcoord(2,ispec) < 0.001d0 .and. xcoord(3,ispec) < 0.001d0) write(22,*) ispec,IRIGHT
+ if(xcoord(3,ispec) < 0.001d0 .and. xcoord(4,ispec) < 0.001d0) write(22,*) ispec,ITOP
+ if(xcoord(4,ispec) < 0.001d0 .and. xcoord(1,ispec) < 0.001d0) write(22,*) ispec,ILEFT
+#endif
+ enddo
+
+ close(22)
+
+! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+! ***
+! *** generer un fichier 'GNUPLOT' pour le controle de la grille ***
+! ***
+
+ if(output_gnuplot_grid) then
+
+ print *
+ print *,' Ecriture de la grille format GNUPLOT...'
+
+ open(unit=20,file='gridfile.gnu',status='unknown')
+
+ do ispec=1,nspec
+
+! draw the four edges of each element (using straight lines to simplify)
+ ia1 = 1
+ ia2 = 2
+ write(20,15) sngl(xcoord(ia1,ispec)),sngl(ycoord(ia1,ispec))
+ write(20,15) sngl(xcoord(ia2,ispec)),sngl(ycoord(ia2,ispec))
+ write(20,10)
+
+ ia1 = 2
+ ia2 = 3
+ write(20,15) sngl(xcoord(ia1,ispec)),sngl(ycoord(ia1,ispec))
+ write(20,15) sngl(xcoord(ia2,ispec)),sngl(ycoord(ia2,ispec))
+ write(20,10)
+
+ ia1 = 3
+ ia2 = 4
+ write(20,15) sngl(xcoord(ia1,ispec)),sngl(ycoord(ia1,ispec))
+ write(20,15) sngl(xcoord(ia2,ispec)),sngl(ycoord(ia2,ispec))
+ write(20,10)
+
+ ia1 = 4
+ ia2 = 1
+ write(20,15) sngl(xcoord(ia1,ispec)),sngl(ycoord(ia1,ispec))
+ write(20,15) sngl(xcoord(ia2,ispec)),sngl(ycoord(ia2,ispec))
+ write(20,10)
+
+ enddo
+
+ close(20)
+
+! cree le script de dessin pour gnuplot
+ open(unit=20,file='plotgrid.gnu',status='unknown')
+ write(20,*) '#set term postscript landscape monochrome solid "Helvetica" 22'
+ write(20,*) '#set output "grille.ps"'
+ write(20,*) 'set term x11'
+ write(20,*) 'set size ratio -1'
+ write(20,*) 'plot "gridfile.gnu" title "Macroblocs mesh" w l'
+ write(20,*) 'pause -1 "Hit any key..."'
+ close(20)
+
+ print *,' Fin ecriture de la grille format GNUPLOT'
+ print *
+
+ 10 format('')
+ 15 format(e12.5,1x,e12.5)
+
+ endif
+
+ end
+
+!-----------------------------------------------------------------------
+
+ subroutine rank(A,IND,N)
+!
+! Use Heap Sort (p 233 Numerical Recipes)
+!
+ implicit none
+
+ integer N
+ double precision A(N)
+ integer IND(N)
+
+ integer i,j,l,ir,indx
+ double precision q
+
+ do J=1,N
+ IND(j)=j
+ enddo
+
+ if(n == 1) return
+ L=n/2+1
+ ir=n
+ 100 continue
+ IF(l > 1) THEN
+ l=l-1
+ indx=ind(l)
+ q=a(indx)
+ ELSE
+ indx=ind(ir)
+ q=a(indx)
+ ind(ir)=ind(1)
+ ir=ir-1
+ if(ir == 1) then
+ ind(1)=indx
+ return
+ endif
+ ENDIF
+ i=l
+ j=l+l
+ 200 continue
+ IF(J <= IR) THEN
+ IF(J < IR) THEN
+ IF(A(IND(j)) < A(IND(j+1))) j=j+1
+ ENDIF
+ IF(q < A(IND(j))) THEN
+ IND(I)=IND(J)
+ I=J
+ J=J+J
+ ELSE
+ J=IR+1
+ ENDIF
+ GOTO 200
+ ENDIF
+ IND(I)=INDX
+ GOTO 100
+
+ end subroutine rank
+
+!-----------------------------------------------------------------------
+
+ subroutine swap(a,w,ind,n)
+!
+! Use IND to sort array A (p 233 Numerical Recipes)
+!
+ implicit none
+
+ integer n
+ double precision A(N),W(N)
+ integer IND(N)
+
+ integer j
+
+ W(:) = A(:)
+
+ do J=1,N
+ A(j) = W(ind(j))
+ enddo
+
+ end subroutine swap
+
+!-----------------------------------------------------------------------
+
+ subroutine iswap(a,w,ind,n)
+!
+! Use IND to sort array A
+!
+ implicit none
+
+ integer n
+ integer A(N),W(N),IND(N)
+
+ integer j
+
+ W(:) = A(:)
+
+ do J=1,N
+ A(j) = W(ind(j))
+ enddo
+
+ end subroutine iswap
+
Added: seismo/2D/SPECFEM2D/trunk/EXAMPLES/global_Earth_ak135f/make_all_check.csh
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/global_Earth_ak135f/make_all_check.csh (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/global_Earth_ak135f/make_all_check.csh 2012-12-15 01:15:45 UTC (rev 21168)
@@ -0,0 +1,8 @@
+#!/bin/csh
+
+rm -f xcreate_mesh_files
+
+ifort -o xcreate_mesh_files -O0 -vec-report0 -std03 -implicitnone -warn truncated_source -warn argument_checking -warn unused -warn declarations -warn alignments -warn ignore_loc -warn usage -check all -align sequence -assume byterecl -ftz -traceback -ftrapuv create_mesh_AK135F_2D_with_central_cube_no_PML.F90
+
+./xcreate_mesh_files
+
Property changes on: seismo/2D/SPECFEM2D/trunk/EXAMPLES/global_Earth_ak135f/make_all_check.csh
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/2D/SPECFEM2D/trunk/EXAMPLES/global_Earth_ak135f/make_all_fast.csh
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/global_Earth_ak135f/make_all_fast.csh (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/global_Earth_ak135f/make_all_fast.csh 2012-12-15 01:15:45 UTC (rev 21168)
@@ -0,0 +1,8 @@
+#!/bin/csh
+
+rm -f xcreate_mesh_files
+
+ifort -o xcreate_mesh_files -O3 -xHost -vec-report0 -std03 -implicitnone -warn truncated_source -warn argument_checking -warn unused -warn declarations -warn alignments -warn ignore_loc -warn usage -check nobounds -mcmodel=medium -shared-intel create_mesh_AK135F_2D_with_central_cube_no_PML.F90
+
+./xcreate_mesh_files
+
Property changes on: seismo/2D/SPECFEM2D/trunk/EXAMPLES/global_Earth_ak135f/make_all_fast.csh
___________________________________________________________________
Name: svn:executable
+ *
More information about the CIG-COMMITS
mailing list