text:
! SELF_Data.F90
!
! Copyright 2020-2022 Fluid Numerics LLC
! Author : Joseph Schoonover (joe@fluidnumerics.com)
! Support : support@fluidnumerics.com
!
! //////////////////////////////////////////////////////////////////////////////////////////////// !

MODULE SELF_Data

  USE SELF_Constants
  USE SELF_Lagrange
  USE SELF_Metadata
  USE FEQParse

  USE ISO_C_BINDING

  IMPLICIT NONE

#include "SELF_Macros.h"

  TYPE,PUBLIC :: SELF_DataObj
  !! The SELF_DataObj class is a base class for all data objects in SELF.
  !! A data object in SELF is a multidimensional array of data, represented
  !! on both host and device, that is associated with an interpolant, metadata,
  !! and (optionally) an equation string.
  !! Type extensions of the SELF_DataObj include scalars, vectors, and tensors
  !! in 1-D, 2-D, and 3-D using the storage patterns that are expected for
  !! derivative and interpolation operations defined in SELF_Lagrange.f90
  !! Additionally, each extended type has the necessary attributes to store
  !! information on element interiors and element boundaries, both of which
  !! are commonly used for spectral element solvers.

    INTEGER :: nVar
    INTEGER :: nElem
    TYPE(Lagrange),POINTER :: interp
    TYPE(Metadata),ALLOCATABLE :: meta(:)
    TYPE(EquationParser),ALLOCATABLE :: eqn(:)

    CONTAINS

    PROCEDURE,PUBLIC :: Init => Init_DataObj
    PROCEDURE,PUBLIC :: Free => Free_DataObj

    ! Procedures for setting metadata for
    PROCEDURE,PUBLIC :: SetName => SetName_DataObj
    PROCEDURE,PUBLIC :: SetDescription => SetDescription_DataObj
    PROCEDURE,PUBLIC :: SetUnits => SetUnits_DataObj
    GENERIC,PUBLIC :: SetEquation => SetEquation_DataObj
    PROCEDURE,PRIVATE :: SetEquation_DataObj

  END TYPE SELF_DataObj

! ---------------------- Scalars ---------------------- !
  TYPE,EXTENDS(SELF_DataObj),PUBLIC :: Scalar1D

    TYPE(hfReal_r3) :: interior
    TYPE(hfReal_r3) :: boundary
    TYPE(hfReal_r3) :: extBoundary
    TYPE(hfReal_r3) :: avgBoundary
    TYPE(hfReal_r3) :: jumpBoundary

  CONTAINS

    PROCEDURE,PUBLIC :: Init => Init_Scalar1D
    PROCEDURE,PUBLIC :: Free => Free_Scalar1D

    PROCEDURE,PUBLIC :: UpdateHost => UpdateHost_Scalar1D
    PROCEDURE,PUBLIC :: UpdateDevice => UpdateDevice_Scalar1D

    PROCEDURE,PUBLIC :: BoundaryInterp => BoundaryInterp_Scalar1D
    PROCEDURE,PUBLIC :: GridInterp => GridInterp_Scalar1D

    GENERIC,PUBLIC :: Derivative => Derivative_Scalar1D
    PROCEDURE,PRIVATE :: Derivative_Scalar1D

    PROCEDURE,PUBLIC :: AbsMaxInterior => AbsMaxInterior_Scalar1D
    PROCEDURE,PUBLIC :: AbsMaxBoundary => AbsMaxBoundary_Scalar1D

    GENERIC,PUBLIC :: ASSIGNMENT(=) => Equals_Scalar1D
    PROCEDURE,PRIVATE,PASS(SELFOut) :: Equals_Scalar1D

  END TYPE Scalar1D

  TYPE,EXTENDS(SELF_DataObj),PUBLIC :: Scalar2D

    TYPE(hfReal_r4) :: interior
    TYPE(hfReal_r4) :: boundary
    TYPE(hfReal_r4) :: extBoundary
    TYPE(hfReal_r4) :: avgBoundary
    TYPE(hfReal_r4) :: jumpBoundary

  CONTAINS

    PROCEDURE,PUBLIC :: Init => Init_Scalar2D
    PROCEDURE,PUBLIC :: Free => Free_Scalar2D
    PROCEDURE,PUBLIC :: UpdateHost => UpdateHost_Scalar2D
    PROCEDURE,PUBLIC :: UpdateDevice => UpdateDevice_Scalar2D
    PROCEDURE,PUBLIC :: BoundaryInterp => BoundaryInterp_Scalar2D
    PROCEDURE,PUBLIC :: GridInterp => GridInterp_Scalar2D

    GENERIC,PUBLIC :: Gradient => Gradient_Scalar2D
    PROCEDURE,PRIVATE :: Gradient_Scalar2D

    PROCEDURE,PUBLIC :: AbsMaxInterior => AbsMaxInterior_Scalar2D
    PROCEDURE,PUBLIC :: AbsMaxBoundary => AbsMaxBoundary_Scalar2D

    GENERIC,PUBLIC :: ASSIGNMENT(=) => Equals_Scalar2D
    PROCEDURE,PRIVATE,PASS(SELFOut) :: Equals_Scalar2D

  END TYPE Scalar2D

  TYPE,EXTENDS(SELF_DataObj),PUBLIC :: Scalar3D

    TYPE(hfReal_r5) :: interior
    TYPE(hfReal_r5) :: boundary
    TYPE(hfReal_r5) :: extBoundary
    TYPE(hfReal_r5) :: avgBoundary
    TYPE(hfReal_r5) :: jumpBoundary

  CONTAINS

    PROCEDURE,PUBLIC :: Init => Init_Scalar3D
    PROCEDURE,PUBLIC :: Free => Free_Scalar3D
    PROCEDURE,PUBLIC :: UpdateHost => UpdateHost_Scalar3D
    PROCEDURE,PUBLIC :: UpdateDevice => UpdateDevice_Scalar3D
    PROCEDURE,PUBLIC :: BoundaryInterp => BoundaryInterp_Scalar3D
    PROCEDURE,PUBLIC :: GridInterp => GridInterp_Scalar3D

    GENERIC,PUBLIC :: Gradient => Gradient_Scalar3D
    PROCEDURE,PRIVATE :: Gradient_Scalar3D

    PROCEDURE,PUBLIC :: AbsMaxInterior => AbsMaxInterior_Scalar3D
    PROCEDURE,PUBLIC :: AbsMaxBoundary => AbsMaxBoundary_Scalar3D

    GENERIC,PUBLIC :: ASSIGNMENT(=) => Equals_Scalar3D
    PROCEDURE,PRIVATE,PASS(SELFOut) :: Equals_Scalar3D

  END TYPE Scalar3D

! ---------------------- Vectors ---------------------- !

  TYPE,EXTENDS(SELF_DataObj),PUBLIC :: Vector2D

    TYPE(hfReal_r5) :: interior
    TYPE(hfReal_r5) :: boundary
    TYPE(hfReal_r5) :: extBoundary
    TYPE(hfReal_r4) :: boundaryNormal

  CONTAINS

    PROCEDURE,PUBLIC :: Init => Init_Vector2D
    PROCEDURE,PUBLIC :: Free => Free_Vector2D
    PROCEDURE,PUBLIC :: UpdateHost => UpdateHost_Vector2D
    PROCEDURE,PUBLIC :: UpdateDevice => UpdateDevice_Vector2D
    PROCEDURE,PUBLIC :: BoundaryInterp => BoundaryInterp_Vector2D
    PROCEDURE,PUBLIC :: GridInterp => GridInterp_Vector2D

    GENERIC,PUBLIC :: Gradient => Gradient_Vector2D
    PROCEDURE,PRIVATE :: Gradient_Vector2D

    GENERIC,PUBLIC :: Divergence => Divergence_Vector2D
    PROCEDURE,PRIVATE :: Divergence_Vector2D

    GENERIC,PUBLIC :: Curl => Curl_Vector2D
    PROCEDURE,PRIVATE :: Curl_Vector2D

    PROCEDURE,PUBLIC :: AbsMaxInterior => AbsMaxInterior_Vector2D
    PROCEDURE,PUBLIC :: AbsMaxBoundary => AbsMaxBoundary_Vector2D

    GENERIC,PUBLIC :: ASSIGNMENT(=) => Equals_Vector2D
    PROCEDURE,PRIVATE,PASS(SELFOut) :: Equals_Vector2D

    GENERIC,PUBLIC :: SetEquation => SetEquation_Vector2D
    PROCEDURE,PRIVATE :: SetEquation_Vector2D

  END TYPE Vector2D

  TYPE,EXTENDS(SELF_DataObj),PUBLIC :: Vector3D

    TYPE(hfReal_r6) :: interior
    TYPE(hfReal_r6) :: boundary
    TYPE(hfReal_r6) :: extBoundary
    TYPE(hfReal_r5) :: boundaryNormal

  CONTAINS

    PROCEDURE,PUBLIC :: Init => Init_Vector3D
    PROCEDURE,PUBLIC :: Free => Free_Vector3D
    PROCEDURE,PUBLIC :: UpdateHost => UpdateHost_Vector3D
    PROCEDURE,PUBLIC :: UpdateDevice => UpdateDevice_Vector3D
    PROCEDURE,PUBLIC :: BoundaryInterp => BoundaryInterp_Vector3D
    PROCEDURE,PUBLIC :: GridInterp => GridInterp_Vector3D

    GENERIC,PUBLIC :: Gradient => Gradient_Vector3D
    PROCEDURE,PRIVATE :: Gradient_Vector3D

    GENERIC,PUBLIC :: Divergence => Divergence_Vector3D
    PROCEDURE,PRIVATE :: Divergence_Vector3D

    GENERIC,PUBLIC :: Curl => Curl_Vector3D
    PROCEDURE,PRIVATE :: Curl_Vector3D

    PROCEDURE,PUBLIC :: AbsMaxInterior => AbsMaxInterior_Vector3D
    PROCEDURE,PUBLIC :: AbsMaxBoundary => AbsMaxBoundary_Vector3D

    GENERIC,PUBLIC :: ASSIGNMENT(=) => Equals_Vector3D
    PROCEDURE,PRIVATE,PASS(SELFOut) :: Equals_Vector3D

    GENERIC,PUBLIC :: SetEquation => SetEquation_Vector3D
    PROCEDURE,PRIVATE :: SetEquation_Vector3D

  END TYPE Vector3D

! ---------------------- Tensors ---------------------- !

  TYPE,EXTENDS(SELF_DataObj),PUBLIC :: Tensor2D

    TYPE(hfReal_r6) :: interior
    TYPE(hfReal_r6) :: boundary
    TYPE(hfReal_r6) :: extBoundary

  CONTAINS

    PROCEDURE,PUBLIC :: Init => Init_Tensor2D
    PROCEDURE,PUBLIC :: Free => Free_Tensor2D
    PROCEDURE,PUBLIC :: UpdateHost => UpdateHost_Tensor2D
    PROCEDURE,PUBLIC :: UpdateDevice => UpdateDevice_Tensor2D
    PROCEDURE,PUBLIC :: BoundaryInterp => BoundaryInterp_Tensor2D
    PROCEDURE,PUBLIC :: GridInterp => GridInterp_Tensor2D

    PROCEDURE,PUBLIC :: Determinant => Determinant_Tensor2D

    GENERIC,PUBLIC :: Divergence => Divergence_Tensor2D
    PROCEDURE,PRIVATE :: Divergence_Tensor2D

    PROCEDURE,PUBLIC :: AbsMaxInterior => AbsMaxInterior_Tensor2D
    PROCEDURE,PUBLIC :: AbsMaxBoundary => AbsMaxBoundary_Tensor2D

    GENERIC,PUBLIC :: ASSIGNMENT(=) => Equals_Tensor2D
    PROCEDURE,PRIVATE,PASS(SELFOut) :: Equals_Tensor2D

    GENERIC,PUBLIC :: SetEquation => SetEquation_Tensor2D
    PROCEDURE,PRIVATE :: SetEquation_Tensor2D

  END TYPE Tensor2D

  TYPE,EXTENDS(SELF_DataObj),PUBLIC :: Tensor3D

    TYPE(hfReal_r7) :: interior
    TYPE(hfReal_r7) :: boundary
    TYPE(hfReal_r7) :: extBoundary

  CONTAINS

    PROCEDURE,PUBLIC :: Init => Init_Tensor3D
    PROCEDURE,PUBLIC :: Free => Free_Tensor3D
    PROCEDURE,PUBLIC :: UpdateHost => UpdateHost_Tensor3D
    PROCEDURE,PUBLIC :: UpdateDevice => UpdateDevice_Tensor3D
    PROCEDURE,PUBLIC :: BoundaryInterp => BoundaryInterp_Tensor3D
    PROCEDURE,PUBLIC :: GridInterp => GridInterp_Tensor3D

    PROCEDURE,PUBLIC :: Determinant => Determinant_Tensor3D

    GENERIC,PUBLIC :: Divergence => Divergence_Tensor3D
    PROCEDURE,PRIVATE :: Divergence_Tensor3D

    PROCEDURE,PUBLIC :: AbsMaxInterior => AbsMaxInterior_Tensor3D
    PROCEDURE,PUBLIC :: AbsMaxBoundary => AbsMaxBoundary_Tensor3D

    GENERIC,PUBLIC :: ASSIGNMENT(=) => Equals_Tensor3D
    PROCEDURE,PRIVATE,PASS(SELFOut) :: Equals_Tensor3D

    GENERIC,PUBLIC :: SetEquation => SetEquation_Tensor3D
    PROCEDURE,PRIVATE :: SetEquation_Tensor3D

  END TYPE Tensor3D

  INTEGER,PARAMETER :: selfStrongForm = 0
  INTEGER,PARAMETER :: selfWeakDGForm = 1
  INTEGER,PARAMETER :: selfWeakCGForm = 2

  INTERFACE
    SUBROUTINE Determinant_Tensor2D_gpu_wrapper(tensor_dev,detTensor_dev,N,nVar,nEl) &
      bind(c,name="Determinant_Tensor2D_gpu_wrapper")
      USE iso_c_binding
      IMPLICIT NONE
      TYPE(c_ptr) :: tensor_dev,detTensor_dev
      INTEGER(C_INT),VALUE :: N,nVar,nEl
    END SUBROUTINE Determinant_Tensor2D_gpu_wrapper
  END INTERFACE

  INTERFACE
    SUBROUTINE Determinant_Tensor3D_gpu_wrapper(tensor_dev,detTensor_dev,N,nVar,nEl) &
      bind(c,name="Determinant_Tensor3D_gpu_wrapper")
      USE iso_c_binding
      IMPLICIT NONE
      TYPE(c_ptr) :: tensor_dev,detTensor_dev
      INTEGER(C_INT),VALUE :: N,nVar,nEl
    END SUBROUTINE Determinant_Tensor3D_gpu_wrapper
  END INTERFACE

CONTAINS

