text:
 MODULE InitializationModule

  USE KindModule, ONLY: &
    DP, Zero, Third, Half, &
    One, Three, &
    Pi, TwoPi, FourPi
  USE UnitsModule, ONLY: &
    Gram, Centimeter, &
    Kilometer, Kelvin, &
    BoltzmannConstant
  USE ProgramHeaderModule, ONLY: &
    ProgramName, &
    iX_B0, iX_E0, &
    iE_B0, iE_E0, &
    nDOF, nDOFX, nDOFE, &
    nZ, nNodes
  USE UtilitiesModule, ONLY: &
    Locate, &
    NodeNumberX, &
    Interpolate1D_Linear
  USE ReferenceElementModuleX, ONLY: &
    NodeNumberTableX
  USE ReferenceElementModule, ONLY: &
    NodeNumbersX, &
    NodeNumberTable, &
    OuterProduct1D3D
  USE MeshModule, ONLY: &
    MeshE, MeshX, &
    NodeCoordinate
  USE GeometryFieldsModule, ONLY: &
    uGF, iGF_Gm_dd_11, iGF_Gm_dd_22, iGF_Gm_dd_33
  USE FluidFieldsModule, ONLY: &
    uPF, iPF_D, iPF_E, iPF_Ne, &
    uAF, iAF_P, iAF_T, iAF_Ye, iAF_S, iAF_E, &
    iAF_Me, iAF_Mp, iAF_Mn, iAF_Xp, iAF_Xn, &
    iAF_Xa, iAF_Xh, iAF_Gm
  USE RadiationFieldsModule, ONLY: &
    nSpecies, &
    uPR, nPR, iPR_D, iPR_I1, iPR_I2, iPR_I3, &
    uCR, nCR, iCR_N, iCR_G1, iCR_G2, iCR_G3
  USE EquationOfStateModule, ONLY: &
    ApplyEquationOfState, &
    ComputeThermodynamicStates_Primitive
  USE EquationOfStateModule_TABLE, ONLY: &
    ApplyEquationOfState_TABLE, &
    ComputeThermodynamicStates_Primitive_TABLE
  USE TwoMoment_UtilitiesModule, ONLY: &
    ComputeConserved_TwoMoment
  USE ProgenitorModule, ONLY: &
    ProgenitorType1D, &
    ReadProgenitor1D

  IMPLICIT NONE
  PRIVATE

  INCLUDE 'mpif.h'

  PUBLIC :: InitializeFields
  PUBLIC :: ComputeError

CONTAINS

  SUBROUTINE InitializeFields &
             ( Direction_Option, SigmaA_Option, SigmaS_Option, Profile_Option )

    CHARACTER(LEN=*), INTENT(in), OPTIONAL :: &
      Direction_Option
    CHARACTER(len=*), INTENT(in), OPTIONAL :: &
      Profile_Option
    REAL(DP),         INTENT(in), OPTIONAL :: &
      SigmaA_Option
    REAL(DP),         INTENT(in), OPTIONAL :: &
      SigmaS_Option
    REAL(DP) :: wTime

    wTime = MPI_WTIME( )

    WRITE(*,*)
    WRITE(*,'(A2,A6,A)') '', 'INFO: ', TRIM( ProgramName )

    SELECT CASE( TRIM( ProgramName ) )

      CASE ( 'SineWaveStreaming' )

        CALL InitializeFields_SineWaveStreaming &
               ( Direction_Option = Direction_Option )

      CASE ( 'SquareWaveStreaming' )

        CALL InitializeFields_SquareWaveStreaming &
               ( Direction_Option = Direction_Option )

      CASE ( 'SineWaveDamping' )

        CALL InitializeFields_SineWaveDamping

      CASE ( 'SineWaveDiffusion' )

        CALL InitializeFields_SineWaveDiffusion &
               ( SigmaS_Option = SigmaS_Option )

      CASE ( 'PackedBeam' )

        CALL InitializeFields_PackedBeam

      CASE ( 'CrossingBeams' )

        CALL InitializeFields_CrossingBeams

      CASE ( 'LineSource' )

        CALL InitializeFields_LineSource

      CASE ( 'RiemannProblem' )

        CALL InitializeFields_RiemannProblem &
               ( Direction_Option = Direction_Option )

      CASE ( 'HomogeneousSphere' )

        CALL InitializeFields_HomogeneousSphere

      CASE ( 'HomogeneousSphere_Spherical' )

        CALL InitializeFields_HomogeneousSphere_Spherical

      CASE ( 'DeleptonizationWave' )

        CALL InitializeFields_DeleptonizationWave( Profile_Option )

      CASE ( 'DeleptonizationWave_Spherical' )

        CALL InitializeFields_DeleptonizationWave_Spherical( Profile_Option )

    END SELECT

    wTime = MPI_WTIME( ) - wTime

