!/*****************************************************************************/
! *
! *  Elmer, A Finite Element Software for Multiphysical Problems
! *
! *  Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland
! * 
! *  This library is free software; you can redistribute it and/or
! *  modify it under the terms of the GNU Lesser General Public
! *  License as published by the Free Software Foundation; either
! *  version 2.1 of the License, or (at your option) any later version.
! *
! *  This library is distributed in the hope that it will be useful,
! *  but WITHOUT ANY WARRANTY; without even the implied warranty of
! *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
! *  Lesser General Public License for more details.
! * 
! *  You should have received a copy of the GNU Lesser General Public
! *  License along with this library (in file ../LGPL-2.1); if not, write 
! *  to the Free Software Foundation, Inc., 51 Franklin Street, 
! *  Fifth Floor, Boston, MA  02110-1301  USA
! *
! *****************************************************************************/
!/******************************************************************************
! *
! *  Subroutine for saving scalar data to files
! *
! ******************************************************************************
! *
! *  Authors: Peter Råback
! *  Email:   Peter.Raback@csc.fi
! *  Web:     http://www.csc.fi/elmer
! *  Address: CSC - IT Center for Science Ltd.
! *           Keilaranta 14
! *           02101 Espoo, Finland 
! *
! *  Original Date: 20 Nov 2001
! *
! *****************************************************************************/