! -- DataObj -- !
  SUBROUTINE Init_DataObj(SELFStorage,interp,nVar,nElem)
    IMPLICIT NONE
    CLASS(SELF_DataObj),INTENT(out)  :: SELFStorage
    TYPE(Lagrange),INTENT(in),TARGET :: interp
    INTEGER,INTENT(in) :: nVar
    INTEGER,INTENT(in) :: nElem

    SELFStorage % nElem = nElem
    SELFStorage % nVar = nVar
    SELFStorage % interp => interp
    ALLOCATE( SELFStorage % meta(1:nVar) )
    ALLOCATE( SELFStorage % eqn(1:nVar) )

  END SUBROUTINE Init_DataObj

  SUBROUTINE Free_DataObj(SELFStorage)
    IMPLICIT NONE
    CLASS(SELF_DataObj),INTENT(inout) :: SELFStorage

    SELFStorage % interp => NULL()
    DEALLOCATE( SELFStorage % meta )
    DEALLOCATE( SELFStorage % eqn )

  END SUBROUTINE Free_DataObj

  SUBROUTINE SetName_DataObj(SELFStorage,ivar,name)
    !! Set the name of the `ivar-th` variable
    IMPLICIT NONE
    CLASS(SELF_DataObj),INTENT(inout) :: SELFStorage
    INTEGER,INTENT(in) :: ivar
    CHARACTER(*),INTENT(in) :: name

    CALL SELFStorage % meta(ivar) % SetName(name)

  END SUBROUTINE SetName_DataObj

  SUBROUTINE SetDescription_DataObj(SELFStorage,ivar,description)
    !! Set the description of the `ivar-th` variable
    IMPLICIT NONE
    CLASS(SELF_DataObj),INTENT(inout) :: SELFStorage
    INTEGER,INTENT(in) :: ivar
    CHARACTER(*),INTENT(in) :: description

    CALL SELFStorage % meta(ivar) % SetDescription(description)

  END SUBROUTINE SetDescription_DataObj

  SUBROUTINE SetUnits_DataObj(SELFStorage,ivar,units)
    !! Set the units of the `ivar-th` variable
    IMPLICIT NONE
    CLASS(SELF_DataObj),INTENT(inout) :: SELFStorage
    INTEGER,INTENT(in) :: ivar
    CHARACTER(*),INTENT(in) :: units

    CALL SELFStorage % meta(ivar) % SetUnits(units)

  END SUBROUTINE SetUnits_DataObj

  SUBROUTINE SetEquation_DataObj(SELFStorage,ivar,eqnChar)
    !! Sets the equation parser for the `ivar-th` variable
    IMPLICIT NONE
    CLASS(SELF_DataObj),INTENT(inout) :: SELFStorage
    INTEGER,INTENT(in) :: ivar
    CHARACTER(*),INTENT(in) :: eqnChar

    SELFStorage % eqn(ivar) = EquationParser( TRIM(eqnChar), &
                                              (/'x','y','z','t'/) )

  END SUBROUTINE SetEquation_DataObj

! -- Scalar1D -- !

  SUBROUTINE Init_Scalar1D(SELFStorage,interp,nVar,nElem)
    IMPLICIT NONE
    CLASS(Scalar1D),INTENT(out) :: SELFStorage
    TYPE(Lagrange),INTENT(in),TARGET :: interp
    INTEGER,INTENT(in) :: nVar
    INTEGER,INTENT(in) :: nElem

    SELFStorage % interp => interp
    SELFStorage % nVar = nVar
    SELFStorage % nElem = nElem

    CALL SELFStorage % interior % Alloc(loBound=(/0,1,1/), &
                                        upBound=(/interp % N,nVar,nElem/))

    CALL SELFStorage % boundary % Alloc(loBound=(/1,1,1/), &
                                        upBound=(/nVar,2,nElem/))

    CALL SELFStorage % extBoundary % Alloc(loBound=(/1,1,1/), &
                                           upBound=(/nVar,2,nElem/))

    CALL SELFStorage % avgBoundary % Alloc(loBound=(/1,1,1/), &
                                           upBound=(/nVar,2,nElem/))

    CALL SELFStorage % jumpBoundary % Alloc(loBound=(/1,1,1/), &
                                           upBound=(/nVar,2,nElem/))

    ALLOCATE( SELFStorage % meta(1:nVar) )
    ALLOCATE( SELFStorage % eqn(1:nVar) )

  END SUBROUTINE Init_Scalar1D

  SUBROUTINE Free_Scalar1D(SELFStorage)
    IMPLICIT NONE
    CLASS(Scalar1D),INTENT(inout) :: SELFStorage

    SELFStorage % interp => NULL()
    CALL SELFStorage % interior % Free()
    CALL SELFStorage % boundary % Free()
    CALL SELFStorage % extBoundary % Free()
    CALL SELFStorage % avgBoundary % Free()
    CALL SELFStorage % jumpBoundary % Free()
    DEALLOCATE( SELFStorage % meta )
    DEALLOCATE( SELFStorage % eqn )

  END SUBROUTINE Free_Scalar1D

  SUBROUTINE UpdateHost_Scalar1D(SELFStorage)
    IMPLICIT NONE
    CLASS(Scalar1D),INTENT(inout) :: SELFStorage

    CALL SELFStorage % interior % UpdateHost()
    CALL SELFStorage % boundary % UpdateHost()
    CALL SELFStorage % extBoundary % UpdateHost()
    CALL SELFStorage % avgBoundary % UpdateHost()
    CALL SELFStorage % jumpBoundary % UpdateHost()

  END SUBROUTINE UpdateHost_Scalar1D

  SUBROUTINE UpdateDevice_Scalar1D(SELFStorage)
    IMPLICIT NONE
    CLASS(Scalar1D),INTENT(inout) :: SELFStorage

    CALL SELFStorage % interior % UpdateDevice()
    CALL SELFStorage % boundary % UpdateDevice()
    CALL SELFStorage % extBoundary % UpdateDevice()
    CALL SELFStorage % avgBoundary % UpdateDevice()
    CALL SELFStorage % jumpBoundary % UpdateDevice()

  END SUBROUTINE UpdateDevice_Scalar1D

  SUBROUTINE BoundaryInterp_Scalar1D(SELFStorage,gpuAccel)
    IMPLICIT NONE
    CLASS(Scalar1D),INTENT(inout) :: SELFStorage
    LOGICAL,INTENT(in) :: gpuAccel

    IF (gpuAccel) THEN
      CALL SELFStorage % interp % ScalarBoundaryInterp_1D(SELFStorage % interior % deviceData, &
                                                          SELFStorage % boundary % deviceData, &
                                                          SELFStorage % nVar, &
                                                          SELFStorage % nElem)
    ELSE
      CALL SELFStorage % interp % ScalarBoundaryInterp_1D(SELFStorage % interior % hostData, &
                                                          SELFStorage % boundary % hostData, &
                                                          SELFStorage % nVar, &
                                                          SELFStorage % nElem)
    END IF

  END SUBROUTINE BoundaryInterp_Scalar1D

  SUBROUTINE GridInterp_Scalar1D(SELFStorage,SELFout,gpuAccel)
    IMPLICIT NONE
    CLASS(Scalar1D),INTENT(in) :: SELFStorage
    TYPE(Scalar1D),INTENT(inout) :: SELFOut
    LOGICAL,INTENT(in) :: gpuAccel

    IF (gpuAccel) THEN
      CALL SELFStorage % interp % ScalarGridInterp_1D(SELFStorage % interior % deviceData, &
                                                      SELFout % interior % deviceData, &
                                                      SELFStorage % nVar, &
                                                      SELFStorage % nElem)
    ELSE
      CALL SELFStorage % interp % ScalarGridInterp_1D(SELFStorage % interior % hostData, &
                                                      SELFout % interior % hostData, &
                                                      SELFStorage % nVar, &
                                                      SELFStorage % nElem)
    END IF

  END SUBROUTINE GridInterp_Scalar1D

  SUBROUTINE Derivative_Scalar1D(SELFStorage,SELFOut,gpuAccel)
    IMPLICIT NONE
    CLASS(Scalar1D),INTENT(in) :: SELFStorage
    TYPE(Scalar1D),INTENT(inout) :: SELFOut
    LOGICAL,INTENT(in) :: gpuAccel

    IF (gpuAccel) THEN
      CALL SELFStorage % interp % Derivative_1D(SELFStorage % interior % deviceData, &
                                                SELFout % interior % deviceData, &
                                                SELFStorage % nVar, &
                                                SELFStorage % nElem)
    ELSE
      CALL SELFStorage % interp % Derivative_1D(SELFStorage % interior % hostData, &
                                                SELFout % interior % hostData, &
                                                SELFStorage % nVar, &
                                                SELFStorage % nElem)
    END IF

  END SUBROUTINE Derivative_Scalar1D

  FUNCTION AbsMaxInterior_Scalar1D(scalar) RESULT(absMax)
    IMPLICIT NONE
    CLASS(Scalar1D) :: scalar
    REAL(prec) :: absMax(1:scalar % nVar)
    ! Local
    INTEGER :: iEl,iVar,i

    absMax = 0.0_prec
    DO iEl = 1,scalar % nElem
      DO iVar = 1,scalar % nVar
        DO i = 0,scalar % interp % N
          absMax(iVar) = MAX(ABS(scalar % interior % hostData(i,iVar,iEl)),absMax(iVar))
        END DO
      END DO
    END DO

  END FUNCTION AbsMaxInterior_Scalar1D

  FUNCTION AbsMaxBoundary_Scalar1D(scalar) RESULT(absMax)
    IMPLICIT NONE
    CLASS(Scalar1D) :: scalar
    REAL(prec) :: absMax(1:scalar % nVar,1:2)
    ! Local
    INTEGER :: iEl,iVar,iSide

    absMax = 0.0_prec
    DO iEl = 1,scalar % nElem
      DO iSide = 1,2
        DO iVar = 1,scalar % nVar
          absMax(iVar,iSide) = MAX(ABS(scalar % boundary % hostData(iVar,iSide,iEl)),absMax(iVar,iSide))
        END DO
      END DO
    END DO

  END FUNCTION AbsMaxBoundary_Scalar1D

  SUBROUTINE Equals_Scalar1D(SELFOut,SELFin)
    IMPLICIT NONE
    CLASS(Scalar1D),INTENT(inout) :: SELFOut
    TYPE(Scalar1D),INTENT(in) :: SELFin

    SELFOut % interior % hostData = SELFin % interior % hostData
    SELFOut % boundary % hostData = SELFin % boundary % hostData

  END SUBROUTINE Equals_Scalar1D

! -- Scalar2D -- !

  SUBROUTINE Init_Scalar2D(SELFStorage,interp,nVar,nElem)
    IMPLICIT NONE
    CLASS(Scalar2D),INTENT(out) :: SELFStorage
    TYPE(Lagrange),INTENT(in),TARGET :: interp
    INTEGER,INTENT(in) :: nVar
    INTEGER,INTENT(in) :: nElem

    SELFStorage % interp => interp
    SELFStorage % nVar = nVar
    SELFStorage % nElem = nElem

    CALL SELFStorage % interior % Alloc(loBound=(/0,0,1,1/), &
                                        upBound=(/interp % N,interp % N,nVar,nElem/))

    CALL SELFStorage % boundary % Alloc(loBound=(/0,1,1,1/), &
                                        upBound=(/interp % N,nVar,4,nElem/))

    CALL SELFStorage % extBoundary % Alloc(loBound=(/0,1,1,1/), &
                                           upBound=(/interp % N,nVar,4,nElem/))

    CALL SELFStorage % avgBoundary % Alloc(loBound=(/0,1,1,1/), &
                                           upBound=(/interp % N,nVar,4,nElem/))

    CALL SELFStorage % jumpBoundary % Alloc(loBound=(/0,1,1,1/), &
                                           upBound=(/interp % N,nVar,4,nElem/))

    ALLOCATE( SELFStorage % meta(1:nVar) )
    ALLOCATE( SELFStorage % eqn(1:nVar) )

  END SUBROUTINE Init_Scalar2D

  SUBROUTINE Free_Scalar2D(SELFStorage)
    IMPLICIT NONE
    CLASS(Scalar2D),INTENT(inout) :: SELFStorage

    SELFStorage % nVar = 0
    SELFStorage % nElem = 0
    SELFStorage % interp => NULL()
    CALL SELFStorage % interior % Free()
    CALL SELFStorage % boundary % Free()
    CALL SELFStorage % extBoundary % Free()
    CALL SELFStorage % avgBoundary % Free()
    CALL SELFStorage % jumpBoundary % Free()

    DEALLOCATE( SELFStorage % meta )
    DEALLOCATE( SELFStorage % eqn )

  END SUBROUTINE Free_Scalar2D

  SUBROUTINE UpdateHost_Scalar2D(SELFStorage)
    IMPLICIT NONE
    CLASS(Scalar2D),INTENT(inout) :: SELFStorage

    CALL SELFStorage % interior % UpdateHost()
    CALL SELFStorage % boundary % UpdateHost()
    CALL SELFStorage % extBoundary % UpdateHost()
    CALL SELFStorage % avgBoundary % UpdateHost()
    CALL SELFStorage % jumpBoundary % UpdateHost()

  END SUBROUTINE UpdateHost_Scalar2D

  SUBROUTINE UpdateDevice_Scalar2D(SELFStorage)
    IMPLICIT NONE
    CLASS(Scalar2D),INTENT(inout) :: SELFStorage

    CALL SELFStorage % interior % UpdateDevice()
    CALL SELFStorage % boundary % UpdateDevice()
    CALL SELFStorage % extBoundary % UpdateDevice()
    CALL SELFStorage % avgBoundary % UpdateDevice()
    CALL SELFStorage % jumpBoundary % UpdateDevice()

  END SUBROUTINE UpdateDevice_Scalar2D

  SUBROUTINE BoundaryInterp_Scalar2D(SELFStorage,gpuAccel)
    IMPLICIT NONE
    CLASS(Scalar2D),INTENT(inout) :: SELFStorage
    LOGICAL,INTENT(in) :: gpuAccel

    IF (gpuAccel) THEN
      CALL SELFStorage % interp % ScalarBoundaryInterp_2D(SELFStorage % interior % deviceData, &
                                                          SELFStorage % boundary % deviceData, &
                                                          SELFStorage % nVar, &
                                                          SELFStorage % nElem)
    ELSE
      CALL SELFStorage % interp % ScalarBoundaryInterp_2D(SELFStorage % interior % hostData, &
                                                          SELFStorage % boundary % hostData, &
                                                          SELFStorage % nVar, &
                                                          SELFStorage % nElem)
    END IF

  END SUBROUTINE BoundaryInterp_Scalar2D

  SUBROUTINE GridInterp_Scalar2D(SELFStorage,SELFOut,gpuAccel)
    IMPLICIT NONE
    CLASS(Scalar2D),INTENT(in) :: SELFStorage
    TYPE(Scalar2D),INTENT(inout) :: SELFOut
    LOGICAL,INTENT(in) :: gpuAccel

    IF (gpuAccel) THEN
      CALL SELFStorage % interp % ScalarGridInterp_2D(SELFStorage % interior % deviceData, &
                                                      SELFout % interior % deviceData, &
                                                      SELFStorage % nVar, &
                                                      SELFStorage % nElem)
    ELSE
      CALL SELFStorage % interp % ScalarGridInterp_2D(SELFStorage % interior % hostData, &
                                                      SELFout % interior % hostData, &
                                                      SELFStorage % nVar, &
                                                      SELFStorage % nElem)
    END IF

  END SUBROUTINE GridInterp_Scalar2D

  SUBROUTINE Gradient_Scalar2D(SELFStorage,SELFOut,gpuAccel)
    IMPLICIT NONE
    CLASS(Scalar2D),INTENT(in) :: SELFStorage
    TYPE(Vector2D),INTENT(inout) :: SELFOut
    LOGICAL,INTENT(in) :: gpuAccel

    IF (gpuAccel) THEN
      CALL SELFStorage % interp % ScalarGradient_2D(SELFStorage % interior % deviceData, &
                                                    SELFout % interior % deviceData, &
                                                    SELFStorage % nVar, &
                                                    SELFStorage % nElem)
    ELSE
      CALL SELFStorage % interp % ScalarGradient_2D(SELFStorage % interior % hostData, &
                                                    SELFout % interior % hostData, &
                                                    SELFStorage % nVar, &
                                                    SELFStorage % nElem)
    END IF

  END SUBROUTINE Gradient_Scalar2D

  FUNCTION AbsMaxInterior_Scalar2D(scalar) RESULT(absMax)
    IMPLICIT NONE
    CLASS(Scalar2D) :: scalar
    REAL(prec) :: absMax(1:scalar % nVar)
    ! Local
    INTEGER :: iEl,iVar,i,j

    absMax = 0.0_prec
    DO iEl = 1,scalar % nElem
      DO iVar = 1,scalar % nVar
        DO j = 0,scalar % interp % N
          DO i = 0,scalar % interp % N
            absMax(iVar) = MAX(ABS(scalar % interior % hostData(i,j,iVar,iEl)),absMax(iVar))
          END DO
        END DO
      END DO
    END DO

  END FUNCTION AbsMaxInterior_Scalar2D

  FUNCTION AbsMaxBoundary_Scalar2D(scalar) RESULT(absMax)
    IMPLICIT NONE
    CLASS(Scalar2D) :: scalar
    REAL(prec) :: absMax(1:scalar % nVar,1:4)
    ! Local
    INTEGER :: iEl,iVar,i,iSide

    absMax = 0.0_prec
    DO iEl = 1,scalar % nElem
      DO iSide = 1,4
        DO iVar = 1,scalar % nVar
          DO i = 0,scalar % interp % N
            absMax(iVar,iSide) = MAX(ABS(scalar % boundary % hostData(i,iVar,iSide,iEl)),absMax(iVar,iSide))
          END DO
        END DO
      END DO
    END DO

  END FUNCTION AbsMaxBoundary_Scalar2D

  SUBROUTINE Equals_Scalar2D(SELFOut,SELFin)
    IMPLICIT NONE
    CLASS(Scalar2D),INTENT(inout) :: SELFOut
    TYPE(Scalar2D),INTENT(in) :: SELFin

    SELFOut % interior % hostData = SELFin % interior % hostData
    SELFOut % boundary % hostData = SELFin % boundary % hostData

  END SUBROUTINE Equals_Scalar2D