!!$    WRITE(*,*)
!!$    WRITE(*,'(A4,A,ES10.4E2)') &
!!$      '', 'InitializeFields: ', wTime
!!$    WRITE(*,*)

  END SUBROUTINE InitializeFields

  SUBROUTINE InitializeFields_SineWaveStreaming( Direction_Option )

    CHARACTER(LEN=*), INTENT(in), OPTIONAL :: &
      Direction_Option

    CHARACTER(8) :: Direction
    INTEGER      :: iE, iX1, iX2, iX3, iS
    INTEGER      :: iNodeX1
    INTEGER      :: iNodeX2
    INTEGER      :: iNodeX3
    INTEGER      :: iNode
    REAL(DP)     :: X1, X2, X3, X, L
    REAL(DP)     :: Gm_dd_11(nDOF)
    REAL(DP)     :: Gm_dd_22(nDOF)
    REAL(DP)     :: Gm_dd_33(nDOF)
    REAL(DP)     :: Ones(nDOFE)

    Direction = 'X'
    IF( PRESENT( Direction_Option ) ) &
      Direction = TRIM( Direction_Option )

    WRITE(*,*)
    WRITE(*,'(A6,A,A)') '', 'Direction = ', TRIM( Direction )

    Ones = 1.0_DP

    DO iS = 1, nSpecies
      DO iX3 = iX_B0(3), iX_E0(3)
        DO iX2 = iX_B0(2), iX_E0(2)
          DO iX1 = iX_B0(1), iX_E0(1)

            Gm_dd_11 &
              = OuterProduct1D3D &
                  ( Ones, nDOFE, uGF(:,iX1,iX2,iX3,iGF_Gm_dd_11), nDOFX )

            Gm_dd_22 &
              = OuterProduct1D3D &
                  ( Ones, nDOFE, uGF(:,iX1,iX2,iX3,iGF_Gm_dd_22), nDOFX )

            Gm_dd_33 &
              = OuterProduct1D3D &
                  ( Ones, nDOFE, uGF(:,iX1,iX2,iX3,iGF_Gm_dd_33), nDOFX )

            DO iE = iE_B0, iE_E0

              SELECT CASE ( TRIM( Direction ) )
              CASE ( 'X' )

                DO iNode = 1, nDOF

                  iNodeX1 = NodeNumberTable(2,iNode)

                  X1 = NodeCoordinate( MeshX(1), iX1, iNodeX1 )

                  uPR(iNode,iE,iX1,iX2,iX3,iPR_D,iS) &
                    = 0.50_DP + 0.49_DP * SIN( TwoPi * X1 )

                  uPR(iNode,iE,iX1,iX2,iX3,iPR_I1,iS) &
                    = uPR(iNode,iE,iX1,iX2,iX3,iPR_D,iS)

                  uPR(iNode,iE,iX1,iX2,iX3,iPR_I2,iS) &
                    = 0.0_DP

                  uPR(iNode,iE,iX1,iX2,iX3,iPR_I3,iS) &
                    = 0.0_DP

                END DO

              CASE ( 'Y' )

                DO iNode = 1, nDOF

                  iNodeX2 = NodeNumberTable(3,iNode)

                  X2 = NodeCoordinate( MeshX(2), iX2, iNodeX2 )

                  uPR(iNode,iE,iX1,iX2,iX3,iPR_D,iS) &
                    = 0.50_DP + 0.49_DP * SIN( TwoPi * X2 )

                  uPR(iNode,iE,iX1,iX2,iX3,iPR_I1,iS) &
                    = 0.0_DP

                  uPR(iNode,iE,iX1,iX2,iX3,iPR_I2,iS) &
                    = uPR(iNode,iE,iX1,iX2,iX3,iPR_D,iS)

                  uPR(iNode,iE,iX1,iX2,iX3,iPR_I3,iS) &
                    = 0.0_DP

                END DO

              CASE ( 'Z' )

                DO iNode = 1, nDOF

                  iNodeX3 = NodeNumberTable(4,iNode)

                  X3 = NodeCoordinate( MeshX(3), iX3, iNodeX3 )

                  uPR(iNode,iE,iX1,iX2,iX3,iPR_D,iS) &
                    = 0.50_DP + 0.49_DP * SIN( TwoPi * X3 )

                  uPR(iNode,iE,iX1,iX2,iX3,iPR_I1,iS) &
                    = 0.0_DP

                  uPR(iNode,iE,iX1,iX2,iX3,iPR_I2,iS) &
                    = 0.0_DP

                  uPR(iNode,iE,iX1,iX2,iX3,iPR_I3,iS) &
                    = uPR(iNode,iE,iX1,iX2,iX3,iPR_D,iS)

                END DO

              CASE ( 'XY' )

                DO iNode = 1, nDOF

                  iNodeX1 = NodeNumberTable(2,iNode)
                  iNodeX2 = NodeNumberTable(3,iNode)

                  X1 = NodeCoordinate( MeshX(1), iX1, iNodeX1 )
                  X2 = NodeCoordinate( MeshX(2), iX2, iNodeX2 )

                  X  = SQRT( 2.0_DP ) * X1 &
                         + SQRT( 2.0_DP ) * X2
                  L  = SQRT( 2.0_DP )

                  uPR(iNode,iE,iX1,iX2,iX3,iPR_D,iS) &
                    = 0.50_DP + 0.49_DP * SIN( TwoPi * X / L )

                  uPR(iNode,iE,iX1,iX2,iX3,iPR_I1,iS) &
                    = SQRT( 2.0_DP ) / 2.0_DP &
                        * uPR(iNode,iE,iX1,iX2,iX3,iPR_D,iS)

                  uPR(iNode,iE,iX1,iX2,iX3,iPR_I2,iS) &
                    = SQRT( 2.0_DP ) / 2.0_DP &
                        * uPR(iNode,iE,iX1,iX2,iX3,iPR_D,iS)

                  uPR(iNode,iE,iX1,iX2,iX3,iPR_I3,iS) &
                    = 0.0_DP

                END DO

              CASE DEFAULT

                WRITE(*,*)
                WRITE(*,'(A6,A,A)') &
                  '', 'Invalid Direction: ', TRIM( Direction )
                WRITE(*,*)
                STOP

              END SELECT

              CALL ComputeConserved_TwoMoment &
                     ( uPR(:,iE,iX1,iX2,iX3,iPR_D, iS), &
                       uPR(:,iE,iX1,iX2,iX3,iPR_I1,iS), &
                       uPR(:,iE,iX1,iX2,iX3,iPR_I2,iS), &
                       uPR(:,iE,iX1,iX2,iX3,iPR_I3,iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_N, iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_G1,iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_G2,iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_G3,iS), &
                       Gm_dd_11(:), Gm_dd_22(:), Gm_dd_33(:) )

            END DO
          END DO
        END DO
      END DO
    END DO

  END SUBROUTINE InitializeFields_SineWaveStreaming

  SUBROUTINE InitializeFields_SquareWaveStreaming( Direction_Option )

    CHARACTER(LEN=*), INTENT(in), OPTIONAL :: &
      Direction_Option

    CHARACTER(8) :: Direction
    INTEGER      :: iE, iX1, iX2, iX3, iS
    INTEGER      :: iNodeX1
    INTEGER      :: iNodeX2
    INTEGER      :: iNodeX3
    INTEGER      :: iNode
    REAL(DP)     :: X1, X2, X3, X, L
    REAL(DP)     :: Gm_dd_11(nDOF)
    REAL(DP)     :: Gm_dd_22(nDOF)
    REAL(DP)     :: Gm_dd_33(nDOF)
    REAL(DP)     :: Ones(nDOFE)

    Direction = 'X'
    IF( PRESENT( Direction_Option ) ) &
      Direction = TRIM( Direction_Option )

    WRITE(*,*)
    WRITE(*,'(A6,A,A)') '', 'Direction = ', TRIM( Direction )

    Ones = 1.0_DP

    DO iS = 1, nSpecies
      DO iX3 = iX_B0(3), iX_E0(3)
        DO iX2 = iX_B0(2), iX_E0(2)
          DO iX1 = iX_B0(1), iX_E0(1)

            Gm_dd_11 &
              = OuterProduct1D3D &
                  ( Ones, nDOFE, uGF(:,iX1,iX2,iX3,iGF_Gm_dd_11), nDOFX )

            Gm_dd_22 &
              = OuterProduct1D3D &
                  ( Ones, nDOFE, uGF(:,iX1,iX2,iX3,iGF_Gm_dd_22), nDOFX )

            Gm_dd_33 &
              = OuterProduct1D3D &
                  ( Ones, nDOFE, uGF(:,iX1,iX2,iX3,iGF_Gm_dd_33), nDOFX )

            DO iE = iE_B0, iE_E0

              SELECT CASE ( TRIM( Direction ) )
              CASE ( 'X' )

                DO iNode = 1, nDOF

                  iNodeX1 = NodeNumberTable(2,iNode)

                  X1 = NodeCoordinate( MeshX(1), iX1, iNodeX1 )

                  IF( (X1 <= 0.40_DP) .AND. (X1 > 0.10_DP) ) THEN
                    uPR(iNode,iE,iX1,iX2,iX3,iPR_D,iS) &
                       = 0.80_DP
                  ELSE
                    uPR(iNode,iE,iX1,iX2,iX3,iPR_D,iS) &
                       = 0.20_DP
                  END IF

                  uPR(iNode,iE,iX1,iX2,iX3,iPR_I1,iS) &
                    = uPR(iNode,iE,iX1,iX2,iX3,iPR_D,iS) &
                      * ( 1.d0 - 1d-5 )

                  uPR(iNode,iE,iX1,iX2,iX3,iPR_I2,iS) &
                    = 0.0_DP

                  uPR(iNode,iE,iX1,iX2,iX3,iPR_I3,iS) &
                    = 0.0_DP

                END DO

              CASE DEFAULT

                WRITE(*,*)
                WRITE(*,'(A6,A,A)') &
                  '', 'Invalid Direction: ', TRIM( Direction )
                WRITE(*,*)
                STOP

              END SELECT

              CALL ComputeConserved_TwoMoment &
                     ( uPR(:,iE,iX1,iX2,iX3,iPR_D, iS), &
                       uPR(:,iE,iX1,iX2,iX3,iPR_I1,iS), &
                       uPR(:,iE,iX1,iX2,iX3,iPR_I2,iS), &
                       uPR(:,iE,iX1,iX2,iX3,iPR_I3,iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_N, iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_G1,iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_G2,iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_G3,iS), &
                       Gm_dd_11(:), Gm_dd_22(:), Gm_dd_33(:) )

            END DO
          END DO
        END DO
      END DO
    END DO

  END SUBROUTINE InitializeFields_SquareWaveStreaming

  SUBROUTINE InitializeFields_SineWaveDamping

    INTEGER      :: iE, iX1, iX2, iX3, iS, iK
    INTEGER      :: iNodeX1
    INTEGER      :: iNode
    REAL(DP)     :: X1
    REAL(DP)     :: Gm_dd_11(nDOF)
    REAL(DP)     :: Gm_dd_22(nDOF)
    REAL(DP)     :: Gm_dd_33(nDOF)
    REAL(DP)     :: Ones(nDOFE)

    Ones = 1.0_DP

    DO iS = 1, nSpecies
      DO iX3 = iX_B0(3), iX_E0(3)
        DO iX2 = iX_B0(2), iX_E0(2)
          DO iX1 = iX_B0(1), iX_E0(1)

            Gm_dd_11 &
              = OuterProduct1D3D &
                  ( Ones, nDOFE, uGF(:,iX1,iX2,iX3,iGF_Gm_dd_11), nDOFX )

            Gm_dd_22 &
              = OuterProduct1D3D &
                  ( Ones, nDOFE, uGF(:,iX1,iX2,iX3,iGF_Gm_dd_22), nDOFX )

            Gm_dd_33 &
              = OuterProduct1D3D &
                  ( Ones, nDOFE, uGF(:,iX1,iX2,iX3,iGF_Gm_dd_33), nDOFX )

            DO iE = iE_B0, iE_E0

              DO iNode = 1, nDOF

                iNodeX1 = NodeNumberTable(2,iNode)

                X1 = NodeCoordinate( MeshX(1), iX1, iNodeX1 )

                uPR(iNode,iE,iX1,iX2,iX3,iPR_D,iS) &
                  = 0.5_DP + 0.49_DP * SIN( TwoPi * X1 )

                uPR(iNode,iE,iX1,iX2,iX3,iPR_I1,iS) &
                  = uPR(iNode,iE,iX1,iX2,iX3,iPR_D,iS)

                uPR(iNode,iE,iX1,iX2,iX3,iPR_I2,iS) &
                  = 0.0_DP

                uPR(iNode,iE,iX1,iX2,iX3,iPR_I3,iS) &
                  = 0.0_DP

              END DO

              CALL ComputeConserved_TwoMoment &
                     ( uPR(:,iE,iX1,iX2,iX3,iPR_D, iS), &
                       uPR(:,iE,iX1,iX2,iX3,iPR_I1,iS), &
                       uPR(:,iE,iX1,iX2,iX3,iPR_I2,iS), &
                       uPR(:,iE,iX1,iX2,iX3,iPR_I3,iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_N, iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_G1,iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_G2,iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_G3,iS), &
                       Gm_dd_11(:), Gm_dd_22(:), Gm_dd_33(:) )

            END DO
          END DO
        END DO
      END DO
    END DO

  END SUBROUTINE InitializeFields_SineWaveDamping

  SUBROUTINE InitializeFields_SineWaveDiffusion( SigmaS_Option )

    REAL(DP), INTENT(in), OPTIONAL :: &
      SigmaS_Option

    INTEGER      :: iE, iX1, iX2, iX3, iS, iK
    INTEGER      :: iNodeX1
    INTEGER      :: iNode
    REAL(DP) :: SigmaS
    REAL(DP)     :: X1
    REAL(DP)     :: Gm_dd_11(nDOF)
    REAL(DP)     :: Gm_dd_22(nDOF)
    REAL(DP)     :: Gm_dd_33(nDOF)
    REAL(DP)     :: Ones(nDOFE)

    SigmaS = 1.0d2
    IF( PRESENT( SigmaS_Option ) ) &
      SigmaS = SigmaS_Option

    Ones = 1.0_DP

    DO iS = 1, nSpecies
      DO iX3 = iX_B0(3), iX_E0(3)
        DO iX2 = iX_B0(2), iX_E0(2)
          DO iX1 = iX_B0(1), iX_E0(1)

            Gm_dd_11 &
              = OuterProduct1D3D &
                  ( Ones, nDOFE, uGF(:,iX1,iX2,iX3,iGF_Gm_dd_11), nDOFX )

            Gm_dd_22 &
              = OuterProduct1D3D &
                  ( Ones, nDOFE, uGF(:,iX1,iX2,iX3,iGF_Gm_dd_22), nDOFX )

            Gm_dd_33 &
              = OuterProduct1D3D &
                  ( Ones, nDOFE, uGF(:,iX1,iX2,iX3,iGF_Gm_dd_33), nDOFX )

            DO iE = iE_B0, iE_E0

              DO iNode = 1, nDOF

                iNodeX1 = NodeNumberTable(2,iNode)

                X1 = NodeCoordinate( MeshX(1), iX1, iNodeX1 )

                uPR(iNode,iE,iX1,iX2,iX3,iPR_D,iS) &
                  = 0.49_DP * SIN( Third * Pi * X1 ) + 0.5_DP

                uPR(iNode,iE,iX1,iX2,iX3,iPR_I1,iS) &
                  = - ( 0.49_DP * Pi / ( 9.0_DP * SigmaS ) ) * COS( Third * Pi * X1 )

                uPR(iNode,iE,iX1,iX2,iX3,iPR_I2,iS) &
                  = 0.0_DP

                uPR(iNode,iE,iX1,iX2,iX3,iPR_I3,iS) &
                  = 0.0_DP

              END DO

              CALL ComputeConserved_TwoMoment &
                     ( uPR(:,iE,iX1,iX2,iX3,iPR_D, iS), &
                       uPR(:,iE,iX1,iX2,iX3,iPR_I1,iS), &
                       uPR(:,iE,iX1,iX2,iX3,iPR_I2,iS), &
                       uPR(:,iE,iX1,iX2,iX3,iPR_I3,iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_N, iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_G1,iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_G2,iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_G3,iS), &
                       Gm_dd_11(:), Gm_dd_22(:), Gm_dd_33(:) )

            END DO
          END DO
        END DO
      END DO
    END DO

  END SUBROUTINE InitializeFields_SineWaveDiffusion

  SUBROUTINE InitializeFields_PackedBeam

    INTEGER   :: iE, iX1, iX2, iX3, iS
    INTEGER   :: iNodeX1
    INTEGER   :: iNode
    REAL(DP)  :: X1, Mu_MP
    REAL(DP)  :: Delta
    REAL(DP)  :: Gm_dd_11(nDOF)
    REAL(DP)  :: Gm_dd_22(nDOF)
    REAL(DP)  :: Gm_dd_33(nDOF)
    REAL(DP)  :: Ones(nDOFE)

    Ones = 1.0_DP

    Mu_MP = 0.0_DP
    Delta = 1.0d-8

    DO iS = 1, nSpecies
      DO iX3 = iX_B0(3), iX_E0(3)
        DO iX2 = iX_B0(2), iX_E0(2)
          DO iX1 = iX_B0(1), iX_E0(1)

            Gm_dd_11 &
              = OuterProduct1D3D &
                  ( Ones, nDOFE, uGF(:,iX1,iX2,iX3,iGF_Gm_dd_11), nDOFX )

            Gm_dd_22 &
              = OuterProduct1D3D &
                  ( Ones, nDOFE, uGF(:,iX1,iX2,iX3,iGF_Gm_dd_22), nDOFX )

            Gm_dd_33 &
              = OuterProduct1D3D &
                  ( Ones, nDOFE, uGF(:,iX1,iX2,iX3,iGF_Gm_dd_33), nDOFX )

            DO iE = iE_B0, iE_E0

              DO iNode = 1, nDOF

                iNodeX1 = NodeNumberTable(2,iNode)

                X1 = NodeCoordinate( MeshX(1), iX1, iNodeX1 )

                IF( X1 .LE. ZERO )THEN

                  uPR(iNode,iE,iX1,iX2,iX3,iPR_D,iS) &
                    = 0.50_DP * ( 1.0_DP + Delta )

                  uPR(iNode,iE,iX1,iX2,iX3,iPR_I1,iS) &
                    = 0.25_DP * ( 1.0_DP - Delta )

                ELSE

                  uPR(iNode,iE,iX1,iX2,iX3,iPR_D,iS) &
                    = Delta

                  uPR(iNode,iE,iX1,iX2,iX3,iPR_I1,iS) &
                    = 0.0d-0

                END IF

                uPR(iNode,iE,iX1,iX2,iX3,iPR_I2,iS) &
                  = 0.0_DP

                uPR(iNode,iE,iX1,iX2,iX3,iPR_I3,iS) &
                  = 0.0_DP

              END DO

              CALL ComputeConserved_TwoMoment &
                     ( uPR(:,iE,iX1,iX2,iX3,iPR_D, iS), &
                       uPR(:,iE,iX1,iX2,iX3,iPR_I1,iS), &
                       uPR(:,iE,iX1,iX2,iX3,iPR_I2,iS), &
                       uPR(:,iE,iX1,iX2,iX3,iPR_I3,iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_N, iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_G1,iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_G2,iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_G3,iS), &
                       Gm_dd_11(:), Gm_dd_22(:), Gm_dd_33(:) )

            END DO
          END DO
        END DO
      END DO
    END DO

  END SUBROUTINE InitializeFields_PackedBeam

  SUBROUTINE InitializeFields_CrossingBeams

    INTEGER  :: iE, iX1, iX2, iX3, iS
    INTEGER  :: iNode
    REAL(DP) :: Ones(nDOF)

    Ones = 1.0_DP

    DO iS = 1, nSpecies
      DO iX3 = iX_B0(3), iX_E0(3)
        DO iX2 = iX_B0(2), iX_E0(2)
          DO iX1 = iX_B0(1), iX_E0(1)
            DO iE = iE_B0, iE_E0

              DO iNode = 1, nDOF

                uPR(iNode,iE,iX1,iX2,iX3,iPR_D,iS) &
                  = 1.0d-12

                uPR(iNode,iE,iX1,iX2,iX3,iPR_I1,iS) &
                  = 0.0_DP

                uPR(iNode,iE,iX1,iX2,iX3,iPR_I2,iS) &
                  = 0.0_DP

                uPR(iNode,iE,iX1,iX2,iX3,iPR_I3,iS) &
                  = 0.0_DP

              END DO

              CALL ComputeConserved_TwoMoment &
                     ( uPR(:,iE,iX1,iX2,iX3,iPR_D, iS), &
                       uPR(:,iE,iX1,iX2,iX3,iPR_I1,iS), &
                       uPR(:,iE,iX1,iX2,iX3,iPR_I2,iS), &
                       uPR(:,iE,iX1,iX2,iX3,iPR_I3,iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_N, iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_G1,iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_G2,iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_G3,iS), &
                       Ones(:), Ones(:), Ones(:) )

            END DO
          END DO
        END DO
      END DO
    END DO

  END SUBROUTINE InitializeFields_CrossingBeams

  SUBROUTINE InitializeFields_LineSource

    INTEGER   :: iE, iX1, iX2, iX3, iS
    INTEGER   :: iNodeX1
    INTEGER   :: iNodeX2
    INTEGER   :: iNode
    REAL(DP)  :: X1, X2, R
    REAL(DP)  :: Gm_dd_11(nDOF)
    REAL(DP)  :: Gm_dd_22(nDOF)
    REAL(DP)  :: Gm_dd_33(nDOF)
    REAL(DP)  :: Ones(nDOFE)

    Ones = 1.0_DP

    DO iS = 1, nSpecies
      DO iX3 = iX_B0(3), iX_E0(3)
        DO iX2 = iX_B0(2), iX_E0(2)
          DO iX1 = iX_B0(1), iX_E0(1)

            Gm_dd_11 &
              = OuterProduct1D3D &
                  ( Ones, nDOFE, uGF(:,iX1,iX2,iX3,iGF_Gm_dd_11), nDOFX )

            Gm_dd_22 &
              = OuterProduct1D3D &
                  ( Ones, nDOFE, uGF(:,iX1,iX2,iX3,iGF_Gm_dd_22), nDOFX )

            Gm_dd_33 &
              = OuterProduct1D3D &
                  ( Ones, nDOFE, uGF(:,iX1,iX2,iX3,iGF_Gm_dd_33), nDOFX )

            DO iE = iE_B0, iE_E0

              DO iNode = 1, nDOF

                iNodeX1 = NodeNumberTable(2,iNode)
                iNodeX2 = NodeNumberTable(3,iNode)

                X1 = NodeCoordinate( MeshX(1), iX1, iNodeX1 )
                X2 = NodeCoordinate( MeshX(2), iX2, iNodeX2 )

                R  = SQRT( X1**2 + X2**2 )

                uPR(iNode,iE,iX1,iX2,iX3,iPR_D,iS) &
                  = 1.0_DP &
                    - MAX( EXP( - 0.5_DP * ( R / 0.03_DP )**2 ), 1.0d-8 )

                uPR(iNode,iE,iX1,iX2,iX3,iPR_I1,iS) &
                  = 0.0_DP

                uPR(iNode,iE,iX1,iX2,iX3,iPR_I2,iS) &
                  = 0.0_DP

                uPR(iNode,iE,iX1,iX2,iX3,iPR_I3,iS) &
                  = 0.0_DP

              END DO

              CALL ComputeConserved_TwoMoment &
                     ( uPR(:,iE,iX1,iX2,iX3,iPR_D, iS), &
                       uPR(:,iE,iX1,iX2,iX3,iPR_I1,iS), &
                       uPR(:,iE,iX1,iX2,iX3,iPR_I2,iS), &
                       uPR(:,iE,iX1,iX2,iX3,iPR_I3,iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_N, iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_G1,iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_G2,iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_G3,iS), &
                       Gm_dd_11(:), Gm_dd_22(:), Gm_dd_33(:) )

            END DO
          END DO
        END DO
      END DO
    END DO

  END SUBROUTINE InitializeFields_LineSource

  SUBROUTINE InitializeFields_RiemannProblem( Direction_Option )

    CHARACTER(LEN=*), INTENT(in), OPTIONAL :: &
      Direction_Option

    CHARACTER(8) :: Direction
    INTEGER      :: iE, iX1, iX2, iX3, iS, iK
    INTEGER      :: iNodeX1, iNodeX2
    INTEGER      :: iNode
    REAL(DP)     :: X1, X2
    REAL(DP)     :: Gm_dd_11(nDOF)
    REAL(DP)     :: Gm_dd_22(nDOF)
    REAL(DP)     :: Gm_dd_33(nDOF)
    REAL(DP)     :: Ones(nDOFE)

    Direction = 'X'
    IF( PRESENT( Direction_Option ) ) &
      Direction = TRIM( Direction_Option )

    WRITE(*,*)
    WRITE(*,'(A6,A,A)') '', 'Direction = ', TRIM( Direction )

    Ones = 1.0_DP

    DO iS = 1, nSpecies
      DO iX3 = iX_B0(3), iX_E0(3)
        DO iX2 = iX_B0(2), iX_E0(2)
          DO iX1 = iX_B0(1), iX_E0(1)

            Gm_dd_11 &
              = OuterProduct1D3D &
                  ( Ones, nDOFE, uGF(:,iX1,iX2,iX3,iGF_Gm_dd_11), nDOFX )

            Gm_dd_22 &
              = OuterProduct1D3D &
                  ( Ones, nDOFE, uGF(:,iX1,iX2,iX3,iGF_Gm_dd_22), nDOFX )

            Gm_dd_33 &
              = OuterProduct1D3D &
                  ( Ones, nDOFE, uGF(:,iX1,iX2,iX3,iGF_Gm_dd_33), nDOFX )

            DO iE = iE_B0, iE_E0

              SELECT CASE ( TRIM( Direction ) )
              CASE ( 'X' )

                DO iNode = 1, nDOF

                  iNodeX1 = NodeNumberTable(2,iNode)

                  X1 = NodeCoordinate( MeshX(1), iX1, iNodeX1 )

                  IF( X1 <= 0.0_DP ) THEN
                    uPR(iNode,iE,iX1,iX2,iX3,iPR_D,iS)  &
                      = 1.0_DP
                    uPR(iNode,iE,iX1,iX2,iX3,iPR_I1,iS) &
                      = 0.9999_DP
                  ELSE
                    uPR(iNode,iE,iX1,iX2,iX3,iPR_D,iS)  &
                      = 0.50_DP
                    uPR(iNode,iE,iX1,iX2,iX3,iPR_I1,iS) &
                      = 0.0_DP
                  END IF

                  uPR(iNode,iE,iX1,iX2,iX3,iPR_I2,iS) &
                    = 0.0_DP

                  uPR(iNode,iE,iX1,iX2,iX3,iPR_I3,iS) &
                    = 0.0_DP

                END DO

              CASE ( 'Y' )

                DO iNode = 1, nDOF

                  iNodeX2 = NodeNumberTable(3,iNode)

                  X2 = NodeCoordinate( MeshX(2), iX2, iNodeX2 )

                  IF( X2 <= 0.0_DP ) THEN
                    uPR(iNode,iE,iX1,iX2,iX3,iPR_D,iS)  &
                      = 1.0_DP
                    uPR(iNode,iE,iX1,iX2,iX3,iPR_I2,iS) &
                      = 0.9999_DP
                  ELSE
                    uPR(iNode,iE,iX1,iX2,iX3,iPR_D,iS)  &
                      = 0.50_DP
                    uPR(iNode,iE,iX1,iX2,iX3,iPR_I2,iS) &
                      = 0.0_DP
                  END IF

                  uPR(iNode,iE,iX1,iX2,iX3,iPR_I1,iS) &
                    = 0.0_DP

                  uPR(iNode,iE,iX1,iX2,iX3,iPR_I3,iS) &
                    = 0.0_DP

                END DO

              END SELECT

              CALL ComputeConserved_TwoMoment &
                     ( uPR(:,iE,iX1,iX2,iX3,iPR_D, iS), &
                       uPR(:,iE,iX1,iX2,iX3,iPR_I1,iS), &
                       uPR(:,iE,iX1,iX2,iX3,iPR_I2,iS), &
                       uPR(:,iE,iX1,iX2,iX3,iPR_I3,iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_N, iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_G1,iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_G2,iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_G3,iS), &
                       Gm_dd_11(:), Gm_dd_22(:), Gm_dd_33(:) )

            END DO
          END DO
        END DO
      END DO
    END DO

  END SUBROUTINE InitializeFields_RiemannProblem

  SUBROUTINE InitializeFields_HomogeneousSphere

    INTEGER   :: iE, iX1, iX2, iX3, iS
    INTEGER   :: iNode
    REAL(DP)  :: Gm_dd_11(nDOF)
    REAL(DP)  :: Gm_dd_22(nDOF)
    REAL(DP)  :: Gm_dd_33(nDOF)
    REAL(DP)  :: Ones(nDOFE)

    Ones = 1.0_DP

    DO iS = 1, nSpecies
      DO iX3 = iX_B0(3), iX_E0(3)
        DO iX2 = iX_B0(2), iX_E0(2)
          DO iX1 = iX_B0(1), iX_E0(1)

            Gm_dd_11 &
              = OuterProduct1D3D &
                  ( Ones, nDOFE, uGF(:,iX1,iX2,iX3,iGF_Gm_dd_11), nDOFX )

            Gm_dd_22 &
              = OuterProduct1D3D &
                  ( Ones, nDOFE, uGF(:,iX1,iX2,iX3,iGF_Gm_dd_22), nDOFX )

            Gm_dd_33 &
              = OuterProduct1D3D &
                  ( Ones, nDOFE, uGF(:,iX1,iX2,iX3,iGF_Gm_dd_33), nDOFX )

            DO iE = iE_B0, iE_E0

              DO iNode = 1, nDOF

                uPR(iNode,iE,iX1,iX2,iX3,iPR_D,iS) &
                  = 1.0d-12

                uPR(iNode,iE,iX1,iX2,iX3,iPR_I1,iS) &
                  = 0.0_DP

                uPR(iNode,iE,iX1,iX2,iX3,iPR_I2,iS) &
                  = 0.0_DP

                uPR(iNode,iE,iX1,iX2,iX3,iPR_I3,iS) &
                  = 0.0_DP

              END DO

              CALL ComputeConserved_TwoMoment &
                     ( uPR(:,iE,iX1,iX2,iX3,iPR_D, iS), &
                       uPR(:,iE,iX1,iX2,iX3,iPR_I1,iS), &
                       uPR(:,iE,iX1,iX2,iX3,iPR_I2,iS), &
                       uPR(:,iE,iX1,iX2,iX3,iPR_I3,iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_N, iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_G1,iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_G2,iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_G3,iS), &
                       Gm_dd_11(:), Gm_dd_22(:), Gm_dd_33(:) )

            END DO
          END DO
        END DO
      END DO
    END DO

  END SUBROUTINE InitializeFields_HomogeneousSphere

  SUBROUTINE InitializeFields_HomogeneousSphere_Spherical

    INTEGER  :: iE, iX1, iX2, iX3, iS
    REAL(DP) :: Gm_dd_11(nDOF)
    REAL(DP) :: Gm_dd_22(nDOF)
    REAL(DP) :: Gm_dd_33(nDOF)

    DO iS = 1, nSpecies
    DO iX3 = iX_B0(3), iX_E0(3)
    DO iX2 = iX_B0(2), iX_E0(2)
    DO iX1 = iX_B0(1), iX_E0(1)

      Gm_dd_11 = uGF(NodeNumbersX,iX1,iX2,iX3,iGF_Gm_dd_11)
      Gm_dd_22 = uGF(NodeNumbersX,iX1,iX2,iX3,iGF_Gm_dd_22)
      Gm_dd_33 = uGF(NodeNumbersX,iX1,iX2,iX3,iGF_Gm_dd_33)

      DO iE = iE_B0, iE_E0

        uPR(:,iE,iX1,iX2,iX3,iPR_D, iS) = 1.0d-6
        uPR(:,iE,iX1,iX2,iX3,iPR_I1,iS) = 0.0_DP
        uPR(:,iE,iX1,iX2,iX3,iPR_I2,iS) = 0.0_DP
        uPR(:,iE,iX1,iX2,iX3,iPR_I3,iS) = 0.0_DP

        CALL ComputeConserved_TwoMoment &
               ( uPR(:,iE,iX1,iX2,iX3,iPR_D, iS), &
                 uPR(:,iE,iX1,iX2,iX3,iPR_I1,iS), &
                 uPR(:,iE,iX1,iX2,iX3,iPR_I2,iS), &
                 uPR(:,iE,iX1,iX2,iX3,iPR_I3,iS), &
                 uCR(:,iE,iX1,iX2,iX3,iCR_N, iS), &
                 uCR(:,iE,iX1,iX2,iX3,iCR_G1,iS), &
                 uCR(:,iE,iX1,iX2,iX3,iCR_G2,iS), &
                 uCR(:,iE,iX1,iX2,iX3,iCR_G3,iS), &
                 Gm_dd_11, Gm_dd_22, Gm_dd_33 )

      END DO

    END DO
    END DO
    END DO
    END DO

  END SUBROUTINE InitializeFields_HomogeneousSphere_Spherical

  SUBROUTINE ComputeError( Time, SigmaA, SigmaS )

    REAL(DP), INTENT(in) :: Time
    REAL(DP), INTENT(in) :: SigmaA
    REAL(DP), INTENT(in) :: SigmaS

    LOGICAL  :: ReportError
    REAL(DP) :: Error_One(nCR)
    REAL(DP) :: Error_Inf(nCR)

    SELECT CASE( TRIM( ProgramName ) )

      CASE ( 'SineWaveStreaming' )

        ReportError = .TRUE.

        CALL ComputeError_SineWaveStreaming &
               ( Time, Error_One, Error_Inf )

      CASE ( 'SineWaveDamping' )

        ReportError = .TRUE.

        CALL ComputeError_SineWaveDamping &
               ( Time, SigmaA, Error_One, Error_Inf )

      CASE ( 'SineWaveDiffusion' )

        ReportError = .TRUE.

        CALL ComputeError_SineWaveDiffusion &
               ( Time, SigmaS, Error_One, Error_Inf )

      CASE DEFAULT

        ReportError = .FALSE.

    END SELECT

    IF( ReportError )THEN

      WRITE(*,*)
      WRITE(*,'(A4,A)') '', '--- Error Analysis ---'
      WRITE(*,*)
      WRITE(*,'(A4,A10,3I5.4)') '', 'nX = ', nZ(2:4)
      WRITE(*,'(A4,A10,1I5.4)') '', 'nE = ', nZ(1)
      WRITE(*,'(A4,A10,1I5.4)') '', 'nNodes = ', nNodes
      WRITE(*,*)
      WRITE(*,'(A4,A12,A4,A12)') '', 'L_1', '', 'L_Inf'
      WRITE(*,*)
      WRITE(*,'(A4,A4,ES12.6E2,A4,ES12.6E2)') &
        '', 'N: ',  Error_One(iCR_N ), '', Error_Inf(iCR_N )
      WRITE(*,'(A4,A4,ES12.6E2,A4,ES12.6E2)') &
        '', 'G1: ', Error_One(iCR_G1), '', Error_Inf(iCR_G1)
      WRITE(*,'(A4,A4,ES12.6E2,A4,ES12.6E2)') &
        '', 'G2: ', Error_One(iCR_G2), '', Error_Inf(iCR_G2)
      WRITE(*,'(A4,A4,ES12.6E2,A4,ES12.6E2)') &
        '', 'G3: ', Error_One(iCR_G3), '', Error_Inf(iCR_G3)
      WRITE(*,*)

    END IF

  END SUBROUTINE ComputeError

  SUBROUTINE ComputeError_SineWaveStreaming( Time, Error_One, Error_Inf )

    REAL(DP), INTENT(in)  :: Time
    REAL(DP), INTENT(out) :: Error_One(nCR)
    REAL(DP), INTENT(out) :: Error_Inf(nCR)

    INTEGER  :: iE, iX1, iX2, iX3, iS
    INTEGER  :: iNodeX1, iNode
    REAL(DP) :: X1, N_A, G1_A, G2_A, G3_A

    Error_One = Zero
    Error_Inf = Zero

    DO iS = 1, nSpecies
      DO iX3 = 1, nZ(4)
        DO iX2 = 1, nZ(3)
          DO iX1 = 1, nZ(2)
            DO iE  = 1, nZ(1)

              DO iNode = 1, nDOF

                iNodeX1 = NodeNumberTable(2,iNode)

                X1 = NodeCoordinate( MeshX(1), iX1, iNodeX1 )

                N_A  = 0.50_DP + 0.49_DP * SIN( TwoPi * ( X1 - Time ) )
                G1_A = N_A
                G2_A = Zero
                G3_A = Zero

                ! --- L1 Error ---

                Error_One(iCR_N) &
                  = Error_One(iCR_N) &
                    + ABS( N_A - uCR(iNode,iE,iX1,iX2,iX3,iCR_N,iS) )

                Error_One(iCR_G1) &
                  = Error_One(iCR_G1) &
                    + ABS( G1_A - uCR(iNode,iE,iX1,iX2,iX3,iCR_G1,iS) )

                Error_One(iCR_G2) &
                  = Error_One(iCR_G2) &
                    + ABS( G2_A - uCR(iNode,iE,iX1,iX2,iX3,iCR_G2,iS) )

                Error_One(iCR_G3) &
                  = Error_One(iCR_G3) &
                    + ABS( G3_A - uCR(iNode,iE,iX1,iX2,iX3,iCR_G3,iS) )

                ! --- Infinity Error ---

                Error_Inf(iCR_N) &
                  = MAX( Error_Inf(iCR_N), &
                         ABS( N_A - uCR(iNode,iE,iX1,iX2,iX3,iCR_N,iS) ) )

                Error_Inf(iCR_G1) &
                  = MAX( Error_Inf(iCR_G1), &
                         ABS( G1_A - uCR(iNode,iE,iX1,iX2,iX3,iCR_G1,iS) ) )

                Error_Inf(iCR_G2) &
                  = MAX( Error_Inf(iCR_G2), &
                         ABS( G2_A - uCR(iNode,iE,iX1,iX2,iX3,iCR_G2,iS) ) )

                Error_Inf(iCR_G3) &
                  = MAX( Error_Inf(iCR_G3), &
                         ABS( G3_A - uCR(iNode,iE,iX1,iX2,iX3,iCR_G3,iS) ) )

              END DO

            END DO
          END DO
        END DO
      END DO
    END DO

    Error_One = Error_One / DBLE( nDOF * PRODUCT( nZ ) * nSpecies )

  END SUBROUTINE ComputeError_SineWaveStreaming

  SUBROUTINE ComputeError_SineWaveDamping( Time, SigmaA, Error_One, Error_Inf )

    REAL(DP), INTENT(in)  :: Time
    REAL(DP), INTENT(in)  :: SigmaA
    REAL(DP), INTENT(out) :: Error_One(nCR)
    REAL(DP), INTENT(out) :: Error_Inf(nCR)

    INTEGER  :: iE, iX1, iX2, iX3, iS
    INTEGER  :: iNodeX1, iNode
    REAL(DP) :: X1, N_A, G1_A, G2_A, G3_A

    Error_One = Zero
    Error_Inf = Zero

    DO iS = 1, nSpecies
      DO iX3 = 1, nZ(4)
        DO iX2 = 1, nZ(3)
          DO iX1 = 1, nZ(2)
            DO iE  = 1, nZ(1)

              DO iNode = 1, nDOF

                iNodeX1 = NodeNumberTable(2,iNode)

                X1 = NodeCoordinate( MeshX(1), iX1, iNodeX1 )

                N_A  = ( 0.50_DP + 0.49_DP * SIN( TwoPi * ( X1 - Time ) ) ) &
                         * EXP( - SigmaA * Time )
                G1_A = N_A
                G2_A = Zero
                G3_A = Zero

                ! --- L1 Error (Relative) ---

                Error_One(iCR_N) &
                  = Error_One(iCR_N) &
                    + ABS( N_A - uCR(iNode,iE,iX1,iX2,iX3,iCR_N,iS) ) / N_A

                Error_One(iCR_G1) &
                  = Error_One(iCR_G1) &
                    + ABS( G1_A - uCR(iNode,iE,iX1,iX2,iX3,iCR_G1,iS) ) / N_A

                Error_One(iCR_G2) &
                  = Error_One(iCR_G2) &
                    + ABS( G2_A - uCR(iNode,iE,iX1,iX2,iX3,iCR_G2,iS) ) / N_A

                Error_One(iCR_G3) &
                  = Error_One(iCR_G3) &
                    + ABS( G3_A - uCR(iNode,iE,iX1,iX2,iX3,iCR_G3,iS) ) / N_A

                ! --- Infinity Error (Relative) ---

                Error_Inf(iCR_N) &
                  = MAX( Error_Inf(iCR_N), &
                         ABS(N_A-uCR(iNode,iE,iX1,iX2,iX3,iCR_N,iS))/N_A )

                Error_Inf(iCR_G1) &
                  = MAX( Error_Inf(iCR_G1), &
                         ABS(G1_A-uCR(iNode,iE,iX1,iX2,iX3,iCR_G1,iS))/N_A )

                Error_Inf(iCR_G2) &
                  = MAX( Error_Inf(iCR_G2), &
                         ABS(G2_A-uCR(iNode,iE,iX1,iX2,iX3,iCR_G2,iS))/N_A )

                Error_Inf(iCR_G3) &
                  = MAX( Error_Inf(iCR_G3), &
                         ABS(G3_A-uCR(iNode,iE,iX1,iX2,iX3,iCR_G3,iS))/N_A )

              END DO

            END DO
          END DO
        END DO
      END DO
    END DO

    Error_One = Error_One / DBLE( nDOF * PRODUCT( nZ ) * nSpecies )

  END SUBROUTINE ComputeError_SineWaveDamping

  SUBROUTINE ComputeError_SineWaveDiffusion &
    ( Time, SigmaS, Error_One, Error_Inf )

    REAL(DP), INTENT(in)  :: Time
    REAL(DP), INTENT(in)  :: SigmaS
    REAL(DP), INTENT(out) :: Error_One(nCR)
    REAL(DP), INTENT(out) :: Error_Inf(nCR)

    INTEGER  :: iE, iX1, iX2, iX3, iS
    INTEGER  :: iNodeX1, iNode
    REAL(DP) :: X1, N_A, G1_A, G2_A, G3_A

    Error_One = Zero
    Error_Inf = Zero

    DO iS = 1, nSpecies
      DO iX3 = 1, nZ(4)
        DO iX2 = 1, nZ(3)
          DO iX1 = 1, nZ(2)
            DO iE  = 1, nZ(1)

              DO iNode = 1, nDOF

                iNodeX1 = NodeNumberTable(2,iNode)

                X1 = NodeCoordinate( MeshX(1), iX1, iNodeX1 )

                N_A  = ( 0.50_DP + 0.49_DP * SIN( Third * Pi * X1 ) &
                         * EXP( - Pi**2 * Time / ( 2.7d1 * SigmaS ) ) )
                G1_A = - ( 0.49_DP * Pi / ( 9.0_DP * SigmaS ) ) &
                         * COS( Third * Pi * X1 ) &
                         * EXP( - Pi**2 * Time / ( 2.7d1 * SigmaS ) )
                G2_A = Zero
                G3_A = Zero

                ! --- L1 Error ---

                Error_One(iCR_N) &
                  = Error_One(iCR_N) &
                    + ABS( N_A - uCR(iNode,iE,iX1,iX2,iX3,iCR_N,iS) )

                Error_One(iCR_G1) &
                  = Error_One(iCR_G1) &
                    + ABS( G1_A - uCR(iNode,iE,iX1,iX2,iX3,iCR_G1,iS) )

                Error_One(iCR_G2) &
                  = Error_One(iCR_G2) &
                    + ABS( G2_A - uCR(iNode,iE,iX1,iX2,iX3,iCR_G2,iS) )

                Error_One(iCR_G3) &
                  = Error_One(iCR_G3) &
                    + ABS( G3_A - uCR(iNode,iE,iX1,iX2,iX3,iCR_G3,iS) )

                ! --- Infinity Error ---

                Error_Inf(iCR_N) &
                  = MAX( Error_Inf(iCR_N), &
                         ABS(N_A-uCR(iNode,iE,iX1,iX2,iX3,iCR_N,iS)) )

                Error_Inf(iCR_G1) &
                  = MAX( Error_Inf(iCR_G1), &
                         ABS(G1_A-uCR(iNode,iE,iX1,iX2,iX3,iCR_G1,iS)) )

                Error_Inf(iCR_G2) &
                  = MAX( Error_Inf(iCR_G2), &
                         ABS(G2_A-uCR(iNode,iE,iX1,iX2,iX3,iCR_G2,iS)) )

                Error_Inf(iCR_G3) &
                  = MAX( Error_Inf(iCR_G3), &
                         ABS(G3_A-uCR(iNode,iE,iX1,iX2,iX3,iCR_G3,iS)) )

              END DO

            END DO
          END DO
        END DO
      END DO
    END DO

    Error_One = Error_One / DBLE( nDOF * PRODUCT( nZ ) * nSpecies )

  END SUBROUTINE ComputeError_SineWaveDiffusion

  SUBROUTINE InitializeFields_DeleptonizationWave( Profile_Option )

    CHARACTER(len=*), INTENT(in), OPTIONAL :: Profile_Option

    IF( PRESENT( Profile_Option ) )THEN

      CALL InitializeFluidFields_DeleptonizationWave_Profile( Profile_Option )

    ELSE

      CALL InitializeFluidFields_DeleptonizationWave

    END IF

    CALL InitializeRadiationFields_DeleptonizationWave

  END SUBROUTINE InitializeFields_DeleptonizationWave

  SUBROUTINE InitializeFields_DeleptonizationWave_Spherical( Profile_Option )

    CHARACTER(len=*), INTENT(in), OPTIONAL :: Profile_Option

    IF( PRESENT( Profile_Option ) )THEN

      CALL InitializeFluidFields_DeleptonizationWave_Spherical_Profile&
             ( Profile_Option )

    ELSE

      CALL InitializeFluidFields_DeleptonizationWave_Spherical

    END IF

    CALL InitializeRadiationFields_DeleptonizationWave

  END SUBROUTINE InitializeFields_DeleptonizationWave_Spherical

  SUBROUTINE InitializeFluidFields_DeleptonizationWave

    ! --- Density Profile ---
    REAL(DP), PARAMETER :: MinD = 1.0d08 * Gram / Centimeter**3
    REAL(DP), PARAMETER :: MaxD = 4.0d14 * Gram / Centimeter**3
    REAL(DP), PARAMETER :: C_D  = 7.5_DP
    REAL(DP), PARAMETER :: R_D  = 2.0d01 * Kilometer
