Harald Klimach avatar Harald Klimach committed 3eda3d0

Updated with developments since october

* Add more functions to the wscript to be used in other projects
* Added q-values as boundary property
* Removed derived quantities subdirectories
* Comments and documentation improvements
* Read Lua configuration on a single process and distribute by MPI_Bcast
* Implement unit tests with waf
* Improved timeControl
* Better stencil handling and communication of element properties
* Bugfixes for multilevel support
* Removed dep-procs array which might have gotten very large
* Initialize random-number generator in init_env
* New randomized spatial function
* Updated Aotus to support quadruple precision
* Bugfix for typed communication patterns
* Proper handling of hg revision string in case of failures

Comments (0)

Files changed (99)

-69e5a947a7561edbb42d5af475855a2f46e6ccbb aotus
+002ceb5996afddce4d44525bf6152544aacf3857 aotus

external/blas/dgemm.f

 *       SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
 * 
 *       .. Scalar Arguments ..
-*       DOUBLE PRECISION ALPHA,BETA
+*       REAL(kind=rk) ALPHA,BETA
 *       INTEGER K,LDA,LDB,LDC,M,N
 *       CHARACTER TRANSA,TRANSB
 *       ..
 *       .. Array Arguments ..
-*       DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
+*       REAL(kind=rk) A(LDA,*),B(LDB,*),C(LDC,*)
 *       ..
 *  
 *
 *>
 *> \param[in] ALPHA
 *> \verbatim
-*>          ALPHA is DOUBLE PRECISION.
+*>          ALPHA is REAL(kind=rk).
 *>           On entry, ALPHA specifies the scalar alpha.
 *> \endverbatim
 *>
 *> \param[in] A
 *> \verbatim
-*>          A is DOUBLE PRECISION array of DIMENSION ( LDA, ka ), where ka is
+*>          A is REAL(kind=rk) array of DIMENSION ( LDA, ka ), where ka is
 *>           k  when  TRANSA = 'N' or 'n',  and is  m  otherwise.
 *>           Before entry with  TRANSA = 'N' or 'n',  the leading  m by k
 *>           part of the array  A  must contain the matrix  A,  otherwise
 *>
 *> \param[in] B
 *> \verbatim
-*>          B is DOUBLE PRECISION array of DIMENSION ( LDB, kb ), where kb is
+*>          B is REAL(kind=rk) array of DIMENSION ( LDB, kb ), where kb is
 *>           n  when  TRANSB = 'N' or 'n',  and is  k  otherwise.
 *>           Before entry with  TRANSB = 'N' or 'n',  the leading  k by n
 *>           part of the array  B  must contain the matrix  B,  otherwise
 *>
 *> \param[in] BETA
 *> \verbatim
-*>          BETA is DOUBLE PRECISION.
+*>          BETA is REAL(kind=rk).
 *>           On entry,  BETA  specifies the scalar  beta.  When  BETA  is
 *>           supplied as zero then C need not be set on input.
 *> \endverbatim
 *>
 *> \param[in,out] C
 *> \verbatim
-*>          C is DOUBLE PRECISION array of DIMENSION ( LDC, n ).
+*>          C is REAL(kind=rk) array of DIMENSION ( LDC, n ).
 *>           Before entry, the leading  m by n  part of the array  C must
 *>           contain the matrix  C,  except when  beta  is zero, in which
 *>           case C need not be set on entry.
 *>
 *  =====================================================================
       SUBROUTINE DGEMM(TRANSA,TRANSB,M,N,K,ALPHA,A,LDA,B,LDB,BETA,C,LDC)
+      use env_module, only: rk
 *
 *  -- Reference BLAS level3 routine (version 3.4.0) --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
 *     November 2011
 *
 *     .. Scalar Arguments ..
-      DOUBLE PRECISION ALPHA,BETA
+      REAL(kind=rk) ALPHA,BETA
       INTEGER K,LDA,LDB,LDC,M,N
       CHARACTER TRANSA,TRANSB
 *     ..
 *     .. Array Arguments ..
-      DOUBLE PRECISION A(LDA,*),B(LDB,*),C(LDC,*)
+      REAL(kind=rk) A(LDA,*),B(LDB,*),C(LDC,*)
 *     ..
 *
 *  =====================================================================
       INTRINSIC MAX
 *     ..
 *     .. Local Scalars ..
-      DOUBLE PRECISION TEMP
+      REAL(kind=rk) TEMP
       INTEGER I,INFO,J,L,NCOLA,NROWA,NROWB
       LOGICAL NOTA,NOTB
 *     ..
 *     .. Parameters ..
-      DOUBLE PRECISION ONE,ZERO
+      REAL(kind=rk) ONE,ZERO
       PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
 *     ..
 *

external/blas/dgemv.f

 *       SUBROUTINE DGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
 * 
 *       .. Scalar Arguments ..
-*       DOUBLE PRECISION ALPHA,BETA
+*       REAL(kind=rk) ALPHA,BETA
 *       INTEGER INCX,INCY,LDA,M,N
 *       CHARACTER TRANS
 *       ..
 *       .. Array Arguments ..
-*       DOUBLE PRECISION A(LDA,*),X(*),Y(*)
+*       REAL(kind=rk) A(LDA,*),X(*),Y(*)
 *       ..
 *  
 *
 *>
 *> \param[in] ALPHA
 *> \verbatim
-*>          ALPHA is DOUBLE PRECISION.
+*>          ALPHA is REAL(kind=rk).
 *>           On entry, ALPHA specifies the scalar alpha.
 *> \endverbatim
 *>
 *> \param[in] A
 *> \verbatim
-*>          A is DOUBLE PRECISION array of DIMENSION ( LDA, n ).
+*>          A is REAL(kind=rk) array of DIMENSION ( LDA, n ).
 *>           Before entry, the leading m by n part of the array A must
 *>           contain the matrix of coefficients.
 *> \endverbatim
 *>
 *> \param[in] X
 *> \verbatim
-*>          X is DOUBLE PRECISION array of DIMENSION at least
+*>          X is REAL(kind=rk) array of DIMENSION at least
 *>           ( 1 + ( n - 1 )*abs( INCX ) ) when TRANS = 'N' or 'n'
 *>           and at least
 *>           ( 1 + ( m - 1 )*abs( INCX ) ) otherwise.
 *>
 *> \param[in] BETA
 *> \verbatim
-*>          BETA is DOUBLE PRECISION.
+*>          BETA is REAL(kind=rk).
 *>           On entry, BETA specifies the scalar beta. When BETA is
 *>           supplied as zero then Y need not be set on input.
 *> \endverbatim
 *>
 *> \param[in,out] Y
 *> \verbatim
-*>          Y is DOUBLE PRECISION array of DIMENSION at least
+*>          Y is REAL(kind=rk) array of DIMENSION at least
 *>           ( 1 + ( m - 1 )*abs( INCY ) ) when TRANS = 'N' or 'n'
 *>           and at least
 *>           ( 1 + ( n - 1 )*abs( INCY ) ) otherwise.
 *>
 *  =====================================================================
       SUBROUTINE DGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY)
+      use env_module, only: rk
 *
 *  -- Reference BLAS level2 routine (version 3.4.0) --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
 *     November 2011
 *
 *     .. Scalar Arguments ..
-      DOUBLE PRECISION ALPHA,BETA
+      REAL(kind=rk) ALPHA,BETA
       INTEGER INCX,INCY,LDA,M,N
       CHARACTER TRANS
 *     ..
 *     .. Array Arguments ..
-      DOUBLE PRECISION A(LDA,*),X(*),Y(*)
+      REAL(kind=rk) A(LDA,*),X(*),Y(*)
 *     ..
 *
 *  =====================================================================
 *
 *     .. Parameters ..
-      DOUBLE PRECISION ONE,ZERO
+      REAL(kind=rk) ONE,ZERO
       PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
 *     ..
 *     .. Local Scalars ..
-      DOUBLE PRECISION TEMP
+      REAL(kind=rk) TEMP
       INTEGER I,INFO,IX,IY,J,JX,JY,KX,KY,LENX,LENY
 *     ..
 *     .. External Functions ..

external/blas/dger.f

 *       SUBROUTINE DGER(M,N,ALPHA,X,INCX,Y,INCY,A,LDA)
 * 
 *       .. Scalar Arguments ..
-*       DOUBLE PRECISION ALPHA
+*       REAL(kind=rk) ALPHA
 *       INTEGER INCX,INCY,LDA,M,N
 *       ..
 *       .. Array Arguments ..
-*       DOUBLE PRECISION A(LDA,*),X(*),Y(*)
+*       REAL(kind=rk) A(LDA,*),X(*),Y(*)
 *       ..
 *  
 *
 *>
 *> \param[in] ALPHA
 *> \verbatim
-*>          ALPHA is DOUBLE PRECISION.
+*>          ALPHA is REAL(kind=rk).
 *>           On entry, ALPHA specifies the scalar alpha.
 *> \endverbatim
 *>
 *> \param[in] X
 *> \verbatim