! -- Scalar3D -- !

  SUBROUTINE Init_Scalar3D(SELFStorage,interp,nVar,nElem)
    IMPLICIT NONE
    CLASS(Scalar3D),INTENT(out) :: SELFStorage
    TYPE(Lagrange),TARGET,INTENT(in) :: interp
    INTEGER,INTENT(in) :: nVar
    INTEGER,INTENT(in) :: nElem
    ! Local
    INTEGER :: N

    SELFStorage % interp => interp
    SELFStorage % nVar = nVar
    SELFStorage % nElem = nElem
    N = interp % N

    CALL SELFStorage % interior % Alloc(loBound=(/0,0,0,1,1/), &
                                        upBound=(/N,N,N,nVar,nElem/))

    CALL SELFStorage % boundary % Alloc(loBound=(/0,0,1,1,1/), &
                                        upBound=(/N,N,nVar,6,nElem/))

    CALL SELFStorage % extBoundary % Alloc(loBound=(/0,0,1,1,1/), &
                                           upBound=(/N,N,nVar,6,nElem/))

    CALL SELFStorage % avgBoundary % Alloc(loBound=(/0,0,1,1,1/), &
                                           upBound=(/N,N,nVar,6,nElem/))

    CALL SELFStorage % jumpBoundary % Alloc(loBound=(/0,0,1,1,1/), &
                                           upBound=(/N,N,nVar,6,nElem/))

    ALLOCATE( SELFStorage % meta(1:nVar) )
    ALLOCATE( SELFStorage % eqn(1:nVar) )

  END SUBROUTINE Init_Scalar3D

  SUBROUTINE Free_Scalar3D(SELFStorage)
    IMPLICIT NONE
    CLASS(Scalar3D),INTENT(inout) :: SELFStorage

    SELFStorage % nVar = 0
    SELFStorage % nElem = 0
    SELFStorage % interp => NULL()
    CALL SELFStorage % interior % Free()
    CALL SELFStorage % boundary % Free()
    CALL SELFStorage % extBoundary % Free()
    CALL SELFStorage % avgBoundary % Free()
    CALL SELFStorage % jumpBoundary % Free()

    DEALLOCATE( SELFStorage % meta )
    DEALLOCATE( SELFStorage % eqn )

  END SUBROUTINE Free_Scalar3D

  SUBROUTINE UpdateHost_Scalar3D(SELFStorage)
    IMPLICIT NONE
    CLASS(Scalar3D),INTENT(inout) :: SELFStorage

    CALL SELFStorage % interior % UpdateHost()
    CALL SELFStorage % boundary % UpdateHost()
    CALL SELFStorage % extBoundary % UpdateHost()
    CALL SELFStorage % avgBoundary % UpdateHost()
    CALL SELFStorage % jumpBoundary % UpdateHost()

  END SUBROUTINE UpdateHost_Scalar3D

  SUBROUTINE UpdateDevice_Scalar3D(SELFStorage)
    IMPLICIT NONE
    CLASS(Scalar3D),INTENT(inout) :: SELFStorage

    CALL SELFStorage % interior % UpdateDevice()
    CALL SELFStorage % boundary % UpdateDevice()
    CALL SELFStorage % extBoundary % UpdateDevice()
    CALL SELFStorage % avgBoundary % UpdateDevice()
    CALL SELFStorage % jumpBoundary % UpdateDevice()

  END SUBROUTINE UpdateDevice_Scalar3D

  SUBROUTINE BoundaryInterp_Scalar3D(SELFStorage,gpuAccel)
    IMPLICIT NONE
    CLASS(Scalar3D),INTENT(inout) :: SELFStorage
    LOGICAL,INTENT(in) :: gpuAccel

    IF (gpuAccel) THEN
      CALL SELFStorage % interp % ScalarBoundaryInterp_3D(SELFStorage % interior % deviceData, &
                                                          SELFStorage % boundary % deviceData, &
                                                          SELFStorage % nVar, &
                                                          SELFStorage % nElem)
    ELSE
      CALL SELFStorage % interp % ScalarBoundaryInterp_3D(SELFStorage % interior % hostData, &
                                                          SELFStorage % boundary % hostData, &
                                                          SELFStorage % nVar, &
                                                          SELFStorage % nElem)
    END IF

  END SUBROUTINE BoundaryInterp_Scalar3D

  SUBROUTINE GridInterp_Scalar3D(SELFStorage,SELFOut,gpuAccel)
    IMPLICIT NONE
    CLASS(Scalar3D),INTENT(in) :: SELFStorage
    TYPE(Scalar3D),INTENT(inout) :: SELFOut
    LOGICAL,INTENT(in) :: gpuAccel

    IF (gpuAccel) THEN
      CALL SELFStorage % interp % ScalarGridInterp_3D(SELFStorage % interior % deviceData, &
                                                      SELFout % interior % deviceData, &
                                                      SELFStorage % nVar, &
                                                      SELFStorage % nElem)
    ELSE
      CALL SELFStorage % interp % ScalarGridInterp_3D(SELFStorage % interior % hostData, &
                                                      SELFout % interior % hostData, &
                                                      SELFStorage % nVar, &
                                                      SELFStorage % nElem)
    END IF

  END SUBROUTINE GridInterp_Scalar3D

  SUBROUTINE Gradient_Scalar3D(SELFStorage,SELFOut,gpuAccel)
    IMPLICIT NONE
    CLASS(Scalar3D),INTENT(in) :: SELFStorage
    TYPE(Vector3D),INTENT(inout) :: SELFOut
    LOGICAL,INTENT(in) :: gpuAccel

    IF (gpuAccel) THEN
      CALL SELFStorage % interp % ScalarGradient_3D(SELFStorage % interior % deviceData, &
                                                    SELFout % interior % deviceData, &
                                                    SELFStorage % nVar, &
                                                    SELFStorage % nElem)
    ELSE
      CALL SELFStorage % interp % ScalarGradient_3D(SELFStorage % interior % hostData, &
                                                    SELFout % interior % hostData, &
                                                    SELFStorage % nVar, &
                                                    SELFStorage % nElem)
    END IF

  END SUBROUTINE Gradient_Scalar3D

  SUBROUTINE Equals_Scalar3D(SELFOut,SELFin)
    IMPLICIT NONE
    CLASS(Scalar3D),INTENT(inout) :: SELFOut
    TYPE(Scalar3D),INTENT(in) :: SELFin

    SELFOut % interior % hostData = SELFin % interior % hostData
    SELFOut % boundary % hostData = SELFin % boundary % hostData

  END SUBROUTINE Equals_Scalar3D

  FUNCTION AbsMaxInterior_Scalar3D(scalar) RESULT(absMax)
    IMPLICIT NONE
    CLASS(Scalar3D) :: scalar
    REAL(prec) :: absMax(1:scalar % nVar)
    ! Local
    INTEGER :: iEl,iVar,i,j,k

    absMax = 0.0_prec
    DO iEl = 1,scalar % nElem
      DO iVar = 1,scalar % nVar
        DO k = 0,scalar % interp % N
          DO j = 0,scalar % interp % N
            DO i = 0,scalar % interp % N
              absMax(iVar) = MAX(ABS(scalar % interior % hostData(i,j,k,iVar,iEl)),absMax(iVar))
            END DO
          END DO
        END DO
      END DO
    END DO

  END FUNCTION AbsMaxInterior_Scalar3D

  FUNCTION AbsMaxBoundary_Scalar3D(scalar) RESULT(absMax)
    IMPLICIT NONE
    CLASS(Scalar3D) :: scalar
    REAL(prec) :: absMax(1:scalar % nVar,1:6)
    ! Local
    INTEGER :: iEl,iVar,i,j,iSide

    absMax = 0.0_prec
    DO iEl = 1,scalar % nElem
      DO iSide = 1,6
        DO iVar = 1,scalar % nVar
          DO j = 0,scalar % interp % N
            DO i = 0,scalar % interp % N
              absMax(iVar,iSide) = MAX(ABS(scalar % boundary % hostData(i,j,iVar,iSide,iEl)),absMax(iVar,iSide))
            END DO
          END DO
        END DO
      END DO
    END DO

  END FUNCTION AbsMaxBoundary_Scalar3D