!    REAL(DP), PARAMETER :: H_D  = 1.0d01 * Kilometer
    REAL(DP), PARAMETER :: H_D  = 5.0d00 * Kilometer
    ! --- Temperature Profile ---
    REAL(DP), PARAMETER :: MinT = 5.0d09 * Kelvin
    REAL(DP), PARAMETER :: MaxT = 1.5d11 * Kelvin
    REAL(DP), PARAMETER :: C_T  = 1.0_DP
    REAL(DP), PARAMETER :: R_T  = 2.5d01 * Kilometer
!    REAL(DP), PARAMETER :: H_T  = 2.0d01 * Kilometer
    REAL(DP), PARAMETER :: H_T  = 5.0d01 * Kilometer
    ! --- Electron Fraction Profile ---
    REAL(DP), PARAMETER :: MinY = 2.5d-1
    REAL(DP), PARAMETER :: MaxY = 4.6d-1
    REAL(DP), PARAMETER :: C_Y  = 1.0_DP
    REAL(DP), PARAMETER :: R_Y  = 4.5d01 * Kilometer
!    REAL(DP), PARAMETER :: H_Y  = 1.0d01 * Kilometer
    REAL(DP), PARAMETER :: H_Y  = 5.0d01 * Kilometer

    INTEGER  :: iX1, iX2, iX3, iNodeX
    INTEGER  :: iNodeX1, iNodeX2, iNodeX3
    REAL(DP) :: X1, X2, X3, R

    DO iX3 = iX_B0(3), iX_E0(3)
      DO iX2 = iX_B0(2), iX_E0(2)
        DO iX1 = iX_B0(1), iX_E0(1)

          DO iNodeX = 1, nDOFX

            iNodeX1 = NodeNumberTableX(1,iNodeX)
            iNodeX2 = NodeNumberTableX(2,iNodeX)
            iNodeX3 = NodeNumberTableX(3,iNodeX)

            X1 = NodeCoordinate( MeshX(1), iX1, iNodeX1 )
            X2 = NodeCoordinate( MeshX(2), iX2, iNodeX2 )
            X3 = NodeCoordinate( MeshX(3), iX3, iNodeX3 )

            R = SQRT( X1**2 + X2**2 + X3**2 )

            uPF(iNodeX,iX1,iX2,iX3,iPF_D) &
