[cig-commits] r22485 - in seismo/3D/SPECFEM3D_GLOBE: branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D trunk/src/specfem3D
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Mon Jul 1 19:26:51 PDT 2013
Author: dkomati1
Date: 2013-07-01 19:26:51 -0700 (Mon, 01 Jul 2013)
New Revision: 22485
Modified:
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_element.F90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/get_attenuation.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_forward_arrays.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_element.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_stacey_crust_mantle.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/get_attenuation.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/get_event_info.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_receivers.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/noise_tomography.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/read_forward_arrays.f90
Log:
more easy merges in src/specfem3D
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_element.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_element.F90 2013-07-02 01:23:11 UTC (rev 22484)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/compute_element.F90 2013-07-02 02:26:51 UTC (rev 22485)
@@ -411,8 +411,6 @@
end subroutine compute_element_iso
-
-
!
!--------------------------------------------------------------------------------------------------
!
@@ -479,7 +477,6 @@
! [alpha,beta,gamma]val reduced to N_SLS to N_SLS*NUM_NODES
real(kind=CUSTOM_REAL), dimension(vx,vy,vz,vnspec) :: one_minus_sum_beta
-
! gravity
double precision, dimension(NRAD_GRAVITY) :: minus_gravity_table,density_table,minus_deriv_gravity_table
@@ -838,7 +835,6 @@
+ rhovpvsq*cosphisq*sinphisq*sinthetafour &
- 0.5_CUSTOM_REAL*eta_aniso*sintwophisq*sinthetafour*(rhovphsq - two_rhovsvsq)
-
! general expression of stress tensor for full Cijkl with 21 coefficients
sigma_xx = c11*duxdxl + c16*duxdyl_plus_duydxl + c12*duydyl + &
c15*duzdxl_plus_duxdzl + c14*duzdyl_plus_duydzl + c13*duzdzl
@@ -1012,7 +1008,6 @@
enddo ! NGLLY
enddo ! NGLLZ
-
end subroutine compute_element_tiso
!
@@ -1038,7 +1033,6 @@
tempz1_att,tempz2_att,tempz3_att, &
epsilondev_loc,rho_s_H,is_backward_field)
-
! this routine is optimized for NGLLX = NGLLY = NGLLZ = 5 using the Deville et al. (2002) inlined matrix-matrix products
implicit none
@@ -1636,7 +1630,6 @@
! IMPROVE we use mu_v here even if there is some anisotropy
! IMPROVE we should probably use an average value instead
-
#ifdef _HANDOPT_ATT
! way 2:
do i_SLS = 1,N_SLS
@@ -1724,238 +1717,12 @@
enddo ! i_SLS
-!daniel: att - debug
-! if( .not. is_backward_field ) then
-!
-! do i_SLS = 1,N_SLS
-! ! reformatted R_memory to handle large factor_common and reduced [alpha,beta,gamma]val
-! if( USE_3D_ATTENUATION_ARRAYS ) then
-! if(ANISOTROPIC_3D_MANTLE_VAL) then
-! factor_common_c44_muv(:,:,:) = factor_common(i_SLS,:,:,:,ispec) * c44store(:,:,:,ispec)
-! else
-! factor_common_c44_muv(:,:,:) = factor_common(i_SLS,:,:,:,ispec) * muvstore(:,:,:,ispec)
-! endif
-! else
-! if(ANISOTROPIC_3D_MANTLE_VAL) then
-! factor_common_c44_muv(:,:,:) = factor_common(i_SLS,1,1,1,ispec) * c44store(:,:,:,ispec)
-! else
-! factor_common_c44_muv(:,:,:) = factor_common(i_SLS,1,1,1,ispec) * muvstore(:,:,:,ispec)
-! endif
-! endif
-!
-!!daniel: att - debug original
-! R_xx(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_xx(i_SLS,:,:,:,ispec) + factor_common_c44_muv(:,:,:) * &
-! (betaval(i_SLS) * epsilondev_xx(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(1,:,:,:))
-!
-! R_yy(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_yy(i_SLS,:,:,:,ispec) + factor_common_c44_muv(:,:,:) * &
-! (betaval(i_SLS) * epsilondev_yy(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(2,:,:,:))
-!
-! R_xy(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_xy(i_SLS,:,:,:,ispec) + factor_common_c44_muv(:,:,:) * &
-! (betaval(i_SLS) * epsilondev_xy(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(3,:,:,:))
-!
-! R_xz(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_xz(i_SLS,:,:,:,ispec) + factor_common_c44_muv(:,:,:) * &
-! (betaval(i_SLS) * epsilondev_xz(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(4,:,:,:))
-!
-! R_yz(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_yz(i_SLS,:,:,:,ispec) + factor_common_c44_muv(:,:,:) * &
-! (betaval(i_SLS) * epsilondev_yz(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(5,:,:,:))
-!
-!!daniel: att - debug runge-kutta
-! if( .false. ) then
-! ! classical RK4: R'(t) = - 1/tau * R(t)
-! kappa = - dble(deltat)/tau_sigma_dble(i_SLS)
-!
-! R_xx(i_SLS,:,:,:,ispec) = R_xx(i_SLS,:,:,:,ispec) * &
-! (1.0d0 + kappa*(1.d0 + 0.5d0*kappa*(1.d0 + 1.0d0/6.0d0*kappa*(1.d0 + 1.0d0/24.0d0*kappa))))
-!
-! R_yy(i_SLS,:,:,:,ispec) = R_yy(i_SLS,:,:,:,ispec) * &
-! (1.0d0 + kappa*(1.d0 + 0.5d0*kappa*(1.d0 + 1.0d0/6.0d0*kappa*(1.d0 + 1.0d0/24.0d0*kappa))))
-!
-! R_xy(i_SLS,:,:,:,ispec) = R_xy(i_SLS,:,:,:,ispec) * &
-! (1.0d0 + kappa*(1.d0 + 0.5d0*kappa*(1.d0 + 1.0d0/6.0d0*kappa*(1.d0 + 1.0d0/24.0d0*kappa))))
-!
-! R_xz(i_SLS,:,:,:,ispec) = R_xz(i_SLS,:,:,:,ispec) * &
-! (1.0d0 + kappa*(1.d0 + 0.5d0*kappa*(1.d0 + 1.0d0/6.0d0*kappa*(1.d0 + 1.0d0/24.0d0*kappa))))
-!
-! R_yz(i_SLS,:,:,:,ispec) = R_yz(i_SLS,:,:,:,ispec) * &
-! (1.0d0 + kappa*(1.d0 + 0.5d0*kappa*(1.d0 + 1.0d0/6.0d0*kappa*(1.d0 + 1.0d0/24.0d0*kappa))))
-!
-! endif
-! if( .false. ) then
-! Butcher RK5:
-! 0 |
-! 1/4 | 1/4
-! 1/4 | 1/8 1/8
-! 1/2 | 0 -1/2 1
-! 3/4 | 3/16 0 0 9/16
-! 1 | -3/7 2/7 12/7 -12/7 8/7
-! -----------------------------------------------------------------------------
-! 7/90 0 32/90 12/90 32/90 7/90
-!
-! k1 = dt * ( -1.0/tau * R_n )
-! k2 = dt * ( -1.0/tau * (R_n + 1./4.*k1) )
-! k3 = dt * ( -1.0/tau * (R_n + 1./8.*k1 + 1./8.*k2) )
-! k4 = dt * ( -1.0/tau * (R_n + 0.*k1 - 1./2.*k2 + 1.* k3) )
-! k5 = dt * ( -1.0/tau * (R_n + 3./16.*k1 + 0.*k2 + 0.*k3 + 9./16.*k4) )
-! k6 = dt * ( -1.0/tau * (R_n - 3./7.*k1 + 2./7.*k2 + 12./7.*k3 - 12./7.*k4 + 8./7.*k5) )
-!
-! R_nplus1 = R_n + 7./90.*k1 + 0.*k2 + 32./90.*k3 + 12./90.*k4 + 32./90.*k5 + 7./90.*k6
-!
-! or:
-! k = - dt/tau
-! R_nplus1 = R_n*(1.0 + k*(1.0 + 0.5*k*(1.0 + 1./6.*k*(1.0 + 1./24.*k*(1.0 + 1./120.*k*(1.0 + 1./640.*k))))))
-!
-! kappa = - dble(deltat)/tau_sigma_dble(i_SLS)
-!
-! R_xx(i_SLS,:,:,:,ispec) = R_xx(i_SLS,:,:,:,ispec) * &
-! (1.d0 + kappa*(1.d0 + 0.5d0*kappa*(1.d0 + 1.d0/6.d0*kappa*(1.d0 + &
-! 1.d0/24.d0*kappa*(1.d0 + 1.d0/120.d0*kappa*(1.d0 + 1.d0/640.d0*kappa))))))
-!
-! R_yy(i_SLS,:,:,:,ispec) = R_yy(i_SLS,:,:,:,ispec) * &
-! (1.d0 + kappa*(1.d0 + 0.5d0*kappa*(1.d0 + 1.d0/6.d0*kappa*(1.d0 + &
-! 1.d0/24.d0*kappa*(1.d0 + 1.d0/120.d0*kappa*(1.d0 + 1.d0/640.d0*kappa))))))
-!
-! R_xy(i_SLS,:,:,:,ispec) = R_xy(i_SLS,:,:,:,ispec) * &
-! (1.d0 + kappa*(1.d0 + 0.5d0*kappa*(1.d0 + 1.d0/6.d0*kappa*(1.d0 + &
-! 1.d0/24.d0*kappa*(1.d0 + 1.d0/120.d0*kappa*(1.d0 + 1.d0/640.d0*kappa))))))
-!
-! R_xz(i_SLS,:,:,:,ispec) = R_xz(i_SLS,:,:,:,ispec) * &
-! (1.d0 + kappa*(1.d0 + 0.5d0*kappa*(1.d0 + 1.d0/6.d0*kappa*(1.d0 + &
-! 1.d0/24.d0*kappa*(1.d0 + 1.d0/120.d0*kappa*(1.d0 + 1.d0/640.d0*kappa))))))
-!
-! R_yz(i_SLS,:,:,:,ispec) = R_yz(i_SLS,:,:,:,ispec) * &
-! (1.d0 + kappa*(1.d0 + 0.5d0*kappa*(1.d0 + 1.d0/6.d0*kappa*(1.d0 + &
-! 1.d0/24.d0*kappa*(1.d0 + 1.d0/120.d0*kappa*(1.d0 + 1.d0/640.d0*kappa))))))
-!
-! endif
-!
-! enddo ! i_SLS
-!
-! else
-!
-! ! backward/reconstruction, re-constructs previous memory variables, note strain arrays are now switched
-!
-! do i_SLS = 1,N_SLS
-! ! reformatted R_memory to handle large factor_common and reduced [alpha,beta,gamma]val
-! if( USE_3D_ATTENUATION_ARRAYS ) then
-! if(ANISOTROPIC_3D_MANTLE_VAL) then
-! factor_common_c44_muv(:,:,:) = factor_common(i_SLS,:,:,:,ispec) * c44store(:,:,:,ispec)
-! else
-! factor_common_c44_muv(:,:,:) = factor_common(i_SLS,:,:,:,ispec) * muvstore(:,:,:,ispec)
-! endif
-! else
-! if(ANISOTROPIC_3D_MANTLE_VAL) then
-! factor_common_c44_muv(:,:,:) = factor_common(i_SLS,1,1,1,ispec) * c44store(:,:,:,ispec)
-! else
-! factor_common_c44_muv(:,:,:) = factor_common(i_SLS,1,1,1,ispec) * muvstore(:,:,:,ispec)
-! endif
-! endif
-!
-!!daniel: att - debug original
-! R_xx(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_xx(i_SLS,:,:,:,ispec) + factor_common_c44_muv(:,:,:) * &
-! (betaval(i_SLS) * epsilondev_xx(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(1,:,:,:))
-!
-! R_yy(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_yy(i_SLS,:,:,:,ispec) + factor_common_c44_muv(:,:,:) * &
-! (betaval(i_SLS) * epsilondev_yy(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(2,:,:,:))
-!
-! R_xy(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_xy(i_SLS,:,:,:,ispec) + factor_common_c44_muv(:,:,:) * &
-! (betaval(i_SLS) * epsilondev_xy(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(3,:,:,:))
-!
-! R_xz(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_xz(i_SLS,:,:,:,ispec) + factor_common_c44_muv(:,:,:) * &
-! (betaval(i_SLS) * epsilondev_xz(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(4,:,:,:))
-!
-! R_yz(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_yz(i_SLS,:,:,:,ispec) + factor_common_c44_muv(:,:,:) * &
-! (betaval(i_SLS) * epsilondev_yz(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(5,:,:,:))
-!
-!!daniel: att - debug runge-kutta
-! if( .false. ) then
-! classical RK 4: R'(t) = - 1/tau * R(t)
-!
-! Butcher RK5:
-! 0 |
-! 1/2 | 1/2
-! 1/2 | 0 1/2
-! 1 | 0 1
-! -----------------------------------------------------------------------------
-! 1/6 1/3 1/3 1/6
-!
-! kappa = - dble(b_deltat)/tau_sigma_dble(i_SLS)
-!
-! R_xx(i_SLS,:,:,:,ispec) = R_xx(i_SLS,:,:,:,ispec) * &
-! (1.0d0 + kappa*(1.d0 + 0.5d0*kappa*(1.d0 + 1.0d0/6.0d0*kappa*(1.d0 + 1.0d0/24.0d0*kappa))))
-!
-! R_yy(i_SLS,:,:,:,ispec) = R_yy(i_SLS,:,:,:,ispec) * &
-! (1.0d0 + kappa*(1.d0 + 0.5d0*kappa*(1.d0 + 1.0d0/6.0d0*kappa*(1.d0 + 1.0d0/24.0d0*kappa))))
-!
-! R_xy(i_SLS,:,:,:,ispec) = R_xy(i_SLS,:,:,:,ispec) * &
-! (1.0d0 + kappa*(1.d0 + 0.5d0*kappa*(1.d0 + 1.0d0/6.0d0*kappa*(1.d0 + 1.0d0/24.0d0*kappa))))
-!
-! R_xz(i_SLS,:,:,:,ispec) = R_xz(i_SLS,:,:,:,ispec) * &
-! (1.0d0 + kappa*(1.d0 + 0.5d0*kappa*(1.d0 + 1.0d0/6.0d0*kappa*(1.d0 + 1.0d0/24.0d0*kappa))))
-!
-! R_yz(i_SLS,:,:,:,ispec) = R_yz(i_SLS,:,:,:,ispec) * &
-! (1.0d0 + kappa*(1.d0 + 0.5d0*kappa*(1.d0 + 1.0d0/6.0d0*kappa*(1.d0 + 1.0d0/24.0d0*kappa))))
-!
-! endif
-!
-!!daniel: att - debug
-! if( .false. ) then
-! Butcher RK5:
-! 0 |
-! 1/4 | 1/4
-! 1/4 | 1/8 1/8
-! 1/2 | 0 -1/2 1
-! 3/4 | 3/16 0 0 9/16
-! 1 | -3/7 2/7 12/7 -12/7 8/7
-! -----------------------------------------------------------------------------
-! 7/90 0 32/90 12/90 32/90 7/90
-!
-! k1 = dt * ( -1.0/tau * R_n )
-! k2 = dt * ( -1.0/tau * (R_n + 1./4.*k1) )
-! k3 = dt * ( -1.0/tau * (R_n + 1./8.*k1 + 1./8.*k2) )
-! k4 = dt * ( -1.0/tau * (R_n + 0.*k1 - 1./2.*k2 + 1.* k3) )
-! k5 = dt * ( -1.0/tau * (R_n + 3./16.*k1 + 0.*k2 + 0.*k3 + 9./16.*k4) )
-! k6 = dt * ( -1.0/tau * (R_n - 3./7.*k1 + 2./7.*k2 + 12./7.*k3 - 12./7.*k4 + 8./7.*k5) )
-!
-! R_nplus1 = R_n + 7./90.*k1 + 0.*k2 + 32./90.*k3 + 12./90.*k4 + 32./90.*k5 + 7./90.*k6
-!
-! or:
-! k = - dt/tau
-! R_nplus1 = R_n*(1.0 + k*(1.0 + 0.5*k*(1.0 + 1./6.*k*(1.0 + 1./24.*k*(1.0 + 1./120.*k*(1.0 + 1./640.*k))))))
-! kappa = - dble(b_deltat)/tau_sigma_dble(i_SLS)
-!
-! R_xx(i_SLS,:,:,:,ispec) = R_xx(i_SLS,:,:,:,ispec) * &
-! (1.d0 + kappa*(1.d0 + 0.5d0*kappa*(1.d0 + 1.d0/6.d0*kappa*(1.d0 + &
-! 1.d0/24.d0*kappa*(1.d0 + 1.d0/120.d0*kappa*(1.d0 + 1.d0/640.d0*kappa))))))
-!
-! R_yy(i_SLS,:,:,:,ispec) = R_yy(i_SLS,:,:,:,ispec) * &
-! (1.d0 + kappa*(1.d0 + 0.5d0*kappa*(1.d0 + 1.d0/6.d0*kappa*(1.d0 + &
-! 1.d0/24.d0*kappa*(1.d0 + 1.d0/120.d0*kappa*(1.d0 + 1.d0/640.d0*kappa))))))
-!
-! R_xy(i_SLS,:,:,:,ispec) = R_xy(i_SLS,:,:,:,ispec) * &
-! (1.d0 + kappa*(1.d0 + 0.5d0*kappa*(1.d0 + 1.d0/6.d0*kappa*(1.d0 + &
-! 1.d0/24.d0*kappa*(1.d0 + 1.d0/120.d0*kappa*(1.d0 + 1.d0/640.d0*kappa))))))
-!
-! R_xz(i_SLS,:,:,:,ispec) = R_xz(i_SLS,:,:,:,ispec) * &
-! (1.d0 + kappa*(1.d0 + 0.5d0*kappa*(1.d0 + 1.d0/6.d0*kappa*(1.d0 + &
-! 1.d0/24.d0*kappa*(1.d0 + 1.d0/120.d0*kappa*(1.d0 + 1.d0/640.d0*kappa))))))
-!
-! R_yz(i_SLS,:,:,:,ispec) = R_yz(i_SLS,:,:,:,ispec) * &
-! (1.d0 + kappa*(1.d0 + 0.5d0*kappa*(1.d0 + 1.d0/6.d0*kappa*(1.d0 + &
-! 1.d0/24.d0*kappa*(1.d0 + 1.d0/120.d0*kappa*(1.d0 + 1.d0/640.d0*kappa))))))
-!
-! endif
-!
-! enddo ! i_SLS
-!
-! endif ! is_backward_field
! dummy to avoid compiler warning
if( is_backward_field ) then
endif
-
#endif
-
end subroutine compute_element_att_memory_cm
!
@@ -2106,101 +1873,12 @@
enddo
-!daniel: att - debug
-! if( .not. is_backward_field ) then
-! do i_SLS = 1,N_SLS
-! ! reformatted R_memory to handle large factor_common and reduced [alpha,beta,gamma]val
-! if( USE_3D_ATTENUATION_ARRAYS ) then
-! factor_common_use(:,:,:) = factor_common(i_SLS,:,:,:,ispec) * muvstore(:,:,:,ispec)
-! else
-! factor_common_use(:,:,:) = factor_common(i_SLS,1,1,1,ispec) * muvstore(:,:,:,ispec)
-! endif
-!
-! do i_memory = 1,5
-! R_memory(i_memory,i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_memory(i_memory,i_SLS,:,:,:,ispec) &
-! + muvstore(:,:,:,ispec) * factor_common_use(:,:,:) * &
-! (betaval(i_SLS) * epsilondev(i_memory,:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(i_memory,:,:,:))
-! enddo
-!
-! R_xx(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_xx(i_SLS,:,:,:,ispec) + factor_common_use(:,:,:) * &
-! (betaval(i_SLS) * epsilondev_xx(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(1,:,:,:))
-!
-! R_yy(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_yy(i_SLS,:,:,:,ispec) + factor_common_use(:,:,:) * &
-! (betaval(i_SLS) * epsilondev_yy(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(2,:,:,:))
-!
-! R_xy(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_xy(i_SLS,:,:,:,ispec) + factor_common_use(:,:,:) * &
-! (betaval(i_SLS) * epsilondev_xy(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(3,:,:,:))
-!
-! R_xz(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_xz(i_SLS,:,:,:,ispec) + factor_common_use(:,:,:) * &
-! (betaval(i_SLS) * epsilondev_xz(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(4,:,:,:))
-!
-! R_yz(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_yz(i_SLS,:,:,:,ispec) + factor_common_use(:,:,:) * &
-! (betaval(i_SLS) * epsilondev_yz(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(5,:,:,:))
-!
-! enddo
-!
-! else
-!
-! ! backward/reconstruction, strain arrays are not switched
-! do i_SLS = 1,N_SLS
-! ! reformatted R_memory to handle large factor_common and reduced [alpha,beta,gamma]val
-! if( USE_3D_ATTENUATION_ARRAYS ) then
-! factor_common_use(:,:,:) = factor_common(i_SLS,:,:,:,ispec) * muvstore(:,:,:,ispec)
-! else
-! factor_common_use(:,:,:) = factor_common(i_SLS,1,1,1,ispec) * muvstore(:,:,:,ispec)
-! endif
-!
-!daniel: att - debug original
-! R_xx(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_xx(i_SLS,:,:,:,ispec) + factor_common_use(:,:,:) * &
-! (betaval(i_SLS) * epsilondev_xx(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(1,:,:,:))
-!
-! R_yy(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_yy(i_SLS,:,:,:,ispec) + factor_common_use(:,:,:) * &
-! (betaval(i_SLS) * epsilondev_yy(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(2,:,:,:))
-!
-! R_xy(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_xy(i_SLS,:,:,:,ispec) + factor_common_use(:,:,:) * &
-! (betaval(i_SLS) * epsilondev_xy(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(3,:,:,:))
-!
-! R_xz(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_xz(i_SLS,:,:,:,ispec) + factor_common_use(:,:,:) * &
-! (betaval(i_SLS) * epsilondev_xz(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(4,:,:,:))
-!
-! R_yz(i_SLS,:,:,:,ispec) = alphaval(i_SLS) * R_yz(i_SLS,:,:,:,ispec) + factor_common_use(:,:,:) * &
-! (betaval(i_SLS) * epsilondev_yz(:,:,:,ispec) + gammaval(i_SLS) * epsilondev_loc(5,:,:,:))
-!
-!daniel: att - debug switched
-! if( .false.) then
-! R_xx(i_SLS,:,:,:,ispec) = 1.0_CUSTOM_REAL/alphaval(i_SLS) * R_xx(i_SLS,:,:,:,ispec) - factor_common_use(:,:,:) * &
-! (gammaval(i_SLS) * epsilondev_xx(:,:,:,ispec) + betaval(i_SLS) * epsilondev_loc(1,:,:,:))
-!
-! R_yy(i_SLS,:,:,:,ispec) = 1.0_CUSTOM_REAL/alphaval(i_SLS) * R_yy(i_SLS,:,:,:,ispec) - factor_common_use(:,:,:) * &
-! (gammaval(i_SLS) * epsilondev_yy(:,:,:,ispec) + betaval(i_SLS) * epsilondev_loc(2,:,:,:))
-!
-! R_xy(i_SLS,:,:,:,ispec) = 1.0_CUSTOM_REAL/alphaval(i_SLS) * R_xy(i_SLS,:,:,:,ispec) - factor_common_use(:,:,:) * &
-! (gammaval(i_SLS) * epsilondev_xy(:,:,:,ispec) + betaval(i_SLS) * epsilondev_loc(3,:,:,:))
-!
-! R_xz(i_SLS,:,:,:,ispec) = 1.0_CUSTOM_REAL/alphaval(i_SLS) * R_xz(i_SLS,:,:,:,ispec) - factor_common_use(:,:,:) * &
-! (gammaval(i_SLS) * epsilondev_xz(:,:,:,ispec) + betaval(i_SLS) * epsilondev_loc(4,:,:,:))
-!
-! R_yz(i_SLS,:,:,:,ispec) = 1.0_CUSTOM_REAL/alphaval(i_SLS) * R_yz(i_SLS,:,:,:,ispec) - factor_common_use(:,:,:) * &
-! (gammaval(i_SLS) * epsilondev_yz(:,:,:,ispec) + betaval(i_SLS) * epsilondev_loc(5,:,:,:))
-! endif
-!
-! enddo
-!
-! endif ! is_backward_field
! dummy to avoid compiler warning
if( is_backward_field ) then
endif
#endif
-
-!daniel: att - debug no att
-! R_xx = 0.0
-! R_yy = 0.0
-! R_xy = 0.0
-! R_xz = 0.0
-! R_yz = 0.0
-
end subroutine compute_element_att_memory_ic
@@ -2220,7 +1898,6 @@
use specfem_par_crustmantle,only: factor_common=>factor_common_crust_mantle
-
implicit none
include "constants.h"
@@ -2253,14 +1930,6 @@
real(kind=CUSTOM_REAL) :: factor_common_c44_muv
integer :: i_SLS
- ! R'(t) = - 1/tau * R(t) + 1/tau * dc:grad(s)
- !
- ! here we add integral contribution for: 1/tau * dc:grad(s)
- !
- ! R_n+1 = R_n+1 + dt * ( 1/tau * dc:grad(s_n+1)
- !
- ! note: s at this point is known at: s(t+dt)
-
if( .not. is_backward_field ) then
dt = dble(deltat)
else
@@ -2295,7 +1964,6 @@
(1.0d0 + kappa*(1.d0 + 0.5d0*kappa*(1.d0 + 1.0d0/6.0d0*kappa*(1.d0 + 1.0d0/24.0d0*kappa))))
endif
-
! reformatted R_memory to handle large factor_common and reduced [alpha,beta,gamma]val
if( USE_3D_ATTENUATION_ARRAYS ) then
factor_common_c44_muv = factor_common(i_SLS,i,j,k,ispec) * c44_muv
@@ -2310,7 +1978,6 @@
R_xz_loc(i_SLS) = R_xz_loc(i_SLS) + 0.5d0 * dt * dble(factor_common_c44_muv) * dble(epsilondev_loc(4))
R_yz_loc(i_SLS) = R_yz_loc(i_SLS) + 0.5d0 * dt * dble(factor_common_c44_muv) * dble(epsilondev_loc(5))
-
enddo ! i_SLS
end subroutine compute_element_att_mem_up_cm
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/get_attenuation.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/get_attenuation.f90 2013-07-02 01:23:11 UTC (rev 22484)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/get_attenuation.f90 2013-07-02 02:26:51 UTC (rev 22485)
@@ -25,7 +25,6 @@
!
!=====================================================================
-
subroutine get_attenuation_model_3D_or_1D(myrank, prname, &
one_minus_sum_beta, &
factor_common, &
@@ -41,8 +40,20 @@
integer :: myrank
integer :: vx,vy,vz,vnspec
+
+!! DK DK to Daniel, Jul 2013
+!! DK DK to Daniel, Jul 2013
+!! DK DK to Daniel, Jul 2013
+!! DK DK to Daniel, Jul 2013
+!! DK DK to Daniel, Jul 2013: BEWARE, declared real(kind=CUSTOM_REAL) in trunk and
+!! DK DK to Daniel, Jul 2013: double precision in branch, let us check which one is right
+!! DK DK to Daniel, Jul 2013
+!! DK DK to Daniel, Jul 2013
+!! DK DK to Daniel, Jul 2013
+!! DK DK to Daniel, Jul 2013
double precision, dimension(vx,vy,vz,vnspec) :: one_minus_sum_beta, scale_factor
double precision, dimension(N_SLS,vx,vy,vz,vnspec) :: factor_common
+
double precision, dimension(N_SLS) :: tau_s
character(len=150) :: prname
@@ -141,10 +152,16 @@
beta(:) = 1.0d0 - tau_e(:) / tau_s(:)
one_minus_sum_beta = 1.0d0
- do i = 1, N_SLS
+ do i = 1,N_SLS
one_minus_sum_beta = one_minus_sum_beta - beta(i)
enddo
+!ZN beware, here the expression differs from the strain used in memory variable equation (6) in D. Komatitsch and J. Tromp 1999,
+!ZN here Brian Savage uses the engineering strain which are epsilon = 1/2*(grad U + (grad U)^T),
+!ZN where U is the displacement vector and grad the gradient operator, i.e. there is a 1/2 factor difference between the two.
+!ZN Both expressions are fine, but we need to keep in mind that if we have put the 1/2 factor there we need to remove it
+!ZN from the expression in which we use the strain here in the code.
+!ZN This is why here Brian Savage multiplies beta(:) * tauinv(:) by 2.0 to compensate for the 1/2 factor used before
factor_common(:) = 2.0d0 * beta(:) * tauinv(:)
end subroutine get_attenuation_property_values
@@ -215,14 +232,12 @@
end subroutine get_attenuation_scale_factor
-
!
!-------------------------------------------------------------------------------------------------
!
+ subroutine get_attenuation_memory_values(tau_s,deltat, alphaval,betaval,gammaval)
- subroutine get_attenuation_memory_values(tau_s, deltat, alphaval,betaval,gammaval)
-
implicit none
include 'constants.h'
@@ -243,517 +258,3 @@
+ deltat**3*tauinv(:)**2 / 24.d0
end subroutine get_attenuation_memory_values
-
-
-!
-!-------------------------------------------------------------------------------------------------
-!
-! not used anymore...
-!
-! subroutine get_attenuation_model_1D(myrank, prname, iregion_code, tau_s, one_minus_sum_beta, &
-! factor_common, scale_factor, vn,vx,vy,vz, AM_V)
-!
-! implicit none
-!
-! include 'mpif.h'
-! include 'constants.h'
-!
-!! model_attenuation_variables
-! type model_attenuation_variables
-! sequence
-! double precision min_period, max_period
-! double precision :: QT_c_source ! Source Frequency
-! double precision, dimension(:), pointer :: Qtau_s ! tau_sigma
-! double precision, dimension(:), pointer :: QrDisc ! Discontinutitues Defined
-! double precision, dimension(:), pointer :: Qr ! Radius
-! integer, dimension(:), pointer :: interval_Q ! Steps
-! double precision, dimension(:), pointer :: Qmu ! Shear Attenuation
-! double precision, dimension(:,:), pointer :: Qtau_e ! tau_epsilon
-! double precision, dimension(:), pointer :: Qomsb, Qomsb2 ! one_minus_sum_beta
-! double precision, dimension(:,:), pointer :: Qfc, Qfc2 ! factor_common
-! double precision, dimension(:), pointer :: Qsf, Qsf2 ! scale_factor
-! integer, dimension(:), pointer :: Qrmin ! Max and Mins of idoubling
-! integer, dimension(:), pointer :: Qrmax ! Max and Mins of idoubling
-! integer :: Qn ! Number of points
-! integer dummy_pad ! padding 4 bytes to align the structure
-! end type model_attenuation_variables
-!
-! type (model_attenuation_variables) AM_V
-!! model_attenuation_variables
-!
-! integer myrank, iregion_code
-! character(len=150) prname
-! integer vn, vx,vy,vz
-! double precision, dimension(N_SLS) :: tau_s
-! double precision, dimension(vx,vy,vz,vn) :: scale_factor, one_minus_sum_beta
-! double precision, dimension(N_SLS, vx,vy,vz,vn) :: factor_common
-!
-! integer i,j,ier,rmax
-! double precision scale_t
-! double precision Qp1, Qpn, radius, fctmp
-! double precision, dimension(:), allocatable :: Qfctmp, Qfc2tmp
-!
-! integer, save :: first_time_called = 1
-!
-! if(myrank == 0 .AND. iregion_code == IREGION_CRUST_MANTLE .AND. first_time_called == 1) then
-! first_time_called = 0
-! open(unit=27, file=prname(1:len_trim(prname))//'1D_Q.bin', status='unknown', form='unformatted')
-! read(27) AM_V%QT_c_source
-! read(27) tau_s
-! read(27) AM_V%Qn
-!
-! allocate(AM_V%Qr(AM_V%Qn))
-! allocate(AM_V%Qmu(AM_V%Qn))
-! allocate(AM_V%Qtau_e(N_SLS,AM_V%Qn))
-!
-! read(27) AM_V%Qr
-! read(27) AM_V%Qmu
-! read(27) AM_V%Qtau_e
-! close(27)
-! endif
-!
-! ! Synch up after the Read
-! call MPI_BCAST(AM_V%QT_c_source,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
-! call MPI_BCAST(tau_s,N_SLS,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
-! call MPI_BCAST(AM_V%Qn,1,MPI_INTEGER,0,MPI_COMM_WORLD,ier)
-!
-! if(myrank /= 0) then
-! allocate(AM_V%Qr(AM_V%Qn))
-! allocate(AM_V%Qmu(AM_V%Qn))
-! allocate(AM_V%Qtau_e(N_SLS,AM_V%Qn))
-! endif
-!
-! call MPI_BCAST(AM_V%Qr,AM_V%Qn,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
-! call MPI_BCAST(AM_V%Qmu,AM_V%Qn,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
-! call MPI_BCAST(AM_V%Qtau_e,AM_V%Qn*N_SLS,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
-!
-! scale_t = ONE/dsqrt(PI*GRAV*RHOAV)
-!
-! ! Scale the Attenuation Values
-! tau_s(:) = tau_s(:) / scale_t
-! AM_V%Qtau_e(:,:) = AM_V%Qtau_e(:,:) / scale_t
-! AM_V%QT_c_source = 1000.0d0 / AM_V%QT_c_source / scale_t
-! AM_V%Qr(:) = AM_V%Qr(:) / R_EARTH
-!
-! allocate(AM_V%Qsf(AM_V%Qn))
-! allocate(AM_V%Qomsb(AM_V%Qn))
-! allocate(AM_V%Qfc(N_SLS,AM_V%Qn))
-!
-! allocate(AM_V%Qsf2(AM_V%Qn))
-! allocate(AM_V%Qomsb2(AM_V%Qn))
-! allocate(AM_V%Qfc2(N_SLS,AM_V%Qn))
-!
-! allocate(AM_V%interval_Q(AM_V%Qn))
-!
-! allocate(Qfctmp(AM_V%Qn))
-! allocate(Qfc2tmp(AM_V%Qn))
-!
-! do i = 1,AM_V%Qn
-! if(AM_V%Qmu(i) == 0.0d0) then
-! AM_V%Qomsb(i) = 0.0d0
-! AM_V%Qfc(:,i) = 0.0d0
-! AM_V%Qsf(i) = 0.0d0
-! else
-! call attenuation_property_values(tau_s, AM_V%Qtau_e(:,i), AM_V%Qfc(:,i), AM_V%Qomsb(i))
-! call attenuation_scale_factor(myrank, AM_V%QT_c_source, AM_V%Qtau_e(:,i), tau_s, AM_V%Qmu(i), AM_V%Qsf(i))
-! endif
-! enddo
-!
-! ! Determine the Spline Coefficients or Second Derivatives
-! call pspline_construction(AM_V%Qr, AM_V%Qsf, AM_V%Qn, Qp1, Qpn, AM_V%Qsf2, AM_V%interval_Q)
-! call pspline_construction(AM_V%Qr, AM_V%Qomsb, AM_V%Qn, Qp1, Qpn, AM_V%Qomsb2, AM_V%interval_Q)
-! do i = 1,N_SLS
-!! copy the sub-arrays to temporary arrays to avoid a warning by some compilers
-!! about temporary arrays being created automatically when using this expression
-!! directly in the call to the subroutine
-! Qfctmp(:) = AM_V%Qfc(i,:)
-! Qfc2tmp(:) = AM_V%Qfc2(i,:)
-! call pspline_construction(AM_V%Qr, Qfctmp, AM_V%Qn, Qp1, Qpn, Qfc2tmp, AM_V%interval_Q)
-!! copy the arrays back to the sub-arrays, since these sub-arrays are used
-!! as input and output
-! AM_V%Qfc(i,:) = Qfctmp(:)
-! AM_V%Qfc2(i,:) = Qfc2tmp(:)
-! enddo
-!
-! radius = 0.0d0
-! rmax = nint(TABLE_ATTENUATION)
-! do i = 1,rmax
-! call attenuation_lookup_value(i, radius)
-! call pspline_evaluation(AM_V%Qr, AM_V%Qsf, AM_V%Qsf2, AM_V%Qn, radius, scale_factor(1,1,1,i), AM_V%interval_Q)
-! call pspline_evaluation(AM_V%Qr, AM_V%Qomsb, AM_V%Qomsb2, AM_V%Qn, radius, one_minus_sum_beta(1,1,1,i), AM_V%interval_Q)
-! do j = 1,N_SLS
-! Qfctmp = AM_V%Qfc(j,:)
-! Qfc2tmp = AM_V%Qfc2(j,:)
-! call pspline_evaluation(AM_V%Qr, Qfctmp, Qfc2tmp, AM_V%Qn, radius, fctmp, AM_V%interval_Q)
-! factor_common(j,1,1,1,i) = fctmp
-! enddo
-! enddo
-! do i = rmax+1,NRAD_ATTENUATION
-! scale_factor(1,1,1,i) = scale_factor(1,1,1,rmax)
-! one_minus_sum_beta(1,1,1,i) = one_minus_sum_beta(1,1,1,rmax)
-! factor_common(1,1,1,1,i) = factor_common(1,1,1,1,rmax)
-! factor_common(2,1,1,1,i) = factor_common(2,1,1,1,rmax)
-! factor_common(3,1,1,1,i) = factor_common(3,1,1,1,rmax)
-! enddo
-!
-! deallocate(AM_V%Qfc2)
-! deallocate(AM_V%Qsf2)
-! deallocate(AM_V%Qomsb2)
-! deallocate(AM_V%Qfc)
-! deallocate(AM_V%Qsf)
-! deallocate(AM_V%Qomsb)
-! deallocate(AM_V%Qtau_e)
-! deallocate(Qfctmp)
-! deallocate(Qfc2tmp)
-!
-! call MPI_BARRIER(MPI_COMM_WORLD, ier)
-!
-! end subroutine get_attenuation_model_1D
-!
-!
-!-------------------------------------------------------------------------------------------------
-!
-!
-!-------------------------------------------------------------------------------------------------
-!
-! not used anymore...
-!
-! Piecewise Continuous Splines
-! - Added Steps which describes the discontinuities
-! - Steps must be repeats in the dependent variable, X
-! - Derivates at the steps are computed using the point
-! at the derivate and the closest point within that piece
-! - A point lying directly on the discontinuity will recieve the
-! value of the first or smallest piece in terms of X
-! - Beginning and Ending points of the Function become beginning
-! and ending points of the first and last splines
-! - A Step with a value of zero is undefined
-! - Works with functions with steps or no steps
-! See the comment below about the ScS bug
-! subroutine pspline_evaluation(xa, ya, y2a, n, x, y, steps)
-!
-! implicit none
-!
-! integer n
-! double precision xa(n),ya(n),y2a(n)
-! integer steps(n)
-! double precision x, y
-!
-! integer i, l, n1, n2
-!
-! do i = 1,n-1,1
-! if(steps(i+1) == 0) return
-! if(x >= xa(steps(i)) .and. x <= xa(steps(i+1))) then
-! call pspline_piece(i,n1,n2,l,n,steps)
-! call spline_evaluation(xa(n1), ya(n1), y2a(n1), l, x, y)
-!! return <-- Commented out to fix ScS bug
-! endif
-! enddo
-!
-! end subroutine pspline_evaluation
-!
-!
-!-------------------------------------------------------------------------------------------------
-!
-! not used anymore...
-!
-! subroutine pspline_piece(i,n1,n2,l,n,s)
-!
-! implicit none
-!
-! integer i, n1, n2, l, n, s(n)
-! n1 = s(i)+1
-! if(i == 1) n1 = s(i)
-! n2 = s(i+1)
-! l = n2 - n1 + 1
-!
-! end subroutine pspline_piece
-!
-!
-!-------------------------------------------------------------------------------------------------
-!
-! not used anymore...
-!
-! subroutine pspline_construction(x, y, n, yp1, ypn, y2, steps)
-!
-! implicit none
-!
-! integer n
-! double precision x(n),y(n),y2(n)
-! double precision yp1, ypn
-! integer steps(n)
-!
-! integer i,r, l, n1,n2
-!
-! steps(:) = 0
-!
-! ! Find steps in x, defining pieces
-! steps(1) = 1
-! r = 2
-! do i = 2,n
-! if(x(i) == x(i-1)) then
-! steps(r) = i-1
-! r = r + 1
-! endif
-! enddo
-! steps(r) = n
-!
-! ! Run spline for each piece
-! do i = 1,r-1
-! call pspline_piece(i,n1,n2,l,n,steps)
-! ! Determine the First Derivates at Begin/End Points
-! yp1 = ( y(n1+1) - y(n1) ) / ( x(n1+1) - x(n1))
-! ypn = ( y(n2) - y(n2-1) ) / ( x(n2) - x(n2-1))
-! call spline_construction(x(n1),y(n1),l,yp1,ypn,y2(n1))
-! enddo
-!
-! end subroutine pspline_construction
-!
-!
-!-------------------------------------------------------------------------------------------------
-!
-!
-! not used anymore...
-!
-! subroutine attenuation_lookup_value(i, r)
-!
-! implicit none
-!
-! include 'constants.h'
-!
-! integer i
-! double precision r
-!
-! r = dble(i) / TABLE_ATTENUATION
-!
-! end subroutine attenuation_lookup_value
-!
-!
-!-------------------------------------------------------------------------------------------------
-!
-! not used anymore...
-!
-! subroutine attenuation_save_arrays(prname, iregion_code, AM_V)
-!
-! implicit none
-!
-! include 'mpif.h'
-! include 'constants.h'
-!
-!! model_attenuation_variables
-! type model_attenuation_variables
-! sequence
-! double precision min_period, max_period
-! double precision :: QT_c_source ! Source Frequency
-! double precision, dimension(:), pointer :: Qtau_s ! tau_sigma
-! double precision, dimension(:), pointer :: QrDisc ! Discontinutitues Defined
-! double precision, dimension(:), pointer :: Qr ! Radius
-! integer, dimension(:), pointer :: interval_Q ! Steps
-! double precision, dimension(:), pointer :: Qmu ! Shear Attenuation
-! double precision, dimension(:,:), pointer :: Qtau_e ! tau_epsilon
-! double precision, dimension(:), pointer :: Qomsb, Qomsb2 ! one_minus_sum_beta
-! double precision, dimension(:,:), pointer :: Qfc, Qfc2 ! factor_common
-! double precision, dimension(:), pointer :: Qsf, Qsf2 ! scale_factor
-! integer, dimension(:), pointer :: Qrmin ! Max and Mins of idoubling
-! integer, dimension(:), pointer :: Qrmax ! Max and Mins of idoubling
-! integer :: Qn ! Number of points
-! end type model_attenuation_variables
-!
-! type (model_attenuation_variables) AM_V
-!! model_attenuation_variables
-!
-! integer iregion_code
-! character(len=150) prname
-! integer ier
-! integer myrank
-! integer, save :: first_time_called = 1
-!
-! call MPI_COMM_RANK(MPI_COMM_WORLD, myrank, ier)
-! if(myrank == 0 .AND. iregion_code == IREGION_CRUST_MANTLE .AND. first_time_called == 1) then
-! first_time_called = 0
-! open(unit=27,file=prname(1:len_trim(prname))//'1D_Q.bin',status='unknown',form='unformatted')
-! write(27) AM_V%QT_c_source
-! write(27) AM_V%Qtau_s
-! write(27) AM_V%Qn
-! write(27) AM_V%Qr
-! write(27) AM_V%Qmu
-! write(27) AM_V%Qtau_e
-! close(27)
-! endif
-!
-! end subroutine attenuation_save_arrays
-!
-!
-!-------------------------------------------------------------------------------------------------
-!
-! not used anymore...
-!
-! subroutine get_attenuation_index(iflag, radius, index, inner_core, AM_V)
-!
-! implicit none
-!
-! include 'constants.h'
-!
-!! model_attenuation_variables
-! type model_attenuation_variables
-! sequence
-! double precision min_period, max_period
-! double precision :: QT_c_source ! Source Frequency
-! double precision, dimension(:), pointer :: Qtau_s ! tau_sigma
-! double precision, dimension(:), pointer :: QrDisc ! Discontinutitues Defined
-! double precision, dimension(:), pointer :: Qr ! Radius
-! integer, dimension(:), pointer :: interval_Q ! Steps
-! double precision, dimension(:), pointer :: Qmu ! Shear Attenuation
-! double precision, dimension(:,:), pointer :: Qtau_e ! tau_epsilon
-! double precision, dimension(:), pointer :: Qomsb, Qomsb2 ! one_minus_sum_beta
-! double precision, dimension(:,:), pointer :: Qfc, Qfc2 ! factor_common
-! double precision, dimension(:), pointer :: Qsf, Qsf2 ! scale_factor
-! integer, dimension(:), pointer :: Qrmin ! Max and Mins of idoubling
-! integer, dimension(:), pointer :: Qrmax ! Max and Mins of idoubling
-! integer :: Qn ! Number of points
-! end type model_attenuation_variables
-!
-! type (model_attenuation_variables) AM_V
-!! model_attenuation_variables
-!
-! integer iflag, iregion, index
-! double precision radius
-!
-! ! Inner Core or not
-! logical inner_core
-!
-! index = nint(radius * TABLE_ATTENUATION)
-!
-!!! DK DK this seems incorrect and is difficult to read anyway
-!!! DK DK therefore let me rewrite it better
-!! if(inner_core) then
-!! if(iflag >= IFLAG_INNER_CORE_NORMAL) then
-!! iregion = IREGION_ATTENUATION_INNER_CORE
-!! else if(iflag >= IFLAG_OUTER_CORE_NORMAL) then
-!! iregion = 6
-!! endif
-!! else
-!! if(iflag >= IFLAG_MANTLE_NORMAL) then
-!! iregion = IREGION_ATTENUATION_CMB_670
-!! else if(iflag == IFLAG_670_220) then
-!! iregion = IREGION_ATTENUATION_670_220
-!! else if(iflag <= IFLAG_220_80) then
-!! iregion = IREGION_ATTENUATION_220_80
-!! else
-!! iregion = IREGION_ATTENUATION_80_SURFACE
-!! endif
-!! endif
-! if(inner_core) then
-!
-! if(iflag == IFLAG_INNER_CORE_NORMAL .or. iflag == IFLAG_MIDDLE_CENTRAL_CUBE .or. &
-! iflag == IFLAG_BOTTOM_CENTRAL_CUBE .or. iflag == IFLAG_TOP_CENTRAL_CUBE .or. &
-! iflag == IFLAG_IN_FICTITIOUS_CUBE) then
-! iregion = IREGION_ATTENUATION_INNER_CORE
-! else
-!! this is fictitious for the outer core, which has no Qmu attenuation since it is fluid
-!! iregion = IREGION_ATTENUATION_80_SURFACE + 1
-! iregion = IREGION_ATTENUATION_UNDEFINED
-! endif
-!
-! else
-!
-! if(iflag == IFLAG_MANTLE_NORMAL) then
-! iregion = IREGION_ATTENUATION_CMB_670
-! else if(iflag == IFLAG_670_220) then
-! iregion = IREGION_ATTENUATION_670_220
-! else if(iflag == IFLAG_220_80) then
-! iregion = IREGION_ATTENUATION_220_80
-! else if(iflag == IFLAG_CRUST .or. iflag == IFLAG_80_MOHO) then
-! iregion = IREGION_ATTENUATION_80_SURFACE
-! else
-!! this is fictitious for the outer core, which has no Qmu attenuation since it is fluid
-!! iregion = IREGION_ATTENUATION_80_SURFACE + 1
-! iregion = IREGION_ATTENUATION_UNDEFINED
-! endif
-!
-! endif
-!
-!! Clamp regions
-! if(index < AM_V%Qrmin(iregion)) index = AM_V%Qrmin(iregion)
-! if(index > AM_V%Qrmax(iregion)) index = AM_V%Qrmax(iregion)
-!
-! end subroutine get_attenuation_index
-!
-!
-!-------------------------------------------------------------------------------------------------
-!
-! not used anymore...
-!
-! subroutine set_attenuation_regions_1D(RICB, RCMB, R670, R220, R80, AM_V)
-!
-! implicit none
-!
-! include 'constants.h'
-!
-!! model_attenuation_variables
-! type model_attenuation_variables
-! sequence
-! double precision min_period, max_period
-! double precision :: QT_c_source ! Source Frequency
-! double precision, dimension(:), pointer :: Qtau_s ! tau_sigma
-! double precision, dimension(:), pointer :: QrDisc ! Discontinutitues Defined
-! double precision, dimension(:), pointer :: Qr ! Radius
-! integer, dimension(:), pointer :: interval_Q ! Steps
-! double precision, dimension(:), pointer :: Qmu ! Shear Attenuation
-! double precision, dimension(:,:), pointer :: Qtau_e ! tau_epsilon
-! double precision, dimension(:), pointer :: Qomsb, Qomsb2 ! one_minus_sum_beta
-! double precision, dimension(:,:), pointer :: Qfc, Qfc2 ! factor_common
-! double precision, dimension(:), pointer :: Qsf, Qsf2 ! scale_factor
-! integer, dimension(:), pointer :: Qrmin ! Max and Mins of idoubling
-! integer, dimension(:), pointer :: Qrmax ! Max and Mins of idoubling
-! integer :: Qn ! Number of points
-! end type model_attenuation_variables
-!
-! type (model_attenuation_variables) AM_V
-!! model_attenuation_variables
-!
-! double precision RICB, RCMB, R670, R220, R80
-! integer i
-!
-! allocate(AM_V%Qrmin(6))
-! allocate(AM_V%Qrmax(6))
-! allocate(AM_V%QrDisc(5))
-!
-! AM_V%QrDisc(1) = RICB
-! AM_V%QrDisc(2) = RCMB
-! AM_V%QrDisc(3) = R670
-! AM_V%QrDisc(4) = R220
-! AM_V%QrDisc(5) = R80
-!
-! ! INNER CORE
-! AM_V%Qrmin(IREGION_ATTENUATION_INNER_CORE) = 1 ! Center of the Earth
-! i = nint(RICB / 100.d0) ! === BOUNDARY === INNER CORE / OUTER CORE
-! AM_V%Qrmax(IREGION_ATTENUATION_INNER_CORE) = i - 1 ! Inner Core Boundary (Inner)
-!
-! ! OUTER_CORE
-! AM_V%Qrmin(6) = i ! Inner Core Boundary (Outer)
-! i = nint(RCMB / 100.d0) ! === BOUNDARY === INNER CORE / OUTER CORE
-! AM_V%Qrmax(6) = i - 1
-!
-! ! LOWER MANTLE
-! AM_V%Qrmin(IREGION_ATTENUATION_CMB_670) = i
-! i = nint(R670 / 100.d0) ! === BOUNDARY === 670 km
-! AM_V%Qrmax(IREGION_ATTENUATION_CMB_670) = i - 1
-!
-! ! UPPER MANTLE
-! AM_V%Qrmin(IREGION_ATTENUATION_670_220) = i
-! i = nint(R220 / 100.d0) ! === BOUNDARY === 220 km
-! AM_V%Qrmax(IREGION_ATTENUATION_670_220) = i - 1
-!
-! ! MANTLE ISH LITHOSPHERE
-! AM_V%Qrmin(IREGION_ATTENUATION_220_80) = i
-! i = nint(R80 / 100.d0) ! === BOUNDARY === 80 km
-! AM_V%Qrmax(IREGION_ATTENUATION_220_80) = i - 1
-!
-! ! CRUST ISH LITHOSPHERE
-! AM_V%Qrmin(IREGION_ATTENUATION_80_SURFACE) = i
-! AM_V%Qrmax(IREGION_ATTENUATION_80_SURFACE) = NRAD_ATTENUATION
-!
-! end subroutine set_attenuation_regions_1D
-!
-
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_forward_arrays.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_forward_arrays.f90 2013-07-02 01:23:11 UTC (rev 22484)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/read_forward_arrays.f90 2013-07-02 02:26:51 UTC (rev 22485)
@@ -33,6 +33,7 @@
use specfem_par_crustmantle
use specfem_par_innercore
use specfem_par_outercore
+
implicit none
! local parameters
@@ -150,7 +151,7 @@
b_B_array_rotation = 0._CUSTOM_REAL
endif
- if(ATTENUATION_VAL) then
+ if (ATTENUATION_VAL) then
b_R_xx_crust_mantle = 0._CUSTOM_REAL
b_R_yy_crust_mantle = 0._CUSTOM_REAL
b_R_xy_crust_mantle = 0._CUSTOM_REAL
@@ -179,9 +180,10 @@
use specfem_par_crustmantle
use specfem_par_innercore
use specfem_par_outercore
+
implicit none
- !local parameters
+ ! local parameters
integer :: ier
character(len=150) outputname
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_element.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_element.F90 2013-07-02 01:23:11 UTC (rev 22484)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_element.F90 2013-07-02 02:26:51 UTC (rev 22485)
@@ -66,13 +66,15 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPECMAX_ISO_MANTLE) :: &
kappavstore,muvstore
+ ! variable sized array variables
+ integer :: vx,vy,vz,vnspec
+
! attenuation
! memory variables for attenuation
! memory variables R_ij are stored at the local rather than global level
! to allow for optimization of cache access by compiler
real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: R_memory
- integer :: vx,vy,vz,vnspec
real(kind=CUSTOM_REAL), dimension(vx, vy, vz, vnspec) :: one_minus_sum_beta
! gravity
@@ -187,7 +189,7 @@
if(ATTENUATION_VAL) mul = mul * one_minus_sum_beta_use
lambdalplus2mul = kappal + FOUR_THIRDS * mul
- lambdal = lambdalplus2mul - 2.0*mul
+ lambdal = lambdalplus2mul - 2.0_CUSTOM_REAL*mul
! compute stress sigma
sigma_xx = lambdalplus2mul*duxdxl + lambdal*duydyl_plus_duzdzl
@@ -337,8 +339,6 @@
end subroutine compute_element_iso
-
-
!
!--------------------------------------------------------------------------------------------------
!
@@ -392,10 +392,11 @@
! to allow for optimization of cache access by compiler
real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: R_memory
- integer vx,vy,vz,vnspec
+ ! variable sized array variables
+ integer :: vx,vy,vz,vnspec
! [alpha,beta,gamma]val reduced to N_SLS to N_SLS*NUM_NODES
- real(kind=CUSTOM_REAL), dimension(vx, vy, vz, vnspec) :: one_minus_sum_beta
+ real(kind=CUSTOM_REAL), dimension(vx,vy,vz,vnspec) :: one_minus_sum_beta
! gravity
double precision, dimension(NRAD_GRAVITY) :: minus_gravity_table,density_table,minus_deriv_gravity_table
@@ -607,87 +608,87 @@
! way 2: reordering operations to facilitate compilation, avoiding divisions, using locality for temporary values
c11 = rhovphsq*sinphifour &
- + 2.0*cosphisq*sinphisq* &
+ + 2.0_CUSTOM_REAL*cosphisq*sinphisq* &
( rhovphsq*costhetasq + sinthetasq*(eta_aniso*rhovphsq + two_rhovsvsq - two_eta_aniso*rhovsvsq) ) &
+ cosphifour*(rhovphsq*costhetafour &
- + 2.0*costhetasq*sinthetasq*(eta_aniso*rhovphsq + two_rhovsvsq - two_eta_aniso*rhovsvsq) &
+ + 2.0_CUSTOM_REAL*costhetasq*sinthetasq*(eta_aniso*rhovphsq + two_rhovsvsq - two_eta_aniso*rhovsvsq) &
+ rhovpvsq*sinthetafour)
- c12 = 0.25*costhetasq*(rhovphsq - two_rhovshsq)*(3.0 + cosfourphi) &
+ c12 = 0.25_CUSTOM_REAL*costhetasq*(rhovphsq - two_rhovshsq)*(3.0_CUSTOM_REAL + cosfourphi) &
- four_rhovshsq*cosphisq*costhetasq*sinphisq &
- + 0.03125*rhovphsq*sintwophisq*(11.0 + cosfourtheta + 4.0*costwotheta) &
+ + 0.03125_CUSTOM_REAL*rhovphsq*sintwophisq*(11.0_CUSTOM_REAL + cosfourtheta + 4.0*costwotheta) &
+ eta_aniso*sinthetasq*(rhovphsq - two_rhovsvsq) &
- *(cosphifour + sinphifour + 2.0*cosphisq*costhetasq*sinphisq) &
+ *(cosphifour + sinphifour + 2.0_CUSTOM_REAL*cosphisq*costhetasq*sinphisq) &
+ rhovpvsq*cosphisq*sinphisq*sinthetafour &
- rhovsvsq*sintwophisq*sinthetafour
- c13 = 0.125*cosphisq*(rhovphsq + six_eta_aniso*rhovphsq + rhovpvsq - four_rhovsvsq &
- - 12.0*eta_aniso*rhovsvsq + cosfourtheta*templ1) &
+ c13 = 0.125_CUSTOM_REAL*cosphisq*(rhovphsq + six_eta_aniso*rhovphsq + rhovpvsq - four_rhovsvsq &
+ - 12.0_CUSTOM_REAL*eta_aniso*rhovsvsq + cosfourtheta*templ1) &
+ sinphisq*(eta_aniso*costhetasq*(rhovphsq - two_rhovsvsq) + sinthetasq*(rhovphsq - two_rhovshsq))
! uses temporary templ1 from c13
c15 = cosphi*costheta*sintheta* &
- ( 0.5*cosphisq* (rhovpvsq - rhovphsq + costwotheta*templ1) &
+ ( 0.5_CUSTOM_REAL*cosphisq* (rhovpvsq - rhovphsq + costwotheta*templ1) &
+ etaminone*sinphisq*(rhovphsq - two_rhovsvsq))
c14 = costheta*sinphi*sintheta* &
- ( 0.5*cosphisq*(templ2_cos + four_rhovshsq - four_rhovsvsq) &
- + sinphisq*(etaminone*rhovphsq + 2.0*(rhovshsq - eta_aniso*rhovsvsq)) )
+ ( 0.5_CUSTOM_REAL*cosphisq*(templ2_cos + four_rhovshsq - four_rhovsvsq) &
+ + sinphisq*(etaminone*rhovphsq + 2.0_CUSTOM_REAL*(rhovshsq - eta_aniso*rhovsvsq)) )
! uses temporary templ2_cos from c14
- c16 = 0.5*cosphi*sinphi*sinthetasq* &
+ c16 = 0.5_CUSTOM_REAL*cosphi*sinphi*sinthetasq* &
( cosphisq*templ2_cos &
- + 2.0*etaminone*sinphisq*(rhovphsq - two_rhovsvsq) )
+ + 2.0_CUSTOM_REAL*etaminone*sinphisq*(rhovphsq - two_rhovsvsq) )
- c22 = rhovphsq*cosphifour + 2.0*cosphisq*sinphisq* &
+ c22 = rhovphsq*cosphifour + 2.0_CUSTOM_REAL*cosphisq*sinphisq* &
(rhovphsq*costhetasq + sinthetasq*(eta_aniso*rhovphsq + two_rhovsvsq - two_eta_aniso*rhovsvsq)) &
+ sinphifour* &
- (rhovphsq*costhetafour + 2.0*costhetasq*sinthetasq*(eta_aniso*rhovphsq &
+ (rhovphsq*costhetafour + 2.0_CUSTOM_REAL*costhetasq*sinthetasq*(eta_aniso*rhovphsq &
+ two_rhovsvsq - two_eta_aniso*rhovsvsq) + rhovpvsq*sinthetafour)
! uses temporary templ1 from c13
- c23 = 0.125*sinphisq*(rhovphsq + six_eta_aniso*rhovphsq &
- + rhovpvsq - four_rhovsvsq - 12.0*eta_aniso*rhovsvsq + cosfourtheta*templ1) &
+ c23 = 0.125_CUSTOM_REAL*sinphisq*(rhovphsq + six_eta_aniso*rhovphsq &
+ + rhovpvsq - four_rhovsvsq - 12.0_CUSTOM_REAL*eta_aniso*rhovsvsq + cosfourtheta*templ1) &
+ cosphisq*(eta_aniso*costhetasq*(rhovphsq - two_rhovsvsq) + sinthetasq*(rhovphsq - two_rhovshsq))
! uses temporary templ1 from c13
c24 = costheta*sinphi*sintheta* &
( etaminone*cosphisq*(rhovphsq - two_rhovsvsq) &
- + 0.5*sinphisq*(rhovpvsq - rhovphsq + costwotheta*templ1) )
+ + 0.5_CUSTOM_REAL*sinphisq*(rhovpvsq - rhovphsq + costwotheta*templ1) )
! uses temporary templ2_cos from c14
c25 = cosphi*costheta*sintheta* &
- ( cosphisq*(etaminone*rhovphsq + 2.0*(rhovshsq - eta_aniso*rhovsvsq)) &
- + 0.5*sinphisq*(templ2_cos + four_rhovshsq - four_rhovsvsq) )
+ ( cosphisq*(etaminone*rhovphsq + 2.0_CUSTOM_REAL*(rhovshsq - eta_aniso*rhovsvsq)) &
+ + 0.5_CUSTOM_REAL*sinphisq*(templ2_cos + four_rhovshsq - four_rhovsvsq) )
! uses temporary templ2_cos from c14
- c26 = 0.5*cosphi*sinphi*sinthetasq* &
- ( 2.0*etaminone*cosphisq*(rhovphsq - two_rhovsvsq) &
+ c26 = 0.5_CUSTOM_REAL*cosphi*sinphi*sinthetasq* &
+ ( 2.0_CUSTOM_REAL*etaminone*cosphisq*(rhovphsq - two_rhovsvsq) &
+ sinphisq*templ2_cos )
c33 = rhovpvsq*costhetafour &
- + 2.0*costhetasq*sinthetasq*(two_rhovsvsq + eta_aniso*(rhovphsq - two_rhovsvsq)) &
+ + 2.0_CUSTOM_REAL*costhetasq*sinthetasq*(two_rhovsvsq + eta_aniso*(rhovphsq - two_rhovsvsq)) &
+ rhovphsq*sinthetafour
! uses temporary templ1_cos from c13
- c34 = - 0.25*sinphi*sintwotheta*templ1_cos
+ c34 = - 0.25_CUSTOM_REAL*sinphi*sintwotheta*templ1_cos
! uses temporary templ1_cos from c34
- c35 = - 0.25*cosphi*sintwotheta*templ1_cos
+ c35 = - 0.25_CUSTOM_REAL*cosphi*sintwotheta*templ1_cos
! uses temporary templ1_cos from c34
- c36 = - 0.25*sintwophi*sinthetasq*(templ1_cos - four_rhovshsq + four_rhovsvsq)
+ c36 = - 0.25_CUSTOM_REAL*sintwophi*sinthetasq*(templ1_cos - four_rhovshsq + four_rhovsvsq)
c44 = cosphisq*(rhovsvsq*costhetasq + rhovshsq*sinthetasq) &
+ sinphisq*(rhovsvsq*costwothetasq + costhetasq*sinthetasq*templ3)
! uses temporary templ3 from c44
c46 = - cosphi*costheta*sintheta* &
- ( cosphisq*(rhovshsq - rhovsvsq) - 0.5*sinphisq*templ3_cos )
+ ( cosphisq*(rhovshsq - rhovsvsq) - 0.5_CUSTOM_REAL*sinphisq*templ3_cos )
! uses templ3 from c46
- c45 = 0.25*sintwophi*sinthetasq* &
- (templ3_two + costwotheta*(rhovphsq + rhovpvsq - two_eta_aniso*rhovphsq + 4.0*etaminone*rhovsvsq))
+ c45 = 0.25_CUSTOM_REAL*sintwophi*sinthetasq* &
+ (templ3_two + costwotheta*(rhovphsq + rhovpvsq - two_eta_aniso*rhovphsq + 4.0_CUSTOM_REAL*etaminone*rhovsvsq))
c55 = sinphisq*(rhovsvsq*costhetasq + rhovshsq*sinthetasq) &
+ cosphisq*(rhovsvsq*costwothetasq &
@@ -695,17 +696,18 @@
! uses temporary templ3_cos from c46
c56 = costheta*sinphi*sintheta* &
- ( 0.5*cosphisq*templ3_cos + sinphisq*(rhovsvsq - rhovshsq) )
+ ( 0.5_CUSTOM_REAL*cosphisq*templ3_cos + sinphisq*(rhovsvsq - rhovshsq) )
c66 = rhovshsq*costwophisq*costhetasq &
- - 2.0*cosphisq*costhetasq*sinphisq*(rhovphsq - two_rhovshsq) &
- + 0.03125*rhovphsq*sintwophisq*(11.0 + 4.0*costwotheta + cosfourtheta) &
- - 0.125*rhovsvsq*sinthetasq*( -6.0 - 2.0*costwotheta - 2.0*cosfourphi &
- + cos(4.0*phi - 2.0*theta) + cos(2.0*(2.0*phi + theta)) ) &
+ - 2.0_CUSTOM_REAL*cosphisq*costhetasq*sinphisq*(rhovphsq - two_rhovshsq) &
+ + 0.03125_CUSTOM_REAL*rhovphsq*sintwophisq*(11.0_CUSTOM_REAL + 4.0_CUSTOM_REAL*costwotheta + cosfourtheta) &
+ - 0.125_CUSTOM_REAL*rhovsvsq*sinthetasq* &
+ ( -6.0_CUSTOM_REAL - 2.0_CUSTOM_REAL*costwotheta - 2.0_CUSTOM_REAL*cosfourphi &
+ + cos(4.0_CUSTOM_REAL*phi - 2.0_CUSTOM_REAL*theta) &
+ + cos(2.0_CUSTOM_REAL*(2.0_CUSTOM_REAL*phi + theta)) ) &
+ rhovpvsq*cosphisq*sinphisq*sinthetafour &
- - 0.5*eta_aniso*sintwophisq*sinthetafour*(rhovphsq - two_rhovsvsq)
+ - 0.5_CUSTOM_REAL*eta_aniso*sintwophisq*sinthetafour*(rhovphsq - two_rhovsvsq)
-
! general expression of stress tensor for full Cijkl with 21 coefficients
sigma_xx = c11*duxdxl + c16*duxdyl_plus_duydxl + c12*duydyl + &
c15*duzdxl_plus_duxdzl + c14*duzdyl_plus_duydzl + c13*duzdzl
@@ -863,7 +865,6 @@
enddo ! NGLLY
enddo ! NGLLZ
-
end subroutine compute_element_tiso
!
@@ -921,10 +922,11 @@
! to allow for optimization of cache access by compiler
real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: R_memory
- integer vx,vy,vz,vnspec
+ ! variable sized array variables
+ integer :: vx,vy,vz,vnspec
! [alpha,beta,gamma]val reduced to N_SLS to N_SLS*NUM_NODES
- real(kind=CUSTOM_REAL), dimension(vx, vy, vz, vnspec) :: one_minus_sum_beta
+ real(kind=CUSTOM_REAL), dimension(vx,vy,vz,vnspec) :: one_minus_sum_beta
! gravity
double precision, dimension(NRAD_GRAVITY) :: minus_gravity_table,density_table,minus_deriv_gravity_table
@@ -1229,8 +1231,8 @@
!
- subroutine compute_element_att_stress( R_memory_loc, &
- sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz)
+ subroutine compute_element_att_stress(R_memory_loc, &
+ sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz)
implicit none
@@ -1304,8 +1306,10 @@
real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: R_memory
+ ! variable sized array variables
integer :: vx,vy,vz,vnspec
- real(kind=CUSTOM_REAL), dimension(N_SLS, vx, vy, vz, vnspec) :: factor_common
+
+ real(kind=CUSTOM_REAL), dimension(N_SLS,vx,vy,vz,vnspec) :: factor_common
real(kind=CUSTOM_REAL), dimension(N_SLS) :: alphaval,betaval,gammaval
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPECMAX_ANISO_MANTLE) :: c44store
@@ -1315,7 +1319,7 @@
real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc
! local parameters
- real(kind=CUSTOM_REAL), dimension(NGLLX, NGLLY, NGLLZ) :: factor_common_c44_muv
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: factor_common_c44_muv
integer :: i_SLS
integer :: i_memory
@@ -1401,8 +1405,10 @@
real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ATTENUATION) :: R_memory
+ ! variable sized array variables
integer :: vx,vy,vz,vnspec
- real(kind=CUSTOM_REAL), dimension(N_SLS, vx, vy, vz, vnspec) :: factor_common
+
+ real(kind=CUSTOM_REAL), dimension(N_SLS,vx,vy,vz,vnspec) :: factor_common
real(kind=CUSTOM_REAL), dimension(N_SLS) :: alphaval,betaval,gammaval
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE) :: muvstore
@@ -1411,7 +1417,7 @@
real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc_nplus1
! local parameters
- real(kind=CUSTOM_REAL), dimension(NGLLX, NGLLY, NGLLZ) :: factor_common_use
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: factor_common_use
integer :: i_SLS
@@ -1831,18 +1837,20 @@
include "constants.h"
-! include values created by the mesher
-! done for performance only using static allocation to allow for loop unrolling
+ ! include values created by the mesher
+ ! done for performance only using static allocation to allow for loop unrolling
include "OUTPUT_FILES/values_from_mesher.h"
- integer ispec,NSPEC,NGLOB
+ ! element id
+ integer :: ispec,i,j,k
+ integer NSPEC,NGLOB
+
! array with derivatives of Lagrange polynomials and precalculated products
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx
real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLY) :: hprime_yy
real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ) :: hprime_zz
-
! arrays with mesh parameters per slice
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: ibool
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC) :: &
@@ -1852,9 +1860,9 @@
real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ) :: epsilondev_loc
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: eps_trace_over_3_loc
-!local parameters
+! local parameters
integer iglob
- integer i,j,k,l
+ integer l
real(kind=CUSTOM_REAL) tempx1l,tempx2l,tempx3l
real(kind=CUSTOM_REAL) tempy1l,tempy2l,tempy3l
@@ -1869,9 +1877,6 @@
real(kind=CUSTOM_REAL) duxdxl_plus_duydyl,duxdxl_plus_duzdzl,duydyl_plus_duzdzl
real(kind=CUSTOM_REAL) duxdyl_plus_duydxl,duzdxl_plus_duxdzl,duzdyl_plus_duydzl
-
-
-
do k=1,NGLLZ
do j=1,NGLLY
do i=1,NGLLX
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_stacey_crust_mantle.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_stacey_crust_mantle.f90 2013-07-02 01:23:11 UTC (rev 22484)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_stacey_crust_mantle.f90 2013-07-02 02:26:51 UTC (rev 22485)
@@ -47,7 +47,6 @@
absorb_xmin_crust_mantle,absorb_xmax_crust_mantle, &
absorb_ymin_crust_mantle,absorb_ymax_crust_mantle)
-
implicit none
include "constants.h"
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/get_attenuation.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/get_attenuation.f90 2013-07-02 01:23:11 UTC (rev 22484)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/get_attenuation.f90 2013-07-02 02:26:51 UTC (rev 22485)
@@ -25,7 +25,6 @@
!
!=====================================================================
-
subroutine get_attenuation_model_3D(myrank, prname, one_minus_sum_beta, &
factor_common, scale_factor, tau_s, vnspec)
@@ -33,16 +32,31 @@
include 'constants.h'
- integer myrank, vnspec
- character(len=150) prname
+ integer :: myrank
+
+ integer :: vnspec
+
+!! DK DK to Daniel, Jul 2013
+!! DK DK to Daniel, Jul 2013
+!! DK DK to Daniel, Jul 2013
+!! DK DK to Daniel, Jul 2013
+!! DK DK to Daniel, Jul 2013: BEWARE, declared real(kind=CUSTOM_REAL) in trunk and
+!! DK DK to Daniel, Jul 2013: double precision in branch, let us check which one is right
+!! DK DK to Daniel, Jul 2013
+!! DK DK to Daniel, Jul 2013
+!! DK DK to Daniel, Jul 2013
+!! DK DK to Daniel, Jul 2013
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,vnspec) :: one_minus_sum_beta, scale_factor
real(kind=CUSTOM_REAL), dimension(N_SLS,NGLLX,NGLLY,NGLLZ,vnspec) :: factor_common
- double precision, dimension(N_SLS) :: tau_s
- integer i,j,k,ispec
+ double precision, dimension(N_SLS) :: tau_s
+ character(len=150) :: prname
+
+ ! local parameters
+ integer :: i,j,k,ispec
double precision, dimension(N_SLS) :: tau_e, fc
- double precision omsb, Q_mu, sf, T_c_source, scale_t
+ double precision :: omsb, Q_mu, sf, T_c_source, scale_t
! All of the following reads use the output parameters as their temporary arrays
! use the filename to determine the actual contents of the read
@@ -93,10 +107,11 @@
include 'constants.h'
- double precision, dimension(N_SLS) :: tau_s, tau_e, beta, factor_common
- double precision one_minus_sum_beta
+ double precision, dimension(N_SLS),intent(in) :: tau_s, tau_e
+ double precision, dimension(N_SLS),intent(out) :: factor_common
+ double precision,intent(out) :: one_minus_sum_beta
- double precision, dimension(N_SLS) :: tauinv
+ double precision, dimension(N_SLS) :: tauinv,beta
integer i
tauinv(:) = -1.0d0 / tau_s(:)
@@ -184,31 +199,29 @@
end subroutine get_attenuation_scale_factor
-
!
!-------------------------------------------------------------------------------------------------
!
+ subroutine get_attenuation_memory_values(tau_s,deltat, alphaval,betaval,gammaval)
- subroutine get_attenuation_memory_values(tau_s, deltat, alphaval,betaval,gammaval)
-
implicit none
include 'constants.h'
- double precision, dimension(N_SLS) :: tau_s, alphaval, betaval,gammaval
- real(kind=CUSTOM_REAL) deltat
+ double precision, dimension(N_SLS), intent(in) :: tau_s
+ double precision, dimension(N_SLS), intent(out) :: alphaval, betaval,gammaval
+ real(kind=CUSTOM_REAL), intent(in) :: deltat
double precision, dimension(N_SLS) :: tauinv
- tauinv(:) = - 1.0 / tau_s(:)
+ tauinv(:) = - 1.d0 / tau_s(:)
- alphaval(:) = 1 + deltat*tauinv(:) + deltat**2*tauinv(:)**2 / 2. + &
- deltat**3*tauinv(:)**3 / 6. + deltat**4*tauinv(:)**4 / 24.
- betaval(:) = deltat / 2. + deltat**2*tauinv(:) / 3. &
- + deltat**3*tauinv(:)**2 / 8. + deltat**4*tauinv(:)**3 / 24.
- gammaval(:) = deltat / 2. + deltat**2*tauinv(:) / 6. &
- + deltat**3*tauinv(:)**2 / 24.0
+ alphaval(:) = 1.d0 + deltat*tauinv(:) + deltat**2*tauinv(:)**2 / 2.d0 + &
+ deltat**3*tauinv(:)**3 / 6.d0 + deltat**4*tauinv(:)**4 / 24.d0
+ betaval(:) = deltat / 2.d0 + deltat**2*tauinv(:) / 3.d0 &
+ + deltat**3*tauinv(:)**2 / 8.d0 + deltat**4*tauinv(:)**3 / 24.d0
+ gammaval(:) = deltat / 2.d0 + deltat**2*tauinv(:) / 6.d0 &
+ + deltat**3*tauinv(:)**2 / 24.d0
end subroutine get_attenuation_memory_values
-
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/get_event_info.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/get_event_info.f90 2013-07-02 01:23:11 UTC (rev 22484)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/get_event_info.f90 2013-07-02 02:26:51 UTC (rev 22485)
@@ -48,27 +48,19 @@
!--- input or output arguments of the subroutine below
integer, intent(in) :: myrank
+ integer, intent(in) :: NSOURCES ! must be given
integer, intent(out) :: yr,jda,ho,mi
+ double precision, intent(out) :: sec
real, intent(out) :: mb
- double precision, intent(out) :: tshift_cmt,elat,elon,depth,cmt_lat,cmt_lon,cmt_depth,cmt_hdur,sec
-
- !character(len=12), intent(out) :: ename
-
- integer, intent(in) :: NSOURCES ! must be given
+ double precision, intent(out) :: tshift_cmt,elat,elon,depth,cmt_lat,cmt_lon,cmt_depth,cmt_hdur
double precision, intent(out) :: t_shift
character(len=20), intent(out) :: event_name
+ ! local parameters
+ integer :: ier
-
-!--- local variables below
-
- integer ier
-
- !integer, parameter :: LENGTH_REGION_NAME = 150
- !character(len=LENGTH_REGION_NAME) region
-
! get event information for SAC header on the master
if(myrank == 0) then
@@ -161,9 +153,7 @@
character(len=5) datasource
character(len=150) string,CMTSOLUTION
- !character(len=150) string,dummystring,CMTSOLUTION
-
!
!---- read hypocenter info
!
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_receivers.f90 2013-07-02 01:23:11 UTC (rev 22484)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_receivers.f90 2013-07-02 02:26:51 UTC (rev 22485)
@@ -334,10 +334,10 @@
r0=r0*(1.0d0-(2.0d0/3.0d0)*ell*p20)
endif
-! subtract station burial depth (in meters)
+ ! subtract station burial depth (in meters)
r0 = r0 - stbur(irec)/R_EARTH
-! compute the Cartesian position of the receiver
+ ! compute the Cartesian position of the receiver
x_target(irec) = r0*sin(theta)*cos(phi)
y_target(irec) = r0*sin(theta)*sin(phi)
z_target(irec) = r0*cos(theta)
@@ -350,8 +350,8 @@
do ispec=1,nspec
-! loop only on points inside the element
-! exclude edges to ensure this point is not shared with other elements
+ ! loop only on points inside the element
+ ! exclude edges to ensure this point is not shared with other elements
do k=2,NGLLZ-1
do j=2,NGLLY-1
do i=2,NGLLX-1
@@ -377,7 +377,7 @@
! end of loop on all the spectral elements in current slice
enddo
-! end of loop on all the stations
+ ! end of loop on all the stations
enddo
! create RECORDHEADER file with usual format for normal-mode codes
@@ -393,62 +393,47 @@
! compute total number of samples for normal modes with 1 sample per second
open(unit=1,file=trim(OUTPUT_FILES)//'/RECORDHEADERS',status='unknown')
nsamp = nint(dble(NSTEP-1)*DT)
+
do irec = 1,nrec
if(stele(irec) >= -999.9999) then
-! write(1,"(a8,1x,a3,6x,f8.4,1x,f9.4,1x,f6.1,1x,f6.1,f6.1,1x,f6.1,1x,f12.4,1x,i7,1x,i4,1x,i3,1x,i2,1x,i2,1x,f6.3)") &
-! station_name(irec),'LHN',stlat(irec),stlon(irec),stele(irec),stbur(irec), &
-! 0.,0.,1.,nsamp,yr,jda,ho,mi,sec
-! write(1,"(a8,1x,a3,6x,f8.4,1x,f9.4,1x,f6.1,1x,f6.1,f6.1,1x,f6.1,1x,f12.4,1x,i7,1x,i4,1x,i3,1x,i2,1x,i2,1x,f6.3)") &
-! station_name(irec),'LHE',stlat(irec),stlon(irec),stele(irec),stbur(irec), &
-! 90.,0.,1.,nsamp,yr,jda,ho,mi,sec
-! write(1,"(a8,1x,a3,6x,f8.4,1x,f9.4,1x,f6.1,1x,f6.1,f6.1,1x,f6.1,1x,f12.4,1x,i7,1x,i4,1x,i3,1x,i2,1x,i2,1x,f6.3)") &
-! station_name(irec),'LHZ',stlat(irec),stlon(irec),stele(irec),stbur(irec), &
-! 0.,-90.,1.,nsamp,yr,jda,ho,mi,sec
- write(1,"(a8,1x,a3,6x,f8.4,1x,f9.4,1x,f6.1,1x,f6.1,f6.1,1x,f6.1,1x,f12.4,1x,i7,1x,i4,1x,i3,1x,i2,1x,i2,1x,f6.3)") &
- station_name(irec),bic(1:2)//'N',stlat(irec),stlon(irec),stele(irec),stbur(irec), &
- 0.,0.,1.,nsamp,yr,jda,ho,mi,sec
- write(1,"(a8,1x,a3,6x,f8.4,1x,f9.4,1x,f6.1,1x,f6.1,f6.1,1x,f6.1,1x,f12.4,1x,i7,1x,i4,1x,i3,1x,i2,1x,i2,1x,f6.3)") &
- station_name(irec),bic(1:2)//'E',stlat(irec),stlon(irec),stele(irec),stbur(irec), &
- 90.,0.,1.,nsamp,yr,jda,ho,mi,sec
- write(1,"(a8,1x,a3,6x,f8.4,1x,f9.4,1x,f6.1,1x,f6.1,f6.1,1x,f6.1,1x,f12.4,1x,i7,1x,i4,1x,i3,1x,i2,1x,i2,1x,f6.3)") &
- station_name(irec),bic(1:2)//'Z',stlat(irec),stlon(irec),stele(irec),stbur(irec), &
- 0.,-90.,1.,nsamp,yr,jda,ho,mi,sec
-
+ write(1,500) station_name(irec),bic(1:2)//'N', &
+ stlat(irec),stlon(irec),stele(irec),stbur(irec), &
+ 0.,0.,1.,nsamp,yr,jda,ho,mi,sec
+ write(1,500) station_name(irec),bic(1:2)//'E', &
+ stlat(irec),stlon(irec),stele(irec),stbur(irec), &
+ 90.,0.,1.,nsamp,yr,jda,ho,mi,sec
+ write(1,500) station_name(irec),bic(1:2)//'Z', &
+ stlat(irec),stlon(irec),stele(irec),stbur(irec), &
+ 0.,-90.,1.,nsamp,yr,jda,ho,mi,sec
else
! very deep ocean-bottom stations such as H2O are not compatible
! with the standard RECORDHEADERS format because of the f6.1 format
! therefore suppress decimals for depth in that case
-! write(1,"(a8,1x,a3,6x,f8.4,1x,f9.4,1x,i6,1x,f6.1,f6.1,1x,f6.1,1x,f12.4,1x,i7,1x,i4,1x,i3,1x,i2,1x,i2,1x,f6.3)") &
-! station_name(irec),'LHN',stlat(irec),stlon(irec),nint(stele(irec)),stbur(irec), &
-! 0.,0.,1.,nsamp,yr,jda,ho,mi,sec
-! write(1,"(a8,1x,a3,6x,f8.4,1x,f9.4,1x,i6,1x,f6.1,f6.1,1x,f6.1,1x,f12.4,1x,i7,1x,i4,1x,i3,1x,i2,1x,i2,1x,f6.3)") &
-! station_name(irec),'LHE',stlat(irec),stlon(irec),nint(stele(irec)),stbur(irec), &
-! 90.,0.,1.,nsamp,yr,jda,ho,mi,sec
-! write(1,"(a8,1x,a3,6x,f8.4,1x,f9.4,1x,i6,1x,f6.1,f6.1,1x,f6.1,1x,f12.4,1x,i7,1x,i4,1x,i3,1x,i2,1x,i2,1x,f6.3)") &
-! station_name(irec),'LHZ',stlat(irec),stlon(irec),nint(stele(irec)),stbur(irec), &
-! 0.,-90.,1.,nsamp,yr,jda,ho,mi,sec
- write(1,"(a8,1x,a3,6x,f8.4,1x,f9.4,1x,i6,1x,f6.1,f6.1,1x,f6.1,1x,f12.4,1x,i7,1x,i4,1x,i3,1x,i2,1x,i2,1x,f6.3)") &
- station_name(irec),bic(1:2)//'N',stlat(irec),stlon(irec),nint(stele(irec)),stbur(irec), &
- 0.,0.,1.,nsamp,yr,jda,ho,mi,sec
- write(1,"(a8,1x,a3,6x,f8.4,1x,f9.4,1x,i6,1x,f6.1,f6.1,1x,f6.1,1x,f12.4,1x,i7,1x,i4,1x,i3,1x,i2,1x,i2,1x,f6.3)") &
- station_name(irec),bic(1:2)//'E',stlat(irec),stlon(irec),nint(stele(irec)),stbur(irec), &
- 90.,0.,1.,nsamp,yr,jda,ho,mi,sec
- write(1,"(a8,1x,a3,6x,f8.4,1x,f9.4,1x,i6,1x,f6.1,f6.1,1x,f6.1,1x,f12.4,1x,i7,1x,i4,1x,i3,1x,i2,1x,i2,1x,f6.3)") &
- station_name(irec),bic(1:2)//'Z',stlat(irec),stlon(irec),nint(stele(irec)),stbur(irec), &
- 0.,-90.,1.,nsamp,yr,jda,ho,mi,sec
-
+ write(1,600) station_name(irec),bic(1:2)//'N', &
+ stlat(irec),stlon(irec),nint(stele(irec)),stbur(irec), &
+ 0.,0.,1.,nsamp,yr,jda,ho,mi,sec
+ write(1,600) station_name(irec),bic(1:2)//'E', &
+ stlat(irec),stlon(irec),nint(stele(irec)),stbur(irec), &
+ 90.,0.,1.,nsamp,yr,jda,ho,mi,sec
+ write(1,600) station_name(irec),bic(1:2)//'Z', &
+ stlat(irec),stlon(irec),nint(stele(irec)),stbur(irec), &
+ 0.,-90.,1.,nsamp,yr,jda,ho,mi,sec
endif
enddo
close(1)
endif
+500 format(a8,1x,a3,6x,f9.4,1x,f9.4,1x,f6.1,1x,f6.1,1x,f6.1,1x,f6.1,1x,f12.4,1x,i7,1x,i4.4,1x,i3.3,1x,i2.2,1x,i2.2,1x,f6.3)
+600 format(a8,1x,a3,6x,f9.4,1x,f9.4,1x,i6,1x,f6.1,f6.1,1x,f6.1,1x,f12.4,1x,i7,1x,i4.4,1x,i3.3,1x,i2.2,1x,i2.2,1x,f6.3)
+
+
! ****************************************
! find the best (xi,eta) for each receiver
! ****************************************
-!! loop on all the receivers to iterate in that slice
+! ! loop on all the receivers to iterate in that slice
! do irec = 1,nrec
! loop on all the receivers
@@ -495,8 +480,7 @@
eta = yigll(iy_initial_guess(irec))
gamma = zigll(iz_initial_guess(irec))
-! define coordinates of the control points of the element
-
+ ! define coordinates of the control points of the element
do ia=1,NGNOD
if(iaddx(ia) == 0) then
@@ -536,17 +520,17 @@
enddo
-! iterate to solve the non linear system
+ ! iterate to solve the non linear system
do iter_loop = 1,NUM_ITER
-! impose receiver exactly at the surface
+ ! impose receiver exactly at the surface
if(.not. RECEIVERS_CAN_BE_BURIED) gamma = 1.d0
-! recompute jacobian for the new point
+ ! recompute jacobian for the new point
call recompute_jacobian(xelm,yelm,zelm,xi,eta,gamma,x,y,z, &
- xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz)
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz)
-! compute distance to target location
+ ! compute distance to target location
dx = - (x - x_target(irec))
dy = - (y - y_target(irec))
dz = - (z - z_target(irec))
@@ -562,11 +546,11 @@
eta = eta + deta
if(RECEIVERS_CAN_BE_BURIED) gamma = gamma + dgamma
-! impose that we stay in that element
-! (useful if user gives a receiver outside the mesh for instance)
-! we can go slightly outside the [1,1] segment since with finite elements
-! the polynomial solution is defined everywhere
-! can be useful for convergence of iterative scheme with distorted elements
+ ! impose that we stay in that element
+ ! (useful if user gives a receiver outside the mesh for instance)
+ ! we can go slightly outside the [1,1] segment since with finite elements
+ ! the polynomial solution is defined everywhere
+ ! can be useful for convergence of iterative scheme with distorted elements
if (xi > 1.10d0) xi = 1.10d0
if (xi < -1.10d0) xi = -1.10d0
if (eta > 1.10d0) eta = 1.10d0
@@ -584,7 +568,7 @@
call recompute_jacobian(xelm,yelm,zelm,xi,eta,gamma,x,y,z, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz)
-! store xi,eta and x,y,z of point found
+ ! store xi,eta and x,y,z of point found
xi_receiver_subset(irec_in_this_subset) = xi
eta_receiver_subset(irec_in_this_subset) = eta
gamma_receiver_subset(irec_in_this_subset) = gamma
@@ -595,7 +579,7 @@
y_found_subset(irec_in_this_subset) = y_found(irec)
z_found_subset(irec_in_this_subset) = z_found(irec)
-! compute final distance between asked and found (converted to km)
+ ! compute final distance between asked and found (converted to km)
final_distance(irec) = dsqrt((x_target(irec)-x_found(irec))**2 + &
(y_target(irec)-y_found(irec))**2 + (z_target(irec)-z_found(irec))**2)*R_EARTH/1000.d0
@@ -603,7 +587,7 @@
enddo ! end of loop on all stations within current subset
-! for MPI version, gather information from all the nodes
+ ! for MPI version, gather information from all the nodes
ispec_selected_rec_all(:,:) = -1
call MPI_GATHER(ispec_selected_rec_subset,nrec_SUBSET_current_size,MPI_INTEGER,ispec_selected_rec_all,nrec_SUBSET_current_size, &
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/noise_tomography.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/noise_tomography.f90 2013-07-02 01:23:11 UTC (rev 22484)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/noise_tomography.f90 2013-07-02 02:26:51 UTC (rev 22485)
@@ -25,10 +25,6 @@
!
!=====================================================================
-! =============================================================================================================
-! =============================================================================================================
-! =============================================================================================================
-
! subroutine for NOISE TOMOGRAPHY
! chracterize noise statistics
! for a given point (xcoord,ycoord,zcoord), specify the noise direction "normal_x/y/z_noise"
@@ -76,6 +72,7 @@
! subroutine for NOISE TOMOGRAPHY
! read parameters
+
subroutine read_parameters_noise(myrank,nrec,NSTEP,nmovie_points, &
islice_selected_rec,xi_receiver,eta_receiver,gamma_receiver,nu, &
noise_sourcearray,xigll,yigll,zigll,nspec_top, &
@@ -83,6 +80,7 @@
xstore_crust_mantle,ystore_crust_mantle,zstore_crust_mantle, &
irec_master_noise,normal_x_noise,normal_y_noise,normal_z_noise,mask_noise)
implicit none
+
include 'mpif.h'
include "precision.h"
include "constants.h"
@@ -418,6 +416,7 @@
ibool_crust_mantle,islice_selected_rec,ispec_selected_rec, &
it,irec_master_noise)
implicit none
+
include "constants.h"
include "OUTPUT_FILES/values_from_mesher.h"
! input parameters
@@ -427,13 +426,14 @@
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NGLLZ,NSTEP) :: noise_sourcearray
real(kind=CUSTOM_REAL),dimension(NDIM,NGLOB_CRUST_MANTLE) :: accel_crust_mantle ! both input and output
! output parameters
+
! local parameters
integer :: i,j,k,iglob,it
! adds noise source (only if this proc carries the noise)
if(myrank == islice_selected_rec(irec_master_noise)) then
- ! adds nosie source contributions
+ ! adds noise source contributions
do k=1,NGLLZ
do j=1,NGLLY
do i=1,NGLLX
@@ -468,6 +468,7 @@
ibelm_top_crust_mantle,ibool_crust_mantle, &
nspec_top,noise_surface_movie,it)
implicit none
+
include "constants.h"
include "OUTPUT_FILES/values_from_mesher.h"
! input parameters
@@ -476,6 +477,7 @@
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: ibool_crust_mantle
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_CRUST_MANTLE) :: displ_crust_mantle
! output parameters
+
! local parameters
integer :: ispec2D,ispec,i,j,k,iglob
real(kind=CUSTOM_REAL),dimension(NDIM,NGLLX,NGLLY,nspec_top) :: noise_surface_movie
@@ -610,17 +612,17 @@
do ispec2D = 1, nspec_top ! NSPEC2D_TOP(IREGION_CRUST_MANTLE)
ispec = ibelm_top_crust_mantle(ispec2D)
- k = NGLLZ
+ k = NGLLZ
- ! loop on all the points inside the element
- do j = 1,NGLLY
- do i = 1,NGLLX
- ipoin = ipoin + 1
- iglob = ibool_crust_mantle(i,j,k,ispec)
+ ! loop on all the points inside the element
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+ ipoin = ipoin + 1
+ iglob = ibool_crust_mantle(i,j,k,ispec)
- eta = noise_surface_movie(1,i,j,ispec2D) * normal_x_noise(ipoin) + &
- noise_surface_movie(2,i,j,ispec2D) * normal_y_noise(ipoin) + &
- noise_surface_movie(3,i,j,ispec2D) * normal_z_noise(ipoin)
+ eta = noise_surface_movie(1,i,j,ispec2D) * normal_x_noise(ipoin) + &
+ noise_surface_movie(2,i,j,ispec2D) * normal_y_noise(ipoin) + &
+ noise_surface_movie(3,i,j,ispec2D) * normal_z_noise(ipoin)
accel_crust_mantle(1,iglob) = accel_crust_mantle(1,iglob) + eta * mask_noise(ipoin) * normal_x_noise(ipoin) &
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/read_forward_arrays.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/read_forward_arrays.f90 2013-07-02 01:23:11 UTC (rev 22484)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/read_forward_arrays.f90 2013-07-02 02:26:51 UTC (rev 22485)
@@ -122,21 +122,26 @@
close(55)
endif
+ ! initializes backward/reconstructed arrays
if (SIMULATION_TYPE == 3) then
- ! initializes
+ ! initializes wavefields
b_displ_crust_mantle = 0._CUSTOM_REAL
b_veloc_crust_mantle = 0._CUSTOM_REAL
b_accel_crust_mantle = 0._CUSTOM_REAL
+
b_displ_inner_core = 0._CUSTOM_REAL
b_veloc_inner_core = 0._CUSTOM_REAL
b_accel_inner_core = 0._CUSTOM_REAL
+
b_displ_outer_core = 0._CUSTOM_REAL
b_veloc_outer_core = 0._CUSTOM_REAL
b_accel_outer_core = 0._CUSTOM_REAL
+
if (ROTATION_VAL) then
b_A_array_rotation = 0._CUSTOM_REAL
b_B_array_rotation = 0._CUSTOM_REAL
endif
+
if (ATTENUATION_VAL) then
b_R_memory_crust_mantle = 0._CUSTOM_REAL
b_R_memory_inner_core = 0._CUSTOM_REAL
@@ -156,7 +161,7 @@
b_R_memory_crust_mantle,b_R_memory_inner_core, &
b_A_array_rotation,b_B_array_rotation,LOCAL_PATH)
-! reads in saved wavefields
+! reads in saved wavefields to reconstruct/backward wavefield
implicit none
@@ -184,7 +189,7 @@
character(len=150) LOCAL_PATH
- !local parameters
+ ! local parameters
character(len=150) outputname
write(outputname,'(a,i6.6,a)') 'proc',myrank,'_save_forward_arrays.bin'
More information about the CIG-COMMITS
mailing list