! -- Vector2D -- !

  SUBROUTINE Init_Vector2D(SELFStorage,interp,nVar,nElem)
    IMPLICIT NONE
    CLASS(Vector2D),INTENT(out) :: SELFStorage
    TYPE(Lagrange),TARGET,INTENT(in) :: interp
    INTEGER,INTENT(in) :: nVar
    INTEGER,INTENT(in) :: nElem
    ! Local
    INTEGER :: N

    SELFStorage % interp => interp
    SELFStorage % nVar = nVar
    SELFStorage % nElem = nElem
    N = interp % N

    CALL SELFStorage % interior % Alloc(loBound=(/1,0,0,1,1/), &
                                        upBound=(/2,N,N,nVar,nElem/))

    CALL SELFStorage % boundary % Alloc(loBound=(/1,0,1,1,1/), &
                                        upBound=(/2,N,nVar,4,nElem/))

    CALL SELFStorage % boundaryNormal % Alloc(loBound=(/0,1,1,1/), &
                                        upBound=(/N,nVar,4,nElem/))

    CALL SELFStorage % extBoundary % Alloc(loBound=(/1,0,1,1,1/), &
                                           upBound=(/2,N,nVar,4,nElem/))

    ALLOCATE( SELFStorage % meta(1:nVar) )
    ALLOCATE( SELFStorage % eqn(1:2*nVar) )

  END SUBROUTINE Init_Vector2D

  SUBROUTINE Free_Vector2D(SELFStorage)
    IMPLICIT NONE
    CLASS(Vector2D),INTENT(inout) :: SELFStorage

    SELFStorage % interp => NULL()
    SELFStorage % nVar = 0
    SELFStorage % nElem = 0
    CALL SELFStorage % interior % Free()
    CALL SELFStorage % boundary % Free()
    CALL SELFStorage % boundaryNormal % Free()
    CALL SELFStorage % extBoundary % Free()

    DEALLOCATE( SELFStorage % meta )
    DEALLOCATE( SELFStorage % eqn )

  END SUBROUTINE Free_Vector2D

  SUBROUTINE SetEquation_Vector2D(SELFStorage,idir,ivar,eqnChar)
    !! Sets the equation parser for the `idir` direction and `ivar-th` variable
    IMPLICIT NONE
    CLASS(Vector2D),INTENT(inout) :: SELFStorage
    INTEGER,INTENT(in) :: idir,ivar
    CHARACTER(*),INTENT(in) :: eqnChar

    SELFStorage % eqn(idir+2*(ivar-1)) = EquationParser( TRIM(eqnChar), &
                                              (/'x','y','z','t'/) )

  END SUBROUTINE SetEquation_Vector2D

  SUBROUTINE UpdateHost_Vector2D(SELFStorage)
    IMPLICIT NONE
    CLASS(Vector2D),INTENT(inout) :: SELFStorage

    CALL SELFStorage % interior % UpdateHost()
    CALL SELFStorage % boundary % UpdateHost()
    CALL SELFStorage % boundaryNormal % UpdateHost()
    CALL SELFStorage % extBoundary % UpdateHost()

  END SUBROUTINE UpdateHost_Vector2D

  SUBROUTINE UpdateDevice_Vector2D(SELFStorage)
    IMPLICIT NONE
    CLASS(Vector2D),INTENT(inout) :: SELFStorage

    CALL SELFStorage % interior % UpdateDevice()
    CALL SELFStorage % boundary % UpdateDevice()
    CALL SELFStorage % boundaryNormal % UpdateDevice()
    CALL SELFStorage % extBoundary % UpdateDevice()

  END SUBROUTINE UpdateDevice_Vector2D

  SUBROUTINE BoundaryInterp_Vector2D(SELFStorage,gpuAccel)
    IMPLICIT NONE
    CLASS(Vector2D),INTENT(inout) :: SELFStorage
    LOGICAL,INTENT(in) :: gpuAccel

    IF (gpuAccel) THEN
      CALL SELFStorage % interp % VectorBoundaryInterp_2D(SELFStorage % interior % deviceData, &
                                                          SELFStorage % boundary % deviceData, &
                                                          SELFStorage % nVar, &
                                                          SELFStorage % nElem)
    ELSE
      CALL SELFStorage % interp % VectorBoundaryInterp_2D(SELFStorage % interior % hostData, &
                                                          SELFStorage % boundary % hostData, &
                                                          SELFStorage % nVar, &
                                                          SELFStorage % nElem)
    END IF

  END SUBROUTINE BoundaryInterp_Vector2D

  SUBROUTINE GridInterp_Vector2D(SELFStorage,SELFOut,gpuAccel)
    IMPLICIT NONE
    CLASS(Vector2D),INTENT(in) :: SELFStorage
    TYPE(Vector2D),INTENT(inout) :: SELFOut
    LOGICAL,INTENT(in) :: gpuAccel

    IF (gpuAccel) THEN
      CALL SELFStorage % interp % VectorGridInterp_2D(SELFStorage % interior % deviceData, &
                                                      SELFout % interior % deviceData, &
                                                      SELFStorage % nVar, &
                                                      SELFStorage % nElem)
    ELSE
      CALL SELFStorage % interp % VectorGridInterp_2D(SELFStorage % interior % hostData, &
                                                      SELFout % interior % hostData, &
                                                      SELFStorage % nVar, &
                                                      SELFStorage % nElem)
    END IF

  END SUBROUTINE GridInterp_Vector2D

  SUBROUTINE Gradient_Vector2D(SELFStorage,SELFOut,gpuAccel)
    IMPLICIT NONE
    CLASS(Vector2D),INTENT(in) :: SELFStorage
    TYPE(Tensor2D),INTENT(inout) :: SELFOut
    LOGICAL,INTENT(in) :: gpuAccel

    IF (gpuAccel) THEN
      CALL SELFStorage % interp % VectorGradient_2D(SELFStorage % interior % deviceData, &
                                                    SELFout % interior % deviceData, &
                                                    SELFStorage % nVar, &
                                                    SELFStorage % nElem)
    ELSE
      CALL SELFStorage % interp % VectorGradient_2D(SELFStorage % interior % hostData, &
                                                    SELFout % interior % hostData, &
                                                    SELFStorage % nVar, &
                                                    SELFStorage % nElem)
    END IF

  END SUBROUTINE Gradient_Vector2D

  SUBROUTINE Divergence_Vector2D(SELFStorage,SELFOut,dForm,gpuAccel)
    IMPLICIT NONE
    CLASS(Vector2D),INTENT(in) :: SELFStorage
    TYPE(Scalar2D),INTENT(inout) :: SELFOut
    INTEGER,INTENT(in) :: dForm
    LOGICAL,INTENT(in) :: gpuAccel

    IF (dForm == selfWeakDGForm) THEN

      IF (gpuAccel) THEN
        CALL SELFStorage % interp % VectorDGDivergence_2D(SELFStorage % interior % deviceData, &
                                                          SELFStorage % boundaryNormal % deviceData, &
                                                          SELFout % interior % deviceData, &
                                                          SELFStorage % nVar, &
                                                          SELFStorage % nElem)
      ELSE
        CALL SELFStorage % interp % VectorDGDivergence_2D(SELFStorage % interior % hostData, &
                                                          SELFStorage % boundaryNormal % hostData, &
                                                          SELFout % interior % hostData, &
                                                          SELFStorage % nVar, &
                                                          SELFStorage % nElem)
      END IF

    ELSE IF (dForm == selfStrongForm) THEN

      IF (gpuAccel) THEN
        CALL SELFStorage % interp % VectorDivergence_2D(SELFStorage % interior % deviceData, &
                                                        SELFout % interior % deviceData, &
                                                        SELFStorage % nVar, &
                                                        SELFStorage % nElem)
      ELSE
        CALL SELFStorage % interp % VectorDivergence_2D(SELFStorage % interior % hostData, &
                                                        SELFout % interior % hostData, &
                                                        SELFStorage % nVar, &
                                                        SELFStorage % nElem)
      END IF

    END IF

  END SUBROUTINE Divergence_Vector2D

  SUBROUTINE Curl_Vector2D(SELFStorage,SELFOut,gpuAccel)
    IMPLICIT NONE
    CLASS(Vector2D),INTENT(in) :: SELFStorage
    TYPE(Scalar2D),INTENT(inout) :: SELFOut
    LOGICAL,INTENT(in) :: gpuAccel

    IF (gpuAccel) THEN
      CALL SELFStorage % interp % VectorCurl_2D(SELFStorage % interior % deviceData, &
                                                SELFout % interior % deviceData, &
                                                SELFStorage % nVar, &
                                                SELFStorage % nElem)
    ELSE
      CALL SELFStorage % interp % VectorCurl_2D(SELFStorage % interior % hostData, &
                                                SELFout % interior % hostData, &
                                                SELFStorage % nVar, &
                                                SELFStorage % nElem)
    END IF

  END SUBROUTINE Curl_Vector2D

  SUBROUTINE Equals_Vector2D(SELFOut,SELFin)
    IMPLICIT NONE
    CLASS(Vector2D),INTENT(inout) :: SELFOut
    TYPE(Vector2D),INTENT(in) :: SELFin

    SELFOut % interior % hostData = SELFin % interior % hostData
    SELFOut % boundary % hostData = SELFin % boundary % hostData

  END SUBROUTINE Equals_Vector2D

  FUNCTION AbsMaxInterior_Vector2D(vector) RESULT(absMax)
    IMPLICIT NONE
    CLASS(Vector2D) :: vector
    REAL(prec) :: absMax(1:vector % nVar)
    ! Local
    INTEGER :: iEl,iVar,i,j,iDir

    absMax = 0.0_prec
    DO iEl = 1,vector % nElem
      DO iVar = 1,vector % nVar
        DO j = 0,vector % interp % N
          DO i = 0,vector % interp % N
            DO iDir = 1,2
              absMax(iVar) = MAX(ABS(vector % interior % hostData(iDir,i,j,iVar,iEl)),absMax(iVar))
            END DO
          END DO
        END DO
      END DO
    END DO

  END FUNCTION AbsMaxInterior_Vector2D

  FUNCTION AbsMaxBoundary_Vector2D(vector) RESULT(absMax)
    IMPLICIT NONE
    CLASS(Vector2D) :: vector
    REAL(prec) :: absMax(1:vector % nVar,1:4)
    ! Local
    INTEGER :: iEl,iVar,i,iDir,iSide

    absMax = 0.0_prec
    DO iEl = 1,vector % nElem
      DO iSide = 1,4
        DO iVar = 1,vector % nVar
          DO i = 0,vector % interp % N
            DO iDir = 1,2
              absMax(iVar,iSide) = MAX(ABS(vector % boundary % hostData(iDir,i,iVar,iSide,iEl)),absMax(iVar,iSide))
            END DO
          END DO
        END DO
      END DO
    END DO

  END FUNCTION AbsMaxBoundary_Vector2D

! -- Vector3D -- !

  SUBROUTINE Init_Vector3D(SELFStorage,interp,nVar,nElem)
    IMPLICIT NONE
    CLASS(Vector3D),INTENT(out) :: SELFStorage
    TYPE(Lagrange),TARGET,INTENT(in) :: interp
    INTEGER,INTENT(in) :: nVar
    INTEGER,INTENT(in) :: nElem
    ! Local
    INTEGER :: N

    SELFStorage % interp => interp
    SELFStorage % nVar = nVar
    SELFStorage % nElem = nElem
    N = interp % N

    CALL SELFStorage % interior % Alloc(loBound=(/1,0,0,0,1,1/), &
                                        upBound=(/3,N,N,N,nVar,nElem/))

    CALL SELFStorage % boundary % Alloc(loBound=(/1,0,0,1,1,1/), &
                                        upBound=(/3,N,N,nVar,6,nElem/))

    CALL SELFStorage % boundaryNormal % Alloc(loBound=(/0,0,1,1,1/), &
                                        upBound=(/N,N,nVar,6,nElem/))

    CALL SELFStorage % extBoundary % Alloc(loBound=(/1,0,0,1,1,1/), &
                                           upBound=(/3,N,N,nVar,6,nElem/))

    ALLOCATE( SELFStorage % meta(1:nVar) )
    ALLOCATE( SELFStorage % eqn(1:3*nVar) )

  END SUBROUTINE Init_Vector3D

  SUBROUTINE Free_Vector3D(SELFStorage)
    IMPLICIT NONE
    CLASS(Vector3D),INTENT(inout) :: SELFStorage

    SELFStorage % interp => NULL()
    SELFStorage % nVar = 0
    SELFStorage % nElem = 0
    CALL SELFStorage % interior % Free()
    CALL SELFStorage % boundary % Free()
    CALL SELFStorage % boundaryNormal % Free()
    CALL SELFStorage % extBoundary % Free()

    DEALLOCATE( SELFStorage % meta )
    DEALLOCATE( SELFStorage % eqn )

  END SUBROUTINE Free_Vector3D

  SUBROUTINE SetEquation_Vector3D(SELFStorage,idir,ivar,eqnChar)
    !! Sets the equation parser for the `idir` direction and `ivar-th` variable
    IMPLICIT NONE
    CLASS(Vector3D),INTENT(inout) :: SELFStorage
    INTEGER,INTENT(in) :: idir,ivar
    CHARACTER(*),INTENT(in) :: eqnChar

    SELFStorage % eqn(idir+3*(ivar-1)) = EquationParser( TRIM(eqnChar), &
                                              (/'x','y','z','t'/) )

  END SUBROUTINE SetEquation_Vector3D

  SUBROUTINE UpdateHost_Vector3D(SELFStorage)
    IMPLICIT NONE
    CLASS(Vector3D),INTENT(inout) :: SELFStorage

    CALL SELFStorage % interior % UpdateHost()
    CALL SELFStorage % boundary % UpdateHost()
    CALL SELFStorage % boundaryNormal % UpdateHost()
    CALL SELFStorage % extBoundary % UpdateHost()

  END SUBROUTINE UpdateHost_Vector3D

  SUBROUTINE UpdateDevice_Vector3D(SELFStorage)
    IMPLICIT NONE
    CLASS(Vector3D),INTENT(inout) :: SELFStorage

    CALL SELFStorage % interior % UpdateDevice()
    CALL SELFStorage % boundary % UpdateDevice()
    CALL SELFStorage % boundaryNormal % UpdateDevice()
    CALL SELFStorage % extBoundary % UpdateDevice()

  END SUBROUTINE UpdateDevice_Vector3D

  SUBROUTINE BoundaryInterp_Vector3D(SELFStorage,gpuAccel)
    IMPLICIT NONE
    CLASS(Vector3D),INTENT(inout) :: SELFStorage
    LOGICAL,INTENT(in) :: gpuAccel

    IF (gpuAccel) THEN
      CALL SELFStorage % interp % VectorBoundaryInterp_3D(SELFStorage % interior % deviceData, &
                                                          SELFStorage % boundary % deviceData, &
                                                          SELFStorage % nVar, &
                                                          SELFStorage % nElem)
    ELSE
      CALL SELFStorage % interp % VectorBoundaryInterp_3D(SELFStorage % interior % hostData, &
                                                          SELFStorage % boundary % hostData, &
                                                          SELFStorage % nVar, &
                                                          SELFStorage % nElem)
    END IF

  END SUBROUTINE BoundaryInterp_Vector3D

  SUBROUTINE GridInterp_Vector3D(SELFStorage,SELFOut,gpuAccel)
    IMPLICIT NONE
    CLASS(Vector3D),INTENT(in) :: SELFStorage
    TYPE(Vector3D),INTENT(inout) :: SELFOut
    LOGICAL,INTENT(in) :: gpuAccel

    IF (gpuAccel) THEN
      CALL SELFStorage % interp % VectorGridInterp_3D(SELFStorage % interior % deviceData, &
                                                      SELFout % interior % deviceData, &
                                                      SELFStorage % nVar, &
                                                      SELFStorage % nElem)
    ELSE
      CALL SELFStorage % interp % VectorGridInterp_3D(SELFStorage % interior % hostData, &
                                                      SELFout % interior % hostData, &
                                                      SELFStorage % nVar, &
                                                      SELFStorage % nElem)
    END IF

  END SUBROUTINE GridInterp_Vector3D

  SUBROUTINE Gradient_Vector3D(SELFStorage,SELFOut,gpuAccel)
    IMPLICIT NONE
    CLASS(Vector3D),INTENT(in) :: SELFStorage
    TYPE(Tensor3D),INTENT(inout) :: SELFOut
    LOGICAL,INTENT(in) :: gpuAccel

    IF (gpuAccel) THEN
      CALL SELFStorage % interp % VectorGradient_3D(SELFStorage % interior % deviceData, &
                                                    SELFout % interior % deviceData, &
                                                    SELFStorage % nVar, &
                                                    SELFStorage % nElem)
    ELSE
      CALL SELFStorage % interp % VectorGradient_3D(SELFStorage % interior % hostData, &
                                                    SELFout % interior % hostData, &
                                                    SELFStorage % nVar, &
                                                    SELFStorage % nElem)
    END IF

  END SUBROUTINE Gradient_Vector3D

  SUBROUTINE Divergence_Vector3D(SELFStorage,SELFOut,gpuAccel)
    IMPLICIT NONE
    CLASS(Vector3D),INTENT(in) :: SELFStorage
    TYPE(Scalar3D),INTENT(inout) :: SELFOut
    LOGICAL,INTENT(in) :: gpuAccel

    IF (gpuAccel) THEN
      CALL SELFStorage % interp % VectorDivergence_3D(SELFStorage % interior % deviceData, &
                                                      SELFout % interior % deviceData, &
                                                      SELFStorage % nVar, &
                                                      SELFStorage % nElem)
    ELSE
      CALL SELFStorage % interp % VectorDivergence_3D(SELFStorage % interior % hostData, &
                                                      SELFout % interior % hostData, &
                                                      SELFStorage % nVar, &
                                                      SELFStorage % nElem)
    END IF

  END SUBROUTINE Divergence_Vector3D

  SUBROUTINE Curl_Vector3D(SELFStorage,SELFOut,gpuAccel)
    IMPLICIT NONE
    CLASS(Vector3D),INTENT(in) :: SELFStorage
    TYPE(Vector3D),INTENT(inout) :: SELFOut
    LOGICAL,INTENT(in) :: gpuAccel

    IF (gpuAccel) THEN
      CALL SELFStorage % interp % VectorCurl_3D(SELFStorage % interior % deviceData, &
                                                SELFout % interior % deviceData, &
                                                SELFStorage % nVar, &
                                                SELFStorage % nElem)
    ELSE
      CALL SELFStorage % interp % VectorCurl_3D(SELFStorage % interior % hostData, &
                                                SELFout % interior % hostData, &
                                                SELFStorage % nVar, &
                                                SELFStorage % nElem)
    END IF

  END SUBROUTINE Curl_Vector3D

  FUNCTION AbsMaxInterior_Vector3D(vector) RESULT(absMax)
    IMPLICIT NONE
    CLASS(Vector3D) :: vector
    REAL(prec) :: absMax(1:vector % nVar)
    ! Local
    INTEGER :: iEl,iVar,i,j,k,iDir

    absMax = 0.0_prec
    DO iEl = 1,vector % nElem
      DO iVar = 1,vector % nVar
        DO k = 0,vector % interp % N
          DO j = 0,vector % interp % N
            DO i = 0,vector % interp % N
              DO iDir = 1,3
                absMax(iVar) = MAX(ABS(vector % interior % hostData(iDir,i,j,k,iVar,iEl)),absMax(iVar))
              END DO
            END DO
          END DO
        END DO
      END DO
    END DO

  END FUNCTION AbsMaxInterior_Vector3D

  FUNCTION AbsMaxBoundary_Vector3D(vector) RESULT(absMax)
    IMPLICIT NONE
    CLASS(Vector3D) :: vector
    REAL(prec) :: absMax(1:vector % nVar,1:6)
    ! Local
    INTEGER :: iEl,iVar,i,j,iSide,iDir

    absMax = 0.0_prec
    DO iEl = 1,vector % nElem
      DO iSide = 1,6
        DO iVar = 1,vector % nVar
          DO j = 0,vector % interp % N
            DO i = 0,vector % interp % N
              DO iDir = 1,3
                absMax(iVar,iSide) = MAX(ABS(vector % boundary % hostData(iDir,i,j,iVar,iSide,iEl)),absMax(iVar,iSide))
              END DO
            END DO
          END DO
        END DO
      END DO
    END DO

  END FUNCTION AbsMaxBoundary_Vector3D

  SUBROUTINE Equals_Vector3D(SELFOut,SELFin)
    IMPLICIT NONE
    CLASS(Vector3D),INTENT(inout) :: SELFOut
    TYPE(Vector3D),INTENT(in) :: SELFin

    SELFOut % interior % hostData = SELFin % interior % hostData
    SELFOut % boundary % hostData = SELFin % boundary % hostData

  END SUBROUTINE Equals_Vector3D

