[cig-commits] commit: improve html documentation with doxygen
Mercurial
hg at geodynamics.org
Tue Sep 20 12:13:14 PDT 2011
changeset: 16:0a07ece6fb5a
user: Sylvain Barbot <sylbar.vainbot at gmail.com>
date: Thu May 12 19:01:34 2011 -0700
files: .bzrignore elastic3d.f90 export.f90 fourier.f90 friction3d.f90 green.f90 relax.f90 viscoelastic3d.f90
description:
improve html documentation with doxygen
diff -r 777cc20187d0 -r 0a07ece6fb5a .bzrignore
--- a/.bzrignore Thu May 12 17:04:47 2011 -0700
+++ b/.bzrignore Thu May 12 19:01:34 2011 -0700
@@ -3,3 +3,4 @@
*~
doc
./relax
+examples/run1_proj.sh
diff -r 777cc20187d0 -r 0a07ece6fb5a elastic3d.f90
--- a/elastic3d.f90 Thu May 12 17:04:47 2011 -0700
+++ b/elastic3d.f90 Thu May 12 19:01:34 2011 -0700
@@ -54,9 +54,9 @@ CONTAINS
CONTAINS
!------------------------------------------------------------
- ! function SIGN
- ! returns the sign of the input -1 for negtive, 0 for zero
- ! and +1 for positive arguments.
+ !> function SIGN
+ !! returns the sign of the input -1 for negtive, 0 for zero
+ !! and +1 for positive arguments.
!------------------------------------------------------------
REAL*8 FUNCTION sign(x)
REAL*8, INTENT(IN) :: x
@@ -73,8 +73,8 @@ CONTAINS
END FUNCTION sign
!------------------------------------------------------------
- ! function fix
- ! returns the closest integer scalar
+ !> function fix
+ !! returns the closest integer scalar
!
! sylvain barbot (08/25/07) - original form
!------------------------------------------------------------
@@ -94,8 +94,8 @@ CONTAINS
END FUNCTION fix
!------------------------------------------------------------
- ! function SINH
- ! computes the hyperbolic sine
+ !> function SINH
+ !! computes the hyperbolic sine
!------------------------------------------------------------
REAL*8 FUNCTION sinh(x)
REAL*8, INTENT(IN) :: x
@@ -108,8 +108,8 @@ CONTAINS
END FUNCTION sinh
!------------------------------------------------------------
- ! function ASINH
- ! computes the inverse hyperbolic sine
+ !> function ASINH
+ !! computes the inverse hyperbolic sine
!------------------------------------------------------------
REAL*8 FUNCTION asinh(x)
REAL*8, INTENT(IN) :: x
@@ -117,14 +117,14 @@ CONTAINS
END FUNCTION asinh
!-----------------------------------------------------------------
- ! subroutine Neighbor
- ! computes the indices of neighbor samples (l points away)
- ! bracketing the current samples location i1,i2,i3 and
- ! assuming periodic boundary condition.
- !
- ! i1m < i1 < i1p
- ! i2m < i2 < i2p
- ! i3m < i3 < i3p
+ !> subroutine Neighbor
+ !! computes the indices of neighbor samples (l points away)
+ !! bracketing the current samples location i1,i2,i3 and
+ !! assuming periodic boundary condition.
+ !!
+ !! i1m < i1 < i1p
+ !! i2m < i2 < i2p
+ !! i3m < i3 < i3p
!-----------------------------------------------------------------
SUBROUTINE neighbor(i1,i2,i3,sx1,sx2,sx3,l,i1m,i1p,i2m,i2p,i3m,i3p)
INTEGER, INTENT(IN) :: i1,i2,i3,sx1,sx2,sx3,l
@@ -140,9 +140,9 @@ CONTAINS
END SUBROUTINE neighbor
!---------------------------------------------------------------
- ! subroutine IsotropicStressStrain
- ! computes in place the isotropic stress tensor from a given
- ! strain tensor using Hooke's law stress-strain relationship.
+ !> subroutine IsotropicStressStrain
+ !! computes in place the isotropic stress tensor from a given
+ !! strain tensor using Hooke's law stress-strain relationship.
!
! sylvain barbot (10/14/07) - original form
!---------------------------------------------------------------
@@ -162,8 +162,8 @@ CONTAINS
END SUBROUTINE isotropicstressstrain
!------------------------------------------------------------
- ! function TensorDiff
- ! computes the difference between two tensors: t=t1-t2
+ !> function TensorDiff
+ !! computes the difference between two tensors: t=t1-t2
!
! sylvain barbot (07/09/07) - original form
!------------------------------------------------------------
@@ -180,8 +180,8 @@ CONTAINS
END FUNCTION tensordiff
!------------------------------------------------------------
- ! function TensorPlus
- ! computes the sum of two tensors: t=t1-t2
+ !> function TensorPlus
+ !! computes the sum of two tensors: t=t1-t2
!
! sylvain barbot (07/09/07) - original form
!------------------------------------------------------------
@@ -198,8 +198,8 @@ CONTAINS
END FUNCTION tensorplus
!------------------------------------------------------------
- ! function TensorScalarProd
- ! multiplies a tensor with a scalar
+ !> function TensorScalarProd
+ !! multiplies a tensor with a scalar
!
! sylvain barbot (07/09/07) - original form
!------------------------------------------------------------
@@ -217,9 +217,9 @@ CONTAINS
END FUNCTION tensorscalarprod
!------------------------------------------------------------
- ! function TensorSymmetricDyadProd
- ! computes the dyadic product of two vectors to obtain a
- ! symmetric second order tensor
+ !> function TensorSymmetricDyadProd
+ !! computes the dyadic product of two vectors to obtain a
+ !! symmetric second order tensor
!
! sylvain barbot (07/09/07) - original form
!------------------------------------------------------------
@@ -238,9 +238,9 @@ CONTAINS
END FUNCTION tensorsymmetricdyadprod
!------------------------------------------------------------
- ! function TensorVectorDotProd
- ! compute the dot product T.v where T is a second-order
- ! tensor and v is a vector.
+ !> function TensorVectorDotProd
+ !! compute the dot product T.v where T is a second-order
+ !! tensor and v is a vector.
!
! sylvain barbot (07/09/07) - original form
!------------------------------------------------------------
@@ -257,9 +257,9 @@ CONTAINS
END FUNCTION tensorvectordotprod
!------------------------------------------------------------
- ! function TensorVectorDotProd
- ! compute the dot product T.v where T is a second-order
- ! tensor and v is a vector.
+ !> function TensorVectorDotProd
+ !! compute the dot product T.v where T is a second-order
+ !! tensor and v is a vector.
!
! sylvain barbot (07/09/07) - original form
!------------------------------------------------------------
@@ -281,8 +281,8 @@ CONTAINS
END FUNCTION tensordeviatoric
!------------------------------------------------------------
- ! function TensorTrace
- ! computes the trace of a second order tensor
+ !> function TensorTrace
+ !! computes the trace of a second order tensor
!
! sylvain barbot (07/09/07) - original form
!------------------------------------------------------------
@@ -294,8 +294,8 @@ CONTAINS
END FUNCTION tensortrace
!------------------------------------------------------------
- ! function TensorNorm
- ! computes the Frobenius norm of a second order tensor
+ !> function TensorNorm
+ !! computes the Frobenius norm of a second order tensor
!
! sylvain barbot (07/09/07) - original form
!------------------------------------------------------------
@@ -310,13 +310,13 @@ CONTAINS
END FUNCTION tensornorm
!------------------------------------------------------------
- ! function TensorDecomposition
- ! writes a tensor t as the product of a norm and a direction
- !
- ! t = gamma * R
- !
- ! where gamma is a scalar, the norm of t, and R is a unitary
- ! tensor. t is assumed to be a deviatoric tensor.
+ !> function TensorDecomposition
+ !! writes a tensor t as the product of a norm and a direction
+ !!
+ !! t = gamma * R
+ !!
+ !! where gamma is a scalar, the norm of t, and R is a unitary
+ !! tensor. t is assumed to be a deviatoric tensor.
!
! sylvain barbot (07/09/07) - original form
!------------------------------------------------------------
@@ -338,8 +338,8 @@ CONTAINS
!------------------------------------------------------------
- ! function TensorForbeniusNorm
- ! computes the Frobenius norm of a second order tensor
+ !> function TensorForbeniusNorm
+ !! computes the Frobenius norm of a second order tensor
!
! sylvain barbot (07/09/07) - original form
!------------------------------------------------------------
@@ -354,8 +354,8 @@ CONTAINS
END FUNCTION tensorfrobeniusnorm
!------------------------------------------------------------
- ! function VectorFieldNormMax
- ! computes the maximum value of the norm of a vector field
+ !> function VectorFieldNormMax
+ !! computes the maximum value of the norm of a vector field
!------------------------------------------------------------
SUBROUTINE vectorfieldnormmax(v1,v2,v3,sx1,sx2,sx3,maximum,location)
INTEGER, INTENT(IN) :: sx1,sx2,sx3
@@ -386,8 +386,8 @@ CONTAINS
END SUBROUTINE vectorfieldnormmax
!------------------------------------------------------------
- ! function TensorMean
- ! computesthe mean of the norm of a tensor field
+ !> function TensorMean
+ !! computesthe mean of the norm of a tensor field
!------------------------------------------------------------
REAL*8 FUNCTION tensormean(t)
TYPE(TENSOR), INTENT(IN), DIMENSION(:,:,:) :: t
@@ -409,8 +409,8 @@ CONTAINS
END FUNCTION tensormean
!------------------------------------------------------------
- ! function TensorAmplitude
- ! computes the integral of the norm of a tensor field
+ !> function TensorAmplitude
+ !! computes the integral of the norm of a tensor field
!------------------------------------------------------------
REAL*8 FUNCTION tensoramplitude(t,dx1,dx2,dx3)
TYPE(TENSOR), INTENT(IN), DIMENSION(:,:,:) :: t
@@ -435,8 +435,8 @@ CONTAINS
END FUNCTION tensoramplitude
!------------------------------------------------------------
- ! function TensorMeanTrace
- ! computesthe mean of the norm of a tensor field
+ !> function TensorMeanTrace
+ !! computesthe mean of the norm of a tensor field
!------------------------------------------------------------
REAL*8 FUNCTION tensormeantrace(t)
TYPE(TENSOR), INTENT(IN), DIMENSION(:,:,:) :: t
@@ -459,8 +459,8 @@ CONTAINS
END FUNCTION tensormeantrace
!------------------------------------------------------------
- ! sinc function
- ! computes sin(pi*x)/(pi*x)
+ !> sinc function
+ !! computes sin(pi*x)/(pi*x)
!
! sylvain barbot (04-14-07) - original form
!------------------------------------------------------------
@@ -475,7 +475,7 @@ CONTAINS
END FUNCTION sinc
!-------------------------------------------------------------------------
- ! function gauss computes the normalized gaussian function
+ !> function gauss computes the normalized gaussian function
!
! Sylvain Barbot (06-29-07)
!-------------------------------------------------------------------------
@@ -487,7 +487,7 @@ CONTAINS
END FUNCTION gauss
!-------------------------------------------------------------------------
- ! function gaussp computes the normalized gaussian derivative
+ !> function gaussp computes the normalized gaussian derivative
!
! Sylvain Barbot (06-29-07)
!-------------------------------------------------------------------------
@@ -499,7 +499,7 @@ CONTAINS
END FUNCTION gaussp
!-------------------------------------------------------------------------
- ! function omega computes raised-cosine taper in the space domain
+ !> function omega computes raised-cosine taper in the space domain
!
! Sylvain Barbot (06-29-07)
!-------------------------------------------------------------------------
@@ -519,8 +519,8 @@ CONTAINS
END FUNCTION omega
!-------------------------------------------------------------------------
- ! function omegap computes raised-cosine taper derivative
- ! in the space domain
+ !> function omegap computes raised-cosine taper derivative
+ !! in the space domain
!
! Sylvain Barbot (06-29-07)
!-------------------------------------------------------------------------
@@ -538,13 +538,13 @@ CONTAINS
END FUNCTION omegap
!-------------------------------------------------------------------------
- ! tapered step function (raised-cosine) of unit area in the Fourier domain
- !
- ! INPUT
- ! k wavenumber
- ! beta roll-off parameter 0<beta<0.5
- ! no smoothing for beta close to 0
- ! string smoothing for beta close to 0.5
+ !> tapered step function (raised-cosine) of unit area in the Fourier domain
+ !!
+ !! INPUT
+ !! @param k wavenumber
+ !! @param beta roll-off parameter 0<beta<0.5
+ !! no smoothing for beta close to 0
+ !! string smoothing for beta close to 0.5
!
! sylvain barbot (04-14-07) - original form
!-------------------------------------------------------------------------
@@ -562,10 +562,10 @@ CONTAINS
END FUNCTION omegak
!----------------------------------------------------------------
- ! subroutine TensorStructure
- ! constructs a vertically-stratified tensor field.
- ! The structure is defined by its interfaces: changes can be
- ! gradual or discontinuous.
+ !> subroutine TensorStructure
+ !! constructs a vertically-stratified tensor field.
+ !! The structure is defined by its interfaces: changes can be
+ !! gradual or discontinuous.
!
! sylvain barbot (10/25/08) - original form
!----------------------------------------------------------------
@@ -643,31 +643,31 @@ CONTAINS
!----------------------------------------------------------------
- ! subroutine ViscoElasticStructure
- ! constructs a vertically-stratified viscoelastic structure.
- ! The structure is defined by its interfaces: changes can be
- ! gradual or discontinuous.
- !
- ! EXAMPLE INPUTS:
- !
- ! 1- elastic plate over linear viscous half-space
- ! 1
- ! 1 1.0 1.0 1.0
- !
- ! 2- elastic plate over powerlaw viscous half-space (n=3)
- ! 1
- ! 1 1.0 1.0 3.0
- !
- ! 3- elastic plate over viscous half-space with depth-dependent
- ! viscosity
- ! 2
- ! 1 01.0 1.0 1.0
- ! 2 10.0 6.0 1.0
- !
- ! in this last example, the grid does not have to reach down
- ! to x3=10.
- !
- ! sylvain barbot (08/07/07) - original form
+ !> subroutine ViscoElasticStructure
+ !! constructs a vertically-stratified viscoelastic structure.
+ !! The structure is defined by its interfaces: changes can be
+ !! gradual or discontinuous.
+ !!
+ !! EXAMPLE INPUTS:
+ !!
+ !! 1- elastic plate over linear viscous half-space
+ !! 1
+ !! 1 1.0 1.0 1.0
+ !!
+ !! 2- elastic plate over powerlaw viscous half-space (n=3)
+ !! 1
+ !! 1 1.0 1.0 3.0
+ !!
+ !! 3- elastic plate over viscous half-space with depth-dependent
+ !! viscosity
+ !! 2
+ !! 1 01.0 1.0 1.0
+ !! 2 10.0 6.0 1.0
+ !!
+ !! in this last example, the grid does not have to reach down
+ !! to x3=10.
+ !!
+ !! \author sylvain barbot (08/07/07) - original form
!----------------------------------------------------------------
SUBROUTINE viscoelasticstructure(vstruct,layers,dx3)
TYPE(LAYER_STRUCT), INTENT(IN), DIMENSION(:) :: layers
@@ -762,16 +762,16 @@ CONTAINS
!------------------------------------------------------------------
- ! function OptimalFilter
- ! load predefined Finite Impulse Response (FIR) filters of various
- ! lengths and select the most appropriate ones based on the
- ! computational grid size. result is filter kernels always smaller
- ! than available computational length.
- ! this is useful in the special cases of infinite faults where
- ! deformation is essentially two-dimensional, despite the actual
- ! three-dimensional computation. in the direction of symmetry,
- ! no strain occurs and high accuracy derivative estimates are not
- ! needed.
+ !> function OptimalFilter
+ !! load predefined Finite Impulse Response (FIR) filters of various
+ !! lengths and select the most appropriate ones based on the
+ !! computational grid size. result is filter kernels always smaller
+ !! than available computational length.
+ !! this is useful in the special cases of infinite faults where
+ !! deformation is essentially two-dimensional, despite the actual
+ !! three-dimensional computation. in the direction of symmetry,
+ !! no strain occurs and high accuracy derivative estimates are not
+ !! needed.
!
! Sylvain Barbot (03/05/08) - original form
!------------------------------------------------------------------
@@ -839,48 +839,51 @@ CONTAINS
END SUBROUTINE optimalfilter
!-----------------------------------------------------------------
- ! subroutine StressUpdate
- ! computes the 3-d stress tensor sigma_ij' from the current
- ! deformation field. Strain is the second order tensor
- !
- ! epsilon_ij = 1/2 ( u_i,j + u_j,i )
- !
- ! The displacement derivatives are approximated numerically by the
- ! application of a differentiator space-domain finite impulse
- ! response filter. Coefficients of the filter can be obtained with
- ! the MATLAB command line
- !
- ! firpm(14, ...
- ! [0 7.0e-1 8.000000e-1 8.500000e-1 9.000000e-1 1.0e+0],...
- ! [0 7.0e-1 5.459372e-1 3.825260e-1 2.433534e-1 0.0e+0]*pi,...
- ! 'differentiator');
- !
- ! The kernel is odd and antisymmetric and only half the numbers
- ! are stored in this code. Kernels of different sizes are readilly
- ! available in the 'kernelX.inc' files. Stress tensor field is
- ! obtained by application of Hooke's law
- !
- ! sigma' = - C' : E
- !
- ! or in indicial notation
- !
- ! sigma_ij' = -lambda'*delta_ij*epsilon_kk - 2*mu'*epsilon_ij
- !
- ! where C' is the heterogeneous elastic moduli tensor and lambda'
- ! and mu' are the inhomogeneous lame parameters
- !
- ! C' = C(x) - C_0
- !
- ! For isotropic materials
- !
- ! mu'(x) = mu(x) - mu_0
- ! lambda'(x) = lambda(x) - lambda_0
- !
- ! Optionally, the surface traction sigma_i3 can be sampled.
- !
- ! sylvain barbot (10/10/07) - original form
- ! - optional sample of normal stress
- ! (02/12/09) - OpemMP parallel implementation
+ !> subroutine StressUpdate
+ !! computes the 3-d stress tensor sigma_ij' from the current
+ !! deformation field. Strain is the second order tensor
+ !!
+ !! \f[ \epsilon_{ij} = \frac{1}{2} ( u_{i,j} + u_{j,i} ) \f]
+ !!
+ !! The displacement derivatives are approximated numerically by the
+ !! application of a differentiator space-domain finite impulse
+ !! response filter. Coefficients of the filter can be obtained with
+ !! the MATLAB command line
+ !!
+ !!\verbatim
+ !! firpm(14, ...
+ !! [0 7.0e-1 8.000000e-1 8.500000e-1 9.000000e-1 1.0e+0],...
+ !! [0 7.0e-1 5.459372e-1 3.825260e-1 2.433534e-1 0.0e+0]*pi,...
+ !! 'differentiator');
+ !!\endverbatim
+ !!
+ !! The kernel is odd and antisymmetric and only half the numbers
+ !! are stored in this code. Kernels of different sizes are readilly
+ !! available in the 'kernelX.inc' files. Stress tensor field is
+ !! obtained by application of Hooke's law
+ !!
+ !! \f[ \sigma' = - C' : E \f]
+ !!
+ !! or in indicial notation
+ !!
+ !!
+ !! \f[ \sigma_{ij}' = -\lambda'*\delta_{ij}*\epsilon_{kk} - 2*\mu'*\epsilon_{ij}\f]
+ !!
+ !! where C' is the heterogeneous elastic moduli tensor and lambda'
+ !! and mu' are the inhomogeneous lame parameters
+ !!
+ !! \f[ C' = C(x) - C_0 \f]
+ !!
+ !! For isotropic materials
+ !!
+ !! \f[ \mu'(x) = \mu(x) - \mu_0 \f]
+ !! \f[ \lambda'(x) = \lambda(x) - \lambda_0 \f]
+ !!
+ !! Optionally, the surface traction sigma_i3 can be sampled.
+ !!
+ !! \author sylvain barbot (10/10/07) - original form
+ !! - optional sample of normal stress
+ !! (02/12/09) - OpemMP parallel implementation
!-----------------------------------------------------------------
SUBROUTINE stressupdate(v1,v2,v3,lambda,mu,dx1,dx2,dx3,sx1,sx2,sx3,sig)
INTEGER, INTENT(IN) :: sx1,sx2,sx3
@@ -948,13 +951,13 @@ CONTAINS
CONTAINS
!---------------------------------------------------------------
- ! LocalStrain_FIR2
- ! implements a finite impulse response filter (FIR) to estimate
- ! derivatives and strain components. the compatibility with the
- ! OpenMP parallel execution requires that all variable be
- ! tractable from the calling routine.
- !
- ! sylvain barbot (10/10/07) - original form
+ !> LocalStrain_FIR2
+ !! implements a finite impulse response filter (FIR) to estimate
+ !! derivatives and strain components. the compatibility with the
+ !! OpenMP parallel execution requires that all variable be
+ !! tractable from the calling routine.
+ !!
+ !! \author sylvain barbot (10/10/07) - original form
! (03/05/08) - implements 3 filters
! (02/12/09) - compatibility with OpenMP (scope)
!---------------------------------------------------------------
@@ -1007,12 +1010,12 @@ CONTAINS
END SUBROUTINE localstrain_fir2
!---------------------------------------------------------------
- ! LocalStrain_FIR
- ! implements a finite impulse response filter (FIR) to estimate
- ! derivatives and strain components.
- !
- ! sylvain barbot (10/10/07) - original form
- ! (03/05/08) - implements 3 filters
+ !> LocalStrain_FIR
+ !! implements a finite impulse response filter (FIR) to estimate
+ !! derivatives and strain components.
+ !!
+ !! \author sylvain barbot (10/10/07) - original form
+ !! (03/05/08) - implements 3 filters
!---------------------------------------------------------------
SUBROUTINE localstrain_fir(e)
TYPE(TENSOR), INTENT(OUT) :: e
@@ -1058,10 +1061,10 @@ CONTAINS
END SUBROUTINE localstrain_fir
!---------------------------------------------------------------
- ! LocalStrain_ANI
- ! implements a different finite impulse response filter (FIR)
- ! in each direction (ANIsotropy) to estimate derivatives and
- ! strain components.
+ !> LocalStrain_ANI
+ !! implements a different finite impulse response filter (FIR)
+ !! in each direction (ANIsotropy) to estimate derivatives and
+ !! strain components.
!
! sylvain barbot (10/10/07) - original form
! (03/05/09) - implements 3 filters
@@ -1109,36 +1112,36 @@ CONTAINS
END SUBROUTINE stressupdate
!-----------------------------------------------------------------
- ! subroutine EquivalentBodyForce
- ! computes and updates the equivalent body-force
- !
- ! f = - div.( C : E^i )
- !
- ! and the equivalent surface traction
- !
- ! t = n . C : E^i
- !
- ! with n = (0,0,-1). In indicial notations
- !
- ! f_i = - (C_ijkl E^i_kl),j
- !
- ! and
- !
- ! t_1 = n_j C_ijkl E^i_kl
- !
- ! where f is the equivalent body-force, t is the equivalent surface
- ! traction, C is the elastic moduli tensor and E^i is the moment
- ! density tensor tensor.
- !
- ! Divergence is computed with a mixed numerical scheme including
- ! centered finite-difference (in the vertical direction) and
- ! finite impulse response differentiator filter for derivatives
- ! estimates. see function 'stress' for further explanations.
- !
- ! sylvain barbot (07/09/07) - original form
- ! (10/09/07) - upgrade the finite difference scheme
- ! to a finite impulse response filter
- ! (02/12/09) - OpenMP parallel implementation
+ !> subroutine EquivalentBodyForce
+ !! computes and updates the equivalent body-force
+ !!
+ !! f = - div.( C : E^i )
+ !!
+ !! and the equivalent surface traction
+ !!
+ !! t = n . C : E^i
+ !!
+ !! with n = (0,0,-1). In indicial notations
+ !!
+ !! f_i = - (C_ijkl E^i_kl),j
+ !!
+ !! and
+ !!
+ !! t_1 = n_j C_ijkl E^i_kl
+ !!
+ !! where f is the equivalent body-force, t is the equivalent surface
+ !! traction, C is the elastic moduli tensor and E^i is the moment
+ !! density tensor tensor.
+ !!
+ !! Divergence is computed with a mixed numerical scheme including
+ !! centered finite-difference (in the vertical direction) and
+ !! finite impulse response differentiator filter for derivatives
+ !! estimates. see function 'stress' for further explanations.
+ !!
+ !! \author sylvain barbot (07/09/07) - original form
+ !! (10/09/07) - upgrade the finite difference scheme
+ !! to a finite impulse response filter
+ !! (02/12/09) - OpenMP parallel implementation
!-----------------------------------------------------------------
SUBROUTINE equivalentbodyforce(sig,dx1,dx2,dx3,sx1,sx2,sx3, &
c1,c2,c3,t1,t2,t3,mask)
@@ -1436,12 +1439,12 @@ CONTAINS
!---------------------------------------------------------------------
- ! function SourceSpectrum
- ! computes the equivalent body-forces for a buried dislocation,
- ! with strike-slip and dip-slip components,
- ! slip s, width W, length L in a rigidity mu
- !
- ! sylvain barbot (06-25-07) - original form
+ !> function SourceSpectrum
+ !! computes the equivalent body-forces for a buried dislocation,
+ !! with strike-slip and dip-slip components,
+ !! slip s, width W, length L in a rigidity mu
+ !!
+ !! \author sylvain barbot (06-25-07) - original form
!---------------------------------------------------------------------
SUBROUTINE sourcespectrum(mu,s,x,y,d, &
L,W,strike,dip,rake,beta,dx1,dx2,dx3,f1,f2,f3)
@@ -1524,12 +1527,12 @@ CONTAINS
!---------------------------------------------------------------------
- ! function SourceSpectrumHalfSpace
- ! computes the equivalent body-forces for a buried dislocation,
- ! with strike-slip and dip-slip components,
- ! slip s, width W, length L in a rigidity mu; sources are not imaged
- !
- ! sylvain barbot (06-25-07) - original form
+ !> function SourceSpectrumHalfSpace
+ !! computes the equivalent body-forces for a buried dislocation,
+ !! with strike-slip and dip-slip components,
+ !! slip s, width W, length L in a rigidity mu; sources are not imaged
+ !!
+ !! \author sylvain barbot (06-25-07) - original form
!---------------------------------------------------------------------
SUBROUTINE sourcespectrumhalfspace(mu,s,x,y,d, &
L,W,strike,dip,rake,beta,dx1,dx2,dx3,f1,f2,f3)
@@ -1600,23 +1603,27 @@ CONTAINS
END SUBROUTINE sourcespectrumhalfspace
!---------------------------------------------------------------------
- ! function Source computes the equivalent body-forces
- ! in the space domain for a buried dislocation with strike-slip
- ! and dip-slip components, slip s, width W, length L in a rigidity mu
- !
- ! Default (strike=0, dip=0, rake=0) is a vertical left-lateral
- ! strike-slip fault along the x2 axis. Default fault slip is
- ! represented with the double-couple equivalent body forces:
- !
- ! x1
- ! |
- ! | ^ f2
- ! | |<-----
- ! +---+------+---- x2
- ! ----->|
- ! v f1
- !
- ! sylvain barbot (06-29-07) - original form
+ !> function Source computes the equivalent body-forces
+ !! in the space domain for a buried dislocation with strike-slip
+ !! and dip-slip components, slip s, width W, length L in a rigidity mu
+ !!
+ !! Default (strike=0, dip=0, rake=0) is a vertical left-lateral
+ !! strike-slip fault along the x2 axis. Default fault slip is
+ !! represented with the double-couple equivalent body forces:
+ !!
+ !!\verbatim
+ !!
+ !! x1
+ !! |
+ !! | ^ f2
+ !! | |<-----
+ !! +---+------+---- x2
+ !! ----->|
+ !! v f1
+ !!
+ !!\endverbatim
+ !!
+ !! \author sylvain barbot (06-29-07) - original form
!---------------------------------------------------------------------
SUBROUTINE source(mu,s,x,y,z,L,W,strike,dip,rake, &
beta,sx1,sx2,sx3,dx1,dx2,dx3,f1,f2,f3,t1,t2,t3)
@@ -1802,53 +1809,65 @@ CONTAINS
END SUBROUTINE source
!---------------------------------------------------------------------
- ! function TensileSource
- ! computes the equivalent body-forces in the space domain for a buried
- ! tensile crack with opening s, width W, length L and Lame parameters
- ! lambda, mu.
- !
- ! Default (strike=0, dip=0) is a vertical opening along the x2 axis.
- ! Default fault opening is represented with the double-couple
- ! equivalent body forces:
- !
- ! x1 f1
- ! | ^^^^^^^
- ! | |||||||
- ! | -f2 <--+-------+--> f2
- ! | |||||||
- ! | vvvvvvv
- ! | -f1
- ! |
- ! +----------------------------- x2
- !
- ! The eigenstrain/potency tensor for a point source is
- !
- ! | 1 0 0 |
- ! E^i = | 0 0 0 |
- ! | 0 0 0 |
- !
- ! and the corresponding moment density for a point source is
- !
- ! | lambda+2*mu 0 0 |
- ! m = C : E^i = | 0 lambda 0 |
- ! | 0 0 lambda |
- !
- ! Moment density is integrated along the planar surface
- !
- ! box(x2) delta (x1) box(x3)
- !
- ! where box(x) and delta(x) are the boxcar and the dirac delta
- ! functions, respectively. We use a tapered boxcar, omega_beta(x) and
- ! approximate the delta function by a small gaussian function.
- ! Finally, the equivalent body force is the divergence of the moment
- ! density tensor
- !
- ! f_i = - ( m_ij ),j
- !
- ! derivatives are performed analytically on the gaussian and
- ! omega_beta functions.
- !
- ! sylvain barbot (05-09-08) - original form
+ !> function TensileSource
+ !! computes the equivalent body-forces in the space domain for a buried
+ !! tensile crack with opening s, width W, length L and Lame parameters
+ !! lambda, mu.
+ !!
+ !! Default (strike=0, dip=0) is a vertical opening along the x2 axis.
+ !! Default fault opening is represented with the double-couple
+ !! equivalent body forces:
+ !!
+ !!\verbatim
+ !!
+ !! x1 f1
+ !! | ^^^^^^^
+ !! | |||||||
+ !! | -f2 <--+-------+--> f2
+ !! | |||||||
+ !! | vvvvvvv
+ !! | -f1
+ !! |
+ !! +----------------------------- x2
+ !!
+ !!\endverbatim
+ !!
+ !! The eigenstrain/potency tensor for a point source is
+ !!
+ !!\verbatim
+ !!
+ !! | 1 0 0 |
+ !! E^i = | 0 0 0 |
+ !! | 0 0 0 |
+ !!
+ !!\endverbatim
+ !!
+ !! and the corresponding moment density for a point source is
+ !!
+ !!\verbatim
+ !!
+ !! | lambda+2*mu 0 0 |
+ !! m = C : E^i = | 0 lambda 0 |
+ !! | 0 0 lambda |
+ !!
+ !!\endverbatim
+ !!
+ !! Moment density is integrated along the planar surface
+ !!
+ !! \f[ box(x2) \delta(x1) box(x3) \f]
+ !!
+ !! where box(x) and delta(x) are the boxcar and the dirac delta
+ !! functions, respectively. We use a tapered boxcar, omega_beta(x) and
+ !! approximate the delta function by a small gaussian function.
+ !! Finally, the equivalent body force is the divergence of the moment
+ !! density tensor
+ !!
+ !! \f[ f_i = - ( m_{ij} )_{,j} \f]
+ !!
+ !! derivatives are performed analytically on the gaussian and
+ !! omega_beta functions.
+ !!
+ !! \author sylvain barbot (05-09-08) - original form
!---------------------------------------------------------------------
SUBROUTINE tensilesource(lambda,mu,s,x,y,z,L,W,strike,dip, &
beta,sx1,sx2,sx3,dx1,dx2,dx3,f1,f2,f3)
@@ -1949,45 +1968,49 @@ CONTAINS
END SUBROUTINE tensilesource
!---------------------------------------------------------------------
- ! function MogiSource
- ! computes the equivalent body-forces in the space domain for a buried
- ! dilatation point source.
- !
- ! The point-source opening o with at position xs in the half space is
- ! associated with eigenstrain
- !
- ! E^i = o 1/3 I delta(x-xs)
- !
- ! where I is the diagonal tensor and delta is the Dirac delta function
- ! (or in index notation E^i_{ij} = o delta_{ij} / 3 delta(xs) ) and
- ! with the moment density
- !
- ! m = C : E^i = K o I delta(x-xs)
- !
- ! The equivalent body-force density is
- !
- ! f = - Nabla . m = K o nabla delta(x-xs)
- !
- ! where nabla is the gradient operator. Default source opening is
- ! represented with the isotropic equivalent body-force density:
- !
- ! x1
- ! | f1
- ! | ^
- ! | f2 | f2
- ! +---<--+-->---- x2
- ! |
- ! v f1
- !
- ! x3
- ! | f3
- ! | ^
- ! | f2 | f2
- ! +---<--+-->---- x2
- ! |
- ! v f3
- !
- ! sylvain barbot (03-24-09) - original form
+ !! function MogiSource
+ !! computes the equivalent body-forces in the space domain for a buried
+ !! dilatation point source.
+ !!
+ !! The point-source opening o with at position xs in the half space is
+ !! associated with eigenstrain
+ !!
+ !! \f[ E^i = o \frac{1}{3} I \delta(x-x_s) \f]
+ !!
+ !! where I is the diagonal tensor and delta is the Dirac delta function
+ !! (or in index notation E^i_{ij} = o delta_{ij} / 3 delta(xs) ) and
+ !! with the moment density
+ !!
+ !! \f[ m = C : E^i = K o I \delta(x-x_s) \f]
+ !!
+ !! The equivalent body-force density is
+ !!
+ !! \f[ f = - \nabla \cdot m = K o \nabla \delta(x-x_s) \f]
+ !!
+ !! where nabla is the gradient operator. Default source opening is
+ !! represented with the isotropic equivalent body-force density:
+ !!
+ !!\verbatim
+ !!
+ !! x1
+ !! | f1
+ !! | ^
+ !! | f2 | f2
+ !! +---<--+-->---- x2
+ !! |
+ !! v f1
+ !!
+ !! x3
+ !! | f3
+ !! | ^
+ !! | f2 | f2
+ !! +---<--+-->---- x2
+ !! |
+ !! v f3
+ !!
+ !!\endverbatim
+ !!
+ !! \author sylvain barbot (03-24-09) - original form
!---------------------------------------------------------------------
SUBROUTINE mogisource(lambda,mu,o,xs,ys,zs,sx1,sx2,sx3,dx1,dx2,dx3,f1,f2,f3)
INTEGER, INTENT(IN) :: sx1,sx2,sx3
@@ -2043,61 +2066,73 @@ CONTAINS
END SUBROUTINE mogisource
!---------------------------------------------------------------------
- ! function MomentDensityShear
- ! computes the inelastic irreversible moment density in the space
- ! domain corresponding to a buried dislocation with strike-slip and
- ! dip-slip components (pure shear). A fault along a surface of normal
- ! n_i with a burger vector s_i, is associated with the eigenstrain
- !
- ! E^i_ij = 1/2 ( n_i s_j + s_i n_j )
- !
- ! In a heterogeneous medium of elastic moduli tensor C_ijkl, the
- ! corresponding moment density tensor is
- !
- ! m_ij = C_ijkl E^i_kl
- !
- ! where C = C(x) is a function of space. Equivalent body forces
- ! representing the set of dislocations can be obtained by evaluating
- ! the divergence of the moment density tensor
- !
- ! f_i = - ( m_ji ),j
- !
- ! using the function "EquivalentBodyForce" in this module.
- !
- ! The default dislocation extends in the x2 direction, with a normal
- ! in the x1 direction. Using the following angular convention,
- !
- ! x1 ! x1
- ! n theta | ! n phi |
- ! \ ____| ! \ ____|
- ! \ | ! \ |
- ! \ | ! \ |
- ! -----\+------ x2 ! -----\+------ x3
- ! (x3 down) ! (x2 up)
- !
- ! where theta is the strike and phi is the dip (internal convention),
- ! and introducting the rotation matrices
- !
- ! | cos(theta) sin(theta) 0 |
- ! R1 = | -sin(theta) cos(theta) 0 |
- ! | 0 0 1 |
- !
- ! | cos(phi) 0 sin(phi) |
- ! R2 = | 0 1 0 |
- ! | -sin(phi) 0 cos(phi) |
- !
- ! a normal vector n of arbitrary orientation and the corresponding
- ! strike-slip and dip-slip vector, s and d respectively, are
- !
- ! | 1 | | 0 | | 0 |
- ! n = R1 R2 | 0 |, s = R1 R2 | 1 |, d = R1 R2 | 0 |
- ! | 0 | | 0 | | 1 |
- !
- ! vector n, s and d are orthogonal and the corresponding moment
- ! density second order tensor is deviatoric. The method of images is
- ! used to avoid tapering of the fault at the surface.
- !
- ! sylvain barbot (03-02-08) - original form
+ !! function MomentDensityShear
+ !! computes the inelastic irreversible moment density in the space
+ !! domain corresponding to a buried dislocation with strike-slip and
+ !! dip-slip components (pure shear). A fault along a surface of normal
+ !! n_i with a burger vector s_i, is associated with the eigenstrain
+ !!
+ !! E^i_ij = 1/2 ( n_i s_j + s_i n_j )
+ !!
+ !! In a heterogeneous medium of elastic moduli tensor C_ijkl, the
+ !! corresponding moment density tensor is
+ !!
+ !! m_ij = C_ijkl E^i_kl
+ !!
+ !! where C = C(x) is a function of space. Equivalent body forces
+ !! representing the set of dislocations can be obtained by evaluating
+ !! the divergence of the moment density tensor
+ !!
+ !! f_i = - ( m_ji ),j
+ !!
+ !! using the function "EquivalentBodyForce" in this module.
+ !!
+ !! The default dislocation extends in the x2 direction, with a normal
+ !! in the x1 direction. Using the following angular convention,
+ !!
+ !!\verbatim
+ !!
+ !! x1 ! x1
+ !! n theta | ! n phi |
+ !! \ ____| ! \ ____|
+ !! \ | ! \ |
+ !! \ | ! \ |
+ !! -----\+------ x2 ! -----\+------ x3
+ !! (x3 down) ! (x2 up)
+ !!
+ !!\endverbatim
+ !!
+ !! where theta is the strike and phi is the dip (internal convention),
+ !! and introducting the rotation matrices
+ !!
+ !!\verbatim
+ !!
+ !! | cos(theta) sin(theta) 0 |
+ !! R1 = | -sin(theta) cos(theta) 0 |
+ !! | 0 0 1 |
+ !!
+ !! | cos(phi) 0 sin(phi) |
+ !! R2 = | 0 1 0 |
+ !! | -sin(phi) 0 cos(phi) |
+ !!
+ !!\endverbatim
+ !!
+ !! a normal vector n of arbitrary orientation and the corresponding
+ !! strike-slip and dip-slip vector, s and d respectively, are
+ !!
+ !!\verbatim
+ !!
+ !! | 1 | | 0 | | 0 |
+ !! n = R1 R2 | 0 |, s = R1 R2 | 1 |, d = R1 R2 | 0 |
+ !! | 0 | | 0 | | 1 |
+ !!
+ !!\endverbatim
+ !!
+ !! vector n, s and d are orthogonal and the corresponding moment
+ !! density second order tensor is deviatoric. The method of images is
+ !! used to avoid tapering of the fault at the surface.
+ !!
+ !! \author sylvain barbot (03-02-08) - original form
!---------------------------------------------------------------------
SUBROUTINE momentdensityshear(mu,slip,x,y,z,L,W,strike,dip,rake, &
beta,sx1,sx2,sx3,dx1,dx2,dx3,sig)
@@ -2191,88 +2226,107 @@ CONTAINS
END SUBROUTINE momentdensityshear
!---------------------------------------------------------------------
- ! function MomentDensityTensile
- ! computes the inelastic irreversible moment density in the space
- ! domain corresponding to a buried dislocation with opening (open
- ! crack). A fault along a surface of normal n_i with a burger vector
- ! s_i, is associated with the eigenstrain
- !
- ! E^i_ij = 1/2 ( n_i s_j + s_i n_j )
- !
- ! The eigenstrain/potency tensor for a point source opening crack is
- !
- ! | 1 0 0 |
- ! E^i = | 0 0 0 |
- ! | 0 0 0 |
-
- !
- ! In a heterogeneous medium of elastic moduli tensor C_ijkl, the
- ! corresponding moment density tensor is
- !
- ! m_ij = C_ijkl E^i_kl = lambda E^i_kk delta_ij + 2 mu E^i_ij
- !
- ! where C = C(x) is a function of space. (We use isotropic elastic
- ! solid, and heterogeneous elastic moduli tensor simplifies to
- ! mu=mu(x) and lambda = lambda(x).) The moment density for a point
- ! source opening crack is
- !
- ! | lambda+2*mu 0 0 |
- ! m(x) = | 0 lambda 0 |
- ! | 0 0 lambda |
- !
- ! Moment density m(x) is integrated along the planar surface
- !
- ! box(x2) delta (x1) box(x3)
- !
- ! where box(x) and delta(x) are the boxcar and the dirac delta
- ! functions, respectively. Equivalent body forces representing the
- ! set of dislocations can be obtained by evaluating the divergence
- ! of the moment density tensor
- !
- ! f_i = - ( m_ji ),j
- !
- ! The corresponding equivalent surface traction is simply
- !
- ! t_i = m_ij n_j
- !
- ! Both equivalent body forces and equivalent surface traction are
- ! computed using the function "EquivalentBodyForce" in this module.
- !
- ! The default dislocation extends in the x2 direction, with a normal
- ! in the x1 direction. Using the following angular convention,
- !
- ! x1 ! x1
- ! n theta | ! n phi |
- ! \ ____| ! \ ____|
- ! \ | ! \ |
- ! \ | ! \ |
- ! -----\+------ x2 ! -----\+------ x3
- ! (x3 down) ! (x2 up)
- !
- ! where theta is the strike and phi is the dip, in internal
- ! convention. (Internal angular convention does not correspond to
- ! usual angular convention of geology and conversion between the two
- ! standard is necessary.) Introducting the rotation matrices,
- !
- ! | cos(theta) sin(theta) 0 |
- ! R1 = | -sin(theta) cos(theta) 0 |
- ! | 0 0 1 |
- !
- ! | cos(phi) 0 sin(phi) |
- ! R2 = | 0 1 0 |
- ! | -sin(phi) 0 cos(phi) |
- !
- ! a normal vector n of arbitrary orientation and the corresponding
- ! slip vector s are
- !
- ! | 1 | | 1 |
- ! n = R1 R2 | 0 |, s = n = R1 R2 | 0 |
- ! | 0 | | 0 |
- !
- ! The method of images is used to avoid tapering of the fault at
- ! the surface.
- !
- ! sylvain barbot (03-02-08) - original form
+ !> function MomentDensityTensile
+ !! computes the inelastic irreversible moment density in the space
+ !! domain corresponding to a buried dislocation with opening (open
+ !! crack). A fault along a surface of normal n_i with a burger vector
+ !! s_i, is associated with the eigenstrain
+ !!
+ !! \f[ E^i_{ij} = \frac{1}{2} ( n_i s_j + s_i n_j ) \f]
+ !!
+ !! The eigenstrain/potency tensor for a point source opening crack is
+ !!
+ !!\verbatim
+ !!
+ !! | 1 0 0 |
+ !! E^i = | 0 0 0 |
+ !! | 0 0 0 |
+ !!
+ !!\endverbatim
+ !!
+ !! In a heterogeneous medium of elastic moduli tensor C_ijkl, the
+ !! corresponding moment density tensor is
+ !!
+ !! \f[ m_{ij} = C_{ijkl} E^i_{kl} = \lambda E^i_{kk} \delta_{ij} + 2 \mu E^i_{ij} \f]
+ !!
+ !! where C = C(x) is a function of space. (We use isotropic elastic
+ !! solid, and heterogeneous elastic moduli tensor simplifies to
+ !! mu=mu(x) and lambda = lambda(x).) The moment density for a point
+ !! source opening crack is
+ !!
+ !!\verbatim
+ !!
+ !! | lambda+2*mu 0 0 |
+ !! m(x) = | 0 lambda 0 |
+ !! | 0 0 lambda |
+ !!
+ !!\endverbatim
+ !!
+ !! Moment density m(x) is integrated along the planar surface
+ !!
+ !! box(x2) delta (x1) box(x3)
+ !!
+ !! where box(x) and delta(x) are the boxcar and the dirac delta
+ !! functions, respectively. Equivalent body forces representing the
+ !! set of dislocations can be obtained by evaluating the divergence
+ !! of the moment density tensor
+ !!
+ !! \f[ f_i = - ( m_{ji} ),j \f]
+ !!
+ !! The corresponding equivalent surface traction is simply
+ !!
+ !! \f[ t_i = m_{ij} n_j \f]
+ !!
+ !! Both equivalent body forces and equivalent surface traction are
+ !! computed using the function "EquivalentBodyForce" in this module.
+ !!
+ !! The default dislocation extends in the x2 direction, with a normal
+ !! in the x1 direction. Using the following angular convention,
+ !!
+ !!\verbatim
+ !!
+ !! x1 ! x1
+ !! n theta | ! n phi |
+ !! \ ____| ! \ ____|
+ !! \ | ! \ |
+ !! \ | ! \ |
+ !! -----\+------ x2 ! -----\+------ x3
+ !! (x3 down) ! (x2 up)
+ !!
+ !!\endverbatim
+ !!
+ !! where theta is the strike and phi is the dip, in internal
+ !! convention. (Internal angular convention does not correspond to
+ !! usual angular convention of geology and conversion between the two
+ !! standard is necessary.) Introducting the rotation matrices,
+ !!
+ !!\verbatim
+ !!
+ !! | cos(theta) sin(theta) 0 |
+ !! R1 = | -sin(theta) cos(theta) 0 |
+ !! | 0 0 1 |
+ !!
+ !! | cos(phi) 0 sin(phi) |
+ !! R2 = | 0 1 0 |
+ !! | -sin(phi) 0 cos(phi) |
+ !!
+ !!\endverbatim
+ !!
+ !! a normal vector n of arbitrary orientation and the corresponding
+ !! slip vector s are
+ !!
+ !!\verbatim
+ !!
+ !! | 1 | | 1 |
+ !! n = R1 R2 | 0 |, s = n = R1 R2 | 0 |
+ !! | 0 | | 0 |
+ !!
+ !!\endverbatim
+ !!
+ !! The method of images is used to avoid tapering of the fault at
+ !! the surface.
+ !!
+ !! \author sylvain barbot (03-02-08) - original form
!---------------------------------------------------------------------
SUBROUTINE momentdensitytensile(lambda,mu,slip,x,y,z,L,W,strike,dip,rake, &
beta,sx1,sx2,sx3,dx1,dx2,dx3,sig)
@@ -2355,27 +2409,27 @@ CONTAINS
END SUBROUTINE momentdensitytensile
!---------------------------------------------------------------------
- ! function MomentDensityMogi
- ! computes the inelastic irreversible moment density in the space
- ! domain corresponding to a buried Mogi source.
- ! The Mogi source is associated with the eigenstrain
- !
- ! E^i_ij = o 1/3 delta_ij
- !
- ! In a heterogeneous medium of elastic moduli tensor C_ijkl, the
- ! corresponding moment density tensor is
- !
- ! m_ij = C_ijkl E^i_kl
- !
- ! where C = C(x) is a function of space. Equivalent body forces
- ! representing the set of dislocations can be obtained by evaluating
- ! the divergence of the moment density tensor
- !
- ! f_i = - ( m_ji ),j
- !
- ! using the function "EquivalentBodyForce" in this module.
- !
- ! sylvain barbot (03-24-09) - original form
+ !! function MomentDensityMogi
+ !! computes the inelastic irreversible moment density in the space
+ !! domain corresponding to a buried Mogi source.
+ !! The Mogi source is associated with the eigenstrain
+ !!
+ !! \f[ E^i_{ij} = o \frac{1}{3} \delta_{ij} \f]
+ !!
+ !! In a heterogeneous medium of elastic moduli tensor C_ijkl, the
+ !! corresponding moment density tensor is
+ !!
+ !! \f[ m_{ij} = C_{ijkl} E^i_{kl} \f]
+ !!
+ !! where C = C(x) is a function of space. Equivalent body forces
+ !! representing the set of dislocations can be obtained by evaluating
+ !! the divergence of the moment density tensor
+ !!
+ !! \f[ f_i = - ( m_{ji} ),j \f]
+ !!
+ !! using the function "EquivalentBodyForce" in this module.
+ !!
+ !! \author sylvain barbot (03-24-09) - original form
!---------------------------------------------------------------------
SUBROUTINE momentdensitymogi(lambda,mu,o,xs,ys,zs,sx1,sx2,sx3,dx1,dx2,dx3,sig)
INTEGER, INTENT(IN) :: sx1,sx2,sx3
@@ -2426,36 +2480,44 @@ CONTAINS
END SUBROUTINE momentdensitymogi
!---------------------------------------------------------------------
- ! function Plane
- ! computes the three components, n1, n2 and n3, of the normal vector
- ! corresponding to a rectangular surface of finite size. The plane
- ! is defined by its orientation (strike and dip) and dimension.
- !
- ! W
- ! +-------------+
- ! | |
- ! L | + | - - - > along strike direction
- ! | (x,y,z) |
- ! +-------------|
- ! |
- ! v
- ! down-dip direction
- !
- ! in the default orientation, for which strike=0 and dip=0, the plane
- ! is vertical along the x2 axis, such as n2(x) = n3(x) = 0 for all x.
- ! internal angular conventions are as follows:
- !
- ! n x1 n x1
- ! \ | \ |
- ! \ | \ |
- ! 90 - strike \ | 90 - dip \ |
- ! ( \| ( \|
- ! ----------+------ x2 ----------+------ x3
- ! (x3 down) (x2 up)
- !
- ! edges of the rectangle are tapered.
- !
- ! sylvain barbot (09-15-07) - original form
+ !> function Plane
+ !! computes the three components, n1, n2 and n3, of the normal vector
+ !! corresponding to a rectangular surface of finite size. The plane
+ !! is defined by its orientation (strike and dip) and dimension.
+ !!
+ !!\verbatim
+ !!
+ !! W
+ !! +-------------+
+ !! | |
+ !! L | + | - - - > along strike direction
+ !! | (x,y,z) |
+ !! +-------------|
+ !! |
+ !! v
+ !! down-dip direction
+ !!
+ !!\endverbatim
+ !!
+ !! in the default orientation, for which strike=0 and dip=0, the plane
+ !! is vertical along the x2 axis, such as n2(x) = n3(x) = 0 for all x.
+ !! internal angular conventions are as follows:
+ !!
+ !!\verbatim
+ !!
+ !! n x1 n x1
+ !! \ | \ |
+ !! \ | \ |
+ !! 90 - strike \ | 90 - dip \ |
+ !! ( \| ( \|
+ !! ----------+------ x2 ----------+------ x3
+ !! (x3 down) (x2 up)
+ !!
+ !!\endverbatim
+ !!
+ !! edges of the rectangle are tapered.
+ !!
+ !! \author sylvain barbot (09-15-07) - original form
!---------------------------------------------------------------------
SUBROUTINE plane(x,y,z,L,W,strike,dip, &
beta,sx1,sx2,sx3,dx1,dx2,dx3,n1,n2,n3)
@@ -2528,10 +2590,10 @@ CONTAINS
END SUBROUTINE plane
!---------------------------------------------------------------------
- ! function MonitorField
- ! samples a scalar field along a specified planar surface.
- !
- ! sylvain barbot (10-16-07) - original form
+ !> function MonitorField
+ !! samples a scalar field along a specified planar surface.
+ !!
+ !! \author sylvain barbot (10-16-07) - original form
!---------------------------------------------------------------------
SUBROUTINE monitorfield(x,y,z,L,W,strike,dip,beta, &
sx1,sx2,sx3,dx1,dx2,dx3,slip,patch)
@@ -2591,13 +2653,13 @@ CONTAINS
CONTAINS
!--------------------------------------------------------------
- ! subroutine sample
- ! interpolates the value of a discretized 3-dimensional field
- ! at a subpixel location. method consists in correlating the
- ! 3D field with a delta function filter. the delta function is
- ! approximated with a narrow normalized gaussian.
- !
- ! sylvain barbot (10-17-07) - original form
+ !> subroutine sample
+ !! interpolates the value of a discretized 3-dimensional field
+ !! at a subpixel location. method consists in correlating the
+ !! 3D field with a delta function filter. the delta function is
+ !! approximated with a narrow normalized gaussian.
+ !!
+ !! \author sylvain barbot (10-17-07) - original form
!--------------------------------------------------------------
SUBROUTINE sample(x1,x2,x3,dx1,dx2,dx3,sx1,sx2,sx3,field,value)
INTEGER, INTENT(IN) :: sx1,sx2,sx3
@@ -2799,13 +2861,13 @@ CONTAINS
END SUBROUTINE sliceadd
!-----------------------------------------------------------------
- ! subroutine TensorFieldAdd
- ! computes the linear combination of two tensor fields
- !
- ! t1 = c1 * t1 + c2 * t2
- !
- ! where t1 and t2 are two tensor fields and c1 and c2 are scalars.
- ! only tensor field t1 is modified.
+ !> subroutine TensorFieldAdd
+ !! computes the linear combination of two tensor fields
+ !!
+ !! t1 = c1 * t1 + c2 * t2
+ !!
+ !! where t1 and t2 are two tensor fields and c1 and c2 are scalars.
+ !! only tensor field t1 is modified.
!
! sylvain barbot (07/27/07) - original form
!-----------------------------------------------------------------
@@ -2900,8 +2962,8 @@ 50 CONTINUE
END SUBROUTINE tensorintegrate
!---------------------------------------------------------------------
- ! subroutine coordinates computes the xi coordinates from the
- ! array index and sampling interval
+ !> subroutine coordinates computes the xi coordinates from the
+ !! array index and sampling interval
!---------------------------------------------------------------------
SUBROUTINE coordinates(i1,i2,i3,sx1,sx2,sx3,dx1,dx2,dx3,x1,x2,x3)
INTEGER, INTENT(IN) :: i1,i2,i3,sx1,sx2,sx3
@@ -2914,11 +2976,11 @@ 50 CONTINUE
END SUBROUTINE coordinates
!---------------------------------------------------------------------
- ! subroutine ShiftedCoordinates
- ! computes the xi coordinates from the array index and sampling
- ! interval assuming data is order like fftshift.
- !
- ! sylvain barbot (07/31/07) - original form
+ !> subroutine ShiftedCoordinates
+ !! computes the xi coordinates from the array index and sampling
+ !! interval assuming data is order like fftshift.
+ !!
+ !! \author sylvain barbot (07/31/07) - original form
!---------------------------------------------------------------------
SUBROUTINE shiftedcoordinates(i1,i2,i3,sx1,sx2,sx3,dx1,dx2,dx3,x1,x2,x3)
INTEGER, INTENT(IN) :: i1,i2,i3,sx1,sx2,sx3
@@ -2944,17 +3006,17 @@ 50 CONTINUE
END SUBROUTINE shiftedcoordinates
!----------------------------------------------------------------------
- ! subroutine ShiftedIndex
- ! returns the integer index corresponding to the specified coordinates
- ! assuming the data are ordered following fftshift. input coordinates
- ! are assumed bounded -sx/2 <= x <= sx/2-1. out of bound input
- ! purposefully triggers a fatal error. in the x3 direction, coordinates
- ! are assumed bounded by 0 <= x3 <= (sx3-1)*dx3
- !
- ! CALLED BY:
- ! monitorfield/sample
- !
- ! sylvain barbot (07/31/07) - original form
+ !> subroutine ShiftedIndex
+ !! returns the integer index corresponding to the specified coordinates
+ !! assuming the data are ordered following fftshift. input coordinates
+ !! are assumed bounded -sx/2 <= x <= sx/2-1. out of bound input
+ !! purposefully triggers a fatal error. in the x3 direction, coordinates
+ !! are assumed bounded by 0 <= x3 <= (sx3-1)*dx3
+ !!
+ !! CALLED BY:
+ !! monitorfield/sample
+ !!
+ !! \author sylvain barbot (07/31/07) - original form
!----------------------------------------------------------------------
SUBROUTINE shiftedindex(x1,x2,x3,sx1,sx2,sx3,dx1,dx2,dx3,i1,i2,i3)
REAL*8, INTENT(IN) :: x1,x2,x3,dx1,dx2,dx3
@@ -3050,11 +3112,11 @@ 50 CONTINUE
END SUBROUTINE exportslice
!-----------------------------------------------------------------
- ! subroutine ExportSpatial
- ! transfer a horizontal layer from array 'data' to smaller array
- ! 'p' and shift center position so that coordinates (0,0) are in
- ! center of array 'p'. optional parameter 'doflip' generates
- ! output compatible with grd binary format.
+ !> subroutine ExportSpatial
+ !! transfer a horizontal layer from array 'data' to smaller array
+ !! 'p' and shift center position so that coordinates (0,0) are in
+ !! center of array 'p'. optional parameter 'doflip' generates
+ !! output compatible with grd binary format.
!
! sylvain barbot (07/09/07) - original form
! (03/19/08) - compatibility with grd output
diff -r 777cc20187d0 -r 0a07ece6fb5a export.f90
--- a/export.f90 Thu May 12 17:04:47 2011 -0700
+++ b/export.f90 Thu May 12 19:01:34 2011 -0700
@@ -171,21 +171,21 @@ CONTAINS
#ifdef PROJ
!------------------------------------------------------------------
- ! subroutine ExportStressPROJ
- ! export a map view of stress with coordinates in
- ! longitude/latitude. Text format output is the GMT-compatible
- ! .xyz file format where data in each file is organized as follows
- !
- ! longitude latitude s11
- ! longitude latitude s12
- ! longitude latitude s13
- ! longitude latitude s22
- ! longitude latitude s23
- ! longitude latitude s33
- !
- ! this is an interface to exportproj.
- !
- ! sylvain barbot (05/22/10) - original form
+ !> subroutine ExportStressPROJ
+ !! export a map view of stress with coordinates in
+ !! longitude/latitude. Text format output is the GMT-compatible
+ !! .xyz file format where data in each file is organized as follows
+ !!
+ !! longitude latitude s11
+ !! longitude latitude s12
+ !! longitude latitude s13
+ !! longitude latitude s22
+ !! longitude latitude s23
+ !! longitude latitude s33
+ !!
+ !! this is an interface to exportproj.
+ !!
+ !! \author sylvain barbot (05/22/10) - original form
!------------------------------------------------------------------
SUBROUTINE exportstressproj(sig,sx1,sx2,sx3,dx1,dx2,dx3,oz, &
x0,y0,lon0,lat0,zone,scale,wdir,index)
@@ -238,20 +238,20 @@ CONTAINS
END SUBROUTINE exportstressproj
!------------------------------------------------------------------
- ! subroutine ExportPROJ
- ! export a map view of displacements with coordinates in
- ! longitude/latitude. Text format output is the GMT-compatible
- ! .xyz file format where data in each file is organized as follows
- !
- ! longitude latitude u1,
- ! longitude latitude u2 and
- ! longitude latitude -u3
- !
- ! for index-geo-north.xyz,
- ! index-geo-east.xyz and
- ! index-geo-up.xyz, respectively.
- !
- ! sylvain barbot (05/22/10) - original form
+ !> subroutine ExportPROJ
+ !! export a map view of displacements with coordinates in
+ !! longitude/latitude. Text format output is the GMT-compatible
+ !! .xyz file format where data in each file is organized as follows
+ !!
+ !! longitude latitude u1,
+ !! longitude latitude u2 and
+ !! longitude latitude -u3
+ !!
+ !! for index-geo-north.xyz,
+ !! index-geo-east.xyz and
+ !! index-geo-up.xyz, respectively.
+ !!
+ !! \author sylvain barbot (05/22/10) - original form
!------------------------------------------------------------------
SUBROUTINE exportproj(c1,c2,c3,sx1,sx2,sx3,dx1,dx2,dx3,oz, &
x0,y0,lon0,lat0,zone,scale,wdir,i,convention)
@@ -335,16 +335,16 @@ CONTAINS
#ifdef XYZ
!------------------------------------------------------------------
- ! subroutine ExportXYZ
- ! export a map view of surface displacement into the GMT-compatible
- ! .xyz file format where data in each file is organized as follows
- !
- ! x1 x2 u1, x1 x2 u2 and x1 x2 -u3
- !
- ! for index-north.xyz, index-east.xyz and index-up.xyz,
- ! respectively.
- !
- ! sylvain barbot (06/10/09) - original form
+ !> subroutine ExportXYZ
+ !! export a map view of surface displacement into the GMT-compatible
+ !! .xyz file format where data in each file is organized as follows
+ !!
+ !! x1 x2 u1, x1 x2 u2 and x1 x2 -u3
+ !!
+ !! for index-north.xyz, index-east.xyz and index-up.xyz,
+ !! respectively.
+ !!
+ !! \author sylvain barbot (06/10/09) - original form
!------------------------------------------------------------------
SUBROUTINE exportxyz(c1,c2,c3,sx1,sx2,sx3,oz,dx1,dx2,dx3,i,wdir)
INTEGER, INTENT(IN) :: i,sx1,sx2,sx3
@@ -435,16 +435,16 @@ CONTAINS
#endif
!------------------------------------------------------------------
- ! subroutine exportpoints
- ! sample a vector field at a series of points for export.
- ! each location is attributed a file in which the time evolution
- ! of the vector value is listed in the format:
- !
- ! t_0 u(t_0) v(t_0) w(t_0)
- ! t_1 u(t_1) v(t_1) w(t_1)
- ! ...
- !
- ! sylvain barbot (11/10/07) - original form
+ !> subroutine exportpoints
+ !! sample a vector field at a series of points for export.
+ !! each location is attributed a file in which the time evolution
+ !! of the vector value is listed in the format:
+ !!
+ !! t_0 u(t_0) v(t_0) w(t_0)
+ !! t_1 u(t_1) v(t_1) w(t_1)
+ !! ...
+ !!
+ !! \author sylvain barbot (11/10/07) - original form
!------------------------------------------------------------------
SUBROUTINE exportpoints(c1,c2,c3,sx1,sx2,sx3,dx1,dx2,dx3, &
opts,ptsname,time,wdir,isnew,x0,y0,rot)
@@ -542,28 +542,28 @@ CONTAINS
END SUBROUTINE exportpoints
!---------------------------------------------------------------------
- ! subroutine exportEigenstrain
- ! samples the value of an input scalar field at the location of
- ! defined plane (position, strike, dip, length and width).
- !
- ! input variables
- ! field - sampled scalar array
- ! nop - number of observation planes
- ! op - structure of observation planes (position, orientation)
- ! x0, y0 - origin position of coordinate system
- ! dx1,2,3 - sampling size
- ! sx1,2,3 - size of the scalar field
- ! wdir - output directory for writing
- ! i - loop index to suffix file names
- !
- ! creates files
- !
- ! wdir/index.s00001.estrain.txt with TXT_EXPORTEIGENSTRAIN option
- !
- ! wdir/index.s00001.estrain.grd with GRD_EXPORTEIGENSTRAIN option
- !
- ! sylvain barbot (01/01/07) - original form
- ! (02/25/10) - output in TXT and GRD formats
+ !> subroutine exportEigenstrain
+ !! samples the value of an input scalar field at the location of
+ !! defined plane (position, strike, dip, length and width).
+ !!
+ !! input variables
+ !! @param field - sampled scalar array
+ !! @param nop - number of observation planes
+ !! @param op - structure of observation planes (position, orientation)
+ !! @param x0, y0 - origin position of coordinate system
+ !! @param dx1,2,3 - sampling size
+ !! @param sx1,2,3 - size of the scalar field
+ !! @param wdir - output directory for writing
+ !! @param i - loop index to suffix file names
+ !!
+ !! creates files
+ !!
+ !! wdir/index.s00001.estrain.txt with TXT_EXPORTEIGENSTRAIN option
+ !!
+ !! wdir/index.s00001.estrain.grd with GRD_EXPORTEIGENSTRAIN option
+ !!
+ !! \author sylvain barbot (01/01/07) - original form
+ ! (02/25/10) - output in TXT and GRD formats
!---------------------------------------------------------------------
SUBROUTINE exporteigenstrain(field,nop,op,x0,y0,dx1,dx2,dx3,sx1,sx2,sx3,wdir,i)
INTEGER, INTENT(IN) :: nop,sx1,sx2,sx3,i
@@ -662,43 +662,43 @@ END SUBROUTINE exporteigenstrain
END SUBROUTINE exporteigenstrain
!---------------------------------------------------------------------
- ! subroutine exportCreep
- ! evaluates the value of creep velocity at the location of
- ! defined plane (position, strike, dip, length and width).
- !
- ! input variables
- ! np - number of frictional planes
- ! n - array of frictional planes (position, orientation)
- ! structure - array of depth-dependent frictional properties
- ! x0, y0 - origin position of coordinate system
- ! dx1,2,3 - sampling size
- ! sx1,2,3 - size of the stress tensor field
- ! beta - smoothing factor controlling the extent of planes
- ! wdir - output directory for writing
- ! i - loop index to suffix file names
- !
- ! creates files
- !
- ! wdir/index.s00001.creep.txt
- !
- ! containing
- !
- ! x,y,z,x',y',sqrt(vx'^2+vy'^2),vx',vy'
- !
- ! with TXT_EXPORTCREEP option and
- !
- ! wdir/index.s00001.creep-north.grd
- ! wdir/index.s00001.creep-east.grd
- ! wdir/index.s00001.creep-up.grd
- !
- ! with GRD_EXPORTCREEP option where the suffix -north stands for
- ! dip slip, -east for strike slip and -up for amplitude of slip.
- !
- ! file wdir/index.s00001.creep.txt is subsampled by a factor "skip"
- ! compared to the grd files.
- !
- ! sylvain barbot (01/01/07) - original form
- ! (02/25/10) - output in TXT and GRD formats
+ !> subroutine exportCreep
+ !! evaluates the value of creep velocity at the location of
+ !! defined plane (position, strike, dip, length and width).
+ !!
+ !! input variables
+ !! @param np - number of frictional planes
+ !! @param n - array of frictional planes (position, orientation)
+ !! @param structure - array of depth-dependent frictional properties
+ !! @param x0, y0 - origin position of coordinate system
+ !! @param dx1,2,3 - sampling size
+ !! @param sx1,2,3 - size of the stress tensor field
+ !! @param beta - smoothing factor controlling the extent of planes
+ !! @param wdir - output directory for writing
+ !! @param i - loop index to suffix file names
+ !!
+ !! creates files
+ !!
+ !! wdir/index.s00001.creep.txt
+ !!
+ !! containing
+ !!
+ !! x,y,z,x',y',sqrt(vx'^2+vy'^2),vx',vy'
+ !!
+ !! with TXT_EXPORTCREEP option and
+ !!
+ !! wdir/index.s00001.creep-north.grd
+ !! wdir/index.s00001.creep-east.grd
+ !! wdir/index.s00001.creep-up.grd
+ !!
+ !! with GRD_EXPORTCREEP option where the suffix -north stands for
+ !! dip slip, -east for strike slip and -up for amplitude of slip.
+ !!
+ !! file wdir/index.s00001.creep.txt is subsampled by a factor "skip"
+ !! compared to the grd files.
+ !!
+ !! \author sylvain barbot (01/01/07) - original form
+ !! (02/25/10) - output in TXT and GRD formats
!---------------------------------------------------------------------
SUBROUTINE exportcreep(np,n,beta,sig,structure, &
sx1,sx2,sx3,dx1,dx2,dx3,x0,y0,wdir,i)
@@ -802,12 +802,12 @@ END SUBROUTINE exportcreep
#ifdef GRD
!------------------------------------------------------------------
- ! subroutine ExportStressGRD
- ! writes the 6 components of deformation in map view in the GMT
- ! (Generic Mapping Tools) GRD binary format. This is an interface
- ! to exportgrd.
- !
- ! sylvain barbot 03/19/08 - original form
+ !> subroutine ExportStressGRD
+ !! writes the 6 components of deformation in map view in the GMT
+ !! (Generic Mapping Tools) GRD binary format. This is an interface
+ !! to exportgrd.
+ !!
+ !! \author sylvain barbot 03/19/08 - original form
!------------------------------------------------------------------
SUBROUTINE exportstressgrd(sig,sx1,sx2,sx3,dx1,dx2,dx3, &
oz,origx,origy,wdir,index)
@@ -861,11 +861,11 @@ END SUBROUTINE exportcreep
!------------------------------------------------------------------
- ! subroutine ExportGRD
- ! writes the 3 components of deformation in map view in the GMT
- ! (Generic Mapping Tools) GRD binary format.
- !
- ! sylvain barbot 03/19/08 - original form
+ !> subroutine ExportGRD
+ !! writes the 3 components of deformation in map view in the GMT
+ !! (Generic Mapping Tools) GRD binary format.
+ !!
+ !! \author sylvain barbot 03/19/08 - original form
!------------------------------------------------------------------
SUBROUTINE exportgrd(c1,c2,c3,sx1,sx2,sx3,dx1,dx2,dx3,oz,origx,origy,&
wdir,i,convention)
@@ -956,11 +956,11 @@ END SUBROUTINE exportcreep
#ifdef VTK
!------------------------------------------------------------------
- ! subroutine ExportVTK_Grid
- ! creates a .vtp file (in the VTK PolyData XML format) containing
- ! the dimension of the computational grid
- !
- ! sylvain barbot 06/24/09 - original form
+ !> subroutine ExportVTK_Grid
+ !! creates a .vtp file (in the VTK PolyData XML format) containing
+ !! the dimension of the computational grid
+ !!
+ !! \author sylvain barbot 06/24/09 - original form
!------------------------------------------------------------------
SUBROUTINE exportvtk_grid(sx1,sx2,sx3,dx1,dx2,dx3,origx,origy,cgfilename)
INTEGER, INTENT(IN) :: sx1,sx2,sx3
@@ -1024,14 +1024,14 @@ END SUBROUTINE exportcreep
END SUBROUTINE exportvtk_grid
!------------------------------------------------------------------
- ! subroutine ExportXY_RFaults
- ! creates a .xy file (in the GMT closed-polygon format) containing
- ! the rectangular faults. Each fault segemnt is described by a
- ! closed polygon (rectangle) associated with a slip amplitude.
- ! use pxzy with the -Cpalette.cpt -L -M options to color rectangles
- ! by slip.
- !
- ! sylvain barbot 03/05/11 - original form
+ !> subroutine ExportXY_RFaults
+ !! creates a .xy file (in the GMT closed-polygon format) containing
+ !! the rectangular faults. Each fault segemnt is described by a
+ !! closed polygon (rectangle) associated with a slip amplitude.
+ !! use pxzy with the -Cpalette.cpt -L -M options to color rectangles
+ !! by slip.
+ !!
+ !! \author sylvain barbot 03/05/11 - original form
!------------------------------------------------------------------
SUBROUTINE exportxy_rfaults(e,rffilename)
TYPE(EVENT_STRUC), INTENT(IN) :: e
@@ -1101,12 +1101,12 @@ END SUBROUTINE exportcreep
END SUBROUTINE exportxy_rfaults
!------------------------------------------------------------------
- ! subroutine ExportVTK_RFaults
- ! creates a .vtp file (in the VTK PolyData XML format) containing
- ! the rectangular faults. The faults are characterized with a set
- ! of subsegments (rectangles) each associated with a slip vector.
- !
- ! sylvain barbot 06/24/09 - original form
+ !> subroutine ExportVTK_RFaults
+ !! creates a .vtp file (in the VTK PolyData XML format) containing
+ !! the rectangular faults. The faults are characterized with a set
+ !! of subsegments (rectangles) each associated with a slip vector.
+ !!
+ !! \author sylvain barbot 06/24/09 - original form
!------------------------------------------------------------------
SUBROUTINE exportvtk_rfaults(e,rffilename)
TYPE(EVENT_STRUC), INTENT(IN) :: e
@@ -1222,11 +1222,11 @@ END SUBROUTINE exportcreep
END SUBROUTINE exportvtk_rfaults
!------------------------------------------------------------------
- ! subroutine ExportVTK_Rectangle
- ! creates a .vtp file (in the VTK PolyData XML format) containing
- ! a rectangle.
- !
- ! sylvain barbot 06/24/09 - original form
+ !> subroutine ExportVTK_Rectangle
+ !! creates a .vtp file (in the VTK PolyData XML format) containing
+ !! a rectangle.
+ !!
+ !! \author sylvain barbot 06/24/09 - original form
!------------------------------------------------------------------
SUBROUTINE exportvtk_rectangle(x1,x2,x3,L,W,strike,dip,filename)
REAL*8 :: x1,x2,x3,L,W,strike,dip
@@ -1310,11 +1310,11 @@ END SUBROUTINE exportcreep
END SUBROUTINE exportvtk_rectangle
!------------------------------------------------------------------
- ! subroutine ExportVTK_Brick
- ! creates a .vtp file (in the VTK PolyData XML format) containing
- ! a brick (3d rectangle, cuboid).
- !
- ! sylvain barbot 06/24/09 - original form
+ !> subroutine ExportVTK_Brick
+ !! creates a .vtp file (in the VTK PolyData XML format) containing
+ !! a brick (3d rectangle, cuboid).
+ !!
+ !! \author sylvain barbot 06/24/09 - original form
!------------------------------------------------------------------
SUBROUTINE exportvtk_brick(x1,x2,x3,L,W,T,strike,dip,filename)
REAL*8 :: x1,x2,x3,L,W,T,strike,dip
@@ -1443,11 +1443,11 @@ END SUBROUTINE exportcreep
END SUBROUTINE exportvtk_brick
!------------------------------------------------------------------
- ! subroutine ExportVTK_Vectors
- ! creates a .vtr file (in the VTK Rectilinear XML format)
- ! containing a vector field.
- !
- ! sylvain barbot 06/25/09 - original form
+ !> subroutine ExportVTK_Vectors
+ !! creates a .vtr file (in the VTK Rectilinear XML format)
+ !! containing a vector field.
+ !!
+ !! \author sylvain barbot 06/25/09 - original form
!------------------------------------------------------------------
SUBROUTINE exportvtk_vectors(u1,u2,u3,sx1,sx2,sx3,dx1,dx2,dx3,j1,j2,j3,vcfilename)
INTEGER, INTENT(IN) :: sx1,sx2,sx3,j1,j2,j3
@@ -1571,11 +1571,11 @@ END SUBROUTINE exportcreep
END SUBROUTINE exportvtk_vectors
!------------------------------------------------------------------
- ! subroutine ExportVTK_Vectors_Slice
- ! creates a .vtr file (in the VTK Rectilinear XML format)
- ! containing a vector field.
- !
- ! sylvain barbot 06/25/09 - original form
+ !> subroutine ExportVTK_Vectors_Slice
+ !! creates a .vtr file (in the VTK Rectilinear XML format)
+ !! containing a vector field.
+ !!
+ !! \author sylvain barbot 06/25/09 - original form
!------------------------------------------------------------------
SUBROUTINE exportvtk_vectors_slice(u1,u2,u3,sx1,sx2,sx3,dx1,dx2,dx3,oz,j1,j2,vcfilename)
INTEGER, INTENT(IN) :: sx1,sx2,sx3,j1,j2
diff -r 777cc20187d0 -r 0a07ece6fb5a fourier.f90
--- a/fourier.f90 Thu May 12 17:04:47 2011 -0700
+++ b/fourier.f90 Thu May 12 19:01:34 2011 -0700
@@ -38,20 +38,28 @@ CONTAINS
CONTAINS
!---------------------------------------------------------------------
- ! subroutine wavenumbers
- ! computes the values of the wavenumbers
- ! in the sequential order required when using subroutine FOURT
- ! to perform forward and backward inverse transforms.
- !
- ! INPUT
- ! i1 i3 running index in the discrete Fourier domain array
- ! sx1 sx3 number of elements in the 2 directions
- ! dx1 dx3 sampling interval in the 2 directions
- !
- ! OUTPUT
- ! k1 k3 wavenumbers in the 2 direction
- !
- ! sylvain barbot (04-14-07) - original form
+ !> subroutine wavenumbers
+ !! computes the values of the wavenumbers
+ !! in the sequential order required when using subroutine FOURT
+ !! to perform forward and backward inverse transforms.
+ !!
+ !! INPUT
+ !! @param i1 running index in the discrete Fourier domain array
+ !! @param i2 running index in the discrete Fourier domain array
+ !! @param i3 running index in the discrete Fourier domain array
+ !! @param sx1 number of elements in the x1-direction
+ !! @param sx2 number of elements in the x2-direction
+ !! @param sx3 number of elements in the x3-direction
+ !! @param dx1 sampling interval in the x1-direction
+ !! @param dx2 sampling interval in the x2-direction
+ !! @param dx3 sampling interval in the x3-direction
+ !!
+ !! OUTPUT
+ !! @param k1 wavenumber in the x1 direction
+ !! @param k2 wavenumber in the x2 direction
+ !! @param k3 wavenumber in the x3 direction
+ !!
+ !! \author sylvain barbot (04-14-07) - original form
!---------------------------------------------------------------------
SUBROUTINE wavenumbers(i1,i2,i3,sx1,sx2,sx3,dx1,dx2,dx3,k1,k2,k3)
INTEGER, INTENT(IN) :: i1, i2, i3, sx1, sx2, sx3
@@ -147,16 +155,17 @@ CONTAINS
END SUBROUTINE fftshift_tf
!----------------------------------------------------------------------
- ! subroutine FFT3 performs normalized forward and
- ! inverse fourier transforms of real 3d data
+ !> subroutine FFT3 performs normalized forward and
+ !! inverse fourier transforms of real 3d data
!
- ! USES
- ! ctfft (Brenner, 1968) by default
- ! fftw3 (Frigo & Jonhson) with preproc FFTW3 flag
- ! scfft (SGI library) with preproc SGI_FFT flag
- !
- ! for real array the fourier transform returns a sx1/2+1 complex array
- ! and the enough space must be reserved
+ !! USES
+ !! ctfft (Brenner, 1968) by default
+ !! fftw3 (Frigo & Jonhson) with preproc FFTW3 flag
+ !! scfft (SGI library) with preproc SGI_FFT flag
+ !! ctfft (Cooley-Tuckey) by default (slowest FFT)
+ !!
+ !! for real array the fourier transform returns a sx1/2+1 complex array
+ !! and the enough space must be reserved
!----------------------------------------------------------------------
#ifdef FFTW3
!--------------------------------------------------------
@@ -323,16 +332,16 @@ CONTAINS
#endif
#endif
!----------------------------------------------------------------------
- ! subroutine FFT2 performs normalized forward and
- ! inverse fourier transforms of real 2d data
- !
- ! USES subroutine FOURT
- ! ctfft(data,n,ndim,isign,iform,work,nwork)
- ! or
- ! fftw3
- !
- ! for real array the fourier transform returns a sx1/2+1 complex array
- ! and the enough space must be reserved
+ !> subroutine FFT2 performs normalized forward and
+ !! inverse fourier transforms of real 2d data
+ !!
+ !! USES subroutine FOURT
+ !! ctfft(data,n,ndim,isign,iform,work,nwork)
+ !! or
+ !! fftw3
+ !!
+ !! for real array the fourier transform returns a sx1/2+1 complex array
+ !! and the enough space must be reserved
!----------------------------------------------------------------------
#ifdef FFTW3
SUBROUTINE fft2(data,sx1,sx2,dx1,dx2,direction)
@@ -479,14 +488,14 @@ CONTAINS
#endif
!-----------------------------------------------------------------
- ! subroutine FFT1
- ! performs a one dimensional complex to complex Fourier
- ! transform
- !
- ! uses complex DFT ctfft (N. Brenner, 1968) by default
- ! or CCFFT (SGI library) with compile flag SGI_FFT
- !
- ! sylvain barbot (05-02-07) - original form
+ !> subroutine FFT1
+ !! performs a one dimensional complex to complex Fourier
+ !! transform
+ !!
+ !! uses complex DFT ctfft (N. Brenner, 1968) by default
+ !! or CCFFT (SGI library) with compile flag SGI_FFT
+ !!
+ !! \author sylvain barbot (05-02-07) - original form
!-----------------------------------------------------------------
#ifdef SGI_FFT
!------------------------------------------------------
diff -r 777cc20187d0 -r 0a07ece6fb5a friction3d.f90
--- a/friction3d.f90 Thu May 12 17:04:47 2011 -0700
+++ b/friction3d.f90 Thu May 12 19:01:34 2011 -0700
@@ -32,9 +32,9 @@ CONTAINS
CONTAINS
!-----------------------------------------------------------------
- ! subroutine FrictionPlaneExpEigenStress
- !
- ! *** this function is deprecated ***
+ !> subroutine FrictionPlaneExpEigenStress
+ !!
+ !! *** this function is deprecated ***
!
! compute the eigen-stress (forcing moment) to be relaxed by
! rate-dependent inelastic deformation in the case of a frictional
@@ -179,51 +179,51 @@ CONTAINS
END SUBROUTINE frictionplaneeigenstress
!-----------------------------------------------------------------
- ! subroutine FrictionEigenStress
- ! compute the eigen-stress (forcing moment) to be relaxed by
- ! rate-dependent inelastic deformation in the case of a frictional
- ! surface:
- !
- ! sigma^i = C:F:sigma
- !
- ! where C is the elastic moduli tensor, F is the heterogeneous
- ! fluidity moduli tensor and sigma is the instantaneous stress
- ! tensor. for a frictional surface, the eigenstrain-rate is given
- ! by
- !
- ! epsilon^i^dot = F:sigma = gamma^dot R
- !
- ! where gamma^dot is the slip rate (a scalar) and R is the
- ! deviatoric, symmetric, and unitary, tensor:
- !
- ! R_ij = 1/2 ( t_i n_j + t_j n_i ) / sqrt( t_j t_j )
- !
- ! where the shear traction t_i is the projection of the traction
- ! vector on the plane surface. the strain amplitude is given by
- !
- ! gamma^dot = H( t_j r_j ) 2 vo sinh( taus / (t_c )
- !
- ! where taus is the effective shear on the fault plane,
- !
- ! taus = tau + mu*sigma
- !
- ! where tau is the shear and sigma the normal stress. tau and sigma
- ! assumed to be the co-seismic change only, not the absolute
- ! stress. vo is a reference slip velocity, and t_c, the critical
- ! stress, corresponds to (a-b)*sigma in the framework of rate-and-
- ! state friction. the effective viscosity eta* and the fluidity
- !
- ! eta* = tau / gamma^dot
- ! fluidity = 1 / eta*
- !
- ! are used to compute the optimal time-step. H( x ) is the
- ! Heaviside function and r_i is the rake vector. I impose
- ! gamma^dot to be zero is t_j r_j < 0. This constraint is
- ! enforced to ensure that no back slip occurs on faults.
- !
- ! sylvain barbot (07/24/07) - original form
- ! (02/28/11) - add constraints on the direction of
- ! afterslip
+ !> subroutine FrictionEigenStress
+ !! compute the eigen-stress (forcing moment) to be relaxed by
+ !! rate-dependent inelastic deformation in the case of a frictional
+ !! surface:
+ !!
+ !! sigma^i = C:F:sigma
+ !!
+ !! where C is the elastic moduli tensor, F is the heterogeneous
+ !! fluidity moduli tensor and sigma is the instantaneous stress
+ !! tensor. for a frictional surface, the eigenstrain-rate is given
+ !! by
+ !!
+ !! epsilon^i^dot = F:sigma = gamma^dot R
+ !!
+ !! where gamma^dot is the slip rate (a scalar) and R is the
+ !! deviatoric, symmetric, and unitary, tensor:
+ !!
+ !! R_ij = 1/2 ( t_i n_j + t_j n_i ) / sqrt( t_j t_j )
+ !!
+ !! where the shear traction t_i is the projection of the traction
+ !! vector on the plane surface. the strain amplitude is given by
+ !!
+ !! gamma^dot = H( t_j r_j ) 2 vo sinh( taus / (t_c )
+ !!
+ !! where taus is the effective shear on the fault plane,
+ !!
+ !! taus = tau + mu*sigma
+ !!
+ !! where tau is the shear and sigma the normal stress. tau and sigma
+ !! assumed to be the co-seismic change only, not the absolute
+ !! stress. vo is a reference slip velocity, and t_c, the critical
+ !! stress, corresponds to (a-b)*sigma in the framework of rate-and-
+ !! state friction. the effective viscosity eta* and the fluidity
+ !!
+ !! eta* = tau / gamma^dot
+ !! fluidity = 1 / eta*
+ !!
+ !! are used to compute the optimal time-step. H( x ) is the
+ !! Heaviside function and r_i is the rake vector. I impose
+ !! gamma^dot to be zero is t_j r_j < 0. This constraint is
+ !! enforced to ensure that no back slip occurs on faults.
+ !!
+ !! \author sylvain barbot (07/24/07) - original form
+ !! (02/28/11) - add constraints on the direction
+ !! of afterslip
!-----------------------------------------------------------------
SUBROUTINE frictioneigenstress(x,y,z,L,W,strike,dip,rake,beta, &
sig,mu,structure,sx1,sx2,sx3,dx1,dx2,dx3,moment,maxwelltime,vel)
@@ -381,24 +381,24 @@ CONTAINS
END SUBROUTINE frictioneigenstress
!---------------------------------------------------------------------
- ! function MonitorFriction
- ! samples a scalar field along a specified planar surface.
- !
- ! input:
- ! x,y,z coordinates of the creeping segment
- ! L dimension of segment in the depth direction
- ! W dimension of segment in the strike direction
- ! beta smoothing factor
- ! sx1,2,3 dimension of the stress tensor array
- ! dx1,2,3 sampling size
- ! sig stress tensor array
- ! structure frictional properties as a function of depth
- !
- ! output:
- ! patch list of strike- and dip-slip as a function of position
- ! on the fault.
- !
- ! sylvain barbot (10-16-07) - original form
+ !> function MonitorFriction
+ !! samples a scalar field along a specified planar surface.
+ !!
+ !! input:
+ !! @param x,y,z coordinates of the creeping segment
+ !! @param L dimension of segment in the depth direction
+ !! @param W dimension of segment in the strike direction
+ !! @param beta smoothing factor
+ !! @param sx1,2,3 dimension of the stress tensor array
+ !! @param dx1,2,3 sampling size
+ !! @param sig stress tensor array
+ !! @param structure frictional properties as a function of depth
+ !!
+ !! output:
+ !! @param patch list of strike- and dip-slip as a function of position
+ !! on the fault.
+ !!
+ !! \author sylvain barbot (10-16-07) - original form
!---------------------------------------------------------------------
SUBROUTINE monitorfriction(x,y,z,L,W,strike,dip,rake,beta, &
sx1,sx2,sx3,dx1,dx2,dx3,sig,structure,patch)
diff -r 777cc20187d0 -r 0a07ece6fb5a green.f90
--- a/green.f90 Thu May 12 17:04:47 2011 -0700
+++ b/green.f90 Thu May 12 19:01:34 2011 -0700
@@ -35,18 +35,18 @@ CONTAINS
CONTAINS
!------------------------------------------------------------------------
- ! Subroutine ElasticResponse
- ! apply the 2d elastic (half-space) transfert function
- ! to the set of body forces.
- !
- ! INPUT:
- ! mu shear modulus
- ! f2 equivalent body-forces in the Fourier domain
- ! sx1, sx3
- !
- ! sylvain barbot (04/14/07) - original form
- ! (02/06/09) - parallel implementation with MPI and OpenMP
- ! (01/06/11) - remove implementation with MPI
+ !> Subroutine ElasticResponse
+ !! apply the 2d elastic (half-space) transfert function
+ !! to the set of body forces.
+ !!
+ !! INPUT:
+ !! @param mu shear modulus
+ !! @param f1,2,3 equivalent body-forces in the Fourier domain
+ !! @param dx1,2,3 sampling size
+ !!
+ !! \author sylvain barbot (04/14/07) - original form
+ !! (02/06/09) - parallel implementation with MPI and OpenMP
+ !! (01/06/11) - remove implementation with MPI
!------------------------------------------------------------------------
SUBROUTINE elasticresponse(lambda,mu,f1,f2,f3,dx1,dx2,dx3)
REAL*8, INTENT(IN) :: lambda,mu,dx1,dx2,dx3
@@ -101,13 +101,13 @@ CONTAINS
END SUBROUTINE elasticresponse
!---------------------------------------------------------------------
- ! subroutine SurfaceNormalTraction
- ! computes the two-dimensional field of surface normal stress
- ! expressed in the Fourier domain.
- ! The surface (x3=0) solution is obtained by integrating over the
- ! wavenumbers in 3-direction in the Fourier domain.
- !
- ! sylvain barbot (05-01-07) - original form
+ !> subroutine SurfaceNormalTraction
+ !! computes the two-dimensional field of surface normal stress
+ !! expressed in the Fourier domain.
+ !! The surface (x3=0) solution is obtained by integrating over the
+ !! wavenumbers in 3-direction in the Fourier domain.
+ !!
+ !! \author sylvain barbot (05-01-07) - original form
!---------------------------------------------------------------------
SUBROUTINE surfacenormaltraction(lambda, mu, u1, u2, u3, dx1, dx2, dx3, p)
REAL*4, INTENT(IN), DIMENSION(:,:,:) :: u1, u2, u3
@@ -147,10 +147,10 @@ CONTAINS
END SUBROUTINE surfacenormaltraction
!---------------------------------------------------------------------
- ! subroutine Boussinesq3D
- ! computes the deformation field in the 3-dimensional grid
- ! due to a normal stress at the surface. Apply the Fourier domain
- ! solution of Steketee [1958].
+ !> subroutine Boussinesq3D
+ !! computes the deformation field in the 3-dimensional grid
+ !! due to a normal stress at the surface. Apply the Fourier domain
+ !! solution of Steketee [1958].
!---------------------------------------------------------------------
SUBROUTINE boussinesq3d(p,lambda,mu,u1,u2,u3,dx1,dx2,dx3)
REAL*4, DIMENSION(:,:), INTENT(IN) :: p
@@ -206,14 +206,14 @@ CONTAINS
CONTAINS
!-----------------------------------------------------------------
- ! subroutine SteketeeSolution
- ! computes the spectrum (two-dimensional Fourier transform)
- ! of the 3 components of the deformation field u1, u2, u3
- ! at wavenumbers k1, k2 and position x3. This is the analytical
- ! solution of [J. A. Steketee, On Volterra's dislocations in a
- ! semi-infinite elastic medium, Canadian Journal of Physics, 1958]
- !
- ! sylvain barbot (05-02-07) - original form
+ !> subroutine SteketeeSolution
+ !! computes the spectrum (two-dimensional Fourier transform)
+ !! of the 3 components of the deformation field u1, u2, u3
+ !! at wavenumbers k1, k2 and position x3. This is the analytical
+ !! solution of [J. A. Steketee, On Volterra's dislocations in a
+ !! semi-infinite elastic medium, Canadian Journal of Physics, 1958]
+ !!
+ !! \author sylvain barbot (05-02-07) - original form
!-----------------------------------------------------------------
SUBROUTINE steketeesolution(p,alpha,u1,u2,u3,k1,k2,x3)
COMPLEX, INTENT(INOUT) :: u1, u2, u3
@@ -245,14 +245,14 @@ CONTAINS
END SUBROUTINE boussinesq3d
!---------------------------------------------------------------------
- ! subroutine SurfaceTraction
- ! computes the two-dimensional field of surface normal stress
- ! expressed in the Fourier domain.
- ! The surface (x3=0) solution is obtained by integrating over the
- ! wavenumbers in 3-direction in the Fourier domain.
- !
- ! sylvain barbot (07-07-07) - original form
- ! (02-09-09) - parallelized with mpi and openmp
+ !> subroutine SurfaceTraction
+ !! computes the two-dimensional field of surface normal stress
+ !! expressed in the Fourier domain.
+ !! The surface (x3=0) solution is obtained by integrating over the
+ !! wavenumbers in 3-direction in the Fourier domain.
+ !!
+ !! \author sylvain barbot (07-07-07) - original form
+ ! (02-09-09) - parallelized with mpi and openmp
!---------------------------------------------------------------------
SUBROUTINE surfacetraction(lambda,mu,u1,u2,u3,dx1,dx2,dx3,p1,p2,p3)
REAL*4, INTENT(IN), DIMENSION(:,:,:) :: u1,u2,u3
@@ -308,21 +308,21 @@ CONTAINS
END SUBROUTINE surfacetraction
!---------------------------------------------------------------------
- ! subroutine SurfaceTractionCowling
- ! computes the two-dimensional field of the resulting traction
- ! expressed in the Fourier domain in the presence of gravity.
- !
- ! The surface solution (x3=0) is obtained from the Fourier domain
- ! array by integrating over the wavenumbers in 3-direction.
- !
- ! The effective traction at x3=0 is
- !
- ! t_1 = sigma_13
- ! t_2 = sigma_23
- ! t_3 = sigma_33 - r g u3
- ! = sigma_33 - 2 mu alpha gamma u3
- !
- ! sylvain barbot (07-07-07) - original form
+ !> subroutine SurfaceTractionCowling
+ !! computes the two-dimensional field of the resulting traction
+ !! expressed in the Fourier domain in the presence of gravity.
+ !!
+ !! The surface solution (x3=0) is obtained from the Fourier domain
+ !! array by integrating over the wavenumbers in 3-direction.
+ !!
+ !! The effective traction at x3=0 is
+ !!
+ !! t_1 = sigma_13
+ !! t_2 = sigma_23
+ !! t_3 = sigma_33 - r g u3
+ !! = sigma_33 - 2 mu alpha gamma u3
+ !!
+ !! \author sylvain barbot (07-07-07) - original form
!---------------------------------------------------------------------
SUBROUTINE surfacetractioncowling(lambda,mu,gamma,u1,u2,u3,dx1,dx2,dx3, &
p1,p2,p3)
@@ -377,11 +377,11 @@ CONTAINS
END SUBROUTINE surfacetractioncowling
!---------------------------------------------------------------------
- ! subroutine Cerruti3D
- ! computes the deformation field in the 3-dimensional grid
- ! due to an arbitrary surface traction.
- !
- ! sylvain barbot (07/07/07) - original form
+ !> subroutine Cerruti3D
+ !! computes the deformation field in the 3-dimensional grid
+ !! due to an arbitrary surface traction.
+ !!
+ !! \author sylvain barbot (07/07/07) - original form
! (02/01/09) - parallelized with MPI and OpenMP
! (01/06/11) - remove parallelized version with MPI
!---------------------------------------------------------------------
@@ -463,15 +463,15 @@ CONTAINS
CONTAINS
!-----------------------------------------------------------------
- ! subroutine CerrutiSolution
- ! computes the general solution for the deformation field in an
- ! elastic half-space due to an arbitrary surface traction.
- ! the 3 components u1, u2, u3 of the deformation field are
- ! expressed in the horizontal Fourier at depth x3.
- ! this combines the solution to the Boussinesq's and the Cerruti's
- ! problem in a half-space.
- !
- ! sylvain barbot (07-07-07) - original form
+ !> subroutine CerrutiSolution
+ !! computes the general solution for the deformation field in an
+ !! elastic half-space due to an arbitrary surface traction.
+ !! the 3 components u1, u2, u3 of the deformation field are
+ !! expressed in the horizontal Fourier at depth x3.
+ !! this combines the solution to the Boussinesq's and the Cerruti's
+ !! problem in a half-space.
+ !!
+ !! \author sylvain barbot (07-07-07) - original form
!-----------------------------------------------------------------
SUBROUTINE cerrutisolution(mu,p1,p2,p3,alpha,u1,u2,u3,k1,k2,x3)
COMPLEX(KIND=4), INTENT(INOUT) :: u1,u2,u3
@@ -518,14 +518,14 @@ CONTAINS
END SUBROUTINE cerruti3d
!---------------------------------------------------------------------
- ! subroutine CerrutiCowling
- ! computes the deformation field in the 3-dimensional grid
- ! due to an arbitrary surface traction.
- !
- ! sylvain barbot - (07/07/07) - original form
- ! (21/11/08) - gravity effect
- ! (02/01/09) - parallelized with MPI and OpenMP
- ! (01/06/11) - remove parallelized version with MPI
+ !> subroutine CerrutiCowling
+ !! computes the deformation field in the 3-dimensional grid
+ !! due to an arbitrary surface traction.
+ !!
+ !! \author sylvain barbot - (07/07/07) - original form
+ !! (21/11/08) - gravity effect
+ !! (02/01/09) - parallelized with MPI and OpenMP
+ !! (01/06/11) - remove parallelized version with MPI
!---------------------------------------------------------------------
SUBROUTINE cerruticowling(p1,p2,p3,lambda,mu,gamma,u1,u2,u3,dx1,dx2,dx3)
REAL*4, DIMENSION(:,:), INTENT(IN) :: p1,p2,p3
@@ -607,16 +607,16 @@ CONTAINS
CONTAINS
!-----------------------------------------------------------------
- ! subroutine CerrutiSolCowling
- ! computes the general solution for the deformation field in an
- ! elastic half-space due to an arbitrary surface traction in the
- ! presence of gravity.
- !
- ! The 3 components u1, u2 and u3 of the deformation field are
- ! expressed in the horizontal Fourier at depth x3.
- !
- ! Combines the solution to the Boussinesq's and the Cerruti's
- ! problem in a half-space with buoyancy boundary conditions.
+ !> subroutine CerrutiSolCowling
+ !! computes the general solution for the deformation field in an
+ !! elastic half-space due to an arbitrary surface traction in the
+ !! presence of gravity.
+ !!
+ !! The 3 components u1, u2 and u3 of the deformation field are
+ !! expressed in the horizontal Fourier at depth x3.
+ !!
+ !! Combines the solution to the Boussinesq's and the Cerruti's
+ !! problem in a half-space with buoyancy boundary conditions.
!
! sylvain barbot (07-07-07) - original form
! (08-30-10) - account for net surface traction
@@ -672,9 +672,9 @@ CONTAINS
END SUBROUTINE cerruticowling
!---------------------------------------------------------------------
- ! subroutine CerrutiCowlingSerial
- ! computes the deformation field in the 3-dimensional grid
- ! due to an arbitrary surface traction. No parallel version.
+ !> subroutine CerrutiCowlingSerial
+ !! computes the deformation field in the 3-dimensional grid
+ !! due to an arbitrary surface traction. No parallel version.
!
! sylvain barbot - 07/07/07 - original form
! 21/11/08 - gravity effect
@@ -734,16 +734,16 @@ CONTAINS
CONTAINS
!-----------------------------------------------------------------
- ! subroutine CerrutiSolCowling
- ! computes the general solution for the deformation field in an
- ! elastic half-space due to an arbitrary surface traction in the
- ! presence of gravity.
- !
- ! The 3 components u1, u2 and u3 of the deformation field are
- ! expressed in the horizontal Fourier at depth x3.
- !
- ! Combines the solution to the Boussinesq's and the Cerruti's
- ! problem in a half-space with buoyancy boundary conditions.
+ !> subroutine CerrutiSolCowling
+ !! computes the general solution for the deformation field in an
+ !! elastic half-space due to an arbitrary surface traction in the
+ !! presence of gravity.
+ !!
+ !! The 3 components u1, u2 and u3 of the deformation field are
+ !! expressed in the horizontal Fourier at depth x3.
+ !!
+ !! Combines the solution to the Boussinesq's and the Cerruti's
+ !! problem in a half-space with buoyancy boundary conditions.
!
! sylvain barbot (07-07-07) - original form
!-----------------------------------------------------------------
@@ -794,17 +794,17 @@ CONTAINS
END SUBROUTINE cerruticowlingserial
!------------------------------------------------------------------
- ! subroutine GreenFunction
- ! computes (inplace) the displacement components due to a set of
- ! 3-D body-forces by application of the semi-analytic Green's
- ! function. The solution satisfies quasi-static Navier's equation
- ! including vanishing of normal traction at the surface.
- !
- ! The surface traction can be made to vanish by application of
- ! 1) method of images + boussinesq problem (grn_method=GRN_IMAGE)
- ! 2) boussinesq's and cerruti's problems (grn_method=GRN_HS)
- ! in the first case, the body-forces are supposed by have been
- ! imaged appropriately.
+ !> subroutine GreenFunction
+ !! computes (inplace) the displacement components due to a set of
+ !! 3-D body-forces by application of the semi-analytic Green's
+ !! function. The solution satisfies quasi-static Navier's equation
+ !! including vanishing of normal traction at the surface.
+ !!
+ !! The surface traction can be made to vanish by application of
+ !! 1) method of images + boussinesq problem (grn_method=GRN_IMAGE)
+ !! 2) boussinesq's and cerruti's problems (grn_method=GRN_HS)
+ !! in the first case, the body-forces are supposed by have been
+ !! imaged appropriately.
!
! sylvain barbot (07/07/07) - original form
!------------------------------------------------------------------
@@ -867,30 +867,30 @@ CONTAINS
END SUBROUTINE greenfunction
!------------------------------------------------------------------
- ! subroutine GreensFunctionCowling
- ! computes (inplace) the displacement components due to a set of
- ! 3-D body-forces by application of the semi-analytic Green's
- ! function. The solution satisfies quasi-static Navier's equation
- ! with buoyancy boundary condition to simulate the effect of
- ! gravity (the Cowling approximation).
- !
- ! the importance of gravity depends upon the density contrast rho
- ! at the surface, the acceleration of gravity g and the value of
- ! shear modulus mu in the crust. effect on the displacement field
- ! is governed by the gradient
- !
- ! gamma = (1 - nu) rho g / mu
- ! = rho g / (2 mu alpha)
- !
- ! where nu is the Poisson's ratio. For a Poisson's solid with
- ! nu = 1/4, with a density contrast rho = 3200 kg/m^3 and a shear
- ! modulus mu = 30 GPa, we have gamma = 0.8e-6 /m.
- !
- ! INPUT:
- ! . c1,c2,c3 is a set of body forces
- ! . dx1,dx2,dx3 are the sampling size
- ! . lambda,mu are the Lame elastic parameters
- ! . gamma is the gravity coefficient
+ !> subroutine GreensFunctionCowling
+ !! computes (inplace) the displacement components due to a set of
+ !! 3-D body-forces by application of the semi-analytic Green's
+ !! function. The solution satisfies quasi-static Navier's equation
+ !! with buoyancy boundary condition to simulate the effect of
+ !! gravity (the Cowling approximation).
+ !!
+ !! the importance of gravity depends upon the density contrast rho
+ !! at the surface, the acceleration of gravity g and the value of
+ !! shear modulus mu in the crust. effect on the displacement field
+ !! is governed by the gradient
+ !!
+ !! gamma = (1 - nu) rho g / mu
+ !! = rho g / (2 mu alpha)
+ !!
+ !! where nu is the Poisson's ratio. For a Poisson's solid with
+ !! nu = 1/4, with a density contrast rho = 3200 kg/m^3 and a shear
+ !! modulus mu = 30 GPa, we have gamma = 0.8e-6 /m.
+ !!
+ !! INPUT:
+ !! @param c1,c2,c3 is a set of body forces
+ !! @param dx1,dx2,dx3 are the sampling size
+ !! @param lambda,mu are the Lame elastic parameters
+ !! @param gamma is the gravity coefficient
!
! sylvain barbot (07/07/07) - original function greenfunction
! (11/21/08) - effect of gravity
diff -r 777cc20187d0 -r 0a07ece6fb5a relax.f90
--- a/relax.f90 Thu May 12 17:04:47 2011 -0700
+++ b/relax.f90 Thu May 12 19:01:34 2011 -0700
@@ -17,147 +17,171 @@
! along with RELAX. If not, see <http://www.gnu.org/licenses/>.
!-----------------------------------------------------------------------
+ !-----------------------------------------------------------------------
+ !> \mainpage
+ !! program relax
+ !! <hr>
+ !! PURPOSE:
+ !! The program RELAX computes nonlinear time-dependent viscoelastic
+ !! deformation with powerlaw rheology and rate-strengthening friction
+ !! in a cubic, periodic grid due to coseismic stress changes, initial
+ !! stress, surface loads, and/or moving faults.
+ !!
+ !! ONLINE DOCUMENTATION:
+ !! generate html documentation from the source directory with the
+ !! doxygen (http://www.stack.nl/~dimitri/doxygen/index.html)
+ !! program with command:
+ !!
+ !! doxygen .doxygen
+ !!
+ !! DESCRIPTION:
+ !! Computation is done semi-analytically inside a cartesian grid.
+ !! The grid is defined by its size sx1*sx2*sx3 and the sampling
+ !! intervals dx1, dx2 and dx3. rule of thumb is to allow for at least
+ !! five samples per fault length or width, and to have the tip of any
+ !! fault at least 10 fault widths away from any edge of the
+ !! computational grid.
+ !!
+ !! Coseismic stress changes and initial coseismic deformation results
+ !! from the presence of dislocations in the brittle layer. Fault
+ !! geometry is prescribed following Okada or Wang's convention, with the
+ !! usual slip, strike, dip and rake and is converted to a double-couple
+ !! equivalent body-force analytically. Current implementation allows
+ !! shear fault (strike slip and dip slip), dykes, Mogi source, and
+ !! surface traction. Faults and dykes can be of arbitrary orientation
+ !! in the half space.
+ !!
+ !! <hr>
+ !!
+ !! METHOD:
+ !! The current implementation is organized to integrate stress/strain-
+ !! rate constitutive laws (rheologies) of the form
+ !! \f[
+ !! \dot{\epsilon} = f(\sigma)
+ !! \f]
+ !! as opposed to epsilon^dot = f(sigma,epsilon) wich would include work-
+ !! hardening (or weakening). The time-stepping implements a second-order
+ !! Runge-Kutta numerical integration scheme with a variable time-step.
+ !! The Runge-Kutta method integrating the ODE y'=f(x,y) can be summarized
+ !! as follows:
+ !! \f[
+ !! y_(n+1) = y_n + k_2
+ !! k_1 = h * f(x_n, y_n)
+ !! k_2 = h * f(x_n + h, y_n + k_1)
+ !! \f]
+ !! where h is the time-step and n is the time-index. The elastic response
+ !! in the computational grid is obtained using elastic Greens functions.
+ !! The Greens functions are applied in the Fourier domain. Strain,
+ !! stress and body-forces are obtained by application of a finite impulse
+ !! response (FIR) differentiator filter in the space domain.
+ !!
+ !! <hr>
+ !!
+ !! INPUT:
+ !! Static dislocation sources are discretized into a series of planar
+ !! segments. Slip patches are defined in terms of position, orientation,
+ !! and slip, as illustrated in the following figure:
+ !!\verbatim
+ !! N (x1)
+ !! /
+ !! /| Strike
+ !! x1,x2,x3 ->@------------------------ (x2)
+ !! |\ p . \ W
+ !! :-\ i . \ i
+ !! | \ l . \ d
+ !! :90 \ S . \ t
+ !! |-Dip\ . \ h
+ !! : \. | Rake \
+ !! | -------------------------
+ !! : L e n g t h
+ !! Z (x3)
+ !!\endverbatim
+ !! Dislocations are converted to double-couple equivalent body-force
+ !! analytically. Solution displacement is obtained by application of
+ !! the Greens functions in the Fourier domain.
+ !!
+ !! For friction faults where slip rates are evaluated from stress and
+ !! a constitutive law, the rake corresponds to the orientation of slip.
+ !! That is, if r_i is the rake vector and v_i is the instantaneous
+ !! velocity vector, then r_j v_j >= 0.
+ !!
+ !! <hr>
+ !!
+ !! OUTPUT:
+ !! The vector-valued deformation is computed everywhere in a cartesian
+ !! grid. The vector field is sampled 1) along a horizontal surface at a
+ !! specified depth and 2) at specific points. Format is always North (x1),
+ !! East (x2) and Down (x3) components, following the right-handed reference
+ !! system convention. North corresponds to x1-direction, East to the
+ !! x2-direction and down to the x3-direction. The Generic Mapping Tool
+ !! output files are labeled explicitely ???-north.grd, ???-east.grd and
+ !! ???-up.grd (or say, ???-geo-up.grd for outputs in geographic
+ !! coordinates), where ??? stands for an output index: 001, 002, ...
+ !!
+ !! The amplitude of the inelastic (irreversible) deformation is also
+ !! tracked and can be output along a plane of arbitrary orientation.
+ !! The inelastic deformation includes the initial, constrained, slip on
+ !! fault surfaces, the time-dependent slip on frictional surfaces and
+ !! the cumulative amplitude of bulk strain in viscoelastic regions.
+ !! Slip is provided as a function of local coordinates along strike and
+ !! dip as well as a function of the Cartesian coordinates for three-
+ !! dimensional display.
+ !!
+ !! Time integration uses adaptive time steps to ensure accuracy but
+ !! results can be output either 1) at specified uniform time intervals
+ !! or 2) at the same intervals as computed. In the later case, output
+ !! intervals is chosen internally depending on instantaneous relaxation
+ !! rates.
+ !!
+ !! <hr>
+ !!
+ !! TECHNICAL ASPECTS:
+ !! Most of the computational burden comes from 1) applying the elastic
+ !! Green function and 2) computing the current strain from a displacement
+ !! field. The convolution of body forces with the Green function is
+ !! performed in the Fourier domain and the efficiency of the computation
+ !! depends essentially upon a choice of the discrete Fourier transform.
+ !! Current implementation is compatible with the Couley-Tuckey, the
+ !! Fast Fourier transform of the West (FFTW), the SGI FFT and the intel
+ !! FFT from the intel MKL library. Among these choices, the MKL FFT is
+ !! the most efficient. The FFTW, SGI FFT and MKL FFT can all be ran
+ !! in parallel on shared-memory computers.
+ !!
+ !! Strain is computed using a Finite Impulse Response differentiator
+ !! filter in the space domain. Use of FIR filter give rise to very
+ !! accurate derivatives but is computationally expensive. The filter
+ !! kernels are provided in the kernel???.inc files. Use of a compact
+ !! kernel may accelerate computation significantly.
+ !!
+ !! Compilation options are defined in the include.f90 file and specify
+ !! for instance the choice of DFT and the kind of output provided.
+ !!
+ !! MODIFICATIONS:
+ !! \author Sylvain Barbot
+ !! (07-06-07) - original form
+ !! (08-28-08) - FFTW/SGI_FFT support, FIR derivatives,
+ !! Runge-Kutta integration, tensile cracks,
+ !! GMT output, comments in input file
+ !! (10-24-08) - interseismic loading, postseismic signal
+ !! output in separate files
+ !! (12-08-09) - slip distribution smoothing
+ !! (05-05-10) - lateral variations in viscous properties
+ !! Intel MKL implementation of the FFT
+ !! (06-04-10) - output in geographic coordinates
+ !! and output components of stress tensor
+ !! (07-19-10) - includes surface tractions initial condition
+ !! output geometry in VTK format for Paraview
+ !! (02-28-11) - add constraints on the broad direction of
+ !! afterslip, export faults to GMT xy format
+ !! and allow scaling of computed time steps.
+ !! (04-26-11) - include command-line arguments
+ !!
+ !! \todo
+ !! - export VTK in binary instead of ascii format
+ !! - homogenize VTK output so that geometry of events match event index
+ !! - evaluate Green function, stress and body forces in GPU
+ !------------------------------------------------------------------------
PROGRAM relax
- !-----------------------------------------------------------------------
- ! PURPOSE:
- ! The program RELAX computes nonlinear time-dependent viscoelastic
- ! deformation with powerlaw rheology and rate-strengthening friction
- ! in a cubic, periodic grid due to coseismic stress changes, initial
- ! stress, surface loads, and/or moving faults.
- !
- ! DESCRIPTION:
- ! Computation is done semi-analytically inside a cartesian grid.
- ! The grid is defined by its size sx1*sx2*sx3 and the sampling
- ! intervals dx1, dx2 and dx3. rule of thumb is to allow for at least
- ! five samples per fault length or width, and to have the tip of any
- ! fault at least 10 fault widths away from any edge of the
- ! computational grid.
- !
- ! Coseismic stress changes and initial coseismic deformation results
- ! from the presence of dislocations in the brittle layer. Fault
- ! geometry is prescribed following Okada or Wang's convention, with the
- ! usual slip, strike, dip and rake and is converted to a double-couple
- ! equivalent body-force analytically. Current implementation allows
- ! shear fault (strike slip and dip slip), dykes, Mogi source, and
- ! surface traction. Faults and dykes can be of arbitrary orientation
- ! in the half space.
- !
- ! METHOD:
- ! The current implementation is organized to integrate stress/strain-
- ! rate constitutive laws (rheologies) of the form
- !
- ! epsilon^dot = f(sigma)
- !
- ! as opposed to epsilon^dot = f(sigma,epsilon) wich would include work-
- ! hardening (or weakening). The time-stepping implements a second-order
- ! Runge-Kutta numerical integration scheme with a variable time-step.
- ! The Runge-Kutta method integrating the ODE y'=f(x,y) can be summarized
- ! as follows:
- !
- ! y_(n+1) = y_n + k_2
- ! k_1 = h * f(x_n, y_n)
- ! k_2 = h * f(x_n + h, y_n + k_1)
- !
- ! where h is the time-step and n is the time-index. The elastic response
- ! in the computational grid is obtained using elastic Greens functions.
- ! The Greens functions are applied in the Fourier domain. Strain,
- ! stress and body-forces are obtained by application of a finite impulse
- ! response (FIR) differentiator filter in the space domain.
- !
- ! INPUT:
- ! Static dislocation sources are discretized into a series of planar
- ! segments. Slip patches are defined in terms of position, orientation,
- ! and slip, as illustrated in the following figure:
- !
- ! N (x1)
- ! /
- ! /| Strike
- ! x1,x2,x3 ->@------------------------ (x2)
- ! |\ p . \ W
- ! :-\ i . \ i
- ! | \ l . \ d
- ! :90 \ S . \ t
- ! |-Dip\ . \ h
- ! : \. | Rake \
- ! | -------------------------
- ! : L e n g t h
- ! Z (x3)
- !
- ! Dislocations are converted to double-couple equivalent body-force
- ! analytically. Solution displacement is obtained by application of
- ! the Greens functions in the Fourier domain.
- !
- ! For friction faults where slip rates are evaluated from stress and
- ! a constitutive law, the rake corresponds to the orientation of slip.
- ! That is, if r_i is the rake vector and v_i is the instantaneous
- ! velocity vector, then r_j v_j >= 0.
- !
- ! OUTPUT:
- ! The vector-valued deformation is computed everywhere in a cartesian
- ! grid. The vector field is sampled 1) along a horizontal surface at a
- ! specified depth and 2) at specific points. Format is always North (x1),
- ! East (x2) and Down (x3) components, following the right-handed reference
- ! system convention. North corresponds to x1-direction, East to the
- ! x2-direction and down to the x3-direction. The Generic Mapping Tool
- ! output files are labeled explicitely ???-north.grd, ???-east.grd and
- ! ???-up.grd (or say, ???-geo-up.grd for outputs in geographic
- ! coordinates), where ??? stands for an output index: 001, 002, ...
- !
- ! The amplitude of the inelastic (irreversible) deformation is also
- ! tracked and can be output along a plane of arbitrary orientation.
- ! The inelastic deformation includes the initial, constrained, slip on
- ! fault surfaces, the time-dependent slip on frictional surfaces and
- ! the cumulative amplitude of bulk strain in viscoelastic regions.
- ! Slip is provided as a function of local coordinates along strike and
- ! dip as well as a function of the Cartesian coordinates for three-
- ! dimensional display.
- !
- ! Time integration uses adaptive time steps to ensure accuracy but
- ! results can be output either 1) at specified uniform time intervals
- ! or 2) at the same intervals as computed. In the later case, output
- ! intervals is chosen internally depending on instantaneous relaxation
- ! rates.
- !
- ! TECHNICAL ASPECTS:
- ! Most of the computational burden comes from 1) applying the elastic
- ! Green function and 2) computing the current strain from a displacement
- ! field. The convolution of body forces with the Green function is
- ! performed in the Fourier domain and the efficiency of the computation
- ! depends essentially upon a choice of the discrete Fourier transform.
- ! Current implementation is compatible with the Couley-Tuckey, the
- ! Fast Fourier transform of the West (FFTW), the SGI FFT and the intel
- ! FFT from the intel MKL library. Among these choices, the MKL FFT is
- ! the most efficient. The FFTW, SGI FFT and MKL FFT can all be ran
- ! in parallel on shared-memory computers.
- !
- ! Strain is computed using a Finite Impulse Response differentiator
- ! filter in the space domain. Use of FIR filter give rise to very
- ! accurate derivatives but is computationally expensive. The filter
- ! kernels are provided in the kernel???.inc files. Use of a compact
- ! kernel may accelerate computation significantly.
- !
- ! Compilation options are defined in the include.f90 file and specify
- ! for instance the choice of DFT and the kind of output provided.
- !
- ! MODIFICATIONS:
- ! sylvain barbot (07-06-07) - original form
- ! (08-28-08) - FFTW/SGI_FFT support, FIR derivatives,
- ! Runge-Kutta integration, tensile cracks,
- ! GMT output, comments in input file
- ! (10-24-08) - interseismic loading, postseismic signal
- ! output in separate files
- ! (12-08-09) - slip distribution smoothing
- ! (05-05-10) - lateral variations in viscous properties
- ! Intel MKL implementation of the FFT
- ! (06-04-10) - output in geographic coordinates
- ! and output components of stress tensor
- ! (07-19-10) - includes surface tractions initial condition
- ! output geometry in VTK format for Paraview
- ! (02-28-11) - add constraints on the broad direction of
- ! afterslip, export faults to GMT xy format
- ! and allow scaling of computed time steps.
- ! (04-26-11) - include command-line arguments
- !-----------------------------------------------------------------------
USE types
USE input
@@ -695,10 +719,10 @@ CONTAINS
CONTAINS
!---------------------------------------------------------------------
- ! subroutine Traction
- ! assigns the traction vector at the surface
- !
- ! sylvain barbot (07-19-07) - original form
+ !> subroutine Traction
+ !! assigns the traction vector at the surface
+ !!
+ !! \author sylvain barbot (07-19-07) - original form
!---------------------------------------------------------------------
SUBROUTINE traction(mu,e,sx1,sx2,dx1,dx2,t3)
TYPE(EVENT_STRUC), INTENT(IN) :: e
@@ -722,11 +746,11 @@ CONTAINS
END SUBROUTINE traction
!--------------------------------------------------------------------
- ! subroutine dislocations
- ! assigns equivalent body forces or moment density to simulate
- ! shear dislocations and fault opening. add the corresponding moment
- ! density in the cumulative relaxed moment so that fault slip does
- ! not reverse in the postseismic time.
+ !> subroutine dislocations
+ !! assigns equivalent body forces or moment density to simulate
+ !! shear dislocations and fault opening. add the corresponding moment
+ !! density in the cumulative relaxed moment so that fault slip does
+ !! not reverse in the postseismic time.
!--------------------------------------------------------------------
SUBROUTINE dislocations(event,lambda,mu,beta,sx1,sx2,sx3,dx1,dx2,dx3, &
v1,v2,v3,t1,t2,t3,tau,factor,eigenstress)
@@ -844,15 +868,15 @@ CONTAINS
END SUBROUTINE dislocations
!--------------------------------------------------------------------
- ! function IsOutput
- ! checks if output should be written based on user choices: if output
- ! time interval (odt) is positive, output is written only if time
- ! is an integer of odt. If odt is negative output is written at times
- ! corresponding to internally chosen time steps.
- !
- ! IsOutput is true only at discrete intervals (skip=0,odt>0),
- ! or at every "skip" computational steps (skip>0,odt<0),
- ! or anytime a coseismic event occurs
+ !> function IsOutput
+ !! checks if output should be written based on user choices: if output
+ !! time interval (odt) is positive, output is written only if time
+ !! is an integer of odt. If odt is negative output is written at times
+ !! corresponding to internally chosen time steps.
+ !!
+ !! @return IsOutput is true only at discrete intervals (skip=0,odt>0),
+ !! or at every "skip" computational steps (skip>0,odt<0),
+ !! or anytime a coseismic event occurs
!
! Sylvain Barbot (07/06/09) - original form
!--------------------------------------------------------------------
@@ -871,46 +895,46 @@ CONTAINS
END FUNCTION isoutput
!--------------------------------------------------------------------
- ! subroutine IntegrationStep
- ! find the time-integration forward step for the predictor-corrector
- ! scheme.
- !
- ! input file line
- !
- ! time interval, (positive dt step) or (negative skip and scaling)
- !
- ! can be filled by either 1)
- !
- ! T, dt
- !
- ! where T is the time interval of the simulation and dt is the
- ! output time step, or 2)
- !
- ! T, -n, t_s
- !
- ! where n indicates the number of computational steps before
- ! outputing results, t_s is a scaling applied to internally
- ! computed time step.
- !
- ! for case 1), an optimal time step is evaluated internally to
- ! ensure stability (t_m/10) of time integration. The actual
- ! time step Dt is chosen as
- !
- ! Dt = min( t_m/10, ((t%odt)+1)*odt-t )
- !
- ! where t is the current time in the simulation. regardless of
- ! time step Dt, results are output if t is a multiple of dt.
- !
- ! for case 2), the time step is chosen internally based on an
- ! estimate of the relaxation time (t_m/10). Results are output
- ! every n steps. The actual time step is chosen as
- !
- ! Dt = min( t_m/10*t_s, t(next event)-t )
- !
- ! where index is the number of computational steps after a coseismic
- ! event and t(next event) is the time of the next coseismic event.
- !
- ! sylvain barbot (01/01/08) - original form
+ !> subroutine IntegrationStep
+ !! find the time-integration forward step for the predictor-corrector
+ !! scheme.
+ !!
+ !! input file line
+ !!
+ !! time interval, (positive dt step) or (negative skip and scaling)
+ !!
+ !! can be filled by either 1)
+ !!
+ !! T, dt
+ !!
+ !! where T is the time interval of the simulation and dt is the
+ !! output time step, or 2)
+ !!
+ !! T, -n, t_s
+ !!
+ !! where n indicates the number of computational steps before
+ !! outputing results, t_s is a scaling applied to internally
+ !! computed time step.
+ !!
+ !! for case 1), an optimal time step is evaluated internally to
+ !! ensure stability (t_m/10) of time integration. The actual
+ !! time step Dt is chosen as
+ !!
+ !! Dt = min( t_m/10, ((t%odt)+1)*odt-t )
+ !!
+ !! where t is the current time in the simulation. regardless of
+ !! time step Dt, results are output if t is a multiple of dt.
+ !!
+ !! for case 2), the time step is chosen internally based on an
+ !! estimate of the relaxation time (t_m/10). Results are output
+ !! every n steps. The actual time step is chosen as
+ !!
+ !! Dt = min( t_m/10*t_s, t(next event)-t )
+ !!
+ !! where index is the number of computational steps after a coseismic
+ !! event and t(next event) is the time of the next coseismic event.
+ !!
+ !! \author sylvain barbot (01/01/08) - original form
!--------------------------------------------------------------------
SUBROUTINE integrationstep(tm,Dt,t,oi,odt,skip,tscale,events,e,ne)
REAL*8, INTENT(INOUT) :: tm,Dt,odt
diff -r 777cc20187d0 -r 0a07ece6fb5a viscoelastic3d.f90
--- a/viscoelastic3d.f90 Thu May 12 17:04:47 2011 -0700
+++ b/viscoelastic3d.f90 Thu May 12 19:01:34 2011 -0700
@@ -32,20 +32,20 @@ CONTAINS
CONTAINS
!-----------------------------------------------------------------
- ! subroutine ViscoElasticDeviatoricStress
- ! computes the instantaneous deviatoric stress tensor sigma_ij'
- !
- ! sigma_ij' = 2*mu*(-delta_ij epsilon_kk/3 + epsilon_ij) - tau_ij
- !
- ! such as
- !
- ! sigma_kk'= 0
- !
- ! where tau_ij is a second-order deviatoric symmetric tensor
- ! that integrates the history of the relaxed stress. strain is
- ! estimated using a centered finite difference derivative.
- !
- ! sylvain barbot (07/07/07) - original form
+ !> subroutine ViscoElasticDeviatoricStress
+ !! computes the instantaneous deviatoric stress tensor sigma_ij'
+ !!
+ !! sigma_ij' = 2*mu*(-delta_ij epsilon_kk/3 + epsilon_ij) - tau_ij
+ !!
+ !! such as
+ !!
+ !! sigma_kk'= 0
+ !!
+ !! where tau_ij is a second-order deviatoric symmetric tensor
+ !! that integrates the history of the relaxed stress. strain is
+ !! estimated using a centered finite difference derivative.
+ !!
+ !! \author sylvain barbot (07/07/07) - original form
!-----------------------------------------------------------------
SUBROUTINE viscoelasticdeviatoricstress(mu,u1,u2,u3,tau,&
dx1,dx2,dx3,sx1,sx2,sx3,sig)
@@ -108,17 +108,17 @@ CONTAINS
END SUBROUTINE viscoelasticdeviatoricstress
!-----------------------------------------------------------------
- ! subroutine ViscousEigenstress
- ! computes the moment density rate due to a layered viscoelastic
- ! structure with powerlaw creep
- !
- ! d Ei / dt = C:F:sigma'
- !
- ! where C is the elastic moduli tensor, F is the heterogeneous
- ! fluidity tensor and sigma' is the instantaneous deviatoric
- ! stress. F is stress dependent (powerlaw creep.)
- !
- ! sylvain barbot (08/30/08) - original form
+ !> subroutine ViscousEigenstress
+ !! computes the moment density rate due to a layered viscoelastic
+ !! structure with powerlaw creep
+ !!
+ !! d Ei / dt = C:F:sigma'
+ !!
+ !! where C is the elastic moduli tensor, F is the heterogeneous
+ !! fluidity tensor and sigma' is the instantaneous deviatoric
+ !! stress. F is stress dependent (powerlaw creep.)
+ !!
+ !! \author sylvain barbot (08/30/08) - original form
!-----------------------------------------------------------------
SUBROUTINE viscouseigenstress(mu,structure,ductilezones,sig,sx1,sx2,sx3, &
dx1,dx2,dx3,moment,beta,maxwelltime,gamma)
@@ -214,12 +214,12 @@ CONTAINS
CONTAINS
!---------------------------------------------------------
- ! function dgammadot0
- ! evaluates the change of fluidity at position x1,x2,x3
- ! due to the presence of weak ductile zones. the extent
- ! and magnitude of ductile zones is tapered (beta).
- !
- ! sylvain barbot (3/29/10) - original form
+ !> function dgammadot0
+ !! evaluates the change of fluidity at position x1,x2,x3
+ !! due to the presence of weak ductile zones. the extent
+ !! and magnitude of ductile zones is tapered (beta).
+ !!
+ !! \author sylvain barbot (3/29/10) - original form
!---------------------------------------------------------
REAL*8 FUNCTION dgammadot0(zones,x1,x2,x3,beta)
TYPE(WEAK_STRUCT), INTENT(IN), DIMENSION(:) :: zones
More information about the CIG-COMMITS
mailing list