[cig-commits] r20196 - in seismo/3D/SPECFEM3D/trunk/utils/unused_routines: . ampuero_implicit_Clayton_now_done_fixed_in_our_official_SVN_code
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Sun May 20 06:54:07 PDT 2012
Author: dkomati1
Date: 2012-05-20 06:54:06 -0700 (Sun, 20 May 2012)
New Revision: 20196
Added:
seismo/3D/SPECFEM3D/trunk/utils/unused_routines/ampuero_implicit_Clayton_now_done_fixed_in_our_official_SVN_code/
seismo/3D/SPECFEM3D/trunk/utils/unused_routines/ampuero_implicit_Clayton_now_done_fixed_in_our_official_SVN_code/ampuero_fix_Stacey_Hughes_1987_implicit_ABC_problem_specfem3D.f90
Removed:
seismo/3D/SPECFEM3D/trunk/utils/unused_routines/ampuero_implicit_Clayton/
seismo/3D/SPECFEM3D/trunk/utils/unused_routines/ampuero_implicit_Clayton_now_done_fixed_in_our_official_SVN_code/ampuero_implicit_ABC_specfem3D.f90
Log:
renamed directory ampuero_implicit_Clayton_now_done_fixed_in_our_official_SVN_code
Copied: seismo/3D/SPECFEM3D/trunk/utils/unused_routines/ampuero_implicit_Clayton_now_done_fixed_in_our_official_SVN_code (from rev 20195, seismo/3D/SPECFEM3D/trunk/utils/unused_routines/ampuero_implicit_Clayton)
Copied: seismo/3D/SPECFEM3D/trunk/utils/unused_routines/ampuero_implicit_Clayton_now_done_fixed_in_our_official_SVN_code/ampuero_fix_Stacey_Hughes_1987_implicit_ABC_problem_specfem3D.f90 (from rev 20195, seismo/3D/SPECFEM3D/trunk/utils/unused_routines/ampuero_implicit_Clayton/ampuero_implicit_ABC_specfem3D.f90)
===================================================================
--- seismo/3D/SPECFEM3D/trunk/utils/unused_routines/ampuero_implicit_Clayton_now_done_fixed_in_our_official_SVN_code/ampuero_fix_Stacey_Hughes_1987_implicit_ABC_problem_specfem3D.f90 (rev 0)
+++ seismo/3D/SPECFEM3D/trunk/utils/unused_routines/ampuero_implicit_Clayton_now_done_fixed_in_our_official_SVN_code/ampuero_fix_Stacey_Hughes_1987_implicit_ABC_problem_specfem3D.f90 2012-05-20 13:54:06 UTC (rev 20196)
@@ -0,0 +1,1853 @@
+ program specfem3D
+!=====================================================================
+!
+! S p e c f e m 3 D B a s i n V e r s i o n 1 . 1
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Princeton University, USA and University of Pau / CNRS / INRIA
+! (c) California Institute of Technology October 2002
+!
+! 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.
+!
+!=====================================================================
+!
+! United States Government Sponsorship Acknowledged.
+!
+!jpampuero Mon Oct 20 23:09:00 EDT 2003
+!jpampuero I am implementing the FAULT boundary conditions in place of one
+!jpampuero of the absorbing boundaries.
+!jpampuero It is a quick and dirty implementation, for SCEC 1st benchmark test.
+!jpampuero No split nodes need to be implemented.
+!jpampuero The model is restricted to be symetric.
+!jpampuero
+!jpampuero All lines containing "!jpampuero" were touched.
+!jpampuero Lines containing "!!jpampuero" were only commented out.
+!jpampuero Most of my variables names have ThisFormat
+
+ implicit none
+
+! standard include of the MPI library
+ include 'mpif.h'
+
+ include "constants.h"
+ include "precision.h"
+
+! include values created by the mesher
+ include "OUTPUT_FILES/values_from_mesher.h"
+
+!jpampuero Fault parameters
+ include "fault_parameters.h" !jpampuero
+
+!=======================================================================!
+! !
+! specfem3D is a 3-D spectral-element solver for a basin. !
+! It uses a mesh generated by program meshfem3D !
+! !
+!=======================================================================!
+!
+! If you use this code for your own research, please send an email
+! to Jeroen Tromp <jtromp at gps.caltech.edu> for information, and cite:
+!
+! @article{KoTr99,
+! author={D. Komatitsch and J. Tromp},
+! year=1999,
+! title={Introduction to the spectral-element method for 3-{D} seismic wave propagation},
+! journal={Geophys. J. Int.},
+! volume=139,
+! pages={806-822}}
+!
+! Evolution of the code:
+! ---------------------
+!
+! MPI v. 1.1 Dimitri Komatitsch, Caltech, October 2002: Zhu's Moho map, scaling
+! of Vs with depth, Hauksson's regional model, attenuation, oceans, movies
+! MPI v. 1.0 Dimitri Komatitsch, Caltech, May 2002: first MPI version
+! based on global code
+
+! memory variables and standard linear solids for attenuation
+ double precision, dimension(N_SLS) :: tau_mu_dble,tau_sigma_dble,beta_dble
+ double precision factor_scale_dble,one_minus_sum_beta_dble
+ real(kind=CUSTOM_REAL), dimension(NUM_REGIONS_ATTENUATION,N_SLS) :: tau_mu,tau_sigma,beta
+ real(kind=CUSTOM_REAL), dimension(NUM_REGIONS_ATTENUATION) :: factor_scale,one_minus_sum_beta
+
+ real(kind=CUSTOM_REAL), dimension(NUM_REGIONS_ATTENUATION,N_SLS) :: tauinv,alphaval,betaval,gammaval,factor_common
+ integer iattenuation
+ double precision scale_factor
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION,N_SLS) :: &
+ R_xx,R_yy,R_xy,R_xz,R_yz
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION) :: &
+ epsilondev_xx,epsilondev_yy,epsilondev_xy,epsilondev_xz,epsilondev_yz
+
+ integer NPOIN2DMAX_XY
+
+ integer, dimension(:), allocatable :: ibelm_xmin,ibelm_xmax, &
+ ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: &
+ jacobian2D_xmin,jacobian2D_xmax,jacobian2D_ymin,jacobian2D_ymax, &
+ jacobian2D_bottom,jacobian2D_top
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: normal_xmin, &
+ normal_xmax,normal_ymin,normal_ymax, &
+ normal_bottom,normal_top
+
+! buffers for send and receive between faces of the slices and the chunks
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: buffer_send_faces,buffer_received_faces
+
+! -----------------
+
+! mesh parameters
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
+
+!! DK DK NOJACOBIAN real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,1) :: &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,jacobian
+ real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: xstore,ystore,zstore
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
+ kappastore,mustore
+
+! Stacey
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
+ rho_vp,rho_vs
+
+! local to global mapping
+ integer, dimension(NSPEC_AB) :: idoubling
+
+! mass matrix
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: rmass !jpampuero ABC
+
+! displacement, velocity, acceleration
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displ,veloc,accel
+
+ real(kind=CUSTOM_REAL) xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
+ real(kind=CUSTOM_REAL) duxdxl,duxdyl,duxdzl,duydxl,duydyl,duydzl,duzdxl,duzdyl,duzdzl
+
+ real(kind=CUSTOM_REAL) duxdxl_plus_duydyl,duxdxl_plus_duzdzl,duydyl_plus_duzdzl
+ real(kind=CUSTOM_REAL) duxdyl_plus_duydxl,duzdxl_plus_duxdzl,duzdyl_plus_duydzl
+
+ real(kind=CUSTOM_REAL) sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz
+
+ real(kind=CUSTOM_REAL) hp1,hp2,hp3
+ real(kind=CUSTOM_REAL) fac1,fac2,fac3
+ real(kind=CUSTOM_REAL) lambdal,kappal,mul,lambdalplus2mul
+
+ real(kind=CUSTOM_REAL) tempx1l,tempx2l,tempx3l
+ real(kind=CUSTOM_REAL) tempy1l,tempy2l,tempy3l
+ real(kind=CUSTOM_REAL) tempz1l,tempz2l,tempz3l
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
+ tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3
+
+! for attenuation
+ real(kind=CUSTOM_REAL) R_xx_val,R_yy_val
+ real(kind=CUSTOM_REAL) factor_loc,alphaval_loc,betaval_loc,gammaval_loc,Sn,Snp1
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: epsilondev_xx_loc, &
+ epsilondev_yy_loc, epsilondev_xy_loc, epsilondev_xz_loc, epsilondev_yz_loc
+ real(kind=CUSTOM_REAL) epsilon_trace_over_3
+
+ integer l
+ integer i_SLS
+ real(kind=CUSTOM_REAL) one_minus_sum_beta_use,minus_sum_beta
+ integer iselected
+
+! --------
+
+! parameters for the source
+ integer it,isource
+ integer, dimension(:), allocatable :: islice_selected_source,ispec_selected_source
+ integer yr,jda,ho,mi
+ real(kind=CUSTOM_REAL) stf_used
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: sourcearray
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: sourcearrays
+ double precision sec,stf
+ double precision, dimension(:), allocatable :: Mxx,Myy,Mzz,Mxy,Mxz,Myz
+ double precision, dimension(:), allocatable :: xi_source,eta_source,gamma_source
+ double precision, dimension(:), allocatable :: t_cmt,hdur
+ double precision, dimension(:), allocatable :: utm_x_source,utm_y_source
+ double precision, external :: comp_source_time_function
+
+
+! time scheme
+ real(kind=CUSTOM_REAL) deltat,deltatover2,deltatsqover2
+
+! receiver information
+ integer nrec,nrec_local,nrec_tot_found
+ integer irec_local
+ integer, allocatable, dimension(:) :: islice_selected_rec,ispec_selected_rec,number_receiver_global
+ double precision, allocatable, dimension(:) :: xi_receiver,eta_receiver
+ double precision hlagrange
+
+! timing information for the stations
+ double precision, allocatable, dimension(:,:,:) :: nu
+ character(len=8), allocatable, dimension(:) :: station_name,network_name
+
+! seismograms
+ double precision uxd,uyd,uzd
+ double precision vxd,vyd,vzd
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: seismograms_d,seismograms_v
+
+ integer i,j,k,ispec,irec,iglob
+
+! Gauss-Lobatto-Legendre points of integration and weights
+ double precision, dimension(NGLLX) :: xigll,wxgll
+ double precision, dimension(NGLLY) :: yigll,wygll
+ double precision, dimension(NGLLZ) :: zigll,wzgll
+
+! array with derivatives of Lagrange polynomials and precalculated products
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx,hprimewgll_xx
+ real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLY) :: hprime_yy,hprimewgll_yy
+ real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ) :: hprime_zz,hprimewgll_zz
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY) :: wgllwgll_xy
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ) :: wgllwgll_xz
+ real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLZ) :: wgllwgll_yz
+
+! Lagrange interpolators at receivers
+ double precision, dimension(:), allocatable :: hxir,hetar,hpxir,hpetar
+ double precision, dimension(:,:), allocatable :: hxir_store,hetar_store
+
+! 2-D addressing and buffers for summation between slices
+ integer, dimension(:), allocatable :: iboolleft_xi, &
+ iboolright_xi,iboolleft_eta,iboolright_eta
+
+! for addressing of the slices
+ integer, dimension(:,:), allocatable :: addressing
+ integer, dimension(:), allocatable :: iproc_xi_slice,iproc_eta_slice
+
+! proc numbers for MPI
+ integer myrank,sizeprocs,ier
+
+ integer npoin2D_xi,npoin2D_eta
+
+ integer iproc_xi,iproc_eta,iproc,iproc_read
+
+! maximum of the norm of the displacement
+ real(kind=CUSTOM_REAL) Usolidnorm,Usolidnorm_all
+
+! timer MPI
+ integer ihours,iminutes,iseconds,int_tCPU
+ double precision time_start,tCPU
+
+! parameters read from parameter file
+ integer NER_TOP_LAYER,NER_MIDDLE_LAYER,NER_BOTTOM_HALFSPACE,NEX_ETA,NEX_XI, &
+ NPROC_ETA,NPROC_XI,NSEIS,NSTEP,UTM_PROJECTION_ZONE
+ integer NSOURCES
+
+ double precision UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK
+ double precision DT,LAT_MIN,LAT_MAX,LONG_MIN,LONG_MAX
+ logical ATTENUATION,STACEY_ABS_CONDITIONS
+
+ character(len=256) LOCAL_PATH,prname
+
+! parameters deduced from parameters read from file
+ integer NPROC,NEX_PER_PROC_XI,NEX_PER_PROC_ETA
+ integer NER
+
+ integer NSPEC2D_A_XI,NSPEC2D_B_XI, &
+ NSPEC2D_A_ETA,NSPEC2D_B_ETA, &
+ NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX, &
+ NSPEC2D_BOTTOM,NSPEC2D_TOP, &
+ NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX
+
+! names of the data files for all the processors in MPI
+ character(len=256) outputname
+
+! Stacey conditions put back
+ integer nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax,ispec2D
+ integer, dimension(:,:), allocatable :: nimin,nimax,njmin,njmax,nkmin_xi,nkmin_eta
+ real(kind=CUSTOM_REAL) vx,vy,vz,nx,ny,nz,tx,ty,tz,vn,weight
+
+!jpampuero Fault variables
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: FaultStrength,FaultTxInit,FaultTx & !jpampuero
+ , FaultTzInit,FaultTz,FaultMuS,FaultMuD,FaultSWR ,FaultB !jpampuero
+ real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: FaultDisplX,FaultDisplZ & !jpampuero
+ ,FaultVelocX,FaultVelocZ & !jpampuero
+ ,FaultStressX,FaultStressZ !jpampuero
+ real(kind=CUSTOM_REAL) :: FaultDVxFree,FaultDVzFree,FaultZ,va,ta,txa,tza,dirx,dirz !jpampuero
+ integer, dimension(FaultOutputNPTS) :: FaultIglobOut,FaultIgllOut,FaultKgllOut,FaultElmtOut !jpampuero
+ character(len=70) FaultOutputName !jpampuero
+ character(len=70), dimension(FaultOutputNPTS) :: FaultOutputPtsName !jpampuero
+ logical :: FaultInThisProc,NotInRuptureRegion,InNucleationRegion !jpampuero
+
+! ************** PROGRAM STARTS HERE **************
+
+! initialize the MPI communicator and start the NPROC MPI processes.
+! sizeprocs returns number of processes started
+! (should be equal to NPROC)
+! myrank is the rank of each process, between 0 and sizeprocs-1.
+! as usual in MPI, process 0 is in charge of coordinating everything
+! and also takes care of the main output
+ call MPI_INIT(ier)
+ call MPI_COMM_SIZE(MPI_COMM_WORLD,sizeprocs,ier)
+ call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ier)
+
+! read the parameter file
+ call read_parameter_file(LAT_MIN,LAT_MAX,LONG_MIN,LONG_MAX, &
+ UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK, &
+ NER_TOP_LAYER,NER_MIDDLE_LAYER,NER_BOTTOM_HALFSPACE, &
+ NEX_ETA,NEX_XI,NPROC_ETA,NPROC_XI,NSEIS,NSTEP,UTM_PROJECTION_ZONE,DT, &
+ ATTENUATION,LOCAL_PATH,NSOURCES,STACEY_ABS_CONDITIONS)
+
+! compute other parameters based upon values read
+ call compute_parameters(NER,NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA, &
+ NPROC,NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
+ NER_BOTTOM_HALFSPACE,NER_MIDDLE_LAYER,NER_TOP_LAYER, &
+ NSPEC_AB,NSPEC2D_A_XI,NSPEC2D_B_XI, &
+ NSPEC2D_A_ETA,NSPEC2D_B_ETA, &
+ NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
+ NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NGLOB_AB)
+
+!jpampuero moved this section (used to be right before time loop)
+!jpampuero need some of these during initialization
+! distinguish whether single or double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ deltat = sngl(DT)
+ else
+ deltat = DT
+ endif
+ deltatover2 = deltat/2.
+ deltatsqover2 = deltat*deltat/2.
+
+! open main output file, only written to by process 0
+ if(myrank == 0 .and. IMAIN /= ISTANDARD_OUTPUT) &
+ open(unit=IMAIN,file='OUTPUT_FILES/output_solver.txt',status='unknown')
+
+ if(myrank == 0) then
+
+ write(IMAIN,*)
+ write(IMAIN,*) '**********************************************'
+ write(IMAIN,*) '**** Specfem 3-D Solver - MPI version f90 ****'
+ write(IMAIN,*) '**********************************************'
+ write(IMAIN,*)
+ write(IMAIN,*)
+
+ if(FIX_UNDERFLOW_PROBLEM) write(IMAIN,*) 'Fixing slow underflow trapping problem using small initial field'
+
+ write(IMAIN,*)
+ write(IMAIN,*) 'There are ',sizeprocs,' MPI processes'
+ write(IMAIN,*) 'Processes are numbered from 0 to ',sizeprocs-1
+ write(IMAIN,*)
+
+ write(IMAIN,*) 'There are ',NEX_XI,' elements along xi'
+ write(IMAIN,*) 'There are ',NEX_ETA,' elements along eta'
+ write(IMAIN,*)
+ write(IMAIN,*) 'There are ',NPROC_XI,' slices along xi'
+ write(IMAIN,*) 'There are ',NPROC_ETA,' slices along eta'
+ write(IMAIN,*) 'There is a total of ',NPROC,' slices'
+
+ write(IMAIN,*)
+ write(IMAIN,*) ' NDIM = ',NDIM
+ write(IMAIN,*)
+ write(IMAIN,*) ' NGLLX = ',NGLLX
+ write(IMAIN,*) ' NGLLY = ',NGLLY
+ write(IMAIN,*) ' NGLLZ = ',NGLLZ
+ write(IMAIN,*)
+
+! write information about precision used for floating-point operations
+ if(CUSTOM_REAL == SIZE_REAL) then
+ write(IMAIN,*) 'using single precision for the calculations'
+ else
+ write(IMAIN,*) 'using double precision for the calculations'
+ endif
+ write(IMAIN,*)
+ write(IMAIN,*) 'smallest and largest possible floating-point numbers are: ',tiny(1._CUSTOM_REAL),huge(1._CUSTOM_REAL)
+ write(IMAIN,*)
+
+ endif
+
+! check that the code is running with the requested nb of processes
+ if(sizeprocs /= NPROC) call exit_MPI(myrank,'wrong number of MPI processes')
+
+! check that we have more than 0 and less than 1000 sources
+ if(NSOURCES < 0 .or. NSOURCES > 999) call exit_MPI(myrank,'invalid number of sources')
+
+! dynamic allocation of arrays
+
+! 2-D addressing and buffers for summation between slices, and point codes
+! use number of elements found in the mantle since it is the largest region
+
+! crust and mantle
+ allocate(iboolleft_xi(NPOIN2DMAX_XMIN_XMAX))
+ allocate(iboolright_xi(NPOIN2DMAX_XMIN_XMAX))
+ allocate(iboolleft_eta(NPOIN2DMAX_YMIN_YMAX))
+ allocate(iboolright_eta(NPOIN2DMAX_YMIN_YMAX))
+
+! for addressing of the slices
+ allocate(addressing(0:NPROC_XI-1,0:NPROC_ETA-1))
+ allocate(iproc_xi_slice(0:NPROC-1))
+ allocate(iproc_eta_slice(0:NPROC-1))
+
+! open file with global slice number addressing
+ open(unit=IIN,file='OUTPUT_FILES/addressing.txt',status='old')
+ do iproc = 0,NPROC-1
+ read(IIN,*) iproc_read,iproc_xi,iproc_eta
+ if(iproc_read /= iproc) call exit_MPI(myrank,'incorrect slice number read')
+ addressing(iproc_xi,iproc_eta) = iproc
+ iproc_xi_slice(iproc) = iproc_xi
+ iproc_eta_slice(iproc) = iproc_eta
+ enddo
+ close(IIN)
+
+! determine local slice coordinates using addressing
+ iproc_xi = iproc_xi_slice(myrank)
+ iproc_eta = iproc_eta_slice(myrank)
+
+ FaultInThisProc = iproc_eta==0 !jpampuero
+
+! define maximum size for message buffers
+ NPOIN2DMAX_XY = max(NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX)
+
+ allocate(buffer_send_faces(NDIM,NPOIN2DMAX_XY))
+ allocate(buffer_received_faces(NDIM,NPOIN2DMAX_XY))
+
+! start reading the databases
+
+! read arrays created by the mesher
+
+ call read_arrays_solver(myrank,xstore,ystore,zstore, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,jacobian, &
+ rho_vp,rho_vs,kappastore,mustore,ibool,idoubling,rmass(1,:),LOCAL_PATH) !jpampuero ABC
+ rmass(2,:) = rmass(1,:) !jpampuero ABC
+ rmass(3,:) = rmass(1,:) !jpampuero ABC
+
+! check that the number of points in this slice is correct
+
+ if(minval(ibool(:,:,:,:)) /= 1 .or. maxval(ibool(:,:,:,:)) /= NGLOB_AB) &
+ call exit_MPI(myrank,'incorrect global numbering: iboolmax does not equal nglob in crust and mantle')
+
+! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+
+if (NSOURCES==0) then !jpampuero
+ allocate(hdur(1))
+ allocate(utm_x_source(1))
+ allocate(utm_y_source(1))
+ hdur(1) = 0.d0
+ utm_x_source(1) = 0.d0
+ utm_y_source(1) = 0.d0
+else
+! allocate arrays for source
+ allocate(sourcearray(NDIM,NGLLX,NGLLY,NGLLZ))
+ allocate(sourcearrays(NSOURCES,NDIM,NGLLX,NGLLY,NGLLZ))
+ allocate(islice_selected_source(NSOURCES))
+ allocate(ispec_selected_source(NSOURCES))
+ allocate(Mxx(NSOURCES))
+ allocate(Myy(NSOURCES))
+ allocate(Mzz(NSOURCES))
+ allocate(Mxy(NSOURCES))
+ allocate(Mxz(NSOURCES))
+ allocate(Myz(NSOURCES))
+ allocate(xi_source(NSOURCES))
+ allocate(eta_source(NSOURCES))
+ allocate(gamma_source(NSOURCES))
+ allocate(t_cmt(NSOURCES))
+ allocate(hdur(NSOURCES))
+ allocate(utm_x_source(NSOURCES))
+ allocate(utm_y_source(NSOURCES))
+
+! locate sources in the mesh
+ do isource = 1,NSOURCES
+ call locate_source(ibool,isource,NSOURCES,myrank,NSPEC_AB,NGLOB_AB,xstore,ystore,zstore, &
+ xigll,yigll,zigll,NPROC, &
+ sec,t_cmt(isource),yr,jda,ho,mi,utm_x_source(isource),utm_y_source(isource), &
+ NSTEP,DT,hdur(isource),Mxx(isource),Myy(isource),Mzz(isource),Mxy(isource),Mxz(isource),Myz(isource), &
+ islice_selected_source(isource),ispec_selected_source(isource), &
+ xi_source(isource),eta_source(isource),gamma_source(isource), &
+ LAT_MIN,LAT_MAX,LONG_MIN,LONG_MAX,Z_DEPTH_BLOCK,UTM_PROJECTION_ZONE)
+ enddo
+
+ if(t_cmt(1) /= 0.) call exit_MPI(myrank,'t_cmt for the first source should be zero')
+ do isource = 2,NSOURCES
+ if(t_cmt(isource) < 0.) call exit_MPI(myrank,'t_cmt should not be less than zero')
+ enddo
+endif !jpampuero
+
+ open(unit=IIN,file='DATA/STATIONS_FILTERED',status='old')
+ read(IIN,*) nrec
+ close(IIN)
+
+ if(myrank == 0) then
+ write(IMAIN,*)
+ write(IMAIN,*) 'Total number of receivers = ',nrec
+ write(IMAIN,*)
+ endif
+
+ if(nrec < 1) call exit_MPI(myrank,'need at least one receiver')
+
+! allocate memory for receiver arrays
+ allocate(islice_selected_rec(nrec))
+ allocate(ispec_selected_rec(nrec))
+ allocate(xi_receiver(nrec))
+ allocate(eta_receiver(nrec))
+ allocate(station_name(nrec))
+ allocate(network_name(nrec))
+ allocate(nu(NDIM,NDIM,nrec))
+
+! locate receivers in the crust in the mesh
+ call locate_receivers(ibool,myrank,NSPEC_AB,NGLOB_AB,idoubling, &
+ xstore,ystore,zstore,xigll,yigll, &
+ nrec,islice_selected_rec,ispec_selected_rec, &
+ xi_receiver,eta_receiver,station_name,network_name,nu, &
+ NPROC,utm_x_source(1),utm_y_source(1),UTM_PROJECTION_ZONE)
+
+! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+
+! read 2-D addressing for summation between slices with MPI
+
+ call read_arrays_buffers_solver(myrank,iboolleft_xi, &
+ iboolright_xi,iboolleft_eta,iboolright_eta, &
+ npoin2D_xi,npoin2D_eta, &
+ NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,LOCAL_PATH)
+
+! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
+
+ if(myrank == 0) then
+ write(IMAIN,*) '******************************************'
+ write(IMAIN,*) 'There are ',NEX_XI,' elements along xi'
+ write(IMAIN,*) 'There are ',NEX_ETA,' elements along eta'
+ write(IMAIN,*)
+ write(IMAIN,*) 'There are ',NPROC_XI,' slices along xi'
+ write(IMAIN,*) 'There are ',NPROC_ETA,' slices along eta'
+ write(IMAIN,*) 'There is a total of ',NPROC,' slices'
+ write(IMAIN,*) '******************************************'
+ write(IMAIN,*)
+ endif
+
+! set up GLL points, weights and derivation matrices
+ call define_derivation_matrices(xigll,yigll,zigll,wxgll,wygll,wzgll, &
+ hprime_xx,hprime_yy,hprime_zz, &
+ hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
+ wgllwgll_xy,wgllwgll_xz,wgllwgll_yz)
+
+! allocate 1-D Lagrange interpolators and derivatives
+ allocate(hxir(NGLLX))
+ allocate(hpxir(NGLLX))
+ allocate(hetar(NGLLY))
+ allocate(hpetar(NGLLY))
+
+! create name of database
+ call create_name_database(prname,myrank,LOCAL_PATH)
+
+! dynamic allocation of arrays
+
+! boundary parameters locator
+ allocate(ibelm_xmin(NSPEC2DMAX_XMIN_XMAX))
+ allocate(ibelm_xmax(NSPEC2DMAX_XMIN_XMAX))
+ allocate(ibelm_ymin(NSPEC2DMAX_YMIN_YMAX))
+ allocate(ibelm_ymax(NSPEC2DMAX_YMIN_YMAX))
+ allocate(ibelm_bottom(NSPEC2D_BOTTOM))
+ allocate(ibelm_top(NSPEC2D_TOP))
+
+ allocate(jacobian2D_xmin(NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX))
+ allocate(jacobian2D_xmax(NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX))
+ allocate(jacobian2D_ymin(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX))
+ allocate(jacobian2D_ymax(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX))
+ allocate(jacobian2D_bottom(NGLLX,NGLLY,NSPEC2D_BOTTOM))
+ allocate(jacobian2D_top(NGLLX,NGLLY,NSPEC2D_TOP))
+
+! Stacey put back
+ allocate(nimin(2,NSPEC2DMAX_YMIN_YMAX))
+ allocate(nimax(2,NSPEC2DMAX_YMIN_YMAX))
+ allocate(njmin(2,NSPEC2DMAX_XMIN_XMAX))
+ allocate(njmax(2,NSPEC2DMAX_XMIN_XMAX))
+ allocate(nkmin_xi(2,NSPEC2DMAX_XMIN_XMAX))
+ allocate(nkmin_eta(2,NSPEC2DMAX_YMIN_YMAX))
+
+! normals
+ allocate(normal_xmin(NDIM,NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX))
+ allocate(normal_xmax(NDIM,NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX))
+!!jpampuero allocate(normal_ymin(NDIM,NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX))
+ allocate(normal_ymax(NDIM,NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX))
+ allocate(normal_bottom(NDIM,NGLLX,NGLLY,NSPEC2D_BOTTOM))
+ allocate(normal_top(NDIM,NGLLX,NGLLY,NSPEC2D_TOP))
+
+!jpampuero Fault
+ if (FaultInThisProc) then
+ allocate(FaultB(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)) !jpampuero
+ allocate(FaultTx(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)) !jpampuero
+ allocate(FaultTxInit(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)) !jpampuero
+ allocate(FaultTz(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)) !jpampuero
+ allocate(FaultTzInit(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)) !jpampuero
+ allocate(FaultStrength(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)) !jpampuero
+ allocate(FaultMuS(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)) !jpampuero
+ allocate(FaultMuD(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)) !jpampuero
+ allocate(FaultSWR(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)) !jpampuero
+ allocate(FaultDisplX(FaultOutputNPTS,NSTEP)) !jpampuero
+ allocate(FaultDisplZ(FaultOutputNPTS,NSTEP)) !jpampuero
+ allocate(FaultVelocX(FaultOutputNPTS,NSTEP)) !jpampuero
+ allocate(FaultVelocZ(FaultOutputNPTS,NSTEP)) !jpampuero
+ allocate(FaultStressX(FaultOutputNPTS,NSTEP)) !jpampuero
+ allocate(FaultStressZ(FaultOutputNPTS,NSTEP)) !jpampuero
+ endif
+
+! boundary parameters
+ open(unit=27,file=prname(1:len_trim(prname))//'ibelm.bin',status='old',form='unformatted')
+ read(27) ibelm_xmin
+ read(27) ibelm_xmax
+ read(27) ibelm_ymin
+ read(27) ibelm_ymax
+ read(27) ibelm_bottom
+ read(27) ibelm_top
+ close(27)
+
+ open(unit=27,file=prname(1:len_trim(prname))//'normal.bin',status='old',form='unformatted')
+ read(27) normal_xmin
+ read(27) normal_xmax
+!!jpampuero read(27) normal_ymin
+ read(27) normal_ymax !jpampuero
+ read(27) normal_ymax
+ read(27) normal_bottom
+ read(27) normal_top
+ close(27)
+
+! Stacey put back
+ open(unit=27,file=prname(1:len_trim(prname))//'nspec2D.bin',status='unknown',form='unformatted')
+ read(27) nspec2D_xmin
+ read(27) nspec2D_xmax
+ read(27) nspec2D_ymin
+ read(27) nspec2D_ymax
+ close(27)
+
+ open(unit=27,file=prname(1:len_trim(prname))//'jacobian2D.bin',status='old',form='unformatted')
+ read(27) jacobian2D_xmin
+ read(27) jacobian2D_xmax
+ read(27) jacobian2D_ymin
+ read(27) jacobian2D_ymax
+ read(27) jacobian2D_bottom
+ read(27) jacobian2D_top
+ close(27)
+
+! Stacey put back
+! read arrays for Stacey conditions
+
+ if(STACEY_ABS_CONDITIONS) then
+ open(unit=27,file=prname(1:len_trim(prname))//'nimin.bin',status='unknown',form='unformatted')
+ read(27) nimin
+ close(27)
+
+ open(unit=27,file=prname(1:len_trim(prname))//'nimax.bin',status='unknown',form='unformatted')
+ read(27) nimax
+ close(27)
+
+ open(unit=27,file=prname(1:len_trim(prname))//'njmin.bin',status='unknown',form='unformatted')
+ read(27) njmin
+ close(27)
+
+ open(unit=27,file=prname(1:len_trim(prname))//'njmax.bin',status='unknown',form='unformatted')
+ read(27) njmax
+ close(27)
+
+ open(unit=27,file=prname(1:len_trim(prname))//'nkmin_xi.bin',status='unknown',form='unformatted')
+ read(27) nkmin_xi
+ close(27)
+
+ open(unit=27,file=prname(1:len_trim(prname))//'nkmin_eta.bin',status='unknown',form='unformatted')
+ read(27) nkmin_eta
+ close(27)
+ endif
+
+!jpampuero ======= BEGIN INITIALIZE FAULT ========================
+ if (FaultInThisProc) then
+
+ accel(:,:) = 0._CUSTOM_REAL
+
+ !OUTPUTS: locate selected fault points
+ call fault_locate_rec(FaultIglobOut,FaultIgllOut,FaultKgllOut,FaultElmtOut,FaultOutputPtsName &
+ ,nspec2D_ymin,ibelm_ymin,NSPEC2DMAX_YMIN_YMAX,ibool,xstore,zstore) !jpampuero
+
+ !Set inital stress and strength for fault
+ FaultTzInit = 0._CUSTOM_REAL
+ FaultTxInit = BackgroundShearStress
+ do ispec2D=1,nspec2D_ymin
+ ispec=ibelm_ymin(ispec2D)
+ do k=1,NGLLZ
+ do i=1,NGLLX
+
+ iglob=ibool(i,1,k,ispec)
+ NotInRuptureRegion = abs(zstore(iglob))>RuptureDepth &
+ .OR. abs(xstore(iglob))>RuptureHalfLength
+ InNucleationRegion = abs(xstore(iglob)-HypoX)<=NucleationHalfSizeX &
+ .AND. abs(zstore(iglob)-HypoZ)<=NucleationHalfSizeZ
+
+ if ( NotInRuptureRegion ) then
+ FaultMuS(i,k,ispec2D) = BarrierStrength
+ FaultMuD(i,k,ispec2D) = BarrierStrength
+ else
+ FaultMuS(i,k,ispec2D) = StaticStrength
+ FaultMuD(i,k,ispec2D) = DynamicStrength
+ if ( InNucleationRegion ) then
+ FaultTxInit(i,k,ispec2D) = NucleationShearStress
+ endif
+ endif
+
+ !WARNING: factor 2 comes from using displacement instead of slip
+ FaultSWR(i,k,ispec2D) = 2._CUSTOM_REAL*(FaultMuS(i,k,ispec2D)-FaultMuD(i,k,ispec2D))/Dc
+
+ !assemble boundary weights, temporarily store in "accel" global array
+ accel(1,iglob) = accel(1,iglob) + jacobian2D_ymin(i,k,ispec2D)*wgllwgll_xz(i,k)
+
+ enddo
+ enddo
+ enddo
+
+ do ispec2D=1,nspec2D_ymin
+ ispec=ibelm_ymin(ispec2D)
+ do k=1,NGLLZ
+ do i=1,NGLLX
+ iglob=ibool(i,1,k,ispec)
+ FaultB(i,k,ispec2D) = accel(1,iglob)
+ enddo
+ enddo
+ enddo
+
+ endif
+!jpampuero ======= END INITIALIZE FAULT ========================
+
+ do isource = 1,NSOURCES
+
+! check that the source slice number is okay
+ if(islice_selected_source(isource) < 0 .or. islice_selected_source(isource) > NPROC-1) &
+ call exit_MPI(myrank,'something is wrong with the source slice number')
+
+! compute source arrays in source slice
+ if(myrank == islice_selected_source(isource)) then
+ call compute_arrays_source(ispec_selected_source(isource), &
+ xi_source(isource),eta_source(isource),gamma_source(isource),sourcearray, &
+ Mxx(isource),Myy(isource),Mzz(isource),Mxy(isource),Mxz(isource),Myz(isource), &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+ xigll,yigll,zigll,NSPEC_AB)
+ sourcearrays(isource,:,:,:,:) = sourcearray(:,:,:,:)
+ endif
+
+ enddo
+
+ if(myrank == 0) then
+ write(IMAIN,*)
+ write(IMAIN,*) 'Total number of samples for seismograms = ',NSTEP
+ write(IMAIN,*)
+ endif
+
+!--- select local receivers
+
+! count number of receivers located in this slice
+ nrec_local = 0
+
+ do irec = 1,nrec
+
+! check that the receiver slice number is okay
+ if(islice_selected_rec(irec) < 0 .or. islice_selected_rec(irec) > NPROC-1) &
+ call exit_MPI(myrank,'something is wrong with the receiver slice number')
+
+! write info about that receiver
+ if(myrank == islice_selected_rec(irec)) then
+
+ nrec_local = nrec_local + 1
+
+ endif
+
+ enddo
+
+! allocate seismogram array
+ allocate(seismograms_d(NDIM,nrec_local,NSTEP))
+ allocate(seismograms_v(NDIM,nrec_local,NSTEP))
+
+! allocate Lagrange interpolators for receivers
+ allocate(hxir_store(nrec_local,NGLLX))
+ allocate(hetar_store(nrec_local,NGLLY))
+
+! define local to global receiver numbering mapping
+ allocate(number_receiver_global(nrec_local))
+ irec_local = 0
+ do irec = 1,nrec
+ if(myrank == islice_selected_rec(irec)) then
+ irec_local = irec_local + 1
+ number_receiver_global(irec_local) = irec
+ endif
+ enddo
+
+! define and store Lagrange interpolators at all the receivers
+ do irec_local = 1,nrec_local
+ irec = number_receiver_global(irec_local)
+ call lagrange_any(xi_receiver(irec),NGLLX,xigll,hxir,hpxir)
+ call lagrange_any(eta_receiver(irec),NGLLY,yigll,hetar,hpetar)
+ hxir_store(irec_local,:) = hxir(:)
+ hetar_store(irec_local,:) = hetar(:)
+ enddo
+
+! check that the sum of the number of receivers in each slice is nrec
+ call MPI_REDUCE(nrec_local,nrec_tot_found,1,MPI_INTEGER,MPI_SUM,0, &
+ MPI_COMM_WORLD,ier)
+ if(myrank == 0) then
+ write(IMAIN,*)
+ write(IMAIN,*) 'found a total of ',nrec_tot_found,' receivers in all the slices'
+ if(nrec_tot_found /= nrec) then
+ call exit_MPI(myrank,'problem when dispatching the receivers')
+ else
+ write(IMAIN,*) 'this total is okay'
+ endif
+ endif
+
+! initialize seismograms
+ seismograms_d(:,:,:) = 0._CUSTOM_REAL
+ seismograms_v(:,:,:) = 0._CUSTOM_REAL
+
+ if(myrank == 0) then
+
+ if(NSOURCES > 1) write(IMAIN,*) 'Using ',NSOURCES,' point sources'
+
+ write(IMAIN,*)
+ if(ATTENUATION) then
+ write(IMAIN,*) 'incorporating attenuation using ',N_SLS,' standard linear solids'
+ else
+ write(IMAIN,*) 'no attenuation'
+ endif
+ write(IMAIN,*)
+
+ endif
+
+! synchronize all the processes before assembling the mass matrix
+! to make sure all the nodes have finished to read their databases
+ call MPI_BARRIER(MPI_COMM_WORLD,ier)
+
+!jpampuero ABC ==== BEGIN modification of mass matrix for ============================
+!jpampuero ABC IMPLICIT implementation of absorbing boundaries
+!jpampuero ABC M --> M+dt/2*C
+!jpampuero ABC where C is related to the boundary term -C*v
+!jpampuero ABC NOTE: assumes straight edges
+ if(STACEY_ABS_CONDITIONS) then
+
+ ! vn=vx*nx+vy*ny+vz*nz
+ ! tx=rho_vp(i,j,k,ispec)*vn*nx+rho_vs(i,j,k,ispec)*(vx-vn*nx)
+ ! ty=rho_vp(i,j,k,ispec)*vn*ny+rho_vs(i,j,k,ispec)*(vy-vn*ny)
+ ! tz=rho_vp(i,j,k,ispec)*vn*nz+rho_vs(i,j,k,ispec)*(vz-vn*nz)
+
+ ! xmin: (nx,ny,nz) = (-1,0,0)
+ ! vn=-vx
+ ! tx=rho_vp(i,j,k,ispec)*vx
+ ! ty=rho_vs(i,j,k,ispec)*vy
+ ! tz=rho_vs(i,j,k,ispec)*vz
+ do ispec2D=1,nspec2D_xmin
+ ispec=ibelm_xmin(ispec2D)
+ ! exclude elements that are not on absorbing edges
+ if(nkmin_xi(1,ispec2D) == 0 .or. njmin(1,ispec2D) == 0) cycle
+ i=1
+ do k=nkmin_xi(1,ispec2D),NGLLZ
+ do j=njmin(1,ispec2D),njmax(1,ispec2D)
+ iglob=ibool(i,j,k,ispec)
+ weight=jacobian2D_xmin(j,k,ispec2D)*wgllwgll_yz(j,k)
+ rmass(1,iglob)=rmass(1,iglob) + deltatover2*weight*rho_vp(i,j,k,ispec)
+ rmass(2,iglob)=rmass(2,iglob) + deltatover2*weight*rho_vs(i,j,k,ispec)
+ rmass(3,iglob)=rmass(3,iglob) + deltatover2*weight*rho_vs(i,j,k,ispec)
+ enddo
+ enddo
+ enddo
+
+
+ ! xmax: (nx,ny,nz) = (1,0,0)
+ ! vn=vx
+ ! tx=rho_vp(i,j,k,ispec)*vx
+ ! ty=rho_vs(i,j,k,ispec)*vy
+ ! tz=rho_vs(i,j,k,ispec)*vz
+ do ispec2D=1,nspec2D_xmax
+ ispec=ibelm_xmax(ispec2D)
+ ! exclude elements that are not on absorbing edges
+ if(nkmin_xi(2,ispec2D) == 0 .or. njmin(2,ispec2D) == 0) cycle
+ i=NGLLX
+ do k=nkmin_xi(2,ispec2D),NGLLZ
+ do j=njmin(2,ispec2D),njmax(2,ispec2D)
+ iglob=ibool(i,j,k,ispec)
+ weight=jacobian2D_xmax(j,k,ispec2D)*wgllwgll_yz(j,k)
+ rmass(1,iglob)=rmass(1,iglob) + deltatover2*weight*rho_vp(i,j,k,ispec)
+ rmass(2,iglob)=rmass(2,iglob) + deltatover2*weight*rho_vs(i,j,k,ispec)
+ rmass(3,iglob)=rmass(3,iglob) + deltatover2*weight*rho_vs(i,j,k,ispec)
+ enddo
+ enddo
+ enddo
+
+ ! ymax: (nx,ny,nz) = (0,1,0)
+ ! vn=vy
+ ! tx=rho_vs(i,j,k,ispec)*vx
+ ! ty=rho_vp(i,j,k,ispec)*vy
+ ! tz=rho_vs(i,j,k,ispec)*vz
+ do ispec2D=1,nspec2D_ymax
+ ispec=ibelm_ymax(ispec2D)
+ ! exclude elements that are not on absorbing edges
+ if(nkmin_eta(2,ispec2D) == 0 .or. nimin(2,ispec2D) == 0) cycle
+ j=NGLLY
+ do k=nkmin_eta(2,ispec2D),NGLLZ
+ do i=nimin(2,ispec2D),nimax(2,ispec2D)
+ iglob=ibool(i,j,k,ispec)
+ weight=jacobian2D_ymax(i,k,ispec2D)*wgllwgll_xz(i,k)
+ rmass(1,iglob)=rmass(1,iglob) + deltatover2*weight*rho_vs(i,j,k,ispec)
+ rmass(2,iglob)=rmass(2,iglob) + deltatover2*weight*rho_vp(i,j,k,ispec)
+ rmass(3,iglob)=rmass(3,iglob) + deltatover2*weight*rho_vs(i,j,k,ispec)
+ enddo
+ enddo
+ enddo
+
+ ! zmin: (nx,ny,nz) = (0,0,-1)
+ ! vn=-vz
+ ! tx=rho_vs(i,j,k,ispec)*vx
+ ! ty=rho_vs(i,j,k,ispec)*vy
+ ! tz=rho_vp(i,j,k,ispec)*vz
+ do ispec2D=1,NSPEC2D_BOTTOM
+ ispec=ibelm_bottom(ispec2D)
+ k=1
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob=ibool(i,j,k,ispec)
+ weight=jacobian2D_bottom(i,j,ispec2D)*wgllwgll_xy(i,j)
+ rmass(1,iglob)=rmass(1,iglob) + deltatover2*weight*rho_vs(i,j,k,ispec)
+ rmass(2,iglob)=rmass(2,iglob) + deltatover2*weight*rho_vs(i,j,k,ispec)
+ rmass(3,iglob)=rmass(3,iglob) + deltatover2*weight*rho_vp(i,j,k,ispec)
+ enddo
+ enddo
+ enddo
+
+ endif
+!jpampuero ABC ======= END modification of mass matrix ====================================
+
+! the mass matrix needs to be assembled with MPI here once and for all
+ call assemble_MPI_vector(rmass,iproc_xi,iproc_eta,addressing, & !jpampuero ABC
+ iboolleft_xi,iboolright_xi,iboolleft_eta,iboolright_eta, &
+ buffer_send_faces,buffer_received_faces,npoin2D_xi,npoin2D_eta, &
+ NPROC_XI,NPROC_ETA,NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NPOIN2DMAX_XY)
+
+ if(myrank == 0) write(IMAIN,*) 'end assembling MPI mass matrix'
+
+! check that mass matrix is positive
+ if(minval(rmass) <= 0.) call exit_MPI(myrank,'negative mass matrix term')
+
+! for efficiency, invert final mass matrix once and for all in each slice
+ rmass = 1.0 / rmass
+
+! if attenuation is on, shift PREM to right frequency
+! rescale mu in PREM to average frequency for attenuation
+
+ if(ATTENUATION) then
+
+! get and store PREM attenuation model
+ do iattenuation = 1,NUM_REGIONS_ATTENUATION
+
+ call get_attenuation_model(myrank,iattenuation,tau_mu_dble, &
+ tau_sigma_dble,beta_dble,one_minus_sum_beta_dble,factor_scale_dble)
+
+! distinguish whether single or double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ tau_mu(iattenuation,:) = sngl(tau_mu_dble(:))
+ tau_sigma(iattenuation,:) = sngl(tau_sigma_dble(:))
+ beta(iattenuation,:) = sngl(beta_dble(:))
+ factor_scale(iattenuation) = sngl(factor_scale_dble)
+ one_minus_sum_beta(iattenuation) = sngl(one_minus_sum_beta_dble)
+ else
+ tau_mu(iattenuation,:) = tau_mu_dble(:)
+ tau_sigma(iattenuation,:) = tau_sigma_dble(:)
+ beta(iattenuation,:) = beta_dble(:)
+ factor_scale(iattenuation) = factor_scale_dble
+ one_minus_sum_beta(iattenuation) = one_minus_sum_beta_dble
+ endif
+ enddo
+
+! rescale shear modulus according to attenuation model
+
+ do ispec = 1,NSPEC_AB
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ scale_factor = factor_scale(IATTENUATION_UNIFORM)
+ mustore(i,j,k,ispec) = mustore(i,j,k,ispec) * scale_factor
+ enddo
+ enddo
+ enddo
+ enddo
+
+ endif
+
+! initialize arrays to zero
+ displ(:,:) = 0._CUSTOM_REAL
+ veloc(:,:) = 0._CUSTOM_REAL
+ accel(:,:) = 0._CUSTOM_REAL
+
+! put negligible initial value to avoid very slow underflow trapping
+ if(FIX_UNDERFLOW_PROBLEM) displ(:,:) = VERYSMALLVAL
+
+! synchronize all processes to make sure everybody is ready to start time loop
+ call MPI_BARRIER(MPI_COMM_WORLD,ier)
+ if(myrank == 0) write(IMAIN,*) 'All processes are synchronized before time loop'
+
+!
+! s t a r t t i m e i t e r a t i o n s
+!
+
+ if(myrank == 0) then
+ write(IMAIN,*)
+ write(IMAIN,*) ' time step: ',sngl(DT),' s'
+ write(IMAIN,*) 'number of time steps: ',NSTEP
+ write(IMAIN,*) 'total simulated time: ',sngl(NSTEP*DT),' seconds'
+ write(IMAIN,*)
+ endif
+
+! precompute Runge-Kutta coefficients if attenuation
+ if(ATTENUATION) then
+ tauinv(:,:) = - 1. / tau_sigma(:,:)
+ factor_common(:,:) = 2. * beta(:,:) * tauinv(:,:)
+ alphaval(:,:) = 1 + deltat*tauinv(:,:) + deltat**2*tauinv(:,:)**2 / 2. + &
+ deltat**3*tauinv(:,:)**3 / 6. + deltat**4*tauinv(:,:)**4 / 24.
+ betaval(:,:) = deltat / 2. + deltat**2*tauinv(:,:) / 3. + deltat**3*tauinv(:,:)**2 / 8. + deltat**4*tauinv(:,:)**3 / 24.
+ gammaval(:,:) = deltat / 2. + deltat**2*tauinv(:,:) / 6. + deltat**3*tauinv(:,:)**2 / 24.
+ endif
+
+ if(myrank == 0) then
+ write(IMAIN,*)
+ write(IMAIN,*) 'Starting time iteration loop...'
+ write(IMAIN,*)
+ endif
+
+! create an empty file to monitor the start of the simulation
+ if(myrank == 0) then
+ open(unit=IOUT,file='OUTPUT_FILES/starttimeloop.txt',status='unknown')
+ write(IOUT,*) 'starting time loop'
+ close(IOUT)
+ endif
+
+! get MPI starting time
+ time_start = MPI_WTIME()
+
+! clear memory variables if attenuation
+ if(ATTENUATION) then
+ R_xx(:,:,:,:,:) = 0._CUSTOM_REAL
+ R_yy(:,:,:,:,:) = 0._CUSTOM_REAL
+ R_xy(:,:,:,:,:) = 0._CUSTOM_REAL
+ R_xz(:,:,:,:,:) = 0._CUSTOM_REAL
+ R_yz(:,:,:,:,:) = 0._CUSTOM_REAL
+
+ if(FIX_UNDERFLOW_PROBLEM) then
+ R_xx(:,:,:,:,:) = VERYSMALLVAL
+ R_yy(:,:,:,:,:) = VERYSMALLVAL
+ R_xy(:,:,:,:,:) = VERYSMALLVAL
+ R_xz(:,:,:,:,:) = VERYSMALLVAL
+ R_yz(:,:,:,:,:) = VERYSMALLVAL
+ endif
+
+ endif
+
+! *********************************************************
+! ************* MAIN LOOP OVER THE TIME STEPS *************
+! *********************************************************
+
+ do it=1,NSTEP
+
+! compute the maximum of the norm of the displacement
+! in all the slices using an MPI reduction
+! and output timestamp file to check that simulation is running fine
+ if(mod(it,ITAFF_TIME_STEPS) == 0 .or. it == 5) then
+
+! compute maximum of norm of displacement in each slice
+ Usolidnorm = sqrt(maxval(displ(1,:)**2 + displ(2,:)**2 + displ(3,:)**2))
+
+! compute the maximum of the maxima for all the slices using an MPI reduction
+ call MPI_REDUCE(Usolidnorm,Usolidnorm_all,1,CUSTOM_MPI_TYPE,MPI_MAX,0, &
+ MPI_COMM_WORLD,ier)
+
+ if(myrank == 0) then
+
+ write(IMAIN,*) 'Time step # ',it
+ write(IMAIN,*) 'Time: ',sngl((it-1)*DT-hdur(1)),' seconds'
+
+! elapsed time since beginning of the simulation
+ tCPU = MPI_WTIME() - time_start
+ int_tCPU = int(tCPU)
+ ihours = int_tCPU / 3600
+ iminutes = (int_tCPU - 3600*ihours) / 60
+ iseconds = int_tCPU - 3600*ihours - 60*iminutes
+ write(IMAIN,*) 'Elapsed time in seconds = ',tCPU
+ write(IMAIN,"(' Elapsed time in hh:mm:ss = ',i4,' h ',i2.2,' m ',i2.2,' s')") ihours,iminutes,iseconds
+ write(IMAIN,*) 'Mean elapsed time per time step in seconds = ',tCPU/dble(it)
+ write(IMAIN,*) 'Max norm displacement vector U in solid in all slices (m) = ',Usolidnorm_all
+ write(IMAIN,*)
+
+! write time stamp file to give information about progression of simulation
+ write(outputname,"('OUTPUT_FILES/timestamp',i6.6)") it
+ open(unit=IOUT,file=outputname,status='unknown')
+ write(IOUT,*) 'Time step # ',it
+ write(IOUT,*) 'Time: ',sngl((it-1)*DT-hdur(1)),' seconds'
+ write(IOUT,*) 'Elapsed time in seconds = ',tCPU
+ write(IOUT,"(' Elapsed time in hh:mm:ss = ',i4,' h ',i2.2,' m ',i2.2,' s')") ihours,iminutes,iseconds
+ write(IOUT,*) 'Mean elapsed time per time step in seconds = ',tCPU/dble(it)
+ write(IOUT,*) 'Max norm displacement vector U in solid in all slices (m) = ',Usolidnorm_all
+ close(IOUT)
+
+! check stability of the code, exit if unstable
+ if(Usolidnorm_all > STABILITY_THRESHOLD) call exit_MPI(myrank,'code became unstable and blew up')
+
+ endif
+ endif
+
+! update displacement using finite difference time scheme
+ do i=1,NGLOB_AB
+ displ(:,i) = displ(:,i) + deltat*veloc(:,i) + deltatsqover2*accel(:,i)
+ veloc(:,i) = veloc(:,i) + deltatover2*accel(:,i)
+ accel(:,i) = 0._CUSTOM_REAL
+ enddo
+
+ do ispec = 1,NSPEC_AB
+
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+
+ tempx1l = 0.
+ tempx2l = 0.
+ tempx3l = 0.
+
+ tempy1l = 0.
+ tempy2l = 0.
+ tempy3l = 0.
+
+ tempz1l = 0.
+ tempz2l = 0.
+ tempz3l = 0.
+
+ do l=1,NGLLX
+ hp1 = hprime_xx(l,i)
+ iglob = ibool(l,j,k,ispec)
+ tempx1l = tempx1l + displ(1,iglob)*hp1
+ tempy1l = tempy1l + displ(2,iglob)*hp1
+ tempz1l = tempz1l + displ(3,iglob)*hp1
+ enddo
+
+ do l=1,NGLLY
+ hp2 = hprime_yy(l,j)
+ iglob = ibool(i,l,k,ispec)
+ tempx2l = tempx2l + displ(1,iglob)*hp2
+ tempy2l = tempy2l + displ(2,iglob)*hp2
+ tempz2l = tempz2l + displ(3,iglob)*hp2
+ enddo
+
+ do l=1,NGLLZ
+ hp3 = hprime_zz(l,k)
+ iglob = ibool(i,j,l,ispec)
+ tempx3l = tempx3l + displ(1,iglob)*hp3
+ tempy3l = tempy3l + displ(2,iglob)*hp3
+ tempz3l = tempz3l + displ(3,iglob)*hp3
+ enddo
+
+! get derivatives of ux, uy and uz with respect to x, y and z
+!! DK DK NOJACOBIAN xixl = xix(i,j,k,ispec)
+!! DK DK NOJACOBIAN xiyl = xiy(i,j,k,ispec)
+!! DK DK NOJACOBIAN xizl = xiz(i,j,k,ispec)
+!! DK DK NOJACOBIAN etaxl = etax(i,j,k,ispec)
+!! DK DK NOJACOBIAN etayl = etay(i,j,k,ispec)
+!! DK DK NOJACOBIAN etazl = etaz(i,j,k,ispec)
+!! DK DK NOJACOBIAN gammaxl = gammax(i,j,k,ispec)
+!! DK DK NOJACOBIAN gammayl = gammay(i,j,k,ispec)
+!! DK DK NOJACOBIAN gammazl = gammaz(i,j,k,ispec)
+!! DK DK NOJACOBIAN jacobianl = jacobian(i,j,k,ispec)
+ xixl = xix(i,j,k,1)
+ xiyl = xiy(i,j,k,1)
+ xizl = xiz(i,j,k,1)
+ etaxl = etax(i,j,k,1)
+ etayl = etay(i,j,k,1)
+ etazl = etaz(i,j,k,1)
+ gammaxl = gammax(i,j,k,1)
+ gammayl = gammay(i,j,k,1)
+ gammazl = gammaz(i,j,k,1)
+ jacobianl = jacobian(i,j,k,1)
+
+ duxdxl = xixl*tempx1l + etaxl*tempx2l + gammaxl*tempx3l
+ duxdyl = xiyl*tempx1l + etayl*tempx2l + gammayl*tempx3l
+ duxdzl = xizl*tempx1l + etazl*tempx2l + gammazl*tempx3l
+
+ duydxl = xixl*tempy1l + etaxl*tempy2l + gammaxl*tempy3l
+ duydyl = xiyl*tempy1l + etayl*tempy2l + gammayl*tempy3l
+ duydzl = xizl*tempy1l + etazl*tempy2l + gammazl*tempy3l
+
+ duzdxl = xixl*tempz1l + etaxl*tempz2l + gammaxl*tempz3l
+ duzdyl = xiyl*tempz1l + etayl*tempz2l + gammayl*tempz3l
+ duzdzl = xizl*tempz1l + etazl*tempz2l + gammazl*tempz3l
+
+! precompute some sums to save CPU time
+ duxdxl_plus_duydyl = duxdxl + duydyl
+ duxdxl_plus_duzdzl = duxdxl + duzdzl
+ duydyl_plus_duzdzl = duydyl + duzdzl
+ duxdyl_plus_duydxl = duxdyl + duydxl
+ duzdxl_plus_duxdzl = duzdxl + duxdzl
+ duzdyl_plus_duydzl = duzdyl + duydzl
+
+! precompute terms for attenuation if needed
+ if(ATTENUATION) then
+
+! compute deviatoric strain
+ epsilon_trace_over_3 = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+ epsilondev_xx_loc(i,j,k) = duxdxl - epsilon_trace_over_3
+ epsilondev_yy_loc(i,j,k) = duydyl - epsilon_trace_over_3
+ epsilondev_xy_loc(i,j,k) = 0.5 * duxdyl_plus_duydxl
+ epsilondev_xz_loc(i,j,k) = 0.5 * duzdxl_plus_duxdzl
+ epsilondev_yz_loc(i,j,k) = 0.5 * duzdyl_plus_duydzl
+
+! distinguish attenuation factors
+ iselected = IATTENUATION_UNIFORM
+ one_minus_sum_beta_use = one_minus_sum_beta(iselected)
+ minus_sum_beta = one_minus_sum_beta_use - 1.
+
+ endif
+
+ kappal = kappastore(i,j,k,ispec)
+ mul = mustore(i,j,k,ispec)
+
+! use unrelaxed parameters if attenuation
+ if(ATTENUATION) mul = mul * one_minus_sum_beta_use
+
+ lambdalplus2mul = kappal + FOUR_THIRDS * mul
+ lambdal = lambdalplus2mul - 2.*mul
+
+! compute stress sigma
+ sigma_xx = lambdalplus2mul*duxdxl + lambdal*duydyl_plus_duzdzl
+ sigma_yy = lambdalplus2mul*duydyl + lambdal*duxdxl_plus_duzdzl
+ sigma_zz = lambdalplus2mul*duzdzl + lambdal*duxdxl_plus_duydyl
+
+ sigma_xy = mul*duxdyl_plus_duydxl
+ sigma_xz = mul*duzdxl_plus_duxdzl
+ sigma_yz = mul*duzdyl_plus_duydzl
+
+! subtract memory variables if attenuation
+ if(ATTENUATION) then
+ do i_sls = 1,N_SLS
+ R_xx_val = R_xx(i,j,k,ispec,i_sls)
+ R_yy_val = R_yy(i,j,k,ispec,i_sls)
+ sigma_xx = sigma_xx - R_xx_val
+ sigma_yy = sigma_yy - R_yy_val
+ sigma_zz = sigma_zz + R_xx_val + R_yy_val
+ sigma_xy = sigma_xy - R_xy(i,j,k,ispec,i_sls)
+ sigma_xz = sigma_xz - R_xz(i,j,k,ispec,i_sls)
+ sigma_yz = sigma_yz - R_yz(i,j,k,ispec,i_sls)
+ enddo
+ endif
+
+! form dot product with test vector, symmetric form
+ tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_xy*xiyl + sigma_xz*xizl)
+ tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_yz*xizl)
+ tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl)
+
+ tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_xy*etayl + sigma_xz*etazl)
+ tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_yz*etazl)
+ tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl)
+
+ tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_xy*gammayl + sigma_xz*gammazl)
+ tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_yz*gammazl)
+ tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl)
+
+ enddo
+ enddo
+ enddo
+
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+
+ tempx1l = 0.
+ tempy1l = 0.
+ tempz1l = 0.
+
+ tempx2l = 0.
+ tempy2l = 0.
+ tempz2l = 0.
+
+ tempx3l = 0.
+ tempy3l = 0.
+ tempz3l = 0.
+
+ do l=1,NGLLX
+ fac1 = hprimewgll_xx(i,l)
+ tempx1l = tempx1l + tempx1(l,j,k)*fac1
+ tempy1l = tempy1l + tempy1(l,j,k)*fac1
+ tempz1l = tempz1l + tempz1(l,j,k)*fac1
+ enddo
+
+ do l=1,NGLLY
+ fac2 = hprimewgll_yy(j,l)
+ tempx2l = tempx2l + tempx2(i,l,k)*fac2
+ tempy2l = tempy2l + tempy2(i,l,k)*fac2
+ tempz2l = tempz2l + tempz2(i,l,k)*fac2
+ enddo
+
+ do l=1,NGLLZ
+ fac3 = hprimewgll_zz(k,l)
+ tempx3l = tempx3l + tempx3(i,j,l)*fac3
+ tempy3l = tempy3l + tempy3(i,j,l)*fac3
+ tempz3l = tempz3l + tempz3(i,j,l)*fac3
+ enddo
+
+ fac1 = wgllwgll_yz(j,k)
+ fac2 = wgllwgll_xz(i,k)
+ fac3 = wgllwgll_xy(i,j)
+
+! sum contributions from each element to the global mesh and add gravity terms
+
+ iglob = ibool(i,j,k,ispec)
+
+ accel(1,iglob) = accel(1,iglob) - (fac1*tempx1l + fac2*tempx2l + fac3*tempx3l)
+ accel(2,iglob) = accel(2,iglob) - (fac1*tempy1l + fac2*tempy2l + fac3*tempy3l)
+ accel(3,iglob) = accel(3,iglob) - (fac1*tempz1l + fac2*tempz2l + fac3*tempz3l)
+
+! update memory variables based upon the Runge-Kutta scheme
+
+ if(ATTENUATION) then
+
+! use Runge-Kutta scheme to march in time
+ do i_sls = 1,N_SLS
+
+! get coefficients for that standard linear solid
+
+ factor_loc = mustore(i,j,k,ispec) * factor_common(iselected,i_sls)
+ alphaval_loc = alphaval(iselected,i_sls)
+ betaval_loc = betaval(iselected,i_sls)
+ gammaval_loc = gammaval(iselected,i_sls)
+
+! term in xx
+ Sn = factor_loc * epsilondev_xx(i,j,k,ispec)
+ Snp1 = factor_loc * epsilondev_xx_loc(i,j,k)
+ R_xx(i,j,k,ispec,i_sls) = alphaval_loc * R_xx(i,j,k,ispec,i_sls) + betaval_loc * Sn + gammaval_loc * Snp1
+
+! term in yy
+ Sn = factor_loc * epsilondev_yy(i,j,k,ispec)
+ Snp1 = factor_loc * epsilondev_yy_loc(i,j,k)
+ R_yy(i,j,k,ispec,i_sls) = alphaval_loc * R_yy(i,j,k,ispec,i_sls) + betaval_loc * Sn + gammaval_loc * Snp1
+
+! term in zz not computed since zero trace
+
+! term in xy
+ Sn = factor_loc * epsilondev_xy(i,j,k,ispec)
+ Snp1 = factor_loc * epsilondev_xy_loc(i,j,k)
+ R_xy(i,j,k,ispec,i_sls) = alphaval_loc * R_xy(i,j,k,ispec,i_sls) + betaval_loc * Sn + gammaval_loc * Snp1
+
+! term in xz
+ Sn = factor_loc * epsilondev_xz(i,j,k,ispec)
+ Snp1 = factor_loc * epsilondev_xz_loc(i,j,k)
+ R_xz(i,j,k,ispec,i_sls) = alphaval_loc * R_xz(i,j,k,ispec,i_sls) + betaval_loc * Sn + gammaval_loc * Snp1
+
+! term in yz
+ Sn = factor_loc * epsilondev_yz(i,j,k,ispec)
+ Snp1 = factor_loc * epsilondev_yz_loc(i,j,k)
+ R_yz(i,j,k,ispec,i_sls) = alphaval_loc * R_yz(i,j,k,ispec,i_sls) + betaval_loc * Sn + gammaval_loc * Snp1
+
+ enddo ! end of loop on memory variables
+
+ endif ! end attenuation
+
+ enddo
+ enddo
+ enddo
+
+! save deviatoric strain for Runge-Kutta scheme
+ if(ATTENUATION) then
+ epsilondev_xx(:,:,:,ispec) = epsilondev_xx_loc(:,:,:)
+ epsilondev_yy(:,:,:,ispec) = epsilondev_yy_loc(:,:,:)
+ epsilondev_xy(:,:,:,ispec) = epsilondev_xy_loc(:,:,:)
+ epsilondev_xz(:,:,:,ispec) = epsilondev_xz_loc(:,:,:)
+ epsilondev_yz(:,:,:,ispec) = epsilondev_yz_loc(:,:,:)
+ endif
+
+ enddo ! spectral element loop
+
+!jpampuero *********** BEGIN ABSORBING BOUNDARY CONDITIONS ******************
+! add Stacey conditions
+
+ if(STACEY_ABS_CONDITIONS) then
+
+! xmin
+ do ispec2D=1,nspec2D_xmin
+
+ ispec=ibelm_xmin(ispec2D)
+
+! exclude elements that are not on absorbing edges
+ if(nkmin_xi(1,ispec2D) == 0 .or. njmin(1,ispec2D) == 0) cycle
+
+ i=1
+ do k=nkmin_xi(1,ispec2D),NGLLZ
+ do j=njmin(1,ispec2D),njmax(1,ispec2D)
+ iglob=ibool(i,j,k,ispec)
+
+ vx=veloc(1,iglob)
+ vy=veloc(2,iglob)
+ vz=veloc(3,iglob)
+
+ nx=normal_xmin(1,j,k,ispec2D)
+ ny=normal_xmin(2,j,k,ispec2D)
+ nz=normal_xmin(3,j,k,ispec2D)
+
+ vn=vx*nx+vy*ny+vz*nz
+
+ tx=rho_vp(i,j,k,ispec)*vn*nx+rho_vs(i,j,k,ispec)*(vx-vn*nx)
+ ty=rho_vp(i,j,k,ispec)*vn*ny+rho_vs(i,j,k,ispec)*(vy-vn*ny)
+ tz=rho_vp(i,j,k,ispec)*vn*nz+rho_vs(i,j,k,ispec)*(vz-vn*nz)
+
+ weight=jacobian2D_xmin(j,k,ispec2D)*wgllwgll_yz(j,k)
+
+ accel(1,iglob)=accel(1,iglob) - tx*weight
+ accel(2,iglob)=accel(2,iglob) - ty*weight
+ accel(3,iglob)=accel(3,iglob) - tz*weight
+
+ enddo
+ enddo
+ enddo
+
+! xmax
+ do ispec2D=1,nspec2D_xmax
+ ispec=ibelm_xmax(ispec2D)
+
+! exclude elements that are not on absorbing edges
+ if(nkmin_xi(2,ispec2D) == 0 .or. njmin(2,ispec2D) == 0) cycle
+
+ i=NGLLX
+ do k=nkmin_xi(2,ispec2D),NGLLZ
+ do j=njmin(2,ispec2D),njmax(2,ispec2D)
+ iglob=ibool(i,j,k,ispec)
+
+ vx=veloc(1,iglob)
+ vy=veloc(2,iglob)
+ vz=veloc(3,iglob)
+
+ nx=normal_xmax(1,j,k,ispec2D)
+ ny=normal_xmax(2,j,k,ispec2D)
+ nz=normal_xmax(3,j,k,ispec2D)
+
+ vn=vx*nx+vy*ny+vz*nz
+
+ tx=rho_vp(i,j,k,ispec)*vn*nx+rho_vs(i,j,k,ispec)*(vx-vn*nx)
+ ty=rho_vp(i,j,k,ispec)*vn*ny+rho_vs(i,j,k,ispec)*(vy-vn*ny)
+ tz=rho_vp(i,j,k,ispec)*vn*nz+rho_vs(i,j,k,ispec)*(vz-vn*nz)
+
+ weight=jacobian2D_xmax(j,k,ispec2D)*wgllwgll_yz(j,k)
+
+ accel(1,iglob)=accel(1,iglob) - tx*weight
+ accel(2,iglob)=accel(2,iglob) - ty*weight
+ accel(3,iglob)=accel(3,iglob) - tz*weight
+
+ enddo
+ enddo
+ enddo
+
+!jpampuero ================ BEGIN DYNAMIC FAULT =========================
+!jpampuero I am implementing the FAULT boundary conditions in place of one
+!jpampuero of the absorbing boundaries.
+!jpampuero It is a quick and dirty implementation, for SCEC 1st benchmark test.
+!jpampuero No split nodes need to be implemented.
+!jpampuero The model is symmetric.
+!jpampuero Assume only ONE proc contains the fault: NPROC_XI=1
+!jpampuero Fault located at y=ymin
+!jpampuero
+ if (FaultInThisProc) then
+ do ispec2D=1,nspec2D_ymin
+
+ ! exclude elements that are not on absorbing edges
+ if(nkmin_eta(1,ispec2D) == 0 .or. nimin(1,ispec2D) == 0) cycle
+
+ !ispec = bulk element associated to the ispec2D'th boundary element
+ ispec=ibelm_ymin(ispec2D)
+
+ !skip element if its center is not in rupture region
+ ! iglob=ibool(3,1,3,ispec)
+ ! if (-zstore(iglob)>RuptureDepth .OR. abs(xstore(iglob))>RuptureHalfLength) cycle
+
+ j=1
+ do k=nkmin_eta(1,ispec2D),NGLLZ
+ do i=nimin(1,ispec2D),nimax(1,ispec2D)
+! do k=1,NGLLZ
+! do i=1,NGLLX
+ iglob=ibool(i,j,k,ispec)
+
+ ! Update strength :
+ ! For symmetry reasons slip=2*displ
+ FaultStrength(i,k,ispec2D) = max(FaultMuD(i,k,ispec2D), &
+ FaultMuS(i,k,ispec2D)-FaultSWR(i,k,ispec2D)*sqrt(displ(1,iglob)**2+displ(3,iglob)**2) )
+
+ ! From Newmark: v_new = v_pred +dt/2*a_new
+ ! so: a_new = (v_new-v_pred)*2/dt = M^(-1) *( -K*d_new -B*t_new )
+
+ !velocity assuming the fault is stress free during this timestep (t_new=0)
+ !v_free = v_pred - dt/2 *M^(-1) *K*d_new
+ FaultDVxFree = veloc(1,iglob)+deltatover2*rmass(1,iglob)*accel(1,iglob)
+ FaultDVzFree = veloc(3,iglob)+deltatover2*rmass(3,iglob)*accel(3,iglob)
+ !NOTE: I assume that only ONE proc contains the fault plane
+ ! else accel=-K*d needs to be MPI assembled prior to compute V_free
+
+ !"stick" or trial stress, assuming v_new = 0 during this timestep
+ !T_trial = B^(-1) * M * 2/dt*V_free
+ weight=FaultB(i,k,ispec2D) ! assembled
+ tx = FaultDVxFree/(deltatover2*weight*rmass(1,iglob))
+ tz = FaultDVzFree/(deltatover2*weight*rmass(3,iglob))
+
+! !VECTORIAL FRICTION VERSION
+! !test yield on magnitude of absolute stress (+initial)
+! txa = tx+FaultTxInit(i,k,ispec2D)
+! tza = tz+FaultTzInit(i,k,ispec2D)
+! ta = sqrt( txa*txa + tza*tza )
+! if ( ta>FaultStrength(i,k,ispec2D) ) then
+! va = sqrt(FaultDVxFree**2+ FaultDVzFree**2)
+! if ( va<VERYSMALLVAL ) then
+! dirx = txa/ta
+! dirz = tza/ta
+! else
+! dirx = FaultDVxFree/va
+! dirz = FaultDVzFree/va
+! endif
+! ta = FaultStrength(i,k,ispec2D)
+! tx = ta*dirx -FaultTxInit(i,k,ispec2D)
+! tz = ta*dirz -FaultTzInit(i,k,ispec2D)
+! endif
+
+! !SCALAR FRICTION VERSION
+ tx = min(tx, FaultStrength(i,k,ispec2D) -FaultTxInit(i,k,ispec2D))
+
+ FaultTx(i,k,ispec2D) = tx !store relative shear stress
+ FaultTz(i,k,ispec2D) = tz !store relative shear stress
+
+
+ enddo
+ enddo
+ enddo
+
+ do ispec2D=1,nspec2D_ymin
+ if(nkmin_eta(1,ispec2D) == 0 .or. nimin(1,ispec2D) == 0) cycle
+ ispec=ibelm_ymin(ispec2D)
+ j=1
+ do k=nkmin_eta(1,ispec2D),NGLLZ
+ do i=nimin(1,ispec2D),nimax(1,ispec2D)
+ iglob=ibool(i,j,k,ispec)
+ weight=jacobian2D_ymin(i,k,ispec2D)*wgllwgll_xz(i,k) ! not assembled
+ accel(1,iglob)=accel(1,iglob) - FaultTx(i,k,ispec2D)*weight
+ !accel(2,iglob): Stress-free on the Y direction (normal to the fault plane): do nothing
+! accel(3,iglob)=accel(3,iglob) - FaultTz(i,k,ispec2D)*weight ! VECTORIAL friction
+ accel(3,iglob)=0._CUSTOM_REAL ! if SCALAR friction
+ enddo
+ enddo
+ enddo
+
+ endif
+
+!jpampuero ================ END DYNAMIC FAULT =========================
+
+! ymax
+ do ispec2D=1,nspec2D_ymax
+
+ ispec=ibelm_ymax(ispec2D)
+
+! exclude elements that are not on absorbing edges
+ if(nkmin_eta(2,ispec2D) == 0 .or. nimin(2,ispec2D) == 0) cycle
+
+ j=NGLLY
+ do k=nkmin_eta(2,ispec2D),NGLLZ
+ do i=nimin(2,ispec2D),nimax(2,ispec2D)
+ iglob=ibool(i,j,k,ispec)
+
+ vx=veloc(1,iglob)
+ vy=veloc(2,iglob)
+ vz=veloc(3,iglob)
+
+ nx=normal_ymax(1,i,k,ispec2D)
+ ny=normal_ymax(2,i,k,ispec2D)
+ nz=normal_ymax(3,i,k,ispec2D)
+
+ vn=vx*nx+vy*ny+vz*nz
+
+ tx=rho_vp(i,j,k,ispec)*vn*nx+rho_vs(i,j,k,ispec)*(vx-vn*nx)
+ ty=rho_vp(i,j,k,ispec)*vn*ny+rho_vs(i,j,k,ispec)*(vy-vn*ny)
+ tz=rho_vp(i,j,k,ispec)*vn*nz+rho_vs(i,j,k,ispec)*(vz-vn*nz)
+
+ weight=jacobian2D_ymax(i,k,ispec2D)*wgllwgll_xz(i,k)
+
+ accel(1,iglob)=accel(1,iglob) - tx*weight
+ accel(2,iglob)=accel(2,iglob) - ty*weight
+ accel(3,iglob)=accel(3,iglob) - tz*weight
+
+ enddo
+ enddo
+ enddo
+
+
+! bottom (zmin)
+ do ispec2D=1,NSPEC2D_BOTTOM
+
+ ispec=ibelm_bottom(ispec2D)
+
+ k=1
+ do j=1,NGLLY
+ do i=1,NGLLX
+
+ iglob=ibool(i,j,k,ispec)
+
+ vx=veloc(1,iglob)
+ vy=veloc(2,iglob)
+ vz=veloc(3,iglob)
+
+ nx=normal_bottom(1,i,j,ispec2D)
+ ny=normal_bottom(2,i,j,ispec2D)
+ nz=normal_bottom(3,i,j,ispec2D)
+
+ vn=vx*nx+vy*ny+vz*nz
+
+ tx=rho_vp(i,j,k,ispec)*vn*nx+rho_vs(i,j,k,ispec)*(vx-vn*nx)
+ ty=rho_vp(i,j,k,ispec)*vn*ny+rho_vs(i,j,k,ispec)*(vy-vn*ny)
+ tz=rho_vp(i,j,k,ispec)*vn*nz+rho_vs(i,j,k,ispec)*(vz-vn*nz)
+
+ weight=jacobian2D_bottom(i,j,ispec2D)*wgllwgll_xy(i,j)
+
+ accel(1,iglob)=accel(1,iglob) - tx*weight
+ accel(2,iglob)=accel(2,iglob) - ty*weight
+ accel(3,iglob)=accel(3,iglob) - tz*weight
+
+ enddo
+ enddo
+ enddo
+
+ endif ! end of Stacey conditions
+!jpampuero ************* END OF ABSORBING BOUNDARY CONDITIONS *********************
+
+ do isource = 1,NSOURCES
+
+! add the source (only if this proc carries the source)
+ if(myrank == islice_selected_source(isource)) then
+
+ stf = comp_source_time_function(dble(it-1)*DT-hdur(isource)-t_cmt(isource),hdur(isource))
+
+! distinguish whether single or double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ stf_used = sngl(stf)
+ else
+ stf_used = stf
+ endif
+
+! add source array
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+ iglob = ibool(i,j,k,ispec_selected_source(isource))
+ accel(:,iglob) = accel(:,iglob) + sourcearrays(isource,:,i,j,k)*stf_used
+ enddo
+ enddo
+ enddo
+
+ endif
+
+ enddo
+
+! assemble all the contributions between slices using MPI
+ call assemble_MPI_vector(accel,iproc_xi,iproc_eta,addressing, &
+ iboolleft_xi,iboolright_xi,iboolleft_eta,iboolright_eta, &
+ buffer_send_faces,buffer_received_faces,npoin2D_xi,npoin2D_eta, &
+ NPROC_XI,NPROC_ETA,NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NPOIN2DMAX_XY)
+
+ do i=1,NGLOB_AB
+ accel(1,i) = accel(1,i)*rmass(1,i) !jpampuero ABC
+ accel(2,i) = accel(2,i)*rmass(2,i) !jpampuero ABC
+ accel(3,i) = accel(3,i)*rmass(3,i) !jpampuero ABC
+
+ veloc(:,i) = veloc(:,i) + deltatover2*accel(:,i)
+ enddo
+
+! write the seismograms with time shift
+
+ do irec_local = 1,nrec_local
+
+! get global number of that receiver
+ irec = number_receiver_global(irec_local)
+
+! perform the general interpolation using Lagrange polynomials
+ uxd = ZERO
+ uyd = ZERO
+ uzd = ZERO
+
+ vxd = ZERO
+ vyd = ZERO
+ vzd = ZERO
+
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+
+! receivers are always located at the surface in the crustal region of the mesh
+ iglob = ibool(i,j,NGLLZ,ispec_selected_rec(irec))
+
+ hlagrange = hxir_store(irec_local,i)*hetar_store(irec_local,j)
+
+! save displacement
+ uxd = uxd + dble(displ(1,iglob))*hlagrange
+ uyd = uyd + dble(displ(2,iglob))*hlagrange
+ uzd = uzd + dble(displ(3,iglob))*hlagrange
+
+! save velocity
+ vxd = vxd + dble(veloc(1,iglob))*hlagrange
+ vyd = vyd + dble(veloc(2,iglob))*hlagrange
+ vzd = vzd + dble(veloc(3,iglob))*hlagrange
+
+ enddo
+ enddo
+
+! store North, East and Vertical components
+
+! distinguish whether single or double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ seismograms_d(:,irec_local,it) = sngl((nu(:,1,irec)*uxd + nu(:,2,irec)*uyd + nu(:,3,irec)*uzd))
+ seismograms_v(:,irec_local,it) = sngl((nu(:,1,irec)*vxd + nu(:,2,irec)*vyd + nu(:,3,irec)*vzd))
+ else
+ seismograms_d(:,irec_local,it) = (nu(:,1,irec)*uxd + nu(:,2,irec)*uyd + nu(:,3,irec)*uzd)
+ seismograms_v(:,irec_local,it) = (nu(:,1,irec)*vxd + nu(:,2,irec)*vyd + nu(:,3,irec)*vzd)
+ endif
+
+ enddo
+
+! write the current seismograms
+ if(mod(it,NSEIS) == 0) then
+ call write_seismograms_d(myrank,seismograms_d,number_receiver_global,station_name, &
+ network_name,nrec,nrec_local,it,DT,NSTEP,hdur(1),LOCAL_PATH)
+ !network_name,nrec,nrec_local,it,DT,NSTEP,hdur(1),'OUTPUT_FILES') !jpampuero
+ call write_seismograms_v(myrank,seismograms_v,number_receiver_global,station_name, &
+ network_name,nrec,nrec_local,it,DT,NSTEP,hdur(1),LOCAL_PATH)
+ !network_name,nrec,nrec_local,it,DT,NSTEP,hdur(1),'OUTPUT_FILES') !jpampuero
+ endif
+
+!jpampuero =========== BEGIN OUTPUT FAULT SNAPSHOTS ================
+ if (FaultInThisProc) then
+
+ !OUTPUT: time series at selected fault points
+ FaultDisplX(:,it) = 2.0*displ(1,FaultIglobOut)
+ FaultDisplZ(:,it) = 2.0*displ(3,FaultIglobOut)
+ FaultVelocX(:,it) = 2.0*veloc(1,FaultIglobOut)
+ FaultVelocZ(:,it) = 2.0*veloc(3,FaultIglobOut)
+ do i=1,FaultOutputNPTS
+ FaultStressX(i,it) = FaultTx(FaultIgllOut(i),FaultKgllOut(i),FaultElmtOut(i))
+ FaultStressZ(i,it) = FaultTz(FaultIgllOut(i),FaultKgllOut(i),FaultElmtOut(i))
+ enddo
+
+ !OUTPUT: slip rate snapshot
+ if (mod(it,FaultSnapshotNDT) == 0) then
+ !write(FaultOutputName,'(A,'/Snapshot',I2.2,'_',I3.3,'.bin')') LOCAL_PATH,it/FaultSnapshotNDT,myrank
+ write(FaultOutputName,"('OUTPUT_FILES/Snapshot',I0,'.bin')") it/FaultSnapshotNDT
+ open(unit=IOUT, file= trim(FaultOutputName), status='replace', form='unformatted')
+ do ispec2D=1,nspec2D_ymin
+ ispec=ibelm_ymin(ispec2D)
+ j=1
+ !skip element if its center is not in rupture region
+ iglob=ibool(3,j,3,ispec)
+ if (-zstore(iglob)>RuptureDepth .OR. abs(xstore(iglob))>RuptureHalfLength) cycle
+
+ do k=1,NGLLZ-1,2 ! watch out: do not repeat a node
+ do i=1,NGLLX-1,2 ! assuming NGLL=5
+ iglob=ibool(i,j,k,ispec)
+ write(IOUT) xstore(iglob),zstore(iglob) &
+ ,2.*veloc(1,iglob),2.*veloc(3,iglob) &
+ ,2.*displ(1,iglob),2.*displ(3,iglob) &
+ ,FaultTx(i,k,ispec2D),FaultTz(i,k,ispec2D)
+ enddo
+ enddo
+ enddo
+ close(IOUT)
+
+ !intermediate output at selected fault points
+ do i=1,FaultOutputNPTS
+ open(IOUT,file='OUTPUT_FILES/'//FaultOutputPtsName(i),status='replace')
+ write(IOUT,*) "# author=Jean-Paul Ampuero (am)"
+ write(IOUT,*) "# date=2004/09/08"
+ write(IOUT,*) "# code=SPECFEM3D (se)"
+ write(IOUT,*) "# code_version=0904"
+ write(IOUT,*) "# element_size=333.33 m (5 nodes per element edge)"
+ write(IOUT,*) "# time_step=",DT
+ write(IOUT,*) "# num_time_steps=",NSTEP
+ write(IOUT,*) "# location=",FaultOutputPtsName(i)
+ write(IOUT,*) "# Time series in 7 column of E14.6"
+ write(IOUT,*) "# Column #1 = Time (s)"
+ write(IOUT,*) "# Column #2 = horizontal right-lateral slip (m)"
+ write(IOUT,*) "# Column #3 = horizontal right-lateral slip rate (m/s)"
+ write(IOUT,*) "# Column #4 = horizontal right-lateral shear stress (MPa)"
+ write(IOUT,*) "# Column #5 = vertical up-dip slip (m)"
+ write(IOUT,*) "# Column #6 = vertical up-dip slip rate (m/s)"
+ write(IOUT,*) "# Column #7 = vertical up-dip shear stress (MPa)"
+ do k=1,NSTEP
+ write(IOUT,'(7(E14.6))') k*DT,FaultDisplX(i,k),FaultVelocX(i,k),FaultStressX(i,k)/1e6 &
+ ,FaultDisplZ(i,k),FaultVelocZ(i,k),FaultStressZ(i,k)/1e6
+ enddo
+ close(IOUT)
+ enddo
+
+ endif
+
+ endif
+!jpampuero =========== END OUTPUT FAULT SNAPSHOTS ================
+
+!
+!---- end of time iteration loop
+!
+ enddo ! end of main time loop
+
+! write the final seismograms
+ call write_seismograms_d(myrank,seismograms_d,number_receiver_global,station_name, &
+ network_name,nrec,nrec_local,it,DT,NSTEP,hdur(1),LOCAL_PATH)
+ !network_name,nrec,nrec_local,it,DT,NSTEP,hdur(1),'OUTPUT_FILES') !jpampuero
+ call write_seismograms_v(myrank,seismograms_v,number_receiver_global,station_name, &
+ network_name,nrec,nrec_local,it,DT,NSTEP,hdur(1),LOCAL_PATH)
+ !network_name,nrec,nrec_local,it,DT,NSTEP,hdur(1),'OUTPUT_FILES') !jpampuero
+
+!jpampuero === BEGIN OUTPUT: time series at selected fault points ===
+ if (FaultInThisProc) then
+ FaultStressX = FaultStressX/1e6
+ FaultStressZ = FaultStressZ/1e6
+ do i=1,FaultOutputNPTS
+ open(IOUT,file='OUTPUT_FILES/'//FaultOutputPtsName(i),status='replace')
+ write(IOUT,*) "# author=Jean-Paul Ampuero (am)"
+ write(IOUT,*) "# date=2004/09/08"
+ write(IOUT,*) "# code=SPECFEM3D (se)"
+ write(IOUT,*) "# code_version=0904"
+ write(IOUT,*) "# element_size=333.33 m (5 nodes per element edge)"
+ write(IOUT,*) "# time_step=",DT
+ write(IOUT,*) "# num_time_steps=",NSTEP
+ write(IOUT,*) "# location=",FaultOutputPtsName(i)
+ write(IOUT,*) "# Time series in 7 column of E14.6"
+ write(IOUT,*) "# Column #1 = Time (s)"
+ write(IOUT,*) "# Column #2 = horizontal right-lateral slip (m)"
+ write(IOUT,*) "# Column #3 = horizontal right-lateral slip rate (m/s)"
+ write(IOUT,*) "# Column #4 = horizontal right-lateral shear stress (MPa)"
+ write(IOUT,*) "# Column #5 = vertical up-dip slip (m)"
+ write(IOUT,*) "# Column #6 = vertical up-dip slip rate (m/s)"
+ write(IOUT,*) "# Column #7 = vertical up-dip shear stress (MPa)"
+ do it=1,NSTEP
+ write(IOUT,'(7(E14.6))') it*DT,FaultDisplX(i,it),FaultVelocX(i,it),FaultStressX(i,it) &
+ ,FaultDisplZ(i,it),FaultVelocZ(i,it),FaultStressZ(i,it)
+ enddo
+ close(IOUT)
+ enddo
+! open(unit=IOUT,file='OUTPUT_FILES/FaultDispl.bin',status='replace',form='unformatted')
+! write(IOUT) FaultDispl
+! close(IOUT)
+! open(unit=IOUT,file='OUTPUT_FILES/FaultVeloc.bin',status='replace',form='unformatted')
+! write(IOUT) FaultVeloc
+! close(IOUT)
+ endif
+!jpampuero === END OUTPUT: time series at selected fault points ===
+
+! close the main output file
+ if(myrank == 0) then
+ write(IMAIN,*)
+ write(IMAIN,*) 'End of the simulation'
+ write(IMAIN,*)
+ close(IMAIN)
+ endif
+
+! synchronize all the processes to make sure everybody has finished
+ call MPI_BARRIER(MPI_COMM_WORLD,ier)
+
+! stop all the MPI processes, and exit
+ call MPI_FINALIZE(ier)
+
+ end program specfem3D
+
Deleted: seismo/3D/SPECFEM3D/trunk/utils/unused_routines/ampuero_implicit_Clayton_now_done_fixed_in_our_official_SVN_code/ampuero_implicit_ABC_specfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D/trunk/utils/unused_routines/ampuero_implicit_Clayton/ampuero_implicit_ABC_specfem3D.f90 2012-05-20 02:25:14 UTC (rev 20195)
+++ seismo/3D/SPECFEM3D/trunk/utils/unused_routines/ampuero_implicit_Clayton_now_done_fixed_in_our_official_SVN_code/ampuero_implicit_ABC_specfem3D.f90 2012-05-20 13:54:06 UTC (rev 20196)
@@ -1,1853 +0,0 @@
- program specfem3D
-!=====================================================================
-!
-! S p e c f e m 3 D B a s i n V e r s i o n 1 . 1
-! --------------------------------------------------
-!
-! Main authors: Dimitri Komatitsch and Jeroen Tromp
-! Princeton University, USA and University of Pau / CNRS / INRIA
-! (c) California Institute of Technology October 2002
-!
-! 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.
-!
-!=====================================================================
-!
-! United States Government Sponsorship Acknowledged.
-!
-!jpampuero Mon Oct 20 23:09:00 EDT 2003
-!jpampuero I am implementing the FAULT boundary conditions in place of one
-!jpampuero of the absorbing boundaries.
-!jpampuero It is a quick and dirty implementation, for SCEC 1st benchmark test.
-!jpampuero No split nodes need to be implemented.
-!jpampuero The model is restricted to be symetric.
-!jpampuero
-!jpampuero All lines containing "!jpampuero" were touched.
-!jpampuero Lines containing "!!jpampuero" were only commented out.
-!jpampuero Most of my variables names have ThisFormat
-
- implicit none
-
-! standard include of the MPI library
- include 'mpif.h'
-
- include "constants.h"
- include "precision.h"
-
-! include values created by the mesher
- include "OUTPUT_FILES/values_from_mesher.h"
-
-!jpampuero Fault parameters
- include "fault_parameters.h" !jpampuero
-
-!=======================================================================!
-! !
-! specfem3D is a 3-D spectral-element solver for a basin. !
-! It uses a mesh generated by program meshfem3D !
-! !
-!=======================================================================!
-!
-! If you use this code for your own research, please send an email
-! to Jeroen Tromp <jtromp at gps.caltech.edu> for information, and cite:
-!
-! @article{KoTr99,
-! author={D. Komatitsch and J. Tromp},
-! year=1999,
-! title={Introduction to the spectral-element method for 3-{D} seismic wave propagation},
-! journal={Geophys. J. Int.},
-! volume=139,
-! pages={806-822}}
-!
-! Evolution of the code:
-! ---------------------
-!
-! MPI v. 1.1 Dimitri Komatitsch, Caltech, October 2002: Zhu's Moho map, scaling
-! of Vs with depth, Hauksson's regional model, attenuation, oceans, movies
-! MPI v. 1.0 Dimitri Komatitsch, Caltech, May 2002: first MPI version
-! based on global code
-
-! memory variables and standard linear solids for attenuation
- double precision, dimension(N_SLS) :: tau_mu_dble,tau_sigma_dble,beta_dble
- double precision factor_scale_dble,one_minus_sum_beta_dble
- real(kind=CUSTOM_REAL), dimension(NUM_REGIONS_ATTENUATION,N_SLS) :: tau_mu,tau_sigma,beta
- real(kind=CUSTOM_REAL), dimension(NUM_REGIONS_ATTENUATION) :: factor_scale,one_minus_sum_beta
-
- real(kind=CUSTOM_REAL), dimension(NUM_REGIONS_ATTENUATION,N_SLS) :: tauinv,alphaval,betaval,gammaval,factor_common
- integer iattenuation
- double precision scale_factor
-
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION,N_SLS) :: &
- R_xx,R_yy,R_xy,R_xz,R_yz
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION) :: &
- epsilondev_xx,epsilondev_yy,epsilondev_xy,epsilondev_xz,epsilondev_yz
-
- integer NPOIN2DMAX_XY
-
- integer, dimension(:), allocatable :: ibelm_xmin,ibelm_xmax, &
- ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top
- real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: &
- jacobian2D_xmin,jacobian2D_xmax,jacobian2D_ymin,jacobian2D_ymax, &
- jacobian2D_bottom,jacobian2D_top
- real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: normal_xmin, &
- normal_xmax,normal_ymin,normal_ymax, &
- normal_bottom,normal_top
-
-! buffers for send and receive between faces of the slices and the chunks
- real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: buffer_send_faces,buffer_received_faces
-
-! -----------------
-
-! mesh parameters
- integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
-
-!! DK DK NOJACOBIAN real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,1) :: &
- xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,jacobian
- real(kind=CUSTOM_REAL), dimension(NGLOB_AB) :: xstore,ystore,zstore
-
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
- kappastore,mustore
-
-! Stacey
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
- rho_vp,rho_vs
-
-! local to global mapping
- integer, dimension(NSPEC_AB) :: idoubling
-
-! mass matrix
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: rmass !jpampuero ABC
-
-! displacement, velocity, acceleration
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displ,veloc,accel
-
- real(kind=CUSTOM_REAL) xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
- real(kind=CUSTOM_REAL) duxdxl,duxdyl,duxdzl,duydxl,duydyl,duydzl,duzdxl,duzdyl,duzdzl
-
- real(kind=CUSTOM_REAL) duxdxl_plus_duydyl,duxdxl_plus_duzdzl,duydyl_plus_duzdzl
- real(kind=CUSTOM_REAL) duxdyl_plus_duydxl,duzdxl_plus_duxdzl,duzdyl_plus_duydzl
-
- real(kind=CUSTOM_REAL) sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz
-
- real(kind=CUSTOM_REAL) hp1,hp2,hp3
- real(kind=CUSTOM_REAL) fac1,fac2,fac3
- real(kind=CUSTOM_REAL) lambdal,kappal,mul,lambdalplus2mul
-
- real(kind=CUSTOM_REAL) tempx1l,tempx2l,tempx3l
- real(kind=CUSTOM_REAL) tempy1l,tempy2l,tempy3l
- real(kind=CUSTOM_REAL) tempz1l,tempz2l,tempz3l
-
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
- tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3
-
-! for attenuation
- real(kind=CUSTOM_REAL) R_xx_val,R_yy_val
- real(kind=CUSTOM_REAL) factor_loc,alphaval_loc,betaval_loc,gammaval_loc,Sn,Snp1
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: epsilondev_xx_loc, &
- epsilondev_yy_loc, epsilondev_xy_loc, epsilondev_xz_loc, epsilondev_yz_loc
- real(kind=CUSTOM_REAL) epsilon_trace_over_3
-
- integer l
- integer i_SLS
- real(kind=CUSTOM_REAL) one_minus_sum_beta_use,minus_sum_beta
- integer iselected
-
-! --------
-
-! parameters for the source
- integer it,isource
- integer, dimension(:), allocatable :: islice_selected_source,ispec_selected_source
- integer yr,jda,ho,mi
- real(kind=CUSTOM_REAL) stf_used
- real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: sourcearray
- real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: sourcearrays
- double precision sec,stf
- double precision, dimension(:), allocatable :: Mxx,Myy,Mzz,Mxy,Mxz,Myz
- double precision, dimension(:), allocatable :: xi_source,eta_source,gamma_source
- double precision, dimension(:), allocatable :: t_cmt,hdur
- double precision, dimension(:), allocatable :: utm_x_source,utm_y_source
- double precision, external :: comp_source_time_function
-
-
-! time scheme
- real(kind=CUSTOM_REAL) deltat,deltatover2,deltatsqover2
-
-! receiver information
- integer nrec,nrec_local,nrec_tot_found
- integer irec_local
- integer, allocatable, dimension(:) :: islice_selected_rec,ispec_selected_rec,number_receiver_global
- double precision, allocatable, dimension(:) :: xi_receiver,eta_receiver
- double precision hlagrange
-
-! timing information for the stations
- double precision, allocatable, dimension(:,:,:) :: nu
- character(len=8), allocatable, dimension(:) :: station_name,network_name
-
-! seismograms
- double precision uxd,uyd,uzd
- double precision vxd,vyd,vzd
- real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: seismograms_d,seismograms_v
-
- integer i,j,k,ispec,irec,iglob
-
-! Gauss-Lobatto-Legendre points of integration and weights
- double precision, dimension(NGLLX) :: xigll,wxgll
- double precision, dimension(NGLLY) :: yigll,wygll
- double precision, dimension(NGLLZ) :: zigll,wzgll
-
-! array with derivatives of Lagrange polynomials and precalculated products
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx,hprimewgll_xx
- real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLY) :: hprime_yy,hprimewgll_yy
- real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ) :: hprime_zz,hprimewgll_zz
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY) :: wgllwgll_xy
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ) :: wgllwgll_xz
- real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLZ) :: wgllwgll_yz
-
-! Lagrange interpolators at receivers
- double precision, dimension(:), allocatable :: hxir,hetar,hpxir,hpetar
- double precision, dimension(:,:), allocatable :: hxir_store,hetar_store
-
-! 2-D addressing and buffers for summation between slices
- integer, dimension(:), allocatable :: iboolleft_xi, &
- iboolright_xi,iboolleft_eta,iboolright_eta
-
-! for addressing of the slices
- integer, dimension(:,:), allocatable :: addressing
- integer, dimension(:), allocatable :: iproc_xi_slice,iproc_eta_slice
-
-! proc numbers for MPI
- integer myrank,sizeprocs,ier
-
- integer npoin2D_xi,npoin2D_eta
-
- integer iproc_xi,iproc_eta,iproc,iproc_read
-
-! maximum of the norm of the displacement
- real(kind=CUSTOM_REAL) Usolidnorm,Usolidnorm_all
-
-! timer MPI
- integer ihours,iminutes,iseconds,int_tCPU
- double precision time_start,tCPU
-
-! parameters read from parameter file
- integer NER_TOP_LAYER,NER_MIDDLE_LAYER,NER_BOTTOM_HALFSPACE,NEX_ETA,NEX_XI, &
- NPROC_ETA,NPROC_XI,NSEIS,NSTEP,UTM_PROJECTION_ZONE
- integer NSOURCES
-
- double precision UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK
- double precision DT,LAT_MIN,LAT_MAX,LONG_MIN,LONG_MAX
- logical ATTENUATION,STACEY_ABS_CONDITIONS
-
- character(len=256) LOCAL_PATH,prname
-
-! parameters deduced from parameters read from file
- integer NPROC,NEX_PER_PROC_XI,NEX_PER_PROC_ETA
- integer NER
-
- integer NSPEC2D_A_XI,NSPEC2D_B_XI, &
- NSPEC2D_A_ETA,NSPEC2D_B_ETA, &
- NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX, &
- NSPEC2D_BOTTOM,NSPEC2D_TOP, &
- NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX
-
-! names of the data files for all the processors in MPI
- character(len=256) outputname
-
-! Stacey conditions put back
- integer nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax,ispec2D
- integer, dimension(:,:), allocatable :: nimin,nimax,njmin,njmax,nkmin_xi,nkmin_eta
- real(kind=CUSTOM_REAL) vx,vy,vz,nx,ny,nz,tx,ty,tz,vn,weight
-
-!jpampuero Fault variables
- real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: FaultStrength,FaultTxInit,FaultTx & !jpampuero
- , FaultTzInit,FaultTz,FaultMuS,FaultMuD,FaultSWR ,FaultB !jpampuero
- real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: FaultDisplX,FaultDisplZ & !jpampuero
- ,FaultVelocX,FaultVelocZ & !jpampuero
- ,FaultStressX,FaultStressZ !jpampuero
- real(kind=CUSTOM_REAL) :: FaultDVxFree,FaultDVzFree,FaultZ,va,ta,txa,tza,dirx,dirz !jpampuero
- integer, dimension(FaultOutputNPTS) :: FaultIglobOut,FaultIgllOut,FaultKgllOut,FaultElmtOut !jpampuero
- character(len=70) FaultOutputName !jpampuero
- character(len=70), dimension(FaultOutputNPTS) :: FaultOutputPtsName !jpampuero
- logical :: FaultInThisProc,NotInRuptureRegion,InNucleationRegion !jpampuero
-
-! ************** PROGRAM STARTS HERE **************
-
-! initialize the MPI communicator and start the NPROC MPI processes.
-! sizeprocs returns number of processes started
-! (should be equal to NPROC)
-! myrank is the rank of each process, between 0 and sizeprocs-1.
-! as usual in MPI, process 0 is in charge of coordinating everything
-! and also takes care of the main output
- call MPI_INIT(ier)
- call MPI_COMM_SIZE(MPI_COMM_WORLD,sizeprocs,ier)
- call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ier)
-
-! read the parameter file
- call read_parameter_file(LAT_MIN,LAT_MAX,LONG_MIN,LONG_MAX, &
- UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK, &
- NER_TOP_LAYER,NER_MIDDLE_LAYER,NER_BOTTOM_HALFSPACE, &
- NEX_ETA,NEX_XI,NPROC_ETA,NPROC_XI,NSEIS,NSTEP,UTM_PROJECTION_ZONE,DT, &
- ATTENUATION,LOCAL_PATH,NSOURCES,STACEY_ABS_CONDITIONS)
-
-! compute other parameters based upon values read
- call compute_parameters(NER,NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA, &
- NPROC,NEX_PER_PROC_XI,NEX_PER_PROC_ETA, &
- NER_BOTTOM_HALFSPACE,NER_MIDDLE_LAYER,NER_TOP_LAYER, &
- NSPEC_AB,NSPEC2D_A_XI,NSPEC2D_B_XI, &
- NSPEC2D_A_ETA,NSPEC2D_B_ETA, &
- NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
- NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NGLOB_AB)
-
-!jpampuero moved this section (used to be right before time loop)
-!jpampuero need some of these during initialization
-! distinguish whether single or double precision for reals
- if(CUSTOM_REAL == SIZE_REAL) then
- deltat = sngl(DT)
- else
- deltat = DT
- endif
- deltatover2 = deltat/2.
- deltatsqover2 = deltat*deltat/2.
-
-! open main output file, only written to by process 0
- if(myrank == 0 .and. IMAIN /= ISTANDARD_OUTPUT) &
- open(unit=IMAIN,file='OUTPUT_FILES/output_solver.txt',status='unknown')
-
- if(myrank == 0) then
-
- write(IMAIN,*)
- write(IMAIN,*) '**********************************************'
- write(IMAIN,*) '**** Specfem 3-D Solver - MPI version f90 ****'
- write(IMAIN,*) '**********************************************'
- write(IMAIN,*)
- write(IMAIN,*)
-
- if(FIX_UNDERFLOW_PROBLEM) write(IMAIN,*) 'Fixing slow underflow trapping problem using small initial field'
-
- write(IMAIN,*)
- write(IMAIN,*) 'There are ',sizeprocs,' MPI processes'
- write(IMAIN,*) 'Processes are numbered from 0 to ',sizeprocs-1
- write(IMAIN,*)
-
- write(IMAIN,*) 'There are ',NEX_XI,' elements along xi'
- write(IMAIN,*) 'There are ',NEX_ETA,' elements along eta'
- write(IMAIN,*)
- write(IMAIN,*) 'There are ',NPROC_XI,' slices along xi'
- write(IMAIN,*) 'There are ',NPROC_ETA,' slices along eta'
- write(IMAIN,*) 'There is a total of ',NPROC,' slices'
-
- write(IMAIN,*)
- write(IMAIN,*) ' NDIM = ',NDIM
- write(IMAIN,*)
- write(IMAIN,*) ' NGLLX = ',NGLLX
- write(IMAIN,*) ' NGLLY = ',NGLLY
- write(IMAIN,*) ' NGLLZ = ',NGLLZ
- write(IMAIN,*)
-
-! write information about precision used for floating-point operations
- if(CUSTOM_REAL == SIZE_REAL) then
- write(IMAIN,*) 'using single precision for the calculations'
- else
- write(IMAIN,*) 'using double precision for the calculations'
- endif
- write(IMAIN,*)
- write(IMAIN,*) 'smallest and largest possible floating-point numbers are: ',tiny(1._CUSTOM_REAL),huge(1._CUSTOM_REAL)
- write(IMAIN,*)
-
- endif
-
-! check that the code is running with the requested nb of processes
- if(sizeprocs /= NPROC) call exit_MPI(myrank,'wrong number of MPI processes')
-
-! check that we have more than 0 and less than 1000 sources
- if(NSOURCES < 0 .or. NSOURCES > 999) call exit_MPI(myrank,'invalid number of sources')
-
-! dynamic allocation of arrays
-
-! 2-D addressing and buffers for summation between slices, and point codes
-! use number of elements found in the mantle since it is the largest region
-
-! crust and mantle
- allocate(iboolleft_xi(NPOIN2DMAX_XMIN_XMAX))
- allocate(iboolright_xi(NPOIN2DMAX_XMIN_XMAX))
- allocate(iboolleft_eta(NPOIN2DMAX_YMIN_YMAX))
- allocate(iboolright_eta(NPOIN2DMAX_YMIN_YMAX))
-
-! for addressing of the slices
- allocate(addressing(0:NPROC_XI-1,0:NPROC_ETA-1))
- allocate(iproc_xi_slice(0:NPROC-1))
- allocate(iproc_eta_slice(0:NPROC-1))
-
-! open file with global slice number addressing
- open(unit=IIN,file='OUTPUT_FILES/addressing.txt',status='old')
- do iproc = 0,NPROC-1
- read(IIN,*) iproc_read,iproc_xi,iproc_eta
- if(iproc_read /= iproc) call exit_MPI(myrank,'incorrect slice number read')
- addressing(iproc_xi,iproc_eta) = iproc
- iproc_xi_slice(iproc) = iproc_xi
- iproc_eta_slice(iproc) = iproc_eta
- enddo
- close(IIN)
-
-! determine local slice coordinates using addressing
- iproc_xi = iproc_xi_slice(myrank)
- iproc_eta = iproc_eta_slice(myrank)
-
- FaultInThisProc = iproc_eta==0 !jpampuero
-
-! define maximum size for message buffers
- NPOIN2DMAX_XY = max(NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX)
-
- allocate(buffer_send_faces(NDIM,NPOIN2DMAX_XY))
- allocate(buffer_received_faces(NDIM,NPOIN2DMAX_XY))
-
-! start reading the databases
-
-! read arrays created by the mesher
-
- call read_arrays_solver(myrank,xstore,ystore,zstore, &
- xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,jacobian, &
- rho_vp,rho_vs,kappastore,mustore,ibool,idoubling,rmass(1,:),LOCAL_PATH) !jpampuero ABC
- rmass(2,:) = rmass(1,:) !jpampuero ABC
- rmass(3,:) = rmass(1,:) !jpampuero ABC
-
-! check that the number of points in this slice is correct
-
- if(minval(ibool(:,:,:,:)) /= 1 .or. maxval(ibool(:,:,:,:)) /= NGLOB_AB) &
- call exit_MPI(myrank,'incorrect global numbering: iboolmax does not equal nglob in crust and mantle')
-
-! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
-
-if (NSOURCES==0) then !jpampuero
- allocate(hdur(1))
- allocate(utm_x_source(1))
- allocate(utm_y_source(1))
- hdur(1) = 0.d0
- utm_x_source(1) = 0.d0
- utm_y_source(1) = 0.d0
-else
-! allocate arrays for source
- allocate(sourcearray(NDIM,NGLLX,NGLLY,NGLLZ))
- allocate(sourcearrays(NSOURCES,NDIM,NGLLX,NGLLY,NGLLZ))
- allocate(islice_selected_source(NSOURCES))
- allocate(ispec_selected_source(NSOURCES))
- allocate(Mxx(NSOURCES))
- allocate(Myy(NSOURCES))
- allocate(Mzz(NSOURCES))
- allocate(Mxy(NSOURCES))
- allocate(Mxz(NSOURCES))
- allocate(Myz(NSOURCES))
- allocate(xi_source(NSOURCES))
- allocate(eta_source(NSOURCES))
- allocate(gamma_source(NSOURCES))
- allocate(t_cmt(NSOURCES))
- allocate(hdur(NSOURCES))
- allocate(utm_x_source(NSOURCES))
- allocate(utm_y_source(NSOURCES))
-
-! locate sources in the mesh
- do isource = 1,NSOURCES
- call locate_source(ibool,isource,NSOURCES,myrank,NSPEC_AB,NGLOB_AB,xstore,ystore,zstore, &
- xigll,yigll,zigll,NPROC, &
- sec,t_cmt(isource),yr,jda,ho,mi,utm_x_source(isource),utm_y_source(isource), &
- NSTEP,DT,hdur(isource),Mxx(isource),Myy(isource),Mzz(isource),Mxy(isource),Mxz(isource),Myz(isource), &
- islice_selected_source(isource),ispec_selected_source(isource), &
- xi_source(isource),eta_source(isource),gamma_source(isource), &
- LAT_MIN,LAT_MAX,LONG_MIN,LONG_MAX,Z_DEPTH_BLOCK,UTM_PROJECTION_ZONE)
- enddo
-
- if(t_cmt(1) /= 0.) call exit_MPI(myrank,'t_cmt for the first source should be zero')
- do isource = 2,NSOURCES
- if(t_cmt(isource) < 0.) call exit_MPI(myrank,'t_cmt should not be less than zero')
- enddo
-endif !jpampuero
-
- open(unit=IIN,file='DATA/STATIONS_FILTERED',status='old')
- read(IIN,*) nrec
- close(IIN)
-
- if(myrank == 0) then
- write(IMAIN,*)
- write(IMAIN,*) 'Total number of receivers = ',nrec
- write(IMAIN,*)
- endif
-
- if(nrec < 1) call exit_MPI(myrank,'need at least one receiver')
-
-! allocate memory for receiver arrays
- allocate(islice_selected_rec(nrec))
- allocate(ispec_selected_rec(nrec))
- allocate(xi_receiver(nrec))
- allocate(eta_receiver(nrec))
- allocate(station_name(nrec))
- allocate(network_name(nrec))
- allocate(nu(NDIM,NDIM,nrec))
-
-! locate receivers in the crust in the mesh
- call locate_receivers(ibool,myrank,NSPEC_AB,NGLOB_AB,idoubling, &
- xstore,ystore,zstore,xigll,yigll, &
- nrec,islice_selected_rec,ispec_selected_rec, &
- xi_receiver,eta_receiver,station_name,network_name,nu, &
- NPROC,utm_x_source(1),utm_y_source(1),UTM_PROJECTION_ZONE)
-
-! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
-
-! read 2-D addressing for summation between slices with MPI
-
- call read_arrays_buffers_solver(myrank,iboolleft_xi, &
- iboolright_xi,iboolleft_eta,iboolright_eta, &
- npoin2D_xi,npoin2D_eta, &
- NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,LOCAL_PATH)
-
-! $$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
-
- if(myrank == 0) then
- write(IMAIN,*) '******************************************'
- write(IMAIN,*) 'There are ',NEX_XI,' elements along xi'
- write(IMAIN,*) 'There are ',NEX_ETA,' elements along eta'
- write(IMAIN,*)
- write(IMAIN,*) 'There are ',NPROC_XI,' slices along xi'
- write(IMAIN,*) 'There are ',NPROC_ETA,' slices along eta'
- write(IMAIN,*) 'There is a total of ',NPROC,' slices'
- write(IMAIN,*) '******************************************'
- write(IMAIN,*)
- endif
-
-! set up GLL points, weights and derivation matrices
- call define_derivation_matrices(xigll,yigll,zigll,wxgll,wygll,wzgll, &
- hprime_xx,hprime_yy,hprime_zz, &
- hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
- wgllwgll_xy,wgllwgll_xz,wgllwgll_yz)
-
-! allocate 1-D Lagrange interpolators and derivatives
- allocate(hxir(NGLLX))
- allocate(hpxir(NGLLX))
- allocate(hetar(NGLLY))
- allocate(hpetar(NGLLY))
-
-! create name of database
- call create_name_database(prname,myrank,LOCAL_PATH)
-
-! dynamic allocation of arrays
-
-! boundary parameters locator
- allocate(ibelm_xmin(NSPEC2DMAX_XMIN_XMAX))
- allocate(ibelm_xmax(NSPEC2DMAX_XMIN_XMAX))
- allocate(ibelm_ymin(NSPEC2DMAX_YMIN_YMAX))
- allocate(ibelm_ymax(NSPEC2DMAX_YMIN_YMAX))
- allocate(ibelm_bottom(NSPEC2D_BOTTOM))
- allocate(ibelm_top(NSPEC2D_TOP))
-
- allocate(jacobian2D_xmin(NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX))
- allocate(jacobian2D_xmax(NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX))
- allocate(jacobian2D_ymin(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX))
- allocate(jacobian2D_ymax(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX))
- allocate(jacobian2D_bottom(NGLLX,NGLLY,NSPEC2D_BOTTOM))
- allocate(jacobian2D_top(NGLLX,NGLLY,NSPEC2D_TOP))
-
-! Stacey put back
- allocate(nimin(2,NSPEC2DMAX_YMIN_YMAX))
- allocate(nimax(2,NSPEC2DMAX_YMIN_YMAX))
- allocate(njmin(2,NSPEC2DMAX_XMIN_XMAX))
- allocate(njmax(2,NSPEC2DMAX_XMIN_XMAX))
- allocate(nkmin_xi(2,NSPEC2DMAX_XMIN_XMAX))
- allocate(nkmin_eta(2,NSPEC2DMAX_YMIN_YMAX))
-
-! normals
- allocate(normal_xmin(NDIM,NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX))
- allocate(normal_xmax(NDIM,NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX))
-!!jpampuero allocate(normal_ymin(NDIM,NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX))
- allocate(normal_ymax(NDIM,NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX))
- allocate(normal_bottom(NDIM,NGLLX,NGLLY,NSPEC2D_BOTTOM))
- allocate(normal_top(NDIM,NGLLX,NGLLY,NSPEC2D_TOP))
-
-!jpampuero Fault
- if (FaultInThisProc) then
- allocate(FaultB(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)) !jpampuero
- allocate(FaultTx(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)) !jpampuero
- allocate(FaultTxInit(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)) !jpampuero
- allocate(FaultTz(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)) !jpampuero
- allocate(FaultTzInit(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)) !jpampuero
- allocate(FaultStrength(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)) !jpampuero
- allocate(FaultMuS(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)) !jpampuero
- allocate(FaultMuD(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)) !jpampuero
- allocate(FaultSWR(NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)) !jpampuero
- allocate(FaultDisplX(FaultOutputNPTS,NSTEP)) !jpampuero
- allocate(FaultDisplZ(FaultOutputNPTS,NSTEP)) !jpampuero
- allocate(FaultVelocX(FaultOutputNPTS,NSTEP)) !jpampuero
- allocate(FaultVelocZ(FaultOutputNPTS,NSTEP)) !jpampuero
- allocate(FaultStressX(FaultOutputNPTS,NSTEP)) !jpampuero
- allocate(FaultStressZ(FaultOutputNPTS,NSTEP)) !jpampuero
- endif
-
-! boundary parameters
- open(unit=27,file=prname(1:len_trim(prname))//'ibelm.bin',status='old',form='unformatted')
- read(27) ibelm_xmin
- read(27) ibelm_xmax
- read(27) ibelm_ymin
- read(27) ibelm_ymax
- read(27) ibelm_bottom
- read(27) ibelm_top
- close(27)
-
- open(unit=27,file=prname(1:len_trim(prname))//'normal.bin',status='old',form='unformatted')
- read(27) normal_xmin
- read(27) normal_xmax
-!!jpampuero read(27) normal_ymin
- read(27) normal_ymax !jpampuero
- read(27) normal_ymax
- read(27) normal_bottom
- read(27) normal_top
- close(27)
-
-! Stacey put back
- open(unit=27,file=prname(1:len_trim(prname))//'nspec2D.bin',status='unknown',form='unformatted')
- read(27) nspec2D_xmin
- read(27) nspec2D_xmax
- read(27) nspec2D_ymin
- read(27) nspec2D_ymax
- close(27)
-
- open(unit=27,file=prname(1:len_trim(prname))//'jacobian2D.bin',status='old',form='unformatted')
- read(27) jacobian2D_xmin
- read(27) jacobian2D_xmax
- read(27) jacobian2D_ymin
- read(27) jacobian2D_ymax
- read(27) jacobian2D_bottom
- read(27) jacobian2D_top
- close(27)
-
-! Stacey put back
-! read arrays for Stacey conditions
-
- if(STACEY_ABS_CONDITIONS) then
- open(unit=27,file=prname(1:len_trim(prname))//'nimin.bin',status='unknown',form='unformatted')
- read(27) nimin
- close(27)
-
- open(unit=27,file=prname(1:len_trim(prname))//'nimax.bin',status='unknown',form='unformatted')
- read(27) nimax
- close(27)
-
- open(unit=27,file=prname(1:len_trim(prname))//'njmin.bin',status='unknown',form='unformatted')
- read(27) njmin
- close(27)
-
- open(unit=27,file=prname(1:len_trim(prname))//'njmax.bin',status='unknown',form='unformatted')
- read(27) njmax
- close(27)
-
- open(unit=27,file=prname(1:len_trim(prname))//'nkmin_xi.bin',status='unknown',form='unformatted')
- read(27) nkmin_xi
- close(27)
-
- open(unit=27,file=prname(1:len_trim(prname))//'nkmin_eta.bin',status='unknown',form='unformatted')
- read(27) nkmin_eta
- close(27)
- endif
-
-!jpampuero ======= BEGIN INITIALIZE FAULT ========================
- if (FaultInThisProc) then
-
- accel(:,:) = 0._CUSTOM_REAL
-
- !OUTPUTS: locate selected fault points
- call fault_locate_rec(FaultIglobOut,FaultIgllOut,FaultKgllOut,FaultElmtOut,FaultOutputPtsName &
- ,nspec2D_ymin,ibelm_ymin,NSPEC2DMAX_YMIN_YMAX,ibool,xstore,zstore) !jpampuero
-
- !Set inital stress and strength for fault
- FaultTzInit = 0._CUSTOM_REAL
- FaultTxInit = BackgroundShearStress
- do ispec2D=1,nspec2D_ymin
- ispec=ibelm_ymin(ispec2D)
- do k=1,NGLLZ
- do i=1,NGLLX
-
- iglob=ibool(i,1,k,ispec)
- NotInRuptureRegion = abs(zstore(iglob))>RuptureDepth &
- .OR. abs(xstore(iglob))>RuptureHalfLength
- InNucleationRegion = abs(xstore(iglob)-HypoX)<=NucleationHalfSizeX &
- .AND. abs(zstore(iglob)-HypoZ)<=NucleationHalfSizeZ
-
- if ( NotInRuptureRegion ) then
- FaultMuS(i,k,ispec2D) = BarrierStrength
- FaultMuD(i,k,ispec2D) = BarrierStrength
- else
- FaultMuS(i,k,ispec2D) = StaticStrength
- FaultMuD(i,k,ispec2D) = DynamicStrength
- if ( InNucleationRegion ) then
- FaultTxInit(i,k,ispec2D) = NucleationShearStress
- endif
- endif
-
- !WARNING: factor 2 comes from using displacement instead of slip
- FaultSWR(i,k,ispec2D) = 2._CUSTOM_REAL*(FaultMuS(i,k,ispec2D)-FaultMuD(i,k,ispec2D))/Dc
-
- !assemble boundary weights, temporarily store in "accel" global array
- accel(1,iglob) = accel(1,iglob) + jacobian2D_ymin(i,k,ispec2D)*wgllwgll_xz(i,k)
-
- enddo
- enddo
- enddo
-
- do ispec2D=1,nspec2D_ymin
- ispec=ibelm_ymin(ispec2D)
- do k=1,NGLLZ
- do i=1,NGLLX
- iglob=ibool(i,1,k,ispec)
- FaultB(i,k,ispec2D) = accel(1,iglob)
- enddo
- enddo
- enddo
-
- endif
-!jpampuero ======= END INITIALIZE FAULT ========================
-
- do isource = 1,NSOURCES
-
-! check that the source slice number is okay
- if(islice_selected_source(isource) < 0 .or. islice_selected_source(isource) > NPROC-1) &
- call exit_MPI(myrank,'something is wrong with the source slice number')
-
-! compute source arrays in source slice
- if(myrank == islice_selected_source(isource)) then
- call compute_arrays_source(ispec_selected_source(isource), &
- xi_source(isource),eta_source(isource),gamma_source(isource),sourcearray, &
- Mxx(isource),Myy(isource),Mzz(isource),Mxy(isource),Mxz(isource),Myz(isource), &
- xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
- xigll,yigll,zigll,NSPEC_AB)
- sourcearrays(isource,:,:,:,:) = sourcearray(:,:,:,:)
- endif
-
- enddo
-
- if(myrank == 0) then
- write(IMAIN,*)
- write(IMAIN,*) 'Total number of samples for seismograms = ',NSTEP
- write(IMAIN,*)
- endif
-
-!--- select local receivers
-
-! count number of receivers located in this slice
- nrec_local = 0
-
- do irec = 1,nrec
-
-! check that the receiver slice number is okay
- if(islice_selected_rec(irec) < 0 .or. islice_selected_rec(irec) > NPROC-1) &
- call exit_MPI(myrank,'something is wrong with the receiver slice number')
-
-! write info about that receiver
- if(myrank == islice_selected_rec(irec)) then
-
- nrec_local = nrec_local + 1
-
- endif
-
- enddo
-
-! allocate seismogram array
- allocate(seismograms_d(NDIM,nrec_local,NSTEP))
- allocate(seismograms_v(NDIM,nrec_local,NSTEP))
-
-! allocate Lagrange interpolators for receivers
- allocate(hxir_store(nrec_local,NGLLX))
- allocate(hetar_store(nrec_local,NGLLY))
-
-! define local to global receiver numbering mapping
- allocate(number_receiver_global(nrec_local))
- irec_local = 0
- do irec = 1,nrec
- if(myrank == islice_selected_rec(irec)) then
- irec_local = irec_local + 1
- number_receiver_global(irec_local) = irec
- endif
- enddo
-
-! define and store Lagrange interpolators at all the receivers
- do irec_local = 1,nrec_local
- irec = number_receiver_global(irec_local)
- call lagrange_any(xi_receiver(irec),NGLLX,xigll,hxir,hpxir)
- call lagrange_any(eta_receiver(irec),NGLLY,yigll,hetar,hpetar)
- hxir_store(irec_local,:) = hxir(:)
- hetar_store(irec_local,:) = hetar(:)
- enddo
-
-! check that the sum of the number of receivers in each slice is nrec
- call MPI_REDUCE(nrec_local,nrec_tot_found,1,MPI_INTEGER,MPI_SUM,0, &
- MPI_COMM_WORLD,ier)
- if(myrank == 0) then
- write(IMAIN,*)
- write(IMAIN,*) 'found a total of ',nrec_tot_found,' receivers in all the slices'
- if(nrec_tot_found /= nrec) then
- call exit_MPI(myrank,'problem when dispatching the receivers')
- else
- write(IMAIN,*) 'this total is okay'
- endif
- endif
-
-! initialize seismograms
- seismograms_d(:,:,:) = 0._CUSTOM_REAL
- seismograms_v(:,:,:) = 0._CUSTOM_REAL
-
- if(myrank == 0) then
-
- if(NSOURCES > 1) write(IMAIN,*) 'Using ',NSOURCES,' point sources'
-
- write(IMAIN,*)
- if(ATTENUATION) then
- write(IMAIN,*) 'incorporating attenuation using ',N_SLS,' standard linear solids'
- else
- write(IMAIN,*) 'no attenuation'
- endif
- write(IMAIN,*)
-
- endif
-
-! synchronize all the processes before assembling the mass matrix
-! to make sure all the nodes have finished to read their databases
- call MPI_BARRIER(MPI_COMM_WORLD,ier)
-
-!jpampuero ABC ==== BEGIN modification of mass matrix for ============================
-!jpampuero ABC IMPLICIT implementation of absorbing boundaries
-!jpampuero ABC M --> M+dt/2*C
-!jpampuero ABC where C is related to the boundary term -C*v
-!jpampuero ABC NOTE: assumes straight edges
- if(STACEY_ABS_CONDITIONS) then
-
- ! vn=vx*nx+vy*ny+vz*nz
- ! tx=rho_vp(i,j,k,ispec)*vn*nx+rho_vs(i,j,k,ispec)*(vx-vn*nx)
- ! ty=rho_vp(i,j,k,ispec)*vn*ny+rho_vs(i,j,k,ispec)*(vy-vn*ny)
- ! tz=rho_vp(i,j,k,ispec)*vn*nz+rho_vs(i,j,k,ispec)*(vz-vn*nz)
-
- ! xmin: (nx,ny,nz) = (-1,0,0)
- ! vn=-vx
- ! tx=rho_vp(i,j,k,ispec)*vx
- ! ty=rho_vs(i,j,k,ispec)*vy
- ! tz=rho_vs(i,j,k,ispec)*vz
- do ispec2D=1,nspec2D_xmin
- ispec=ibelm_xmin(ispec2D)
- ! exclude elements that are not on absorbing edges
- if(nkmin_xi(1,ispec2D) == 0 .or. njmin(1,ispec2D) == 0) cycle
- i=1
- do k=nkmin_xi(1,ispec2D),NGLLZ
- do j=njmin(1,ispec2D),njmax(1,ispec2D)
- iglob=ibool(i,j,k,ispec)
- weight=jacobian2D_xmin(j,k,ispec2D)*wgllwgll_yz(j,k)
- rmass(1,iglob)=rmass(1,iglob) + deltatover2*weight*rho_vp(i,j,k,ispec)
- rmass(2,iglob)=rmass(2,iglob) + deltatover2*weight*rho_vs(i,j,k,ispec)
- rmass(3,iglob)=rmass(3,iglob) + deltatover2*weight*rho_vs(i,j,k,ispec)
- enddo
- enddo
- enddo
-
-
- ! xmax: (nx,ny,nz) = (1,0,0)
- ! vn=vx
- ! tx=rho_vp(i,j,k,ispec)*vx
- ! ty=rho_vs(i,j,k,ispec)*vy
- ! tz=rho_vs(i,j,k,ispec)*vz
- do ispec2D=1,nspec2D_xmax
- ispec=ibelm_xmax(ispec2D)
- ! exclude elements that are not on absorbing edges
- if(nkmin_xi(2,ispec2D) == 0 .or. njmin(2,ispec2D) == 0) cycle
- i=NGLLX
- do k=nkmin_xi(2,ispec2D),NGLLZ
- do j=njmin(2,ispec2D),njmax(2,ispec2D)
- iglob=ibool(i,j,k,ispec)
- weight=jacobian2D_xmax(j,k,ispec2D)*wgllwgll_yz(j,k)
- rmass(1,iglob)=rmass(1,iglob) + deltatover2*weight*rho_vp(i,j,k,ispec)
- rmass(2,iglob)=rmass(2,iglob) + deltatover2*weight*rho_vs(i,j,k,ispec)
- rmass(3,iglob)=rmass(3,iglob) + deltatover2*weight*rho_vs(i,j,k,ispec)
- enddo
- enddo
- enddo
-
- ! ymax: (nx,ny,nz) = (0,1,0)
- ! vn=vy
- ! tx=rho_vs(i,j,k,ispec)*vx
- ! ty=rho_vp(i,j,k,ispec)*vy
- ! tz=rho_vs(i,j,k,ispec)*vz
- do ispec2D=1,nspec2D_ymax
- ispec=ibelm_ymax(ispec2D)
- ! exclude elements that are not on absorbing edges
- if(nkmin_eta(2,ispec2D) == 0 .or. nimin(2,ispec2D) == 0) cycle
- j=NGLLY
- do k=nkmin_eta(2,ispec2D),NGLLZ
- do i=nimin(2,ispec2D),nimax(2,ispec2D)
- iglob=ibool(i,j,k,ispec)
- weight=jacobian2D_ymax(i,k,ispec2D)*wgllwgll_xz(i,k)
- rmass(1,iglob)=rmass(1,iglob) + deltatover2*weight*rho_vs(i,j,k,ispec)
- rmass(2,iglob)=rmass(2,iglob) + deltatover2*weight*rho_vp(i,j,k,ispec)
- rmass(3,iglob)=rmass(3,iglob) + deltatover2*weight*rho_vs(i,j,k,ispec)
- enddo
- enddo
- enddo
-
- ! zmin: (nx,ny,nz) = (0,0,-1)
- ! vn=-vz
- ! tx=rho_vs(i,j,k,ispec)*vx
- ! ty=rho_vs(i,j,k,ispec)*vy
- ! tz=rho_vp(i,j,k,ispec)*vz
- do ispec2D=1,NSPEC2D_BOTTOM
- ispec=ibelm_bottom(ispec2D)
- k=1
- do j=1,NGLLY
- do i=1,NGLLX
- iglob=ibool(i,j,k,ispec)
- weight=jacobian2D_bottom(i,j,ispec2D)*wgllwgll_xy(i,j)
- rmass(1,iglob)=rmass(1,iglob) + deltatover2*weight*rho_vs(i,j,k,ispec)
- rmass(2,iglob)=rmass(2,iglob) + deltatover2*weight*rho_vs(i,j,k,ispec)
- rmass(3,iglob)=rmass(3,iglob) + deltatover2*weight*rho_vp(i,j,k,ispec)
- enddo
- enddo
- enddo
-
- endif
-!jpampuero ABC ======= END modification of mass matrix ====================================
-
-! the mass matrix needs to be assembled with MPI here once and for all
- call assemble_MPI_vector(rmass,iproc_xi,iproc_eta,addressing, & !jpampuero ABC
- iboolleft_xi,iboolright_xi,iboolleft_eta,iboolright_eta, &
- buffer_send_faces,buffer_received_faces,npoin2D_xi,npoin2D_eta, &
- NPROC_XI,NPROC_ETA,NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NPOIN2DMAX_XY)
-
- if(myrank == 0) write(IMAIN,*) 'end assembling MPI mass matrix'
-
-! check that mass matrix is positive
- if(minval(rmass) <= 0.) call exit_MPI(myrank,'negative mass matrix term')
-
-! for efficiency, invert final mass matrix once and for all in each slice
- rmass = 1.0 / rmass
-
-! if attenuation is on, shift PREM to right frequency
-! rescale mu in PREM to average frequency for attenuation
-
- if(ATTENUATION) then
-
-! get and store PREM attenuation model
- do iattenuation = 1,NUM_REGIONS_ATTENUATION
-
- call get_attenuation_model(myrank,iattenuation,tau_mu_dble, &
- tau_sigma_dble,beta_dble,one_minus_sum_beta_dble,factor_scale_dble)
-
-! distinguish whether single or double precision for reals
- if(CUSTOM_REAL == SIZE_REAL) then
- tau_mu(iattenuation,:) = sngl(tau_mu_dble(:))
- tau_sigma(iattenuation,:) = sngl(tau_sigma_dble(:))
- beta(iattenuation,:) = sngl(beta_dble(:))
- factor_scale(iattenuation) = sngl(factor_scale_dble)
- one_minus_sum_beta(iattenuation) = sngl(one_minus_sum_beta_dble)
- else
- tau_mu(iattenuation,:) = tau_mu_dble(:)
- tau_sigma(iattenuation,:) = tau_sigma_dble(:)
- beta(iattenuation,:) = beta_dble(:)
- factor_scale(iattenuation) = factor_scale_dble
- one_minus_sum_beta(iattenuation) = one_minus_sum_beta_dble
- endif
- enddo
-
-! rescale shear modulus according to attenuation model
-
- do ispec = 1,NSPEC_AB
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- scale_factor = factor_scale(IATTENUATION_UNIFORM)
- mustore(i,j,k,ispec) = mustore(i,j,k,ispec) * scale_factor
- enddo
- enddo
- enddo
- enddo
-
- endif
-
-! initialize arrays to zero
- displ(:,:) = 0._CUSTOM_REAL
- veloc(:,:) = 0._CUSTOM_REAL
- accel(:,:) = 0._CUSTOM_REAL
-
-! put negligible initial value to avoid very slow underflow trapping
- if(FIX_UNDERFLOW_PROBLEM) displ(:,:) = VERYSMALLVAL
-
-! synchronize all processes to make sure everybody is ready to start time loop
- call MPI_BARRIER(MPI_COMM_WORLD,ier)
- if(myrank == 0) write(IMAIN,*) 'All processes are synchronized before time loop'
-
-!
-! s t a r t t i m e i t e r a t i o n s
-!
-
- if(myrank == 0) then
- write(IMAIN,*)
- write(IMAIN,*) ' time step: ',sngl(DT),' s'
- write(IMAIN,*) 'number of time steps: ',NSTEP
- write(IMAIN,*) 'total simulated time: ',sngl(NSTEP*DT),' seconds'
- write(IMAIN,*)
- endif
-
-! precompute Runge-Kutta coefficients if attenuation
- if(ATTENUATION) then
- tauinv(:,:) = - 1. / tau_sigma(:,:)
- factor_common(:,:) = 2. * beta(:,:) * tauinv(:,:)
- alphaval(:,:) = 1 + deltat*tauinv(:,:) + deltat**2*tauinv(:,:)**2 / 2. + &
- deltat**3*tauinv(:,:)**3 / 6. + deltat**4*tauinv(:,:)**4 / 24.
- betaval(:,:) = deltat / 2. + deltat**2*tauinv(:,:) / 3. + deltat**3*tauinv(:,:)**2 / 8. + deltat**4*tauinv(:,:)**3 / 24.
- gammaval(:,:) = deltat / 2. + deltat**2*tauinv(:,:) / 6. + deltat**3*tauinv(:,:)**2 / 24.
- endif
-
- if(myrank == 0) then
- write(IMAIN,*)
- write(IMAIN,*) 'Starting time iteration loop...'
- write(IMAIN,*)
- endif
-
-! create an empty file to monitor the start of the simulation
- if(myrank == 0) then
- open(unit=IOUT,file='OUTPUT_FILES/starttimeloop.txt',status='unknown')
- write(IOUT,*) 'starting time loop'
- close(IOUT)
- endif
-
-! get MPI starting time
- time_start = MPI_WTIME()
-
-! clear memory variables if attenuation
- if(ATTENUATION) then
- R_xx(:,:,:,:,:) = 0._CUSTOM_REAL
- R_yy(:,:,:,:,:) = 0._CUSTOM_REAL
- R_xy(:,:,:,:,:) = 0._CUSTOM_REAL
- R_xz(:,:,:,:,:) = 0._CUSTOM_REAL
- R_yz(:,:,:,:,:) = 0._CUSTOM_REAL
-
- if(FIX_UNDERFLOW_PROBLEM) then
- R_xx(:,:,:,:,:) = VERYSMALLVAL
- R_yy(:,:,:,:,:) = VERYSMALLVAL
- R_xy(:,:,:,:,:) = VERYSMALLVAL
- R_xz(:,:,:,:,:) = VERYSMALLVAL
- R_yz(:,:,:,:,:) = VERYSMALLVAL
- endif
-
- endif
-
-! *********************************************************
-! ************* MAIN LOOP OVER THE TIME STEPS *************
-! *********************************************************
-
- do it=1,NSTEP
-
-! compute the maximum of the norm of the displacement
-! in all the slices using an MPI reduction
-! and output timestamp file to check that simulation is running fine
- if(mod(it,ITAFF_TIME_STEPS) == 0 .or. it == 5) then
-
-! compute maximum of norm of displacement in each slice
- Usolidnorm = sqrt(maxval(displ(1,:)**2 + displ(2,:)**2 + displ(3,:)**2))
-
-! compute the maximum of the maxima for all the slices using an MPI reduction
- call MPI_REDUCE(Usolidnorm,Usolidnorm_all,1,CUSTOM_MPI_TYPE,MPI_MAX,0, &
- MPI_COMM_WORLD,ier)
-
- if(myrank == 0) then
-
- write(IMAIN,*) 'Time step # ',it
- write(IMAIN,*) 'Time: ',sngl((it-1)*DT-hdur(1)),' seconds'
-
-! elapsed time since beginning of the simulation
- tCPU = MPI_WTIME() - time_start
- int_tCPU = int(tCPU)
- ihours = int_tCPU / 3600
- iminutes = (int_tCPU - 3600*ihours) / 60
- iseconds = int_tCPU - 3600*ihours - 60*iminutes
- write(IMAIN,*) 'Elapsed time in seconds = ',tCPU
- write(IMAIN,"(' Elapsed time in hh:mm:ss = ',i4,' h ',i2.2,' m ',i2.2,' s')") ihours,iminutes,iseconds
- write(IMAIN,*) 'Mean elapsed time per time step in seconds = ',tCPU/dble(it)
- write(IMAIN,*) 'Max norm displacement vector U in solid in all slices (m) = ',Usolidnorm_all
- write(IMAIN,*)
-
-! write time stamp file to give information about progression of simulation
- write(outputname,"('OUTPUT_FILES/timestamp',i6.6)") it
- open(unit=IOUT,file=outputname,status='unknown')
- write(IOUT,*) 'Time step # ',it
- write(IOUT,*) 'Time: ',sngl((it-1)*DT-hdur(1)),' seconds'
- write(IOUT,*) 'Elapsed time in seconds = ',tCPU
- write(IOUT,"(' Elapsed time in hh:mm:ss = ',i4,' h ',i2.2,' m ',i2.2,' s')") ihours,iminutes,iseconds
- write(IOUT,*) 'Mean elapsed time per time step in seconds = ',tCPU/dble(it)
- write(IOUT,*) 'Max norm displacement vector U in solid in all slices (m) = ',Usolidnorm_all
- close(IOUT)
-
-! check stability of the code, exit if unstable
- if(Usolidnorm_all > STABILITY_THRESHOLD) call exit_MPI(myrank,'code became unstable and blew up')
-
- endif
- endif
-
-! update displacement using finite difference time scheme
- do i=1,NGLOB_AB
- displ(:,i) = displ(:,i) + deltat*veloc(:,i) + deltatsqover2*accel(:,i)
- veloc(:,i) = veloc(:,i) + deltatover2*accel(:,i)
- accel(:,i) = 0._CUSTOM_REAL
- enddo
-
- do ispec = 1,NSPEC_AB
-
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
-
- tempx1l = 0.
- tempx2l = 0.
- tempx3l = 0.
-
- tempy1l = 0.
- tempy2l = 0.
- tempy3l = 0.
-
- tempz1l = 0.
- tempz2l = 0.
- tempz3l = 0.
-
- do l=1,NGLLX
- hp1 = hprime_xx(l,i)
- iglob = ibool(l,j,k,ispec)
- tempx1l = tempx1l + displ(1,iglob)*hp1
- tempy1l = tempy1l + displ(2,iglob)*hp1
- tempz1l = tempz1l + displ(3,iglob)*hp1
- enddo
-
- do l=1,NGLLY
- hp2 = hprime_yy(l,j)
- iglob = ibool(i,l,k,ispec)
- tempx2l = tempx2l + displ(1,iglob)*hp2
- tempy2l = tempy2l + displ(2,iglob)*hp2
- tempz2l = tempz2l + displ(3,iglob)*hp2
- enddo
-
- do l=1,NGLLZ
- hp3 = hprime_zz(l,k)
- iglob = ibool(i,j,l,ispec)
- tempx3l = tempx3l + displ(1,iglob)*hp3
- tempy3l = tempy3l + displ(2,iglob)*hp3
- tempz3l = tempz3l + displ(3,iglob)*hp3
- enddo
-
-! get derivatives of ux, uy and uz with respect to x, y and z
-!! DK DK NOJACOBIAN xixl = xix(i,j,k,ispec)
-!! DK DK NOJACOBIAN xiyl = xiy(i,j,k,ispec)
-!! DK DK NOJACOBIAN xizl = xiz(i,j,k,ispec)
-!! DK DK NOJACOBIAN etaxl = etax(i,j,k,ispec)
-!! DK DK NOJACOBIAN etayl = etay(i,j,k,ispec)
-!! DK DK NOJACOBIAN etazl = etaz(i,j,k,ispec)
-!! DK DK NOJACOBIAN gammaxl = gammax(i,j,k,ispec)
-!! DK DK NOJACOBIAN gammayl = gammay(i,j,k,ispec)
-!! DK DK NOJACOBIAN gammazl = gammaz(i,j,k,ispec)
-!! DK DK NOJACOBIAN jacobianl = jacobian(i,j,k,ispec)
- xixl = xix(i,j,k,1)
- xiyl = xiy(i,j,k,1)
- xizl = xiz(i,j,k,1)
- etaxl = etax(i,j,k,1)
- etayl = etay(i,j,k,1)
- etazl = etaz(i,j,k,1)
- gammaxl = gammax(i,j,k,1)
- gammayl = gammay(i,j,k,1)
- gammazl = gammaz(i,j,k,1)
- jacobianl = jacobian(i,j,k,1)
-
- duxdxl = xixl*tempx1l + etaxl*tempx2l + gammaxl*tempx3l
- duxdyl = xiyl*tempx1l + etayl*tempx2l + gammayl*tempx3l
- duxdzl = xizl*tempx1l + etazl*tempx2l + gammazl*tempx3l
-
- duydxl = xixl*tempy1l + etaxl*tempy2l + gammaxl*tempy3l
- duydyl = xiyl*tempy1l + etayl*tempy2l + gammayl*tempy3l
- duydzl = xizl*tempy1l + etazl*tempy2l + gammazl*tempy3l
-
- duzdxl = xixl*tempz1l + etaxl*tempz2l + gammaxl*tempz3l
- duzdyl = xiyl*tempz1l + etayl*tempz2l + gammayl*tempz3l
- duzdzl = xizl*tempz1l + etazl*tempz2l + gammazl*tempz3l
-
-! precompute some sums to save CPU time
- duxdxl_plus_duydyl = duxdxl + duydyl
- duxdxl_plus_duzdzl = duxdxl + duzdzl
- duydyl_plus_duzdzl = duydyl + duzdzl
- duxdyl_plus_duydxl = duxdyl + duydxl
- duzdxl_plus_duxdzl = duzdxl + duxdzl
- duzdyl_plus_duydzl = duzdyl + duydzl
-
-! precompute terms for attenuation if needed
- if(ATTENUATION) then
-
-! compute deviatoric strain
- epsilon_trace_over_3 = ONE_THIRD * (duxdxl + duydyl + duzdzl)
- epsilondev_xx_loc(i,j,k) = duxdxl - epsilon_trace_over_3
- epsilondev_yy_loc(i,j,k) = duydyl - epsilon_trace_over_3
- epsilondev_xy_loc(i,j,k) = 0.5 * duxdyl_plus_duydxl
- epsilondev_xz_loc(i,j,k) = 0.5 * duzdxl_plus_duxdzl
- epsilondev_yz_loc(i,j,k) = 0.5 * duzdyl_plus_duydzl
-
-! distinguish attenuation factors
- iselected = IATTENUATION_UNIFORM
- one_minus_sum_beta_use = one_minus_sum_beta(iselected)
- minus_sum_beta = one_minus_sum_beta_use - 1.
-
- endif
-
- kappal = kappastore(i,j,k,ispec)
- mul = mustore(i,j,k,ispec)
-
-! use unrelaxed parameters if attenuation
- if(ATTENUATION) mul = mul * one_minus_sum_beta_use
-
- lambdalplus2mul = kappal + FOUR_THIRDS * mul
- lambdal = lambdalplus2mul - 2.*mul
-
-! compute stress sigma
- sigma_xx = lambdalplus2mul*duxdxl + lambdal*duydyl_plus_duzdzl
- sigma_yy = lambdalplus2mul*duydyl + lambdal*duxdxl_plus_duzdzl
- sigma_zz = lambdalplus2mul*duzdzl + lambdal*duxdxl_plus_duydyl
-
- sigma_xy = mul*duxdyl_plus_duydxl
- sigma_xz = mul*duzdxl_plus_duxdzl
- sigma_yz = mul*duzdyl_plus_duydzl
-
-! subtract memory variables if attenuation
- if(ATTENUATION) then
- do i_sls = 1,N_SLS
- R_xx_val = R_xx(i,j,k,ispec,i_sls)
- R_yy_val = R_yy(i,j,k,ispec,i_sls)
- sigma_xx = sigma_xx - R_xx_val
- sigma_yy = sigma_yy - R_yy_val
- sigma_zz = sigma_zz + R_xx_val + R_yy_val
- sigma_xy = sigma_xy - R_xy(i,j,k,ispec,i_sls)
- sigma_xz = sigma_xz - R_xz(i,j,k,ispec,i_sls)
- sigma_yz = sigma_yz - R_yz(i,j,k,ispec,i_sls)
- enddo
- endif
-
-! form dot product with test vector, symmetric form
- tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_xy*xiyl + sigma_xz*xizl)
- tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_yz*xizl)
- tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl)
-
- tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_xy*etayl + sigma_xz*etazl)
- tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_yz*etazl)
- tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl)
-
- tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_xy*gammayl + sigma_xz*gammazl)
- tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_yz*gammazl)
- tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl)
-
- enddo
- enddo
- enddo
-
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
-
- tempx1l = 0.
- tempy1l = 0.
- tempz1l = 0.
-
- tempx2l = 0.
- tempy2l = 0.
- tempz2l = 0.
-
- tempx3l = 0.
- tempy3l = 0.
- tempz3l = 0.
-
- do l=1,NGLLX
- fac1 = hprimewgll_xx(i,l)
- tempx1l = tempx1l + tempx1(l,j,k)*fac1
- tempy1l = tempy1l + tempy1(l,j,k)*fac1
- tempz1l = tempz1l + tempz1(l,j,k)*fac1
- enddo
-
- do l=1,NGLLY
- fac2 = hprimewgll_yy(j,l)
- tempx2l = tempx2l + tempx2(i,l,k)*fac2
- tempy2l = tempy2l + tempy2(i,l,k)*fac2
- tempz2l = tempz2l + tempz2(i,l,k)*fac2
- enddo
-
- do l=1,NGLLZ
- fac3 = hprimewgll_zz(k,l)
- tempx3l = tempx3l + tempx3(i,j,l)*fac3
- tempy3l = tempy3l + tempy3(i,j,l)*fac3
- tempz3l = tempz3l + tempz3(i,j,l)*fac3
- enddo
-
- fac1 = wgllwgll_yz(j,k)
- fac2 = wgllwgll_xz(i,k)
- fac3 = wgllwgll_xy(i,j)
-
-! sum contributions from each element to the global mesh and add gravity terms
-
- iglob = ibool(i,j,k,ispec)
-
- accel(1,iglob) = accel(1,iglob) - (fac1*tempx1l + fac2*tempx2l + fac3*tempx3l)
- accel(2,iglob) = accel(2,iglob) - (fac1*tempy1l + fac2*tempy2l + fac3*tempy3l)
- accel(3,iglob) = accel(3,iglob) - (fac1*tempz1l + fac2*tempz2l + fac3*tempz3l)
-
-! update memory variables based upon the Runge-Kutta scheme
-
- if(ATTENUATION) then
-
-! use Runge-Kutta scheme to march in time
- do i_sls = 1,N_SLS
-
-! get coefficients for that standard linear solid
-
- factor_loc = mustore(i,j,k,ispec) * factor_common(iselected,i_sls)
- alphaval_loc = alphaval(iselected,i_sls)
- betaval_loc = betaval(iselected,i_sls)
- gammaval_loc = gammaval(iselected,i_sls)
-
-! term in xx
- Sn = factor_loc * epsilondev_xx(i,j,k,ispec)
- Snp1 = factor_loc * epsilondev_xx_loc(i,j,k)
- R_xx(i,j,k,ispec,i_sls) = alphaval_loc * R_xx(i,j,k,ispec,i_sls) + betaval_loc * Sn + gammaval_loc * Snp1
-
-! term in yy
- Sn = factor_loc * epsilondev_yy(i,j,k,ispec)
- Snp1 = factor_loc * epsilondev_yy_loc(i,j,k)
- R_yy(i,j,k,ispec,i_sls) = alphaval_loc * R_yy(i,j,k,ispec,i_sls) + betaval_loc * Sn + gammaval_loc * Snp1
-
-! term in zz not computed since zero trace
-
-! term in xy
- Sn = factor_loc * epsilondev_xy(i,j,k,ispec)
- Snp1 = factor_loc * epsilondev_xy_loc(i,j,k)
- R_xy(i,j,k,ispec,i_sls) = alphaval_loc * R_xy(i,j,k,ispec,i_sls) + betaval_loc * Sn + gammaval_loc * Snp1
-
-! term in xz
- Sn = factor_loc * epsilondev_xz(i,j,k,ispec)
- Snp1 = factor_loc * epsilondev_xz_loc(i,j,k)
- R_xz(i,j,k,ispec,i_sls) = alphaval_loc * R_xz(i,j,k,ispec,i_sls) + betaval_loc * Sn + gammaval_loc * Snp1
-
-! term in yz
- Sn = factor_loc * epsilondev_yz(i,j,k,ispec)
- Snp1 = factor_loc * epsilondev_yz_loc(i,j,k)
- R_yz(i,j,k,ispec,i_sls) = alphaval_loc * R_yz(i,j,k,ispec,i_sls) + betaval_loc * Sn + gammaval_loc * Snp1
-
- enddo ! end of loop on memory variables
-
- endif ! end attenuation
-
- enddo
- enddo
- enddo
-
-! save deviatoric strain for Runge-Kutta scheme
- if(ATTENUATION) then
- epsilondev_xx(:,:,:,ispec) = epsilondev_xx_loc(:,:,:)
- epsilondev_yy(:,:,:,ispec) = epsilondev_yy_loc(:,:,:)
- epsilondev_xy(:,:,:,ispec) = epsilondev_xy_loc(:,:,:)
- epsilondev_xz(:,:,:,ispec) = epsilondev_xz_loc(:,:,:)
- epsilondev_yz(:,:,:,ispec) = epsilondev_yz_loc(:,:,:)
- endif
-
- enddo ! spectral element loop
-
-!jpampuero *********** BEGIN ABSORBING BOUNDARY CONDITIONS ******************
-! add Stacey conditions
-
- if(STACEY_ABS_CONDITIONS) then
-
-! xmin
- do ispec2D=1,nspec2D_xmin
-
- ispec=ibelm_xmin(ispec2D)
-
-! exclude elements that are not on absorbing edges
- if(nkmin_xi(1,ispec2D) == 0 .or. njmin(1,ispec2D) == 0) cycle
-
- i=1
- do k=nkmin_xi(1,ispec2D),NGLLZ
- do j=njmin(1,ispec2D),njmax(1,ispec2D)
- iglob=ibool(i,j,k,ispec)
-
- vx=veloc(1,iglob)
- vy=veloc(2,iglob)
- vz=veloc(3,iglob)
-
- nx=normal_xmin(1,j,k,ispec2D)
- ny=normal_xmin(2,j,k,ispec2D)
- nz=normal_xmin(3,j,k,ispec2D)
-
- vn=vx*nx+vy*ny+vz*nz
-
- tx=rho_vp(i,j,k,ispec)*vn*nx+rho_vs(i,j,k,ispec)*(vx-vn*nx)
- ty=rho_vp(i,j,k,ispec)*vn*ny+rho_vs(i,j,k,ispec)*(vy-vn*ny)
- tz=rho_vp(i,j,k,ispec)*vn*nz+rho_vs(i,j,k,ispec)*(vz-vn*nz)
-
- weight=jacobian2D_xmin(j,k,ispec2D)*wgllwgll_yz(j,k)
-
- accel(1,iglob)=accel(1,iglob) - tx*weight
- accel(2,iglob)=accel(2,iglob) - ty*weight
- accel(3,iglob)=accel(3,iglob) - tz*weight
-
- enddo
- enddo
- enddo
-
-! xmax
- do ispec2D=1,nspec2D_xmax
- ispec=ibelm_xmax(ispec2D)
-
-! exclude elements that are not on absorbing edges
- if(nkmin_xi(2,ispec2D) == 0 .or. njmin(2,ispec2D) == 0) cycle
-
- i=NGLLX
- do k=nkmin_xi(2,ispec2D),NGLLZ
- do j=njmin(2,ispec2D),njmax(2,ispec2D)
- iglob=ibool(i,j,k,ispec)
-
- vx=veloc(1,iglob)
- vy=veloc(2,iglob)
- vz=veloc(3,iglob)
-
- nx=normal_xmax(1,j,k,ispec2D)
- ny=normal_xmax(2,j,k,ispec2D)
- nz=normal_xmax(3,j,k,ispec2D)
-
- vn=vx*nx+vy*ny+vz*nz
-
- tx=rho_vp(i,j,k,ispec)*vn*nx+rho_vs(i,j,k,ispec)*(vx-vn*nx)
- ty=rho_vp(i,j,k,ispec)*vn*ny+rho_vs(i,j,k,ispec)*(vy-vn*ny)
- tz=rho_vp(i,j,k,ispec)*vn*nz+rho_vs(i,j,k,ispec)*(vz-vn*nz)
-
- weight=jacobian2D_xmax(j,k,ispec2D)*wgllwgll_yz(j,k)
-
- accel(1,iglob)=accel(1,iglob) - tx*weight
- accel(2,iglob)=accel(2,iglob) - ty*weight
- accel(3,iglob)=accel(3,iglob) - tz*weight
-
- enddo
- enddo
- enddo
-
-!jpampuero ================ BEGIN DYNAMIC FAULT =========================
-!jpampuero I am implementing the FAULT boundary conditions in place of one
-!jpampuero of the absorbing boundaries.
-!jpampuero It is a quick and dirty implementation, for SCEC 1st benchmark test.
-!jpampuero No split nodes need to be implemented.
-!jpampuero The model is symmetric.
-!jpampuero Assume only ONE proc contains the fault: NPROC_XI=1
-!jpampuero Fault located at y=ymin
-!jpampuero
- if (FaultInThisProc) then
- do ispec2D=1,nspec2D_ymin
-
- ! exclude elements that are not on absorbing edges
- if(nkmin_eta(1,ispec2D) == 0 .or. nimin(1,ispec2D) == 0) cycle
-
- !ispec = bulk element associated to the ispec2D'th boundary element
- ispec=ibelm_ymin(ispec2D)
-
- !skip element if its center is not in rupture region
- ! iglob=ibool(3,1,3,ispec)
- ! if (-zstore(iglob)>RuptureDepth .OR. abs(xstore(iglob))>RuptureHalfLength) cycle
-
- j=1
- do k=nkmin_eta(1,ispec2D),NGLLZ
- do i=nimin(1,ispec2D),nimax(1,ispec2D)
-! do k=1,NGLLZ
-! do i=1,NGLLX
- iglob=ibool(i,j,k,ispec)
-
- ! Update strength :
- ! For symmetry reasons slip=2*displ
- FaultStrength(i,k,ispec2D) = max(FaultMuD(i,k,ispec2D), &
- FaultMuS(i,k,ispec2D)-FaultSWR(i,k,ispec2D)*sqrt(displ(1,iglob)**2+displ(3,iglob)**2) )
-
- ! From Newmark: v_new = v_pred +dt/2*a_new
- ! so: a_new = (v_new-v_pred)*2/dt = M^(-1) *( -K*d_new -B*t_new )
-
- !velocity assuming the fault is stress free during this timestep (t_new=0)
- !v_free = v_pred - dt/2 *M^(-1) *K*d_new
- FaultDVxFree = veloc(1,iglob)+deltatover2*rmass(1,iglob)*accel(1,iglob)
- FaultDVzFree = veloc(3,iglob)+deltatover2*rmass(3,iglob)*accel(3,iglob)
- !NOTE: I assume that only ONE proc contains the fault plane
- ! else accel=-K*d needs to be MPI assembled prior to compute V_free
-
- !"stick" or trial stress, assuming v_new = 0 during this timestep
- !T_trial = B^(-1) * M * 2/dt*V_free
- weight=FaultB(i,k,ispec2D) ! assembled
- tx = FaultDVxFree/(deltatover2*weight*rmass(1,iglob))
- tz = FaultDVzFree/(deltatover2*weight*rmass(3,iglob))
-
-! !VECTORIAL FRICTION VERSION
-! !test yield on magnitude of absolute stress (+initial)
-! txa = tx+FaultTxInit(i,k,ispec2D)
-! tza = tz+FaultTzInit(i,k,ispec2D)
-! ta = sqrt( txa*txa + tza*tza )
-! if ( ta>FaultStrength(i,k,ispec2D) ) then
-! va = sqrt(FaultDVxFree**2+ FaultDVzFree**2)
-! if ( va<VERYSMALLVAL ) then
-! dirx = txa/ta
-! dirz = tza/ta
-! else
-! dirx = FaultDVxFree/va
-! dirz = FaultDVzFree/va
-! endif
-! ta = FaultStrength(i,k,ispec2D)
-! tx = ta*dirx -FaultTxInit(i,k,ispec2D)
-! tz = ta*dirz -FaultTzInit(i,k,ispec2D)
-! endif
-
-! !SCALAR FRICTION VERSION
- tx = min(tx, FaultStrength(i,k,ispec2D) -FaultTxInit(i,k,ispec2D))
-
- FaultTx(i,k,ispec2D) = tx !store relative shear stress
- FaultTz(i,k,ispec2D) = tz !store relative shear stress
-
-
- enddo
- enddo
- enddo
-
- do ispec2D=1,nspec2D_ymin
- if(nkmin_eta(1,ispec2D) == 0 .or. nimin(1,ispec2D) == 0) cycle
- ispec=ibelm_ymin(ispec2D)
- j=1
- do k=nkmin_eta(1,ispec2D),NGLLZ
- do i=nimin(1,ispec2D),nimax(1,ispec2D)
- iglob=ibool(i,j,k,ispec)
- weight=jacobian2D_ymin(i,k,ispec2D)*wgllwgll_xz(i,k) ! not assembled
- accel(1,iglob)=accel(1,iglob) - FaultTx(i,k,ispec2D)*weight
- !accel(2,iglob): Stress-free on the Y direction (normal to the fault plane): do nothing
-! accel(3,iglob)=accel(3,iglob) - FaultTz(i,k,ispec2D)*weight ! VECTORIAL friction
- accel(3,iglob)=0._CUSTOM_REAL ! if SCALAR friction
- enddo
- enddo
- enddo
-
- endif
-
-!jpampuero ================ END DYNAMIC FAULT =========================
-
-! ymax
- do ispec2D=1,nspec2D_ymax
-
- ispec=ibelm_ymax(ispec2D)
-
-! exclude elements that are not on absorbing edges
- if(nkmin_eta(2,ispec2D) == 0 .or. nimin(2,ispec2D) == 0) cycle
-
- j=NGLLY
- do k=nkmin_eta(2,ispec2D),NGLLZ
- do i=nimin(2,ispec2D),nimax(2,ispec2D)
- iglob=ibool(i,j,k,ispec)
-
- vx=veloc(1,iglob)
- vy=veloc(2,iglob)
- vz=veloc(3,iglob)
-
- nx=normal_ymax(1,i,k,ispec2D)
- ny=normal_ymax(2,i,k,ispec2D)
- nz=normal_ymax(3,i,k,ispec2D)
-
- vn=vx*nx+vy*ny+vz*nz
-
- tx=rho_vp(i,j,k,ispec)*vn*nx+rho_vs(i,j,k,ispec)*(vx-vn*nx)
- ty=rho_vp(i,j,k,ispec)*vn*ny+rho_vs(i,j,k,ispec)*(vy-vn*ny)
- tz=rho_vp(i,j,k,ispec)*vn*nz+rho_vs(i,j,k,ispec)*(vz-vn*nz)
-
- weight=jacobian2D_ymax(i,k,ispec2D)*wgllwgll_xz(i,k)
-
- accel(1,iglob)=accel(1,iglob) - tx*weight
- accel(2,iglob)=accel(2,iglob) - ty*weight
- accel(3,iglob)=accel(3,iglob) - tz*weight
-
- enddo
- enddo
- enddo
-
-
-! bottom (zmin)
- do ispec2D=1,NSPEC2D_BOTTOM
-
- ispec=ibelm_bottom(ispec2D)
-
- k=1
- do j=1,NGLLY
- do i=1,NGLLX
-
- iglob=ibool(i,j,k,ispec)
-
- vx=veloc(1,iglob)
- vy=veloc(2,iglob)
- vz=veloc(3,iglob)
-
- nx=normal_bottom(1,i,j,ispec2D)
- ny=normal_bottom(2,i,j,ispec2D)
- nz=normal_bottom(3,i,j,ispec2D)
-
- vn=vx*nx+vy*ny+vz*nz
-
- tx=rho_vp(i,j,k,ispec)*vn*nx+rho_vs(i,j,k,ispec)*(vx-vn*nx)
- ty=rho_vp(i,j,k,ispec)*vn*ny+rho_vs(i,j,k,ispec)*(vy-vn*ny)
- tz=rho_vp(i,j,k,ispec)*vn*nz+rho_vs(i,j,k,ispec)*(vz-vn*nz)
-
- weight=jacobian2D_bottom(i,j,ispec2D)*wgllwgll_xy(i,j)
-
- accel(1,iglob)=accel(1,iglob) - tx*weight
- accel(2,iglob)=accel(2,iglob) - ty*weight
- accel(3,iglob)=accel(3,iglob) - tz*weight
-
- enddo
- enddo
- enddo
-
- endif ! end of Stacey conditions
-!jpampuero ************* END OF ABSORBING BOUNDARY CONDITIONS *********************
-
- do isource = 1,NSOURCES
-
-! add the source (only if this proc carries the source)
- if(myrank == islice_selected_source(isource)) then
-
- stf = comp_source_time_function(dble(it-1)*DT-hdur(isource)-t_cmt(isource),hdur(isource))
-
-! distinguish whether single or double precision for reals
- if(CUSTOM_REAL == SIZE_REAL) then
- stf_used = sngl(stf)
- else
- stf_used = stf
- endif
-
-! add source array
- do k=1,NGLLZ
- do j=1,NGLLY
- do i=1,NGLLX
- iglob = ibool(i,j,k,ispec_selected_source(isource))
- accel(:,iglob) = accel(:,iglob) + sourcearrays(isource,:,i,j,k)*stf_used
- enddo
- enddo
- enddo
-
- endif
-
- enddo
-
-! assemble all the contributions between slices using MPI
- call assemble_MPI_vector(accel,iproc_xi,iproc_eta,addressing, &
- iboolleft_xi,iboolright_xi,iboolleft_eta,iboolright_eta, &
- buffer_send_faces,buffer_received_faces,npoin2D_xi,npoin2D_eta, &
- NPROC_XI,NPROC_ETA,NPOIN2DMAX_XMIN_XMAX,NPOIN2DMAX_YMIN_YMAX,NPOIN2DMAX_XY)
-
- do i=1,NGLOB_AB
- accel(1,i) = accel(1,i)*rmass(1,i) !jpampuero ABC
- accel(2,i) = accel(2,i)*rmass(2,i) !jpampuero ABC
- accel(3,i) = accel(3,i)*rmass(3,i) !jpampuero ABC
-
- veloc(:,i) = veloc(:,i) + deltatover2*accel(:,i)
- enddo
-
-! write the seismograms with time shift
-
- do irec_local = 1,nrec_local
-
-! get global number of that receiver
- irec = number_receiver_global(irec_local)
-
-! perform the general interpolation using Lagrange polynomials
- uxd = ZERO
- uyd = ZERO
- uzd = ZERO
-
- vxd = ZERO
- vyd = ZERO
- vzd = ZERO
-
- do j = 1,NGLLY
- do i = 1,NGLLX
-
-! receivers are always located at the surface in the crustal region of the mesh
- iglob = ibool(i,j,NGLLZ,ispec_selected_rec(irec))
-
- hlagrange = hxir_store(irec_local,i)*hetar_store(irec_local,j)
-
-! save displacement
- uxd = uxd + dble(displ(1,iglob))*hlagrange
- uyd = uyd + dble(displ(2,iglob))*hlagrange
- uzd = uzd + dble(displ(3,iglob))*hlagrange
-
-! save velocity
- vxd = vxd + dble(veloc(1,iglob))*hlagrange
- vyd = vyd + dble(veloc(2,iglob))*hlagrange
- vzd = vzd + dble(veloc(3,iglob))*hlagrange
-
- enddo
- enddo
-
-! store North, East and Vertical components
-
-! distinguish whether single or double precision for reals
- if(CUSTOM_REAL == SIZE_REAL) then
- seismograms_d(:,irec_local,it) = sngl((nu(:,1,irec)*uxd + nu(:,2,irec)*uyd + nu(:,3,irec)*uzd))
- seismograms_v(:,irec_local,it) = sngl((nu(:,1,irec)*vxd + nu(:,2,irec)*vyd + nu(:,3,irec)*vzd))
- else
- seismograms_d(:,irec_local,it) = (nu(:,1,irec)*uxd + nu(:,2,irec)*uyd + nu(:,3,irec)*uzd)
- seismograms_v(:,irec_local,it) = (nu(:,1,irec)*vxd + nu(:,2,irec)*vyd + nu(:,3,irec)*vzd)
- endif
-
- enddo
-
-! write the current seismograms
- if(mod(it,NSEIS) == 0) then
- call write_seismograms_d(myrank,seismograms_d,number_receiver_global,station_name, &
- network_name,nrec,nrec_local,it,DT,NSTEP,hdur(1),LOCAL_PATH)
- !network_name,nrec,nrec_local,it,DT,NSTEP,hdur(1),'OUTPUT_FILES') !jpampuero
- call write_seismograms_v(myrank,seismograms_v,number_receiver_global,station_name, &
- network_name,nrec,nrec_local,it,DT,NSTEP,hdur(1),LOCAL_PATH)
- !network_name,nrec,nrec_local,it,DT,NSTEP,hdur(1),'OUTPUT_FILES') !jpampuero
- endif
-
-!jpampuero =========== BEGIN OUTPUT FAULT SNAPSHOTS ================
- if (FaultInThisProc) then
-
- !OUTPUT: time series at selected fault points
- FaultDisplX(:,it) = 2.0*displ(1,FaultIglobOut)
- FaultDisplZ(:,it) = 2.0*displ(3,FaultIglobOut)
- FaultVelocX(:,it) = 2.0*veloc(1,FaultIglobOut)
- FaultVelocZ(:,it) = 2.0*veloc(3,FaultIglobOut)
- do i=1,FaultOutputNPTS
- FaultStressX(i,it) = FaultTx(FaultIgllOut(i),FaultKgllOut(i),FaultElmtOut(i))
- FaultStressZ(i,it) = FaultTz(FaultIgllOut(i),FaultKgllOut(i),FaultElmtOut(i))
- enddo
-
- !OUTPUT: slip rate snapshot
- if (mod(it,FaultSnapshotNDT) == 0) then
- !write(FaultOutputName,'(A,'/Snapshot',I2.2,'_',I3.3,'.bin')') LOCAL_PATH,it/FaultSnapshotNDT,myrank
- write(FaultOutputName,"('OUTPUT_FILES/Snapshot',I0,'.bin')") it/FaultSnapshotNDT
- open(unit=IOUT, file= trim(FaultOutputName), status='replace', form='unformatted')
- do ispec2D=1,nspec2D_ymin
- ispec=ibelm_ymin(ispec2D)
- j=1
- !skip element if its center is not in rupture region
- iglob=ibool(3,j,3,ispec)
- if (-zstore(iglob)>RuptureDepth .OR. abs(xstore(iglob))>RuptureHalfLength) cycle
-
- do k=1,NGLLZ-1,2 ! watch out: do not repeat a node
- do i=1,NGLLX-1,2 ! assuming NGLL=5
- iglob=ibool(i,j,k,ispec)
- write(IOUT) xstore(iglob),zstore(iglob) &
- ,2.*veloc(1,iglob),2.*veloc(3,iglob) &
- ,2.*displ(1,iglob),2.*displ(3,iglob) &
- ,FaultTx(i,k,ispec2D),FaultTz(i,k,ispec2D)
- enddo
- enddo
- enddo
- close(IOUT)
-
- !intermediate output at selected fault points
- do i=1,FaultOutputNPTS
- open(IOUT,file='OUTPUT_FILES/'//FaultOutputPtsName(i),status='replace')
- write(IOUT,*) "# author=Jean-Paul Ampuero (am)"
- write(IOUT,*) "# date=2004/09/08"
- write(IOUT,*) "# code=SPECFEM3D (se)"
- write(IOUT,*) "# code_version=0904"
- write(IOUT,*) "# element_size=333.33 m (5 nodes per element edge)"
- write(IOUT,*) "# time_step=",DT
- write(IOUT,*) "# num_time_steps=",NSTEP
- write(IOUT,*) "# location=",FaultOutputPtsName(i)
- write(IOUT,*) "# Time series in 7 column of E14.6"
- write(IOUT,*) "# Column #1 = Time (s)"
- write(IOUT,*) "# Column #2 = horizontal right-lateral slip (m)"
- write(IOUT,*) "# Column #3 = horizontal right-lateral slip rate (m/s)"
- write(IOUT,*) "# Column #4 = horizontal right-lateral shear stress (MPa)"
- write(IOUT,*) "# Column #5 = vertical up-dip slip (m)"
- write(IOUT,*) "# Column #6 = vertical up-dip slip rate (m/s)"
- write(IOUT,*) "# Column #7 = vertical up-dip shear stress (MPa)"
- do k=1,NSTEP
- write(IOUT,'(7(E14.6))') k*DT,FaultDisplX(i,k),FaultVelocX(i,k),FaultStressX(i,k)/1e6 &
- ,FaultDisplZ(i,k),FaultVelocZ(i,k),FaultStressZ(i,k)/1e6
- enddo
- close(IOUT)
- enddo
-
- endif
-
- endif
-!jpampuero =========== END OUTPUT FAULT SNAPSHOTS ================
-
-!
-!---- end of time iteration loop
-!
- enddo ! end of main time loop
-
-! write the final seismograms
- call write_seismograms_d(myrank,seismograms_d,number_receiver_global,station_name, &
- network_name,nrec,nrec_local,it,DT,NSTEP,hdur(1),LOCAL_PATH)
- !network_name,nrec,nrec_local,it,DT,NSTEP,hdur(1),'OUTPUT_FILES') !jpampuero
- call write_seismograms_v(myrank,seismograms_v,number_receiver_global,station_name, &
- network_name,nrec,nrec_local,it,DT,NSTEP,hdur(1),LOCAL_PATH)
- !network_name,nrec,nrec_local,it,DT,NSTEP,hdur(1),'OUTPUT_FILES') !jpampuero
-
-!jpampuero === BEGIN OUTPUT: time series at selected fault points ===
- if (FaultInThisProc) then
- FaultStressX = FaultStressX/1e6
- FaultStressZ = FaultStressZ/1e6
- do i=1,FaultOutputNPTS
- open(IOUT,file='OUTPUT_FILES/'//FaultOutputPtsName(i),status='replace')
- write(IOUT,*) "# author=Jean-Paul Ampuero (am)"
- write(IOUT,*) "# date=2004/09/08"
- write(IOUT,*) "# code=SPECFEM3D (se)"
- write(IOUT,*) "# code_version=0904"
- write(IOUT,*) "# element_size=333.33 m (5 nodes per element edge)"
- write(IOUT,*) "# time_step=",DT
- write(IOUT,*) "# num_time_steps=",NSTEP
- write(IOUT,*) "# location=",FaultOutputPtsName(i)
- write(IOUT,*) "# Time series in 7 column of E14.6"
- write(IOUT,*) "# Column #1 = Time (s)"
- write(IOUT,*) "# Column #2 = horizontal right-lateral slip (m)"
- write(IOUT,*) "# Column #3 = horizontal right-lateral slip rate (m/s)"
- write(IOUT,*) "# Column #4 = horizontal right-lateral shear stress (MPa)"
- write(IOUT,*) "# Column #5 = vertical up-dip slip (m)"
- write(IOUT,*) "# Column #6 = vertical up-dip slip rate (m/s)"
- write(IOUT,*) "# Column #7 = vertical up-dip shear stress (MPa)"
- do it=1,NSTEP
- write(IOUT,'(7(E14.6))') it*DT,FaultDisplX(i,it),FaultVelocX(i,it),FaultStressX(i,it) &
- ,FaultDisplZ(i,it),FaultVelocZ(i,it),FaultStressZ(i,it)
- enddo
- close(IOUT)
- enddo
-! open(unit=IOUT,file='OUTPUT_FILES/FaultDispl.bin',status='replace',form='unformatted')
-! write(IOUT) FaultDispl
-! close(IOUT)
-! open(unit=IOUT,file='OUTPUT_FILES/FaultVeloc.bin',status='replace',form='unformatted')
-! write(IOUT) FaultVeloc
-! close(IOUT)
- endif
-!jpampuero === END OUTPUT: time series at selected fault points ===
-
-! close the main output file
- if(myrank == 0) then
- write(IMAIN,*)
- write(IMAIN,*) 'End of the simulation'
- write(IMAIN,*)
- close(IMAIN)
- endif
-
-! synchronize all the processes to make sure everybody has finished
- call MPI_BARRIER(MPI_COMM_WORLD,ier)
-
-! stop all the MPI processes, and exit
- call MPI_FINALIZE(ier)
-
- end program specfem3D
-
More information about the CIG-COMMITS
mailing list