-*>          X is DOUBLE PRECISION array of dimension at least
+*>          X is REAL(kind=rk) array of dimension at least
 *>           ( 1 + ( m - 1 )*abs( INCX ) ).
 *>           Before entry, the incremented array X must contain the m
 *>           element vector x.
 *>
 *> \param[in] Y
 *> \verbatim
-*>          Y is DOUBLE PRECISION array of dimension at least
+*>          Y is REAL(kind=rk) array of dimension at least
 *>           ( 1 + ( n - 1 )*abs( INCY ) ).
 *>           Before entry, the incremented array Y must contain the n
 *>           element vector y.
 *>
 *> \param[in,out] A
 *> \verbatim
-*>          A is DOUBLE PRECISION array of DIMENSION ( LDA, n ).
+*>          A is REAL(kind=rk) array of DIMENSION ( LDA, n ).
 *>           Before entry, the leading m by n part of the array A must
 *>           contain the matrix of coefficients. On exit, A is
 *>           overwritten by the updated matrix.
 *>
 *  =====================================================================
       SUBROUTINE DGER(M,N,ALPHA,X,INCX,Y,INCY,A,LDA)
+      use env_module, only: rk
 *
 *  -- Reference BLAS level2 routine (version 3.4.0) --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
 *     November 2011
 *
 *     .. Scalar Arguments ..
-      DOUBLE PRECISION ALPHA
+      REAL(kind=rk) ALPHA
       INTEGER INCX,INCY,LDA,M,N
 *     ..
 *     .. Array Arguments ..
-      DOUBLE PRECISION A(LDA,*),X(*),Y(*)
+      REAL(kind=rk) A(LDA,*),X(*),Y(*)
 *     ..
 *
 *  =====================================================================
 *
 *     .. Parameters ..
-      DOUBLE PRECISION ZERO
+      REAL(kind=rk) ZERO
       PARAMETER (ZERO=0.0D+0)
 *     ..
 *     .. Local Scalars ..
-      DOUBLE PRECISION TEMP
+      REAL(kind=rk) TEMP
       INTEGER I,INFO,IX,J,JY,KX
 *     ..
 *     .. External Subroutines ..

external/blas/dscal.f

 *       SUBROUTINE DSCAL(N,DA,DX,INCX)
 * 
 *       .. Scalar Arguments ..
-*       DOUBLE PRECISION DA
+*       REAL(kind=rk) DA
 *       INTEGER INCX,N
 *       ..
 *       .. Array Arguments ..
-*       DOUBLE PRECISION DX(*)
+*       REAL(kind=rk) DX(*)
 *       ..
 *  
 *
 *>
 *  =====================================================================
       SUBROUTINE DSCAL(N,DA,DX,INCX)
+      use env_module, only: rk
 *
 *  -- Reference BLAS level1 routine (version 3.4.0) --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
 *     November 2011
 *
 *     .. Scalar Arguments ..
-      DOUBLE PRECISION DA
+      REAL(kind=rk) DA
       INTEGER INCX,N
 *     ..
 *     .. Array Arguments ..
-      DOUBLE PRECISION DX(*)
+      REAL(kind=rk) DX(*)
 *     ..
 *
 *  =====================================================================

external/blas/dswap.f

 *       INTEGER INCX,INCY,N
 *       ..
 *       .. Array Arguments ..
-*       DOUBLE PRECISION DX(*),DY(*)
+*       REAL(kind=rk) DX(*),DY(*)
 *       ..
 *  
 *
 *>
 *  =====================================================================
       SUBROUTINE DSWAP(N,DX,INCX,DY,INCY)
+      use env_module, only: rk
 *
 *  -- Reference BLAS level1 routine (version 3.4.0) --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
       INTEGER INCX,INCY,N
 *     ..
 *     .. Array Arguments ..
-      DOUBLE PRECISION DX(*),DY(*)
+      REAL(kind=rk) DX(*),DY(*)
 *     ..
 *
 *  =====================================================================
 *
 *     .. Local Scalars ..
-      DOUBLE PRECISION DTEMP
+      REAL(kind=rk) DTEMP
       INTEGER I,IX,IY,M,MP1
 *     ..
 *     .. Intrinsic Functions ..

external/blas/dtrmm.f

 *       SUBROUTINE DTRMM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
 * 
 *       .. Scalar Arguments ..
-*       DOUBLE PRECISION ALPHA
+*       REAL(kind=rk) ALPHA
 *       INTEGER LDA,LDB,M,N
 *       CHARACTER DIAG,SIDE,TRANSA,UPLO
 *       ..
 *       .. Array Arguments ..
-*       DOUBLE PRECISION A(LDA,*),B(LDB,*)
+*       REAL(kind=rk) A(LDA,*),B(LDB,*)
 *       ..
 *  
 *
 *>
 *> \param[in] ALPHA
 *> \verbatim
-*>          ALPHA is DOUBLE PRECISION.
+*>          ALPHA is REAL(kind=rk).
 *>           On entry,  ALPHA specifies the scalar  alpha. When  alpha is
 *>           zero then  A is not referenced and  B need not be set before
 *>           entry.
 *>
 *> \param[in] A
 *> \verbatim
-*>           A is DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m
+*>           A is REAL(kind=rk) array of DIMENSION ( LDA, k ), where k is m
 *>           when  SIDE = 'L' or 'l'  and is  n  when  SIDE = 'R' or 'r'.
 *>           Before entry  with  UPLO = 'U' or 'u',  the  leading  k by k
 *>           upper triangular part of the array  A must contain the upper
 *>
 *> \param[in,out] B
 *> \verbatim
-*>          B is DOUBLE PRECISION array of DIMENSION ( LDB, n ).
+*>          B is REAL(kind=rk) array of DIMENSION ( LDB, n ).
 *>           Before entry,  the leading  m by n part of the array  B must
 *>           contain the matrix  B,  and  on exit  is overwritten  by the
 *>           transformed matrix.
 *>
 *  =====================================================================
       SUBROUTINE DTRMM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
+      use env_module, only: rk
 *
 *  -- Reference BLAS level3 routine (version 3.4.0) --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
 *     November 2011
 *
 *     .. Scalar Arguments ..
-      DOUBLE PRECISION ALPHA
+      REAL(kind=rk) ALPHA
       INTEGER LDA,LDB,M,N
       CHARACTER DIAG,SIDE,TRANSA,UPLO
 *     ..
 *     .. Array Arguments ..
-      DOUBLE PRECISION A(LDA,*),B(LDB,*)
+      REAL(kind=rk) A(LDA,*),B(LDB,*)
 *     ..
 *
 *  =====================================================================
       INTRINSIC MAX
 *     ..
 *     .. Local Scalars ..
-      DOUBLE PRECISION TEMP
+      REAL(kind=rk) TEMP
       INTEGER I,INFO,J,K,NROWA
       LOGICAL LSIDE,NOUNIT,UPPER
 *     ..
 *     .. Parameters ..
-      DOUBLE PRECISION ONE,ZERO
+      REAL(kind=rk) ONE,ZERO
       PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
 *     ..
 *

external/blas/dtrmv.f

 *       CHARACTER DIAG,TRANS,UPLO
 *       ..
 *       .. Array Arguments ..
-*       DOUBLE PRECISION A(LDA,*),X(*)
+*       REAL(kind=rk) A(LDA,*),X(*)
 *       ..
 *  
 *
 *>
 *> \param[in] A
 *> \verbatim
-*>          A is DOUBLE PRECISION array of DIMENSION ( LDA, n ).
+*>          A is REAL(kind=rk) array of DIMENSION ( LDA, n ).
 *>           Before entry with  UPLO = 'U' or 'u', the leading n by n
 *>           upper triangular part of the array A must contain the upper
 *>           triangular matrix and the strictly lower triangular part of
 *>
 *> \param[in,out] X
 *> \verbatim
-*>          X is DOUBLE PRECISION array of dimension at least
+*>          X is REAL(kind=rk) array of dimension at least
 *>           ( 1 + ( n - 1 )*abs( INCX ) ).
 *>           Before entry, the incremented array X must contain the n
 *>           element vector x. On exit, X is overwritten with the
 *>
 *  =====================================================================
       SUBROUTINE DTRMV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX)
+      use env_module, only: rk
 *
 *  -- Reference BLAS level2 routine (version 3.4.0) --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
       CHARACTER DIAG,TRANS,UPLO
 *     ..
 *     .. Array Arguments ..
-      DOUBLE PRECISION A(LDA,*),X(*)
+      REAL(kind=rk) A(LDA,*),X(*)
 *     ..
 *
 *  =====================================================================
 *
 *     .. Parameters ..
-      DOUBLE PRECISION ZERO
+      REAL(kind=rk) ZERO
       PARAMETER (ZERO=0.0D+0)
 *     ..
 *     .. Local Scalars ..
-      DOUBLE PRECISION TEMP
+      REAL(kind=rk) TEMP
       INTEGER I,INFO,IX,J,JX,KX
       LOGICAL NOUNIT
 *     ..

external/blas/dtrsm.f

 *       SUBROUTINE DTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
 * 
 *       .. Scalar Arguments ..
