! copyright 2005,  J. E. Akin, all rights reserved.
! file: scp_lib.f  Super Convergent Patch Library

SUBROUTINE TURN_OFF_SCP
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
!  FORCE NO SUPER CONVERGENT PATCHES, DEFAULTS TO ON IF N_QP > 0
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
 Use SCP_Type_Data
! U_SCPR     = UNIT FOR SUPER_CONVERGENT PATCH RECOVERY, IF > 0
   U_SCPR = 0   ! SET TO INACTIVE
   PRINT *, 'NOTE: SUPER_CONVERGENT PATCH RECOVERY TURNED OFF'
END SUBROUTINE TURN_OFF_SCP

SUBROUTINE TURN_ON_SCP
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
!  FORCE USE OF SUPER CONVERGENT PATCHES
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
 Use SCP_Type_Data
! U_SCPR     = UNIT FOR SUPER_CONVERGENT PATCH RECOVERY, IF > 0
   U_SCPR = 29   ! SET TO ACTIVE
   PRINT *, 'NOTE: SUPER_CONVERGENT PATCH RECOVERY TURNED ON'
END SUBROUTINE TURN_ON_SCP

SUBROUTINE INQUIRE_ABOUT_SCP_UNIT ()
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
!  DEBUG AID FOR SCP RECOVERY DIRECT ACCESS UNIT, U_SCPR
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
  Use SCP_Type_Data
  IMPLICIT NONE
    CHARACTER(LEN=11) :: ACCESS, ACTION, DIRECT, FORM, POS
    LOGICAL           :: EXIST,  OPEN   
    INTEGER           :: IOSTAT, NEXTREC, RECL

! U_SCPR     = UNIT FOR SUPER_CONVERGENT PATCH RECOVERY, IF > 0
      
      INQUIRE (U_SCPR, ACCESS = ACCESS, ACTION = ACTION,          &
               DIRECT = DIRECT, EXIST = EXIST, FORM = FORM,       &
               IOSTAT = IOSTAT, NEXTREC = NEXTREC, OPENED = OPEN, &
               POSITION = POS,  RECL = RECL )
      PRINT *,''
      PRINT *,'** SCP RANDOM ACCESS UNIT STATUS **'
      PRINT *,'UNIT     = ', U_SCPR
      PRINT *,'ACCESS   = ', ACCESS
      PRINT *,'ACTION   = ', ACTION
      PRINT *,'DIRECT   = ', DIRECT
      PRINT *,'EXIST    = ', EXIST
      PRINT *,'FORM     = ', FORM
      PRINT *,'IOSTAT   = ', IOSTAT
      PRINT *,'NEXTREC  = ', NEXTREC
      PRINT *,'OPENED   = ', OPEN
      PRINT *,'POSITION = ', POS
      PRINT *,'RECL     = ', RECL
      PRINT *,''
END SUBROUTINE INQUIRE_ABOUT_SCP_UNIT 

SUBROUTINE CALC_SCP_AVE_NODE_FLUXES  (NODES, X, L_NEIGH, &
                                      SCP_AVERAGES)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
!  CALCULATE THE SUPER_CONVERGENCE_PATCH AVERAGE FLUXES   
!              AT ALL NODES IN THE MESH
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
! (USES A LEAST SQUARES FIT THROUGH THE GAUSS POINT VALUES, AND
! SOLVES THAT LOCAL PROBLEM BY THE SINGULAR VALUE DECOMPOSITION 
! SOLUTION METHOD)
  Use System_Constants ! for MAX_NP, NEIGH_P, N_ELEMS, N_PATCH, N_QP, 
                       ! SCP_FIT, U_SCPR, ON_BOUNDARY
  Use Sys_Properties_Data  ! for DELETED
  Use Elem_Type_Data       ! for LAST_LT, LT_*
  Use SCP_Type_Data        ! for SCP_H (SCP_N), SCP_DLH (SCP_PARM,SCP_N)
  Use Interface_Header
   IMPLICIT NONE
   INTEGER,  INTENT (IN)  :: NODES        (L_S_TOT, NOD_PER_EL)
   INTEGER,  INTENT (IN)  :: L_NEIGH      (NEIGH_P, N_PATCH)  
   REAL(DP), INTENT (IN)  :: X            (MAX_NP, N_SPACE)
   REAL(DP), INTENT (OUT) :: SCP_AVERAGES (MAX_NP, SCP_FIT)

   INTEGER :: MEMBERS     (NEIGH_P + 1)   ! ELEMENTS IN PATCH  
   INTEGER :: SCP_COUNTS  (MAX_NP)        ! PATCH HITS PER NODE

   REAL(DP), ALLOCATABLE :: PATCH_SQ  (:, :)  ! SCRATCH ARRAY
   REAL(DP), ALLOCATABLE :: PATCH_DAT (:, :)  ! SCRATCH ARRAY
   REAL(DP), ALLOCATABLE :: PATCH_P   (:, :)  ! SCRATCH ARRAY
   REAL(DP), ALLOCATABLE :: PATCH_WRK (:)     ! SCRATCH ARRAY

   REAL(DP), ALLOCATABLE :: PATCH_FIT (:, :)         ! ANSWERS

   REAL(DP) :: XYZ_MIN (N_SPACE), XYZ_MAX (N_SPACE)  ! BOUNDS
   REAL(DP) :: XYZ     (N_SPACE), FLUX    (N_R_B)    ! POINT & FLUX 
   REAL(DP) :: POINT   (N_SPACE)                     ! SCP POINT

   INTEGER  :: POINTS                ! SIZES OF CURRENT PATCH MATRICES
   INTEGER  :: L_IN_PATCH            ! NUMBER OF ELEMS IN CURRENT PATCH
   INTEGER  :: FIT, IL, IP, IQ       ! LOOPS
   INTEGER  :: ROW                   ! IN PATCH ARRAYS
   INTEGER  :: LM, LT, REC_LM_IQ, SCP_STAT
   INTEGER, PARAMETER :: ONE = 1

!                   VARIABLES:
! DELETED    = TRUE IF THE ELEMENT HAS BEEN DELETED
! E          = CONSTITUTIVE MATRIX
! ELEM_NODES = THE NOD_PER_EL INCIDENCES OF THE ELEMENT
! INDEX      = SYSTEM DOF NUMBERS ASSOCIATED WITH ELEMENT
! L_HOMO     = 1, IF ELEMENT PROPERTIES ARE HOMOGENEOUS
! L_SHAPE    = SHAPE FLAG FOR QUADRATURE RULE SELECTION
! L_S_TOT    = TOTAL NUMBER OF ELEMENTS & THEIR SEGMENTS
! L_NEIGH    = ELEMENTS FORMING THE PATCH    !b
! L_TYPE     = ELEMENT TYPE NUMBER
! LT         = ELEMENT TYPE NUMBER (IF USED)
! LT_DATA    = ARRAY OF CONSTANTS FOR EACH ELEMENT TYPE
! LT_QP      = NUMBER OF QUADRATURE PTS FOR ELEMENT TYPE
! LT_USER    = OPTIONAL USER CODE ASSOCIATED WITH TYPE
! MAX_NP     = NUMBER OF SYSTEM NODES
! MEMBERS    = ELEMENT NUMBERS MACKING UP A SCP 
! NEIGH_P    = MAXIMUM NUMBER OF ELEMENTS IN A PATCH
! NOD_PER_EL = NUMBER OF NODES PER ELEMENT
! NODES      = NODAL INCIDENCES OF ALL ELEMENTS
! N_ELEMS    = NUMBER OF ELEMENTS IN SYSTEM
! N_PATCH    = NUMBER OF PATCHES, MAX_NP OR N_ELEMS
! N_QP       = MAXIMUN NUMBER OF QUADRATURE POINTS, >= LT_QP
! N_R_B      = NUMBER OF ROWS IN B AND E MATRICES
! N_SPACE    = DIMENSION OF SPACE
! PATCH_FIT  = LOCAL PATCH VALUES FOR FLUX AT ITS NODES
! POINT      = LOCAL POINT IN PATCH INTERPOLATION SPACE
! POINTS     = TOTAL NUMBER OF QUADRATURE POINTS IN CURRENT PATCH
! SCP_AVERAGES = AVERAGED FLUXES AT ALL NODES IN MESH
! SCP_FIT    = NUMBER IF TERMS BEING FIT, OR AVERAGED: = N_R_B MIN    
! SCP_INC    = INCREASE PATCH POLYNOMIAL DEGREE BY THIS (1 or 2)
! SCP_RECORD_NUMBER = SCP DIRECT ACCESS RECORD LOCATOR: (N_ELEMS, N_QP)
! U_SCPR     = UNIT FOR SUPER_CONVERGENT PATCH RECOVERY, IF > 0
! X          = COORDINATES OF SYSTEM NODES
! XYZ        = SPACE COORDINATES AT A POINT
! XYZ_MAX    = UPPER BOUNDS FOR SCP GEOMETRY
! XYZ_MIN    = LOWER BOUNDS FOR SCP GEOMETRY

    !b IF ( I_BUG > 0 ) CALL INQUIRE_ABOUT_SCP_UNIT 
    LAST_LT = 0 ; SCP_COUNTS = 0 ; SCP_AVERAGES = 0.d0  ! INITIALIZE
    IF ( U_SCPR  == 0 ) RETURN                          ! NO DATA 
    IF ( N_QP    == 0 ) RETURN                          ! NO DATA   

! Note: Replace next line with logical function SCP_DATA_STORED to
! allow for other random access data structures
    IF ( MAXVAL ( SCP_RECORD_NUMBER ) < 1 ) RETURN      ! NO DATA
! End Note

    IF ( N_ELEMS == 1 .AND. N_QP < NOD_PER_EL ) THEN
      PRINT *, 'ERROR, a single elem patch requires N_Q_P > NOD_PER_EL'
      PRINT *, 'Single element patch error, CALC_SCP_AVE_NODE_FLUXES'
      RETURN ! assume expert user
     !b STOP 'Single element patch error, CALC_SCP_AVE_NODE_FLUXES'
    END IF ! SINGLE ELEMENT MESH

    LT = 1                           ! INITIALIZE
    IF ( N_PATCH == 0 ) THEN ! No patch type selected
      PRINT *,'WARNING: CALC_SCP_AVE_NODE_FLUXES, NO PATCH TYPE GIVEN'
      PRINT *,'PLACE KEYWORD OPTION scp_neigh_el OR scp_neigh_face OR'
      PRINT *,'scp_neigh_pt IN THE keyword CONTROLS.'
      PRINT *,'Skipping SCP averaging and error estimates.'
      N_WARN = N_WARN + 1 ; RETURN
    END IF
    IF ( SCP_INC > 0 ) PRINT *, 'NOTE: PATCH DEGREE INCREASED BY ', &
                                 SCP_INC

    DO IP = 1, N_PATCH               ! LOOP OVER EACH PATCH  
      MEMBERS  = 0 ; POINTS = 0                   ! INITIALIZE

!       ELEMENT OR NODAL CENTERED PATCH TYPE
      IF ( .NOT. SCP_NEIGH_PT ) THEN ! ELEMENT BASED ELEMENTS
        IF ( IP == 1 ) PRINT *,'NOTE: ELEMENT BASED PATCHES '
        IF ( DELETED (IP) ) CYCLE      ! TO AN ACTIVE ELEMENT

!        GET ELEMENT NEIGHBORS TO DEFINE THE PATCH
        MEMBERS = (/ IP, L_NEIGH (:, IP) /)      ! TO SIMPLIFY LOOP 
      ELSE ! NODAL BASED PATCH OF ELEMENTS
        IF ( IP == 1 ) PRINT *,'NOTE: NODAL BASED PATCHES'

!        GET ELEMENT NEIGHBORS TO DEFINE THE PATCH
        !b MEMBERS    = L_NEIGH (:, IP) 
        MEMBERS (1:NEIGH_P) = L_NEIGH (:, IP) ! bug fix
      END IF ! PATCH BASIS
      L_IN_PATCH = COUNT ( MEMBERS > 0 )
      IF ( L_IN_PATCH <= 1 ) CYCLE ! TO AN ACTIVE PATCH

!      FIND TYPE OF SCP FOR THIS GROUP OF ELEMENTS, VERIFY GEOMETRY
      IF ( DEBUG_SCP ) THEN
        PRINT *, 'IP, L_IN_PATCH ', IP, L_IN_PATCH
        PRINT *, 'SIZE(MEMBERS) ', SIZE(MEMBERS)
        PRINT *, 'MEMBERS ', MEMBERS
      END IF
      !b CALL DETERMINE_SCP_BOUNDS (L_IN_PATCH, MEMBERS, NODES, X,    &
                                 !b XYZ_MIN, XYZ_MAX, POINTS)
      CALL DETERMINE_SCP_BOUNDS (L_IN_PATCH, MEMBERS(1:L_IN_PATCH),   &
                                 NODES, X, XYZ_MIN, XYZ_MAX, POINTS)
!    Note, DETERMINE_SCP_BOUNDS now includes 
!    CALL SET_SCP_TYPE_DATA  (LT) ! CONTROLS FOR SCP TYPE
!    which could be broken out and placed here

      IF ( PATCH_ALLOC_STATUS ) THEN
!         DEALLOCATE LOCAL PATCH RELATED ARRAYS
        DEALLOCATE ( PATCH_WRK )
        DEALLOCATE ( PATCH_SQ  )
        DEALLOCATE ( PATCH_FIT )
        DEALLOCATE ( PATCH_DAT )
        DEALLOCATE ( PATCH_P   )
        PATCH_ALLOC_STATUS = .FALSE.
      END IF ! STATUS CHECK

!       ALLOCATE LOCAL PATCH RELATED ARRAYS
      ALLOCATE ( PATCH_P   (POINTS, SCP_N  ) )
      ALLOCATE ( PATCH_DAT (POINTS, SCP_FIT) )
      ALLOCATE ( PATCH_FIT (SCP_N , N_R_B  ) ) 
      ALLOCATE ( PATCH_SQ  (SCP_N , SCP_N  ) )
      ALLOCATE ( PATCH_WRK (SCP_N )          )
      PATCH_ALLOC_STATUS = .TRUE. 

!      ZERO PATCH WORKSPACE AND RESULTS ARRYS
      PATCH_P   = 0.d0 ; PATCH_DAT = 0.d0 ; PATCH_FIT = 0.d0     
      PATCH_SQ  = 0.d0 ; PATCH_WRK = 0.d0             

      IF ( DEBUG_EL_TYPE .OR. DEBUG_SCP ) THEN
        PRINT *,'DEBUG_ Elem type data ' 
        CALL IPRINT (LT_DATA, MAX_DAT,  N_L_TYPE) 
        DO IL = 1, N_L_TYPE 
          CALL GET_ELEM_TYPE_DATA  (IL) 
          CALL LIST_ELEM_TYPE_DATA (IL) 
        END DO 
      END IF ! Debug

!      PREPARE LEAST SQUARES FIT MATRICES
      LAST_LT = 0                           ! INITIALIZE, AGAIN
      ROW = 0                               ! INITIALIZE
      DO IL = 1, L_IN_PATCH                 ! PATCH ELEMENT MEMBER LOOP
        LM = MEMBERS (IL)                   ! ELEMENT IN PATCH
! XXX died in pt based patch , Subscript 1 of DELETED (value 0) !b turn off
        !b IF ( DELETED (LM)     ) CYCLE       ! TO ACTIVE ELEMENT ! XXX

!b!         OMIT FROM FLUX AVERAGES IF A BOUNDARY NODE
!bIF ( IL == 1 .AND. ON_BOUNDARY (LM) ) THEN 
!b  PRINT *,'NOTE: SKIPPING LM ON BOUNDARY PATCH ', IL 
!b  CYCLE       ! OMIT IN FIT 
!bEND IF 

        IF ( N_L_TYPE > 1) LT = L_TYPE (LM) ! GET ELEMENT TYPE NUMBER
        IF ( LT /= LAST_LT ) THEN           ! NEW ELEMENT TYPE  
          CALL GET_ELEM_TYPE_DATA (LT)      ! GET CONTROLS FOR ELEM TYPE
          IF ( TYPE_APLY_ALLOC ) CALL DEALLOCATE_TYPE_APPLICATION
          CALL ALLOCATE_TYPE_APPLICATION     ! for ELEM_NODES (LT_N)
          LAST_LT = LT
        END IF ! a new element type and scp type

        DO IQ = 1, LT_QP                    ! LOOP OVER GAUSS POINTS
          ROW = ROW + 1                           ! UPDATE LOCATION

! Note: The following lines should be replaced with a generic
! function call so other data structures, like arrays, could be
! used for the random access recovery:
! RANDOM_GET_EL_QP_FLUX_DATA (LM, IQ, XYZ, FLUX)
! INTEGER,  INTENT (IN)  :: LM, IQ
! REAL(DP), INTENT (OUT) :: XYZ (N_SPACE), FLUX (N_R_B) 

!           RECOVER EACH FLUX VECTOR TO BE USED IN LEAST SQ FIX
          REC_LM_IQ = SCP_RECORD_NUMBER (LM, IQ)  ! GET RECORD NUMBER
          IF ( REC_LM_IQ < 0 ) THEN               ! INVALID RECORD
            N_WARN = N_WARN + 1
            PRINT *,'WARNING, IS KEYWORD gauss < EL_TYPE RULE ?'
            PRINT *, 'ERROR: INVALID SCP_RECORD_NUMBER IN PATCH ', IP
            PRINT *, 'MEMBER ELEMENT = ', LM, ' AT POINT ', IQ
            STOP 'BAD SCP_RECORD_NUMBER, CALC_SCP_AVE_NODE_FLUXES   '
          END IF ! BAD RECORD

!          GET GAUSS PT COORD & FLUX SAVED IN LIST_ELEM_FLUXES

          READ (U_SCPR, REC = REC_LM_IQ, IOSTAT = SCP_STAT) XYZ, FLUX
          SELECT CASE ( SCP_STAT ) 
            CASE (:-1) 
              PRINT *,'EOR OR EOF, CALC_SCP_AVE_NODE_FLUXES   ' 
              STOP    'EOR OR EOF, CALC_SCP_AVE_NODE_FLUXES   ' 
            CASE (1:)
              PRINT *,'READ ERROR FOR SUPER_CONVERGENT PATCH'
              PRINT *,'MEMBER ELEMENT = ', LM, ' AT POINT ', IQ
              PRINT *,'REC_LM_IQ      = ', REC_LM_IQ
              PRINT *,'BAD SCP_STAT, CALC_SCP_AVE_NODE_FLUXES   '
              STOP    'BAD SCP_STAT, CALC_SCP_AVE_NODE_FLUXES   '
            CASE DEFAULT ! NO READ ERROR
          END SELECT ! RANDOM ACCESS READ ERROR
! End Note for RANDOM_GET_EL_QP_FLUX_DATA (LM, IQ, XYZ, FLUX)

!            CONVERT IQ XYZ TO LOCAL PATCH POINT
          POINT = GET_SCP_PT_AT_XYZ (XYZ, XYZ_MIN, XYZ_MAX)

!             EVALUATE PATCH INTERPOLATION AT LOCAL POINT
          IF ( .NOT. SCP_SCAL_ALLOC ) CALL ALLOCATE_SCP_INTERPOLATIONS
          CALL GEN_ELEM_SHAPE (POINT, SCP_H, SCP_N, N_SPACE, ONE)

!             INSERT FLUX & INTERPOLATIONS INTO PATCH MATRICES
          PATCH_DAT (ROW, 1:N_R_B) = FLUX  (:) 
          PATCH_P   (ROW, :)       = SCP_H (:)
        END DO ! FOR EACH IQ FLUX VECTOR

      END DO ! FOR EACH IL PATCH MEMBER
!   ASSEMBLY OF PATCH COMPLETED

!       VALIDATE CURRENT PATCH
      IF ( POINTS < SCP_N ) THEN 
        PRINT *,'WARNING: CALC_SCP_AVE_NODE_FLUXES SKIPPING PATCH ', &
                 IP, ' WITH ONLY ', POINTS, ' EQUATIONS'
        N_WARN = N_WARN + 1 ! INCREMENT WARNING 
        CYCLE ! TO NEXT PATCH 
      END IF ! INSUFFICIENT DATA 

!        USE SINGULAR VALUE DECOMPOSITION SOLUTION METHOD
      CALL SVDC_FACTOR (PATCH_P, POINTS, SCP_N, PATCH_WRK, PATCH_SQ)

      WHERE ( PATCH_WRK < EPSILON(1.d0) ) PATCH_WRK = 0.D0

      DO FIT = 1, N_R_B   ! LOOP FOR EACH FLUX COMPONENT TO FIT 
        CALL SVDC_BACK_SUBST (PATCH_P, PATCH_WRK, PATCH_SQ,        &
                              POINTS, SCP_N, PATCH_DAT (:, FIT),   &
                              PATCH_FIT (:, FIT))
      END DO ! FOR FLUXES   COMPONENTS

!       INTERPOLATE AVERAGES TO ALL NODES IN THE PATCH,
!       SCATTER PATCH NODAL AVERAGES TO SYSTEM NODES, INCREMENT COUNTS

      IF ( DEBUG_SCP .AND. IP == 1) PRINT *,'AVERAGING ELEMENT FLUXES'
      CALL EVAL_SCP_FIT_AT_PATCH_NODES (IP, NODES, X, L_IN_PATCH,  &
                             MEMBERS, XYZ_MIN, XYZ_MAX, PATCH_FIT, &
                             SCP_AVERAGES, SCP_COUNTS)

! NOTE: Copyright June 25, 1999, J. E. Akin
!  At this point one can use the LS smoothed gradients in this patch to
!  get a locally improved solution value estimate at all of the patch's
!  interior nodes. Algorithm is:
! a) read or reform element matrices, K_e & C_e, here 
! b) use a higher order quadrature rule to form the equilibrating vector
!    V_e = integral (B * (Stress_scp - Stress_e))
! c) assemble the elements in the patch into a local linear system:
!    K * U_new = C - V
! d) apply the previous solution as essential boundary conditions at all
!    nodes on the patch boundary.
! e) solve for the new interior node values (always a small system but
!    especially small for a patch of elements around a single node as in
!    original ZZ patch paper). call them u*_e.
! f) compute norm of u_e - u*_e to use in the final error estimator.
! [This is a form of the Loubignac iterative process.]

!        DEALLOCATE LOCAL PATCH RELATED ARRAYS
      IF ( SCP_SCAL_ALLOC ) CALL DEALLOCATE_SCP_INTERPOLATIONS
      DEALLOCATE ( PATCH_WRK )
      DEALLOCATE ( PATCH_SQ  )
      DEALLOCATE ( PATCH_FIT )
      DEALLOCATE ( PATCH_DAT )
      DEALLOCATE ( PATCH_P   )
      PATCH_ALLOC_STATUS = .FALSE.

    END DO ! FOR EACH IP PATCH IN MESH

!     FINALLY, AVERAGE FLUXES    FOR EACH NODAL HIT COUNT
    DO FIT = 1, MAX_NP ! FOR ALL NODES
      IF ( SCP_COUNTS (FIT) /= 0 ) THEN ! ACTIVE NODE
        SCP_AVERAGES (FIT, :) = SCP_AVERAGES (FIT, :) / SCP_COUNTS (FIT)
      ELSE  ! could skip since initialized
        SCP_AVERAGES (FIT, :) = 0.D0
      END IF   
    END DO ! FOR (AN UNWEIGHTED) AVERAGE 
!b xxx set logic flag assuring average has been done here

!      REPORT THE AVERAGED MAX & MIN VALUES AND LOCATIONS 
    IF ( F95 ) THEN ! Enhanced MAXLOC avilable
      !b CALL MAX_AND_MIN_SCP_AVE_F95 (SCP_AVERAGES) 
      CALL MAX_AND_MIN_SCP_AVE_F90 (SCP_AVERAGES) !b
    ELSE
      CALL MAX_AND_MIN_SCP_AVE_F90 (SCP_AVERAGES)
    END IF ! Enhanced MAXLOC avilable

END SUBROUTINE CALC_SCP_AVE_NODE_FLUXES    

SUBROUTINE EVAL_PT_FLUX_IN_SCP_PATCH (NP_SYS, X, XYZ_MIN, XYZ_MAX, &
                       PATCH_FIT, DIAGONAL_JACOBIAN, SCP_AVERAGES, &
                       SCP_COUNTS, NODE_DONE)          
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! CALCULATE THE SUPER_CONVERGENCE_PATCH FLUXES AT A NODE INSIDE
!    THE PATCH; OPTIONALLY EVALUATE 2ND OR 3RD DERIVATIVES    
!        (AVERAGE VALUES FROM ALL PATCHES LATER)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants ! for MAX_NP, NEIGH_L, N_ELEMS, N_QP,
                     ! SCP_FIT, U_SCPR,
Use Elem_Type_Data
Use SCP_Type_Data       ! for SCP_H (SCP_N), SCP_DLH (SCP_PARM, SCP_N)
Use Interface_Header
 IMPLICIT NONE
 INTEGER,  INTENT (IN)    :: NP_SYS       ! CURRENT GLOBAL NODE 
 REAL(DP), INTENT (IN)    :: X            (MAX_NP, N_SPACE)
 REAL(DP), INTENT (IN)    :: XYZ_MIN      (N_SPACE), XYZ_MAX (N_SPACE)
 REAL(DP), INTENT (IN)    :: PATCH_FIT    (SCP_N , N_R_B  )    !2b
 REAL(DP), INTENT (IN)    :: DIAGONAL_JACOBIAN (N_SPACE) 
 REAL(DP), INTENT (INOUT) :: SCP_AVERAGES (MAX_NP, SCP_FIT)
 INTEGER,  INTENT (INOUT) :: SCP_COUNTS   (MAX_NP)      
 LOGICAL,  INTENT (INOUT) :: NODE_DONE    (MAX_NP) 

 REAL(DP) :: XYZ     (N_SPACE), POINT   (N_SPACE)  ! SCP POINT
 REAL(DP) :: FLUX_PT (N_R_B)                       ! PT FLUXEX
 REAL(DP) :: SCP_DGH (N_SPACE, SCP_N)    ! GLOBAL FIRST DERIVATIVES
 REAL(DP) :: DERIV_1_FLUX (N_R_B * N_SPACE)   ! 2ND DERIVATIVES 

 INTEGER,      SAVE :: HIGHER, STORE (1) 
 INTEGER            :: IS              ! LOOPS
 INTEGER, PARAMETER :: ONE = 1