! -- Tensor2D -- !

  SUBROUTINE Init_Tensor2D(SELFStorage,interp,nVar,nElem)
    IMPLICIT NONE
    CLASS(Tensor2D),INTENT(out) :: SELFStorage
    TYPE(Lagrange),TARGET,INTENT(in) :: interp
    INTEGER,INTENT(in) :: nVar
    INTEGER,INTENT(in) :: nElem
    ! Local
    INTEGER :: N

    SELFStorage % interp => interp
    SELFStorage % nVar = nVar
    SELFStorage % nElem = nElem
    N = interp % N

    CALL SELFStorage % interior % Alloc(loBound=(/1,1,0,0,1,1/), &
                                        upBound=(/2,2,N,N,nVar,nElem/))

    CALL SELFStorage % boundary % Alloc(loBound=(/1,1,0,1,1,1/), &
                                        upBound=(/2,2,N,nVar,4,nElem/))

    CALL SELFStorage % extBoundary % Alloc(loBound=(/1,1,0,1,1,1/), &
                                           upBound=(/2,2,N,nVar,4,nElem/))

    ALLOCATE( SELFStorage % meta(1:nVar) )
    ALLOCATE( SELFStorage % eqn(1:4*nVar) )

  END SUBROUTINE Init_Tensor2D

  SUBROUTINE Free_Tensor2D(SELFStorage)
    IMPLICIT NONE
    CLASS(Tensor2D),INTENT(inout) :: SELFStorage

    SELFStorage % interp => NULL()
    SELFStorage % nVar = 0
    SELFStorage % nElem = 0
    CALL SELFStorage % interior % Free()
    CALL SELFStorage % boundary % Free()
    CALL SELFStorage % extBoundary % Free()

    DEALLOCATE( SELFStorage % meta )
    DEALLOCATE( SELFStorage % eqn )

  END SUBROUTINE Free_Tensor2D

  SUBROUTINE SetEquation_Tensor2D(SELFStorage,row,col,ivar,eqnChar)
    !! Sets the equation parser for row, col  of the ivar-th tensor
    IMPLICIT NONE
    CLASS(Tensor2D),INTENT(inout) :: SELFStorage
    INTEGER,INTENT(in) :: row,col,ivar
    CHARACTER(*),INTENT(in) :: eqnChar
    ! Local
    INTEGER :: ind

    ind = row+2*(col-1+2*(ivar-1))
    SELFStorage % eqn(ind) = EquationParser( TRIM(eqnChar), &
                                              (/'x','y','z','t'/) )

  END SUBROUTINE SetEquation_Tensor2D

  SUBROUTINE UpdateHost_Tensor2D(SELFStorage)
    IMPLICIT NONE
    CLASS(Tensor2D),INTENT(inout) :: SELFStorage

    CALL SELFStorage % interior % UpdateHost()
    CALL SELFStorage % boundary % UpdateHost()
    CALL SELFStorage % extBoundary % UpdateHost()

  END SUBROUTINE UpdateHost_Tensor2D

  SUBROUTINE UpdateDevice_Tensor2D(SELFStorage)
    IMPLICIT NONE
    CLASS(Tensor2D),INTENT(inout) :: SELFStorage

    CALL SELFStorage % interior % UpdateDevice()
    CALL SELFStorage % boundary % UpdateDevice()
    CALL SELFStorage % extBoundary % UpdateDevice()

  END SUBROUTINE UpdateDevice_Tensor2D

  SUBROUTINE BoundaryInterp_Tensor2D(SELFStorage,gpuAccel)
    IMPLICIT NONE
    CLASS(Tensor2D),INTENT(inout) :: SELFStorage
    LOGICAL,INTENT(in) :: gpuAccel

    IF (gpuAccel) THEN
      CALL SELFStorage % interp % TensorBoundaryInterp_2D(SELFStorage % interior % deviceData, &
                                                          SELFStorage % boundary % deviceData, &
                                                          SELFStorage % nVar, &
                                                          SELFStorage % nElem)
    ELSE
      CALL SELFStorage % interp % TensorBoundaryInterp_2D(SELFStorage % interior % hostData, &
                                                          SELFStorage % boundary % hostData, &
                                                          SELFStorage % nVar, &
                                                          SELFStorage % nElem)
    END IF

  END SUBROUTINE BoundaryInterp_Tensor2D

  SUBROUTINE GridInterp_Tensor2D(SELFStorage,SELFOut,gpuAccel)
    IMPLICIT NONE
    CLASS(Tensor2D),INTENT(in) :: SELFStorage
    TYPE(Tensor2D),INTENT(inout) :: SELFOut
    LOGICAL,INTENT(in) :: gpuAccel

    IF (gpuAccel) THEN
      CALL SELFStorage % interp % TensorGridInterp_2D(SELFStorage % interior % deviceData, &
                                                      SELFout % interior % deviceData, &
                                                      SELFStorage % nVar, &
                                                      SELFStorage % nElem)
    ELSE
      CALL SELFStorage % interp % TensorGridInterp_2D(SELFStorage % interior % hostData, &
                                                      SELFout % interior % hostData, &
                                                      SELFStorage % nVar, &
                                                      SELFStorage % nElem)
    END IF

  END SUBROUTINE GridInterp_Tensor2D

  SUBROUTINE Divergence_Tensor2D(SELFStorage,SELFOut,gpuAccel)
    IMPLICIT NONE
    CLASS(Tensor2D),INTENT(in) :: SELFStorage
    TYPE(Vector2D),INTENT(inout) :: SELFOut
    LOGICAL,INTENT(in) :: gpuAccel

    IF (gpuAccel) THEN
      CALL SELFStorage % interp % TensorDivergence_2D(SELFStorage % interior % deviceData, &
                                                      SELFout % interior % deviceData, &
                                                      SELFStorage % nVar, &
                                                      SELFStorage % nElem)
    ELSE
      CALL SELFStorage % interp % TensorDivergence_2D(SELFStorage % interior % hostData, &
                                                      SELFout % interior % hostData, &
                                                      SELFStorage % nVar, &
                                                      SELFStorage % nElem)
    END IF

  END SUBROUTINE Divergence_Tensor2D

  SUBROUTINE Determinant_Tensor2D(SELFStorage,SELFout,gpuAccel)
#undef __FUNC__
#define __FUNC__ "Determinant_Tensor2D"
    IMPLICIT NONE
    CLASS(Tensor2D),INTENT(in) :: SELFStorage
    TYPE(Scalar2D),INTENT(inout) :: SELFOut
    LOGICAL,INTENT(in) :: gpuAccel
    ! Local
    INTEGER :: iEl,iVar,i,j

    IF (gpuAccel) THEN

      CALL Determinant_Tensor2D_gpu_wrapper(SELFStorage % interior % deviceData, &
                                            SELFOut % interior % deviceData, &
                                            SELFStorage % interp % N, &
                                            SELFStorage % nVar, &
                                            SELFStorage % nElem)

    ELSE

      DO iEl = 1,SELFStorage % nElem
        DO iVar = 1,SELFStorage % nVar
          DO j = 0,SELFStorage % interp % N
            DO i = 0,SELFStorage % interp % N

              SELFOut % interior % hostData(i,j,iVar,iEl) = SELFStorage % interior % hostData(1,1,i,j,iVar,iEl)* &
                                                            SELFStorage % interior % hostData(2,2,i,j,iVar,iEl) - &
                                                            SELFStorage % interior % hostData(1,2,i,j,iVar,iEl)* &
                                                            SELFStorage % interior % hostData(2,1,i,j,iVar,iEl)

            END DO
          END DO
        END DO
      END DO

    END IF

  END SUBROUTINE Determinant_Tensor2D

  FUNCTION AbsMaxInterior_Tensor2D(tensor) RESULT(absMax)
    IMPLICIT NONE
    CLASS(Tensor2D) :: tensor
    REAL(prec) :: absMax(1:tensor % nVar)
    ! Local
    INTEGER :: iEl,iVar,i,j,row,col

    absMax = 0.0_prec
    DO iEl = 1,tensor % nElem
      DO iVar = 1,tensor % nVar
        DO j = 0,tensor % interp % N
          DO i = 0,tensor % interp % N
            DO col = 1,2
              DO row = 1,2
                absMax(iVar) = MAX(ABS(tensor % interior % hostData(row,col,i,j,iVar,iEl)),absMax(iVar))
              END DO
            END DO
          END DO
        END DO
      END DO
    END DO

  END FUNCTION AbsMaxInterior_Tensor2D

  FUNCTION AbsMaxBoundary_Tensor2D(tensor) RESULT(absMax)
    IMPLICIT NONE
    CLASS(Tensor2D) :: tensor
    REAL(prec) :: absMax(1:tensor % nVar,1:4)
    ! Local
    INTEGER :: iEl,iVar,i,iSide,row,col

    absMax = 0.0_prec
    DO iEl = 1,tensor % nElem
      DO iSide = 1,4
        DO iVar = 1,tensor % nVar
          DO i = 0,tensor % interp % N
            DO col = 1,2
              DO row = 1,2
                absMax(iVar,iSide) = MAX(ABS(tensor % boundary % hostData(row,col,i,iVar,iSide,iEl)),absMax(iVar,iSide))
              END DO
            END DO
          END DO
        END DO
      END DO
    END DO

  END FUNCTION AbsMaxBoundary_Tensor2D

  SUBROUTINE Equals_Tensor2D(SELFOut,SELFin)
    IMPLICIT NONE
    CLASS(Tensor2D),INTENT(inout) :: SELFOut
    TYPE(Tensor2D),INTENT(in) :: SELFin

    SELFOut % interior % hostData = SELFin % interior % hostData
    SELFOut % boundary % hostData = SELFin % boundary % hostData

  END SUBROUTINE Equals_Tensor2D