!!$              = Half * ( MaxD * ( One - TANH( (R-R_D)/H_D ) ) &
!!$                         + MinD * ( One - TANH( (R_D-R)/H_D ) ) )
              = MaxD * C_D / ( C_D + ( R / H_D )**4 )

            uAF(iNodeX,iX1,iX2,iX3,iAF_T) &
!!$              = Half * ( MaxT * ( One - TANH( (R-R_T)/H_T ) ) &
!!$                         + MinT * ( One - TANH( (R_T-R)/H_T ) ) )
              = MAX( MaxT * C_T / ( C_T + ( R / H_T )**2 ), 1.0d10 * Kelvin )

            uAF(iNodeX,iX1,iX2,iX3,iAF_Ye) &
!!$              = Half * ( MinY * ( One - TANH( (R-R_Y)/H_Y ) ) &
!!$                         + MaxY * ( One - TANH( (R_Y-R)/H_Y ) ) )
              = MinY * ( One + C_Y / ( C_Y + ( R / H_Y )**(-12) ) )

          END DO

          CALL ComputeThermodynamicStates_Primitive_TABLE &
                 ( uPF(:,iX1,iX2,iX3,iPF_D),  uAF(:,iX1,iX2,iX3,iAF_T), &
                   uAF(:,iX1,iX2,iX3,iAF_Ye), uPF(:,iX1,iX2,iX3,iPF_E), &
                   uAF(:,iX1,iX2,iX3,iAF_E),  uPF(:,iX1,iX2,iX3,iPF_Ne) )

          CALL ApplyEquationOfState_TABLE &
                 ( uPF(:,iX1,iX2,iX3,iPF_D ), uAF(:,iX1,iX2,iX3,iAF_T ), &
                   uAF(:,iX1,iX2,iX3,iAF_Ye), uAF(:,iX1,iX2,iX3,iAF_P ), &
                   uAF(:,iX1,iX2,iX3,iAF_S ), uAF(:,iX1,iX2,iX3,iAF_E ), &
                   uAF(:,iX1,iX2,iX3,iAF_Me), uAF(:,iX1,iX2,iX3,iAF_Mp), &
                   uAF(:,iX1,iX2,iX3,iAF_Mn), uAF(:,iX1,iX2,iX3,iAF_Xp), &
                   uAF(:,iX1,iX2,iX3,iAF_Xn), uAF(:,iX1,iX2,iX3,iAF_Xa), &
                   uAF(:,iX1,iX2,iX3,iAF_Xh), uAF(:,iX1,iX2,iX3,iAF_Gm) )

        END DO
      END DO
    END DO

  END SUBROUTINE InitializeFluidFields_DeleptonizationWave

  SUBROUTINE InitializeFluidFields_DeleptonizationWave_Profile( FileName )

    CHARACTER(len=*), INTENT(in) :: FileName

    TYPE(ProgenitorType1D) :: P1D
    CHARACTER(LEN=100)                  :: Format1, Format2, Format3
    CHARACTER(LEN=30)                   :: a
    INTEGER                             :: ii, datasize
    REAL(dp), DIMENSION(:,:), ALLOCATABLE :: database

    INTEGER  :: iX1, iX2, iX3, iNodeX
    INTEGER  :: iNodeX1, iNodeX2, iNodeX3
    REAL(DP) :: X1, X2, X3, R

    CALL ReadProgenitor1D( FileName, P1D )

    ASSOCIATE &
      ( R1D => P1D % Radius, &
        D1D => P1D % MassDensity, &
        T1D => P1D % Temperature, &
        Y1D => P1D % ElectronFraction )

    DO iX3 = iX_B0(3), iX_E0(3)
      DO iX2 = iX_B0(2), iX_E0(2)
        DO iX1 = iX_B0(1), iX_E0(1)

          DO iNodeX = 1, nDOFX

            iNodeX1 = NodeNumberTableX(1,iNodeX)
            iNodeX2 = NodeNumberTableX(2,iNodeX)
            iNodeX3 = NodeNumberTableX(3,iNodeX)

            X1 = NodeCoordinate( MeshX(1), iX1, iNodeX1 )
            X2 = NodeCoordinate( MeshX(2), iX2, iNodeX2 )
            X3 = NodeCoordinate( MeshX(3), iX3, iNodeX3 )

            R = SQRT( X1**2 + X2**2 + X3**2 )

            uPF(iNodeX,iX1,iX2,iX3,iPF_D) &
              = Interpolate1D( R1D, D1D, SIZE( R1D ), R )

            uAF(iNodeX,iX1,iX2,iX3,iAF_T) &
              = Interpolate1D( R1D, T1D, SIZE( R1D ), R )

            uAF(iNodeX,iX1,iX2,iX3,iAF_Ye) &
              = Interpolate1D( R1D, Y1D, SIZE( R1D ), R )

          END DO

          CALL ComputeThermodynamicStates_Primitive_TABLE &
                 ( uPF(:,iX1,iX2,iX3,iPF_D),  uAF(:,iX1,iX2,iX3,iAF_T), &
                   uAF(:,iX1,iX2,iX3,iAF_Ye), uPF(:,iX1,iX2,iX3,iPF_E), &
                   uAF(:,iX1,iX2,iX3,iAF_E),  uPF(:,iX1,iX2,iX3,iPF_Ne) )

          CALL ApplyEquationOfState_TABLE &
                 ( uPF(:,iX1,iX2,iX3,iPF_D ), uAF(:,iX1,iX2,iX3,iAF_T ), &
                   uAF(:,iX1,iX2,iX3,iAF_Ye), uAF(:,iX1,iX2,iX3,iAF_P ), &
                   uAF(:,iX1,iX2,iX3,iAF_S ), uAF(:,iX1,iX2,iX3,iAF_E ), &
                   uAF(:,iX1,iX2,iX3,iAF_Me), uAF(:,iX1,iX2,iX3,iAF_Mp), &
                   uAF(:,iX1,iX2,iX3,iAF_Mn), uAF(:,iX1,iX2,iX3,iAF_Xp), &
                   uAF(:,iX1,iX2,iX3,iAF_Xn), uAF(:,iX1,iX2,iX3,iAF_Xa), &
                   uAF(:,iX1,iX2,iX3,iAF_Xh), uAF(:,iX1,iX2,iX3,iAF_Gm) )

        END DO
      END DO
    END DO

    END ASSOCIATE

  END SUBROUTINE InitializeFluidFields_DeleptonizationWave_Profile

  SUBROUTINE InitializeFluidFields_DeleptonizationWave_Spherical

    ! --- Density Profile ---
    REAL(DP), PARAMETER :: MinD = 1.0d08 * Gram / Centimeter**3
    REAL(DP), PARAMETER :: MaxD = 4.0d14 * Gram / Centimeter**3
    REAL(DP), PARAMETER :: C_D  = 7.5_DP
    REAL(DP), PARAMETER :: R_D  = 2.0d01 * Kilometer