! DIAGONAL_JACOBIAN = CONSTANT GEOMETRY JACOBIAN
! FLUX_PT    = FLUX AT THE NODE INTERPOLATED FROM A SCP 
! HIGHER     = NUMBER OF HIGHER DEGREE DERIVATIVE FLUX TERMS
! INDEX      = SYSTEM DOF NUMBERS ASSOCIATED WITH ELEMENT
! L_SHAPE    = SHAPE FLAG FOR QUADRATURE RULE SELECTION
! L_TYPE     = ELEMENT TYPE NUMBER
! LT         = ELEMENT TYPE NUMBER (IF USED)
! LT_DATA    = ARRAY OF CONSTANTS FOR EACH ELEMENT TYPE
! LT_QP      = NUMBER OF QUADRATURE PTS FOR ELEMENT TYPE
! LT_USER    = OPTIONAL USER CODE ASSOCIATED WITH TYPE
! MAX_NP     = NUMBER OF SYSTEM NODES
! NODE_DONE  = TRUE, IF IT HAS ONE CONTRIBUTION FROM PATCH
! NP_SYS     = CURRENT GLOBAL NODE
! N_QP       = MAXIMUM NUMBER OF QUADRATURE POINTS, >= LT_QP
! N_SPACE    = DIMENSION OF SPACE
! PATCH_FIT  = LOCAL PATCH VALUES FOR FLUX AT ITS NODES
! POINT      = LOCAL POINT IN PATCH INTERPOLATION SPACE
! SCP_AVERAGES = AVERAGED FLUXES AT ALL NODES IN MESH
! SCP_CENTER_ONLY USES "CENTRAL" ELEMENT OR POINT VALUES
! SCP_FIT    = NUMBER IF TERMS BEING AVERAGED, = N_R_B MIN    
! SCP_ONLY_ONCE COMPUTES A NODE FLUX ONLY ONCE PER PATCH
! STORE      = ARRAY TO VECTOR RESHAPE CONTROL
! X          = COORDINATES OF SYSTEM NODES
! XYZ        = SPACE COORDINATES AT A POINT
! XYZ_MAX    = UPPER BOUNDS FOR SCP GEOMETRY
! XYZ_MIN    = LOWER BOUNDS FOR SCP GEOMETRY

   IF ( NP_SYS < 1 ) RETURN        ! INACTIVE NODE
   IF ( NP_SYS > MAX_NP ) STOP           &
     'EVAL_PT_FLUX_IN_SCP_PATCH: NP_SYS > MAX_NP' 
   IF ( SCP_ONLY_ONCE .AND. NODE_DONE (NP_SYS) ) then !b
     print *, 'WARNING NODE_DONE should not be true'
     N_WARN = N_WARN + 1 ; RETURN ; end if ! impossible
   STORE (1) = N_R_B * N_SPACE     ! INITIALIZE
   HIGHER    = N_R_B * N_SPACE     ! INITIALIZE

!         CONVERT IQ XYZ TO LOCAL PATCH POINT
   XYZ   = X (NP_SYS, :)
   POINT = GET_SCP_PT_AT_XYZ (XYZ, XYZ_MIN, XYZ_MAX)

!     EVALUATE PATCH INTERPOLATION AT LOCAL POINT
   CALL GEN_ELEM_SHAPE (POINT, SCP_H, SCP_N, N_SPACE, ONE)

!     INTERPOLATE PATCH NODE FIT TO PHYSICAL NODE
   FLUX_PT (1:N_R_B) = MATMUL (SCP_H (:), PATCH_FIT (:, 1:N_R_B))

!      SCATTER NODE VALUES TO SYSTEM & INCREMENT COUNTS
   SCP_COUNTS   (NP_SYS)          = SCP_COUNTS   (NP_SYS) + 1
   SCP_AVERAGES (NP_SYS, 1:N_R_B) = SCP_AVERAGES (NP_SYS, 1:N_R_B)&
                                  + FLUX_PT      (1:N_R_B)
   NODE_DONE    (NP_SYS)          = .TRUE. 

!      IF 2ND OR 3RD DERIVATIVES TO BE FOUND  
   IF ( SCP_2ND_DERIV .OR. SCP_3RD_DERIV ) THEN 

!       GET PATCH LOCAL & GLOBAL DERIVATIVES AT THIS MESH NODE 
     CALL GEN_ELEM_DERIV (POINT, SCP_DLH (1:N_SPACE, :), &
                          SCP_N, N_SPACE, ONE) 
!       MULTIPLY BY INVERSE OF PATCH'S DIAGONAL_JACOBIAN 
     DO IS = 1, N_SPACE   
       SCP_DGH (IS, :) = SCP_DLH (IS, :) / DIAGONAL_JACOBIAN (IS)
     END DO ! GLOBAL DERIVATIVES 

!       GET GRADIENT OF AVERAGED FLUX COMPONENTS
!     AND STORE N_R_B BY N_SPACE ARRAY AS A VECTOR
     DERIV_1_FLUX = RESHAPE ( TRANSPOSE ( MATMUL (         &
                              SCP_DGH, PATCH_FIT) ), STORE )

!       ADD ADDITIONAL DATA TO COLUMNS TO BE AVERAGED 
     SCP_AVERAGES (NP_SYS, (N_R_B + 1):(N_R_B + HIGHER)) =    &
          SCP_AVERAGES (NP_SYS, (N_R_B + 1):(N_R_B + HIGHER)) &
          +  DERIV_1_FLUX (1:HIGHER) 
 
!b        SAVING FIRST ELEMENT IN EACH PATCH GIVES A SEQUENTIAL
!b         ELEMENT LIST OF THE GRADIENT OF THE NODAL GRADIENTS !b
      !bIF ( LP == 1 .AND. U_PLT10 > 0 )         & !b
      !bCALL SAVE_NP_GRAD_GRAD (XYZ, AVE_2ND_DERIVS) !b

   END IF ! NEED HIGHER DERIVATIVES 
END SUBROUTINE EVAL_PT_FLUX_IN_SCP_PATCH 

SUBROUTINE EVAL_SCP_FIT_AT_PATCH_NODES (IP, NODES, X, L_IN_PATCH,   &
             MEMBERS, XYZ_MIN, XYZ_MAX, PATCH_FIT, SCP_AVERAGES,    &
             SCP_COUNTS)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!  CALCULATE THE SUPER_CONVERGENCE_PATCH AVERAGE FLUXES   
!         AT ALL ELEMENT NODES INSIDE THE PATCH
!       OPTIONALLY EVALUATE 2ND OR 3RD DERIVATIVES    
!        (AVERAGE VALUES FROM ALL PATCHES LATER)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants ! for MAX_NP, NEIGH_L, N_ELEMS, N_QP,
                     ! SCP_FIT, U_SCPR,
Use Sys_Properties_Data ! for DELETED
Use Elem_Type_Data      ! for LT_*
Use SCP_Type_Data       ! for SCP_H (SCP_N), SCP_DLH (SCP_PARM, SCP_N)
Use Interface_Header
 IMPLICIT NONE
 INTEGER,  INTENT (IN)    :: IP           ! CURRENT PATCH NUMBER
 INTEGER,  INTENT (IN)    :: NODES        (L_S_TOT, NOD_PER_EL)
 REAL(DP), INTENT (IN)    :: X            (MAX_NP, N_SPACE)
 INTEGER,  INTENT (IN)    :: L_IN_PATCH
 INTEGER,  INTENT (IN)    :: MEMBERS      (L_IN_PATCH)   
 REAL(DP), INTENT (IN)    :: XYZ_MIN      (N_SPACE), XYZ_MAX (N_SPACE)
!2bREAL(DP), INTENT (IN)    :: PATCH_FIT    (SCP_N , SCP_FIT)    
 REAL(DP), INTENT (IN)    :: PATCH_FIT    (SCP_N , N_R_B  )    !2b
 REAL(DP), INTENT (INOUT) :: SCP_AVERAGES (MAX_NP, SCP_FIT)
 INTEGER,  INTENT (INOUT) :: SCP_COUNTS   (MAX_NP)      


 REAL(DP) :: FLUX_LT (N_R_B  , NOD_PER_EL)         ! NODAL FLUXES 
 REAL(DP) :: DIAGONAL_JACOBIAN (N_SPACE)    ! PATCH GEOM JACOBIAN

 INTEGER, SAVE  :: HIGHER, STORE (1) 
 INTEGER  :: IN, LM, LP, LT, NP_SYS              ! LOOPS
 INTEGER, PARAMETER :: ONE = 1

 LOGICAL :: NODE_DONE (MAX_NP) 

! DELETED    = TRUE IF THE ELEMENT HAS BEEN DELETED
! ELEM_NODES = THE NOD_PER_EL INCIDENCES OF THE ELEMENT
! FLUX_LT    = FLUX AT ELEMENT NODES INTERPOLATED FROM A SCP 
! HIGHER     = NUMBER OF HIGHER DEGREE DERIVATIVE FLUX TERMS
! INDEX      = SYSTEM DOF NUMBERS ASSOCIATED WITH ELEMENT
! IP         = CURRENT PATCH NUMBER
! L_HOMO     = 1, IF ELEMENT PROPERTIES ARE HOMOGENEOUS
! L_IN_PATCH = NUMBER OF ELEMENTS IN CURRENT SCP
! L_SHAPE    = SHAPE FLAG FOR QUADRATURE RULE SELECTION
! L_S_TOT    = TOTAL NUMBER OF ELEMENTS & THEIR SEGMENTS
! L_TYPE     = ELEMENT TYPE NUMBER
! LT         = ELEMENT TYPE NUMBER (IF USED)
! LT_DATA    = ARRAY OF CONSTANTS FOR EACH ELEMENT TYPE
! LT_QP      = NUMBER OF QUADRATURE PTS FOR ELEMENT TYPE
! LT_USER    = OPTIONAL USER CODE ASSOCIATED WITH TYPE
! MAX_NP     = NUMBER OF SYSTEM NODES
! MEMBERS    = ELEMENT NUMBERS MACKING UP A SCP
! NOD_PER_EL = MAXIMUM NUMBER OF NODES PER ELEMENT
! NODES      = NODAL INCIDENCES OF ALL ELEMENTS
! NODE_DONE  = TRUE, IF IT HAS ONE CONTRIBUTION FROM PATCH
! N_QP       = MAXIMUM NUMBER OF QUADRATURE POINTS, >= LT_QP
! N_SPACE    = DIMENSION OF SPACE
! PATCH_FIT  = LOCAL PATCH VALUES FOR FLUX AT ITS NODES
! SCP_AVERAGES = AVERAGED FLUXES AT ALL NODES IN MESH
! SCP_CENTER_ONLY USES "CENTRAL" ELEMENT OR POINT VALUES
! SCP_FIT    = NUMBER IF TERMS BEING AVERAGED, = N_R_B MIN    
! SCP_ONLY_ONCE COMPUTES A NODE FLUX ONLY ONCE PER PATCH
! STORE      = ARRAY TO VECTOR RESHAPE CONTROL
! X          = COORDINATES OF SYSTEM NODES
! XYZ_MAX    = UPPER BOUNDS FOR SCP GEOMETRY
! XYZ_MIN    = LOWER BOUNDS FOR SCP GEOMETRY

   IF ( DEBUG_SCP ) THEN
     PRINT *,'SCP_FIT ', SCP_FIT
     PRINT *,'IP, L_IN_PATCH ', IP, L_IN_PATCH
     PRINT *,'EVAL_SCP_FIT_AT_PATCH_NODES, N_2_DER, SCP_FIT ', &
              N_2_DER, SCP_FIT ; END IF ! Debug

   STORE (1) = N_R_B * N_SPACE           ! INITIALIZE
   HIGHER    = N_R_B * N_SPACE           ! INITIALIZE
   LT = 1                                ! INITIALIZE
   NODE_DONE = .FALSE.                   ! INITIALIZE  
   IF ( SCP_PARM /= N_SPACE ) THEN 
     PRINT *,'WARNING EVAL_SCP_FIT_AT_PATCH_NODES: SCP_PARM/=N_SPACE'
     N_WARN = N_WARN + 1 ; END IF

!    IF 2ND OR 3RD DERIVATIVES TO BE FOUND EVALUATE THE PATCH'S
!            CONSTANT DIAGONAL JACOBIAN, ITS PRODUCTS
   IF ( SCP_2ND_DERIV .OR. SCP_3RD_DERIV ) THEN ! NEED JACOBIAN
     CALL PATCH_CONSTANT_JACOBIAN (XYZ_MIN, XYZ_MAX, &
                                   DIAGONAL_JACOBIAN)
     IF ( DEBUG_SCP .AND. IP == 1) PRINT *,'DIAGONAL JACOBIAN', &
        DIAGONAL_JACOBIAN
   END IF ! PATCH JACOBIAN NEEDED

   IF ( SCP_NEIGH_PT .AND. SCP_CENTER_ONLY ) THEN ! use ctr pt only
     IF ( N_PATCH /= MAX_NP ) THEN
      PRINT *,'WARNING, REMOVE KEYWORDS scp_neigh_pt OR scp_center_only'
      PRINT *,'PATCH COUNT OF ', N_PATCH, ' IS INVALID'
      N_WARN = N_WARN + 1
     END IF ! debug pt base patch options'
     CALL  EVAL_PT_FLUX_IN_SCP_PATCH (IP, X, XYZ_MIN, XYZ_MAX,  &
                    PATCH_FIT, DIAGONAL_JACOBIAN, SCP_AVERAGES, &
                    SCP_COUNTS, NODE_DONE)
     RETURN
   END IF ! USE CENTER NODE ONLY IN NODE BASED PATCH

!    INTERPOLATE AVERAGES TO CENTER OR ALL NODES IN THE PATCH
   DO LP = 1, L_IN_PATCH                 ! PATCH ELEMENT MEMBER LOOP

     IF ( SCP_CENTER_ONLY .AND. (.NOT. SCP_NEIGH_PT) .AND.   &  
          LP > 1 )  EXIT ! THIS LOOP SINCE ONLY THE "CENTER" 
                         ! ELEMENT RESULT IS WANTED FOR AVERAGE 

     LM = MEMBERS (LP)                   ! ELEMENT IN PATCH
     IF ( DELETED (LM) )  CYCLE          ! TO ACTIVE ELEMENT
     IF ( N_L_TYPE > 1) LT = L_TYPE (LM) ! GET ELEMENT TYPE NUMBER
     IF ( LT /= LAST_LT ) THEN           ! NEW ELEMENT TYPE
       CALL GET_ELEM_TYPE_DATA (LT)      ! GET CONTROLS FOR ELEM TYPE
       IF ( TYPE_APLY_ALLOC ) CALL DEALLOCATE_TYPE_APPLICATION
       CALL ALLOCATE_TYPE_APPLICATION !  for ELEM_NODES (LT_N)
       LAST_LT = LT 
     END IF ! a new element type

!-->        EXTRACT ELEMENT NODE NUMBERS
if ( LM < 1 ) stop 'impossible el 1'
     ELEM_NODES = GET_ELEM_NODES (LM, LT_N, NODES)   
     IF ( DEBUG_SCP .AND. LP == 1 ) PRINT *, 'ELEM_NODES ', &
                                              ELEM_NODES  

!       LOOP OVER ELEMENT'S NODE, INTERPOLATE FLUX
     DO IN = 1, LT_N
       NP_SYS = ELEM_NODES (IN)   ! SYSTEM NODE NUMBER
       IF ( NP_SYS < 1 ) CYCLE ! TO ACTIVE NODE
       IF ( NP_SYS > MAX_NP ) THEN ! IMPOSSIBLE              
         PRINT *, 'NP_SYS, LT_N ', NP_SYS, LT_N
         CALL IPRINT (ELEM_NODES, 1, LT_N)                   
         PRINT *,'L_S_TOT, NOD_PER_EL ', L_S_TOT, NOD_PER_EL 
         CALL IPRINT (NODES, L_S_TOT, NOD_PER_EL)
         STOP 'EVAL_SCP_FIT_AT_PATCH_NODES: NP_SYS > MAX_NP' 
       END IF ! INVALID NODE
       IF ( SCP_ONLY_ONCE .AND. NP_SYS > 0 ) THEN ! USE NODE ONLY ONCE
         IF ( NODE_DONE (NP_SYS) ) CYCLE ! TO UNUSED NODE
       END IF

!        EVALUATE ITEMS AT THIS NODE
       CALL  EVAL_PT_FLUX_IN_SCP_PATCH (NP_SYS, X, XYZ_MIN, XYZ_MAX, &
                         PATCH_FIT, DIAGONAL_JACOBIAN, SCP_AVERAGES, &
                         SCP_COUNTS, NODE_DONE)
     END DO ! LOCAL NODES

!         LIST AVERAGED NODAL FLUX VALUES FOR (FIRST) ELEMENT
     IF ( DEBUG_SCP .AND. LP == 1 ) THEN
       PRINT *, '       ELEMENT ', LM
       PRINT *, 'NODE, ', N_R_B, ' FLUX COMPONENTS' 
       DO IN = 1, LT_N
         PRINT *, ELEM_NODES (IN), FLUX_LT (1:N_R_B, IN)
       END DO ! LOCAL NODES
     END IF ! FIRST PATCH MEMBER
   END DO ! OVER PATCH MEMBERS
END SUBROUTINE EVAL_SCP_FIT_AT_PATCH_NODES 

!2b SUBROUTINE AVERAGE_SECOND_CROSS_DERIVATIVES (DERIV_2_SQ, AVE_2ND_DERIVS)
!2b ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!2b !           AVERAGE SECOND DERIVATIVE COMPONENTS
!2b ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
  !2b Use System_Constants ! for MAX_2ND, N_SPACE, N_R_B, N_2_DER
   !2b IMPLICIT NONE
   !2b REAL(DP), INTENT (IN)  :: DERIV_2_SQ (N_R_B, N_R_B) 
   !2b REAL(DP), INTENT (OUT) :: AVE_2ND_DERIVS (N_2_DER)
!2b 
!2b ! AVE_2ND_DERIVS = SECOND CROSS DERIVATIVES AVERAGES
!2b ! DERIV_2_SQ = SYMMETRIC SECOND DERIVATIVES ARRAY
!2b ! N_R_B      = NUMBER OF ROWS IN B AND E MATRICES
!2b ! N_2_DER    = NUMBER OF SECOND DERIVATIVES, 1, 3, OR 6
!2b 
!2b print *, 'AVERAGE_SECOND_CROSS_DERIVATIVES' !b
   !2b SELECT CASE (N_2_DER)
     !2b CASE (1) ! 1_D XX
      !2b AVE_2ND_DERIVS (1) = DERIV_2_SQ (1, 1)
     !2b CASE (2) ! 2_D XX, XY, YY
      !2b AVE_2ND_DERIVS (1) =  DERIV_2_SQ (1, 1)
      !2b AVE_2ND_DERIVS (2) = (DERIV_2_SQ (2, 1) + DERIV_2_SQ (1, 2))*0.5d0
      !2b AVE_2ND_DERIVS (3) =  DERIV_2_SQ (2, 2)
     !2b CASE (3) ! 3_D XX, XY, XZ, YY, YZ, ZZ
      !2b AVE_2ND_DERIVS (1) =  DERIV_2_SQ (1, 1)
      !2b AVE_2ND_DERIVS (2) = (DERIV_2_SQ (2, 1) + DERIV_2_SQ (1, 2))*0.5d0
      !2b AVE_2ND_DERIVS (3) = (DERIV_2_SQ (3, 1) + DERIV_2_SQ (1, 3))*0.5d0
      !2b AVE_2ND_DERIVS (4) =  DERIV_2_SQ (2, 2)
      !2b AVE_2ND_DERIVS (3) = (DERIV_2_SQ (3, 2) + DERIV_2_SQ (2, 3))*0.5d0
      !2b AVE_2ND_DERIVS (6) =  DERIV_2_SQ (3, 3)
     !2b CASE DEFAULT ! UN_PLANNED
      !2b PRINT *, 'AVERAGE_SECOND_CROSS_DERIVATIVES: UNKNOWN N_2_DER', N_2_DER
      !2b STOP     'AVERAGE_SECOND_CROSS_DERIVATIVES: UNKNOWN N_2_DER'
   !2b END SELECT
!2b END SUBROUTINE AVERAGE_SECOND_CROSS_DERIVATIVES 

!SUBROUTINE SAVE_FE_NP_GRAD_FLUX (XYZ, DERIV_1_FLUX) !b
! Use System_Constants ! for N_SPACE, N_R_B, U2_PLT1
!  IMPLICIT NONE
!  REAL(DP), INTENT (IN)  :: XYZ          (N_SPACE)
!  REAL(DP), INTENT (IN)  :: DERIV_1_FLUX (N_R_B * N_SPACE) ! 2nd derivs
!
!!b print *, "SAVE_FE_NP_GRAD_FLUX" !b
!   IF ( U2_PLT1 > 0 )                                &
!     WRITE (U2_PLT1, '( 10(1PE11.3) )') XYZ, DERIV_1_FLUX 
!END SUBROUTINE SAVE_FE_NP_GRAD_FLUX 

SUBROUTINE PATCH_CONSTANT_JACOBIAN (XYZ_MIN, XYZ_MAX, DIAGONAL_JACOBIAN)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!            PATCH JACOBIAN FOR EST OF 2ND DERIV FROM 1ST DERIV
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
  Use System_Constants ! for N_SPACE
  Use SCP_Type_Data    ! for SCP_SHAP      
   IMPLICIT NONE
   REAL(DP), INTENT (IN)   :: XYZ_MIN (N_SPACE), XYZ_MAX (N_SPACE)
   REAL(DP), INTENT (OUT)  :: DIAGONAL_JACOBIAN (N_SPACE)

! DIAGONAL_JACOBIAN = CONSTANT PATCH JACOBIAN
! N_SPACE           = DIMENSION OF SOLUTION SPACE
! SCP_SHAP          = PATCH SHAPE FLAG NUMBER
! XYZ_MAX           = UPPER BOUNDS FOR SCP GEOMETRY 
! XYZ_MIN           = LOWER BOUNDS FOR SCP GEOMETRY

     DIAGONAL_JACOBIAN = XYZ_MAX - XYZ_MIN

!      SCALE FOR PARAMETRIC PATCH TYPE
     SELECT CASE (SCP_SHAP)
       CASE (1) ! LINE IN NATURAL COORDINATES
         DIAGONAL_JACOBIAN = DIAGONAL_JACOBIAN * 0.5d0
       CASE (2) ! "CENTERED" UNIT TRIANGLE
         DIAGONAL_JACOBIAN = DIAGONAL_JACOBIAN * 2.d0
       CASE (3) ! QUADRILATERAL IN NATURAL COORDINATES
         DIAGONAL_JACOBIAN = DIAGONAL_JACOBIAN * 0.5d0
       CASE (4) ! HEXAHEDRON IN NATURAL COORDINATES
         DIAGONAL_JACOBIAN = DIAGONAL_JACOBIAN * 0.5d0
       CASE (5) ! "CENTERED" UNIT TETRAHEDRON
         DIAGONAL_JACOBIAN = DIAGONAL_JACOBIAN * 3.d0
       CASE (6) ! "CENTERED" UNIT-NATURAL WEDGE
         ! UNIT TRIANGLE IN X_Y PLANE TIME QUADS IN Z
         DIAGONAL_JACOBIAN (1:2) = DIAGONAL_JACOBIAN (1:2) * 2.d0 
         DIAGONAL_JACOBIAN ( 3 ) = DIAGONAL_JACOBIAN ( 3 ) * 0.5d0
       CASE DEFAULT
         PRINT *,'PATCH_CONSTANT_JACOBIAN: INVALID SHAPE ', SCP_SHAP
         STOP    'PATCH_CONSTANT_JACOBIAN: INVALID SHAPE'
     END SELECT 
     IF ( ANY (DIAGONAL_JACOBIAN <= 0.d0) ) THEN
       PRINT *,'PATCH_CONSTANT_JACOBIAN: DIAGONAL_JACOBIAN <= 0'
       STOP    'PATCH_CONSTANT_JACOBIAN: DIAGONAL_JACOBIAN <= 0'
     END IF ! FATAL ERROR
END SUBROUTINE PATCH_CONSTANT_JACOBIAN 

SUBROUTINE LIST_SCP_AVE_GRADIENTS (X, SCP_AVERAGES)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
!   LIST AVERAGED NODAL GRADIENTS AT ALL NODES IN THE MESH
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
 Use System_Constants
 IMPLICIT NONE
 REAL(DP), INTENT(IN) :: X            (MAX_NP, N_SPACE)
 REAL(DP), INTENT(IN) :: SCP_AVERAGES (MAX_NP, SCP_FIT)
 INTEGER :: IP

! L_HOMO       = 1, IF ELEMENT PROPERTIES ARE HOMOGENEOUS
! MAX_NP       = NUMBER OF SYSTEM NODES
! N_SPACE      = DIMENSION OF SPACE
! SCP_AVERAGES = AVERAGED GRADIENTS AT ALL NODES 
! SCP_FIT      = NUMBER IF TERMS BEING FIT, = N_R_B MIN       
! X            = COORDINATES OF ALL SYSTEM NODES

  PRINT *, ''
  PRINT *, '** SUPER_CONVERGENT AVERAGED NODAL GRADIENTS **'
  PRINT *,'NODE, ', N_SPACE, ' COORDINATES, ', SCP_FIT, ' GRADIENTS'

  IF ( NULL_SCP_AVE ) THEN ! Nothing to print   !b
    PRINT *, 'NOTE: ALL AVERAGES ARE ZERO, SKIP LISTING'   !b
  ELSE   !b
!b print *, 'debug NULL_SCP_AVE = ', NULL_SCP_AVE !b
    DO IP = 1, MAX_NP
      WRITE (N_PRT, 5) IP, X (IP, :), SCP_AVERAGES (IP, :)
      5 FORMAT (I6, (6(1X, 1PE12.5)) )
    END DO ! FOR ALL NODES
  END IF   !b
END SUBROUTINE LIST_SCP_AVE_GRADIENTS 

SUBROUTINE LIST_SCP_AVE_FLUXES  (X, SCP_AVERAGES)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
!   LIST AVERAGED NODAL FLUXES AT ALL NODES IN THE MESH
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
 Use System_Constants ! for MAX_NP, N_R_B, N_SPACE, SCP_FIT, 
                      ! NULL_SCP_AVE
 IMPLICIT NONE
 REAL(DP), INTENT(IN) :: X            (MAX_NP, N_SPACE)
 REAL(DP), INTENT(IN) :: SCP_AVERAGES (MAX_NP, SCP_FIT)
 INTEGER :: IP