! -- Tensor3D -- !

  SUBROUTINE Init_Tensor3D(SELFStorage,interp,nVar,nElem)
    IMPLICIT NONE
    CLASS(Tensor3D),INTENT(out) :: SELFStorage
    TYPE(Lagrange),TARGET,INTENT(in) :: interp
    INTEGER,INTENT(in) :: nVar
    INTEGER,INTENT(in) :: nElem
    ! Local
    INTEGER :: N

    SELFStorage % interp => interp
    SELFStorage % nVar = nVar
    SELFStorage % nElem = nElem
    N = interp % N

    CALL SELFStorage % interior % Alloc(loBound=(/1,1,0,0,0,1,1/), &
                                        upBound=(/3,3,N,N,N,nVar,nElem/))

    CALL SELFStorage % boundary % Alloc(loBound=(/1,1,0,0,1,1,1/), &
                                        upBound=(/3,3,N,N,nVar,6,nElem/))

    CALL SELFStorage % extBoundary % Alloc(loBound=(/1,1,0,0,1,1,1/), &
                                           upBound=(/3,3,N,N,nVar,6,nElem/))

    ALLOCATE( SELFStorage % meta(1:nVar) )
    ALLOCATE( SELFStorage % eqn(1:9*nVar) )

  END SUBROUTINE Init_Tensor3D

  SUBROUTINE Free_Tensor3D(SELFStorage)
    IMPLICIT NONE
    CLASS(Tensor3D),INTENT(inout) :: SELFStorage

    SELFStorage % interp => NULL()
    SELFStorage % nVar = 0
    SELFStorage % nElem = 0
    CALL SELFStorage % interior % Free()
    CALL SELFStorage % boundary % Free()
    CALL SELFStorage % extBoundary % Free()

    DEALLOCATE( SELFStorage % meta )
    DEALLOCATE( SELFStorage % eqn )

  END SUBROUTINE Free_Tensor3D

  SUBROUTINE SetEquation_Tensor3D(SELFStorage,row,col,ivar,eqnChar)
    !! Sets the equation parser for row, col  of the ivar-th tensor
    IMPLICIT NONE
    CLASS(Tensor3D),INTENT(inout) :: SELFStorage
    INTEGER,INTENT(in) :: row,col,ivar
    CHARACTER(*),INTENT(in) :: eqnChar
    ! Local
    INTEGER :: ind

    ind = row+3*(col-1+3*(ivar-1))
    SELFStorage % eqn(ind) = EquationParser( TRIM(eqnChar), &
                                              (/'x','y','z','t'/) )

  END SUBROUTINE SetEquation_Tensor3D

  SUBROUTINE UpdateHost_Tensor3D(SELFStorage)
    IMPLICIT NONE
    CLASS(Tensor3D),INTENT(inout) :: SELFStorage

    CALL SELFStorage % interior % UpdateHost()
    CALL SELFStorage % boundary % UpdateHost()
    CALL SELFStorage % extBoundary % UpdateHost()

  END SUBROUTINE UpdateHost_Tensor3D

  SUBROUTINE UpdateDevice_Tensor3D(SELFStorage)
    IMPLICIT NONE
    CLASS(Tensor3D),INTENT(inout) :: SELFStorage

    CALL SELFStorage % interior % UpdateDevice()
    CALL SELFStorage % boundary % UpdateDevice()
    CALL SELFStorage % extBoundary % UpdateHost()

  END SUBROUTINE UpdateDevice_Tensor3D

  SUBROUTINE BoundaryInterp_Tensor3D(SELFStorage,gpuAccel)
    IMPLICIT NONE
    CLASS(Tensor3D),INTENT(inout) :: SELFStorage
    LOGICAL,INTENT(in) :: gpuAccel

    IF (gpuAccel) THEN
      CALL SELFStorage % interp % TensorBoundaryInterp_3D(SELFStorage % interior % deviceData, &
                                                          SELFStorage % boundary % deviceData, &
                                                          SELFStorage % nVar, &
                                                          SELFStorage % nElem)
    ELSE
      CALL SELFStorage % interp % TensorBoundaryInterp_3D(SELFStorage % interior % hostData, &
                                                          SELFStorage % boundary % hostData, &
                                                          SELFStorage % nVar, &
                                                          SELFStorage % nElem)
    END IF

  END SUBROUTINE BoundaryInterp_Tensor3D

  SUBROUTINE GridInterp_Tensor3D(SELFStorage,SELFOut,gpuAccel)
    IMPLICIT NONE
    CLASS(Tensor3D),INTENT(in) :: SELFStorage
    TYPE(Tensor3D),INTENT(inout) :: SELFOut
    LOGICAL,INTENT(in) :: gpuAccel

    IF (gpuAccel) THEN
      CALL SELFStorage % interp % TensorGridInterp_3D(SELFStorage % interior % deviceData, &
                                                      SELFout % interior % deviceData, &
                                                      SELFStorage % nVar, &
                                                      SELFStorage % nElem)
    ELSE
      CALL SELFStorage % interp % TensorGridInterp_3D(SELFStorage % interior % hostData, &
                                                      SELFout % interior % hostData, &
                                                      SELFStorage % nVar, &
                                                      SELFStorage % nElem)
    END IF

  END SUBROUTINE GridInterp_Tensor3D

  SUBROUTINE Divergence_Tensor3D(SELFStorage,SELFOut,gpuAccel)
    IMPLICIT NONE
    CLASS(Tensor3D),INTENT(in) :: SELFStorage
    TYPE(Vector3D),INTENT(inout) :: SELFOut
    LOGICAL,INTENT(in) :: gpuAccel

    IF (gpuAccel) THEN
      CALL SELFStorage % interp % TensorDivergence_3D(SELFStorage % interior % deviceData, &
                                                      SELFout % interior % deviceData, &
                                                      SELFStorage % nVar, &
                                                      SELFStorage % nElem)
    ELSE
      CALL SELFStorage % interp % TensorDivergence_3D(SELFStorage % interior % hostData, &
                                                      SELFout % interior % hostData, &
                                                      SELFStorage % nVar, &
                                                      SELFStorage % nElem)
    END IF

  END SUBROUTINE Divergence_Tensor3D

  SUBROUTINE Determinant_Tensor3D(SELFStorage,SELFOut,gpuAccel)
#undef __FUNC__
#define __FUNC__ "Determinant_Tensor3D"
    IMPLICIT NONE
    CLASS(Tensor3D),INTENT(in) :: SELFStorage
    TYPE(Scalar3D),INTENT(inout) :: SELFOut
    LOGICAL,INTENT(in) :: gpuAccel
    ! Local
    INTEGER :: iEl,iVar,i,j,k

    IF (gpuAccel) THEN

      CALL Determinant_Tensor3D_gpu_wrapper(SELFStorage % interior % deviceData, &
                                            SELFOut % interior % deviceData, &
                                            SELFStorage % interp % N, &
                                            SELFStorage % nVar, &
                                            SELFStorage % nElem)

    ELSE

      DO iEl = 1,SELFStorage % nElem
        DO iVar = 1,SELFStorage % nVar
          DO k = 0,SELFStorage % interp % N
            DO j = 0,SELFStorage % interp % N
              DO i = 0,SELFStorage % interp % N

                SELFOut % interior % hostData(i,j,k,iVar,iEl) = &
                  SELFStorage % interior % hostData(1,1,i,j,k,iVar,iEl)* &
                  (SELFStorage % interior % hostData(2,2,i,j,k,iVar,iEl)* &
                   SELFStorage % interior % hostData(3,3,i,j,k,iVar,iEl) - &
                   SELFStorage % interior % hostData(2,3,i,j,k,iVar,iEl)* &
                   SELFStorage % interior % hostData(3,2,i,j,k,iVar,iEl)) - &
                  SELFStorage % interior % hostData(2,1,i,j,k,iVar,iEl)* &
                  (SELFStorage % interior % hostData(1,2,i,j,k,iVar,iEl)* &
                   SELFStorage % interior % hostData(3,3,i,j,k,iVar,iEl) - &
                   SELFStorage % interior % hostData(1,3,i,j,k,iVar,iEl)* &
                   SELFStorage % interior % hostData(3,2,i,j,k,iVar,iEl)) + &
                  SELFStorage % interior % hostData(3,1,i,j,k,iVar,iEl)* &
                  (SELFStorage % interior % hostData(1,2,i,j,k,iVar,iEl)* &
                   SELFStorage % interior % hostData(2,3,i,j,k,iVar,iEl) - &
                   SELFStorage % interior % hostData(1,3,i,j,k,iVar,iEl)* &
                   SELFStorage % interior % hostData(2,2,i,j,k,iVar,iEl))

              END DO
            END DO
          END DO
        END DO
      END DO

    END IF

  END SUBROUTINE Determinant_Tensor3D

  FUNCTION AbsMaxInterior_Tensor3D(tensor) RESULT(absMax)
    IMPLICIT NONE
    CLASS(Tensor3D) :: tensor
    REAL(prec) :: absMax(1:tensor % nVar)
    ! Local
    INTEGER :: iEl,iVar,i,j,k,row,col

    absMax = 0.0_prec
    DO iEl = 1,tensor % nElem
      DO iVar = 1,tensor % nVar
        DO k = 0,tensor % interp % N
          DO j = 0,tensor % interp % N
            DO i = 0,tensor % interp % N
              DO col = 1,3
                DO row = 1,3
                  absMax(iVar) = MAX(ABS(tensor % interior % hostData(row,col,i,j,k,iVar,iEl)),absMax(iVar))
                END DO
              END DO
            END DO
          END DO
        END DO
      END DO
    END DO

  END FUNCTION AbsMaxInterior_Tensor3D

  FUNCTION AbsMaxBoundary_Tensor3D(tensor) RESULT(absMax)
    IMPLICIT NONE
    CLASS(Tensor3D) :: tensor
    REAL(prec) :: absMax(1:tensor % nVar,1:6)
    ! Local
    INTEGER :: iEl,iVar,i,j,iSide,row,col

    absMax = 0.0_prec
    DO iEl = 1,tensor % nElem
      DO iSide = 1,6
        DO iVar = 1,tensor % nVar
          DO j = 0,tensor % interp % N
            DO i = 0,tensor % interp % N
              DO col = 1,3
                DO row = 1,3
              absMax(iVar,iSide) = MAX(ABS(tensor % boundary % hostData(row,col,i,j,iVar,iSide,iEl)),absMax(iVar,iSide))
                END DO
              END DO
            END DO
          END DO
        END DO
      END DO
    END DO

  END FUNCTION AbsMaxBoundary_Tensor3D

  SUBROUTINE Equals_Tensor3D(SELFOut,SELFin)
    IMPLICIT NONE
    CLASS(Tensor3D),INTENT(inout) :: SELFOut
    TYPE(Tensor3D),INTENT(in) :: SELFin

    SELFOut % interior % hostData = SELFin % interior % hostData
    SELFOut % boundary % hostData = SELFin % boundary % hostData

  END SUBROUTINE Equals_Tensor3D

END MODULE SELF_Data

max_stars_repo_name: FluidNumerics/self-fluids

max_issues_repo_name: FluidNumerics/self-fluids

max_forks_repo_name: FluidNumerics/SELF-Fluids


text:
There is an increasing interest for language immersion programs growing in the United States.  As people around the world become more easily connected by better technology, knowledge of a single language is not enough.  In the past five years, a number of private and public schools in the San Francisco Bay Area now offers ChineseMandarin immersion programs.

As someone whose ethnic roots are from China, but struggled to learn Chinese all my life, I am very interested in these ChineseImmersion Programs for my two young children.

I would love it if Davis had a Chinese Immersion program in one of its elementary schools. I have no idea how to go about getting one established. I just wanted to create this page to see if anybody else felt the same way and has any thoughts about this.

20101216 21:41:11 nbsp Definitely interested! Users/PeggyEnderle

20110418 12:21:33 nbsp If youre interested in an immersion program for your kids, this is the best one in town: Dees Wonderland Users/wryezome

20120308 15:55:28 nbsp Id love it too!  The thought never occurred to me until, in college, I met some American School graduates from Taiwan/HK.  These kids essentially attended (American) Englishimmersion schools while living in Asia.  I thought they had the best of both worlds... fluent in both English and Chinese, as well as culturally comfortable in both.   This site about efforts in the Bay Area might be useful: http://miparentscouncil.org/ Users/ChonT

20120308 18:37:23 nbsp Sacramento now has Chinese immersion and their waitlist is years long. The Sac Bee did an article about it last year. It would be awesome if this was more widely available. Unfortunately, theres shortage of teachers who would be qualified to teach and would need a credential for it. My husband is a teacher and looked into this since he speaks Chinese fluently and only a couple colleges in California offer a Chinese credential. Users/jsbmeb

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
// Generated Files
${PROJ_DIR}/lpif/rtl/lpif_txrx/lpif_txrx_x16_f2_slave_top.sv
${PROJ_DIR}/lpif/rtl/lpif_txrx/lpif_txrx_x16_f2_slave_concat.sv
${PROJ_DIR}/lpif/rtl/lpif_txrx/lpif_txrx_x16_f2_slave_name.sv

// Logic Link files
-f ${PROJ_DIR}/llink/rtl/llink.f

// Common Files
-f ${PROJ_DIR}/common/rtl/common.f

max_stars_repo_name: chipsalliance/aib-protocols

max_issues_repo_name: chipsalliance/aib-protocols

max_forks_repo_name: chipsalliance/aib-protocols


text:
! Set initial conditions for q.
subroutine qinit(meqn, mbc, mx, my, mz, xlower, ylower, zlower, dx, dy, dz, q, maux, aux)

    implicit none

    ! Input/Output
    integer, intent(in) :: meqn, mbc, mx, my, mz, maux
    real(kind=8), intent(in) :: xlower, ylower, zlower, dx, dy, dz
    real(kind=8), intent(inout) :: q(meqn, 1-mbc:mx+mbc, 1-mbc:my+mbc, 1-mbc:mz+mbc)
    real(kind=8), intent(in) :: aux(maux, 1-mbc:mx+mbc, 1-mbc:my+mbc, 1-mbc:mz+mbc)

    ! Local storage
    integer :: i, j, k
    real(kind=8) :: x, y, z, r, rho

    do i = 1, mx
        x = xlower + (i - 0.5d0) * dx
        do j = 1, my
            y = ylower + (j - 0.5d0) * dy
            do k = 1, mz
                z = zlower + (k - 0.5d0) * dz
                r = sqrt(x**2 + y**2 + z**2)
                rho = 1.d0 + 10.d0 * exp(-20.d0 * (r - 0.d0)**2)
                q(1, i, j, k) = rho
                q(2, i, j, k) = 0.d0
                q(3, i, j, k) = 0.d0
                q(4, i, j, k) = 0.d0
                q(5, i, j, k) = 0.d0
            end do
        end do
    end do
end subroutine qinit

max_stars_repo_name: rjleveque/amrclaw

max_issues_repo_name: clawpack/amrclaw