!    REAL(DP), PARAMETER :: H_D  = 1.0d01 * Kilometer
    REAL(DP), PARAMETER :: H_D  = 5.0d00 * Kilometer
    ! --- Temperature Profile ---
    REAL(DP), PARAMETER :: MinT = 5.0d09 * Kelvin
    REAL(DP), PARAMETER :: MaxT = 1.5d11 * Kelvin
    REAL(DP), PARAMETER :: C_T  = 1.0_DP
    REAL(DP), PARAMETER :: R_T  = 2.5d01 * Kilometer
!    REAL(DP), PARAMETER :: H_T  = 2.0d01 * Kilometer
    REAL(DP), PARAMETER :: H_T  = 5.0d01 * Kilometer
    ! --- Electron Fraction Profile ---
    REAL(DP), PARAMETER :: MinY = 2.5d-1
    REAL(DP), PARAMETER :: MaxY = 4.6d-1
    REAL(DP), PARAMETER :: C_Y  = 1.0_DP
    REAL(DP), PARAMETER :: R_Y  = 4.5d01 * Kilometer
!    REAL(DP), PARAMETER :: H_Y  = 1.0d01 * Kilometer
    REAL(DP), PARAMETER :: H_Y  = 5.0d01 * Kilometer

    INTEGER  :: iX1, iX2, iX3, iNodeX
    INTEGER  :: iNodeX1, iNodeX2, iNodeX3
    REAL(DP) :: X1, X2, X3, R

    DO iX3 = iX_B0(3), iX_E0(3)
      DO iX2 = iX_B0(2), iX_E0(2)
        DO iX1 = iX_B0(1), iX_E0(1)

          DO iNodeX = 1, nDOFX

            iNodeX1 = NodeNumberTableX(1,iNodeX)
            iNodeX2 = NodeNumberTableX(2,iNodeX)
            iNodeX3 = NodeNumberTableX(3,iNodeX)

            R = NodeCoordinate( MeshX(1), iX1, iNodeX1 )

            uPF(iNodeX,iX1,iX2,iX3,iPF_D) &
!!$              = Half * ( MaxD * ( One - TANH( (R-R_D)/H_D ) ) &
!!$                         + MinD * ( One - TANH( (R_D-R)/H_D ) ) )
              = MaxD * C_D / ( C_D + ( R / H_D )**4 )

            uAF(iNodeX,iX1,iX2,iX3,iAF_T) &
!!$              = Half * ( MaxT * ( One - TANH( (R-R_T)/H_T ) ) &
!!$                         + MinT * ( One - TANH( (R_T-R)/H_T ) ) )
              = MAX( MaxT * C_T / ( C_T + ( R / H_T )**2 ), 1.0d10 * Kelvin )

            uAF(iNodeX,iX1,iX2,iX3,iAF_Ye) &
!!$              = Half * ( MinY * ( One - TANH( (R-R_Y)/H_Y ) ) &
!!$                         + MaxY * ( One - TANH( (R_Y-R)/H_Y ) ) )
              = MinY * ( One + C_Y / ( C_Y + ( R / H_Y )**(-12) ) )

          END DO

          CALL ComputeThermodynamicStates_Primitive_TABLE &
                 ( uPF(:,iX1,iX2,iX3,iPF_D),  uAF(:,iX1,iX2,iX3,iAF_T), &
                   uAF(:,iX1,iX2,iX3,iAF_Ye), uPF(:,iX1,iX2,iX3,iPF_E), &
                   uAF(:,iX1,iX2,iX3,iAF_E),  uPF(:,iX1,iX2,iX3,iPF_Ne) )

          CALL ApplyEquationOfState_TABLE &
                 ( uPF(:,iX1,iX2,iX3,iPF_D ), uAF(:,iX1,iX2,iX3,iAF_T ), &
                   uAF(:,iX1,iX2,iX3,iAF_Ye), uAF(:,iX1,iX2,iX3,iAF_P ), &
                   uAF(:,iX1,iX2,iX3,iAF_S ), uAF(:,iX1,iX2,iX3,iAF_E ), &
                   uAF(:,iX1,iX2,iX3,iAF_Me), uAF(:,iX1,iX2,iX3,iAF_Mp), &
                   uAF(:,iX1,iX2,iX3,iAF_Mn), uAF(:,iX1,iX2,iX3,iAF_Xp), &
                   uAF(:,iX1,iX2,iX3,iAF_Xn), uAF(:,iX1,iX2,iX3,iAF_Xa), &
                   uAF(:,iX1,iX2,iX3,iAF_Xh), uAF(:,iX1,iX2,iX3,iAF_Gm) )

        END DO
      END DO
    END DO

  END SUBROUTINE InitializeFluidFields_DeleptonizationWave_Spherical

  SUBROUTINE InitializeFluidFields_DeleptonizationWave_Spherical_Profile( FileName )

    CHARACTER(len=*), INTENT(in) :: FileName

    TYPE(ProgenitorType1D) :: P1D
    CHARACTER(LEN=100)                  :: Format1, Format2, Format3
    CHARACTER(LEN=30)                   :: a
    INTEGER                             :: ii, datasize
    REAL(dp), DIMENSION(:,:), ALLOCATABLE :: database

    INTEGER  :: iX1, iX2, iX3, iNodeX
    INTEGER  :: iNodeX1, iNodeX2, iNodeX3
    REAL(DP) :: X1, X2, X3, R

    CALL ReadProgenitor1D( FileName, P1D )

    ASSOCIATE &
      ( R1D => P1D % Radius, &
        D1D => P1D % MassDensity, &
        T1D => P1D % Temperature, &
        Y1D => P1D % ElectronFraction )

    DO iX3 = iX_B0(3), iX_E0(3)
      DO iX2 = iX_B0(2), iX_E0(2)
        DO iX1 = iX_B0(1), iX_E0(1)

          DO iNodeX = 1, nDOFX

            iNodeX1 = NodeNumberTableX(1,iNodeX)
            iNodeX2 = NodeNumberTableX(2,iNodeX)
            iNodeX3 = NodeNumberTableX(3,iNodeX)

            X1 = NodeCoordinate( MeshX(1), iX1, iNodeX1 )

            R = X1

            uPF(iNodeX,iX1,iX2,iX3,iPF_D) &
              = Interpolate1D( R1D, D1D, SIZE( R1D ), R )

            uAF(iNodeX,iX1,iX2,iX3,iAF_T) &
              = Interpolate1D( R1D, T1D, SIZE( R1D ), R )

            uAF(iNodeX,iX1,iX2,iX3,iAF_Ye) &
              = Interpolate1D( R1D, Y1D, SIZE( R1D ), R )

          END DO

          CALL ComputeThermodynamicStates_Primitive_TABLE &
                 ( uPF(:,iX1,iX2,iX3,iPF_D),  uAF(:,iX1,iX2,iX3,iAF_T), &
                   uAF(:,iX1,iX2,iX3,iAF_Ye), uPF(:,iX1,iX2,iX3,iPF_E), &
                   uAF(:,iX1,iX2,iX3,iAF_E),  uPF(:,iX1,iX2,iX3,iPF_Ne) )

          CALL ApplyEquationOfState_TABLE &
                 ( uPF(:,iX1,iX2,iX3,iPF_D ), uAF(:,iX1,iX2,iX3,iAF_T ), &
                   uAF(:,iX1,iX2,iX3,iAF_Ye), uAF(:,iX1,iX2,iX3,iAF_P ), &
                   uAF(:,iX1,iX2,iX3,iAF_S ), uAF(:,iX1,iX2,iX3,iAF_E ), &
                   uAF(:,iX1,iX2,iX3,iAF_Me), uAF(:,iX1,iX2,iX3,iAF_Mp), &
                   uAF(:,iX1,iX2,iX3,iAF_Mn), uAF(:,iX1,iX2,iX3,iAF_Xp), &
                   uAF(:,iX1,iX2,iX3,iAF_Xn), uAF(:,iX1,iX2,iX3,iAF_Xa), &
                   uAF(:,iX1,iX2,iX3,iAF_Xh), uAF(:,iX1,iX2,iX3,iAF_Gm) )

        END DO
      END DO
    END DO

    END ASSOCIATE

  END SUBROUTINE InitializeFluidFields_DeleptonizationWave_Spherical_Profile

  SUBROUTINE InitializeRadiationFields_DeleptonizationWave

    INTEGER  :: iE, iX1, iX2, iX3, iS
    INTEGER  :: iNode, iNodeE
    REAL(DP) :: kT(nDOF)
    REAL(DP) :: Mnu(nDOF), E
    REAL(DP) :: Gm_dd_11(nDOF)
    REAL(DP) :: Gm_dd_22(nDOF)
    REAL(DP) :: Gm_dd_33(nDOF)

    DO iS = 1, nSpecies

      DO iX3 = iX_B0(3), iX_E0(3)
        DO iX2 = iX_B0(2), iX_E0(2)
          DO iX1 = iX_B0(1), iX_E0(1)

            Gm_dd_11 &
              = uGF(NodeNumbersX,iX1,iX2,iX3,iGF_Gm_dd_11)

            Gm_dd_22 &
              = uGF(NodeNumbersX,iX1,iX2,iX3,iGF_Gm_dd_22)

            Gm_dd_33 &
              = uGF(NodeNumbersX,iX1,iX2,iX3,iGF_Gm_dd_33)

            kT = BoltzmannConstant &
                 * uAF(NodeNumbersX,iX1,iX2,iX3,iAF_T)
!!!
            Mnu = uAF(NodeNumbersX,iX1,iX2,iX3,iAF_Me) &
                  + uAF(NodeNumbersX,iX1,iX2,iX3,iAF_Mp) &
                  - uAF(NodeNumbersX,iX1,iX2,iX3,iAF_Mn)

            DO iE = iE_B0, iE_E0

              DO iNode = 1, nDOF

                iNodeE = NodeNumberTable(1,iNode)

                E = NodeCoordinate( MeshE, iE, iNodeE )

                uPR(iNode,iE,iX1,iX2,iX3,iPR_D,iS) &
                  = MAX( 1.0d-32, &
                    1.0d0 / ( EXP( (E-Mnu(iNode))/kT(iNode) ) + 1.0_DP ) )
                  != 1.0d-99

                uPR(iNode,iE,iX1,iX2,iX3,iPR_I1,iS) &
                  = Zero

                uPR(iNode,iE,iX1,iX2,iX3,iPR_I2,iS) &
                  = Zero

                uPR(iNode,iE,iX1,iX2,iX3,iPR_I3,iS) &
                  = Zero

              END DO

              CALL ComputeConserved_TwoMoment &
                     ( uPR(:,iE,iX1,iX2,iX3,iPR_D, iS), &
                       uPR(:,iE,iX1,iX2,iX3,iPR_I1,iS), &
                       uPR(:,iE,iX1,iX2,iX3,iPR_I2,iS), &
                       uPR(:,iE,iX1,iX2,iX3,iPR_I3,iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_N, iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_G1,iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_G2,iS), &
                       uCR(:,iE,iX1,iX2,iX3,iCR_G3,iS), &
                       Gm_dd_11, Gm_dd_22, Gm_dd_33 )

            END DO

          END DO
        END DO
      END DO

    END DO

  END SUBROUTINE InitializeRadiationFields_DeleptonizationWave

  REAL(DP) FUNCTION Interpolate1D( x, y, n, xq )

    INTEGER,                INTENT(in) :: n
    REAL(DP), DIMENSION(n), INTENT(in) :: x, y
    REAL(DP),               INTENT(in) :: xq

    INTEGER :: i

    i = Locate( xq, x, n )

    IF( i == 0 )THEN

      ! --- Extrapolate Left ---

      Interpolate1D &
        = Interpolate1D_Linear( xq, x(1), x(2), y(1), y(2) )

    ELSE IF( i == n )THEN

      ! --- Extrapolate Right ---

      Interpolate1D &
        = Interpolate1D_Linear( xq, x(n-1), x(n), y(n-1), y(n) )

    ELSE

      Interpolate1D &
        = Interpolate1D_Linear( xq, x(i), x(i+1), y(i), y(i+1) )

    END IF

    RETURN

  END FUNCTION Interpolate1D