! L_HOMO       = 1, IF ELEMENT PROPERTIES ARE HOMOGENEOUS
! MAX_NP       = NUMBER OF SYSTEM NODES
! N_SPACE      = DIMENSION OF SPACE
! N_R_B        = NUMBER OF ROWS IN B AND E MATRICES
! SCP_AVERAGES = AVERAGED FLUXES AT ALL NODES 
! SCP_FIT      = NUMBER IF TERMS BEING FIT, = N_R_B MIN       
! X            = COORDINATES OF ALL SYSTEM NODES

  PRINT *, ''
  PRINT *, '*** SUPER_CONVERGENT AVERAGED NODAL FLUXES ***'
  !b PRINT *,'NODE, ', N_SPACE, ' COORDINATES, ', N_R_B, ' FLUXES' 
  WRITE (N_PRT, 4) XYZ_NAME (1:N_SPACE), FLUX_NAME (1:N_R_B)
  4 FORMAT ('NODE,  ', 6A12, (6X, 6A12) )

  IF ( NULL_SCP_AVE ) THEN ! Nothing to print   !b
    PRINT *, 'NOTE: ALL AVERAGES ARE ZERO, SKIP LISTING'   !b
  ELSE   !b
!b print *, 'debug NULL_SCP_AVE = ', NULL_SCP_AVE !b
    DO IP = 1, MAX_NP
      WRITE (N_PRT, 5) IP, X (IP, :), SCP_AVERAGES (IP, 1:N_R_B)
      !b 5 FORMAT (I6, (6(1X, 1PE12.5)) )
      5 FORMAT (I6, (6(1X, 1PE11.4)) )
    END DO ! FOR ALL NODES
  END IF   !b
END SUBROUTINE LIST_SCP_AVE_FLUXES    

SUBROUTINE LIST_SCP_AVE_AND_EXACT_FLUXES  (X, SCP_AVERAGES)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
! LIST AVERAGED NODAL AND EXACT FLUXES AT ALL NODES IN THE MESH
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
 Use System_Constants ! for MAX_NP, N_R_B, N_SPACE, SCP_FIT,
                      ! NULL_SCP_AVE
 Use Select_Source
 IMPLICIT NONE
 REAL(DP), INTENT(IN) :: X            (MAX_NP, N_SPACE)
 REAL(DP), INTENT(IN) :: SCP_AVERAGES (MAX_NP, SCP_FIT)
 REAL(DP) :: FLUX (N_R_B)
 INTEGER  :: IP

! FLUX         = EXACT FLUX COMPONENTS AT A POINT
! L_HOMO       = 1, IF ELEMENT PROPERTIES ARE HOMOGENEOUS
! MAX_NP       = NUMBER OF SYSTEM NODES
! N_R_B        = NUMBER OF ROWS IN B AND E MATRICES
! N_SPACE      = DIMENSION OF SPACE
! SCP_AVERAGES = AVERAGED FLUXES AT ALL NODES 
! SCP_FIT      = NUMBER IF TERMS BEING FIT, = N_R_B MIN       
! X            = COORDINATES OF ALL SYSTEM NODES

  PRINT *, ''
  PRINT *,'** SUPER_CONVERGENT AVERAGED NODAL FLUXES & EXACT FLUXES **'
  !b PRINT *,'NODE, ', N_SPACE, ' COORDINATES, ', N_R_B, ' FLUXES, ', & 
           !b N_R_B, ' EXACT FLUXES'
  WRITE (N_PRT, 4) XYZ_NAME (1:N_SPACE), FLUX_NAME (1:N_R_B), &
                   EXACT_NAME (1:N_R_B)
  4 FORMAT ('NODE,  ', 10A12, (6X, 10A12) )
  
  DO IP = 1, MAX_NP
    CALL SELECT_EXACT_FLUX (X (IP, :), FLUX)
    WRITE (N_PRT, 5) IP, X (IP, :), SCP_AVERAGES (IP, 1:N_R_B), FLUX
    5 FORMAT (I6, (10(1X, 1PE10.3)) )
  END DO ! FOR ALL NODES
END SUBROUTINE LIST_SCP_AVE_AND_EXACT_FLUXES    

SUBROUTINE LIST_SCP_AVE_EX_FLUX_GRAD  (X, SCP_AVERAGES)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
! LIST AVERAGED NODAL AND EXACT FLUX GRADIENTS (2ND DERIVATIVES)
!              AT ALL NODES IN THE MESH
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
 Use System_Constants ! for MAX_NP, N_R_B, N_SPACE, SCP_FIT
 Use Select_Source
 IMPLICIT NONE
 REAL(DP), INTENT(IN) :: X            (MAX_NP, N_SPACE)
 REAL(DP), INTENT(IN) :: SCP_AVERAGES (MAX_NP, SCP_FIT)
 REAL(DP) :: FLUX_GRAD (N_R_B * N_SPACE)
 INTEGER  :: IP, HIGHER

! FLUX_GRAD    = EXACT FLUX GRAD (2ND DERIV) COMPONENTS AT A POINT
! HIGHER       = NUMBER OF HIGHER DEGREE DERIVATIVE FLUX TERMS
! L_HOMO       = 1, IF ELEMENT PROPERTIES ARE HOMOGENEOUS
! MAX_NP       = NUMBER OF SYSTEM NODES
! N_R_B        = NUMBER OF ROWS IN B AND E MATRICES
! N_SPACE      = DIMENSION OF SPACE
! SCP_AVERAGES = AVERAGED FLUXES AT ALL NODES 
! SCP_FIT      = NUMBER IF TERMS BEING FIT, = N_R_B MIN       
! X            = COORDINATES OF ALL SYSTEM NODES

  HIGHER    = N_R_B * N_SPACE
  PRINT *, ''
  PRINT *,'** SUPER_CONVERGENT AVERAGE & EXACT NODE FLUX GRAD **'
  IF ( N_SPACE == 1 ) THEN  ! use one line
    WRITE (N_PRT, 7) XYZ_NAME (1), DERIV_2_NAME (1:HIGHER), &
                                   EXACT_NAME   (1:HIGHER)
  ELSE
    WRITE (N_PRT, 7) XYZ_NAME (1:N_SPACE), DERIV_2_NAME (1:HIGHER)
    WRITE (N_PRT, 8) NO_NAME  (1:N_SPACE), EXACT_NAME   (1:HIGHER)
    7 FORMAT (' Node, ', (10A12))
    8 FORMAT ('       ', (10A12))
  END IF ! 1_D
  
  DO IP = 1, MAX_NP
    5 FORMAT (I6, (10( 1PE12.4 )) )
    CALL SELECT_EXACT_FLUX_GRAD (X (IP, :), FLUX_GRAD (1:N_2_DER)) !b
    SELECT CASE (N_SPACE) ! FOR PRETTY FORMAT
      CASE (1) ! U,xx
        WRITE (N_PRT, 5) IP, X (IP, :),                            &
          SCP_AVERAGES (IP, (N_R_B+1):(N_R_B+HIGHER)), FLUX_GRAD (1)
      CASE (2) 
        WRITE (N_PRT, 5) IP, X (IP, :),                &
          SCP_AVERAGES (IP, (N_R_B+1):(N_R_B+HIGHER))
        ! EXACT 3 (U,j),k -> U,xx U,yx=U,xy U,yy
        ! FEA   4 (U,j),k -> U,xx U,yx U,xy U,yy
        WRITE (N_PRT, 2) FLUX_GRAD (1:2), FLUX_GRAD (2:3)
      CASE (3) 
        WRITE (N_PRT, 5) IP, X (IP, :),                &
          SCP_AVERAGES (IP, (N_R_B+1):(N_R_B+HIGHER))
       ! EXACT 6 (U,j),k -> U,xx U,yx=U,xy U,zx=U,xz U,yy U,zy=U,yz U,zz
       ! FEA   9 (U,j),k -> U,xx U,yx U,zx U,xy U,yy U,zy U,xz U,yz U,zz
        WRITE (N_PRT, 3)    FLUX_GRAD !b xxx update consistent notation
      CASE DEFAULT
                 WRITE (N_PRT, 4) FLUX_GRAD (1:N_2_DER)
    END SELECT
  END DO ! FOR ALL NODES
  1 FORMAT (18X, (10(  1PE12.4 )) ) ! SPACE DEPENDENT FORMATS
  2 FORMAT (30X, (10(  1PE12.4 )) )
  3 FORMAT (42X, (10(  1PE12.4 )) )
  4 FORMAT ( 6X, (10(  1PE12.4 )) )
END SUBROUTINE LIST_SCP_AVE_EX_FLUX_GRAD    

SUBROUTINE SAVE_SCP_AVE_AND_EXACT_FLUXES (X, SCP_AVERAGES)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
!       SAVE AVERAGE and EXACT FLUXES TO A FILE
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
 Use System_Constants ! for MAX_NP, N_R_B, N_SPACE, SCP_FIT
 Use Select_Source
 IMPLICIT NONE
 REAL(DP), INTENT(IN) :: X            (MAX_NP, N_SPACE)
 REAL(DP), INTENT(IN) :: SCP_AVERAGES (MAX_NP, SCP_FIT)
 REAL(DP) :: FLUX (N_R_B)
 INTEGER  :: IP, IO_1

! FLUX         = EXACT FLUX COMPONENTS AT A POINT
! MAX_NP       = NUMBER OF SYSTEM NODES
! N_R_B        = NUMBER OF ROWS IN B AND E MATRICES
! SCP_AVERAGES = AVERAGED FLUXES AT ALL NODES 
! SCP_FIT      = NUMBER IF TERMS BEING FIT, = N_R_B MIN       
! X            = COORDINATES OF ALL SYSTEM NODES
! U_PLT5       = UNIT TO OUTPUT NODAL GRADIENTS, IF U_SCPR>0

  IF ( U_SCPR < 1 ) RETURN  ! recovery turned off
  IF ( U_PLT5 < 1 ) RETURN  ! output turned off
  OPEN (U_PLT5, FILE='scp_node_ave_fluxes.tmp', ACTION='WRITE', &
        STATUS='REPLACE', IOSTAT = IO_1)
  IF ( IO_1 > 0 ) THEN
    PRINT *, 'OPEN FAILED FOR scp_node_averages.tmp'
    RETURN
  END IF

  DO IP = 1, MAX_NP
    CALL SELECT_EXACT_FLUX (X (IP, :), FLUX)
    WRITE (U_PLT5, 5)  SCP_AVERAGES (IP, 1:N_R_B), FLUX  
    5 FORMAT ( (8(1PE13.5)) )
  END DO ! FOR ALL NODES
  CLOSE (U_PLT5)
  PRINT *, 'NOTE: SCP AVE & EXACT FLUX SAVED TO scp_node_ave_fluxes.tmp'
END SUBROUTINE SAVE_SCP_AVE_AND_EXACT_FLUXES 

SUBROUTINE SAVE_SCP_AVE_EX_FLUX_GRAD (X, SCP_AVERAGES)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
!   SAVE NODAL AVERAGE and EXACT GRADIENT OF FLUXES TO A FILE
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
 Use System_Constants ! for MAX_NP, N_R_B, N_SPACE, SCP_FIT
 Use Select_Source
 IMPLICIT NONE
 REAL(DP), INTENT(IN) :: X            (MAX_NP, N_SPACE)
 REAL(DP), INTENT(IN) :: SCP_AVERAGES (MAX_NP, SCP_FIT)

 REAL(DP) :: FLUX_GRAD (N_R_B * N_SPACE)
 INTEGER  :: IO_1, IP, HIGHER

! FLUX_GRAD    = EXACT FLUX GRAD (2ND DERIV) COMPONENTS AT A POINT
! HIGHER       = NUMBER OF HIGHER DEGREE DERIVATIVE FLUX TERMS
! L_HOMO       = 1, IF ELEMENT PROPERTIES ARE HOMOGENEOUS
! MAX_NP       = NUMBER OF SYSTEM NODES
! N_R_B        = NUMBER OF ROWS IN B AND E MATRICES
! N_SPACE      = DIMENSION OF SPACE
! SCP_AVERAGES = AVERAGED FLUXES AT ALL NODES 
! SCP_FIT      = NUMBER IF TERMS BEING FIT, = N_R_B MIN       
! X            = COORDINATES OF ALL SYSTEM NODES

  IF ( U_SCPR  < 1 ) RETURN  ! recovery turned off

!   OPEN FILES TO SAVE 2ND DERIVATIVE RESULTS
  IF ( U2_PLT1 < 1 ) RETURN  ! output turned off
    OPEN (U2_PLT1, FILE='pt_ave_grad_flux.tmp', ACTION='WRITE', &
          STATUS='REPLACE', IOSTAT = IO_1)
    IF ( IO_1 > 0 ) THEN
      PRINT *, 'WARNING: SAVE_SCP_AVE_EX_FLUX_GRAD'
      PRINT *, 'OPEN FAILED FOR pt_ave_grad_flux.tmp'
      PRINT *, 'PUT U2_PLT1 51 IN KEYWORD CONTROL'
      N_WARN = N_WARN + 1
      RETURN
    END IF
  IF ( U2_PLT2 < 1 ) RETURN  ! output turned off
    OPEN (U2_PLT2, FILE='pt_ex_grad_flux.tmp', ACTION='WRITE', &
          STATUS='REPLACE', IOSTAT = IO_1)
    IF ( IO_1 > 0 ) THEN
      PRINT *, 'WARNING: SAVE_SCP_AVE_EX_FLUX_GRAD'
      PRINT *, 'OPEN FAILED FOR pt_ex_grad_flux.tmp'
      PRINT *, 'PUT U2_PLT2 52 IN KEYWORD CONTROL'
      N_WARN = N_WARN + 1
      RETURN
    END IF

!    SAVE FEA (AND EXACT) GRADIENTS OF FLUXES
  HIGHER = N_R_B * N_SPACE
  DO IP = 1, MAX_NP
    WRITE (U2_PLT1, '( 10(1PE11.3) )')          &
      SCP_AVERAGES (IP, (N_R_B+1):(N_R_B+HIGHER)) 
    !b check that exact values available first here
    CALL SELECT_EXACT_FLUX_GRAD (X (IP, :), FLUX_GRAD)
    WRITE (U2_PLT2, '( 10(1PE11.3) )') FLUX_GRAD
  END DO ! FOR ALL NODES

  CLOSE (U2_PLT1) ; CLOSE (U2_PLT2)
  PRINT *, 'NOTE: SCP AVE FLUX GRADS SAVED TO ', &
           'pt_ave_grad_flux.tmp'
  PRINT *, 'NOTE: SCP EXACT FLUX GRADS SAVED TO ', &
           'pt_ex_grad_flux.tmp'
END SUBROUTINE SAVE_SCP_AVE_EX_FLUX_GRAD 

SUBROUTINE SAVE_SCP_NODE_AVERAGES (SCP_AVERAGES)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
! SAVE AVERAGE GRADIENTS (OR FLUXES) TO A FILE
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
 Use System_Constants ! for MAX_NP, N_R_B, N_SPACE, SCP_FIT
 IMPLICIT NONE
 REAL(DP), INTENT(IN) :: SCP_AVERAGES (MAX_NP, SCP_FIT)
 INTEGER              :: IP, IO_1

! MAX_NP       = NUMBER OF SYSTEM NODES
! N_R_B        = NUMBER OF ROWS IN B AND E MATRICES
! SCP_AVERAGES = AVERAGED FLUXES AT ALL NODES 
! SCP_FIT      = NUMBER IF TERMS BEING FIT, = N_R_B MIN       
! U_PLT5       = UNIT TO OUTPUT NODAL GRADIENTS, IF U_SCPR>0

  IF ( U_SCPR < 1 ) RETURN  ! recovery turned off
  IF ( U_PLT5 < 1 ) RETURN  ! output turned off
  OPEN (U_PLT5, FILE='scp_node_ave_fluxes.tmp', ACTION='WRITE', &
        STATUS='REPLACE', IOSTAT = IO_1)
  IF ( IO_1 > 0 ) THEN
    PRINT *, 'OPEN FAILED FOR scp_node_averages.tmp'
    RETURN
  END IF

  DO IP = 1, MAX_NP
    WRITE (U_PLT5, 5)  SCP_AVERAGES (IP, 1:N_R_B) 
    5 FORMAT ( (6(1PE13.5)) )
  END DO ! FOR ALL NODES
  CLOSE (U_PLT5)
  PRINT *, 'NOTE: SCP AVERAGES SAVED TO scp_node_ave_fluxes.tmp'
END SUBROUTINE SAVE_SCP_NODE_AVERAGES 

SUBROUTINE RESTORE_SCP_NODE_AVERAGES (SCP_AVERAGES)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
! RESTORE AVERAGE GRADIENTS (OR FLUXES) FROM A FILE
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
 Use System_Constants ! for MAX_NP, N_R_B, N_SPACE, SCP_FIT
 IMPLICIT NONE
 REAL(DP), INTENT(OUT) :: SCP_AVERAGES (MAX_NP, SCP_FIT)
 INTEGER               :: IP, IO_1

! MAX_NP       = NUMBER OF SYSTEM NODES
! N_R_B        = NUMBER OF ROWS IN B AND E MATRICES
! SCP_AVERAGES = AVERAGED FLUXES AT ALL NODES 
! SCP_FIT      = NUMBER IF TERMS BEING FIT, = N_R_B MIN       
! U_PLT5       = UNIT TO OUTPUT NODAL GRADIENTS, IF U_SCPR>0

  IF ( U_SCPR < 1 ) RETURN  ! recovery turned off
  IF ( U_PLT5 < 1 ) RETURN  ! output turned off
  OPEN (U_PLT5, FILE='scp_node_ave_fluxes.tmp', ACTION='READ', &
        STATUS='KEEP', IOSTAT = IO_1)
  IF ( IO_1 > 0 ) THEN
    PRINT *, 'OPEN FAILED FOR scp_node_averages.tmp RESTORE'
    RETURN
  END IF

  DO IP = 1, MAX_NP
    READ (U_PLT5, 5)  SCP_AVERAGES (IP, 1:N_R_B) 
    5 FORMAT ( (6(1PE13.5)) )
  END DO ! FOR ALL NODES
  CLOSE (U_PLT5)
  PRINT *, 'NOTE: SCP AVERAGES RESTORED FROM scp_node_ave_fluxes'
END SUBROUTINE RESTORE_SCP_NODE_AVERAGES 

SUBROUTINE SCP_ERROR_ESTIMATES (NODES, X, SCP_AVERAGES, DOF_SYS, &
                                ELEM_ERROR_ENERGY,               &
                                ELEM_REFINEMENT, ERR_MAX)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
!  USE INTEGRAL OF DIFFERENCE BETWEEN THE RECOVERED AVERAGE NODAL 
!  FLUXES AND ORIGINAL ELEMENT FLUXES TO ESTIMATE ELEMENT ERROR
!                    IN THE ENERGY NORM
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 

!                        References: 
! Zienkiewicz, O. C., Zhu, J. Z., "Superconvergence derivative Recovery
! Techniques and a posteriori Error Estimation in the F.E.M.",
! I.J. Num. Meth. Engr., 33, 1992, pp. 1331-1382.
! Onate, E., Bugeda, G.,"Mesh Optimality Criteria for Adaptive F. E.
! Computations", Math. of F.E. and Applications, J. R. Whiteman (ed), 
! John Wiley, 1994, pp. 121-135.
! Wiberg, N. E., "Superconvergent Patch Recovery", Adv. in Engr. 
! Software, 28, 1997, pp. 85-95.

  Use System_Constants ! for MAX_NP, NEIGH_L, N_ELEMS, N_QP, 
                       ! SCP_FIT, U_SCPR, SKIP_ERROR
  Use Select_Source
  Use Sys_Properties_Data  ! for DELETED (L_S_TOT)
  Use Elem_Type_Data  ! for: 
    ! LT_FREE, LT_GEOM, LT_N, LT_PARM, LT_QP, ELEM_NODES (LT_N),&
    ! COORD (LT_N, N_SPACE), GEOMETRY (LT_GEOM, N_SPACE),       &
    ! C (LT_FREE), D (LT_FREE), S (LT_FREE, LT_FREE),           &
    ! DLG (LT_PARM, LT_GEOM), DLG_QP (LT_PARM, LT_GEOM, LT_QP)  &
    ! DLH (LT_PARM, LT_N),    DLH_QP (LT_PARM, LT_N,    LT_QP), &
    ! DLV (LT_PARM, LT_FREE), DLV_QP (LT_PARM, LT_FREE, LT_QP), &
    ! G (LT_GEOM), G_QP (LT_GEOM, LT_QP),                       &
    ! H (LT_N),    H_QP (LT_N,    LT_QP),                       &
    ! V (LT_FREE), V_QP (LT_FREE, LT_QP),                       &
    ! PT (LT_PARM, LT_QP), WT (LT_QP), D2LH (N_2_DER, LT_N)
  Use SCP_Type_Data
  Use Interface_Header
  Use Geometric_Properties
   IMPLICIT NONE
   INTEGER,  INTENT (IN)  :: NODES             (L_S_TOT, NOD_PER_EL)
   REAL(DP), INTENT (IN)  :: X                 (MAX_NP,  N_SPACE)
   REAL(DP), INTENT (IN)  :: SCP_AVERAGES      (MAX_NP,  SCP_FIT)
   REAL(DP), INTENT (IN)  :: DOF_SYS           (N_D_FRE)
   REAL(DP), INTENT (OUT) :: ELEM_ERROR_ENERGY (N_ELEMS)
   REAL(DP), INTENT (OUT) :: ELEM_REFINEMENT   (N_ELEMS)
   REAL(DP), INTENT (OUT) :: ERR_MAX

   REAL (DP), PARAMETER :: ZERO = 0.d0

   !b INTEGER   :: IE, IN, IQ, LT, QP_LT, SCP_QP, LOC_MAX (1)
   INTEGER   :: IE, IN, IQ, LT, QP_LT, LOC_MAX (1)  !b
   INTEGER   :: I_ERROR ! NON ZERO IFF INVERSION OF E FAILS
   INTEGER   :: HIGHER  ! N_R_B * N_SPACE  

   REAL (DP) :: GLOBAL_ERROR_ENERGY
   REAL (DP) :: GLOBAL_FLUX_NORM, GLOBAL_FLUX_ERROR   
   REAL (DP) :: GLOBAL_FLUX_RMS,  GLOBAL_SOLUTION_L2
   REAL (DP) :: GLOBAL_SOLUTION_ERR
   REAL (DP) :: GLOBAL_STRAIN_ENERGY, STRAIN_ENERGY_NORM
   REAL (DP) :: ELEM_STRAIN_ENERGY,   ALLOWED_ERROR
   REAL (DP) :: ALLOWED_ERR_DENSITY,  ALLOWED_ERR_PER_EL
   REAL (DP) :: ELEM_FLUX_NORM,       ELEM_FLUX_ERROR   
   REAL (DP) :: ELEM_FLUX_RMS,        ELEM_SOLUTION_L2
   REAL (DP) :: ELEM_SOLUTION_ERR,    EL_ERR_ENERGY

   REAL (DP) :: GLOBAL_H1_NORM,  GLOBAL_H2_NORM   
   REAL (DP) :: GLOBAL_H1_ERROR, GLOBAL_H2_ERROR   
   REAL (DP) :: ELEM_H1_NORM,    ELEM_H2_NORM   
   REAL (DP) :: ELEM_H1_ERROR,   ELEM_H2_ERROR   
   REAL (DP) :: EXACT_H1_NORM,   EXACT_H2_NORM   
   REAL (DP) :: EXACT_H1_ERROR,  EXACT_H2_ERROR   
   REAL (DP) :: CROSS_D2_ERROR,  CROSS_D3_ERROR   

   REAL (DP) :: EXACT_2ND_DERIV_NORM,  GLOBAL_2ND_DERIV_NORM 
   REAL (DP) :: EXACT_2ND_DERIV_ERROR, GLOBAL_2ND_DERIV_ERROR

   REAL (DP) :: DET, DET_WT, TEMP, TEST, SCP_VOLUME, VOL

   REAL (DP) :: SOLUTION     (N_G_DOF)        ! PT SOLUTION     
   REAL (DP) :: EXACT_SOL    (N_G_DOF)        ! PT EXACT VALUE  
   REAL (DP) :: SOLUTION_ERR (N_G_DOF)        ! PT ERROR RESULT 
   REAL (DP) :: EXACT_SOL_L2

   REAL (DP) :: EXACT_FLUX_ERROR,    EXACT_FLUX_NORM 
   REAL (DP) :: EXACT_STRAIN_ENERGY, EXACT_ERR_ENERGY, EX_ERR_ENERGY

   REAL (DP), ALLOCATABLE :: DOF_EL (:, :)  ! D RECTANGULAR   

!         Automatic Arrays
   REAL (DP) :: AJ        (N_SPACE, N_SPACE)    ! JACOBIAN      
   REAL (DP) :: AJ_INV    (N_SPACE, N_SPACE)    ! JACOBIAN INVERSE
   REAL (DP) :: P_AJ      (N_2_DER, N_2_DER)    ! AJ PRODUCTS  
   REAL (DP) :: P_AJ_INV  (N_2_DER, N_2_DER)    ! PRODUCTS INVERSE  
   REAL (DP) :: B         (N_R_B,   N_EL_FRE)   ! PROB DIFFERENTIAL OP
   REAL (DP) :: E         (N_R_B,   N_R_B)      ! CONSTITUTIVE MATRIX
   REAL (DP) :: E_INVERSE (N_R_B,   N_R_B)      ! CONSTITUTIVE INVERSE
   REAL (DP) :: DGH       (N_SPACE, NOD_PER_EL) ! FE GRADIENT OF H
   REAL (DP) :: D2GH      (N_2_DER, NOD_PER_EL) ! FE 2ND DERIV OF H 
   REAL (DP) :: FLUX_LT   (N_R_B,   NOD_PER_EL) ! NODAL FLUXES
   REAL (DP) :: XYZ       (N_SPACE)             ! POINT IN SPACE
   REAL (DP) :: SIGMA_SCP (N_R_B)               ! SCP FLUX
   REAL (DP) :: SIGMA_HAT (N_R_B)               ! FEA FLUX
   REAL (DP) :: DIFF      (N_R_B)               ! FLUX DIFFERENCE
   REAL (DP) :: DERIV2_HAT(N_2_DER)             ! EL 2ND DERIVS 
   REAL (DP) :: DERIV2_SCP(N_2_DER)             ! SCP 2ND DERIVS 
   REAL (DP) :: DERIV2_LT (N_R_B * N_SPACE, NOD_PER_EL) ! SCP VALUES 
   REAL (DP) :: DERIV2_AVE(N_2_DER, NOD_PER_EL) ! SCP AVERAGES 
   REAL (DP) :: DIFF_2    (N_2_DER)             ! 2ND DERIV DIFFERENCE