max_forks_repo_name: rjleveque/amrclaw


      SUBROUTINE DGATHR(N, X, INCX, INDEX, INCI, Y, INCY)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  DOUBLE
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:04:22  11/16/01
!...Switches:
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      INTEGER , INTENT(IN) :: N
      INTEGER  :: INCX
      INTEGER  :: INCI
      INTEGER , INTENT(IN) :: INCY
      INTEGER , INTENT(IN) :: INDEX(*)
      REAL(DOUBLE) , INTENT(IN) :: X(*)
      REAL(DOUBLE) , INTENT(OUT) :: Y(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: IY, I
!-----------------------------------------------

!*******
!     Assume incx=inci=1
!*******

      IY = 1
      Y(IY:(N-1)*INCY+IY:INCY) = X(INDEX(:N))
      RETURN
      END SUBROUTINE DGATHR


text:
PROGRAM mpi
    include 'mpif.h'

    integer process_Rank, size_Of_Cluster, ierror, tag

    call MPI_INIT(ierror)
    call MPI_COMM_SIZE(MPI_COMM_WORLD, size_Of_Cluster, ierror)
    call MPI_COMM_RANK(MPI_COMM_WORLD, process_Rank, ierror)

    print *, 'Hello World from process: ', process_Rank, 'of ', size_Of_Cluster

    call MPI_FINALIZE(ierror)
END PROGRAM

max_stars_repo_name: okhlybov/tclbuildtest

max_issues_repo_name: okhlybov/tclbuildtest

max_forks_repo_name: okhlybov/tclbuildtest


text:
c EQUISUBMU evalua la presion en equilibrio hidrostatico
c_______________________________________________________________
c equisubmu rutina que evalua la presion en equilibrio hidrostatico
c teniendo en cuenta el numero de electrones para el calculo de mu
c Luis Bellot y Basilio Ruiz 27/7/95
c Cambio del metodo de integracion Basilio Ruiz 3/9/96
c _______________________________________________________________
	subroutine equisubmu(ntau,tau1,t,pe,pg,z,ro)

	implicit real*4 (a-h,o-z)

	include 'PARAMETER'   !solo por kt, que era igual a 1000
	parameter (nex=28,cgases=83145100.)
        real*4 tau1(kt)
	real*4 tau(kt),t(*),pe(*),pg(*),kac,d2,x(kt),kappa(kt),taue(kt)
        real*4 z1(kt),z(*),ro(*),y(kt)
	real*4 mu,pcontorno
	integer*4 nmaxitera

	real*4 wgt,abu,ei1,ei2,pp(10),tsi,psi,d1(10)
        common/constantes/g,avog	!gravedad,n. avogadro/pmu
	common/mu/cth                   !esto esta deshabilitado (entra 1)
        common/precisoitera/precitera
        common/anguloheliocent/mu
        common/nmaxitera/nmaxitera
c        common/pinicial/pcontorno

        precitera=1.e-5
        nmaxitera=250

	g=mu*2.7414e+4		!gravedad cm/s^2 en fotosfera solar
       	avog=6.023e23
       	cth=1.
       	mu=1.

        do i=1,ntau
           tau(i)=tau1(i)+alog10(cth)  !this is the optical depth in
        end do                        !the VERTICAL direction for
	                              !hydrostatic equilib. computations
        tau0=1.
        imin=1
        do i=2,ntau
           if(abs(tau(i)).lt.tau0)then
              imin=i         !zero of z is at log tau (vertical) =0
              tau0=tau(i)
           end if
        end do

        do i=1,10
           d1(i)=0
        end do
c calculamos el peso molecular medio pmu
	pmusum=0.0
        asum=0.0
	do i=1,92
  	   ii=i
	   call neldatb(ii,0.,wgt,abu,ei1,ei2)
	   pmusum=pmusum+wgt*abu
           asum=asum+abu
	end do

	do i=1,ntau
           taue(i)=10.**tau(i)
        end do

	do i=1,ntau-1
            x(i+1)=g*(tau(i)-tau(i+1))*2.30259
        end do

c inicializamos
	tsi=t(ntau)
        psi=pe(ntau)
c        pg(ntau)=pcontorno
        call gasc(tsi,psi,pg(ntau),pp)
	call kappach(5.e-5,tsi,psi,pp,d1,d1,kac,d2,d2)
        pesomedio=pmusum/(asum+pp(8)) !peso molec. medio
        kappa(ntau)=kac*avog/pmusum
        ro(ntau)=pesomedio*pg(ntau)/tsi/cgases
        y(ntau)=taue(ntau)/kappa(ntau)/ro(ntau)

c integramos
        do i=ntau-1,1,-1
           pgpr=pg(i+1)+x(i+1)*taue(i+1)/kappa(i+1)
	   tsi=t(i)

           call pefrompg11(tsi,pgpr,psi)
           pe(i)=psi

           call gasc(tsi,psi,pgpr,pp)
	   call kappach(5.e-5,tsi,psi,pp,d1,d1,kac,d2,d2)
           pesomedio=pmusum/(asum+pp(8)) !peso molec. medio
           kappa(i)=kac*avog/pmusum

           pg(i)=pg(i+1)+x(i+1)*(taue(i+1)/kappa(i+1)+taue(i)/kappa(i))/2.

           nit=1
           dif=2.*abs((pg(i)-pgpr))/(pg(i)+pgpr)

           do while (nit.lt.20.and.dif.gt.1.e-5)
              nit=nit+1
              pgpr=pg(i)
              call pefrompg11(tsi,pgpr,psi)
              pe(i)=psi

              call gasc(tsi,psi,pgpr,pp)
	      call kappach(5.e-5,tsi,psi,pp,d1,d1,kac,d2,d2)
              pesomedio=pmusum/(asum+pp(8)) !peso molec. medio
              kappa(i)=kac*avog/pmusum
              pg(i)=pg(i+1)+x(i+1)*(taue(i+1)/kappa(i+1)+taue(i)/kappa(i))/2.
              dif=2.*abs((pg(i)-pgpr))/(pg(i)+pgpr)
            end do

            if(dif.gt.0.1)print*,'WARNING:
     & Hydrostatic equilibrium results in inaccurate electron pressures '
            call pefrompg11(tsi,pg(i),psi)
            pe(i)=psi
            ro(i)=pesomedio*pg(i)/tsi/cgases
            y(i)=taue(i)/kappa(i)/ro(i)
	end do
        z1(1)=0.
        do i=2,ntau  !this is the z scale along the vertical direction
           z1(i)=z1(i-1)+x(i)/g*(y(i-1)+y(i))/2.
        end do
        do i=1,ntau
           z(i)=(z1(i)-z1(imin))*1.e-5
        end do

c        do i=1,ntau   !we output geometrical heights ALONG the LOS
c          z(i)=z(i)/mu
c        enddo

       return
       end

c____________________________________________________________________
c PEFROMPG10 (incuida en EQUISUBMU) evalua la presion electonica  desde la
c temperatura y la presion gaseosa
c pefrompg10 evalua la presion electonica p correspondiente a t1 y pg
c
c	subroutine pefrompg10(t,pg,p)
c
c	implicit real*4 (a-h,o-z)
c        common/preciso/prec
c
c        dif=1.
c        n2=0
c        p1=p
c        do while (dif.gt.prec.and.n2.lt.50)
c	     n2=n2+1
c	     p=(p+p1)/2.
c             p1=p
c             call pe_pg10(t,p,pg)
c             dif=abs(abs( (p-p1)/(p+p1) )-1.0)
c         end do
c	return
c	end
c____________________________________________________________________
c PEFROMPG11 (incuida en EQUISUBMU) evalua la presion electonica  desde la
c temperatura y la presion gaseosa
c pefrompg11 evalua la presion electonica p correspondiente a t1 y pg

	subroutine pefrompg11(t,pg,p)

	implicit real*4 (a-h,o-z)
        common/precisoitera/prec
        common/nmaxitera/nmaxitera

        dif=1.
        n2=0

         p0=p
         call pe_pg10(t,p,pg)
         dif=abs( (p-p0)/p )
         if(dif .gt. 1)call inicia_pefrompgt(t,pg,p)

        p1=p
        do while (dif.gt.prec.and.n2.lt.nmaxitera)
	     n2=n2+1
	     p=(p+p1)/2.
             p1=p
             call pe_pg10(t,p,pg)
             dif=abs( (p-p1)/p )
         end do
c         print*,'pefrompg11 n2=',n2,dif,p,p1
	return
	end
c____________________________________________________________________
c PE_PG10 (incuida en EQUISUBMU) evalua la presion electonica  desde la
c la presion gaseosa y una estimacion de la presion electronica
c calcula la presion electronica a partir de la pg y de una estimacion de la pe

      subroutine pe_pg10(t,pe,pg)

      parameter (ncontr=28)
      real*4 cmol(91),alfai(ncontr),chi1(ncontr),chi2(ncontr),
     *u0(ncontr),u1(ncontr),u2(ncontr)

      real*4 du0,du1,du2,dcmol(91)

      if(t.lt.500)then
         print*,'pe_pg10: temperature < 500 K '
	 print*,'temperature = 500 K'
	 t=500.
      end if
      theta=5040./t
      g4=0.
      g5=0.
      if(pe.le.0)then
         pe=1.e-15
         g4=0.
         g5=0.
      else
         call molecb(theta,cmol,dcmol)
         do i=1,2
            call acota(cmol(i),-30.,30.)
         end do
         g4=pe*10.**(cmol(1))
         g5=pe*10.**(cmol(2))
      end if
c ahora calculo los niveles u0,u1,u2 y sus derivadas
      do 5 i=1,ncontr
         iii=i
5     	 call neldatb(iii,0.,weight,alfai(i),chi1(i),chi2(i))
6     do 4 i=1,ncontr
      	 iii=i
4     	 call nelfctb(iii,t,u0(iii),u1(iii),u2(iii),du0,du1,du2)

      g2=saha(theta,chi1(1),u0(1),u1(1),pe)   ! p(h+)/p(h)

      g3=saha(theta,0.754,1.,u0(1),pe)        ! p(h)/p(h-)

      call acota(g3,1.e-30,1.e30)
      g3=1.d0/g3                              ! p(h-)/p(h)

      g1=0.
      do 1 i=2,ncontr
        a=saha(theta,chi1(i),u0(i),u1(i),pe)
        b=saha(theta,chi2(i),u1(i),u2(i),pe)
	c=1.+a*(1.+b)
1	g1=g1+alfai(i)/c*a*(1.+2.*b)

        a=1.+g2+g3
        b=2.*(1.+g2/g5*g4)
        c=g5
        d=g2-g3
        e=g2/g5*g4

	call acotasig(a,1.e-15,1.e15)
	call acotasig(d,1.e-15,1.e15)

        c1=c*b**2+a*d*b-e*a**2

        c2=2.*a*e-d*b+a*b*g1

        c3=-(e+b*g1)

        f1=0.5*c2/c1

        f1=-f1+sign(1.,c1)*sqrt(f1**2-c3/c1)
        f5=(1.-a*f1)/b
        f4=e*f5
        f3=g3*f1
        f2=g2*f1
        fe=f2-f3+f4+g1

        call acota(fe,1.e-30,1.e30)
	phtot=pe/fe

        if(f5.gt.1.e-4) goto 2
          const6=g5/pe*f1**2
          const7=f2-f3+g1
	  do 3 i=1,5
      	       f5=phtot*const6
	       f4=e*f5
	       fe=const7+f4
	       phtot=pe/fe
3	  continue

2       pe=pg/(1.+(f1+f2+f3+f4+f5+0.1014)/fe)
        if(pe.le.0)pe=1.e-15

      return
      end
c ____________________________________________________________________________

c esta rutina interpola en ttau el modelo tau,t,pe y da la salida en
c ttau,tt,ppe por medio de un polinomio de grado ngrado

	subroutine interpolatp(ngrado,n,tau,t,pe,nnew,ttau,tt,ppe)

	real tau(*),t(*),pe(*),ttau(*),tt(*),ppe(*)
	real xa(11),ya(11)

c interpolaremos las presiones en logaritmos neperianos
	num=n
	do i=1,num
	   pe(i)=alog(pe(i))
	end do

c interpolamos
	n2=int(ngrado/2)

	do i=1,nnew
	   CALL LOCATE(TAU,NUM,TTAU(I),J)
	   n3=j-n2-1
           if(n3.lt.0)n3=0
           if(n3+ngrado+1.gt.num)n3=num-ngrado-1
	   do k=1,ngrado+1
	      xa(k)=tau(n3+k)
	   end do
	   do k=1,ngrado+1
	      ya(k)=t(n3+k)
	   end do
	   CALL POLINT(XA,YA,NGRADO+1,TTAU(I),TT(I),ERROR)

	   do k=1,ngrado+1
	      ya(k)=pe(n3+k)
	   end do
	   CALL POLINT(XA,YA,NGRADO+1,TTAU(I),ppe(I),ERROR)
        end do

	do i=1,nnew
           ppe(i)=exp(ppe(i))
        end do
        do i=1,n
           pe(i)=exp(pe(i))
	end do

	return
	end
c ____________________________________________________________________________

	subroutine acota(x,x0,x1)

	if(x.lt.x0)x=x0
        if(x.gt.x1)x=x1

        return
        end
c ____________________________________________________________________________

	subroutine acotasig(x,x0,x1)

        if(x.lt.0)then
           x=-x
           call acota(x,x0,x1)
           x=-x
         else
           call acota(x,x0,x1)
         end if
        return
        end
c ____________________________________________________________________________


text: #redirect Users/Robnum1


text:
Two for One Sessions Now Available

Touchwork Body Therapy is committed to helping you achieve your optimum health and happiness. Make Touchwork a regular part of your selfcare system.

RELAXATION MASSAGE on a regular basis can help you stay in touch with your inner self, relieve issues with muscle, bone, and tissue, and slow down your nervous system.

A blend of Swedish and Asian techniques, with mobility movements and stretching gives you amazing results.

◾Sessions are 75 minutes and customized to your needs.
◾Cost is $81.00.

REFLEXOLOGY focuses on the feet, hands, and ears, and provides for ultimate relaxation, improved circulation, a stimulated nervous system, and increased energy.

◾Sessions are 75 minutes and customized to your needs.
◾Cost is $81.00.

MIND/BODY THERAPY allows us to go deeper into the structures of your life by resolving blockages in both mind and body. Touchwork at the centers of energy concentration in the body can make a profound difference on the nonconscious levels and the autonomic nervous system.

These sessions integrate massage and energy work, while engaging the power of our minds to optimize healing.
◾Sessions are 75 minutes and customized to your needs.
◾Cost is $108.00.

We welcome your inquiries.


text:
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!
!!     AUTHOR: Kenneth Leiter (kenneth.leiter@arl.army.mil)
!!
!!     Use the Xdmf Fortran Bindings to write out a simple mesh consisting of
!!     two hexahedrons.  Link against the XdmfFortran library to compile.
!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

PROGRAM XdmfFortranExample
IMPLICIT NONE

  INTEGER*8 obj
  character*256 filename
  REAL*4 myPoints(3,3,4)
  INTEGER myConnections(8,2)
  REAL*8 myCellAttribute(2), myNodeAttribute(3,4)

  filename = 'my_output'//CHAR(0)

  myPoints(1,1,1) = 0
  myPoints(2,1,1) = 0
  myPoints(3,1,1) = 1
  myPoints(1,2,1) = 1
  myPoints(2,2,1) = 0
  myPoints(3,2,1) = 1
  myPoints(1,3,1) = 3
  myPoints(2,3,1) = 0
  myPoints(3,3,1) = 2
  myPoints(1,1,2) = 0
  myPoints(2,1,2) = 1
  myPoints(3,1,2) = 1
  myPoints(1,2,2) = 1
  myPoints(2,2,2) = 1
  myPoints(3,2,2) = 1
  myPoints(1,3,2) = 3
  myPoints(2,3,2) = 2
  myPoints(3,3,2) = 2
  myPoints(1,1,3) = 0
  myPoints(2,1,3) = 0
  myPoints(3,1,3) = -1
  myPoints(1,2,3) = 1
  myPoints(2,2,3) = 0
  myPoints(3,2,3) = -1
  myPoints(1,3,3) = 3
  myPoints(2,3,3) = 0
  myPoints(3,3,3) = -2
  myPoints(1,1,4) = 0
  myPoints(2,1,4) = 1
  myPoints(3,1,4) = -1
  myPoints(1,2,4) = 1
  myPoints(2,2,4) = 1
  myPoints(3,2,4) = -1
  myPoints(1,3,4) = 3
  myPoints(2,3,4) = 2
  myPoints(3,3,4) = -2

  myConnections(1,1) = 0
  myConnections(2,1) = 1
  myConnections(3,1) = 7
  myConnections(4,1) = 6
  myConnections(5,1) = 3
  myConnections(6,1) = 4
  myConnections(7,1) = 10
  myConnections(8,1) = 9
  myConnections(1,2) = 1
  myConnections(2,2) = 2
  myConnections(3,2) = 8
  myConnections(4,2) = 7
  myConnections(5,2) = 4
  myConnections(6,2) = 5
  myConnections(7,2) = 11
  myConnections(8,2) = 10

  myNodeAttribute(1,1) = 100
  myNodeAttribute(1,2) = 300
  myNodeAttribute(1,3) = 300
  myNodeAttribute(1,4) = 500
  myNodeAttribute(2,1) = 200
  myNodeAttribute(2,2) = 400
  myNodeAttribute(2,3) = 400
  myNodeAttribute(2,4) = 600
  myNodeAttribute(3,1) = 300
  myNodeAttribute(3,2) = 500
  myNodeAttribute(3,3) = 500
  myNodeAttribute(3,4) = 700

  myCellAttribute(1) = 100
  myCellAttribute(2) = 200

  CALL XDMFINIT(obj, filename)
!  CALL XDMFSETGRIDTOPOLOGY(obj, 'Hexahedron'//CHAR(0), 2, myConnections)\
!  CALL XDMFSETGEOMETRY(obj, foo, 36, XDMF_ARRAY_TYPE_FLOAT64, myPoints)
!  CALL XDMFADDGRIDATTRIBUTE(obj, 'NodeValues'//CHAR(0),'XDMF_FLOAT64_TYPE'//CHAR(0), 'NODE'//CHAR(0), &
!       'SCALAR'//CHAR(0), 12, myNodeAttribute)
!  CALL XDMFADDGRIDATTRIBUTE(obj, 'CellValues'//CHAR(0),'XDMF_FLOAT64_TYPE'//CHAR(0), 'CELL'//CHAR(0), &
!       'SCALAR'//CHAR(0), 2, myCellAttribute)
  CALL XDMFADDGRID(obj, 'TestGrid'//CHAR(0))
!  CALL XDMFWRITETOFILE(obj)
!  CALL XDMFSERIALIZE(obj)
  CALL XDMFCLOSE(obj)

END PROGRAM XdmfFortranExample

max_stars_repo_name: scottwedge/xdmf

max_issues_repo_name: scottwedge/xdmf

max_forks_repo_name: scottwedge/xdmf


text:
      subroutine d9upak (x, y, n)
c august 1980 portable edition.  w. fullerton, los alamos scientific lab
c
c unpack floating point number x so that x = y * 2.0**n, where
c 0.5 .le. abs(y) .lt. 1.0 .
c
      double precision x, y, absx
c
      absx = dabs(x)
      n = 0
      y = 0.0d0
      if (x.eq.0.0d0) return
c
 10   if (absx.ge.0.5d0) go to 20
      n = n - 1
      absx = absx*2.0d0
      go to 10
c
 20   if (absx.lt.1.0d0) go to 30
      n = n + 1
      absx = absx*0.5d0
      go to 20
c
 30   y = dsign (absx, x)
      return
c
      end

max_stars_repo_name: nilqed/spadlib

max_issues_repo_name: nilqed/spadlib

max_forks_repo_name: nilqed/spadlib


text:
!==============================================================================!
  subroutine Control_Mod_Turbulent_Schmidt_Number(val, verbose)
!------------------------------------------------------------------------------!
  implicit none
!---------------------------------[Arguments]----------------------------------!
  real              :: val
  logical, optional :: verbose
!==============================================================================!

  call Control_Mod_Read_Real_Item('TURBULENT_SCHMIDT_NUMBER', 0.9,  &
                                   val, verbose)

  end subroutine

max_issues_repo_name: EdinSmartLab/T-Flows

max_forks_repo_name: EdinSmartLab/T-Flows


text:
      SUBROUTINE TGLOC00(NTARGS,KTYPE,TARGNAME,NPARMS,PARMS,
     *     KSYS,KORIGIN,TARGVEC,
     *     KTIME,TSEC50, KPOS,SCPOS, KVEL,SCVEL, KATT,ATT, KTOP,TOPCEN,
     *     LUERR,IERR)
      IMPLICIT REAL*8 (A-H,O-Z)
C
C THIS ROUTINE IS PART OF THE TOSS TARGET LOCATION PACKAGE, TGLOC.
C IT CALLS APPROPRIATE ROUTINES TO COMPUTE TARGET LOCATIONS IN THE
C 'NATURAL' SYSTEM AND, IF NECESSARY, CALLS A TRANSFORMATION PACKAGE
C TO PUT THEM IN THE DESIRED OUTPUT SYSTEM.  ON ANY ONE CALL TO TGLOC00,
C TARGETS ARE ALL OF THE SAME TYPE.
C
C VAR     DIM    TYPE   I/O  DESCRIPTION
C ---     ---    ----   ---  -----------
C
C -->> ARGUMENTS HAVE THE SAME DESCRIPTION AS THOSE IN TGLOC.  NOTE,
C      HOWEVER THAT TGLOC CALLS TGLOC00 WITH INDEXED ARRAYS SO THAT,
C      FOR EXAMPLE, TARGNAME(1) IN THIS ROUTINE MAY CORRESPOND TO
C      TARGNAME(23) IN TGLOC, THE ROUTINE THAT CALLS THIS ONE.
C
C***********************************************************************
C
C BY C PETRUZZO, 8/85.
C     MODIFIED....
C
C***********************************************************************
C
C
C CALLING SEQUENCE ARRAYS:
C
      REAL*8 PARMS(NPARMS,1)              ! ACTUALLY, (NPARMS,NTARGS)
      REAL*8 TARGVEC(3,1)                 ! ACTUALLY, (3,NTARGS)
      CHARACTER*16 TARGNAME(1)            ! ACTUALLY, (NTARGS)
      INTEGER KTYPE(1)                    ! ACTUALLY, (NTARGS)
      REAL*8 SCPOS(3),SCVEL(3),ATT(4),TOPCEN(4)
C
C INTERNAL VARIABLES:
      REAL*8 IDENT(3,3)/1.D0,3*0.D0,1.D0,3*0.D0,1.D0/
      REAL*8 TEMP6(6),DUMTOP(5)/5*9999.D0/,TOPTEMP(5)
      INTEGER KSPEC1(3)/999,999,1/,KSPEC2(3)/999,999,1/
      INTEGER KOUNTERR/0/
      LOGICAL NEEDPACK
C
      IERR = 0
      IF(NTARGS.LE.0) GO TO 9999
C
      IBUG = 0
      LUBUG = 19
C
      IF(IBUG.NE.0) THEN
        WRITE(LUBUG,8502) NTARGS,KSYS,KORIGIN,NPARMS
 8502   FORMAT(' TGLOC00 DEBUG 8502. NTARGS=',I4,
     *            '   KSYS,KORIGIN=',2I3,'  NPARMS=',I3)
        NP = MIN(NPARMS,5)
        IF(NTARGS.GT.0) WRITE(LUBUG,8501)
     *      (I,KTYPE(I),TARGNAME(I),(PARMS(J,I),J=1,NP),I=1,NTARGS)
 8501   FORMAT(<NTARGS>('    ITARG=',I3,'  KTYPE=',I3,'  TARGNAME=',A/,
     *                  '       PARMS(1-..)=',<NP>G13.5/))
        END IF
C
C CHECK FOR ERRONEOUS CALL. ALL KTYPE(I)'S SHOULD BE THE SAME. IF NOT
C AN ERROR HAS CREPT IN.
C
      KTEST = KTYPE(1)
      DO I=1,NTARGS
        IF(KTYPE(I).NE.KTEST) THEN
          TYPE *,' TGLOC00. ERROR END. PROGRAMMER ERROR. SEE SOURCE.'
          STOP   ' TGLOC00. ERROR END. PROGRAMMER ERROR. SEE SOURCE.'
          END IF
        END DO
C
C
C ***************************************
C *  ERROR CHECKS FOR THIS TARGET TYPE  *
C ***************************************
C
      CALL TGLOC00A(KTYPE(1),NPARMS, KTIME,KPOS,KVEL,KATT,KTOP,ATT(1),
     *        KSYS,KORIGIN, LUERR,IERR1)
      IF(IERR1.NE.0) THEN
        IERR = 1
        DO I=1,NTARGS
          TARGVEC(1,I) = 'ERROR'
          TARGVEC(2,I) = 'ERROR'
          TARGVEC(3,I) = 'ERROR'
          END DO
        GO TO 9999
        END IF
C
C
C  *********************************************
C  *  GET TARGET LOCATION IN 'NATURAL' SYSTEM  *
C  *********************************************
C
      GO TO (100,200,300,400,500,600,700,800),KTYPE(1)
C
C---> MOVING CELESTIAL TARGET.
  100 CONTINUE
      CALL TGLOC01(NTARGS,TARGNAME,TSEC50,KSYSNAT,KCENNAT,TARGVEC,
     *       LUERR,IERR1)
      GO TO 5000
C
C---> EARTH-FIXED TARGET.
  200 CONTINUE
      CALL TGLOC02(NTARGS,TARGNAME,NPARMS,PARMS,KSYSNAT,KCENNAT,TARGVEC,
     *        LUERR,IERR1)
      GO TO 5000
C
C---> FIXED CELESTIAL TARGET.
  300 CONTINUE
      CALL TGLOC03(NTARGS,TARGNAME,NPARMS,PARMS,KSYSNAT,KCENNAT,TARGVEC,
     *        LUERR,IERR1)
      GO TO 5000
C
C---> FIXED TO LOCAL ORBITAL COORDINATES.
  400 CONTINUE
      CALL TGLOC04(NTARGS,TARGNAME,NPARMS,PARMS,KSYSNAT,KCENNAT,TARGVEC,
     *        LUERR,IERR1)
      GO TO 5000
C
C---> NOON OR MIDNIGHT ZENITH.
  500 CONTINUE
      CALL TGLOC05(NTARGS,KSYSNAT,KCENNAT,TARGNAME,
     *       TSEC50,SCPOS,SCVEL,TARGVEC,LUERR,IERR1)
      GO TO 5000
C
C---> HORIZON.
  600 CONTINUE
      CALL TGLOC06(NTARGS,TARGNAME,NPARMS,PARMS,SCPOS,
     *       KSYSNAT,KCENNAT,TARGVEC,LUERR,IERR1)
      GO TO 5000
C
C---> A SATELLITE.
  700 CONTINUE
      CALL TGLOC07(NTARGS,TARGNAME,NPARMS,PARMS,TSEC50,
     *         KSYSNAT,KCENNAT,TARGVEC,LUERR,IERR1)
      GO TO 5000
C
C---> NON-SPECIFIC.
  800 CONTINUE
      CALL TGLOC08(NTARGS,KSYSNAT,KCENNAT,TARGVEC)
      GO TO 5000
C
C
 5000 CONTINUE
C
      IF(IERR1.NE.0) THEN
        IERR=1
        IF(LUERR.GT.0) WRITE(LUERR,5001) KTYPE(1)
 5001   FORMAT(/,' TGLOC00. ERROR CONDITION. TARGET TYPE=',I4)
        GO TO 9999
        END IF
C
C
C  ******************************************
C  *  GET TARGET LOCATION IN OUTPUT SYSTEM  *
C  ******************************************
C
      IF(KSYSNAT.EQ.9999) KSYSNAT = KSYS    ! NATRL SYSTEM UNIMPORTANT
      IF(KCENNAT.EQ.9999) KCENNAT = KORIGIN ! NATRL ORIGIN UNIMPORTANT
C
C SET FLAG INDICATING NEED TO CALL CTPAK.
      NEEDPACK = KSYSNAT.NE.KSYS .OR. KCENNAT.NE.KORIGIN
C
      IF(IBUG.NE.0)  WRITE(LUBUG,8507)
     *   KTYPE(1),KSYSNAT,KCENNAT,KSYS,KORIGIN,NEEDPACK,
     *      (TARGNAME(J),(TARGVEC(I,J),I=1,3),J=1,NTARGS)
 8507 FORMAT(' TGLOC DEBUG 8507. KTYPE=',I3,
     *               ' NATL KSYSNAT,KCENNAT=',2I3/,
     *       '    CALLER KSYS,KORIGIN=',2I3,'   NEEDPACK=',L2/,
     *       '    NATURAL TARGVEC='/,(7X,A,3G16.8) )
C
C
      IF(NEEDPACK) THEN
        KSPEC1(1)=KSYSNAT
        KSPEC1(2)=KCENNAT
        KSPEC2(1)=KSYS
        KSPEC2(2)=KORIGIN
        IF(KTOP.NE.0) THEN
          TOPTEMP(1) = TOPCEN(1)
          TOPTEMP(2) = TOPCEN(2)
          TOPTEMP(3) = TOPCEN(3)
          TOPTEMP(4) = TOPCEN(4)
          IF(TOPTEMP(4).NE.0.D0) TOPTEMP(4) = CONST(59)
          TOPTEMP(5) = CONST(53)
          END IF
        CALL CTPAK(NTARGS,
     *    TARGVEC, KSPEC1,
     *     KTIME,TSEC50, KATT,ATT, KPOS,SCPOS, KVEL,SCVEL, KTOP,TOPTEMP,
     *    TARGVEC, KSPEC2,
     *     KTIME,TSEC50, KATT,ATT, KPOS,SCPOS, KVEL,SCVEL, KTOP,TOPTEMP,
     *    0,DUM,DUM, LUERR,IERR1)
        IF(IERR1.NE.0) THEN
          IERR=1
          DO I=1,NTARGS
            TARGVEC(1,I) = 'ERROR'
            TARGVEC(2,I) = 'ERROR'
            TARGVEC(3,I) = 'ERROR'
            END DO
          IF(LUERR.GT.0) WRITE(LUERR,9015)
 9015     FORMAT(' TGLOC00. ERROR RETURN FROM CTPACK.')
          END IF
        END IF
C
 9999 CONTINUE
C
      RETURN
      END

ext: for


text:
!!# PARAMETERS MODULE: <PAR_NumericStrings>
MODULE PAR_NumericStrings

!!## PURPOSE
!! Define the characters which can appear in
!! numeric strings.

!!## EXTERNAL KINDS
USE KND_IntrinsicTypes,ONLY: KIND_S !!((01-A-KND_IntrinsicTypes.f90))

!!## HISTORY
! 1.[waw] Author   = William A. Wieselquist
!         Modified = 104.2006
!         Contact  = william.wieselquist AT gmail.com

!!## DEFAULT IMPLICIT
IMPLICIT NONE

!!## DEFAULT ACCESS
PUBLIC

!!## EXTERNAL PARAMETERS
CHARACTER(LEN=*,KIND=KIND_S),PARAMETER :: decimal_     = "."
CHARACTER(LEN=*,KIND=KIND_S),PARAMETER :: digits_      = "0123456789"
CHARACTER(LEN=*,KIND=KIND_S),PARAMETER :: exponents_   = "Ee"
CHARACTER(LEN=*,KIND=KIND_S),PARAMETER :: signs_       = "+-"
CHARACTER(LEN=*,KIND=KIND_S),PARAMETER :: bool_        = "TF"
CHARACTER(LEN=*,KIND=KIND_S),PARAMETER :: orderedpair_ = "(,)"
CHARACTER(LEN=*,KIND=KIND_S),PARAMETER :: float_       = signs_//decimal_//digits_
CHARACTER(LEN=*,KIND=KIND_S),PARAMETER :: scientific_  = signs_//decimal_//digits_//exponents_
CHARACTER(LEN=*,KIND=KIND_S),PARAMETER :: cfloat_      = signs_//decimal_//digits_//orderedpair_
CHARACTER(LEN=*,KIND=KIND_S),PARAMETER :: cscientific_ = signs_//decimal_//digits_//exponents_//orderedpair_
CHARACTER(LEN=*,KIND=KIND_S),PARAMETER :: signedint_   = signs_//digits_
CHARACTER(LEN=*,KIND=KIND_S),PARAMETER :: unsignedint_ = digits_

END MODULE

max_stars_repo_name: wawiesel/tapack

max_issues_repo_name: wawiesel/tapack

max_forks_repo_name: wawiesel/tapack