END MODULE InitializationModule

max_stars_repo_name: srichers/thornado

max_issues_repo_name: srichers/thornado

max_forks_repo_name: srichers/thornado


subroutine read_write

  USE variables
  implicit none

  integer i,j,k,l,m,n,o
  real*8,  dimension(:,:), allocatable  :: Ov, HOv
  real*8,  dimension(:,:,:,:), allocatable  :: EEOv
  character st1,st2,st3,st4
  integer nn,ii

  open(unit = 100, file = 'molecule.dat', status = 'old', action = 'read')
      read(100,*) st1, n_el
      read(100,*) st2, sp_deg
      read(100,*) st3, n_bas
      read(100,*) st4, nuc_repul
      allocate ( Ov(n_bas,n_bas),HOv(n_bas,n_bas),EEOv(n_bas,n_bas,n_bas,n_bas),STAT = AllocateStatus)

     close(100)
  open(unit = 101, file = 'ov_integral.dat', status = 'old', action = 'read')

      do i=1,n_bas
        read(101,*) Ov(i,:)
      end do
    close(101)

  open(unit = 102, file = 'ham_integral.dat', status = 'old', action = 'read')

      do i=1,n_bas
        read(102,*) HOv(i,:)
      end do

     close(102)
  open(unit = 103, file = '2e_integral.dat', status = 'old', action = 'read')
      read(103,*) nn
      do ii=1,nn
        read(103,*) i,j,k,l,EEOv(i,k,j,l)
        EEOv(j,k,i,l) = EEOv(i,k,j,l)
        EEOv(i,l,j,k) = EEOv(i,k,j,l)
        EEOv(j,l,i,k) = EEOv(i,k,j,l)
        EEOv(k,i,l,j) = EEOv(i,k,j,l)
        EEOv(k,j,l,i) = EEOv(i,k,j,l)
        EEOv(l,i,k,j) = EEOv(i,k,j,l)
        EEOv(l,j,k,i) = EEOv(i,k,j,l)
      end do
     close(103)
   call orthogonal(Ov,HOv,EEOv)
!   write(6,*) "Pass Read"
end subroutine read_write


