[cig-commits] r11742 - seismo/2D/SPECFEM2D/trunk
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Thu Apr 3 04:52:06 PDT 2008
Author: dkomati1
Date: 2008-04-03 04:52:06 -0700 (Thu, 03 Apr 2008)
New Revision: 11742
Added:
seismo/2D/SPECFEM2D/trunk/compute_curl_one_element.f90
seismo/2D/SPECFEM2D/trunk/paco_beyond_critical.f90
seismo/2D/SPECFEM2D/trunk/paco_convolve_fft.f90
Modified:
seismo/2D/SPECFEM2D/trunk/Makefile
seismo/2D/SPECFEM2D/trunk/checkgrid.F90
seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90
seismo/2D/SPECFEM2D/trunk/specfem2D.F90
seismo/2D/SPECFEM2D/trunk/write_seismograms.F90
Log:
- added Ronan Madec's code to use an incident plane wave (analytical expression in the time domain if below the critical angle, Paco Sanchez-Sesma's analytical expression in the frequency domain and then an inverse FFT to convert to the time domain if above the critical angle)
- added Ronan Madec's code to save curl seismograms to compute rotations
- added Paco Sanchez-Sesma's routines for a plane wave
Modified: seismo/2D/SPECFEM2D/trunk/Makefile
===================================================================
--- seismo/2D/SPECFEM2D/trunk/Makefile 2008-04-03 11:35:00 UTC (rev 11741)
+++ seismo/2D/SPECFEM2D/trunk/Makefile 2008-04-03 11:52:06 UTC (rev 11742)
@@ -88,8 +88,9 @@
$O/specfem2D.o $O/write_seismograms.o $O/define_external_model.o $O/createnum_fast.o $O/createnum_slow.o\
$O/define_shape_functions.o $O/attenuation_model.o $O/create_color_image.o $O/compute_vector_field.o $O/compute_pressure.o\
$O/recompute_jacobian.o $O/compute_arrays_source.o $O/locate_source_moment_tensor.o $O/netlib_specfun_erf.o\
- $O/construct_acoustic_surface.o $O/assemble_MPI.o $O/compute_energy.o\
- $O/attenuation_compute_param.o $O/compute_Bielak_conditions.o
+ $O/construct_acoustic_surface.o $O/assemble_MPI.o $O/compute_energy.o $O/compute_curl_one_element.o\
+ $O/attenuation_compute_param.o $O/compute_Bielak_conditions.o $O/paco_beyond_critical.o\
+ $O/paco_convolve_fft.o
default: clean meshfem2D specfem2D convolve_source_timefunction
@@ -194,6 +195,9 @@
$O/compute_pressure.o: compute_pressure.f90 constants.h
${F90} $(FLAGS_CHECK) -c -o $O/compute_pressure.o compute_pressure.f90
+$O/compute_curl_one_element.o: compute_curl_one_element.f90 constants.h
+ ${F90} $(FLAGS_CHECK) -c -o $O/compute_curl_one_element.o compute_curl_one_element.f90
+
$O/compute_Bielak_conditions.o: compute_Bielak_conditions.f90 constants.h
${F90} $(FLAGS_CHECK) -c -o $O/compute_Bielak_conditions.o compute_Bielak_conditions.f90
@@ -226,3 +230,10 @@
$O/attenuation_compute_param.o: attenuation_compute_param.c
${CC} -c -o $O/attenuation_compute_param.o attenuation_compute_param.c
+
+$O/paco_beyond_critical.o: paco_beyond_critical.f90 constants.h
+ ${F90} $(FLAGS_CHECK) -c -o $O/paco_beyond_critical.o paco_beyond_critical.f90
+
+$O/paco_convolve_fft.o: paco_convolve_fft.f90 constants.h
+ ${F90} $(FLAGS_CHECK) -c -o $O/paco_convolve_fft.o paco_convolve_fft.f90
+
Modified: seismo/2D/SPECFEM2D/trunk/checkgrid.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/checkgrid.F90 2008-04-03 11:35:00 UTC (rev 11741)
+++ seismo/2D/SPECFEM2D/trunk/checkgrid.F90 2008-04-03 11:52:06 UTC (rev 11742)
@@ -120,7 +120,6 @@
double precision coorg(NDIM,npgeo)
-
! title of the plot
character(len=60) simulation_title
Added: seismo/2D/SPECFEM2D/trunk/compute_curl_one_element.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_curl_one_element.f90 (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/compute_curl_one_element.f90 2008-04-03 11:52:06 UTC (rev 11742)
@@ -0,0 +1,122 @@
+
+!========================================================================
+!
+! S P E C F E M 2 D Version 5.2
+! ------------------------------
+!
+! Copyright Universite de Pau et des Pays de l'Adour and CNRS, France.
+! Contributors: Dimitri Komatitsch, dimitri DOT komatitsch aT univ-pau DOT fr
+! Nicolas Le Goff, nicolas DOT legoff aT univ-pau DOT fr
+! Roland Martin, roland DOT martin aT univ-pau DOT fr
+!
+! This software is a computer program whose purpose is to solve
+! the two-dimensional viscoelastic anisotropic wave equation
+! using a spectral-element method (SEM).
+!
+! This software is governed by the CeCILL license under French law and
+! abiding by the rules of distribution of free software. You can use,
+! modify and/or redistribute the software under the terms of the CeCILL
+! license as circulated by CEA, CNRS and INRIA at the following URL
+! "http://www.cecill.info".
+!
+! As a counterpart to the access to the source code and rights to copy,
+! modify and redistribute granted by the license, users are provided only
+! with a limited warranty and the software's author, the holder of the
+! economic rights, and the successive licensors have only limited
+! liability.
+!
+! In this respect, the user's attention is drawn to the risks associated
+! with loading, using, modifying and/or developing or reproducing the
+! software by the user in light of its specific status of free software,
+! that may mean that it is complicated to manipulate, and that also
+! therefore means that it is reserved for developers and experienced
+! professionals having in-depth computer knowledge. Users are therefore
+! encouraged to load and test the software's suitability as regards their
+! requirements in conditions enabling the security of their systems and/or
+! data to be ensured and, more generally, to use and operate it in the
+! same conditions as regards security.
+!
+! The full text of the license is available in file "LICENSE".
+!
+!========================================================================
+
+subroutine compute_curl_one_element(curl_element,displ_elastic,elastic, &
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin, ispec)
+
+ ! compute curl in elastic elements (for rotational study)
+
+ implicit none
+
+ include "constants.h"
+
+ integer nspec,npoin,ispec
+
+ integer, dimension(NGLLX,NGLLX,nspec) :: ibool
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz
+
+ ! curl in this element
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: curl_element
+
+ logical, dimension(nspec) :: elastic
+ real(kind=CUSTOM_REAL), dimension(NDIM,npoin) :: displ_elastic
+
+ ! array with derivatives of Lagrange polynomials
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx
+ real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ) :: hprime_zz
+
+ ! local variables
+ integer :: i,j,k,iglob
+
+ ! jacobian
+ real(kind=CUSTOM_REAL) :: xixl,xizl,gammaxl,gammazl
+
+ ! spatial derivatives
+ real(kind=CUSTOM_REAL) :: dux_dxi,dux_dgamma,duz_dxi,duz_dgamma
+ real(kind=CUSTOM_REAL) :: duz_dxl,dux_dzl
+
+
+ if(elastic(ispec)) then
+
+ do j = 1,NGLLZ
+ do i = 1,NGLLX
+
+ ! derivative along x and along z
+ dux_dxi = ZERO
+ duz_dxi = ZERO
+
+ dux_dgamma = ZERO
+ duz_dgamma = ZERO
+
+ ! first double loop over GLL points to compute and store gradients
+ ! we can merge the two loops because NGLLX == NGLLZ
+ do k = 1,NGLLX
+ dux_dxi = dux_dxi + displ_elastic(1,ibool(k,j,ispec))*hprime_xx(i,k)
+ duz_dxi = duz_dxi + displ_elastic(2,ibool(k,j,ispec))*hprime_xx(i,k)
+ dux_dgamma = dux_dgamma + displ_elastic(1,ibool(i,k,ispec))*hprime_zz(j,k)
+ duz_dgamma = duz_dgamma + displ_elastic(2,ibool(i,k,ispec))*hprime_zz(j,k)
+ enddo
+
+ xixl = xix(i,j,ispec)
+ xizl = xiz(i,j,ispec)
+ gammaxl = gammax(i,j,ispec)
+ gammazl = gammaz(i,j,ispec)
+
+ ! derivatives of displacement
+ dux_dzl = dux_dxi*xizl + dux_dgamma*gammazl
+ duz_dxl = duz_dxi*xixl + duz_dgamma*gammaxl
+
+ ! store pressure
+ curl_element(i,j) = - 0.5d0 * (dux_dzl - duz_dxl)
+
+ enddo
+ enddo
+
+ else
+
+ call exit_MPI('no curl in acoustic')
+
+ endif ! end of test if acoustic or elastic element
+
+end subroutine compute_curl_one_element
+
Property changes on: seismo/2D/SPECFEM2D/trunk/compute_curl_one_element.f90
___________________________________________________________________
Name: svn:executable
+ *
Modified: seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90 2008-04-03 11:35:00 UTC (rev 11741)
+++ seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90 2008-04-03 11:52:06 UTC (rev 11742)
@@ -50,7 +50,9 @@
dux_dxl_np1,duz_dzl_np1,duz_dxl_np1,dux_dzl_np1,hprime_xx,hprimewgll_xx, &
hprime_zz,hprimewgll_zz,wxgll,wzgll,inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2,Mu_nu1,Mu_nu2,N_SLS, &
nspec_inner_outer,ispec_inner_outer_to_glob,num_phase_inner_outer,deltat,coord,add_Bielak_conditions, &
- x0_source, z0_source, A_plane, B_plane, C_plane, angleforce_refl, c_inc, c_refl, time_offset,f0)
+ x0_source, z0_source, A_plane, B_plane, C_plane, angleforce_refl, c_inc, c_refl, time_offset,f0, &
+ v0x_left,v0z_left,v0x_right,v0z_right,v0x_bot,v0z_bot,t0x_left,t0z_left,t0x_right,t0z_right,t0x_bot,t0z_bot,&
+ nleft,nright,nbot,over_critical_angle)
! compute forces for the elastic elements
@@ -131,6 +133,13 @@
double precision, dimension(NDIM,npoin), intent(in) :: coord
double precision x0_source, z0_source, angleforce_refl, c_inc, c_refl, time_offset, f0
double precision, dimension(NDIM) :: A_plane, B_plane, C_plane
+!over critical angle
+ logical :: over_critical_angle
+ integer :: nleft, nright, nbot
+ double precision, dimension(nleft) :: v0x_left,v0z_left,t0x_left,t0z_left
+ double precision, dimension(nright) :: v0x_right,v0z_right,t0x_right,t0z_right
+ double precision, dimension(nbot) :: v0x_bot,v0z_bot,t0x_bot,t0z_bot
+ integer count_left,count_right,count_bot
! only for the first call to compute_forces_elastic (during computation on outer elements)
if ( num_phase_inner_outer ) then
@@ -296,6 +305,10 @@
!
if(anyabs) then
+ count_left=1
+ count_right=1
+ count_bot=1
+
do ispecabs = 1,nelemabs
ispec = numabs(ispecabs)
@@ -310,6 +323,7 @@
!--- left absorbing boundary
if(codeabs(ILEFT,ispecabs)) then
+!!$ if(.false.) then
i = 1
@@ -320,11 +334,19 @@
! for analytical initial plane wave for Bielak's conditions
! left or right edge, horizontal normal vector
if(add_Bielak_conditions .and. initialfield) then
- call compute_Bielak_conditions(coord,iglob,npoin,it,deltat,dxUx,dxUz,dzUx,dzUz,veloc_horiz,veloc_vert, &
- x0_source, z0_source, A_plane, B_plane, C_plane, angleforce, angleforce_refl, &
- c_inc, c_refl, time_offset,f0)
- traction_x_t0 = (lambdal_relaxed+2*mul_relaxed)*dxUx + lambdal_relaxed*dzUz
- traction_z_t0 = mul_relaxed*(dxUz + dzUx)
+ if (.not.over_critical_angle) then
+ call compute_Bielak_conditions(coord,iglob,npoin,it,deltat,dxUx,dxUz,dzUx,dzUz,veloc_horiz,veloc_vert, &
+ x0_source, z0_source, A_plane, B_plane, C_plane, angleforce, angleforce_refl, &
+ c_inc, c_refl, time_offset,f0)
+ traction_x_t0 = (lambdal_relaxed+2*mul_relaxed)*dxUx + lambdal_relaxed*dzUz
+ traction_z_t0 = mul_relaxed*(dxUz + dzUx)
+ else
+ veloc_horiz=v0x_left(count_left)
+ veloc_vert=v0z_left(count_left)
+ traction_x_t0=t0x_left(count_left)
+ traction_z_t0=t0z_left(count_left)
+ count_left=count_left+1
+ end if
else
veloc_horiz = 0
veloc_vert = 0
@@ -370,6 +392,7 @@
!--- right absorbing boundary
if(codeabs(IRIGHT,ispecabs)) then
+!!$ if(.false.) then
i = NGLLX
@@ -380,11 +403,19 @@
! for analytical initial plane wave for Bielak's conditions
! left or right edge, horizontal normal vector
if(add_Bielak_conditions .and. initialfield) then
- call compute_Bielak_conditions(coord,iglob,npoin,it,deltat,dxUx,dxUz,dzUx,dzUz,veloc_horiz,veloc_vert, &
- x0_source, z0_source, A_plane, B_plane, C_plane, angleforce, angleforce_refl, &
- c_inc, c_refl, time_offset,f0)
- traction_x_t0 = (lambdal_relaxed+2*mul_relaxed)*dxUx + lambdal_relaxed*dzUz
- traction_z_t0 = mul_relaxed*(dxUz + dzUx)
+ if (.not.over_critical_angle) then
+ call compute_Bielak_conditions(coord,iglob,npoin,it,deltat,dxUx,dxUz,dzUx,dzUz,veloc_horiz,veloc_vert, &
+ x0_source, z0_source, A_plane, B_plane, C_plane, angleforce, angleforce_refl, &
+ c_inc, c_refl, time_offset,f0)
+ traction_x_t0 = (lambdal_relaxed+2*mul_relaxed)*dxUx + lambdal_relaxed*dzUz
+ traction_z_t0 = mul_relaxed*(dxUz + dzUx)
+ else
+ veloc_horiz=v0x_right(count_right)
+ veloc_vert=v0z_right(count_right)
+ traction_x_t0=t0x_right(count_right)
+ traction_z_t0=t0z_right(count_right)
+ count_right=count_right+1
+ end if
else
veloc_horiz = 0
veloc_vert = 0
@@ -446,11 +477,19 @@
! for analytical initial plane wave for Bielak's conditions
! top or bottom edge, vertical normal vector
if(add_Bielak_conditions .and. initialfield) then
- call compute_Bielak_conditions(coord,iglob,npoin,it,deltat,dxUx,dxUz,dzUx,dzUz,veloc_horiz,veloc_vert, &
- x0_source, z0_source, A_plane, B_plane, C_plane, angleforce, angleforce_refl, &
- c_inc, c_refl, time_offset,f0)
- traction_x_t0 = mul_relaxed*(dxUz + dzUx)
- traction_z_t0 = lambdal_relaxed*dxUx + (lambdal_relaxed+2*mul_relaxed)*dzUz
+ if (.not.over_critical_angle) then
+ call compute_Bielak_conditions(coord,iglob,npoin,it,deltat,dxUx,dxUz,dzUx,dzUz,veloc_horiz,veloc_vert, &
+ x0_source, z0_source, A_plane, B_plane, C_plane, angleforce, angleforce_refl, &
+ c_inc, c_refl, time_offset,f0)
+ traction_x_t0 = mul_relaxed*(dxUz + dzUx)
+ traction_z_t0 = lambdal_relaxed*dxUx + (lambdal_relaxed+2*mul_relaxed)*dzUz
+ else
+ veloc_horiz=v0x_bot(count_bot)
+ veloc_vert=v0z_bot(count_bot)
+ traction_x_t0=t0x_bot(count_bot)
+ traction_z_t0=t0z_bot(count_bot)
+ count_bot=count_bot+1
+ end if
else
veloc_horiz = 0
veloc_vert = 0
Added: seismo/2D/SPECFEM2D/trunk/paco_beyond_critical.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/paco_beyond_critical.f90 (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/paco_beyond_critical.f90 2008-04-03 11:52:06 UTC (rev 11742)
@@ -0,0 +1,620 @@
+!
+! This subroutine was written by Paco Sanchez-Sesma and his colleagues
+! from the Autonomous University of Mexico (UNAM), Mexico City, Mexico
+!
+! original name : DISTRAFF.f
+!
+! CALCULO DE DESPLAZAMIENTOS (UX, UZ) y TRACCIONES (TX, TZ) DE CAMPO LIBRE
+! EN UN SEMIESPACIO ELÁSTICO Y EN LA VECINDAD DE LA SUPERFICIE
+!
+! INCIDENCIA DE ONDAS P, SV Y DE RAYLEIGH
+!
+! 7 de febrero de 2007
+!
+! modified by Dimitri Komatitsch and Ronan Madec in March 2008
+! in particular, converted to Fortran90 and to double precision
+
+subroutine paco_beyond_critical(coord,npoin,deltat,NSTEP_global,angleforce,&
+ f0,cp_local,cs_local,INCLUDE_ATTENUATION,QD,source_type,v0x_left,v0z_left,v0x_right,v0z_right,&
+ v0x_bot,v0z_bot,t0x_left,t0z_left,t0x_right,t0z_right,t0x_bot,t0z_bot,left_bound,right_bound,&
+ bot_bound,nleft,nright,nbot,displ_elastic,veloc_elastic,accel_elastic)
+
+ implicit none
+
+ include "constants.h"
+
+ double precision :: f0,cp_local,cs_local,deltat,dt,TP,angleforce,QD,delta_in_period
+ logical :: INCLUDE_ATTENUATION
+ integer :: npt,NSTEP_global,source_type,nleft,nright,nbot,npoin
+
+ integer, dimension(nleft) :: left_bound
+ integer, dimension(nright) :: right_bound
+ integer, dimension(nbot) :: bot_bound
+
+ double precision, dimension(nleft,NSTEP_global) :: v0x_left,v0z_left, t0x_left,t0z_left
+ double precision, dimension(nright,NSTEP_global) :: v0x_right,v0z_right, t0x_right,t0z_right
+ double precision, dimension(nbot,NSTEP_global) :: v0x_bot,v0z_bot, t0x_bot,t0z_bot
+
+ double precision, dimension(2,npoin) :: coord
+ double precision, dimension(2,npoin) :: displ_elastic
+ double precision, dimension(2,npoin) :: veloc_elastic
+ double precision, dimension(2,npoin) :: accel_elastic
+
+ integer, dimension(:),allocatable :: local_pt
+
+ double precision, dimension(:), allocatable :: temp_field
+
+ integer :: J, indice, NSTEP_local, FLAG, N, NFREC, NFREC1
+
+ double precision :: ANU,BEALF,ALFBE,RLM,VNX,VNZ,A1,B1,TOTO,FJ,AKA,AQA,GAMR
+
+!location point
+ double precision :: X, Z, xmin, xmax, zmin, zmax
+ integer :: inode
+
+ double complex CAKA,CAQA,UI,UR
+ double complex UX,UZ,SX,SZ,SXZ,A2,B2,AL,AK,AM
+
+ double complex :: TX,TZ
+
+ double complex, dimension(:),allocatable::Field_Ux,Field_Uz,Field_Tx,Field_Tz
+
+ double precision :: TS
+
+!! move it to move the place where the wave reflects on free surface (offset too)
+ double precision :: offset
+
+!properties of the model
+ xmin=minval(coord(1,:))
+ xmax=maxval(coord(1,:))
+ zmin=minval(coord(2,:))
+ zmax=maxval(coord(2,:))
+
+!! DK DK offset of the origin of time of the Ricker (equivalent to t0 in SPECFEM2D)
+ offset=4.d0*(xmax-xmin)/5.d0
+ TS=2.d0/f0
+
+!! DK DK dominant period of the Ricker ((equivalent to 1/f0 in SPECFEM2D)
+ TP=1.d0/f0
+
+
+!!!find optimal period
+!!!if period is too small, you should see several initial plane wave on your initial field
+ delta_in_period=2.d0
+ do while(delta_in_period<1.5*abs(xmax-xmin)/cs_local)
+ delta_in_period=2.d0*delta_in_period
+ end do
+
+!!!testing of Deltat compatibility
+ DT=256.d0
+ do while(DT>deltat)
+ DT=DT/2.d0
+ end do
+ if (abs(DT-deltat)>1.0d-13) then
+ print *, "you must take a deltat as a power of two (power can be negative)"
+ print *, "for example you can take", DT
+ stop "can't go further, restart with new deltat"
+ end if
+
+ DT=deltat/2.d0
+
+ N=2
+ do while(N<2*NSTEP_global+1)
+ N=2.d0*N
+ end do
+
+ do while(DT<(delta_in_period/N))
+ N=2.d0*N
+ end do
+
+ print *, 'N found to do frequency calcul :', N
+ print *,'number of discrete frequencies = ',N/2
+ print *,'delta in period (seconds) = ',delta_in_period
+ print *,'delta in frequency (Hz) = ',1.d0/delta_in_period
+ print *,'dt (here we need deltat/2) = ', DT
+
+ NFREC=N/2
+ NFREC1=NFREC+1
+
+
+ !
+ ! FDT: FUNCION DE TRASFERENCIA
+ !
+
+!! calcul of Poisson's ratio
+ ANU = (cp_local*cp_local-2.d0*cs_local*cs_local)/(2.d0*(cp_local*cp_local-cs_local*cs_local))
+ print *,"Poisson's ratio = ",ANU
+
+ UI=(0.0d0, 1.0d0)
+ UR=(1.0d0, 0.0d0)
+
+!! DK DK convert angle to radians
+ GAMR = angleforce
+
+ BEALF=SQRT((1.0d0-2.0d0*ANU)/(2.0d0*(1.0d0-ANU)))
+ ALFBE=1.0d0/BEALF
+ RLM=ALFBE**2-2.0d0
+
+ !! RM flags: interior=0, left=1, right=2, bottom=3
+ do FLAG=0,3
+
+ if (FLAG==0) then
+ print *, "calcul of the initial field for every point of the mesh"
+ npt=npoin
+ allocate(local_pt(npt))
+ do inode=1,npt
+ local_pt(inode)=inode
+ end do
+ NSTEP_local=1
+ else if(FLAG==1) then
+ print *, "calcul of every time step on the left absorbing boundary"
+ npt=nleft
+ allocate(local_pt(npt))
+ local_pt=left_bound
+ NSTEP_local=NSTEP_global
+ else if(FLAG==2) then
+ print *, "calcul of every time step on the right absorbing boundary"
+ npt=nright
+ allocate(local_pt(npt))
+ local_pt=right_bound
+ NSTEP_local=NSTEP_global
+ else if(FLAG==3) then
+ print *, "calcul of every time step on the bottom absorbing boundary"
+ npt=nbot
+ allocate(local_pt(npt))
+ local_pt=bot_bound
+ NSTEP_local=NSTEP_global
+ end if
+
+
+ !to distinguish all model case and boundary case
+ allocate(temp_field(NSTEP_local))
+
+ allocate(Field_Ux(NFREC1))
+ allocate(Field_Uz(NFREC1))
+ allocate(Field_Tx(NFREC1))
+ allocate(Field_Tz(NFREC1))
+
+
+ if(mod(N,2) /= 0) stop 'N must be a multiple of 2'
+
+!! DK DK normal vector to the edge at this grid point
+!! DK DK therefore corners between two grid edges must be computed twice
+!! DK DK because the normal will change
+ if (FLAG==1) then
+ VNZ = 0.d0
+ VNX = 1.d0
+ else if (FLAG==2) then
+ VNZ = 0.d0
+ VNX = 1.d0
+ else if (FLAG==3) then
+ VNZ = 1.d0
+ VNX = 0.d0
+ else
+ VNZ = 0.d0
+ VNX = 0.d0
+ end if
+
+
+ do indice=1,npt
+
+ if (FLAG==0) then
+ inode=indice
+ X=coord(1,indice)-offset
+ !!!specfem is implemented from bottom to top whereas for this code
+ !!!we need from top to bottom
+ Z=zmax-coord(2,indice)
+ else
+ inode=local_pt(indice)
+ X=coord(1,inode)-offset
+ !!!specfem is implemented from bottom to top whereas for this code
+ !!!we need from top to bottom
+ Z=zmax-coord(2,inode)
+
+ end if
+
+ if (mod(indice,500)==0) then
+ print *, indice, "points have been treated on ",npt," total points"
+ end if
+
+ !
+ !! DK DK first handle the particular case of zero frequency
+ !
+ TOTO=0.01d0
+ IF (source_type==1) CALL ONDASP(GAMR,0.01d0*BEALF,A1,B1,A2,B2,AL,AK,AM,ANU,BEALF)
+ IF (source_type==2) CALL ONDASS(GAMR,TOTO,0.01d0*BEALF,A1,B1,A2,B2,AL,AK,AM,ANU,BEALF)
+ IF (source_type==3) CALL ONDASR(0.01d0*BEALF,A1,B1,A2,B2,AL,AK,AM,ANU,BEALF)
+
+
+ TOTO=0.0d0
+ CALL DESFXY(TOTO,TOTO,source_type,UX,UZ,SX,SZ,SXZ,A1,B1,A2,B2,AL,AK,AM,RLM,ANU)
+
+ !! DK DK write the frequency seismograms
+ TX = SX *VNX+SXZ*VNZ
+ TZ = SXZ*VNX+SZ *VNZ
+
+ Field_Ux(1)=UX
+ Field_Uz(1)=UZ
+ if (FLAG/=0) then
+ Field_Tx(1)=TX
+ Field_Tz(1)=TZ
+ end if
+
+
+ !
+ ! then loop on all the other discrete frequencies
+ !
+ do J=1,N/2
+
+ !! DK DK compute the value of the frequency (= index * delta in frequency = index * 1/delta in period)
+ FJ = dble(J) * 1.d0 / delta_in_period
+
+ !! DK DK pulsation (= 2 * PI * frequency)
+ AKA=2.0d0*PI*FJ
+
+ AQA=AKA*BEALF
+
+ !! DK DK exclude attenuation completely if needed
+ if(INCLUDE_ATTENUATION) then
+ CAKA=CMPLX(AKA,-AKA/(2.0d0*QD))
+ CAQA=CMPLX(AQA,-AQA/(2.0d0*QD))
+ else
+ CAKA=CMPLX(AKA,0)
+ CAQA=CMPLX(AQA,0)
+ endif
+
+ IF (source_type==1) CALL ONDASP(GAMR,AQA,A1,B1,A2,B2,AL,AK,AM,ANU,BEALF)
+ IF (source_type==2) CALL ONDASS(GAMR,AKA,AQA,A1,B1,A2,B2,AL,AK,AM,ANU,BEALF)
+ IF (source_type==3) CALL ONDASR(AQA,A1,B1,A2,B2,AL,AK,AM,ANU,BEALF)
+
+ CALL DESFXY(X,Z,source_type,UX,UZ,SX,SZ,SXZ,A1,B1,A2,B2,AL,AK,AM,RLM,ANU)
+
+ !! DK DK write the frequency seismograms
+ TX = SX *VNX+SXZ*VNZ
+ TZ = SXZ*VNX+SZ *VNZ
+
+ Field_Ux(J+1)=UX
+ Field_Uz(J+1)=UZ
+ if (FLAG/=0) then
+ Field_Tx(J+1)=TX
+ Field_Tz(J+1)=TZ
+ end if
+
+ enddo
+
+ !to convert frequency field in time field
+ !(number at the end are unit numbers for writing in the good file,
+ !in the case of the traction we fill only one file per call)
+
+ !global model case for initial field
+ if (FLAG==0) then
+ call paco_convolve_fft(Field_Ux,1,NSTEP_local,dt,NFREC,temp_field,TP,TS)
+ displ_elastic(1,indice)=temp_field(1)
+ call paco_convolve_fft(Field_Uz,1,NSTEP_local,dt,NFREC,temp_field,TP,TS)
+ displ_elastic(2,indice)=temp_field(1)
+ call paco_convolve_fft(Field_Ux,2,NSTEP_local,dt,NFREC,temp_field,TP,TS)
+ veloc_elastic(1,indice)=temp_field(1)
+ call paco_convolve_fft(Field_Uz,2,NSTEP_local,dt,NFREC,temp_field,TP,TS)
+ veloc_elastic(2,indice)=temp_field(1)
+ call paco_convolve_fft(Field_Ux,3,NSTEP_local,dt,NFREC,temp_field,TP,TS)
+ accel_elastic(1,indice)=temp_field(1)
+ call paco_convolve_fft(Field_Uz,3,NSTEP_local,dt,NFREC,temp_field,TP,TS)
+ accel_elastic(2,indice)=temp_field(1)
+
+ !absorbing boundaries...
+ !left case
+ else if (FLAG==1) then
+ call paco_convolve_fft(Field_Ux,2,NSTEP_local,dt,NFREC,temp_field,TP,TS)
+ v0x_left(indice,:)=temp_field(:)
+ call paco_convolve_fft(Field_Uz,2,NSTEP_local,dt,NFREC,temp_field,TP,TS)
+ v0z_left(indice,:)=temp_field(:)
+ call paco_convolve_fft(Field_Tx,4,NSTEP_local,dt,NFREC,temp_field,TP,TS)
+ t0x_left(indice,:)=temp_field(:)
+ call paco_convolve_fft(Field_Tz,4,NSTEP_local,dt,NFREC,temp_field,TP,TS)
+ t0z_left(indice,:)=temp_field(:)
+
+ !right case
+ else if (FLAG==2) then
+ call paco_convolve_fft(Field_Ux,2,NSTEP_local,dt,NFREC,temp_field,TP,TS)
+ v0x_right(indice,:)=temp_field(:)
+ call paco_convolve_fft(Field_Uz,2,NSTEP_local,dt,NFREC,temp_field,TP,TS)
+ v0z_right(indice,:)=temp_field(:)
+ call paco_convolve_fft(Field_Tx,4,NSTEP_local,dt,NFREC,temp_field,TP,TS)
+ t0x_right(indice,:)=temp_field(:)
+ call paco_convolve_fft(Field_Tz,4,NSTEP_local,dt,NFREC,temp_field,TP,TS)
+ t0z_right(indice,:)=temp_field(:)
+
+ !bottom case
+ else if (FLAG==3) then
+ call paco_convolve_fft(Field_Ux,2,NSTEP_local,dt,NFREC,temp_field,TP,TS)
+ v0x_bot(indice,:)=temp_field(:)
+ call paco_convolve_fft(Field_Uz,2,NSTEP_local,dt,NFREC,temp_field,TP,TS)
+ v0z_bot(indice,:)=temp_field(:)
+ call paco_convolve_fft(Field_Tx,4,NSTEP_local,dt,NFREC,temp_field,TP,TS)
+ t0x_bot(indice,:)=temp_field(:)
+ call paco_convolve_fft(Field_Tz,4,NSTEP_local,dt,NFREC,temp_field,TP,TS)
+ t0z_bot(indice,:)=temp_field(:)
+ end if
+ enddo
+
+ deallocate(temp_field)
+ deallocate(local_pt)
+
+ deallocate(Field_Ux)
+ deallocate(Field_Uz)
+ deallocate(Field_Tx)
+ deallocate(Field_Tz)
+
+ end do
+
+end subroutine paco_beyond_critical
+
+!---
+
+SUBROUTINE DESFXY(X,Z,ICAS,UX,UZ,SX,SZ,SXZ,A1,B1,A2,B2,AL,AK,AM,RLM,ANU)
+
+ implicit none
+
+ double precision A1,B1,RLM,ANU,X,Z
+ integer ICAS
+ double complex UX,UZ,SX,SZ,SXZ,A2,B2,AL,AK,AM
+ double complex UI,FAC
+ double complex AUX1,AUX2,FI1,FI2,PS1,PS2
+
+
+ UI=(0.0d0,1.0d0)
+ if (A1/=0.0d0) then
+ AUX1=A1*CDEXP(UI*(AM*Z-AL*X)) !campo P incidente
+ else
+ AUX1=CMPLX(0.0d0)
+ end if
+ if (A2/=0.0d0) then
+ AUX2=A2*CDEXP(-UI*(AM*Z+AL*X)) *1.0d0 !campo P reflejado
+ else
+ AUX2=CMPLX(0.0d0)
+ end if
+ FI1=AUX1+AUX2
+ FI2=AUX1-AUX2
+ if (B1/=0.0d0) then
+ AUX1=B1*CDEXP(UI*(AK*Z-AL*X)) !campo S incidente
+ else
+ AUX1=CMPLX(0.0d0)
+ end if
+ if (B2/=0.0d0) then
+ AUX2=B2*CDEXP(-UI*(AK*Z+AL*X)) *1.0d0 !campo S reflejado
+ else
+ AUX2=CMPLX(0.0d0)
+ end if
+ PS1=AUX1+AUX2
+ PS2=AUX1-AUX2
+
+!
+! FAC ES PARA TENER CONSISTENCIA CON AKI & RICHARDS
+!
+ FAC=UI
+ IF (ICAS==2)FAC=-UI
+
+ UX=(-UI*AL*FI1+UI*AK*PS2)*FAC
+
+ UZ=(UI*AM*FI2+UI*AL*PS1)*FAC
+!!!! DK DK Paco's convention is inverted
+ UZ = - UZ
+
+ AUX1=AL*AL+AM*AM
+ SX=(-RLM*AUX1*FI1-2.0d0*AL*(AL*FI1-AK*PS2))*FAC
+ SZ=(-RLM*AUX1*FI1-2.0d0*(AM*AM*FI1+AK*AL*PS2))*FAC
+
+ SXZ=(2.0d0*AM*AL*FI2+(AL*AL-AK*AK)*PS1)*FAC
+!!!! DK DK Paco's convention is inverted
+ SXZ = - SXZ
+
+END SUBROUTINE DESFXY
+
+SUBROUTINE FAFB(CA,CB,FA,FB)
+
+ implicit none
+
+ double precision CA,CB,A,B
+ double complex FA,FB,ZER,UI
+
+ ZER=(0.0d0,0.0d0)
+ UI=(0.0d0,1.0d0)
+ A=CA*CA-1.0d0
+ B=CB*CB-1.0d0
+
+ IF (CA<1.0d0) then
+ FA=-UI*SQRT(-A)
+ else
+ FA=SQRT(A)+ZER
+ end IF
+
+ IF (CB<1.0d0) then
+ FB=-UI*SQRT(-B)
+ else
+ FB=CMPLX(SQRT(B),0.0d0)
+ end IF
+
+ RETURN
+
+END SUBROUTINE FAFB
+
+SUBROUTINE A2B2(FA,FB,A2,B2)
+
+ implicit none
+
+ double complex FA,FB,A2,B2,DEN,AUX
+
+ AUX=FB*FB-1.0d0
+ DEN=4.0d0*FA*FB+AUX*AUX
+ A2=(4.0d0*FA*FB-AUX*AUX)/DEN
+ B2=4.0d0*FA*AUX/DEN
+
+END SUBROUTINE A2B2
+
+!! DK DK calculation of P waves
+SUBROUTINE ONDASP(GP,AQB,A1,B1,A2,B2,AL,AK,AM,ANU,BEALF)
+
+ implicit none
+
+ double precision A1,B1,ANU,CA,CB,GP,AQB,BEALF
+ double complex A2,B2,FA,FB,ZER,AL,AK,AM
+
+ ZER=(0.0d0,0.0d0)
+ BEALF=SQRT((1.0d0-2.0d0*ANU)/2.0d0/(1.0d0-ANU))
+ A1=1.0d0/AQB
+ B1=0.0d0
+
+ IF (GP==0.0d0) then
+
+ AL=ZER
+ AK=ZER
+ AM=AQB+ZER
+ A2=(-1.0d0+ZER)/AQB
+ B2=ZER
+ RETURN
+
+ end IF
+
+ CA=1.0d0/SIN(GP)
+ CB=CA/BEALF
+ AL=AQB/CA+ZER
+ CALL FAFB(CA,CB,FA,FB)
+ AK=AL*FB
+ AM=AL*FA
+ CALL A2B2(FA,FB,A2,B2)
+ A2=A2/AQB
+ B2=B2/AQB
+
+END SUBROUTINE ONDASP
+
+!! DK DK calculation of S waves
+SUBROUTINE ONDASS(GS,AKB,AQB,A1,B1,A2,B2,AL,AK,AM,ANU,BEALF)
+
+ implicit none
+
+ double precision A1,B1,ANU,CA,CB,GS,AQB,BEALF,AKB
+ double complex A2,B2,FA,FB,ZER,AL,AK,AM
+
+ ZER=(0.0d0,0.0d0)
+ BEALF=SQRT((1.0d0-2.0d0*ANU)/2.0d0/(1.0d0-ANU))
+ A1=0.0d0
+ B1=1.0d0/AKB
+
+ IF (GS==0.0d0) then
+
+ AL=ZER
+ AK=AKB+ZER
+ AM=ZER
+ A2=ZER
+ B2=(-1.0d0+ZER)/AKB
+ return
+
+ end IF
+
+ CB=1.0d0/SIN(GS)
+ CA=CB*BEALF
+
+!
+!! DK DK case of the critical angle
+!
+!! DK DK changed that to be safe IF (CA==1.0d0)THEN
+ IF (CA==1.d0) then!>= 0.9999d0)THEN
+ AL=AQB+ZER
+ AM=ZER
+ CALL FAFB(CA,CB,FA,FB)
+ AK=AL*FB
+ B2=-B1
+ A2=-4.0d0*COS(GS)*B1/(1./BEALF-2.*BEALF)
+
+!! DK DK case of an angle that is not critical
+ ELSE
+ AL=AQB/CA+ZER
+ CALL FAFB(CA,CB,FA,FB)
+ AK=AL*FB
+ AM=AL*FA
+ CALL A2B2(FA,FB,B2,A2)
+ A2=-A2*FB/FA
+ A2=A2/AKB
+ B2=B2/AKB
+ endif
+
+ RETURN
+
+END SUBROUTINE ONDASS
+
+!! DK DK calculation of Rayleigh waves
+SUBROUTINE ONDASR(AQB,A1,B1,A2,B2,AL,AK,AM,ANU,BEALF)
+
+ implicit none
+
+ double precision A1,B1,ANU,CA,CB,AQB,BEALF,ba2
+ double complex A2,B2,FA,FB,ZER,AL,AK,AM
+
+ double precision, external :: crb
+
+ ZER=(0.0d0,0.0d0)
+ A1=0.0d0
+ B1=0.0d0
+ B2=1.0d0+ZER
+ BEALF=SQRT((1.0d0-2.0d0*ANU)/2.0d0/(1.0d0-ANU))
+ BA2=BEALF*BEALF
+ CB=CRB(BEALF)
+ CA=CB*BEALF
+ AL=AQB/CA+ZER
+
+ CALL FAFB(CA,CB,FA,FB)
+
+ AK=AL*FB
+ AM=AL*FA
+ A2=2.0d0*FB/(FB*FB-1.0d0)*B2
+ B2=B2/(AL*A2+AK)
+ A2=A2*B2
+
+END SUBROUTINE ONDASR
+
+FUNCTION CRB(BEALF)
+
+ implicit none
+
+ include "constants.h"
+
+ double precision U3,BA2,P,Q,FIND,F1,F2,F12,FACT,CRB,BEALF
+
+ U3=1.0d0/3.0d0
+ BA2=BEALF*BEALF
+ P=8.0d0/3.0d0-16.0d0*BA2
+ Q=272.0d0/27.0d0-80.0d0/3.0d0*BA2
+ FIND=Q*Q/4.0d0+P*P*P/27.0d0
+ IF (FIND>=0.0d0) then
+ F1=SQRT(FIND)-Q/2.0d0
+ IF (F1>0.0d0) then
+ F1=F1**U3
+ else
+ F1=-(-F1)**U3
+ end IF
+ F2=-SQRT(FIND)-Q/2.0d0
+ IF (F2>0.0d0) then
+ F2=F2**U3
+ else
+ F2=-(-F2)**U3
+ end IF
+ FACT=F1+F2+8.0d0/3.0d0
+ CRB=SQRT(FACT)
+ return
+ else
+ F1=-27.0d0*Q*Q/(4.0d0*P*P*P)
+ F1=SQRT(F1)
+ IF (Q<0.0d0) then
+ F1=COS((PI-ACOS(F1))/3.0d0)
+ else
+ F1=COS(ACOS(F1)/3.0d0)
+ end IF
+ F2=-P/3.0d0
+ F2=SQRT(F2)
+ F12=-2.0d0*F1*F2+8.0d0/3.0d0
+ CRB=SQRT(F12)
+ RETURN
+ end IF
+
+END FUNCTION CRB
+
Property changes on: seismo/2D/SPECFEM2D/trunk/paco_beyond_critical.f90
___________________________________________________________________
Name: svn:executable
+ *
Added: seismo/2D/SPECFEM2D/trunk/paco_convolve_fft.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/paco_convolve_fft.f90 (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/paco_convolve_fft.f90 2008-04-03 11:52:06 UTC (rev 11742)
@@ -0,0 +1,237 @@
+!
+! This subroutine was written by Paco Sanchez-Sesma and his colleagues
+! from the Autonomous University of Mexico (UNAM), Mexico City, Mexico
+!
+! PROGRAMA PARA CALCULAR SISMOGRAMAS SINTETICOS DADA LA
+! FUNCION DE TRANSFERENCIA PARA COMPONENTES Ux, Uz, R2
+! Tx y Tz SOLUCION DE CAMPO LIBRE Caso P-SV, RAYLEIGH
+!
+! modified by Dimitri Komatitsch and Ronan Madec in March 2008
+! in particular, converted to Fortran90 and to double precision
+
+subroutine paco_convolve_fft(Field,label,NSTEP,dt,NFREC,output_field,tp,ts)
+
+ implicit none
+
+ integer :: NFREC,N,NSTEP
+
+ double complex, dimension(NFREC+1) :: Field
+
+ double complex CR(2*NFREC)
+
+ double precision, dimension(NSTEP) :: output_field
+
+ integer :: J,i,label
+
+ integer :: beg_UV, end_UV
+
+ double precision :: AN,FUN,RAIZ,X,Z,dt,tp,ts
+
+ double precision, external :: RIC, deRIC, de2RIC
+
+ integer inode
+
+ N=2*NFREC
+
+ AN = N
+
+ !!!
+ !!!label=1 <=> champ U en entree =>convolution par un ricker pour U
+ !!!label=2 <=> champ U en entree =>convolution par la derivee de ricker pour V
+ !!!label=3 <=> champ U en entree =>convolution par la derivee seconde de ricker pour A
+ !!!label=4 <=> champ T en entree =>convolution par un ricker
+ !!!
+ !!!flag=0 on a besoin de U, V et A (pas T)
+ !!!flag/=0 on a besoin de T et V (pas U ni A)
+ !!!
+ !!!NSTEP==1 <=> FLAG==0 (flags: interior=0, left=1, right=2, bottom=3)
+ !!!
+
+
+ do j=1,N
+ if (label==1 .or. label==4) FUN=ric(j,tp,ts,dt)
+ if (label==2) FUN=deric(j,tp,ts,dt)
+ if (label==3) FUN=de2ric(j,tp,ts,dt)
+ CR(j)=CMPLX(FUN,0.0d0)
+ enddo
+
+ CALL Transform(N,CR,-1.0d0)
+
+ RAIZ = SQRT(AN)
+
+ CALL SINTER(Field,output_field,NSTEP,CR,RAIZ,NFREC,label,dt)
+
+END subroutine paco_convolve_fft
+
+SUBROUTINE SINTER(V,output_field,NSTEP,CR,RAIZ,NFREC,label,dt)
+
+ implicit none
+
+ integer NSTEP, j,jn,nar,N,jstep,jbegin,jend,label,nfrec,mult,delay
+
+ double precision :: RAIZ
+
+ double complex VC
+
+ double precision VT(2*NFREC)
+
+ double precision :: filt,eta,dt
+
+
+ double precision, dimension(NSTEP) :: output_field
+
+ double complex, dimension(NFREC+1) :: V
+
+ double complex CY(2*NFREC),CR(2*NFREC)
+
+ N=2*NFREC
+
+ CY(1) = CR(1) * V(1) * RAIZ * dt
+
+ DO J=2,N/2+1
+ FILT = 1.0d0
+ VC = V(J)
+ CY(J)= CR(J)*VC * RAIZ * dt/ FILT
+ JN = N-J+2
+ CY(JN)=CONJG(CY(J))
+ enddo
+
+ CALL Transform(N,CY,1.0d0)
+
+ if (label==1.or.label==3.or.(label==2.and.NSTEP==1)) then
+ !coefs to take wanted time steps (t=0 : first time step)
+ mult=1
+ delay=0
+ else if(label==2.and.NSTEP>1) then
+ !coefs to take wanted time steps (t=i*deltat+1/2 : one step on two starting at 1/2)
+ mult=2
+ delay=0
+ else if(label==4) then
+ !coefs to take wanted time steps (t=i*deltat+1 : one step on two starting at 1)
+ mult=2
+ delay=1
+ end if
+
+
+ do J=1,NSTEP
+ CY(mult*J+delay)=CY(mult*J+delay)/RAIZ/dt
+ VT(mult*J+delay)=REAL(CY(mult*J+delay))
+ output_field(J)=VT(mult*J+delay)
+ enddo
+
+END SUBROUTINE SINTER
+
+!
+!! DK DK Ricker time function
+!
+FUNCTION RIC(J,tp,ts,dt)
+
+ implicit none
+
+ include "constants.h"
+
+ double precision :: A,RIC,tp,ts,dt
+
+ integer j
+
+ A=PI*(dt*(J-1)-ts)/tp
+ A=A*A
+ RIC=0.0d0
+ IF(A>30.0d0) RETURN
+ RIC=(A-0.5)*EXP(-A)
+
+END FUNCTION RIC
+
+!
+!! RM time derivative Ricker time function
+!
+FUNCTION deRIC(J,tp,ts,dt)
+
+ implicit none
+
+ include "constants.h"
+
+ double precision :: A,A_dot,deRIC,tp,ts,dt
+ integer :: label,j
+
+ A=PI*(dt*(J-1)-ts)/tp
+ A=A*A
+ A_dot=2*(PI/tp)**2*(dt*(J-1)-ts)
+ deRIC=0.0d0
+ IF(A>30.0d0) RETURN
+ deRIC=A_dot*(1.5-A)*EXP(-A)
+
+END FUNCTION deRIC
+
+!
+!! RM second time derivative Ricker time function
+!
+FUNCTION de2RIC(J,tp,ts,dt)
+
+ implicit none
+
+ include "constants.h"
+
+ double precision :: A,A_dot,A_dot_dot,de2RIC,tp,ts,dt
+ integer j
+
+ A=PI*(dt*(J-1)-ts)/tp
+ A=A*A
+ A_dot=2*(PI/tp)**2*(dt*(J-1)-ts)
+ A_dot_dot=2*(PI/tp)**2
+ de2RIC=0.0d0
+ IF(A>30.0d0) RETURN
+ de2RIC=(A_dot_dot*(1.5-A)-A_dot*A_dot-A_dot*(1.5-A)*A_dot)*EXP(-A)
+
+END FUNCTION de2RIC
+
+
+!! DK DK Fourier transform
+SUBROUTINE Transform(LX,CX,SIGNI)
+
+ implicit none
+
+ include "constants.h"
+
+ integer LX,i,j,l,istep,m
+
+ double precision SC
+
+ double complex CX(LX),CARG,CW,CTEMP
+
+ double precision SIGNI
+
+ J=1
+ SC=SQRT(1.0d0/LX)
+ DO I=1,LX
+ IF (I<=J) then
+ CTEMP=CX(J)*SC
+ CX(J)=CX(I)*SC
+ CX(I)=CTEMP
+ end IF
+ M=LX/2
+ do while (M>=1 .and. M<J)
+ J=J-M
+ M=M/2
+ end do
+ J=J+M
+ end DO
+ L=1
+
+ do while(L<LX)
+ ISTEP=2*L
+ DO M=1,L
+ CARG=(0.0d0,1.0d0)*(PI*SIGNI*(M-1))/L
+ CW=CDEXP(CARG)
+ DO I=M,LX,ISTEP
+ CTEMP=CW*CX(I+L)
+ CX(I+L)=CX(I)-CTEMP
+ CX(I)=CX(I)+CTEMP
+ end DO
+ end DO
+
+ L=ISTEP
+ end do
+
+END SUBROUTINE Transform
+
Property changes on: seismo/2D/SPECFEM2D/trunk/paco_convolve_fft.f90
___________________________________________________________________
Name: svn:executable
+ *
Modified: seismo/2D/SPECFEM2D/trunk/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2008-04-03 11:35:00 UTC (rev 11741)
+++ seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2008-04-03 11:52:06 UTC (rev 11742)
@@ -151,7 +151,7 @@
character(len=150) dummystring
! for seismograms
- double precision, dimension(:,:), allocatable :: sisux,sisuz
+ double precision, dimension(:,:), allocatable :: sisux,sisuz,siscurl
integer :: seismo_offset, seismo_current
! vector field in an element
@@ -160,9 +160,12 @@
! pressure in an element
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: pressure_element
+! curl in an element
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: curl_element
+
integer :: i,j,k,l,it,irec,ipoin,ip,id,nbpoin,inump,n,ispec,npoin,npgeo,iglob
logical :: anyabs
- double precision :: dxd,dzd,valux,valuz,hlagrange,rhol,cosrot,sinrot,xi,gamma,x,z
+ double precision :: dxd,dzd,dcurld,valux,valuz,valcurl,hlagrange,rhol,cosrot,sinrot,xi,gamma,x,z
! coefficients of the explicit Newmark time scheme
integer NSTEP
@@ -359,6 +362,15 @@
double precision, dimension(2) :: A_plane, B_plane, C_plane
double precision :: PP, PS, SP, SS, z0_source, x0_source, xmax, xmin, zmax, zmin, time_offset
+!over critical angle
+ integer , dimension(:), allocatable :: left_bound,right_bound,bot_bound
+ double precision , dimension(:,:), allocatable :: v0x_left,v0z_left,v0x_right,v0z_right,v0x_bot,v0z_bot
+ double precision , dimension(:,:), allocatable :: t0x_left,t0z_left,t0x_right,t0z_right,t0x_bot,t0z_bot
+ integer count_left,count_right,count_bot,itime,ibegin,iend
+!!$ double precision :: a1,b1,a2,b2,a3,b3,sin0,cos0,Zr,Zr_dot,Zr_dot_dot,Zi,Zi_dot,Zi_dot_dot,Cte
+ double precision :: X_temp,Z_temp,ignore
+ logical :: over_critical_angle
+
!***********************************************************************
!
! i n i t i a l i z a t i o n p h a s e
@@ -445,7 +457,7 @@
read(IIN,"(a80)") datlin
read(IIN,*) seismotype,imagetype
- if(seismotype < 1 .or. seismotype > 4) call exit_MPI('Wrong type for seismogram output')
+ if(seismotype < 1 .or. seismotype > 5) call exit_MPI('Wrong type for seismogram output')
if(imagetype < 1 .or. imagetype > 4) call exit_MPI('Wrong type for snapshots')
read(IIN,"(a80)") datlin
@@ -680,6 +692,7 @@
!
!---- read absorbing boundary data
!
+
read(IIN,"(a80)") datlin
if(anyabs) then
do inum = 1,nelemabs
@@ -888,17 +901,19 @@
!--- save the grid of points in a file
!
if(outputgrid) then
- write(IOUT,*)
- write(IOUT,*) 'Saving the grid in a text file...'
- write(IOUT,*)
- open(unit=55,file='OUTPUT_FILES/grid_points_and_model.txt',status='unknown')
- write(55,*) npoin
- do n = 1,npoin
- write(55,*) (coord(i,n), i=1,NDIM)
- enddo
- close(55)
+ write(IOUT,*)
+ write(IOUT,*) 'Saving the grid in a text file...'
+ write(IOUT,*)
+ open(unit=55,file='OUTPUT_FILES/grid_points_and_model.txt',status='unknown')
+ zmax=maxval(coord(2,:))
+ write(55,*) npoin
+ do n = 1,npoin
+ write(55,*) (coord(i,n), i=1,NDIM)
+ enddo
+ close(55)
endif
+
!
!----- plot the GLL mesh in a Gnuplot file
!
@@ -1010,7 +1025,7 @@
call compute_arrays_source(ispec_selected_source,xi_source,gamma_source,sourcearray, &
Mxx,Mzz,Mxz,xix,xiz,gammax,gammaz,xigll,zigll,nspec)
- else
+ else if(.not.initialfield) then
call exit_MPI('incorrect source type')
endif
@@ -1022,6 +1037,7 @@
! allocate seismogram arrays
allocate(sisux(NTSTEP_BETWEEN_OUTPUT_SEISMO,nrecloc))
allocate(sisuz(NTSTEP_BETWEEN_OUTPUT_SEISMO,nrecloc))
+ allocate(siscurl(NTSTEP_BETWEEN_OUTPUT_SEISMO,nrecloc))
! check if acoustic receiver is exactly on the free surface because pressure is zero there
do ispec_acoustic_surface = 1,nelem_acoustic_surface
@@ -1248,6 +1264,7 @@
if(any_elastic) rmass_inverse_elastic(:) = 1._CUSTOM_REAL / rmass_inverse_elastic(:)
if(any_acoustic) rmass_inverse_acoustic(:) = 1._CUSTOM_REAL / rmass_inverse_acoustic(:)
+
! check the mesh, stability and number of points per wavelength
call checkgrid(vpext,vsext,rhoext,density,elastcoef,ibool,kmato,coord,npoin,vpmin,vpmax, &
assign_external_model,nspec,numat,deltat,f0,t0,initialfield,time_function_type, &
@@ -1466,34 +1483,17 @@
!
!---- read initial fields from external file if needed
!
+
+! if we are looking a plane wave over critical angle we use other method
+ over_critical_angle=.false.
+
if(initialfield) then
write(IOUT,*)
write(IOUT,*) 'Reading initial fields from external file...'
write(IOUT,*)
if(any_acoustic) call exit_MPI('initial field currently implemented for purely elastic simulation only')
- if(.not. add_Bielak_conditions) then
- open(unit=55,file='OUTPUT_FILES/wavefields.txt',status='unknown')
- read(55,*) nbpoin
- if(nbpoin /= npoin) call exit_MPI('Wrong number of points in input file')
- allocate(displread(NDIM))
- allocate(velocread(NDIM))
- allocate(accelread(NDIM))
- do n = 1,npoin
- read(55,*) inump, (displread(i), i=1,NDIM), (velocread(i), i=1,NDIM), (accelread(i), i=1,NDIM)
- if(inump<1 .or. inump>npoin) call exit_MPI('Wrong point number')
- displ_elastic(:,inump) = displread
- veloc_elastic(:,inump) = velocread
- accel_elastic(:,inump) = accelread
- enddo
- deallocate(displread)
- deallocate(velocread)
- deallocate(accelread)
- close(55)
-
- else
-
!!$! compute analytical initial plane wave field
!!$! the analytical expression below is specific to an SV wave at 30 degrees and Poisson = 0.3333
!!$ print *,'computing analytical initial plane wave field for SV wave at 30 degrees and Poisson = 0.3333'
@@ -1535,16 +1535,23 @@
!
!=======================================================================
- print *,'Number of grid points: ',npoin
+ print *,'Number of grid points: ',npoin,'\n'
print *,'*** calculation of initial plane wave ***'
if (source_type == 1) then
print *,'initial P wave of', angleforce*180.d0/pi, 'degrees introduced...'
else if (source_type == 2) then
print *,'initial SV wave of', angleforce*180.d0/pi, ' degrees introduced...'
+
+ else if (source_type == 3) then
+ print *,'Rayleigh wave introduced...'
else
call exit_MPI('Unrecognized source_type: should be 1 for plane P waves, 2 for plane SV waves!')
endif
+ if ((angleforce<0.0d0.or.angleforce>=pi/2.d0).and.source_type/=3) then
+ call exit_MPI("bad angleforce: 0<=angleforce<90")
+ end if
+
! only implemented for homogeneous media therefore only 1 material supported
if (numat==1) then
@@ -1562,7 +1569,7 @@
c_inc = cploc
c_refl = csloc
- angleforce_refl = asin(p*csloc)
+ angleforce_refl = asin(p*c_refl)
! from formulas (5.26) and (5.27) p 140 in Aki & Richards (1980)
PP = (- cos(2.d0*angleforce_refl)**2/csloc**3 + 4.d0*p**2*cos(angleforce)*cos(angleforce_refl)/cploc) / &
@@ -1581,15 +1588,15 @@
C_plane(1) = PS * cos(angleforce_refl); C_plane(2) = PS * sin(angleforce_refl)
! SV wave case
- else
+ else if (source_type == 2) then
p=sin(angleforce)/csloc
c_inc = csloc
c_refl = cploc
! if this coefficient is greater than 1, we are beyond the critical SV wave angle and there cannot be a converted P wave
- if (p*cploc<=1.d0) then
- angleforce_refl = asin(p*cploc)
+ if (p*c_refl<=1.d0) then
+ angleforce_refl = asin(p*c_refl)
! from formulas (5.30) and (5.31) p 140 in Aki & Richards (1980)
SS = (cos(2.d0*angleforce_refl)**2/csloc**3 - 4.d0*p**2*cos(angleforce)*cos(angleforce_refl)/cploc) / &
@@ -1600,8 +1607,16 @@
print *,'reflected convert plane wave angle: ', angleforce_refl*180.d0/pi, '\n'
+ !SV45 degree incident plane wave is a particular case
+ else if (angleforce>pi/4.d0-1.0d-11 .and. angleforce<pi/4.d0+1.0d-11) then
+ angleforce_refl = 0.d0
+ SS = -1.0d0
+ SP = 0.d0
else
- call exit_MPI('cannot be included for now: SV angle too high, beyond critical angle')
+ over_critical_angle=.true.
+ angleforce_refl = 0.d0
+ SS = 0.0d0
+ SP = 0.d0
endif
! from Table 5.1 p141 in Aki & Richards (1980)
@@ -1610,10 +1625,15 @@
B_plane(1) = SS * cos(angleforce); B_plane(2) = SS * sin(angleforce)
C_plane(1) = SP * sin(angleforce_refl); C_plane(2) = - SP * cos(angleforce_refl)
+ !Rayleigh case
+ else if (source_type == 3) then
+ over_critical_angle=.true.
+ A_plane(1)=0.d0;A_plane(2)=0.d0
+ B_plane(1)=0.d0;B_plane(2)=0.d0
+ C_plane(1)=0.d0;C_plane(2)=0.d0;
endif
-
else
- call exit_MPI('not possible for now to have several materials with a plane wave (but could be done one day)')
+ call exit_MPI('not possible to have several materials with a plane wave')
endif
! get minimum and maximum values of mesh coordinates
@@ -1622,55 +1642,155 @@
xmax = maxval(coord(1,:))
zmax = maxval(coord(2,:))
- ! initialize the time offset to put the plane wave not too close to the free surface topography
- if (abs(angleforce)<20.d0*pi/180.d0) then
- time_offset=-1.d0*zmax/3.d0/c_inc
+ !initializing of the time offset to put the planewave not to close of the irregularity on the free surface
+ if (abs(angleforce)<10.d0*pi/180.d0 .and. source_type/=3) then
+ time_offset=-1.d0*zmax/2.d0/c_inc
+
+!!$ elseif (abs(angleforce)<31.d0*pi/180.d0 .and. source_type==1) then
+!!$ time_offset=-1.d0*zmax/3.d0/c_inc
+!!$ elseif (abs(angleforce)<29.d0*pi/180.d0 .and. source_type==2) then
+!!$ time_offset=-1.d0*zmax/3.d0/c_inc
+
else
time_offset=0.d0
endif
! to correctly center the initial plane wave in the mesh
z0_source=zmax
- x0_source=xmin + 1.d0*(xmax-xmin)/3.d0
+!!$ if (abs(angleforce)<1.d0*pi/180.d0) then
+!!$ x0_source=xmin + 1.d0*(xmax-xmin)/2.d0
+!!$ else
+!!$ x0_source=xmin + 1.d0*(xmax-xmin)/3.d0
+!!$
+!!$ endif
- do i = 1,npoin
+!!$ x0_source=x_source
+ x0_source=xmin + 1.d0*(xmax-xmin)/4.d0
- x = coord(1,i)
- z = coord(2,i)
+ if (.not.over_critical_angle) then
- ! z is from bottom to top therefore we take -z to make parallel with Aki & Richards
- z = z0_source - z
- x = x - x0_source
+ do i = 1,npoin
- t = 0.d0 + time_offset
+ x = coord(1,i)
+ z = coord(2,i)
- ! formulas for the initial displacement for a plane wave from Aki & Richards (1980)
- displ_elastic(1,i) = A_plane(1) * ricker_Bielak_displ(t - sin(angleforce)*x/c_inc + cos(angleforce)*z/c_inc,f0) &
- + B_plane(1) * ricker_Bielak_displ(t - sin(angleforce)*x/c_inc - cos(angleforce)*z/c_inc,f0) &
- + C_plane(1) * ricker_Bielak_displ(t - sin(angleforce_refl)*x/c_refl - cos(angleforce_refl)*z/c_refl,f0)
- displ_elastic(2,i) = A_plane(2) * ricker_Bielak_displ(t - sin(angleforce)*x/c_inc + cos(angleforce)*z/c_inc,f0) &
- + B_plane(2) * ricker_Bielak_displ(t - sin(angleforce)*x/c_inc - cos(angleforce)*z/c_inc,f0) &
- + C_plane(2) * ricker_Bielak_displ(t - sin(angleforce_refl)*x/c_refl - cos(angleforce_refl)*z/c_refl,f0)
+ ! z is from bottom to top therefore we take -z to make parallel with Aki & Richards
+ z = z0_source - z
+ x = x - x0_source
- ! formulas for the initial velocity for a plane wave (first derivative in time of the displacement)
- veloc_elastic(1,i) = A_plane(1) * ricker_Bielak_veloc(t - sin(angleforce)*x/c_inc + cos(angleforce)*z/c_inc,f0) &
- + B_plane(1) * ricker_Bielak_veloc(t - sin(angleforce)*x/c_inc - cos(angleforce)*z/c_inc,f0) &
- + C_plane(1) * ricker_Bielak_veloc(t - sin(angleforce_refl)*x/c_refl - cos(angleforce_refl)*z/c_refl,f0)
- veloc_elastic(2,i) = A_plane(2) * ricker_Bielak_veloc(t - sin(angleforce)*x/c_inc + cos(angleforce)*z/c_inc,f0) &
- + B_plane(2) * ricker_Bielak_veloc(t - sin(angleforce)*x/c_inc - cos(angleforce)*z/c_inc,f0) &
- + C_plane(2) * ricker_Bielak_veloc(t - sin(angleforce_refl)*x/c_refl - cos(angleforce_refl)*z/c_refl,f0)
+ t = 0.d0 + time_offset
- ! formulas for the initial acceleration for a plane wave (second derivative in time of the displacement)
- accel_elastic(1,i) = A_plane(1) * ricker_Bielak_accel(t - sin(angleforce)*x/c_inc + cos(angleforce)*z/c_inc,f0) &
- + B_plane(1) * ricker_Bielak_accel(t - sin(angleforce)*x/c_inc - cos(angleforce)*z/c_inc,f0) &
- + C_plane(1) * ricker_Bielak_accel(t - sin(angleforce_refl)*x/c_refl - cos(angleforce_refl)*z/c_refl,f0)
- accel_elastic(2,i) = A_plane(2) * ricker_Bielak_accel(t - sin(angleforce)*x/c_inc + cos(angleforce)*z/c_inc,f0) &
- + B_plane(2) * ricker_Bielak_accel(t - sin(angleforce)*x/c_inc - cos(angleforce)*z/c_inc,f0) &
- + C_plane(2) * ricker_Bielak_accel(t - sin(angleforce_refl)*x/c_refl - cos(angleforce_refl)*z/c_refl,f0)
- enddo
- endif ! add_Bielak
+ ! formulas for the initial displacement for a plane wave from Aki & Richards (1980)
+ displ_elastic(1,i) = A_plane(1) * ricker_Bielak_displ(t - sin(angleforce)*x/c_inc + cos(angleforce)*z/c_inc,f0) &
+ + B_plane(1) * ricker_Bielak_displ(t - sin(angleforce)*x/c_inc - cos(angleforce)*z/c_inc,f0) &
+ + C_plane(1) * ricker_Bielak_displ(t - sin(angleforce_refl)*x/c_refl - cos(angleforce_refl)*z/c_refl,f0)
+ displ_elastic(2,i) = A_plane(2) * ricker_Bielak_displ(t - sin(angleforce)*x/c_inc + cos(angleforce)*z/c_inc,f0) &
+ + B_plane(2) * ricker_Bielak_displ(t - sin(angleforce)*x/c_inc - cos(angleforce)*z/c_inc,f0) &
+ + C_plane(2) * ricker_Bielak_displ(t - sin(angleforce_refl)*x/c_refl - cos(angleforce_refl)*z/c_refl,f0)
+ ! formulas for the initial velocity for a plane wave (first derivative in time of the displacement)
+ veloc_elastic(1,i) = A_plane(1) * ricker_Bielak_veloc(t - sin(angleforce)*x/c_inc + cos(angleforce)*z/c_inc,f0) &
+ + B_plane(1) * ricker_Bielak_veloc(t - sin(angleforce)*x/c_inc - cos(angleforce)*z/c_inc,f0) &
+ + C_plane(1) * ricker_Bielak_veloc(t - sin(angleforce_refl)*x/c_refl - cos(angleforce_refl)*z/c_refl,f0)
+ veloc_elastic(2,i) = A_plane(2) * ricker_Bielak_veloc(t - sin(angleforce)*x/c_inc + cos(angleforce)*z/c_inc,f0) &
+ + B_plane(2) * ricker_Bielak_veloc(t - sin(angleforce)*x/c_inc - cos(angleforce)*z/c_inc,f0) &
+ + C_plane(2) * ricker_Bielak_veloc(t - sin(angleforce_refl)*x/c_refl - cos(angleforce_refl)*z/c_refl,f0)
+
+ ! formulas for the initial acceleration for a plane wave (second derivative in time of the displacement)
+ accel_elastic(1,i) = A_plane(1) * ricker_Bielak_accel(t - sin(angleforce)*x/c_inc + cos(angleforce)*z/c_inc,f0) &
+ + B_plane(1) * ricker_Bielak_accel(t - sin(angleforce)*x/c_inc - cos(angleforce)*z/c_inc,f0) &
+ + C_plane(1) * ricker_Bielak_accel(t - sin(angleforce_refl)*x/c_refl - cos(angleforce_refl)*z/c_refl,f0)
+ accel_elastic(2,i) = A_plane(2) * ricker_Bielak_accel(t - sin(angleforce)*x/c_inc + cos(angleforce)*z/c_inc,f0) &
+ + B_plane(2) * ricker_Bielak_accel(t - sin(angleforce)*x/c_inc - cos(angleforce)*z/c_inc,f0) &
+ + C_plane(2) * ricker_Bielak_accel(t - sin(angleforce_refl)*x/c_refl - cos(angleforce_refl)*z/c_refl,f0)
+
+ end do
+
+ else !over critical angle
+
+ if (source_type/=3) then
+ print *, '\n You are over critical angle ( > ',asin(c_inc/c_refl)*180d0/pi,')'
+ end if
+
+ print *, '\n *************'
+ print *, 'We have to calcul initial field in the frequency domain'
+ print *, 'and then pass it in the time domain(can be long... be patient...)'
+ print *, '*************\n'
+
+!!$ print *, "if first reflection is not at a good place, change TS and offset in paco_beyond_critical.f90 and recompute\n"
+
+ allocate(left_bound(nelemabs*NGLLX))
+ allocate(right_bound(nelemabs*NGLLX))
+ allocate(bot_bound(nelemabs*NGLLZ))
+
+ count_bot=0
+ count_left=0
+ count_right=0
+ do ispecabs=1,nelemabs
+ ispec=numabs(ispecabs)
+ if(codeabs(ILEFT,ispecabs)) then
+ i = 1
+ do j = 1,NGLLZ
+ count_left=count_left+1
+ iglob = ibool(i,j,ispec)
+ left_bound(count_left)=iglob
+ end do
+ end if
+ if(codeabs(IRIGHT,ispecabs)) then
+ i = NGLLX
+ do j = 1,NGLLZ
+ count_right=count_right+1
+ iglob = ibool(i,j,ispec)
+ right_bound(count_right)=iglob
+ end do
+ end if
+ if(codeabs(IBOTTOM,ispecabs)) then
+ j = 1
+ ! exclude corners to make sure there is no contradiction on the normal
+ ibegin = 1
+ iend = NGLLX
+ if(codeabs(ILEFT,ispecabs)) ibegin = 2
+ if(codeabs(IRIGHT,ispecabs)) iend = NGLLX-1
+ do i = ibegin,iend
+ count_bot=count_bot+1
+ iglob = ibool(i,j,ispec)
+ bot_bound(count_bot)=iglob
+ end do
+ end if
+ end do
+
+ allocate(v0x_left(count_left,NSTEP))
+ allocate(v0z_left(count_left,NSTEP))
+ allocate(t0x_left(count_left,NSTEP))
+ allocate(t0z_left(count_left,NSTEP))
+
+ allocate(v0x_right(count_right,NSTEP))
+ allocate(v0z_right(count_right,NSTEP))
+ allocate(t0x_right(count_right,NSTEP))
+ allocate(t0z_right(count_right,NSTEP))
+
+ allocate(v0x_bot(count_bot,NSTEP))
+ allocate(v0z_bot(count_bot,NSTEP))
+ allocate(t0x_bot(count_bot,NSTEP))
+ allocate(t0z_bot(count_bot,NSTEP))
+
+ !call paco routine to do calcul in frequency and pass in time by fourier transform
+ call paco_beyond_critical(coord,npoin,deltat,NSTEP,angleforce,&
+ f0,cploc,csloc,TURN_ATTENUATION_ON,Qp_attenuation,source_type,v0x_left,v0z_left,&
+ v0x_right,v0z_right,v0x_bot,v0z_bot,t0x_left,t0z_left,t0x_right,t0z_right,&
+ t0x_bot,t0z_bot,left_bound(:count_left),right_bound(:count_right),bot_bound(:count_bot)&
+ ,count_left,count_right,count_bot,displ_elastic,veloc_elastic,accel_elastic)
+
+
+
+ print *, '\n ***********'
+ print *, 'initial field calculated, starting propagation'
+ print *, '***********\n\n'
+
+
+ endif !over critical angle
+
write(IOUT,*) 'Max norm of initial elastic displacement = ',maxval(sqrt(displ_elastic(1,:)**2 + displ_elastic(2,:)**2))
endif ! initialfield
@@ -1903,7 +2023,7 @@
! if acoustic absorbing element and acoustic/elastic coupled element is the same
if(ispec_acoustic == ispec) then
- if(iedge_acoustic == IBOTTOM) then
+ if(iedge_acoustic == IBOTTOM) then
jbegin_left(ispecabs) = 2
jbegin_right(ispecabs) = 2
endif
@@ -2212,7 +2332,10 @@
dux_dxl_np1,duz_dzl_np1,duz_dxl_np1,dux_dzl_np1,hprime_xx,hprimewgll_xx, &
hprime_zz,hprimewgll_zz,wxgll,wzgll,inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2,Mu_nu1,Mu_nu2,N_SLS, &
nspec_outer, ispec_outer_to_glob,.true.,deltat,coord,add_Bielak_conditions, x0_source, z0_source, &
- A_plane, B_plane, C_plane, angleforce_refl, c_inc, c_refl, time_offset, f0)
+ A_plane, B_plane, C_plane, angleforce_refl, c_inc, c_refl, time_offset, f0,&
+ v0x_left(:,it),v0z_left(:,it),v0x_right(:,it),v0z_right(:,it),v0x_bot(:,it),v0z_bot(:,it), &
+ t0x_left(:,it),t0z_left(:,it),t0x_right(:,it),t0z_right(:,it),t0x_bot(:,it),t0z_bot(:,it), &
+ count_left,count_right,count_bot,over_critical_angle)
! *********************************************************
! ************* add coupling with the acoustic side
@@ -2304,8 +2427,12 @@
dux_dxl_np1,duz_dzl_np1,duz_dxl_np1,dux_dzl_np1,hprime_xx,hprimewgll_xx, &
hprime_zz,hprimewgll_zz,wxgll,wzgll,inv_tau_sigma_nu1,phi_nu1,inv_tau_sigma_nu2,phi_nu2,Mu_nu1,Mu_nu2,N_SLS, &
nspec_inner, ispec_inner_to_glob,.false.,deltat,coord,add_Bielak_conditions, x0_source, z0_source, &
- A_plane, B_plane, C_plane, angleforce_refl, c_inc, c_refl, time_offset, f0)
+ A_plane, B_plane, C_plane, angleforce_refl, c_inc, c_refl, time_offset, f0,&
+ v0x_left(:,it),v0z_left(:,it),v0x_right(:,it),v0z_right(:,it),v0x_bot(:,it),v0z_bot(:,it), &
+ t0x_left(:,it),t0z_left(:,it),t0x_right(:,it),t0z_right(:,it),t0x_bot(:,it),t0z_bot(:,it), &
+ count_left,count_right,count_bot,over_critical_angle)
+
! assembling accel_elastic for elastic elements (receive)
#ifdef USE_MPI
if ( nproc > 1 .and. any_elastic .and. ninterface_elastic > 0) then
@@ -2420,14 +2547,14 @@
! loop on all the receivers to compute and store the seismograms
do irecloc = 1,nrecloc
- irec = recloc(irecloc)
+ irec = recloc(irecloc)
ispec = ispec_selected_rec(irec)
! compute pressure in this element if needed
if(seismotype == 4) then
- call compute_pressure_one_element(pressure_element,potential_dot_dot_acoustic,displ_elastic,elastic, &
+ call compute_pressure_one_element(pressure_element,potential_dot_dot_acoustic,displ_elastic,elastic, &
xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,assign_external_model, &
numat,kmato,elastcoef,vpext,vsext,rhoext,ispec,e1,e11, &
TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,Mu_nu1,Mu_nu2,N_SLS)
@@ -2435,22 +2562,26 @@
else if(.not. elastic(ispec)) then
! for acoustic medium, compute vector field from gradient of potential for seismograms
- if(seismotype == 1) then
- call compute_vector_one_element(vector_field_element,potential_acoustic,displ_elastic,elastic, &
+ if(seismotype == 1) then
+ call compute_vector_one_element(vector_field_element,potential_acoustic,displ_elastic,elastic, &
xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,ispec,numat,kmato,density,rhoext,assign_external_model)
- else if(seismotype == 2) then
- call compute_vector_one_element(vector_field_element,potential_dot_acoustic,veloc_elastic,elastic, &
+ else if(seismotype == 2) then
+ call compute_vector_one_element(vector_field_element,potential_dot_acoustic,veloc_elastic,elastic, &
xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,ispec,numat,kmato,density,rhoext,assign_external_model)
- else if(seismotype == 3) then
- call compute_vector_one_element(vector_field_element,potential_dot_dot_acoustic,accel_elastic,elastic, &
+ else if(seismotype == 3) then
+ call compute_vector_one_element(vector_field_element,potential_dot_dot_acoustic,accel_elastic,elastic, &
xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin,ispec,numat,kmato,density,rhoext,assign_external_model)
- endif
+ endif
+ else if(seismotype == 5) then
+ call compute_curl_one_element(curl_element,displ_elastic,elastic, &
+ xix,xiz,gammax,gammaz,ibool,hprime_xx,hprime_zz,nspec,npoin, ispec)
endif
! perform the general interpolation using Lagrange polynomials
valux = ZERO
valuz = ZERO
+ valcurl = ZERO
do j = 1,NGLLZ
do i = 1,NGLLX
@@ -2459,6 +2590,8 @@
hlagrange = hxir_store(irec,i)*hgammar_store(irec,j)
+ dcurld=ZERO
+
if(seismotype == 4) then
dxd = pressure_element(i,j)
@@ -2484,11 +2617,18 @@
dxd = accel_elastic(1,iglob)
dzd = accel_elastic(2,iglob)
+ else if(seismotype == 5) then
+
+ dxd = displ_elastic(1,iglob)
+ dzd = displ_elastic(2,iglob)
+ dcurld = curl_element(i,j)
+
endif
! compute interpolated field
valux = valux + dxd*hlagrange
valuz = valuz + dzd*hlagrange
+ valcurl = valcurl + dcurld*hlagrange
enddo
enddo
@@ -2501,8 +2641,9 @@
sisux(seismo_current,irecloc) = valux
sisuz(seismo_current,irecloc) = ZERO
endif
+ siscurl(seismo_current,irecloc) = valcurl
- enddo
+ enddo
!
!---- display results at given time steps
@@ -2699,7 +2840,7 @@
endif
!---- save temporary or final seismograms
- call write_seismograms(sisux,sisuz,station_name,network_name,NSTEP, &
+ call write_seismograms(sisux,sisuz,siscurl,station_name,network_name,NSTEP, &
nrecloc,which_proc_receiver,nrec,myrank,deltat,seismotype,st_xval,t0, &
NTSTEP_BETWEEN_OUTPUT_SEISMO,seismo_offset,seismo_current &
)
@@ -2734,6 +2875,27 @@
enddo ! end of the main time loop
+
+ if (over_critical_angle) then
+ deallocate(v0x_left)
+ deallocate(v0z_left)
+ deallocate(t0x_left)
+ deallocate(t0z_left)
+ deallocate(left_bound)
+
+ deallocate(v0x_right)
+ deallocate(v0z_right)
+ deallocate(t0x_right)
+ deallocate(t0z_right)
+ deallocate(right_bound)
+
+ deallocate(v0x_bot)
+ deallocate(v0z_bot)
+ deallocate(t0x_bot)
+ deallocate(t0z_bot)
+ deallocate(bot_bound)
+ end if
+
!---- close energy file and create a gnuplot script to display it
if(OUTPUT_ENERGY) then
close(IENERGY)
Modified: seismo/2D/SPECFEM2D/trunk/write_seismograms.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/write_seismograms.F90 2008-04-03 11:35:00 UTC (rev 11741)
+++ seismo/2D/SPECFEM2D/trunk/write_seismograms.F90 2008-04-03 11:52:06 UTC (rev 11742)
@@ -42,7 +42,7 @@
! write seismograms to text files
- subroutine write_seismograms(sisux,sisuz,station_name,network_name, &
+ subroutine write_seismograms(sisux,sisuz,siscurl,station_name,network_name, &
NSTEP,nrecloc,which_proc_receiver,nrec,myrank,deltat,seismotype,st_xval,t0, &
NTSTEP_BETWEEN_OUTPUT_SEISMO,seismo_offset,seismo_current &
)
@@ -61,7 +61,7 @@
integer, intent(in) :: nrecloc,myrank
integer, dimension(nrec),intent(in) :: which_proc_receiver
- double precision, dimension(NTSTEP_BETWEEN_OUTPUT_SEISMO,nrecloc), intent(in) :: sisux,sisuz
+ double precision, dimension(NTSTEP_BETWEEN_OUTPUT_SEISMO,nrecloc), intent(in) :: sisux,sisuz,siscurl
double precision st_xval(nrec)
@@ -101,6 +101,8 @@
component = 'a'
else if(seismotype == 4) then
component = 'p'
+ else if(seismotype == 5) then
+ component = 'c'
else
call exit_MPI('wrong component to save for seismograms')
endif
@@ -109,6 +111,8 @@
! only one seismogram if pressurs
if(seismotype == 4) then
number_of_components = 1
+ else if(seismotype == 5) then
+ number_of_components = NDIM+1
else
number_of_components = NDIM
endif
@@ -137,6 +141,12 @@
open(unit=11,file='OUTPUT_FILES/Uz_file_double.bin',status='unknown')
close(11,status='delete')
+ open(unit=11,file='OUTPUT_FILES/Curl_file_single.bin',status='unknown')
+ close(11,status='delete')
+
+ open(unit=11,file='OUTPUT_FILES/Curl_file_double.bin',status='unknown')
+ close(11,status='delete')
+
endif
if ( myrank == 0 ) then
@@ -154,13 +164,20 @@
open(unit=13,file='OUTPUT_FILES/Ux_file_double.bin',status='unknown',access='direct',recl=8)
endif
-! no Z component seismogram if pressurs
+! no Z component seismogram if pressure
if(seismotype /= 4) then
open(unit=14,file='OUTPUT_FILES/Uz_file_single.bin',status='unknown',access='direct',recl=4)
open(unit=15,file='OUTPUT_FILES/Uz_file_double.bin',status='unknown',access='direct',recl=8)
end if
+! curl output
+ if(seismotype == 5) then
+ open(unit=16,file='OUTPUT_FILES/Curl_file_single.bin',status='unknown',access='direct',recl=4)
+ open(unit=17,file='OUTPUT_FILES/Curl_file_double.bin',status='unknown',access='direct',recl=8)
+
+ end if
+
end if
@@ -174,6 +191,9 @@
buffer_binary(:,1) = sisux(:,irecloc)
if ( number_of_components == 2 ) then
buffer_binary(:,2) = sisuz(:,irecloc)
+ else if ( number_of_components == 3 ) then
+ buffer_binary(:,2) = sisuz(:,irecloc)
+ buffer_binary(:,3) = siscurl(:,irecloc)
end if
#ifdef USE_MPI
@@ -184,6 +204,12 @@
call MPI_RECV(buffer_binary(1,2),NTSTEP_BETWEEN_OUTPUT_SEISMO,MPI_DOUBLE_PRECISION,&
which_proc_receiver(irec),irec,MPI_COMM_WORLD,status,ierror)
end if
+ if ( number_of_components == 3 ) then
+ call MPI_RECV(buffer_binary(1,2),NTSTEP_BETWEEN_OUTPUT_SEISMO,MPI_DOUBLE_PRECISION,&
+ which_proc_receiver(irec),irec,MPI_COMM_WORLD,status,ierror)
+ call MPI_RECV(buffer_binary(1,3),NTSTEP_BETWEEN_OUTPUT_SEISMO,MPI_DOUBLE_PRECISION,&
+ which_proc_receiver(irec),irec,MPI_COMM_WORLD,status,ierror)
+ end if
#endif
@@ -196,6 +222,8 @@
chn = 'BHX'
else if(iorientation == 2) then
chn = 'BHZ'
+ else if(iorientation == 3) then
+ chn = 'cur'
else
call exit_MPI('incorrect channel value')
endif
@@ -251,6 +279,10 @@
write(14,rec=(irec-1)*NSTEP+seismo_offset+isample) sngl(buffer_binary(isample,2))
write(15,rec=(irec-1)*NSTEP+seismo_offset+isample) buffer_binary(isample,2)
end if
+ if ( seismotype == 5 ) then
+ write(16,rec=(irec-1)*NSTEP+seismo_offset+isample) sngl(buffer_binary(isample,3))
+ write(17,rec=(irec-1)*NSTEP+seismo_offset+isample) buffer_binary(isample,3)
+ end if
enddo
#ifdef USE_MPI
@@ -258,9 +290,12 @@
if ( which_proc_receiver(irec) == myrank ) then
irecloc = irecloc + 1
call MPI_SEND(sisux(1,irecloc),NTSTEP_BETWEEN_OUTPUT_SEISMO,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ierror)
- if ( number_of_components == 2 ) then
+ if ( number_of_components >= 2 ) then
call MPI_SEND(sisuz(1,irecloc),NTSTEP_BETWEEN_OUTPUT_SEISMO,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ierror)
end if
+ if ( number_of_components == 3 ) then
+ call MPI_SEND(siscurl(1,irecloc),NTSTEP_BETWEEN_OUTPUT_SEISMO,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ierror)
+ end if
end if
#endif
@@ -275,6 +310,10 @@
close(14)
close(15)
end if
+ if ( seismotype == 5 ) then
+ close(16)
+ close(17)
+ end if
!----
More information about the cig-commits
mailing list