!b REAL (DP) :: FLUX_GRAD (N_R_B * N_SPACE)     ! 2ND DERIV EXACT     
   REAL (DP) :: FLUX_GRAD (N_2_DER)             ! 2ND DERIV EXACT     
!b REAL (DP) :: VALUES    (N_NP_FLO)            ! PT PROPERTIES   !b
   REAL (DP) :: MEASURE   (N_ELEMS)             ! ELEMENT MEASURE
   REAL (DP) :: ELEM_ERROR_DENSITY (N_ELEMS)    ! ER*SQRT(MEAS/SCP_VOL)
   REAL (DP) :: EXAC_ERROR_ENERGY (N_ELEMS)     ! EXACT ERROR ENERGY
   REAL (DP) :: NEAR_ZERO

   LOGICAL :: ANY_FLUX_NORM = .TRUE.  !b move to System_Constants
!b LOGICAL :: CONSTANT_J    = .FALSE. !b move to System_Constants !b
   LOGICAL :: DEFINE_TERMS  = .TRUE.  !b move to System_Constants
   LOGICAL :: ENERGY_ERROR  = .TRUE.  !b move to System_Constants
   LOGICAL :: FLUX_L2       = .TRUE.  !b move to System_Constants
   LOGICAL :: FLUX_ERROR    = .TRUE.  !b move to System_Constants
   LOGICAL :: FLUX_RMS      = .TRUE.  !b move to System_Constants
   LOGICAL :: WARN_VP       = .TRUE.  !b move to System_Constants

!b  LT_N unknown at this point, set up as allocatable !b--
   !b REAL(DP) :: PRT_L_PT  (LT_N, N_NP_FLO)       ! NODE PROPERTIES  !b
   !b INTEGER  :: L_PT_PROP (LT_N, N_NP_FIX)       ! NODE PROPERTIES  !b

! ANY_FLUX_NORM = TRUE IF ANY ENERGY_ERROR, FLUX_L2, FLUX_ERROR TRUE
! B          = GRADIENT VERSUS DOF MATRIX 
! DEFINE_TERMS = TRUE TO LIST ERROR TERM DEFINITIONS
! DELETED    = TRUE IF THE ELEMENT HAS BEEN DELETED
! DGH        = GLOBAL DERIVS OF SCALAR FUNCTIONS H
! DOF_EL     = DEGREES OF FREEDOM OF THE ELEMENT
! DOF_SYS    = DEGREES OF FREEDOM OF THE SYSTEM
! E          = CONSTITUTIVE MATRIX, INVERSE IS E_INVERSE
! ELEM_ERROR_ENERGY  = ESTIMATED ELEMENT ERROR, % OF ENERGY NORM 
! ELEM_ERROR_DENSITY = ESTIMATED ELEMENT ERROR / SQRT(MEASURE)
! ELEM_NODES         = THE NOD_PER_EL INCIDENCES OF THE ELEMENT
! ELEM_REFINEMENT    = INDICATOR, >1 REFINE, <1 DE-REFINE, 1 NO CHANGE
! ENERGY_ERROR       = TRUE IF THE ENERGY NORM IS USED (DEFAULT)
! EXAC_ERROR_ENERGY  = ENERGY IN ERROR FROM EXACT SOLUTION, IF ANY
! FLUX_L2    = TRUE, IF L2 FLUX NORM IS TO BE PRINTED
! FLUX_ERROR = TRUE, IF L2 FLUX ERROR NORM IS TO BE PRINTED
! FLUX_RMS   = TRUE, IF L2 FLUX ERROR NORM DENSITY IS TO BE PRINTED
! GRAD_NORM  = TRUE, IF L2 GRADIENT NORM IS TO BE PRINTED
! INDEX      = SYSTEM DOF NUMBERS ASSOCIATED WITH ELEMENT
! L_HOMO     = 1, IF ELEMENT PROPERTIES ARE HOMOGENEOUS
! L_IN_PATCH = NUMBER OF ELEMENTS IN CURRENT SCP
! L_PT_PROP  = INTEGER PROPERTIES AT EACH ELEMENT NODE
! L_SHAPE    = SHAPE FLAG FOR QUADRATURE RULE SELECTION
! L_S_TOT    = TOTAL NUMBER OF ELEMENTS & THEIR SEGMENTS
! L_TYPE     = ELEMENT TYPE NUMBER
! LT         = ELEMENT TYPE NUMBER (IF USED)
! LT_DATA    = ARRAY OF CONSTANTS FOR EACH ELEMENT TYPE
! LT_QP      = NUMBER OF QUADRATURE PTS FOR ELEMENT TYPE
! LT_USER    = OPTIONAL USER CODE ASSOCIATED WITH TYPE
! MAX_NP     = NUMBER OF SYSTEM NODES
! MEASURE    = ELEMENT MEASURE (GENERALIZED VOLUME)
! MEMBERS    = ELEMENT NUMBERS MACKING UP A SCP
! NOD_PER_EL = MAXIMUM NUMBER OF NODES PER ELEMENT
! NODES      = NODAL INCIDENCES OF ALL ELEMENTS
! N_NP_FLO   = NUMBER OF FLOATING POINT PROPERTIES PER NODE
! N_QP       = MAXIMUM NUMBER OF QUADRATURE POINTS, >= LT_QP
! N_R_B      = NUMBER OF ROWS IN B AND E MATRICES
! N_SPACE    = DIMENSION OF SPACE
! PATCH_FIT  = LOCAL PATCH VALUES FOR FLUX AT ITS NODES
! POINT      = LOCAL POINT IN PATCH INTERPOLATION SPACE
! POINTS     = TOTAL NUMBER OF QUADRATURE POINTS IN CURRENT PATCH
! PRT_L_PT   = REAL PROPERTIES AT ELEMENT NODES
! SCP_AVERAGES = AVERAGED FLUXES AT ALL NODES IN MESH
! SCP_FIT    = NUMBER IF TERMS BEING FIT, OR AVERAGED: = N_R_B MIN    
! SCP_VOLUME = SCP VOLUME USED IN GETTING RMS VALUES
! SIGMA_HAT  = FLUX COMPONENTS AT PT FROM ORIGINAL ELEMENT 
! SIGMA_SCP  = FLUX COMPONENTS AT PT FROM SMOOTHED SCP
! SOLUTION_ERR_L2 = TRUE, IF L2 SOLUTION ERROR NORM IS TO BE PRINTED
! SOLUTION_ERR = DIFFERENCE BETWEEN SOLUTION AND EXACT_SOL 
! EXACT_SOL = VALUE FROM USER SUPPLIED ROUTINE 
! VALUES     = LOCAL VALUES OF PROPERTIES, AT XYZ
! X          = COORDINATES OF SYSTEM NODES
! XYZ        = SPACE COORDINATES AT A POINT
! XYZ_MAX    = UPPER BOUNDS FOR SCP GEOMETRY
! XYZ_MIN    = LOWER BOUNDS FOR SCP GEOMETRY

    WRITE (N_PRT, "(/, '*** BEGINNING ELEMENT ERROR ESTIMATES ***')")
    ELEM_ERROR_ENERGY = 0.d0 ; ERR_MAX = 0.d0              ! INITIALIZE
    ELEM_REFINEMENT   = 0.d0 ; EXAC_ERROR_ENERGY = 0.d0    ! INITIALIZE
    HIGHER            = N_R_B * N_SPACE                    ! CONSTANT
    NEAR_ZERO         = TINY (1.d0)                        ! CONSTANT

    IF ( SKIP_ERROR ) THEN
      PRINT *, 'NOTE: SCP_ERROR_ESTIMATES, SKIPPING ERROR ESTIMATES'
      IF ( NULL_SCP_AVE ) PRINT *, 'NOTE: SCP AVERAGES WERE ZERO'
      RETURN
    END IF ! Omit error estimates
    IF ( GRAD_BASE_ERROR ) PRINT *, 'NOTE: ERROR ESTIMATES ARE ', &
        'BASED ON GRADIENTS ONLY'

    IF ( LT_PARM /= N_SPACE ) THEN ! INVALID DATA
      PRINT *, 'RETURN. LT_PARM /= N_SPACE, SCP_ERROR_ESTIMATES'
      RETURN ! INVALID DATA
    END IF ! CONTINUUM MODEL

    IF ( SCP_2ND_DERIV .AND. N_SPACE > 2 ) THEN
      PRINT *, 'WARNING: SCP_ERROR_ESTIMATES, 2ND DERIVATIVES ', &
               'NOT CHECKED FOR THREE-DIMENSIONAL SPACE'
      N_WARN = N_WARN + 1
    END IF

! INITIALIZE
    GLOBAL_ERROR_ENERGY  = 0.d0                            
    GLOBAL_FLUX_NORM  = 0.d0 ; GLOBAL_FLUX_ERROR    = 0.d0 
    GLOBAL_FLUX_RMS   = 0.d0 ; GLOBAL_SOLUTION_L2   = 0.d0 
    SCP_VOLUME        = 0.d0 ; GLOBAL_SOLUTION_ERR  = 0.d0 
    MEASURE           = 0.d0 ; GLOBAL_STRAIN_ENERGY = 0.d0 

    EXACT_FLUX_ERROR    = 0.d0 ; EXACT_FLUX_NORM    = 0.d0 
    EXACT_STRAIN_ENERGY = 0.d0 ; EXACT_ERR_ENERGY   = 0.d0 
    EXACT_SOL           = 0.d0 ; EXACT_SOL_L2       = 0.d0 
    EXAC_ERROR_ENERGY   = 0.d0 !b

    GLOBAL_H1_NORM  = 0.d0 ;  GLOBAL_H2_NORM  = 0.d0   
    GLOBAL_H1_ERROR = 0.d0 ;  GLOBAL_H2_ERROR = 0.d0  
    ELEM_H1_NORM    = 0.d0 ;  ELEM_H2_NORM    = 0.d0 
    ELEM_H1_ERROR   = 0.d0 ;  ELEM_H2_ERROR   = 0.d0   
    EXACT_H1_NORM   = 0.d0 ;  EXACT_H2_NORM   = 0.d0   
    EXACT_H1_ERROR  = 0.d0 ;  EXACT_H2_ERROR  = 0.d0   
    CROSS_D2_ERROR  = 0.d0 ;  CROSS_D3_ERROR  = 0.d0   

    DIFF = 0.d0 ; DIFF_2 = 0.d0

!        CHECK A SPECIAL CASE
    !b IF ( NO_PROPERTIES ) CALL REAL_IDENTITY (N_R_B, E) ! INITIALIZE
    IF ( NO_PROPERTIES .OR. GRAD_BASE_ERROR ) &
         CALL REAL_IDENTITY (N_R_B, E) ! INITIALIZE
    XYZ = 0.d0 ; IE = 1                                ! DUMMY DEBUG 
  
!        CHECK FOR FLUX NORM (DEFAULT) OR GRADIENT NORM
    ANY_FLUX_NORM = (ENERGY_ERROR .OR. FLUX_L2 .OR. FLUX_ERROR)  

    IF ( ANY_FLUX_NORM ) THEN  
!        CHECK FOR POSSIBLE CONSTANT CONSTITUTIVE MATRIX
      IF ( L_HOMO == 1 ) THEN ! HOMOGENEOUS  
          XYZ = 0.d0 ; IE = 1 ! DUMMY ARGUMENTS
          !b CALL SELECT_APPLICATION_E_MATRIX (IE, XYZ, E) !b
          IF ( .NOT. GRAD_BASE_ERROR ) CALL &
                SELECT_APPLICATION_E_MATRIX (IE, XYZ, E) !b
          IF ( ENERGY_ERROR ) THEN ! NEED E INVERSE
            CALL INV_SMALL_MAT (N_R_B, E, E_INVERSE, I_ERROR)  
            IF ( I_ERROR /= 0 ) THEN ! INVERSE FAILED
              PRINT *, 'BAD HOMOGENEOUS CONSTITUTIVE MATRIX, SET TO AN'
              PRINT *, 'IDENTITY MATRIX IN SCP_ERROR_ESTIMATES'
              CALL RPRINT (E, N_R_B, N_R_B, 1)
              E_INVERSE = GET_REAL_IDENTITY (N_R_B)
            END IF ! ERROR
!b if ( ie == 1) then !b
!b print *,'E, E_INV 0 ' !b
!b call rprint(E        , N_R_B, N_R_B, 1) !b
!b call rprint(E_INVERSE, N_R_B, N_R_B, 1) !b
!b end if !b
        END IF ! INVERSE NEEDED  
      END IF ! HOMOGENEOUS MATERIAL
   END IF ! ANY FLUX NORM USED  
    
!b if ( ie == 1) then !b
!b print *,'E, E_INV 1 ' !b
!b call rprint(E        , N_R_B, N_R_B, 1) !b
!b call rprint(E_INVERSE, N_R_B, N_R_B, 1) !b
!b end if !b

!     LOOP OVER ALL STANDARD ELEMENTS
    LT = 1 ; LAST_LT = 0
!be IE = 0 ; print *,'IE, EXACT_ERR_ENERGY ', IE, EXACT_ERR_ENERGY, LT_QP !b
    DO  IE = 1, N_ELEMS  ! LOOP OVER ALL ELEMENTS  ---------------------
      IF ( DELETED (IE) ) CYCLE ! TO ACTIVE ELEMENT

!-->     GET ELEMENT TYPE NUMBER
      IF ( N_L_TYPE > 1) LT = L_TYPE (IE)  ! SAME AS LAST TYPE ?
      IF ( LT /= LAST_LT ) THEN            ! this is a new type
        CALL GET_ELEM_TYPE_DATA (LT)       ! CONTROLS FOR THIS TYPE
!x      IF ( TYPE_APLY_ALLOC ) CALL DEALLOCATE_TYPE_APPLICATION
!x      CALL ALLOCATE_TYPE_APPLICATION     ! for ELEM_NODES (LT_N)
!x      IF ( TYPE_NTRP_ALLOC ) CALL DEALLOCATE_TYPE_INTERPOLATIONS
!x      CALL ALLOCATE_TYPE_INTERPOLATIONS  ! for DLH (LT_PARM, LT_N)
        LAST_LT = LT

!-->       GET UPGRADED QUADRATURE RULE FOR PATCH "ELEMENT" TYPE 
        CALL GET_PATCH_QUADRATURE_ORDER (LT_SHAP, LT_QP, SCP_QP)

!        SINCE SCP_QP >= LT_QP MUST REALLOCATE SOME ARRAYS
        QP_LT = LT_QP   ! Copy to prevent overwrite
        LT_QP = SCP_QP  ! Remember to return to original value below
        IF ( TYPE_APLY_ALLOC ) CALL DEALLOCATE_TYPE_APPLICATION
        CALL ALLOCATE_TYPE_APPLICATION     ! for ELEM_NODES (LT_N)
        IF ( TYPE_NTRP_ALLOC ) CALL DEALLOCATE_TYPE_INTERPOLATIONS
        CALL ALLOCATE_TYPE_INTERPOLATIONS  ! for DLH (LT_PARM, LT_N)
        IF ( ALLOCATED (DOF_EL) ) DEALLOCATE (DOF_EL) !b
        ALLOCATE (DOF_EL (N_G_DOF, LT_N)) !b

 !-->     GET QUADRATURE RULE FOR PATCH "ELEMENT" TYPE
        IF ( LT_QP > 0 ) THEN
           CALL GET_ELEM_QUADRATURES
           CALL FILL_TYPE_INTERPOLATIONS
        END IF ! quadratures used
      END IF ! a new element type

!-->    EXTRACT ELEMENT NODE NUMBERS
if ( IE < 1 ) stop 'impossible el 2'
      ELEM_NODES = GET_ELEM_NODES (IE, LT_N, NODES)    

!-->    EXTRACT NODAL COORDINATES
      CALL ELEM_COORD (LT_N, N_SPACE, X, COORD, ELEM_NODES)
!be if (IE==1) print *,'IE, LT_N, LT_QP ', IE, LT_N, LT_QP !b

!-->    CALCULATE DEGREE OF FREEDOM NUMBERS
      INDEX = GET_ELEM_INDEX (LT_N, ELEM_NODES)

!-->    EXTRACT NODAL PARAMETERS OF THE ELEMENT
      D      = GET_ELEM_DOF (DOF_SYS)
      DOF_EL = RESHAPE (D, (/ N_G_DOF, LT_N /)) 
!be if (IE==1) print *,'IE, D ', D                         !b
!b if ( ie == 1) then !b
!b print *, 'INDEX ', INDEX !b 
!b print *, 'D ', D !b
!b call rprint (DOF_EL, N_G_DOF, LT_N, 1) !b
!b endif !b

      IF ( ANY_FLUX_NORM ) THEN  
!-->       CHECK FOR FIXED CONSTITUTIVE MATRIX AT ALL QP
        IF ( L_HOMO == 0 .AND. N_NP_FLO == 0 ) THEN ! CONSTANT IN IE
          XYZ = 0.d0  ! DUMMY ARGUMENT           
!b xxx ? if homo then e e-inverse exist here & are constants ? 
          !b CALL SELECT_APPLICATION_E_MATRIX (IE, XYZ, E) 
          IF ( .NOT. GRAD_BASE_ERROR ) CALL &
                SELECT_APPLICATION_E_MATRIX (IE, XYZ, E) !b
!b if ( ie == 1) then !b
!b print *,'E, E_INV 2 ' !b
!b call rprint(E        , N_R_B, N_R_B, 1) !b
!b call rprint(E_INVERSE, N_R_B, N_R_B, 1) !b
!b end if !b
          IF ( ENERGY_ERROR ) THEN ! NEED E INVERSE
            CALL INV_SMALL_MAT (N_R_B, E, E_INVERSE, I_ERROR)  
            IF ( I_ERROR /= 0 ) THEN ! INVERSE FAILED
              PRINT *, 'BAD CONSTITUTIVE MATRIX, AT ELEMENT ', IE
              PRINT *, 'SET TO IDENTITY MATRIX. SCP_ERROR_ESTIMATES,1'
              CALL RPRINT (E, N_R_B, N_R_B, 1)
              E_INVERSE = GET_REAL_IDENTITY (N_R_B)
            END IF ! ERROR
          END IF ! INVERSE NEEDED  
        END IF ! CONSTITUTIVE LAW TYPE
      END IF ! ANY FLUX NORM USED  

!-->    EXTRACT SCP NODAL FLUXES (NOW GATHER_LT_SCP_AVERAGES)
      DO IN = 1, LT_N  ! OVER NODES OF ELEMENT  ........................
        IF ( ELEM_NODES (IN) < 1 ) CYCLE ! TO VALID NODE
          FLUX_LT (1:N_R_B, IN) = SCP_AVERAGES(ELEM_NODES (IN), 1:N_R_B)
          IF ( SCP_2ND_DERIV ) THEN  ! (SECOND DERIVATIVE INFO)
!           GET NODAL AVERAGE FLUX GRADIENTS 
            DERIV2_LT (1:HIGHER, IN) = SCP_AVERAGES (ELEM_NODES (IN), &
                                           (N_R_B+1):(N_R_B+HIGHER)) 
!             AVERAGE CROSS TERMS 
            SELECT CASE (N_SPACE) 
              CASE (1) ; DERIV2_AVE = DERIV2_LT 
              CASE (2) 
                !          1    2         3
                ! EXACT 3: U,xx U,yx=U,xy U,yy
                ! FE    4: U,xx U,yx U,xy U,yy 
                !          1    2    3    4
                DERIV2_AVE (1, IN) = DERIV2_LT (1, IN) 
                DERIV2_AVE (2, IN) = SUM (DERIV2_LT (2:3, IN)) * 0.5d0
                DERIV2_AVE (3, IN) = DERIV2_LT (4, IN) 
              CASE (3)
                !          1    2         3         4    5         6
                ! EXACT 6: U,xx U,yx=U,xy U,zx=U,xz U,yy U,zy=U,yz U,zz
                ! FE    9: U,xx U,yx U,zx U,xy U,yy U,zy U,xz U,yz U,zz
                !          1    2    3    4    5    6    7    8    9
                DERIV2_AVE (1, IN) =  DERIV2_LT (1, IN)
                DERIV2_AVE (2, IN) = (DERIV2_LT (2, IN)        &
                                   +  DERIV2_LT (4, IN)) * 0.5d0
                DERIV2_AVE (3, IN) = (DERIV2_LT (3, IN)        & 
                                   +  DERIV2_LT (7, IN)) * 0.5d0
                DERIV2_AVE (4, IN) =  DERIV2_LT (5, IN)
                DERIV2_AVE (5, IN) = (DERIV2_LT (6, IN)        &
                                   +  DERIV2_LT (8, IN)) * 0.5d0
                DERIV2_AVE (6, IN) =  DERIV2_LT (9, IN)
              CASE DEFAULT 
                STOP 'SCP_ERROR_ESTIMATES, DERIV2_AVE UNDEFINED' 
            END SELECT  ! SPACE
          END IF  ! GET 2ND DERIVATIVE DATA
      END DO ! FOR NODES ON ELEMENT  ........................
           
      ELEM_FLUX_NORM         = 0.d0  ! INITIALIZE NORMS
      ELEM_FLUX_ERROR        = 0.d0  ! INITIALIZE NORMS
      ELEM_FLUX_RMS          = 0.d0  ! INITIALIZE NORMS
      ELEM_SOLUTION_L2       = 0.d0  ! INITIALIZE NORMS
      ELEM_SOLUTION_ERR      = 0.d0  ! INITIALIZE NORMS
      ELEM_STRAIN_ENERGY     = 0.d0  ! INITIALIZE 
      ELEM_ERROR_ENERGY (IE) = 0.d0  ! INITIALIZE
      EL_ERR_ENERGY          = 0.d0  ! INITIALIZE 
      EX_ERR_ENERGY          = 0.d0  ! INITIALIZE 
      VOL                    = 0.d0  ! INITIALIZE

!       LOOP OVER ERROR ESTIMATOR QUADRATURE POINTS
      DO IQ = 1, LT_QP    ! q q q q q q q q q q q q q q q q q q

!            EVALUATE INTERPOLATION FUNCTIONS
        H = GET_H_AT_QP (IQ) 

!           FIND GLOBAL COORD, XYZ = H*COORD (ISOPARAMETRIC)
        XYZ = MATMUL (H, COORD) 

!            FIND LOCAL DERIVATIVES
        DLH = GET_DLH_AT_QP (IQ) 

!           FIND JACOBIAN AT THE PT, INVERSE AND DETERMINANT
        AJ = MATMUL (DLH (1:N_SPACE, :), COORD)
        CALL INVERT_JACOBIAN (AJ, AJ_INV, DET, N_SPACE)
        IF ( DET <= ZERO ) STOP 'BAD DET, SCP_ERROR_ESTIMATES'
        DET_WT     = DET * WT(IQ)

!           UPDATE ELEMENT VOLUME MEASURE
        IF ( AXISYMMETRIC ) THEN ! X IS RADIUS
          DET_WT = DET_WT * XYZ (1) * TWO_PI
        END IF ! AXISYMMETRIC
        VOL = VOL + DET_WT

!         EVALUATE SOLUTION L2 NORM (ASSUMING C_0 FOR EACH N_G_DOF)
!b     IF ( SOLUTION_L2 ) THEN ! currently incomplete
         SOLUTION         = MATMUL (DOF_EL, H) 
         ELEM_SOLUTION_L2 = ELEM_SOLUTION_L2 + DET_WT      & 
                          * DOT_PRODUCT (SOLUTION, SOLUTION) 

!          EVALUATE APPLICATION DEPENDENT EXACT VALUE & ERROR HERE 
         IF ( USE_EXACT )  THEN
           CALL SELECT_EXACT_SOLUTION (XYZ, EXACT_SOL)
           EXACT_SOL_L2      = EXACT_SOL_L2 + DET_WT            &
                             * DOT_PRODUCT (EXACT_SOL, EXACT_SOL)
           SOLUTION_ERR      = EXACT_SOL - SOLUTION
           ELEM_SOLUTION_ERR = ELEM_SOLUTION_ERR + DET_WT            &
                            * DOT_PRODUCT (SOLUTION_ERR, SOLUTION_ERR)
           IF ( DEBUG_ERROR .OR. DEBUG_EXACT .AND. IE == 1 ) THEN
                print *, 'IQ, DET_WT, ELEM_SOLUTION_ERR ', &
                          IQ, DET_WT, ELEM_SOLUTION_ERR 
                print *, 'SOLUTION, EXACT_SOL, SOLUTION_ERR ', &
                          SOLUTION, EXACT_SOL, SOLUTION_ERR 
           END IF ! debug
         END IF ! EXACT SOLUTION GIVEN
!b     END IF ! SOLUTION_L2

!         EVALUATE GLOBAL DERIVATIVES & STANDARD GRADIENT ESTIMATE
!b if ( ie == 1) then !b
!b print *,'E, E_INV B- ' !b
!b call rprint(E        , N_R_B, N_R_B, 1) !b
!b end if !b
!be if (IE==1 .and. IQ==1) print *,'L '; call rprint(DLH,LT_PARM, LT_N,1) !b
!b print *,'IE, LT_PARM, LT_N, SIZE (DLH)', IE, LT_PARM, LT_N, SIZE (DLH)
        !b DGH  = MATMUL (AJ_INV, DLH)
        DGH (:, 1:LT_N)  = MATMUL (AJ_INV, DLH) !b
!be if (IE==1 .and. IQ==1) print *,'G '; call rprint(DGH,LT_PARM, LT_N,1) !b
        CALL SELECT_APPLICATION_B_MATRIX (DGH, XYZ, B (:, 1:LT_FREE)) !b
!be if (IE==1 .and. IQ==1) print *,'X ', XYZ  !b
!be if (IE==1 .and. IQ==1) print *,'b '; call rprint(B,N_R_B,LT_FREE,1) !b
!b if ( ie == 1) then !b
!b print *,'E, E_INV B+ ' !b
!b call rprint(E        , N_R_B, N_R_B, 1) !b
!b print *,'DLH ' !b
!b call rprint (DLH, LT_PARM, LT_N, 1) !b
!b print *,'DGH ' !b
!b call rprint (DGH, N_SPACE, NOD_PER_EL, 1) !b
!b print *,'B   ' !b
!b call rprint (B, N_R_B, N_EL_FRE, 1) !b
!b print *, 'D ', D !b
!b end if !b

!         GET LOCAL STRAINS (STORE IN SIGMA_HAT)
        SIGMA_HAT (1:N_R_B) = MATMUL (B (:, 1:LT_FREE), D (1:LT_FREE)) 
