[cig-commits] r17885 - in seismo/3D/ADJOINT_TOMO/iterate_adj/pangu: model_vp_vs/src smooth/src
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Wed Feb 16 13:49:17 PST 2011
Author: danielpeter
Date: 2011-02-16 13:49:17 -0800 (Wed, 16 Feb 2011)
New Revision: 17885
Removed:
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/mpif.h
Modified:
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/add_model_globe_iso.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/add_model_globe_tiso.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/add_model_globe_tiso_iso.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/src/smooth_sem_globe.f90
Log:
updates smoothing and model update routines for SPECFEM3D_GLOBE models in iterate_adj/pangu/
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/add_model_globe_iso.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/add_model_globe_iso.f90 2011-02-16 15:53:48 UTC (rev 17884)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/add_model_globe_iso.f90 2011-02-16 21:49:17 UTC (rev 17885)
@@ -1,37 +1,37 @@
! add_model_globe_iso
!
-! this program can be used to update ISOTROPIC model files with
+! this program can be used to update ISOTROPIC model files with
! (smoothed & summed) event kernels.
! the kernels are given for isotropic parameters (alpha,beta,rho) or ( bulk_c,beta,rho).
!
-! the algorithm uses a steepest descent method with a step length
+! the algorithm uses a steepest descent method with a step length
! determined by the given maximum update percentage.
-!
+!
! input:
-! - step_fac : step length to update the models, f.e. 0.03 for plusminus 3%
+! - step_fac : step length to update the models, f.e. 0.03 for plusminus 3%
!
! setup:
!
!- INPUT_MODEL/ contains:
-! proc000***_reg1_vs.bin &
-! proc000***_reg1_vp.bin &
+! proc000***_reg1_vs.bin &
+! proc000***_reg1_vp.bin &
! proc000***_reg1_rho.bin
!
-!- INPUT_GRADIENT/ contains:
-! proc000***_reg1_bulk_c_kernel_smooth.bin &
+!- INPUT_GRADIENT/ contains:
+! proc000***_reg1_bulk_c_kernel_smooth.bin &
! proc000***_reg1_bulk_beta_kernel_smooth.bin &
-! proc000***_reg1_rho_kernel_smooth.bin
+! proc000***_reg1_rho_kernel_smooth.bin
! or
-! proc000***_reg1_alpha_kernel_smooth.bin &
+! proc000***_reg1_alpha_kernel_smooth.bin &
! proc000***_reg1_beta_kernel_smooth.bin &
-! proc000***_reg1_rho_kernel_smooth.bin
-!
+! proc000***_reg1_rho_kernel_smooth.bin
+!
!- topo/ contains:
! proc000***_reg1_solver_data_1.bin
!
! new models are stored in
-!- OUTPUT_MODEL/ as
-! proc000***_reg1_vp_new.bin and
+!- OUTPUT_MODEL/ as
+! proc000***_reg1_vp_new.bin and
! proc000***_reg1_vs_new.bin and
! proc000***_reg1_rho_new.bin and
!
@@ -47,9 +47,9 @@
include 'values_from_mesher_globe.h'
! ======================================================
-
+
! USER PARAMETERS
-
+
! by default, this algorithm uses (bulk,bulk_beta,rho) kernels to update vp,vs,rho
! if you prefer using (alpha,beta,rho) kernels, set this flag to true
logical, parameter :: USE_ALPHA_BETA_RHO = .false.
@@ -60,23 +60,23 @@
! in case of rho scaling, specifies density scaling factor with shear perturbations
! see e.g. Montagner & Anderson (1989), Panning & Romanowicz (2006)
real(kind=CUSTOM_REAL),parameter :: RHO_SCALING = 0.33_CUSTOM_REAL
-
+
! ======================================================
integer, parameter :: NSPEC = NSPEC_CRUST_MANTLE
integer, parameter :: NGLOB = NGLOB_CRUST_MANTLE
character(len=150) :: sline, m_file, fname
-
+
! model values, model updates and kernels
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: &
model_vp, model_vs, model_rho, &
model_dA, model_dB, model_dR, &
kernel_a, kernel_b, kernel_rho, total_model
- ! steepest descent lengths
+ ! steepest descent lengths
real(kind=CUSTOM_REAL) :: step_fac,step_length
-
+
integer :: nfile, myrank, sizeprocs, ier
integer :: i, j, k,ispec, iglob, ishell, n, it, j1, ib, npts_sem, ios
@@ -89,8 +89,8 @@
real(kind=CUSTOM_REAL), dimension(NGLOB) :: x, y, z
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: ibool
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: volume_spec
- real(kind=CUSTOM_REAL), dimension(NGLOB) :: volume_glob,volume_glob_sum
-
+ real(kind=CUSTOM_REAL), dimension(NGLOB) :: volume_glob,volume_glob_sum
+
! Gauss-Lobatto-Legendre points of integration and weights
double precision, dimension(NGLLX) :: xigll, wxgll
double precision, dimension(NGLLY) :: yigll, wygll
@@ -98,29 +98,29 @@
! array with all the weights in the cube
double precision, dimension(NGLLX,NGLLY,NGLLZ) :: wgll_cube
-
+
! jacobian
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: jacobian
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
real(kind=CUSTOM_REAL) xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl, &
jacobianl,volumel
-
+
! calculates volume associated with (global) point
logical,parameter:: VOLUME_GLOBALPOINT = .true.
real(kind=CUSTOM_REAL) :: kernel_integral_alpha,kernel_integral_beta,kernel_integral_rho
real(kind=CUSTOM_REAL) :: integral_alpha_sum,integral_beta_sum,integral_rho_sum
-
+
! ============ program starts here =====================
-
+
! initialize the MPI communicator and start the NPROCTOT MPI processes
call MPI_INIT(ier)
call MPI_COMM_SIZE(MPI_COMM_WORLD,sizeprocs,ier)
call MPI_COMM_RANK(MPI_COMM_WORLD,myrank,ier)
if( sizeprocs /= nchunks_val*nproc_xi_val*nproc_eta_val ) then
- print*,'sizeprocs:',sizeprocs,nchunks_val,nproc_xi_val,nproc_eta_val
- call exit_mpi('error number sizeprocs')
+ print*,'sizeprocs:',sizeprocs,nchunks_val,nproc_xi_val,nproc_eta_val
+ call exit_mpi(myrank,'error number sizeprocs')
endif
model_vp = 0.0_CUSTOM_REAL
@@ -132,10 +132,10 @@
kernel_a = 0.0_CUSTOM_REAL
kernel_b = 0.0_CUSTOM_REAL
kernel_rho = 0.0_CUSTOM_REAL
-
+
!--------------------------------------------------------
- ! subjective step length to multiply to the gradient
+ ! subjective step length to multiply to the gradient
! e.g. step_fac = 0.03
call getarg(1,s_step_fac)
@@ -150,11 +150,11 @@
print*,'error: step factor ',step_fac
call exit_MPI(myrank,'error step factor')
endif
-
+
if (myrank == 0) then
print*,'defaults'
- print*,' NPROC_XI , NPROC_ETA: ',nproc_xi_val,nproc_eta_val
- print*,' NCHUNKS: ',nchunks_val
+ print*,' NPROC_XI , NPROC_ETA: ',nproc_xi_val,nproc_eta_val
+ print*,' NCHUNKS: ',nchunks_val
print*
print*,'model update for vs & vp & rho'
print*,' step_fac = ',step_fac
@@ -169,21 +169,21 @@
print*,'scaling rho perturbations'
print*
endif
-
+
open(20,file='step_fac',status='unknown',action='write')
write(20,'(1e24.12)') step_fac
close(20)
endif
!-----------------------------------------------------
- ! reads in current vp & vs & rho model
+ ! reads in current vp & vs & rho model
! vp model
write(m_file,'(a,i6.6,a)') 'INPUT_MODEL/proc',myrank,'_reg1_vp.bin'
open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) then
print*,'error opening: ',trim(m_file)
- call exit_mpi('file not found')
- endif
+ call exit_mpi(myrank,'file not found')
+ endif
read(12) model_vp(:,:,:,1:nspec)
close(12)
@@ -192,8 +192,8 @@
open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) then
print*,'error opening: ',trim(m_file)
- call exit_mpi('file not found')
- endif
+ call exit_mpi(myrank,'file not found')
+ endif
read(12) model_vs(:,:,:,1:nspec)
close(12)
@@ -202,21 +202,21 @@
open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) then
print*,'error opening: ',trim(m_file)
- call exit_mpi('file not found')
- endif
+ call exit_mpi(myrank,'file not found')
+ endif
read(12) model_rho(:,:,:,1:nspec)
close(12)
-
+
! statistics
call mpi_reduce(minval(model_vp),min_vp,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
call mpi_reduce(maxval(model_vp),max_vp,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
-
+
call mpi_reduce(minval(model_vs),min_vs,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
call mpi_reduce(maxval(model_vs),max_vs,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
-
+
call mpi_reduce(minval(model_rho),min_rho,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
call mpi_reduce(maxval(model_rho),max_rho,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
-
+
if( myrank == 0 ) then
print*,'initial models:'
print*,' vs min/max: ',min_vs,max_vs
@@ -237,11 +237,11 @@
open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) then
print*,'error opening: ',trim(m_file)
- call exit_mpi('file not found')
- endif
+ call exit_mpi(myrank,'file not found')
+ endif
read(12) kernel_a(:,:,:,1:nspec)
close(12)
-
+
! beta kernel
if( USE_ALPHA_BETA_RHO ) then
! reads in beta kernel
@@ -254,29 +254,29 @@
open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) then
print*,'error opening: ',trim(m_file)
- call exit_mpi('file not found')
- endif
+ call exit_mpi(myrank,'file not found')
+ endif
read(12) kernel_b(:,:,:,1:nspec)
close(12)
! rho kernel
if( USE_RHO_SCALING ) then
-
+
! uses scaling relation with shear perturbations
kernel_rho(:,:,:,:) = RHO_SCALING * kernel_b(:,:,:,:)
-
+
else
-
- ! uses rho kernel
+
+ ! uses rho kernel
write(m_file,'(a,i6.6,a)') 'INPUT_GRADIENT/proc',myrank,'_reg1_rho_kernel_smooth.bin'
open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) then
print*,'error opening: ',trim(m_file)
- call exit_mpi('file not found')
- endif
+ call exit_mpi(myrank,'file not found')
+ endif
read(12) kernel_rho(:,:,:,1:nspec)
close(12)
- endif
+ endif
! statistics
call mpi_reduce(minval(kernel_a),min_vp,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
@@ -295,10 +295,10 @@
print*,' rho min/max: ',min_rho,max_rho
print*
endif
-
-
+
+
!--------------------------------------------------------
-
+
if( VOLUME_GLOBALPOINT ) then
! computes volume element associated with points
! GLL points
@@ -319,20 +319,20 @@
open(11,file=trim(m_file),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) then
print*,'error opening: ',trim(m_file)
- call exit_mpi('file not found')
+ call exit_mpi(myrank,'file not found')
endif
read(11) x(1:nglob)
read(11) y(1:nglob)
read(11) z(1:nglob)
read(11) ibool(:,:,:,1:nspec)
close(11)
-
- ! builds jacobian
+
+ ! builds jacobian
write(m_file,'(a,i6.6,a)') 'topo/proc',myrank,'_reg1_solver_data_1.bin'
open(11,file=trim(m_file),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) then
print*,'error opening: ',trim(m_file)
- call exit_mpi('file not found')
+ call exit_mpi(myrank,'file not found')
endif
read(11) xix
read(11) xiy
@@ -342,14 +342,14 @@
read(11) etaz
read(11) gammax
read(11) gammay
- read(11) gammaz
+ read(11) gammaz
close(11)
-
+
jacobian = 0.0
do ispec = 1, NSPEC
do k = 1, NGLLZ
do j = 1, NGLLY
- do i = 1, NGLLX
+ do i = 1, NGLLX
! gets derivatives of ux, uy and uz with respect to x, y and z
xixl = xix(i,j,k,ispec)
xiyl = xiy(i,j,k,ispec)
@@ -370,7 +370,7 @@
enddo
enddo
enddo
-
+
! volume associated with global point
volume_glob = 0.0
kernel_integral_alpha = 0._CUSTOM_REAL
@@ -379,7 +379,7 @@
do ispec = 1, NSPEC
do k = 1, NGLLZ
do j = 1, NGLLY
- do i = 1, NGLLX
+ do i = 1, NGLLX
iglob = ibool(i,j,k,ispec)
if( iglob == 0 ) then
print*,'iglob zero',i,j,k,ispec
@@ -387,13 +387,13 @@
print*,'ibool:',ispec
print*,ibool(:,:,:,ispec)
print*
- call exit_MPI('error ibool')
+ call exit_MPI(myrank,'error ibool')
endif
-
- ! volume associated with GLL point
- volumel = jacobian(i,j,k,ispec)*wgll_cube(i,j,k)
-
- ! kernel integration: for each element
+
+ ! volume associated with GLL point
+ volumel = jacobian(i,j,k,ispec)*wgll_cube(i,j,k)
+
+ ! kernel integration: for each element
kernel_integral_alpha = kernel_integral_alpha &
+ volumel * kernel_a(i,j,k,ispec)
@@ -402,12 +402,12 @@
kernel_integral_rho = kernel_integral_rho &
+ volumel * kernel_rho(i,j,k,ispec)
-
+
enddo
enddo
enddo
enddo
-
+
! statistics
! kernel integration: for whole volume
call mpi_reduce(kernel_integral_alpha,integral_alpha_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
@@ -421,26 +421,26 @@
print*,' rho : ',integral_rho_sum
print*
endif
-
+
endif
-
+
! calculates gradient
do ispec = 1, NSPEC
do k = 1, NGLLZ
do j = 1, NGLLY
- do i = 1, NGLLX
+ do i = 1, NGLLX
! gradient in direction of steepest descent
-
+
! for vp
model_dA(i,j,k,ispec) = - kernel_a(i,j,k,ispec)
-
+
! for shear
model_dB(i,j,k,ispec) = - kernel_b(i,j,k,ispec)
! for rho
model_dR(i,j,k,ispec) = - kernel_rho(i,j,k,ispec)
-
+
enddo
enddo
enddo
@@ -463,15 +463,15 @@
print*,' rho min/max: ',min_rho,max_rho
print*
endif
-
+
! master determines step length based on maximum gradient value (either vp or vs)
if( myrank == 0 ) then
minmax(1) = abs(min_vs)
minmax(2) = abs(max_vs)
-
+
minmax(3) = abs(min_vp)
minmax(4) = abs(max_vp)
-
+
max = maxval(minmax)
step_length = step_fac/max
@@ -481,18 +481,18 @@
call mpi_bcast(step_length,1,CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
- ! gradient length
- max_vp = sum( model_dA * model_dA )
+ ! gradient length
+ max_vp = sum( model_dA * model_dA )
max_vs = sum( model_dB * model_dB )
max_rho = sum( model_dR * model_dR )
max_vp = sqrt(max_vp)
max_vs = sqrt(max_vs)
max_rho = sqrt(max_rho)
- call mpi_reduce(max_vp,vp_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(max_vp,vp_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
call mpi_reduce(max_vs,vs_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
call mpi_reduce(max_rho,rho_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
-
+
if( myrank == 0 ) then
print*,' initial beta length : ',vs_sum
print*,' initial a length: ',vp_sum
@@ -502,7 +502,7 @@
! model updates:
! multiply model updates by a subjective factor that will change the step
-
+
model_dA(:,:,:,:) = step_length * model_dA(:,:,:,:)
model_dB(:,:,:,:) = step_length * model_dB(:,:,:,:)
model_dR(:,:,:,:) = step_length * model_dR(:,:,:,:)
@@ -530,9 +530,9 @@
! new model:
! computes new model values for alpha, beta and rho
! and stores new model files
-
+
! S wavespeed model
- total_model = 0._CUSTOM_REAL
+ total_model = 0._CUSTOM_REAL
total_model = model_vs * exp( model_dB )
call mpi_reduce(maxval(total_model),max_vs,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
@@ -545,9 +545,9 @@
close(12)
! Rho density model
- total_model = 0._CUSTOM_REAL
+ total_model = 0._CUSTOM_REAL
total_model = model_rho * exp( model_dR )
-
+
call mpi_reduce(maxval(total_model),max_rho,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
call mpi_reduce(minval(total_model),min_rho,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
@@ -558,12 +558,12 @@
close(12)
! P wavespeed model
- total_model = 0._CUSTOM_REAL
+ total_model = 0._CUSTOM_REAL
if( USE_ALPHA_BETA_RHO ) then
! new vp values use alpha model update
total_model = model_vp * exp( model_dA )
else
- ! new vp values use bulk model update:
+ ! new vp values use bulk model update:
! this is based on vp_new = sqrt( bulk_new**2 + 4/3 vs_new**2 )
total_model = sqrt( model_vp**2 * exp(2.0*model_dA) + FOUR_THIRDS * model_vs**2 *( &
exp(2.0*model_dB) - exp(2.0*model_dA) ) )
@@ -596,16 +596,16 @@
write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_dvpvp.bin'
open(12,file=trim(m_file),form='unformatted',action='write')
write(12) total_model(:,:,:,1:nspec)
- close(12)
+ close(12)
call mpi_reduce(maxval(total_model),max_vp,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
call mpi_reduce(minval(total_model),min_vp,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
! relative S model perturbations
- total_model = 0._CUSTOM_REAL
+ total_model = 0._CUSTOM_REAL
do ispec = 1, NSPEC
do k = 1, NGLLZ
do j = 1, NGLLY
- do i = 1, NGLLX
+ do i = 1, NGLLX
if ( model_vs(i,j,k,ispec) > 1.e-10 ) then
total_model(i,j,k,ispec) = ( model_vs(i,j,k,ispec)*exp(model_dB(i,j,k,ispec)) &
- model_vs(i,j,k,ispec) ) / model_vs(i,j,k,ispec)
@@ -617,16 +617,16 @@
write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_dvsvs.bin'
open(12,file=trim(m_file),form='unformatted',action='write')
write(12) total_model(:,:,:,1:nspec)
- close(12)
+ close(12)
call mpi_reduce(maxval(total_model),max_vs,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
call mpi_reduce(minval(total_model),min_vs,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
! relative rho model perturbations
- total_model = 0._CUSTOM_REAL
+ total_model = 0._CUSTOM_REAL
do ispec = 1, NSPEC
do k = 1, NGLLZ
do j = 1, NGLLY
- do i = 1, NGLLX
+ do i = 1, NGLLX
if ( model_rho(i,j,k,ispec) > 1.e-10 ) then
total_model(i,j,k,ispec) = ( model_rho(i,j,k,ispec)*exp(model_dR(i,j,k,ispec)) &
- model_rho(i,j,k,ispec) ) / model_rho(i,j,k,ispec)
@@ -641,8 +641,8 @@
close(12)
call mpi_reduce(maxval(total_model),max_rho,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
call mpi_reduce(minval(total_model),min_rho,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
-
+
if( myrank == 0 ) then
print*,'relative Vs & Vp update:'
print*,' dvs/vs min/max: ',min_vs,max_vs
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/add_model_globe_tiso.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/add_model_globe_tiso.f90 2011-02-16 15:53:48 UTC (rev 17884)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/add_model_globe_tiso.f90 2011-02-16 21:49:17 UTC (rev 17885)
@@ -1,12 +1,12 @@
! add_model_globe_tiso
!
-! this program can be used to update TRANSVERSE ISOTROPIC model files
+! this program can be used to update TRANSVERSE ISOTROPIC model files
! based on smoothed event kernels.
! the kernels are given for tranverse isotropic parameters (bulk_c,bulk_betav,bulk_betah,eta).
!
-! the algorithm uses a steepest descent method with a step length
+! the algorithm uses a steepest descent method with a step length
! determined by the given maximum update percentage.
-!
+!
! input:
! - step_fac : step length to update the models, f.e. 0.03 for plusminus 3%
!
@@ -24,7 +24,7 @@
! proc000***_reg1_bulk_c_kernel_smooth.bin &
! proc000***_reg1_bulk_betav_kernel_smooth.bin &
! proc000***_reg1_bulk_betah_kernel_smooth.bin &
-! proc000***_reg1_eta_kernel_smooth.bin
+! proc000***_reg1_eta_kernel_smooth.bin
!
!- topo/ contains:
! proc000***_reg1_solver_data_1.bin
@@ -40,7 +40,7 @@
!
! USAGE: ./add_model_globe_tiso 0.3
-module model_update
+module model_update_tiso
include 'mpif.h'
include 'constants_globe.h'
@@ -48,7 +48,7 @@
include 'values_from_mesher_globe.h'
! ======================================================
-
+
! density scaling factor with shear perturbations
! see e.g. Montagner & Anderson (1989), Panning & Romanowicz (2006)
real(kind=CUSTOM_REAL),parameter :: RHO_SCALING = 0.33_CUSTOM_REAL
@@ -75,24 +75,12 @@
! kernels
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: &
kernel_bulk,kernel_betav,kernel_betah,kernel_eta
-
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: total_model
! volume
real(kind=CUSTOM_REAL), dimension(NGLOB) :: x, y, z
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: volume_spec
- real(kind=CUSTOM_REAL), dimension(NGLOB) :: volume_glob,volume_glob_sum
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: ibool
integer, dimension(NSPEC) :: idoubling
- ! Gauss-Lobatto-Legendre points of integration and weights
- double precision, dimension(NGLLX) :: xigll, wxgll
- double precision, dimension(NGLLY) :: yigll, wygll
- double precision, dimension(NGLLZ) :: zigll, wzgll
-
- ! array with all the weights in the cube
- double precision, dimension(NGLLX,NGLLY,NGLLZ) :: wgll_cube
-
! gradient vector norm ( v^T * v )
real(kind=CUSTOM_REAL) :: norm_bulk,norm_betav,norm_betah,norm_eta
real(kind=CUSTOM_REAL) :: norm_bulk_sum,norm_betav_sum, &
@@ -104,17 +92,17 @@
real(kind=CUSTOM_REAL) :: min_vpv,min_vph,min_vsv,min_vsh, &
max_vpv,max_vph,max_vsv,max_vsh,min_eta,max_eta,min_bulk,max_bulk, &
min_rho,max_rho,max,minmax(4)
-
+
real(kind=CUSTOM_REAL) :: betav1,betah1,betav0,betah0,rho1,rho0, &
betaiso1,betaiso0,eta1,eta0,alphav1,alphav0,alphah1,alphah0
real(kind=CUSTOM_REAL) :: dbetaiso,dbulk
-
+
integer :: nfile, myrank, sizeprocs, ier
integer :: i, j, k,ispec, iglob, ishell, n, it, j1, ib, npts_sem, ios
character(len=150) :: sline, m_file, fname
-end module model_update
+end module model_update_tiso
!
!-------------------------------------------------------------------------------------------------
@@ -122,7 +110,7 @@
program add_model
- use model_update
+ use model_update_tiso
implicit none
@@ -134,9 +122,9 @@
! reads in parameters needed
call read_parameters()
- ! reads in current transverse isotropic model files: vpv.. & vsv.. & eta & rho
+ ! reads in current transverse isotropic model files: vpv.. & vsv.. & eta & rho
call read_model()
-
+
! reads in smoothed kernels: bulk, betav, betah, eta
call read_kernels()
@@ -146,8 +134,8 @@
! calculates gradient
! steepest descent method
call get_gradient()
-
- ! compute new model in terms of alpha, beta, eta and rho
+
+ ! compute new model in terms of alpha, beta, eta and rho
! (see also Carl's Latex notes)
! model update:
@@ -156,8 +144,8 @@
do ispec = 1, NSPEC
do k = 1, NGLLZ
do j = 1, NGLLY
- do i = 1, NGLLX
-
+ do i = 1, NGLLX
+
! initial model values
eta0 = model_eta(i,j,k,ispec)
betav0 = model_vsv(i,j,k,ispec)
@@ -172,47 +160,47 @@
rho1 = 0._CUSTOM_REAL
alphav1 = 0._CUSTOM_REAL
alphah1 = 0._CUSTOM_REAL
-
+
! do not use transverse isotropy except if element is between d220 and Moho
if(.not. ( idoubling(ispec)==IFLAG_220_80 .or. idoubling(ispec)==IFLAG_80_MOHO) ) then
-
+
! isotropic model update
-
+
! no eta perturbation, since eta = 1 in isotropic media
eta1 = eta0
-
+
! shear values
- ! isotropic kernel K_beta = K_betav + K_betah
+ ! isotropic kernel K_beta = K_betav + K_betah
! with same scaling step_length the model update dbeta_iso = dbetav + dbetah
- ! note:
+ ! note:
! this step length can be twice as big as that given by the input
dbetaiso = model_dbetav(i,j,k,ispec) + model_dbetah(i,j,k,ispec)
betav1 = betav0 * exp( dbetaiso )
- betah1 = betah0 * exp( dbetaiso )
- ! note: betah is probably not really used in isotropic layers
- ! (see SPECFEM3D_GLOBE/get_model.f90)
-
+ betah1 = betah0 * exp( dbetaiso )
+ ! note: betah is probably not really used in isotropic layers
+ ! (see SPECFEM3D_GLOBE/get_model.f90)
+
! density: uses scaling relation with isotropic shear perturbations
! dln rho = RHO_SCALING * dln betaiso
- rho1 = rho0 * exp( RHO_SCALING * dbetaiso )
-
+ rho1 = rho0 * exp( RHO_SCALING * dbetaiso )
+
! alpha values
dbulk = model_dbulk(i,j,k,ispec)
alphav1 = sqrt( alphav0**2 * exp(2.0*dbulk) + FOUR_THIRDS * betav0**2 * ( &
- exp(2.0*dbetaiso) - exp(2.0*dbulk) ) )
+ exp(2.0*dbetaiso) - exp(2.0*dbulk) ) )
alphah1 = sqrt( alphah0**2 * exp(2.0*dbulk) + FOUR_THIRDS * betah0**2 * ( &
- exp(2.0*dbetaiso) - exp(2.0*dbulk) ) )
- ! note: alphah probably not used in SPECFEM3D_GLOBE
-
+ exp(2.0*dbetaiso) - exp(2.0*dbulk) ) )
+ ! note: alphah probably not used in SPECFEM3D_GLOBE
+
else
-
+
! transverse isotropic model update
-
+
! eta value : limits updated values for eta range constraint
- eta1 = eta0 * exp( model_deta(i,j,k,ispec) )
+ eta1 = eta0 * exp( model_deta(i,j,k,ispec) )
if( eta1 < LIMIT_ETA_MIN ) eta1 = LIMIT_ETA_MIN
if( eta1 > LIMIT_ETA_MAX ) eta1 = LIMIT_ETA_MAX
-
+
! shear values
betav1 = betav0 * exp( model_dbetav(i,j,k,ispec) )
betah1 = betah0 * exp( model_dbetah(i,j,k,ispec) )
@@ -222,16 +210,16 @@
betaiso1 = sqrt( ( 2.0 * betav1**2 + betah1**2 ) / 3.0 )
dbetaiso = log( betaiso1 / betaiso0 )
rho1 = rho0 * exp( RHO_SCALING * dbetaiso )
-
+
! alpha values
- dbulk = model_dbulk(i,j,k,ispec)
+ dbulk = model_dbulk(i,j,k,ispec)
alphav1 = sqrt( alphav0**2 * exp(2.0*dbulk) &
+ FOUR_THIRDS * betav0**2 * ( &
- exp(2.0*model_dbetav(i,j,k,ispec)) - exp(2.0*dbulk) ) )
+ exp(2.0*model_dbetav(i,j,k,ispec)) - exp(2.0*dbulk) ) )
alphah1 = sqrt( alphah0**2 * exp(2.0*dbulk) &
+ FOUR_THIRDS * betah0**2 * ( &
- exp(2.0*model_dbetah(i,j,k,ispec)) - exp(2.0*dbulk) ) )
-
+ exp(2.0*model_dbetah(i,j,k,ispec)) - exp(2.0*dbulk) ) )
+
endif
@@ -242,7 +230,7 @@
model_vsh_new(i,j,k,ispec) = betah1
model_eta_new(i,j,k,ispec) = eta1
model_rho_new(i,j,k,ispec) = rho1
-
+
enddo
enddo
enddo
@@ -253,7 +241,7 @@
! stores relative model perturbations
call store_perturbations()
-
+
! stop all the MPI processes, and exit
call MPI_FINALIZE(ier)
@@ -267,8 +255,8 @@
! initializes arrays
- use model_update
- implicit none
+ use model_update_tiso
+ implicit none
! initialize the MPI communicator and start the NPROCTOT MPI processes
call MPI_INIT(ier)
@@ -277,7 +265,7 @@
if( sizeprocs /= nchunks_val*nproc_xi_val*nproc_eta_val ) then
print*,'sizeprocs:',sizeprocs,nchunks_val,nproc_xi_val,nproc_eta_val
- call exit_mpi('error number sizeprocs')
+ call exit_mpi(myrank,'error number sizeprocs')
endif
! model
@@ -317,10 +305,10 @@
! reads in parameters needed
- use model_update
+ use model_update_tiso
implicit none
character(len=150) :: s_step_fac
-
+
! subjective step length to multiply to the gradient
!step_fac = 0.03
@@ -357,9 +345,9 @@
subroutine read_model()
-! reads in current transverse isotropic model: vpv.. & vsv.. & eta & rho
+! reads in current transverse isotropic model: vpv.. & vsv.. & eta & rho
- use model_update
+ use model_update_tiso
implicit none
! vpv model
@@ -367,7 +355,7 @@
open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) then
print*,'error opening: ',trim(m_file)
- call exit_mpi('file not found')
+ call exit_mpi(myrank,'file not found')
endif
read(12) model_vpv(:,:,:,1:nspec)
close(12)
@@ -377,7 +365,7 @@
open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) then
print*,'error opening: ',trim(m_file)
- call exit_mpi('file not found')
+ call exit_mpi(myrank,'file not found')
endif
read(12) model_vph(:,:,:,1:nspec)
close(12)
@@ -387,7 +375,7 @@
open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) then
print*,'error opening: ',trim(m_file)
- call exit_mpi('file not found')
+ call exit_mpi(myrank,'file not found')
endif
read(12) model_vsv(:,:,:,1:nspec)
close(12)
@@ -397,7 +385,7 @@
open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) then
print*,'error opening: ',trim(m_file)
- call exit_mpi('file not found')
+ call exit_mpi(myrank,'file not found')
endif
read(12) model_vsh(:,:,:,1:nspec)
close(12)
@@ -407,7 +395,7 @@
open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) then
print*,'error opening: ',trim(m_file)
- call exit_mpi('file not found')
+ call exit_mpi(myrank,'file not found')
endif
read(12) model_eta(:,:,:,1:nspec)
close(12)
@@ -417,7 +405,7 @@
open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) then
print*,'error opening: ',trim(m_file)
- call exit_mpi('file not found')
+ call exit_mpi(myrank,'file not found')
endif
read(12) model_rho(:,:,:,1:nspec)
close(12)
@@ -462,7 +450,7 @@
! reads in smoothed kernels: bulk, betav, betah, eta
- use model_update
+ use model_update_tiso
implicit none
! bulk kernel
@@ -470,7 +458,7 @@
open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) then
print*,'error opening: ',trim(m_file)
- call exit_mpi('file not found')
+ call exit_mpi(myrank,'file not found')
endif
read(12) kernel_bulk(:,:,:,1:nspec)
close(12)
@@ -480,7 +468,7 @@
open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) then
print*,'error opening: ',trim(m_file)
- call exit_mpi('file not found')
+ call exit_mpi(myrank,'file not found')
endif
read(12) kernel_betav(:,:,:,1:nspec)
close(12)
@@ -490,7 +478,7 @@
open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) then
print*,'error opening: ',trim(m_file)
- call exit_mpi('file not found')
+ call exit_mpi(myrank,'file not found')
endif
read(12) kernel_betah(:,:,:,1:nspec)
close(12)
@@ -500,7 +488,7 @@
open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) then
print*,'error opening: ',trim(m_file)
- call exit_mpi('file not found')
+ call exit_mpi(myrank,'file not found')
endif
read(12) kernel_eta(:,:,:,1:nspec)
close(12)
@@ -538,7 +526,7 @@
! computes volume element associated with points
- use model_update
+ use model_update_tiso
implicit none
! jacobian
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: jacobian
@@ -549,11 +537,18 @@
! integration values
real(kind=CUSTOM_REAL) :: integral_bulk_sum,integral_betav_sum, &
integral_betah_sum,integral_eta_sum
- real(kind=CUSTOM_REAL) :: kernel_integral_bulk,kernel_integral_betav,&
- kernel_integral_betah,kernel_integral_eta
+ real(kind=CUSTOM_REAL) :: integral_bulk,integral_betav,&
+ integral_betah,integral_eta
+ real(kind=CUSTOM_REAL) :: volume_glob,volume_glob_sum
+ ! Gauss-Lobatto-Legendre points of integration and weights
+ double precision, dimension(NGLLX) :: xigll, wxgll
+ double precision, dimension(NGLLY) :: yigll, wygll
+ double precision, dimension(NGLLZ) :: zigll, wzgll
+ ! array with all the weights in the cube
+ double precision, dimension(NGLLX,NGLLY,NGLLZ) :: wgll_cube
! GLL points
- wgll_cube = 0.0
+ wgll_cube = 0.0d0
call zwgljd(xigll,wxgll,NGLLX,GAUSSALPHA,GAUSSBETA)
call zwgljd(yigll,wygll,NGLLY,GAUSSALPHA,GAUSSBETA)
call zwgljd(zigll,wzgll,NGLLZ,GAUSSALPHA,GAUSSBETA)
@@ -570,7 +565,7 @@
open(11,file=trim(m_file),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) then
print*,'error opening: ',trim(m_file)
- call exit_mpi('file not found')
+ call exit_mpi(myrank,'file not found')
endif
read(11) x(1:nglob)
read(11) y(1:nglob)
@@ -584,7 +579,7 @@
open(11,file=trim(m_file),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) then
print*,'error opening: ',trim(m_file)
- call exit_mpi('file not found')
+ call exit_mpi(myrank,'file not found')
endif
read(11) xix
read(11) xiy
@@ -629,10 +624,14 @@
! volume associated with global point
volume_glob = 0.0
- kernel_integral_bulk = 0._CUSTOM_REAL
- kernel_integral_betav = 0._CUSTOM_REAL
- kernel_integral_betah = 0._CUSTOM_REAL
- kernel_integral_eta = 0._CUSTOM_REAL
+ integral_bulk = 0._CUSTOM_REAL
+ integral_betav = 0._CUSTOM_REAL
+ integral_betah = 0._CUSTOM_REAL
+ integral_eta = 0._CUSTOM_REAL
+ norm_bulk = 0._CUSTOM_REAL
+ norm_betav = 0._CUSTOM_REAL
+ norm_betah = 0._CUSTOM_REAL
+ norm_eta = 0._CUSTOM_REAL
do ispec = 1, NSPEC
do k = 1, NGLLZ
do j = 1, NGLLY
@@ -644,31 +643,41 @@
print*,'ibool:',ispec
print*,ibool(:,:,:,ispec)
print*
- call exit_MPI('error ibool')
+ call exit_MPI(myrank,'error ibool')
endif
! volume associated with GLL point
volumel = jacobian(i,j,k,ispec)*wgll_cube(i,j,k)
+ volume_glob = volume_glob + volumel
! kernel integration: for each element
- kernel_integral_bulk = kernel_integral_bulk &
+ integral_bulk = integral_bulk &
+ volumel * kernel_bulk(i,j,k,ispec)
- kernel_integral_betav = kernel_integral_betav &
+ integral_betav = integral_betav &
+ volumel * kernel_betav(i,j,k,ispec)
- kernel_integral_betah = kernel_integral_betah &
+ integral_betah = integral_betah &
+ volumel * kernel_betah(i,j,k,ispec)
- kernel_integral_eta = kernel_integral_eta &
+ integral_eta = integral_eta &
+ volumel * kernel_eta(i,j,k,ispec)
! gradient vector norm sqrt( v^T * v )
- norm_bulk = norm_bulk + kernel_bulk(i,j,k,ispec)**2
- norm_betav = norm_betav + kernel_betav(i,j,k,ispec)**2
- norm_betah = norm_betah + kernel_betah(i,j,k,ispec)**2
- norm_eta = norm_eta + kernel_eta(i,j,k,ispec)**2
+ norm_bulk = norm_bulk + kernel_bulk(i,j,k,ispec)*kernel_bulk(i,j,k,ispec)
+ norm_betav = norm_betav + kernel_betav(i,j,k,ispec)*kernel_betav(i,j,k,ispec)
+ norm_betah = norm_betah + kernel_betah(i,j,k,ispec)*kernel_betah(i,j,k,ispec)
+ norm_eta = norm_eta + kernel_eta(i,j,k,ispec)*kernel_eta(i,j,k,ispec)
+ ! checks number
+ if( isNaN(integral_bulk) ) then
+ print*,'error NaN: ',integral_bulk
+ print*,'rank:',myrank
+ print*,'i,j,k,ispec:',i,j,k,ispec
+ print*,'volumel: ',volumel,'kernel_bulk:',kernel_bulk(i,j,k,ispec)
+ call exit_MPI(myrank,'error NaN')
+ endif
+
enddo
enddo
enddo
@@ -676,10 +685,11 @@
! statistics
! kernel integration: for whole volume
- call mpi_reduce(kernel_integral_bulk,integral_bulk_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
- call mpi_reduce(kernel_integral_betav,integral_betav_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
- call mpi_reduce(kernel_integral_betah,integral_betah_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
- call mpi_reduce(kernel_integral_eta,integral_eta_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(integral_bulk,integral_bulk_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(integral_betav,integral_betav_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(integral_betah,integral_betah_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(integral_eta,integral_eta_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(volume_glob,volume_glob_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
if( myrank == 0 ) then
print*,'integral kernels:'
@@ -688,6 +698,8 @@
print*,' betah : ',integral_betah_sum
print*,' eta : ',integral_eta_sum
print*
+ print*,' total volume:',volume_glob_sum
+ print*
endif
! norms: for whole volume
@@ -701,7 +713,6 @@
norm_betah = sqrt(norm_betah_sum)
norm_eta = sqrt(norm_eta_sum)
-
if( myrank == 0 ) then
print*,'norm kernels:'
print*,' bulk : ',norm_bulk
@@ -721,9 +732,21 @@
! calculates gradient by steepest descent method
- use model_update
+ use model_update_tiso
implicit none
+ ! local parameters
+ ! ------------------------------------------------------------------------
+ ! sets maximum update in this depth range
+ logical,parameter :: use_depth_maximum = .true.
+ ! normalized radii
+ real(kind=CUSTOM_REAL),parameter :: R_top = (6371.0 - 50.0 ) / R_EARTH_KM ! shallow depth
+ real(kind=CUSTOM_REAL),parameter :: R_bottom = (6371.0 - 100.0 ) / R_EARTH_KM ! deep depth
+ real(kind=CUSTOM_REAL):: r,depth_max
+ ! ------------------------------------------------------------------------
+ ! initializes kernel maximum
+ max = 0._CUSTOM_REAL
+
! gradient in negative direction for steepest descent
do ispec = 1, NSPEC
do k = 1, NGLLZ
@@ -740,6 +763,23 @@
! for eta
model_deta(i,j,k,ispec) = - kernel_eta(i,j,k,ispec)
+ ! determines maximum kernel betav value within given radius
+ if( use_depth_maximum ) then
+ ! get radius of point
+ iglob = ibool(i,j,k,ispec)
+ r = sqrt( x(iglob)*x(iglob) + y(iglob)*y(iglob) + z(iglob)*z(iglob) )
+
+ ! stores maximum kernel betav value in this depth slice, since betav is most likely dominating
+ if( r < R_top .and. r > R_bottom ) then
+ ! kernel betav value
+ max_vsv = abs( kernel_betav(i,j,k,ispec) )
+ if( max < max_vsv ) then
+ max = max_vsv
+ depth_max = r
+ endif
+ endif
+ endif
+
enddo
enddo
enddo
@@ -769,21 +809,43 @@
print*,' eta min/max : ',min_eta,max_eta
print*
endif
-
- ! determines step length based on maximum gradient value (either vsv or vsh)
+
+ ! determines maximum kernel betav value within given radius
+ if( use_depth_maximum ) then
+ ! maximum of all processes stored in max_vsv
+ call mpi_reduce(max,max_vsv,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
+ max = max_vsv
+ depth_max = 6371.0 *( 1.0 - depth_max )
+ endif
+
+ ! determines step length
+ ! based on maximum gradient value (either vsv or vsh)
if( myrank == 0 ) then
- ! maximum gradient values
- minmax(1) = abs(min_vsv)
- minmax(2) = abs(max_vsv)
- minmax(3) = abs(min_vsh)
- minmax(4) = abs(max_vsh)
+ ! determines maximum kernel betav value within given radius
+ if( use_depth_maximum ) then
+ print*,' using depth maximum between 50km and 100km: ',max
+ print*,' approximate depth maximum: ',depth_max
+ print*
+ else
+ ! maximum gradient values
+ minmax(1) = abs(min_vsv)
+ minmax(2) = abs(max_vsv)
+ minmax(3) = abs(min_vsh)
+ minmax(4) = abs(max_vsh)
+
+ ! maximum value of all kernel maxima
+ max = maxval(minmax)
+ print*,' using maximum: ',max
+ print*
+ endif
+
! chooses step length such that it becomes the desired, given step factor as inputted
- max = maxval(minmax)
step_length = step_fac/max
- print*,' step length : ',step_length,max
+ print*,' step length : ',step_length
print*
+
endif
call mpi_bcast(step_length,1,CUSTOM_MPI_TYPE,0,MPI_COMM_WORLD,ier)
@@ -850,33 +912,33 @@
subroutine store_kernel_updates()
-! file output for new model
+! file output for new model
- use model_update
- implicit none
+ use model_update_tiso
+ implicit none
- ! kernel updates
+ ! kernel updates
fname = 'dbulk_c'
write(m_file,'(a,i6.6,a)') 'INPUT_GRADIENT/proc',myrank,'_reg1_'//trim(fname)//'.bin'
- open(12,file=trim(m_file),form='unformatted')
+ open(12,file=trim(m_file),form='unformatted',action='write')
write(12) model_dbulk
close(12)
fname = 'dbetav'
write(m_file,'(a,i6.6,a)') 'INPUT_GRADIENT/proc',myrank,'_reg1_'//trim(fname)//'.bin'
- open(12,file=trim(m_file),form='unformatted')
+ open(12,file=trim(m_file),form='unformatted',action='write')
write(12) model_dbetav
close(12)
fname = 'dbetah'
write(m_file,'(a,i6.6,a)') 'INPUT_GRADIENT/proc',myrank,'_reg1_'//trim(fname)//'.bin'
- open(12,file=trim(m_file),form='unformatted')
+ open(12,file=trim(m_file),form='unformatted',action='write')
write(12) model_dbetah
close(12)
fname = 'deta'
write(m_file,'(a,i6.6,a)') 'INPUT_GRADIENT/proc',myrank,'_reg1_'//trim(fname)//'.bin'
- open(12,file=trim(m_file),form='unformatted')
+ open(12,file=trim(m_file),form='unformatted',action='write')
write(12) model_deta
close(12)
@@ -888,62 +950,62 @@
subroutine store_new_model()
-! file output for new model
+! file output for new model
- use model_update
- implicit none
+ use model_update_tiso
+ implicit none
- ! vpv model
+ ! vpv model
call mpi_reduce(maxval(model_vpv_new),max_vpv,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
call mpi_reduce(minval(model_vpv_new),min_vpv,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
fname = 'vpv_new'
write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_'//trim(fname)//'.bin'
- open(12,file=trim(m_file),form='unformatted')
+ open(12,file=trim(m_file),form='unformatted',action='write')
write(12) model_vpv_new
close(12)
- ! vph model
+ ! vph model
call mpi_reduce(maxval(model_vph_new),max_vph,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
call mpi_reduce(minval(model_vph_new),min_vph,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
fname = 'vph_new'
write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_'//trim(fname)//'.bin'
- open(12,file=trim(m_file),form='unformatted')
+ open(12,file=trim(m_file),form='unformatted',action='write')
write(12) model_vph_new
close(12)
- ! vsv model
+ ! vsv model
call mpi_reduce(maxval(model_vsv_new),max_vsv,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
call mpi_reduce(minval(model_vsv_new),min_vsv,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
fname = 'vsv_new'
write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_'//trim(fname)//'.bin'
- open(12,file=trim(m_file),form='unformatted')
+ open(12,file=trim(m_file),form='unformatted',action='write')
write(12) model_vsv_new
close(12)
- ! vsh model
+ ! vsh model
call mpi_reduce(maxval(model_vsh_new),max_vsh,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
call mpi_reduce(minval(model_vsh_new),min_vsh,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
fname = 'vsh_new'
write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_'//trim(fname)//'.bin'
- open(12,file=trim(m_file),form='unformatted')
+ open(12,file=trim(m_file),form='unformatted',action='write')
write(12) model_vsh_new
close(12)
- ! eta model
+ ! eta model
call mpi_reduce(maxval(model_eta_new),max_eta,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
call mpi_reduce(minval(model_eta_new),min_eta,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
fname = 'eta_new'
write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_'//trim(fname)//'.bin'
- open(12,file=trim(m_file),form='unformatted')
+ open(12,file=trim(m_file),form='unformatted',action='write')
write(12) model_eta_new
close(12)
- ! rho model
+ ! rho model
call mpi_reduce(maxval(model_rho_new),max_rho,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
call mpi_reduce(minval(model_rho_new),min_rho,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
fname = 'rho_new'
write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_'//trim(fname)//'.bin'
- open(12,file=trim(m_file),form='unformatted')
+ open(12,file=trim(m_file),form='unformatted',action='write')
write(12) model_rho_new
close(12)
@@ -968,12 +1030,13 @@
subroutine store_perturbations()
-! file output for new model
+! file output for new model
- use model_update
- implicit none
+ use model_update_tiso
+ implicit none
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: total_model
- ! vpv relative perturbations
+ ! vpv relative perturbations
! logarithmic perturbation: log( v_new) - log( v_old) = log( v_new / v_old )
total_model = 0.0_CUSTOM_REAL
where( model_vpv /= 0.0 ) total_model = log( model_vpv_new / model_vpv)
@@ -990,7 +1053,7 @@
! vph relative perturbations
total_model = 0.0_CUSTOM_REAL
- where( model_vph /= 0.0 ) total_model = log( model_vph_new / model_vph)
+ where( model_vph /= 0.0 ) total_model = log( model_vph_new / model_vph)
write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_dvphvph.bin'
open(12,file=trim(m_file),form='unformatted',action='write')
write(12) total_model
@@ -1000,7 +1063,7 @@
! vsv relative perturbations
total_model = 0.0_CUSTOM_REAL
- where( model_vsv /= 0.0 ) total_model = log( model_vsv_new / model_vsv)
+ where( model_vsv /= 0.0 ) total_model = log( model_vsv_new / model_vsv)
write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_dvsvvsv.bin'
open(12,file=trim(m_file),form='unformatted',action='write')
write(12) total_model
@@ -1030,7 +1093,7 @@
! rho relative model perturbations
total_model = 0.0_CUSTOM_REAL
- where( model_rho /= 0.0 ) total_model = log( model_rho_new / model_rho)
+ where( model_rho /= 0.0 ) total_model = log( model_rho_new / model_rho)
write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_drhorho.bin'
open(12,file=trim(m_file),form='unformatted',action='write')
write(12) total_model
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/add_model_globe_tiso_iso.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/add_model_globe_tiso_iso.f90 2011-02-16 15:53:48 UTC (rev 17884)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/add_model_globe_tiso_iso.f90 2011-02-16 21:49:17 UTC (rev 17885)
@@ -1,12 +1,12 @@
! add_model_globe_tiso_iso
!
-! this program can be used to update TRANSVERSE ISOTROPIC model files
+! this program can be used to update TRANSVERSE ISOTROPIC model files
! based on smoothed, summed, ISOTROPIC event kernels.
! the kernels are given for isotropic parameters (bulk_c,bulk_beta,rho) or (alpha,beta,rho).
!
-! the algorithm uses a steepest descent method with a step length
+! the algorithm uses a steepest descent method with a step length
! determined by the given maximum update percentage.
-!
+!
! input:
! - step_fac : step length to update the models, f.e. 0.03 for plusminus 3%
!
@@ -23,11 +23,11 @@
!- INPUT_GRADIENT/ contains:
! proc000***_reg1_bulk_c_kernel_smooth.bin &
! proc000***_reg1_bulk_beta_kernel_smooth.bin &
-! proc000***_reg1_rho_kernel_smooth.bin
+! proc000***_reg1_rho_kernel_smooth.bin
! or
! proc000***_reg1_alpha_kernel_smooth.bin &
! proc000***_reg1_beta_kernel_smooth.bin &
-! proc000***_reg1_rho_kernel_smooth.bin
+! proc000***_reg1_rho_kernel_smooth.bin
!
!- topo/ contains:
! proc000***_reg1_solver_data_1.bin
@@ -53,7 +53,7 @@
! ======================================================
! USER PARAMETERS
-
+
! by default, this algorithm uses (bulk_c,bulk_beta,rho) kernels to update vpv,vph,vsv,vsh,rho,eta
! if you prefer using (alpha,beta,rho) kernels, set this flag to true
logical, parameter :: USE_ALPHA_BETA_RHO = .false.
@@ -83,7 +83,7 @@
! kernels
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: &
kernel_a,kernel_b,kernel_rho
-
+
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: total_model
! volume
@@ -111,11 +111,11 @@
real(kind=CUSTOM_REAL) :: min_vpv,min_vph,min_vsv,min_vsh, &
max_vpv,max_vph,max_vsv,max_vsh,min_b,max_b,min_a,max_a, &
min_rho,max_rho,min_eta,max_eta,max,minmax(4)
-
+
real(kind=CUSTOM_REAL) :: betav1,betah1,betav0,betah0,rho1,rho0, &
betaiso1,betaiso0,eta1,eta0,alphav1,alphav0,alphah1,alphah0
real(kind=CUSTOM_REAL) :: dbetaiso,dbulk
-
+
integer :: nfile, myrank, sizeprocs, ier
integer :: i, j, k,ispec, iglob, ishell, n, it, j1, ib, npts_sem, ios
character(len=150) :: sline, m_file, fname
@@ -141,9 +141,9 @@
! reads in parameters needed
call read_parameters()
- ! reads in current transverse isotropic model files: vpv.. & vsv.. & eta & rho
+ ! reads in current transverse isotropic model files: vpv.. & vsv.. & eta & rho
call read_model()
-
+
! reads in smoothed kernels: bulk, beta, rho
call read_kernels()
@@ -153,8 +153,8 @@
! calculates gradient
! steepest descent method
call get_gradient()
-
- ! compute new model in terms of alpha, beta, eta and rho
+
+ ! compute new model in terms of alpha, beta, eta and rho
! (see also Carl's Latex notes)
! model update:
@@ -163,8 +163,8 @@
do ispec = 1, NSPEC
do k = 1, NGLLZ
do j = 1, NGLLY
- do i = 1, NGLLX
-
+ do i = 1, NGLLX
+
! initial model values
eta0 = model_eta(i,j,k,ispec)
betav0 = model_vsv(i,j,k,ispec)
@@ -179,29 +179,29 @@
rho1 = 0._CUSTOM_REAL
alphav1 = 0._CUSTOM_REAL
alphah1 = 0._CUSTOM_REAL
-
+
! isotropic model update
-
+
! no eta perturbation for isotropic updates, takes the same as initial values
eta1 = eta0
-
+
! shear values
dbetaiso = model_dB(i,j,k,ispec)
betav1 = betav0 * exp( dbetaiso )
- betah1 = betah0 * exp( dbetaiso )
- ! note: betah is probably not really used in isotropic layers
- ! (see SPECFEM3D_GLOBE/get_model.f90)
-
+ betah1 = betah0 * exp( dbetaiso )
+ ! note: betah is probably not really used in isotropic layers
+ ! (see SPECFEM3D_GLOBE/get_model.f90)
+
! density
- rho1 = rho0 * exp( model_dR(i,j,k,ispec) )
-
+ rho1 = rho0 * exp( model_dR(i,j,k,ispec) )
+
! alpha values
dbulk = model_dA(i,j,k,ispec)
alphav1 = sqrt( alphav0**2 * exp(2.0*dbulk) + FOUR_THIRDS * betav0**2 * ( &
- exp(2.0*dbetaiso) - exp(2.0*dbulk) ) )
+ exp(2.0*dbetaiso) - exp(2.0*dbulk) ) )
alphah1 = sqrt( alphah0**2 * exp(2.0*dbulk) + FOUR_THIRDS * betah0**2 * ( &
- exp(2.0*dbetaiso) - exp(2.0*dbulk) ) )
- ! note: alphah probably not used in SPECFEM3D_GLOBE
+ exp(2.0*dbetaiso) - exp(2.0*dbulk) ) )
+ ! note: alphah probably not used in SPECFEM3D_GLOBE
! stores new model values
model_vpv_new(i,j,k,ispec) = alphav1
@@ -210,7 +210,7 @@
model_vsh_new(i,j,k,ispec) = betah1
model_eta_new(i,j,k,ispec) = eta1
model_rho_new(i,j,k,ispec) = rho1
-
+
enddo
enddo
enddo
@@ -221,7 +221,7 @@
! stores relative model perturbations
call store_perturbations()
-
+
! stop all the MPI processes, and exit
call MPI_FINALIZE(ier)
@@ -236,7 +236,7 @@
! initializes arrays
use model_update_iso
- implicit none
+ implicit none
! initialize the MPI communicator and start the NPROCTOT MPI processes
call MPI_INIT(ier)
@@ -245,7 +245,7 @@
if( sizeprocs /= nchunks_val*nproc_xi_val*nproc_eta_val ) then
print*,'sizeprocs:',sizeprocs,nchunks_val,nproc_xi_val,nproc_eta_val
- call exit_mpi('error number sizeprocs')
+ call exit_mpi(myrank,'error number sizeprocs')
endif
! model
@@ -286,7 +286,7 @@
use model_update_iso
implicit none
character(len=150) :: s_step_fac
-
+
! subjective step length to multiply to the gradient
!step_fac = 0.03
@@ -333,7 +333,7 @@
subroutine read_model()
-! reads in current transverse isotropic model: vpv.. & vsv.. & eta & rho
+! reads in current transverse isotropic model: vpv.. & vsv.. & eta & rho
use model_update_iso
implicit none
@@ -343,7 +343,7 @@
open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) then
print*,'error opening: ',trim(m_file)
- call exit_mpi('file not found')
+ call exit_mpi(myrank,'file not found')
endif
read(12) model_vpv(:,:,:,1:nspec)
close(12)
@@ -353,7 +353,7 @@
open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) then
print*,'error opening: ',trim(m_file)
- call exit_mpi('file not found')
+ call exit_mpi(myrank,'file not found')
endif
read(12) model_vph(:,:,:,1:nspec)
close(12)
@@ -363,7 +363,7 @@
open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) then
print*,'error opening: ',trim(m_file)
- call exit_mpi('file not found')
+ call exit_mpi(myrank,'file not found')
endif
read(12) model_vsv(:,:,:,1:nspec)
close(12)
@@ -373,7 +373,7 @@
open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) then
print*,'error opening: ',trim(m_file)
- call exit_mpi('file not found')
+ call exit_mpi(myrank,'file not found')
endif
read(12) model_vsh(:,:,:,1:nspec)
close(12)
@@ -383,7 +383,7 @@
open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) then
print*,'error opening: ',trim(m_file)
- call exit_mpi('file not found')
+ call exit_mpi(myrank,'file not found')
endif
read(12) model_eta(:,:,:,1:nspec)
close(12)
@@ -393,7 +393,7 @@
open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) then
print*,'error opening: ',trim(m_file)
- call exit_mpi('file not found')
+ call exit_mpi(myrank,'file not found')
endif
read(12) model_rho(:,:,:,1:nspec)
close(12)
@@ -453,7 +453,7 @@
open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) then
print*,'error opening: ',trim(m_file)
- call exit_mpi('file not found')
+ call exit_mpi(myrank,'file not found')
endif
read(12) kernel_a(:,:,:,1:nspec)
close(12)
@@ -470,29 +470,29 @@
open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) then
print*,'error opening: ',trim(m_file)
- call exit_mpi('file not found')
+ call exit_mpi(myrank,'file not found')
endif
read(12) kernel_b(:,:,:,1:nspec)
close(12)
! rho kernel
if( USE_RHO_SCALING ) then
-
+
! uses scaling relation with shear perturbations
kernel_rho(:,:,:,:) = RHO_SCALING * kernel_b(:,:,:,:)
-
+
else
-
- ! uses rho kernel
+
+ ! uses rho kernel
write(m_file,'(a,i6.6,a)') 'INPUT_GRADIENT/proc',myrank,'_reg1_rho_kernel_smooth.bin'
open(12,file=trim(m_file),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) then
print*,'error opening: ',trim(m_file)
- call exit_mpi('file not found')
- endif
+ call exit_mpi(myrank,'file not found')
+ endif
read(12) kernel_rho(:,:,:,1:nspec)
close(12)
- endif
+ endif
! statistics
@@ -553,7 +553,7 @@
open(11,file=trim(m_file),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) then
print*,'error opening: ',trim(m_file)
- call exit_mpi('file not found')
+ call exit_mpi(myrank,'file not found')
endif
read(11) x(1:nglob)
read(11) y(1:nglob)
@@ -567,7 +567,7 @@
open(11,file=trim(m_file),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) then
print*,'error opening: ',trim(m_file)
- call exit_mpi('file not found')
+ call exit_mpi(myrank,'file not found')
endif
read(11) xix
read(11) xiy
@@ -615,6 +615,9 @@
kernel_integral_a = 0._CUSTOM_REAL
kernel_integral_b = 0._CUSTOM_REAL
kernel_integral_rho = 0._CUSTOM_REAL
+ norm_a = 0._CUSTOM_REAL
+ norm_beta = 0._CUSTOM_REAL
+ norm_rho = 0._CUSTOM_REAL
do ispec = 1, NSPEC
do k = 1, NGLLZ
do j = 1, NGLLY
@@ -626,11 +629,12 @@
print*,'ibool:',ispec
print*,ibool(:,:,:,ispec)
print*
- call exit_MPI('error ibool')
+ call exit_MPI(myrank,'error ibool')
endif
! volume associated with GLL point
volumel = jacobian(i,j,k,ispec)*wgll_cube(i,j,k)
+ volume_glob = volume_glob + volumel
! kernel integration: for each element
kernel_integral_a = kernel_integral_a &
@@ -657,6 +661,7 @@
call mpi_reduce(kernel_integral_a,integral_a_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
call mpi_reduce(kernel_integral_b,integral_b_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
call mpi_reduce(kernel_integral_rho,integral_rho_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
+ call mpi_reduce(volume_glob,volume_glob_sum,1,CUSTOM_MPI_TYPE,MPI_SUM,0,MPI_COMM_WORLD,ier)
if( myrank == 0 ) then
print*,'integral kernels:'
@@ -664,6 +669,8 @@
print*,' beta : ',integral_b_sum
print*,' rho : ',integral_rho_sum
print*
+ print*,' total volume:',volume_glob_sum
+ print*
endif
! norms: for whole volume
@@ -675,7 +682,6 @@
norm_beta = sqrt(norm_beta_sum)
norm_rho = sqrt(norm_rho_sum)
-
if( myrank == 0 ) then
print*,'norm kernels:'
print*,' a : ',norm_a
@@ -737,7 +743,7 @@
print*,' rho min/max : ',min_rho,max_rho
print*
endif
-
+
! determines step length based on maximum gradient value (either shear or bulk)
if( myrank == 0 ) then
! maximum gradient values
@@ -809,12 +815,12 @@
subroutine store_kernel_updates()
-! file output for new model
+! file output for new model
use model_update_iso
- implicit none
+ implicit none
- ! kernel updates
+ ! kernel updates
fname = 'dA'
write(m_file,'(a,i6.6,a)') 'INPUT_GRADIENT/proc',myrank,'_reg1_'//trim(fname)//'.bin'
open(12,file=trim(m_file),form='unformatted')
@@ -841,12 +847,12 @@
subroutine store_new_model()
-! file output for new model
+! file output for new model
use model_update_iso
- implicit none
+ implicit none
- ! vpv model
+ ! vpv model
call mpi_reduce(maxval(model_vpv_new),max_vpv,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
call mpi_reduce(minval(model_vpv_new),min_vpv,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
fname = 'vpv_new'
@@ -855,7 +861,7 @@
write(12) model_vpv_new
close(12)
- ! vph model
+ ! vph model
call mpi_reduce(maxval(model_vph_new),max_vph,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
call mpi_reduce(minval(model_vph_new),min_vph,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
fname = 'vph_new'
@@ -864,7 +870,7 @@
write(12) model_vph_new
close(12)
- ! vsv model
+ ! vsv model
call mpi_reduce(maxval(model_vsv_new),max_vsv,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
call mpi_reduce(minval(model_vsv_new),min_vsv,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
fname = 'vsv_new'
@@ -873,7 +879,7 @@
write(12) model_vsv_new
close(12)
- ! vsh model
+ ! vsh model
call mpi_reduce(maxval(model_vsh_new),max_vsh,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
call mpi_reduce(minval(model_vsh_new),min_vsh,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
fname = 'vsh_new'
@@ -882,7 +888,7 @@
write(12) model_vsh_new
close(12)
- ! eta model
+ ! eta model
call mpi_reduce(maxval(model_eta_new),max_eta,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
call mpi_reduce(minval(model_eta_new),min_eta,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
fname = 'eta_new'
@@ -891,7 +897,7 @@
write(12) model_eta_new
close(12)
- ! rho model
+ ! rho model
call mpi_reduce(maxval(model_rho_new),max_rho,1,CUSTOM_MPI_TYPE,MPI_MAX,0,MPI_COMM_WORLD,ier)
call mpi_reduce(minval(model_rho_new),min_rho,1,CUSTOM_MPI_TYPE,MPI_MIN,0,MPI_COMM_WORLD,ier)
fname = 'rho_new'
@@ -921,12 +927,12 @@
subroutine store_perturbations()
-! file output for new model
+! file output for new model
use model_update_iso
- implicit none
+ implicit none
- ! vpv relative perturbations
+ ! vpv relative perturbations
! logarithmic perturbation: log( v_new) - log( v_old) = log( v_new / v_old )
total_model = 0.0_CUSTOM_REAL
where( model_vpv /= 0.0 ) total_model = log( model_vpv_new / model_vpv)
@@ -943,7 +949,7 @@
! vph relative perturbations
total_model = 0.0_CUSTOM_REAL
- where( model_vph /= 0.0 ) total_model = log( model_vph_new / model_vph)
+ where( model_vph /= 0.0 ) total_model = log( model_vph_new / model_vph)
write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_dvphvph.bin'
open(12,file=trim(m_file),form='unformatted',action='write')
write(12) total_model
@@ -953,7 +959,7 @@
! vsv relative perturbations
total_model = 0.0_CUSTOM_REAL
- where( model_vsv /= 0.0 ) total_model = log( model_vsv_new / model_vsv)
+ where( model_vsv /= 0.0 ) total_model = log( model_vsv_new / model_vsv)
write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_dvsvvsv.bin'
open(12,file=trim(m_file),form='unformatted',action='write')
write(12) total_model
@@ -983,7 +989,7 @@
! rho relative model perturbations
total_model = 0.0_CUSTOM_REAL
- where( model_rho /= 0.0 ) total_model = log( model_rho_new / model_rho)
+ where( model_rho /= 0.0 ) total_model = log( model_rho_new / model_rho)
write(m_file,'(a,i6.6,a)') 'OUTPUT_MODEL/proc',myrank,'_reg1_drhorho.bin'
open(12,file=trim(m_file),form='unformatted',action='write')
write(12) total_model
Deleted: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/mpif.h
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/mpif.h 2011-02-16 15:53:48 UTC (rev 17884)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/model_vp_vs/src/mpif.h 2011-02-16 21:49:17 UTC (rev 17885)
@@ -1,265 +0,0 @@
-! mpif.h. Generated from /opt/mpich/myrinet/intel//include/mpif.h by configure.
-
-!
-!
-! (C) 1993 by Argonne National Laboratory and Mississipi State University.
-! All rights reserved. See COPYRIGHT in top-level directory.
-!
-!
-! user include file for MPI programs, with no dependencies
-!
-! It really is not possible to make a perfect include file that can
-! be used by both F77 and F90 compilers, but this is close. We have removed
-! continuation lines (allows free form input in F90); systems whose
-! Fortran compilers support ! instead of just C or * for comments can
-! globally replace a C in the first column with !; the resulting file
-! should work for both Fortran 77 and Fortran 90.
-!
-! If your Fortran compiler supports ! for comments, you can run this
-! through sed with
-! sed -e 's/^C/\!/g'
-!
-! We have also removed the use of contractions (involving the single quote)
-! character because some users use .F instead of .f files (to invoke the
-! cpp preprocessor) and further, their preprocessor is determined to find
-! matching single quote pairs (and probably double quotes; given the
-! different rules in C and Fortran, this sounds like a disaster). Rather than
-! take the position that the poor users should get a better system, we
-! have removed the text that caused problems. Of course, the users SHOULD
-! get a better system...
-!
-! return codes
- INTEGER MPI_SUCCESS,MPI_ERR_BUFFER,MPI_ERR_COUNT,MPI_ERR_TYPE
- INTEGER MPI_ERR_TAG,MPI_ERR_COMM,MPI_ERR_RANK,MPI_ERR_ROOT
- INTEGER MPI_ERR_GROUP
- INTEGER MPI_ERR_OP,MPI_ERR_TOPOLOGY,MPI_ERR_DIMS,MPI_ERR_ARG
- INTEGER MPI_ERR_UNKNOWN,MPI_ERR_TRUNCATE,MPI_ERR_OTHER
- INTEGER MPI_ERR_INTERN,MPI_ERR_IN_STATUS,MPI_ERR_PENDING
- INTEGER MPI_ERR_REQUEST, MPI_ERR_LASTCODE
- PARAMETER (MPI_SUCCESS=0,MPI_ERR_BUFFER=1,MPI_ERR_COUNT=2)
- PARAMETER (MPI_ERR_TYPE=3,MPI_ERR_TAG=4,MPI_ERR_COMM=5)
- PARAMETER (MPI_ERR_RANK=6,MPI_ERR_ROOT=7,MPI_ERR_GROUP=8)
- PARAMETER (MPI_ERR_OP=9,MPI_ERR_TOPOLOGY=10,MPI_ERR_DIMS=11)
- PARAMETER (MPI_ERR_ARG=12,MPI_ERR_UNKNOWN=13)
- PARAMETER (MPI_ERR_TRUNCATE=14,MPI_ERR_OTHER=15)
- PARAMETER (MPI_ERR_INTERN=16,MPI_ERR_IN_STATUS=17)
- PARAMETER (MPI_ERR_PENDING=18,MPI_ERR_REQUEST=19)
- PARAMETER (MPI_ERR_LASTCODE=1073741823)
-!
- INTEGER MPI_UNDEFINED
- parameter (MPI_UNDEFINED = (-32766))
-!
- INTEGER MPI_GRAPH, MPI_CART
- PARAMETER (MPI_GRAPH = 1, MPI_CART = 2)
- INTEGER MPI_PROC_NULL
- PARAMETER ( MPI_PROC_NULL = (-1) )
-!
- INTEGER MPI_BSEND_OVERHEAD
- PARAMETER ( MPI_BSEND_OVERHEAD = 512 )
-
- INTEGER MPI_SOURCE, MPI_TAG, MPI_ERROR
- PARAMETER(MPI_SOURCE=2, MPI_TAG=3, MPI_ERROR=4)
- INTEGER MPI_STATUS_SIZE
- PARAMETER (MPI_STATUS_SIZE=4)
- INTEGER MPI_MAX_PROCESSOR_NAME, MPI_MAX_ERROR_STRING
- PARAMETER (MPI_MAX_PROCESSOR_NAME=256)
- PARAMETER (MPI_MAX_ERROR_STRING=512)
- INTEGER MPI_MAX_NAME_STRING
- PARAMETER (MPI_MAX_NAME_STRING=63)
- INTEGER MPI_MAX_PORT_NAME
- PARAMETER (MPI_MAX_PORT_NAME=256)
-!
- INTEGER MPI_COMM_NULL
- PARAMETER (MPI_COMM_NULL=0)
-!
- INTEGER MPI_DATATYPE_NULL
- PARAMETER (MPI_DATATYPE_NULL = 0)
-
- INTEGER MPI_ERRHANDLER_NULL
- PARAMETER (MPI_ERRHANDLER_NULL = 0)
-
- INTEGER MPI_GROUP_NULL
- PARAMETER (MPI_GROUP_NULL = 0)
-
- INTEGER MPI_KEYVAL_INVALID
- PARAMETER (MPI_KEYVAL_INVALID = 0)
-
- INTEGER MPI_REQUEST_NULL
- PARAMETER (MPI_REQUEST_NULL = 0)
-!
- INTEGER MPI_IDENT, MPI_CONGRUENT, MPI_SIMILAR, MPI_UNEQUAL
- PARAMETER (MPI_IDENT=0, MPI_CONGRUENT=1, MPI_SIMILAR=2)
- PARAMETER (MPI_UNEQUAL=3)
-!
-! MPI_BOTTOM needs to be a known address; here we put it at the
-! beginning of the common block. The point-to-point and collective
-! routines know about MPI_BOTTOM, but MPI_TYPE_STRUCT as yet does not.
-!
-! MPI_STATUS_IGNORE and MPI_STATUSES_IGNORE are similar objects
-! Until the underlying MPI library implements the C version of these
-! (a null pointer), these are declared as arrays of MPI_STATUS_SIZE
-!
-! The types MPI_INTEGER1,2,4 and MPI_REAL4,8 are OPTIONAL.
-! Their values are zero if they are not available. Note that
-! using these reduces the portability of code (though may enhance
-! portability between Crays and other systems)
-!
- INTEGER MPI_TAG_UB, MPI_HOST, MPI_IO
- INTEGER MPI_BOTTOM
- INTEGER MPI_STATUS_IGNORE(MPI_STATUS_SIZE)
- INTEGER MPI_STATUSES_IGNORE(MPI_STATUS_SIZE)
- INTEGER MPI_INTEGER, MPI_REAL, MPI_DOUBLE_PRECISION
- INTEGER MPI_COMPLEX, MPI_DOUBLE_COMPLEX,MPI_LOGICAL
- INTEGER MPI_CHARACTER, MPI_BYTE, MPI_2INTEGER, MPI_2REAL
- INTEGER MPI_2DOUBLE_PRECISION, MPI_2COMPLEX, MPI_2DOUBLE_COMPLEX
- INTEGER MPI_UB, MPI_LB
- INTEGER MPI_PACKED, MPI_WTIME_IS_GLOBAL
- INTEGER MPI_COMM_WORLD, MPI_COMM_SELF, MPI_GROUP_EMPTY
- INTEGER MPI_SUM, MPI_MAX, MPI_MIN, MPI_PROD, MPI_LAND, MPI_BAND
- INTEGER MPI_LOR, MPI_BOR, MPI_LXOR, MPI_BXOR, MPI_MINLOC
- INTEGER MPI_MAXLOC
- INTEGER MPI_OP_NULL
- INTEGER MPI_ERRORS_ARE_FATAL, MPI_ERRORS_RETURN
-!
- PARAMETER (MPI_ERRORS_ARE_FATAL=119)
- PARAMETER (MPI_ERRORS_RETURN=120)
-!
- PARAMETER (MPI_COMPLEX=23,MPI_DOUBLE_COMPLEX=24,MPI_LOGICAL=25)
- PARAMETER (MPI_REAL=26,MPI_DOUBLE_PRECISION=27,MPI_INTEGER=28)
- PARAMETER (MPI_2INTEGER=29,MPI_2COMPLEX=30,MPI_2DOUBLE_COMPLEX=31)
- PARAMETER (MPI_2REAL=32,MPI_2DOUBLE_PRECISION=33,MPI_CHARACTER=1)
- PARAMETER (MPI_BYTE=3,MPI_UB=16,MPI_LB=15,MPI_PACKED=14)
-
- INTEGER MPI_ORDER_C, MPI_ORDER_FORTRAN
- PARAMETER (MPI_ORDER_C=56, MPI_ORDER_FORTRAN=57)
- INTEGER MPI_DISTRIBUTE_BLOCK, MPI_DISTRIBUTE_CYCLIC
- INTEGER MPI_DISTRIBUTE_NONE, MPI_DISTRIBUTE_DFLT_DARG
- PARAMETER (MPI_DISTRIBUTE_BLOCK=121, MPI_DISTRIBUTE_CYCLIC=122)
- PARAMETER (MPI_DISTRIBUTE_NONE=123)
- PARAMETER (MPI_DISTRIBUTE_DFLT_DARG=-49767)
- INTEGER MPI_MAX_INFO_KEY, MPI_MAX_INFO_VAL
- PARAMETER (MPI_MAX_INFO_KEY=255, MPI_MAX_INFO_VAL=1024)
- INTEGER MPI_INFO_NULL
- PARAMETER (MPI_INFO_NULL=0)
-
-!
-! Optional Fortran Types. Configure attempts to determine these.
-!
- INTEGER MPI_INTEGER1, MPI_INTEGER2, MPI_INTEGER4, MPI_INTEGER8
- INTEGER MPI_INTEGER16
- INTEGER MPI_REAL4, MPI_REAL8, MPI_REAL16
- INTEGER MPI_COMPLEX8, MPI_COMPLEX16, MPI_COMPLEX32
- PARAMETER (MPI_INTEGER1=1,MPI_INTEGER2=4)
- PARAMETER (MPI_INTEGER4=6)
- PARAMETER (MPI_INTEGER8=8)
- PARAMETER (MPI_INTEGER16=0)
- PARAMETER (MPI_REAL4=10)
- PARAMETER (MPI_REAL8=11)
- PARAMETER (MPI_REAL16=12)
- PARAMETER (MPI_COMPLEX8=23)
- PARAMETER (MPI_COMPLEX16=24)
- PARAMETER (MPI_COMPLEX32=0)
-!
-! This is now handled with either the "pointer" extension or this same
-! code, appended at the end.
-! COMMON /MPIPRIV/ MPI_BOTTOM,MPI_STATUS_IGNORE,MPI_STATUSES_IGNORE
-!C
-!C Without this save, some Fortran implementations may make the common
-!C dynamic!
-!C
-!C For a Fortran90 module, we might replace /MPIPRIV/ with a simple
-!C SAVE MPI_BOTTOM
-!C
-! SAVE /MPIPRIV/
-!
- PARAMETER (MPI_MAX=100,MPI_MIN=101,MPI_SUM=102,MPI_PROD=103)
- PARAMETER (MPI_LAND=104,MPI_BAND=105,MPI_LOR=106,MPI_BOR=107)
- PARAMETER (MPI_LXOR=108,MPI_BXOR=109,MPI_MINLOC=110)
- PARAMETER (MPI_MAXLOC=111, MPI_OP_NULL=0)
-!
- PARAMETER (MPI_GROUP_EMPTY=90,MPI_COMM_WORLD=91,MPI_COMM_SELF=92)
- PARAMETER (MPI_TAG_UB=80,MPI_HOST=82,MPI_IO=84)
- PARAMETER (MPI_WTIME_IS_GLOBAL=86)
-!
- INTEGER MPI_ANY_SOURCE
- PARAMETER (MPI_ANY_SOURCE = (-2))
- INTEGER MPI_ANY_TAG
- PARAMETER (MPI_ANY_TAG = (-1))
-!
- INTEGER MPI_VERSION, MPI_SUBVERSION
- PARAMETER (MPI_VERSION = 1, MPI_SUBVERSION = 2)
-!
-! There are additional MPI-2 constants
- INTEGER MPI_ADDRESS_KIND, MPI_OFFSET_KIND
- PARAMETER (MPI_ADDRESS_KIND=8)
- PARAMETER (MPI_OFFSET_KIND=8)
-!
-! All other MPI routines are subroutines
-! This may cause some Fortran compilers to complain about defined and
-! not used. Such compilers should be improved.
-!
-! Some Fortran compilers will not link programs that contain
-! external statements to routines that are not provided, even if
-! the routine is never called. Remove PMPI_WTIME and PMPI_WTICK
-! if you have trouble with them.
-!
- DOUBLE PRECISION MPI_WTIME, MPI_WTICK,PMPI_WTIME,PMPI_WTICK
- EXTERNAL MPI_WTIME, MPI_WTICK,PMPI_WTIME,PMPI_WTICK
-!
-! The attribute copy/delete subroutines are symbols that can be passed
-! to MPI routines
-!
- EXTERNAL MPI_NULL_COPY_FN, MPI_NULL_DELETE_FN, MPI_DUP_FN
- COMMON /MPIPRIV/ MPI_BOTTOM,MPI_STATUS_IGNORE,MPI_STATUSES_IGNORE
-!
-! Without this save, some Fortran implementations may make the common
-! dynamic!
-!
-! For a Fortran90 module, we might replace /MPIPRIV/ with a simple
-! SAVE MPI_BOTTOM
-!
- SAVE /MPIPRIV/
-!
-! $Id: mpiof.h.in,v 1.3 1999/08/06 18:33:09 thakur Exp $
-!
-! Copyright (C) 1997 University of Chicago.
-! See COPYRIGHT notice in top-level directory.
-!
-!
-! user include file for Fortran MPI-IO programs
-!
- INTEGER MPI_MODE_RDONLY, MPI_MODE_RDWR, MPI_MODE_WRONLY
- INTEGER MPI_MODE_DELETE_ON_CLOSE, MPI_MODE_UNIQUE_OPEN
- INTEGER MPI_MODE_CREATE, MPI_MODE_EXCL
- INTEGER MPI_MODE_APPEND, MPI_MODE_SEQUENTIAL
- PARAMETER (MPI_MODE_RDONLY=2, MPI_MODE_RDWR=8, MPI_MODE_WRONLY=4)
- PARAMETER (MPI_MODE_CREATE=1, MPI_MODE_DELETE_ON_CLOSE=16)
- PARAMETER (MPI_MODE_UNIQUE_OPEN=32, MPI_MODE_EXCL=64)
- PARAMETER (MPI_MODE_APPEND=128, MPI_MODE_SEQUENTIAL=256)
-!
- INTEGER MPI_FILE_NULL
- PARAMETER (MPI_FILE_NULL=0)
-!
- INTEGER MPI_MAX_DATAREP_STRING
- PARAMETER (MPI_MAX_DATAREP_STRING=128)
-!
- INTEGER MPI_SEEK_SET, MPI_SEEK_CUR, MPI_SEEK_END
- PARAMETER (MPI_SEEK_SET=600, MPI_SEEK_CUR=602, MPI_SEEK_END=604)
-!
- INTEGER MPIO_REQUEST_NULL
- PARAMETER (MPIO_REQUEST_NULL=0)
-!
-!
-!
-
-
-
-
-
-
-
-!
-!
-!
-!
-!
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/src/smooth_sem_globe.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/src/smooth_sem_globe.f90 2011-02-16 15:53:48 UTC (rev 17884)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/pangu/smooth/src/smooth_sem_globe.f90 2011-02-16 21:49:17 UTC (rev 17885)
@@ -3,33 +3,33 @@
! this program can be used for smoothing a (summed) event kernel,
! where it smooths files with a given input kernel name:
!
-! Usage:
+! Usage:
! ./smooth_sem_globe sigma_h(km) sigma_v(km) kernel_file_name scratch_file_dir scratch_topo_dir
! e.g.
-! ./smooth_sem_globe 160 10 bulk_c_kernel OUTPUT_SUM/ topo/
-!
+! ./smooth_sem_globe 160 10 bulk_c_kernel OUTPUT_SUM/ topo/
+!
! where:
! sigma_h - gaussian width for horizontal smoothing (in km)
! sigma_v - gaussian width for vertical smoothing (in km)
-! kernel_file_name - takes file with this kernel name,
+! kernel_file_name - takes file with this kernel name,
! e.g. "bulk_c_kernel"
-! scratch_file_dir - directory containing kernel files,
+! scratch_file_dir - directory containing kernel files,
! e.g. proc***_reg1_bulk_c_kernel.bin
-! topo_dir - directory containing mesh topo files:
+! topo_dir - directory containing mesh topo files:
! proc***_solver_data1.bin, proc***_solver_data2.bin
-! outputs:
+! outputs:
! puts the resulting, smoothed kernel files into the same directory as scratch_file_dir/
! with a file ending "proc***_kernel_smooth.bin"
program smooth_sem_globe
-! this is the embarassingly-parallel program that smooths any specfem function (primarily the kernels)
+! this is the embarassingly-parallel program that smooths any specfem function (primarily the kernels)
! that has the dimension of (NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT)
!
! notice that it uses the constants_globe.h and precision_globe.h files
-! from the original SPECFEM3D_GLOBE package, and the
+! from the original SPECFEM3D_GLOBE package, and the
! values_from_mesher_globe.h file from the output of the mesher (or create_header_file),
-! therefore, you need to compile it for your specific case
+! therefore, you need to compile it for your specific case
!
! NOTE: smoothing can be different in radial & horizontal directions; mesh is in spherical geometry.
! algorithm uses vector components in radial/horizontal direction
@@ -42,13 +42,13 @@
! ======================================================
! USER PARAMETERS
-
+
! copy from Par_file values
integer, parameter :: NPROC_XI = 12
integer, parameter :: NPROC_ETA = 12
integer, parameter :: NCHUNKS = 1
- ! taken from values_from_mesher.h:
+ ! taken from values_from_mesher.h:
! average size of a spectral element in km = ...
! e.g. nproc 12x12, nex 192: element_size = 52.122262
real(kind=CUSTOM_REAL),parameter:: element_size = 52.12262
@@ -58,7 +58,7 @@
!takes region 1 kernels
integer, parameter :: NSPEC_MAX = NSPEC_CRUST_MANTLE_ADJOINT
integer, parameter :: NGLOB_MAX = NGLOB_CRUST_MANTLE
-
+
! only include the neighboring 3 x 3 slices
integer, parameter :: NSLICES = 3
integer ,parameter :: NSLICES2 = NSLICES * NSLICES
@@ -76,7 +76,7 @@
character(len=256) :: ks_file, reg_name
character(len=256), dimension(NSLICES2) :: x_file, y_file, z_file, j_file, k_file, i_file
character(len=256), dimension(NSLICES2) :: solver1_file,solver2_file
-
+
logical :: global_code
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_MAX) :: ibool
@@ -86,7 +86,7 @@
real(kind=CUSTOM_REAL) xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_MAX) :: &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
-
+
real(kind=CUSTOM_REAL), dimension(NGLOB_MAX) :: x, y, z
real(kind=CUSTOM_REAL), dimension(NSPEC_MAX) :: cx0, cy0, cz0, cx, cy, cz
@@ -100,9 +100,9 @@
! array with all the weights in the cube
double precision, dimension(NGLLX,NGLLY,NGLLZ) :: wgll_cube
- real(kind=CUSTOM_REAL) :: dist_h,dist_v
+ real(kind=CUSTOM_REAL) :: dist_h,dist_v
real(kind=CUSTOM_REAL) :: r1,theta1
-
+
! ============ program starts here =====================
! initialize the MPI communicator and start the NPROCTOT MPI processes
@@ -112,7 +112,7 @@
if (myrank == 0) print*,"smooth:"
call mpi_barrier(MPI_COMM_WORLD,ier)
-
+
! arguments
call getarg(1,s_sigma_h)
call getarg(2,s_sigma_v)
@@ -144,12 +144,12 @@
! user output
if (myrank == 0) then
print*,"defaults:"
- print*," NPROC_XI , NPROC_ETA: ",NPROC_XI,NPROC_ETA
- print*," NCHUNKS : ",NCHUNKS
- print*," element size on surface(km): ",element_size
+ print*," NPROC_XI , NPROC_ETA: ",NPROC_XI,NPROC_ETA
+ print*," NCHUNKS : ",NCHUNKS
+ print*," element size on surface(km): ",element_size
print*," smoothing sigma_h , sigma_v: ",sigma_h,sigma_v
endif
- ! synchronizes
+ ! synchronizes
call mpi_barrier(MPI_COMM_WORLD,ier)
! initializes lengths
@@ -157,18 +157,18 @@
if (global_code) element_size_m = element_size_m/R_EARTH
sigma_h = sigma_h * 1000.0 ! m
- if (global_code) sigma_h = sigma_h / R_EARTH ! scale
+ if (global_code) sigma_h = sigma_h / R_EARTH ! scale
sigma_v = sigma_v * 1000.0 ! m
if (global_code) sigma_v = sigma_v / R_EARTH ! scale
-
+
sigma_h2 = 2.0 * sigma_h ** 2 ! factor two for gaussian distribution with standard variance sigma
sigma_v2 = 2.0 * sigma_v ** 2
! search radius
- sigma_h3 = 3.0 * sigma_h + element_size_m
- sigma_v3 = 3.0 * sigma_v + element_size_m
+ sigma_h3 = 3.0 * sigma_h + element_size_m
+ sigma_v3 = 3.0 * sigma_v + element_size_m
- ! theoretic normal value
+ ! theoretic normal value
! (see integral over -inf to +inf of exp[- x*x/(2*sigma) ] = sigma * sqrt(2*pi) )
! note: smoothing is using a gaussian (ellipsoid for sigma_h /= sigma_v),
@@ -213,7 +213,7 @@
islice(j) = islice0(i)
endif
enddo
- nums = j
+ nums = j
if( myrank == 0 ) then
print *,'slices:',nums
@@ -238,17 +238,17 @@
! point locations
open(11,file=solver2_file(1),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error opening solver2 file')
-
+
read(11) x(1:nglob(1))
read(11) y(1:nglob(1))
read(11) z(1:nglob(1))
read(11) ibool(:,:,:,1:nspec(1))
close(11)
-
+
! jacobian
open(11,file=solver1_file(1),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error opening solver1 file')
-
+
read(11) xix
read(11) xiy
read(11) xiz
@@ -257,7 +257,7 @@
read(11) etaz
read(11) gammax
read(11) gammay
- read(11) gammaz
+ read(11) gammaz
close(11)
! get the location of the center of the elements
@@ -269,8 +269,8 @@
xl(i,j,k,ispec) = x(iglob)
yl(i,j,k,ispec) = y(iglob)
zl(i,j,k,ispec) = z(iglob)
-
- ! build jacobian
+
+ ! build jacobian
! get derivatives of ux, uy and uz with respect to x, y and z
xixl = xix(i,j,k,ispec)
xiyl = xiy(i,j,k,ispec)
@@ -287,9 +287,9 @@
+ xizl*(etaxl*gammayl-etayl*gammaxl))
jacobian(i,j,k,ispec) = jacobianl
-
-
-
+
+
+
enddo
enddo
enddo
@@ -310,7 +310,7 @@
! loop over all the slices
do iproc = 1, nums
-
+
! read in the topology, kernel files, calculate center of elements
! point locations
! given in cartesian coordinates
@@ -325,7 +325,7 @@
open(11,file=solver1_file(iproc),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error opening slices: solver1 file')
-
+
read(11) xix
read(11) xiy
read(11) xiz
@@ -334,7 +334,7 @@
read(11) etaz
read(11) gammax
read(11) gammay
- read(11) gammaz
+ read(11) gammaz
close(11)
! get the location of the center of the elements
@@ -342,7 +342,7 @@
do k = 1, NGLLZ
do j = 1, NGLLY
do i = 1, NGLLX
- ! build jacobian
+ ! build jacobian
! get derivatives of ux, uy and uz with respect to x, y and z
xixl = xix(i,j,k,ispec)
xiyl = xiy(i,j,k,ispec)
@@ -357,7 +357,7 @@
jacobianl = 1._CUSTOM_REAL / (xixl*(etayl*gammazl-etazl*gammayl) &
- xiyl*(etaxl*gammazl-etazl*gammaxl) &
+ xizl*(etaxl*gammayl-etayl*gammaxl))
- jacobian(i,j,k,ispec) = jacobianl
+ jacobian(i,j,k,ispec) = jacobianl
enddo
enddo
enddo
@@ -366,7 +366,7 @@
! kernel file
open(11,file=k_file(iproc),status='old',form='unformatted',iostat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error opening kernel file')
-
+
read(11) kernel(:,:,:,1:nspec(iproc))
close(11)
@@ -385,7 +385,7 @@
iglob = ibool(i,j,k,ispec2)
xx(i,j,k,ispec2) = x(iglob)
yy(i,j,k,ispec2) = y(iglob)
- zz(i,j,k,ispec2) = z(iglob)
+ zz(i,j,k,ispec2) = z(iglob)
enddo
enddo
enddo
@@ -395,50 +395,64 @@
enddo
! loop over elements to be smoothed in the current slice
- do ispec = 1, nspec(1)
+ do ispec = 1, nspec(1)
! --- only double loop over the elements in the search radius ---
do ispec2 = 1, nspec(iproc)
-
+
! calculates horizontal and vertical distance between two element centers
-
+
! vector approximation
call get_distance_vec(dist_h,dist_v,cx0(ispec),cy0(ispec),cz0(ispec),&
cx(ispec2),cy(ispec2),cz(ispec2))
-
+
! note: distances and sigmah, sigmav are normalized by R_EARTH
-
- ! checks distance between centers of elements
+
+ ! checks distance between centers of elements
if ( dist_h > sigma_h3 .or. abs(dist_v) > sigma_v3 ) cycle
- ! integration factors
- factor(:,:,:) = jacobian(:,:,:,ispec2) * wgll_cube(:,:,:)
+ ! integration factors:
+ ! uses volume assigned to GLL points
+ factor(:,:,:) = jacobian(:,:,:,ispec2) * wgll_cube(:,:,:)
+ ! no volume
+ !factor(:,:,:) = 1.0_CUSTOM_REAL
-
! loop over GLL points of the elements in current slice (ispec)
- do k = 1, NGLLZ
+ do k = 1, NGLLZ
do j = 1, NGLLY
do i = 1, NGLLX
-
+
+ ! reference location
! current point (i,j,k,ispec) location, cartesian coordinates
- x0 = xl(i,j,k,ispec)
- y0 = yl(i,j,k,ispec)
- z0 = zl(i,j,k,ispec)
+ x0 = xl(i,j,k,ispec)
+ y0 = yl(i,j,k,ispec)
+ z0 = zl(i,j,k,ispec)
! calculate weights based on gaussian smoothing
call smoothing_weights_vec(x0,y0,z0,ispec2,sigma_h2,sigma_v2,exp_val,&
xx(:,:,:,ispec2),yy(:,:,:,ispec2),zz(:,:,:,ispec2))
-
+
! adds GLL integration weights
exp_val(:,:,:) = exp_val(:,:,:) * factor(:,:,:)
! adds contribution of element ispec2 to smoothed kernel values
tk(i,j,k,ispec) = tk(i,j,k,ispec) + sum(exp_val(:,:,:) * kernel(:,:,:,ispec2))
-
+
! normalization, integrated values of gaussian smoothing function
bk(i,j,k,ispec) = bk(i,j,k,ispec) + sum(exp_val(:,:,:))
- enddo
+ ! checks number
+ !if( isNaN(tk(i,j,k,ispec)) ) then
+ ! print*,'error tk NaN: ',tk(i,j,k,ispec)
+ ! print*,'rank:',myrank
+ ! print*,'i,j,k,ispec:',i,j,k,ispec
+ ! print*,'tk: ',tk(i,j,k,ispec),'bk:',bk(i,j,k,ispec)
+ ! print*,'sum exp_val: ',sum(exp_val(:,:,:)),'sum factor:',sum(factor(:,:,:))
+ ! print*,'sum kernel:',sum(kernel(:,:,:,ispec2))
+ ! call exit_MPI('error NaN')
+ !endif
+
+ enddo
enddo
enddo ! (i,j,k)
enddo ! (ispec2)
@@ -454,30 +468,39 @@
do i = 1, NGLLX
! checks the normalization criterion
- ! e.g. sigma_h 160km, sigma_v 40km:
+ ! e.g. sigma_h 160km, sigma_v 40km:
! norm (not squared sigma_h ) ~ 0.001
! norm ( squared sigma_h) ~ 6.23 * e-5
- if (abs(bk(i,j,k,ispec) - norm) > 1.e-4 ) then
+ if (abs(bk(i,j,k,ispec) - norm) > 1.e-4 ) then
print *, 'Problem norm here --- ', myrank, ispec, i, j, k, bk(i,j,k,ispec), norm
!call exit_mpi(myrank, 'Error computing Gaussian function on the grid')
endif
-
- ! normalizes smoothed kernel values by integral value of gaussian weighting
- kernel_smooth(i,j,k,ispec) = tk(i,j,k,ispec) / bk(i,j,k,ispec)
-
+
+ ! normalizes smoothed kernel values by integral value of gaussian weighting
+ kernel_smooth(i,j,k,ispec) = tk(i,j,k,ispec) / bk(i,j,k,ispec)
+
+ ! checks number
+ if( isNaN(kernel_smooth(i,j,k,ispec)) ) then
+ print*,'error kernel_smooth NaN: ',kernel_smooth(i,j,k,ispec)
+ print*,'rank:',myrank
+ print*,'i,j,k,ispec:',i,j,k,ispec
+ print*,'tk: ',tk(i,j,k,ispec),'bk:',bk(i,j,k,ispec)
+ call exit_MPI('error NaN')
+ endif
+
enddo
enddo
enddo
enddo
if (myrank == 0) write(*,*) ' norm: ',norm
-
+
! file output
open(11,file=trim(ks_file),status='unknown',form='unformatted',iostat=ier)
if( ier /= 0 ) call exit_mpi(myrank,'error opening smoothed kernel file')
-
+
! Note: output the following instead of kernel_smooth(:,:,:,1:nspec(1)) to create files of the same sizes
write(11) kernel_smooth(:,:,:,:)
- close(11)
+ close(11)
if (myrank == 0) print *,' written:',trim(ks_file)
@@ -506,9 +529,9 @@
implicit none
include "constants_globe.h"
-
+
real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ),intent(out) :: exp_val
- real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ),intent(in) :: xx_elem, yy_elem, zz_elem
+ real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ),intent(in) :: xx_elem, yy_elem, zz_elem
real(kind=CUSTOM_REAL),intent(in) :: x0,y0,z0,sigma_h2,sigma_v2
integer,intent(in) :: ispec2
@@ -528,24 +551,34 @@
! -(yy(:,:,:,ispec2)-y0)**2/(sigma_h2) &
! -(zz(:,:,:,ispec2)-z0)**2/(sigma_v2) ) * factor(:,:,:)
! >>>>>
-
- do kk = 1, NGLLZ
+
+ do kk = 1, NGLLZ
do jj = 1, NGLLY
do ii = 1, NGLLX
- ! point in second slice
-
+ ! point in second slice
+
! vector approximation:
call get_distance_vec(dist_h,dist_v,x0,y0,z0, &
xx_elem(ii,jj,kk),yy_elem(ii,jj,kk),zz_elem(ii,jj,kk))
-
+
! gaussian function
- exp_val(ii,jj,kk) = exp( - dist_h*dist_h/sigma_h2 &
- - dist_v*dist_v/sigma_v2 ) ! * factor(ii,jj,kk)
+ exp_val(ii,jj,kk) = exp( - (dist_h*dist_h)/sigma_h2 &
+ - (dist_v*dist_v)/sigma_v2 ) ! * factor(ii,jj,kk)
+
+ ! checks number
+ !if( isNaN(exp_val(ii,jj,kk)) ) then
+ ! print*,'error exp_val NaN: ',exp_val(ii,jj,kk)
+ ! print*,'i,j,k:',ii,jj,kk
+ ! print*,'dist_h: ',dist_h,'dist_v:',dist_v
+ ! print*,'sigma_h2:',sigma_h2,'sigma_v2:',sigma_v2
+ ! call exit_MPI('error NaN')
+ !endif
+
enddo
enddo
enddo
-
+
end subroutine smoothing_weights_vec
@@ -555,41 +588,59 @@
subroutine get_distance_vec(dist_h,dist_v,x0,y0,z0,x1,y1,z1)
-! returns vector lengths as distances in radial and horizontal direction
+! returns vector lengths as distances in radial and horizontal direction
implicit none
include "constants_globe.h"
-
+
real(kind=CUSTOM_REAL),intent(out) :: dist_h,dist_v
real(kind=CUSTOM_REAL),intent(in) :: x0,y0,z0,x1,y1,z1
-
+
! local parameters
- real(kind=CUSTOM_REAL) :: r0,r1,alpha
- real(kind=CUSTOM_REAL) :: vx,vy,vz
-
- ! vertical distance
+ real(kind=CUSTOM_REAL) :: r0,r1
+ real(kind=CUSTOM_REAL) :: theta,ratio
+ !real(kind=CUSTOM_REAL) :: vx,vy,vz,alpha
+
+ ! vertical distance
r0 = sqrt( x0*x0 + y0*y0 + z0*z0 ) ! length of first position vector
- r1 = sqrt( x1*x1 + y1*y1 + z1*z1 )
- dist_v = r1 - r0
+ r1 = sqrt( x1*x1 + y1*y1 + z1*z1 )
+ dist_v = r1 - r0
! only for flat earth with z in depth: dist_v = sqrt( (cz(ispec2)-cz0(ispec))** 2)
-
- ! horizontal distance
- ! length of vector from point 0 to point 1
+
+ ! epicentral distance
+ ! (accounting for spherical curvature)
+ ! calculates distance of circular segment
+ ! angle between r0 and r1 in radian
+ ! given by dot-product of two vectors
+ ratio = (x0*x1 + y0*y1 + z0*z1)/(r0 * r1)
+
+ ! checks boundaries of ratio (due to numerical inaccuracies)
+ if( ratio > 1.0_CUSTOM_REAL ) ratio = 1.0_CUSTOM_REAL
+ if( ratio < -1.0_CUSTOM_REAL ) ratio = -1.0_CUSTOM_REAL
+
+ theta = acos( ratio )
+
+ ! segment length at heigth of r1
+ dist_h = r1 * theta
+
+ ! vector approximation (fast computation): neglects curvature
+ ! horizontal distance
+ ! length of vector from point 0 to point 1
! assuming small earth curvature (since only for neighboring elements)
-
+
! scales r0 to have same length as r1
- alpha = r1 / r0
- vx = alpha * x0
- vy = alpha * y0
- vz = alpha * z0
-
+ !alpha = r1 / r0
+ !vx = alpha * x0
+ !vy = alpha * y0
+ !vz = alpha * z0
+
! vector in horizontal between new r0 and r1
- vx = x1 - vx
- vy = y1 - vy
- vz = z1 - vz
-
- ! distance is vector length
- dist_h = sqrt( vx*vx + vy*vy + vz*vz )
-
+ !vx = x1 - vx
+ !vy = y1 - vy
+ !vz = z1 - vz
+
+ ! distance is vector length
+ !dist_h = sqrt( vx*vx + vy*vy + vz*vz )
+
end subroutine get_distance_vec
More information about the CIG-COMMITS
mailing list