-*       DOUBLE PRECISION ALPHA
+*       REAL(kind=rk) ALPHA
 *       INTEGER LDA,LDB,M,N
 *       CHARACTER DIAG,SIDE,TRANSA,UPLO
 *       ..
 *       .. Array Arguments ..
-*       DOUBLE PRECISION A(LDA,*),B(LDB,*)
+*       REAL(kind=rk) A(LDA,*),B(LDB,*)
 *       ..
 *  
 *
 *>
 *> \param[in] ALPHA
 *> \verbatim
-*>          ALPHA is DOUBLE PRECISION.
+*>          ALPHA is REAL(kind=rk).
 *>           On entry,  ALPHA specifies the scalar  alpha. When  alpha is
 *>           zero then  A is not referenced and  B need not be set before
 *>           entry.
 *>
 *> \param[in] A
 *> \verbatim
-*>          A is DOUBLE PRECISION array of DIMENSION ( LDA, k ),
+*>          A is REAL(kind=rk) array of DIMENSION ( LDA, k ),
 *>           where k is m when SIDE = 'L' or 'l'  
 *>             and k is n when SIDE = 'R' or 'r'.
 *>           Before entry  with  UPLO = 'U' or 'u',  the  leading  k by k
 *>
 *> \param[in,out] B
 *> \verbatim
-*>          B is DOUBLE PRECISION array of DIMENSION ( LDB, n ).
+*>          B is REAL(kind=rk) array of DIMENSION ( LDB, n ).
 *>           Before entry,  the leading  m by n part of the array  B must
 *>           contain  the  right-hand  side  matrix  B,  and  on exit  is
 *>           overwritten by the solution matrix  X.
 *>
 *  =====================================================================
       SUBROUTINE DTRSM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
+      use env_module, only: rk
 *
 *  -- Reference BLAS level3 routine (version 3.4.0) --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
 *     November 2011
 *
 *     .. Scalar Arguments ..
-      DOUBLE PRECISION ALPHA
+      REAL(kind=rk) ALPHA
       INTEGER LDA,LDB,M,N
       CHARACTER DIAG,SIDE,TRANSA,UPLO
 *     ..
 *     .. Array Arguments ..
-      DOUBLE PRECISION A(LDA,*),B(LDB,*)
+      REAL(kind=rk) A(LDA,*),B(LDB,*)
 *     ..
 *
 *  =====================================================================
       INTRINSIC MAX
 *     ..
 *     .. Local Scalars ..
-      DOUBLE PRECISION TEMP
+      REAL(kind=rk) TEMP
       INTEGER I,INFO,J,K,NROWA
       LOGICAL LSIDE,NOUNIT,UPPER
 *     ..
 *     .. Parameters ..
-      DOUBLE PRECISION ONE,ZERO
+      REAL(kind=rk) ONE,ZERO
       PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
 *     ..
 *

external/blas/idamax.f

 *       INTEGER INCX,N
 *       ..
 *       .. Array Arguments ..
-*       DOUBLE PRECISION DX(*)
+*       REAL(kind=rk) DX(*)
 *       ..
 *  
 *
 *>
 *  =====================================================================
       INTEGER FUNCTION IDAMAX(N,DX,INCX)
+      use env_module, only: rk
 *
 *  -- Reference BLAS level1 routine (version 3.4.0) --
 *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
       INTEGER INCX,N
 *     ..
 *     .. Array Arguments ..
-      DOUBLE PRECISION DX(*)
+      REAL(kind=rk) DX(*)
 *     ..
+*      INTRINSIC DABS
 *
 *  =====================================================================
 *
 *     .. Local Scalars ..
-      DOUBLE PRECISION DMAX
+      REAL(kind=rk) DMAX
       INTEGER I,IX
 *     ..
-*     .. Intrinsic Functions ..
-      INTRINSIC DABS
 *     ..
       IDAMAX = 0
       IF (N.LT.1 .OR. INCX.LE.0) RETURN
 *
 *        code for increment equal to 1
 *
-         DMAX = DABS(DX(1))
+         DMAX = ABS(DX(1))
          DO I = 2,N
-            IF (DABS(DX(I)).GT.DMAX) THEN
+            IF (ABS(DX(I)).GT.DMAX) THEN
                IDAMAX = I
-               DMAX = DABS(DX(I))
+               DMAX = ABS(DX(I))
             END IF
          END DO
       ELSE
 *        code for increment not equal to 1
 *
          IX = 1
-         DMAX = DABS(DX(1))
+         DMAX = ABS(DX(1))
          IX = IX + INCX
          DO I = 2,N
-            IF (DABS(DX(IX)).GT.DMAX) THEN
+            IF (ABS(DX(IX)).GT.DMAX) THEN
                IDAMAX = I
-               DMAX = DABS(DX(IX))
+               DMAX = ABS(DX(IX))
             END IF
             IX = IX + INCX
          END DO

external/lapack/dgetf2.f

 *
 *  =====================================================================
       SUBROUTINE DGETF2( M, N, A, LDA, IPIV, INFO )