!be if (IE==1 .and. IQ==1) print *,'D ', LT_FREE, D  !b
!be if (IE==1 .and. IQ==1) print *,'B,  SIGMA_HAT ', IQ, SIGMA_HAT   !b
!b if ( ie == 1) then !b
!b print *,'hat 1 ' !b
!b print *,SIGMA_HAT (1:N_R_B) !b
!b end if !b

        IF ( ANY_FLUX_NORM ) THEN  
!         CHECK CONSTITUTIVE MATRIX AT QUADRATURE POINT    
          !b IF ( N_NP_FLO > 0 ) THEN ! PROBABLY VARIABLE     
          IF ( VARY_E) THEN ! VARIABLE E MATRIX
            !b IF ( WARN_VP ) THEN ; WARN_VP = .FALSE.
            IF ( ( N_NP_FLO > 0 ) .AND. WARN_VP ) THEN 
              WARN_VP = .FALSE.
              PRINT *, 'WARNING SCP_ERROR_ESTIMATES: NO VARIABLE ', &
                       ' NODAL PROPERTY CONSIDERED'
              N_WARN = N_WARN + 1 ! INCREMENT WARNING
            END IF ! Variable property waring

!             EXTRACT NODAL POINT PROPERTIES    !b
!b          IF ( N_NP_FLO > 0 ) THEN ! PROBABLY VARIABLE E
!b            CALL ELEM_NODE_PROPERTIES (LT_N, PRT_L_PT, L_PT_PROP, &  !b
!b                                       ELEM_NODES)    !b
!b
!b!           INTERPOLATE FOR REAL PROPERTIES AT QUADRATURE POINT    !b
!b            VALUES = MATMUL ( H, PRT_L_PT )    !b
!b          END IF ! VARIABLE NODAL POINT PROPERTIES

!             USE LOCAL VALUES IN CONSTITUTIVE EQUATIONS    
            !b CALL SELECT_APPLICATION_E_MATRIX (IE, XYZ, E) 
            IF ( .NOT. GRAD_BASE_ERROR ) CALL &
                  SELECT_APPLICATION_E_MATRIX (IE, XYZ, E) !b
            IF ( ENERGY_ERROR ) THEN ! NEED E INVERSE
              CALL INV_SMALL_MAT (N_R_B, E, E_INVERSE, I_ERROR)  
              IF ( I_ERROR /= 0 ) THEN ! INVERSE FAILED
                PRINT *, 'BAD VARIABLE CONSTITUTIVE LAW AT ELEMENT ', IE
                PRINT *, 'GAUSS POINT ', IQ, ' LOCATED AT ', XYZ
                PRINT *, 'SET TO IDENTITY MATRIX, SCP_ERROR_ESTIMATES,2'
                E_INVERSE = GET_REAL_IDENTITY (N_R_B)
              END IF ! ERROR
            END IF ! INVERSE NEEDED  
          END IF ! VARIABLE CONSTITUTIVE LAW TYPE    
        END IF ! ANY FLUX NORM USED  

!           APPLY CONSTITUTIVE RELATION (TO STRAINS IN SIGMA_SCP)       
        SIGMA_HAT (1:N_R_B) = MATMUL (E, SIGMA_HAT (1:N_R_B)) 
!be if (IE==1 .and. IQ==1) print *,'E,  SIGMA_HAT ', IQ, SIGMA_HAT   !b

!            GET SCP FLUX ESTIMATES & DIFFERENCE
        SIGMA_SCP (:) = MATMUL (H, TRANSPOSE( FLUX_LT (:, 1:LT_N)))
!b XXX Here FLUX_LT is a flux not a grad, so it contains E*grad
!b XXX already, I think, so skip following matmul  JEA 4/6/02
        !b SIGMA_SCP (1:N_R_B) = MATMUL (E, SIGMA_SCP (1:N_R_B))  !bm
        DIFF          = SIGMA_SCP - SIGMA_HAT     ! SIGMA ERROR EST
!b if ( ie == 1) then !b
!b print *,'E, E_INV 3 ' !b
!b call rprint(E        , N_R_B, N_R_B, 1) !b
!b call rprint(E_INVERSE, N_R_B, N_R_B, 1) !bm
!b print *,'hat, star, dif' !bm
!b print *,SIGMA_HAT (1:N_R_B) !bm
!b print *,SIGMA_SCP (1:N_R_B) !bm
!b print *,DIFF      (1:N_R_B) !bm
!b end if !b
        IF ( FLUX_L2    ) ELEM_FLUX_NORM  = ELEM_FLUX_NORM            &
                          + DET_WT * DOT_PRODUCT (SIGMA_SCP, SIGMA_SCP)
        IF ( FLUX_ERROR ) ELEM_FLUX_ERROR = ELEM_FLUX_ERROR &
                          + DET_WT * DOT_PRODUCT (DIFF, DIFF)

!            INCREMENT STRAIN ENERGY & ENERGY IN THE ERROR  
        TEST = DOT_PRODUCT (SIGMA_SCP, MATMUL (E_INVERSE, SIGMA_SCP))
        IF ( TEST < 0.d0 ) THEN ! PHYSICALLY IMPOSSIBLE ENERGY
          PRINT *,'SCP_ERROR_ESTIMATES: FATAL CONSTITUTIVE MATRIX (1)'
          PRINT *,'INVALID SIGMA_SCP PRODUCT, ELEM ', IE, ' POINT ', IQ
          PRINT *,'PRODUCT VALUE = ', TEST
          PRINT *,'SIGMA_SCP ARRAY: ', SIGMA_SCP
          PRINT *, 'E scp_lib' ;CALL RPRINT (E        , N_R_B, N_R_B, 1)
          PRINT *, 'E_INVERSE' ;CALL RPRINT (E_INVERSE, N_R_B, N_R_B, 1)
          STOP    'SCP_ERROR_ESTIMATES: FATAL CONSTITUTIVE MATRIX (1)'
        END IF
        ELEM_STRAIN_ENERGY = ELEM_STRAIN_ENERGY + DET_WT * TEST

        TEST = DOT_PRODUCT (DIFF, MATMUL (E_INVERSE, DIFF))
        IF ( TEST < 0.d0 ) THEN ! PHYSICALLY IMPOSSIBLE ENERGY
          PRINT *,'SCP_ERROR_ESTIMATES: FATAL CONSTITUTIVE MATRIX (2)'
          PRINT *,'INVALID DIFF PRODUCT, ELEMENT ', IE, ' POINT ', IQ
          PRINT *,'PRODUCT VALUE = ', TEST
          PRINT *,'SIGMA_SCP ARRAY: ', SIGMA_SCP
          PRINT *,'SIGMA_HAT ARRAY: ', SIGMA_HAT
          PRINT *,'DIFF      ARRAY: ', DIFF
          PRINT *, 'E scp_lib' ;CALL RPRINT (E        , N_R_B, N_R_B, 1)
          PRINT *, 'E_INVERSE' ;CALL RPRINT (E_INVERSE, N_R_B, N_R_B, 1)
          STOP    'SCP_ERROR_ESTIMATES: FATAL CONSTITUTIVE MATRIX (2)'
        END IF
        EL_ERR_ENERGY = EL_ERR_ENERGY + DET_WT * TEST
        IF (EL_ERR_ENERGY < NEAR_ZERO ) EL_ERR_ENERGY = 0.d0
!b if (EL_ERR_ENERGY < 0.d0 ) then
!b print *,'IE, EL_ERR_ENERGY ', IE, EL_ERR_ENERGY
!b end if

!           EVALUATE SECOND DERIVATIVE CONTRIBUTIONS, IF ACTIVE   
        IF ( SCP_2ND_DERIV ) THEN   
          CALL JACOBIAN_PRODUCTS (AJ, P_AJ) ! ALWAYS NEEDED   
          ! NOTE: D2GH     = P_AJ_INV (D2LH - HESSIAN * DGH)   
          !       HESSIAN  = D2LH * COORD   
          CALL INV_SMALL_MAT (N_2_DER, P_AJ, P_AJ_INV, I_ERROR)   
          IF ( I_ERROR > 0 ) THEN   
            PRINT *, 'ERROR: SCP_ERROR_ESTIMATES, NO PRODUCT INVERSE' 
            STOP     'ERROR: SCP_ERROR_ESTIMATES, NO PRODUCT INVERSE' 
          END IF ! INVERSION ERROR   
          !b D2LH = GET_D2LH_AT_QP (IQ)       ! 2ND LOCAL DERIVATIVES   
          CALL SCALAR_2ND_DERIVS (PT(:, IQ), D2LH) ! some elem only
!b if ( IQ == 1 .AND. IE == 1 ) THEN
  !b call rprint (D2LH, N_2_DER, LT_N, 1) !b
  !b call DCHECK (D2LH, LT_N, N_2_DER) !b
!b end if
          IF ( CONSTANT_J ) THEN ! HESSIAN = 0   
            D2GH = MATMUL (P_AJ_INV, D2LH)   
          ELSE  
            D2GH = MATMUL (P_AJ_INV, (D2LH            &   
                 - MATMUL (MATMUL (D2LH, COORD), DGH)))   
          END IF ! ZERO HESSIAN  
!b if(IQ==1 .AND. IE==N_ELEMS) call rprint (D2GH,N_2_DER,LT_N,1)!b

!           ACTUAL SOLUTION SECOND DERIVATIVES AT THE POINT
          DERIV2_HAT = MATMUL ( D2GH, D (1:LT_N) ) 

!            GET SCP FLUX GRADIENT ESTIMATES & DIFFERENCE 
        DERIV2_SCP (:) = MATMUL (H, TRANSPOSE( DERIV2_AVE(:, 1:LT_N)))
        DIFF_2         = DERIV2_SCP - DERIV2_HAT ! DERIV2 ERROR EST 
        ELEM_H2_NORM   = ELEM_H2_NORM                                &
                       + DET_WT * DOT_PRODUCT (DERIV2_SCP, DERIV2_SCP)
        ELEM_H2_ERROR  = ELEM_H2_ERROR                       &
                       + DET_WT * DOT_PRODUCT (DIFF_2, DIFF_2)

!               GET EXACT RESULT IF AVALIABLE
          !b CALL EXACT_FLUX_GRAD (XYZ, FLUX_GRAD) ! U,1,2 = U,2,1
          CALL SELECT_EXACT_FLUX_GRAD (XYZ, FLUX_GRAD) ! U,1,2 = U,2,1
          EXACT_H2_NORM = EXACT_H2_NORM                              &
                        + DOT_PRODUCT (FLUX_GRAD (1:N_2_DER),        &
                                       FLUX_GRAD (1:N_2_DER)) * DET_WT
!b  (/ FLUX_GRAD (1:2), FLUX_GRAD (4) /),        &
!b  (/ FLUX_GRAD (1:2), FLUX_GRAD (4) /)) * DET_WT
!b if ( IQ == 1 .AND. IE == 1 ) print *, FLUX_GRAD

!           TRUE DIFFERENCE FOR NORM
          DIFF_2 (1:N_2_DER) = FLUX_GRAD  (1:N_2_DER) &
                             - DERIV2_SCP (1:N_2_DER)
!b DIFF_2 (1:2)   = FLUX_GRAD (1:2) - DERIV2_SCP (1:2)  
!b DIFF_2 (3)     = FLUX_GRAD (4)   - DERIV2_SCP (3)  
          EXACT_H2_ERROR = EXACT_H2_ERROR                      &
                         + DET_WT * DOT_PRODUCT (DIFF_2, DIFF_2) 

!           CALCULATE CROSS DERIVATIVE ERROR NORM
            SELECT CASE (N_SPACE) 
              CASE (1) ; CROSS_D2_ERROR = 0.d0 ! NO CROSS DERIVATIVE 
              CASE (2)  
                CROSS_D2_ERROR = CROSS_D2_ERROR + DET_WT * DOT_PRODUCT &
                  (H (1:LT_N), (DERIV2_LT (2, 1:LT_N)      &
                              - DERIV2_LT (3, 1:LT_N))) **2
              CASE (3)
                ! EXACT 6: U,xx U,yx=U,xy U,zx=U,xz U,yy U,zy=U,yz U,zz
                ! FE    9: U,xx U,yx U,zx U,xy U,yy U,zy U,xz U,yz U,zz
                STOP 'SCP_ERROR_ESTIMATES, CROSS_D2 3D NOT YET DEFINED'
              CASE DEFAULT 
                STOP 'SCP_ERROR_ESTIMATES, CROSS_D2 UNDEFINED' 
            END SELECT 
        END IF ! 2ND DERIVS  

!              GET EXACT VALUES, WHEN AVAILABLE  
        IF ( USE_EXACT_FLUX )  THEN
          CALL SELECT_EXACT_FLUX (XYZ, SIGMA_SCP (1:N_R_B))  
          DIFF             = SIGMA_SCP - SIGMA_HAT ! EXACT FLUX ERROR 
          EXACT_FLUX_NORM  = EXACT_FLUX_NORM                           &
                           + DET_WT * DOT_PRODUCT (SIGMA_SCP, SIGMA_SCP)
          EXACT_FLUX_ERROR = EXACT_FLUX_ERROR                &  
                           + DET_WT * DOT_PRODUCT (DIFF, DIFF)  
          TEST = DOT_PRODUCT (SIGMA_SCP, MATMUL (E_INVERSE, SIGMA_SCP))
          IF ( TEST < 0.d0 ) THEN  
            STOP  'SCP_ERROR_ESTIMATES: NEGATIVE EXACT_STRAIN_ENERGY'  
          END IF  
          EXACT_STRAIN_ENERGY = EXACT_STRAIN_ENERGY + DET_WT * TEST  
          TEST = DOT_PRODUCT (DIFF, MATMUL (E_INVERSE, DIFF))
!be if (IE==1 .and. IQ==1) print *,'IQ, SIGMA_SCP ', IQ, SIGMA_SCP    !b
!be if (IE==1 .and. IQ==1) print *,'IQ, SIGMA_HAT ', IQ, SIGMA_HAT    !b
!be if (IE==1 .and. IQ==1) print *,'IQ, DIFF ', IQ, DIFF        !b
!be if (IE==1 .and. IQ==1) print *,'IQ, TEST, DET_WT ',IQ,TEST,DET_WT !b
          IF ( TEST < 0.d0 ) THEN  
            STOP  'SCP_ERROR_ESTIMATES: FATAL EXACT_ERR_ENERGY'  
          END IF  
          EX_ERR_ENERGY    = EX_ERR_ENERGY    + DET_WT * TEST !b
          EXACT_ERR_ENERGY = EXACT_ERR_ENERGY + DET_WT * TEST 
        END IF ! EXACT FLUXES GIVEN

!be if (IE==1 .and. IQ==1) print *,'IE, EXACT_ERR_ENERGY ', IE, EXACT_ERR_ENERGY, IQ !b
      END DO ! OVER ERROR EST QP  q q q q q q q q q q q q q q q q q q

!b if (EL_ERR_ENERGY < 0.d0 ) then
!b print *,'IE, EL_ERR_ENERGY ', IE, EL_ERR_ENERGY !bm
!b end if
      EXAC_ERROR_ENERGY (IE) = EX_ERR_ENERGY + NEAR_ZERO !b
      ELEM_ERROR_ENERGY (IE) = EL_ERR_ENERGY + NEAR_ZERO
      MEASURE (IE)           = VOL
      SCP_VOLUME             = SCP_VOLUME + VOL   

!           COMBINE AND NORMALIZE ERROR TERMS
      GLOBAL_STRAIN_ENERGY = GLOBAL_STRAIN_ENERGY + ELEM_STRAIN_ENERGY
      GLOBAL_ERROR_ENERGY  = GLOBAL_ERROR_ENERGY  + EL_ERR_ENERGY
      GLOBAL_FLUX_NORM     = GLOBAL_FLUX_NORM     + ELEM_FLUX_NORM
      GLOBAL_FLUX_ERROR    = GLOBAL_FLUX_ERROR    + ELEM_FLUX_ERROR
      GLOBAL_SOLUTION_L2   = GLOBAL_SOLUTION_L2   + ELEM_SOLUTION_L2
      GLOBAL_SOLUTION_ERR  = GLOBAL_SOLUTION_ERR  + ELEM_SOLUTION_ERR

      GLOBAL_H2_NORM       = GLOBAL_H2_NORM       + ELEM_H2_NORM
      GLOBAL_H2_ERROR      = GLOBAL_H2_ERROR      + ELEM_H2_ERROR

!b if ( ELEM_ERROR_ENERGY (IE) < 0.d0 ) then
!b print *,'IE, ELEM_ERROR_ENERGY ', IE, ELEM_ERROR_ENERGY (IE) !bm
!b STOP '1463'
!b end if
      ELEM_ERROR_ENERGY (IE) = SQRT (ELEM_ERROR_ENERGY (IE))
      EXAC_ERROR_ENERGY (IE) = SQRT (EXAC_ERROR_ENERGY (IE)) !b
 
    END DO ! OVER ALL ELEMENTS  ---------------------------------------
!       RESET LT_QP TO ITS TRUE VALUE
    LT_QP = QP_LT

!        CHECK VOLUME AGAINST GEOMETRIC PROPERTIES
    IF ( ABS(SCP_VOLUME - VOLUME) > (SCP_VOLUME / 100.d0) ) THEN 
      WRITE (N_PRT, 5) SCP_VOLUME, VOLUME ! MORE THAN 1% ERROR
      5 FORMAT ('WARNING: SCP VOLUME OF          ', 1PE12.5, /, &
      &         'DOES NOT MATCH SYSTEM VOLUME OF ', 1PE12.5, / )
      N_WARN = N_WARN + 1 ! INCREMENT WARNING
    END IF

!           FINAL GLOBAL COMBINATIONS
    GLOBAL_STRAIN_ENERGY   = GLOBAL_STRAIN_ENERGY + GLOBAL_ERROR_ENERGY
    EXACT_2ND_DERIV_NORM   = SQRT (EXACT_H2_NORM)
    GLOBAL_2ND_DERIV_NORM  = SQRT (GLOBAL_H2_NORM)
    EXACT_2ND_DERIV_ERROR  = SQRT (EXACT_H2_ERROR)
    GLOBAL_2ND_DERIV_ERROR = SQRT (GLOBAL_H2_ERROR)

    EXACT_H2_NORM  = EXACT_H2_NORM  + EXACT_FLUX_ERROR &
                   + EXACT_SOL_L2
    GLOBAL_H2_NORM = GLOBAL_H2_NORM + GLOBAL_FLUX_NORM &
                   + GLOBAL_SOLUTION_L2 

!b if (SCP_VOLUME           < 0.d0 ) STOP '1483'
    STRAIN_ENERGY_NORM = SQRT (GLOBAL_STRAIN_ENERGY)
    ALLOWED_ERROR      = STRAIN_ENERGY_NORM * (PERCENT_ERR_MAX/100.d0)
    ALLOWED_ERR_DENSITY= ALLOWED_ERROR / SQRT (SCP_VOLUME)
    ALLOWED_ERR_PER_EL = ALLOWED_ERROR / SQRT (FLOAT(N_ELEMS))

!         AVOID DIVISION BY ZERO IF THE ERROR IS ZERO
    ALLOWED_ERROR        = ALLOWED_ERROR       + NEAR_ZERO
    ALLOWED_ERR_DENSITY  = ALLOWED_ERR_DENSITY + NEAR_ZERO
    ALLOWED_ERR_PER_EL   = ALLOWED_ERR_PER_EL  + NEAR_ZERO
    ERR_MAX              = ALLOWED_ERROR
    EXACT_FLUX_ERROR     = EXACT_FLUX_ERROR    + NEAR_ZERO
    EXACT_H1_ERROR       = EXACT_H1_ERROR      + NEAR_ZERO
    !b EXACT_H2_ERROR       = EXACT_H2_ERROR      + NEAR_ZERO
    GLOBAL_ERROR_ENERGY  = GLOBAL_ERROR_ENERGY + NEAR_ZERO
    GLOBAL_FLUX_ERROR    = GLOBAL_FLUX_ERROR   + NEAR_ZERO
    GLOBAL_H1_ERROR      = GLOBAL_H1_ERROR     + NEAR_ZERO
    !b GLOBAL_H2_ERROR      = GLOBAL_H2_ERROR     + NEAR_ZERO
    GLOBAL_SOLUTION_ERR  = GLOBAL_SOLUTION_ERR + NEAR_ZERO

    GLOBAL_ERROR_ENERGY = SQRT (GLOBAL_ERROR_ENERGY)
    GLOBAL_FLUX_NORM    = SQRT (GLOBAL_FLUX_NORM)
    GLOBAL_SOLUTION_L2  = SQRT (GLOBAL_SOLUTION_L2)
    EXACT_SOL_L2        = SQRT (EXACT_SOL_L2)
    GLOBAL_SOLUTION_ERR = SQRT (GLOBAL_SOLUTION_ERR)
    IF ( FLUX_RMS .AND. SCP_VOLUME > 0.d0 )   &
      GLOBAL_FLUX_RMS   = SQRT (GLOBAL_FLUX_ERROR / SCP_VOLUME)
    GLOBAL_FLUX_ERROR   = SQRT (GLOBAL_FLUX_ERROR)

!        GET EXACT VALUES, WHEN AVAILABLE  
    EXACT_STRAIN_ENERGY = EXACT_STRAIN_ENERGY + EXACT_ERR_ENERGY  
    EXACT_STRAIN_ENERGY = SQRT (EXACT_STRAIN_ENERGY)  
    EXACT_FLUX_NORM     = SQRT (EXACT_FLUX_NORM)  
!b if (EXACT_FLUX_ERROR  < 0.d0 ) stop '1523'
!b if (EXACT_H2_ERROR    < 0.d0 ) stop '1523'
    EXACT_FLUX_ERROR    = SQRT (EXACT_FLUX_ERROR)  

    EXACT_H2_NORM       = SQRT (EXACT_H2_NORM)   
    GLOBAL_H2_NORM      = SQRT (GLOBAL_H2_NORM)   
    EXACT_H2_ERROR      = SQRT (EXACT_H2_ERROR) 
    GLOBAL_H2_ERROR     = SQRT (GLOBAL_H2_ERROR)   
    CROSS_D2_ERROR      = SQRT (CROSS_D2_ERROR)      

    PRINT *," " 
    PRINT *,"***  S_C_P  ENERGY  NORM  ERROR  ESTIMATE  DATA ***"
    PRINT *," "

!        LIST DEFINITIONS, IF REQUESTED
    IF ( DEFINE_TERMS ) CALL LIST_ERROR_TERMS

    PRINT *, " "
    PRINT *, "DOMAIN MEASURE ................", SCP_VOLUME
    PRINT *, "AVERAGE ELEMENT MEASURE .......", SCP_VOLUME / N_ELEMS
    PRINT *, "GLOBAL_SOLUTION_L2 ............", GLOBAL_SOLUTION_L2
    IF ( USE_EXACT ) THEN
      PRINT *, "EXACT_SOLUTION_L2 .............", EXACT_SOL_L2
      PRINT *, "GLOBAL_SOLUTION_ERR............", GLOBAL_SOLUTION_ERR
    END IF ! EXACT SOLUTION GIVEN
    PRINT *, " "
    PRINT *, "STRAIN_ENERGY_NORM ............", STRAIN_ENERGY_NORM
    IF ( USE_EXACT_FLUX ) PRINT *,                                  &
             "EXACT_STRAIN_ENERGY_NORM ......", EXACT_STRAIN_ENERGY 
    PRINT *, "ALLOWED_PER_CENT_ERROR ........", PERCENT_ERR_MAX
    PRINT *, "ALLOWED_GLOBAL_ERROR ..........", ALLOWED_ERROR 
    PRINT *, "ALLOWED_ERROR_DENSITY .........", ALLOWED_ERR_DENSITY
    PRINT *, "ALLOWED_ERROR_PER_ELEM ........", ALLOWED_ERR_PER_EL
    PRINT *, " "
    PRINT *, "GLOBAL_ERROR_ENERGY ...........", GLOBAL_ERROR_ENERGY
    PRINT *, "GLOBAL_ERROR_PARAMETER ........", &
              GLOBAL_ERROR_ENERGY / ALLOWED_ERROR
    PRINT *, " "
    PRINT *, "GLOBAL_FLUX_ERROR .............", GLOBAL_FLUX_ERROR
    IF ( USE_EXACT_FLUX ) PRINT *,                                 &
             "EXACT_FLUX_ERROR ..............", EXACT_FLUX_ERROR 
    PRINT *, "GLOBAL_FLUX_NORM ..............", GLOBAL_FLUX_NORM 
    IF ( USE_EXACT_FLUX ) PRINT *,                                 &
             "EXACT_FLUX_NORM ...............", EXACT_FLUX_NORM 
    PRINT *, "GLOBAL_FLUX_RMS ...............", GLOBAL_FLUX_RMS

    IF ( SCP_2ND_DERIV ) THEN
      PRINT *," "
      PRINT *,"EXACT_2ND_DERIV_NORM .........", EXACT_2ND_DERIV_NORM
      PRINT *,"GLOBAL_2ND_DERIV_NORM ........", GLOBAL_2ND_DERIV_NORM
      PRINT *,"EXACT_2ND_DERIV_ERROR ........", EXACT_2ND_DERIV_ERROR
      PRINT *,"GLOBAL_2ND_DERIV_ERROR .......", GLOBAL_2ND_DERIV_ERROR
      PRINT *,"EXACT_H2_NORM ................", EXACT_H2_NORM       
      PRINT *,"GLOBAL_H2_NORM ...............", GLOBAL_H2_NORM       
      PRINT *,"EXACT_H2_ERROR ...............", EXACT_H2_ERROR   
      PRINT *,"GLOBAL_H2_ERROR ..............", GLOBAL_H2_ERROR       
      PRINT *,"CROSS_D2_ERROR ...............", CROSS_D2_ERROR       
    END IF ! 2nd derivatives

!     CONVERT TO ELEMENT ERROR DENSITY
    WHERE ( MEASURE > 0.d0 )
      ELEM_ERROR_DENSITY = ELEM_ERROR_ENERGY / SQRT (MEASURE)
    ELSEWHERE
      ELEM_ERROR_DENSITY = 0.d0
    END WHERE

!       LIST AVERAGE TOTAL ERROR
    TEMP = STRAIN_ENERGY_NORM / 100.d0
    TEST    = SUM ( ELEM_ERROR_ENERGY ) / TEMP
    WRITE (N_PRT,  '("TOTAL % ERROR IN ENERGY NORM = ", 1PE8.2)') TEST
    !b WRITE (U_HIST, '("% ERROR IN ENERGY NORM = ", 1PE9.3)') TEST
    !b WRITE (U_LOG , '("% ERROR IN ENERGY NORM = ", 1PE9.3)') TEST
    PRINT *, " "

