[cig-commits] r12642 - in seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta: setup src
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Thu Aug 14 20:35:20 PDT 2008
Author: dkomati1
Date: 2008-08-14 20:35:19 -0700 (Thu, 14 Aug 2008)
New Revision: 12642
Added:
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/combine_AVS_DX.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/create_name_database.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/write_AVS_DX_global_chunks_data.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/write_AVS_DX_global_data.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/write_AVS_DX_global_faces_data.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/write_AVS_DX_surface_data.f90
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/setup/constants.h
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/create_header_file.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/create_regions_mesh.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/main_program.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/meshfem3D.F90
Log:
added codes to visualise AVS/OpenDX meshes;
fixed a small bug when calling attenuation;
added code to debug second sorting of ibool;
added MESHER_ONLY option
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/setup/constants.h
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/setup/constants.h 2008-08-15 00:52:09 UTC (rev 12641)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/setup/constants.h 2008-08-15 03:35:19 UTC (rev 12642)
@@ -171,6 +171,9 @@
logical, parameter :: ACTUALLY_COUPLE_FLUID_CMB = .true.
logical, parameter :: ACTUALLY_COUPLE_FLUID_ICB = .true.
+! flag to only create the mesh but not start the solver (for instance to check the mesh obtained)
+ logical, parameter :: MESHER_ONLY = .true.
+
!------------------------------------------------------
!----------- do not modify anything below -------------
!------------------------------------------------------
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/combine_AVS_DX.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/combine_AVS_DX.F90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/combine_AVS_DX.F90 2008-08-15 03:35:19 UTC (rev 12642)
@@ -0,0 +1,1207 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 4 . 1
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory, California Institute of Technology, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) California Institute of Technology and University of Pau / CNRS / INRIA
+! August 2008
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+! combine AVS or DX global data files to check the mesh
+! this is done in postprocessing after running the mesh generator
+
+ program combine_AVS_DX
+
+ implicit none
+
+ include "constants.h"
+
+! threshold for number of points per wavelength displayed
+! otherwise the scale is too large and we cannot see the small values
+! all values above this threshold are truncated
+ double precision, parameter :: THRESHOLD_GRIDPOINTS = 12.
+
+! non-linear scaling factor for elevation if topography for Earth model
+ double precision, parameter :: SCALE_NON_LINEAR = 0.3
+
+! maximum polynomial degree for which we can compute the stability condition
+ integer, parameter :: NGLL_MAX_STABILITY = 15
+
+ integer iproc,nspec,npoin
+ integer ispec
+ integer iglob1,iglob2,iglob3,iglob4
+ integer ipoin,numpoin,numpoin2,iglobpointoffset,ntotpoin,ntotspec
+ integer numelem,numelem2,iglobelemoffset,idoubling,maxdoubling
+ integer iformat,ivalue,icolor,itarget_doubling
+ integer imaterial,imatprop,ispec_scale_AVS_DX
+ integer nrec,ir,iregion_code
+ integer ntotpoinAVS_DX,ntotspecAVS_DX
+
+ real(kind=CUSTOM_REAL) vmin,vmax,deltavp,deltavs
+ double precision xval,yval,zval
+ double precision val_color,rnorm_factor
+
+ logical threshold_used
+ logical USE_OPENDX
+
+! for source location
+ integer yr,jda,ho,mi
+ double precision x_target_source,y_target_source,z_target_source
+ double precision r_target_source
+ double precision x_source_trgl1,y_source_trgl1,z_source_trgl1
+ double precision x_source_trgl2,y_source_trgl2,z_source_trgl2
+ double precision x_source_trgl3,y_source_trgl3,z_source_trgl3
+ double precision theta,phi,delta_trgl
+ double precision sec,t_cmt,hdur
+ double precision lat,long,depth
+ double precision moment_tensor(6)
+
+! for receiver location
+ integer irec,ios
+ double precision r_target
+ double precision, allocatable, dimension(:) :: stlat,stlon,stele,stbur
+ character(len=MAX_LENGTH_STATION_NAME), allocatable, dimension(:) :: station_name
+ character(len=MAX_LENGTH_NETWORK_NAME), allocatable, dimension(:) :: network_name
+ character(len=150) dummystring
+
+ double precision, allocatable, dimension(:) :: x_target,y_target,z_target
+
+! for the reference ellipsoid
+ double precision reference,radius_dummy,theta_s,phi_s
+
+! processor identification
+ character(len=150) prname
+
+! small offset for source and receiver line in AVS_DX
+! (small compared to normalized radius of the Earth)
+
+! for full Earth
+ double precision, parameter :: small_offset_source_earth = 0.025d0
+ double precision, parameter :: small_offset_receiver_earth = 0.0125d0
+
+! for oceans only
+ logical OCEANS_ONLY
+ integer ioceans
+ integer above_zero,below_zero
+
+! for stability condition
+ double precision, dimension (:), allocatable :: stability_value,gridpoints_per_wavelength,elevation_sphere
+ double precision, dimension (:), allocatable :: dvp,dvs
+ double precision, dimension (:), allocatable :: xcoord,ycoord,zcoord,vmincoord,vmaxcoord
+ double precision stability_value_min,stability_value_max
+ double precision gridpoints_per_wavelength_min,gridpoints_per_wavelength_max
+ integer iloop_corners,istab,jstab
+ integer ipointnumber1_horiz,ipointnumber2_horiz
+ integer ipointnumber1_vert,ipointnumber2_vert
+ double precision distance_horiz,distance_vert
+ double precision stabmax,gridmin,scale_factor
+ integer NGLL_current_horiz,NGLL_current_vert
+ double precision :: percent_GLL(NGLL_MAX_STABILITY)
+
+! for chunk numbering
+ integer :: iproc_eta,iproc_xi,iprocnum,ichunk
+
+ integer, dimension(:), allocatable :: ichunk_slice
+
+! parameters read from parameter file
+ integer MIN_ATTENUATION_PERIOD,MAX_ATTENUATION_PERIOD,NER_CRUST, &
+ NER_80_MOHO,NER_220_80,NER_400_220,NER_600_400,NER_670_600,NER_771_670, &
+ NER_TOPDDOUBLEPRIME_771,NER_CMB_TOPDDOUBLEPRIME,NER_OUTER_CORE, &
+ NER_TOP_CENTRAL_CUBE_ICB,NEX_XI,NEX_ETA,RMOHO_FICTITIOUS_IN_MESHER, &
+ NPROC_XI,NPROC_ETA,NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
+ NTSTEP_BETWEEN_READ_ADJSRC,NSTEP,NTSTEP_BETWEEN_FRAMES, &
+ NTSTEP_BETWEEN_OUTPUT_INFO,NUMBER_OF_RUNS,NUMBER_OF_THIS_RUN,NCHUNKS,SIMULATION_TYPE, &
+ REFERENCE_1D_MODEL,THREE_D_MODEL,MOVIE_VOLUME_TYPE,MOVIE_START,MOVIE_STOP, &
+ ifirst_layer_aniso,ilast_layer_aniso
+
+ double precision DT,ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES,CENTER_LONGITUDE_IN_DEGREES, &
+ CENTER_LATITUDE_IN_DEGREES,GAMMA_ROTATION_AZIMUTH,ROCEAN,RMIDDLE_CRUST, &
+ RMOHO,R80,R120,R220,R400,R600,R670,R771,RTOPDDOUBLEPRIME,RCMB,RICB, &
+ R_CENTRAL_CUBE,RHO_TOP_OC,RHO_BOTTOM_OC,RHO_OCEANS,HDUR_MOVIE, &
+ MOVIE_TOP,MOVIE_BOTTOM,MOVIE_WEST,MOVIE_EAST,MOVIE_NORTH,MOVIE_SOUTH
+
+ logical TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE, &
+ CRUSTAL,ELLIPTICITY,GRAVITY,ONE_CRUST,ROTATION,ISOTROPIC_3D_MANTLE, &
+ TOPOGRAPHY,OCEANS,MOVIE_SURFACE,MOVIE_VOLUME,MOVIE_COARSE,ATTENUATION_3D, &
+ RECEIVERS_CAN_BE_BURIED,PRINT_SOURCE_TIME_FUNCTION, &
+ SAVE_MESH_FILES,ATTENUATION, &
+ ABSORBING_CONDITIONS,INCLUDE_CENTRAL_CUBE,INFLATE_CENTRAL_CUBE,SAVE_FORWARD,CASE_3D, &
+ OUTPUT_SEISMOS_ASCII_TEXT,OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY, &
+ ROTATE_SEISMOGRAMS_RT,HONOR_1D_SPHERICAL_MOHO,WRITE_SEISMOGRAMS_BY_MASTER,&
+ SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE
+
+ character(len=150) OUTPUT_FILES,MODEL
+
+! parameters deduced from parameters read from file
+ integer NPROC,NPROCTOT,NEX_PER_PROC_XI,NEX_PER_PROC_ETA,ratio_divide_central_cube
+
+! for all the regions
+ integer, dimension(MAX_NUM_REGIONS) :: NSPEC2D_XI,NSPEC2D_ETA, &
+ NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX, &
+ NSPEC2D_BOTTOM,NSPEC2D_TOP, &
+ NSPEC1D_RADIAL,NGLOB1D_RADIAL, &
+ NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX, &
+ nglob
+
+ integer region_min,region_max
+
+ double precision small_offset_source,small_offset_receiver
+
+ integer proc_p1,proc_p2
+
+! computed in read_compute_parameters
+ integer, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: ner,ratio_sampling_array
+ integer, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: doubling_index
+ double precision, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: r_bottom,r_top
+ logical, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: this_layer_has_a_doubling
+ double precision, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: rmins,rmaxs
+
+
+ logical :: CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA
+ integer, dimension(NB_SQUARE_CORNERS,NB_CUT_CASE) :: DIFF_NSPEC1D_RADIAL
+ integer, dimension(NB_SQUARE_EDGES_ONEDIR,NB_CUT_CASE) :: DIFF_NSPEC2D_XI,DIFF_NSPEC2D_ETA
+
+! ************** PROGRAM STARTS HERE **************
+
+ print *
+ print *,'Recombining all AVS or DX files for slices'
+ print *
+
+ print *
+ print *,'reading parameter file'
+ print *
+
+! read the parameter file and compute additional parameters
+ call read_compute_parameters(MIN_ATTENUATION_PERIOD,MAX_ATTENUATION_PERIOD,NER_CRUST, &
+ NER_80_MOHO,NER_220_80,NER_400_220,NER_600_400,NER_670_600,NER_771_670, &
+ NER_TOPDDOUBLEPRIME_771,NER_CMB_TOPDDOUBLEPRIME,NER_OUTER_CORE, &
+ NER_TOP_CENTRAL_CUBE_ICB,NEX_XI,NEX_ETA,RMOHO_FICTITIOUS_IN_MESHER, &
+ NPROC_XI,NPROC_ETA,NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
+ NTSTEP_BETWEEN_READ_ADJSRC,NSTEP,NTSTEP_BETWEEN_FRAMES, &
+ NTSTEP_BETWEEN_OUTPUT_INFO,NUMBER_OF_RUNS,NUMBER_OF_THIS_RUN,NCHUNKS,DT, &
+ ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES,CENTER_LONGITUDE_IN_DEGREES, &
+ CENTER_LATITUDE_IN_DEGREES,GAMMA_ROTATION_AZIMUTH,ROCEAN,RMIDDLE_CRUST, &
+ RMOHO,R80,R120,R220,R400,R600,R670,R771,RTOPDDOUBLEPRIME,RCMB,RICB, &
+ R_CENTRAL_CUBE,RHO_TOP_OC,RHO_BOTTOM_OC,RHO_OCEANS,HDUR_MOVIE,MOVIE_VOLUME_TYPE, &
+ MOVIE_TOP,MOVIE_BOTTOM,MOVIE_WEST,MOVIE_EAST,MOVIE_NORTH,MOVIE_SOUTH,MOVIE_START,MOVIE_STOP, &
+ TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE, &
+ ANISOTROPIC_INNER_CORE,CRUSTAL,ELLIPTICITY,GRAVITY,ONE_CRUST, &
+ ROTATION,ISOTROPIC_3D_MANTLE,TOPOGRAPHY,OCEANS,MOVIE_SURFACE, &
+ MOVIE_VOLUME,MOVIE_COARSE,ATTENUATION_3D,RECEIVERS_CAN_BE_BURIED, &
+ PRINT_SOURCE_TIME_FUNCTION,SAVE_MESH_FILES, &
+ ATTENUATION,REFERENCE_1D_MODEL,THREE_D_MODEL,ABSORBING_CONDITIONS, &
+ INCLUDE_CENTRAL_CUBE,INFLATE_CENTRAL_CUBE,MODEL,SIMULATION_TYPE,SAVE_FORWARD, &
+ NPROC,NPROCTOT,NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
+ NSPEC, &
+ NSPEC2D_XI, &
+ NSPEC2D_ETA, &
+ NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
+ NSPEC1D_RADIAL,NGLOB1D_RADIAL, &
+ NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX,NGLOB, &
+ ratio_sampling_array, ner, doubling_index,r_bottom,r_top,this_layer_has_a_doubling,rmins,rmaxs,CASE_3D, &
+ OUTPUT_SEISMOS_ASCII_TEXT,OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY, &
+ ROTATE_SEISMOGRAMS_RT,ratio_divide_central_cube,HONOR_1D_SPHERICAL_MOHO,CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA,&
+ DIFF_NSPEC1D_RADIAL,DIFF_NSPEC2D_XI,DIFF_NSPEC2D_ETA,&
+ WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE, &
+ USE_BINARY_FOR_LARGE_FILE,ifirst_layer_aniso,ilast_layer_aniso,.false.)
+
+ if(.not. SAVE_MESH_FILES) stop 'AVS or DX files were not saved by the mesher'
+
+! get the base pathname for output files
+ call get_value_string(OUTPUT_FILES, 'OUTPUT_FILES', 'OUTPUT_FILES')
+
+ print *,'1 = create files in OpenDX format'
+ print *,'2 = create files in AVS UCD format'
+ print *,'any other value = exit'
+ print *
+ print *,'enter value:'
+ read(5,*) iformat
+ if(iformat<1 .or. iformat>2) stop 'exiting...'
+ if(iformat == 1) then
+ USE_OPENDX = .true.
+ else
+ USE_OPENDX = .false.
+ endif
+
+ print *
+ print *,'1 = edges of all the slices only'
+ print *,'2 = edges of the chunks only'
+ print *,'3 = surface of the model only'
+ print *,'any other value = exit'
+ print *
+ print *,'enter value:'
+ read(5,*) ivalue
+ if(ivalue<1 .or. ivalue>3) stop 'exiting...'
+
+! warning if surface elevation
+ if(ivalue == 3) then
+ print *,'******************************************'
+ print *,'*** option 7 to color using topography ***'
+ print *,'******************************************'
+ endif
+
+ print *
+ print *,'1 = color by doubling flag'
+ print *,'2 = by slice number'
+ print *,'3 = by stability value'
+ print *,'4 = by gridpoints per wavelength'
+ print *,'5 = dvp/vp'
+ print *,'6 = dvs/vs'
+ print *,'7 = elevation of Earth model'
+ print *,'8 = by region number'
+ print *,'9 = focus on one doubling flag only'
+ print *,'any other value=exit'
+ print *
+ print *,'enter value:'
+ read(5,*) icolor
+ if(icolor<1 .or. icolor >9) stop 'exiting...'
+ if((icolor == 3 .or. icolor == 4) .and. ivalue /= 2) &
+ stop 'need chunks only to represent stability or gridpoints per wavelength'
+
+ if(icolor == 9) then
+ print *
+ print *,'enter value of target doubling flag:'
+ read(5,*) itarget_doubling
+ endif
+
+! for oceans only
+ OCEANS_ONLY = .false.
+ if(ivalue == 3 .and. icolor == 7) then
+ print *
+ print *,'1 = represent full topography (topo + oceans)'
+ print *,'2 = represent oceans only'
+ print *
+ read(5,*) ioceans
+ if(ioceans == 1) then
+ OCEANS_ONLY = .false.
+ else if(ioceans == 2) then
+ OCEANS_ONLY = .true.
+ else
+ stop 'incorrect option for the oceans'
+ endif
+ endif
+
+ print *
+ print *,'1 = material property by doubling flag'
+ print *,'2 = by slice number'
+ print *,'3 = by region number'
+ print *,'4 = by chunk number'
+ print *,'any other value=exit'
+ print *
+ print *,'enter value:'
+ read(5,*) imaterial
+ if(imaterial < 1 .or. imaterial > 4) stop 'exiting...'
+
+! user can specify a range of processors here
+#ifdef USE_MPI
+ print *
+ print *,'enter first proc (proc numbers start at 0) = '
+ read(5,*) proc_p1
+ if(proc_p1 < 0) proc_p1 = 0
+ if(proc_p1 > NPROCTOT-1) proc_p1 = NPROCTOT-1
+
+ print *,'enter last proc (enter -1 for all procs) = '
+ read(5,*) proc_p2
+ if(proc_p2 == -1) proc_p2 = NPROCTOT-1
+ if(proc_p2 < 0) proc_p2 = 0
+
+#else
+ proc_p1 = 0
+ proc_p2 = 0
+#endif
+
+ if(proc_p2 > NPROCTOT-1) proc_p2 = NPROCTOT-1
+
+! set interval to maximum if user input is not correct
+ if(proc_p1 <= 0) proc_p1 = 0
+ if(proc_p2 < 0) proc_p2 = NPROCTOT - 1
+
+ print *
+ print *,'There are ',NPROCTOT,' slices numbered from 0 to ',NPROCTOT-1
+ print *
+
+! loop on all the chunks to create global slice addressing for solver
+ write(*,*) 'creating global slice addressing'
+ write(*,*)
+
+ allocate(ichunk_slice(0:NPROCTOT-1))
+ ichunk_slice(:) = 0
+
+ do ichunk = 1,NCHUNKS
+ do iproc_eta = 0,NPROC_ETA-1
+ do iproc_xi = 0,NPROC_XI-1
+ iprocnum = (ichunk-1)*NPROC + iproc_eta * NPROC_XI + iproc_xi
+ ichunk_slice(iprocnum) = ichunk
+ enddo
+ enddo
+ enddo
+
+! define percentage of smallest distance between GLL points for NGLL points
+! percentages were computed by calling the GLL points routine for each degree
+ percent_GLL(2) = 100.d0
+ percent_GLL(3) = 50.d0
+ percent_GLL(4) = 27.639320225002102d0
+ percent_GLL(5) = 17.267316464601141d0
+ percent_GLL(6) = 11.747233803526763d0
+ percent_GLL(7) = 8.4888051860716516d0
+ percent_GLL(8) = 6.4129925745196719d0
+ percent_GLL(9) = 5.0121002294269914d0
+ percent_GLL(10) = 4.0233045916770571d0
+ percent_GLL(11) = 3.2999284795970416d0
+ percent_GLL(12) = 2.7550363888558858d0
+ percent_GLL(13) = 2.3345076678918053d0
+ percent_GLL(14) = 2.0032477366369594d0
+ percent_GLL(15) = 1.7377036748080721d0
+
+! convert to real percentage
+ percent_GLL(:) = percent_GLL(:) / 100.d0
+
+! clear flag to detect if threshold used
+ threshold_used = .false.
+
+! set length of segments for source and receiver representation
+ small_offset_source = small_offset_source_earth
+ small_offset_receiver = small_offset_receiver_earth
+
+! set total number of points and elements to zero
+ ntotpoin = 0
+ ntotspec = 0
+
+ region_min = 1
+ region_max = MAX_NUM_REGIONS
+
+! if representing surface elements, only one region
+ if(ivalue == 3) then
+ region_min = IREGION_CRUST_MANTLE
+ region_max = IREGION_CRUST_MANTLE
+ endif
+
+ do iregion_code = region_min,region_max
+
+! loop on the selected range of processors
+ do iproc = proc_p1,proc_p2
+
+ print *,'Reading slice ',iproc,' in region ',iregion_code
+
+! create the name for the database of the current slide
+ call create_name_database(prname,iproc,iregion_code)
+
+ if(ivalue == 1) then
+ open(unit=10,file='OUTPUT_FILES/'//prname(1:len_trim(prname))//'AVS_DXpointsfaces.txt',status='old',action='read')
+ else if(ivalue == 2) then
+ open(unit=10,file='OUTPUT_FILES/'//prname(1:len_trim(prname))//'AVS_DXpointschunks.txt',status='old',action='read')
+ else if(ivalue == 3) then
+ open(unit=10,file='OUTPUT_FILES/'//prname(1:len_trim(prname))//'AVS_DXpointssurface.txt',status='old',action='read')
+ endif
+
+ read(10,*) npoin
+ print *,'There are ',npoin,' global AVS or DX points in the slice'
+ ntotpoin = ntotpoin + npoin
+ close(10)
+
+ if(ivalue == 1) then
+ open(unit=10,file='OUTPUT_FILES/'//prname(1:len_trim(prname))//'AVS_DXelementsfaces.txt',status='old',action='read')
+ else if(ivalue == 2) then
+ open(unit=10,file='OUTPUT_FILES/'//prname(1:len_trim(prname))//'AVS_DXelementschunks.txt',status='old',action='read')
+ else if(ivalue == 3) then
+ open(unit=10,file='OUTPUT_FILES/'//prname(1:len_trim(prname))//'AVS_DXelementssurface.txt',status='old',action='read')
+ endif
+
+ read(10,*) nspec
+ print *,'There are ',nspec,' AVS or DX elements in the slice'
+ ntotspec = ntotspec + nspec
+ close(10)
+
+ enddo
+ enddo
+
+ print *
+ print *,'There is a total of ',ntotspec,' elements in all the slices'
+ print *,'There is a total of ',ntotpoin,' points in all the slices'
+ print *
+
+ ntotpoinAVS_DX = ntotpoin
+ ntotspecAVS_DX = ntotspec
+
+! write AVS or DX header with element data
+ if(USE_OPENDX) then
+ open(unit=11,file=trim(OUTPUT_FILES)//'/DX_fullmesh.dx',status='unknown')
+ write(11,*) 'object 1 class array type float rank 1 shape 3 items ',ntotpoinAVS_DX,' data follows'
+ else
+ open(unit=11,file=trim(OUTPUT_FILES)//'/AVS_fullmesh.inp',status='unknown')
+ write(11,*) ntotpoinAVS_DX,' ',ntotspecAVS_DX,' 0 1 0'
+ endif
+
+! allocate array for stability condition
+ allocate(stability_value(ntotspecAVS_DX))
+ allocate(gridpoints_per_wavelength(ntotspecAVS_DX))
+ allocate(elevation_sphere(ntotspecAVS_DX))
+ allocate(dvp(ntotspecAVS_DX))
+ allocate(dvs(ntotspecAVS_DX))
+ allocate(xcoord(ntotpoinAVS_DX))
+ allocate(ycoord(ntotpoinAVS_DX))
+ allocate(zcoord(ntotpoinAVS_DX))
+ allocate(vmincoord(ntotpoinAVS_DX))
+ allocate(vmaxcoord(ntotpoinAVS_DX))
+
+! ************* generate points ******************
+
+! set global point offset to zero
+ iglobpointoffset = 0
+
+ do iregion_code = region_min,region_max
+
+! loop on the selected range of processors
+ do iproc=proc_p1,proc_p2
+
+ print *,'Reading slice ',iproc,' in region ',iregion_code
+
+! create the name for the database of the current slide
+ call create_name_database(prname,iproc,iregion_code)
+
+ if(ivalue == 1) then
+ open(unit=10,file='OUTPUT_FILES/'//prname(1:len_trim(prname))//'AVS_DXpointsfaces.txt',status='old',action='read')
+ else if(ivalue == 2) then
+ open(unit=10,file='OUTPUT_FILES/'//prname(1:len_trim(prname))//'AVS_DXpointschunks.txt',status='old',action='read')
+ open(unit=12,file='OUTPUT_FILES/'//prname(1:len_trim(prname))//'AVS_DXpointschunks_stability.txt',status='old',action='read')
+ else if(ivalue == 3) then
+ open(unit=10,file='OUTPUT_FILES/'//prname(1:len_trim(prname))//'AVS_DXpointssurface.txt',status='old',action='read')
+ endif
+
+ read(10,*) npoin
+ print *,'There are ',npoin,' global AVS or DX points in the slice'
+
+! read local points in this slice and output global AVS or DX points
+ do ipoin=1,npoin
+ read(10,*) numpoin,xval,yval,zval
+ if(ivalue == 2) then
+ read(12,*) numpoin2,vmin,vmax
+ else
+ numpoin2 = 0
+ vmin = 0.
+ vmax = 0.
+ endif
+ if(numpoin /= ipoin) stop 'incorrect point number'
+ if(ivalue == 2 .and. numpoin2 /= ipoin) stop 'incorrect point number'
+! write to AVS or DX global file with correct offset
+ if(USE_OPENDX) then
+ write(11,"(f10.7,1x,f10.7,1x,f10.7)") xval,yval,zval
+ else
+ write(11,"(i6,1x,f10.7,1x,f10.7,1x,f10.7)") numpoin + iglobpointoffset,xval,yval,zval
+ endif
+
+! save coordinates in global array of points for stability condition
+ xcoord(numpoin + iglobpointoffset) = xval
+ ycoord(numpoin + iglobpointoffset) = yval
+ zcoord(numpoin + iglobpointoffset) = zval
+ vmincoord(numpoin + iglobpointoffset) = vmin
+ vmaxcoord(numpoin + iglobpointoffset) = vmax
+
+ enddo
+
+ iglobpointoffset = iglobpointoffset + npoin
+
+ close(10)
+ if(ivalue == 2) close(12)
+
+ enddo
+ enddo
+
+! ************* generate elements ******************
+
+! get source information for frequency for number of points per lambda
+ print *,'reading source duration from the CMTSOLUTION file'
+ call get_cmt(yr,jda,ho,mi,sec,t_cmt,hdur,lat,long,depth,moment_tensor,DT,1)
+
+! set global element and point offsets to zero
+ iglobpointoffset = 0
+ iglobelemoffset = 0
+ maxdoubling = -1
+ above_zero = 0
+ below_zero = 0
+
+ if(USE_OPENDX) write(11,*) 'object 2 class array type int rank 1 shape 4 items ',ntotspecAVS_DX,' data follows'
+
+ do iregion_code = region_min,region_max
+
+! loop on the selected range of processors
+ do iproc=proc_p1,proc_p2
+
+ print *,'Reading slice ',iproc,' in region ',iregion_code
+
+! create the name for the database of the current slide
+ call create_name_database(prname,iproc,iregion_code)
+
+ if(ivalue == 1) then
+ open(unit=10,file='OUTPUT_FILES/'//prname(1:len_trim(prname))//'AVS_DXelementsfaces.txt',status='old',action='read')
+ open(unit=12,file='OUTPUT_FILES/'//prname(1:len_trim(prname))//'AVS_DXpointsfaces.txt',status='old',action='read')
+ else if(ivalue == 2) then
+ open(unit=10,file='OUTPUT_FILES/'//prname(1:len_trim(prname))//'AVS_DXelementschunks.txt',status='old',action='read')
+ if(icolor == 5 .or. icolor == 6) &
+ open(unit=13,file='OUTPUT_FILES/'//prname(1:len_trim(prname))//'AVS_DXelementschunks_dvp_dvs.txt',status='old',action='read')
+ open(unit=12,file='OUTPUT_FILES/'//prname(1:len_trim(prname))//'AVS_DXpointschunks.txt',status='old',action='read')
+ else if(ivalue == 3) then
+ open(unit=10,file='OUTPUT_FILES/'//prname(1:len_trim(prname))//'AVS_DXelementssurface.txt',status='old',action='read')
+ open(unit=12,file='OUTPUT_FILES/'//prname(1:len_trim(prname))//'AVS_DXpointssurface.txt',status='old',action='read')
+ endif
+
+ read(10,*) nspec
+ read(12,*) npoin
+ print *,'There are ',npoin,' global AVS or DX points in the slice'
+ print *,'There are ',nspec,' AVS or DX elements in the slice'
+
+! read local elements in this slice and output global AVS or DX elements
+ do ispec=1,nspec
+ read(10,*) numelem,idoubling,iglob1,iglob2,iglob3,iglob4
+ if(icolor == 5 .or. icolor == 6) then
+ read(13,*) numelem2,deltavp,deltavs
+ dvp(numelem + iglobelemoffset) = deltavp
+ dvs(numelem + iglobelemoffset) = deltavs
+ else
+ numelem2 = 0
+ endif
+ if(numelem /= ispec) stop 'incorrect element number'
+ if((icolor == 5 .or. icolor == 6) .and. numelem2 /= ispec) stop 'incorrect element number'
+! compute max of the doubling flag
+ maxdoubling = max(maxdoubling,idoubling)
+
+! assign material property (which can be filtered later in AVS_DX)
+ if(imaterial == 1) then
+ imatprop = idoubling
+ else if(imaterial == 2) then
+ imatprop = iproc
+ else if(imaterial == 3) then
+ imatprop = iregion_code
+ else if(imaterial == 4) then
+ imatprop = ichunk_slice(iproc)
+ else
+ stop 'invalid code for material property'
+ endif
+
+! write to AVS or DX global file with correct offset
+
+! quadrangles (2-D)
+ iglob1 = iglob1 + iglobpointoffset
+ iglob2 = iglob2 + iglobpointoffset
+ iglob3 = iglob3 + iglobpointoffset
+ iglob4 = iglob4 + iglobpointoffset
+
+! in the case of OpenDX, node numbers start at zero
+! in the case of AVS, node numbers start at one
+ if(USE_OPENDX) then
+! point order in OpenDX is 1,4,2,3 *not* 1,2,3,4 as in AVS
+ write(11,"(i6,1x,i6,1x,i6,1x,i6)") iglob1-1,iglob4-1,iglob2-1,iglob3-1
+ else
+ write(11,"(i6,1x,i3,' quad ',i6,1x,i6,1x,i6,1x,i6)") numelem + iglobelemoffset,imatprop,iglob1,iglob2,iglob3,iglob4
+ endif
+
+! get number of GLL points in current element
+ NGLL_current_horiz = NGLLX
+ NGLL_current_vert = NGLLZ
+
+! check that the degree is not above the threshold for list of percentages
+ if(NGLL_current_horiz > NGLL_MAX_STABILITY .or. &
+ NGLL_current_vert > NGLL_MAX_STABILITY) &
+ stop 'degree too high to compute stability value'
+
+! scaling factor to compute real value of stability condition
+ scale_factor = dsqrt(PI*GRAV*RHOAV)
+
+! compute stability value
+ stabmax = -1.d0
+ gridmin = HUGEVAL
+
+ if(idoubling == IFLAG_CRUST) then
+
+! distinguish between horizontal and vertical directions in crust
+! because we have a different polynomial degree in each direction
+! this works because the mesher always creates the 2D surfaces starting
+! from the lower-left corner, continuing to the lower-right corner and so on
+ do iloop_corners = 1,2
+
+ select case(iloop_corners)
+
+ case(1)
+ ipointnumber1_horiz = iglob1
+ ipointnumber2_horiz = iglob2
+
+ ipointnumber1_vert = iglob1
+ ipointnumber2_vert = iglob4
+
+ case(2)
+ ipointnumber1_horiz = iglob4
+ ipointnumber2_horiz = iglob3
+
+ ipointnumber1_vert = iglob2
+ ipointnumber2_vert = iglob3
+
+ end select
+
+ distance_horiz = &
+ dsqrt((xcoord(ipointnumber2_horiz)-xcoord(ipointnumber1_horiz))**2 &
+ + (ycoord(ipointnumber2_horiz)-ycoord(ipointnumber1_horiz))**2 &
+ + (zcoord(ipointnumber2_horiz)-zcoord(ipointnumber1_horiz))**2)
+
+ distance_vert = &
+ dsqrt((xcoord(ipointnumber2_vert)-xcoord(ipointnumber1_vert))**2 &
+ + (ycoord(ipointnumber2_vert)-ycoord(ipointnumber1_vert))**2 &
+ + (zcoord(ipointnumber2_vert)-zcoord(ipointnumber1_vert))**2)
+
+! compute stability value using the scaled interval
+ stabmax = dmax1(scale_factor*DT*vmaxcoord(ipointnumber1_horiz)/(distance_horiz*percent_GLL(NGLL_current_horiz)),stabmax)
+ stabmax = dmax1(scale_factor*DT*vmaxcoord(ipointnumber1_vert)/(distance_vert*percent_GLL(NGLL_current_vert)),stabmax)
+
+! compute number of points per wavelength
+ gridmin = dmin1(scale_factor*hdur*vmincoord(ipointnumber1_horiz)*dble(NGLL_current_horiz)/distance_horiz,gridmin)
+ gridmin = dmin1(scale_factor*hdur*vmincoord(ipointnumber1_vert)*dble(NGLL_current_vert)/distance_vert,gridmin)
+
+ enddo
+
+! regular regions with same polynomial degree everywhere
+
+ else
+
+ do istab = 1,4
+ do jstab = 1,4
+ if(jstab /= istab) then
+
+ if(istab == 1) then
+ ipointnumber1_vert = iglob1
+ else if(istab == 2) then
+ ipointnumber1_vert = iglob2
+ else if(istab == 3) then
+ ipointnumber1_vert = iglob3
+ else if(istab == 4) then
+ ipointnumber1_vert = iglob4
+ endif
+
+ if(jstab == 1) then
+ ipointnumber2_vert = iglob1
+ else if(jstab == 2) then
+ ipointnumber2_vert = iglob2
+ else if(jstab == 3) then
+ ipointnumber2_vert = iglob3
+ else if(jstab == 4) then
+ ipointnumber2_vert = iglob4
+ endif
+
+ distance_vert = &
+ dsqrt((xcoord(ipointnumber2_vert)-xcoord(ipointnumber1_vert))**2 &
+ + (ycoord(ipointnumber2_vert)-ycoord(ipointnumber1_vert))**2 &
+ + (zcoord(ipointnumber2_vert)-zcoord(ipointnumber1_vert))**2)
+
+! compute stability value using the scaled interval
+ stabmax = dmax1(scale_factor*DT*vmaxcoord(ipointnumber1_vert)/(distance_vert*percent_GLL(NGLL_current_vert)),stabmax)
+
+! compute number of points per wavelength
+ gridmin = dmin1(scale_factor*hdur*vmincoord(ipointnumber1_vert)*dble(NGLL_current_vert)/distance_vert,gridmin)
+
+ endif
+ enddo
+ enddo
+
+ endif
+
+ stability_value(numelem + iglobelemoffset) = stabmax
+ gridpoints_per_wavelength(numelem + iglobelemoffset) = gridmin
+
+! compute elevation to represent ellipticity or topography at the surface
+! use point iglob1 for instance and subtract reference
+
+! get colatitude and longitude of current point
+ xval = xcoord(iglob1)
+ yval = ycoord(iglob1)
+ zval = zcoord(iglob1)
+
+ call xyz_2_rthetaphi_dble(xval,yval,zval,radius_dummy,theta_s,phi_s)
+ call reduce(theta_s,phi_s)
+
+! if topography then subtract reference ellipsoid or sphere for color code
+! if ellipticity then subtract reference sphere for color code
+! otherwise subtract nothing
+ if(TOPOGRAPHY .or. CRUSTAL) then
+ if(ELLIPTICITY) then
+ reference = 1.d0 - (3.d0*dcos(theta_s)**2 - 1.d0)/3.d0/299.8d0
+ else
+ reference = R_UNIT_SPHERE
+ endif
+ else if(ELLIPTICITY) then
+ reference = R_UNIT_SPHERE
+ else
+ reference = 0.
+ endif
+
+! compute elevation
+ elevation_sphere(numelem + iglobelemoffset) = &
+ (dsqrt(xval**2 + yval**2 + zval**2) - reference)
+
+ enddo
+
+ iglobelemoffset = iglobelemoffset + nspec
+ iglobpointoffset = iglobpointoffset + npoin
+
+ close(10)
+ close(12)
+ if(icolor == 5 .or. icolor == 6) close(13)
+
+ enddo
+ enddo
+
+! saturate color scale for elevation since small values
+! apply non linear scaling if topography to enhance regions around sea level
+
+ if(TOPOGRAPHY .or. CRUSTAL) then
+
+! compute absolute maximum
+ rnorm_factor = maxval(dabs(elevation_sphere(:)))
+
+! map to [-1,1]
+ elevation_sphere(:) = elevation_sphere(:) / rnorm_factor
+
+! apply non-linear scaling
+ do ispec_scale_AVS_DX = 1,ntotspecAVS_DX
+
+ xval = elevation_sphere(ispec_scale_AVS_DX)
+
+! compute total area consisting of oceans
+! and suppress areas that are not considered oceans if needed
+! use arbitrary threshold to suppress artefacts in ETOPO5 model
+ if(xval >= -0.018) then
+ if(OCEANS_ONLY) xval = 0.
+ above_zero = above_zero + 1
+ else
+ below_zero = below_zero + 1
+ endif
+
+ if(xval >= 0.) then
+ if(.not. OCEANS_ONLY) then
+ elevation_sphere(ispec_scale_AVS_DX) = xval ** SCALE_NON_LINEAR
+ else
+ elevation_sphere(ispec_scale_AVS_DX) = 0.
+ endif
+ else
+ elevation_sphere(ispec_scale_AVS_DX) = - dabs(xval) ** SCALE_NON_LINEAR
+ endif
+
+ enddo
+
+ else
+
+! regular scaling to real distance if no topography
+ elevation_sphere(:) = R_EARTH * elevation_sphere(:)
+
+ endif
+
+ if(ISOTROPIC_3D_MANTLE) then
+
+! compute absolute maximum for dvp
+ rnorm_factor = maxval(dabs(dvp(:)))
+
+! map to [-1,1]
+ dvp(:) = dvp(:) / rnorm_factor
+
+! apply non-linear scaling
+ do ispec_scale_AVS_DX = 1,ntotspecAVS_DX
+ xval = dvp(ispec_scale_AVS_DX)
+ if(xval >= 0.) then
+ dvp(ispec_scale_AVS_DX) = xval ** SCALE_NON_LINEAR
+ else
+ dvp(ispec_scale_AVS_DX) = - dabs(xval) ** SCALE_NON_LINEAR
+ endif
+ enddo
+
+! compute absolute maximum for dvs
+ rnorm_factor = maxval(dabs(dvs(:)))
+
+! map to [-1,1]
+ dvs(:) = dvs(:) / rnorm_factor
+
+! apply non-linear scaling
+ do ispec_scale_AVS_DX = 1,ntotspecAVS_DX
+ xval = dvs(ispec_scale_AVS_DX)
+ if(xval >= 0.) then
+ dvs(ispec_scale_AVS_DX) = xval ** SCALE_NON_LINEAR
+ else
+ dvs(ispec_scale_AVS_DX) = - dabs(xval) ** SCALE_NON_LINEAR
+ endif
+ enddo
+
+ endif
+
+! ************* generate element data values ******************
+
+! output AVS or DX header for data
+ if(USE_OPENDX) then
+ write(11,*) 'attribute "element type" string "quads"'
+ write(11,*) 'attribute "ref" string "positions"'
+ write(11,*) 'object 3 class array type float rank 0 items ',ntotspecAVS_DX,' data follows'
+ else
+ write(11,*) '1 1'
+ write(11,*) 'Zcoord, meters'
+ endif
+
+! set global element and point offsets to zero
+ iglobelemoffset = 0
+
+ do iregion_code = region_min,region_max
+
+! loop on the selected range of processors
+ do iproc=proc_p1,proc_p2
+
+ print *,'Reading slice ',iproc,' in region ',iregion_code
+
+! create the name for the database of the current slide
+ call create_name_database(prname,iproc,iregion_code)
+
+ if(ivalue == 1) then
+ open(unit=10,file='OUTPUT_FILES/'//prname(1:len_trim(prname))//'AVS_DXelementsfaces.txt',status='old',action='read')
+ else if(ivalue == 2) then
+ open(unit=10,file='OUTPUT_FILES/'//prname(1:len_trim(prname))//'AVS_DXelementschunks.txt',status='old',action='read')
+ else if(ivalue == 3) then
+ open(unit=10,file='OUTPUT_FILES/'//prname(1:len_trim(prname))//'AVS_DXelementssurface.txt',status='old',action='read')
+ endif
+
+ read(10,*) nspec
+ print *,'There are ',nspec,' AVS or DX elements in the slice'
+
+! read local elements in this slice and output global AVS or DX elements
+ do ispec=1,nspec
+ read(10,*) numelem,idoubling,iglob1,iglob2,iglob3,iglob4
+ if(numelem /= ispec) stop 'incorrect element number'
+
+! data is either the slice number or the mesh doubling region flag
+ if(icolor == 1) then
+ val_color = dble(idoubling)
+ else if(icolor == 2) then
+ val_color = dble(iproc)
+ else if(icolor == 3) then
+ val_color = stability_value(numelem + iglobelemoffset)
+ else if(icolor == 4) then
+ val_color = gridpoints_per_wavelength(numelem + iglobelemoffset)
+! put a threshold for number of points per wavelength displayed
+! otherwise the scale is too large and we cannot see the small values
+ if(val_color > THRESHOLD_GRIDPOINTS) then
+ val_color = THRESHOLD_GRIDPOINTS
+ threshold_used = .true.
+ endif
+ else if(icolor == 5) then
+! minus sign to get the color scheme right: blue is fast (+) and red is slow (-)
+ val_color = -dvp(numelem + iglobelemoffset)
+ else if(icolor == 6) then
+! minus sign to get the color scheme right: blue is fast (+) and red is slow (-)
+ val_color = -dvs(numelem + iglobelemoffset)
+ else if(icolor == 7) then
+ val_color = elevation_sphere(numelem + iglobelemoffset)
+
+ else if(icolor == 8) then
+ val_color = iregion_code
+ else if(icolor == 9) then
+ if(idoubling == itarget_doubling) then
+ val_color = dble(iregion_code)
+ else
+ val_color = dble(IFLAG_DUMMY)
+ endif
+ else
+ stop 'incorrect coloring code'
+ endif
+
+! write to AVS or DX global file with correct offset
+ if(USE_OPENDX) then
+ write(11,*) sngl(val_color)
+ else
+ write(11,*) numelem + iglobelemoffset,' ',sngl(val_color)
+ endif
+ enddo
+
+ iglobelemoffset = iglobelemoffset + nspec
+
+ close(10)
+
+ enddo
+ enddo
+
+! define OpenDX field
+ if(USE_OPENDX) then
+ write(11,*) 'attribute "dep" string "connections"'
+ write(11,*) 'object "irregular positions irregular connections" class field'
+ write(11,*) 'component "positions" value 1'
+ write(11,*) 'component "connections" value 2'
+ write(11,*) 'component "data" value 3'
+ write(11,*) 'end'
+ endif
+
+ close(11)
+
+ print *
+ print *,'maximum value of doubling flag in all slices = ',maxdoubling
+ print *
+
+! print min and max of stability and points per lambda
+
+ if(ivalue == 2) then
+
+! compute minimum and maximum of stability value and points per wavelength
+
+ stability_value_min = minval(stability_value)
+ stability_value_max = maxval(stability_value)
+
+ gridpoints_per_wavelength_min = minval(gridpoints_per_wavelength)
+ gridpoints_per_wavelength_max = maxval(gridpoints_per_wavelength)
+
+ print *
+ print *,'stability value min, max, ratio = ', &
+ stability_value_min,stability_value_max,stability_value_max / stability_value_min
+
+ print *
+ print *,'number of points per wavelength min, max, ratio = ', &
+ gridpoints_per_wavelength_min,gridpoints_per_wavelength_max, &
+ gridpoints_per_wavelength_max / gridpoints_per_wavelength_min
+
+ print *
+ print *,'half duration of ',sngl(hdur),' s used for points per wavelength'
+ print *
+
+ if(hdur < 5.*DT) then
+ print *,'***************************************************************'
+ print *,'Source time function is a Heaviside, points per wavelength meaningless'
+ print *,'***************************************************************'
+ print *
+ endif
+
+ if(icolor == 4 .and. threshold_used) then
+ print *,'***************************************************************'
+ print *,'the number of points per wavelength have been cut above a threshold'
+ print *,'of ',THRESHOLD_GRIDPOINTS,' to avoid saturation of color scale'
+ print *,'***************************************************************'
+ print *
+ endif
+ endif
+
+! if we have the surface for the Earth, print min and max of elevation
+
+ if(ivalue == 3) then
+ print *
+ print *,'elevation min, max = ',minval(elevation_sphere),maxval(elevation_sphere)
+ if(TOPOGRAPHY .or. CRUSTAL) print *,'elevation has been normalized for topography'
+ print *
+
+! print percentage of oceans at surface of the globe
+ print *
+ print *,'the oceans represent ',100. * below_zero / (above_zero + below_zero),' % of the surface of the mesh'
+ print *
+
+ endif
+
+!
+! create an AVS or DX file with the source and the receivers as well
+!
+
+! get source information
+ print *,'reading position of the source from the CMTSOLUTION file'
+ call get_cmt(yr,jda,ho,mi,sec,t_cmt,hdur,lat,long,depth,moment_tensor,DT,1)
+
+! convert geographic latitude lat (degrees)
+! to geocentric colatitude theta (radians)
+ theta=PI/2.0d0-atan(0.99329534d0*tan(dble(lat)*PI/180.0d0))
+ phi=dble(long)*PI/180.0d0
+ call reduce(theta,phi)
+
+! compute Cartesian position of the source (ignore ellipticity for AVS_DX)
+! the point for the source is put at the surface for clarity (depth ignored)
+! even slightly above to superimpose to real surface
+ r_target_source = 1.02d0
+ x_target_source = r_target_source*sin(theta)*cos(phi)
+ y_target_source = r_target_source*sin(theta)*sin(phi)
+ z_target_source = r_target_source*cos(theta)
+
+! save triangle for AVS or DX representation of epicenter
+ r_target_source = 1.05d0
+ delta_trgl = 1.8 * pi / 180.
+ x_source_trgl1 = r_target_source*sin(theta+delta_trgl)*cos(phi-delta_trgl)
+ y_source_trgl1 = r_target_source*sin(theta+delta_trgl)*sin(phi-delta_trgl)
+ z_source_trgl1 = r_target_source*cos(theta+delta_trgl)
+
+ x_source_trgl2 = r_target_source*sin(theta+delta_trgl)*cos(phi+delta_trgl)
+ y_source_trgl2 = r_target_source*sin(theta+delta_trgl)*sin(phi+delta_trgl)
+ z_source_trgl2 = r_target_source*cos(theta+delta_trgl)
+
+ x_source_trgl3 = r_target_source*sin(theta-delta_trgl)*cos(phi)
+ y_source_trgl3 = r_target_source*sin(theta-delta_trgl)*sin(phi)
+ z_source_trgl3 = r_target_source*cos(theta-delta_trgl)
+
+ ntotpoinAVS_DX = 2
+ ntotspecAVS_DX = 1
+
+ print *
+ print *,'reading position of the receivers'
+
+! get number of stations from receiver file
+ open(unit=11,file='DATA/STATIONS',iostat=ios,status='old',action='read')
+ nrec = 0
+ do while(ios == 0)
+ read(11,"(a)",iostat=ios) dummystring
+ if(ios == 0) nrec = nrec + 1
+ enddo
+ close(11)
+
+ print *,'There are ',nrec,' three-component stations'
+ print *
+ if(nrec < 1) stop 'incorrect number of stations read - need at least one'
+
+ allocate(station_name(nrec))
+ allocate(network_name(nrec))
+ allocate(stlat(nrec))
+ allocate(stlon(nrec))
+ allocate(stele(nrec))
+ allocate(stbur(nrec))
+
+ allocate(x_target(nrec))
+ allocate(y_target(nrec))
+ allocate(z_target(nrec))
+
+! loop on all the stations
+ open(unit=11,file='DATA/STATIONS',status='old',action='read')
+ do irec=1,nrec
+ read(11,*) station_name(irec),network_name(irec),stlat(irec),stlon(irec),stele(irec),stbur(irec)
+
+! convert geographic latitude stlat (degrees)
+! to geocentric colatitude theta (radians)
+ theta=PI/2.0d0-atan(0.99329534d0*dtan(stlat(irec)*PI/180.0d0))
+ phi=stlon(irec)*PI/180.0d0
+ call reduce(theta,phi)
+
+! compute the Cartesian position of the receiver (ignore ellipticity for AVS_DX)
+! points for the receivers are put at the surface for clarity (depth ignored)
+ r_target=1.0d0
+ x_target(irec) = r_target*dsin(theta)*dcos(phi)
+ y_target(irec) = r_target*dsin(theta)*dsin(phi)
+ z_target(irec) = r_target*dcos(theta)
+
+ enddo
+
+ close(11)
+
+! duplicate source to have right color normalization in AVS_DX
+ ntotpoinAVS_DX = ntotpoinAVS_DX + 2*nrec + 1
+ ntotspecAVS_DX = ntotspecAVS_DX + nrec + 1
+
+! write AVS or DX header with element data
+! add source and receivers (small AVS or DX lines)
+! duplicate source to have right color normalization in AVS_DX
+ if(USE_OPENDX) then
+ open(unit=11,file=trim(OUTPUT_FILES)//'/DX_source_receivers.dx',status='unknown')
+ write(11,*) 'object 1 class array type float rank 1 shape 3 items ',ntotpoinAVS_DX,' data follows'
+ write(11,*) sngl(x_target_source),' ',sngl(y_target_source),' ',sngl(z_target_source)
+ write(11,*) sngl(x_target_source+0.1*small_offset_source),' ', &
+ sngl(y_target_source+0.1*small_offset_source),' ',sngl(z_target_source+0.1*small_offset_source)
+ write(11,*) sngl(x_target_source+1.3*small_offset_source),' ', &
+ sngl(y_target_source+1.3*small_offset_source),' ',sngl(z_target_source+1.3*small_offset_source)
+ do ir=1,nrec
+ write(11,*) sngl(x_target(ir)),' ',sngl(y_target(ir)),' ',sngl(z_target(ir))
+ write(11,*) sngl(x_target(ir)+small_offset_receiver),' ', &
+ sngl(y_target(ir)+small_offset_receiver),' ',sngl(z_target(ir)+small_offset_receiver)
+ enddo
+ write(11,*) 'object 2 class array type int rank 1 shape 2 items ',ntotspecAVS_DX,' data follows'
+ write(11,*) '0 1'
+ do ir=1,nrec
+ write(11,*) 4+2*(ir-1)-1,' ',4+2*(ir-1)
+ enddo
+ write(11,*) '0 2'
+ write(11,*) 'attribute "element type" string "lines"'
+ write(11,*) 'attribute "ref" string "positions"'
+ write(11,*) 'object 3 class array type float rank 0 items ',ntotspecAVS_DX,' data follows'
+ write(11,*) '1.'
+ do ir=1,nrec
+ write(11,*) ' 255.'
+ enddo
+ write(11,*) ' 120.'
+ write(11,*) 'attribute "dep" string "connections"'
+ write(11,*) 'object "irregular connections irregular positions" class field'
+ write(11,*) 'component "positions" value 1'
+ write(11,*) 'component "connections" value 2'
+ write(11,*) 'component "data" value 3'
+ write(11,*) 'end'
+ close(11)
+ else
+ open(unit=11,file=trim(OUTPUT_FILES)//'/AVS_source_receivers.inp',status='unknown')
+ write(11,*) ntotpoinAVS_DX,' ',ntotspecAVS_DX,' 0 1 0'
+ write(11,*) '1 ',sngl(x_target_source),' ',sngl(y_target_source),' ',sngl(z_target_source)
+ write(11,*) '2 ',sngl(x_target_source+0.1*small_offset_source),' ', &
+ sngl(y_target_source+0.1*small_offset_source),' ',sngl(z_target_source+0.1*small_offset_source)
+ write(11,*) '3 ',sngl(x_target_source+1.3*small_offset_source),' ', &
+ sngl(y_target_source+1.3*small_offset_source),' ',sngl(z_target_source+1.3*small_offset_source)
+ do ir=1,nrec
+ write(11,*) 4+2*(ir-1),' ',sngl(x_target(ir)),' ',sngl(y_target(ir)),' ',sngl(z_target(ir))
+ write(11,*) 4+2*(ir-1)+1,' ',sngl(x_target(ir)+small_offset_receiver),' ', &
+ sngl(y_target(ir)+small_offset_receiver),' ',sngl(z_target(ir)+small_offset_receiver)
+ enddo
+ write(11,*) '1 1 line 1 2'
+ do ir=1,nrec
+ write(11,*) ir+1,' 1 line ',4+2*(ir-1),' ',4+2*(ir-1)+1
+ enddo
+ write(11,*) ir+1,' 1 line 1 3'
+ write(11,*) '1 1'
+ write(11,*) 'Zcoord, meters'
+ write(11,*) '1 1.'
+ do ir=1,nrec
+ write(11,*) ir+1,' 255.'
+ enddo
+ write(11,*) ir+1,' 120.'
+ close(11)
+ endif
+
+! create a file with the epicenter only, represented as a triangle
+
+! write AVS or DX header with element data
+ if(USE_OPENDX) then
+ open(unit=11,file=trim(OUTPUT_FILES)//'/DX_epicenter.dx',status='unknown')
+ write(11,*) 'object 1 class array type float rank 1 shape 3 items 3 data follows'
+ write(11,*) sngl(x_source_trgl1),' ',sngl(y_source_trgl1),' ',sngl(z_source_trgl1)
+ write(11,*) sngl(x_source_trgl2),' ',sngl(y_source_trgl2),' ',sngl(z_source_trgl2)
+ write(11,*) sngl(x_source_trgl3),' ',sngl(y_source_trgl3),' ',sngl(z_source_trgl3)
+ write(11,*) 'object 2 class array type int rank 1 shape 3 items 1 data follows'
+ write(11,*) '0 1 2'
+ write(11,*) 'attribute "element type" string "triangles"'
+ write(11,*) 'attribute "ref" string "positions"'
+ write(11,*) 'object 3 class array type float rank 0 items 1 data follows'
+ write(11,*) '1.'
+ write(11,*) 'attribute "dep" string "connections"'
+ write(11,*) 'object "irregular connections irregular positions" class field'
+ write(11,*) 'component "positions" value 1'
+ write(11,*) 'component "connections" value 2'
+ write(11,*) 'component "data" value 3'
+ write(11,*) 'end'
+ close(11)
+ else
+ open(unit=11,file=trim(OUTPUT_FILES)//'/AVS_epicenter.inp',status='unknown')
+ write(11,*) '3 1 0 1 0'
+ write(11,*) '1 ',sngl(x_source_trgl1),' ',sngl(y_source_trgl1),' ',sngl(z_source_trgl1)
+ write(11,*) '2 ',sngl(x_source_trgl2),' ',sngl(y_source_trgl2),' ',sngl(z_source_trgl2)
+ write(11,*) '3 ',sngl(x_source_trgl3),' ',sngl(y_source_trgl3),' ',sngl(z_source_trgl3)
+ write(11,*) '1 1 tri 1 2 3'
+ write(11,*) '1 1'
+ write(11,*) 'Zcoord, meters'
+ write(11,*) '1 1.'
+ close(11)
+ endif
+
+ end program combine_AVS_DX
+
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/create_header_file.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/create_header_file.f90 2008-08-15 00:52:09 UTC (rev 12641)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/create_header_file.f90 2008-08-15 03:35:19 UTC (rev 12642)
@@ -108,7 +108,7 @@
call get_value_string(HEADER_FILE, 'solver.HEADER_FILE', 'OUTPUT_FILES/values_from_mesher.h')
print *
- print *,'creating file ', trim(HEADER_FILE), ' to compile solver with correct values'
+ print *,'creating file ', trim(HEADER_FILE), ' to compile the code with correct static values'
! read the parameter file and compute additional parameters
call read_compute_parameters(MIN_ATTENUATION_PERIOD,MAX_ATTENUATION_PERIOD,NER_CRUST, &
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/create_name_database.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/create_name_database.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/create_name_database.f90 2008-08-15 03:35:19 UTC (rev 12642)
@@ -0,0 +1,43 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 4 . 1
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory, California Institute of Technology, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) California Institute of Technology and University of Pau / CNRS / INRIA
+! August 2008
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+ subroutine create_name_database(prname,myrank,iregion_code)
+
+! create the name of the database for a given processor and given mesh region
+
+ implicit none
+
+ integer myrank,iregion_code
+
+! name of the database file
+ character(len=150) prname
+
+! create the name for the database of the current slide and region
+ write(prname,"('proc',i6.6,'_reg',i1,'_')") myrank,iregion_code
+
+ end subroutine create_name_database
+
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/create_regions_mesh.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/create_regions_mesh.F90 2008-08-15 00:52:09 UTC (rev 12641)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/create_regions_mesh.F90 2008-08-15 03:35:19 UTC (rev 12642)
@@ -50,7 +50,7 @@
xread1D_leftxi_righteta,xread1D_rightxi_righteta,yread1D_leftxi_lefteta,yread1D_rightxi_lefteta,yread1D_leftxi_righteta, &
yread1D_rightxi_righteta,zread1D_leftxi_lefteta,zread1D_rightxi_lefteta,zread1D_leftxi_righteta,zread1D_rightxi_righteta, &
#endif
- rho_vp,rho_vs,Qmu_store,tau_e_store,ifirst_layer_aniso,ilast_layer_aniso)
+ rho_vp,rho_vs,Qmu_store,tau_e_store,ifirst_layer_aniso,ilast_layer_aniso,SAVE_MESH_FILES)
! create the different regions of the mesh
@@ -373,7 +373,7 @@
integer npointot
- logical ELLIPTICITY,TOPOGRAPHY
+ logical ELLIPTICITY,TOPOGRAPHY,SAVE_MESH_FILES
logical TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE,ISOTROPIC_3D_MANTLE,CRUSTAL,ONE_CRUST,OCEANS
logical ATTENUATION,ATTENUATION_3D,INCLUDE_CENTRAL_CUBE,ABSORBING_CONDITIONS
@@ -381,7 +381,7 @@
double precision R_CENTRAL_CUBE,RICB,RHO_OCEANS,RCMB,R670,RMOHO, &
RTOPDDOUBLEPRIME,R600,R220,R771,R400,R120,R80,RMIDDLE_CRUST,ROCEAN
- character(len=150) errmsg
+ character(len=150) errmsg,prname
! use integer array to store values
integer ibathy_topo(NX_BATHY,NY_BATHY)
@@ -1291,7 +1291,12 @@
enddo
enddo
- if(minval(ibool) /= 1 .or. maxval(ibool) /= nglob_theor) call exit_MPI(myrank,'incorrect global numbering after sorting')
+ if(minval(ibool) /= 1 .or. maxval(ibool) /= nglob_theor) then
+ print *,'Incorrect global numbering after sorting: mask_ibool ',minval(mask_ibool),maxval(mask_ibool)
+ print *,'Incorrect global numbering after sorting: ibool ',minval(ibool),maxval(ibool)
+ print *,'Incorrect global numbering after sorting: expected min/max = 1, ',nglob_theor
+ call exit_MPI(myrank,'incorrect global numbering after sorting')
+ endif
! create MPI buffers
! arrays locval(npointot) and ifseg(npointot) used to save memory
@@ -1315,6 +1320,28 @@
iregion_code)
#endif
+! create AVS or DX mesh data for the slices
+ if(SAVE_MESH_FILES) then
+ call create_name_database(prname,myrank,iregion_code)
+
+ call write_AVS_DX_global_data(myrank,prname,nspec,ibool,idoubling,xstore,ystore,zstore,locval,ifseg,npointot,iregion_code)
+
+ call write_AVS_DX_global_faces_data(myrank,prname,nspec,iMPIcut_xi,iMPIcut_eta,ibool, &
+ idoubling,xstore,ystore,zstore,locval,ifseg,npointot,iregion_code)
+
+ call write_AVS_DX_global_chunks_data(myrank,prname,nspec,iboun,ibool, &
+ idoubling,xstore,ystore,zstore,locval,ifseg,npointot, &
+!! DK DK changed for now because array rhostore is not available in v4.1 anymore
+!! DK DK rhostore,kappavstore,muvstore,nspl,rspl,espl,espl2, &
+ nspl,rspl,espl,espl2, &
+ ELLIPTICITY,ISOTROPIC_3D_MANTLE,CRUSTAL,ONE_CRUST,REFERENCE_1D_MODEL, &
+ RICB,RCMB,RTOPDDOUBLEPRIME,R600,R670,R220,R771,R400,R120,R80,RMOHO, &
+ RMIDDLE_CRUST,ROCEAN,M1066a_V,Mak135_V,Mref_V,SEA1DM_V,iregion_code)
+
+ call write_AVS_DX_surface_data(myrank,prname,nspec,iboun,ibool, &
+ idoubling,xstore,ystore,zstore,locval,ifseg,npointot,iregion_code)
+ endif
+
! only create mass matrix and save all the final arrays in the second pass
else if(ipass == 2) then
@@ -1430,7 +1457,7 @@
!! DK DK save Brian's attenuation files to a shared disk
!! DK DK obviously we should do this with MPI or with subroutine arguments
!! DK DK shared by the mesher and the solver subroutines at some point
- call attenuation_save_arrays(iregion_code, AM_V)
+ if(ATTENUATION_VAL) call attenuation_save_arrays(iregion_code, AM_V)
! compute volume, bottom and top area of that part of the slice
volume_local = ZERO
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/main_program.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/main_program.F90 2008-08-15 00:52:09 UTC (rev 12641)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/main_program.F90 2008-08-15 03:35:19 UTC (rev 12642)
@@ -416,6 +416,8 @@
! YYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYYY
+ if(.not. MESHER_ONLY) then
+
!! DK DK for merged version, temporary patch for David's code to cut the superbrick
!! DK DK which I have not fully ported to the merged version yet: I do not
!! DK DK yet distinguish the two values of each array, therefore let me set them
@@ -482,11 +484,15 @@
#endif
AM_V)
+! synchronize all the processes to make sure everybody has finished
#ifdef USE_MPI
-! synchronize all the processes to make sure everybody has finished
call MPI_BARRIER(MPI_COMM_WORLD,ier)
+#endif
+ endif ! end of test if running the mesher only but not the solver
+
! stop all the MPI processes, and exit
+#ifdef USE_MPI
call MPI_FINALIZE(ier)
#endif
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/meshfem3D.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/meshfem3D.F90 2008-08-15 00:52:09 UTC (rev 12641)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/meshfem3D.F90 2008-08-15 03:35:19 UTC (rev 12642)
@@ -902,8 +902,8 @@
write(IMAIN,*)
endif
do ichunk = 1,NCHUNKS
- do iproc_eta=0,NPROC_ETA-1
- do iproc_xi=0,NPROC_XI-1
+ do iproc_eta = 0,NPROC_ETA-1
+ do iproc_xi = 0,NPROC_XI-1
iprocnum = (ichunk-1)*NPROC + iproc_eta * NPROC_XI + iproc_xi
addressing(ichunk,iproc_xi,iproc_eta) = iprocnum
ichunk_slice(iprocnum) = ichunk
@@ -1455,7 +1455,7 @@
xread1D_leftxi_righteta,xread1D_rightxi_righteta,yread1D_leftxi_lefteta,yread1D_rightxi_lefteta,yread1D_leftxi_righteta, &
yread1D_rightxi_righteta,zread1D_leftxi_lefteta,zread1D_rightxi_lefteta,zread1D_leftxi_righteta,zread1D_rightxi_righteta, &
#endif
- rho_vp,rho_vs,Qmu_store,tau_e_store,ifirst_layer_aniso,ilast_layer_aniso)
+ rho_vp,rho_vs,Qmu_store,tau_e_store,ifirst_layer_aniso,ilast_layer_aniso,SAVE_MESH_FILES)
else if(iregion_code == IREGION_OUTER_CORE) then
! outer_core
@@ -1490,7 +1490,7 @@
xread1D_leftxi_righteta,xread1D_rightxi_righteta,yread1D_leftxi_lefteta,yread1D_rightxi_lefteta,yread1D_leftxi_righteta, &
yread1D_rightxi_righteta,zread1D_leftxi_lefteta,zread1D_rightxi_lefteta,zread1D_leftxi_righteta,zread1D_rightxi_righteta, &
#endif
- rho_vp,rho_vs,Qmu_store,tau_e_store,ifirst_layer_aniso,ilast_layer_aniso)
+ rho_vp,rho_vs,Qmu_store,tau_e_store,ifirst_layer_aniso,ilast_layer_aniso,SAVE_MESH_FILES)
else if(iregion_code == IREGION_INNER_CORE) then
! inner_core
@@ -1525,7 +1525,7 @@
xread1D_leftxi_righteta,xread1D_rightxi_righteta,yread1D_leftxi_lefteta,yread1D_rightxi_lefteta,yread1D_leftxi_righteta, &
yread1D_rightxi_righteta,zread1D_leftxi_lefteta,zread1D_rightxi_lefteta,zread1D_leftxi_righteta,zread1D_rightxi_righteta, &
#endif
- rho_vp,rho_vs,Qmu_store,tau_e_store,ifirst_layer_aniso,ilast_layer_aniso)
+ rho_vp,rho_vs,Qmu_store,tau_e_store,ifirst_layer_aniso,ilast_layer_aniso,SAVE_MESH_FILES)
else
stop 'DK DK incorrect region in merged code'
@@ -1754,15 +1754,24 @@
if(myrank == 0) then
#ifdef USE_MPI
tCPU = MPI_WTIME() - time_start
+#else
+ tCPU = 0
+#endif
+
write(IMAIN,*)
+#ifdef USE_MPI
write(IMAIN,*) 'Elapsed time for mesh generation and buffer creation in seconds = ',tCPU
+#endif
write(IMAIN,*) 'End of mesh generation'
write(IMAIN,*)
-#else
- tCPU = 0
-#endif
+
+ if(MESHER_ONLY) then
+ write(IMAIN,*) 'Generation of the mesh only, not starting the solver because MESHER_ONLY is true'
+ write(IMAIN,*)
+ endif
+
! close main output file
- if(myrank == 0 .and. IMAIN /= ISTANDARD_OUTPUT) close(IMAIN)
+ if(IMAIN /= ISTANDARD_OUTPUT) close(IMAIN)
endif
end subroutine meshfem3D
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/write_AVS_DX_global_chunks_data.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/write_AVS_DX_global_chunks_data.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/write_AVS_DX_global_chunks_data.f90 2008-08-15 03:35:19 UTC (rev 12642)
@@ -0,0 +1,729 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 4 . 1
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory, California Institute of Technology, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) California Institute of Technology and University of Pau / CNRS / INRIA
+! August 2008
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+! create AVS or DX 2D data for the faces of the global chunks,
+! to be recombined in postprocessing
+ subroutine write_AVS_DX_global_chunks_data(myrank,prname,nspec,iboun, &
+ ibool,idoubling,xstore,ystore,zstore,num_ibool_AVS_DX,mask_ibool, &
+!! DK DK changed for now because array rhostore is not available in v4.1 anymore
+!! DK DK npointot,rhostore,kappavstore,muvstore,nspl,rspl,espl,espl2, &
+ npointot,nspl,rspl,espl,espl2, &
+ ELLIPTICITY,ISOTROPIC_3D_MANTLE,CRUSTAL,ONE_CRUST,REFERENCE_1D_MODEL, &
+ RICB,RCMB,RTOPDDOUBLEPRIME,R600,R670,R220,R771,R400,R120,R80,RMOHO, &
+ RMIDDLE_CRUST,ROCEAN,M1066a_V,Mak135_V,Mref_V,SEA1DM_V,iregion_code)
+
+ implicit none
+
+ include "constants.h"
+
+ integer nspec,myrank,REFERENCE_1D_MODEL,iregion_code
+ integer ibool(NGLLX,NGLLY,NGLLZ,nspec)
+
+ integer idoubling(nspec)
+
+ logical iboun(6,nspec),ELLIPTICITY,ISOTROPIC_3D_MANTLE,CRUSTAL,ONE_CRUST
+
+ double precision RICB,RCMB,RTOPDDOUBLEPRIME,R600,R670,R220,R771,R400,R120,R80,RMOHO,RMIDDLE_CRUST,ROCEAN
+
+ double precision xstore(NGLLX,NGLLY,NGLLZ,nspec)
+ double precision ystore(NGLLX,NGLLY,NGLLZ,nspec)
+ double precision zstore(NGLLX,NGLLY,NGLLZ,nspec)
+
+!! DK DK changed for now because array rhostore is not available in v4.1 anymore
+!! DK DK real(kind=CUSTOM_REAL) rhostore(NGLLX,NGLLY,NGLLZ,nspec)
+!! DK DK real(kind=CUSTOM_REAL) kappavstore(NGLLX,NGLLY,NGLLZ,nspec)
+!! DK DK real(kind=CUSTOM_REAL) muvstore(NGLLX,NGLLY,NGLLZ,nspec)
+
+! logical mask used to output global points only once
+ integer npointot
+ logical mask_ibool(npointot)
+
+! numbering of global AVS or DX points
+ integer num_ibool_AVS_DX(npointot)
+
+ integer ispec,idoubling_to_use
+ integer i,j,k,np
+ integer, dimension(8) :: iglobval
+ integer npoin,numpoin,nspecface,ispecface
+
+ real(kind=CUSTOM_REAL) vmin,vmax
+
+ double precision r,rho,vp,vs,Qkappa,Qmu
+ double precision vpv,vph,vsv,vsh,eta_aniso
+ double precision x,y,z,theta,phi_dummy,cost,p20,ell,factor
+ real(kind=CUSTOM_REAL) dvp,dvs
+
+! for ellipticity
+ integer nspl
+ double precision rspl(NR),espl(NR),espl2(NR)
+
+! processor identification
+ character(len=150) prname
+
+! model_1066a_variables
+ type model_1066a_variables
+ sequence
+ double precision, dimension(NR_1066A) :: radius_1066a
+ double precision, dimension(NR_1066A) :: density_1066a
+ double precision, dimension(NR_1066A) :: vp_1066a
+ double precision, dimension(NR_1066A) :: vs_1066a
+ double precision, dimension(NR_1066A) :: Qkappa_1066a
+ double precision, dimension(NR_1066A) :: Qmu_1066a
+ end type model_1066a_variables
+
+ type (model_1066a_variables) M1066a_V
+! model_1066a_variables
+
+! model_ak135_variables
+ type model_ak135_variables
+ sequence
+ double precision, dimension(NR_AK135) :: radius_ak135
+ double precision, dimension(NR_AK135) :: density_ak135
+ double precision, dimension(NR_AK135) :: vp_ak135
+ double precision, dimension(NR_AK135) :: vs_ak135
+ double precision, dimension(NR_AK135) :: Qkappa_ak135
+ double precision, dimension(NR_AK135) :: Qmu_ak135
+ end type model_ak135_variables
+
+ type (model_ak135_variables) Mak135_V
+! model_ak135_variables
+
+! model_ref_variables
+ type model_ref_variables
+ sequence
+ double precision, dimension(NR_REF) :: radius_ref
+ double precision, dimension(NR_REF) :: density_ref
+ double precision, dimension(NR_REF) :: vpv_ref
+ double precision, dimension(NR_REF) :: vph_ref
+ double precision, dimension(NR_REF) :: vsv_ref
+ double precision, dimension(NR_REF) :: vsh_ref
+ double precision, dimension(NR_REF) :: eta_ref
+ double precision, dimension(NR_REF) :: Qkappa_ref
+ double precision, dimension(NR_REF) :: Qmu_ref
+ end type model_ref_variables
+
+ type (model_ref_variables) Mref_V
+! model_ref_variables
+
+! sea1d_model_variables
+ type sea1d_model_variables
+ sequence
+ double precision, dimension(NR_SEA1D) :: radius_sea1d
+ double precision, dimension(NR_SEA1D) :: density_sea1d
+ double precision, dimension(NR_SEA1D) :: vp_sea1d
+ double precision, dimension(NR_SEA1D) :: vs_sea1d
+ double precision, dimension(NR_SEA1D) :: Qkappa_sea1d
+ double precision, dimension(NR_SEA1D) :: Qmu_sea1d
+ end type sea1d_model_variables
+
+ type (sea1d_model_variables) SEA1DM_V
+! sea1d_model_variables
+
+! writing points
+ open(unit=10,file='OUTPUT_FILES/'//prname(1:len_trim(prname))//'AVS_DXpointschunks.txt',status='unknown')
+ open(unit=11,file='OUTPUT_FILES/'//prname(1:len_trim(prname))//'AVS_DXpointschunks_stability.txt',status='unknown')
+
+! erase the logical mask used to mark points already found
+ mask_ibool(:) = .false.
+
+ nspecface = 0
+
+! mark global AVS or DX points
+ do ispec=1,nspec
+! only if on face
+ if(iboun(1,ispec) .or. iboun(2,ispec) .or. &
+ iboun(3,ispec) .or. iboun(4,ispec)) then
+ iglobval(1)=ibool(1,1,1,ispec)
+ iglobval(2)=ibool(NGLLX,1,1,ispec)
+ iglobval(3)=ibool(NGLLX,NGLLY,1,ispec)
+ iglobval(4)=ibool(1,NGLLY,1,ispec)
+ iglobval(5)=ibool(1,1,NGLLZ,ispec)
+ iglobval(6)=ibool(NGLLX,1,NGLLZ,ispec)
+ iglobval(7)=ibool(NGLLX,NGLLY,NGLLZ,ispec)
+ iglobval(8)=ibool(1,NGLLY,NGLLZ,ispec)
+
+! face xi = xi_min
+ if(iboun(1,ispec)) then
+ nspecface = nspecface + 1
+ mask_ibool(iglobval(1)) = .true.
+ mask_ibool(iglobval(4)) = .true.
+ mask_ibool(iglobval(8)) = .true.
+ mask_ibool(iglobval(5)) = .true.
+ endif
+
+! face xi = xi_max
+ if(iboun(2,ispec)) then
+ nspecface = nspecface + 1
+ mask_ibool(iglobval(2)) = .true.
+ mask_ibool(iglobval(3)) = .true.
+ mask_ibool(iglobval(7)) = .true.
+ mask_ibool(iglobval(6)) = .true.
+ endif
+
+! face eta = eta_min
+ if(iboun(3,ispec)) then
+ nspecface = nspecface + 1
+ mask_ibool(iglobval(1)) = .true.
+ mask_ibool(iglobval(2)) = .true.
+ mask_ibool(iglobval(6)) = .true.
+ mask_ibool(iglobval(5)) = .true.
+ endif
+
+! face eta = eta_max
+ if(iboun(4,ispec)) then
+ nspecface = nspecface + 1
+ mask_ibool(iglobval(4)) = .true.
+ mask_ibool(iglobval(3)) = .true.
+ mask_ibool(iglobval(7)) = .true.
+ mask_ibool(iglobval(8)) = .true.
+ endif
+
+ endif
+ enddo
+
+! count global number of AVS or DX points
+ npoin = count(mask_ibool(:))
+
+! number of points in AVS or DX file
+ write(10,*) npoin
+
+! erase the logical mask used to mark points already found
+ mask_ibool(:) = .false.
+
+! output global AVS or DX points
+ numpoin = 0
+
+ do ispec = 1,nspec
+
+! only if on face
+ if(iboun(1,ispec) .or. iboun(2,ispec) .or. iboun(3,ispec) .or. iboun(4,ispec)) then
+
+ iglobval(1)=ibool(1,1,1,ispec)
+ iglobval(2)=ibool(NGLLX,1,1,ispec)
+ iglobval(3)=ibool(NGLLX,NGLLY,1,ispec)
+ iglobval(4)=ibool(1,NGLLY,1,ispec)
+ iglobval(5)=ibool(1,1,NGLLZ,ispec)
+ iglobval(6)=ibool(NGLLX,1,NGLLZ,ispec)
+ iglobval(7)=ibool(NGLLX,NGLLY,NGLLZ,ispec)
+ iglobval(8)=ibool(1,NGLLY,NGLLZ,ispec)
+
+ if(iregion_code == IREGION_OUTER_CORE) then
+! doubling array is not created in the outer core to save memory because it is useless
+! therefore impose the doubling flag manually if needed
+ idoubling_to_use = IFLAG_OUTER_CORE_NORMAL
+ else
+ idoubling_to_use = idoubling(ispec)
+ endif
+
+! face xi = xi_min
+ if(iboun(1,ispec)) then
+
+ if(.not. mask_ibool(iglobval(1))) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglobval(1)) = numpoin
+ write(10,*) numpoin,sngl(xstore(1,1,1,ispec)), &
+ sngl(ystore(1,1,1,ispec)),sngl(zstore(1,1,1,ispec))
+ vmax = 333 !! sqrt((kappavstore(1,1,1,ispec)+4.*muvstore(1,1,1,ispec)/3.)/rhostore(1,1,1,ispec))
+ vmin = 333 !! sqrt(muvstore(1,1,1,ispec)/rhostore(1,1,1,ispec))
+! particular case of the outer core (muvstore contains 1/rho)
+ if(iregion_code == IREGION_OUTER_CORE) then
+ r = dsqrt(xstore(1,1,1,ispec)**2 + ystore(1,1,1,ispec)**2 + zstore(1,1,1,ispec)**2)
+ call prem_display_outer_core(myrank,r,rho,vp,vs,Qkappa,Qmu,idoubling_to_use)
+ vmax = vp
+ vmin = vp
+ endif
+ if(vmin == 0.0) vmin=vmax
+ write(11,*) numpoin,vmin,vmax
+ endif
+
+ if(.not. mask_ibool(iglobval(4))) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglobval(4)) = numpoin
+ write(10,*) numpoin,sngl(xstore(1,NGLLY,1,ispec)), &
+ sngl(ystore(1,NGLLY,1,ispec)),sngl(zstore(1,NGLLY,1,ispec))
+ vmax = 333 !! sqrt((kappavstore(1,NGLLY,1,ispec)+4.*muvstore(1,NGLLY,1,ispec)/3.)/rhostore(1,NGLLY,1,ispec))
+ vmin = 333 !! sqrt(muvstore(1,NGLLY,1,ispec)/rhostore(1,NGLLY,1,ispec))
+! particular case of the outer core (muvstore contains 1/rho)
+ if(iregion_code == IREGION_OUTER_CORE) then
+ r = dsqrt(xstore(1,NGLLY,1,ispec)**2 + ystore(1,NGLLY,1,ispec)**2 + zstore(1,NGLLY,1,ispec)**2)
+ call prem_display_outer_core(myrank,r,rho,vp,vs,Qkappa,Qmu,idoubling_to_use)
+ vmax = vp
+ vmin = vp
+ endif
+ if(vmin == 0.0) vmin=vmax
+ write(11,*) numpoin,vmin,vmax
+ endif
+
+ if(.not. mask_ibool(iglobval(8))) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglobval(8)) = numpoin
+ write(10,*) numpoin,sngl(xstore(1,NGLLY,NGLLZ,ispec)), &
+ sngl(ystore(1,NGLLY,NGLLZ,ispec)),sngl(zstore(1,NGLLY,NGLLZ,ispec))
+ vmax = 333 !! sqrt((kappavstore(1,NGLLY,NGLLZ,ispec)+4.*muvstore(1,NGLLY,NGLLZ,ispec)/3.)/rhostore(1,NGLLY,NGLLZ,ispec))
+ vmin = 333 !! sqrt(muvstore(1,NGLLY,NGLLZ,ispec)/rhostore(1,NGLLY,NGLLZ,ispec))
+! particular case of the outer core (muvstore contains 1/rho)
+ if(iregion_code == IREGION_OUTER_CORE) then
+ r = dsqrt(xstore(1,NGLLY,NGLLZ,ispec)**2 + ystore(1,NGLLY,NGLLZ,ispec)**2 + zstore(1,NGLLY,NGLLZ,ispec)**2)
+ call prem_display_outer_core(myrank,r,rho,vp,vs,Qkappa,Qmu,idoubling_to_use)
+ vmax = vp
+ vmin = vp
+ endif
+ if(vmin == 0.0) vmin=vmax
+ write(11,*) numpoin,vmin,vmax
+ endif
+
+ if(.not. mask_ibool(iglobval(5))) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglobval(5)) = numpoin
+ write(10,*) numpoin,sngl(xstore(1,1,NGLLZ,ispec)), &
+ sngl(ystore(1,1,NGLLZ,ispec)),sngl(zstore(1,1,NGLLZ,ispec))
+ vmax = 333 !! sqrt((kappavstore(1,1,NGLLZ,ispec)+4.*muvstore(1,1,NGLLZ,ispec)/3.)/rhostore(1,1,NGLLZ,ispec))
+ vmin = 333 !! sqrt(muvstore(1,1,NGLLZ,ispec)/rhostore(1,1,NGLLZ,ispec))
+! particular case of the outer core (muvstore contains 1/rho)
+ if(iregion_code == IREGION_OUTER_CORE) then
+ r = dsqrt(xstore(1,1,NGLLZ,ispec)**2 + ystore(1,1,NGLLZ,ispec)**2 + zstore(1,1,NGLLZ,ispec)**2)
+ call prem_display_outer_core(myrank,r,rho,vp,vs,Qkappa,Qmu,idoubling_to_use)
+ vmax = vp
+ vmin = vp
+ endif
+ if(vmin == 0.0) vmin=vmax
+ write(11,*) numpoin,vmin,vmax
+ endif
+
+ mask_ibool(iglobval(1)) = .true.
+ mask_ibool(iglobval(4)) = .true.
+ mask_ibool(iglobval(8)) = .true.
+ mask_ibool(iglobval(5)) = .true.
+ endif
+
+! face xi = xi_max
+ if(iboun(2,ispec)) then
+
+ if(.not. mask_ibool(iglobval(2))) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglobval(2)) = numpoin
+ write(10,*) numpoin,sngl(xstore(NGLLX,1,1,ispec)), &
+ sngl(ystore(NGLLX,1,1,ispec)),sngl(zstore(NGLLX,1,1,ispec))
+ vmax = 333 !! sqrt((kappavstore(NGLLX,1,1,ispec)+4.*muvstore(NGLLX,1,1,ispec)/3.)/rhostore(NGLLX,1,1,ispec))
+ vmin = 333 !! sqrt(muvstore(NGLLX,1,1,ispec)/rhostore(NGLLX,1,1,ispec))
+! particular case of the outer core (muvstore contains 1/rho)
+ if(iregion_code == IREGION_OUTER_CORE) then
+ r = dsqrt(xstore(NGLLX,1,1,ispec)**2 + ystore(NGLLX,1,1,ispec)**2 + zstore(NGLLX,1,1,ispec)**2)
+ call prem_display_outer_core(myrank,r,rho,vp,vs,Qkappa,Qmu,idoubling_to_use)
+ vmax = vp
+ vmin = vp
+ endif
+ if(vmin == 0.0) vmin=vmax
+ write(11,*) numpoin,vmin,vmax
+ endif
+
+ if(.not. mask_ibool(iglobval(3))) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglobval(3)) = numpoin
+ write(10,*) numpoin,sngl(xstore(NGLLX,NGLLY,1,ispec)), &
+ sngl(ystore(NGLLX,NGLLY,1,ispec)),sngl(zstore(NGLLX,NGLLY,1,ispec))
+ vmax = 333 !! sqrt((kappavstore(NGLLX,NGLLY,1,ispec)+4.*muvstore(NGLLX,NGLLY,1,ispec)/3.)/rhostore(NGLLX,NGLLY,1,ispec))
+ vmin = 333 !! sqrt(muvstore(NGLLX,NGLLY,1,ispec)/rhostore(NGLLX,NGLLY,1,ispec))
+! particular case of the outer core (muvstore contains 1/rho)
+ if(iregion_code == IREGION_OUTER_CORE) then
+ r = dsqrt(xstore(NGLLX,NGLLY,1,ispec)**2 + ystore(NGLLX,NGLLY,1,ispec)**2 + zstore(NGLLX,NGLLY,1,ispec)**2)
+ call prem_display_outer_core(myrank,r,rho,vp,vs,Qkappa,Qmu,idoubling_to_use)
+ vmax = vp
+ vmin = vp
+ endif
+ if(vmin == 0.0) vmin=vmax
+ write(11,*) numpoin,vmin,vmax
+ endif
+
+ if(.not. mask_ibool(iglobval(7))) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglobval(7)) = numpoin
+ write(10,*) numpoin,sngl(xstore(NGLLX,NGLLY,NGLLZ,ispec)), &
+ sngl(ystore(NGLLX,NGLLY,NGLLZ,ispec)),sngl(zstore(NGLLX,NGLLY,NGLLZ,ispec))
+ vmax = 333 !! sqrt((kappavstore(NGLLX,NGLLY,NGLLZ,ispec)+ &
+!! 4.*muvstore(NGLLX,NGLLY,NGLLZ,ispec)/3.)/rhostore(NGLLX,NGLLY,NGLLZ,ispec))
+ vmin = 333 !! sqrt(muvstore(NGLLX,NGLLY,NGLLZ,ispec)/rhostore(NGLLX,NGLLY,NGLLZ,ispec))
+! particular case of the outer core (muvstore contains 1/rho)
+ if(iregion_code == IREGION_OUTER_CORE) then
+ r = dsqrt(xstore(NGLLX,NGLLY,NGLLZ,ispec)**2 + ystore(NGLLX,NGLLY,NGLLZ,ispec)**2 + zstore(NGLLX,NGLLY,NGLLZ,ispec)**2)
+ call prem_display_outer_core(myrank,r,rho,vp,vs,Qkappa,Qmu,idoubling_to_use)
+ vmax = vp
+ vmin = vp
+ endif
+ if(vmin == 0.0) vmin=vmax
+ write(11,*) numpoin,vmin,vmax
+ endif
+
+ if(.not. mask_ibool(iglobval(6))) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglobval(6)) = numpoin
+ write(10,*) numpoin,sngl(xstore(NGLLX,1,NGLLZ,ispec)), &
+ sngl(ystore(NGLLX,1,NGLLZ,ispec)),sngl(zstore(NGLLX,1,NGLLZ,ispec))
+ vmax = 333 !! sqrt((kappavstore(NGLLX,1,NGLLZ,ispec)+4.*muvstore(NGLLX,1,NGLLZ,ispec)/3.)/rhostore(NGLLX,1,NGLLZ,ispec))
+ vmin = 333 !! sqrt(muvstore(NGLLX,1,NGLLZ,ispec)/rhostore(NGLLX,1,NGLLZ,ispec))
+! particular case of the outer core (muvstore contains 1/rho)
+ if(iregion_code == IREGION_OUTER_CORE) then
+ r = dsqrt(xstore(NGLLX,1,NGLLZ,ispec)**2 + ystore(NGLLX,1,NGLLZ,ispec)**2 + zstore(NGLLX,1,NGLLZ,ispec)**2)
+ call prem_display_outer_core(myrank,r,rho,vp,vs,Qkappa,Qmu,idoubling_to_use)
+ vmax = vp
+ vmin = vp
+ endif
+ if(vmin == 0.0) vmin=vmax
+ write(11,*) numpoin,vmin,vmax
+ endif
+
+ mask_ibool(iglobval(2)) = .true.
+ mask_ibool(iglobval(3)) = .true.
+ mask_ibool(iglobval(7)) = .true.
+ mask_ibool(iglobval(6)) = .true.
+ endif
+
+! face eta = eta_min
+ if(iboun(3,ispec)) then
+
+ if(.not. mask_ibool(iglobval(1))) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglobval(1)) = numpoin
+ write(10,*) numpoin,sngl(xstore(1,1,1,ispec)), &
+ sngl(ystore(1,1,1,ispec)),sngl(zstore(1,1,1,ispec))
+ vmax = 333 !! sqrt((kappavstore(1,1,1,ispec)+4.*muvstore(1,1,1,ispec)/3.)/rhostore(1,1,1,ispec))
+ vmin = 333 !! sqrt(muvstore(1,1,1,ispec)/rhostore(1,1,1,ispec))
+! particular case of the outer core (muvstore contains 1/rho)
+ if(iregion_code == IREGION_OUTER_CORE) then
+ r = dsqrt(xstore(1,1,1,ispec)**2 + ystore(1,1,1,ispec)**2 + zstore(1,1,1,ispec)**2)
+ call prem_display_outer_core(myrank,r,rho,vp,vs,Qkappa,Qmu,idoubling_to_use)
+ vmax = vp
+ vmin = vp
+ endif
+ if(vmin == 0.0) vmin=vmax
+ write(11,*) numpoin,vmin,vmax
+ endif
+
+ if(.not. mask_ibool(iglobval(2))) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglobval(2)) = numpoin
+ write(10,*) numpoin,sngl(xstore(NGLLX,1,1,ispec)), &
+ sngl(ystore(NGLLX,1,1,ispec)),sngl(zstore(NGLLX,1,1,ispec))
+ vmax = 333 !! sqrt((kappavstore(NGLLX,1,1,ispec)+4.*muvstore(NGLLX,1,1,ispec)/3.)/rhostore(NGLLX,1,1,ispec))
+ vmin = 333 !! sqrt(muvstore(NGLLX,1,1,ispec)/rhostore(NGLLX,1,1,ispec))
+! particular case of the outer core (muvstore contains 1/rho)
+ if(iregion_code == IREGION_OUTER_CORE) then
+ r = dsqrt(xstore(NGLLX,1,1,ispec)**2 + ystore(NGLLX,1,1,ispec)**2 + zstore(NGLLX,1,1,ispec)**2)
+ call prem_display_outer_core(myrank,r,rho,vp,vs,Qkappa,Qmu,idoubling_to_use)
+ vmax = vp
+ vmin = vp
+ endif
+ if(vmin == 0.0) vmin=vmax
+ write(11,*) numpoin,vmin,vmax
+ endif
+
+ if(.not. mask_ibool(iglobval(6))) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglobval(6)) = numpoin
+ write(10,*) numpoin,sngl(xstore(NGLLX,1,NGLLZ,ispec)), &
+ sngl(ystore(NGLLX,1,NGLLZ,ispec)),sngl(zstore(NGLLX,1,NGLLZ,ispec))
+ vmax = 333 !! sqrt((kappavstore(NGLLX,1,NGLLZ,ispec)+4.*muvstore(NGLLX,1,NGLLZ,ispec)/3.)/rhostore(NGLLX,1,NGLLZ,ispec))
+ vmin = 333 !! sqrt(muvstore(NGLLX,1,NGLLZ,ispec)/rhostore(NGLLX,1,NGLLZ,ispec))
+! particular case of the outer core (muvstore contains 1/rho)
+ if(iregion_code == IREGION_OUTER_CORE) then
+ r = dsqrt(xstore(NGLLX,1,NGLLZ,ispec)**2 + ystore(NGLLX,1,NGLLZ,ispec)**2 + zstore(NGLLX,1,NGLLZ,ispec)**2)
+ call prem_display_outer_core(myrank,r,rho,vp,vs,Qkappa,Qmu,idoubling_to_use)
+ vmax = vp
+ vmin = vp
+ endif
+ if(vmin == 0.0) vmin=vmax
+ write(11,*) numpoin,vmin,vmax
+ endif
+
+ if(.not. mask_ibool(iglobval(5))) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglobval(5)) = numpoin
+ write(10,*) numpoin,sngl(xstore(1,1,NGLLZ,ispec)), &
+ sngl(ystore(1,1,NGLLZ,ispec)),sngl(zstore(1,1,NGLLZ,ispec))
+ vmax = 333 !! sqrt((kappavstore(1,1,NGLLZ,ispec)+4.*muvstore(1,1,NGLLZ,ispec)/3.)/rhostore(1,1,NGLLZ,ispec))
+ vmin = 333 !! sqrt(muvstore(1,1,NGLLZ,ispec)/rhostore(1,1,NGLLZ,ispec))
+! particular case of the outer core (muvstore contains 1/rho)
+ if(iregion_code == IREGION_OUTER_CORE) then
+ r = dsqrt(xstore(1,1,NGLLZ,ispec)**2 + ystore(1,1,NGLLZ,ispec)**2 + zstore(1,1,NGLLZ,ispec)**2)
+ call prem_display_outer_core(myrank,r,rho,vp,vs,Qkappa,Qmu,idoubling_to_use)
+ vmax = vp
+ vmin = vp
+ endif
+ if(vmin == 0.0) vmin=vmax
+ write(11,*) numpoin,vmin,vmax
+ endif
+
+ mask_ibool(iglobval(1)) = .true.
+ mask_ibool(iglobval(2)) = .true.
+ mask_ibool(iglobval(6)) = .true.
+ mask_ibool(iglobval(5)) = .true.
+ endif
+
+! face eta = eta_max
+ if(iboun(4,ispec)) then
+
+ if(.not. mask_ibool(iglobval(4))) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglobval(4)) = numpoin
+ write(10,*) numpoin,sngl(xstore(1,NGLLY,1,ispec)), &
+ sngl(ystore(1,NGLLY,1,ispec)),sngl(zstore(1,NGLLY,1,ispec))
+ vmax = 333 !! sqrt((kappavstore(1,NGLLY,1,ispec)+4.*muvstore(1,NGLLY,1,ispec)/3.)/rhostore(1,NGLLY,1,ispec))
+ vmin = 333 !! sqrt(muvstore(1,NGLLY,1,ispec)/rhostore(1,NGLLY,1,ispec))
+! particular case of the outer core (muvstore contains 1/rho)
+ if(iregion_code == IREGION_OUTER_CORE) then
+ r = dsqrt(xstore(1,NGLLY,1,ispec)**2 + ystore(1,NGLLY,1,ispec)**2 + zstore(1,NGLLY,1,ispec)**2)
+ call prem_display_outer_core(myrank,r,rho,vp,vs,Qkappa,Qmu,idoubling_to_use)
+ vmax = vp
+ vmin = vp
+ endif
+ if(vmin == 0.0) vmin=vmax
+ write(11,*) numpoin,vmin,vmax
+ endif
+
+ if(.not. mask_ibool(iglobval(3))) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglobval(3)) = numpoin
+ write(10,*) numpoin,sngl(xstore(NGLLX,NGLLY,1,ispec)), &
+ sngl(ystore(NGLLX,NGLLY,1,ispec)),sngl(zstore(NGLLX,NGLLY,1,ispec))
+ vmax = 333 !! sqrt((kappavstore(NGLLX,NGLLY,1,ispec)+4.*muvstore(NGLLX,NGLLY,1,ispec)/3.)/rhostore(NGLLX,NGLLY,1,ispec))
+ vmin = 333 !! sqrt(muvstore(NGLLX,NGLLY,1,ispec)/rhostore(NGLLX,NGLLY,1,ispec))
+! particular case of the outer core (muvstore contains 1/rho)
+ if(iregion_code == IREGION_OUTER_CORE) then
+ r = dsqrt(xstore(NGLLX,NGLLY,1,ispec)**2 + ystore(NGLLX,NGLLY,1,ispec)**2 + zstore(NGLLX,NGLLY,1,ispec)**2)
+ call prem_display_outer_core(myrank,r,rho,vp,vs,Qkappa,Qmu,idoubling_to_use)
+ vmax = vp
+ vmin = vp
+ endif
+ if(vmin == 0.0) vmin=vmax
+ write(11,*) numpoin,vmin,vmax
+ endif
+
+ if(.not. mask_ibool(iglobval(7))) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglobval(7)) = numpoin
+ write(10,*) numpoin,sngl(xstore(NGLLX,NGLLY,NGLLZ,ispec)), &
+ sngl(ystore(NGLLX,NGLLY,NGLLZ,ispec)),sngl(zstore(NGLLX,NGLLY,NGLLZ,ispec))
+ vmax = 333 !! sqrt((kappavstore(NGLLX,NGLLY,NGLLZ,ispec)+ &
+!! 4.*muvstore(NGLLX,NGLLY,NGLLZ,ispec)/3.)/rhostore(NGLLX,NGLLY,NGLLZ,ispec))
+ vmin = 333 !! sqrt(muvstore(NGLLX,NGLLY,NGLLZ,ispec)/rhostore(NGLLX,NGLLY,NGLLZ,ispec))
+! particular case of the outer core (muvstore contains 1/rho)
+ if(iregion_code == IREGION_OUTER_CORE) then
+ r = dsqrt(xstore(NGLLX,NGLLY,NGLLZ,ispec)**2 + ystore(NGLLX,NGLLY,NGLLZ,ispec)**2 + zstore(NGLLX,NGLLY,NGLLZ,ispec)**2)
+ call prem_display_outer_core(myrank,r,rho,vp,vs,Qkappa,Qmu,idoubling_to_use)
+ vmax = vp
+ vmin = vp
+ endif
+ if(vmin == 0.0) vmin=vmax
+ write(11,*) numpoin,vmin,vmax
+ endif
+
+ if(.not. mask_ibool(iglobval(8))) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglobval(8)) = numpoin
+ write(10,*) numpoin,sngl(xstore(1,NGLLY,NGLLZ,ispec)), &
+ sngl(ystore(1,NGLLY,NGLLZ,ispec)),sngl(zstore(1,NGLLY,NGLLZ,ispec))
+ vmax = 333 !! sqrt((kappavstore(1,NGLLY,NGLLZ,ispec)+4.*muvstore(1,NGLLY,NGLLZ,ispec)/3.)/rhostore(1,NGLLY,NGLLZ,ispec))
+ vmin = 333 !! sqrt(muvstore(1,NGLLY,NGLLZ,ispec)/rhostore(1,NGLLY,NGLLZ,ispec))
+! particular case of the outer core (muvstore contains 1/rho)
+ if(iregion_code == IREGION_OUTER_CORE) then
+ r = dsqrt(xstore(1,NGLLY,NGLLZ,ispec)**2 + ystore(1,NGLLY,NGLLZ,ispec)**2 + zstore(1,NGLLY,NGLLZ,ispec)**2)
+ call prem_display_outer_core(myrank,r,rho,vp,vs,Qkappa,Qmu,idoubling_to_use)
+ vmax = vp
+ vmin = vp
+ endif
+ if(vmin == 0.0) vmin=vmax
+ write(11,*) numpoin,vmin,vmax
+ endif
+
+ mask_ibool(iglobval(4)) = .true.
+ mask_ibool(iglobval(3)) = .true.
+ mask_ibool(iglobval(7)) = .true.
+ mask_ibool(iglobval(8)) = .true.
+ endif
+
+ endif
+ enddo
+
+! check that number of global points output is okay
+ if(numpoin /= npoin) &
+ call exit_MPI(myrank,'incorrect number of global points in AVS or DX file creation')
+
+ close(10)
+ close(11)
+
+! output global AVS or DX elements
+
+! writing elements
+ open(unit=10,file='OUTPUT_FILES/'//prname(1:len_trim(prname))//'AVS_DXelementschunks.txt',status='unknown')
+ if(ISOTROPIC_3D_MANTLE) &
+ open(unit=11,file='OUTPUT_FILES/'//prname(1:len_trim(prname))//'AVS_DXelementschunks_dvp_dvs.txt',status='unknown')
+
+! number of elements in AVS or DX file
+ write(10,*) nspecface
+
+ ispecface = 0
+ do ispec=1,nspec
+! only if on face
+ if(iboun(1,ispec) .or. iboun(2,ispec) .or. &
+ iboun(3,ispec) .or. iboun(4,ispec)) then
+ iglobval(1)=ibool(1,1,1,ispec)
+ iglobval(2)=ibool(NGLLX,1,1,ispec)
+ iglobval(3)=ibool(NGLLX,NGLLY,1,ispec)
+ iglobval(4)=ibool(1,NGLLY,1,ispec)
+ iglobval(5)=ibool(1,1,NGLLZ,ispec)
+ iglobval(6)=ibool(NGLLX,1,NGLLZ,ispec)
+ iglobval(7)=ibool(NGLLX,NGLLY,NGLLZ,ispec)
+ iglobval(8)=ibool(1,NGLLY,NGLLZ,ispec)
+
+! include lateral variations if needed
+
+ if(ISOTROPIC_3D_MANTLE) then
+! pick a point within the element and get its radius
+ r=dsqrt(xstore(2,2,2,ispec)**2+ystore(2,2,2,ispec)**2+zstore(2,2,2,ispec)**2)
+
+ if(r > RCMB/R_EARTH .and. r < R_UNIT_SPHERE) then
+! average over the element
+ dvp = 0.0
+ dvs = 0.0
+ np =0
+ do k=2,NGLLZ-1
+ do j=2,NGLLY-1
+ do i=2,NGLLX-1
+ np=np+1
+ x=xstore(i,j,k,ispec)
+ y=ystore(i,j,k,ispec)
+ z=zstore(i,j,k,ispec)
+ r=dsqrt(x*x+y*y+z*z)
+! take out ellipticity
+ if(ELLIPTICITY) then
+ call xyz_2_rthetaphi_dble(x,y,z,r,theta,phi_dummy)
+ cost=dcos(theta)
+ p20=0.5d0*(3.0d0*cost*cost-1.0d0)
+ call spline_evaluation(rspl,espl,espl2,nspl,r,ell)
+ factor=ONE-(TWO/3.0d0)*ell*p20
+ r=r/factor
+ endif
+
+ if(iregion_code == IREGION_OUTER_CORE) then
+! doubling array is not created in the outer core to save memory because it is useless
+! therefore impose the doubling flag manually if needed
+ idoubling_to_use = IFLAG_OUTER_CORE_NORMAL
+ else
+ idoubling_to_use = idoubling(ispec)
+ endif
+
+ if(REFERENCE_1D_MODEL == REFERENCE_MODEL_IASP91) then
+ call model_iasp91(myrank,r,rho,vp,vs,Qkappa,Qmu,idoubling_to_use,ONE_CRUST, &
+ .true.,RICB,RCMB,RTOPDDOUBLEPRIME,R771,R670,R400,R220,R120,RMOHO,RMIDDLE_CRUST)
+
+ else if(REFERENCE_1D_MODEL == REFERENCE_MODEL_PREM) then
+ call prem_iso(myrank,r,rho,vp,vs,Qkappa,Qmu,idoubling_to_use, &
+ CRUSTAL,ONE_CRUST,.true.,RICB,RCMB,RTOPDDOUBLEPRIME, &
+ R600,R670,R220,R771,R400,R80,RMOHO,RMIDDLE_CRUST,ROCEAN)
+
+ else if(REFERENCE_1D_MODEL == REFERENCE_MODEL_1066A) then
+ call model_1066a(r,rho,vp,vs,Qkappa,Qmu,idoubling_to_use,M1066a_V)
+
+ else if(REFERENCE_1D_MODEL == REFERENCE_MODEL_AK135) then
+ call model_ak135(r,rho,vp,vs,Qkappa,Qmu,idoubling_to_use,Mak135_V)
+
+ else if(REFERENCE_1D_MODEL == REFERENCE_MODEL_REF) then
+ call model_ref(r,rho,vpv,vph,vsv,vsh,eta_aniso,Qkappa,Qmu,idoubling_to_use,CRUSTAL,Mref_V)
+ vp = vpv
+ vs = vsv
+ else if(REFERENCE_1D_MODEL == REFERENCE_MODEL_JP1D) then
+ call model_jp1d(myrank,r,rho,vp,vs,Qkappa,Qmu,idoubling_to_use, &
+ .true.,RICB,RCMB,RTOPDDOUBLEPRIME, &
+ R670,R220,R771,R400,R80,RMOHO,RMIDDLE_CRUST)
+ else if(REFERENCE_1D_MODEL == REFERENCE_MODEL_SEA1D) then
+ call model_sea1d(r,rho,vp,vs,Qkappa,Qmu,idoubling_to_use,SEA1DM_V)
+ else
+ call exit_MPI(myrank,'unknown 1D reference Earth model in writing of AVS/DX data')
+ endif
+
+ dvp = 333 !! dvp + (sqrt((kappavstore(i,j,k,ispec)+ &
+!! 4.*muvstore(i,j,k,ispec)/3.)/rhostore(i,j,k,ispec)) - sngl(vp))/sngl(vp)
+ dvs = 333 !! dvs + (sqrt(muvstore(i,j,k,ispec)/rhostore(i,j,k,ispec)) - sngl(vs))/sngl(vs)
+ enddo
+ enddo
+ enddo
+ dvp = dvp / np
+ dvs = dvs / np
+ else
+ dvp = 0.0
+ dvs = 0.0
+ endif
+ endif
+
+! face xi = xi_min
+ if(iboun(1,ispec)) then
+ ispecface = ispecface + 1
+ write(10,*) ispecface,idoubling_to_use,num_ibool_AVS_DX(iglobval(1)), &
+ num_ibool_AVS_DX(iglobval(4)),num_ibool_AVS_DX(iglobval(8)), &
+ num_ibool_AVS_DX(iglobval(5))
+ if(ISOTROPIC_3D_MANTLE) write(11,*) ispecface,dvp,dvs
+ endif
+
+! face xi = xi_max
+ if(iboun(2,ispec)) then
+ ispecface = ispecface + 1
+ write(10,*) ispecface,idoubling_to_use,num_ibool_AVS_DX(iglobval(2)), &
+ num_ibool_AVS_DX(iglobval(3)),num_ibool_AVS_DX(iglobval(7)), &
+ num_ibool_AVS_DX(iglobval(6))
+ if(ISOTROPIC_3D_MANTLE) write(11,*) ispecface,dvp,dvs
+ endif
+
+! face eta = eta_min
+ if(iboun(3,ispec)) then
+ ispecface = ispecface + 1
+ write(10,*) ispecface,idoubling_to_use,num_ibool_AVS_DX(iglobval(1)), &
+ num_ibool_AVS_DX(iglobval(2)),num_ibool_AVS_DX(iglobval(6)), &
+ num_ibool_AVS_DX(iglobval(5))
+ if(ISOTROPIC_3D_MANTLE) write(11,*) ispecface,dvp,dvs
+ endif
+
+! face eta = eta_max
+ if(iboun(4,ispec)) then
+ ispecface = ispecface + 1
+ write(10,*) ispecface,idoubling_to_use,num_ibool_AVS_DX(iglobval(4)), &
+ num_ibool_AVS_DX(iglobval(3)),num_ibool_AVS_DX(iglobval(7)), &
+ num_ibool_AVS_DX(iglobval(8))
+ if(ISOTROPIC_3D_MANTLE) write(11,*) ispecface,dvp,dvs
+ endif
+
+ endif
+ enddo
+
+! check that number of surface elements output is okay
+ if(ispecface /= nspecface) &
+ call exit_MPI(myrank,'incorrect number of surface elements in AVS or DX file creation')
+
+ close(10)
+ if(ISOTROPIC_3D_MANTLE) close(11)
+
+ end subroutine write_AVS_DX_global_chunks_data
+
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/write_AVS_DX_global_data.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/write_AVS_DX_global_data.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/write_AVS_DX_global_data.f90 2008-08-15 03:35:19 UTC (rev 12642)
@@ -0,0 +1,205 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 4 . 1
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory, California Institute of Technology, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) California Institute of Technology and University of Pau / CNRS / INRIA
+! August 2008
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+! create AVS or DX 3D data for the slice, to be recombined in postprocessing
+ subroutine write_AVS_DX_global_data(myrank,prname,nspec,ibool,idoubling, &
+ xstore,ystore,zstore,num_ibool_AVS_DX,mask_ibool,npointot,iregion_code)
+
+ implicit none
+
+ include "constants.h"
+
+ integer nspec,myrank
+ integer ibool(NGLLX,NGLLY,NGLLZ,nspec)
+
+ integer idoubling(nspec)
+
+ double precision xstore(NGLLX,NGLLY,NGLLZ,nspec)
+ double precision ystore(NGLLX,NGLLY,NGLLZ,nspec)
+ double precision zstore(NGLLX,NGLLY,NGLLZ,nspec)
+
+! logical mask used to output global points only once
+ integer npointot,iregion_code
+ logical mask_ibool(npointot)
+
+! numbering of global AVS or DX points
+ integer num_ibool_AVS_DX(npointot)
+
+ integer ispec,idoubling_to_use
+ integer iglob1,iglob2,iglob3,iglob4,iglob5,iglob6,iglob7,iglob8
+ integer npoin,numpoin
+
+! processor identification
+ character(len=150) prname
+
+! writing points
+ open(unit=10,file='OUTPUT_FILES/'//prname(1:len_trim(prname))//'AVS_DXpoints.txt',status='unknown')
+
+! erase the logical mask used to mark points already found
+ mask_ibool(:) = .false.
+
+! mark global AVS or DX points
+ do ispec=1,nspec
+ iglob1=ibool(1,1,1,ispec)
+ iglob2=ibool(NGLLX,1,1,ispec)
+ iglob3=ibool(NGLLX,NGLLY,1,ispec)
+ iglob4=ibool(1,NGLLY,1,ispec)
+ iglob5=ibool(1,1,NGLLZ,ispec)
+ iglob6=ibool(NGLLX,1,NGLLZ,ispec)
+ iglob7=ibool(NGLLX,NGLLY,NGLLZ,ispec)
+ iglob8=ibool(1,NGLLY,NGLLZ,ispec)
+ mask_ibool(iglob1) = .true.
+ mask_ibool(iglob2) = .true.
+ mask_ibool(iglob3) = .true.
+ mask_ibool(iglob4) = .true.
+ mask_ibool(iglob5) = .true.
+ mask_ibool(iglob6) = .true.
+ mask_ibool(iglob7) = .true.
+ mask_ibool(iglob8) = .true.
+ enddo
+
+! count global number of AVS or DX points
+ npoin = count(mask_ibool(:))
+
+! number of points in AVS or DX file
+ write(10,*) npoin
+
+! erase the logical mask used to mark points already found
+ mask_ibool(:) = .false.
+
+! output global AVS or DX points
+ numpoin = 0
+ do ispec=1,nspec
+ iglob1=ibool(1,1,1,ispec)
+ iglob2=ibool(NGLLX,1,1,ispec)
+ iglob3=ibool(NGLLX,NGLLY,1,ispec)
+ iglob4=ibool(1,NGLLY,1,ispec)
+ iglob5=ibool(1,1,NGLLZ,ispec)
+ iglob6=ibool(NGLLX,1,NGLLZ,ispec)
+ iglob7=ibool(NGLLX,NGLLY,NGLLZ,ispec)
+ iglob8=ibool(1,NGLLY,NGLLZ,ispec)
+ if(.not. mask_ibool(iglob1)) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglob1) = numpoin
+ write(10,*) numpoin,sngl(xstore(1,1,1,ispec)), &
+ sngl(ystore(1,1,1,ispec)),sngl(zstore(1,1,1,ispec))
+ endif
+ if(.not. mask_ibool(iglob2)) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglob2) = numpoin
+ write(10,*) numpoin,sngl(xstore(NGLLX,1,1,ispec)), &
+ sngl(ystore(NGLLX,1,1,ispec)),sngl(zstore(NGLLX,1,1,ispec))
+ endif
+ if(.not. mask_ibool(iglob3)) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglob3) = numpoin
+ write(10,*) numpoin,sngl(xstore(NGLLX,NGLLY,1,ispec)), &
+ sngl(ystore(NGLLX,NGLLY,1,ispec)),sngl(zstore(NGLLX,NGLLY,1,ispec))
+ endif
+ if(.not. mask_ibool(iglob4)) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglob4) = numpoin
+ write(10,*) numpoin,sngl(xstore(1,NGLLY,1,ispec)), &
+ sngl(ystore(1,NGLLY,1,ispec)),sngl(zstore(1,NGLLY,1,ispec))
+ endif
+ if(.not. mask_ibool(iglob5)) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglob5) = numpoin
+ write(10,*) numpoin,sngl(xstore(1,1,NGLLZ,ispec)), &
+ sngl(ystore(1,1,NGLLZ,ispec)),sngl(zstore(1,1,NGLLZ,ispec))
+ endif
+ if(.not. mask_ibool(iglob6)) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglob6) = numpoin
+ write(10,*) numpoin,sngl(xstore(NGLLX,1,NGLLZ,ispec)), &
+ sngl(ystore(NGLLX,1,NGLLZ,ispec)),sngl(zstore(NGLLX,1,NGLLZ,ispec))
+ endif
+ if(.not. mask_ibool(iglob7)) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglob7) = numpoin
+ write(10,*) numpoin,sngl(xstore(NGLLX,NGLLY,NGLLZ,ispec)), &
+ sngl(ystore(NGLLX,NGLLY,NGLLZ,ispec)),sngl(zstore(NGLLX,NGLLY,NGLLZ,ispec))
+ endif
+ if(.not. mask_ibool(iglob8)) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglob8) = numpoin
+ write(10,*) numpoin,sngl(xstore(1,NGLLY,NGLLZ,ispec)), &
+ sngl(ystore(1,NGLLY,NGLLZ,ispec)),sngl(zstore(1,NGLLY,NGLLZ,ispec))
+ endif
+ mask_ibool(iglob1) = .true.
+ mask_ibool(iglob2) = .true.
+ mask_ibool(iglob3) = .true.
+ mask_ibool(iglob4) = .true.
+ mask_ibool(iglob5) = .true.
+ mask_ibool(iglob6) = .true.
+ mask_ibool(iglob7) = .true.
+ mask_ibool(iglob8) = .true.
+ enddo
+
+! check that number of global points output is okay
+ if(numpoin /= npoin) &
+ call exit_MPI(myrank,'incorrect number of global points in AVS or DX file creation')
+
+ close(10)
+
+! writing elements
+ open(unit=10,file='OUTPUT_FILES/'//prname(1:len_trim(prname))//'AVS_DXelements.txt',status='unknown')
+
+! number of elements in AVS or DX file
+ write(10,*) nspec
+
+! output global AVS or DX elements
+ do ispec = 1,nspec
+
+ iglob1=ibool(1,1,1,ispec)
+ iglob2=ibool(NGLLX,1,1,ispec)
+ iglob3=ibool(NGLLX,NGLLY,1,ispec)
+ iglob4=ibool(1,NGLLY,1,ispec)
+ iglob5=ibool(1,1,NGLLZ,ispec)
+ iglob6=ibool(NGLLX,1,NGLLZ,ispec)
+ iglob7=ibool(NGLLX,NGLLY,NGLLZ,ispec)
+ iglob8=ibool(1,NGLLY,NGLLZ,ispec)
+
+ if(iregion_code == IREGION_OUTER_CORE) then
+! doubling array is not created in the outer core to save memory because it is useless
+! therefore impose the doubling flag manually if needed
+ idoubling_to_use = IFLAG_OUTER_CORE_NORMAL
+ else
+ idoubling_to_use = idoubling(ispec)
+ endif
+
+ write(10,*) ispec,idoubling_to_use,num_ibool_AVS_DX(iglob1), &
+ num_ibool_AVS_DX(iglob2),num_ibool_AVS_DX(iglob3), &
+ num_ibool_AVS_DX(iglob4),num_ibool_AVS_DX(iglob5), &
+ num_ibool_AVS_DX(iglob6),num_ibool_AVS_DX(iglob7), &
+ num_ibool_AVS_DX(iglob8)
+ enddo
+
+ close(10)
+
+ end subroutine write_AVS_DX_global_data
+
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/write_AVS_DX_global_faces_data.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/write_AVS_DX_global_faces_data.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/write_AVS_DX_global_faces_data.f90 2008-08-15 03:35:19 UTC (rev 12642)
@@ -0,0 +1,362 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 4 . 1
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory, California Institute of Technology, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) California Institute of Technology and University of Pau / CNRS / INRIA
+! August 2008
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+! create AVS or DX 2D data for the faces of the slice,
+! to be recombined in postprocessing
+
+ subroutine write_AVS_DX_global_faces_data(myrank,prname,nspec,iMPIcut_xi,iMPIcut_eta, &
+ ibool,idoubling,xstore,ystore,zstore,num_ibool_AVS_DX,mask_ibool, &
+ npointot,iregion_code)
+
+ implicit none
+
+ include "constants.h"
+
+ integer nspec,myrank,iregion_code
+ integer ibool(NGLLX,NGLLY,NGLLZ,nspec)
+
+ integer idoubling(nspec)
+
+ logical iMPIcut_xi(2,nspec)
+ logical iMPIcut_eta(2,nspec)
+
+ double precision xstore(NGLLX,NGLLY,NGLLZ,nspec)
+ double precision ystore(NGLLX,NGLLY,NGLLZ,nspec)
+ double precision zstore(NGLLX,NGLLY,NGLLZ,nspec)
+
+! logical mask used to output global points only once
+ integer npointot
+ logical mask_ibool(npointot)
+
+! numbering of global AVS or DX points
+ integer num_ibool_AVS_DX(npointot)
+
+ integer ispec,idoubling_to_use
+ integer iglob1,iglob2,iglob3,iglob4,iglob5,iglob6,iglob7,iglob8
+ integer npoin,numpoin,nspecface,ispecface
+
+! processor identification
+ character(len=150) prname
+
+! writing points
+ open(unit=10,file='OUTPUT_FILES/'//prname(1:len_trim(prname))//'AVS_DXpointsfaces.txt',status='unknown')
+
+! erase the logical mask used to mark points already found
+ mask_ibool(:) = .false.
+
+ nspecface = 0
+
+! mark global AVS or DX points
+ do ispec=1,nspec
+! only if on face
+ if(iMPIcut_xi(1,ispec) .or. iMPIcut_xi(2,ispec) .or. &
+ iMPIcut_eta(1,ispec) .or. iMPIcut_eta(2,ispec)) then
+ iglob1=ibool(1,1,1,ispec)
+ iglob2=ibool(NGLLX,1,1,ispec)
+ iglob3=ibool(NGLLX,NGLLY,1,ispec)
+ iglob4=ibool(1,NGLLY,1,ispec)
+ iglob5=ibool(1,1,NGLLZ,ispec)
+ iglob6=ibool(NGLLX,1,NGLLZ,ispec)
+ iglob7=ibool(NGLLX,NGLLY,NGLLZ,ispec)
+ iglob8=ibool(1,NGLLY,NGLLZ,ispec)
+
+! face xi = xi_min
+ if(iMPIcut_xi(1,ispec)) then
+ nspecface = nspecface + 1
+ mask_ibool(iglob1) = .true.
+ mask_ibool(iglob4) = .true.
+ mask_ibool(iglob8) = .true.
+ mask_ibool(iglob5) = .true.
+ endif
+
+! face xi = xi_max
+ if(iMPIcut_xi(2,ispec)) then
+ nspecface = nspecface + 1
+ mask_ibool(iglob2) = .true.
+ mask_ibool(iglob3) = .true.
+ mask_ibool(iglob7) = .true.
+ mask_ibool(iglob6) = .true.
+ endif
+
+! face eta = eta_min
+ if(iMPIcut_eta(1,ispec)) then
+ nspecface = nspecface + 1
+ mask_ibool(iglob1) = .true.
+ mask_ibool(iglob2) = .true.
+ mask_ibool(iglob6) = .true.
+ mask_ibool(iglob5) = .true.
+ endif
+
+! face eta = eta_max
+ if(iMPIcut_eta(2,ispec)) then
+ nspecface = nspecface + 1
+ mask_ibool(iglob4) = .true.
+ mask_ibool(iglob3) = .true.
+ mask_ibool(iglob7) = .true.
+ mask_ibool(iglob8) = .true.
+ endif
+
+ endif
+ enddo
+
+! count global number of AVS or DX points
+ npoin = count(mask_ibool(:))
+
+! number of points in AVS or DX file
+ write(10,*) npoin
+
+! erase the logical mask used to mark points already found
+ mask_ibool(:) = .false.
+
+! output global AVS or DX points
+ numpoin = 0
+ do ispec=1,nspec
+! only if on face
+ if(iMPIcut_xi(1,ispec) .or. iMPIcut_xi(2,ispec) .or. &
+ iMPIcut_eta(1,ispec) .or. iMPIcut_eta(2,ispec)) then
+ iglob1=ibool(1,1,1,ispec)
+ iglob2=ibool(NGLLX,1,1,ispec)
+ iglob3=ibool(NGLLX,NGLLY,1,ispec)
+ iglob4=ibool(1,NGLLY,1,ispec)
+ iglob5=ibool(1,1,NGLLZ,ispec)
+ iglob6=ibool(NGLLX,1,NGLLZ,ispec)
+ iglob7=ibool(NGLLX,NGLLY,NGLLZ,ispec)
+ iglob8=ibool(1,NGLLY,NGLLZ,ispec)
+
+! face xi = xi_min
+ if(iMPIcut_xi(1,ispec)) then
+ if(.not. mask_ibool(iglob1)) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglob1) = numpoin
+ write(10,*) numpoin,sngl(xstore(1,1,1,ispec)), &
+ sngl(ystore(1,1,1,ispec)),sngl(zstore(1,1,1,ispec))
+ endif
+ if(.not. mask_ibool(iglob4)) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglob4) = numpoin
+ write(10,*) numpoin,sngl(xstore(1,NGLLY,1,ispec)), &
+ sngl(ystore(1,NGLLY,1,ispec)),sngl(zstore(1,NGLLY,1,ispec))
+ endif
+ if(.not. mask_ibool(iglob8)) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglob8) = numpoin
+ write(10,*) numpoin,sngl(xstore(1,NGLLY,NGLLZ,ispec)), &
+ sngl(ystore(1,NGLLY,NGLLZ,ispec)),sngl(zstore(1,NGLLY,NGLLZ,ispec))
+ endif
+ if(.not. mask_ibool(iglob5)) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglob5) = numpoin
+ write(10,*) numpoin,sngl(xstore(1,1,NGLLZ,ispec)), &
+ sngl(ystore(1,1,NGLLZ,ispec)),sngl(zstore(1,1,NGLLZ,ispec))
+ endif
+ mask_ibool(iglob1) = .true.
+ mask_ibool(iglob4) = .true.
+ mask_ibool(iglob8) = .true.
+ mask_ibool(iglob5) = .true.
+ endif
+
+! face xi = xi_max
+ if(iMPIcut_xi(2,ispec)) then
+ if(.not. mask_ibool(iglob2)) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglob2) = numpoin
+ write(10,*) numpoin,sngl(xstore(NGLLX,1,1,ispec)), &
+ sngl(ystore(NGLLX,1,1,ispec)),sngl(zstore(NGLLX,1,1,ispec))
+ endif
+ if(.not. mask_ibool(iglob3)) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglob3) = numpoin
+ write(10,*) numpoin,sngl(xstore(NGLLX,NGLLY,1,ispec)), &
+ sngl(ystore(NGLLX,NGLLY,1,ispec)),sngl(zstore(NGLLX,NGLLY,1,ispec))
+ endif
+ if(.not. mask_ibool(iglob7)) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglob7) = numpoin
+ write(10,*) numpoin,sngl(xstore(NGLLX,NGLLY,NGLLZ,ispec)), &
+ sngl(ystore(NGLLX,NGLLY,NGLLZ,ispec)),sngl(zstore(NGLLX,NGLLY,NGLLZ,ispec))
+ endif
+ if(.not. mask_ibool(iglob6)) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglob6) = numpoin
+ write(10,*) numpoin,sngl(xstore(NGLLX,1,NGLLZ,ispec)), &
+ sngl(ystore(NGLLX,1,NGLLZ,ispec)),sngl(zstore(NGLLX,1,NGLLZ,ispec))
+ endif
+ mask_ibool(iglob2) = .true.
+ mask_ibool(iglob3) = .true.
+ mask_ibool(iglob7) = .true.
+ mask_ibool(iglob6) = .true.
+ endif
+
+! face eta = eta_min
+ if(iMPIcut_eta(1,ispec)) then
+ if(.not. mask_ibool(iglob1)) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglob1) = numpoin
+ write(10,*) numpoin,sngl(xstore(1,1,1,ispec)), &
+ sngl(ystore(1,1,1,ispec)),sngl(zstore(1,1,1,ispec))
+ endif
+ if(.not. mask_ibool(iglob2)) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglob2) = numpoin
+ write(10,*) numpoin,sngl(xstore(NGLLX,1,1,ispec)), &
+ sngl(ystore(NGLLX,1,1,ispec)),sngl(zstore(NGLLX,1,1,ispec))
+ endif
+ if(.not. mask_ibool(iglob6)) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglob6) = numpoin
+ write(10,*) numpoin,sngl(xstore(NGLLX,1,NGLLZ,ispec)), &
+ sngl(ystore(NGLLX,1,NGLLZ,ispec)),sngl(zstore(NGLLX,1,NGLLZ,ispec))
+ endif
+ if(.not. mask_ibool(iglob5)) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglob5) = numpoin
+ write(10,*) numpoin,sngl(xstore(1,1,NGLLZ,ispec)), &
+ sngl(ystore(1,1,NGLLZ,ispec)),sngl(zstore(1,1,NGLLZ,ispec))
+ endif
+ mask_ibool(iglob1) = .true.
+ mask_ibool(iglob2) = .true.
+ mask_ibool(iglob6) = .true.
+ mask_ibool(iglob5) = .true.
+ endif
+
+! face eta = eta_max
+ if(iMPIcut_eta(2,ispec)) then
+ if(.not. mask_ibool(iglob4)) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglob4) = numpoin
+ write(10,*) numpoin,sngl(xstore(1,NGLLY,1,ispec)), &
+ sngl(ystore(1,NGLLY,1,ispec)),sngl(zstore(1,NGLLY,1,ispec))
+ endif
+ if(.not. mask_ibool(iglob3)) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglob3) = numpoin
+ write(10,*) numpoin,sngl(xstore(NGLLX,NGLLY,1,ispec)), &
+ sngl(ystore(NGLLX,NGLLY,1,ispec)),sngl(zstore(NGLLX,NGLLY,1,ispec))
+ endif
+ if(.not. mask_ibool(iglob7)) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglob7) = numpoin
+ write(10,*) numpoin,sngl(xstore(NGLLX,NGLLY,NGLLZ,ispec)), &
+ sngl(ystore(NGLLX,NGLLY,NGLLZ,ispec)),sngl(zstore(NGLLX,NGLLY,NGLLZ,ispec))
+ endif
+ if(.not. mask_ibool(iglob8)) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglob8) = numpoin
+ write(10,*) numpoin,sngl(xstore(1,NGLLY,NGLLZ,ispec)), &
+ sngl(ystore(1,NGLLY,NGLLZ,ispec)),sngl(zstore(1,NGLLY,NGLLZ,ispec))
+ endif
+ mask_ibool(iglob4) = .true.
+ mask_ibool(iglob3) = .true.
+ mask_ibool(iglob7) = .true.
+ mask_ibool(iglob8) = .true.
+ endif
+
+ endif
+ enddo
+
+! check that number of global points output is okay
+ if(numpoin /= npoin) &
+ call exit_MPI(myrank,'incorrect number of global points in AVS or DX file creation')
+
+ close(10)
+
+! output global AVS or DX elements
+
+! writing elements
+ open(unit=10,file='OUTPUT_FILES/'//prname(1:len_trim(prname))//'AVS_DXelementsfaces.txt',status='unknown')
+
+! number of elements in AVS or DX file
+ write(10,*) nspecface
+
+ ispecface = 0
+
+ do ispec = 1,nspec
+
+! only if on face
+ if(iMPIcut_xi(1,ispec) .or. iMPIcut_xi(2,ispec) .or. &
+ iMPIcut_eta(1,ispec) .or. iMPIcut_eta(2,ispec)) then
+
+ iglob1=ibool(1,1,1,ispec)
+ iglob2=ibool(NGLLX,1,1,ispec)
+ iglob3=ibool(NGLLX,NGLLY,1,ispec)
+ iglob4=ibool(1,NGLLY,1,ispec)
+ iglob5=ibool(1,1,NGLLZ,ispec)
+ iglob6=ibool(NGLLX,1,NGLLZ,ispec)
+ iglob7=ibool(NGLLX,NGLLY,NGLLZ,ispec)
+ iglob8=ibool(1,NGLLY,NGLLZ,ispec)
+
+ if(iregion_code == IREGION_OUTER_CORE) then
+! doubling array is not created in the outer core to save memory because it is useless
+! therefore impose the doubling flag manually if needed
+ idoubling_to_use = IFLAG_OUTER_CORE_NORMAL
+ else
+ idoubling_to_use = idoubling(ispec)
+ endif
+
+! face xi = xi_min
+ if(iMPIcut_xi(1,ispec)) then
+ ispecface = ispecface + 1
+ write(10,*) ispecface,idoubling_to_use,num_ibool_AVS_DX(iglob1), &
+ num_ibool_AVS_DX(iglob4),num_ibool_AVS_DX(iglob8), &
+ num_ibool_AVS_DX(iglob5)
+ endif
+
+! face xi = xi_max
+ if(iMPIcut_xi(2,ispec)) then
+ ispecface = ispecface + 1
+ write(10,*) ispecface,idoubling_to_use,num_ibool_AVS_DX(iglob2), &
+ num_ibool_AVS_DX(iglob3),num_ibool_AVS_DX(iglob7), &
+ num_ibool_AVS_DX(iglob6)
+ endif
+
+! face eta = eta_min
+ if(iMPIcut_eta(1,ispec)) then
+ ispecface = ispecface + 1
+ write(10,*) ispecface,idoubling_to_use,num_ibool_AVS_DX(iglob1), &
+ num_ibool_AVS_DX(iglob2),num_ibool_AVS_DX(iglob6), &
+ num_ibool_AVS_DX(iglob5)
+ endif
+
+! face eta = eta_max
+ if(iMPIcut_eta(2,ispec)) then
+ ispecface = ispecface + 1
+ write(10,*) ispecface,idoubling_to_use,num_ibool_AVS_DX(iglob4), &
+ num_ibool_AVS_DX(iglob3),num_ibool_AVS_DX(iglob7), &
+ num_ibool_AVS_DX(iglob8)
+ endif
+
+ endif
+ enddo
+
+! check that number of surface elements output is okay
+ if(ispecface /= nspecface) &
+ call exit_MPI(myrank,'incorrect number of surface elements in AVS or DX file creation')
+
+ close(10)
+
+ end subroutine write_AVS_DX_global_faces_data
+
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/write_AVS_DX_surface_data.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/write_AVS_DX_surface_data.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/write_AVS_DX_surface_data.f90 2008-08-15 03:35:19 UTC (rev 12642)
@@ -0,0 +1,201 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 4 . 1
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory, California Institute of Technology, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) California Institute of Technology and University of Pau / CNRS / INRIA
+! August 2008
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+! create AVS or DX 2D data for the surface of the model
+! to be recombined in postprocessing
+ subroutine write_AVS_DX_surface_data(myrank,prname,nspec,iboun, &
+ ibool,idoubling,xstore,ystore,zstore,num_ibool_AVS_DX,mask_ibool,npointot,iregion_code)
+
+ implicit none
+
+ include "constants.h"
+
+ integer nspec,myrank,iregion_code
+ integer ibool(NGLLX,NGLLY,NGLLZ,nspec)
+
+ integer idoubling(nspec)
+
+ logical iboun(6,nspec)
+
+ double precision xstore(NGLLX,NGLLY,NGLLZ,nspec)
+ double precision ystore(NGLLX,NGLLY,NGLLZ,nspec)
+ double precision zstore(NGLLX,NGLLY,NGLLZ,nspec)
+
+! logical mask used to output global points only once
+ integer npointot
+ logical mask_ibool(npointot)
+
+! numbering of global AVS or DX points
+ integer num_ibool_AVS_DX(npointot)
+
+ integer ispec,idoubling_to_use
+ integer, dimension(8) :: iglobval
+ integer npoin,numpoin,nspecface,ispecface
+
+! processor identification
+ character(len=150) prname
+
+! writing points
+ open(unit=10,file='OUTPUT_FILES/'//prname(1:len_trim(prname))//'AVS_DXpointssurface.txt',status='unknown')
+
+! erase the logical mask used to mark points already found
+ mask_ibool(:) = .false.
+
+ nspecface = 0
+
+! mark global AVS or DX points
+ do ispec=1,nspec
+! only if at the surface (top plane)
+ if(iboun(6,ispec)) then
+
+ iglobval(5)=ibool(1,1,NGLLZ,ispec)
+ iglobval(6)=ibool(NGLLX,1,NGLLZ,ispec)
+ iglobval(7)=ibool(NGLLX,NGLLY,NGLLZ,ispec)
+ iglobval(8)=ibool(1,NGLLY,NGLLZ,ispec)
+
+! element is at the surface
+ nspecface = nspecface + 1
+ mask_ibool(iglobval(5)) = .true.
+ mask_ibool(iglobval(6)) = .true.
+ mask_ibool(iglobval(7)) = .true.
+ mask_ibool(iglobval(8)) = .true.
+
+ endif
+ enddo
+
+! count global number of AVS or DX points
+ npoin = count(mask_ibool(:))
+
+! number of points in AVS or DX file
+ write(10,*) npoin
+
+! erase the logical mask used to mark points already found
+ mask_ibool(:) = .false.
+
+! output global AVS or DX points
+ numpoin = 0
+ do ispec=1,nspec
+! only if at the surface
+ if(iboun(6,ispec)) then
+
+ iglobval(5)=ibool(1,1,NGLLZ,ispec)
+ iglobval(6)=ibool(NGLLX,1,NGLLZ,ispec)
+ iglobval(7)=ibool(NGLLX,NGLLY,NGLLZ,ispec)
+ iglobval(8)=ibool(1,NGLLY,NGLLZ,ispec)
+
+! top face
+ if(iboun(6,ispec)) then
+
+ if(.not. mask_ibool(iglobval(5))) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglobval(5)) = numpoin
+ write(10,*) numpoin,sngl(xstore(1,1,NGLLZ,ispec)), &
+ sngl(ystore(1,1,NGLLZ,ispec)),sngl(zstore(1,1,NGLLZ,ispec))
+ endif
+
+ if(.not. mask_ibool(iglobval(6))) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglobval(6)) = numpoin
+ write(10,*) numpoin,sngl(xstore(NGLLX,1,NGLLZ,ispec)), &
+ sngl(ystore(NGLLX,1,NGLLZ,ispec)),sngl(zstore(NGLLX,1,NGLLZ,ispec))
+ endif
+
+ if(.not. mask_ibool(iglobval(7))) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglobval(7)) = numpoin
+ write(10,*) numpoin,sngl(xstore(NGLLX,NGLLY,NGLLZ,ispec)), &
+ sngl(ystore(NGLLX,NGLLY,NGLLZ,ispec)),sngl(zstore(NGLLX,NGLLY,NGLLZ,ispec))
+ endif
+
+ if(.not. mask_ibool(iglobval(8))) then
+ numpoin = numpoin + 1
+ num_ibool_AVS_DX(iglobval(8)) = numpoin
+ write(10,*) numpoin,sngl(xstore(1,NGLLY,NGLLZ,ispec)), &
+ sngl(ystore(1,NGLLY,NGLLZ,ispec)),sngl(zstore(1,NGLLY,NGLLZ,ispec))
+ endif
+
+ mask_ibool(iglobval(5)) = .true.
+ mask_ibool(iglobval(6)) = .true.
+ mask_ibool(iglobval(7)) = .true.
+ mask_ibool(iglobval(8)) = .true.
+
+ endif
+
+ endif
+ enddo
+
+! check that number of global points output is okay
+ if(numpoin /= npoin) &
+ call exit_MPI(myrank,'incorrect number of global points in AVS or DX file creation')
+
+ close(10)
+
+! output global AVS or DX elements
+
+! writing elements
+ open(unit=10,file='OUTPUT_FILES/'//prname(1:len_trim(prname))//'AVS_DXelementssurface.txt',status='unknown')
+
+! number of elements in AVS or DX file
+ write(10,*) nspecface
+
+ ispecface = 0
+
+ do ispec = 1,nspec
+
+! only if at the surface
+ if(iboun(6,ispec)) then
+
+ iglobval(5)=ibool(1,1,NGLLZ,ispec)
+ iglobval(6)=ibool(NGLLX,1,NGLLZ,ispec)
+ iglobval(7)=ibool(NGLLX,NGLLY,NGLLZ,ispec)
+ iglobval(8)=ibool(1,NGLLY,NGLLZ,ispec)
+
+ if(iregion_code == IREGION_OUTER_CORE) then
+! doubling array is not created in the outer core to save memory because it is useless
+! therefore impose the doubling flag manually if needed
+ idoubling_to_use = IFLAG_OUTER_CORE_NORMAL
+ else
+ idoubling_to_use = idoubling(ispec)
+ endif
+
+! top face
+ ispecface = ispecface + 1
+ write(10,*) ispecface,idoubling_to_use,num_ibool_AVS_DX(iglobval(5)), &
+ num_ibool_AVS_DX(iglobval(6)),num_ibool_AVS_DX(iglobval(7)), &
+ num_ibool_AVS_DX(iglobval(8))
+
+ endif
+ enddo
+
+! check that number of surface elements output is okay
+ if(ispecface /= nspecface) &
+ call exit_MPI(myrank,'incorrect number of surface elements in AVS or DX file creation')
+
+ close(10)
+
+ end subroutine write_AVS_DX_surface_data
+
More information about the cig-commits
mailing list