[cig-commits] r19028 - in seismo/2D/SPECFEM2D/trunk: EXAMPLES/noise_uniform src/specfem2D
rmodrak at geodynamics.org
rmodrak at geodynamics.org
Wed Oct 5 23:56:18 PDT 2011
Author: rmodrak
Date: 2011-10-05 23:56:17 -0700 (Wed, 05 Oct 2011)
New Revision: 19028
Modified:
seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/adj_cc.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/check_stability.F90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/noise_tomography.f90
seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
Log:
switch to direct access i/o for writing noise wavefields
Modified: seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/adj_cc.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/adj_cc.f90 2011-10-06 03:31:24 UTC (rev 19027)
+++ seismo/2D/SPECFEM2D/trunk/EXAMPLES/noise_uniform/adj_cc.f90 2011-10-06 06:56:17 UTC (rev 19028)
@@ -9,12 +9,12 @@
logical, parameter :: use_filtering = .true.
!choose exactly one of the following window options
-logical, parameter :: use_negative_branch = .true.
-logical, parameter :: use_positive_branch = .false.
+logical, parameter :: use_negative_branch = .false.
+logical, parameter :: use_positive_branch = .true.
logical, parameter :: use_custom_window = .false.
!choose whether to time reverse, carried out subsequent to all other processing
-logical, parameter :: time_reverse = .true.
+logical, parameter :: time_reverse = .false.
@@ -135,17 +135,25 @@
it_end = ceiling((t_end - t(1))/dt)
if (it_begin < 1) it_begin = 1
if (it_end > nt) it_end = nt
+
elseif (use_positive_branch) then
write(*,*) 'Choosing positive branch'
it_begin = nthalf+1
it_end = nt
+ if (it_begin < nthalf) it_begin = nthalf
+ if (it_end > nt) it_end = nt
+
elseif (use_negative_branch) then
write(*,*) 'Choosing negative branch'
it_begin = 1
it_end = nthalf
+ if (it_begin < 1) it_begin = 1
+ if (it_end > nthalf) it_end = nthalf
+
else
write(*,*) 'Must select one of the following: positive_branch, &
negative_branch, custom_window.'
+
endif
write(*,'(a,2f10.3)') ' Time range: ', t(1), t(nt)
@@ -153,6 +161,7 @@
write(*,'(a,f10.3,f10.3)') ' Filtering: ', 1./freq_high, 1./freq_low
!! Tukey taper
+alpha = w_tukey
k=0
do it = it_begin,it_end
k=k+1
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/check_stability.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/check_stability.F90 2011-10-06 03:31:24 UTC (rev 19027)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/check_stability.F90 2011-10-06 06:56:17 UTC (rev 19028)
@@ -43,7 +43,7 @@
!========================================================================
- subroutine check_stability(myrank,time,it,NSTEP, &
+ subroutine check_stability(myrank,time,it,NSTEP,NOISE_TOMOGRAPHY, &
nglob_acoustic,nglob_elastic,nglob_poroelastic, &
any_elastic_glob,any_elastic,displ_elastic, &
any_poroelastic_glob,any_poroelastic, &
@@ -59,7 +59,7 @@
include "mpif.h"
#endif
- integer :: myrank,it,NSTEP
+ integer :: myrank,it,NSTEP,NOISE_TOMOGRAPHY
double precision :: time
@@ -130,8 +130,13 @@
MPI_MAX, MPI_COMM_WORLD, ier)
#endif
+ if (NOISE_TOMOGRAPHY /= 0) then
+ if (myrank == 0) write(*,*) 'Noise simulation ', NOISE_TOMOGRAPHY, ' of 3'
+ endif
+
+
if (myrank == 0) &
- write(IOUT,*) 'Max norm of vector field in solid (elastic) = ',displnorm_all_glob
+ write(IOUT,*) 'Max norm of vector field in solid (elastic) = ', displnorm_all_glob
! check stability of the code in solid, exit if unstable
! negative values can occur with some compilers when the unstable value is greater
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/noise_tomography.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/noise_tomography.f90 2011-10-06 03:31:24 UTC (rev 19027)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/noise_tomography.f90 2011-10-06 06:56:17 UTC (rev 19028)
@@ -181,17 +181,17 @@
implicit none
include "constants.h"
- !input parameters
+ !input
logical :: p_sv
integer NSTEP, ispec_noise,nglob
integer, dimension(NGLLX,NGLLZ,nglob) :: ibool
real(kind=CUSTOM_REAL) :: deltat, xi_noise, gamma_noise
- !output parameters
+ !output
real(kind=CUSTOM_REAL), dimension(NSTEP) :: time_function_noise
real(kind=CUSTOM_REAL), dimension(3,NGLLX,NGLLZ,NSTEP) :: source_array_noise
- !local parameters
+ !local
integer :: it,i,j,iglob
real(kind=CUSTOM_REAL) :: t
real(kind=CUSTOM_REAL), dimension(NGLLX) :: xigll
@@ -281,7 +281,7 @@
else
- call exit_MPI('unknown noise time function')
+ call exit_MPI('Bad time_function_type in compute_source_array_noise.')
endif
@@ -317,7 +317,7 @@
implicit none
include "constants.h"
- !input parameters
+ !input
logical :: p_sv
integer :: it, NSTEP
integer :: ispec_noise, nglob
@@ -326,7 +326,7 @@
real(kind=CUSTOM_REAL) :: angle_noise
real(kind=CUSTOM_REAL), dimension(3,NGLLX,NGLLZ,NSTEP) :: source_array_noise
- !local variables
+ !local
integer :: i,j,iglob
if(p_sv) then ! P-SV calculation
@@ -354,15 +354,16 @@
! =============================================================================================================
! 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,nspec,nglob,ibool, &
+ subroutine add_surface_movie_noise(p_sv,NOISE_TOMOGRAPHY,it,NSTEP,nspec,nglob,ibool, &
accel_elastic,surface_movie_x_noise,surface_movie_y_noise, &
surface_movie_z_noise,mask_noise,jacobian,wxgll,wzgll)
implicit none
include "constants.h"
- !input parameters
+ !input
logical :: p_sv
+ integer :: NOISE_TOMOGRAPHY
integer :: it, NSTEP
integer :: nspec, nglob
integer :: ibool(NGLLX,NGLLZ,nspec)
@@ -373,22 +374,27 @@
real(kind=CUSTOM_REAL), dimension(3,nglob) :: accel_elastic
real(kind=CUSTOM_REAL) :: wxgll(NGLLX), wzgll(NGLLZ), jacobian(NGLLX,NGLLZ,nspec)
- !local variables
+ !local
integer :: ios, i, j, ispec, iglob
- character(len=60) :: file_in_noise
- write(file_in_noise,"('eta_',i6.6)") NSTEP-it+1
- open(unit=500,file='OUTPUT_FILES/NOISE_TOMOGRAPHY/'//trim(file_in_noise),status='old',form='unformatted',action='read',iostat=ios)
- if( ios /= 0) write(*,*) 'Error retrieving generating wavefield.'
+ if (it==1) then
+ open(unit=500,file='OUTPUT_FILES/NOISE_TOMOGRAPHY/eta',access='direct', &
+ recl=nglob*CUSTOM_REAL,action='read',iostat=ios)
+ if( ios /= 0) call exit_mpi('Error retrieving generating wavefield.')
+ endif
+
if(p_sv) then
- read(500) surface_movie_x_noise
- read(500) surface_movie_z_noise
+ call exit_mpi('P-SV case not yet implemented.')
else
- read(500) surface_movie_y_noise
+ if (NOISE_TOMOGRAPHY==2) read(unit=500,rec=NSTEP-it+1) surface_movie_y_noise
+ if (NOISE_TOMOGRAPHY==3) read(unit=500,rec=it) surface_movie_y_noise
endif
- close(500)
+ if (it==NSTEP) then
+ !close(500)
+ endif
+
do ispec = 1, nspec
do j = 1, NGLLZ
do i = 1, NGLLX
@@ -413,39 +419,51 @@
! =============================================================================================================
! save a snapshot of the "generating wavefield" eta that will be used to drive
! the "ensemble forward wavefield"
- subroutine save_surface_movie_noise(NOISE_TOMOGRAPHY,p_sv,it,nglob,displ_elastic)
+ subroutine save_surface_movie_noise(NOISE_TOMOGRAPHY,p_sv,it,NSTEP,nglob,displ_elastic)
implicit none
include "constants.h"
- !input paramters
+ !input
integer :: NOISE_TOMOGRAPHY
logical :: p_sv
- integer :: it, nglob
+ integer :: it, NSTEP, nglob
real(kind=CUSTOM_REAL), dimension(3,nglob) :: displ_elastic
- !local parameters
- character(len=60) file_out_noise
-
- if (NOISE_TOMOGRAPHY == 1) then
- write(file_out_noise,"('OUTPUT_FILES/NOISE_TOMOGRAPHY/eta_',i6.6)") it
+ !local
+ integer :: ios
- elseif (NOISE_TOMOGRAPHY == 2) then
- write(file_out_noise,"('OUTPUT_FILES/NOISE_TOMOGRAPHY/phi_',i6.6)") it
+ if (it==1) then
- else
- call exit_mpi('Bad value of NOISE_TOMOGRAPHY in save_surface_movie_noise.')
+ if (NOISE_TOMOGRAPHY == 1) then
- endif
-
- open(unit=500,file=trim(file_out_noise),status='unknown',form='unformatted',action='write')
+ open(unit=500,file='OUTPUT_FILES/NOISE_TOMOGRAPHY/eta',access='direct', &
+ recl=nglob*CUSTOM_REAL,action='write',iostat=ios)
+ if( ios /= 0) call exit_mpi('Error saving generating wavefield.')
+
+ elseif (NOISE_TOMOGRAPHY == 2) then
+
+ open(unit=500,file='OUTPUT_FILES/NOISE_TOMOGRAPHY/phi',access='direct', &
+ recl=nglob*CUSTOM_REAL,action='write',iostat=ios)
+ if( ios /= 0) call exit_mpi('Error saving ensemble forward wavefield.')
+
+ else
+ call exit_mpi('Bad value of NOISE_TOMOGRAPHY in save_surface_movie_noise.')
+
+ endif
+
+ endif ! (it==1)
+
if(p_sv) then
- write(500) displ_elastic(1,:)
- write(500) displ_elastic(3,:)
+ call exit_mpi('P-SV case not yet implemented.')
else
- write(500) displ_elastic(2,:)
+ write(unit=500,rec=it) displ_elastic(2,:)
endif
+ if (it==NSTEP) then
+ !close(500)
+ endif
+
end subroutine save_surface_movie_noise
! =============================================================================================================
@@ -482,29 +500,29 @@
! =============================================================================================================
! auxillary routine
- subroutine elem_to_glob(nspec,nglob,ibool,array_elem,array_glob)
+ subroutine spec2glob(nspec,nglob,ibool,array_spec,array_glob)
implicit none
include "constants.h"
- !input paramters
+ !input
integer :: nspec, nglob
integer :: ibool(NGLLX,NGLLZ,nspec)
real(kind=CUSTOM_REAL), dimension(nglob) :: array_glob
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: array_elem
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: array_spec
- !local parameters
+ !local
integer :: i,j,iglob,ispec
do ispec = 1, nspec
do j = 1, NGLLZ
do i = 1, NGLLX
iglob = ibool(i,j,ispec)
- array_glob(iglob) = array_elem(i,j,ispec)
+ array_glob(iglob) = array_spec(i,j,ispec)
enddo
enddo
enddo
- end subroutine elem_to_glob
+ end subroutine spec2glob
Modified: seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2011-10-06 03:31:24 UTC (rev 19027)
+++ seismo/2D/SPECFEM2D/trunk/src/specfem2D/specfem2D.F90 2011-10-06 06:56:17 UTC (rev 19028)
@@ -839,9 +839,10 @@
! For writing noise wavefields
logical :: output_wavefields_noise = .true.
+ logical :: ex, od
integer :: noise_output_ncol
real(kind=CUSTOM_REAL), dimension(:,:), allocatable :: noise_output_array
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: noise_output_dim_5
+ real(kind=CUSTOM_REAL), dimension(:), allocatable :: noise_output_rhokl
character(len=512) :: noise_output_file
! For noise tomography only - specify whether to reconstruct ensemble forward
@@ -3695,7 +3696,7 @@
!prepare array that will hold wavefield snapshots
noise_output_ncol = 5
allocate(noise_output_array(noise_output_ncol,nglob))
- allocate(noise_output_dim_5(nglob))
+ allocate(noise_output_rhokl(nglob))
endif
endif
@@ -5061,13 +5062,13 @@
accel_elastic,angle_noise,source_array_noise)
elseif (NOISE_TOMOGRAPHY == 2) then
- call add_surface_movie_noise(p_sv,it,NSTEP,nspec,nglob,ibool,accel_elastic, &
+ call add_surface_movie_noise(p_sv,NOISE_TOMOGRAPHY,it,NSTEP,nspec,nglob,ibool,accel_elastic, &
surface_movie_x_noise,surface_movie_y_noise, &
surface_movie_z_noise,mask_noise,jacobian,wxgll,wzgll)
elseif (NOISE_TOMOGRAPHY == 3) then
if (.not. save_everywhere) then
- call add_surface_movie_noise(p_sv,NSTEP-it+1,NSTEP,nspec,nglob,ibool,b_accel_elastic, &
+ call add_surface_movie_noise(p_sv,NOISE_TOMOGRAPHY,it,NSTEP,nspec,nglob,ibool,b_accel_elastic, &
surface_movie_x_noise,surface_movie_y_noise, &
surface_movie_z_noise,mask_noise,jacobian,wxgll,wzgll)
endif
@@ -5938,23 +5939,21 @@
!<NOISE_TOMOGRAPHY
if ( NOISE_TOMOGRAPHY == 1 ) then
- call save_surface_movie_noise(NOISE_TOMOGRAPHY,p_sv,it,nglob,displ_elastic)
+ call save_surface_movie_noise(NOISE_TOMOGRAPHY,p_sv,it,NSTEP,nglob,displ_elastic)
elseif ( NOISE_TOMOGRAPHY == 2 .and. save_everywhere ) then
- call save_surface_movie_noise(NOISE_TOMOGRAPHY,p_sv,it,nglob,displ_elastic)
+ call save_surface_movie_noise(NOISE_TOMOGRAPHY,p_sv,it,NSTEP,nglob,displ_elastic)
elseif ( NOISE_TOMOGRAPHY == 3 .and. save_everywhere ) then
- write(noise_output_file,"('phi_',i6.6)") NSTEP-it+1
- open(unit=500,file='OUTPUT_FILES/NOISE_TOMOGRAPHY/'//trim(noise_output_file), &
- status='old',form='unformatted',action='read',iostat=ios)
+ if (it==1) &
+ open(unit=500,file='OUTPUT_FILES/NOISE_TOMOGRAPHY/phi',access='direct', &
+ recl=nglob*CUSTOM_REAL,action='write',iostat=ios)
if( ios /= 0) write(*,*) 'Error retrieving ensemble forward wavefield.'
if(p_sv) then
- read(500) b_displ_elastic(1,:)
- read(500) b_displ_elastic(3,:)
+ call exit_mpi('P-SV case not yet implemented.')
else
- read(500) b_displ_elastic(2,:)
+ read(unit=500,rec=NSTEP-it+1) b_displ_elastic(2,:)
endif
- close(500)
endif
@@ -6011,7 +6010,7 @@
!---- display time step and max of norm of displacement
if(mod(it,NTSTEP_BETWEEN_OUTPUT_INFO) == 0 .or. it == 5 .or. it == NSTEP) then
- call check_stability(myrank,time,it,NSTEP, &
+ call check_stability(myrank,time,it,NSTEP,NOISE_TOMOGRAPHY, &
nglob_acoustic,nglob_elastic,nglob_poroelastic, &
any_elastic_glob,any_elastic,displ_elastic, &
any_poroelastic_glob,any_poroelastic, &
@@ -6572,22 +6571,21 @@
if ( NOISE_TOMOGRAPHY == 3 .and. output_wavefields_noise ) then
!load ensemble foward source
- write(noise_output_file,"('eta_',i6.6)") it
- open(unit=500,file='OUTPUT_FILES/NOISE_TOMOGRAPHY/'//trim(noise_output_file), &
- status='old',form='unformatted',action='read',iostat=ios)
- if( ios /= 0) write(*,*) 'Error preparing noise output.'
- read(500) surface_movie_y_noise
- close(500)
+ inquire(unit=500,exist=ex,opened=od)
+ if (.not. od) &
+ open(unit=500,file='OUTPUT_FILES/NOISE_TOMOGRAPHY/eta',access='direct', &
+ recl=nglob*CUSTOM_REAL,action='write',iostat=ios)
+ read(unit=500,rec=it) surface_movie_y_noise
!load product of fwd, adj wavefields
- call elem_to_glob(nspec,nglob,ibool,rho_kl,noise_output_dim_5)
+ call spec2glob(nspec,nglob,ibool,rho_kl,noise_output_rhokl)
!write text file
noise_output_array(1,:) = surface_movie_y_noise(:) * mask_noise(:)
noise_output_array(2,:) = b_displ_elastic(2,:)
noise_output_array(3,:) = accel_elastic(2,:)
noise_output_array(4,:) = rho_k(:)
- noise_output_array(5,:) = noise_output_dim_5(:)
+ noise_output_array(5,:) = noise_output_rhokl(:)
write(noise_output_file,"('OUTPUT_FILES/snapshot_all_',i6.6)") it
call snapshots_noise(noise_output_ncol,nglob,noise_output_file,noise_output_array)
More information about the CIG-COMMITS
mailing list