!        LIST MAXIMUMS
    TEST    = MAXVAL ( ELEM_ERROR_ENERGY )
    LOC_MAX = MAXLOC ( ELEM_ERROR_ENERGY )
    PRINT *, " "
    WRITE (N_PRT, '("MAXIMUM ELEMENT ENERGY ERROR OF ", 1PE8.2)') TEST
    WRITE (N_PRT, '("OCCURS IN ELEMENT ", I6)') LOC_MAX (1)
    PRINT *, " "
    TEST    = MAXVAL ( ELEM_ERROR_DENSITY )
    LOC_MAX = MAXLOC ( ELEM_ERROR_DENSITY )
    WRITE (N_PRT, '("MAXIMUM ENERGY ERROR DENSITY OF ", 1PE8.2)') TEST
    WRITE (N_PRT, '("OCCURS IN ELEMENT ", I6)') LOC_MAX (1)
    
!        FINALLY, CONVERT REFINEMENT TO TRUE REFINEMENT PARAMETER
    ELEM_REFINEMENT =  ELEM_ERROR_DENSITY * SQRT ( SCP_VOLUME ) &
                    / ALLOWED_ERROR
    WRITE (N_PRT, '("WITH A REFINEMENT PARAMETER OF ", 1PE8.2)')  &
           TEST / ALLOWED_ERR_DENSITY

    PRINT *, " "
    PRINT *, "--------------------------------------------------------"
    PRINT *, "             ERROR IN        % ERROR IN      REFINEMENT"
    PRINT *, "ELEMENT,     ENERGY_NORM,    ENERGY_NORM,    PARAMETER"  
    PRINT *, "--------------------------------------------------------"
    
    DO IE = 1, N_ELEMS ! LOOP OVER ALL ELEMENTS -----------------------
      IF ( DELETED (IE) ) CYCLE ! TO ACTIVE ELEMENT
        IF ( ELEM_ERROR_ENERGY  (IE) > ALLOWED_ERR_PER_EL .OR.  &
             ELEM_ERROR_DENSITY (IE) > ALLOWED_ERR_DENSITY ) THEN
           WRITE (N_PRT,"(I8, 3(1PE16.4),6X,A)")   &
            IE, ELEM_ERROR_ENERGY (IE),            &
            ELEM_ERROR_ENERGY (IE) / TEMP,         &
            ELEM_REFINEMENT (IE), "Refine"
        ELSE
           WRITE (N_PRT,"(I8, 3(1PE16.4))")        &
            IE, ELEM_ERROR_ENERGY (IE),            &
            ELEM_ERROR_ENERGY (IE) / TEMP,         &
            ELEM_REFINEMENT (IE)
        END IF
    END DO ! FOR ALL ELEMENTS ----------------------------------------

!     CONVERT TO % ERROR IN ENERGY NORM * 100
    ELEM_ERROR_ENERGY = ELEM_ERROR_ENERGY / TEMP !b

    WRITE (N_PRT, "(/, 'ELEMENT ERROR ESTIMATES COMPLETED', /)")
!b  following caused errors if eigen error analysis, move to main
    !b IF ( TYPE_APLY_ALLOC ) CALL DEALLOCATE_TYPE_APPLICATION
    !b IF ( TYPE_NTRP_ALLOC ) CALL DEALLOCATE_TYPE_INTERPOLATIONS

    IF ( SAVE_NEW_MESH ) CALL FORM_NEW_EL_SIZES (NODES, X, MEASURE, &
                                                 ELEM_REFINEMENT)   
!   ARE ENERGY ERRORS TO BE SAVED !b
    IF ( .NOT. USE_EXACT_FLUX ) RETURN 
      EXAC_ERROR_ENERGY = 100 * EXAC_ERROR_ENERGY / EXACT_STRAIN_ENERGY
      IF ( U_PLT6 > 0 ) THEN !b
        CALL SAVE_EXAC_ERROR_ESTIMATES (EXAC_ERROR_ENERGY) !b
      END IF ! SAVE
      IF ( U_PLT7 > 0 ) CALL SAVE_PT_AVE_EXACT_ERROR (NODES, &
                               EXAC_ERROR_ENERGY) !b
END SUBROUTINE SCP_ERROR_ESTIMATES 

SUBROUTINE LIST_ERROR_TERMS	
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
!      DEFINE TERMS IN SCP ERROR ESTIMATE OUTPUT
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
  PRINT *,"    ---------------------------------------"
  PRINT *,"    DEFINITIONS OF ERROR MEASURES & REPORTS"
  PRINT *,"    ---------------------------------------"
  PRINT *,"MEASURE    = GENERALIZED VOLUME"
  PRINT *,"R_S_L_I_Q  = SQRT( SUM ( EL_INTEGRAL (Q) ))"
  PRINT *,"L2         = SQRT( SUM ( EL_INTEGRAL (Q^2) ))"
  PRINT *,"SIGMA_HAT  = FLUX COMPONENTS FROM ORIGINAL ELEMENT"
  PRINT *,"SIGMA_SCP  = FLUX COMPONENTS FROM SMOOTHED SCP"
  PRINT *,"FLUX_ERR   = SIGMA_SCP - SIGMA_HAT"
  PRINT *,"E_INVERSE  = INVERSE OF CONSTITUTIVE MATRIX"
  PRINT *,"ERROR_ENERGY_DENSITY   = FLUX_ERR'  * E_INVERSE * FLUX_ERR"
  PRINT *,"STRAIN_ENERGY_DENSITY  = SIGMA_SCP' * E_INVERSE * SIGMA_SCP"
  PRINT *,"GLOBAL_SOLUTION_L2     = L2 (SOLUTION)"
  PRINT *,"GLOBAL_ENERGY_ERROR    = L2 (ERROR_ENERGY_DENSITY)"
  PRINT *,"GLOBAL_FLUX_NORM       = L2 (STRAIN_ENERGY_DENSITY) "
  PRINT *,"STRAIN_ENERGY_NORM     = R_S_L_I_Q (STRAIN_ENERGY_DENSITY"
  PRINT *,"                       + ERROR_ENERGY_DENSITY)"
  PRINT *,"ALLOWED_GLOBAL_ERROR   = STRAIN_ENERGY_NORM * ALLOWED_%"
  PRINT *,"ALLOWED_ERROR_DENSITY  = (ALLOWED_GLOBAL_ERROR "
  PRINT *,"                       / SQRT (DOMAIN_MEASURE))"
  PRINT *,"ALLOWED_ERROR_PER_ELEM = (ALLOWED_GLOBAL_ERROR "    
  PRINT *,"                       / SQRT (NUMBER_OF_ELEMS))"
  PRINT *,"GLOBAL_ERROR_PARAMETER = GLOBAL_ENERGY_ERROR"
  PRINT *,"                       / (ALLOWED_% STRAIN_ENERGY_NORM)"
  PRINT *,"REFINEMENT_PARAMETER   = L2 (ELEM ERROR_ENERGY_DENSITY)"
  PRINT *,"           / (ALLOWED_ERROR_DENSITY * SQRT(ELEM MEASURE))"
  PRINT *,"           <1 DE-REFINE, >1 REFINE, =1 NO CHANGE"
END SUBROUTINE LIST_ERROR_TERMS	

SUBROUTINE SAVE_ELEM_ERROR_ESTIMATES (ELEM_ERROR_ENERGY)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
! Save element error estimate for each element for plotting
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
  Use System_Constants ! for N_ELEMS, U_PLT6
   IMPLICIT NONE
   REAL(DP), INTENT (IN) :: ELEM_ERROR_ENERGY (N_ELEMS)
   !b REAL(DP)              :: TEMP 
   INTEGER               :: IE, IO_1

! ELEM_ERROR_ENERGY = ESTIMATED ERROR, % OF ENERGY NORM * 100
! N_ELEMS           = NUMBER OF ELEMENTS IN SYSTEM
! U_PLT6            = UNIT TO OUTPUT ELEMENT ERRORS FOR PLOTS

  IF ( U_PLT6 < 1 ) RETURN  ! output turned off
  OPEN (U_PLT6, FILE='el_error_est.tmp', ACTION='WRITE', &
        STATUS='REPLACE', IOSTAT = IO_1)
  IF ( IO_1 > 0 ) THEN
    PRINT *, 'WARNING: OPEN FAILED FOR el_error_est.tmp'
    N_WARN = N_WARN + 1 ; RETURN ; END IF

    DO IE = 1, N_ELEMS
      WRITE (U_PLT6, '( 1PE13.4 )') ELEM_ERROR_ENERGY (IE) 
    END DO ! FOR ALL ELEMENTS
    CLOSE (U_PLT6)
    PRINT *, 'NOTE: ELEM ERROR ESTIMATES SAVED TO el_error_est.tmp'
END SUBROUTINE SAVE_ELEM_ERROR_ESTIMATES 

SUBROUTINE SAVE_EXAC_ERROR_ESTIMATES (EXAC_ERROR_ENERGY)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
! Save element exact error for each element for plotting
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
  Use System_Constants ! for N_ELEMS, U_ERR_EX
   IMPLICIT NONE
   REAL(DP), INTENT (IN) :: EXAC_ERROR_ENERGY (N_ELEMS)
   INTEGER :: IE, IO_1, U_ERR_EX, GET_NEXT_IO_UNIT

! EXAC_ERROR_ENERGY = EXACT ERROR, % OF ENERGY NORM * 100
! N_ELEMS           = NUMBER OF ELEMENTS IN SYSTEM
! U_ERR_EX          = UNIT TO OUTPUT ELEMENT ERRORS FOR PLOTS

  U_ERR_EX = GET_NEXT_IO_UNIT ()              
  OPEN (U_ERR_EX, FILE='ex_error_est.tmp', ACTION='WRITE', &
        STATUS='REPLACE', IOSTAT = IO_1)
  IF ( IO_1 > 0 ) THEN
    PRINT *, 'OPEN FAILED FOR ex_error_est.tmp'
    N_WARN = N_WARN + 1 ; RETURN ; END IF

    DO IE = 1, N_ELEMS
      WRITE (U_ERR_EX, '( 1PE13.4 )') EXAC_ERROR_ENERGY (IE) 
    END DO ! FOR ALL ELEMENTS
    CLOSE (U_ERR_EX)
    PRINT *, 'NOTE: EXACT ERROR ESTIMATES SAVED TO ex_error_est.tmp'
END SUBROUTINE SAVE_EXAC_ERROR_ESTIMATES 

SUBROUTINE SAVE_PT_AVE_ERROR_EST (NODES, ELEM_ERROR_ENERGY)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
! Save nodal averages of the element error estimates for 
!            optional plotting as a surface
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
  Use System_Constants ! for L_S_TOT, MAX_NP, NOD_PER_EL, N_ELEMS 
  Use Sys_Properties_Data  ! for DELETED
   IMPLICIT NONE
   INTEGER,  INTENT (IN) :: NODES             (L_S_TOT, NOD_PER_EL)
   REAL(DP), INTENT (IN) :: ELEM_ERROR_ENERGY (N_ELEMS)

   INTEGER   :: IE, IN, LN, IO_1
   REAL (DP) :: ELEM_ERROR   
!         Automatic Arrays
   REAL (DP) :: PT_AVERAGE (MAX_NP), PT_COUNT (MAX_NP)          
    
! DELETED    = TRUE IF ELEMENT HAS BEEN DELETED FROM MESH
! ELEM_ERROR_ENERGY = ESTIMATED ERROR, % OF ENERGY NORM
! L_S_TOT    = N_ELEMS + N_MIXED + N_F_SEG, TOTAL GENERALIZED ELEMENTS
! MAX_NP     = NUMBER OF SYSTEM NODES
! NODES      = NODE INCIDENCES OF ALL ELEMENTS
! NOD_PER_EL = MAXIMUM NUMBER OF NODES PER ELEMENT
! N_ELEMS    = NUMBER OF ELEMENTS IN SYSTEM
! PT_AVERAGE =  AVE NODAL ESTIMATED ERROR, % OF ENERGY NORM
! U_PLT7     = UNIT FOR ELEMENT ERRORS AVERAGED AT NODES, IF U_SCPR>0

    IF ( U_PLT7 < 1 ) THEN
      RETURN  ! output turned off
    END IF
    OPEN (U_PLT7, FILE='pt_ave_error_est.tmp', ACTION='WRITE', &
          STATUS='REPLACE', IOSTAT = IO_1)
    IF ( IO_1 > 0 ) THEN
      PRINT *, 'OPEN FAILED FOR pt_ave_error_est.tmp'
      N_WARN = N_WARN + 1 ; RETURN ; END IF

    PT_AVERAGE = 0.d0 ; PT_COUNT = 0.d0 ! INITIALIZE 

    DO  IE = 1, N_ELEMS  ! LOOP OVER ALL ELEMENTS
      IF ( DELETED (IE) ) CYCLE ! TO ACTIVE ELEMENT

!        ADD TO LOCAL NODES
        ELEM_ERROR = ELEM_ERROR_ENERGY (IE)
        DO IN = 1, NOD_PER_EL
          IF ( NODES (IE, IN) < 1 ) CYCLE ! TO A VALID NODE
            LN                = NODES (IE, IN)
            PT_AVERAGE ( LN ) = PT_AVERAGE ( LN ) + ELEM_ERROR
            PT_COUNT   ( LN ) = PT_COUNT   ( LN ) + 1.d0      
        END DO ! OVER NODES
    END DO ! FOR ALL ELEMENTS

!      GET NODAL AVERAGE OF ELEMENT ERROR ESTIMATES 
    WHERE ( PT_COUNT > 0.d0 ) PT_AVERAGE = PT_AVERAGE / PT_COUNT

!      SAVE NODAL AVERAGES OF ESTIMATES 
    IF ( IO_1 == 0 ) THEN
      DO IN = 1, MAX_NP 
        WRITE (U_PLT7, '( 1PE13.5 )') PT_AVERAGE (IN)
      END DO ! OVER ALL NODES
      CLOSE (U_PLT7)
      PRINT *,'NOTE: AVE NODE ERROR EST SAVED TO pt_ave_error_est.tmp'
    END IF ! FILE WAS OPEN
END SUBROUTINE SAVE_PT_AVE_ERROR_EST 

SUBROUTINE SAVE_PT_AVE_EXACT_ERROR (NODES, EXAC_ERROR_ENERGY)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
! Save nodal averages of the exact element error for 
!            optional plotting as a surface
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
  Use System_Constants ! for L_S_TOT, MAX_NP, NOD_PER_EL, N_ELEMS 
  Use Sys_Properties_Data  ! for DELETED
   IMPLICIT NONE
   INTEGER,  INTENT (IN) :: NODES             (L_S_TOT, NOD_PER_EL)
   REAL(DP), INTENT (IN) :: EXAC_ERROR_ENERGY (N_ELEMS)

   INTEGER   :: IE, IN, LN, IO_1, U_EXACT, GET_NEXT_IO_UNIT
   REAL (DP) :: EXAC_ERROR   
!         Automatic Arrays
   REAL (DP) :: PT_AVERAGE (MAX_NP), PT_COUNT (MAX_NP)          
    
! DELETED    = TRUE IF ELEMENT HAS BEEN DELETED FROM MESH
! EXAC_ERROR_ENERGY = EXACT ERROR, % OF ENERGY NORM
! L_S_TOT    = N_ELEMS + N_MIXED + N_F_SEG, TOTAL GENERALIZED ELEMENTS
! MAX_NP     = NUMBER OF SYSTEM NODES
! NODES      = NODE INCIDENCES OF ALL ELEMENTS
! NOD_PER_EL = MAXIMUM NUMBER OF NODES PER ELEMENT
! N_ELEMS    = NUMBER OF ELEMENTS IN SYSTEM
! PT_AVERAGE =  AVE NODAL ESTIMATED ERROR, % OF ENERGY NORM
! U_EXACT    = UNIT FOR ELEMENT ERRORS AVERAGED AT NODES

    U_EXACT = GET_NEXT_IO_UNIT ()
    OPEN (U_EXACT, FILE='pt_ave_ex_error.tmp', ACTION='WRITE', &
          STATUS='REPLACE', IOSTAT = IO_1)
    IF ( IO_1 > 0 ) THEN
      PRINT *, 'OPEN FAILED FOR pt_ave_ex_error.tmp ON ', U_EXACT
      N_WARN = N_WARN + 1 ; RETURN
    END IF

    PT_AVERAGE = 0.d0 ; PT_COUNT = 0.d0 ! INITIALIZE 

    DO  IE = 1, N_ELEMS  ! LOOP OVER ALL ELEMENTS
      IF ( DELETED (IE) ) CYCLE ! TO ACTIVE ELEMENT

!        ADD TO LOCAL NODES
        EXAC_ERROR = EXAC_ERROR_ENERGY (IE)
        DO IN = 1, NOD_PER_EL
          IF ( NODES (IE, IN) < 1 ) CYCLE ! TO A VALID NODE
            LN                = NODES (IE, IN)
            PT_AVERAGE ( LN ) = PT_AVERAGE ( LN ) + EXAC_ERROR
            PT_COUNT   ( LN ) = PT_COUNT   ( LN ) + 1.d0      
        END DO ! OVER NODES
    END DO ! FOR ALL ELEMENTS

!      GET NODAL AVERAGE OF ELEMENT ERROR ESTIMATES 
    WHERE ( PT_COUNT > 0.d0 ) PT_AVERAGE = PT_AVERAGE / PT_COUNT

!      SAVE NODAL AVERAGES OF ESTIMATES 
    IF ( IO_1 == 0 ) THEN
      DO IN = 1, MAX_NP 
        WRITE (U_EXACT, '( 1PE13.5 )') PT_AVERAGE (IN)
      END DO ! OVER ALL NODES
      CLOSE (U_EXACT)
      PRINT *,'NOTE: AVE EXACT ERROR SAVED TO pt_ave_ex_error.tmp'
    END IF ! FILE WAS OPEN
END SUBROUTINE SAVE_PT_AVE_EXACT_ERROR 

SUBROUTINE GET_PATCH_QUADRATURE_ORDER (LT_SHAP, LT_QP, SCP_QP)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
!   GET UPGRADED QUADRATURE RULE FOR PATCH "ELEMENT" TYPE 
! NOTE: ANY PATCH OF ELEMENTS SURROUNDING A SUBJECT ELEMENT
!       HAS THE SAME SHAPE AS THE SUBJECT ELEMENT, LT_SHAP, 
!       BUT IT REQUIRES A HIGHER ORDER RULE THAN LT_QP
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
 IMPLICIT NONE        
 INTEGER, INTENT (IN)  :: LT_SHAP, LT_QP
 INTEGER, INTENT (OUT) :: SCP_QP

! LT_SHAP    = ELEMENT TYPE SHAPE FLAG NUMBER
!              1=LINE 2=TRI 3=QUAD 4=HEX 5=TET 6=WEDGE
! LT_QP      = NUMBER OF QUADRATURE POINTS FOR ELEMENT TYPE
! SCP_QP     = NUMBER OF QUADRATURE POINTS NEEDED IN A SCP

   SELECT CASE (LT_SHAP)  ! BRANCH ON PATCH SHAPE
     CASE (1); SCP_QP = LT_QP + 1  ! LINE
     CASE (2)                      ! TRIANGLE
!    FOR  N_QP = 1,3,4,6,7,12,13,16,19,25,27,33,37,42,48,52,61
!       DEGREE = 1,2,3,4,5, 6, 7, 8, 9,10,11,12,13,14,15,16,17
       SELECT CASE (LT_QP)
         CASE ( 1) ; SCP_QP = 3
         CASE ( 3) ; SCP_QP = 4
         CASE ( 4) ; SCP_QP = 6
         CASE ( 6) ; SCP_QP = 7
         CASE ( 7) ; SCP_QP = 12
         CASE (12) ; SCP_QP = 13
         CASE (13) ; SCP_QP = 16
         CASE (16) ; SCP_QP = 19
         CASE (19) ; SCP_QP = 25
         CASE (25) ; SCP_QP = 27
         CASE (27) ; SCP_QP = 33
         CASE (33) ; SCP_QP = 37
         CASE (37) ; SCP_QP = 42
         CASE (42) ; SCP_QP = 48
         CASE (48) ; SCP_QP = 52
         CASE (52) ; SCP_QP = 61
         CASE DEFAULT
           STOP 'GET_PATCH_QUADRATURE_ORDER: DUNAVANT QUADRATURE LIMIT'
       END SELECT
     CASE (3)                      ! QUADRILATERAL
       SCP_QP = NINT (( SQRT (FLOAT(LT_QP)) + 1) **2)
     CASE (4)                      ! HEXAHEDRAL
       SCP_QP = NINT (( (FLOAT(LT_QP))**1./3. + 1) **3)
     CASE (5); SCP_QP = LT_QP + 4  ! TETRAHEDRAL
     CASE (6); SCP_QP = LT_QP + 6  ! WEDGE
     CASE DEFAULT 
       STOP 'INVALID LT_SHAP, CALC_SCP_AVE_NODE_GRADIENTS'
   END SELECT ! BASED ON PATCH SHAPE
END SUBROUTINE GET_PATCH_QUADRATURE_ORDER 

SUBROUTINE SAVE_QP_FLUX_FOR_SCP (IE, IQ, XYZ, E, STRAIN)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! SAVE ELEMENT FLUXES AT QUADRATURE POINTS TO SCP FILE      
! Usually LIST_ELEM_FLUXES is used instead
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants ! for GRAD_BASE_ERROR
Use Elem_Type_Data   
 IMPLICIT NONE
!    ALWAYS USED
 INTEGER,  INTENT(IN) :: IE, IQ  ! ELEMENT & QUADRATURE PT
 REAL(DP), INTENT(IN) :: XYZ (N_SPACE), E (N_R_B, N_R_B), &
                         STRAIN (N_R_B + 2)
 REAL(DP) :: STRESS (N_R_B)  ! Automatic array

!                 VARIABLES:
! E         = ELEMENT CONSTITUTIVE MATRIX AT GAUSS POINT
! IE        = ELEMENT NUMBER
! IQ        = QUADRATURE POINT NUMBER IN IE
! N_R_B     = NUMBER OF ROWS IN B AND E MATRICES
! N_SPACE   = DIMENSION OF SPACE
! STRAIN    = STRAIN OR GRADIENT VECTOR
! STRESS    = FLUX COMPONENTS, = E * STRAIN
! XYZ       = SPACE COORDINATES AT A POINT

    !b STRESS = MATMUL (E, STRAIN (1:N_R_B))   4/6/02 XXX
    !b currently noy called (see LIST_ELEM_FLUXES)
    IF ( GRAD_BASE_ERROR ) THEN
      STRESS (1:N_R_B) = STRAIN (1:N_R_B)
    ELSE
      STRESS (1:N_R_B) = MATMUL (E, STRAIN (1:N_R_B))
    END IF ! GRAD VS FLUX
    RECORD_NUMBER = RECORD_NUMBER + 1
    SCP_RECORD_NUMBER (IE, IQ) = RECORD_NUMBER
    WRITE (U_SCPR, REC = RECORD_NUMBER) XYZ, STRESS (1:N_R_B) 
    !b WRITE (U_SCPR, REC = RECORD_NUMBER) XYZ, E, STRESS (1:N_R_B) 
END SUBROUTINE SAVE_QP_FLUX_FOR_SCP

SUBROUTINE FORM_NEW_EL_SIZES (NODES, X, MEASURE, ELEM_REFINEMENT)   
! ---------------------------------------------------------------
!   AVERAGE NEW ELEMENT SIZES FOR ADAPTIVE MESH GENERATOR
!    H_NEW = H_OLD / (ELEM_REFINEMENT ** (1. / EL_DEG))
!       USUALLY CALLED WHEN SAVE_NEW_MESH IS TRUE
! ---------------------------------------------------------------
Use System_Constants     ! for MAX_NP, NOD_PER_EL, N_ELEMS, N_SPACE
                         ! L_SHAPE, AVE_H_WT
Use Elem_Type_Data       ! for: LT_N, LT_SHA, ELEM_NODES (LT_N)
Use Geometric_Properties ! for VOLUME
Use Sys_Properties_Data  ! for DELETED
IMPLICIT NONE
  INTEGER,  INTENT (IN) :: NODES            (L_S_TOT, NOD_PER_EL)
  REAL(DP), INTENT (IN) :: X                (MAX_NP,  N_SPACE)
  REAL(DP), INTENT (IN) :: MEASURE          (N_ELEMS)          
  REAL(DP), INTENT (IN) :: ELEM_REFINEMENT  (N_ELEMS)

  INTEGER  :: NODE_AVE_COUNT   (MAX_NP) 
  REAL(DP) :: NODE_AVE_EL_SIZE (MAX_NP), H_SCALE (N_ELEMS), LOW_R
  
  INTEGER  :: IE, U_NEW_H, IO_1, LT, J, K, U_NEW_S
  INTEGER  :: EL_DEG, GET_SCALAR_DEGREE, GET_NEXT_IO_UNIT

  REAL(DP), PARAMETER :: TRI_WEDGE = 2.309d0, TET = 7.994d0, &
                         AVE_TO_LIMIT = 1.d-2
  REAL(DP) :: H_OLD, H_OLD_AVE, H_OLD_MAX, H_OLD_MID, H_OLD_MIN
  REAL(DP) :: H_NEW, H_NEW_AVE, H_NEW_MAX, H_NEW_MID, H_NEW_MIN
  REAL(DP) :: H_LIMIT, REFINE, WEIGHT
  REAL(DP) :: R_MAX, R_MIN, R_MID, R_AVE
  REAL(DP) :: AVE_VOL, EL_VOL
  
! AVE_H_WT = (KEYWORD) WEIGHT NEW SIZES. 1.0-REDUCE AVERAGE AND
!            LARGER SIZES, 0.0-REDUCE ALL ELEMENT SIZES
! H_LIMIT  = SMALLEST NEW SIZE, = H_OLD_AVE * AVE_TO_LIMIT
! H_SCALE  = ELEMENT SCALE, NEW SIZE = OLD SIZE * H_SCALE
! LT_SHAP  = ELEMENT TYPE SHAPE FLAG NUMBER
! L_SHAPE  = ELEM SHAPE, 1=LINE 2=TRI 3=QUAD 4=HEX 5=TET 6=WEDGE
! NODE_AVE_EL_SIZE = NEW ELEMENT SIZES AVERAGED AT ALL NODES
! REFINE   = FINAL SCALED, WEIGHTED REFINEMENT VALUE