text:
HANSARD REVISE * NUMERO 194
Le jeudi 11 mars 1999
LES TRAVAUX DE LA CHAMBRE
LES COMITES DE LA CHAMBRE
LOI SUR LE SYSTEME DE JUSTICE PENALE POUR LES ADOLESCENTS
Presentation et premiere lecture
LOI INTERDISANT L'EXPORTATION D'EAU HORS DU CANADA
Presentation et premiere lecture
LA LOI SUR LA REVISION DES LIMITES DES CIRCONSCRIPTIONS
Presentation et premiere lecture
LOI SUR LA JOURNEE NATIONALE DE LUTTE CONTRE LA CONDUITE
LES COMITES DE LA CHAMBRE
LOI INSTITUANT UN SYMBOLE DE L'UNITE CANADIENNE
L'Accord multilateral sur l'investissement
LA LOI SUR LES ALLOCATIONS AUX ANCIENS COMBATTANTS
L'hon. Harbance Singh Dhaliwal
L'hon. Harbance Singh Dhaliwal
LA DEPUTEE D'EDMONTON-NORD
LES SERVICES A L'ENFANCE ET A LA FAMILLE
LA SEMAINE NATIONALE DE LA SECURITE A LA FERME
LA DEPUTEE D'EDMONTON-NORD
Mme Rose-Marie Ur
M. Gurbax Singh Malhi
LA DEPUTEE D'EDMONTON-NORD
M. Svend J. Robinson
Mme Caroline St-Hilaire
LA DEPUTEE D'EDMONTON-NORD
L'hon. Arthur C. Eggleton
L'hon. Arthur C. Eggleton
L'hon. Pierre S. Pettigrew
L'hon. Pierre S. Pettigrew
Mme Karen Kraft Sloan
L'hon. Arthur C. Eggleton
Mme Diane St-Jacques
L'hon. Arthur C. Eggleton
Mme Jocelyne Girard-Bujold
Mme Diane St-Jacques
Certains propos tenus au cours de la periode des questions
LES TRAVAUX DE LA CHAMBRE
M. Rey D. Pagtakhan
M. Svend J. Robinson
LES TRAVAUX DE LA CHAMBRE
Demande et report des votes
Le Programme d'adaptation pour les travailleurs ages
Mme Jocelyne Girard-Bujold
HANSARD REVISE * NUMERO 194
Le jeudi 11 mars 1999
La seance est ouverte a 10 heures.
LES TRAVAUX DE LA CHAMBRE
L'hon. Don Boudria (leader du gouvernement a la Chambre des communes, Lib.):
Monsieur le President, j'invoque le Reglement.
Les leaders de tous les partis a la Chambre ont discute de la motion suivante:
LES COMITES DE LA CHAMBRE
M. Charles Hubbard (Miramichi, Lib.):
Monsieur le President, j'invoque le Reglement.
Le president suppleant (M. McClelland):
Le document est depose.
M. Randy White (Langley-Abbotsford, Ref.):
Le president suppleant (M. McClelland):
A ce stade-ci, on ne fait que deposer le document.
A ce moment-la, on l'etudiera.
Monsieur le President, je voudrais des eclaircissements la-dessus.
M. Bob Kilger (Stormont-Dundas, Lib.):
Le president suppleant (M. McClelland):
Il se peut que nous mettions la charrue devant les boeufs.
Nous en sommes seulement au depot du rapport.
La Chambre n'est saisie d'aucune motion.
Le president suppleant (M. McClelland):
Nous allons nous arreter ici et passer aux projets de loi d'initiative ministerielle.
LOI SUR LE SYSTEME DE JUSTICE PENALE POUR LES ADOLESCENTS
LOI INTERDISANT L'EXPORTATION D'EAU HORS DU CANADA
L'eau sera certainement le facteur strategique du siecle prochain et des siecles a venir.
Les Canadiens desirent plus que toute chose proteger leur eau contre l'exportation.
C'est ce que le projet de loi vise a faire.
LA LOI SUR LA REVISION DES LIMITES DES CIRCONSCRIPTIONS ELECTORALES
LOI SUR LA JOURNEE NATIONALE DE LUTTE CONTRE LA CONDUITE AVEC DES FACULTES AFFAIBLIES
Ce projet de loi devrait etre surnomme la loi Cindy.
M. Randy White (Langley-Abbotsford, Ref.):
Le president suppleant (M. McClelland):
LES COMITES DE LA CHAMBRE
LOI INSTITUANT UN SYMBOLE DE L'UNITE CANADIENNE
M. Lynn Myers (Waterloo-Wellington, Lib.):
Le president suppleant (M. McClelland):
Y a-t-il consentement unanime?
(Les ordres sont annules, et les projets de loi sont retires.)
M. Paul Szabo (Mississauga-Sud, Lib.):
M. Randy White (Langley-Abbotsford, Ref.):
D'autres noms vont s'ajouter a cette petition.
Je demande instamment a la Chambre d'y donner suite.
M. Peter Goldring (Edmonton-Est, Ref.):
M. Stephane Bergeron (Vercheres-Les-Patriotes, BQ):
L'ACCORD MULTILATERAL SUR L'INVESTISSEMENT
M. John Solomon (Regina-Lumsden-Lake Centre, NPD):
M. Richard Marceau (Charlesbourg, BQ):
Cette petition se joint a celles presentees par mes autres collegues du Bloc quebecois.
Monsieur le President, je suggere que toutes les questions soient reservees.
Le president suppleant (M. McClelland):
LA LOI SUR LES ALLOCATIONS AUX ANCIENS COMBATTANTS
Le president suppleant (M. McClelland):
Quand le projet de loi sera-t-il lu pour la troisieme fois?
M. Bob Wood (secretaire parlementaire du ministre des Anciens combattants, Lib.):
C'est precisement l'objet du projet de loi.
Il ne fera pas la manchette.
Que fait le projet de loi C-61 pour les anciens combattants?
Plus de 35 000 veuves sont ainsi visees.
J'exhorte tous les deputes a appuyer ce projet de loi.
M. Peter Goldring (Edmonton-Est, Ref.):
Ces braves Canadiens et Terre-Neuviens ont paye cher leur devouement envers le service.
La premiere victime de guerre du service appartenait a la marine marchande.
La marine marchande du Canada etait tres peu imposante aux premiers jours de la guerre.
A cette epoque-la, elle ne comprenait que 38 navires oceaniques.
Une nuit de pleine lune, calme et claire, provoquait de l'insomnie.
En 1942, la flotte de sous-marins allemands est passee de 91 a 212 unites.
Cela a fait empirer la situation des navires marchands.
La bataille de l'Atlantique s'est poursuivie implacablement durant toute l'annee 1942.
Les sous-marins allemands infestaient les eaux canadiennes.
Plusieurs navires ont ete coules dans le Saint-Laurent.
La contribution canadienne est devenue d'autant plus importante a mesure que durait la guerre.
Ce fut vraiment une guerre d'approvisionnement en materiel.
Tout ce materiel a ete transporte par notre marine marchande.
Certains etaient meme des laquiers.
On trouvait une aussi grande variete dans les equipages.
Il n'y avait donc ni navire type, ni marin type.
Beaucoup avaient des familles, des enfants et des petits-enfants.
Des 1 500 marins marchands qui sont morts, huit etaient des femmes.
Elle transportait le materiel necessaire a la reconstruction et au reapprovisionnement de l'Europe.
La marine marchande a donc perdu en tout 72 navires.
De nombreuses personnes devouees ont travaille pour faire avancer la cause de la marine marchande.
On se souviendra de leur contribution a cette longue bataille pour l'egalite.
Je suis heureux que nous puissions adopter ce projet de loi sans retard deraisonnable.
Il ne reste que la question de l'indemnisation a regler.
M. Maurice Godin (Chateauguay, BQ):
Ce projet de loi repond a des demandes de certaines associations.
Pourquoi ce gouvernement nie-t-il aujourd'hui ce qu'il souhaitait hier?
Je suis personnellement au courant des changements.
Cela se passait a l'automne de 1991.
On s'apercoit que leurs revendications avancent a petits pas.
Cependant, il faut tenir compte que la moyenne d'age est de 75 ans.
Il manque 50 annees de benefices sociaux auxquels ces marins ont ete prives.
Ils ont souffert toute leur vie de ce refus.
Pourtant, ces marins etaient sur la ligne de feu.
C'etait des conditions de travail horribles et de lourdes pertes de vies humaines.
Pourtant, aujourd'hui, il demeure de glace devant cette grande souffrance profonde et surtout mentale.
Seul le gouvernement liberal demeure intraitable.
Outre les avantages sociaux et une pension d'invalidite, il faut un regime d'indemnisation.
Le gouvernement a-t-il les sommes necessaires?
Au total, 887 960 424 $ non utilises.
La demande de la marine marchande est d'environ 40 millions de dollars.
Ce projet de loi aurait pu regenerer ces sommes d'argent non depensees.
Et j'en passe.
J'ai assiste a ces pelerinages.
Les marins marchands ne manifestaient aucune amertume.
Pour quelles raisons refuse-t-on encore, apres 50 ans, de reconnaitre leurs droits acquis?
Cela lui donne une belle image.
Les Europeens et les Asiatiques l'ont subie et s'en souviennent.
La plupart du temps, une poignee de spectateurs du milieu assistent a ces ceremonies.
Pourquoi ne pas instaurer une veritable Journee nationale du souvenir au Canada?
Ce sont les Canadiennes et les Canadiens que nous devons informer et faire participer.
Faire une veritable journee de respect, de remerciement et de commemoration.
En comite, j'ai presente un amendement qui fut juge non recevable.
Pourtant, cet amendement, incorpore a cette loi, aurait solutionne ce probleme une fois pour toutes.
Il se lisait comme suit:
M. Gordon Earle (Halifax-Ouest, NPD):
Ces braves Canadiens ont joue un role nevralgique dans les efforts de guerre du Canada.
Quel sort a ete reserve aux marins marchands canadiens a leur retour au pays?
Ce paiement a ete promis en decembre dernier.
Cela me semble au moins un effort d'en arriver a un reglement juste.
Que dire des anciens combattants autochtones du Canada?
Le gouvernement liberal a echoue miserablement la ou bien d'autres gouvernements ont reussi.
Mme Elsie Wayne (Saint John, PC):
Fort heureusement, mes deux freres sont rentres au pays sains et saufs.
Je n'oublierai jamais ce jour-la, je devais avoir sept ans.
Nous allons a la gare chercher les garcons qui reviennent.
Bref, les anciens combattants sont tres chers a mon coeur.
Comment vous sentiriez-vous?
Vous vous sentiriez blesse, rejete.
Je n'etais pas la.
Ils viennent parfois me voir a mon bureau, les larmes aux yeux.
Leurs epouses aussi communiquent avec moi.
Ils se battent depuis 54 ans pour avoir l'egalite.
Je crois personnellement que la plupart des Canadiens repondraient oui a cette question.
Pourquoi seraient-ils de cet avis?
Parce que de nos jours, les gens sont bien informes.
Ils savent que ces anciens combattants ne recevaient pas un traitement tres eleve.
On peut voir les hommes dans l'eau, dont certains sont en train de mourir.
Un marin sur huit a peri en mer.
Je suis tres heureuse qu'ils aient pu beneficier de ces avantages.
Mais qu'a-t-on offert aux membres de la marine marchande?
On a prefere les occuper a nettoyer les mers.
Ces hommes sont tous tres ages maintenant.
Ils ne peuvent plus pelleter la neige qui tombe dans leur entree.
Ils ne peuvent plus s'occuper du nettoyage courant.
Ils ont du mal a accomplir les taches menageres quotidiennes.
Ils ont besoin de cette aide ponctuelle.
Nous ne voulons pas qu'ils se retrouvent seuls.
Je suis fiere d'etre Canadienne.
Je ne crois pas que personne en doute.
Il y a de l'espoir toutefois.
Je tiens a etre bien claire.
Le maximum prevu est de 30 000 $ je crois.
Un tres petit nombre d'entre eux y auraient droit, cependant.
Il en reste environ 2 000.
L'argent est la.
Nous pouvons rendre leur dignite a ces hommes.
Je sais qu'un grand nombre de deputes ministeriels sont d'accord avec moi.
C'est ce dont nous avons besoin.
C'est le genre de dialogue qu'il faut pour se comprendre.
J'en profite aujourd'hui pour le remercier sincerement.
Le president suppleant (M. McClelland):
La Chambre est-elle prete a se prononcer?
Le president suppleant (M. McClelland):
Plait-il a la Chambre d'adopter la motion?
Le president suppleant (M. McClelland):
Je declare la motion adoptee.
(Le projet de loi est lu pour la troisieme fois et adopte.)
Mme Carolyn Parrish (secretaire parlementaire du ministre des Travaux publics et des Services gouvernementaux, Lib.):
Les realisations de la SCHL ont grandement profite a notre pays.
Bon nombre de familles canadiennes et de localites canadiennes ont recu de l'aide.
Les fonds aident a reparer les logis insalubres et dangereux.
La mesure a l'etude presente un triple avantage.
Troisiemement, la SCHL sera en mesure d'offrir un meilleur service a tous les Canadiens.
Il etait accompagne de 30 representants de l'industrie du logement.
Les deputes de Kelowna et de Quebec ont egalement participe a cette mission commerciale.
J'ai mentionne plus tot que ces modifications refletent les priorites de notre gouvernement.
Cela se voit de plusieurs facons.
Cependant, nous voulons faire davantage.
Des millions de familles canadiennes ont achete leur premiere maison grace au programme d'assurance.
En fait, un acheteur de maison sur trois a recu cette aide de la SCHL.
Ainsi, elle a montre au gouvernement une nouvelle facon d'assurer les services.
M. Werner Schmidt (Kelowna, Ref.):
Il apporte aussi un amendement correlatif a une autre loi.
Quels en seront les couts financiers?
Quels en seront les couts economiques?
Qu'en est-il de la qualite de vie?
Voila les questions sur lesquelles je vais m'arreter.
N'est-ce pas une facon interessante de dire qu'on delaisse cette affaire?
Le gouvernement n'a jamais consulte les provinces.
Je reviendrai la-dessus tout a l'heure.
Ou en est la SCHL a l'heure actuelle?
Je voudrais entrer ici dans les details.
Ce qui suit est aussi tire du plan d'entreprise.
Pour ce qui est de l'assurance-pret hypothecaire, voici ce que dit le plan:
C'est un merveilleux programme.
Voila ce qu'il fait.
Il ne pouvait depasser les 150 milliards de dollars.
Autrement dit, le plafond est passe a 200 milliards de dollars.
Il convient de souligner que le projet de loi ne modifie pas cela.
C'est une evolution interessante.
La societe fait la pluie et le beau temps.
Cette diminution en 1997 est attribuable a une augmentation des depenses liees aux reglements.
En 1996, la societe a enregistre des benefices.
Cela s'applique aux nouvelles activites commerciales d'assurance hypothecaire entreprises apres 1996.
Ce qui commence a se produire a la SCHL est tres interessant.
Il ne s'agit plus d'une simple societe d'Etat.
Legalement, c'en est une, mais elle entre en concurrence directe avec les institutions financieres.
J'expliquerai cela un peu plus tard.
C'est un secteur.
La SCHL exerce egalement son activite dans le secteur des titres hypothecaires.
Passons maintenant a d'autres secteurs.
La SCHL veut exercer son activite dans le secteur de l'aide au logement.
En fait, c'est un peu moins de deux milliards de dollars.
C'est tres interessant.
Ce theme prendra de l'importance a mesure que nous avancerons.
Ce qui se passe n'est pas sans importance.
Aujourd'hui, la SCHL est aussi active sur la scene internationale.
J'etais l'un de ces participants a cette entreprise commerciale.
Ce fut une bonne mission.
La mission commerciale n'a beneficie d'aucune subvention gouvernementale.
J'ai vraiment ete impressionne par l'organisation de cette mission.
C'est une question tout a fait differente.
Il y a aussi d'autres initiatives.
Le compte d'entreprise est un autre secteur.
C'est la situation qui existe maintenant.
Cette loi devait permettre de construire, jusqu'en 1938, presque 5 000 logements.
C'etait la le debut.
C'est la raison pour laquelle la societe a ete creee.
Pourquoi l'institution financiere n'augmenterait-elle pas le montant d'argent accessible?
Il est exprime tres brievement.
En 1946, le gouvernement federal a cree la Societe centrale d'hypotheques et de logement.
Il existe egalement un partenariat historique entre le gouvernement federal et les provinces.
Le ministre avait annonce que ce programme allait prendre fin.
L'administration de ces programmes ete dans la plupart des cas confiee aux provinces.
C'est vraiment la question.
On leur a donne l'initiative dans ce domaine.
Elles sont tres capables de le faire.
Il concerne les cooperatives de logement.
Il a dit qu'il voulait se decharger de cette responsabilite.
Finalement, les pressions ont ete telles que le ministre a change d'avis.
Il a dit qu'il continuerait d'assumer la responsabilite des cooperatives de logement federales.
Elle donne un sentiment de propriete.
Il existe deux types de cooperatives d'habitation.
C'est ce que nous souhaitons.
Je sais que leur nombre est partage a peu pres a parts egales.
An action plan for Toronto.
Voila qui devient plutot precis, mais il en est de meme dans les faits.
On a enregistre une augmentation considerable de leur nombre.
Certains refuges sont maintenant occupes en grande partie par des familles.
Que devient le volet sante dans tout cela?
On compte quelque 110 recommandations a la fin de ce rapport.
Nous devons nous attaquer a ce probleme.
Un gouvernement ne peut pas y arriver seul.
En feuilletant ce plan general, on note avec interet une citation de la SCHL.
Dans le plan de 1998-2002, la SCHL affirme qu'elle:
De combien de perspectives differentes pouvons-nous examiner ce probleme?
On pourrait croire que la ville de Toronto y travaillait toute seule de son cote.
Ce groupe de travail a collabore avec toutes sortes de specialistes.
A-t-il obtenu de l'aide simplement de Toronto?
Il l'a obtenue recue de villes canadiennes comme Calgary, Montreal, Ottawa, Vancouver et Winnipeg.
Ce groupe de travail n'est pas forme d'amateurs.
Ce ne sont pas des gens qui se fichent de la politique.
La SCHL dit maintenant qu'elle va organiser un forum.
Nous avons les renseignements dont nous avons besoin.
Le groupe de travail sait tout cela.
En 1935, cette intervention s'etait revelee tres utile.
La participation de multiples ordres de gouvernement mene essentiellement a l'inefficacite.
Les intervenants en viennent aux tiraillements.
Nous pourrions nous en passer.
En fait, la situation pourrait facilement etre simplifiee.
Celui-ci prefere les solutions rapides et ponctuelles.
Passons aux actes, maintenant.
Ce projet de loi ne s'attaque pas aux problemes de ce genre.
Il ne fait que legaliser l'etat actuel des choses.
Je crois que la SCHL s'est ecartee du chemin a un autre egard.
On peut aussi lire ce qui suit:
Voici les tactiques cles:
La SCHL prevoit ameliorer sa gamme de produits d'assurance hypothecaire.
La Societe cherchera a creer des blocs de creances hypothecaires de plusieurs preteurs.
On elabore actuellement une formule de paiement annuel au gouvernement federal.
C'est exactement ce que fait cette loi.
Elle entre tout simplement dans ce secteur en concurrence directe.
Cela me parait mauvais en principe.
Cela n'est pas du tout equitable.
Qu'est-ce que c'est qu'un pret hypothecaire inverse?
Peut-etre que ce marche grandira.
C'est une autre affaire.
M. Ghislain Lebel (Chambly, BQ):
La presidente suppleante (Mme Thibeault):
Je m'excuse d'interrompre l'honorable depute.
Le parti que je represente ici, a la Chambre, reconnait tout cela.
Donc, il ne s'agit pas d'un cadeau venu de Dieu.
Le createur ne nous a pas envoye un cheque pour mettre la SCHL en place.
Bien des choses ont change depuis son arrivee, a ce semeur de zizanie.
C'est un semeur de zizanie.
Il est permis de se poser des questions.
Ce n'est pas nous.
Au Quebec, on n'a jamais ete la.
On le voit dans plusieurs domaines.
L'esprit du projet de loi C-66 commence a refleter cela.
C'est pour cela qu'il s'est donne une Societe d'habitation du Quebec.
On s'aligne vers une politique nouvelle et, comme toujours, ce gouvernement est inconsequent.
C'est arrive dans mon comte.
C'est frequent, ce n'est pas le premier cas.
On a appris a vivre avec cela.
On vient s'immiscer dans un domaine qui n'est pas le notre.
Comme depute du Bloc quebecois, je ne peux pas accepter cela.
Une fois de plus, on pietine une competence qui releve des provinces.
Les depenses gouvernementales comme telles n'ont pas tellement baisse depuis les cinq dernieres annees.
On ne se laisse pas leurrer.
Pour ma part, ce projet de loi ne m'interesse pas du tout.
La Societe canadienne d'hypotheques et de logement prete le reste.
Cela n'a pas grand bon sens.
Ce n'est pas le gouvernement qui perd.
Je termine donc la-dessus.
M. Real Menard (Hochelaga-Maisonneuve, BQ):
Cela etant dit, nous aurions aime etre en faveur de ce projet de loi.
Parce que tous, autant que nous sommes, nous sommes des gens positifs.
Je pense que c'est connu.
Mais il y a une difficulte.
J'aimerais que quelqu'un me l'explique en cours de debat.
Ce qui pose probleme, ce sont deux articles.
Celui-ci m'a repondu: Oui, c'est quelque chose qui est possible.
Il prevoit, et je cite:
J'aurai l'occasion d'y revenir.
Je vois que le ministre des Affaires intergouvernementales opine du bonnet.
C'est une bonne nouvelle.
On attendait cela depuis longtemps.
Sauf que les montants en question sont completement ridicules.
Or, quel est le poids demographique du Quebec a l'interieur de la federation canadienne?
Or, nous recevons 18,7 p. 100 des depenses federales en matiere d'habitation.
Il y a un probleme, et cette tendance, tous les gouvernements l'ont deploree.
La presidente suppleante (Mme Thibeault):
L'honorable depute demande la permission de la Chambre pour deposer un document.
Y a-t-il consentement unanime?
Madame la Presidente, je pense que les deputes vont en sortir enrichis.
Faites-le circuler de ce cote-la, par exemple.
J'ai la repartition.
Ce n'est que le bon sens qui dicte cela.
Mais il ne veut pas le faire avec des pertes.
Ce n'est pas rien.
Il doit s'inscrire dans une intervention de securite du revenu.
Il doit s'inscrire dans une intervention de lutte a la pauvrete.
C'est evidemment le gouvernement du Quebec.
Il y a un travail de proximite avec les citoyens.
Je connais sa preoccupation pour ce probleme.
Qu'est-ce qu'une habitation a loyer modique?
C'est un projet de loi qui envoie un signal contradictoire.
Nous souhaitons que tout cela soit inscrit dans le projet de loi.
Je demande aux deputes ministeriels de s'inspirer du paragraphe 88(2).
C'est la condition qui nous permettrait d'appuyer le projet de loi.
M. Preston Manning (chef de l'opposition, Ref.):
C'est la que la prevention et la readaptation peuvent se faire.
Les Canadiens exigent un changement pour les jeunes de 10 a 15 ans.
M. Chuck Cadman (Surrey-Nord, Ref.):
Le gouvernement a refuse en invoquant des motifs economiques et politiques.
Quelles considerations politiques ont amene la ministre a ne pas respecter les souhaits des Canadiens?
M. Antoine Dube (Levis-et-Chutes-de-la-Chaudiere, BQ):
Il faut frapper a la racine du probleme.
Monsieur le President, rien de tout cela.
L'hon. Arthur C. Eggleton (ministre de la Defense nationale, Lib.):
M. John Reynolds (West Vancouver-Sunshine Coast, Ref.):
Monsieur le President, ma question s'adresse a la ministre de la Justice.
Le depute de West Vancouver-Sunshine Coast a la parole.
M. John Reynolds (West Vancouver-Sunshine Coast, Ref.):
L'hon. Anne McLellan (ministre de la Justice et procureur general du Canada, Lib.):
Monsieur le President, le droit penal est le meme.
Je pense que notre reforme est orientee dans le sens de l'interet des femmes.
Mme Christiane Gagnon (Quebec, BQ):
Pour celles qui travaillent moins de 15 heures certaines semaines, manifestement, c'est le cas.
M. Monte Solberg (Medicine Hat, Ref.):
Mme Elinor Caplan (secretaire parlementaire du ministre de la Sante, Lib.):
M. Grant Hill (Macleod, Ref.):
J'ai l'arrete ici devant moi.
Pourquoi le gouvernement refuse-t-il de l'admettre?
Le premier ministre a fait construire une nouvelle residence.
L'hon. Anne McLellan (ministre de la Justice et procureur general du Canada, Lib.):
Je prie la ministre de la Justice de repondre a la question.
L'hon. Anne McLellan:
Mme Karen Kraft Sloan (York-Nord, Lib.):
L'hon. Arthur C. Eggleton (ministre de la Defense nationale, Lib.):
Nous faisons et nous continuerons de faire tout ce que nous pouvons.
Monsieur le President, le depute fait erreur sur toute la ligne.
Le Canada est toujours pret a se defendre.
Monsieur le President, ma question s'adresse au secretaire d'Etat a l'Agriculture.
C'est ce que nous faisons actuellement.
C'est pourquoi nous travaillons si energiquement avec les gouvernements provinciaux et territoriaux.
L'hon. Herb Gray (vice-premier ministre, Lib.):
Monsieur le President, nos priorites sont les bonnes.
C'est le depute qui est dans l'erreur.
Nous n'avons pas a nous excuser de defendre le Canada.
Le depute devrait se reveiller et savoir qui sont ses amis.
LE BOGUE DE L'AN 2000
M. Sarkis Assadourian (Brampton-Centre, Lib.):
L'hon. Christine Stewart (ministre de l'Environnement, Lib.):
Nous avons pris les mesures qui s'imposaient.
Mme Libby Davies (Vancouver-Est, NPD):
L'hon. Herb Gray (vice-premier ministre, Lib.):
M. Murray Calder (Dufferin-Peel-Wellington-Grey, Lib.):
Monsieur le President, ma question s'adresse au ministre du Commerce international.
Nous devrions tous etre tres fiers de le faire savoir au monde entier.
Mme Diane St-Jacques (Shefford, PC):
M. Gurmant Grewal (Surrey-Centre, Ref.):
Le debat sur ce projet de loi prendra fin demain.
La Chambre reprend l'etude de la motion:
Toutes les formes de logement sont liees.
Dans ma circonscription, nous avons des cooperatives d'habitation federales.
Reglera-t-il tous les problemes?
Reglera-t-il completement le probleme des sans-abri?
Mme Libby Davies (Vancouver-Est, NPD):
C'est une vraie honte.
Il est garanti dans la Declaration universelle des droits de l'homme.
Nous n'avons pas vu un sou pour les logements sociaux.
C'est une honte!
M. Paul Szabo (Mississauga-Sud, Lib.):
Cependant, supposons qu'ils existent et qu'ils soient adequats.
Ils participaient a des programmes pour les jeunes et avaient droit a du counselling.
Je constate cela sans cesse.
Je ne suis pas necessairement en desaccord avec ce que le depute dit.
Les autres programmes sont tous essentiels.
Il est impossible d'aller travailler.
Ils sont enclins a exagerer.
J'espere que oui.
Voila la question que je voudrais poser au depute.
C'est pour cela qu'on voit plus de gens dans les rues.
Or, je suis d'avis que le logement est une question de sante.
J'invite le depute a le faire.
M. Reed Elley (Nanaimo-Cowichan, Ref.):
La question du logement s'etend a une foule de sujets.
Voila qui donne une toute autre signification a l'expression maison de ses reves.
J'apprecie les observations du depute et je vais y reflechir.
M. Werner Schmidt (Kelowna, Ref.):
Si on adopte une solution unique pour tous, est-ce que cela aide vraiment?
C'est l'attitude des gouvernements, de nos jours.
M. Pat Martin (Winnipeg-Centre, NPD):
Or, on ne construit plus de ces logements sociaux.
C'est la un autre point sur lequel nous manquons le coche.
En tant que charpentier, j'ai construit beaucoup de maisons.
Je vous lis un court passage de la page titre du rapport.
J'ai deja souleve la question.
Il y a un groupe a Winnipeg appele Rossbrooke House.
Cette maison est dirigee par deux religieuses catholiques, soeur Leslie et soeur Bernadette.
Elles font un excellent travail.
Je connais l'importance de l'industrie.
Personne ne faisait de cadeaux.
La dette est remboursee en totalite.
Il s'est rendu compte qu'au fil des ans il l'avait neglige.
Les partisans du libre-echange n'ont pas fourni les unites de logement necessaires.
De 1995-1996 a nos jours, les choses ont encore empire.
Les sans-abri sont des personnes comme les autres.
Nous, ici, au Canada, que fait-on pour les humains qui ont plus de difficulte?
Il a essaye de presenter les choses a sa facon.
Ils n'exigeaient meme pas necessairement des subventions.
Voila ce qui rendrait les projets realisables.
Je voudrais parler de quelques-uns de ces amendements.
Ces Canadiens vivent dans des conditions tres inconfortables.
Elles sont privees de la commodite qu'offre un four a micro-ondes, par exemple.
Leurs enfants ont faim.
Ce serait une bonne blague si ce n'etait pas vrai.
Environ 200 000 Canadiens sont consideres comme etant des sans-abri.
Cela devrait avoir pour effet de stimuler l'emploi au Canada.
Personnellement, je n'en suis pas certain, mais je suis pret a l'accepter.
En parlant d'efficacite, le projet de loi est muet a propos de son application.
Je n'y vois aucune disposition concernant la facon dont il sera applique.
Nous avons besoin de connaitre plus de details a propos du projet de loi.
M. Svend J. Robinson (Burnaby-Douglas, NPD):
J'ai une question a poser au depute reformiste.
Lorsqu'ils grandissent, ils se demandent s'ils trouveront du travail.
LA LOI NATIONALE SUR L'HABITATION
M. Gurmant Grewal (Surrey-Centre, Ref.):
C'est ca la solution.
Mais pour acheter une maison, il faut avoir un emploi.
M. Gordon Earle (Halifax-Ouest, NPD):
Cela m'inquiete beaucoup aussi.
Il y a aussi la commercialisation de l'assurance hypothecaire.
Ce genre de situation se produit a la grandeur du pays.
Plus de 100 000 Canadiens sont des itinerants.
Le gouvernement, lui, en fait.
Il a supprime toute aide financiere aux logements sociaux.
M. Keith Martin (Esquimalt-Juan de Fuca, Ref.):
C'est un grave probleme dans notre pays aujourd'hui.
Le plus triste dans tout cela, c'est ce que nous perdons.
Mais il n'y a pas que celles-la.
Elle n'a pas eu la vie facile, passant d'un refuge a un autre.
Un de ces hommes vit sous un pont depuis 17 mois.
Nous pouvons certainement faire mieux.
Je veux parler des gens des premieres nations.
C'est tout ce qu'elles possedaient.
M. John Bryden (Wentworth-Burlington, Lib.):
Le projet de loi est une amorce de solution.
M. Keith Martin (Esquimalt-Juan de Fuca, Ref.):
Cela s'est avere encore et encore.
LE PROGRAMME D'ADAPTATION POUR LES TRAVAILLEURS AGES
Mme Jocelyne Girard-Bujold (Jonquiere, BQ):
La reponse du ministre ne m'a pas satisfaite du tout.
Le monde du travail evolue rapidement.
Il est en pleine mutation et il exige de nouvelles connaissances et de nouvelles competences.
Le marche du travail ne regle pas, de lui-meme, les inegalites et les injustices.
Le ministre lui-meme l'admet dans son livre.
L'Etat doit corriger les iniquites du marche du travail.
On hesite a leur donner un emploi a cause de leur age.
Ils ont un long chemin a parcourir pour se recycler.
On ne leur ouvre pas toutes grandes les portes du marche du travail.
Monsieur le President, je suis heureux de repondre au nom du ministre.
Le gouvernement entend bien aider les travailleurs sans emploi, y compris les travailleurs ages.
Ainsi, nos efforts visent a aider ces travailleurs ages qui ont besoin de notre aide.
Je remercie la deputee de sa question et de sa preoccupation.
La motion portant que la Chambre s'ajourne maintenant est reputee adoptee.
(La seance est levee a 17 h 37.)