+      use env_module, only: rk
 *
 *  -- LAPACK computational routine (version 3.4.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *     ..
 *     .. Array Arguments ..
       INTEGER            IPIV( * )
-      DOUBLE PRECISION   A( LDA, * )
+      REAL(kind=rk)   A( LDA, * )
 *     ..
 *
 *  =====================================================================
 *
 *     .. Parameters ..
-      DOUBLE PRECISION   ONE, ZERO
+      REAL(kind=rk)   ONE, ZERO
       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 *     ..
 *     .. Local Scalars ..
-      DOUBLE PRECISION   SFMIN 
+      REAL(kind=rk)   SFMIN 
       INTEGER            I, J, JP
 *     ..
 *     .. External Functions ..
-      DOUBLE PRECISION   DLAMCH      
+      REAL(kind=rk)   DLAMCH      
       INTEGER            IDAMAX
       EXTERNAL           DLAMCH, IDAMAX
 *     ..

external/lapack/dgetrf.f

 *
 *  =====================================================================
       SUBROUTINE DGETRF( M, N, A, LDA, IPIV, INFO )
+      use env_module, only: rk
 *
 *  -- LAPACK computational routine (version 3.4.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *     ..
 *     .. Array Arguments ..
       INTEGER            IPIV( * )
-      DOUBLE PRECISION   A( LDA, * )
+      REAL(kind=rk)   A( LDA, * )
 *     ..
 *
 *  =====================================================================
 *
 *     .. Parameters ..
-      DOUBLE PRECISION   ONE
+      REAL(kind=rk)   ONE
       PARAMETER          ( ONE = 1.0D+0 )
 *     ..
 *     .. Local Scalars ..

external/lapack/dgetri.f

 *
 *  =====================================================================
       SUBROUTINE DGETRI( N, A, LDA, IPIV, WORK, LWORK, INFO )
+      use env_module, only: rk
 *
 *  -- LAPACK computational routine (version 3.4.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *     ..
 *     .. Array Arguments ..
       INTEGER            IPIV( * )
-      DOUBLE PRECISION   A( LDA, * ), WORK( * )
+      REAL(kind=rk)   A( LDA, * ), WORK( * )
 *     ..
 *
 *  =====================================================================
 *
 *     .. Parameters ..
-      DOUBLE PRECISION   ZERO, ONE
+      REAL(kind=rk)   ZERO, ONE
       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
 *     ..
 *     .. Local Scalars ..

external/lapack/dlamch.f

 *  Definition:
 *  ===========
 *
-*      DOUBLE PRECISION FUNCTION DLAMCH( CMACH )
+*      REAL(kind=rk) FUNCTION DLAMCH( CMACH )
 *  
 *
 *> \par Purpose:
 *> \ingroup auxOTHERauxiliary
 *
 *  =====================================================================
-      DOUBLE PRECISION FUNCTION DLAMCH( CMACH )
+      FUNCTION DLAMCH( CMACH ) result( res_dlamch )
+      use env_module, only:rk
 *
 *  -- LAPACK auxiliary routine (version 3.4.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *
 *     .. Scalar Arguments ..
       CHARACTER          CMACH
+      REAL(kind=rk)      res_dlamch
 *     ..
 *
 * =====================================================================
 *
 *     .. Parameters ..
-      DOUBLE PRECISION   ONE, ZERO
+      REAL(kind=rk)   ONE, ZERO
       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 *     ..
 *     .. Local Scalars ..
-      DOUBLE PRECISION   RND, EPS, SFMIN, SMALL, RMACH
+      REAL(kind=rk)   RND, EPS, SFMIN, SMALL, RMACH
 *     ..
 *     .. External Functions ..
       LOGICAL            LSAME
          RMACH = ZERO
       END IF
 *
-      DLAMCH = RMACH
+      res_dlamch = RMACH
       RETURN
 *
 *     End of DLAMCH
 *>
 *> \param[in] A
 *> \verbatim
-*>          A is a DOUBLE PRECISION
+*>          A is a REAL(kind=rk)
 *> \endverbatim
 *>
 *> \param[in] B
 *> \verbatim
-*>          B is a DOUBLE PRECISION
+*>          B is a REAL(kind=rk)
 *>          The values A and B.
 *> \endverbatim
 *>
-      DOUBLE PRECISION FUNCTION DLAMC3( A, B )
+      FUNCTION DLAMC3( A, B ) result ( res_dlamc3 )
+      use env_module, only:rk
 *
 *  -- LAPACK auxiliary routine (version 3.4.0) --
 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 *     November 2010
 *
 *     .. Scalar Arguments ..
-      DOUBLE PRECISION   A, B
+      REAL(kind=rk)   A, B
+      REAL(kind=rk)   res_dlamc3
 *     ..
 * =====================================================================
 *
 *     .. Executable Statements ..
 *
-      DLAMC3 = A + B
+      res_dlamc3 = A + B
 *
       RETURN
 *

external/lapack/dlaswp.f

 *       ..
 *       .. Array Arguments ..
 *       INTEGER            IPIV( * )
-*       DOUBLE PRECISION   A( LDA, * )
+*       REAL(kind=rk)   A( LDA, * )
 *       ..
 *  
 *
 *>
 *> \param[in,out] A
 *> \verbatim
-*>          A is DOUBLE PRECISION array, dimension (LDA,N)
+*>          A is REAL(kind=rk) array, dimension (LDA,N)
 *>          On entry, the matrix of column dimension N to which the row
 *>          interchanges will be applied.
 *>          On exit, the permuted matrix.
 *>
 *  =====================================================================
       SUBROUTINE DLASWP( N, A, LDA, K1, K2, IPIV, INCX )
+      use env_module, only: rk
 *
 *  -- LAPACK auxiliary routine (version 3.4.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *     ..
 *     .. Array Arguments ..
       INTEGER            IPIV( * )
-      DOUBLE PRECISION   A( LDA, * )
+      REAL(kind=rk)   A( LDA, * )
 *     ..
 *
 * =====================================================================
 *
 *     .. Local Scalars ..
       INTEGER            I, I1, I2, INC, IP, IX, IX0, J, K, N32
-      DOUBLE PRECISION   TEMP
+      REAL(kind=rk)   TEMP
 *     ..
 *     .. Executable Statements ..
 *

external/lapack/dtrti2.f

 *       INTEGER            INFO, LDA, N
 *       ..
 *       .. Array Arguments ..
-*       DOUBLE PRECISION   A( LDA, * )
+*       REAL(kind=rk)   A( LDA, * )
 *       ..
 *  
 *
 *>
 *> \param[in,out] A
 *> \verbatim
-*>          A is DOUBLE PRECISION array, dimension (LDA,N)
+*>          A is REAL(kind=rk) array, dimension (LDA,N)
 *>          On entry, the triangular matrix A.  If UPLO = 'U', the
 *>          leading n by n upper triangular part of the array A contains
 *>          the upper triangular matrix, and the strictly lower
 *
 *  =====================================================================
       SUBROUTINE DTRTI2( UPLO, DIAG, N, A, LDA, INFO )
+      use env_module, only: rk
 *
 *  -- LAPACK computational routine (version 3.4.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
       INTEGER            INFO, LDA, N
 *     ..
 *     .. Array Arguments ..
-      DOUBLE PRECISION   A( LDA, * )
+      REAL(kind=rk)   A( LDA, * )
 *     ..
 *
 *  =====================================================================
 *
 *     .. Parameters ..
-      DOUBLE PRECISION   ONE
+      REAL(kind=rk)   ONE
       PARAMETER          ( ONE = 1.0D+0 )
 *     ..
 *     .. Local Scalars ..
       LOGICAL            NOUNIT, UPPER
       INTEGER            J
-      DOUBLE PRECISION   AJJ
+      REAL(kind=rk)   AJJ
 *     ..
 *     .. External Functions ..
       LOGICAL            LSAME

external/lapack/dtrtri.f

 *       INTEGER            INFO, LDA, N
 *       ..
 *       .. Array Arguments ..
-*       DOUBLE PRECISION   A( LDA, * )
+*       REAL(kind=rk)   A( LDA, * )
 *       ..
 *  
 *
 *>
 *> \param[in,out] A
 *> \verbatim
-*>          A is DOUBLE PRECISION array, dimension (LDA,N)
+*>          A is REAL(kind=rk) array, dimension (LDA,N)
 *>          On entry, the triangular matrix A.  If UPLO = 'U', the
 *>          leading N-by-N upper triangular part of the array A contains
 *>          the upper triangular matrix, and the strictly lower
 *
 *  =====================================================================
       SUBROUTINE DTRTRI( UPLO, DIAG, N, A, LDA, INFO )
+      use env_module, only:rk
 *
 *  -- LAPACK computational routine (version 3.4.0) --
 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
       INTEGER            INFO, LDA, N
 *     ..
 *     .. Array Arguments ..
-      DOUBLE PRECISION   A( LDA, * )
+      REAL(kind=rk)   A( LDA, * )
 *     ..
 *
 *  =====================================================================
 *
 *     .. Parameters ..
-      DOUBLE PRECISION   ONE, ZERO
+      REAL(kind=rk)   ONE, ZERO
       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 *     ..
 *     .. Local Scalars ..

source/arrayMacros.inc

     type(grw_?tname?Array_type), intent(out) :: me !< Dynamic Array to init
     integer, intent(in), optional :: length !< Initial length of the container
 
-    if (present(length)) me%ContainerSize = length
-
+    if (present(length)) then
+      me%ContainerSize = length
+    else
+      me%ContainerSize = zeroLength
+    end if
+    ! Deallocate ...
+    if( allocated( me%Val ))     &
+      deallocate(me%Val)
+    ! ... and reallocate
     allocate(me%Val(me%ContainerSize))
 
   end subroutine init_GA_?tname?
 
     allocate(swpval(me%nVals))
     swpval(:me%nVals) = me%Val(:me%nVals)
+    if( allocated(me%Val)) &
     deallocate(me%Val)
     if( present( pos ) ) then
       me%ContainerSize = max(minLength, me%ContainerSize*2, pos)
     integer, intent(in), optional :: length !< Initial length of the container
     !-----------------------------------------------------------------
 
-    if (present(length)) me%ContainerSize = length
+    if (present(length)) then
+      me%ContainerSize = length
+    else
+      me%ContainerSize = zeroLength
+    end if
 
     me%Unique = unique
+
+    ! Deallocate ...
+    if( allocated( me%Val ))     &
+      deallocate(me%Val)
+    if( allocated( me%sorted ))     &
+      deallocate(me%sorted)
+    ! ... and reallocate
     allocate(me%Val(me%ContainerSize))
     allocate(me%sorted(me%ContainerSize))
+
   end subroutine init_DA_?tname?
 
   !> Destruction of a dynamic array
     !------------------------------------------------------------------------
     type(dyn_?tname?Array_type) :: me !< Array to sortTruncate
     !------------------------------------------------------------------------
-    type(dyn_?tname?Array_type) :: tArray !< temporary Array
+    type(dyn_?tname?Array_type) :: tArray ! temporary Array
     integer :: iVal
     integer :: dPos
     !------------------------------------------------------------------------
     !------------------------------------------------------------------------
     type(dyn_?tname?Array_type) :: me !< Array to sortTruncate
     !------------------------------------------------------------------------
-    type(dyn_?tname?Array_type) :: tArray !< temporary Array
+    type(dyn_?tname?Array_type) :: tArray ! temporary Array
     integer :: iVal
     integer :: dPos
     !------------------------------------------------------------------------
 
     allocate(swpsort(me%nVals))
     swpsort = me%sorted(:me%nVals)
-    deallocate(me%sorted)
+    if( allocated( me%sorted)) &
+      deallocate(me%sorted)
     allocate(me%sorted(me%ContainerSize))
     me%sorted(1:me%nVals) = swpsort
     deallocate(swpsort)

source/atl_derived/README

-This equation directory contains the descriptions for the various supported
-equation systems of the ateles solver.
-The general structure for the equation system in the solver is given by
-atl_equation_module, which is supported by a set of routines acting on this
-data structure given in atl_eqn_helpers_module.
-All the modules are supposed to have the atl_eqn_ prefix to easily identify
-the equation specific things.
-Each specific equation system should get its own descriptive module and in
-addition some supporting helpers module.
-The variables of the different equation types are defined in
-atl_eqn_*_var_module.
-Routines to calculate derived quantities of the different equation types are
-defined in atl_eqn_*_derive_module.

source/atl_derived/atl_eqn_advdiff_module.f90

-! See copyright notice in the ../COPYRIGHT file.
-module atl_eqn_advdiff_module
-  use env_module, only: rk
-  use aotus_module, only: flu_State, aot_get_val
-
-  implicit none
-
-  private
-
-  public :: AdvectionDiffusion_type
-  public :: init_AdvectionDiffusion
-
-  !> The Advection-Diffusion equation properties  are stored here
-  type AdvectionDiffusion_type
-    real(kind=rk) :: velX   !< transport propery in x
-    real(kind=rk) :: velY   !< transport propery in y
-    real(kind=rk) :: velZ   !< transport propery in z
-    real(kind=rk) :: dfsn  !< diffusion coefficient
-    real(kind=rk) :: sigma  !< stability
-    character(len=128) :: TimeStep
-    integer :: TimeStepSize
-  end type AdvectionDiffusion_type
-
-contains
-
-  subroutine init_AdvectionDiffusion(AdvectionDiffusion, conf, ic_table)
-    !--------------------------------------------------------------------------!
-    type(AdvectionDiffusion_type), intent(out) :: AdvectionDiffusion
-    type(flu_State)                            :: conf
-    integer, intent(in)                        :: ic_table
-    !--------------------------------------------------------------------------!
-    integer                                    :: iError
-    !--------------------------------------------------------------------------!
-
-    !HK: was severly broken and opened initial condition table and scheme table?
-    !HK: Now just reading all these values from the table, given...
-    !HK: Propably does not work right now !!
-    !read the data from the equation table of the lua file
-    call aot_get_val(L = conf, thandle = ic_table, &
-      &                key = 'velocityX', &
-      &                val = AdvectionDiffusion%velX, &
-      &                ErrCode = iError)
-    call aot_get_val(L = conf, thandle = ic_table, &
-      &                key = 'velocityY', &
-      &                val = AdvectionDiffusion%velY, &
-      &                ErrCode = iError)
-    call aot_get_val(L = conf, thandle = ic_table, &
-      &                key = 'velocityZ', &
-      &                val = AdvectionDiffusion%velZ, &
-      &                ErrCode = iError)
-    call aot_get_val(L = conf, thandle = ic_table, &
-      &                key = 'diffusion', &
-      &                val = AdvectionDiffusion%dfsn, &
-      &                ErrCode = iError)
-    call aot_get_val(L = conf, thandle = ic_table, &
-      &                key = 'sigma', &
-      &                val = AdvectionDiffusion%sigma, &
-      &                ErrCode = iError)
-    call aot_get_val(L = conf, thandle = ic_table, &
-      &                key = 'TimeStep', &
-      &                val = AdvectionDiffusion%TimeStep, &
-      &                ErrCode = iError)
-    call aot_get_val(L = conf, thandle = ic_table, &
-      &                key = 'TimeStepSize', &
-      &                val = AdvectionDiffusion%TimeStepSize, &
-      &                ErrCode = iError)
-
-  end subroutine init_AdvectionDiffusion
-
-
-end module atl_eqn_advdiff_module

source/atl_derived/atl_eqn_bbm_module.f90

-! See copyright notice in the ../COPYRIGHT file.
-!> Module describing a simple, basic membrane model
-module atl_eqn_bbm_module
-  use env_module, only: rk
-  use aotus_module, only: flu_State, aot_get_val
-
-  use tem_tools_module, only: tem_write
-
-  implicit none
-
-  private
-
-  public :: BBMEM_type
-  public :: init_BBMEM
-
-  !> Parameters describing the model
-  type BBMEM_type
-    real(kind=rk) :: tna
-    real(kind=rk) :: tcl
-    real(kind=rk) :: kappa
-  end type BBMEM_type
-
-contains
-
-  !> \brief subroutine to intialize BBM
-  !! material parameters are to be filled here.
-  subroutine init_BBMEM(BBM, conf, eq_table)
-    !--------------------------------------------------------------------------!
-    type(BBMEM_type), intent(out) :: BBM
-    type(flu_State)             :: conf
-    integer, intent(in)         :: eq_table
-    !--------------------------------------------------------------------------!
-    integer :: iError
-    !--------------------------------------------------------------------------!
-
-    call tem_write('Initializing BBM')
-    call aot_get_val(L = conf, thandle = eq_table, key = 'tna', &
-      &              val = BBM%tna, &
-      &              ErrCode = iError)
-    call tem_write('  tna  : ', BBM%tna)
-
-    call aot_get_val(L = conf, thandle = eq_table, key = 'tcl', &
-      &              val = BBM%tcl, &
-      &              ErrCode = iError)
-    call tem_write('  tcl  : ', BBM%tcl)
-
-    call aot_get_val(L = conf, thandle = eq_table, key = 'kappa', &
-      &              val = BBM%kappa, &
-      &              ErrCode = iError)
-    call tem_write('  kappa: ', BBM%kappa)
-
-  end subroutine init_BBMEM
-
-end module atl_eqn_bbm_module

source/atl_derived/atl_eqn_euler_derive_module.f90

-! See copyright notice in the ../COPYRIGHT file.
-!> Routines to derive quantities from the state in the Euler equation system.
-module atl_eqn_euler_derive_module
-  use env_module,             only: rk
-  use aotus_module,           only: flu_State, aot_get_val
-  use aot_table_module,       only: aot_table_open, aot_table_close
-
-  use tem_aux_module,         only: tem_abort
-  use tem_tools_module,       only: tem_write
-  use tem_var_system_module,  only: tem_var_system_type, &
-    &                               tem_var_append
-  use tem_stencil_module,     only: tem_stencilHeader_type
-  use tem_time_module,        only: tem_time_type
-  use tem_variable_module,    only: tem_variable_type
-  use treelmesh_module, only: treelmesh_type
-
-  use atl_eqn_euler_module,   only: euler_type, init_euler
-  use atl_equation_module,    only: equations_type
-
-  implicit none
-
-  private
-
-  public :: deriveVelocityFromConsVar
-  public :: derivePressureFromConsVar
-  public :: deriveTemperatureFromConsVar
-  public :: deriveMachNumberFromConsVar
-  public :: deriveMachVectorFromConsVar
-  public :: deriveKineticEnergyFromConsVar
-
-  public :: eqn_euler_cons2prim
-  public :: eqn_euler_prim2cons
-
-contains
-
-
-  !> Calculate the velocity from conservative Euler variables.
-  !!
-  !! The interface has to comply to the abstract interface
-  !! tem_var_system_module#derive.
-  !! See there for a more detailed description of the parameters.
-  subroutine deriveVelocityFromConsVar( input, nElems, tree, requVar, globSys, &
-    &                                stencil, res, conf, bary, minIDbound, time)
-    !------------------------------------------------------------------------
-    real(kind=rk), target, intent(inout)         :: input(:)
-    integer, intent(in)                          :: nElems
-    !> tree to work on (e.g. global or local from tracking)
-    type(treelmesh_type), intent(in) :: tree
-    type(tem_variable_type), intent(in)          :: requVar
-    type(tem_var_system_type), intent(in)        :: globSys
-    type(tem_stencilHeader_type)                       :: stencil
-    real(kind=rk), optional, target, intent(out) :: res(:)
-    type(flu_state)                              :: conf
-    real(kind=rk), optional, intent(in)          :: bary(:,:)
-    !> starting bound for the barycenters and the treeID list in the tree in
-    !! case of a chunk usage (needed for spacetime functions and meshInfo)
-    integer, optional, intent(in) :: minIDbound
-    type(tem_time_type), optional, intent(in)    :: time
-    !------------------------------------------------------------------------
-    !local variable
-    integer :: iElem
-    integer :: offset, offsetRes
-    integer :: iDensity, iMomentum(3)
-    integer :: pos_Density, pos_Momentum(3)
-    real(kind=rk), pointer, dimension(:) :: local_res => NULL()
-    !------------------------------------------------------------------------
-
-    local_res => input
-    ! If res present, no in-place computation.
-    if (present(res)) local_res => res
-
-    call get_cvarPos(varsys       = globSys%variable%val, &
-      &              requVar      = requVar%mapping, &
-      &              pos_Density  = pos_Density, &
-      &              pos_Momentum = pos_Momentum)
-
-    ! Offsets in the serialized data array
-    offset = 0 ! offset in the input data
-    offsetRes = 0 ! offset in the output data
-    do iElem = 1, nElems
-      iDensity  = offset + pos_Density
-      iMomentum = offset + pos_Momentum
-
-      ! Compute velocity as momentum/density
-      local_res(offsetRes+1:offsetRes+3) = input(iMomentum) / input(iDensity)
-
-      ! Advance the input array by the number of scalars in that system.
-      offset = offset + globSys%nScalars
-
-      ! Advance the output offset by the 3 components written with this
-      ! velocity computation.
-      offsetRes = offsetRes + 3
-    end do
-
-  end subroutine deriveVelocityFromConsVar
-
-
-  !> Calculate the pressure from conservative input
-  !!
-  !! The interface has to comply to the abstract interface
-  !! tem_var_system_module#derive.
-  !! See there for a more detailed description of the parameters.
-  subroutine derivePressureFromConsVar( input, nElems, tree, requVar, globSys, &
-    &                               stencil, res, conf, bary, minIDbound, time )
-    !------------------------------------------------------------------------
-    real(kind=rk), target, intent(inout)         :: input(:)
-    integer, intent(in)                          :: nElems
-    !> tree to work on (e.g. global or local from tracking)
-    type(treelmesh_type), intent(in) :: tree
-    type(tem_variable_type), intent(in)          :: requVar
-    type(tem_var_system_type), intent(in)        :: globSys
-    type(tem_stencilHeader_type)                       :: stencil
-    real(kind=rk), optional, target, intent(out) :: res(:)
-    type(flu_state)                              :: conf
-    real(kind=rk), optional, intent(in)          :: bary(:,:)
-    !> starting bound for the barycenters and the treeID list in the tree in
-    !! case of a chunk usage (needed for spacetime functions and meshInfo)
-    integer, optional, intent(in) :: minIDbound
-    type(tem_time_type), optional, intent(in)    :: time
-    !------------------------------------------------------------------------
-    !local variable
-    integer :: iElem
-    integer :: offset
-    integer :: iDensity, iMomentum(3), iEnergy
-    integer :: pos_Density, pos_Momentum(3), pos_Energy
-    real(kind=rk), pointer, dimension(:) :: local_res => NULL()
-    type(Euler_type) :: euler
-    !------------------------------------------------------------------------
-
-    local_res => input
-    ! If res present, no in-place computation.
-    if (present(res)) local_res => res
-
-    call get_EulerParams(conf, euler)
-
-    call get_cvarPos(varsys       = globSys%variable%val, &
-      &              requVar      = requVar%mapping, &
-      &              pos_Density  = pos_Density, &
-      &              pos_Momentum = pos_Momentum, &
-      &              pos_Energy   = pos_Energy)
-
-    ! Offsets in the serialized data array
-    offset = 0 ! Offset in the input
-    ! No need to track offset of output, as pressure is scalar and the
-    ! resulting offset is equal to iElem.
-    do iElem = 1, nElems
-      iDensity  = offset + pos_Density
-      iMomentum = offset + pos_Momentum
-      iEnergy   = offset + pos_Energy
-
-      local_res(iElem) = (euler%isen_coef - 1.0_rk) * &
-        &            (input(iEnergy) - 0.5_rk / input(iDensity) &
-        &                               * sum(input(iMomentum)**2))
-
-      offset = offset + globSys%nScalars
-    end do
-
-  end subroutine derivePressureFromConsVar
-
-
-  !> Calculate the temperature from conservative input
-  !!
-  !! The interface has to comply to the abstract interface
-  !! tem_var_system_module#derive.
-  !! See there for a more detailed description of the parameters.
-  subroutine deriveTemperatureFromConsVar( input, nElems, tree, requVar,       &
-    &                                      globSys, stencil, res, conf, bary,  &
-    &                                      minIDbound, time )
-    !------------------------------------------------------------------------
-    real(kind=rk), target, intent(inout)         :: input(:)
-    integer, intent(in)                          :: nElems
-    !> tree to work on (e.g. global or local from tracking)
-    type(treelmesh_type), intent(in) :: tree
-    type(tem_variable_type), intent(in)          :: requVar
-    type(tem_var_system_type), intent(in)        :: globSys
-    type(tem_stencilHeader_type)                       :: stencil
-    real(kind=rk), optional, target, intent(out) :: res(:)
-    type(flu_state)                              :: conf
-    real(kind=rk), optional, intent(in)          :: bary(:,:)
-    !> starting bound for the barycenters and the treeID list in the tree in
-    !! case of a chunk usage (needed for spacetime functions and meshInfo)
-    integer, optional, intent(in) :: minIDbound
-    type(tem_time_type), optional, intent(in)    :: time
-    !------------------------------------------------------------------------
-    !local variable
-    integer :: iElem
-    integer :: offset
-    integer :: iDensity, iMomentum(3), iEnergy
-    integer :: pos_Density, pos_Momentum(3), pos_Energy
-    real(kind=rk), pointer, dimension(:) :: local_res => NULL()
-    real(kind=rk) :: pressure
-    real(kind=rk) :: R_q
-    type(Euler_type) :: euler
-    !------------------------------------------------------------------------
-
-    local_res => input
-    ! If res present, no in-place computation.
-    if (present(res)) local_res => res
-
-    call get_EulerParams(conf, euler)
-
-    call get_cvarPos(varsys = globSys%variable%val, &
-      &              requVar = requVar%mapping, &
-      &              pos_Density = pos_Density, &
-      &              pos_Momentum = pos_Momentum, &
-      &              pos_Energy = pos_Energy)
-
-    R_q = 1.0_rk / euler%R
-
-    ! Offsets in the serialized data array
-    offset = 0 ! Offset in the input
-    ! No need to track offset of output, as temperature is scalar and the
-    ! resulting offset is equal to iElem.
-    do iElem = 1, nElems
-      iDensity  = offset + pos_Density
-      iMomentum = offset + pos_Momentum
-      iEnergy   = offset + pos_Energy
-
-      pressure = (euler%isen_coef - 1.0_rk) * &
-        &            (input(iEnergy) - 0.5_rk / input(iDensity) &
-        &                                     * sum(input(iMomentum)**2))
-
-      local_res(iElem) = R_q * pressure / input(iDensity)
-
-      offset = offset + globSys%nScalars
-    end do
-
-  end subroutine deriveTemperatureFromConsVar
-
-
-  !> Calculate the speed of sound from conservative input
-  !!
-  !! The interface has to comply to the abstract interface
-  !! tem_var_system_module#derive.
-  !! See there for a more detailed description of the parameters.
-  subroutine deriveSpeedOfSoundFromConsVar(input, nElems, tree, requVar, &
-    &                                      globSys, stencil, res, conf, bary, &
-    &                                      minIDbound, time)
-    !------------------------------------------------------------------------
-    real(kind=rk), target, intent(inout)         :: input(:)
-    integer, intent(in)                          :: nElems
-    !> tree to work on (e.g. global or local from tracking)
-    type(treelmesh_type), intent(in) :: tree
-    type(tem_variable_type), intent(in)          :: requVar
-    type(tem_var_system_type), intent(in)        :: globSys
-    type(tem_stencilHeader_type)                       :: stencil
-    real(kind=rk), optional, target, intent(out) :: res(:)
-    type(flu_state)                              :: conf
-    real(kind=rk), optional, intent(in)          :: bary(:,:)
-    !> starting bound for the barycenters and the treeID list in the tree in
-    !! case of a chunk usage (needed for spacetime functions and meshInfo)
-    integer, optional, intent(in) :: minIDbound
-    type(tem_time_type), optional, intent(in)    :: time
-    !------------------------------------------------------------------------
-    !local variable
-    integer :: iElem
-    integer :: offset
-    integer :: iDensity, iMomentum(3), iEnergy
-    integer :: pos_Density, pos_Momentum(3), pos_Energy
-    real(kind=rk), pointer, dimension(:) :: local_res => NULL()
-    real(kind=rk) :: pressure
-    type(Euler_type) :: euler
-    !------------------------------------------------------------------------
-
-    local_res => input
-    ! If res present, no in-place computation.
-    if (present(res)) local_res => res
-
-    call get_EulerParams(conf, euler)
-
-    call get_cvarPos(varsys       = globSys%variable%val, &
-      &              requVar      = requVar%mapping, &
-      &              pos_Density  = pos_Density, &
-      &              pos_Momentum = pos_Momentum, &
-      &              pos_Energy   = pos_Energy)
-
-    ! Offsets in the serialized data array
-    offset = 0 ! Offset in the input
-    ! No need to track offset of output, as speed of sound is scalar and the
-    ! resulting offset is equal to iElem.
-    do iElem = 1, nElems
-      iDensity  = offset + pos_Density
-      iMomentum = offset + pos_Momentum
-      iEnergy   = offset + pos_Energy
-
-      pressure = (euler%isen_coef - 1.0_rk) * &
-        &            (input(iEnergy) - 0.5_rk / input(iDensity) &
-        &                                     * sum(input(iMomentum)**2))
-
-      local_res(iElem) = sqrt(pressure/input(iDensity)*euler%isen_coef)
-
-      offset = offset + globSys%nScalars
-    end do
-
-  end subroutine deriveSpeedOfSoundFromConsVar
-
-
-  !> Calculate the Mach number from conservative input
-  !!
-  !! The interface has to comply to the abstract interface
-  !! tem_var_system_module#derive.
-  !! See there for a more detailed description of the parameters.
-  subroutine deriveMachNumberFromConsVar(input, nElems, tree, requVar, &
-    &                                    globSys, stencil, res, conf, bary, &
-    &                                    minIDbound, time)
-    !------------------------------------------------------------------------
-    real(kind=rk), target, intent(inout)         :: input(:)
-    integer, intent(in)                          :: nElems
-    !> tree to work on (e.g. global or local from tracking)
-    type(treelmesh_type), intent(in) :: tree
-    type(tem_variable_type), intent(in)          :: requVar
-    type(tem_var_system_type), intent(in)        :: globSys
-    type(tem_stencilHeader_type)                       :: stencil
-    real(kind=rk), optional, target, intent(out) :: res(:)
-    type(flu_state)                              :: conf
-    real(kind=rk), optional, intent(in)          :: bary(:,:)
-    !> starting bound for the barycenters and the treeID list in the tree in
-    !! case of a chunk usage (needed for spacetime functions and meshInfo)
-    integer, optional, intent(in) :: minIDbound
-    type(tem_time_type), optional, intent(in)    :: time
-    !------------------------------------------------------------------------
-    !local variable
-    integer :: iElem
-    integer :: offset
-    integer :: iDensity, iMomentum(3), iEnergy
-    integer :: pos_Density, pos_Momentum(3), pos_Energy
-    real(kind=rk), pointer, dimension(:) :: local_res => NULL()
-    real(kind=rk) :: pressure, speedOfSound
-    type(Euler_type) :: euler
-    !------------------------------------------------------------------------
-
-    local_res => input
-    ! If res present, no in-place computation.
-    if (present(res)) local_res => res
-
-    call get_EulerParams(conf, euler)
-
-    call get_cvarPos(varsys       = globSys%variable%val, &
-      &              requVar      = requVar%mapping, &
-      &              pos_Density  = pos_Density, &
-      &              pos_Momentum = pos_Momentum, &
-      &              pos_Energy   = pos_Energy)
-
-    ! Offsets in the serialized data array
-    offset = 0 ! Offset in the input
-    ! No need to track offset of output, as speed of sound is scalar and the
-    ! resulting offset is equal to iElem.
-    do iElem = 1, nElems
-      iDensity  = offset + pos_Density
-      iMomentum = offset + pos_Momentum
-      iEnergy   = offset + pos_Energy
-
-      pressure = (euler%isen_coef - 1.0_rk) * &
-        &            (input(iEnergy) - 0.5_rk / input(iDensity) &
-        &                                     * sum(input(iMomentum)**2))
-
-      speedOfSound = sqrt(pressure/input(iDensity)*euler%isen_coef)
-
-      local_res(iElem) = sqrt(sum((input(iMomentum)/input(iDensity))**2)) &
-        &           / speedOfSound
-
-      offset = offset + globSys%nScalars
-    end do
-
-  end subroutine deriveMachNumberFromConsVar
-
-
-  !> Calculate the Mach vector (v/c) from conservative input
-  !!
-  !! The interface has to comply to the abstract interface
-  !! tem_var_system_module#derive.
-  !! See there for a more detailed description of the parameters.
-  subroutine deriveMachVectorFromConsVar(input, nElems, tree, requVar, &
-    &                                    globSys, stencil, res, conf, bary, &
-    &                                    minIDbound, time)
-    !------------------------------------------------------------------------
-    real(kind=rk), target, intent(inout)         :: input(:)
-    integer, intent(in)                          :: nElems
-    !> tree to work on (e.g. global or local from tracking)
-    type(treelmesh_type), intent(in) :: tree
-    type(tem_variable_type), intent(in)          :: requVar
-    type(tem_var_system_type), intent(in)        :: globSys
-    type(tem_stencilHeader_type)                       :: stencil
-    real(kind=rk), optional, target, intent(out) :: res(:)
-    type(flu_state)                              :: conf
-    real(kind=rk), optional, intent(in)          :: bary(:,:)
-    !> starting bound for the barycenters and the treeID list in the tree in
-    !! case of a chunk usage (needed for spacetime functions and meshInfo)
-    integer, optional, intent(in) :: minIDbound
-    type(tem_time_type), optional, intent(in)    :: time
-    !------------------------------------------------------------------------
-    !local variable
-    integer :: iElem
-    integer :: offset, resOffset
-    integer :: iDensity, iMomentum(3), iEnergy
-    integer :: pos_Density, pos_Momentum(3), pos_Energy
-    real(kind=rk), pointer, dimension(:) :: local_res => NULL()
-    real(kind=rk) :: pressure, speedOfSound
-    type(Euler_type) :: euler
-    !------------------------------------------------------------------------
-
-    local_res => input
-    ! If res present, no in-place computation.
-    if (present(res)) local_res => res
-
-    call get_EulerParams(conf, euler)
-
-    call get_cvarPos(varsys       = globSys%variable%val, &
-      &              requVar      = requVar%mapping, &
-      &              pos_Density  = pos_Density, &
-      &              pos_Momentum = pos_Momentum, &
-      &              pos_Energy   = pos_Energy)
-
-    offset = 0
-    resOffset = 0
-    do iElem = 1, nElems
-      iDensity  = offset + pos_Density
-      iMomentum = offset + pos_Momentum
-      iEnergy   = offset + pos_Energy
-
-      pressure = (euler%isen_coef - 1.0_rk) * &
-        &            (input(iEnergy) - 0.5_rk / input(iDensity) &
-        &                                     * sum(input(iMomentum)**2))
-
-      speedOfSound = sqrt(pressure/input(iDensity)*euler%isen_coef)
-
-      local_res(resOffset+1:resOffset+3) = input(iMomentum)/input(iDensity) &
-        &           / speedOfSound
-
-      offset = offset + globSys%nScalars
-    end do
-
-  end subroutine deriveMachVectorFromConsVar
-
-
-  !> Calculate the kinetic energy from conservative input
-  !!
-  !! The interface has to comply to the abstract interface
-  !! tem_var_system_module#derive.
-  !! See there for a more detailed description of the parameters.
-  subroutine deriveKineticEnergyFromConsVar( input, nElems, tree, requVar, &
-    &        globSys, stencil, res, conf, bary, minIDbound, time )
-    !------------------------------------------------------------------------
-    real(kind=rk), target, intent(inout)         :: input(:)
-    integer, intent(in)                          :: nElems
-    !> tree to work on (e.g. global or local from tracking)
-    type(treelmesh_type), intent(in) :: tree
-    type(tem_variable_type), intent(in)          :: requVar
-    type(tem_var_system_type), intent(in)        :: globSys
-    type(tem_stencilHeader_type)                       :: stencil
-    real(kind=rk), optional, target, intent(out) :: res(:)
-    type(flu_state)                              :: conf
-    real(kind=rk), optional, intent(in)          :: bary(:,:)
-    !> starting bound for the barycenters and the treeID list in the tree in
-    !! case of a chunk usage (needed for spacetime functions and meshInfo)
-    integer, optional, intent(in) :: minIDbound
-    type(tem_time_type), optional, intent(in)    :: time
-    !------------------------------------------------------------------------
-    !local variable
-    integer :: iElem
-    integer :: offset
-    integer :: iDensity, iMomentum(3)
-    integer :: pos_Density, pos_Momentum(3)
-    real(kind=rk), pointer, dimension(:) :: local_res => NULL()
-    !------------------------------------------------------------------------
-
-    local_res => input
-    ! If res present, no in-place computation.
-    if (present(res)) local_res => res
-
-    call get_cvarPos(varsys       = globSys%variable%val, &
-      &              requVar      = requVar%mapping, &
-      &              pos_Density  = pos_Density, &
-      &              pos_Momentum = pos_Momentum)
-
-    offset = 0
-    do iElem = 1, nElems
-      iDensity  = offset + pos_Density
-      iMomentum = offset + pos_Momentum
-
-      local_res(iElem) = 0.5_rk*sum(input(iMomentum)**2) / input(iDensity)
-
-      offset = offset + globSys%nScalars
-    end do
-
-  end subroutine deriveKineticEnergyFromConsVar
-
-
-  !> Convert primitive varibales to conservative variables.
-  !!
-  !! The interface has to comply to the abstract interface
-  !! atl_equation_module#eqn_var_trafo "eqn_var_trafo".
-  subroutine eqn_euler_prim2cons(equation, instate, outstate, nElems)
-    !--------------------------------------------------------------------------!
-    !> Description of the equation system.
-    class(equations_type), intent(in) :: equation
-
-    !> Primitive variables to convert.
-    real(kind=rk), intent(in)         :: instate(:,:,:)
-
-    !> Converted variables.
-    real(kind=rk), intent(out)        :: outstate(:,:,:)
-
-    !> Number of elements to act on (first index in the state arrays).
-    integer, intent(in)               :: nElems
-    !--------------------------------------------------------------------------!
-
-    outstate(1:nElems,:,1) = instate(1:nElems,:,1)
-    outstate(1:nElems,:,2) = instate(1:nElems,:,1)*instate(1:nElems,:,2)
-    outstate(1:nElems,:,3) = instate(1:nElems,:,1)*instate(1:nElems,:,3)
-    outstate(1:nElems,:,4) = instate(1:nElems,:,1)*instate(1:nElems,:,4)
-    outstate(1:nElems,:,5) = instate(1:nElems,:,5) &
-      &                      / (equation%euler%isen_coef-1.0_rk) &
-      &                    + 0.5_rk*instate(1:nElems,:,1) &
-      &                            * ( instate(1:nElems,:,2)**2 &
-      &                               + instate(1:nElems,:,3)**2 &
-      &                               + instate(1:nElems,:,4)**2 )
-  end subroutine eqn_euler_prim2cons
-
-
-  !> Convert conservative to primitive variables.
-  !!
-  !! The interface has to comply to the abstract interface
-  !! atl_equation_module#eqn_var_trafo "eqn_var_trafo".
-  subroutine eqn_euler_cons2prim(equation, instate, outstate, nElems)
-    !--------------------------------------------------------------------------!
-    !> Description of the equation system.
-    class(equations_type), intent(in) :: equation
-
-    !> Primitive variables to convert.
-    real(kind=rk), intent(in)         :: instate(:,:,:)
-
-    !> Converted variables.
-    real(kind=rk), intent(out)        :: outstate(:,:,:)
-
-    !> Number of elements to act on (first index in the state arrays).
-    integer, intent(in)               :: nElems
-    !--------------------------------------------------------------------------!
-
-    outstate(1:nElems,:,1) = instate(1:nElems,:,1)
-    outstate(1:nElems,:,2) = instate(1:nElems,:,2)/instate(1:nElems,:,1)
-    outstate(1:nElems,:,3) = instate(1:nElems,:,3)/instate(1:nElems,:,1)
-    outstate(1:nElems,:,4) = instate(1:nElems,:,4)/instate(1:nElems,:,1)
-    outstate(1:nElems,:,5) = (equation%euler%isen_coef - 1.0_rk) &
-      &           * ( instate(1:nElems,:,5) - 0.5_rk * instate(1:nElems,:,1) &
-      &                                     * (outstate(1:nElems,:,2)**2 &
-      &                                        + outstate(1:nElems,:,3)**2 &
-      &                                        + outstate(1:nElems,:,4)**2) &
-      &             )
-  end subroutine eqn_euler_cons2prim
-
-
-  !> Internal helping subroutine to obtain positions of the variables in the
-  !! conservative state of the system.
-  !!
-  !! This routine exploits the implicit knowledge, that most derived quantities
-  !! use the same ordering of the variables they depend on:
-  !! 1. Density
-  !! 2. Momentum
-  !! 3. Energy
-  !! Take care, not to use this routine, if the requested variable does not
-  !! have this ordering.
-  !!
-  !! \todo HK: Not quite sure, this effort is actually needed here, we can
-  !!           savely assume, the conservative variables to be properly ordered
-  !!           directly access them, instead of going over the dependencies.
-  !!           In my opinion it would be more useful, to have some integer
-  !!           variables declared in the euler system, that point to the correct
-  !!           positions after initialization and use those directly.
-  !!           The dependency stuff is more geared towards the usage of various
-  !!           derived quantities to find yet another one.
-  !!           However, the way it is right now, provides the most generic
-  !!           implementation and should work, no matter where the equation
-  !!           system description comes from. Thus, this small subroutine to
-  !!           avoid too much code duplication.
-  subroutine get_cvarPos(varsys, requVar, pos_Density, pos_Momentum, pos_Energy)
-    !> List of all variables in the available system.
-    type(tem_variable_type), intent(in) :: varsys(:)
-
-    !> Index of the requested variable.
-    integer, intent(in) :: requVar
-
-    !> Position of the density in the serialized variable system.
-    integer, optional, intent(out) :: pos_Density
-
-    !> Position of the 3 momentum components in the serialized variable system.
-    integer, optional, intent(out) :: pos_Momentum(3)
-
-    !> Position of the energy in the serialized variable system.
-    integer, optional, intent(out) :: pos_Energy
-
-    integer :: iDep
-
-    ! First dependency holds the density.
-    if (present(pos_Density)) then
-      iDep = varsys(requVar)%varDep(1)
-      pos_Density = varsys(iDep)%varPos(1)
-    end if
-
-    ! Second dependency holds the momentum.
-    if (present(pos_Momentum)) then
-      iDep = varsys(requVar)%varDep(2)
-      pos_Momentum = varsys(iDep)%varPos(1:3)
-    end if
-
-    ! Third dependency holds the energy.
-    if (present(pos_Energy)) then
-      iDep = varsys(requVar)%varDep(3)
-      pos_Energy = varsys(iDep)%varPos(1)
-    end if
-
-  end subroutine get_cvarPos
-
-
-  !> Small internal helper routine, to read Euler system parameters from a
-  !! Lua script.
-  !!
-  !! If the reading of parameters fails, the routine aborts the execution!
-  subroutine get_EulerParams(conf, euler)
-    !> Handle to the Lua script.
-    type(flu_State) :: conf
-
-    !> Description of Euler parameters, obtained from the Lua script.
-    type(Euler_type), intent(out) :: euler
-
-    integer :: eq_table
-
-    ! Try to open the equation table and abort if it fails.
-    call aot_table_open(L=conf, thandle=eq_table, key='equation')
-    if (eq_table.eq.0) then
-      call tem_write('ERROR in derive-routine: no equation table in ' &
-        &            // 'lua configuration file found, stopping...' )
-      call tem_abort()
-    end if
-
-    call init_euler(Euler, conf, eq_table)
-    call aot_table_close(L=conf, thandle=eq_table)
-
-  end subroutine get_EulerParams
-
-end module atl_eqn_euler_derive_module

source/atl_derived/atl_eqn_euler_module.f90

-! See copyright notice in the ../COPYRIGHT file.
-!> A module describing the Euler equation system.
-module atl_eqn_euler_module
-  use env_module, only: rk
-  use aotus_module, only: flu_State, aot_get_val
-
-  implicit none
-
-  private
-
-  public :: euler_type
-  public :: init_euler
-
-  !> The Euler equation properties are stored here
-  type Euler_type
-    real(kind=rk) :: isen_coef  !< isentropic coefficient
-    real(kind=rk) :: r          !< ideal gas constant
-    real(kind=rk) :: therm_cond !< thermal conductivity
-    real(kind=rk) :: cv         !< heat capacity
-  end type Euler_type
-
-contains
-
-  !> subroutine to initialize an equation of type euler equation
-  !! as defined in the configuration file
-  subroutine init_euler(euler, conf, eq_table)
-    !--------------------------------------------------------------------------!
-    !> Resulting description of the Euler equation parameters.
-    type(Euler_type), intent(out) :: euler
-
-    !> Handle to the configuration script, to load the parameters from.
-    type(flu_State)               :: conf
-
-    !> Handle to the table containing the description for the equation
-    !! system.
-    integer, intent(in)           :: eq_table
-    !--------------------------------------------------------------------------!
-    integer                       :: iError
-    !--------------------------------------------------------------------------!
-
-    !read the data from the equation table of the lua file
-    call aot_get_val(L = conf, thandle = eq_table, key = 'isen_coef', &
-      &              val = Euler%isen_coef, &
-      &              ErrCode = iError)
-    call aot_get_val(L = conf, thandle = eq_table, key = 'therm_cond', &
-      &              val = Euler%therm_cond, &
-      &              ErrCode = iError)
-    call aot_get_val(L = conf, thandle = eq_table, key = 'r', &
-      &              val = Euler%r, &
-      &              ErrCode = iError)
-    call aot_get_val(L = conf, thandle = eq_table, key = 'cv', &
-      &              val = Euler%cv, &
-      &              ErrCode = iError)
-
-  end subroutine init_euler
-
-end module atl_eqn_euler_module

source/atl_derived/atl_eqn_euler_var_module.f90

-! See copyright notice in the ../COPYRIGHT file.
-!> Module to configure the variables of the Euler equations.
-module atl_eqn_euler_var_module
-  use env_module,                  only: rk, newUnit
-  use tem_tools_module,            only: tem_write
-  use tem_var_system_module,       only: tem_var_system_type, &
-    &                                    tem_var_sys_init, tem_var_sys_mark, &
-    &                                    tem_var_append
-  use tem_dyn_array_module,        only: init, append
-  use atl_eqn_euler_derive_module, only: deriveVelocityFromConsVar, &
-    &                                    derivePressureFromConsVar, &
-    &                                    deriveTemperatureFromConsVar, &
-    &                                    deriveMachNumberFromConsVar, &
-    &                                    deriveMachVectorFromConsVar, &
-    &                                    deriveKineticEnergyFromConsVar
-  use atl_eqn_euler_module,        only: euler_type
-  use atl_eqn_source_module,       only: eval_source_spongeLayer
-  use atl_equation_module,         only: eqn_sourceMap_type
-
-
-  implicit none
-
-  private
-
-  public :: init_euler_vars
-  public :: append_euler_consVars
-  public :: append_euler_primVars
-  public :: append_euler_derivedVars
-  public :: init_euler_sourceTerms
-
-  contains
-
-  !> Init the variable system for Euler (inviscid) flow simulations.
-  !!
-  !! The variable system describes, which variables are to be used and how
-  !! they are organized in the memory.
-  !! The first few variables up to the sys_mark are those, describing the
-  !! state, and thus describe the output for regular restart files.
-  !! Here these are the conservative variables density, momentum and energy.
-  !! After the mark, there are additional values described that can be deduced
-  !! from the state variables.
-  subroutine init_euler_vars(varSys, hasPrimitiveVariables, primVar)
-    !--------------------------------------------------------------------------
-    !> The resulting variable system used in the Euler equations
-    type(tem_var_system_type), intent(out)   :: varSys
-
-    !> A flag indicating, that this system has primitive variables.
-    logical, intent(out)                     :: hasPrimitiveVariables
-
-    !> Indices to the primitive variables.
-    integer, allocatable, intent(out)        :: primVar(:)
-    !--------------------------------------------------------------------------
-    !--------------------------------------------------------------------------
-
-    ! Initialize variable system
-    call tem_var_sys_init( sys = varSys )
-    varSys%solSpecUnit = newUnit()
-
-    ! Append the conservative variables
-    call append_euler_consVars(varSys)
-
-    ! Set mark, to separate state variables from the rest
-    call tem_var_sys_mark(sys = varSys, systemName = 'euler_conservative')
-