[cig-commits] r19323 - in seismo/2D/SPECFEM2D/trunk: EXAMPLES EXAMPLES/Rayleigh_wave EXAMPLES/attenuation setup src/specfem2D

dkomati1 at geodynamics.org dkomati1 at geodynamics.org
Wed Dec 28 06:36:38 PST 2011


Author: dkomati1
Date: 2011-12-28 06:36:37 -0800 (Wed, 28 Dec 2011)
New Revision: 19323

Added:
   seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave/
   seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave/Par_file_Rayleigh_2D
   seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave/README
   seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave/SOURCE_Rayleigh_2D
   seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave/interfaces_Rayleigh_flat.dat
   seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave/process.sh
Modified:
   seismo/2D/SPECFEM2D/trunk/EXAMPLES/attenuation/README
   seismo/2D/SPECFEM2D/trunk/setup/constants.h.in
   seismo/2D/SPECFEM2D/trunk/src/specfem2D/paco_beyond_critical.f90
   seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
Log:
added EXAMPLES/Rayleigh_wave and fixed a small problem in the horizontal offset of Rayleigh sources


Added: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave/Par_file_Rayleigh_2D
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave/Par_file_Rayleigh_2D	                        (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave/Par_file_Rayleigh_2D	2011-12-28 14:36:37 UTC (rev 19323)
@@ -0,0 +1,133 @@
+# title of job
+title                           = Rayleigh wave in a homogeneous 2D medium
+
+# forward or adjoint simulation
+SIMULATION_TYPE                 = 1   # 1 = forward, 2 = adjoint + kernels
+NOISE_TOMOGRAPHY                = 0   # 0 = earthquake simulation, 1/2/3 = noise simulation
+SAVE_FORWARD                    = .false.  # save the last frame, needed for adjoint simulation
+
+# parameters concerning partitioning
+nproc                           = 1              # number of processes
+partitioning_method             = 3              # SCOTCH = 3, ascending order (very bad idea) = 1
+PERFORM_CUTHILL_MCKEE           = .true.         # perform inverse Cuthill-McKee (1969) optimization/permutation for mesh numbering
+
+ngnod                           = 4              # number of control nodes per element (4 or 9)
+initialfield                    = .true.         # use a plane wave as source or not
+add_Bielak_conditions           = .true.         # add Bielak conditions or not if initial plane wave
+assign_external_model           = .false.        # define external earth model or not
+READ_EXTERNAL_SEP_FILE          = .false.        # Read external SEP file from DATA/model_velocity.dat_input, or use routine
+TURN_ATTENUATION_ON             = .false.        # turn attenuation on or off for solid medium
+TURN_VISCATTENUATION_ON         = .false.        # turn viscous attenuation on or off
+Q0                              =  1             # quality factor for viscous attenuation
+freq0                           =  10            # frequency for viscous attenuation
+p_sv                            = .true.         # set the type of calculation (P-SV or SH/membrane waves)
+
+# time step parameters
+nt                              = 5000           # total number of time steps
+deltat                          = 3.90625e-3     # duration of a time step
+USER_T0                         = 0.0d0          # use this t0 as earliest starting time rather than the automatically calculated one
+
+# source parameters
+NSOURCES                        = 1              # number of sources [source info read in CMTSOLUTION file]
+force_normal_to_surface         = .false.        # angleforce normal to surface (external mesh and curve file needed)
+
+# constants for attenuation
+N_SLS                           = 2                      # number of standard linear solids for attenuation
+f0_attenuation                  = 5.196152422706633      # (Hz) relevant only if source is a Dirac or a Heaviside, else it is f0
+
+# receiver line parameters for seismograms
+seismotype                      = 1              # record 1=displ 2=veloc 3=accel 4=pressure
+generate_STATIONS               = .true.         # creates a STATION file in ./DATA
+nreceiverlines                  = 1              # number of receiver lines
+anglerec                        = 0.d0           # angle to rotate components at receivers
+rec_normal_to_surface           = .false.        # base anglerec normal to surface (external mesh and curve file needed)
+SU_FORMAT                       = .false.        # output seismograms in Seismic Unix format (adjoint traces will be read in the same format)
+
+# first receiver line
+nrec                            = 5              # number of receivers
+xdeb                            = 14.5           # first receiver x in meters
+zdeb                            = 9.             # first receiver z in meters
+xfin                            = 17.5           # last receiver x in meters (ignored if onlyone receiver)
+zfin                            = 9.             # last receiver z in meters (ignored if onlyone receiver)
+enreg_surf_same_vertical        = .false.        # receivers inside the medium or at the surface
+
+# display parameters
+NTSTEP_BETWEEN_OUTPUT_INFO      = 200            # display frequency in time steps
+output_postscript_snapshot      = .true.         # output Postscript snapshot of the results
+output_color_image              = .true.         # output color image of the results
+imagetype                       = 1              # display 1=displ 2=veloc 3=accel 4=pressure
+cutsnaps                        = 1.             # minimum amplitude in % for snapshots
+meshvect                        = .true.         # display mesh on vector plots or not
+modelvect                       = .false.        # display velocity model on vector plots
+boundvect                       = .true.         # display boundary conditions on plots
+interpol                        = .true.         # interpolation of the display or not
+pointsdisp                      = 6              # points for interpolation of display (set to 1 for lower-left corner only)
+subsamp_postscript              = 1              # subsampling of background velocity model in PostScript snapshots
+factor_subsample_image          = 1              # factor to subsample color images output by the code (useful for very large models)
+POWER_DISPLAY_COLOR             = 0.30d0         # non linear display to enhance small amplitudes in color images
+DRAW_WATER_CONSTANT_BLUE_IN_JPG = .true.         # display acoustic layers as constant blue in JPEG images, because they likely correspond to water
+sizemax_arrows                  = 1.d0           # maximum size of arrows on vector plots in cm
+US_LETTER                       = .false.        # US letter paper or European A4
+USE_SNAPSHOT_NUMBER_IN_FILENAME = .false.        # use snapshot number in the file name of JPEG color snapshots instead of the time step
+gnuplot                         = .false.        # generate a GNUPLOT file for the grid
+output_grid                     = .false.        # save the grid in a text file or not
+output_energy                   = .false.        # compute and output acoustic and elastic energy (slows down the code significantly)
+output_wavefield_snapshot       = .false.        # output Ux,Uy,Uz text file for each output time (big files)
+
+# velocity and density models
+nbmodels                        = 1              # nb of different models
+# define models as
+# I:   (model_number 1 rho Vp Vs 0 0 QKappa Qmu 0 0 0 0 0 0) or
+# II:  (model_number 2 rho c11 c13 c15 c33 c35 c55 0 0 0 0 0 0) or
+# III: (model_number 3 rhos rhof phi c kxx kxz kzz Ks Kf Kfr etaf mufr Qmu).
+# For istropic elastic/acoustic material use I and set Vs to zero to make a given model acoustic, for anisotropic elastic use II,
+# and for isotropic poroelastic material use III. The mesh can contain acoustic, elastic, and poroelastic models simultaneously.
+1 1 1.d0 2.d0 1.d0 0 0 9999.d0 9999.d0 0 0 0 0 0 0
+
+# external mesh or not
+read_external_mesh              = .false.
+
+# absorbing boundary active or not
+absorbing_conditions            = .true.
+
+# for horizontal periodic conditions: detect common points between left and right edges
+ADD_PERIODIC_CONDITIONS         = .false.
+
+# horizontal periodicity distance for periodic conditions
+PERIODIC_horiz_dist             = 0.3597d0
+
+# grid point detection tolerance for periodic conditions
+PERIODIC_DETECT_TOL             = 3.3334d-6
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR EXTERNAL MESHING
+
+# data concerning mesh, when generated using third-party app (more info in README)
+# (see also absorbing_conditions above)
+mesh_file                       = ./DATA/mesh/canyon_mesh_file   # file containing the mesh
+nodes_coords_file               = ./DATA/mesh/canyon_nodes_coords_file   # file containing the nodes coordinates
+materials_file                  = ./DATA/mesh/canyon_materials_file   # file containing the material number for each element
+free_surface_file               = ./DATA/mesh/canyon_free_surface_file   # file containing the free surface
+absorbing_surface_file          = ./DATA/mesh/canyon_absorbing_surface_file   # file containing the absorbing surface
+tangential_detection_curve_file = ./DATA/courbe_eros_nodes # file containing the curve delimiting the velocity model
+
+#-----------------------------------------------------------------------------
+# PARAMETERS FOR INTERNAL MESHING
+
+# file containing interfaces for internal mesh
+interfacesfile                  = interfaces_Rayleigh_flat.dat
+
+# geometry of the model (origin lower-left corner = 0,0) and mesh description
+xmin                            = 0.d0           # abscissa of left side of the model
+xmax                            = 19.d0          # abscissa of right side of the model
+nx                              = 60             # number of elements along X
+
+# absorbing boundary parameters (see absorbing_conditions above)
+absorbbottom                    = .true.
+absorbright                     = .true.
+absorbtop                       = .false.
+absorbleft                      = .true.
+
+# define the different regions of the model in the (nx,nz) spectral element mes
+nbregions                       = 1              # nb of regions and model number for each
+1 60  1 28 1

Added: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave/README
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave/README	                        (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave/README	2011-12-28 14:36:37 UTC (rev 19323)
@@ -0,0 +1,29 @@
+----------------------------------------------------------------------
+README
+----------------------------------------------------------------------
+
+This example propagates an initial (incident) Rayleigh wave along the surface of a 2D homogeneous medium, similar to Figure 3 of Komatitsch and Tromp (1999) but without the canyon.
+
+TO RUN:
+
+0. Read the user manual in SPECFEM2D/doc/manual_SPECFEM2D.pdf
+
+1. in the SPECFEM2D root directory, configure, e.g.,
+   ./configure FC=gfortran
+
+2. compile:
+   make all
+
+3. cd EXAMPLES/Rayleigh_wave
+
+4. execute script to run mesher and solver for the PSV case:
+   ./process.sh
+
+5. check out the output files in the local directory OUTPUT_FILES
+
+References:
+-----------
+
+Dimitri Komatitsch, Jean-Pierre Vilotte, Rossana Vai, José M. Castillo-Covarrubias and Francisco J. Sánchez-Sesma, The Spectral Element method for elastic wave equations: application to 2D and 3D seismic problems, International Journal for Numerical Methods in Engineering, vol. 45, p. 1139-1164 (1999).
+
+----------------------------------------------------------------------

Added: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave/SOURCE_Rayleigh_2D
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave/SOURCE_Rayleigh_2D	                        (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave/SOURCE_Rayleigh_2D	2011-12-28 14:36:37 UTC (rev 19323)
@@ -0,0 +1,13 @@
+# source parameters
+source_surf                     = .false.        # source inside the medium or at the surface
+xs                              = 3.5            # source location x in meters
+zs                              = 9.             # source location z in meters
+source_type                     = 3              # 1 for plane P waves, 2 for plane SV waves, 3 for Rayleigh wave
+time_function_type              = 1              # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4, Heaviside = 5
+f0                              = 1.0            # dominant source frequency (Hz) if not Dirac or Heaviside
+tshift                          = 0.0            # time shift when multi sources (if one source, must be zero)
+angleforce                      = 0.             # angle of the source (for a force only)
+Mxx                             = 1.             # Mxx component (for a moment tensor source only)
+Mzz                             = 1.             # Mzz component (for a moment tensor source only)
+Mxz                             = 0.             # Mxz component (for a moment tensor source only)
+factor                          = 1.d10          # amplification factor

Added: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave/interfaces_Rayleigh_flat.dat
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave/interfaces_Rayleigh_flat.dat	                        (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave/interfaces_Rayleigh_flat.dat	2011-12-28 14:36:37 UTC (rev 19323)
@@ -0,0 +1,26 @@
+#
+# number of interfaces
+#
+ 2
+#
+# for each interface below, we give the number of points and then x,y for each point
+#
+#
+# interface number 1 (bottom of the mesh)
+#
+ 2
+ 0 0
+ 19 0
+#
+# interface number 2 (topography, top of the mesh)
+#
+ 2
+    0 9
+ 19 9
+#
+# for each layer, we give the number of spectral elements in the vertical direction
+#
+#
+# layer number 1 (bottom layer)
+#
+ 28

Added: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave/process.sh
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave/process.sh	                        (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave/process.sh	2011-12-28 14:36:37 UTC (rev 19323)
@@ -0,0 +1,69 @@
+#!/bin/bash
+#
+# script runs mesher and solver (in serial)
+# using this example setup
+#
+
+echo "running example: `date`"
+currentdir=`pwd`
+
+echo
+echo "(will take a few minutes)"
+echo
+
+# sets up directory structure in current example directoy
+echo
+echo "   setting up example..."
+echo
+
+mkdir -p OUTPUT_FILES
+mkdir -p DATA
+
+# sets up local DATA/ directory
+cd DATA/
+cp ../Par_file_Rayleigh_2D Par_file
+cp ../interfaces_Rayleigh_flat.dat .
+cp ../SOURCE_Rayleigh_2D SOURCE
+cd ../
+
+# cleans output files
+rm -rf OUTPUT_FILES/*
+
+# compiles executables in root directory
+cd ../../
+make > tmp.log
+cd $currentdir
+
+# links executables
+rm -f xmeshfem2D xspecfem2D
+ln -s ../../bin/xmeshfem2D
+ln -s ../../bin/xspecfem2D
+
+# stores setup
+cp DATA/Par_file OUTPUT_FILES/
+cp DATA/SOURCE OUTPUT_FILES/
+
+# runs database generation
+echo
+echo "  running mesher..."
+echo
+./xmeshfem2D > OUTPUT_FILES/output_mesher.txt
+./xmeshfem2D
+
+# runs simulation
+echo
+echo "  running solver..."
+echo
+#./xspecfem2D > OUTPUT_FILES/output_solver.txt
+./xspecfem2D
+
+# stores output
+cp DATA/SOURCE_xz.dat OUTPUT_FILES/
+cp DATA/STATIONS OUTPUT_FILES/
+cp DATA/STATIONS_target OUTPUT_FILES/
+
+echo
+echo "see results in directory: OUTPUT_FILES/"
+echo
+echo "done"
+echo `date`


Property changes on: seismo/2D/SPECFEM2D/trunk/EXAMPLES/Rayleigh_wave/process.sh
___________________________________________________________________
Name: svn:executable
   + *

Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/attenuation/README
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/attenuation/README	2011-12-28 03:32:05 UTC (rev 19322)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/attenuation/README	2011-12-28 14:36:37 UTC (rev 19323)
@@ -2,7 +2,7 @@
 README
 ----------------------------------------------------------------------
 
-This default example creates the 2D attenuation example of Komatitsch and Tromp (1999, Figure 16).
+This example creates the 2D attenuation benchmark of Komatitsch and Tromp (1999, Figure 16).
 
 TO RUN:
 
@@ -14,7 +14,7 @@
 2. compile:
    make all
 
-3. cd EXAMPLES/attenuation/
+3. cd EXAMPLES/attenuation
 
 4. execute script to run mesher and solver for the PSV case:
    ./process.sh

Modified: seismo/2D/SPECFEM2D/trunk/setup/constants.h.in
===================================================================
--- seismo/2D/SPECFEM2D/trunk/setup/constants.h.in	2011-12-28 03:32:05 UTC (rev 19322)
+++ seismo/2D/SPECFEM2D/trunk/setup/constants.h.in	2011-12-28 14:36:37 UTC (rev 19323)
@@ -48,6 +48,9 @@
   logical, parameter :: ACTUALLY_IMPLEMENT_PERM_INN = .false.
   logical, parameter :: ACTUALLY_IMPLEMENT_PERM_WHOLE = .true.
 
+! create file DATA/model_velocity.dat_output or not (can be huge and slow, thus off by default)
+  logical, parameter :: OUTPUT_MODEL_VELOCITY_FILE = .false.
+
 ! in PostScript files about mesh quality, show all elements above this threshold
   double precision, parameter :: THRESHOLD_POSTSCRIPT = 95.d0 / 100.d0
 

Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/paco_beyond_critical.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/paco_beyond_critical.f90	2011-12-28 03:32:05 UTC (rev 19322)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/paco_beyond_critical.f90	2011-12-28 14:36:37 UTC (rev 19323)
@@ -17,7 +17,7 @@
 subroutine paco_beyond_critical(coord,nglob,deltat,NSTEP_global,angleforce,&
      f0,cp_local,cs_local,INCLUDE_ATTENUATION,QD,source_type,v0x_left,v0z_left,v0x_right,v0z_right,&
      v0x_bot,v0z_bot,t0x_left,t0z_left,t0x_right,t0z_right,t0x_bot,t0z_bot,left_bound,right_bound,&
-     bot_bound,nleft,nright,nbot,displ_elastic,veloc_elastic,accel_elastic)
+     bot_bound,nleft,nright,nbot,displ_elastic,veloc_elastic,accel_elastic,x_source)
 
   implicit none
 
@@ -59,24 +59,22 @@
 
   complex(selected_real_kind(15,300)), dimension(:),allocatable::Field_Ux,Field_Uz,Field_Tx,Field_Tz
 
-  double precision :: TS
+  double precision :: TS,offset,x_source
 
-! to move the place where the wave reflects on free surface (offset too)
-  double precision :: offset
-
 ! size of the model
   xmin=minval(coord(1,:))
   xmax=maxval(coord(1,:))
   zmin=minval(coord(2,:))
   zmax=maxval(coord(2,:))
 
-! offset of the origin of time of the Ricker (equivalent to t0 in SPECFEM2D)
-  offset=4.d0*(xmax-xmin)/5.d0
-  TS=2.d0/f0
+  TS=1.2d0/f0
 
-! dominant period of the Ricker (equivalent to 1/f0 in SPECFEM2D)
+! dominant period of the Ricker
   TP=1.d0/f0
 
+! offset to move the initial location of the source in the horizontal direction of the mesh
+  offset = x_source
+
 ! find optimal period
 ! if period is too small, you should see several initial plane wave on your initial field
   delta_in_period=2.d0
@@ -93,7 +91,7 @@
      print *, "you must take a deltat that is a power of two (power can be negative)"
      print *, "for example you can take", DT
      stop "cannot go further, restart with new deltat"
-  end if
+  endif
 
   DT=deltat/2.d0
 
@@ -163,7 +161,7 @@
         allocate(local_pt(npt))
         local_pt=bot_bound
         NSTEP_local=NSTEP_global
-     end if
+     endif
 
 ! to distinguish all model case and boundary case
      allocate(temp_field(NSTEP_local))
@@ -191,28 +189,28 @@
      else
         VNZ = 0.d0
         VNX = 0.d0
-     end if
+     endif
 
 
      do indice=1,npt
 
         if (FLAG==0) then
            inode=indice
-           X=coord(1,indice)-offset
+           X=coord(1,indice) - offset
 ! specfem coordinate axes are implemented from bottom to top whereas for this code
 ! we need from top to bottom
            Z=zmax-coord(2,indice)
         else
            inode=local_pt(indice)
-           X=coord(1,inode)-offset
+           X=coord(1,inode) - offset
 ! specfem coordinate axes are implemented from bottom to top whereas for this code
 ! we need from top to bottom
            Z=zmax-coord(2,inode)
-        end if
+        endif
 
         if (mod(indice,500)==0) then
            print *,indice,"points have been computed out of ",npt
-        end if
+        endif
 
 !
 ! first handle the particular case of zero frequency
@@ -235,7 +233,7 @@
         if (FLAG/=0) then
            Field_Tx(1)=TX
            Field_Tz(1)=TZ
-        end if
+        endif
 
 !
 ! then loop on all the other discrete frequencies
@@ -274,7 +272,7 @@
            if (FLAG/=0) then
               Field_Tx(J+1)=TX
               Field_Tz(J+1)=TZ
-           end if
+           endif
 
         enddo
 
@@ -331,7 +329,7 @@
            t0x_bot(indice,:)=temp_field(:)
            call paco_convolve_fft(Field_Tz,4,NSTEP_local,dt,NFREC,temp_field,TP,TS)
            t0z_bot(indice,:)=temp_field(:)
-        end if
+        endif
      enddo
 
      deallocate(temp_field)
@@ -363,24 +361,24 @@
      AUX1=A1*EXP(UI*(AM*Z-AL*X))         ! campo P incidente
   else
      AUX1=CMPLX(0.0d0)
-  end if
+  endif
   if (A2/=0.0d0) then
      AUX2=A2*EXP(-UI*(AM*Z+AL*X)) *1.0d0      ! campo P reflejado
   else
      AUX2=CMPLX(0.0d0)
-  end if
+  endif
   FI1=AUX1+AUX2
   FI2=AUX1-AUX2
   if (B1/=0.0d0) then
      AUX1=B1*EXP(UI*(AK*Z-AL*X))            ! campo S incidente
   else
      AUX1=CMPLX(0.0d0)
-  end if
+  endif
   if (B2/=0.0d0) then
      AUX2=B2*EXP(-UI*(AK*Z+AL*X)) *1.0d0      ! campo S reflejado
   else
      AUX2=CMPLX(0.0d0)
-  end if
+  endif
   PS1=AUX1+AUX2
   PS2=AUX1-AUX2
 

Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90	2011-12-28 03:32:05 UTC (rev 19322)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90	2011-12-28 14:36:37 UTC (rev 19323)
@@ -1,4 +1,6 @@
+
   program specfem2D
+
 !========================================================================
 !
 !                   S P E C F E M 2 D  Version 6 . 2
@@ -3038,7 +3040,7 @@
               f0(1),cploc,csloc,TURN_ATTENUATION_ON,QKappa_attenuation(1),source_type(1),v0x_left,v0z_left, &
               v0x_right,v0z_right,v0x_bot,v0z_bot,t0x_left,t0z_left,t0x_right,t0z_right, &
               t0x_bot,t0z_bot,left_bound(1:count_left),right_bound(1:count_right),bot_bound(1:count_bottom), &
-              count_left,count_right,count_bottom,displ_elastic,veloc_elastic,accel_elastic)
+              count_left,count_right,count_bottom,displ_elastic,veloc_elastic,accel_elastic,x_source(1))
 
       deallocate(left_bound)
       deallocate(right_bound)
@@ -7096,16 +7098,12 @@
 !----  close energy file
   if(output_energy .and. myrank == 0) close(IOUT_ENERGY)
 
-  if (.not. any_poroelastic) then
+  if (OUTPUT_MODEL_VELOCITY_FILE .and. .not. any_poroelastic) then
     open(unit=1001,file='DATA/model_velocity.dat_output',status='unknown')
     if ( .NOT. assign_external_model) then
       allocate(rho_local(ngllx,ngllz,nspec)); rho_local=0.
       allocate(vp_local(ngllx,ngllz,nspec)); vp_local=0.
       allocate(vs_local(ngllx,ngllz,nspec)); vs_local=0.
-!!      write(1001,*) nglob
-!!      do iglob = 1,nglob
-!!         write(1001,*) coord(1,iglob),coord(2,iglob),rho_global(iglob),vp_global(iglob),vs_global(iglob)
-!!      end do
       do ispec = 1,nspec
         do j = 1,NGLLZ
           do i = 1,NGLLX
@@ -7119,10 +7117,6 @@
         end do
       end do
     else
-!!     write(1001,*) nglob
-!!  do iglob = 1,nglob
-!!     write(1001,*) coord(1,iglob),coord(2,iglob),rhoext_global(iglob),vpext_global(iglob),vsext_global(iglob)
-!!  end do
       do ispec = 1,nspec
         do j = 1,NGLLZ
           do i = 1,NGLLX



More information about the CIG-COMMITS mailing list