ext: f


text:
      subroutine tred2(nm,n,a,d,e,z)
c
      integer i,j,k,l,n,ii,nm,jp1
      double precision a(nm,n),d(n),e(n),z(nm,n)
      double precision f,g,h,hh,scale
c
c     this subroutine is a translation of the algol procedure tred2,
c     num. math. 11, 181-195(1968) by martin, reinsch, and wilkinson.
c     handbook for auto. comp., vol.ii-linear algebra, 212-226(1971).
c
c     this subroutine reduces a real symmetric matrix to a
c     symmetric tridiagonal matrix using and accumulating
c     orthogonal similarity transformations.
c
c     on input
c
c        nm must be set to the row dimension of two-dimensional
c          array parameters as declared in the calling program
c          dimension statement.
c
c        n is the order of the matrix.
c
c        a contains the real symmetric input matrix.  only the
c          lower triangle of the matrix need be supplied.
c
c     on output
c
c        d contains the diagonal elements of the tridiagonal matrix.
c
c        e contains the subdiagonal elements of the tridiagonal
c          matrix in its last n-1 positions.  e(1) is set to zero.
c
c        z contains the orthogonal transformation matrix
c          produced in the reduction.
c
c        a and z may coincide.  if distinct, a is unaltered.
c
c     questions and comments should be directed to burton s. garbow,
c     mathematics and computer science div, argonne national laboratory
c
c     this version dated august 1983.
c
c     ------------------------------------------------------------------
c
      do 100 i = 1, n
c
         do 80 j = i, n
   80    z(j,i) = a(j,i)
c
         d(i) = a(n,i)
  100 continue
c
      if (n .eq. 1) go to 510
c     .......... for i=n step -1 until 2 do -- ..........
      do 300 ii = 2, n
         i = n + 2 - ii
         l = i - 1
         h = 0.0d0
         scale = 0.0d0
         if (l .lt. 2) go to 130
c     .......... scale row (algol tol then not needed) ..........
         do 120 k = 1, l
  120    scale = scale + dabs(d(k))
c
         if (scale .ne. 0.0d0) go to 140
  130    e(i) = d(l)
c
         do 135 j = 1, l
            d(j) = z(l,j)
            z(i,j) = 0.0d0
            z(j,i) = 0.0d0
  135    continue
c
         go to 290
c
  140    do 150 k = 1, l
            d(k) = d(k) / scale
            h = h + d(k) * d(k)
  150    continue
c
         f = d(l)
         g = -dsign(dsqrt(h),f)
         e(i) = scale * g
         h = h - f * g
         d(l) = f - g
c     .......... form a*u ..........
         do 170 j = 1, l
  170    e(j) = 0.0d0
c
         do 240 j = 1, l
            f = d(j)
            z(j,i) = f
            g = e(j) + z(j,j) * f
            jp1 = j + 1
            if (l .lt. jp1) go to 220
c
            do 200 k = jp1, l
               g = g + z(k,j) * d(k)
               e(k) = e(k) + z(k,j) * f
  200       continue
c
  220       e(j) = g
  240    continue
c     .......... form p ..........
         f = 0.0d0
c
         do 245 j = 1, l
            e(j) = e(j) / h
            f = f + e(j) * d(j)
  245    continue
c
         hh = f / (h + h)
c     .......... form q ..........
         do 250 j = 1, l
  250    e(j) = e(j) - hh * d(j)
c     .......... form reduced a ..........
         do 280 j = 1, l
            f = d(j)
            g = e(j)
c
            do 260 k = j, l
  260       z(k,j) = z(k,j) - f * e(k) - g * d(k)
c
            d(j) = z(l,j)
            z(i,j) = 0.0d0
  280    continue
c
  290    d(i) = h
  300 continue
c     .......... accumulation of transformation matrices ..........
      do 500 i = 2, n
         l = i - 1
         z(n,l) = z(l,l)
         z(l,l) = 1.0d0
         h = d(i)
         if (h .eq. 0.0d0) go to 380
c
         do 330 k = 1, l
  330    d(k) = z(k,i) / h
c
         do 360 j = 1, l
            g = 0.0d0
c
            do 340 k = 1, l
  340       g = g + z(k,i) * z(k,j)
c
            do 360 k = 1, l
               z(k,j) = z(k,j) - g * d(k)
  360    continue
c
  380    do 400 k = 1, l
  400    z(k,i) = 0.0d0
c
  500 continue
c
  510 do 520 i = 1, n
         d(i) = z(n,i)
         z(n,i) = 0.0d0
  520 continue
c
      z(n,n) = 1.0d0
      e(1) = 0.0d0
      return
      end


text:
program double_loop_1

    integer,parameter :: n=139

    real*8 A(n, n)

    !$xev loop_tag
    do i=1,n
        do j=1,n
            A(i,j)=A(i-1,j)+A(i,j-1)
        end do
    end do

end program double_loop_1

max_stars_repo_name: xevolver/xevxml

max_issues_repo_name: xevolver/xevxml

max_forks_repo_name: xevolver/xevxml


      PROGRAM LA_SGBSVX_EXAMPLE

!  -- LAPACK95 EXAMPLE DRIVER ROUTINE (VERSION 1.0) --
!     UNI-C, DENMARK
!     DECEMBER, 1999
!
!  .. "Use Statements"
      USE LA_PRECISION, ONLY: WP => SP
      USE F95_LAPACK, ONLY: LA_GBSVX
!  .. "Implicit Statement" ..
      IMPLICIT NONE
!  .. "Local Scalars" ..
      INTEGER :: K, KL, KU, I, J, N, NRHS
      CHARACTER(LEN=1) :: EQUED
!  .. "Local Arrays" ..
      REAL(WP), ALLOCATABLE :: AB(:,:), B(:,:), X(:,:), R(:), C(:), &
                               FERR(:), BERR(:)
!  .. "Executable Statements" ..
      WRITE (*,*) 'LA_GBSVX Example Program Results.'
      N = 6; KL = 2; KU = 1; NRHS = 2
!     ALLOCATE ( AB(2*KL+KU+1,N), B(N,NRHS), X(N,NRHS), R(N), C(N), &
      ALLOCATE ( AB(  KL+KU+1,N), B(N,NRHS), X(N,NRHS), R(N), C(N), &
                 FERR(NRHS),BERR(NRHS) )

      OPEN(UNIT=21,FILE='gbsv.ma',STATUS='UNKNOWN')
!     DO I=KL+1,2*KL+KU+1
      DO I=   1,  KL+KU+1
      DO J=1,N
         READ(21,'(F2.0)') AB(I,J)
      ENDDO
      ENDDO
      CLOSE(21)

      AB(:,1:1)=1e-6*AB(:,1:1)
      DO J = 1, N
        DO K = 1, 1+KU-MIN(J-1,KU)-1; AB(K,J) = HUGE(1.0_WP); END DO
        DO K = 1+ KU+MIN(N-J,KL)+KL+1, 1+KU+KL; AB(K,J) = HUGE(1.0_WP); END DO
      END DO
      WRITE(*,*) 'The array AB:'
      DO I=1,N; WRITE(*,*) AB(I,:); ENDDO

      B = 0.0_WP
      DO I = 1, NRHS;
        DO J = 1, N
          DO K = 1+KU-MIN(J-1,KU), 1+ KU+MIN(N-J,KL)
            B(J,I) = AB(K,J)+B(J,I);
          ENDDO
        ENDDO;
        B(:,I) = B(:,I)*I;
      ENDDO

      WRITE(*,*) 'The RHS matrix B:'
      DO I=1,N; WRITE(*,"(1(I3,1X),I3,1X)") INT(B(I,:)); ENDDO

      WRITE(*,*) "CALL LA_GBSVX( AB, B, X, 2, TRANS='T', EQUED=EQUED,", &
                 " R=R, C=C, FERR=FERR, BERR=BERR )"

      CALL LA_GBSVX( AB, B, X, 2, FACT='E',TRANS='T',EQUED=EQUED, &
                     R=R, C=C, FERR=FERR, BERR=BERR )

      WRITE(*,*)'X on exit: '
      DO I=1,N; WRITE(*,"(1(E12.6,1X),E12.6,1X)") X(I,:); ENDDO
      WRITE(*,*)'EQUED = ', EQUED
      WRITE(*,*)'R = ', R
      WRITE(*,*)'C = ', C
      WRITE(*,*)'FERR = ', FERR
      WRITE(*,*)'BERR = ', BERR

      END PROGRAM LA_SGBSVX_EXAMPLE