!   SAVE NEW MESH SCALING AND DENSITY RESULTS
  U_NEW_H = GET_NEXT_IO_UNIT ()
  OPEN (U_NEW_H, FILE='msh_new_el_size.tmp', ACTION='WRITE', &
    STATUS='REPLACE', IOSTAT = IO_1)
  IF ( IO_1 > 0 ) THEN
    PRINT *, 'WARNING: OPEN FAILED FOR msh_new_el_size.tmp'
    N_WARN = N_WARN + 1 ; RETURN
  END IF
  U_NEW_S = GET_NEXT_IO_UNIT ()
  OPEN (U_NEW_S, FILE='msh_new_el_scale.tmp', ACTION='WRITE', &
    STATUS='REPLACE', IOSTAT = IO_1)
  IF ( IO_1 > 0 ) THEN
    PRINT *, 'WARNING: OPEN FAILED FOR msh_new_el_scale.tmp'
    N_WARN = N_WARN + 1 ; RETURN
  END IF

! INITIALIZE REFINEMENT VALUES. REFINE IF > 1
  NODE_AVE_EL_SIZE = 0.d0 ; NODE_AVE_COUNT = 0 ; H_SCALE = 1.d0
  LOW_R = 0.7  !b parameter xxx
  !b IF ( .NOT. DE_REFINE ) LOW_R = 1.0 !B add to global
  R_MAX = MAXVAL (ELEM_REFINEMENT)
  R_MIN = MINVAL (ELEM_REFINEMENT)
  R_MID = (R_MIN + R_MAX) * 0.5d0
  R_AVE = SUM (ELEM_REFINEMENT) / N_ELEMS
!b print *,'R_ave, R_max, R_mid, R_min '
!b print *, R_ave, R_max, R_mid, R_min
  IF ( R_AVE > 1.d0 ) THEN ! SCALE DOWN TO SLOW MESH GROWTH
    R_MAX = R_MAX / R_AVE
    R_MID = R_MID / R_AVE
    R_MIN = R_MIN / R_AVE
!b print *,'Scaled: R_max, R_mid, R_min '
!b print *, R_max, R_mid, R_min
  END IF ! Scaling

  AVE_VOL = VOLUME / N_ELEMS
  IF ( L_SHAPE == 2 .OR. L_SHAPE == 6 ) AVE_VOL = AVE_VOL * TRI_WEDGE
  IF ( L_SHAPE == 5 )                   AVE_VOL = AVE_VOL * TET      
  SELECT CASE (N_SPACE)          ! GET AVERAGE MESH SIZE
    CASE (1) ; H_OLD_AVE = AVE_VOL
    CASE (2) ; H_OLD_AVE = SQRT (AVE_VOL)
    CASE DEFAULT ; H_OLD_AVE = AVE_VOL **(1.d0 / N_SPACE)
  END SELECT ! for ave size
  H_OLD_MIN = HUGE (H_OLD_MIN) ; H_OLD_MAX = -HUGE (H_OLD_MAX) 
  H_NEW_MIN = HUGE (H_NEW_MIN) ; H_NEW_MAX = -HUGE (H_NEW_MAX) 
  H_LIMIT   = H_OLD_AVE * AVE_TO_LIMIT
  !b print *,'H_old_ave, H_LIMIT ', H_old_ave, H_LIMIT
  H_OLD_AVE = 0.d0 ; H_NEW_AVE = 0.d0

!     LOOP OVER ALL STANDARD ELEMENTS
  LT = 1 ; LAST_LT = 0
  DO IE = 1, N_ELEMS  ! LOOP OVER ALL ELEMENTS  -----------------------
    IF ( DELETED (IE) ) CYCLE ! TO ACTIVE ELEMENT
    CALL SET_THIS_ELEMENT_NUMBER (IE)

!-->     GET ELEMENT TYPE NUMBER
    IF ( N_L_TYPE > 1) LT = L_TYPE (IE)  ! SAME AS LAST TYPE ?
    IF ( LT /= LAST_LT ) THEN            ! this is a new type
      CALL GET_ELEM_TYPE_DATA (LT)       ! CONTROLS FOR THIS TYPE
      LAST_LT = LT
      IF ( TYPE_APLY_ALLOC ) CALL DEALLOCATE_TYPE_APPLICATION
      CALL ALLOCATE_TYPE_APPLICATION ! for ELEM_NODES (LT_N)
      EL_DEG = GET_SCALAR_DEGREE ()
    END IF ! a new element type
     
!     GET AVERAGE ELEMENT SIZE
    EL_VOL = MEASURE (IE)
    IF ( LT_SHAP == 2 .OR. LT_SHAP == 6 ) EL_VOL = EL_VOL * TRI_WEDGE
    IF ( LT_SHAP == 5 )                   EL_VOL = EL_VOL * TET      
    SELECT CASE (N_SPACE)
      CASE (1) ; H_OLD = EL_VOL
      CASE (2) ; H_OLD = SQRT (EL_VOL)
      CASE DEFAULT ; H_OLD = EL_VOL **(1.d0 / N_SPACE)
    END SELECT ! for old size
    H_OLD_AVE = H_OLD_AVE + H_OLD
    IF ( H_OLD > H_OLD_MAX ) H_OLD_MAX = H_OLD
    IF ( H_OLD < H_OLD_MIN ) H_OLD_MIN = H_OLD
     
!    GET NEW ELEMENT SIZE, WEIGHT TO CONTROL GROWTH
    WEIGHT = R_MIN * (1.d0 - AVE_H_WT) + R_AVE * AVE_H_WT 
    REFINE = ELEM_REFINEMENT (IE) / WEIGHT 
    IF ( REFINE > LOW_R ) THEN
      SELECT CASE (EL_DEG)
        CASE (1) ; H_SCALE (IE) = 1.d0 / REFINE
        CASE (2) ; H_SCALE (IE) = 1.d0 / SQRT (REFINE)
        CASE DEFAULT ; H_SCALE (IE) = 1.d0 / REFINE **(1.d0 / EL_DEG)
      END SELECT
    END IF ! reduce size (no increases)
    H_NEW = H_OLD * H_SCALE (IE)
    IF ( H_NEW < H_LIMIT ) H_NEW = H_LIMIT
    IF ( H_NEW > H_NEW_MAX ) H_NEW_MAX = H_NEW
    IF ( H_NEW < H_NEW_MIN ) H_NEW_MIN = H_NEW
    H_NEW_AVE = H_NEW_AVE + H_NEW

!     SEND SIZES FOR AVERAGING AT NODES
if ( IE < 1 ) stop 'impossible el 3'
    ELEM_NODES = GET_ELEM_NODES (IE, LT_N, NODES)
    DO J = 1, LT_N
      K = ELEM_NODES (J)
      IF ( K < 1 ) CYCLE ! to a valid node
        NODE_AVE_EL_SIZE (K) = NODE_AVE_EL_SIZE (K) + H_NEW
        NODE_AVE_COUNT   (K) = NODE_AVE_COUNT   (K) + 1     
    END DO ! over element nodes
  END DO ! over elements  ---------------------------------------------
  H_OLD_AVE = H_OLD_AVE / N_ELEMS
  H_NEW_AVE = H_NEW_AVE / N_ELEMS
  H_NEW_MID = (H_NEW_MIN + H_NEW_MAX) * 0.5d0
  H_OLD_MID = (H_OLD_MIN + H_OLD_MAX) * 0.5d0
  !b print *,'H_old_ave, H_old_max, H_old_mid, H_old_min '
  !b print *, H_old_ave, H_old_max, H_old_mid, H_old_min
  !b print *,'H_new_ave, H_new_max, H_new_mid, H_new_min '
  !b print *, H_new_ave, H_new_max, H_new_mid, H_new_min

!  AVERAGE SIZES AT NODES FOR MESH GENERATOR
  WHERE ( NODE_AVE_COUNT > 0 )
    NODE_AVE_EL_SIZE = NODE_AVE_EL_SIZE / NODE_AVE_COUNT
  ELSEWHERE
    NODE_AVE_EL_SIZE = H_NEW_AVE
  ENDWHERE

!    SAVE RESULTS FOR PLOTS AND/OR MESH GENERATOR
  WRITE (U_NEW_H, '(5I12)') MAX_NP, (K, K = 1, N_SPACE) ! for matlab
  DO J = 1, MAX_NP  ! save averaged sizes at nodes
    WRITE (U_NEW_H, '(3(1PE12.4))') X (J, :), NODE_AVE_EL_SIZE (J)
  END DO ! over nodes
  IF (IO_1 == 0) PRINT *,'NOTE: NEW SIZES SAVED TO msh_new_el_size.tmp'
  DO J = 1, N_ELEMS  ! save scaling of old elements
    WRITE (U_NEW_S, '(1PE12.4)') H_SCALE (J)
  END DO ! over elements
  IF (IO_1 == 0) PRINT *,'NOTE: EL SCALES SAVED TO msh_new_el_scale.tmp'
  CLOSE (U_NEW_H) ; CLOSE (U_NEW_S)
END SUBROUTINE FORM_NEW_EL_SIZES 

SUBROUTINE GATHER_LT_SCP_FLUX_2ND (SCP_AVERAGES, DERIV2_LT, &
                                   DERIV2_AVE, FLUX_LT)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!              Gather flux values to an element
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants ! for MAX_NP, NEIGH_L, N_ELEMS, N_QP,
                     ! SCP_FIT, U_SCPR, SKIP_ERROR
Use Sys_Properties_Data  ! for DELETED (L_S_TOT)
Use Elem_Type_Data  ! for:
  ! LT_FREE, LT_GEOM, LT_N, LT_PARM, LT_QP, ELEM_NODES (LT_N),&
  ! COORD (LT_N, N_SPACE), GEOMETRY (LT_GEOM, N_SPACE),       &
  ! C (LT_FREE), D (LT_FREE), S (LT_FREE, LT_FREE),           &
  ! DLG (LT_PARM, LT_GEOM), DLG_QP (LT_PARM, LT_GEOM, LT_QP)  &
  ! DLH (LT_PARM, LT_N),    DLH_QP (LT_PARM, LT_N,    LT_QP), &
  ! DLV (LT_PARM, LT_FREE), DLV_QP (LT_PARM, LT_FREE, LT_QP), &
  ! G (LT_GEOM), G_QP (LT_GEOM, LT_QP),                       &
  ! H (LT_N),    H_QP (LT_N,    LT_QP),                       &
  ! V (LT_FREE), V_QP (LT_FREE, LT_QP),                       &
  ! PT (LT_PARM, LT_QP), WT (LT_QP), D2LH (N_2_DER, LT_N)
Use SCP_Type_Data
Use Interface_Header
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: SCP_AVERAGES (MAX_NP,  SCP_FIT)     ! GLOBAL
 REAL(DP), INTENT(OUT) :: DERIV2_LT (N_R_B*N_SPACE, NOD_PER_EL) ! LONG
 REAL(DP), INTENT(OUT) :: DERIV2_AVE(N_2_DER, NOD_PER_EL) ! SCP AVERAGES
 REAL(DP), INTENT(OUT) :: FLUX_LT   (N_R_B,   NOD_PER_EL) ! NODAL FLUXES

 INTEGER :: HIGHER  ! N_R_B * N_SPACE
 INTEGER :: IN

!                   VARIABLES:
! DELETED    = TRUE IF THE ELEMENT HAS BEEN DELETED
! DERIV2_AVE = CALCULUS SECOND DERIVATIVES (CROSS AVERAGED)
! DERIV2_LT  = CALCULUS SECOND DERIVATIVES (NOT CROSS AVERAGED)
!              D2 u/Dx Dy /= D2 u/Dy Dx, currently
! ELEM_NODES = THE NOD_PER_EL INCIDENCES OF THE ELEMENT
! FLUX_LT    = FLUX AT ELEMENT NODES INTERPOLATED FROM A SCP
! HIGHER     = NUMBER OF HIGHER DEGREE DERIVATIVE FLUX TERMS
! NOD_PER_EL = MAXIMUM NUMBER OF NODES PER ELEMENT
! N_R_B      = NUMBER OF ROWS IN B AND E MATRICES
! N_SPACE    = DIMENSION OF SPACE
! N_2_DER    = NUMBER OF CALCULUS SECOND DERIVATIVES = 1, 3, or 6
! SCP_AVERAGES = AVERAGED FLUXES AT ALL NODES IN MESH
! SCP_FIT    = NUMBER IF TERMS BEING AVERAGED, = N_R_B MIN

!-->    EXTRACT SCP NODAL FLUXES
  HIGHER            = N_R_B * N_SPACE                    ! CONSTANT
  DO IN = 1, LT_N  ! OVER NODES OF ELEMENT
    IF ( ELEM_NODES (IN) < 1 ) CYCLE ! TO VALID NODE
      FLUX_LT (1:N_R_B, IN) = SCP_AVERAGES (ELEM_NODES (IN), 1:N_R_B)
      IF ( SCP_2ND_DERIV ) THEN  ! (SECOND DERIVATIVE INFO)
!         GET NODAL AVERAGE FLUX GRADIENTS 
        DERIV2_LT (1:HIGHER, IN) = SCP_AVERAGES (ELEM_NODES (IN), &
                                       (N_R_B+1):(N_R_B+HIGHER)) 
!           AVERAGE CROSS TERMS 
        SELECT CASE (N_SPACE) 
          CASE (1) ; DERIV2_AVE = DERIV2_LT 
          CASE (2) 
            !          1    2         3
            ! EXACT 3: U,xx U,yx=U,xy U,yy
            ! FE    4: U,xx U,yx U,xy U,yy 
            !          1    2    3    4
            DERIV2_AVE (1, IN) = DERIV2_LT (1, IN) 
            DERIV2_AVE (2, IN) = SUM (DERIV2_LT (2:3, IN)) * 0.5d0
            DERIV2_AVE (3, IN) = DERIV2_LT (4, IN) 
          CASE (3)
            !          1    2         3         4    5         6
            ! EXACT 6: U,xx U,yx=U,xy U,zx=U,xz U,yy U,zy=U,yz U,zz
            ! FE    9: U,xx U,yx U,zx U,xy U,yy U,zy U,xz U,yz U,zz
            !          1    2    3    4    5    6    7    8    9
            DERIV2_AVE (1, IN) =  DERIV2_LT (1, IN)
            DERIV2_AVE (2, IN) = (DERIV2_LT (2, IN)        &
                               +  DERIV2_LT (4, IN)) * 0.5d0
            DERIV2_AVE (3, IN) = (DERIV2_LT (3, IN)        & 
                               +  DERIV2_LT (7, IN)) * 0.5d0
            DERIV2_AVE (4, IN) =  DERIV2_LT (5, IN)
            DERIV2_AVE (5, IN) = (DERIV2_LT (6, IN)        &
                               +  DERIV2_LT (8, IN)) * 0.5d0
            DERIV2_AVE (6, IN) =  DERIV2_LT (9, IN)
          CASE DEFAULT 
            STOP 'GATHER_LT_SCP_FLUX_2ND, DERIV2_AVE UNDEFINED' 
        END SELECT  ! SPACE
      END IF  ! GET 2ND DERIVATIVE DATA
  END DO ! FOR NODES ON ELEMENT
END SUBROUTINE GATHER_LT_SCP_FLUX_2ND

SUBROUTINE GATHER_LT_SCP_FLUXES (SCP_AVERAGES, FLUX_LT)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!              Gather flux values to an element
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants ! for MAX_NP, NEIGH_L, N_ELEMS, N_QP,
                     ! SCP_FIT, U_SCPR, SKIP_ERROR
Use Sys_Properties_Data  ! for DELETED (L_S_TOT)
Use Elem_Type_Data  ! for:
  ! LT_FREE, LT_GEOM, LT_N, LT_PARM, LT_QP, ELEM_NODES (LT_N),&
  ! COORD (LT_N, N_SPACE), GEOMETRY (LT_GEOM, N_SPACE),       &
  ! C (LT_FREE), D (LT_FREE), S (LT_FREE, LT_FREE),           &
  ! DLG (LT_PARM, LT_GEOM), DLG_QP (LT_PARM, LT_GEOM, LT_QP)  &
  ! DLH (LT_PARM, LT_N),    DLH_QP (LT_PARM, LT_N,    LT_QP), &
  ! DLV (LT_PARM, LT_FREE), DLV_QP (LT_PARM, LT_FREE, LT_QP), &
  ! G (LT_GEOM), G_QP (LT_GEOM, LT_QP),                       &
  ! H (LT_N),    H_QP (LT_N,    LT_QP),                       &
  ! V (LT_FREE), V_QP (LT_FREE, LT_QP),                       &
  ! PT (LT_PARM, LT_QP), WT (LT_QP), D2LH (N_2_DER, LT_N)
Use SCP_Type_Data
!b Use Interface_Header
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: SCP_AVERAGES (MAX_NP,  SCP_FIT)     ! GLOBAL
 REAL(DP), INTENT(OUT) :: FLUX_LT   (N_R_B,   NOD_PER_EL) ! NODAL FLUXES

 INTEGER :: IN  ! loops

! DELETED    = TRUE IF THE ELEMENT HAS BEEN DELETED
! ELEM_NODES = THE NOD_PER_EL INCIDENCES OF THE ELEMENT
! FLUX_LT    = FLUX AT ELEMENT NODES INTERPOLATED FROM A SCP
! MAX_NP     = NUMBER OF SYSTEM NODES
! NOD_PER_EL = MAXIMUM NUMBER OF NODES PER ELEMENT
! N_R_B      = NUMBER OF ROWS IN B AND E MATRICES
! SCP_AVERAGES = AVERAGED FLUXES AT ALL NODES IN MESH
! SCP_FIT    = NUMBER IF TERMS BEING AVERAGED, = N_R_B MIN

!-->    EXTRACT SCP NODAL FLUXES
  DO IN = 1, LT_N  ! OVER NODES OF ELEMENT
    IF ( ELEM_NODES (IN) < 1 ) CYCLE ! TO VALID NODE
      FLUX_LT (1:N_R_B, IN) = SCP_AVERAGES (ELEM_NODES (IN), 1:N_R_B)
  END DO ! FOR NODES ON ELEMENT
END SUBROUTINE GATHER_LT_SCP_FLUXES

SUBROUTINE SURF_SCP_AVE_NODE_FLUXES  (NODES, X, L_NEIGH, &
                                      DD, SCP_AVERAGES)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
!  CALCULATE THE SUPER_CONVERGENCE_PATCH AVERAGE FLUXES   
!  AT ALL NODES ON A PARAMETRIC SURFACE MESH (OR PLANE)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
! (USES A LEAST SQUARES FIT OF THE DISPLACEMENTS THROUGH THE 
! MESH NODES. IT SOLVES THAT LOCAL PROBLEM BY THE SINGULAR VALUE 
! DECOMPOSITION SOLUTION METHOD AND THEN RECOVERS LOCAL FLUXES)
  Use System_Constants ! for MAX_NP, NEIGH_P, N_ELEMS, N_PATCH, N_QP, 
                       ! SCP_FIT, U_SCPR, ON_BOUNDARY
  Use Sys_Properties_Data  ! for DELETED
  Use Elem_Type_Data       ! for LAST_LT, LT_*
  Use SCP_Type_Data        ! for SCP_H (SCP_N), SCP_DLH (SCP_PARM,SCP_N)
  Use Interface_Header
  Use Select_Source
   IMPLICIT NONE
   INTEGER,  INTENT (IN)  :: NODES        (L_S_TOT, NOD_PER_EL)
   INTEGER,  INTENT (IN)  :: L_NEIGH      (NEIGH_P, N_PATCH)  
   REAL(DP), INTENT (IN)  :: X            (MAX_NP, N_SPACE)
   REAL(DP), INTENT (IN)  :: DD           (N_D_FRE)        
   REAL(DP), INTENT (OUT) :: SCP_AVERAGES (MAX_NP, SCP_FIT)

   INTEGER :: NODES_PER   (N_PATCH)       ! NODES PER PATCH
   INTEGER :: MEMBERS     (NEIGH_P + 1)   ! ELEMENTS IN PATCH  
   INTEGER :: FIRST       (NOD_PER_EL)    ! NODES OF FIRST PATCH ELEMENTS !b
   INTEGER :: SCP_COUNTS  (MAX_NP)        ! PATCH HITS PER NODE

   INTEGER,  ALLOCATABLE :: NP_SURF (:)       ! SYS NODES ON SURFACE PATCH

   REAL(DP), ALLOCATABLE :: PATCH_SQ  (:, :)  ! SCRATCH ARRAY
   REAL(DP), ALLOCATABLE :: PATCH_DAT (:, :)  ! SCRATCH ARRAY
   REAL(DP), ALLOCATABLE :: PATCH_P   (:, :)  ! SCRATCH ARRAY
   REAL(DP), ALLOCATABLE :: PATCH_WRK (:)     ! SCRATCH ARRAY
   REAL(DP), ALLOCATABLE :: PATCH_FIT (:, :)  ! ANSWERS

   REAL(DP), ALLOCATABLE :: SURF_XY (:, :)  ! NODES IN PARMETRIC COORDS !b
   REAL(DP), ALLOCATABLE :: SCP_DGH (:, :)  ! GLOBAL DERIVATIVES OF SCP_H !b

   REAL(DP) :: EL_MIN  (N_SPACE), EL_MAX  (N_SPACE)  ! BOUNDS !b
   REAL(DP) :: XYZ_MIN (N_SPACE), XYZ_MAX (N_SPACE)  ! BOUNDS
   REAL(DP) :: XYZ     (N_SPACE), FLUX    (N_R_B)    ! POINT & FLUX 
   REAL(DP) :: BORDER  (N_SPACE)                     ! MAPPING SAFETY
   REAL(DP) :: POINT   (N_PARM)                      ! SCP POINT !b
   REAL(DP) :: GRAD_PT (N_R_B)                       ! PT GRAD !b
   REAL(DP) :: FLUX_PT (N_R_B)                       ! PT FLUX !b
   REAL(DP) :: E       (N_R_B, N_R_B)                ! CONSTITUTIVE !b
   REAL(DP) :: DIAGONAL_JACOBIAN (N_SPACE) ! XXX flat only
   REAL(DP) :: GRAD_DOF (2, N_G_DOF) ! GLOBAL DERIVATIVES OF DOF !b

   INTEGER  :: NODES_IN_IP           ! NUMBER OF MESH NODES IN PATCH  !b
   INTEGER  :: FIRST_EL, FIRST_N     ! FIRST EL NUMBER, NODES COUNT !b
   INTEGER  :: POINTS                ! SIZES OF CURRENT PATCH MATRICES
   INTEGER  :: L_IN_PATCH            ! NUMBER OF ELEMS IN CURRENT PATCH
   INTEGER  :: FIT, IL, IP, IQ, JP, IS ! LOOPS !b
   INTEGER  :: ROW                   ! IN PATCH ARRAYS
   INTEGER  :: LM, LT, REC_LM_IQ, SCP_STAT, SYS, unit
   INTEGER  :: COUNT_4 = 0, COUNT_9 = 0, COUNT_16 = 0, COUNT_17 = 0
   INTEGER, PARAMETER :: ONE = 1

!b interface
  !b function get_next_io_unit () result (next)
    !b integer :: next
  !b end function get_next_io_unit
!b end interface

!                   VARIABLES:
! DD         = SYSTEM LIST OF NODAL PARAMETERS
! DELETED    = TRUE IF THE ELEMENT HAS BEEN DELETED
! E          = CONSTITUTIVE MATRIX
! ELEM_NODES = THE NOD_PER_EL INCIDENCES OF THE ELEMENT
! FIRST      = LIST OF NODES ON FIRST ELEMENT IN A PATCH
! FIRST_EL   = FIRST ELEMENT NUMBER IN A PATCH
! FIRST_N    = NUMBER OF NODES ON FIRST ELEMENT IN A PATCH
! INDEX      = SYSTEM DOF NUMBERS ASSOCIATED WITH ELEMENT
! L_HOMO     = 1, IF ELEMENT PROPERTIES ARE HOMOGENEOUS
! L_SHAPE    = SHAPE FLAG FOR QUADRATURE RULE SELECTION
! L_S_TOT    = TOTAL NUMBER OF ELEMENTS & THEIR SEGMENTS
! L_NEIGH    = ELEMENTS FORMING THE PATCH    !b
! L_TYPE     = ELEMENT TYPE NUMBER
! LT         = ELEMENT TYPE NUMBER (IF USED)
! LT_DATA    = ARRAY OF CONSTANTS FOR EACH ELEMENT TYPE
! LT_QP      = NUMBER OF QUADRATURE PTS FOR ELEMENT TYPE
! LT_USER    = OPTIONAL USER CODE ASSOCIATED WITH TYPE
! MAX_NP     = NUMBER OF SYSTEM NODES
! MEMBERS    = ELEMENT NUMBERS MAKING UP A SCP 
! NEIGH_P    = MAXIMUM NUMBER OF ELEMENTS IN A PATCH
! NOD_PER_EL = NUMBER OF NODES PER ELEMENT
! NODES      = NODAL INCIDENCES OF ALL ELEMENTS
! N_ELEMS    = NUMBER OF ELEMENTS IN SYSTEM
! N_PATCH    = NUMBER OF PATCHES, MAX_NP OR N_ELEMS
! N_QP       = MAXIMUN NUMBER OF QUADRATURE POINTS, >= LT_QP
! N_R_B      = NUMBER OF ROWS IN B AND E MATRICES
! N_SPACE    = DIMENSION OF SPACE
! PATCH_FIT  = LOCAL PATCH VALUES FOR FLUX AT ITS NODES
! POINT      = LOCAL POINT IN PATCH INTERPOLATION SPACE
! POINTS     = TOTAL NUMBER OF QUADRATURE POINTS IN CURRENT PATCH
! SCP_AVERAGES = AVERAGED FLUXES AT ALL NODES IN MESH
! SCP_FIT    = NUMBER IF TERMS BEING FIT, OR AVERAGED: = N_R_B MIN    
! SCP_INC    = INCREASE PATCH POLYNOMIAL DEGREE BY THIS (1 or 2)
! SYS        = SYSTEM NODE NUMBER OF A NODE IN THE PATCH
! U_SCPR     = UNIT FOR SUPER_CONVERGENT PATCH RECOVERY, IF > 0
! X          = COORDINATES OF SYSTEM NODES
! XYZ        = SPACE COORDINATES AT A POINT
! XYZ_MAX    = UPPER BOUNDS FOR SCP GEOMETRY
! XYZ_MIN    = LOWER BOUNDS FOR SCP GEOMETRY

    LAST_LT = 0 ; SCP_COUNTS = 0 ; SCP_AVERAGES = 0.d0  ! INITIALIZE
    IF ( N_ELEMS == 1 ) THEN
      PRINT *, 'Single element patch error, SURF_SCP_AVE_NODE_FLUXES'
      RETURN ! assume expert user
    END IF ! SINGLE ELEMENT MESH

    LT = 1                           ! INITIALIZE
    IF ( N_PATCH == 0 ) THEN ! No patch type selected
      PRINT *,'WARNING: SURF_SCP_AVE_NODE_FLUXES, NO PATCH TYPE GIVEN'
      PRINT *,'PLACE KEYWORD OPTION scp_neigh_el OR scp_neigh_face OR'
      PRINT *,'scp_neigh_pt IN THE keyword CONTROLS.'
      PRINT *,'Skipping SCP averaging and error estimates.'
      N_WARN = N_WARN + 1 ; RETURN
    END IF
    IF ( SCP_INC > 0 ) PRINT *, 'NOTE: PATCH DEGREE INCREASED BY ', &
                                 SCP_INC

    DO IP = 1, N_PATCH ! LOOP OVER EACH PATCH ====>> ====>> ====>> ====>>
      MEMBERS = 0 ; POINTS = 0 ; FIRST = -1      ! INITIALIZE !b

