[cig-commits] r21700 - in seismo/3D/SPECFEM3D/branches/full_aniso_kernels_matharu/src: shared specfem3D
gian at geodynamics.org
gian at geodynamics.org
Tue Apr 2 17:54:30 PDT 2013
Author: gian
Date: 2013-04-02 17:54:29 -0700 (Tue, 02 Apr 2013)
New Revision: 21700
Modified:
seismo/3D/SPECFEM3D/branches/full_aniso_kernels_matharu/src/shared/constants.h.in
seismo/3D/SPECFEM3D/branches/full_aniso_kernels_matharu/src/specfem3D/compute_kernels.f90
seismo/3D/SPECFEM3D/branches/full_aniso_kernels_matharu/src/specfem3D/prepare_timerun.F90
seismo/3D/SPECFEM3D/branches/full_aniso_kernels_matharu/src/specfem3D/read_mesh_databases.f90
seismo/3D/SPECFEM3D/branches/full_aniso_kernels_matharu/src/specfem3D/save_adjoint_kernels.f90
seismo/3D/SPECFEM3D/branches/full_aniso_kernels_matharu/src/specfem3D/specfem3D_par.f90
Log:
Added implementation of anisotropic kernels to custom branch
Modified: seismo/3D/SPECFEM3D/branches/full_aniso_kernels_matharu/src/shared/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D/branches/full_aniso_kernels_matharu/src/shared/constants.h.in 2013-04-02 22:52:02 UTC (rev 21699)
+++ seismo/3D/SPECFEM3D/branches/full_aniso_kernels_matharu/src/shared/constants.h.in 2013-04-03 00:54:29 UTC (rev 21700)
@@ -259,6 +259,21 @@
! postprocess the colors to balance them if Droux (1993) algorithm is not used
logical, parameter :: BALANCE_COLORS_SIMPLE_ALGO = .false.
+!!-----------------------------------------------------------
+!!
+!! adjoint kernel outputs
+!!
+!!-----------------------------------------------------------
+
+! this parameter must be set to .true. to compute anisotropic kernels
+! in crust and mantle (related to the 21 Cij in geographical coordinates)
+! default is .false. to compute isotropic kernels (related to alpha and beta)
+ logical, parameter :: ANISOTROPIC_KL = .false.
+
+! compute transverse isotropic kernels (alpha_v,alpha_h,beta_v,beta_h,eta,rho)
+! rather than fully anisotropic kernels in case ANISOTROPIC_KL is set to .true.
+ logical, parameter :: SAVE_TRANSVERSE_KL = .false.
+
!------------------------------------------------------
!----------- do not modify anything below -------------
!------------------------------------------------------
Modified: seismo/3D/SPECFEM3D/branches/full_aniso_kernels_matharu/src/specfem3D/compute_kernels.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/full_aniso_kernels_matharu/src/specfem3D/compute_kernels.f90 2013-04-02 22:52:02 UTC (rev 21699)
+++ seismo/3D/SPECFEM3D/branches/full_aniso_kernels_matharu/src/specfem3D/compute_kernels.f90 2013-04-03 00:54:29 UTC (rev 21700)
@@ -74,6 +74,7 @@
implicit none
! local parameters
integer :: i,j,k,ispec,iglob
+ real(kind=CUSTOM_REAL),dimension(21) :: prod !, cijkl_kl_local
real(kind=CUSTOM_REAL), dimension(5) :: epsilondev_loc,b_epsilondev_loc
if( .not. GPU_MODE ) then
@@ -87,22 +88,7 @@
do j = 1, NGLLY
do i = 1, NGLLX
iglob = ibool(i,j,k,ispec)
-
- ! isotropic kernels
- ! note: takes displacement from backward/reconstructed (forward) field b_displ
- ! and acceleration from adjoint field accel (containing adjoint sources)
- !
- ! note: : time integral summation uses deltat
- !
- ! compare with Tromp et al. (2005), eq. (14), which takes adjoint displacement
- ! and forward acceleration, that is the symmetric form of what is calculated here
- ! however, this kernel expression is symmetric with regards
- ! to interchange adjoint - forward field
- rho_kl(i,j,k,ispec) = rho_kl(i,j,k,ispec) &
- + deltat * dot_product(accel(:,iglob), b_displ(:,iglob))
-
- ! kernel for shear modulus, see e.g. Tromp et al. (2005), equation (17)
- ! note: multiplication with 2*mu(x) will be done after the time loop
+
epsilondev_loc(1) = epsilondev_xx(i,j,k,ispec)
epsilondev_loc(2) = epsilondev_yy(i,j,k,ispec)
epsilondev_loc(3) = epsilondev_xy(i,j,k,ispec)
@@ -115,18 +101,46 @@
b_epsilondev_loc(4) = b_epsilondev_xz(i,j,k,ispec)
b_epsilondev_loc(5) = b_epsilondev_yz(i,j,k,ispec)
- mu_kl(i,j,k,ispec) = mu_kl(i,j,k,ispec) &
- + deltat * (epsilondev_loc(1)*b_epsilondev_loc(1) + epsilondev_loc(2)*b_epsilondev_loc(2) &
- + (epsilondev_loc(1)+epsilondev_loc(2)) * (b_epsilondev_loc(1)+b_epsilondev_loc(2)) &
- + 2 * (epsilondev_loc(3)*b_epsilondev_loc(3) + epsilondev_loc(4)*b_epsilondev_loc(4) + &
- epsilondev_loc(5)*b_epsilondev_loc(5)) )
+ rho_kl(i,j,k,ispec) = rho_kl(i,j,k,ispec) &
+ + deltat * dot_product(accel(:,iglob), b_displ(:,iglob))
- ! kernel for bulk modulus, see e.g. Tromp et al. (2005), equation (18)
- ! note: multiplication with kappa(x) will be done after the time loop
- kappa_kl(i,j,k,ispec) = kappa_kl(i,j,k,ispec) &
- + deltat * (9 * epsilon_trace_over_3(i,j,k,ispec) &
- * b_epsilon_trace_over_3(i,j,k,ispec))
+ ! For anisotropic kernels
+ if (ANISOTROPIC_KL) then
+ call compute_strain_product(prod,epsilon_trace_over_3(i,j,k,ispec),epsilondev_loc, &
+ b_epsilon_trace_over_3(i,j,k,ispec),b_epsilondev_loc)
+ cijkl_kl(:,i,j,k,ispec) = cijkl_kl(:,i,j,k,ispec) + deltat * prod(:)
+
+ else
+
+ ! isotropic kernels
+ ! note: takes displacement from backward/reconstructed (forward) field b_displ
+ ! and acceleration from adjoint field accel (containing adjoint sources)
+ !
+ ! and acceleration from adjoint field accel (containing adjoint sources)
+ !
+ ! note: : time integral summation uses deltat
+ !
+ ! compare with Tromp et al. (2005), eq. (14), which takes adjoint displacement
+ ! and forward acceleration, that is the symmetric form of what is calculated here
+ ! however, this kernel expression is symmetric with regards
+ ! to interchange adjoint - forward field
+
+ ! kernel for shear modulus, see e.g. Tromp et al. (2005), equation (17)
+ ! note: multiplication with 2*mu(x) will be done after the time loop
+
+ mu_kl(i,j,k,ispec) = mu_kl(i,j,k,ispec) &
+ + deltat * (epsilondev_loc(1)*b_epsilondev_loc(1) + epsilondev_loc(2)*b_epsilondev_loc(2) &
+ + (epsilondev_loc(1)+epsilondev_loc(2)) * (b_epsilondev_loc(1)+b_epsilondev_loc(2)) &
+ + 2 * (epsilondev_loc(3)*b_epsilondev_loc(3) + epsilondev_loc(4)*b_epsilondev_loc(4) + &
+ epsilondev_loc(5)*b_epsilondev_loc(5)) )
+
+ ! kernel for bulk modulus, see e.g. Tromp et al. (2005), equation (18)
+ ! note: multiplication with kappa(x) will be done after the time loop
+ kappa_kl(i,j,k,ispec) = kappa_kl(i,j,k,ispec) &
+ + deltat * (9 * epsilon_trace_over_3(i,j,k,ispec) &
+ * b_epsilon_trace_over_3(i,j,k,ispec))
+ endif
enddo
enddo
enddo
@@ -426,4 +440,65 @@
end subroutine compute_kernels_hessian
+!-------------------------------------------------------------------------------------------------
+!
+! Subroutine to compute the kernels for the 21 elastic coefficients
+! Last modified 19/04/2007
+!-------------------------------------------------------------------
+ subroutine compute_strain_product(prod,eps_trace_over_3,epsdev,&
+ b_eps_trace_over_3,b_epsdev)
+
+ ! Purpose : compute the 21 strain products at a grid point
+ ! (ispec,i,j,k fixed) and at a time t to compute then the kernels cij_kl (Voigt notation)
+ ! (eq. 15 of Tromp et al., 2005)
+ ! prod(1)=eps11*eps11 -> c11, prod(2)=eps11eps22 -> c12, prod(3)=eps11eps33 -> c13, ...
+ ! prod(7)=eps22*eps22 -> c22, prod(8)=eps22eps33 -> c23, prod(9)=eps22eps23 -> c24, ...
+ ! prod(19)=eps13*eps13 -> c55, prod(20)=eps13eps12 -> c56, prod(21)=eps12eps12 -> c66
+ ! This then gives how the 21 kernels are organized
+ ! For crust_mantle
+
+ ! Modif 09/11/2005
+
+ implicit none
+ include "constants.h"
+
+ real(kind=CUSTOM_REAL),dimension(21) :: prod
+ real(kind=CUSTOM_REAL) :: eps_trace_over_3,b_eps_trace_over_3
+ real(kind=CUSTOM_REAL),dimension(5) :: epsdev,b_epsdev
+ real(kind=CUSTOM_REAL), dimension(6) :: eps,b_eps
+ integer :: p,i,j
+
+ ! Building of the local matrix of the strain tensor
+ ! for the adjoint field and the regular backward field
+ eps(1:2)=epsdev(1:2)+eps_trace_over_3 !eps11 et eps22
+ eps(3)=-(eps(1)+eps(2))+3*eps_trace_over_3 !eps33
+ eps(4)=epsdev(5) !eps23
+ eps(5)=epsdev(4) !eps13
+ eps(6)=epsdev(3) !eps12
+
+ b_eps(1:2)=b_epsdev(1:2)+b_eps_trace_over_3
+ b_eps(3)=-(b_eps(1)+b_eps(2))+3*b_eps_trace_over_3
+ b_eps(4)=b_epsdev(5)
+ b_eps(5)=b_epsdev(4)
+ b_eps(6)=b_epsdev(3)
+
+ ! Computing the 21 strain products without assuming eps(i)*b_eps(j) = eps(j)*b_eps(i)
+ p=1
+ do i=1,6
+ do j=i,6
+ prod(p)=eps(i)*b_eps(j)
+ if(j>i) then
+ prod(p)=prod(p)+eps(j)*b_eps(i)
+ if(j>3 .and. i<4) prod(p)=prod(p)*2
+ endif
+ if(i>3) prod(p)=prod(p)*4
+ p=p+1
+ enddo
+ enddo
+
+ end subroutine compute_strain_product
+
+!
+!-------------------------------------------------------------------------------------------------
+
Modified: seismo/3D/SPECFEM3D/branches/full_aniso_kernels_matharu/src/specfem3D/prepare_timerun.F90
===================================================================
--- seismo/3D/SPECFEM3D/branches/full_aniso_kernels_matharu/src/specfem3D/prepare_timerun.F90 2013-04-02 22:52:02 UTC (rev 21699)
+++ seismo/3D/SPECFEM3D/branches/full_aniso_kernels_matharu/src/specfem3D/prepare_timerun.F90 2013-04-03 00:54:29 UTC (rev 21700)
@@ -833,6 +833,7 @@
rho_kl(:,:,:,:) = 0._CUSTOM_REAL
mu_kl(:,:,:,:) = 0._CUSTOM_REAL
kappa_kl(:,:,:,:) = 0._CUSTOM_REAL
+ cijkl_kl(:,:,:,:,:) = 0._CUSTOM_REAL
if ( APPROXIMATE_HESS_KL ) &
hess_kl(:,:,:,:) = 0._CUSTOM_REAL
Modified: seismo/3D/SPECFEM3D/branches/full_aniso_kernels_matharu/src/specfem3D/read_mesh_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/full_aniso_kernels_matharu/src/specfem3D/read_mesh_databases.f90 2013-04-02 22:52:02 UTC (rev 21699)
+++ seismo/3D/SPECFEM3D/branches/full_aniso_kernels_matharu/src/specfem3D/read_mesh_databases.f90 2013-04-03 00:54:29 UTC (rev 21700)
@@ -711,6 +711,10 @@
allocate(kappa_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT),stat=ier)
if( ier /= 0 ) stop 'error allocating array kappa_kl'
+ ! anisotropic kernels
+ allocate(cijkl_kl(21,NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT),stat=ier)
+ if( ier /= 0 ) stop 'error allocating array cijkl_kl'
+
! derived kernels
! density prime kernel
allocate(rhop_kl(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT),stat=ier)
Modified: seismo/3D/SPECFEM3D/branches/full_aniso_kernels_matharu/src/specfem3D/save_adjoint_kernels.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/full_aniso_kernels_matharu/src/specfem3D/save_adjoint_kernels.f90 2013-04-02 22:52:02 UTC (rev 21699)
+++ seismo/3D/SPECFEM3D/branches/full_aniso_kernels_matharu/src/specfem3D/save_adjoint_kernels.f90 2013-04-03 00:54:29 UTC (rev 21700)
@@ -159,7 +159,15 @@
! local parameters
integer:: ispec,i,j,k,iglob,ier
real(kind=CUSTOM_REAL) :: rhol,mul,kappal
+ real(kind=CUSTOM_REAL),dimension(21) :: cijkl_kl_local
+ ! Transverse isotropic paramters
+ real(kind=CUSTOM_REAL) :: A,N,C,L,F,eta
+ real(kind=CUSTOM_REAL), dimension(21) :: an_kl
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ADJOINT):: &
+ alphav_kl,alphah_kl,betav_kl,betah_kl, &
+ eta_kl
+
! finalizes calculation of rhop, beta, alpha kernels
do ispec = 1, NSPEC_AB
@@ -170,41 +178,103 @@
do j = 1, NGLLY
do i = 1, NGLLX
iglob = ibool(i,j,k,ispec)
+
+ ! Store local material values
+ rhol = rho_vs(i,j,k,ispec)*rho_vs(i,j,k,ispec) / mustore(i,j,k,ispec)
+ mul = mustore(i,j,k,ispec)
+ kappal = kappastore(i,j,k,ispec)
- ! isotropic adjoint kernels (see e.g. Tromp et al. 2005)
- rhol = rho_vs(i,j,k,ispec)*rho_vs(i,j,k,ispec) / mustore(i,j,k,ispec)
- mul = mustore(i,j,k,ispec)
- kappal = kappastore(i,j,k,ispec)
+ if (ANISOTROPIC_KL) then
+
+ cijkl_kl_local = - cijkl_kl(:,i,j,k,ispec)
+
+ if (SAVE_TRANSVERSE_KL) then
+
+ ! Computes parameters for an isotropic model
+ A = kappal + FOUR_THIRDS * mul
+ C = A
+ L = mul
+ N = mul
+ F = kappal - 2._CUSTOM_REAL/3._CUSTOM_REAL * mul
+ eta = 1._CUSTOM_REAL
+
+ ! note: cijkl_kl_local() is fully anisotropic C_ij kernel components (non-dimensionalized)
+ ! for GLL point at (i,j,k,ispec)
+
+ ! Purpose : compute the kernels for the An coeffs (an_kl)
+ ! from the kernels for Cij (cijkl_kl_local)
+
+ ! Definition of the input array cij_kl :
+ ! cij_kl(1) = C11 ; cij_kl(2) = C12 ; cij_kl(3) = C13
+ ! cij_kl(4) = C14 ; cij_kl(5) = C15 ; cij_kl(6) = C16
+ ! cij_kl(7) = C22 ; cij_kl(8) = C23 ; cij_kl(9) = C24
+ ! cij_kl(10) = C25 ; cij_kl(11) = C26 ; cij_kl(12) = C33
+ ! cij_kl(13) = C34 ; cij_kl(14) = C35 ; cij_kl(15) = C36
+ ! cij_kl(16) = C44 ; cij_kl(17) = C45 ; cij_kl(18) = C46
+ ! cij_kl(19) = C55 ; cij_kl(20) = C56 ; cij_kl(21) = C66
+ ! where the Cij (Voigt's notation) are defined as function of
+ ! the components of the elastic tensor in spherical coordinates
+ ! by eq. (A.1) of Chen & Tromp, GJI 168 (2007)
+
+ ! From the relations giving Cij in function of An
+ ! Checked with Min Chen's results (routine build_cij)
- ! for a parameterization: (rho,mu,kappa) "primary" kernels
- ! density kernel
- ! multiplies with rho
- rho_kl(i,j,k,ispec) = - rhol * rho_kl(i,j,k,ispec)
+ an_kl(1) = cijkl_kl_local(1)+cijkl_kl_local(2)+cijkl_kl_local(7) !A
+ an_kl(2) = cijkl_kl_local(12) !C
+ an_kl(3) = -2*cijkl_kl_local(2)+cijkl_kl_local(21) !N
+ an_kl(4) = cijkl_kl_local(16)+cijkl_kl_local(19) !L
+ an_kl(5) = cijkl_kl_local(3)+cijkl_kl_local(8) !F
+
+ ! for parameterization: ( alpha_v, alpha_h, beta_v, beta_h, eta, rho )
+ ! K_alpha_v
+ alphav_kl(i,j,k,ispec) = 2*C*an_kl(2)
+ ! K_alpha_h
+ alphah_kl(i,j,k,ispec) = 2*A*an_kl(1) + 2*A*eta*an_kl(5)
+ ! K_beta_v
+ betav_kl(i,j,k,ispec) = 2*L*an_kl(4) - 4*L*eta*an_kl(5)
+ ! K_beta_h
+ betah_kl(i,j,k,ispec) = 2*N*an_kl(3)
+ ! K_eta
+ eta_kl(i,j,k,ispec) = F*an_kl(5)
- ! shear modulus kernel
- mu_kl(i,j,k,ispec) = - 2._CUSTOM_REAL * mul * mu_kl(i,j,k,ispec)
+ ! to check: isotropic kernels from transverse isotropic ones
+ alpha_kl(i,j,k,ispec) = alphav_kl(i,j,k,ispec) &
+ + alphah_kl(i,j,k,ispec)
+ beta_kl(i,j,k,ispec) = betav_kl(i,j,k,ispec) &
+ + betah_kl(i,j,k,ispec)
- ! bulk modulus kernel
- kappa_kl(i,j,k,ispec) = - kappal * kappa_kl(i,j,k,ispec)
+ endif ! SAVE_TRANSVERSE_KL
- ! for a parameterization: (rho,alpha,beta)
- ! density prime kernel
- rhop_kl(i,j,k,ispec) = rho_kl(i,j,k,ispec) + kappa_kl(i,j,k,ispec) + mu_kl(i,j,k,ispec)
+ else
+
+ ! isotropic kernels
- ! vs kernel
- beta_kl(i,j,k,ispec) = 2._CUSTOM_REAL * (mu_kl(i,j,k,ispec) &
- - 4._CUSTOM_REAL * mul / (3._CUSTOM_REAL * kappal) * kappa_kl(i,j,k,ispec))
+ ! isotropic adjoint kernels (see e.g. Tromp et al. 2005)
+ ! for a parameterization: (rho,mu,kappa) "primary" kernels
+ ! density kernel
+ ! multiplies with rho
+ rho_kl(i,j,k,ispec) = - rhol * rho_kl(i,j,k,ispec)
- ! vp kernel
- alpha_kl(i,j,k,ispec) = 2._CUSTOM_REAL * (1._CUSTOM_REAL &
- + 4._CUSTOM_REAL * mul / (3._CUSTOM_REAL * kappal) ) * kappa_kl(i,j,k,ispec)
+ ! shear modulus kernel
+ mu_kl(i,j,k,ispec) = - 2._CUSTOM_REAL * mul * mu_kl(i,j,k,ispec)
- ! for a parameterization: (rho,bulk, beta)
- ! where bulk wave speed is c = sqrt( kappa / rho)
- ! note: rhoprime is the same as for (rho,alpha,beta) parameterization
- !bulk_c_kl_crust_mantle(i,j,k,ispec) = 2._CUSTOM_REAL * kappa_kl(i,j,k,ispec)
- !bulk_beta_kl_crust_mantle(i,j,k,ispec ) = 2._CUSTOM_REAL * mu_kl(i,j,k,ispec)
+ ! bulk modulus kernel
+ kappa_kl(i,j,k,ispec) = - kappal * kappa_kl(i,j,k,ispec)
+ ! for a parameterization: (rho,alpha,beta)
+ ! density prime kernel
+ rhop_kl(i,j,k,ispec) = rho_kl(i,j,k,ispec) + kappa_kl(i,j,k,ispec) + mu_kl(i,j,k,ispec)
+
+ ! vs kernel
+ beta_kl(i,j,k,ispec) = 2._CUSTOM_REAL * (mu_kl(i,j,k,ispec) &
+ - 4._CUSTOM_REAL * mul / (3._CUSTOM_REAL * kappal) * kappa_kl(i,j,k,ispec))
+
+ ! vp kernel
+ alpha_kl(i,j,k,ispec) = 2._CUSTOM_REAL * (1._CUSTOM_REAL &
+ + 4._CUSTOM_REAL * mul / (3._CUSTOM_REAL * kappal) ) * kappa_kl(i,j,k,ispec)
+
+ endif
+
enddo
enddo
enddo
@@ -212,38 +282,85 @@
endif ! elastic
enddo
+
+ if (ANISOTROPIC_KL) then
- ! save kernels to binary files
- open(unit=27,file=prname(1:len_trim(prname))//'rho_kernel.bin',status='unknown',form='unformatted',iostat=ier)
- if( ier /= 0 ) stop 'error opening file rho_kernel.bin'
- write(27) rho_kl
- close(27)
+ ! outputs transverse isotropic kernels only
+ if (SAVE_TRANSVERSE_KL) then
+ ! transverse isotropic kernels
+ ! (alpha_v, alpha_h, beta_v, beta_h, eta, rho ) parameterization
+ open(unit=27,file=trim(prname)//'alphav_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) alphav_kl
+ close(27)
+ open(unit=27,file=trim(prname)//'alphah_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) alphah_kl
+ close(27)
+ open(unit=27,file=trim(prname)//'betav_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) betav_kl
+ close(27)
+ open(unit=27,file=trim(prname)//'betah_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) betah_kl
+ close(27)
+ open(unit=27,file=trim(prname)//'eta_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) eta_kl
+ close(27)
- open(unit=27,file=prname(1:len_trim(prname))//'mu_kernel.bin',status='unknown',form='unformatted',iostat=ier)
- if( ier /= 0 ) stop 'error opening file mu_kernel.bin'
- write(27) mu_kl
- close(27)
+ ! transverse isotropic test kernels
+ open(unit=27,file=trim(prname)//'alpha_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) alpha_kl
+ close(27)
+ open(unit=27,file=trim(prname)//'beta_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) beta_kl
+ close(27)
- open(unit=27,file=prname(1:len_trim(prname))//'kappa_kernel.bin',status='unknown',form='unformatted',iostat=ier)
- if( ier /= 0 ) stop 'error opening file kappa_kernel.bin'
- write(27) kappa_kl
- close(27)
+ else
+ ! fully anisotropic kernels
+ ! note: the C_ij and density kernels are not for relative perturbations (delta ln( m_i) = delta m_i / m_i),
+ ! but absolute perturbations (delta m_i = m_i - m_0).
+ ! Kappa and mu are for absolute perturbations, can be used to check with purely isotropic versions.
+ open(unit=27,file=trim(prname)//'rho_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) - rho_kl
+ close(27)
+ open(unit=27,file=trim(prname)//'cijkl_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) - cijkl_kl
+ close(27)
- open(unit=27,file=prname(1:len_trim(prname))//'rhop_kernel.bin',status='unknown',form='unformatted',iostat=ier)
- if( ier /= 0 ) stop 'error opening file rhop_kernel.bin'
- write(27) rhop_kl
- close(27)
+ endif
- open(unit=27,file=prname(1:len_trim(prname))//'beta_kernel.bin',status='unknown',form='unformatted',iostat=ier)
- if( ier /= 0 ) stop 'error opening file beta_kernel.bin'
- write(27) beta_kl
- close(27)
+ else
+
+ ! save kernels to binary files
+ open(unit=27,file=prname(1:len_trim(prname))//'rho_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file rho_kernel.bin'
+ write(27) rho_kl
+ close(27)
- open(unit=27,file=prname(1:len_trim(prname))//'alpha_kernel.bin',status='unknown',form='unformatted',iostat=ier)
- if( ier /= 0 ) stop 'error opening file alpha_kernel.bin'
- write(27) alpha_kl
- close(27)
+ open(unit=27,file=prname(1:len_trim(prname))//'mu_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file mu_kernel.bin'
+ write(27) mu_kl
+ close(27)
+ open(unit=27,file=prname(1:len_trim(prname))//'kappa_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file kappa_kernel.bin'
+ write(27) kappa_kl
+ close(27)
+
+ open(unit=27,file=prname(1:len_trim(prname))//'rhop_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file rhop_kernel.bin'
+ write(27) rhop_kl
+ close(27)
+
+ open(unit=27,file=prname(1:len_trim(prname))//'beta_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file beta_kernel.bin'
+ write(27) beta_kl
+ close(27)
+
+ open(unit=27,file=prname(1:len_trim(prname))//'alpha_kernel.bin',status='unknown',form='unformatted',iostat=ier)
+ if( ier /= 0 ) stop 'error opening file alpha_kernel.bin'
+ write(27) alpha_kl
+ close(27)
+ endif
+
if (SAVE_MOHO_MESH) then
open(unit=27,file=prname(1:len_trim(prname))//'moho_kernel.bin',status='unknown',form='unformatted',iostat=ier)
if( ier /= 0 ) stop 'error opening file moho_kernel.bin'
Modified: seismo/3D/SPECFEM3D/branches/full_aniso_kernels_matharu/src/specfem3D/specfem3D_par.f90
===================================================================
--- seismo/3D/SPECFEM3D/branches/full_aniso_kernels_matharu/src/specfem3D/specfem3D_par.f90 2013-04-02 22:52:02 UTC (rev 21699)
+++ seismo/3D/SPECFEM3D/branches/full_aniso_kernels_matharu/src/specfem3D/specfem3D_par.f90 2013-04-03 00:54:29 UTC (rev 21700)
@@ -361,7 +361,10 @@
! adjoint kernels
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: rho_kl, mu_kl, kappa_kl, &
rhop_kl, beta_kl, alpha_kl
-
+
+ ! anisotropic kernels
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: cijkl_kl
+
! approximate hessian
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: hess_kl
More information about the CIG-COMMITS
mailing list