[cig-commits] r18969 - in seismo/2D/SPECFEM2D/trunk: EXAMPLES/noise_uniform src/specfem2D
rmodrak at geodynamics.org
rmodrak at geodynamics.org
Fri Sep 23 13:24:20 PDT 2011
Author: rmodrak
Date: 2011-09-23 13:24:19 -0700 (Fri, 23 Sep 2011)
New Revision: 18969
Added:
seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/ellipse_large
Modified:
seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_1
seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_2
seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_3
seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/adj_cc.f90
seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/clean
seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/contribution_1
seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/contribution_2
seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/process.sh
seismo/2D/SPECFEM2D/trunk/src/specfem2D/noise_tomography.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
Log:
for noise tomography, fixes weighting on ensemble forward wavefield source
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_1
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_1 2011-09-23 20:00:18 UTC (rev 18968)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_1 2011-09-23 20:24:19 UTC (rev 18969)
@@ -22,7 +22,7 @@
p_sv = .false. # set the type of calculation (P-SV or SH/membrane waves)
# time step parameters
-nt = 4000 # total number of time steps
+nt = 2000 # total number of time steps
deltat = 4.d-2 # duration of a time step
# source parameters
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_2
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_2 2011-09-23 20:00:18 UTC (rev 18968)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_2 2011-09-23 20:24:19 UTC (rev 18969)
@@ -22,7 +22,7 @@
p_sv = .false. # set the type of calculation (P-SV or SH/membrane waves)
# time step parameters
-nt = 4000 # total number of time steps
+nt = 2000 # total number of time steps
deltat = 4.d-2 # duration of a time step
# source parameters
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_3
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_3 2011-09-23 20:00:18 UTC (rev 18968)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/Par_file_noise_3 2011-09-23 20:24:19 UTC (rev 18969)
@@ -22,7 +22,7 @@
p_sv = .false. # set the type of calculation (P-SV or SH/membrane waves)
# time step parameters
-nt = 4000 # total number of time steps
+nt = 2000 # total number of time steps
deltat = 4.d-2 # duration of a time step
# source parameters
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/adj_cc.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/adj_cc.f90 2011-09-23 20:00:18 UTC (rev 18968)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/adj_cc.f90 2011-09-23 20:24:19 UTC (rev 18969)
@@ -13,8 +13,8 @@
logical, parameter :: use_positive_branch = .false.
logical, parameter :: use_custom_window = .false.
-!choose whether to time reverse (carried out subsequent to all other processing)
-logical, parameter :: reverse = .true.
+!choose whether to time reverse, carried out subsequent to all other processing
+logical, parameter :: time_reverse = .true.
@@ -183,8 +183,8 @@
write(*,*) 'Writing to file: '//trim(file_in)//'.adj'
do it = 1,nt
- if (.not. reverse) write(1002,'(f16.12,1pe18.8)'), t(it), seismo_adj(it)
- if (reverse) write(1002,'(f16.12,1pe18.8)'), t(it), seismo_adj(nt+1-it)
+ if (.not. time_reverse) write(1002,'(f16.12,1pe18.8)'), t(it), seismo_adj(it)
+ if (time_reverse) write(1002,'(f16.12,1pe18.8)'), t(it), seismo_adj(nt+1-it)
end do
close(1002)
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/clean
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/clean 2011-09-23 20:00:18 UTC (rev 18968)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/clean 2011-09-23 20:24:19 UTC (rev 18969)
@@ -5,7 +5,7 @@
echo "Cleaning up."
-rm -rf DATA SEM NOISE_TOMOGRAPHY OUTPUT_FILES OUTPUT_ALL SNAPSHOTS &> /dev/null
+rm -rf DATA SEM OUTPUT_FILES OUTPUT_ALL SNAPSHOTS &> /dev/null
rm xmeshfem2D xspecfem2D &> /dev/null
rm adj_run &> /dev/null
rm log.* &> /dev/null
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/contribution_1
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/contribution_1 2011-09-23 20:00:18 UTC (rev 18968)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/contribution_1 2011-09-23 20:24:19 UTC (rev 18969)
@@ -5,8 +5,8 @@
#prepare directories
-rm -rf DATA SEM NOISE_TOMOGRAPHY OUTPUT_FILES OUTPUT_ALL SNAPSHOTS
-mkdir -p DATA SEM NOISE_TOMOGRAPHY OUTPUT_FILES OUTPUT_ALL SNAPSHOTS DATA/NOISE_TOMOGRAPHY OUTPUT_FILES/NOISE_TOMOGRAPHY
+rm -rf DATA SEM OUTPUT_FILES OUTPUT_ALL
+mkdir -p DATA SEM OUTPUT_FILES OUTPUT_ALL DATA/NOISE_TOMOGRAPHY OUTPUT_FILES/NOISE_TOMOGRAPHY
#prepare files
@@ -15,6 +15,7 @@
cp uniform.dat DATA/
echo 1 > DATA/NOISE_TOMOGRAPHY/irec_master
if [ -f S_squared ]; then cp S_squared DATA/NOISE_TOMOGRAPHY/; fi
+if [ -f model_velocity.dat_input ]; then cp model_velocity.dat_input DATA/; fi
#compile
@@ -31,25 +32,30 @@
cp Par_file_noise_1 DATA/Par_file
./xmeshfem2D; ./xspecfem2D
mkdir OUTPUT_ALL/step_1
-mv OUTPUT_FILES/image* OUTPUT_ALL/step_1
+#mv OUTPUT_FILES/image* OUTPUT_ALL/step_1
mv OUTPUT_FILES/*.semd OUTPUT_ALL/step_1
mv DATA/Par_file OUTPUT_ALL/step_1
-#mv OUTPUT_FILES/snapshot* SNAPSHOTS/
+mv OUTPUT_FILES/mask_noise OUTPUT_ALL/
+#mv OUTPUT_FILES/snapshot* OUTPUT_ALL/
#simulation 2
cp Par_file_noise_2 DATA/Par_file
./xmeshfem2D; ./xspecfem2D
mkdir OUTPUT_ALL/step_2
-mv OUTPUT_FILES/image* OUTPUT_ALL/step_2
+#mv OUTPUT_FILES/image* OUTPUT_ALL/step_2
mv OUTPUT_FILES/*.semd OUTPUT_ALL/step_2
mv DATA/Par_file OUTPUT_ALL/step_2
-#mv OUTPUT_FILES/snapshot* SNAPSHOTS/
+#mv OUTPUT_FILES/snapshot* OUTPUT_ALL/
#prepare adjoint source
ADJ_CODE=adj_cc.f90
-sed -i'.bak' 's/reverse = .[a-z]*./reverse = .false./' $ADJ_CODE
+sed -i'.bak' 's/use_positive_branch = .[a-z]*./use_positive_branch = .true./' $ADJ_CODE
+sed -i'.bak' 's/use_negative_branch = .[a-z]*./use_negative_branch = .false./' $ADJ_CODE
+sed -i'.bak' 's/use_custom_window = .[a-z]*./use_custom_window = .false./' $ADJ_CODE
+sed -i'.bak' 's/time_reverse = .[a-z]*./time_reverse = .true./' $ADJ_CODE
+
gfortran $ADJ_CODE -o adj_run
cp OUTPUT_ALL/step_2/*.semd SEM/
./adj_run SEM/S0003.AA.BXY.semd
@@ -67,9 +73,9 @@
cp Par_file_noise_3 DATA/Par_file
./xmeshfem2D; ./xspecfem2D
mkdir OUTPUT_ALL/step_3
-mv OUTPUT_FILES/image* OUTPUT_ALL/step_3
+#mv OUTPUT_FILES/image* OUTPUT_ALL/step_3
mv SEM/*Y.adj OUTPUT_ALL/step_3
mv OUTPUT_FILES/proc* OUTPUT_ALL/step_3
mv DATA/Par_file OUTPUT_ALL/step_3
-mv OUTPUT_FILES/snapshot* SNAPSHOTS/
+mv OUTPUT_FILES/snapshot* OUTPUT_ALL/
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/contribution_2
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/contribution_2 2011-09-23 20:00:18 UTC (rev 18968)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/contribution_2 2011-09-23 20:24:19 UTC (rev 18969)
@@ -5,8 +5,8 @@
#prepare directories
-rm -rf DATA SEM NOISE_TOMOGRAPHY OUTPUT_FILES OUTPUT_ALL SNAPSHOTS
-mkdir -p DATA SEM NOISE_TOMOGRAPHY OUTPUT_FILES OUTPUT_ALL SNAPSHOTS DATA/NOISE_TOMOGRAPHY OUTPUT_FILES/NOISE_TOMOGRAPHY
+rm -rf DATA SEM OUTPUT_FILES OUTPUT_ALL
+mkdir -p DATA SEM OUTPUT_FILES OUTPUT_ALL DATA/NOISE_TOMOGRAPHY OUTPUT_FILES/NOISE_TOMOGRAPHY
#prepare files
@@ -15,6 +15,7 @@
cp uniform.dat DATA/
echo 3 > DATA/NOISE_TOMOGRAPHY/irec_master
if [ -f S_squared ]; then cp S_squared DATA/NOISE_TOMOGRAPHY/; fi
+if [ -f model_velocity.dat_input ]; then cp model_velocity.dat_input DATA/; fi
#compile
@@ -31,25 +32,30 @@
cp Par_file_noise_1 DATA/Par_file
./xmeshfem2D; ./xspecfem2D
mkdir OUTPUT_ALL/step_1
-mv OUTPUT_FILES/image* OUTPUT_ALL/step_1
+#mv OUTPUT_FILES/image* OUTPUT_ALL/step_1
mv OUTPUT_FILES/*.semd OUTPUT_ALL/step_1
mv DATA/Par_file OUTPUT_ALL/step_1
-#mv OUTPUT_FILES/snapshot* SNAPSHOTS/
+mv OUTPUT_FILES/mask_noise OUTPUT_ALL/
+#mv OUTPUT_FILES/snapshot* OUTPUT_ALL/
#simulation 2
cp Par_file_noise_2 DATA/Par_file
./xmeshfem2D; ./xspecfem2D
mkdir OUTPUT_ALL/step_2
-mv OUTPUT_FILES/image* OUTPUT_ALL/step_2
+#mv OUTPUT_FILES/image* OUTPUT_ALL/step_2
mv OUTPUT_FILES/*.semd OUTPUT_ALL/step_2
mv DATA/Par_file OUTPUT_ALL/step_2
-#mv OUTPUT_FILES/snapshot* SNAPSHOTS/
+#mv OUTPUT_FILES/snapshot* OUTPUT_ALL/
#prepare adjoint source
ADJ_CODE=adj_cc.f90
-sed -i'.bak' 's/reverse = .[a-z]*./reverse = .true./' $ADJ_CODE
+sed -i'.bak' 's/use_positive_branch = .[a-z]*./use_positive_branch = .false./' $ADJ_CODE
+sed -i'.bak' 's/use_negative_branch = .[a-z]*./use_negative_branch = .true./' $ADJ_CODE
+sed -i'.bak' 's/use_custom_window = .[a-z]*./use_custom_window = .false./' $ADJ_CODE
+sed -i'.bak' 's/time_reverse = .[a-z]*./time_reverse = .true./' $ADJ_CODE
+
gfortran $ADJ_CODE -o adj_run
cp OUTPUT_ALL/step_2/*.semd SEM/
./adj_run SEM/S0001.AA.BXY.semd
@@ -67,9 +73,9 @@
cp Par_file_noise_3 DATA/Par_file
./xmeshfem2D; ./xspecfem2D
mkdir OUTPUT_ALL/step_3
-mv OUTPUT_FILES/image* OUTPUT_ALL/step_3
+#mv OUTPUT_FILES/image* OUTPUT_ALL/step_3
mv SEM/*Y.adj OUTPUT_ALL/step_3
mv OUTPUT_FILES/proc* OUTPUT_ALL/step_3
mv DATA/Par_file OUTPUT_ALL/step_3
-mv OUTPUT_FILES/snapshot* SNAPSHOTS/
+mv OUTPUT_FILES/snapshot* OUTPUT_ALL/
Added: seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/ellipse_large
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/ellipse_large (rev 0)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/ellipse_large 2011-09-23 20:24:19 UTC (rev 18969)
@@ -0,0 +1,42 @@
+! =============================================================================================================
+! specify spatial distribution of microseismic noise sources
+! USERS need to modify this subroutine to suit their own needs
+ subroutine create_mask_noise(nglob,coord,mask_noise)
+
+ implicit none
+ include "constants.h"
+
+ !input
+ integer :: nglob
+ real(kind=CUSTOM_REAL), dimension(2,nglob) :: coord
+
+ !output
+ real(kind=CUSTOM_REAL), dimension(nglob) :: mask_noise
+
+ !local
+ integer :: iglob
+ real(kind=CUSTOM_REAL) :: xx,zz,aa,phi
+
+ aa = PI/8.
+
+ !specify distribution of noise sources as a function of xx, zz
+ do iglob = 1,nglob
+ xx = coord(1,iglob)
+ zz = coord(2,iglob)
+
+ phi = (((xx-105.e3)*cos(aa)+(zz-92.e3)*sin(aa))/28.e3)**2. + &
+ (((xx-105.e3)*sin(aa)-(zz-92.e3)*cos(aa))/20.e3)**2.
+
+ if (phi < 0.9) then
+ mask_noise(iglob) = 1.
+ elseif (phi < 1.1) then
+ mask_noise(iglob) = 0.5*(1.+cos(5.*PI*(phi/1.1-0.8)))
+ else
+ mask_noise(iglob) = 0.
+ endif
+
+
+ enddo
+
+ end subroutine create_mask_noise
+
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/process.sh
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/process.sh 2011-09-23 20:00:18 UTC (rev 18968)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/process.sh 2011-09-23 20:24:19 UTC (rev 18969)
@@ -2,9 +2,7 @@
./contribution_1
mv OUTPUT_ALL OUTPUT_ALL_1
-mv SNAPSHOT SNAPSHOT_1
./contribution_2
mv OUTPUT_ALL OUTPUT_ALL_2
-mv SNAPSHOT SNAPSHOT_2
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/noise_tomography.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/noise_tomography.f90 2011-09-23 20:00:18 UTC (rev 18968)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/noise_tomography.f90 2011-09-23 20:24:19 UTC (rev 18969)
@@ -43,10 +43,9 @@
!NOISE TOMOGRAPHY TO DO LIST
-!1. Replace ad hoc scaling in add_surface_movie_noise by expression involving Jacobian
-!2. Use separate STATIONS_ADJOINT file
-!3. Add exploration test case under EXAMPLES
-!4. Update manual
+!1. Use separate STATIONS_ADJOINT file
+!2. Add exploration test case under EXAMPLES
+!3. Update manual
! =============================================================================================================
! specify spatial distribution of microseismic noise sources
@@ -355,9 +354,9 @@
! =============================================================================================================
! read in and inject the "source" that drives the "enemble forward wavefield"
! (recall that the ensemble forward wavefield has a spatially distributed source)
- subroutine add_surface_movie_noise(p_sv,it,NSTEP,nglob, &
+ subroutine add_surface_movie_noise(p_sv,it,NSTEP,nspec,nglob,ibool, &
accel_elastic,surface_movie_x_noise,surface_movie_y_noise, &
- surface_movie_z_noise,mask_noise)
+ surface_movie_z_noise,mask_noise,jacobian,wxgll,wzgll)
implicit none
include "constants.h"
@@ -365,16 +364,17 @@
!input parameters
logical :: p_sv
integer :: it, NSTEP
- integer :: nglob
+ integer :: nspec, nglob
+ integer :: ibool(NGLLX,NGLLZ,nspec)
real(kind=CUSTOM_REAL), dimension(nglob) :: mask_noise
real(kind=CUSTOM_REAL), dimension(nglob) :: &
surface_movie_x_noise, surface_movie_y_noise, surface_movie_z_noise
real(kind=CUSTOM_REAL), dimension(3,nglob) :: accel_elastic
+ real(kind=CUSTOM_REAL) :: wxgll(NGLLX), wzgll(NGLLZ), jacobian(NGLLX,NGLLZ,nspec)
!local variables
- integer :: ios
- real(kind=CUSTOM_REAL) :: factor_noise
+ integer :: ios, i, j, ispec, iglob
character(len=60) :: file_in_noise
@@ -389,30 +389,25 @@
endif
close(500)
- factor_noise = 1.d6
+ do ispec = 1, nspec
+ do j = 1, NGLLZ
+ do i = 1, NGLLX
+ if (p_sv) then ! P-SV calculation
+ iglob = ibool(i,j,ispec)
+ accel_elastic(1,iglob) = accel_elastic(1,iglob) + surface_movie_x_noise(iglob) * &
+ mask_noise(iglob) * wxgll(i)*wzgll(j)*jacobian(i,j,ispec)
+ accel_elastic(3,iglob) = accel_elastic(3,iglob) + surface_movie_z_noise(iglob) * &
+ mask_noise(iglob) * wxgll(i)*wzgll(j)*jacobian(i,j,ispec)
- if(p_sv) then ! P-SV calculation
- accel_elastic(1,:) = accel_elastic(1,:) + factor_noise * mask_noise(:) * surface_movie_x_noise(:)
- accel_elastic(3,:) = accel_elastic(3,:) + factor_noise * mask_noise(:) * surface_movie_z_noise(:)
-! note: above implementation does not yet include non-vertical noise
-! excitation
+ else ! SH (membrane) calculation
+ iglob = ibool(i,j,ispec)
+ accel_elastic(2,iglob) = accel_elastic(2,iglob) + surface_movie_y_noise(iglob) * &
+ mask_noise(iglob) * wxgll(i)*wzgll(j)*jacobian(i,j,ispec)
+ endif
+ enddo
+ enddo
+ enddo
- else ! SH (membrane) calculation
- accel_elastic(2,:) = accel_elastic(2,:) + factor_noise * mask_noise(:) * surface_movie_y_noise(:)
-
-! note: above implementation includes ad hoc scaling by a number
-! "factor_noise", meant to approximate the correct jacobian and gll weighting
-! correct implementation should look something like:
-! do iglob = 1,nglob
-! hlagrange = weight(iglob) * jacobian(iglob)
-! accel_elastic(2,iglob) = accel_elastic(2,iglob) &
-! + hlagrange * surface_movie_y_noise(iglob)
-! enddo
-! enddo
-! enddo
-
- endif
-
end subroutine add_surface_movie_noise
! =============================================================================================================
@@ -458,20 +453,15 @@
!local parameters
integer :: i,iglob
- real(kind=CUSTOM_REAL), dimension(ncol) :: row
open(unit=504,file=filename,status='unknown',action='write')
do iglob = 1,nglob
- row(:) = array_all(:,iglob)
-
do i = 1,ncol-1
- if (abs(row(i)) < 1.e-32) row(i) = 0.
- write(unit=504,fmt='(1pe15.5e3)',advance='no') row(i)
+ write(unit=504,fmt='(1pe12.3e3)',advance='no') array_all(i,iglob)
enddo
- if (abs(row(ncol)) < 1.e-32) row(ncol) = 0.
- write(unit=504,fmt='(1pe15.5e3)') row(ncol)
+ write(unit=504,fmt='(1pe13.3e3)') array_all(ncol,iglob)
enddo
@@ -509,36 +499,3 @@
end subroutine elem_to_glob
-
-! =============================================================================================================
-! auxillary routine
- subroutine snapshot_glob(nglob,coord,filename,array1)
-
- implicit none
- include "constants.h"
-
- !input paramters
- integer :: nglob
- character(len=100) filename
-
- real(kind=CUSTOM_REAL), dimension(2,nglob) :: coord
- real(kind=CUSTOM_REAL), dimension(nglob) :: array1
-
- !local parameters
- integer :: iglob
- real(kind=CUSTOM_REAL) :: xx,zz
-
-
- open(unit=504,file=filename,status='unknown',action='write')
-
- do iglob = 1, nglob
- xx = coord(1,iglob)
- zz = coord(2,iglob)
- write(504,*) xx, zz, array1(iglob)
- enddo
-
- close(504)
-
- end subroutine snapshot_glob
-
-
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2011-09-23 20:00:18 UTC (rev 18968)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2011-09-23 20:24:19 UTC (rev 18969)
@@ -828,7 +828,7 @@
integer :: ncol
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rho_klglob
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: noise_all
- character(len=100) :: fnoise
+ character(len=100) :: noise_file
@@ -3591,11 +3591,6 @@
allocate(surface_movie_y_noise(nglob))
allocate(surface_movie_z_noise(nglob))
- ncol = 6
- allocate(rho_klglob(nglob))
- allocate(noise_all(ncol,nglob))
-
-
!read in parameters for noise tomography
call read_parameters_noise(NOISE_TOMOGRAPHY,SIMULATION_TYPE,SAVE_FORWARD, &
any_acoustic,any_poroelastic,p_sv,irec_master, &
@@ -3610,13 +3605,70 @@
call compute_source_array_noise(p_sv,NSTEP,deltat,nglob,ibool,ispec_noise, &
xi_noise,gamma_noise,xigll,zigll, &
time_function_noise,source_array_noise)
+ !write out local coordinates of mesh
+ open(unit=504,file='OUTPUT_FILES/elem',status='unknown',action='write')
+ do ispec = 1, nspec
+ do j = 1, NGLLZ
+ do i = 1, NGLLX
+ iglob = ibool(i,j,ispec)
+ write(504,'(1pe11.3,1pe11.3,2i3,i7)') &
+ coord(1,iglob), coord(2,iglob), i, j, ispec
+ enddo
+ enddo
+ enddo
+ close(504)
+ !write out spatial distribution of noise sources
+ call create_mask_noise(nglob,coord,mask_noise)
+ open(unit=504,file='OUTPUT_FILES/mask_noise',status='unknown',action='write')
+ do iglob = 1, nglob
+ write(504,'(1pe11.3,1pe11.3,1pe11.3)') &
+ coord(1,iglob), coord(2,iglob), mask_noise(iglob)
+ enddo
+ close(504)
+
+ !write out velocity model
+ if(assign_external_model) then
+ open(unit=504,file='OUTPUT_FILES/model_rho_vp_vs',status='unknown',action='write')
+ do ispec = 1, nspec
+ do j = 1, NGLLZ
+ do i = 1, NGLLX
+ iglob = ibool(i,j,ispec)
+ write(504,'(1pe11.3,1pe11.3,1pe11.3,1pe11.3,1pe11.3)') &
+ coord(1,iglob), coord(2,iglob), &
+ rhoext(i,j,ispec), vpext(i,j,ispec), vsext(i,j,ispec)
+ enddo
+ enddo
+ enddo
+ close(504)
+ else
+ open(unit=504,file='OUTPUT_FILES/model_rho_kappa_mu',status='unknown',action='write')
+ do ispec = 1, nspec
+ do j = 1, NGLLZ
+ do i = 1, NGLLX
+ iglob = ibool(i,j,ispec)
+ write(504,'(1pe11.3,1pe11.3,1pe11.3,1pe11.3,1pe11.3)') &
+ coord(1,iglob), coord(2,iglob), density(1,kmato(ispec)), &
+ poroelastcoef(1,1,kmato(ispec)) + 2.d0/3.d0*poroelastcoef(2,1,kmato(ispec)), &
+ poroelastcoef(2,1,kmato(ispec))
+
+ enddo
+ enddo
+ enddo
+ close(504)
+ endif
+
elseif (NOISE_TOMOGRAPHY == 2) then
call create_mask_noise(nglob,coord,mask_noise)
elseif (NOISE_TOMOGRAPHY == 3) then
call create_mask_noise(nglob,coord,mask_noise)
+ !prepare array that will hold wavefield snapshots
+ ncol = 5
+ allocate(noise_all(ncol,nglob))
+ allocate(rho_klglob(nglob))
+
endif
!>NOISE_TOMOGRAPHY
@@ -4980,14 +5032,14 @@
accel_elastic,angle_noise,source_array_noise)
elseif (NOISE_TOMOGRAPHY == 2) then
- call add_surface_movie_noise(p_sv,it,NSTEP,nglob,accel_elastic, &
+ call add_surface_movie_noise(p_sv,it,NSTEP,nspec,nglob,ibool,accel_elastic, &
surface_movie_x_noise,surface_movie_y_noise, &
- surface_movie_z_noise,mask_noise)
+ surface_movie_z_noise,mask_noise,jacobian,wxgll,wzgll)
elseif (NOISE_TOMOGRAPHY == 3) then
- call add_surface_movie_noise(p_sv,it,NSTEP,nglob,b_accel_elastic, &
+ call add_surface_movie_noise(p_sv,it,NSTEP,nspec,nglob,ibool,b_accel_elastic, &
surface_movie_x_noise,surface_movie_y_noise, &
- surface_movie_z_noise,mask_noise)
+ surface_movie_z_noise,mask_noise,jacobian,wxgll,wzgll)
endif
@@ -6467,25 +6519,20 @@
!<NOISE_TOMOGRAPHY
if (NOISE_TOMOGRAPHY == 1) then
- !write(fnoise,"('OUTPUT_FILES/snapshot_eta_',i6.6)") it
- !call snapshot_glob(nglob,coord,fnoise,accel_elastic(2,:))
elseif (NOISE_TOMOGRAPHY == 2) then
- !write(fnoise,"('OUTPUT_FILES/snapshot_phi_',i6.6)") it
- !call snapshot_glob(nglob,coord,fnoise,displ_elastic(2,:))
elseif (NOISE_TOMOGRAPHY == 3) then
call elem_to_glob(nspec,nglob,ibool,rho_kl,rho_klglob)
- noise_all(1,:) = coord(1,:)
- noise_all(2,:) = coord(2,:)
- noise_all(3,:) = b_displ_elastic(2,:)
- noise_all(4,:) = accel_elastic(2,:)
- noise_all(5,:) = rho_k
- noise_all(6,:) = rho_klglob
+ noise_all(1,:) = surface_movie_y_noise
+ noise_all(2,:) = b_displ_elastic(2,:)
+ noise_all(3,:) = accel_elastic(2,:)
+ noise_all(4,:) = rho_k
+ noise_all(5,:) = rho_klglob
- write(fnoise,"('OUTPUT_FILES/snapshot_all_',i6.6)") it
- call snapshot_all(ncol,nglob,fnoise,noise_all)
+ write(noise_file,"('OUTPUT_FILES/snapshot_all_',i6.6)") it
+ call snapshot_all(ncol,nglob,noise_file,noise_all)
endif
More information about the CIG-COMMITS
mailing list