!       ELEMENT OR NODAL CENTERED PATCH TYPE
      IF ( .NOT. SCP_NEIGH_PT ) THEN ! ELEMENT BASED ELEMENTS
        IF ( IP == 1 ) PRINT *,'NOTE: ELEMENT BASED SURFACE PATCHES'
        IF ( DELETED (IP) ) CYCLE ! TO AN ACTIVE ELEMENT ====>> ====>> 

!        GET ELEMENT NEIGHBORS TO DEFINE THE PATCH
        MEMBERS = (/ IP, L_NEIGH (:, IP) /)      ! TO SIMPLIFY LOOP 
 
      ELSE ! NODAL BASED PATCH OF ELEMENTS
        IF ( IP == 1 ) PRINT *,'NOTE: NODAL BASED SURFACE PATCHES'

!        GET ELEMENT NEIGHBORS TO DEFINE THE PATCH
        !b MEMBERS    = L_NEIGH (:, IP)       ! bug fix
        MEMBERS (1:NEIGH_P) = L_NEIGH (:, IP) ! bug fix
      END IF ! PATCH BASIS

      L_IN_PATCH = COUNT ( MEMBERS > 0 )
      IF ( L_IN_PATCH <= 1 ) CYCLE ! TO AN ACTIVE PATCH ====>> ====>> ====>>

      IF ( DEBUG_SCP ) THEN
        PRINT *, 'NEIGH_P ', NEIGH_P
        PRINT *, L_NEIGH (1:NEIGH_P, IP)
        PRINT *, 'IP, L_IN_PATCH ', IP, L_IN_PATCH
        PRINT *, 'SIZE(MEMBERS) ', SIZE(MEMBERS)
        PRINT *, 'MEMBERS ', MEMBERS
      END IF

!   COUNT MESH NODES IN PATCH !b
      LAST_LT = 0                           ! INITIALIZE, AGAIN !b
      DO IL = 1, L_IN_PATCH ! PATCH ELEMENT MEMBER LOOP ====> ====> ====> !b
        LM = MEMBERS (IL)                   ! ELEMENT IN PATCH !b
if ( LM < 1 ) then
  print *, 'impossible el 4'
  print *,'IP, IL, L_IN_PATCH, LM', IP, IL, L_IN_PATCH, LM
  print *,'members ', members (1:L_IN_PATCH)
end if

        IF ( N_L_TYPE > 1) LT = L_TYPE (LM) ! GET ELEMENT TYPE NUMBER !b
        IF ( LT /= LAST_LT ) THEN           ! NEW ELEMENT TYPE   !b
          CALL GET_ELEM_TYPE_DATA (LT)      ! GET CONTROLS FOR ELEM TYPE !b
          IF ( TYPE_APLY_ALLOC ) CALL DEALLOCATE_TYPE_APPLICATION !b
          CALL ALLOCATE_TYPE_APPLICATION    ! for ELEM_NODES (LT_N) !b
          LAST_LT = LT !b
        END IF ! a new element type and scp type !b

if ( LM < 1 ) stop 'impossible el 4'
        ELEM_NODES (1:LT_N) = GET_ELEM_NODES (LM, LT_N, NODES)     !b

!   UPDATE PHYSICAL SPACE BOUNDS
        CALL ELEM_COORD (LT_N, N_SPACE, X, COORD, ELEM_NODES)
        EL_MAX (:) = MAXVAL ( COORD (1:LT_N, :), DIM = 1 )
        EL_MIN (:) = MINVAL ( COORD (1:LT_N, :), DIM = 1 )
        IF ( IL == 1 ) THEN !b
          XYZ_MAX = EL_MAX ; XYZ_MIN = EL_MIN !b
        END IF !b
        WHERE ( EL_MAX > XYZ_MAX ) XYZ_MAX = EL_MAX
        WHERE ( EL_MIN < XYZ_MIN ) XYZ_MIN = EL_MIN

        IF ( IL == 1 ) THEN ! first element in patch !b
          FIRST (1:LT_N) = ELEM_NODES  ! first connectivity list !b
          FIRST_EL       = LM          ! first element number !b
          FIRST_N        = LT_N        ! number of nodes on first el !b
          NODES_IN_IP    = LT_N        ! start the count !b
          CYCLE                        ! to next IL in patch ====> ====> !b
        END IF ! first el !b
        !b INDEX = GET_ELEM_INDEX (LT_N, ELEM_NODES) !b

        ! count number of nodes occuring since first element !b
        DO JP = 1, LT_N  !b
          SYS = ELEM_NODES (JP)
          IF ( IL > 1 .AND. ANY (FIRST == SYS) ) CYCLE ! to JP ---> ---> 
          !b IF ( ANY ( FIRST /= ELEM_NODES (JP) ) ) &  !b
             NODES_IN_IP = NODES_IN_IP + 1  !b
        END DO ! over JP nodes of this element  !b
      END DO ! over IL elements in this patch <==== <==== <==== <====  !b
      NODES_PER (IP) = NODES_IN_IP  ! for planning
      POINTS         = NODES_IN_IP  ! redundant !b

!     ADD A SMALL BORDER TO PATCH
      BORDER  = ABS ( XYZ_MAX - XYZ_MIN ) * 0.01
      XYZ_MAX = XYZ_MAX + BORDER
      XYZ_MIN = XYZ_MIN - BORDER

!   IS THIS A VALID RECTANGULAR PATCH !b
      IF ( NODES_IN_IP < 4 ) THEN !b
        PRINT *,'WARNING: SKIPPING SINGULAR PATCH NUMBER ', IP !b
        N_WARN = N_WARN + 1 !b
        CYCLE ! to next IP patch ====>> ====>> ====>> ====>> ====>> ====>> !b
      END IF ! singular patch !b
                                 
!      VALIDATE PATCH SHAPE  !b
      IF ( ANY( (XYZ_MAX - XYZ_MIN) == 0.D0 )) THEN  !b
        PRINT *,'WARNING, PATCH ', IP, ' HAS A ZERO DIMENSION'  !b
        PRINT *,'MAX BOUNDS ARE: ', XYZ_MAX  !b
        PRINT *,'MIN BOUNDS ARE: ', XYZ_MIN  !b
        N_WARN = N_WARN + 1  !b
      END IF ! FINITE SHAPED PATCH  !b

!   SELECT SCP_N VALUE !b
      IF ( SCP_NEIGH_PT ) THEN ! point centered
        SELECT CASE (NODES_IN_IP) !b
          CASE ( 4:8)  ; SCP_N = 4  ! Use Q4 patch !b
                         COUNT_4 = COUNT_4 + 1
          CASE ( 9:15) ; SCP_N = 9  ! Use Q9 patch !b
                         COUNT_9 = COUNT_9 + 1
          CASE (16:30) ; SCP_N = 16 ! Use Q16 patch !b
                         COUNT_16 = COUNT_16 + 1
          CASE (31:)   ; SCP_N = 17 ! Use Q17 patch !b
                         COUNT_17 = COUNT_17 + 1
          ! also change summary report at the if case ranges changed
          CASE DEFAULT 
            PRINT *, 'PATCH NUMBER ', IP, ' TYPE UNDEFINED FOR'
            PRINT *, 'THE NUMBER OF MESH NODES = ', NODES_IN_IP
            PRINT *, 'ERROR SURF_SCP_AVE_NODE_FLUXES STOPPING'
            STOP 'SELECT SCP_N VALUE FAILED' !b
            !XXX also need scp qp_n here
        END SELECT !b
      ELSE ! ELEMENT CENTER
        SELECT CASE (NODES_IN_IP) !b
          CASE ( 4:10) ; SCP_N = 4  ! Use Q4 patch !b
                         COUNT_4 = COUNT_4 + 1
          CASE (11:20) ; SCP_N = 9  ! Use Q9 patch !b
                         COUNT_9 = COUNT_9 + 1
          CASE (21:30) ; SCP_N = 16 ! Use Q16 patch !b
                         COUNT_16 = COUNT_16 + 1
          CASE (31:)   ; SCP_N = 17 ! Use Q17 patch !b
                         COUNT_17 = COUNT_17 + 1
          ! also change summary report at the if case ranges changed
          CASE DEFAULT 
            PRINT *, 'PATCH NUMBER ', IP, ' TYPE UNDEFINED FOR'
            PRINT *, 'THE NUMBER OF MESH NODES = ', NODES_IN_IP
            PRINT *, 'ERROR SURF_SCP_AVE_NODE_FLUXES STOPPING'
            STOP 'SELECT SCP_N VALUE FAILED' !b
            !XXX also need scp qp_n here
        END SELECT !b
      END IF ! patch type
      SCP_SHAP = 3 ; SCP_PARM = 2 ; SCP_LT = 1 !b
      SCP_GEOM = SCP_N ! or analytic !b
      IF ( DEBUG_SCP ) THEN
        PRINT *,'SCP_GEOM, SCP_LT, SCP_N, SCP_PARM, SCP_SHAP', &
                 SCP_GEOM, SCP_LT, SCP_N, SCP_PARM, SCP_SHAP
      END IF

!XXX flat surface only
      CALL PATCH_CONSTANT_JACOBIAN (XYZ_MIN, XYZ_MAX, DIAGONAL_JACOBIAN)!b

      IF ( PATCH_ALLOC_STATUS ) THEN
!         DEALLOCATE LOCAL PATCH RELATED ARRAYS
        DEALLOCATE ( NP_SURF   )
        DEALLOCATE ( SCP_DGH   )
        DEALLOCATE ( SURF_XY   )
        DEALLOCATE ( PATCH_WRK )
        DEALLOCATE ( PATCH_SQ  )
        DEALLOCATE ( PATCH_FIT )
        DEALLOCATE ( PATCH_DAT )
        DEALLOCATE ( PATCH_P   )
        PATCH_ALLOC_STATUS = .FALSE.
      END IF ! STATUS CHECK

!       ALLOCATE LOCAL PATCH RELATED ARRAYS
      ALLOCATE ( PATCH_P   (POINTS, SCP_N  ) )
      !b ALLOCATE ( PATCH_DAT (POINTS, SCP_FIT) )
      ALLOCATE ( PATCH_DAT (POINTS, N_G_DOF) )  !b
      !b ALLOCATE ( PATCH_FIT (SCP_N , N_R_B  ) ) 
      ALLOCATE ( PATCH_FIT (SCP_N , N_G_DOF) )   !b
      ALLOCATE ( PATCH_SQ  (SCP_N , SCP_N  ) )
      ALLOCATE ( PATCH_WRK (SCP_N )          )
      ALLOCATE ( SURF_XY   (2     , NODES_IN_IP) )
      ALLOCATE ( SCP_DGH   (2     , SCP_N) )
      ALLOCATE ( NP_SURF   (NODES_IN_IP) )
      PATCH_ALLOC_STATUS = .TRUE. 

!      ZERO PATCH WORKSPACE AND RESULTS ARRYS
      PATCH_P   = 0.d0 ; PATCH_DAT = 0.d0 ; PATCH_FIT = 0.d0     
      PATCH_SQ  = 0.d0 ; PATCH_WRK = 0.d0             

      IF ( DEBUG_EL_TYPE .OR. DEBUG_SCP ) THEN
        PRINT *,'DEBUG_ Elem type data ' 
        CALL IPRINT (LT_DATA, MAX_DAT,  N_L_TYPE) 
        DO IL = 1, N_L_TYPE 
          CALL GET_ELEM_TYPE_DATA  (IL) 
          CALL LIST_ELEM_TYPE_DATA (IL) 
        END DO 
      END IF ! Debug

!      PREPARE LEAST SQUARES FIT MATRICES FOR NODE DOF !b
      LAST_LT = 0                           ! INITIALIZE, AGAIN
      ROW = 0                               ! INITIALIZE
!b print *,'FIRST ', FIRST(1:FIRST_N)
!b print *,'IP, NODES_IN_IP',IP, NODES_IN_IP
      DO IL = 1, L_IN_PATCH ! PATCH ELEMENT MEMBER LOOP ====> ====> ====>
        LM = MEMBERS (IL)                   ! ELEMENT IN PATCH
        IF ( N_L_TYPE > 1) LT = L_TYPE (LM) ! GET ELEMENT TYPE NUMBER
        IF ( LT /= LAST_LT ) THEN           ! NEW ELEMENT TYPE  
          CALL GET_ELEM_TYPE_DATA (LT)      ! GET CONTROLS FOR ELEM TYPE
          IF ( TYPE_APLY_ALLOC ) CALL DEALLOCATE_TYPE_APPLICATION
          CALL ALLOCATE_TYPE_APPLICATION     ! for ELEM_NODES (LT_N)
          LAST_LT = LT
        END IF ! a new element type and scp type
if ( LM < 1 ) stop 'impossible el 5'
        ELEM_NODES (1:LT_N) = GET_ELEM_NODES (LM, LT_N, NODES)     !b
!b print *,'IL, LM, LT_N, ELEM_NODES (1:LT_N)', &
         !b IL, LM, LT_N, ELEM_NODES (1:LT_N)

        IF ( .NOT. SCP_SCAL_ALLOC ) CALL ALLOCATE_SCP_INTERPOLATIONS
        DO JP = 1, LT_N        ! LOOP OVER NODE POINTS ---> ---> ---> --->
          SYS = ELEM_NODES (JP) ! GLOBAL NODE NUMBER
!b print *, 'JP, SYS, ANY_IN', JP, SYS, ANY (FIRST == SYS)
!b if ( sys < 1 .OR. sys > MAX_NP ) stop 'impossible node'
! verify a new point, get dof, parmetric position (saved)     
          IF ( IL > 1 .AND. ANY (FIRST == SYS) ) CYCLE ! to JP ---> ---> !b
            ROW = ROW + 1                    ! UPDATE ROW LOCATION
            NP_SURF (ROW) = SYS              ! DEBUG AID

!             CONVERT JP XYZ TO LOCAL PATCH SURFACE POINT
!b print *,'ROW, SYS, N_SPACE, X(1,1)', ROW, SYS, N_SPACE, X(1,1)
            XYZ (:) = X (SYS, :) !b
! XXX   following for planar only, convert to shallow shell
            POINT = GET_SCP_PT_AT_XYZ (XYZ, XYZ_MIN, XYZ_MAX)
            SURF_XY (1:2, ROW) = POINT (1:2)

!             EVALUATE PATCH INTERPOLATION AT LOCAL POINT
            !b IF ( .NOT. SCP_SCAL_ALLOC ) CALL ALLOCATE_SCP_INTERPOLATIONS
            !b CALL GEN_ELEM_SHAPE (POINT, SCP_H, SCP_N, N_SPACE, ONE)
            CALL GEN_ELEM_SHAPE (POINT, SCP_H, SCP_N, SCP_PARM, ONE)

!             INSERT FLUX & INTERPOLATIONS INTO PATCH MATRICES
            INDEX (1:N_G_DOF) = GET_INDEX_AT_PT (SYS)
            PATCH_DAT (ROW, 1:N_G_DOF) = DD (INDEX(1:N_G_DOF)) !b
            PATCH_P   (ROW, 1:SCP_N)   = SCP_H (1:SCP_N)
!b if ( IP == 1 ) THEN
 !b print *,'IP, ROW, SYS, PATCH_DAT (ROW, 1:N_G_DOF)'
 !b print *, IP, ROW, SYS, PATCH_DAT (ROW, 1:N_G_DOF)
!b END IF
        END DO ! FOR EACH JP NODAL DOF VALUES     <--- <--- <--- <--- <---
!b print *,'IL, NODES_IN_IP, ROW, POINTS', IL, NODES_IN_IP, ROW, POINTS

      END DO ! FOR EACH IL PATCH MEMBER      <==== <==== <==== <==== <====

!   ASSEMBLY OF PATCH COMPLETED
!b if (ROW /= POINTS) STOP 'ROW /= POINTS, as required' !b

!       VALIDATE CURRENT PATCH
      IF ( POINTS < SCP_N ) THEN 
        PRINT *,'WARNING: SURF_SCP_AVE_NODE_FLUXES SKIPPING PATCH ', &
                 IP, ' WITH ONLY ', POINTS, ' EQUATIONS'
        N_WARN = N_WARN + 1 ! INCREMENT WARNING 
        CYCLE ! TO NEXT PATCH    ====>> ====>> ====>> ====>> ====>> ====>>
      END IF ! INSUFFICIENT DATA 

!        USE SINGULAR VALUE DECOMPOSITION SOLUTION METHOD
      CALL SVDC_FACTOR (PATCH_P, POINTS, SCP_N, PATCH_WRK, PATCH_SQ)

      WHERE ( PATCH_WRK < EPSILON(1.d0) ) PATCH_WRK = 0.D0

!b    DO FIT = 1, N_R_B   ! LOOP FOR EACH FLUX COMPONENT TO FIT 
      DO FIT = 1, N_G_DOF ! LOOP FOR EACH NODAL DOF TO FIT 
        CALL SVDC_BACK_SUBST (PATCH_P, PATCH_WRK, PATCH_SQ,        &
                              POINTS, SCP_N, PATCH_DAT (:, FIT),   &
                              PATCH_FIT (:, FIT))
      END DO ! FOR DOF COMPONENTS
!b IF ( IP == 1 ) THEN
  !b print *, 'FIT DISPLACEMENTS IP=1'; call RPRINT (PATCH_FIT,SCP_N,N_G_DOF,1)
!b END IF

!  LOOP OVER NODES OF PATCH AND GET INTERPOLATED GRADIENTS
      XYZ = 0.d0 !b fake for E recovery
      CALL SELECT_APPLICATION_E_MATRIX (FIRST_EL, XYZ, E) ! constitutive law
!b if (  ip <= 2 ) then
  !b print *,'E matrix'; call rprint (E, N_R_B, N_R_B, 1)
!b end if

      DO JP = 1, NODES_IN_IP ! --> -->  -->  -->  -->  -->  -->  -->  -->
        IF ( SCP_NEIGH_PT .AND. SCP_CENTER_ONLY .AND. JP > 1 ) CYCLE !-->
        IF ( SCP_CENTER_ONLY .AND. JP > FIRST_N ) CYCLE !   -->  -->  -->
          POINT (1:2) = SURF_XY (1:2, JP) ! recover parametric point
          SYS         = NP_SURF (JP)      ! recover system node
          CALL GEN_ELEM_DERIV (POINT, SCP_DLH (1:SCP_PARM, 1:SCP_N), &
                               SCP_N, SCP_PARM, ONE)
IF ( IP<=2 .AND. JP<=LT_N) THEN  !b
! print node displacement as check  !b
  CALL GEN_ELEM_SHAPE (POINT, SCP_H, SCP_N, SCP_PARM, ONE)  !b
  print *,'node / dof ', NP_SURF (JP), MATMUL (SCP_H, PATCH_FIT)  !b
  INDEX (1:N_G_DOF) = GET_INDEX_AT_PT (NP_SURF (JP))  !b
  print *,'node / true', NP_SURF (JP), DD (INDEX(1:N_G_DOF)) !b
END IF  !b

!       MULTIPLY BY INVERSE OF PATCH'S DIAGONAL_JACOBIAN !XXX
!b print *,'IP, JP, size(SCP_DLH), size(SCP_DGH)', &
!b          IP, JP, size(SCP_DLH), size(SCP_DGH)
          DO IS = 1, N_SPACE !XXX
            SCP_DGH (IS, :) = SCP_DLH (IS, :) / DIAGONAL_JACOBIAN (IS) !XXX
          END DO ! GLOBAL DERIVATIVES !XXX

!       DOT WITH PATCH FITTED DOF FOR SURFACE GRADIENTS
          GRAD_DOF = MATMUL (SCP_DGH (1:SCP_PARM, 1:SCP_N), & !b
                             PATCH_FIT (1:SCP_N, 1:N_G_DOF)) !b
!       RESHAPE TO VECTOR FORM
          GRAD_PT = (/ GRAD_DOF (1,1),  GRAD_DOF (2,2),  &
                      (GRAD_DOF (2,1) + GRAD_DOF (1,2)) /)

!       OPTIONALLY (FOR SHALLOW SHELL) GET SURFACE DISPLACEMENT !b
          !b CALL GEN_ELEM_SHAPE (POINT, SCP_H, SCP_N, SCP_PARM, ONE) !b
          !b would need surf normal: w = normal.dof at node JP !b

!       CONVERT GRADIENTS TO STRAINS OR STRESSES (ASSUME E CONSTANT) !b
          FLUX_PT = MATMUL (E, GRAD_PT) !b

!b IF ( IP==2 .AND. JP<=4) THEN
  !b !b print *,'SCP_DLH 1 ', SUM ( SCP_DLH (1, :) )
  !b !b print *,'SCP_DLH 2 ', SUM ( SCP_DLH (2, :) )
  !b !b print *,'SCP_DLH '; call rprint (SCP_DLH,N_SPACE,SCP_N,1)
  !b !b print *,'SCP_DGH '; call rprint (SCP_DGH,N_SPACE,SCP_N,1)
  !b print *,'GRAD_DOF  '; call rprint (GRAD_DOF, SCP_PARM, N_G_DOF,1)
  !b print *,'GRAD_PT   ', GRAD_PT
  !b print *,'FLUX_PT   ', FLUX_PT
!b END IF 

!       SCATTER POINT VALUES TO SYSTEM NODES, INCREMENT COUNTS !b
          SCP_COUNTS   (SYS)          = SCP_COUNTS   (SYS) + 1
          SCP_AVERAGES (SYS, 1:N_R_B) = SCP_AVERAGES (SYS, 1:N_R_B) & !b
                                      + FLUX_PT      (1:N_R_B) !b
          !b NODE_DONE    (SYS)          = .TRUE.
      END DO ! for mesh nodes in patch  <--  <--  <--  <--  <--  <--  <--

      !b IF ( DEBUG_SCP .AND. IP == 1) PRINT *,'AVERAGING NODAL DOF'     

!        DEALLOCATE LOCAL PATCH RELATED ARRAYS
      IF ( SCP_SCAL_ALLOC ) CALL DEALLOCATE_SCP_INTERPOLATIONS
      DEALLOCATE ( NP_SURF   )
      DEALLOCATE ( SCP_DGH   )
      DEALLOCATE ( SURF_XY   )
      DEALLOCATE ( PATCH_WRK )
      DEALLOCATE ( PATCH_SQ  )
      DEALLOCATE ( PATCH_FIT )
      DEALLOCATE ( PATCH_DAT )
      DEALLOCATE ( PATCH_P   )
      PATCH_ALLOC_STATUS = .FALSE.
    END DO ! FOR EACH IP PATCH IN MESH <<==== <<==== <<==== <<==== <<====

    PRINT *,'' ; PRINT *,'PATCH SUMMARIES'
    PRINT *,'MINIMUM NODES PER PATCH = ', MINVAL (NODES_PER)
    PRINT *,'MAXIMUM NODES PER PATCH = ', MAXVAL (NODES_PER)
    PRINT *,'AVERAGE NODES PER PATCH = ', SUM (NODES_PER) / N_PATCH
    PRINT *,'FIRST 10 ', NODES_PER (1:10)
    PRINT *,'LAST  10 ', NODES_PER (N_PATCH-9:N_PATCH)
    PRINT *,'Q4  PATCH COUNT = ', COUNT_4
    PRINT *,'Q9  PATCH COUNT = ', COUNT_9
    PRINT *,'Q16 PATCH COUNT = ', COUNT_16
    PRINT *,'Q17 PATCH COUNT = ', COUNT_17
    PRINT *,'TOTAL COUNT =     ', COUNT_4+COUNT_9+COUNT_16+COUNT_17

!     FINALLY, AVERAGE FLUXES    FOR EACH NODAL HIT COUNT
!b unit = get_next_io_unit ()  !b
!b open (unit, file = "scp_ave_dof.tmp", ACTION='WRITE', &  !b
      !b STATUS='REPLACE')  !b
    DO FIT = 1, MAX_NP ! FOR ALL NODES
      IF ( SCP_COUNTS (FIT) /= 0 ) THEN ! ACTIVE NODE
        SCP_AVERAGES (FIT, :) = SCP_AVERAGES (FIT, :) / SCP_COUNTS (FIT)
      ELSE  ! could skip since initialized
        SCP_AVERAGES (FIT, :) = 0.D0
      END IF   
!b write (unit, 5) SCP_AVERAGES (FIT, 1:N_G_DOF)  !b
!b 5 FORMAT ( (10(1X, 1PE12.5)) )  !b
    END DO ! FOR (AN UNWEIGHTED) AVERAGE 
!b xxx set logic flag assuring average has been done here

!      REPORT THE AVERAGED MAX & MIN VALUES AND LOCATIONS 
    !b IF ( F95 ) THEN ! Enhanced MAXLOC avilable
      !b !b CALL MAX_AND_MIN_SCP_AVE_F95 (SCP_AVERAGES) 
      !b CALL MAX_AND_MIN_SCP_AVE_F90 (SCP_AVERAGES) !b
    !b ELSE
      CALL MAX_AND_MIN_SCP_AVE_F90 (SCP_AVERAGES)
    !b END IF ! Enhanced MAXLOC avilable

!b print *,'first 5 node averages'
!b do JP=1,5
  !b print *, JP, SCP_AVERAGES (JP, :)
!b end do
!b PRINT *,'SCP_GEOM, SCP_LT, SCP_N, SCP_PARM, SCP_SHAP, SCP_FIT', &
         !b SCP_GEOM, SCP_LT, SCP_N, SCP_PARM, SCP_SHAP, SCP_FIT
END SUBROUTINE SURF_SCP_AVE_NODE_FLUXES    