!> \ingroup Solvers
!> \{

SUBROUTINE SaveScalars_init( Model,Solver,dt,TransientSimulation )
!------------------------------------------------------------------------------
  USE DefUtils

  IMPLICIT NONE
!------------------------------------------------------------------------------
  TYPE(Solver_t), TARGET :: Solver
  TYPE(Model_t) :: Model
  REAL(KIND=dp) :: dt
  LOGICAL :: TransientSimulation
!------------------------------------------------------------------------------
! Local variables
!------------------------------------------------------------------------------
  INTEGER :: LineInd, MarkerUnit, i
  LOGICAL :: GotIt, MarkFailed, AvoidFailed
  CHARACTER(LEN=MAX_NAME_LEN) :: Name
  
  
  ! If we want to show a pseudonorm add a variable for which the norm
  ! is associated with.
  IF ( ListCheckPresent( Solver % Values,'Show Norm Index') .OR. &
      ListCheckPresent( Solver % Values,'Show Norm Name') .OR. &
      ListCheckPresent( Solver % Values,'Reference Values') ) THEN
    Name = ListGetString( Solver % Values, 'Equation',GotIt)
    IF(.NOT. GotIt) Name = "SaveScalars"
    IF( .NOT. ListCheckPresent( Solver % Values,'Variable') ) THEN
      CALL ListAddString( Solver % Values,'Variable',&
          '-nooutput -global '//TRIM(Name)//'_var')
    END IF
  END IF

  IF( ParEnv % MyPe == 0 ) THEN
    MarkFailed = ListGetLogical( Solver % Values,'Mark Failed Strategy',GotIt)
    AvoidFailed = ListGetLogical( Solver % Values,'Avoid Failed Strategy',GotIt)
    IF(.NOT. GotIt) AvoidFailed = MarkFailed
    
    IF( MarkFailed .OR. AvoidFailed ) THEN
      LineInd = ListGetInteger( Solver % Values,'Line Marker',GotIt)
      IF(.NOT. GotIt) THEN
        CALL Fatal('SaveScalars_init','Failed strategy marked requires > Line Marker <')
      END IF
      Name = 'FINISHED_MARKER_'//I2S(LineInd)
    END IF

    IF( AvoidFailed ) THEN
      INQUIRE(FILE=TRIM(Name),EXIST=GotIt)
      IF( GotIt ) THEN     
        OPEN(NEWUNIT=MarkerUnit, FILE=Name)
        READ(MarkerUnit,*) i
        IF( i == 0 ) THEN
          CALL Fatal('SaveScalars_init','Strategy already failed before!')
        END IF
        CLOSE(MarkerUnit)
      END IF
    END IF
    
    ! Save a negative status during the execution such that if the
    ! program terminates the negative status will prevail
    IF( MarkFailed ) THEN
      CALL Info('SaveScalars_init','Saving False marker at start')
      i = 0
      OPEN(NEWUNIT=MarkerUnit,FILE=Name,STATUS='Unknown')
      WRITE(MarkerUnit,'(I0)') i
      CLOSE(MarkerUnit)
    END IF
  END IF


  
END SUBROUTINE SaveScalars_init



!------------------------------------------------------------------------------
!> This subroutine saves scalar values in ascii format to an external file.
!> This is a dynamically loaded solver with a standard interface.
!------------------------------------------------------------------------------
SUBROUTINE SaveScalars( Model,Solver,dt,TransientSimulation )
!------------------------------------------------------------------------------
  USE DefUtils
  USE Interpolation
  USE ElementUtils
  USE MeshUtils, ONLY : ClosestElementInMesh, ClosestNodeInMesh, NodeToDGIndex
  USE SaveUtils
  
  IMPLICIT NONE
!------------------------------------------------------------------------------
  TYPE(Solver_t), TARGET :: Solver
  TYPE(Model_t) :: Model
  REAL(KIND=dp) :: dt
  LOGICAL :: TransientSimulation
!------------------------------------------------------------------------------
! Local variables
!------------------------------------------------------------------------------
  TYPE(Solver_t), POINTER :: ParSolver
  TYPE(ValueList_t), POINTER :: Params
  TYPE(ValueListEntry_t), POINTER :: Lst
  TYPE(Variable_t), POINTER :: Var, OldVar, Var2, Var3
  TYPE(Mesh_t), POINTER :: Mesh
  TYPE(Nodes_t) :: ElementNodes
  TYPE(Element_t),POINTER :: Element

  LOGICAL :: MovingMesh, GotCoeff, &
      GotIt, GotOper, GotParOper, GotVar, GotOldVar, ExactCoordinates, VariablesExist, &
      ComplexEigenVectors, ComplexEigenValues, IsParallel, ParallelWrite, SaveCVS, &
      FileAppend, SaveEigenValue, SaveEigenFreq, IsInteger, ParallelReduce, WriteCore, &
      Hit, SaveToFile, EchoValues, BodyOper, BodyForceOper, &
      MaterialOper, MaskOper, GotMaskName, GotOldOper, ElementalVar, ComponentVar, &
      Numbering, NodalOper, GotNodalOper, SaveFluxRange, PosOper, NegOper, SaveJson, &
      StartNewFile, SimulationVisited

  LOGICAL, ALLOCATABLE :: ValuesInteger(:)

  LOGICAL, ALLOCATABLE :: ActiveBC(:)

  REAL (KIND=DP) :: Minimum, Maximum, AbsMinimum, AbsMaximum, &
      Mean, Variance, MinDist, x, y, z, Vol, Intmean, intvar, &
      KineticEnergy, PotentialEnergy, &
      Coords(3), LocalCoords(3), TempCoordinates(3), Val, Val2, &
      Change = 0._dp, Norm = 0.0_dp, PrevNorm, ParallelHits, ParallelCands

  REAL (KIND=DP), ALLOCATABLE :: Values(:), &
      CoordinateBasis(:), ElementValues(:), BoundaryFluxes(:),BoundaryAreas(:)

  REAL (KIND=DP), POINTER :: PointCoordinates(:,:), LineCoordinates(:,:), WrkPntr(:)
  INTEGER, ALLOCATABLE :: BoundaryHits(:)
  INTEGER, POINTER :: PointIndex(:), NodeIndexes(:), SaveIndex(:)
  INTEGER, ALLOCATABLE, TARGET :: ClosestIndex(:)
  CHARACTER(LEN=MAX_NAME_LEN), ALLOCATABLE :: ValueNames(:)
  CHARACTER(LEN=MAX_PATH_LEN) :: ScalarsFile, ScalarNamesFile, ScalarParFile
  CHARACTER(LEN=MAX_NAME_LEN) ::  DateStr, VariableName, OldVariableName, ResultPrefix, &
      Suffix, Oper, Oper0, OldOper0, ParOper, Name, CoefficientName, &
      MinOper, MaxOper, MaskName, OldMaskName, SaveName
  CHARACTER(:), ALLOCATABLE :: OutputDirectory
  INTEGER :: i,j,k,l,lpar,q,n,ierr,No,NoPoints,NoCoordinates,NoLines,NumberOfVars,&
      NoDims, NoDofs, NoOper, NoElements, NoVar, NoValues, PrevNoValues, DIM, &
      MaxVars, NoEigenValues, Ind, EigenDofs, LineInd, NormInd, CostInd, istat, nlen, &
      jsonpos, PassiveCoordinate
  INTEGER :: IntVal, FirstInd, LastInd, ScalarsUnit, MarkerUnit, NamesUnit, RunInd, PrevRunInd=-1
  LOGICAL, ALLOCATABLE :: NodeMask(:)
  INTEGER, ALLOCATABLE :: DGIndex(:)
  REAL (KIND=DP) :: CT, RT  
  LOGICAL :: SlicesReduce, TimesReduce, DoIt
  INTEGER :: PrevComm, CommRank, CommSize, nSlices, nTimes
  REAL(KIND=dp) :: Vals(100)
  LOGICAL :: GotEigen, GotEdge
  INTEGER :: NoVals
  CHARACTER(*), PARAMETER :: Caller = 'SaveScalars'
  
  SAVE :: jsonpos, PrevRunInd
  
!------------------------------------------------------------------------------

  CALL Info(Caller, '-----------------------------------------', Level=4 )
  CALL Info(Caller,'Saving scalar values of various kinds',Level=4)

  
  Mesh => GetMesh()
  DIM = CoordinateSystemDimension()
  Params => GetSolverParams()	

  MovingMesh = ListGetLogical(Params,'Moving Mesh',GotIt )

  FileAppend = ListGetLogical( Params,'File Append',GotIt)  

  SaveFluxRange = ListGetLogical( Params,'Save Flux Range',GotIt)
  IF(.NOT. GotIt) SaveFluxRange = .TRUE.

  RunInd = 0
  Var => VariableGet( Model % Variables,'run')
  IF( ASSOCIATED( Var ) ) RunInd = NINT( Var % Values(1) )

  SimulationVisited = .FALSE.
  
  ScalarsFile = ListGetString(Params,'Filename',SaveToFile )
  IF( SaveToFile ) THEN    
    ! Optionally number files by the number of partitions
    ! This makes the benchmarking more convenient since each case 
    ! may use the same command file
    IF(ListGetLogical(Params,'Partition Numbering',GotIt)) THEN
      i = INDEX( ScalarsFile,'.',.TRUE. )
      j = LEN_TRIM(ScalarsFile)
      IF(i > 0) THEN
        Suffix = ScalarsFile(i:j)
        WRITE( ScalarsFile,'(A,I0,A)') &
            ScalarsFile(1:i-1)//'_np',ParEnv % PEs,Suffix(1:j-i+1)
      ELSE
        WRITE( ScalarsFile,'(A,I0)') &
            ScalarsFile(1:j)//'_np',ParEnv % PEs 
      END IF
    END IF

    ! Optionally number files by the run control loop.
    IF(ListGetLogical(Params,'Run Control Numbering',GotIt)) THEN
      i = INDEX( ScalarsFile,'.',.TRUE. )
      j = LEN_TRIM(ScalarsFile)
      IF(i > 0) THEN
        Suffix = ScalarsFile(i:j)
        WRITE( ScalarsFile,'(A,I0,A)') &
            ScalarsFile(1:i-1)//'_run',RunInd,Suffix(1:j-i+1)
      ELSE
        WRITE( ScalarsFile,'(A,I0)') &
            ScalarsFile(1:j)//'_run',RunInd
      END IF
      IF( RunInd /= PrevRunInd ) Solver % TimesVisited = 0      
    END IF
    
    CALL SolverOutputDirectory( Solver, ScalarsFile, OutputDirectory )
    ! Make parallel reduction to ensure that the output directory has been created
    IF( Solver % TimesVisited == 0 ) THEN
      i = 1; i = ParallelReduction(i)
    END IF
    
    ScalarsFile = TRIM(OutputDirectory)// '/' //TRIM(ScalarsFile)
    
    Numbering = ListGetLogical(Params,'Filename Numbering',GotIt)

    IF( Numbering  ) THEN
      IF( Solver % TimesVisited > 0  ) THEN
        ScalarsFile = NextFreeFilename( ScalarsFile, LastExisting = .TRUE. ) 
      ELSE
        ScalarsFile = NextFreeFilename( ScalarsFile ) 
      END IF
    END IF
      
    ScalarNamesFile = TRIM(ScalarsFile) // TRIM(".names")
    SaveCVS = ListGetLogical(Params,'Live Graph',GotIt)
    IF(.NOT. GotIt) SaveCVS = ListGetLogical(Params,'CVS Format',GotIt)
    SaveJson = ListGetLogical(Params,'Json Format',GotIt)
  END IF

  EchoValues = ListGetLogical( Params,'Echo Values',GotIt)
  IF(.NOT. GotIt) EchoValues = .NOT. SaveToFile

  ResultPrefix = ListGetString(Params,'Scalars Prefix',GotIt )
  IF(.NOT. gotIt) ResultPrefix = 'res:'


  IsParallel = .FALSE.
  WriteCore = .TRUE.
  ParallelWrite = .FALSE.
  ParallelReduce = .FALSE.
  SlicesReduce = .FALSE.
  TimesReduce = .FALSE.
  PrevComm = 0
  CommRank = 0
  CommSize = ParEnv % PEs
  
  IF( CommSize > 1 ) THEN
    IsParallel = .TRUE.
    ParallelReduce = GetLogical( Params,'Parallel Reduce',GotIt)
    
    SlicesReduce = GetLogical( Params,'Slices Reduce',GotIt)
    TimesReduce = GetLogical( Params,'Times Reduce',GotIt)

    nSlices = 0
    nTimes = 0
    IF( SlicesReduce .OR. TimesReduce ) THEN
      ParallelReduce = .TRUE.
      nSlices = ListGetInteger( Model % Simulation,'Number Of Slices',GotIt)
      nTimes = ListGetInteger( Model % Simulation,'Number Of Times',GotIt)
      IF( nSlices > 1 .AND. nTimes > 1 ) THEN
        PrevComm = ParEnv % ActiveComm 
        IF( SlicesReduce ) ParEnv % ActiveComm = ParallelSlicesComm()
        IF( TimesReduce ) ParEnv % ActiveComm = ParallelTimesComm()
      END IF
      CommRank = ParallelPieceRank(ParEnv % ActiveComm)
      CommSize = ParallelPieceSize(ParEnv % ActiveComm)
    ELSE
      CommRank = ParEnv % MyPe 
    END IF

    IF( CommRank > 0 ) EchoValues = .FALSE.
    
    IF( ParallelReduce ) THEN
      WriteCore = ( CommRank == 0 )
      ParallelWrite = ( CommSize < ParEnv % PEs )
    ELSE
      ParallelWrite = .TRUE.      
    END IF
  END IF

  NoLines = 0
  LineCoordinates => ListGetConstRealArray(Params,'Polyline Coordinates',gotIt)
  IF(gotIt) THEN
    NoLines = SIZE(LineCoordinates,1) / 2
    NoDims = SIZE(LineCoordinates,2)
  END IF

  NoPoints = 0
  PointIndex => ListGetIntegerArray( Params,'Save Points',GotIt)
  IF ( gotIt ) THEN    
    NoPoints = SIZE(PointIndex)
  END IF

  NoCoordinates = 0
  NoElements = 0
  PointCoordinates => ListGetConstRealArray(Params,'Save Coordinates',gotIt)
  IF(gotIt) THEN
    NoDims = SIZE(PointCoordinates,2)
    ExactCoordinates = ListGetLogical(Params,'Exact Coordinates',GotIt )      

    IF( ParallelReduce .AND. .NOT. ExactCoordinates) THEN
      CALL Warn(Caller,'Only Exact Save Coordinates works in parallel, enforcing...')
      ExactCoordinates = .TRUE.
    END IF
    
    IF(ExactCoordinates) THEN            
      ! Look for the value at the given coordinate point really.
      NoElements = SIZE(PointCoordinates,1)
      GotIt = .FALSE.
      
      IF( .NOT. MovingMesh ) THEN
        SaveIndex => ListGetIntegerArray( Params,'Save Coordinate Elements',GotIt )
      END IF
      IF(.NOT. GotIt ) THEN
        CALL Info(Caller,'Searching for elements containing save coordinates',Level=8)

        ALLOCATE(ClosestIndex(NoElements), STAT=istat)
        IF( istat /= 0 ) CALL Fatal(Caller,'Memory allocation error for CoordinateElemNo')         
        DO j=1,NoElements
          Coords(1:NoDims) = PointCoordinates(j,1:NoDims)
          IF(NoDims < 3 ) Coords(NoDims+1:3) = 0.0_dp
          ClosestIndex(j) = ClosestElementInMesh( Mesh, Coords )
        END DO

        SaveIndex => ClosestIndex
        IF( .NOT. MovingMesh ) THEN
          CALL ListAddIntegerArray( Params,'Save Coordinate Elements',&
              NoElements,ClosestIndex )
        END IF
      END IF
    ELSE      
      ! Find the indexes of minimum distances
      NoCoordinates = SIZE(PointCoordinates,1)
      GotIt = .FALSE.
      IF( .NOT. MovingMesh ) THEN
        SaveIndex => ListGetIntegerArray( Params,'Save Coordinate Indexes',GotIt )
      END IF
      IF( .NOT. GotIt ) THEN
        CALL Info(Caller,'Searching for closest nodes to coordinates',Level=8)
        ALLOCATE(ClosestIndex(NoCoordinates), STAT=istat)
        IF( istat /= 0 ) CALL Fatal(Caller,'Memory allocation error for SaveIndex') 
        DO j=1,NoCoordinates 
          Coords(1:NoDims) = PointCoordinates(j,1:NoDims)
          IF(NoDims < 3 ) Coords(NoDims+1:3) = 0.0_dp
          ClosestIndex(j) = ClosestNodeInMesh( Mesh, Coords, DoParallel = .TRUE. )
        END DO
        
        SaveIndex => ClosestIndex
        IF( .NOT. MovingMesh ) THEN
          CALL ListAddIntegerArray( Params,'Save Coordinate Indexes',&
              NoCoordinates,ClosestIndex )
        END IF
      END IF
    END IF
  END IF
  
  n = NoPoints + NoCoordinates
  IF( n > 0 ) THEN
    IF( ASSOCIATED( Mesh % Elements(1) % DGIndexes ) ) THEN
      ALLOCATE( DGIndex(n) )       
      DGIndex = 0
      DO i=1,n
        IF( i<= NoPoints ) THEN
          DGIndex(i) = NodeToDGIndex(Mesh,PointIndex(i))
        ELSE
          j = SaveIndex(i-NoPoints)
          IF(j>0) DGIndex(i) = NodeToDGIndex(Mesh,j)
        END IF
      END DO
    END IF
  END IF
    
  
!------------------------------------------------------------------------------

  n = Mesh % MaxElementNodes 
  ALLOCATE( ElementNodes % x(n), ElementNodes % y(n), ElementNodes % z(n), &
      ElementValues( n ), CoordinateBasis(n), STAT=istat)
  IF( istat /= 0 ) CALL Fatal(Caller,'Memory allocation error 1') 	

  n = MAX( Model % NumberOfBodies, MAX(Model % NumberOfBCs, NoLines))
  ALLOCATE( BoundaryFluxes(n), BoundaryAreas(n), BoundaryHits(n), STAT=istat )
  IF( istat /= 0 ) CALL Fatal(Caller,'Memory allocation error 2') 	
  
  ALLOCATE( ActiveBC( Model % NumberOfBCs ), STAT=istat )
  IF( istat /= 0 ) CALL Fatal(Caller,'Memory allocation error 3') 	

  ComplexEigenVectors = ListGetLogical(Params,'Complex Eigen Vectors',GotIt)
  
      
  !------------------------------------------------------------------------------
  ! Go through the variables and compute the desired statistical data
  !------------------------------------------------------------------------------
  NoValues = 0
  GotVar  = .TRUE.
  GotOper = .FALSE.
  GotOldOper = .FALSE.
  NULLIFY(OldVar)
  NoVar = 0
  MinOper = 'min'
  MaxOper = 'max'
  GotNodalOper = .FALSE.
  OldMaskName = 'save scalars'
  PosOper = .FALSE.
  NegOper = .FALSE.
  
  
  DO WHILE(GotVar .OR. GotOper)
    
    GotOper = .FALSE.
    NULLIFY(Var)
    
    NoVar = NoVar + 1
    WRITE (Name,'(A,I0)') 'Variable ',NoVar
  
    VariableName = ListGetString( Params, TRIM(Name), GotVar )


    IF(TRIM(VariableName) == 'cpu time' .OR. TRIM(VariableName) == 'cpu memory') THEN
      CALL Warn(Caller,'This variable should now be invoked as an operator: '//TRIM(VariableName))
      CYCLE
    END IF
    
    GotOldVar = .FALSE.

    IF(GotVar) THEN
      Var => VariableGet( Model % Variables, TRIM(VariableName) )
      IF ( .NOT. ASSOCIATED( Var ) )  THEN
        Var => VariableGet( Model % Variables, TRIM(VariableName)//' 1' )
        IF( ASSOCIATED( Var ) ) THEN
          CALL Info(Caller,'Treating a component variable: '//TRIM(VariableName),Level=8)
          ComponentVar = .TRUE.
          Var2 => VariableGet( Model % Variables, TRIM(VariableName)//' 2' )
          Var3 => VariableGet( Model % Variables, TRIM(VariableName)//' 3' )          
        ELSE
          val = ListGetConstReal( Model % Simulation, TRIM(VariableName), GotIt )
          IF( GotIt ) THEN
            CALL Info(Caller,'Adding keyword from simulation section as variable: '&
                //TRIM(VariableName),Level=10)
            CALL AddToSaveList(TRIM(VariableName), val, .FALSE.)
            SimulationVisited = .TRUE.
            CYCLE
          ELSE
            CALL Warn(Caller,'Requested variable does not exist: '//TRIM(VariableName))
            CYCLE
          END IF
        END IF
      ELSE
        ComponentVar = .FALSE.
      END IF    
      OldVar => Var
      OldVariableName = VariableName
      
      ! A 0D variable cannot really be much operated, hence save it as is
      !-------------------------------------------------------------------
      IF(SIZE(Var % Values) == Var % Dofs) THEN
        IsInteger = .FALSE.
        IF( VariableName == 'timestep' ) IsInteger = .TRUE.
        IF( VariableName == 'nonlin iter' ) IsInteger = .TRUE.
        IF( VariableName == 'coupled iter' ) IsInteger = .TRUE.
        IF( VariableName == 'run' ) IsInteger = .TRUE.
        
        IF( Var % Dofs == 1 ) THEN
          CALL AddToSaveList('value: '//TRIM(VariableName)//' scalar variable', &
                              Var % Values(1), IsInteger )
        ELSE
          DO j=1,Var % DOfs
            CALL AddToSaveList('value: '//ComponentName(VariableName,j)//' scalar variable', Var % Values(j))          
          END DO
        END IF
        CYCLE
      END IF

      !EdgeBasis = ( Var % TYPE == variable_on_edges )
      !DGVar = ( Var % TYPE == variable_on_nodes_on_elements ) 
      !IpVar = ( Var % TYPE == variable_on_gauss_points )
      !ElemVar = ( Var % TYPE == Variable_on_elements )       
      !PiolaVersion = .FALSE.
      !pElem = .FALSE.
      
      WRITE (Name,'(A,I0)') 'Nodal Variable ',NoVar
      NodalOper = ListGetLogical(Params,TRIM(Name),GotNodalOper)   
    ELSE
      IF(ASSOCIATED(OldVar)) THEN
        Var => OldVar
        VariableName = OldVariableName
	GotOldVar = .TRUE.
      END IF
    END IF

    WRITE (Name,'(A,I0)') 'Nodal Variable ',NoVar
    IF( ListCheckPresent( Params,TRIM(Name) ) ) THEN
      NodalOper = ListGetLogical(Params,TRIM(Name),GotNodalOper)   
    END IF
      
    NoOper = NoVar     
    MaskOper = .FALSE.
    WRITE (Name,'(A,I0)') 'Operator ',NoOper
    Oper0 = ListGetString(Params,TRIM(Name),GotOper)
    IF(.NOT. GotOper .AND. GotOldOper ) Oper0 = OldOper0


    WRITE (Name,'(A,I0)') 'Mask Name ',NoOper
    MaskName = ListGetString(Params,TRIM(Name),GotMaskName)
    IF( GotMaskName ) THEN
      OldMaskName = MaskName
    ELSE
      MaskName = OldMaskName
    END IF

    
    IF(.NOT. (GotOper .OR. GotVar .OR. GotMaskName ) ) CYCLE


    IF( ASSOCIATED( Var ) ) THEN
      CALL Info(Caller,'Treating variable: '//TRIM(VariableName),Level=12)
      ElementalVar = ( Var % TYPE == Variable_on_nodes_on_elements ) 
    END IF

    IF( GotOper ) THEN
      CALL Info(Caller,'Treating operator: '//TRIM(Oper0),Level=12)
      OldOper0 = Oper0
      GotOldOper = .TRUE.
    ELSE IF( GotOldOper ) THEN
      Oper0 = OldOper0
      GotOper = GotOldOper
    ELSE
      CALL Info(Caller,'No operator given for variable: '//TRIM(VariableName))
      CYCLE
    END IF

    Oper = ""
    BodyOper = .FALSE.
    BodyForceOper = .FALSE.      
    MaterialOper = .FALSE.
    nlen = LEN_TRIM(Oper0) 
    IF( Oper0(1:11) == 'body force ') THEN
      BodyForceOper = .TRUE.
      j = 11
    ELSE IF( Oper0(1:5) == 'body ') THEN
      BodyOper = .TRUE.
      j = 5
    ELSE IF( Oper0(1:9) == 'material ') THEN
      MaterialOper = .TRUE.
      j = 9
    ELSE
      j = 0
    END IF
    Oper(1:nlen-j) = Oper0(1+j:nlen)
    
    MaskOper = ( BodyForceOper .OR. BodyOper .OR. MaterialOper )
    IF( MaskOper ) THEN
      CALL Info(Caller,'Operator to be masked: '//TRIM(Oper),Level=12)
    END IF

    nlen = LEN_TRIM(Oper) 
    PosOper = .FALSE.
    NegOper = .FALSE.
    j = 0
    IF( Oper(1:9) == 'positive ') THEN
      PosOper = .TRUE.
      j = 9
    ELSE IF( Oper(1:9) == 'negative ') THEN
      NegOper = .TRUE.
      j = 9 
    END IF
    IF(j>0) Oper(1:nlen-j) = Oper(1+j:nlen)
    
    ! We may want to do integrals over projected surfaces
    PassiveCoordinate = ListGetInteger( Params,'Passive Coordinate',GotIt )
    IF(.NOT. GotIt ) THEN
      PassiveCoordinate = ListGetInteger( Params,'Passive Coordinate '//I2S(NoOper), GotIt )
    END IF

    WRITE (Name,'(A,I0)') 'Coefficient ',NoOper
    CoefficientName = ListGetString(Params,TRIM(Name),GotCoeff )

    WRITE (Name,'(A,I0)') 'Parallel Operator ',NoOper
    ParOper = ListGetString(Params,TRIM(Name),GotParOper)
    IF(.NOT. GotParOper) ParOper = OperToParOperMap(Oper)

    IF( MaskOper ) THEN
      GotIt = .FALSE.
      IF( BodyOper ) THEN
        GotIt = ListGetLogicalAnyBody( Model, MaskName )
      ELSE IF( BodyForceOper ) THEN
        GotIt = ListGetLogicalAnyBodyForce( Model, MaskName )
      ELSE IF( MaterialOper ) THEN
        GotIt = ListGetLogicalAnyMaterial( Model, MaskName )
      END IF
      IF(.NOT. GotIt ) THEN
        CALL Warn(Caller,'Masked operators require mask: '//TRIM(MaskName))
      END IF
    END IF
    
    ActiveBC = .FALSE.
    DO j=1,Model % NumberOfBCs
      ActiveBC(j) =  &
          ListGetLogical(Model % BCs(j) % Values,'Flux Integrate',gotIt) .OR. &
          ListGetLogical(Model % BCs(j) % Values, MaskName, gotIt)
    END DO
    
    IF ( GotOper ) THEN
      
      SELECT CASE( Oper ) 
        
      CASE ('partitions')
      CASE ('partition checksum')
      CASE ('partition neighbours checksum')
      CASE ('cpu time')
      CASE ('wall time')
      CASE ('cpu memory') 
      CASE ('nodes')
      CASE ('elements')
      CASE ('bounding box')
      CASE ('area')
      CASE ('volume')
      CASE ('threads')

      CASE DEFAULT
        IF( .NOT. (GotVar .OR. GotOldVar ) ) THEN
          CALL Error(Caller,'Operator > '//TRIM(Oper)//' < requires variable')
          CYCLE
        END IF
      END SELECT


      ! Set default name for saving 
      IF(GotVar .OR. GotOldVar ) THEN
        SaveName = TRIM(Oper0)//': '//TRIM(VariableName)
        IF( GotMaskName ) THEN
          SaveName = TRIM(SaveName)//' mask '//TRIM(MaskName)
        END IF
        IF( GotNodalOper ) THEN
          IF( NodalOper ) THEN
            SaveName = TRIM(SaveName)//' nodal'
          ELSE
            SaveName = TRIM(SaveName)//' non-nodal'            
          END IF
        END IF
      END IF
        
      
      SELECT CASE(Oper)

      CASE ('partitions')
        Val = 1.0_dp * ParEnv % PEs 
        SaveName = 'value: number of partitions'
        CALL AddToSaveList(SaveName,Val,.TRUE.,ParOper)

      CASE ('threads')
        Val = 1.0_dp * ParEnv % NumberOfThreads 
        SaveName = 'value: number of threads'
        CALL AddToSaveList(SaveName,Val,.TRUE.,ParOper)

      CASE ('partition checksum')
        Val = 0.0_dp
        IF( IsParallel ) THEN
          Val = 1.0_dp * SUM( 1.0_dp * Mesh % ParallelInfo % GlobalDOFS ) 
          ! Give different partition different weight to create something like a checksum
          Val = ( ParEnv % MyPe + 1 ) * Val
        END IF
        SaveName = 'value: partition checksum'
        ! Don't use integer as type because it can exceed the bounds, long should be used
        CALL AddToSaveList(SaveName,Val,.FALSE.,ParOper)

      CASE ('partition neighbours checksum')
        Val = 0.0_dp
        IF( IsParallel ) THEN
          DO j=1,Mesh % NumberOfNodes 
            Val = Val + 1.0_dp * SUM( Mesh % ParallelInfo % NeighbourList(j) % Neighbours )
          END DO
          Val = ( ParEnv % MyPe + 1 ) * Val
        END IF
        SaveName = 'value: partition neighbours checksum'
        ! Don't use integer as type because it can exceed the bounds, long should be used
        CALL AddToSaveList(SaveName,Val,.FALSE.,ParOper)

      CASE ('cpu time')
        Val = CPUTime()
        SaveName = 'value: cpu time (s)'
        CALL AddToSaveList(SaveName,Val,.FALSE.,ParOper)
      
      CASE ('wall time')
        Val = RealTime()
        SaveName = 'value: real time (s)'
        CALL AddToSaveList(SaveName,Val,.FALSE.,ParOper)
      
      CASE ('cpu memory') 
        Val = CPUMemory()
        SaveName = 'value: maximum memory usage (kb)'
        CALL AddToSaveList(SaveName,Val,.FALSE.,ParOper)

      CASE ('nodes')
        Val = 1.0_dp * Solver % Mesh % NumberOfNodes
        SaveName = TRIM(Oper)//': '//TRIM(VariableName)
        CALL AddToSaveList(SaveName,Val,.TRUE.,ParOper)

      CASE ('elements')
        Val = 1.0_dp * Solver % Mesh % NumberOfBulkElements
        SaveName = TRIM(Oper)//': '//TRIM(VariableName)
        CALL AddToSaveList(SaveName,Val,.TRUE.,ParOper)

      CASE ('bounding box')
        Val = MINVAL( Solver % Mesh % Nodes % x ) 
        CALL AddToSaveList(TRIM(Oper)//' min x',Val,ParallelOperator=MinOper)
        Val = MAXVAL( Solver % Mesh % Nodes % x ) 
        CALL AddToSaveList(TRIM(Oper)//' max x',Val,ParallelOperator=MaxOper)
        Val = MINVAL( Solver % Mesh % Nodes % y ) 
        CALL AddToSaveList(TRIM(Oper)//' min y',Val,ParallelOperator=MinOper)
        Val = MAXVAL( Solver % Mesh % Nodes % y ) 
        CALL AddToSaveList(TRIM(Oper)//' max y',Val,ParallelOperator=MaxOper)
        IF( Solver % Mesh % MeshDim > 2 ) THEN
          Val = MINVAL( Solver % Mesh % Nodes % z ) 
          CALL AddToSaveList(TRIM(Oper)//' min z',Val,ParallelOperator=MinOper)
          Val = MAXVAL( Solver % Mesh % Nodes % z )
          CALL AddToSaveList(TRIM(Oper)//' max z',Val,ParallelOperator=MaxOper)
        END IF

      ! Operators that require a variable (except for 'area' and 'volume')
      !-----------------------------------------------------------------
      CASE ('norm')
        Val = Var % Norm
        CALL AddToSaveList(SaveName,Val,.FALSE.,ParOper)
        
      CASE ('nonlin change')
        Val = Var % NonlinChange
        CALL AddToSaveList(SaveName,Val,.FALSE.,ParOper)
        
      CASE ('steady state change')
        Val = Var % SteadyChange
        CALL AddToSaveList(SaveName,Val,.FALSE.,ParOper)
        
      CASE ('nonlin iter')
        Val = Var % NonlinIter
        CALL AddToSaveList(SaveName,Val,.TRUE.,ParOper)
        
      CASE ('nonlin converged')
        Val = 1.0_dp * Var % NonlinConverged
        CALL AddToSaveList(SaveName,Val,.TRUE.,ParOper)
        
      CASE ('steady converged')
        Val = 1.0_dp * Var % SteadyConverged
        CALL AddToSaveList(SaveName,Val,.TRUE.,ParOper)
        
      CASE ('dofs')
        Val = 1.0_dp * SIZE(Var % Values)
        CALL AddToSaveList(SaveName,Val,.TRUE.,ParOper)

      CASE ('nans')  ! number of NaN:s in vector
        GotIt = .FALSE.
        j = 0
        DO i=1,SIZE(Var % Values)
          IF( Var % Values(i) /= Var % Values(i) ) THEN
            j = j + 1
          END IF
        END DO
        Val = 1.0_dp * j
        CALL AddToSaveList(SaveName,Val,.TRUE.,ParOper)
       
      CASE ('sum','sum abs','mean abs','max','max abs','min','min abs',&
          'mean','variance','range','sum square','mean square','rms')
        IF( MaskOper ) CALL CreateNodeMask()
        IF( GotNodalOper ) THEN
          Val = VectorStatistics(Var,Oper,NodalOper)
        ELSE
          Val = VectorStatistics(Var,Oper)
        END IF
        CALL AddToSaveList(SaveName,Val,.FALSE.,ParOper)
        
      CASE ('deviation')
        IF( MaskOper ) CALL CreateNodeMask()
        Val = VectorMeanDeviation(Var,Oper)
        CALL AddToSaveList(SaveName, Val,.FALSE.,ParOper)
        
      CASE ('int','int mean','int square','int square mean','int rms','int abs','int abs mean',&
          'int variance','volume','potential energy','diffusive energy','convective energy')
        
        IF( MaskOper ) CALL CreateNodeMask()
        Val = BulkIntegrals(Var, Oper, GotCoeff, CoefficientName)
        IF(GotCoeff) THEN
          SaveName = TRIM(SaveName)//' with '//TRIM(CoefficientName)
        END IF
        CALL AddToSaveList(SaveName, Val,.FALSE.,ParOper)
        
      CASE('boundary sum','boundary dofs','boundary max','boundary max abs','boundary min',&
          'boundary min abs','boundary mean')
        
        IF( .NOT. ANY( ActiveBC ) ) THEN
          CALL Error(Caller,'No flag > '//TRIM(MaskName)// &
              ' < active for operator: '// TRIM(Oper))
          CYCLE
        END IF

        ! Statistical operators do not really like being split among bundaries.
        ! So if you want split the results use different mask names.
        BoundaryHits = 0
        BoundaryFluxes = 0.0_dp

        CALL BoundaryStatistics(Var, Oper, GotCoeff, &
            CoefficientName, BoundaryFluxes, BoundaryHits)
        
        IF( TRIM(Oper) == 'boundary mean' ) THEN
          IF( IsParallel .AND. ParallelReduce ) THEN
            CALL Warn(Caller,'Operator > boundary mean < not implemented in parallel!')
          ELSE IF( BoundaryHits(1) > 0 ) THEN
            BoundaryFluxes(1) = BoundaryFluxes(1) / BoundaryHits(1)
          END IF
        END IF
        WRITE (Name,'(A,A,A,A,I0)') TRIM(Oper),': ',TRIM(VariableName),' over bc '//TRIM(MaskName)
        CALL AddToSaveList( TRIM(Name), BoundaryFluxes(1),.FALSE.,ParOper)
        
      CASE ('boundary int','boundary int mean','area','diffusive flux','convective flux')
        
        IF( .NOT. ANY( ActiveBC ) ) THEN
          CALL Fatal(Caller,'No flag > '//TRIM(MaskName)// &
              '< active for operator: '// TRIM(Oper))
        ELSE
          BoundaryHits = 0
          BoundaryFluxes = 0.0_dp
          BoundaryAreas = 0.0_dp
          IF( SaveFluxRange ) THEN
            Minimum = HUGE(Minimum) 
            Maximum = -HUGE(Maximum)          
          END IF
            
          CALL BoundaryIntegrals(Var, Oper, GotCoeff, CoefficientName,&
              BoundaryFluxes,BoundaryAreas,BoundaryHits)
          
          DO j=1,Model % NumberOfBCs
            IF( ActiveBC(j) ) THEN
              IF( TRIM(Oper) == 'boundary int mean' ) THEN
                IF( IsParallel .AND. ParallelReduce ) THEN
                  CALL Warn(Caller,'Operator > boundary int mean < not implemented in parallel!')
                ELSE IF( BoundaryAreas(j) > 0.0 ) THEN
                  BoundaryFluxes(j) = BoundaryFluxes(j) / BoundaryAreas(j)
                END IF
              END IF
              WRITE (Name,'(A,A,A,A,I0)') TRIM(Oper),': ',TRIM(VariableName),' over bc ',j
              CALL AddToSaveList( TRIM(Name), BoundaryFluxes(j),.FALSE.,ParOper)
            END IF
          END DO

          IF( SaveFluxRange ) THEN
            IF(TRIM(Oper) == 'diffusive flux' .OR. TRIM(Oper) == 'convective flux') THEN
              ParOper = 'min'
              WRITE (Name,'(A,A,A,A)') 'min ',TRIM(Oper),': ',TRIM(VariableName)
              CALL AddToSaveList( TRIM(Name), Minimum,.FALSE.,ParOper)
              
              ParOper = 'max'
              WRITE (Name,'(A,A,A,A)') 'max ',TRIM(Oper),': ',TRIM(VariableName)
              CALL AddToSaveList( TRIM(Name), Maximum,.FALSE.,ParOper)
            END IF
          END IF
        END IF


        IF(NoLines > 0) THEN
          BoundaryHits = 0
          BoundaryFluxes = 0.0_dp
          BoundaryAreas = 0.0_dp
          
          CALL PolylineIntegrals(Var, Oper, GotCoeff, CoefficientName,&
              BoundaryFluxes,BoundaryAreas, BoundaryHits)
          
          DO j=1,NoLines
            IF( TRIM(Oper) == 'boundary int mean' ) THEN
              IF( BoundaryHits(j) > 0 ) THEN 
                BoundaryFluxes(j) = BoundaryFluxes(j) / BoundaryAreas(j)
              END IF
            END IF
            WRITE (Name,'(A,A,A,A,I0)') TRIM(Oper),': ',TRIM(VariableName),' over polyline ',j
            CALL AddToSaveList( TRIM(Name), BoundaryFluxes(j),.FALSE.,ParOper)
          END DO
        END IF
        
      CASE DEFAULT 
        
        WRITE (Message,'(A,A)') 'Unknown operator: ',TRIM(Oper)
        CALL WARN(Caller,Message)
        
      END SELECT

    END IF
  END DO

  IF( NoVar > 0 ) THEN
    WRITE (Message,'(A,I0,A)') 'Performed ',NoVar-1,' reduction operations'
    CALL Info(Caller,Message,Level=7)
  END IF


  !------------------------------------------------------------------------------
  !   Add eigenvalues on the list if told to
  !------------------------------------------------------------------------------

  SaveEigenValue = ListGetLogical( Params, 'Save Eigenvalues', GotIt )
  IF(.NOT. GotIt) & 
      SaveEigenValue = ListGetLogical( Params, 'Save Eigen values', GotIt )

  SaveEigenFreq = ListGetLogical( Params, 'Save Eigenfrequencies', GotIt )
  IF(.NOT. GotIt) &
      SaveEigenFreq = ListGetLogical( Params, 'Save Eigen Frequencies', GotIt ) 

  IF ( SaveEigenValue .OR. SaveEigenFreq ) THEN
    ComplexEigenValues = ListGetLogical(Params,'Complex Eigen Values',GotIt)
    IF(.NOT. GotIt) &
        ComplexEigenValues = ListGetLogical(Params,'Complex EigenValues',GotIt)
    
    l = 0
    DO i = 1, Model % NumberOfSolvers       

      ! When we have strong coupling then two solvers will have exactly same eigenvalues.
      ! Let's skip the one which is passive.
      IF( ListGetLogical( Model % Solvers(i) % Values,'Linear System Solver Disabled',GotIt ) ) CYCLE

      IF ( Model % Solvers(i) % NOFEigenValues > 0 ) THEN
        DO k = 1, Model % Solvers(i) % NOFEigenValues
          
          Val = REAL( Model % Solvers(i) % Variable % EigenValues(k) )
          l = l + 1

          IF ( SaveEigenValue ) THEN
            WRITE( Name, '("eigen: Re value ", I0)' ) k
            CALL AddToSaveList(TRIM(Name), Val)
            IF( ComplexEigenValues ) THEN
              Val2 = AIMAG( Model % Solvers(i) % Variable % EigenValues(k) )
              WRITE( Name, '("eigen: Im value ", I0)' ) k
              CALL AddToSaveList(TRIM(Name), Val2)
            END IF
          END IF 

          IF ( SaveEigenFreq ) THEN
            WRITE( Name, '("eigen: frequency ", I0, " [Hz]")' ) k
            IF ( Val >= 0.0 ) THEN
              Val = SQRT(Val) / (2*PI)
            ELSE
              ! If the eigenvalue is negative take take take a square root of its absolute value and
              ! return a negative frequency. 
              Val = -SQRT(-Val) / (2*PI)
            END IF
            CALL AddToSaveList(TRIM(Name), Val)
          END IF

        END DO
      END IF
    END DO
    WRITE (Message,'(A,I0,A)') 'Found ',l,' Eigenvalues'
    CALL Info(Caller,Message)
  END IF
    
  !------------------------------------------------------------------------------
  ! Get the info at given node points
  !------------------------------------------------------------------------------
  DO k=1,NoPoints+NoCoordinates

    IF(k <= NoPoints) THEN
      l = PointIndex(k)
    ELSE
      l = SaveIndex(k-NoPoints)
    END IF

    lpar = ParallelReduction(l, 2 ) 
    IF( lpar == 0 ) CYCLE
        
    Var => Model % Variables

    DO WHILE( ASSOCIATED( Var ) )
      
      IF ( .NOT. Var % Output .OR. SIZE(Var % Values) == Var % DOFs ) THEN
        
        CONTINUE

      ELSE IF (ASSOCIATED (Var % EigenVectors)) THEN

        NoEigenValues = SIZE(Var % EigenValues) 
        EigenDofs = SIZE( Var % EigenVectors(1,:) ) / SIZE( Var % Perm )

        IF(EigenDofs == Var % DOFs) THEN
          DO j=1,NoEigenValues
            DO i=1,Var % DOFs              
              Ind = 0
              IF( l>0 ) Ind = Var % Perm(l)              

              NoVals = 0
              IF( Ind > 0 ) THEN
                Vals(1) = REAL( Var % EigenVectors(j, Var%Dofs*(Ind-1)+i) )
                NoVals = 1
                IF(ComplexEigenVectors) THEN
                  Vals(2) = AIMAG( Var % EigenVectors(j, Var%Dofs*(Ind-1)+i) )
                  NoVals = 2
                END IF
              END IF
                
              NoVals = ParallelReduction(NoVals)
              IF( NoVals > 0 ) THEN              
                IF(Var % DOFs == 1) THEN
                  WRITE(Name,'("value: Re Eigen ",I0," ",A," at node ",I0)') j,TRIM(Var % Name),lpar
                ELSE
                  WRITE(Name,'("value: Re Eigen ",I0," ",A,I2," at node ",I0)') j,TRIM(Var % Name),i,lpar
                END IF
                CALL AddToSaveList( TRIM(Name), Vals(1))                               
                IF(ComplexEigenVectors) THEN
                  IF(Var % DOFs == 1) THEN
                    WRITE(Name,'("value: Im Eigen ",I0," ",A," at node ",I0)') j,TRIM(Var % Name),lpar
                  ELSE
                    WRITE(Name,'("value: Im Eigen ",I0," ",A,I2," at node ",I0)') j,TRIM(Var % Name),i,lpar
                  END IF
                  CALL AddToSaveList( TRIM(Name), Vals(2))
                END IF
              END IF
            END DO
          END DO
        END IF

      ELSE IF( Var % Dofs == 1) THEN          
        ! The variables exist always also as scalars, therefore omit vectors.         
        Ind = l
        NoVals = 0
        IF(l>0) THEN
          IF(ASSOCIATED(Var % Perm)) THEN
            IF( Var % TYPE == variable_on_nodes_on_elements ) THEN
              Ind = Var % Perm(DGIndex(k))
            ELSE
              IF(ASSOCIATED(Var % Perm)) Ind = Var % Perm(l)
            END IF
          END IF
          val = 0.0_dp
          IF(Ind > 0) THEN
            val = Var % Values(Ind)
            NoVals = 1
          END IF
        ELSE
          val = -HUGE(val)
        END IF

        NoVals = ParallelReduction(NoVals,2)        
        IF(NoVals > 0) THEN
          val = ParallelReduction(val,2)
          WRITE(Name,'("value: ",A," at node ",I0)') TRIM( Var % Name ), lpar
          CALL AddToSaveList( TRIM(Name), val)        
        END IF
      END IF

      Var => Var % Next      
    END DO
  END DO
  
  IF( NoPoints + NoCoordinates > 0 ) THEN
    WRITE (Message,'(A,I0,A)') 'Tabulated all field values at ',NoPoints+NoCoordinates,' points'
    CALL Info(Caller,Message)
  END IF
  
  !------------------------------------------------------------------------------
  ! Get the info at exact coordinates within elements
  !------------------------------------------------------------------------------

  !------------------------------------------------------------------------------
  ! For parallel cases set the default to -HUGE and take the max of the values
  ParOper = 'max'

  DO k=1,NoElements        
    l = SaveIndex(k)

    ! In parallel only one partition should have found the element.
    lpar = l
    IF( l > 0 ) THEN
      Element => Mesh % Elements(l)
      IF( IsParallel ) lpar = Element % GElementIndex
    END IF

    lpar = ParallelReduction(lpar, 2 ) 
    IF( lpar == 0 ) CYCLE

    ! Only try to find the data in the active element/partition. 
    IF( l > 0 ) THEN
      n = Element % TYPE % NumberOfNodes

      NodeIndexes => Element % NodeIndexes

      Coords(1:NoDims) = PointCoordinates(k,1:NoDims)
      IF(NoDims < 3 ) Coords(NoDims+1:3) = 0.0_dp

      CALL CopyElementNodesFromMesh( ElementNodes, Mesh, n, NodeIndexes)
      
      Hit = PointInElement( Element, ElementNodes, &
          Coords, LocalCoords, GlobalEps = 1.0_dp, LocalEps=0.1_dp )	          

      ElementValues(1:n) = 0.0d0          
      CoordinateBasis = 0.0_dp
      DO q=1,N
        ElementValues(q) = 1.0d0
        CoordinateBasis(q) = InterpolateInElement( Element, ElementValues, &
            LocalCoords(1), LocalCoords(2), LocalCoords(3) )
        ElementValues(q) = 0.0d0
      END DO
    END IF


    GotEigen = .FALSE.
    GotEdge  = .FALSE.

    Var => Model % Variables
    DO WHILE( ASSOCIATED( Var ) )

      IF ( .NOT. Var % Output .OR. SIZE(Var % Values) == Var % DOFs) THEN
        CONTINUE 
      ELSE IF( Var % Dofs == 1 ) THEN
        NoVals = 0

        ! The active partition knows the number of variables. 
        IF(l>0) THEN
          CALL EvaluateVariableAtGivenPoint(NoVals,Vals,Mesh,Var,Element=Element,&
              LocalCoord=LocalCoords, GotEigen=GotEigen, GotEdge=GotEdge)
        END IF
        NoVals = ParallelReduction(NoVals,2)

        DO i = 1, NoVals
          IF(l>0) THEN
            val = Vals(i)
          ELSE
            val = -HUGE(val)
          END IF
          val = ParallelReduction(val,2)

          IF( GotEigen ) THEN 
            Name = "value: Eigen "//TRIM(I2S(i))//" "//TRIM(Var % Name)
          ELSE IF( GotEdge ) THEN
            IF( i <= 3 ) THEN
              Name = "value: "//TRIM(Var % Name)//' {e} '//TRIM(I2S(i))
            ELSE IF(i == 4 ) THEN
              Name = "value: "//TRIM(Var % Name)//' nodal'
            END IF
          ELSE 
            IF( NoVals == 1 ) THEN
              Name = "value: "//TRIM(Var % Name)
            ELSE
              Name = "value: "//TRIM(Var % Name)
            END IF
          END IF
          Name = TRIM(Name)//' in element '//TRIM(I2S(lpar))
          CALL AddToSaveList(TRIM(Name), val,.FALSE.,ParOper)
        END DO
      END IF

      Var => Var % Next      
    END DO
  END DO

  IF( NoElements > 0 ) THEN
    WRITE (Message,'(A,I0,A)') 'Tabulated points within ',NoElements,' elements'
    CALL Info(Caller,Message)
  END IF

  !------------------------------------------------------------------------------
  ! Update time elapsed from start if requested
  !------------------------------------------------------------------------------
  IF( ListGetLogical( Model % Simulation,'Simulation Timing',GotIt ) ) THEN
    CALL Info(Caller,'Adding information on simulation timing',Level=10)
    CT = CPUTime() - ListGetConstReal( Model % Simulation,'cputime0',GotIt)
    RT = RealTime() - ListGetConstReal( Model % Simulation,'realtime0',GotIt)
    CALL AddToSaveList('simulation: cpu time (s)',CT )
    CALL AddToSaveList('simulation: real time (s)',RT )
  END IF

  !------------------------------------------------------------------------------
  ! Read the scalars defined in other modules
  !------------------------------------------------------------------------------
  Lst => ListHead(Model % Simulation)
  l = 0
  DO WHILE( ASSOCIATED( Lst ) )    
    IF ( Lst % Name(1:4) == TRIM(ResultPrefix) ) THEN
      IF ( ASSOCIATED(Lst % Fvalues) ) THEN
        CALL AddToSaveList(Lst % Name, Lst % Fvalues(1,1,1), &
            CheckForDuplicates = SimulationVisited )
        l = l + 1
      END IF
    END IF
    Lst => Lst % Next
  END DO
  IF(l > 0 ) THEN
    WRITE (Message,'(A,I0,A)') 'Found ',l,' result scalars in simulation section'
    CALL Info(Caller,Message,Level=7)
  END IF

  !------------------------------------------------------------------------------
  !   Add free form expressions with GetCReal
  !------------------------------------------------------------------------------
  NoVar = 0
  GotVar = .TRUE.
  DO WHILE( GotVar )
    NoVar = NoVar + 1    
    WRITE (Name,'(A,I0)') 'Expression ',NoVar
    Val = ListGetCReal( Params, TRIM(Name), GotVar )
    IF( GotVar ) CALL AddToSaveList(TRIM(Name),Val)
  END DO

  !------------------------------------------------------------------------------
  ! Add results in Components
  !------------------------------------------------------------------------------
  DoIt = ListGetLogical( Params,'Save Component Results',GotIt )
  IF(.NOT. GotIt) DoIt = (Model % NumberOfComponents > 0)
  IF(DoIt) THEN
    CALL Info(Caller,'Saving results from component - if any',Level=10)
    l = 0
    DO i = 1, Model % NumberOfComponents
      Lst => ListHead( Model % Components(i) % Values )
      DO WHILE( ASSOCIATED( Lst ) )    
        IF ( Lst % Name(1:4) == TRIM(ResultPrefix) ) THEN
          IF ( ASSOCIATED(Lst % Fvalues) ) THEN
            CALL AddToSaveList('component '//I2S(i)//': '//TRIM(Lst % Name), Lst % Fvalues(1,1,1))
            l = l + 1
           END IF
        END IF
        Lst => Lst % Next
      END DO
    END DO
    IF(l > 0 ) THEN
      WRITE (Message,'(A,I0,A)') 'Found ',l,' result scalars in components section'
      CALL Info(Caller,Message,Level=7)
    END IF
  END IF
  

  !------------------------------------------------------------------------------
  ! If there are no values 
  !------------------------------------------------------------------------------
  IF( NoValues == 0 ) THEN
    CALL Warn(Caller,'Found no values to save')
    RETURN
  ELSE
    CALL Info(Caller,'Found '//I2S(NoValues)//' values to save in total',Level=6)
  END IF


  
  !------------------------------------------------------------------------------
  ! Finally save all the scalars into a file 
  !------------------------------------------------------------------------------
  IF( SaveToFile ) THEN

    LineInd = ListGetInteger( Params,'Line Marker',GotIt)
    PrevNoValues = ListGetInteger( Params,'Save Scalars Dofs',GotIt) 

    IF(WriteCore .AND. ( NoValues /= PrevNoValues .OR. RunInd /= PrevRunInd ) ) THEN 
      CALL ListAddInteger( Params,'Save Scalars Dofs',NoValues )

      WRITE( Message, '(A)' ) 'Saving names of values to file: '//TRIM(ScalarNamesFile)
      CALL Info( Caller, Message, Level=4 )
      
      IF( Solver % TimesVisited > 0 .AND. NoValues /= PrevNoValues ) THEN
        WRITE ( Message,'(A,I0,A,I0)') 'Number of scalar values differ from previous time: ',&
            NoValues,' vs. ',PrevNoValues
        CALL Warn(Caller,Message)
      END IF
      
      IF(ParallelWrite) CALL Info(Caller,'Parallel data is written into separate files',Level=6)
      IF(ParallelReduce) CALL Info(Caller,'Parallel data is reduced into one file',Level=6)
      IF(FileAppend) CALL Info(Caller,'Data is appended to existing file',Level=6)
      
      OPEN(NEWUNIT=NamesUnit, FILE=ScalarNamesFile,IOSTAT=istat)
      IF(istat /= 0) THEN
        CALL Fatal(Caller,'Could not open fie for saving: '//TRIM(ScalarNamesFile))
      END IF
      
      Message = ListGetString(Model % Simulation,'Comment',GotIt)
      IF( GotIt ) THEN
        WRITE(NamesUnit,'(A)') TRIM(Message)
      END IF

      Message = ListGetString(Params,'Comment',GotIt)
      IF( GotIt ) THEN
        WRITE(NamesUnit,'(A)') TRIM(Message)
      END IF

      WRITE(NamesUnit,'(A,A)') 'Metadata for SaveScalars file: ',TRIM(ScalarsFile)

      DateStr = GetVersion()
      WRITE( NamesUnit,'(A)') 'Elmer version: '//TRIM(DateStr)     

      DateStr = GetRevision( GotIt )
      IF( GotIt ) THEN
        WRITE( NamesUnit,'(A)') 'Elmer revision: '//TRIM(DateStr)
      END IF        

      DateStr = GetCompilationDate( GotIt )
      IF( GotIt ) THEN
        WRITE( NamesUnit,'(A)') 'Elmer compilation date: '//TRIM(DateStr)
      END IF

      DateStr = GetSifName( GotIt ) 
      IF( GotIt ) THEN
        WRITE( NamesUnit,'(A)') 'Solver input file: '//TRIM(DateStr)
      END IF
            
      DateStr = FormatDate()      
      WRITE( NamesUnit,'(A)') 'File started at: '//TRIM(DateStr)

      WRITE(NamesUnit,'(A)') ' '
      WRITE(NamesUnit,'(A)') 'Variables in columns of matrix: '
      IF( LineInd /= 0 ) THEN
        i = 1
        WRITE(NamesUnit,'(I4,": ",A)') 1,'Line Marker'
      ELSE
        i = 0
      END IF
      DO No=1,NoValues 
        WRITE(NamesUnit,'(I4,": ",A)') No+i,TRIM(ValueNames(No))
      END DO
      CLOSE(NamesUnit)
    END IF
    
    !------------------------------------------------------------------------------
    ! In parallel case save the data in different files
    !------------------------------------------------------------------------------
    
    WRITE( Message,'(A)' ) 'Saving values to file: '// TRIM(ScalarsFile)
    CALL Info( Caller, Message, Level=4 )
    
    IF ( ParallelWrite ) THEN      
      IF(WriteCore) WRITE( ScalarParFile, '(A,i0)' ) TRIM(ScalarsFile)//'.', ParEnv % MyPe      
      IF( Solver % TimesVisited > 0 .OR. FileAppend) THEN 
        OPEN(NEWUNIT=ScalarsUnit, FILE=ScalarParFile,POSITION='APPEND',IOStat=istat)
      ELSE 
        OPEN(NEWUNIT=ScalarsUnit, FILE=ScalarParFile,IOSTAT=istat)
      END IF
    ELSE IF( WriteCore ) THEN 
      IF( Solver % TimesVisited > 0 .OR. FileAppend) THEN 
        OPEN(NEWUNIT=ScalarsUnit, FILE=ScalarsFile,POSITION='APPEND',IOStat=istat)
      ELSE 
        OPEN(NEWUNIT=ScalarsUnit, FILE=ScalarsFile,IOStat=istat)
      END IF
    END IF
    IF( istat /= 0) THEN
      CALL Fatal(Caller,'Could not open file for saving: '//TRIM(ScalarsFile))
    END IF
    

    IF( WriteCore ) THEN
      ! If there are multiple lines it may be a good idea to mark each by an index
      IF( LineInd /= 0) THEN
        WRITE (ScalarsUnit,'(I10)',advance='no') LineInd
      END IF
      DO No=1,NoValues-1
        IF( ValuesInteger(No) ) THEN
          IntVal = NINT( Values(No) )
          WRITE (ScalarsUnit,'(I10)',advance='no') IntVal
        ELSE
          WRITE (ScalarsUnit,'(ES22.12E3)',advance='no') Values(No)
        END IF
      END DO
      IF( ValuesInteger(NoValues) ) THEN
        IntVal = NINT( Values(No) )
        WRITE (ScalarsUnit,'(I10)') IntVal
      ELSE
        WRITE (ScalarsUnit,'(ES22.12E3)') Values(NoValues)
      END IF
      CLOSE(ScalarsUnit)
    
      !------------------------------------------------------------------------------
      ! Save comments by line in a metadata file
      !------------------------------------------------------------------------------
      IF( Solver % TimesVisited == 0 .AND. FileAppend .AND. LineInd /= 0 ) THEN      
        Message = ListGetString(Params,'Comment',GotIt)
        Name = TRIM(ScalarsFile) // '.' // TRIM("marker")
        IF( GotIt ) THEN
          OPEN(NEWUNIT=ScalarsUnit, FILE=Name,POSITION='APPEND',IOstat=istat)
          WRITE(ScalarsUnit,'(I6,A,A)') LineInd,': ',TRIM(Message)
          CLOSE(ScalarsUnit)
        END IF
      END IF
    
      !------------------------------------------------------------------------------
      ! Save data in CVS format for, e.g. Livegraph 
      !------------------------------------------------------------------------------
    
      IF(SaveCVS .AND. WriteCore ) THEN
        ! Save data as comma-separated-values (cvs-file)
        IF( Solver % TimesVisited > 0 .OR. FileAppend) THEN 
          OPEN(NEWUNIT=ScalarsUnit, FILE=TRIM(ScalarsFile)//'.csv',POSITION='APPEND')      
        ELSE 
          OPEN(NEWUNIT=ScalarsUnit, FILE=TRIM(ScalarsFile)//'.csv')
          DO No=1,NoValues-1
            WRITE (ScalarsUnit,'(A)',advance='no') TRIM(ValueNames(No))//","
          END DO
          WRITE (ScalarsUnit,'(A)') TRIM(ValueNames(NoValues))
        END IF
      
        DO No=1,NoValues-1
          WRITE (ScalarsUnit,'(ES22.12E3,A)',advance='no') Values(No),","
        END DO
        WRITE (ScalarsUnit,'(ES22.12E3)') Values(NoValues)
        CLOSE(ScalarsUnit)
      END IF

      IF(SaveJson .AND. WriteCore ) THEN
        CALL Info(Caller,'Saving data in jason format o file: '&
            //TRIM(ScalarsFile)//'.json',Level=6)

        IF(ParallelWrite) CALL Warn(Caller,'Use "Parallel Reduce=True" with JSON format!')
        IF(FileAppend) CALL Info(Caller,'Data is appended to existing file',Level=6)
    
        ! Save data in JSON format
        StartNewFile = ( Solver % TimesVisited == 0 .AND. .NOT. FileAppend)

        IF( StartNewFile ) THEN 
          OPEN(NEWUNIT=ScalarsUnit, ACCESS='stream',FORM='formatted',&
              FILE=TRIM(ScalarsFile)//'.json')

          WRITE( ScalarsUnit, '(A)' ) '{'

          DateStr = GetVersion()
          WRITE( ScalarsUnit,'(A)') '  "elmerver": "'//TRIM(DateStr)//'",'     

          DateStr = GetRevision( GotIt )
          IF( GotIt ) THEN
            WRITE( ScalarsUnit,'(A)') '  "elmerrev": "'//TRIM(DateStr)//'",'
          END IF

          DateStr = GetCompilationDate( GotIt )
          IF( GotIt ) THEN
            WRITE( ScalarsUnit,'(A)') '  "elmercompiled": "'//TRIM(DateStr)//'",'
          END IF

          DateStr = GetSifName( GotIt ) 
          IF( GotIt ) THEN
            WRITE( ScalarsUnit,'(A)') '  "siffile": "'//TRIM(DateStr)//'",'
          END IF
            
          DateStr = FormatDate()      
          WRITE( ScalarsUnit,'(A)') '  "starttime": "'//TRIM(DateStr)//'",'

          WRITE( ScalarsUnit,'(A)') '  "columns": '//I2S(NoValues)//','
                    
          WRITE( ScalarsUnit, '(A)',ADVANCE='no' ) '  "names":['
 
          DO No=1,NoValues 
            WRITE(ScalarsUnit,'(A)',ADVANCE='no') '"'//TRIM(ValueNames(No))//'"'
            IF(No<NoValues) WRITE(ScalarsUnit,'(A)',ADVANCE='no') ','
          END DO
          WRITE( ScalarsUnit, '(A)' ) '],'
          WRITE( ScalarsUnit, '(A)' ) '  "values":'
          WRITE( ScalarsUnit, '(A)' ) '  ['
          WRITE( ScalarsUnit, '(A)', ADVANCE='no' ) '    ['
        ELSE
          OPEN(NEWUNIT=ScalarsUnit, ACCESS='stream',FORM='formatted',&
              STATUS='old', POSITION='append',&
              FILE=TRIM(ScalarsFile)//'.json')

          ! Jump to the start position from previous visit
          WRITE( ScalarsUnit,'(A)', POS=jsonpos, ADVANCE='no' ) '   ,['
        END IF
                 
        DO No=1,NoValues-1
          WRITE (ScalarsUnit,'(ES22.12E3,A)',advance='no') Values(No),","
        END DO
        WRITE (ScalarsUnit,'(ES22.12E3,A)') Values(No),"]"

        ! Mark the start position for next visit
        INQUIRE(ScalarsUnit,POS=jsonpos)    
        
        WRITE( ScalarsUnit, '(A)' ) '  ]'
        WRITE( ScalarsUnit, '(A)' ) '}'

        CLOSE(ScalarsUnit)
      END IF

      
    END IF
  END IF ! of SaveFile
  
  !------------------------------------------------------------------------------
  ! Echo values if requested, this is the default if no output to file
  !------------------------------------------------------------------------------
  IF( EchoValues ) THEN
    CALL Info( Caller,'Showing computed results:')
    DO No=1,NoValues
      WRITE(Message,'(I4,": ",A,ES22.12E3)') No,TRIM(ValueNames(No)),Values(No)
      CALL Info(Caller,Message)
    END DO
  END IF

  !------------------------------------------------------------------------------
  ! Add data to the Simulation block for possible future use
  !------------------------------------------------------------------------------
  DO No=1,NoValues
    CALL ListAddConstReal(Model % Simulation,TRIM(ValueNames(No)),Values(No))
  END DO
 
  !------------------------------------------------------------------------------
  ! For consistency checks one may print out a value imitating ComputeChange
  !------------------------------------------------------------------------------
  CALL SetPseudoNorm()

  !------------------------------------------------------------------------------
  ! One may also export desired data as a cost function for FindOptimum solver
  !------------------------------------------------------------------------------
  CostInd = ListGetInteger( Params,'Cost Function Index',GotIt)
  IF( CostInd > 0 .AND. CostInd <= NoValues ) THEN
    CALL ListAddConstReal( Model % Simulation,'Cost Function',Values(CostInd) )
  END IF


  IF( ParEnv % MyPe == 0 ) THEN
    IF( ListGetLogical( Params,'Mark Failed Strategy',GotIt) ) THEN
      LineInd = ListGetInteger( Params,'Line Marker',GotIt)
      IF(.NOT. GotIt) THEN
        CALL Fatal(Caller,'Failed strategy marked requires > Line Marker <')
      END IF

      CALL Info(Caller,'Saving True marker at end')
      Name = 'FINISHED_MARKER_'//I2S(LineInd)
      i = 1
      OPEN(NEWUNIT=MarkerUnit,FILE=Name,STATUS='Unknown')
      WRITE(MarkerUnit,'(I0)') i
      CLOSE(MarkerUnit)
    END IF
  END IF

  
  DEALLOCATE( ElementNodes % x, ElementNodes % y, ElementNodes % z )

  n = 1
  n = ParallelReduction(n) 

  PrevRunInd = RunInd

  IF( PrevComm > 0 ) ParEnv % ActiveComm = PrevComm

  
  CALL Info(Caller, '-----------------------------------------', Level=7 )

  
!------------------------------------------------------------------------------

CONTAINS


  FUNCTION OperToParOperMap( LocalOper ) RESULT ( ParOper )

    CHARACTER(LEN=MAX_NAME_LEN) :: LocalOper, ParOper


    ParOper = 'none'
    IF(.NOT. ParallelReduce ) RETURN

    SELECT CASE(LocalOper)
      
    CASE('nodes','elements','dofs','sum','sum square','sum abs','int','int square','int abs','volume',&
        'potential energy', 'convective energy','diffusive energy','boundary sum','boundary dofs',&
        'boundary int','area','diffusive flux','convective flux','nans','partition checksum',&
        'partition neighbours checksum')
      ParOper = 'sum'
            
    CASE('max','max abs','boundary max','boundary max abs')
      ParOper = 'max'
      
    CASE('min','min abs','boundary min','boundary min abs')
      ParOper = 'min'

      ! These operators should already be more of less parallel
    CASE('partitions','cpu time','wall time','cpu memory','norm','nonlin change','steady state change',&
	'nonlin iter','nonlin converged','steady converged','threads',&
        'mean','mean abs','mean square','variance','deviation','int mean','int square mean', &
        'int abs mean','int variance')
      ParOper = 'none'

    CASE DEFAULT
      ParOper = ListGetString( Params,'Default Parallel Operator',GotIt)
      IF( .NOT. GotIt ) THEN
        ParOper = 'none'
        CALL Warn(Caller,'Reduction not implemented in parallel:'//TRIM(LocalOper))
      END IF      

    END SELECT

    CALL Info('OperToParOperMap',TRIM(LocalOper)//' -> '//TRIM(ParOper),Level=12)

  END FUNCTION OperToParOperMap




  SUBROUTINE AddToSaveList(Name, Val, ValueIsInteger, ParallelOperator, CheckForDuplicates )

    CHARACTER(LEN=*) :: Name
    REAL(KIND=dp) :: Val
    LOGICAL, OPTIONAL :: ValueIsInteger
    CHARACTER(LEN=*), OPTIONAL :: ParallelOperator
    LOGICAL, OPTIONAL :: CheckForDuplicates
    !------------------------------------------------------------------------
    INTEGER :: i,n
    REAL(KIND=dp) :: ParVal
    CHARACTER(LEN=MAX_NAME_LEN) :: Str, ParOper

    LOGICAL, ALLOCATABLE :: TmpValuesInteger(:)     
    REAL(KIND=dp), ALLOCATABLE :: TmpValues(:)
    CHARACTER(LEN=MAX_NAME_LEN), ALLOCATABLE :: TmpValueNames(:)

    TYPE(Variable_t), POINTER :: TargetVar
    INTEGER :: MPIOper
    INTEGER :: istat
    LOGICAL :: GotParOper, GotIt

    SAVE TmpValues, TmpValueNames

    ! For the first time allocate some space
    IF(.NOT. ALLOCATED(Values)) THEN
      n = 20
      ALLOCATE( Values(n), ValueNames(n), ValuesInteger(n), STAT=istat )
      IF( istat /= 0 ) CALL Fatal(Caller,'Memory allocation error 6') 	
      Values = 0._dp
    END IF

    n = NoValues

    ! We save always in lower case as we may be doing string comparisons etc.
    nlen = StringToLowerCase(str,Name)

    IF( PRESENT( CheckForDuplicates ) ) THEN      
      IF( CheckForDuplicates ) THEN
        DO i=1,n
          IF( TRIM(ValueNames(i)) == str(1:nlen) ) EXIT
        END DO
        IF( i<=n) THEN
          CALL Info(Caller,'Found duplicate at '//I2S(i)//' for: '//TRIM(str),Level=12)
          RETURN
        END IF
      END IF
    END IF
    
    ! If vectors are too small make some more room in a rather dummy way
    IF(n >= SIZE(Values) ) THEN
      ALLOCATE(TmpValues(n), TmpValueNames(n), TmpValuesInteger(n),STAT=istat)
      IF( istat /= 0 ) CALL Fatal(Caller,'Memory allocation error 8') 		
	
      TmpValues(1:n) = Values(1:n)
      TmpValuesInteger(1:n) = ValuesInteger(1:n)
      TmpValueNames(1:n) = ValueNames(1:n)
      DEALLOCATE(Values,ValueNames,ValuesInteger)

      ALLOCATE(Values(n+10), ValueNames(n+10), ValuesInteger(n+10), STAT=istat )
      IF( istat /= 0 ) CALL Fatal(Caller,'Memory allocation error 9') 		
      Values = 0._dp

      Values(1:n) = TmpValues(1:n)
      ValuesInteger(1:n) = TmpValuesInteger(1:n)
      ValueNames(1:n) = TmpValueNames(1:n)
      DEALLOCATE(TmpValues,TmpValueNames,TmpValuesInteger)
    END IF
    
    n = n + 1
    Values(n) = Val

    ! If we have positive or negative operator then revert that back in the save name
    IF( PosOper ) THEN
      ValueNames(n) = 'positive '//TRIM(str)
    ELSE IF( NegOper ) THEN
     ValueNames(n) = 'negative '//TRIM(str)
    ELSE
      ValueNames(n) = TRIM(str)
    END IF
      
    IF( PRESENT( ValueIsInteger ) ) THEN
      ValuesInteger(n) = ValueIsInteger
    ELSE
      ValuesInteger(n) = .FALSE.
    END IF
    NoValues = n

    IF( ParallelReduce ) THEN
      GotParOper = .FALSE.
      IF( PRESENT(ParallelOperator)) THEN
	ParOper = ParallelOperator
        GotParOper = .TRUE.
      ELSE

!------------------------------------------------------------------------------------------------
! This is to ensure that parallel operators may be applied to also other scalars than those 
! computed within this solver with the serial operators. The indexing may not be known in advance
! but may be checked by a trial run. Note that the indexing does not often coinsice with the 
! normal operators but then this loop is never reached.
!------------------------------------------------------------------------------------------------

        WRITE (Str,'(A,I0)') 'Parallel Operator ',n 
        ParOper = ListGetString(Params,TRIM(Str),GotParOper)
      END IF

      ! Note: the parallel operators for ParallelReduction have different numbering
      ! conventions that MPI_MIN, MPI_MAX, etc. 
      
      IF( GotParOper ) THEN
        SELECT CASE( ParOper )
         
        CASE('sum' )
          MPIOper = 0 
        CASE('min' )
          MPIOper = 1 
        CASE('max' )
          MPIOper = 2 
        CASE('mean' )
          MPIOper = 3
        CASE('none' )
          MPIOper = -1
          
        CASE DEFAULT
          CALL Warn(Caller,'Unknown parallel operator: '//TRIM(ParOper))
          MPIOper = -1
          
        END SELECT
        
        IF( MPIOper >= 0 ) THEN
          ParVal = ParallelReduction( Val, MPIOper )
          Values(n) = ParVal
          WRITE( ValueNames(n),'(A)') TRIM( ValueNames(n) )//' : mpi_'//TRIM(ParOper)
        END IF
        
      END IF
    END IF
  
    !------------------------------------------------------------------------------
    ! If requested, create variable of the result
    ! This is performed already here so the variable can be used 
    ! in subsequent definitions within SaveScalars.
    !------------------------------------------------------------------------------
    WRITE (Str,'(A,I0)') 'Target Variable ',n
    VariableName = ListGetString( Params, TRIM(Str), GotIt )
    IF( GotIt ) THEN
      TargetVar => VariableGet( Model % Variables, TRIM(VariableName) )
      IF(.NOT. ASSOCIATED(TargetVar)) THEN
        WrkPntr => Null()
        ALLOCATE(WrkPntr(1),STAT=istat)
	IF( istat /= 0 ) CALL Fatal(Caller,'Memory allocation error 5') 	
 
        CALL VariableAdd( Model % Variables, Mesh, Solver, &
            TRIM(VariableName), 1, WrkPntr )
        TargetVar => VariableGet( Model % Variables, TRIM(VariableName) )       
      END IF
      TargetVar % Values(1) = Values(n)
      CALL Info(Caller,'Defining: '//TRIM(VariableName)//' = '//TRIM(ValueNames(n)),Level=8)
    END IF
    
  END SUBROUTINE AddToSaveList


  ! Create table for masking nodes in statistical operators.
  !-------------------------------------------------------------------------
  SUBROUTINE CreateNodeMask( ) 

    INTEGER :: t
    TYPE(Element_t), POINTER :: Element
    TYPE(ValueList_t), POINTER :: MaskList

    IF(.NOT. MaskOper ) RETURN

    CALL Info(Caller,'Creating mask for: '//TRIM(MaskName),Level=10)
    
    IF( ALLOCATED( NodeMask ) ) THEN
      IF( SIZE( NodeMask ) /= SIZE( Var % Perm ) ) DEALLOCATE( NodeMask ) 
    END IF
   
    IF(.NOT. ALLOCATED( NodeMask ) ) THEN
      ALLOCATE( NodeMask( SIZE( Var % Perm ) ) ) 
    END IF
    NodeMask = .FALSE.
    
    DO t = 1, Mesh % NumberOfBulkElements 
      Element => Mesh % Elements(t)
      n = Element % TYPE % NumberOfNodes

      IF( ElementalVar ) THEN
        NodeIndexes => Element % DgIndexes
      ELSE        
        NodeIndexes => Element % NodeIndexes
      END IF

      ! If we are masking operators with correct (body, body force, or material) then do it here
      IF( BodyOper ) THEN        
        MaskList => GetBodyParams( Element ) 
        IF(.NOT. ASSOCIATED(MaskList)) CYCLE
      ELSE IF( BodyForceOper ) THEN
        MaskList => GetBodyForce( Element, GotIt )
        IF( .NOT. GotIt ) CYCLE
      ELSE IF( MaterialOper ) THEN
        MaskList => GetMaterial( Element, GotIt ) 
        IF(.NOT. GotIt ) CYCLE
      ELSE
        CALL Fatal('MaskNodes','Unknown mask strategy')
      END IF
      IF( ListGetLogical( MaskList, MaskName, GotIt ) ) THEN
        NodeMask(NodeIndexes(1:n)) = .TRUE.
      END IF
    END DO
    
    t = COUNT( NodeMask )    
    CALL Info(Caller,'Created mask of size: '&
        //I2S(t),Level=12)

  END SUBROUTINE CreateNodeMask




  FUNCTION VectorStatistics(Var,OperName,NodalOper) RESULT (operx)

    TYPE(Variable_t), POINTER :: Var
    CHARACTER(LEN=MAX_NAME_LEN) :: OperName
    LOGICAL, OPTIONAL :: NodalOper
    REAL(KIND=dp) :: operx

    REAL(KIND=dp) :: Minimum, Maximum, AbsMinimum, AbsMaximum, &
        Mean, sumx, sumxx, sumabsx, x, Variance2
    INTEGER :: Nonodes, i, j, k, l, NoDofs, sumi
    TYPE(NeighbourList_t), POINTER :: nlist(:)
    INTEGER, POINTER :: PPerm(:)
    
    CALL Info(Caller,'Computing operator: '//TRIM(OperName),Level=12)

    operx = 0.0_dp

    sumi = 0
    sumx = 0.0
    sumxx = 0.0
    sumabsx = 0.0

    Maximum = -HUGE(x)
    Minimum = HUGE(x)
    AbsMaximum = 0.0_dp
    AbsMinimum = HUGE(x)

    PPerm => Var % Perm
    IF( Var % TYPE == Variable_on_gauss_points .OR. &
        Var % TYPE == Variable_on_elements ) THEN
      NULLIFY( PPerm )
    END IF

    
    NoDofs = Var % Dofs
    IF(ASSOCIATED (PPerm)) THEN
      Nonodes = SIZE(PPerm) 
    ELSE
      Nonodes = SIZE(Var % Values) / NoDofs
    END IF

    
    IF( MaskOper ) THEN
      IF( NoNodes > SIZE(NodeMask) ) THEN
        CALL Info(Caller,'Decreasing operator range to size of mask: '&
            //I2S(NoNodes)//' vs. '//I2S(SIZE(NodeMask)), Level=8)
        NoNodes = SIZE(NodeMask)
      END IF
    END IF

    nlist => NULL()
    IF( IsParallel ) THEN
      IF(ASSOCIATED(Var % Solver)) THEN
        IF(ASSOCIATED(Var % Solver % Matrix)) THEN
          IF(ASSOCIATED(Var % Solver % Matrix % ParallelInfo)) THEN
            nlist => Var % Solver % Matrix % ParallelInfo % NeighbourList
          END IF
        END IF
      END IF
    END IF

    IF( PRESENT( NodalOper ) ) THEN
      IF( NodalOper ) THEN
        FirstInd = 1
        LastInd = Mesh % NumberOfNodes
      ELSE
        FirstInd = Mesh % NumberOfNodes + 1
        LastInd = SIZE( PPerm ) 
      END IF
    ELSE
      FirstInd = 1
      LastInd = NoNodes
    END IF

    
    DO i=FirstInd,LastInd 
      IF( MaskOper ) THEN
        IF( .NOT. NodeMask(i) ) CYCLE
      END IF

      j = i
      IF(ASSOCIATED(PPerm)) j = Var % Perm(i)

      IF(j > 0) THEN
        IF( IsParallel .AND. ASSOCIATED(nlist) ) THEN
          IF( ASSOCIATED( nlist(j) % Neighbours ) ) THEN
            IF( SIZE( nlist ) >= j ) THEN
              IF( nlist(j) % Neighbours(1) /= ParEnv % MyPE ) CYCLE
            ELSE
              PRINT *,'Nlist too small:',SIZE(nlist), j
            END IF
          END IF
        END IF

        IF(NoDofs <= 1) THEN
          x = Var % Values(j)
        ELSE
          x = 0.0d0
          DO l=1,NoDofs
            x = x + Var % Values(NoDofs*(j-1)+l) ** 2
          END DO
          x = SQRT(x)
        END IF

        IF( PosOper .AND. x < 0 ) CYCLE
        IF( NegOper .AND. x > 0 ) CYCLE
        
        sumi = sumi + 1
        sumx = sumx + x
        sumxx = sumxx + x*x
        sumabsx = sumabsx + ABS( x )

        Maximum = MAX(x,Maximum)
        Minimum = MIN(x,Minimum)
        IF(ABS(x) > ABS(AbsMaximum) ) AbsMaximum = x
        IF(ABS(x) < ABS(AbsMinimum) ) AbsMinimum = x
      END IF
    END DO

    ! These operators cannot be done in parallel afterwards, since they are more complex
    ! than 'min','max', or 'sum'
    IF( ParallelReduce ) THEN
      IF( OperName == 'mean' .OR. OperName == 'mean abs' .OR. &
          OperName == 'mean square' .OR. OperName == 'variance' ) THEN
        sumi = ParallelReduction(sumi)
      END IF
    END IF
      
    ! If there are no dofs avoid division by zero
    IF(sumi == 0) RETURN


    SELECT CASE(OperName)
      
    CASE ('sum')
      operx = sumx

    CASE ('sum square')
      operx = sumxx 
      
    CASE ('sum abs')
      operx = sumabsx      
      
    CASE ('max')
      operx = Maximum
      
    CASE ('min')
      operx = Minimum
      
    CASE ('range')
      operx = Maximum - Minimum

    CASE ('max abs')
      operx = AbsMaximum
      
    CASE ('min abs')
      operx = AbsMinimum
      
    CASE ('mean')
      IF(ParallelReduce) sumx = ParallelReduction(sumx)
      operx = sumx / sumi 

    CASE ('mean square')
      IF(ParallelReduce) sumxx = ParallelReduction(sumxx)
      operx = sumxx / sumi 

    CASE ('rms')
      IF(ParallelReduce) sumxx = ParallelReduction(sumxx)
      operx = SQRT( sumxx / sumi )
      
    CASE ('mean abs')
      IF(ParallelReduce) sumabsx = ParallelReduction(sumabsx)
      operx = sumabsx / sumi

    CASE('variance')      
      IF(ParallelReduce) THEN
        sumx = ParallelReduction(sumx)
        sumxx = ParallelReduction(sumxx)
      END IF
      mean = sumx / sumi
      Variance2 = sumxx/sumi-mean*mean
      IF(Variance2 > 0.0d0) operx = SQRT(Variance2) 
      
    CASE DEFAULT 
      CALL Warn(Caller,'Unknown statistical operator!')

    END SELECT
          
    CALL Info(Caller,'Finished computing operator',Level=12)

  END FUNCTION VectorStatistics



  
!------------------------------------------------------------------------------

  FUNCTION VectorMeanDeviation(Var,OperName) RESULT (Deviation)

    TYPE(Variable_t), POINTER :: Var
    CHARACTER(LEN=MAX_NAME_LEN) :: OperName
    REAL(KIND=dp) :: Mean, Deviation
    REAL(KIND=dp) :: sumx, sumdx, x, dx
    INTEGER :: Nonodes, i, j, k, NoDofs, sumi
    INTEGER, POINTER :: PPerm(:)
    
    NoDofs = Var % Dofs

    PPerm => Var % Perm
    IF( Var % TYPE == Variable_on_gauss_points .OR. &
        Var % TYPE == Variable_on_elements ) THEN
      NULLIFY( PPerm )
    END IF

    IF(ASSOCIATED (PPerm)) THEN
      Nonodes = SIZE(PPerm) 
    ELSE
      Nonodes = SIZE(Var % Values) / NoDofs
    END IF

    IF( MaskOper ) THEN
      IF( NoNodes > SIZE(NodeMask) ) THEN
        CALL Info(Caller,'Decreasing operator range to size of mask: '&
            //I2S(NoNodes)//' vs. '//I2S(SIZE(NodeMask)) )
        NoNodes = SIZE(NodeMask)
      END IF
    END IF

    sumi = 0
    sumx = 0.0
    Deviation = 0.0

    DO i=1,Nonodes
      IF( MaskOper ) THEN
        IF( .NOT. NodeMask(i) ) CYCLE
      END IF

      j = i
      IF( IsParallel ) THEN
        IF( Mesh % ParallelInfo % NeighbourList(j) % Neighbours(1) /= ParEnv % MyPE ) CYCLE
      END IF

      IF(ASSOCIATED(PPerm)) j = PPerm(i)
      IF(j > 0) THEN
        IF(NoDofs <= 1) THEN
          x = Var % Values(j)
        ELSE
          x = 0.0d0
          DO k=1,NoDofs
            x = x + Var % Values(NoDofs*(j-1)+k) ** 2
          END DO
          x = SQRT(x)
        END IF

        IF( PosOper .AND. x < 0 ) CYCLE
        IF( NegOper .AND. x > 0 ) CYCLE

        sumi = sumi + 1
        sumx = sumx + x
      END IF
    END DO

    IF( ParallelReduce ) sumi = ParallelReduction(sumi)   
    
    IF(sumi == 0) RETURN

    IF( ParallelReduce ) sumx = ParallelReduction(sumx)
    Mean = sumx / sumi
    sumi = 0
    sumdx = 0.0
    DO i=1,Nonodes
      IF( MaskOper ) THEN
        IF( .NOT. NodeMask(i) ) CYCLE
      END IF

      j = i
      IF( IsParallel ) THEN
        IF( Mesh % ParallelInfo % NeighbourList(j) % Neighbours(1) /= ParEnv % MyPE ) CYCLE
      END IF

      IF(ASSOCIATED(PPerm)) j = PPerm(i)
      IF(j > 0) THEN
        IF(NoDofs <= 1) THEN
          x = Var % Values(j)
        ELSE
          x = 0.0d0
          DO k=1,NoDofs
            x = x + Var % Values(NoDofs*(j-1)+k) ** 2
          END DO
          x = SQRT(x)
        END IF

        IF( PosOper .AND. x < 0 ) CYCLE
        IF( NegOper .AND. x > 0 ) CYCLE

        dx = ABS(x-Mean)
        sumi = sumi + 1
        sumdx = sumdx + dx
      END IF
    END DO
    IF(ParallelReduce) sumdx = ParallelReduction(sumdx)
    Deviation = sumdx / sumi

  END FUNCTION VectorMeanDeviation

!------------------------------------------------------------------------------

  FUNCTION BulkIntegrals(Var, OperName, GotCoeff, CoeffName) RESULT (operx)
    TYPE(Variable_t), POINTER :: Var
    CHARACTER(LEN=MAX_NAME_LEN) :: OperName, CoeffName
    LOGICAL :: GotCoeff
    REAL(KIND=dp) :: operx, vol
    
    INTEGER :: t, hits
    TYPE(Element_t), POINTER :: Element
    INTEGER, POINTER :: NodeIndexes(:), PermIndexes(:)
    REAL(KIND=dp) :: Basis(Model % MaxElementNodes)
    REAL(KIND=dp) :: EnergyTensor(3,3,Model % MaxElementNodes),&
        EnergyCoeff(Model % MaxElementNodes), ElemVals(Model % MaxElementNodes) 
    REAL(KIND=dp) :: DetJ,U,V,W,S,A,L,C(3,3),x,y,z,Vals(3),uvw(3)
    REAL(KIND=dp) :: func, coeff, integral1, integral2, Grad(3), CoeffGrad(3)
    REAL(KIND=DP), POINTER :: Pwrk(:,:,:) => Null()
    LOGICAL :: Stat
    TYPE(ValueList_t), POINTER :: MaskList, Material

    INTEGER :: i,j,k,p,q,DIM,NoDofs,No
    
    TYPE(GaussIntegrationPoints_t) :: IntegStuff
    LOGICAL :: DiffEnergy
    
    hits = 0
    integral1 = 0._dp
    integral2 = 0._dp
    vol = 0._dp
    EnergyCoeff = 1.0d0

    NoDofs = Var % Dofs

    DIM = CoordinateSystemDimension()

    DiffEnergy = ( OperName == 'diffusive energy' )  

    
    DO t = 1, Mesh % NumberOfBulkElements + Mesh % NumberOfBoundaryElements

      IF(t == Mesh % NumberOfBulkElements + 1 .AND. hits > 0) GOTO 10

      Element => Mesh % Elements(t)
      Model % CurrentElement => Mesh % Elements(t)
      n = Element % TYPE % NumberOfNodes

      IF ( Element % TYPE % ElementCode == 101 ) CYCLE
      
      IF( ElementalVar ) THEN
        PermIndexes => Element % DgIndexes
      ELSE
        PermIndexes => Element % NodeIndexes
      END IF

      IF(ASSOCIATED(Var % Perm)) THEN
        IF ( ANY(Var % Perm(PermIndexes) == 0 ) ) CYCLE      
      END IF
      hits = hits + 1
        
      NodeIndexes => Element % NodeIndexes 
      CALL CopyElementNodesFromMesh( ElementNodes, Mesh, n, NodeIndexes)

      IF( PassiveCoordinate > 0 ) THEN
        SELECT CASE( PassiveCoordinate ) 
        CASE(1) 
          ElementNodes % x(1:n) = 0.0_dp 
        CASE(2)           
          ElementNodes % y(1:n) = 0.0_dp
        CASE(3) 
          ElementNodes % z(1:n) = 0.0_dp
        END SELECT
      END IF

      ! If we are masking operators with correct (body, body force, or material) then do it here
      IF( BodyOper ) THEN        
        MaskList => GetBodyParams( Element ) 
      ELSE IF( BodyForceOper ) THEN
        MaskList => GetBodyForce( Element, GotIt )
        IF( .NOT. GotIt ) CYCLE
      ELSE IF( MaterialOper ) THEN
        MaskList => GetMaterial( Element, GotIt ) 
        IF(.NOT. GotIt ) CYCLE
      END IF
      IF( MaskOper ) THEN
        IF( .NOT. ListGetLogical( MaskList, MaskName, GotIt ) ) CYCLE
      END IF

      IF( GotCoeff ) Material => GetMaterial( Element, GotIt ) 

      IF( DiffEnergy ) THEN
        EnergyTensor = 0.0d0
        IF(GotCoeff) THEN
          CALL ListGetRealArray( Material, TRIM(CoeffName), Pwrk, n, NodeIndexes )

          IF ( SIZE(Pwrk,1) == 1 ) THEN
            DO i=1,3
              EnergyTensor( i,i,1:n ) = Pwrk( 1,1,1:n )
            END DO
          ELSE IF ( SIZE(Pwrk,2) == 1 ) THEN
            DO i=1,MIN(3,SIZE(Pwrk,1))
              EnergyTensor(i,i,1:n) = Pwrk(i,1,1:n)
            END DO
          ELSE
            DO i=1,MIN(3,SIZE(Pwrk,1))
              DO j=1,MIN(3,SIZE(Pwrk,2))
                EnergyTensor( i,j,1:n ) = Pwrk(i,j,1:n)
              END DO
            END DO
          END IF
        ELSE 
          DO i=1,3          
            EnergyTensor(i,i,1:n) = 1.0d0
          END DO
        END IF
      ELSE
        IF(GotCoeff) THEN
          EnergyCoeff = ListGetReal( Material, TRIM(CoeffName), n, NodeIndexes(1:n) )
        ELSE
          EnergyCoeff(1:n) = 1.0d0
        END IF
      END IF

!------------------------------------------------------------------------------
!    Numerical integration
!------------------------------------------------------------------------------
      IntegStuff = GaussPoints( Element )
            
      DO i=1,IntegStuff % n
        U = IntegStuff % u(i)
        V = IntegStuff % v(i)
        W = IntegStuff % w(i)
!------------------------------------------------------------------------------
!        Basis function values & derivatives at the integration point
!------------------------------------------------------------------------------
        stat = ElementInfo( Element,ElementNodes,U,V,W,DetJ,Basis)
!------------------------------------------------------------------------------
!      Coordinatesystem dependent info
!------------------------------------------------------------------------------
        
        s = DetJ * IntegStuff % s(i)
        IF ( CurrentCoordinateSystem() /= Cartesian ) THEN
          x = SUM( ElementNodes % x(1:n)*Basis(1:n) )
          s = 2._dp * PI * x * s 
        END IF

        coeff = SUM( EnergyCoeff(1:n) * Basis(1:n))
        vol =  coeff * vol + S

        uvw(1)=u;uvw(2)=v;uvw(3)=w;

        No = 0
        IF( DiffEnergy ) THEN
          CALL EvaluateVariableAtGivenPoint(No,Vals,Mesh,Var,Element=Element,LocalCoord=uvw,&
              DoGrad=.TRUE.)
          Grad(1:dim) = Vals(1:dim)
        ELSE
          CALL EvaluateVariableAtGivenPoint(No,Vals,Mesh,Var,Element=Element,LocalCoord=uvw)
        END IF

        IF(No > 1) THEN
          func = SQRT(SUM(Vals(1:No)**2))
        ELSE
          func = Vals(1)
        END IF
        
          
        SELECT CASE(OperName)

        CASE ('volume')
          integral1 = integral1 + coeff * S

        CASE ('int','int mean')
          !          func = SUM( ElemVals(1:n) * Basis(1:n) )
          IF( PosOper ) func = MAX( 0.0_dp, func )
          IF( NegOper ) func = MIN( 0.0_dp, func ) 
          integral1 = integral1 + S * coeff * func 

        CASE ('int square','int square mean','int rms')
          !          func = SUM( ElemVals(1:n) * Basis(1:n) )
          IF( PosOper ) func = MAX( 0.0_dp, func )
          IF( NegOper ) func = MIN( 0.0_dp, func ) 
          integral1 = integral1 + S * coeff * func**2 

        CASE ('int abs','int abs mean')
          !          func = ABS( SUM( ElemVals(1:n) * Basis(1:n) ) )
          func = ABS(func)
          IF( PosOper ) func = MAX( 0.0_dp, func )
          IF( NegOper ) func = MIN( 0.0_dp, func ) 
          integral1 = integral1 + S * coeff * func 

        CASE ('int variance')
          !          func = SUM( ElemVals(1:n) * Basis(1:n) )
          IF( PosOper ) func = MAX( 0.0_dp, func )
          IF( NegOper ) func = MIN( 0.0_dp, func ) 
          integral1 = integral1 + S * coeff * func 
          integral2 = integral2 + S * coeff * func**2 

        CASE ('diffusive energy')
          CoeffGrad = 0.0d0
          DO j = 1, DIM
            DO k = 1, DIM
              CoeffGrad(j) = CoeffGrad(j) + SUM( EnergyTensor(j,k,1:n) * Basis(1:n) ) * Grad(k)
            END DO
          END DO

          integral1 = integral1 + s * SUM( Grad(1:DIM) * CoeffGrad(1:DIM) )

        CASE ('convective energy')
          IF( PosOper ) func = MAX( 0.0_dp, func )
          IF( NegOper ) func = MIN( 0.0_dp, func ) 

          IF(NoDofs == 1) THEN
            integral1 = integral1 + s * coeff * func**2
          ELSE
            j=MIN(DIM,NoDofs)
            integral1 = integral1 + s * coeff * SUM(Grad(1:j)**2)
          END IF

        CASE ('potential energy')

          IF( PosOper ) func = MAX( 0.0_dp, func )
          IF( NegOper ) func = MIN( 0.0_dp, func ) 

          integral1 = integral1 + s * coeff * func

        CASE DEFAULT 
          CALL Warn(Caller,'Unknown statistical operator!')

        END SELECT

      END DO

    END DO
    
10  CONTINUE 
    
    operx = 0.0d0
    IF(ParallelReduce) THEN
      IF( OperName == 'int mean' .OR. OperName == 'int square mean' .OR. &
          OperName == 'int abs mean' .OR. OperName == 'int variance' ) THEN
        hits = ParallelReduction(hits)
      END IF
    END IF

    
    IF(hits == 0) RETURN

    SELECT CASE(OperName)

    CASE ('volume','int','int abs','int square')
      operx = integral1

    CASE ('int mean','int square mean','int abs mean')
      IF( ParallelReduce ) THEN
        integral1 = ParallelReduction(integral1)
        vol = ParallelReduction(vol)
      END IF
      operx = integral1 / vol        

    CASE ('int rms')
      IF( ParallelReduce ) THEN
        integral1 = ParallelReduction(integral1)
        vol = ParallelReduction(vol)
      END IF
      operx = SQRT( integral1 / vol ) 

    CASE ('int variance')
      IF( ParallelReduce ) THEN
        integral1 = ParallelReduction(integral1)
        integral2 = ParallelReduction(integral2)
        vol = ParallelReduction(vol)
      END IF
      operx = SQRT(integral2/vol-(integral1/vol)**2)

    CASE ('diffusive energy','convective energy')
      operx = 0.5d0 * integral1

    CASE ('potential energy')
      operx = integral1

    END SELECT


  END FUNCTION BulkIntegrals
!------------------------------------------------------------------------------


  SUBROUTINE BoundaryIntegrals(Var, OperName, GotCoeff, &
      CoeffName, fluxes, areas, fluxescomputed)

    TYPE(Variable_t), POINTER :: Var
    CHARACTER(LEN=MAX_NAME_LEN) :: OperName, CoeffName
    LOGICAL :: GotCoeff
    REAL(KIND=dp) :: fluxes(:), areas(:)
    INTEGER :: fluxescomputed(:)
    
    TYPE(Variable_t), POINTER :: VeloVar
    INTEGER :: t, FluxBody, LBody, RBody, NActive
    TYPE(Element_t), POINTER :: Element, Parent    
    TYPE(ValueList_t), POINTER :: Material, BCVal
    REAL(KIND=dp) :: Basis(Model % MaxElementNodes),ParentBasis(Model % MaxElementNodes),&
        EnergyTensor(3,3,Model % MaxElementNodes),EnergyCoeff(Model % MaxElementNodes)
    REAL(KIND=dp) :: DetJ,U,V,W,up,vp,wp,S,A,L,C(3,3),x,y,z
    REAL(KIND=dp) :: func, coeff, Normal(3), Flow(3), Grad(3), flux
    REAL(KIND=DP), POINTER :: Pwrk(:,:,:) => Null()
    INTEGER, POINTER :: ParentIndexes(:), PermIndexes(:)
    REAL(KIND=dp) :: LocalVectorSolution(3,35), LocalVeloSolution(3,35)

    LOGICAL :: Stat, Permutated    
    INTEGER :: i,j,k,p,q,DIM,bc,NoDofs,pn,hits,istat
    TYPE(GaussIntegrationPoints_t) :: IntegStuff
    TYPE(Nodes_t) :: ParentNodes

    n = Model % MaxElementNodes

    DIM = CoordinateSystemDimension()

    ALLOCATE(ParentNodes % x(n), ParentNodes % y(n), ParentNodes % z(n), PermIndexes(n), STAT=istat )
    IF( istat /= 0 ) CALL Fatal('BoundaryIntegrals','Memory allocation error') 


    IF( SaveFluxRange ) THEN
      Minimum = HUGE(minimum)
      Maximum = -HUGE(maximum)
    END IF
      
    IF( OperName == 'area' ) THEN
      NoDofs = 1
      Permutated = .FALSE.
    ELSE
      NoDofs = Var % Dofs
      Permutated = ASSOCIATED(Var % Perm)
    END IF

    DO i=1,3          
      EnergyTensor(i,i,:) = 1.0d0
    END DO

    SELECT CASE(OperName)
      
      CASE('diffusive flux') 
      IF(NoDofs /= 1) THEN
        CALL Fatal(Caller,'diffusive flux & NoDofs /= 1?')
      END IF

      CASE ('convective flux')
      IF( NoDofs < DIM) THEN
        IF( NoDofs == 1 ) THEN
          VeloVar => VariableGet( Mesh % Variables,'Flow Solution')
          IF( .NOT. ASSOCIATED( VeloVar ) ) THEN
            CALL Fatal(Caller,'For convective flux of a scalar we need the convective field!')
          END IF
        ELSE
          CALL Warn(Caller,'convective flux & NoDofs < DIM?')
          GOTO 100
        END IF
      END IF
      
      CASE ('area','boundary int','boundary int mean')

    CASE DEFAULT 
      CALL Warn(Caller,'Unknown statistical OPERATOR')

    END SELECT

    fluxes = 0.0_dp
    areas = 0.0_dp
    coeff = 1.0_dp


    NActive = GetNOFBoundaryElements()

    DO t = 1, NActive 

      Element => GetBoundaryElement(t, Var % Solver)
      BCVal => GetBC()
      IF (.NOT. ASSOCIATED(BCVal) ) CYCLE

      IF (NoDOFs > 1) THEN
        CALL GetVectorLocalSolution(LocalVectorSolution, &
            UElement=Element, USolver=Var % Solver, UVariable=Var)
      ELSE
        CALL GetScalarLocalSolution(LocalVectorSolution(1,:), &
            UElement=Element, USolver=Var % Solver, UVariable=Var)
        IF( ComponentVar ) THEN
          IF( ASSOCIATED( Var2 ) ) THEN
            CALL GetScalarLocalSolution(LocalVectorSolution(2,:), &
                UElement=Element, USolver=Var % Solver, UVariable=Var2)
            NoDofs = 2
          ELSE
            LocalVectorSolution(2,:) = 0.0_dp
          END IF
          IF( ASSOCIATED( Var3 ) ) THEN
            CALL GetScalarLocalSolution(LocalVectorSolution(3,:), &
                UElement=Element, USolver=Var % Solver, UVariable=Var3)
            NoDofs = 3
          ELSE
            LocalVectorSolution(3,:) = 0.0_dp
          END IF
        END IF
      END IF

      Model % CurrentElement => Element

      IF ( Element % TYPE % ElementCode == 101 ) CYCLE

      n = Element % TYPE % NumberOfNodes
      NodeIndexes => Element % NodeIndexes

      IF(Permutated) THEN
        PermIndexes(1:n) = Var % Perm(NodeIndexes(1:n))
        IF (ANY( PermIndexes(1:n) == 0)) CYCLE
      ELSE
        PermIndexes(1:n) = NodeIndexes(1:n)        
      END IF


      DO bc=1, Model % NumberOfBCs

        IF ( Model % BCs(bc) % Tag /= Element % BoundaryInfo % Constraint ) CYCLE

        IF(.NOT. ListGetLogical(Model % BCs(bc) % Values,'Flux Integrate',gotIt ) .AND. &
            .NOT. ListGetLogical(Model % BCs(bc) % Values, MaskName, gotIt ) ) CYCLE
                 
        CALL CopyElementNodesFromMesh( ElementNodes, Mesh, n, NodeIndexes)
                 
        IF( PassiveCoordinate > 0 ) THEN
          SELECT CASE( PassiveCoordinate ) 
          CASE(1) 
            ElementNodes % x(1:n) = 0.0_dp 
          CASE(2)           
            ElementNodes % y(1:n) = 0.0_dp
          CASE(3) 
            ElementNodes % z(1:n) = 0.0_dp
          END SELECT
        END IF

        
        SELECT CASE(OperName)
          
        CASE('diffusive flux') 
          
          FluxBody = ListGetInteger( Model % BCs(bc) % Values,'Flux Integrate Body', gotIt ) 
          IF ( GotIt ) THEN
            IF ( ASSOCIATED( Element  % BoundaryInfo % Left ) ) &
                Lbody = Element % BoundaryInfo % Left % BodyId
            IF ( ASSOCIATED( Element  % BoundaryInfo % Right ) ) &
                Rbody = Element % BoundaryInfo % Right % BodyId
            
            IF ( Lbody == FluxBody ) THEN
              Parent => Element % BoundaryInfo % Left
            ELSEIF ( Rbody == FluxBody ) THEN
              Parent => Element % BoundaryInfo % Right
            ELSE
              CALL Fatal( Caller,'No such flux integrate body on bc '&
                  //I2S(Element % BoundaryInfo % Constraint))
            END IF
          ELSE        
            Parent => ELement % BoundaryInfo % Left
            stat = ASSOCIATED( Parent )
            
            IF(Permutated) THEN
              IF(stat) stat = ALL(Var % Perm(Parent % NodeIndexes) > 0)              
              IF ( .NOT. stat ) THEN
                Parent => ELement % BoundaryInfo % Right              
                stat = ASSOCIATED( Parent )
                IF(stat) stat = ALL(Var % Perm(Parent % NodeIndexes) > 0)
              END IF
            END IF
            IF ( .NOT. stat )  CALL Fatal( Caller,&
                'No solution available for specified boundary' )
          END IF                   
          
          ParentIndexes => Parent % NodeIndexes
          i = ListGetInteger( Model % Bodies(Parent % BodyId) % Values, 'Material', &
              minv=1, maxv=Model % NumberOFMaterials )
          Material => Model % Materials(i) % Values

          EnergyTensor = 0._dp

          IF(GotCoeff) THEN
            CALL ListGetRealArray( Material, CoeffName, Pwrk, n, Element % NodeIndexes )

            IF ( SIZE(Pwrk,1) == 1 ) THEN
              DO i=1,3
                EnergyTensor( i,i,1:n ) = Pwrk( 1,1,1:n )
              END DO
            ELSE IF ( SIZE(Pwrk,2) == 1 ) THEN
              DO i=1,MIN(3,SIZE(Pwrk,1))
                EnergyTensor(i,i,1:n) = Pwrk(i,1,1:n)
              END DO
            ELSE
              DO i=1,MIN(3,SIZE(Pwrk,1))
                DO j=1,MIN(3,SIZE(Pwrk,2))
                  EnergyTensor( i,j,1:n ) = Pwrk(i,j,1:n)
                END DO
              END DO
            END IF
          ELSE
            DO i=1,3          
              EnergyTensor(i,i,1:n) = 1.0_dp
            END DO
          END IF
          EnergyCoeff(1:n) = 1.0_dp
          fluxescomputed(bc) = fluxescomputed(bc) + 1


        CASE ('convective flux')
          
          FluxBody = ListGetInteger( Model % BCs(bc) % Values, &
              'Flux Integrate Body', gotIt ) 
          IF ( GotIt ) THEN
            IF ( ASSOCIATED( Element  % BoundaryInfo % Left ) ) &
                Lbody = Element % BoundaryInfo % Left % BodyId
            IF ( ASSOCIATED( Element  % BoundaryInfo % Right ) ) &
                Rbody = Element % BoundaryInfo % Right % BodyId
            
            IF ( Lbody == FluxBody ) THEN
              Parent => Element % BoundaryInfo % Left
            ELSEIF ( Rbody == FluxBody ) THEN
              Parent => Element % BoundaryInfo % Right
            ELSE
              WRITE( Message, * ) 'No such flux integrate body on bc ', &
                  Element % BoundaryInfo % Constraint
              CALL Fatal( Caller, Message )
            END IF
          ELSE        
            Parent => ELement % BoundaryInfo % Left
            stat = ASSOCIATED( Parent )
            
            IF(Permutated) THEN
              IF(stat) stat = ALL(Var % Perm(Parent % NodeIndexes) > 0)
              
              IF ( .NOT. stat ) THEN
                Parent => ELement % BoundaryInfo % Right              
                stat = ASSOCIATED( Parent )
                IF(stat) stat = ALL(Var % Perm(Parent % NodeIndexes) > 0)
              END IF
            END IF
            IF ( .NOT. stat )  CALL Fatal( Caller,&
                'No solution available for specified boundary' )
          END IF                   

          IF( NoDofs == 1 ) THEN
            CALL GetVectorLocalSolution(LocalVeloSolution, &
                UElement=Element, USolver=VeloVar % Solver, UVariable=VeloVar)

          END IF
          
          pn = Parent % TYPE % NumberOfNodes
          ParentIndexes => Parent % NodeIndexes
          i = ListGetInteger( Model % Bodies(Parent % BodyId) % Values, 'Material', &
              minv=1, maxv=Model % NumberOFMaterials )
          Material => Model % Materials(i) % Values

          CALL CopyElementNodesFromMesh( ParentNodes, Mesh, pn, ParentIndexes)

          IF(GotCoeff) THEN
            EnergyCoeff(1:n) = ListGetReal( Material, CoeffName, n, NodeIndexes )
          ELSE
            EnergyCoeff(1:n) = 1.0d0
          END IF          
          fluxescomputed(bc) = fluxescomputed(bc) + 1

        CASE ('area','boundary int','boundary int mean')
           IF(GotCoeff) THEN
             i = ListGetInteger( Model % Bodies(Element % BodyId) % Values, 'Material', &
                 minv=1, maxv=Model % NumberOFMaterials )
             Material => Model % Materials(i) % Values
             EnergyCoeff(1:n) = ListGetReal( Material, CoeffName, n, NodeIndexes )
          ELSE
            EnergyCoeff(1:n) = 1.0d0
          END IF
          fluxescomputed(bc) = fluxescomputed(bc) + 1


        END SELECT

!------------------------------------------------------------------------------
!    Numerical integration
!------------------------------------------------------------------------------
        IntegStuff = GaussPoints( Element )
        
        DO i=1,IntegStuff % n
          U = IntegStuff % u(i)
          V = IntegStuff % v(i)
          W = IntegStuff % w(i)
!------------------------------------------------------------------------------
!        Basis function values & derivatives at the integration point
!------------------------------------------------------------------------------
          stat = ElementInfo( Element,ElementNodes,U,V,W,DetJ,Basis)
!------------------------------------------------------------------------------
!      Coordinatesystem dependent info
!------------------------------------------------------------------------------
          s = DetJ * IntegStuff % s(i)
          IF ( CurrentCoordinateSystem() /= Cartesian ) THEN
            x = SUM( ElementNodes % x(1:n)*Basis(1:n) )
            s = 2._dp * PI * x * s
          END IF
          
          
          SELECT CASE(OperName)
            
          CASE ('diffusive flux')

            BLOCK
              REAL(KIND=dp) :: uvw(3)
              uvw(1) = u
              uvw(2) = v
              uvw(3) = w
              Grad = 0.0_dp
              k = 0
              ! This can use p-element basis etc. when evaluating the gradient in the parent element.
              CALL EvaluateVariableAtGivenPoint(k,Grad,Mesh,Var,Element=Element,LocalCoord=uvw,&
                  LocalBasis=Basis,DoGrad=.TRUE.,Parent=Parent)              
            END BLOCK
                        
            Flow = 0.0d0
            DO j = 1, DIM
              DO k = 1, DIM
                Flow(j) = Flow(j) + SUM( EnergyTensor(j,k,1:n) * Basis(1:n) ) * Grad(k)
              END DO
            END DO
            
            Normal = NormalVector( Element,ElementNodes,U,V,.TRUE. )
            flux = SUM(Normal(1:DIM) * Flow(1:DIM))

            IF( SaveFluxRange ) THEN
              Minimum = MIN(flux,Minimum)
              Maximum = MAX(flux,Maximum)
            END IF

            IF( PosOper ) flux = MAX( flux, 0.0_dp )
            IF( NegOper ) flux = MIN( flux, 0.0_dp ) 
                         
            fluxes(bc) = fluxes(bc) + s * flux
            
          CASE ('convective flux')
            
            Normal = NormalVector( Element,ElementNodes,u,v,.TRUE. )
            coeff = SUM( EnergyCoeff(1:n) * Basis(1:n))
            
            IF(NoDofs == 1) THEN              
              DO j=1,DIM
                Flow(j) = coeff * SUM(LocalVectorSolution(1,1:n) * LocalVeloSolution(j,1:n)*Basis(1:n))
              END DO             
            ELSE
              DO j=1,DIM
                Flow(j) = coeff * SUM(LocalVectorSolution(j,1:n)*Basis(1:n))
              END DO
            END IF

            flux = SUM( Normal(1:dim) * Flow(1:dim) ) 
            
            IF( SaveFluxRange ) THEN
              Minimum = MIN(flux,Minimum)
              Maximum = MAX(flux,Maximum)
            END IF

            ! Enable positive and negative flux integrals
            IF( PosOper ) flux = MAX( flux, 0.0_dp )
            IF( NegOper ) flux = MIN( flux, 0.0_dp ) 
          
            fluxes(bc) = fluxes(bc) + s * flux
              
          CASE ('boundary int','boundary int mean')
            
            coeff = SUM( EnergyCoeff(1:n) * Basis(1:n))
            
            IF(NoDofs == 1) THEN
              func = SUM( LocalVectorSolution(1,1:n) * Basis(1:n) )
              flux = coeff * func 
              IF( PosOper ) flux = MAX( flux, 0.0_dp )
              IF( NegOper ) flux = MIN( flux, 0.0_dp ) 
              fluxes(bc) = fluxes(bc) + s * flux
            ELSE 
              flux = 0.0_dp
              DO j=1,NoDofs
                flux = flux + SUM( Var % Values(NoDofs*(PermIndexes(1:n)-1)+j) * Basis(1:n) )**2
              END DO
              flux = coeff * SQRT( flux )
              IF( PosOper ) flux = MAX( flux, 0.0_dp )
              IF( NegOper ) flux = MIN( flux, 0.0_dp ) 
              fluxes(bc) = fluxes(bc) + s * flux
            END IF
     
          CASE ('area')
            coeff = SUM( EnergyCoeff(1:n) * Basis(1:n))
            fluxes(bc) = fluxes(bc) + s * coeff 
            
         END SELECT
          
         areas(bc) = areas(bc) + s * coeff

        END DO

      END DO

    END DO

100 CONTINUE

    DEALLOCATE( ParentNodes % x, ParentNodes % y, ParentNodes % z, PermIndexes )

  END SUBROUTINE BoundaryIntegrals


!------------------------------------------------------------------------------
! Take nodal sum over given boundaries
!------------------------------------------------------------------------------

  SUBROUTINE BoundaryStatistics(Var, OperName, GotCoeff, &
      CoeffName, fluxes, fluxescomputed)

    TYPE(Variable_t), POINTER :: Var
    CHARACTER(LEN=MAX_NAME_LEN) :: OperName, CoeffName
    LOGICAL :: GotCoeff, FindMinMax
    REAL(KIND=dp) :: fluxes(:), val
    INTEGER :: fluxescomputed(:)
    LOGICAL, ALLOCATABLE :: nodescomputed(:)    
    TYPE(Element_t), POINTER :: Element, Parent    
    TYPE(ValueList_t), POINTER :: Material
    LOGICAL :: Stat, Permutated, NodalVar    
    INTEGER :: i,j,j2,k,p,q,t,DIM,bc,n,nd,hits,istat
    INTEGER :: tmpDofs
    INTEGER, TARGET :: Indexes(100)

    IF( ASSOCIATED( Var % Perm ) ) THEN
      n = SIZE( Var % Perm )
    ELSE
      n = Mesh % NumberOfNodes
    END IF
    ALLOCATE(NodesComputed(n),STAT=istat)
    IF( istat /= 0 ) CALL Fatal('BoundaryStatistics','Memory allocation error') 
	
    NodesComputed = .FALSE.
    hits = 0
    
    NoDofs = Var % Dofs
    Permutated = ASSOCIATED(Var % Perm)
    FindMinMax = .FALSE.
    
    SELECT CASE(OperName)
      
    CASE('boundary sum','boundary dofs','boundary mean') 
      fluxes = 0.0_dp

    CASE('boundary min','boundary min abs')
      fluxes = HUGE( val )
      FindMinMax = .TRUE.
      
    CASE('boundary max','boundary max abs') 
      fluxes = -HUGE( val )
      FindMinMax = .TRUE.
      
    CASE DEFAULT 
      CALL Warn(Caller,'Unknown statistical operator')
      
    END SELECT


    DO t = Mesh % NumberOfBulkElements+1, Mesh % NumberOfBulkElements &
        + Mesh % NumberOfBoundaryElements

      Element => Mesh % Elements(t)
      Model % CurrentElement => Mesh % Elements(t)

      n = Element % TYPE % NumberOfNodes
      IF( Var % TYPE == variable_on_nodes_on_elements )  THEN
        NodeIndexes => Element % NodeIndexes
        nd = n
        NodalVar = .TRUE.
      ELSE
        ! Do we have p-elements
        IF( ASSOCIATED( Var % Solver ) ) THEN
          nd = mGetElementDOFs( Indexes, Element, Var % Solver)
        ELSE
          nd = mGetElementDOFs( Indexes, Element )
        END IF
        NodeIndexes => Indexes
        NodalVar = .FALSE.
      END IF
        
      DO bc=1, Model % NumberOfBCs

        IF ( Model % BCs(bc) % Tag /= Element % BoundaryInfo % Constraint ) CYCLE
        IF(.NOT. ListGetLogical(Model % BCs(bc) % Values, MaskName, gotIt ) ) CYCLE

        hits = hits + 1
        
        DO i=1,nd
          j = NodeIndexes(i)

          j2 = j
          IF( Permutated ) j2 = Var % Perm(j)
          IF( j2 == 0) CYCLE            
          IF( nodescomputed(j2) ) CYCLE
          

          IF( .NOT. FindMinMax .AND. IsParallel ) THEN
            IF(ASSOCIATED( Var % Solver) ) THEN
              IF( ASSOCIATED( Var % SOlver % Matrix ) ) THEN
                tmpDofs=Var % Solver % Variable % Dofs
                IF( Var % Solver % Matrix % ParallelInfo % NeighbourList(tmpDofs*(j2-1)+1) % Neighbours(1) /= ParEnv % MyPE ) CYCLE
              END IF
            END IF
          END IF



          IF( NoDofs == 1 ) THEN
            val = Var % Values(j2)
          ELSE
            val = 0.0_dp
            DO k=1,NoDofs
              val = val + Var % Values(NoDofs*(j2-1)+k)
            END DO
            val = SQRT( val )
          END IF
            
          
          IF(FindMinMax) THEN
            SELECT CASE(OperName)              

            CASE('boundary min')
              fluxes(1) = MIN( val, fluxes(1) )

            CASE('boundary max')
              fluxes(1) = MAX( val, fluxes(1) )

            CASE('boundary min abs')
              IF(ABS(fluxes(1)) < ABS( val )) fluxes(1) = val

            CASE('boundary max abs')
              IF(ABS(fluxes(1)) > ABS( val )) fluxes(1) = val

            END SELECT
          ELSE
            fluxes(1) = fluxes(1) + val
          END IF

          nodescomputed(j2) = .TRUE.         
          fluxescomputed(1) = fluxescomputed(1) + 1          
        END DO        
      END DO
    END DO


    SELECT CASE(OperName)
      
    CASE('boundary dofs') 
      fluxes = 1.0 * fluxescomputed
      
    END SELECT

  END SUBROUTINE BoundaryStatistics
  

!------------------------------------------------------------------------------


  SUBROUTINE PolylineIntegrals(Var, OperName, GotCoeff, &
      CoeffName, fluxes, areas, fluxescomputed)

    TYPE(Variable_t), POINTER :: Var
    CHARACTER(LEN=MAX_NAME_LEN) :: OperName, CoeffName
    LOGICAL :: GotCoeff
    REAL(KIND=dp) :: fluxes(:),areas(:)
    INTEGER :: fluxescomputed(:)
    
    INTEGER :: t
    TYPE(Element_t), TARGET :: SideElement
    TYPE(Element_t), POINTER :: Element, Parent    
    TYPE(ValueList_t), POINTER :: Material
    REAL(KIND=dp) :: LocalCoordinates(3),Point(3)
    REAL(KIND=dp) :: Basis(Model % MaxElementNodes),ParentBasis(Model % MaxElementNodes),&
        EnergyTensor(3,3,Model % MaxElementNodes),&
        EnergyCoeff(Model % MaxElementNodes) 
    REAL(KIND=dp) :: DetJ,U,V,W,up,vp,wp,S,A,L,C(3,3),x,y,z,dx,dy,dz,ds,dsmax,Grad(3)
    REAL(KIND=dp) :: func, coeff, Normal(3), Flow(3), x0, y0, z0, pos(2), flux
    REAL(KIND=DP), POINTER :: Pwrk(:,:,:) => Null()
    INTEGER, POINTER :: ParentIndexes(:), PermIndexes(:), SideIndexes(:), OnLine(:,:)

    LOGICAL :: Stat, Permutated, Inside    
    INTEGER :: i,j,k,p,q,DIM,bc,NoDofs,pn,Line, NoSides, Side, NodeNumber, LineNode(2), istat
    TYPE(GaussIntegrationPoints_t) :: IntegStuff
    TYPE(Nodes_t) :: ParentNodes, LineNodes, SideNodes


    n = Model % MaxElementNodes
    DIM = CoordinateSystemDimension()

    ALLOCATE(ParentNodes % x(n), ParentNodes % y(n), ParentNodes % z(n), PermIndexes(n), &
	SideNodes % x(n), SideNodes % y(n), SideNodes % z(n), SideIndexes(n), &
	LineNodes % x(n), LineNodes % y(n), LineNodes % z(n), &
	OnLine(Mesh % NumberOfNodes,2), STAT=istat )
    IF( istat /= 0 ) CALL Fatal('PolylineIntegrals','Memory allocation error') 


    CALL Info('PolylineIntegrals','OPERATOR: '//TRIM(OperName))

    areas = 0.0_dp
    coeff = 1.0_dp   

    NoDofs = Var % Dofs
    IF( OperName == 'area' ) THEN
      Permutated = .FALSE.
    ELSE
      Permutated = ASSOCIATED(Var % Perm)
    END IF

    DO i=1,3          
      EnergyTensor(i,i,:) = 1.0d0
    END DO

    SELECT CASE(OperName)
      
      CASE('diffusive flux') 
      IF(NoDofs /= 1) THEN
        CALL Warn(Caller,'diffusive flux & NoDofs /= 1?')
        RETURN
      END IF

      CASE ('convective flux')
      IF(NoDofs /= 1 .AND. NoDofs < DIM) THEN
        CALL Warn(Caller,'convective flux & NoDofs < DIM?')
        RETURN
      END IF
      
      CASE ('area','boundary int','boundary int mean')

    CASE DEFAULT 
      CALL Warn(Caller,'Unknown physical operator')

    END SELECT

    DIM = CoordinateSystemDimension()
    IF(DIM < 3) THEN
      LineNodes % z = 0.0d0
      SideNodes % z(1:2) = 0.0d0
    END IF

    fluxes = 0.0d0
    SideElement % TYPE => GetElementType( 202, .FALSE.)
    SideElement % Bdofs = 0
    Element => SideElement

!   /* Go through the line segments */
    DO Line = 1,NoLines 

      LineNodes % x(1:2) = LineCoordinates(2*Line-1:2*Line,1) 
      LineNodes % y(1:2) = LineCoordinates(2*Line-1:2*Line,2) 
      IF(DIM == 3) LineNodes % z(1:2) = LineCoordinates(2*Line-1:2*Line,3) 
      OnLine = 0

      DO t = 1, Mesh % NumberOfBulkElements 
        
        Parent => Mesh % Elements(t)
        Model % CurrentElement => Mesh % Elements(t)

        NoSides = Parent % TYPE % ElementCode / 100  
        IF(NoSides < 3 .OR. NoSides > 4) CYCLE
        
        pn = Parent % TYPE % NumberOfNodes
        ParentIndexes => Parent % NodeIndexes

        CALL CopyElementNodesFromMesh( ParentNodes, Mesh, pn, ParentIndexes)

        IF(Permutated) THEN
          PermIndexes(1:pn) = Var % Perm(ParentIndexes(1:pn))
          IF (ANY( PermIndexes(1:pn) == 0)) CYCLE
        ELSE
          PermIndexes(1:pn) = ParentIndexes(1:pn)        
        END IF

        NodeNumber = 0

        DO Side = 1, NoSides

          SideIndexes(1) = ParentIndexes(Side)
          SideIndexes(2) = ParentIndexes(MOD(Side,NoSides)+1)
          
          SideNodes % x(1:2) = Mesh % Nodes % x(SideIndexes(1:2))
          SideNodes % y(1:2) = Mesh % Nodes % y(SideIndexes(1:2))
          IF(DIM == 3) SideNodes % z(1:2) = Mesh % Nodes % z(SideIndexes(1:2))
  
          CALL LineIntersectionCoords(SideNodes,LineNodes,Inside,x0,y0,z0,u)

          IF(.NOT. Inside) CYCLE

          NodeNumber = NodeNumber + 1        
          ElementNodes % x(NodeNumber) = x0
          ElementNodes % y(NodeNumber) = y0
          ElementNodes % z(NodeNumber) = z0
          pos(NodeNumber) = u

        END DO

        IF(NodeNumber == 0) CYCLE

        IF(NodeNumber > 2) THEN
          CALL Warn('PolylineIntergrals','There should not be more than 2 intersections!')
          CYCLE
        END IF

        !---------------------------------------------------------------------------
        ! If there is only one intersection the other end of the node must lie
        ! inside the element. Assuming that the line is long compared to the 
        ! element the correct end of the line segment may be easily deduced.
        !---------------------------------------------------------------------------
        IF(NodeNumber == 1) THEN
          IF(pos(1) < 0.5d0) THEN
            i = 1
            pos(2) = 0.0
          ELSE
            i = 2
            pos(2) = 1.0
          END IF
          x0 = LineNodes % x(i)
          y0 = LineNodes % y(i)
          z0 = LineNodes % z(i)            

          ElementNodes % x(2) = LineNodes % x(i)
          ElementNodes % y(2) = LineNodes % y(i)
          ElementNodes % z(2) = LineNodes % z(i)            
        END IF

        IF(ABS(pos(1)-pos(2)) < 1.0d-8) CYCLE

        !-----------------------------------------------------------------------------
        ! Change the order of nodes so that the normal always points to the same direction          
        !-----------------------------------------------------------------------------
        IF(pos(1) < pos(2)) THEN
          ElementNodes % x(2) = ElementNodes % x(1)
          ElementNodes % y(2) = ElementNodes % y(1)
          ElementNodes % z(2) = ElementNodes % z(1)
          ElementNodes % x(1) = x0
          ElementNodes % y(1) = y0
          ElementNodes % z(1) = z0           
        END IF
        
        !--------------------------------------------------------------------------------
        ! The following avoids the cases where the line goes exactly at the element 
        ! interface and therefore the flux would be computed twice
        !--------------------------------------------------------------------------------
        dx = ElementNodes % x(1) - ElementNodes % x(2)
        dy = ElementNodes % y(1) - ElementNodes % y(2)
        dsmax = SQRT(dx*dx+dy*dy)
        LineNode = 0


        DO i=1,Parent % TYPE % ElementCode / 100 
          DO j=1,2
            dx = ParentNodes % x(i) - ElementNodes % x(j)
            dy = ParentNodes % y(i) - ElementNodes % y(j)
            ds = SQRT(dx*dx+dy*dy)
            IF(ds < 1.0d-4 * dsmax) LineNode(j) = ParentIndexes(i)
          END DO
        END DO

        IF(ALL(LineNode(1:2) > 0)) THEN
          IF(ANY(OnLine(LineNode(1),:) == LineNode(2))) CYCLE
          IF(ANY(OnLine(LineNode(2),:) == LineNode(1))) CYCLE

          IF(OnLine(LineNode(1),1) == 0) THEN
            OnLine(LineNode(1),1) = LineNode(2)
          ELSE IF(OnLine(LineNode(1),2) == 0) THEN
            OnLine(LineNode(1),2) = LineNode(2)
          ELSE
            CALL Warn('PolylineIntegrate','This should never happen')
          END IF
          
          IF(OnLine(LineNode(2),1) == 0) THEN
            OnLine(LineNode(2),1) = LineNode(1)
          ELSE IF(OnLine(LineNode(2),2) == 0) THEN
            OnLine(LineNode(2),2) = LineNode(1)
          ELSE
            CALL Warn('PolylineIntegrate','This should never happen')
          END IF
        END IF
        
        i = ListGetInteger( Model % Bodies(Parent % BodyId) % Values, 'Material', &
            minv=1, maxv=Model % NumberOFMaterials )
        Material => Model % Materials(i) % Values
        fluxescomputed(Line) = fluxescomputed(Line) + 1
        
        
        SELECT CASE(OperName)
          
          CASE('diffusive flux') 
          
            IF(GotCoeff) THEN
              CALL ListGetRealArray( Material, CoeffName, Pwrk, &
                  pn, ParentIndexes, gotIt )

              EnergyTensor = 0._dp
              IF(GotIt) THEN
                IF ( SIZE(Pwrk,1) == 1 ) THEN
                  DO i=1,3
                    EnergyTensor( i,i,1:pn ) = Pwrk( 1,1,1:pn )
                  END DO
                ELSE IF ( SIZE(Pwrk,2) == 1 ) THEN
                  DO i=1,MIN(3,SIZE(Pwrk,1))
                    EnergyTensor(i,i,1:pn) = Pwrk(i,1,1:pn)
                  END DO
                ELSE
                  DO i=1,MIN(3,SIZE(Pwrk,1))
                    DO j=1,MIN(3,SIZE(Pwrk,2))
                      EnergyTensor( i,j,1:pn ) = Pwrk(i,j,1:pn)
                    END DO
                  END DO
                END IF
              ELSE 
                DO i=1,3          
                  EnergyTensor(i,i,1:pn) = 1.0d0
                END DO
              END IF
            END IF

          CASE ('convective flux','area')
            EnergyCoeff(1:n) = ListGetReal( Material, CoeffName, pn, ParentIndexes, gotIt )
            IF(.NOT. GotIt) EnergyCoeff(1:pn) = 1.0d0
          
        END SELECT

!------------------------------------------------------------------------------
!    Numerical integration
!------------------------------------------------------------------------------

        IntegStuff = GaussPoints( Element, 2 )
        
        DO i=1,IntegStuff % n
          U = IntegStuff % u(i)
          V = IntegStuff % v(i)
          W = IntegStuff % w(i)

!------------------------------------------------------------------------------
!        Basis function values & derivatives at the integration point
!------------------------------------------------------------------------------
          stat = ElementInfo( Element,ElementNodes,U,V,W,DetJ,Basis)
          
!------------------------------------------------------------------------------
!      Coordinatesystem dependent info
!------------------------------------------------------------------------------
          s = DetJ * IntegStuff % s(i)            
          IF(CurrentCoordinateSystem() /= Cartesian ) THEN
            x = SUM( ElementNodes % x(1:n) * Basis(1:n) )
            s = 2._dp * PI * x * s
          END IF
            
          Normal = NormalVector( Element,ElementNodes,U,V,.FALSE. )

!------------------------------------------------------------------------------
!      Because the intersection nodes do not really exist the field variables 
!      must be evaluated using the nodes of the parent element.
!------------------------------------------------------------------------------
            
          Point(1) = x
          Point(2) = y
          Point(3) = z
          
          IF( .NOT. PointInElement( Parent, ParentNodes, Point, LocalCoordinates ) ) THEN
            CALL Warn('PolylineIntegrals','The node should be in the element by construction!')
            CYCLE
          END IF
          
          Up = LocalCoordinates(1)
          Vp = LocalCoordinates(2)
          Wp = LocalCoordinates(3)
                    
          stat = ElementInfo( Parent,ParentNodes,Up,Vp,Wp,DetJ,ParentBasis)
  
          SELECT CASE(OperName)
            
            CASE ('diffusive flux')           

              Grad = 0.0_dp
              k = 0
              ! This can use p-element basis etc. when evaluating the gradient in the parent element.
              CALL EvaluateVariableAtGivenPoint(k,Grad,Mesh,Var,Element=Parent,LocalCoord=LocalCoordinates,&
                  LocalBasis=ParentBasis,DoGrad=.TRUE.)
              
              Flow = 0.0d0
              DO j = 1, DIM
                DO k = 1, DIM
                  Flow(j) = Flow(j) + SUM( EnergyTensor(j,k,1:n) * Basis(1:n) ) * Grad(k)
                END DO
              END DO              
              fluxes(Line) = fluxes(Line) + s * SUM(Normal(1:DIM) * Flow(1:DIM))
                        
            CASE ('convective flux')            
          
              coeff = SUM( EnergyCoeff(1:pn) * ParentBasis(1:pn))              
              IF(NoDofs == 1) THEN
                func = SUM( Var % Values(PermIndexes(1:pn)) * ParentBasis(1:pn) )                
                fluxes(Line) = fluxes(Line) + s * coeff * func
              ELSE 
                DO j=1,DIM
                  Flow(j) = coeff * &
                      SUM( Var % Values(NoDofs*(PermIndexes(1:pn)-1)+j) * ParentBasis(1:pn) )
                END DO
                fluxes(Line) = fluxes(Line) + s * coeff * SUM(Normal * Flow)
              END IF
 
           CASE ('boundary int','boundary int mean')            
              coeff = SUM( EnergyCoeff(1:pn) * ParentBasis(1:pn))              
              IF(NoDofs == 1) THEN
                func = SUM( Var % Values(PermIndexes(1:pn)) * ParentBasis(1:pn) )                
                fluxes(Line) = fluxes(Line) + s * coeff * func
              ELSE 
                flux = 0.0_dp
                DO j=1,NoDofs
                  flux = flux + &
                      SUM( Var % Values(NoDofs*(PermIndexes(1:pn)-1)+j) * ParentBasis(1:pn) )**2
                END DO
                flux = coeff * SQRT(flux) 
                fluxes(Line) = fluxes(Line) + s * flux
              END IF
                         
            CASE ('area')                        
              coeff = SUM( EnergyCoeff(1:pn) * ParentBasis(1:pn))
              fluxes(Line) = fluxes(Line) + s * coeff 
              
            END SELECT

            areas(Line) = areas(Line) + s * coeff

        END DO

      END DO

    END DO

    DEALLOCATE( ParentNodes % x, ParentNodes % y, ParentNodes % z, PermIndexes, &
        SideNodes % x, SideNodes % y, SideNodes % z, SideIndexes, &
        LineNodes % x, LineNodes % y, LineNodes % z, OnLine)

  END SUBROUTINE PolylineIntegrals
!------------------------------------------------------------------------------


!------------------------------------------------------------------------------
!> This subroutine tests whether the line segment goes through the current
!> face of the element. 
!------------------------------------------------------------------------------
  SUBROUTINE LineIntersectionCoords(Plane,Line,Inside,x0,y0,z0,frac)

    TYPE(Nodes_t) :: Plane, Line
    LOGICAL :: Inside
    REAL (KIND=dp) :: x0, y0, z0, frac

    REAL (KIND=dp) :: A(3,3),B(3),C(3),eps=1.0d-6,detA,absA

    Inside = .FALSE.
    
    ! In 2D the intersection is between two lines
    A(1,1) = Line % x(2) - Line % x(1)
    A(2,1) = Line % y(2) - Line % y(1)
    A(1,2) = Plane % x(1) - Plane % x(2)
    A(2,2) = Plane % y(1) - Plane % y(2)

    detA = A(1,1)*A(2,2)-A(1,2)*A(2,1)
    absA = SUM(ABS(A(1,1:2))) * SUM(ABS(A(2,1:2)))
    
    IF(ABS(detA) <= eps * absA + 1.0d-20) RETURN
    
    B(1) = Plane % x(1) - Line % x(1) 
    B(2) = Plane % y(1) - Line % y(1) 
    
    CALL InvertMatrix( A,2 )
    C(1:2) = MATMUL(A(1:2,1:2),B(1:2))
    
    IF(ANY(C(1:2) < 0.0) .OR. ANY(C(1:2) > 1.0d0)) RETURN
    
    Inside = .TRUE.
    frac = C(1)
    X0 = Line % x(1) + C(1) * (Line % x(2) - Line % x(1))
    Y0 = Line % y(1) + C(1) * (Line % y(2) - Line % y(1))
    Z0 = Line % z(1) + C(1) * (Line % z(2) - Line % z(1))
    
  END SUBROUTINE LineIntersectionCoords


  SUBROUTINE SetPseudoNorm()
    LOGICAL :: GotNorm
    REAL(KIND=dp) :: c
    REAL(KIND=dp), POINTER :: RefVals(:,:)

    GotNorm = .FALSE.
    
    NormInd = 0
    Name = ListGetString( Params,'Show Norm Name',GotIt)
    IF(GotIt) THEN
      DO No=1,NoValues
        IF( TRIM(Name) == TRIM(ValueNames(No)) ) THEN
          NormInd = No
          CALL Info(Caller,'Associating scalar '//I2S(No)//' to norm!',Level=8)
          EXIT
        END IF
      END DO
      IF(NormInd == 0 ) THEN
        CALL Warn(Caller,'Could not find scalar for norm: '//TRIM(Name))
      END IF
    ELSE
      NormInd = ListGetInteger( Params,'Show Norm Index',GotIt)
    END IF

    IF( NormInd > 0 .AND. NormInd <= NoValues ) THEN
      Norm = Values( NormInd )
      GotNorm = .TRUE.

      Name = ListGetString( Params, 'Equation',GotIt)
      IF(.NOT. GotIt) Name = Caller

    ELSE IF( ListCheckPresent( Params,'Reference Values') ) THEN
      RefVals => ListGetConstRealArray( Params,'Reference Values') 
      n = SIZE( RefVals, 1 )
      IF( n > NoValues ) THEN
        CALL Fatal('SaveScalars','Size of "Reference Values" bigger than computed values!')
      END IF      
      Norm = 0.0_dp
      DO i=1,n
        IF( ABS(RefVals(i,1)) > EPSILON(c) ) THEN
          c = Values(i) / RefVals(i,1)
          c = MAX( c, 1.0_dp /c ) 
        ELSE
          c = 1.0_dp + ABS(Values(i))
        END IF
        Norm = Norm + c
      END DO
      Norm = Norm / n
      GotNorm = .TRUE.
      
      ! By construction the reference norm is 1.
      CALL ListAddNewConstReal( Params,'Reference Norm',1.0_dp)

      IF( InfoActive(20 ) ) THEN
        PRINT *,'Reference  for SaveScalars:'
        PRINT *,'ThisResults:',Values(1:n)
        PRINT *,'RefResults:',RefVals(1:n,1)
      END IF
    END IF

    IF( GotNorm ) THEN
      Solver % Variable % Values = Norm 
      Solver % Variable % Norm = ABS( Norm ) 
      
      ! Here the name is ComputeChange in order to get the change also to ElmerGUI
      ! albeit in a very dirty style. One could also edit ElmerGUI....
      WRITE( Message, '(a,g15.8,g15.8,a)') &
          'SS (ITER=1) (NRM,RELC): (',Norm, Change,&
          ' ) :: '//TRIM( Name )
      CALL Info( 'ComputeChange', Message, Level=3 )
    END IF
    
  END SUBROUTINE SetPseudoNorm
    
  
!------------------------------------------------------------------------------
END SUBROUTINE SaveScalars
!------------------------------------------------------------------------------

!> \}
