! File: model_lib.f
! copyright 2005,  J. E. Akin, all rights reserved.
! ---------------------------------------------------------------
!  Alphabetical listings follow
! ---------------------------------------------------------------

SUBROUTINE B_MATRIX_EQUAL_DGH (DGH, B)
! * * * * * * * * * * * * * * * * * * * * * * * * * * *
!  B == DGH FOR MOST SCALAR PROBLEMS, SO JUST COPY IT
! * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data  ! for LT_N, LT_FREE  
 IMPLICIT NONE
 REAL (DP), INTENT(IN)  :: DGH (N_SPACE, LT_N)
 REAL (DP), INTENT(OUT) :: B   (N_R_B,   LT_FREE)

   IF ( N_SPACE == N_R_B .AND. LT_N == LT_FREE )  THEN
     B (1:N_SPACE, 1:LT_N) = DGH
   ELSE ! INVALID CALL
     PRINT *,'ERROR: B_MATRIX_EQUAL_DGH, INCOMPATABLE MATRIX SIZES'
     PRINT *,'EXPECTED SIZE: ',   N_R_B, ' BY ', LT_FREE
     PRINT *,'RECEIVED SIZE: ', N_SPACE, ' BY ', LT_N    
     STOP 'INCOMPATABLE MATRIX SIZES, B_MATRIX_EQUAL_DGH'
   END IF
END SUBROUTINE B_MATRIX_EQUAL_DGH

!SUBROUTINE BOUNDARY_FLUX (FLUX, COORD_SEG, C_SEG, S_SEG, I_OPT, &
!                          NSQP, NSPARM, ISEG)
!! * * * * * * * * * * * * * * * * * * * * * * * * * *
!!    PROBLEM DEPENDENT BOUNDARY FLUX CONTRIBUTIONS
!!    >>> Outdated version, inactive in F90 version
!! * * * * * * * * * * * * * * * * * * * * * * * * * *
!Use System_Constants
!Use Sys_Properties_Data
!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),              &
!  ! EL_M (LT_FREE, LT_FREE), DIAG (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)  
!Use Interface_Header     ! for functions
! IMPLICIT NONE
!!    ALWAYS USED
! INTEGER,  INTENT(IN)  :: NSQP, NSPARM, ISEG
! REAL(DP), INTENT(IN)  :: COORD_SEG (L_B_N, N_SPACE)
! REAL(DP), INTENT(IN)  :: FLUX (L_B_N, N_G_FLUX)
! REAL(DP), INTENT(OUT) :: C_SEG (N_D_FLUX), S_SEG (N_D_FLUX, N_D_FLUX)
! INTEGER,  INTENT(OUT) :: I_OPT
!!          Automatic Arrays
! !b REAL(DP) :: XYZ (N_SPACE), AJ (LT_PARM, LT_PARM),              &
!             !b DGH (N_SPACE, LT_N), DGV (N_SPACE, LT_FREE),       &
!             !b AJ_INV (LT_PARM, LT_PARM), GPT (LT_QP), GWT (LT_QP)
!
!! AJ         = JACOBIAN
!! AJ_INV     = JACOBIAN INVERSE
!! C_SEG      = BOUNDARY FLUX COLUMN MATRIX CONTRIBUTIONS
!! COORD_SEG  = SPATIAL COORDINATES OF SEGMENT NODES
!! DGH        = GLOBAL DERIVATIVES SCALAR INTERPOLATION FUNCTIONS
!! DGV        = GLOBAL DERIVATIVES VECTOR INTERPOLATION FUNCTIONS
!! DLG        = LOCAL DERIVATIVES GEOMETRIC INTERPOLATION
!! DLH        = LOCAL DERIVATIVES SCALAR INTERPOLATION FUNCTIONS
!! DLV        = LOCAL DERIVATIVES VECTOR INTERPOLATION FUNCTIONS
!! FLT_SP     = REAL PROPERTIES ON THE SEGMENTS
!! FLUX       = SPECIFIED BOUNDARY FLUX COMPONENTS
!! G          = GEOMETRIC INTERPOLATION FUNCTIONS
!! H          = SOLUTION INTERPOLATION FUNCTIONS
!! I_OPT      = PROBLEM MATRIX REQUIREMENTS, MUST BE SET.
!!            = 1, CALCULATE C ONLY
!!            = 2, CALCULATE S ONLY
!!            = 3, CALCULATE BOTH C AND S
!! ISEG       = SEGMENT NUMBER
!! L_B_N      = NO. OF NODES ON AN ELEMENT BOUNDARY SEGMENT
!! L_HOMO     = 1 IF SEGMENT PROPERTIES ARE HOMOGENEOUS
!! NOD_PER_EL = NUMBER OF SOLUTION NODES, NOD_PER_EL = L_B_N USUALLY
!! N_BS_FIX   = NUMBER OF INTEGER PROPERTIES PER SEGMENT
!! N_BS_FLO   = NUMBER OF REAL PROPERTIES PER SEGMENT
!! SP_FIX     = INTEGER PROPERTIES ON THE SEGMENTS
!! N_D_FLUX   = L_B_N*N_G_DOF = MAXIMUM NUMBER OF FLUX CONTRIBUTIONS
!! N_G_DOF    = NUMBER OF PARAMETERS PER NODE POINT
!! N_PARM     = PARAMETRIC GEOMETRY NODES, = N_SPACE USUALLY
!! N_QP       = NUMBER OF QUADRATURE POINTS
!! N_F_SEG      = MAX NUMBER OF SEGMENTS
!! N_SPACE    = DIMENSION OF SOLUTION SPACE
!! PT         = QUADRATURE COORDINATES
!! S_SEG      = BOUNDARY FLUX SQUARE MATRIX
!! V          = VECTOR INTERPOLATION FUNCTIONS
!! WT         = QUADRATURE WEIGHTS
!! XYZ        = SPACE COORDINATES AT A POINT
!! ....................................................
!!   ** BOUNDARY_FLUX PROBLEM DEPENDENT STATEMENTS FOLLOW **
!! ....................................................
!! I_OPT       = PROBLEM MATRIX REQUIREMENTS, MUST BE SET.
!  S_SEG = 0.d0 ; C_SEG = 0.d0
!  I_OPT = 0 ! or REPLACE with: include 'my_flux_inc'
!! use print to supress compiler warnings
!if ( I_BUG > 0 ) print *,  NSQP, NSPARM, ISEG  !b
!if ( I_BUG > 0 ) print *,  COORD_SEG  !b
!if ( I_BUG > 0 ) print *,  FLUX  !b
!!       include 'my_flux_inc'
!END SUBROUTINE BOUNDARY_FLUX 

SUBROUTINE CALC_PROP_AT_PT (H, PRT_L_PT, VALUES)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!   CALCULATE N_NP_FLO PROPERTIES AT A LOCAL PT USING
!    ELEMENT'S NODAL PROPERTIES, PRT_L_PT, AND THE N
!      INTERPOLATION FUNCTIONS, H, AT THE POINT
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: H (NOD_PER_EL)
 REAL(DP), INTENT(IN)  :: PRT_L_PT (NOD_PER_EL, N_NP_FLO)
 REAL(DP), INTENT(OUT) :: VALUES (N_NP_FLO)

! NOD_PER_EL = NUMBER OF NODES PER ELEMENT
! N_NP_FLO   = NO. OF FLOATING POINT NODAL PROPERTIES
! H          = C^0 INTERPOLATION FUNCTIONS FOR AN ELEMENT
! PRT_L_PT   = FLOATING PT PROPS OF ELEMENT'S NODES
! VALUES     = LOCAL VALUES OF PROPERTIES

  IF ( N_NP_FLO < 1 ) STOP 'pt_real = 0, CALC_PROP_AT_PT'
  VALUES = MATMUL (H, PRT_L_PT)
END SUBROUTINE CALC_PROP_AT_PT
                              
SUBROUTINE CHANGE_IN_SOLUTION (DD, DD_OLD, OLD_NORM, NEW_NORM, &
                               DIFF_NORM, RATIO)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! CALCULATE THE MEAN CHANGE IN NODAL PARAMETERS FROM
!                  THE LAST ITERATION
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: DD (N_D_FRE),  DD_OLD(N_D_FRE)
 REAL(DP), INTENT(OUT) :: OLD_NORM,  NEW_NORM
 REAL(DP), INTENT(OUT) :: DIFF_NORM, RATIO
 REAL(DP), PARAMETER   :: EPS = EPSILON (1.d0)

! * CHANGE SHOULD BE CALLED BEFORE CORRECT_ITERATION *
! RATIO     = DIFF_NORM/NEW_NORM
! DIFF_NORM = SQRT(SUM OF (DD(I) - DD_OLD(I))**2)
! OLD_NORM  = SQRT(SUM OF DD_OLD(I)**2)
! NEW_NORM  = SQRT(SUM OF DD(I)**2)
! DD_OLD    = NODAL PARAMETER LIST FROM LAST ITERATION
! DD        = NODAL PARAMETERS FROM CURRENT ITERATION
! N_D_FRE   = TOTAL NO OF DEGREES OF FREEDOM IN SYS

  NEW_NORM  = SQRT ( SUM (DD**2) )
  OLD_NORM  = SQRT ( SUM (DD_OLD**2) )
  DIFF_NORM = SQRT ( SUM ((DD - DD_OLD) **2) )
  RATIO = DIFF_NORM / (NEW_NORM + EPS)
  WRITE (N_PRT, 5) OLD_NORM, NEW_NORM, DIFF_NORM, RATIO 
   5 FORMAT (                                               &
  '*** NODAL DOF FOR CURRENT AND PREVIOUS ITERATIONS ***',/,&
  'ROOT MEAN SQ OF PREVIOUS VALUES . . . ',1PE13.5,/,       &
  'ROOT MEAN SQ OF CURRENT VALUES  . . . ',1PE13.5,/,       &
  'ROOT MEAN SQ OF DIFFERENCES . . . . . ',1PE13.5,/,       &
  'RATIO OF LAST TWO QUANTITIES  . . . . ',1PE13.5,/)
END SUBROUTINE CHANGE_IN_SOLUTION
                                 
SUBROUTINE CHECK_EL_SHAPE
! * * * * * * * * * * * * * * * * * * * * * * * * * * *
!  CHECK L_SHAPE DATA FOR FREQUENT USER ERRORS
! * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 LOGICAL :: WARN 
! L_B_N      = NUMBER OF BOUNDARY NODES ON ELEM
! L_SHAPE    = ELEMENT SHAPE CODE
! NOD_PER_EL = NUMBER OF NODES PER ELEMENT
! N_SPACE    = DIMENSION OF ELEMENT SPACE

! >>> Note the CONTINUITY constant is now available here too
  WARN = .FALSE.

!   FORCE THE INPUT VALUE IF NEGATIVE
  IF ( L_SHAPE < 0 ) THEN
    L_SHAPE = IABS (L_SHAPE)
    PRINT *, 'WARNING, L_SHAPE SIGN CORRECTED'
    N_WARN = N_WARN + 1 ! INCREMENT WARNING
    RETURN
  END IF

!   LINE
  IF ( NOD_PER_EL == 2 .AND. L_SHAPE /= 1 ) THEN
    L_SHAPE = 1 ; WARN = .TRUE.               
  END IF
!   TRIANGLE
  IF ( N_SPACE > 1 .AND. NOD_PER_EL == 3 .AND. L_SHAPE /= 2 ) THEN
    L_SHAPE = 2 ; WARN = .TRUE.               
  END IF
!  QUADRILATERAL
  IF ( N_SPACE > 1 .AND. NOD_PER_EL == 4 .AND. L_SHAPE /= 3 ) THEN
    L_SHAPE = 3 ; WARN = .TRUE.               
  END IF
!   TRI OR WEDGE
  IF ( NOD_PER_EL == 6 .AND. L_SHAPE /= 2 ) THEN
    IF ( N_SPACE == 2 ) THEN
      L_SHAPE = 2 ; WARN = .TRUE.               
    END IF
  END IF
!   QUAD OR HEX
  IF ( NOD_PER_EL == 8 .AND. L_SHAPE /= 3 ) THEN
    IF ( N_SPACE /= 3 ) THEN
      L_SHAPE = 3 ; WARN = .TRUE.               
    END IF
  END IF
  IF ( WARN ) THEN
    WRITE (N_PRT, 1) L_SHAPE
    1 FORMAT ('********************************',/,        &
              '* WARNING L_SHAPE CHANGED TO', I2, ' *', /, &
              '********************************')
    N_WARN = N_WARN + 1 ! INCREMENT WARNING
  END IF

  WARN = .FALSE.
  IF ( L_B_N > 0 ) THEN
!    CHECK EDGE_FACE_VOLUME DATA (INCOMPLETE, NEEDS L_SHAPE ?) !bb xx
    IF (NOD_PER_EL == 2 .AND. L_B_N /= 1) WARN = .TRUE.              
    IF (NOD_PER_EL == 3 .AND. L_B_N /= 2) WARN = .TRUE.              
    IF (NOD_PER_EL == 4 .AND. L_B_N /= 2) WARN = .TRUE.              
    IF (NOD_PER_EL == 6 .AND. L_B_N /= 3) WARN = .TRUE.              
    IF (NOD_PER_EL == 8 .AND. L_B_N >  4) WARN = .TRUE.              
    IF ( WARN ) THEN
      PRINT *, 'WARNING: CHECK_EL_SHAPE, el_segment VALUE ', L_B_N
      PRINT *, 'UNEXPECTED FOR el_nodes VALUE ', NOD_PER_EL
      N_WARN = N_WARN + 1 ! INCREMENT WARNING
    END IF
  END IF

  WARN = .FALSE.
  IF ( M_B_N > 0 ) THEN
!    CHECK EDGE_FACE_VOLUME DATA (INCOMPLETE, NEEDS L_SHAPE ?) !bb xx
    IF (NOD_PER_EL == 2 .AND. M_B_N /= 1) WARN = .TRUE.              
    IF (NOD_PER_EL == 3 .AND. M_B_N /= 2) WARN = .TRUE.              
    IF (NOD_PER_EL == 4 .AND. M_B_N /= 2) WARN = .TRUE.              
    IF (NOD_PER_EL == 6 .AND. M_B_N /= 3) WARN = .TRUE.              
    IF (NOD_PER_EL == 8 .AND. M_B_N >  4) WARN = .TRUE.              
    IF ( WARN ) THEN
      PRINT *, 'WARNING: CHECK_EL_SHAPE, MIXED NODES VALUE ', M_B_N
      PRINT *, 'UNEXPECTED FOR el_nodes VALUE ', NOD_PER_EL
      N_WARN = N_WARN + 1 ! INCREMENT WARNING
    END IF
  END IF
END SUBROUTINE CHECK_EL_SHAPE
                             
SUBROUTINE CONSTRAINT_COUNT (N_RES_EQ, I_BC)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! CALCULATE NUMBER OF CONSTRAINT FLAGS OF EACH TYPE
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 INTEGER, INTENT(IN)    :: I_BC (MAX_NP)
 INTEGER, INTENT(INOUT) :: N_RES_EQ (MAX_TYP)
 INTEGER                :: I, ITEST, J, K
!           Automatic Arrays
 !b INTEGER :: KODES (N_G_DOF)
 INTEGER :: KODES (N_G_DOF), USED (N_G_DOF)

! MAX_NP   = TOTAL NUMBER OF SYSTEM NODES
! N_G_DOF  = NO. OF PARAMETERS (DOF) PER NODE
! I_BC     = NODAL POINT BOUNDARY RESTRAINT INDICATOR
! KODES    = LIST OF RESTRAINT INDICATORS AT A NODE
! N_RES_EQ = NUMBER OF CONSTR EQS ON EXIT, N_RES_EQ
! MAX_TYP  = MAX NO OF DIFFERENT CONSTRAINT TYPES
! MAX_ACT  = ACTIVE NO OF TYPES
! USED     = DATA CHECKING WORK ARRAY

  USED     = 0 ! INITIALIZATION
  N_RES_EQ = 0 ! INITIALIZATION
  DO I = 1, MAX_NP
!    DOES NODE I HAVE A NODAL PARAMETER CONSTRAINT
    ITEST = IABS (I_BC (I) )
    IF ( ITEST > 0 ) THEN
!      EXTRACT PARAMETER CODES
      CALL PT_CODE (I, ITEST, KODES)
      DO J = 1, N_G_DOF
        K = KODES (J)
!        UPDATE CONSTRAINT COUNTERS
        !b IF ( K > 0) N_RES_EQ (K) = N_RES_EQ (K) + 1
        IF ( K > 0 ) THEN
          N_RES_EQ (K) = N_RES_EQ (K) + 1
          USED     (J) = 1
        END IF ! USED
      END DO
    END IF
  END DO

!  CONVERT TO EQUATION COUNTERS
  N_CEQ = 0
  MAX_ACT = 0
  WRITE (N_PRT, 5)
  5 FORMAT ( /,'***  NODAL PARAMETER CONSTRAINT LIST  ***', /,&
               'TYPE      EQUATIONS')
  DO I = 1, MAX_TYP
    K = N_RES_EQ (I)
    IF ( K > 0) MAX_ACT = I
    IF ( ( (K / I) *I)  < K ) THEN
      WRITE (N_BUG,  * ) 'INVALID DATA FOR CONSTRAINT TYPE', I
      N_WARN = N_WARN + 1
    END IF
    N_RES_EQ (I) = N_RES_EQ (I) / I
    IF ( N_RES_EQ (I) > 0 ) WRITE (N_PRT, 5020) I, N_RES_EQ (I)
      5020 FORMAT ( I4, I10 )
    N_CEQ = N_CEQ + N_RES_EQ (I)
  END DO

!  VERIFY ALL DOF USED !b
  DO J = 1, N_G_DOF !b
    IF ( USED (J) <= 0 ) THEN !b
      PRINT *,'WARNING: PARAMETER ', J, ' HAS NO BOUNDARY CONDITION' !b
      N_WARN = N_WARN + 1 !b
    END IF ! UNUSED IN BC !b
  END DO ! VALIDATION !b
  !b IF ( ABS(N_CEQ) > N_D_FRE ) STOP 'impossible constraints'
  IF ( ABS(N_CEQ) > N_D_FRE ) THEN
    PRINT *,'ERROR: CONSTRAINT_COUNT EXCEEDS NUMBER OF EQUATIONS'
    STOP    'ERROR: CONSTRAINT_COUNT EXCEEDS NUMBER OF EQUATIONS'
  END IF ! IMPOSSIBLE DATA
END SUBROUTINE CONSTRAINT_COUNT

FUNCTION COPY_DGH_INTO_B_MATRIX (DGH) RESULT (B)
! * * * * * * * * * * * * * * * * * * * * * * * * * * *
!  B == DGH FOR MOST SCALAR PROBLEMS, SO JUST COPY IT
! * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data  ! for LT_N, LT_FREE  
 IMPLICIT NONE
 REAL (DP), INTENT(IN)  :: DGH (N_SPACE, LT_N)
 REAL (DP)              :: B   (N_R_B,   LT_FREE)

   IF ( N_SPACE == N_R_B .AND. LT_N == LT_FREE )  THEN
     B (1:N_SPACE, 1:LT_N) = DGH
   ELSE ! INVALID CALL
     STOP 'INCOMPATABLE MATRIX SIZES, COPY_DGH_INTO_B_MATRIX'
   END IF
END FUNCTION COPY_DGH_INTO_B_MATRIX
                               
SUBROUTINE PUT_DGH_INTO_B_MATRIX (DGH, B)
! * * * * * * * * * * * * * * * * * * * * * * * * * * *
!  B == DGH FOR MOST SCALAR PROBLEMS, SO JUST COPY IT
! * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data  ! for LT_N, LT_FREE  
 IMPLICIT NONE
 REAL (DP), INTENT(IN)  :: DGH (N_SPACE, LT_N)
 REAL (DP)              :: B   (N_R_B,   LT_FREE)

   IF ( N_SPACE == N_R_B .AND. LT_N == LT_FREE )  THEN
     B (1:N_SPACE, 1:LT_N) = DGH
   ELSE ! INVALID CALL
     STOP 'INCOMPATABLE MATRIX SIZES, PUT_DGH_INTO_B_MATRIX'
   END IF
END SUBROUTINE PUT_DGH_INTO_B_MATRIX
                               
SUBROUTINE CORRECT_ITERATION (DD, DD_OLD)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! CALCULATE NEW STARTING VALUES FOR NEXT ITERATION
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!    OVER RELAXATION METHOD
Use System_Constants
 IMPLICIT NONE
 REAL(DP), INTENT(IN)    :: DD (N_D_FRE)
 REAL(DP), INTENT(INOUT) :: DD_OLD (N_D_FRE)
 REAL(DP) :: OMEGA (N_G_DOF) 
 INTEGER  :: I, J, K

! DD      = CALCULATED DOF FROM LAST ITERATION
! DD_OLD  = DOF TO BE USED TO START NEXT ITERATION
! N_D_FRE = TOTAL NO OF SYS DEGREES OF FREEDOM
! N_G_DOF = NUMBER OF PARAMETERS (DOF) PER NODE
! OMEGA   = RELAXATION FACTOR FOR EACH NODAL DOF
  !b OMEGA = 1.25_DP    ! TYPICAL VALUE

 IF ( SPACE_TIME .AND. (.NOT. S_T_CONTINUOUS) ) THEN
   IF ( MAX_NP /= (MAX_NP/2)*2 ) THEN ! not space time
     PRINT *,'WARNING, CORRECT_ITERATION: Odd max node number'
     N_WARN = N_WARN + 1
   END IF
!  move t_(n+1) answers to t_(n) for weak BC arrays
     I = N_D_FRE / 2
     DD_OLD (1:I) = DD_OLD (I+1:I+N_D_FRE/2)
 ELSE ! standard iteration
  OMEGA = RELAXATION ! TYPICAL VALUE
  I = 0 
  DO J = 1, MAX_NP 
    DO K = 1, N_G_DOF 
      I = I + 1
      DD_OLD (I) = DD_OLD (I) + OMEGA (K) * (DD (I) - DD_OLD (I) )
    END DO
  END DO
 END IF
END SUBROUTINE CORRECT_ITERATION
                                
FUNCTION COUNT_EL_TYPE_NODES  (LT) RESULT (COUNT)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!   GET DATA BASED ON ELEMENT TYPE (SEE CONTROL_INPUT )
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data     ! for LT_DATA (MAX_DAT, N_L_TYPE)
 IMPLICIT NONE
 INTEGER, INTENT(IN)  :: LT
 INTEGER              :: COUNT
! LT       = ELEMENT TYPE NUMBER
! N_L_TYPE = NUMBER OF ELEMENT TYPES
! LT_N     = NUMBER OF NODES PER ELEMENT
    LT_N  = LT_DATA (1, LT) 
    COUNT = LT_N
END FUNCTION COUNT_EL_TYPE_NODES
                                
!SUBROUTINE DEGPAR (I_PT, J_PARM, N_G_DOF, INDEX) 
!! * * * * * * * * * * * * * * * * * * * * * * *
!!  DETERMINE THE DEGREE OF FREEDOM NUMBER
!!  OF NODAL PARAMETER J_PARM AT NODE POINT I_PT
!! * * * * * * * * * * * * * * * * * * * * * * *
!Use Precision_Module
! IMPLICIT NONE
! INTEGER, INTENT(IN)  :: I_PT, J_PARM, N_G_DOF
! INTEGER, INTENT(OUT) :: INDEX
!
!! I_PT    = A SPECIFIC NODE NUMBER
!! J_PARM  = THE PARAMETER AT THAT NODE
!! N_G_DOF = NUMBER OF PARAMETERS PER NODE
!
!  !b IF ( DOF_VARY ) STOP 'DEGPAR: INVALID FOR DOF_VARY'
!
!  INDEX = N_G_DOF * (I_PT - 1) + J_PARM
!END SUBROUTINE DEGPAR
                     
SUBROUTINE ELEM_COORD (NOD_ON_EL, L_SPACE, X, COORD, ELEM_NODES)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! DETERMINE COORDINATES OF NODES ON ELEMENT
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 INTEGER,  INTENT(IN)  :: NOD_ON_EL, L_SPACE
 REAL(DP), INTENT(IN)  :: X (MAX_NP, L_SPACE)
 INTEGER,  INTENT(IN)  :: ELEM_NODES (NOD_ON_EL)
 REAL(DP), INTENT(OUT) :: COORD (NOD_ON_EL, L_SPACE)
 INTEGER :: I

! MAX_NP     = NUMBER OF NODES IN SYSTEM
! L_SPACE    = DIMENSION OF SPACE
! NOD_ON_EL  = NUMBER OF NODES PER ELEMENT
! X          = COORDINATES OF SYSTEM NODES
! COORD      = COORDINATES OF ELEMENT NODES
! ELEM_NODES = NOD_ON_EL INCIDENCES OF ELEMENT

  DO I = 1, NOD_ON_EL
!     ALLOW FOR OMITTED NODES
    IF ( ELEM_NODES (I) > 0) &
         COORD (I, 1:L_SPACE) = X (ELEM_NODES (I), 1:L_SPACE)
  END DO
END SUBROUTINE ELEM_COORD
                     
FUNCTION GET_ELEM_COORD (X, ELEM_NODES) RESULT (COORD)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! DETERMINE COORDINATES OF NODES ON ELEMENT
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 REAL(DP), INTENT(IN) :: X (MAX_NP, N_SPACE)
 INTEGER, INTENT(IN)  :: ELEM_NODES (NOD_PER_EL)
 REAL(DP)             :: COORD (NOD_PER_EL, N_SPACE)
 INTEGER :: I

! MAX_NP     = NUMBER OF NODES IN SYSTEM
! N_SPACE    = DIMENSION OF SPACE
! NOD_PER_EL = NUMBER OF NODES PER ELEMENT
! X          = COORDINATES OF SYSTEM NODES
! COORD      = COORDINATES OF ELEMENT NODES
! ELEM_NODES = NOD_PER_EL INCIDENCES OF ELEMENT

  DO I = 1, NOD_PER_EL
!      ALLOW FOR OMITTED NODES
    IF ( ELEM_NODES (I) > 0) &   
         COORD (I, 1:N_SPACE) = X (ELEM_NODES (I), 1:N_SPACE)
  END DO
END FUNCTION GET_ELEM_COORD
                     
SUBROUTINE ELEMENT_NODES (L_ID, LT_N, NODES, ELEM_NODES)
! * * * * * * * * * * * * * * * * * * * * * * * *
!   EXTRACT NODES ASSOCIATED WITH ELEMENT L_ID
! * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 INTEGER, INTENT(IN)  :: NODES (L_S_TOT, NOD_PER_EL), L_ID, LT_N
 INTEGER, INTENT(OUT) :: ELEM_NODES (LT_N)

! N_ELEMS    = NUMBER OF ELEMENTS IN SYSTEM
! NOD_PER_EL = MAXIMUN NUMBER OF NODES PER ELEMENT
! LT_N       = NUMBER OF NODES PER ELEMENT
! L_ID       = ELEMENT NUMBER
! NODES      = NODAL INCIDENCES OF ALL ELEMENTS
! ELEM_NODES = THE NOD_PER_EL NODAL INCIDENCES OF THE ELEMENT

  IF ( LT_N > NOD_PER_EL ) STOP 'DATA ERROR ELEMENT_NODES'
  ELEM_NODES (1:LT_N) = NODES (L_ID, 1:LT_N)
END SUBROUTINE ELEMENT_NODES
                            
SUBROUTINE EL_FRE (N_D_FRE, N_EL_FRE, D, DD, INDEX)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! EXTRACT ELEMENT DEGREES OF FREEDOM FROM SYSTEM DOF
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 IMPLICIT NONE
 INTEGER,  INTENT(IN)  :: N_D_FRE, N_EL_FRE, INDEX (N_EL_FRE)
 REAL(DP), INTENT(IN)  :: DD (N_D_FRE)
 REAL(DP), INTENT(OUT) :: D  (N_EL_FRE)

! D        = NODAL PARAMETERS ASSOCIATED WITH ELEMENT
! DD       = SYSTEM ARRAY OF NODAL PARAMETERS
! INDEX    = ARRAY OF SYSTEM DEGREE OF FREEDOM NUMBERS
! N_EL_FRE = NUMBER OF DEGREES OF FREEDOM PER ELEMENT
! N_D_FRE  = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM

  WHERE ( INDEX > 0 ) 
    D = DD (INDEX)
  ELSEWHERE
    D = 0.d0
  END WHERE
END SUBROUTINE EL_FRE
                    
FUNCTION GET_ELEM_DOF (DD) RESULT (D_LT)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! EXTRACT ELEMENT DEGREES OF FREEDOM FROM SYSTEM DOF
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants    ! for N_D_FRE
Use Elem_Type_Data      ! for LT_FREE, INDEX (LT_FREE)
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: DD    (N_D_FRE)
 REAL(DP)              :: D_LT  (LT_FREE)

! D_LT     = NODAL PARAMETERS ASSOCIATED WITH ELEMENT
! DD       = SYSTEM ARRAY OF NODAL PARAMETERS
! INDEX    = ARRAY OF SYSTEM DEGREE OF FREEDOM NUMBERS
! LT_FREE  = NUMBER OF DEGREES OF FREEDOM FOR ELEMENT TYPE
! N_D_FRE  = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM

  WHERE ( INDEX > 0 ) 
    D_LT = DD (INDEX)
  ELSEWHERE
    D_LT = 0.d0
  END WHERE
END FUNCTION GET_ELEM_DOF
                    
SUBROUTINE GET_TYPE_DOF (LT_D, DD, LT_INDEX)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! GET ELEMENT TYPE DEGREES OF FREEDOM FROM SYSTEM DOF
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants    ! for N_D_FREE
Use Elem_Type_Data      ! for LT_FREE
 IMPLICIT NONE
 INTEGER,  INTENT(IN)  :: LT_INDEX (LT_FREE)
 REAL(DP), INTENT(IN)  :: DD       (N_D_FRE)
 REAL(DP), INTENT(OUT) :: LT_D     (LT_FREE)

! LT_D     = NODAL PARAMETERS ASSOCIATED WITH ELEMENT
! DD       = SYSTEM ARRAY OF NODAL PARAMETERS
! LT_INDEX = ARRAY OF SYSTEM DEGREE OF FREEDOM NUMBERS
! LT_FREE  = NUMBER OF DEGREES OF FREEDOM FOR ELEMENT TYPE
! N_D_FRE  = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM

  WHERE ( LT_INDEX > 0 ) 
    LT_D = DD (LT_INDEX)
  ELSEWHERE
    LT_D = 0.d0
  END WHERE
END SUBROUTINE GET_TYPE_DOF
                    
SUBROUTINE ELEM_PROPERTIES (L_ID, L_PROP, EL_PROP)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! EXTRACT PROPERTIES OF ELEMENT, L_ID, FROM TOTAL 
! PROPERTIES ARRAYS  >>> now inactive, but correct
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Sys_Properties_Data
 IMPLICIT NONE
 INTEGER,  INTENT(IN)  :: L_ID
 REAL(DP), INTENT(OUT) :: EL_PROP (N_LP_FLO)
 INTEGER,  INTENT(OUT) :: L_PROP  (N_LP_FIX)
 INTEGER :: I

! L_ID     = ELEMENT NUMBER
! LP_FIX   = SYSTEM ARRAY OF FIXED PT ELEM PROPERTIES
! L_PROP   = ELEM FIXED PT PROPERTIES ARRAY
! FLT_EL   = SYS ARRAY OF FLOATING PT NODAL PROP
! EL_PROP  = ELEM FLOATING PT PROPERTIES ARRAY
! L_HOMO   = 1, IF PROPERTIES ARE SAME IN ALL ELEMENTS
! N_LP_FIX = NUMBER OF INTEGER ELEMENT PROPERTIES
! N_LP_FLO = NUMBER OF REAL ELEMENT PROPERTIES

  I = L_ID
  IF ( L_HOMO == 1 ) I = 1
!  FLOATING POINT PROPERTIES
  EL_PROP (1:N_LP_FLO) = FLT_EL (I, 1:N_LP_FLO)
!  FIXED POINT PROPERTIES
  L_PROP (1:N_LP_FIX) = LP_FIX (I, 1:N_LP_FIX)
END SUBROUTINE ELEM_PROPERTIES

SUBROUTINE GENERATE_EL_MATRICES_NEW (IE, X, DD_OLD)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!      GENERATE ELEMENT MATRICES AND POST SOLUTION DATA
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data       ! for PT (LT_PARM, LT_QP), WT (LT_QP),
                         ! G (LT_GEOM, LT_QP), H (LT_N   ),
                         ! DLG (LT_PARM, LT_GEOM, LT_QP),
                         ! DLH (LT_PARM, LT_N   , LT_QP),
                         ! EL_M (LT_FREE, LT_FREE), DIAG_M (LT_FREE),
                         ! S (LT_FREE, LT_FREE), C (LT_FREE),
                         ! ELEM_NODES (LT_N), 
                         ! COORD (LT_N, N_SPACE),
                         ! GEOMETRY (LT_GEOM, N_SPACE)
                         ! INDEX (LT_FREE), D (LT_FREE)
Use Sys_Properties_Data  ! for FLT_MISC, MISC_FIX, DELETED
Use Select_Source
Use Interface_Header
 IMPLICIT NONE
 INTEGER,  INTENT(IN)  :: IE
!     SYSTEM DATA
 REAL(DP), INTENT(IN) :: X (MAX_NP, N_SPACE), DD_OLD (N_D_FRE)
!          Automatic Arrays
 REAL(DP) :: E (N_R_B, N_R_B) 
 REAL(DP) :: H_INTG (LT_N   )
 REAL(DP) :: PRT_L_PT (LT_N, N_NP_FLO)
 REAL(DP) :: PRT_MAT (MISC_FL) 
 REAL(DP) :: FLUX (L_B_N, N_G_FLUX)  !b
 INTEGER  :: L_PT_PROP (LT_N, N_NP_FIX)

!                 VARIABLES:
! COORD      = SPATIAL COORDINATES OF ELEMENT'S NODES
! D          = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT
! DD_OLD     = SYSTEM NODAL PARAMETERS FROM LAST ITERATION
! DGH        = GLOBAL DERIVATIVES INTERPOLATION FUNCTIONS
! DLG        = LOCAL DERIVATIVES GEOMETRIC INTERPOLATION
! DLH        = LOCAL DERIVATIVES INTERPOLATION FUNCTIONS
! E          = CONSTITUTIVE MATRIX
! FLT_EL     = REAL PROPERTIES OF ELEMENTS
! FLT_MISC   = MISCELLANEOUS REAL PROPERTIES OF SYSTEM
! FLT_NP     = REAL PROPERTIES OF SYSTEM NODES
! FLT_SP     = REAL PROPERTIES OF SYSTEM SEGMENTS
! G          = GEOMETRIC INTERPOLATION FUNCTIONS
! GEOMETRY   = COORDINATES OF ELEM GEOMETRIC NODES
! H          = SOLUTION INTERPOLATION FUNCTIONS
! H_INTG     = INTEGRAL OF INTERPOLATION FUNCTIONS
! IE         = CURRENT ELEMENT NUMBER
! INDEX      = SYSTEM DOF NUMBERS ASSOCIATED WITH ELEMENT
! L_HOMO     = 1, IF ELEMENT PROPERTIES ARE HOMOGENEOUS
! ELEM_NODES = THE NOD_PER_EL INCIDENCES OF THE ELEMENT
! LP_FIX     = SYSTEM ARRAY OF INTEGER ELEM PROPERTIES
! L_PT_PROP  = INTEGER PROPERTIES AT EACH ELEMENT NODE
! LP_TEST    > 0, IF ELEMENT PROPERTIES HAVE BEEN DEFINED
! LT         = ELEMENT TYPE NUMBER (IF USED)
! LT_USER    = OPTIONAL USER CODE ASSOCIATED WITH TYPE
! MAX_NP     = NUMBER OF SYSTEM NODES
! MISC_FIX   = MISCELLANEOUS INTEGER SYSTEM PROPERTIES
! LT_N       = NUMBER OF NODES PER ELEMENT
! N_D_FRE    = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM
! N_ELEMS    = NUMBER OF ELEMENTS
! LT_FREE    = NUMBER OF DEGREES OF FREEDOM PER ELEMENT
! LT_GEOM    = NUMBER OF GEOMETRY NODES
! N_HOMO     = 1, IF NODAL PROPERTIES ARE HOMOGENEOUS
! N_ITER     = NO. OF ITERATIONS TO BE RUN (USUALLY 1)
! N_MAT      = MATERIAL TYPE NUMBER
! LT_PARM    = DIMENSION OF PARAMETRIC SPACE
! NP_FIX     = INTEGER PROPERTIES AT ALL NODES
! LT_QP      = NUMBER OF QUADRATURE POINTS
! N_R_B      = NUMBER OF ROWS IN B AND E MATRICES
! N_SPACE    = DIMENSION OF SPACE
! N_FILE1    = UNIT FOR POST SOLUTION MATRICES STORAGE
! N_FILE2    = OPTIONAL UNIT FOR USER (USED WHEN > 0)
! N_FILE#    = OPTIONAL UNITS FOR USER (USED WHEN > 0)
! NUL_COL    = 1, IF ELEMENT COLUMN MATRIX IS ALWAYS ZERO
! PRT_L_PT   = REAL PROPERTIES AT ELEMENT NODES
! PRT_MAT    = REAL ELEM PROPERTIES BASED ON MATERIAL NUMBER
! PT         = QUADRATURE COORDINATES
! S          = ELEMENT SQUARE MATRIX
! SP_FIX     = INTEGER PROPERTIES AT ALL SEGMENTS
! WT         = QUADRATURE WEIGHTS
! X          = COORDINATES OF SYSTEM NODES

  IF ( I_BUG > 0 .AND. (N_ITER > 0 .OR. TRANSIENT .OR. DYNAMIC)) THEN
    PRINT *,'Enter GENERATE_EL_MATRICES_NEW with DD_OLD'
    !b CALL RPRINT (DD_OLD, 1, N_D_FRE, 1)
  END IF

!  GET ELEMENT GROUP: STANDARD, MIXED BC, SEGMENT, USER DEFINED
  CALL SET_TO_MIXED_OR_SEG_OR_STD (IE) ! For IS_STD, IS_MIX, IS_SEG

!-->   EXTRACT ANALYSIS AND GEOMETRIC NODAL COORDINATES
  CALL ELEM_COORD (LT_N,    N_SPACE, X, COORD,    ELEM_NODES)
  CALL ELEM_COORD (LT_GEOM, N_SPACE, X, GEOMETRY, ELEM_NODES(1:LT_GEOM))

!-->    EXTRACT NODAL PARAMETERS FROM LAST ITERATION (IF ANY)
  IF ( N_ITER > 1 .OR. TRANSIENT ) D = GET_ELEM_DOF (DD_OLD)
  !b call at(693); print *,'DD_OLD ', DD_OLD ; print *, 'D ', THIS_EL, D

!-->    EXTRACT NODAL POINT PROPERTIES (IF ANY)
  IF ( N_NP_FLO > 0) CALL ELEM_NODE_PROPERTIES (LT_N, PRT_L_PT, &
                                                L_PT_PROP, ELEM_NODES)

!b call at(640)
!-->   EXTRACT REAL MATERIAL PROPERTIES (IF ANY)
  IF ( N_MAT > 0) CALL MATERIAL_PROP (N_MAT, PRT_MAT)

!-->   GENERATE ELEMENT SQUARE AND COLUMN MATRICES
  E      = GET_REAL_IDENTITY (N_R_B) ! DEFAULT TO IDENTITY MATRIX
  H_INTG = 0.d0                      ! ZERO INTEGRAL
                                         
  IF (IS_ELEMENT) THEN ! Get sq and maybe column matrices
    CALL SELECT_ELEM_SQ_MATRIX (E, H_INTG, PRT_L_PT, PRT_MAT,     &
                                L_PT_PROP, IE, X) !b for space frame
    IF ( DEBUG_EL_COL .AND. IE <= 3 )  CALL RPRINT &!b
      (C, 1, LT_FREE, 1) !b
    IF ( NUL_COL == 0 ) THEN ! Need column matrix
      CALL SELECT_ELEM_COL_MATRIX (E, H_INTG, PRT_L_PT,  &
                                   PRT_MAT, L_PT_PROP, IE)
    END IF ! Non-zero column matrix
    IF ( DEBUG_EL_COL .AND. IE <= 3 )  CALL RPRINT &!b
      (C, 1, LT_FREE, 1) !b
  END IF ! Standard element

  IF (IS_MIXED_BC) CALL SELECT_MIXED_SQ_MATRIX (E, H_INTG, PRT_L_PT,&
                                               PRT_MAT, L_PT_PROP, IE)
  IF (IS_FLUX_SEG) THEN ! Get and use boundary flux input
!-->   EXTRACT FLUX IF BOUNDARY SEGMENT ELEMENT
    IF ( IS_SEG > 0 ) FLUX (:, :) = S_FLUX(:, :, IS_SEG) !b make call
    CALL SELECT_SEG_COL_MATRIX  (E, H_INTG, PRT_L_PT, PRT_MAT, &
                                 L_PT_PROP, IE, FLUX)
  END IF ! Boundary flux source

!-->   STORE DATA FOR POST SOLUTION CALCULATIONS (IF ANY)
  IF ( N_FILE1 > 0 .OR. N_FILE2 > 0 )                         &
    CALL SELECT_ELEM_POST_DATA (E, H_INTG, PRT_L_PT, PRT_MAT, &
                                L_PT_PROP, IE, FLUX)
!  NOTE:  SYSTEM PROPERTIES UPDATE COULD BE DONE HERE
END SUBROUTINE GENERATE_EL_MATRICES_NEW
                              
SUBROUTINE GENERATE_EL_MATRICES (IE, X, DD_OLD)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!      GENERATE ELEMENT MATRICES AND POST SOLUTION DATA
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data       ! for PT (LT_PARM, LT_QP), WT (LT_QP),
                         ! G (LT_GEOM, LT_QP), H (LT_N   ),
                         ! DLG (LT_PARM, LT_GEOM, LT_QP),
                         ! DLH (LT_PARM, LT_N   , LT_QP),
                         ! S (LT_FREE, LT_FREE), C (LT_FREE),
                         ! ELEM_NODES (LT_N), 
                         ! COORD (LT_N,    N_SPACE),
                         ! GEOMETRY (LT_GEOM, N_SPACE)
                         ! INDEX (LT_FREE), D (LT_FREE)
Use Sys_Properties_Data  ! for FLT_MISC, MISC_FIX, DELETED
Use Select_Source
Use Interface_Header
 IMPLICIT NONE
 INTEGER,  INTENT(IN)  :: IE
!     SYSTEM DATA
 REAL(DP), INTENT(IN) :: X (MAX_NP, N_SPACE), DD_OLD (N_D_FRE)
!          Automatic Arrays
 REAL(DP) :: E (N_R_B, N_R_B) 
 REAL(DP) :: H_INTG (LT_N   )
 !b REAL(DP) :: PRT_L_PT (LT_N, N_NP_FLO), EL_PROP (N_LP_FLO)
 REAL(DP) :: PRT_L_PT (LT_N, N_NP_FLO)
 REAL(DP) :: PRT_MAT (MISC_FL) 
 REAL(DP) :: FLUX (L_B_N, N_G_FLUX)  !b
 !b INTEGER  :: L_PROP (N_LP_FIX), L_PT_PROP (LT_N, N_NP_FIX)
 INTEGER  :: L_PT_PROP (LT_N, N_NP_FIX)
 !b INTEGER  :: I_SEG !b

!                 VARIABLES:
! COORD      = SPATIAL COORDINATES OF ELEMENT'S NODES
! D          = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT
! DD_OLD     = SYSTEM NODAL PARAMETERS FROM LAST ITERATION
! DGH        = GLOBAL DERIVATIVES INTERPOLATION FUNCTIONS
! DLG        = LOCAL DERIVATIVES GEOMETRIC INTERPOLATION
! DLH        = LOCAL DERIVATIVES INTERPOLATION FUNCTIONS
! E          = CONSTITUTIVE MATRIX
! FLT_EL     = REAL PROPERTIES OF ELEMENTS
! FLT_MISC   = MISCELLANEOUS REAL PROPERTIES OF SYSTEM
! FLT_NP     = REAL PROPERTIES OF SYSTEM NODES
! FLT_SP     = REAL PROPERTIES OF SYSTEM SEGMENTS
! G          = GEOMETRIC INTERPOLATION FUNCTIONS
! GEOMETRY   = COORDINATES OF ELEM GEOMETRIC NODES
! H          = SOLUTION INTERPOLATION FUNCTIONS
! H_INTG     = INTEGRAL OF INTERPOLATION FUNCTIONS
! IE         = CURRENT ELEMENT NUMBER
! INDEX      = SYSTEM DOF NUMBERS ASSOCIATED WITH ELEMENT
! L_HOMO     = 1, IF ELEMENT PROPERTIES ARE HOMOGENEOUS
! ELEM_NODES = THE NOD_PER_EL INCIDENCES OF THE ELEMENT
! LP_FIX     = SYSTEM ARRAY OF INTEGER ELEM PROPERTIES
! L_PT_PROP  = INTEGER PROPERTIES AT EACH ELEMENT NODE
! LP_TEST    > 0, IF ELEMENT PROPERTIES HAVE BEEN DEFINED
! LT         = ELEMENT TYPE NUMBER (IF USED)
! LT_USER    = OPTIONAL USER CODE ASSOCIATED WITH TYPE
! MAX_NP     = NUMBER OF SYSTEM NODES
! MISC_FIX   = MISCELLANEOUS INTEGER SYSTEM PROPERTIES
! LT_N       = NUMBER OF NODES PER ELEMENT
! N_D_FRE    = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM
! N_ELEMS    = NUMBER OF ELEMENTS
! LT_FREE    = NUMBER OF DEGREES OF FREEDOM PER ELEMENT
! LT_GEOM    = NUMBER OF GEOMETRY NODES
! N_HOMO     = 1, IF NODAL PROPERTIES ARE HOMOGENEOUS
! N_ITER     = NO. OF ITERATIONS TO BE RUN (USUALLY 1)
! N_MAT      = MATERIAL TYPE NUMBER
! LT_PARM    = DIMENSION OF PARAMETRIC SPACE
! NP_FIX     = INTEGER PROPERTIES AT ALL NODES
! LT_QP      = NUMBER OF QUADRATURE POINTS
! N_R_B      = NUMBER OF ROWS IN B AND E MATRICES
! N_SPACE    = DIMENSION OF SPACE
! N_FILE1    = UNIT FOR POST SOLUTION MATRICES STORAGE
! N_FILE2    = OPTIONAL UNIT FOR USER (USED WHEN > 0)
! N_FILE#    = OPTIONAL UNITS FOR USER (USED WHEN > 0)
! NUL_COL    = 1, IF ELEMENT COLUMN MATRIX IS ALWAYS ZERO
! PRT_L_PT   = REAL PROPERTIES AT ELEMENT NODES
! PRT_MAT    = REAL ELEM PROPERTIES BASED ON MATERIAL NUMBER
! PT         = QUADRATURE COORDINATES
! S          = ELEMENT SQUARE MATRIX
! SP_FIX     = INTEGER PROPERTIES AT ALL SEGMENTS
! WT         = QUADRATURE WEIGHTS
! X          = COORDINATES OF SYSTEM NODES
 
  IF ( DEBUG_ASSEMBLY ) THEN
    PRINT *,'Enter GENERATE_EL_MATRICES with DD_OLD'
    CALL RPRINT (DD_OLD, 1, N_D_FRE, 1)
  END IF

!  GET ELEMENT GROUP: STANDARD, MIXED BC, SEGMENT, USER DEFINED
  CALL SET_TO_MIXED_OR_SEG_OR_STD (IE) ! For IS_STD, IS_MIX, IS_SEG

!-->   EXTRACT ANALYSIS AND GEOMETRIC NODAL COORDINATES
  CALL ELEM_COORD (LT_N,    N_SPACE, X, COORD,    ELEM_NODES)
  CALL ELEM_COORD (LT_GEOM, N_SPACE, X, GEOMETRY, ELEM_NODES(1:LT_GEOM))

!-->    EXTRACT NODAL PARAMETERS FROM LAST ITERATION (IF ANY)
  IF ( TRANSIENT .OR. N_ITER > 1) D = GET_ELEM_DOF (DD_OLD)
  !b call at(829); print *,'DD_OLD ', DD_OLD ; print *, 'D ', THIS_EL, D

!-->    EXTRACT NODAL POINT PROPERTIES (IF ANY)
  IF ( N_NP_FLO > 0) CALL ELEM_NODE_PROPERTIES (LT_N, PRT_L_PT, &
                                                L_PT_PROP, ELEM_NODES)

!-->   EXTRACT REAL MATERIAL PROPERTIES (IF ANY)
  IF ( N_MAT > 0) CALL MATERIAL_PROP (N_MAT, PRT_MAT)

!-->   GENERATE ELEMENT SQUARE AND COLUMN MATRICES
  E      = GET_REAL_IDENTITY (N_R_B) ! DEFAULT TO IDENTITY MATRIX
  H_INTG = 0.d0                      ! ZERO INTEGRAL
  IF (IS_ELEMENT ) THEN ! Get sq and maybe column matrices
     CALL ELEM_SQ_MATRIX (E, H_INTG, PRT_L_PT, PRT_MAT, &
                          L_PT_PROP, IE, X) !b for space frame
                          !b L_PT_PROP, IE)
     IF ( NUL_COL == 0 ) CALL ELEM_COL_MATRIX (E, H_INTG, PRT_L_PT,  &
                                               PRT_MAT, L_PT_PROP, IE)
  END IF ! Standard element

  !b IF (IS_MIXED_BC) CALL MIXED_SQ_MATRIX (E, H_INTG, PRT_L_PT,    &
  IF (IS_MIXED_BC) CALL SELECT_MIXED_SQ_MATRIX (E, H_INTG, PRT_L_PT,&
                                         PRT_MAT, L_PT_PROP, IE)
  IF (IS_FLUX_SEG) THEN
!-->   EXTRACT FLUX IF BOUNDARY SEGMENT ELEMENT
    IF ( IS_SEG > 0 ) FLUX (:, :) = S_FLUX(:, :, IS_SEG)
    CALL SEG_COL_MATRIX  (E, H_INTG, PRT_L_PT, PRT_MAT, L_PT_PROP, &
                          IE, FLUX)
  END IF

!-->   STORE DATA FOR POST SOLUTION CALCULATIONS (IF ANY)
  IF ( N_FILE1 > 0 .OR. N_FILE2 > 0 )                           &
    CALL ELEM_POST_DATA (E, H_INTG, PRT_L_PT, PRT_MAT, &
                         L_PT_PROP, IE, FLUX)
!  NOTE:  SYSTEM PROPERTIES UPDATE COULD BE DONE HERE
END SUBROUTINE GENERATE_EL_MATRICES
                                   
FUNCTION GET_PT_DOF_INDEX (I_PT, J_PARM) RESULT (INDEX)
! * * * * * * * * * * * * * * * * * * * * * * *
!  DETERMINE THE DEGREE OF FREEDOM NUMBER
!  OF NODAL PARAMETER J_PARM AT NODE POINT I_PT
! * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants    ! for N_G_DOF
 IMPLICIT NONE
 INTEGER, INTENT(IN) :: I_PT, J_PARM
 INTEGER             :: INDEX

! DOF_PT     = NUMBER OF DOF AT EACH NODE:  (0:MAX_NP)
! DOF_PT_SUM = TOTAL NUMBER OF DOF, INCLUDING THIS NODE: (-1:MAX_NP)
! I_PT       = A SPECIFIC NODE NUMBER
! J_PARM     = THE PARAMETER AT THAT NODE
! N_G_DOF    = NUMBER OF PARAMETERS PER NODE

  IF ( DOF_VARY ) THEN ! VARIABLE DOF PER NODE
    INDEX = DOF_PT_SUM (I_PT-1) + J_PARM
    IF ( J_PARM > DOF_PT (I_PT) ) THEN
      PRINT *,'WARNING: INVALID DOF ', J_PARM, ' AT NODE ', I_PT
      N_WARN = N_WARN + 1
      INDEX = 0
    END IF
  ELSE ! ALL NODES HAVE N_G_DOF
    INDEX = N_G_DOF * (I_PT - 1) + J_PARM
  END IF
END FUNCTION GET_PT_DOF_INDEX
                     
FUNCTION GET_INDEX_AT_PT (I_PT) RESULT (INDEX)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!  DETERMINE DEGREES OF FREEDOM NUMBERS AT A NODE
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants ! for N_G_DOF, DOF_PT DOF_PT_SUM
 IMPLICIT NONE
 INTEGER, INTENT(IN)  :: I_PT
 INTEGER              :: INDEX (N_G_DOF)
 INTEGER :: EQ_PT, J ! implied loop

! DOF_PT     = NUMBER OF DOF AT EACH NODE:  (0:MAX_NP)
! DOF_PT_SUM = TOTAL NUMBER OF DOF, INCLUDING THIS NODE: (-1:MAX_NP)
! N_G_DOF    = NUMBER OF PARAMETERS (DOF) PER NODE
! I_PT       = SYSTEM NODE NUMBER
! INDEX      = SYSTEM DOF NOS OF NODAL DOF

  IF ( DOF_VARY ) THEN ! VARIABLE DOF PER NODE

! EQ_IJ = DOF_PT_SUM (I-1) + J, J <= DOF_PT(I), I > 0
    INDEX = 0 ! DEFAULT MAXIMUM ARRAY
    EQ_PT = DOF_PT (I_PT)      ! DOF AT THIS NODE
    IF ( EQ_PT < 1 ) RETURN    ! GEOMETRIC ONLY
      DO J = 1, EQ_PT
        INDEX (J) = DOF_PT_SUM (I_PT-1) + J
      END DO ! ON J

  ELSE ! ALL NODES HAVE N_G_DOF

! INDEX (J) = N_G_DOF * (I_PT - 1) + J
    INDEX = (/ (N_G_DOF * (I_PT - 1) + J, J = 1, N_G_DOF) /)
  END IF
END FUNCTION GET_INDEX_AT_PT

FUNCTION GET_ELEM_INDEX (LT_N, ELEM_NODES) RESULT(INDEX)
! * * * * * * * * * * * * * * * * * * * * * * * * * * *
!   DETERMINE DEGREES OF FREEDOM NUMBERS OF ELEMENT
! * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants ! for N_G_DOF, DOF_PT DOF_PT_SUM
 IMPLICIT NONE
 INTEGER, INTENT(IN) :: LT_N, ELEM_NODES (LT_N)
 INTEGER             :: INDEX (LT_N * N_G_DOF)   ! OUT
 INTEGER :: EQ_ELEM, EQ_PT, EQ_SYS, IG, K, SYS_K ! LOOPS
! INTEGER :: DOF_PT (0:MAX_NP), DOF_PT_SUM (-1:MAX_NP)

! DOF_PT     = NUMBER OF DOF AT EACH NODE:  (0:MAX_NP)
! DOF_PT_SUM = TOTAL NUMBER OF DOF, INCLUDING THIS NODE: (-1:MAX_NP)
! N_G_DOF    = NUMBER OF GENERAL PARAMETERS (DOF) PER NODE
! LT_N       = NUMBER OF NODES PER ELEMENT
! ELEM_NODES = NODAL INCIDENCES OF THE ELEMENT
! EQ_ELEM    = LOCAL  EQUATION NUMBER
! EQ_SYS     = SYSTEM EQUATION NUMBER
! INDEX      = SYSTEM DOF NUMBERS OF ELEMENT DOF NUMBERS

  IF ( DOF_VARY ) THEN ! VARIABLE DOF PER ELEMENT NODE

! EQ_IJ = DOF_PT_SUM (I-1) + J, J <= DOF_PT(I), I > 0
    INDEX = 0 ! DEFAULT MAXIMUM ARRAY
    EQ_ELEM = 0
    DO K = 1, LT_N               ! LOOP OVER NODES OF ELEMENT
      SYS_K = ELEM_NODES (K)     ! SYSTEM NODE NUMBER
      IF ( SYS_K < 1 ) CYCLE     ! TO A REAL NODE
      EQ_PT = DOF_PT (SYS_K)     ! DOF AT THIS NODE
      IF ( EQ_PT < 1 ) CYCLE     ! TO A NON_GEOMETRIC NODE
      DO IG = 1, EQ_PT           ! LOOP OVER GENERALIZED DOF
        EQ_ELEM = EQ_ELEM + 1                     ! LOCAL  EQ
        EQ_SYS  = DOF_PT_SUM (SYS_K - 1) + IG     ! SYSTEM EQ
        INDEX (EQ_ELEM) = EQ_SYS
      END DO ! OVER DOF
    END DO ! OVER LOCAL NODES

  ELSE ! ALL NODES HAVE N_G_DOF

! INDEX (N_G_DOF*(K-1)+IG) = N_G_DOF*(ELEM_NODES(K)-1) + IG
    DO K = 1, LT_N               ! LOOP OVER NODES OF ELEMENT
      SYS_K = ELEM_NODES (K)     ! SYSTEM NODE NUMBER
      DO IG = 1, N_G_DOF         ! LOOP OVER GENERALIZED DOF
        EQ_ELEM = IG + N_G_DOF * (K     - 1)      ! LOCAL  EQ
        EQ_SYS  = IG + N_G_DOF * (SYS_K - 1)      ! SYSTEM EQ
        IF ( SYS_K > 0 ) THEN           ! VALID NODE
          INDEX (EQ_ELEM) = EQ_SYS
        ELSE                            ! ALLOW MISSING NODE
          INDEX (EQ_ELEM) = 0
        END IF ! MISSING NODE
      END DO ! OVER DOF
    END DO ! OVER LOCAL NODES

  END IF ! VARIABLE DOF
END FUNCTION GET_ELEM_INDEX

!FUNCTION GET_ELEM_INDEX (LT_N, ELEM_NODES) RESULT(INDEX)
!! * * * * * * * * * * * * * * * * * * * * * * * * * *
!! DETERMINE DEGREES OF FREEDOM NUMBERS OF ELEMENT
!! * * * * * * * * * * * * * * * * * * * * * * * * * *
!Use System_Constants    ! for N_G_DOF
! IMPLICIT NONE
! INTEGER, INTENT(IN) :: LT_N, ELEM_NODES (LT_N)
! INTEGER             :: INDEX (LT_N * N_G_DOF)
! INTEGER             :: IDOF, IELM, IG, K, CONST
!
!! LT_N       = NUMBER OF NODES PER ELEMENT
!! N_G_DOF    = NUMBER OF PARAMETERS (DOF) PER NODE
!! ELEM_NODES = NODAL INCIDENCES OF THE ELEMENT
!! INDEX      = SYSTEM DOF NOS OF ELEMENT DOF
!
!  DO K = 1, LT_N           !  LOOP OVER NODES OF ELEMENT
!    IDOF = - N_G_DOF
!    IF ( ELEM_NODES (K) > 0) IDOF = IDOF + N_G_DOF * ELEM_NODES (K)
!    CONST = N_G_DOF * (K - 1)
!!    LOOP OVER GENERALIZED DEGREES OF FREEDOM
!    DO IG = 1, N_G_DOF
!      IELM = CONST + IG
!!     INDEX (N_G_DOF * (K-1)+IG) = N_G_DOF * (ELEM_NODES(K)-1) + IG
!      INDEX (IELM) = IDOF + IG
!    END DO
!  END DO
!END FUNCTION GET_ELEM_INDEX
                           
FUNCTION GET_ELEM_NODES (L_ID, LT_N, NODES) RESULT (ELEM_NODES)
! * * * * * * * * * * * * * * * * * * * * * * * *
!      EXTRACT NODES ASSOCIATED WITH ELEMENT L_ID
! * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 INTEGER, INTENT(IN) :: L_ID, LT_N
 INTEGER, INTENT(IN) :: NODES      (L_S_TOT, NOD_PER_EL)
 INTEGER             :: ELEM_NODES (LT_N)

! ELEM_NODES = THE NOD_PER_EL NODAL INCIDENCES OF THE ELEMENT
! L_ID       = ELEMENT NUMBER
! LT_N       = NUMBER OF NODES FOR THIS ELEMENT TYPE
! L_S_TOT    = N_ELEMS + N_MIXED + N_F_SEG, TOTAL GENERALIZED ELEMENTS
! N_ELEMS    = NUMBER OF ELEMENTS IN SYSTEM
! NODES      = NODAL INCIDENCES OF ALL ELEMENTS
! NOD_PER_EL = NUMBER OF NODES PER ELEMENT

  ELEM_NODES (1:LT_N) = NODES (L_ID, 1:LT_N)
  if ( minval(elem_nodes) < 0      .or.        &
       maxval(elem_nodes) > MAX_NP)            & 
       stop 'invalid elem_nodes, GET_ELEM_NODES' 
END FUNCTION GET_ELEM_NODES
                           
SUBROUTINE GET_ELEM_QUADRATURES 
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!  RECOVER QUADRATURE DATA BASED ON ELEMENT SHAPE
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
Use Elem_Type_Data  ! for: 
  ! LT_PARM, LT_QP, PT (LT_PARM, LT_QP), WT (LT_QP)  
 IMPLICIT NONE
 INTEGER :: N_GP, NIP ! work, 1D rule Number of Integration Pts

! LT_SHAP  = ELEMENT SHAPE CODE
! LT_QP    = NUMBER OF QUADRATURE POINTS
! LT_PARM  = DIMENSION OF SPACE
! PT       = QUADRATURE COORDINATES
! WT       = QUADRATURE WEIGHTS

  IF ( .NOT. TYPE_GAUS_ALLOC ) STOP 'ALLOC ERROR, GET_ELEM_QUADRATURES'
  N_GP = LT_QP
  SELECT CASE ( LT_SHAP )
    CASE (0:1) !    LINE ELEMENT
      CALL GAUSS_1D (N_GP, N_GP, PT, WT)
    CASE (2) !    TRIANGULAR (UNIT COORDINATES)
      CALL DUNAVANT_UNIT_TRIANGLE_RULE (N_GP, PT, WT)
    CASE (3) !    QUADRILATERAL 
      SELECT CASE ( LT_QP ) ! integer squared only
        CASE (0:1) ; NIP = 1
        CASE (4)   ; NIP = 2
        CASE (9)   ; NIP = 3
        CASE (16)  ; NIP = 4
        CASE DEFAULT
          NIP = SQRT (FLOAT (N_GP) ) + 0.1   !9 ceiling
      END SELECT
      CALL GAUSS_2D (NIP, N_GP, PT, WT)
    CASE (4) !    HEXAHEDRA
      SELECT CASE ( LT_QP ) ! integer cubed only
        CASE (0:1) ; NIP = 1
        CASE (8)   ; NIP = 2
        CASE (27)  ; NIP = 3
        CASE (64)  ; NIP = 4
        CASE DEFAULT
          NIP = (FLOAT (N_GP) ) ** (1. / 3.) + 0.1   !9 ceiling
      END SELECT
      CALL GAUSS_3D (NIP, N_GP, PT, WT)
    CASE (5) !    TETRAHEDRA (UNIT COORDINATES)
      CALL KEAST_UNIT_TET_RULE (N_GP, PT, WT)
    CASE (6) !    WEDGE
      STOP 'WEDGE NOT IN GET_ELEM_QUADRATURES'
    !b CASE (7) !    USER SUPPLIED QUADRATURE SUBROUTINE
    CASE DEFAULT !    UNSUPPORTED
      STOP 'INVALID ELEMENT SHAPE OPTION, GET_ELEM_QUADRATURES'
  END SELECT

  IF ( LT_QP == N_GP) RETURN
    WRITE (6, * ) 'LT_SHAP, LT_QP, NIP, N_GP', LT_SHAP, &
    LT_QP, NIP, N_GP, ' FATAL ERROR, GET_ELEM_QUADRATURES'
    STOP 'FATAL ERROR, GET_ELEM_QUADRATURES'
END SUBROUTINE GET_ELEM_QUADRATURES
                                   
SUBROUTINE GET_NUMB_QUADRATURES (PT_N, WT_N, N_GP)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!     RECOVER QUADRATURE BASED ON NUMBER OF POINTS    
!        SEE: FUNCTION GET_NUMB_INERTIA_PTS
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
Use Elem_Type_Data ! for LT_PARM, LT_SHAP
 IMPLICIT NONE
 INTEGER,  INTENT(IN)  :: N_GP
 REAL(DP), INTENT(OUT) :: PT_N (LT_PARM, N_GP), WT_N (N_GP)
 INTEGER  :: NIP
! LT_SHAP = ELEMENT SHAPE CODE
! N_GP    = NUMBER OF QUADRATURE POINTS
! LT_PARM = DIMENSION OF SPACE
! PT_N    = QUADRATURE COORDINATES
! WT_N    = QUADRATURE WEIGHTS
!
  SELECT CASE ( LT_SHAP )
    CASE (0:1) !    LINE ELEMENT
      CALL GAUSS_1D (N_GP, N_GP, PT_N, WT_N)
    CASE (2) !    TRIANGULAR (UNIT COORDINATES)
      CALL DUNAVANT_UNIT_TRIANGLE_RULE (N_GP, PT_N, WT_N)
    CASE (3) !    QUADRILATERAL
      NIP = SQRT (FLOAT (N_GP) ) + 0.1
      CALL GAUSS_2D (NIP, N_GP, PT_N, WT_N)
    CASE (4) !    HEXAHEDRA
      NIP = (FLOAT (N_GP) ) ** (1. / 3.) + 0.1
      CALL GAUSS_3D (NIP, N_GP, PT_N, WT_N)
    CASE (5) !    TETRAHEDRA (UNIT COORDINATES)
      CALL KEAST_UNIT_TET_RULE (N_GP, PT_N, WT_N)
    CASE (6) !    WEDGE
      STOP 'WEDGE NOT IN GET_NUMB_QUADRATURES'
    !b CASE (7) !    USER SUPPLIED QUADRATURE SUBROUTINE
    CASE DEFAULT !    UNSUPPORTED
      STOP 'INVALID ELEMENT SHAPE OPTION, GET_NUMB_QUADRATURES'
  END SELECT
END SUBROUTINE GET_NUMB_QUADRATURES
                                   
SUBROUTINE GET_ONE_PT_RULE (ONE_PT, ONE_WT)   
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!     RECOVER ONE POINT QUADRATURE RULE
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
Use Elem_Type_Data ! for LT_PARM, LT_SHAP
 IMPLICIT NONE
 REAL(DP), INTENT(OUT) :: ONE_PT (LT_PARM), ONE_WT
 INTEGER,  PARAMETER   :: N_GP = 1
 REAL(DP) :: PT_N (LT_PARM, N_GP), WT_N (N_GP)
 
! LT_SHAP = ELEMENT SHAPE CODE
! N_GP    = NUMBER OF QUADRATURE POINTS
! LT_PARM = DIMENSION OF SPACE
! PT_N    = QUADRATURE COORDINATES
! WT_N    = QUADRATURE WEIGHTS
!
  SELECT CASE ( LT_SHAP )
    CASE (0:1) !    LINE ELEMENT
      CALL GAUSS_1D (N_GP, N_GP, PT_N, WT_N)
    CASE (2) !    TRIANGULAR (UNIT COORDINATES)
      CALL DUNAVANT_UNIT_TRIANGLE_RULE (N_GP, PT_N, WT_N)
    CASE (3) !    QUADRILATERAL
      CALL GAUSS_2D (N_GP, N_GP, PT_N, WT_N)
    CASE (4) !    HEXAHEDRA
      CALL GAUSS_3D (N_GP, N_GP, PT_N, WT_N)
    CASE (5) !    TETRAHEDRA (UNIT COORDINATES)
      CALL KEAST_UNIT_TET_RULE (N_GP, PT_N, WT_N)
    CASE (6) !    WEDGE
      STOP 'WEDGE NOT IN GET_ONE_PT_RULE'
    !b CASE (7) !    USER SUPPLIED QUADRATURE SUBROUTINE
    CASE DEFAULT !    UNSUPPORTED
      STOP 'INVALID ELEMENT SHAPE OPTION, GET_ONE_PT_RULE'
  END SELECT
  ONE_PT (:) = PT_N (:, 1) ; ONE_WT = WT_N (1) ! copy
END SUBROUTINE GET_ONE_PT_RULE
                                   
!b FUNCTION GET_NUMB_INERTIA_PTS ()  RESULT (N_GP)
!b ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!b !   RECOVER NUMBER OF POINTS TO INTEGRATE INERTIA TENSOR
!b ! * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!b Use Precision_Module
!b Use Elem_Type_Data
 !b IMPLICIT NONE
 !b INTEGER :: N_GP
!b 
!b ! LT_SHAP = ELEMENT SHAPE CODE
!b ! N_GP    = NUMBER OF QUADRATURE POINTS
!b ! LT_PARM = DIMENSION OF SPACE
!b !
  !b SELECT CASE ( LT_SHAP )
    !b CASE (0:1) ; N_GP = 2  !    LINE ELEMENT
    !b CASE (2)   ; N_GP = 3  !    TRIANGULAR ELEMENT
    !b CASE (3)   ; N_GP = 4  !    QUADRILATERAL ELEMENT
    !b CASE (4)   ; N_GP = 8  !    HEXAHEDRAL ELEMENT
    !b CASE (5)   ; N_GP = 4  !    TETRAHEDRAL ELEMENT
    !b CASE (6)   ; N_GP = 6  !    WEDGE ELEMENT
    !b CASE (7)               !    USER SUPPLIED 
      !b ! CALL USER_INERTIA_PT (N_GP, LT_PARM, PT, WT)
      !b !b STOP 'USER OPTION MISSING, GET_NUMB_INERTIA_PTS'
    !b CASE DEFAULT 
      !b STOP 'INVALID OPTION, GET_NUMB_INERTIA_PTS'
  !b END SELECT
!b END FUNCTION GET_NUMB_INERTIA_PTS
                                   
SUBROUTINE GET_ELEM_TYPE_DATA (LT)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!   GET DATA BASED ON ELEMENT TYPE (SEE CONTROL_INPUT )
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data       ! for LAST_LT, N_L_TYPE
 IMPLICIT NONE
 INTEGER, INTENT(IN)  :: LT

! LAST_LT  = LAST ACTIVE ELEMENT TYPE NUMBER
! LT       = ELEMENT TYPE NUMBER
! LT_N     = NUMBER OF NODES PER ELEMENT
! LT_QP    = NUMBER OF QUADRATURE POINTS
! LT_GEOM  = NUMBER OF ELEMENT GEOMERTY NODES
! LT_PARM  = NUMBER OF PARAMETRIC SPACES FOR ELEMENT
! LT_FACES = NUMBER OF (FACING) NEIGHBORS
! LT_SONS  = NUMBER OF SONS IF DIVIDED
! LT_SHAP  = ELEMENT SHAPE FLAG NUMBER
! LT_USER  = APPLICATION DEPENDENT OPTIONAL USER ITEM
! N_L_TYPE = NUMBER OF ELEMENT TYPES

  IF ( LT < 1 .OR. LT > N_L_TYPE ) STOP &
       'ELEMENT TYPE IMPOSSIBLE, GET_ELEM_TYPE_DATA'

!      GET CONTROLS FOR THIS TYPE
  LT_N     = LT_DATA (1, LT)
  LT_QP    = LT_DATA (2, LT)
  LT_GEOM  = LT_DATA (3, LT)
  LT_PARM  = LT_DATA (4, LT)
  LT_SHAP  = LT_DATA (5, LT)
  LT_FACES = LT_DATA (6, LT)
  LT_SONS  = LT_DATA (7, LT)
  LT_USER  = LT_DATA (8, LT)
  LT_FREE  = LT_N * N_G_DOF
END SUBROUTINE GET_ELEM_TYPE_DATA
                                 
SUBROUTINE LIST_ELEM_TYPE_DATA (LT)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!   LIST DATA BASED ON ELEMENT TYPE (SEE CONTROL_INPUT )
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data       ! for LAST_LT, N_L_TYPE
 IMPLICIT NONE
 INTEGER, INTENT(IN)  :: LT

! LAST_LT  = LAST ACTIVE ELEMENT TYPE NUMBER
! LT       = ELEMENT TYPE NUMBER
! LT_N     = NUMBER OF NODES PER ELEMENT
! LT_QP    = NUMBER OF QUADRATURE POINTS
! LT_GEOM  = NUMBER OF ELEMENT GEOMERTY NODES
! LT_PARM  = NUMBER OF PARAMETRIC SPACES FOR ELEMENT
! LT_FACES = NUMBER OF (FACING) NEIGHBORS
! LT_SONS  = NUMBER OF SONS IF DIVIDED
! LT_SHAP  = ELEMENT SHAPE FLAG NUMBER
! LT_USER  = APPLICATION DEPENDENT OPTIONAL USER ITEM
! N_L_TYPE = NUMBER OF ELEMENT TYPES

  IF ( LT < 1 .OR. LT > N_L_TYPE ) THEN
    PRINT *, 'WARNING: ELEMENT TYPE ', &
    LT, ' IS IMPOSSIBLE, LIST_ELEM_TYPE_DATA'
    N_WARN  = N_WARN + 1 ! INCREMENT WARNING
  END IF
  LT_FREE = LT_N * N_G_DOF

!     LIST CONTROLS FOR THIS TYPE
  PRINT *, ' ' ;
  PRINT *, '*** LIST OF DATA FOR ELEMENT TYPE ', LT, ' ***'
  PRINT *, 'LT_N     = ', LT_N      
  PRINT *, 'LT_QP    = ', LT_QP     
  PRINT *, 'LT_GEOM  = ', LT_GEOM   
  PRINT *, 'LT_PARM  = ', LT_PARM   
  PRINT *, 'LT_SHAP  = ', LT_SHAP   
  PRINT *, 'LT_FACES = ', LT_FACES  
  PRINT *, 'LT_SONS  = ', LT_SONS   
  PRINT *, 'LT_USER  = ', LT_USER   
  PRINT *, 'LT_FREE  = ', LT_FREE  
END SUBROUTINE LIST_ELEM_TYPE_DATA
                                 
SUBROUTINE SET_ELEM_TYPE_INFO (LT)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!   SET DATA BASED ON ELEMENT TYPE (SEE CONTROL_INPUT )
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data  ! for: MAX_DAT, N_L_TYPE,                         &
  ! LT_DATA (MAX_DAT, N_L_TYPE),                                      &
  ! 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)  
 IMPLICIT NONE
 INTEGER, INTENT(IN)  :: LT

! LAST_LT  = LAST ACTIVE ELEMENT TYPE NUMBER
! LT       = ELEMENT TYPE NUMBER
! LT_DATA  = ARRAY OF CONSTANTS FOR EACH ELEMENT TYPE
!    1, LT_N     = NUMBER OF NODES PER ELEMENT
!    2, LT_QP    = NUMBER OF QUADRATURE POINTS
!    3, LT_GEOM  = NUMBER OF ELEMENT GEOMERTY NODES
!    4, LT_PARM  = NUMBER OF PARAMETRIC SPACES FOR ELEMENT
!    5, LT_FACES = NUMBER OF (FACING) NEIGHBORS
!    6, LT_SONS  = NUMBER OF SONS IF DIVIDED
!    7, LT_SHAP  = ELEMENT SHAPE FLAG NUMBER
!    8, LT_USER  = APPLICATION DEPENDENT OPTIONAL USER ITEM
! LT_FREE  = LT_N * N_G_DOF
! N_L_TYPE = NUMBER OF ELEMENT TYPES

  IF ( LT < 1 .OR. LT > N_L_TYPE ) STOP &
       'ELEMENT TYPE WRONG, SET_ELEM_TYPE_INFO'

!      GET CONTROLS FOR THIS TYPE:
!      LT_N,    LT_QP,    LT_GEOM, LT_PARM, 
!      LT_SHAP, LT_FACES, LT_SONS, LT_USER, LT_FREE
  CALL GET_ELEM_TYPE_DATA (LT)

!      RESET LAST, ALLOCATE ARRAY FOR AN ELEMENT TYPE
  IF ( LT /= LAST_LT )  THEN
    LAST_LT = LT

    IF ( TYPE_APLY_ALLOC ) CALL DEALLOCATE_TYPE_APPLICATION
    CALL   ALLOCATE_TYPE_APPLICATION
    IF ( TYPE_NTRP_ALLOC ) CALL DEALLOCATE_TYPE_INTERPOLATIONS
    CALL ALLOCATE_TYPE_INTERPOLATIONS

 !-->     GET QUADRATURE RULE FOR ELEMENT TYPE AND SHAPE
    IF ( LT_QP > 0 ) THEN
!b print *, 'SET_ELEM_TYPE_INFO LT, LT_QP ', LT, LT_QP !b
      CALL GET_ELEM_QUADRATURES
      CALL FILL_TYPE_INTERPOLATIONS
    END IF ! quadratures used
  END IF ! new element type
END SUBROUTINE SET_ELEM_TYPE_INFO
                                 
SUBROUTINE GET_REACTIONS (DD)
! * * * * * * * * * * * * * * * * * * * * * * * * *
!-->     COMPUTE REACTIONS IN FROM STORED SYSTEM ROWS
! * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 REAL(DP), INTENT(IN) :: DD (N_D_FRE)
!          Automatic Arrays
 REAL(DP) :: TOTAL (N_G_DOF), SUM_POS (N_G_DOF), SUM_NEG (N_G_DOF)
 INTEGER  :: IE, NODE, IG, J1, J2, J
 REAL(DP) :: R, CI, SIJ

! DD         = COMPUTED SOLUTION VECTOR, SS*DD=CC
! MODE       = STORAGE FLAG, 0=SYMMETRIC, OR 1=UNSYMMETRIC SKYLINE
! NOD_PER_EL = DOF NUMBER OF ESSENTIAL B.C.
! N_G_DOF    = NUMBER OF DOF PER NODE
! N_REACT    = SEQUENTIAL UNIT TO STORE REACTION DATA
! N_D_FRE    = TOTAL NUMBER OF EQUATIONS
! TOTAL      = SUM OF REACTIONS FOR N_G_DOF PARAMETERS

  IF ( N_REACT > 0 ) THEN
    REWIND N_REACT
    WRITE (6, 5)
    5 FORMAT ( /, '*** REACTION RECOVERY ***', /,  &
               '   NODE, PARAMETER,     REACTION, EQUATION')
    TOTAL = 0.d0 ; SUM_POS = 0.d0 ; SUM_NEG = 0.d0
  ELSE
    PRINT *,'WARNING: NO REACTION FILE, BUT GET_REACTIONS WAS CALLED'
    N_WARN = N_WARN + 1
    RETURN
  END IF

!    READ NODE, PARAMETER, RANGE OF NON-ZERO TERMS
10 READ (N_REACT, END = 30) NODE, IG, J1, J2
   R = 0.d0
   DO J = J1, J2
     READ (N_REACT) SIJ
     R = R + SIJ * DD (J)
   END DO
   READ (N_REACT) CI
   R = R - CI
   IE = N_G_DOF * (NODE-1) + IG
   WRITE (6, 5010) NODE, DOF_NAME (IG), R, IE
   5010 FORMAT ( I7, ', ', A12, 1X, 1PE11.4, I6 )
   !b TOTAL (IG) = TOTAL (IG) + R
   IF ( R >= 0.d0 ) THEN   !b
     SUM_POS (IG) = SUM_POS (IG) + R   !b
   ELSE   !b
     SUM_NEG (IG) = SUM_NEG (IG) + R   !b
   END IF   !b
   GO TO 10 !b replace with iostat in infinite loop !xxx

30 WRITE (6, 5020)
   5020 FORMAT ('REACTION RESULTANTS', /, &
                'PARAMETER,    SUM          POSITIVE     NEGATIVE')
   TOTAL = SUM_POS + SUM_NEG !b
   DO J = 1, N_G_DOF
     WRITE (6, 5030) DOF_NAME (J), TOTAL (J), SUM_POS (J), SUM_NEG (J)
     5030 FORMAT ( A12, 1PE11.4, 2X, 1PE11.4, 2X, 1PE11.4 )
   END DO
END SUBROUTINE GET_REACTIONS
                            
!SUBROUTINE INDXEL (NOD_PER_EL, N_EL_FRE, N_G_DOF, ELEM_NODES, INDEX)
!! * * * * * * * * * * * * * * * * * * * * * * * * * *
!! DETERMINE DEGREES OF FREEDOM NUMBERS OF ELEMENT
!! * * * * * * * * * * * * * * * * * * * * * * * * * *
!Use Precision_Module
! IMPLICIT NONE
! INTEGER, INTENT(IN)  :: NOD_PER_EL, N_EL_FRE, N_G_DOF
! INTEGER, INTENT(IN)  :: ELEM_NODES (NOD_PER_EL)
! INTEGER, INTENT(OUT) :: INDEX (N_EL_FRE)
! INTEGER :: CONST, IDOF, IG, IELM, K
!
!! NOD_PER_EL = NUMBER OF NODES PER ELEMENT
!! N_G_DOF    = NUMBER OF PARAMETERS (DOF) PER NODE
!! N_EL_FRE   = NOD_PER_EL*N_G_DOF = NUMBER OF DOF PER ELEMENT
!! ELEM_NODES = NODAL INCIDENCES OF THE ELEMENT
!! INDEX      = SYSTEM DOF NOS OF ELEMENT DOF
!
!  !b IF ( DOF_VARY ) STOP 'INDXEL: INVALID FOR DOF_VARY'
!
!!  LOOP OVER NODES OF ELEMENT
!  DO K = 1, NOD_PER_EL
!    IDOF = - N_G_DOF
!    IF ( ELEM_NODES (K) > 0) IDOF = IDOF + N_G_DOF * ELEM_NODES (K)
!    CONST = N_G_DOF * (K - 1)
!!    LOOP OVER GENERALIZED DEGREES OF FREEDOM
!    DO IG = 1, N_G_DOF
!      IELM = CONST + IG
!!     INDEX (N_G_DOF * (K-1)+IG) = N_G_DOF * (ELEM_NODES(K)-1) + IG
!      INDEX (IELM) = IDOF + IG
!    END DO
!  END DO
!END SUBROUTINE INDXEL
!                     
!SUBROUTINE INDXPT (I_PT, N_G_DOF, INDEX)
!! * * * * * * * * * * * * * * * * * * * * * * * * * *
!! DETERMINE DEGREES OF FREEDOM NUMBERS AT A NODE
!! * * * * * * * * * * * * * * * * * * * * * * * * * *
! IMPLICIT NONE
! INTEGER, INTENT(IN)  :: I_PT, N_G_DOF
! INTEGER, INTENT(OUT) :: INDEX (N_G_DOF)
! INTEGER :: J, CONST
!
!! I_PT     = SYSTEM NODE NUMBER
!! N_G_DOF = NUMBER OF PARAMETERS (DOF) PER NODE
!! INDEX   = SYSTEM DOF NOS OF NODAL DOF
!
!  !b IF ( DOF_VARY ) STOP 'INDXPT: INVALID FOR DOF_VARY'
!
!  CONST = N_G_DOF * (I_PT - 1)
!  DO J = 1, N_G_DOF
!!   INDEX (J) = N_G_DOF*(I_PT - 1) + J
!    INDEX (J) = CONST + J
!  END DO
!END SUBROUTINE INDXPT
                     
SUBROUTINE ELEM_NODE_PROPERTIES (LT_N, PRT_L_PT, L_PT_PROP, &
                                 ELEM_NODES)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! EXTRACT FLOATING POINT PROPERTIES AT NODAL POINTS
!                  OF AN ELEMENT
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Sys_Properties_Data
 IMPLICIT NONE
 INTEGER,  INTENT(IN)  :: LT_N
 INTEGER,  INTENT(IN)  :: ELEM_NODES (LT_N)
 REAL(DP), INTENT(OUT) :: PRT_L_PT  (NOD_PER_EL, N_NP_FLO)
 INTEGER,  INTENT(OUT) :: L_PT_PROP (NOD_PER_EL, N_NP_FIX)
 INTEGER :: I, IROW, J

! FLT_NP     = FLOATING POINT PROP ARRAY OF SYSTEM NODES
! PRT_L_PT   = FLOATING POINT PROP ARRAY OF ELEMENT NODES
! NP_FIX     = INTEGER PROPERTY ARRAY OF SYSTEM NODES
! L_PT_PROP  = INTEGER PROPERTY ARRAY OF ELEMENT NODES
! ELEM_NODES = ELEMENT INCIDENCES ARRAY OF THE ELEMENT
! MAX_NP     = NUMBER OF SYSTEM NODES
! LT_N       = NUMBER OF NODES PER ELEMENT
! N_NP_FIX   = NUMBER OF INTEGER PROPERTIES PER NODE
! N_NP_FLO   = NUMBER OF REAL PROPERTIES PER NODE
! N_HOMO     = 1, IF PROPERTIES ARE SAME AT EACH NODE

  DO I = 1, LT_N
    IROW = ELEM_NODES (I)
!    ALLOW FOR OMITTED NODES          !9 use WHERE
    IF ( IROW > 0 ) THEN
      IF ( N_HOMO == 1) IROW = 1
      IF ( N_NP_FLO > 0 ) THEN
        DO J = 1, N_NP_FLO
          PRT_L_PT (I, J) = FLT_NP (IROW, J)
        END DO
      END IF
      IF ( N_NP_FIX > 0 ) THEN
        DO J = 1, N_NP_FIX
          L_PT_PROP (I, J) = NP_FIX (IROW, J)
        END DO
      END IF
    END IF
  END DO
END SUBROUTINE ELEM_NODE_PROPERTIES
                                   
SUBROUTINE MATERIAL_PROP (NUM, PRT_MAT)  !b not right
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!    EXTRACT REAL MATERIAL PROPERTIES OF MATERIAL
!   NUM FROM MISCELLANEOUS REAL SYSTEM PROPERTIES
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Sys_Properties_Data
 IMPLICIT NONE
 INTEGER,  INTENT(IN)  :: NUM
 REAL(DP), INTENT(OUT) :: PRT_MAT (MISC_FL) 
 INTEGER               :: I, ISTART
! NUM    = MATERIAL NUMBER
! N_MAT   = MAXIMUM ALLOWABLE MATERIAL NUMBER
! N_LP_FLO = NUMBER OF REAL ELEMENT PROPERTIES
! MISC_FL = NO OF MISC REAL SYSTEM PROPERTIES
! FLT_MISC = SYSTEM STORAGE FOR MISC REAL PROP
! PRT_MAT = REAL PROPERTY ARRAY FOR MATERIAL NUM
!9 see if the following is still true
!  PROPERTIES ARE STORED IN FLT_MISC IN ORDER OF MAT NO

if ( num >= 0 ) stop 'MATERIAL_PROP invalid code' !b
  IF ( N_LP_FLO < 1 ) STOP 'BAD N_LP_FLO, MATERIAL_PROP'
  IF ( NUM > N_MAT ) STOP 'DIMENSIONS EXCEEDED, MATERIAL_PROP'
  ISTART = N_LP_FLO * (NUM - 1)
  DO I = 1, N_LP_FLO
    PRT_MAT (I) = FLT_MISC (ISTART + I)
  END DO
END SUBROUTINE MATERIAL_PROP
                            
SUBROUTINE PENLTY (NP_FREE, COEF_C_EQ, C_Pen, S_Pen, WT)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! DEFINE CONSTRAINT PENALTY SQUARE AND COLUMN MATRICES
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 INTEGER, INTENT(IN)   :: NP_FREE
 REAL(DP), INTENT(IN)  :: COEF_C_EQ  (NP_FREE), WT
 REAL(DP), INTENT(OUT) :: C_Pen (NP_FREE), S_Pen (NP_FREE, NP_FREE)

! NP_FREE         = NO DOF IN CONSTRAINT EQUATION
! COEF_C_EQ (I) = CONSTR EQ COEFFICIENT I+1
! C_Pen         = CONSTRAINT COLUMN MATRIX
! S_Pen         = CONSTRAINT SQUARE MATRIX
! WT            = PENALTY WEIGHT FACTOR

!  INITIAL CALCULATIONS
  C_Pen (1) = 1.0
  TEMP = COEF_C_EQ (NP_FREE)
  IF ( NP_FREE > 1 ) THEN
    DO I = 2, NP_FREE
      C_Pen (I) = COEF_C_EQ (I - 1)
    END DO
  END IF
!  CALCULATE LEAST SQ CONSTRAINT FORMS
  DO I = 1, NP_FREE
    DO J = 1, NP_FREE
      S_Pen (J, I) = WT * C_Pen (I) * C_Pen (J)
    END DO
  END DO
  DO I = 1, NP_FREE
    C_Pen (I) = C_Pen (I) * TEMP * WT
  END DO
END SUBROUTINE PENLTY
                     
SUBROUTINE POST_PROCESS_APPLICATION (NODES, DD, DD_OLD, ITER, X)
!b SUBROUTINE POST_PROCESS_APPLICATION (NODES, DD, ITER, X)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!       ELEMENT LEVEL POST-SOLUTION CALCULATIONS
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data       ! for PT (LT_PARM, LT_QP), WT (LT_QP),
                         ! G (LT_GEOM, LT_QP), H (LT_N   ),
                         ! DLG (LT_PARM, LT_GEOM, LT_QP),
                         ! DLH (LT_PARM, LT_N   , LT_QP),
                         ! S (LT_FREE, LT_FREE), C (LT_FREE),
                         ! ELEM_NODES (LT_N)
Use Sys_Properties_Data  ! for FLT_MISC, MISC_FIX
Use Select_Source        
Use Interface_Header
 IMPLICIT NONE
! Material option for future uses. Add to arguments, etc.
!    ALWAYS USED
 REAL(DP), INTENT(INOUT) :: DD (N_D_FRE), DD_OLD (N_D_FRE)   
 REAL(DP), INTENT(IN) :: X (MAX_NP, N_SPACE)
 INTEGER,  INTENT(IN) :: NODES (L_S_TOT, NOD_PER_EL), ITER

!    Automatic Arrays
 !b REAL(DP) :: EL_PROP (N_LP_FLO), PRT_MAT (MISC_FL) 
 !b INTEGER  :: L_PROP (N_LP_FIX), L_PT_PROP (N_NP_FIX)  !b ??

 !B REAL(DP) :: FIT_LT (LT_N, SCP_FIT) !b
 REAL(DP) :: PRT_MAT (MISC_FL) 
 INTEGER  :: L_PT_PROP (N_NP_FIX)  !b ??
 INTEGER  :: IE, LT
 INTEGER, SAVE :: NOTE = 0
!          Allocatable Arrays
 REAL(DP), ALLOCATABLE :: PRT_L_PT (:, :)

!                 VARIABLES:
! COORD      = COORDINATES OF ALL NODES ON ELEMENT
! D          = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT
! DD         = ARRAY OF SYSTEM DEGREES OF FREEDOM
! DGH        = GLOBAL DERIVATIVES INTERPOLATION FUNCTIONS
! DLG        = LOCAL DERIVATIVES GEOMETRIC INTERPOLATION
! DLH        = LOCAL DERIVATIVES INTERPOLATION FUNCTIONS
! E          = CONSTITUTIVE MATRIX
! EB         = PRODUCT OF E*B
! FLT_MISC   = SYSTEM STORAGE OF FLOATING PT MISC PROP
! FLUX_LT    = FLUX AT ELEMENT NODES FROM A SCP: (SCP_FIT, LT_N)
! G          = GEOMETRIC INTERPOLATION FUNCTIONS
! H          = SOLUTION INTERPOLATION FUNCTIONS
! H_INTG     = INTEGRAL OF INTERPOLATION FUNCTIONS
! INDEX      = SYSTEM DOF NOS ASSOCIATED WITH ELEMENT
! ITER       = CURRENT ITERATION NUMBER
! ELEM_NODES = THE NOD_PER_EL INCIDENCES OF THE ELEMENT
! L_PT_PROP  = INTEGER PROPERTIES AT EACH ELEMENT NODE
! MAX_NP     = TOTAL NUMBER OF NODES
! MISC_FIX   = MISCELLANEOUS INTEGER SYSTEM PROPERTIES
! NOD_PER_EL = NUMBER OF NODES PER ELEMENT
! N_D_FRE    = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM
! N_ELEMS    = NUMBER OF ELEMENTS IN SYSTEM
! N_EL_FRE   = NUMBER OF DEGREES OF FREEDOM PER ELEMENT
! N_G_DOF    = NUMBERS PARAMETERS PER NODE
! N_GEOM     = NUMBER OF GEOMETRY NODES
! N_ITER     = MAX NUMBER OF ITERATIONS
! N_MAT      = NUMBER OF MATERIAL TYPES
! NODES      = ELEMENT INCIDENCES OF ALL ELEMENTS
! N_PARM     = DIMENSION OF PARAMETRIC SPACE
! N_QP       = NUMBER OF QUADRATURE POINTS
! N_R_B      = NUMBER OF ROWS IN B AND E MATRICES
! N_SPACE    = DIMENSION OF SPACE
! N_FILE1    = UNIT FOR POST SOLUTION MATRICES STORAGE
! N_FILE2    = UNIT FOR POST SOLUTION MATRICES STORAGE
! N_FILE3-5  = OPTIONAL UNITS FOR USER (USED WHEN > 0)
! PRT_L_PT   = REAL PROPERTIES AT ELEMENT NODES
! PRT_MAT    = REAL ELEM PROPERTIES BASED ON MATERIAL NUMBER
! PT         = QUADRATURE COORDINATES
! S          = ELEMENT SQUARE MATRIX
! SCP_AVERAGES = AVERAGED FLUX ITEMS AT ALL NODES
! SCP_FIT    = NUMBER IF TERMS BEING FIT,   = N_R_B USUALLY
! STRAIN     = STRAIN OR GRADIENT VECTOR
! STRAIN_0   = INITIAL STRAIN OR GRADIENT VECTOR
! STRESS     = STRESS VECTOR
! WT         = QUADRATURE WEIGHTS
! XYZ        = SPACE COORDINATES AT A POINT
                                        
!  N_FILE1 OR 2 MUST BE > 0
  IF ( N_FILE1 ==  0 .AND. N_FILE2 == 0 ) THEN
    STOP 'NO N_FILE1 OR NFILE2 IN POST_PROCESS_APPLICATION'
  ELSE
    IF ( N_FILE1 > 0 ) REWIND (N_FILE1)
    IF ( N_FILE2 > 0 ) REWIND (N_FILE2)
  END IF
  IF ( U_FLUX  > 0 ) THEN 
    REWIND (U_FLUX )  
  ELSE 
    PRINT *, 'NOTE: UNIT U_FLUX NOT AVAILABLE FOR APPLICATION'
  END IF 

  WRITE (N_PRT, "(/,'** BEGIN ELEMENT APPLICATION POST PROCESSING **')")
  LT = 1 ; LAST_LT = 0               ! INITIALIZE ELEMENT TYPES 
  !b STRAIN = 0.d0 ; STRAIN_0 = 0.d0 ; STRESS = 0.d0 ! INITRIALIZE !b

!-->   LOOP OVER ELEMENTS
  DO IE = 1, L_S_TOT ! for elements and any boundary segments
    CALL SET_TO_MIXED_OR_SEG_OR_STD (IE) ! For IS_STD, IS_MIX, IS_SEG
    IF ( IE == (N_ELEMS + 1) ) WRITE (N_PRT, &
       "('** BEGIN APPLICATION BOUNDARY SEGMENTS POST PROCESSING **')")

!-->     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
!      GET CONTROLS FOR THIS TYPE
      CALL SET_ELEM_TYPE_INFO (LT)

!-->        Allocate Arrays for this type
      IF ( ALLOCATED (PRT_L_PT) ) DEALLOCATE ( PRT_L_PT ) !ba
      ALLOCATE ( PRT_L_PT ( LT_N, N_NP_FLO ) )
!b  need the same for integer prop  !b
      IF ( I_BUG > 0 ) PRINT *, ALLOCATED ( PRT_L_PT )

!-->       GET QUADRATURE RULE FOR ELEMENT TYPE AND SHAPE
      IF ( LT_QP > 0 ) CALL GET_ELEM_QUADRATURES 
    END IF ! a new element type

!        INITIALIZE AND CHECK IF ACTIVE
    IF ( DELETED (IE) ) CYCLE ! to skip deleted element

!-->     EXTRACT ELEMENT NODE NUMBERS
      ELEM_NODES = GET_ELEM_NODES (IE, LT_N, NODES)    !9

!-->     EXTRACT NODAL COORDINATES
      CALL ELEM_COORD (LT_N, N_SPACE, X, COORD, ELEM_NODES)

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

!-->     EXTRACT NODAL PARAMETERS OF THE ELEMENT
      D = GET_ELEM_DOF (DD)

!-->     EXTRACT NODAL POINT PROPERTIES (IF ANY)
      IF ( N_NP_FLO > 0) CALL ELEM_NODE_PROPERTIES (LT_N, &
                              PRT_L_PT, L_PT_PROP, ELEM_NODES)

!b !-->     EXTRACT ELEMENT PROPERTIES (IF ANY)
      !b IF ( LP_TEST > 0) CALL ELEM_PROPERTIES (IE, L_PROP, EL_PROP)

!B !-->     EXTRACT ELEMENT NODAL FLUXES AT ELEM_NODES (IF ANY) !b
!B     IF ( SCP_RECORDS_SAVED ) THEN !b
!B       FIT_LT = GET_SCP_ELEM_FIT (SCP_AVERAGES) !b
!B     ELSE ; FIT_LT = 0.d0 ; END IF !b

!-->     EXTRACT MATERIAL PROPERTIES (IF ANY)
      IF ( N_MAT > 0) CALL MATERIAL_PROP (N_MAT, PRT_MAT)
!
!-->     PERFORM PROBLEM DEPENDENT CALCULATIONS AND OUTPUT
!
      !b IF ( POST_EL .AND. IS_ELEMENT ) THEN  !b
      IF ( IE <= N_ELEMS ) THEN  !b
        IF ( DEBUG_POST_EL ) PRINT *,'Debug IS_STD ', IS_STD
        IF ( POST_EL .AND. IS_ELEMENT ) CALL                 &
             SELECT_POST_PROCESS_ELEM (PRT_L_PT, PRT_MAT, &  
                                       L_PT_PROP, ITER, IE, &
                                       DD, DD_OLD)  !b
                                       !b L_PT_PROP, ITER, IE)  !b
      ELSE IF ( POST_MIXED .AND. IS_MIXED_BC ) THEN  !b
        IF ( DEBUG_MIX_SQ ) PRINT *,'Debug IS_MIX ', IS_MIX
        CALL SELECT_POST_PROCESS_MIXED (PRT_L_PT, PRT_MAT,     &  !b
                                        L_PT_PROP, ITER, IS_MIX)  !b
      ELSE  !b
        IF ( NOTE == 0 ) THEN ; NOTE = 1;
          PRINT *,'WARNING, SELECT_POST_PROCESS UNKNOWN OPTION'  !b
          PRINT *, THIS_EL, 'NOT A STANDARD OR MIXED SEGMENT ELEMENT'
          PRINT *,'IS_STD, IS_MIX, IS_SEG', IS_STD, IS_MIX, IS_SEG
          PRINT *,'POST_EL, IS_ELEMENT', POST_EL, IS_ELEMENT
          PRINT *,'POST_MIXED, IS_MIXED_BC', POST_MIXED, IS_MIXED_BC
          N_WARN = N_WARN + 1  !b
        END IF ! Warn only once
      END IF ! Element class  !b
      IF ( IE == N_ELEMS ) WRITE (N_PRT, &
                "('ELEMENT POST PROCESSING COMPLETE.')")
      IF ( IE == L_S_TOT .AND. (N_MIXED > 0 .OR. N_F_SEG > 0) )   &
        WRITE (N_PRT, "('POST PROCESSING COMPLETE.', /)")
    END DO ! over all elements and boundary segments

!-->    POST-PROCESSING COMPLETE
  !b CALL DEALLOCATE_TYPE_INTERPOLATIONS
  IF ( ALLOCATED (PRT_L_PT) ) DEALLOCATE ( PRT_L_PT ) 
END SUBROUTINE POST_PROCESS_APPLICATION

SUBROUTINE POST_PROCESS_GRADS (NODES, DD, ITER)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!    SAVE ELEMENT GRADIENTS AS SCP INPUT RECORDS
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants     ! for L_S_TOT, N_D_FRE, NOD_PER_EL
    ! N_L_TYPE, THIS_EL, U_FLUX
Use Elem_Type_Data       ! for PT (LT_PARM, LT_QP), WT (LT_QP),
                         ! G (LT_GEOM, LT_QP), H (LT_N   ),
                         ! DLG (LT_PARM, LT_GEOM, LT_QP),
                         ! DLH (LT_PARM, LT_N   , LT_QP),
                         ! S (LT_FREE, LT_FREE), C (LT_FREE),
                         ! ELEM_NODES (LT_N), D (LT_FREE)
!b Use Sys_Properties_Data  ! for FLT_MISC, MISC_FIX
Use Interface_Header
 IMPLICIT NONE

 REAL(DP), INTENT(IN) :: DD (N_D_FRE)
 INTEGER,  INTENT(IN) :: NODES (L_S_TOT, NOD_PER_EL), ITER

 INTEGER  :: IE, LT

!                 VARIABLES:
! D          = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT
! DD         = ARRAY OF SYSTEM DEGREES OF FREEDOM
! DGH        = GLOBAL DERIVATIVES INTERPOLATION FUNCTIONS
! E          = CONSTITUTIVE MATRIX
! INDEX      = SYSTEM DOF NOS ASSOCIATED WITH ELEMENT
! ITER       = CURRENT ITERATION NUMBER
! ELEM_NODES = THE NOD_PER_EL INCIDENCES OF THE ELEMENT
! NOD_PER_EL = NUMBER OF NODES PER ELEMENT
! N_D_FRE    = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM
! N_L_TYPE   = NUMBER OF DIFFERENT ELEMENT TYPES USED
! N_ELEMS    = NUMBER OF ELEMENTS IN SYSTEM
! NODES      = ELEMENT INCIDENCES OF ALL ELEMENTS
! N_QP       = NUMBER OF QUADRATURE POINTS
! N_R_B      = NUMBER OF ROWS IN B AND E MATRICES
! U_FLUX     = BINARY UNIT TO STORE GRADIENTS OR FLUXES
                                        
!  U_FLUX  MUST BE > 0
  IF ( U_FLUX  <=  0 ) THEN
    PRINT *, 'NO U_FLUX  IN POST_PROCESS_GRADS'
  ELSE
    REWIND (U_FLUX )
  END IF
                                        
  WRITE (N_PRT, "(/, 'BEGIN ELEMENT FLUX SAVE, ITER =', I4)") ITER
  LT = 1 ! INITIALIZE ELEMENT TYPES 
  !b STRAIN = 0.d0 ; STRAIN_0 = 0.d0 ; STRESS = 0.d0 ! INITRIALIZE !b
  !b LAST_LT = 0 !b
!b call inquire_about_file ("GRAD.BIN") !b
!b call inquire_about_file ("N_FILE1.BIN") !b

!-->   LOOP OVER ELEMENTS
  DO IE = 1, N_ELEMS ! for elements and any boundary segments
    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
!      GET CONTROLS FOR THIS TYPE
      CALL SET_ELEM_TYPE_INFO (LT) !B set or get ??? !b
    END IF ! a new element type

!        INITIALIZE AND CHECK IF ACTIVE
    !b IF ( DELETED (IE) ) CYCLE ! to skip deleted element !XXX

!-->     EXTRACT ELEMENT NODE NUMBERS
      ELEM_NODES = GET_ELEM_NODES (IE, LT_N, NODES)    !9

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

!-->     EXTRACT NODAL PARAMETERS OF THE ELEMENT
      D = GET_ELEM_DOF (DD)

!-->  RECOVER FLUXES, LIST, SAVE FOR SCP
!b print *, 'IE, LT, LT_QP, U_FLUX ', IE, LT, LT_QP, U_FLUX 
!b call LIST_SYSTEM_UNIT_NUMBERS !b
      IF ( USE_EXACT_FLUX ) THEN
        CALL LIST_ELEM_AND_EXACT_FLUXES (U_FLUX, IE) !b
      ELSE
        CALL LIST_ELEM_FLUXES (U_FLUX, IE) !b
      END IF
                              
      IF ( IE == N_ELEMS ) WRITE (N_PRT, &
           "('ELEMENT FLUX SAVE COMPLETE (TO UNIT U_FLUX).')")
    END DO ! over all elements 
END SUBROUTINE POST_PROCESS_GRADS

SUBROUTINE PT_COORD (I_PT, P_SPACE, X, COORD)
! * * * * * * * * * * * * * * * * * * * * * * * * *
! EXTRACT COORDINATES OF POINT NUMBER I_PT
! * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 INTEGER,  INTENT(IN)  :: I_PT, P_SPACE
 REAL(DP), INTENT(IN)  :: X (MAX_NP, P_SPACE)
 REAL(DP), INTENT(OUT) :: COORD (1, P_SPACE)

! X      = SPATIAL COORDINATES OF ALL SYSTEM NODES
! COORD  = SPATIAL COORDINATES OF THE NODE
! MAX_NP = TOTAL NUMBER OF NODES IN SYSTEM
! NOD_PER_EL      = NUMBER OF NODES PER ELEMENT
! P_SPACE = DIMENSION OF THE SPACE
    COORD (1, 1:P_SPACE) = X (I_PT, 1:P_SPACE)
END SUBROUTINE PT_COORD
                     
FUNCTION GET_PT_COORD (I_PT, X) RESULT (COORD_PT)
! * * * * * * * * * * * * * * * * * * * * * * * * *
! EXTRACT COORDINATES OF POINT NUMBER I_PT
! * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 INTEGER,  INTENT(IN) :: I_PT
 REAL(DP), INTENT(IN) :: X (MAX_NP, N_SPACE)
 REAL(DP)             :: COORD_PT (N_SPACE)

! X          = SPATIAL COORDINATES OF ALL SYSTEM NODES
! COORD_PT   = SPATIAL COORDINATES OF THE NODE
! MAX_NP     = TOTAL NUMBER OF NODES IN SYSTEM
! NOD_PER_EL = NUMBER OF NODES PER ELEMENT
! N_SPACE    = DIMENSION OF THE SPACE
   COORD_PT (1:N_SPACE) = X (I_PT, 1:N_SPACE)
END FUNCTION GET_PT_COORD
                     
SUBROUTINE START_ITERATION (X, DD)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!    INITIALIZE SYSTEM DOF FOR ITERATIVE SOLUTION
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Interface_Header
IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: X (MAX_NP, N_SPACE)
 REAL(DP), INTENT(OUT) :: DD (N_D_FRE)
 REAL(DP) :: START_DOF_VALUE
 INTEGER  :: I, INDX, J, IPRINT = 1
!           Automatic Arrays
 REAL(DP) :: XYZ_P (N_SPACE)
 INTEGER  :: INDEX (N_G_DOF)
!
! MAX_NP  = NUMBER OF NODES IN SYSTEM
! N_G_DOF = NUMBER OF PARAMETERS (DOF) PER NODE
! N_SPACE = DIMENSION OF SPACE
! N_D_FRE = TOTAL NUMBER OF SYSTEM DOF
! INDEX   = SYSTEM DOF NUMBERS FOR DOF AT NODE
! X       = COORDINATES OF SYSTEM NODES
! XYZ_P   = SPATIAL COORDINATE ARRAY OF A NODE
! DD      = SYSTEM ARRAY OF DEGREES OF FREEDOM
! IPRINT > 0, PRINT THE STARTING VALUES

  IF ( ECHO_START ) WRITE (N_PRT, 5) N_G_DOF ; 5 FORMAT ( /, &
    '** STARTING VALUES FOR ITERATIVE SOLUTION **',/,        &
    'NODE ', I2,' PARAMETER & VALUE PAIRS')
  DO I = 1, MAX_NP
!    FIND PT COORDS AND DOF NOS
    INDEX = GET_INDEX_AT_PT (I) 
    XYZ_P = GET_PT_COORD (I, X) 
    IF ( ECHO_START ) WRITE (N_PRT, '(I5)', ADVANCE="NO") I 
    DO J = 1, N_G_DOF
      INDX = INDEX (J)
      DD (INDX) = START_DOF_VALUE (J, XYZ_P)
!     START IS A FUNCTION TO DEFINE INITIAL VALUES
!     OF THE SYSTEM DEGREES OF FREEDOM
      IF ( ECHO_START ) WRITE (N_PRT, 510, ADVANCE="NO") J, DD (INDX) 
        510 FORMAT ( I4, 1PE13.5 ) 
    END DO
    IF ( ECHO_START ) WRITE (N_PRT, '(" ")') 
  END DO
END SUBROUTINE START_ITERATION
                              
SUBROUTINE SET_INITIAL_CONDITION (X, DD)
! * * * * * * * * * * * * * * * * * * * * * * * * * * *
! INITIALIZE SYSTEM DOF FOR TIME DEPENDENT INTEGRATION
!        Keyword initial_function should be true
! * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Interface_Header
IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: X (MAX_NP, N_SPACE)
 REAL(DP), INTENT(OUT) :: DD (N_D_FRE)
 REAL(DP) :: START_DOF_VALUE
 INTEGER  :: I, INDX, J, IPRINT = 1
!           Automatic Arrays
 REAL(DP) :: XYZ_P (N_SPACE)
 INTEGER  :: INDEX (N_G_DOF)
!
! MAX_NP  = NUMBER OF NODES IN SYSTEM
! N_G_DOF = NUMBER OF PARAMETERS (DOF) PER NODE
! N_SPACE = DIMENSION OF SPACE
! N_D_FRE = TOTAL NUMBER OF SYSTEM DOF
! INDEX   = SYSTEM DOF NUMBERS FOR DOF AT NODE
! X       = COORDINATES OF SYSTEM NODES
! XYZ_P   = SPATIAL COORDINATE ARRAY OF A NODE
! DD      = SYSTEM ARRAY OF DEGREES OF FREEDOM
! IPRINT > 0, PRINT THE STARTING VALUES

  IF ( IPRINT > 0 ) WRITE (N_PRT, 5) N_G_DOF 
    5 FORMAT ( /, '** SETTING DOF INITIAL CONDITIONS **',/, &
    'NODE ', I2,' PARAMETER & VALUE PAIRS')
  DO I = 1, MAX_NP
!    FIND PT COORDS AND DOF NOS
    INDEX = GET_INDEX_AT_PT (I) 
    XYZ_P = GET_PT_COORD (I, X) 
    IF ( IPRINT > 0 ) WRITE (N_PRT, '(I5)', ADVANCE="NO") I 
    DO J = 1, N_G_DOF
      INDX = INDEX (J)
      DD (INDX) = START_DOF_VALUE   (J, XYZ_P)
   !b DD (INDX) = SET_INITIAL_VALUE (J, XYZ_P) !b
!     START IS A FUNCTION TO DEFINE INITIAL VALUES
!     OF THE SYSTEM DEGREES OF FREEDOM
      IF ( IPRINT > 0 ) WRITE (N_PRT, 510, ADVANCE="NO") J, DD (INDX) 
        510 FORMAT ( I4, 1PE13.5 ) 
    END DO
    IF ( IPRINT > 0 ) WRITE (N_PRT, '(" ")') 
  END DO
END SUBROUTINE SET_INITIAL_CONDITION
                              
SUBROUTINE STORE_COLUMN (LT_FREE, INDEX, C, CC)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! STORE ELEMENT COLUMN MATRIX IN SYSTEM COLUMN MATRIX
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 INTEGER,  INTENT(IN)    :: LT_FREE, INDEX (LT_FREE)
 REAL(DP), INTENT(IN)    :: C (LT_FREE)
 REAL(DP), INTENT(INOUT) :: CC (N_D_FRE)
 INTEGER :: I, J
! INDEX   = SYSTEM DOF NOS OF THE ELEMENT DOF
! C       = ELEMENT COLUMN MATRIX
! CC      = SYSTEM COLUMN MATRIX
! N_D_FRE = NO DEGREES OF FREEDOM IN THE SYSTEM
! LT_FREE = NUMBER OF DEGREES OF FREEDOM PER ELEMENT
  DO I = 1, LT_FREE
    J = INDEX (I)
    IF ( J > 0) CC (J) = CC (J) + C (I)
  END DO
END SUBROUTINE STORE_COLUMN
                           
SUBROUTINE SUM_RHS (CC)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!      SUM INPUT VALUES IN FORCING VECTOR, CC
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 REAL(DP), INTENT(IN) :: CC (N_D_FRE)
!          Automatic Arrays
 REAL(DP) :: TOTAL (N_G_DOF), SUM_POS (N_G_DOF), SUM_NEG (N_G_DOF) !b
 INTEGER  :: I, J, IJ_EQ

! N_D_FRE = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM
! N_G_DOF = NUMBER OF PARAMETERS PER NODE
! CC      = SYSTEM EQUATIONS COLUMN MATRIX
! MAX_NP  = TOTAL NUMBER OF NODES IN THE SYSTEM

  TOTAL = 0.d0 ; SUM_POS = 0.d0 ; SUM_NEG = 0.d0 !b
  WRITE (6, 5030)
  5030 FORMAT (/, '*** INPUT SOURCE RESULTANTS ***',/,               &
                  'ITEM       SUM           POSITIVE      NEGATIVE')
  DO I = 1, MAX_NP
    DO J = 1, N_G_DOF
      IJ_EQ = N_G_DOF * (I - 1) + J
      !b TOTAL (J) = TOTAL (J) + CC (IJ_EQ)
      IF ( CC (IJ_EQ) >= 0.d0 ) THEN   !b
        SUM_POS (J) = SUM_POS (J) + CC (IJ_EQ)   !b
      ELSE   !b
        SUM_NEG (J) = SUM_NEG (J) + CC (IJ_EQ)   !b
      END IF   !b
    END DO
  END DO
  TOTAL = SUM_POS + SUM_NEG
  DO J = 1, N_G_DOF
    !b WRITE (6, '( I5, 2X, 1PE12.4 )') J, TOTAL (J)
    WRITE (6, '( I5, 3(2X, 1PE12.4) )')    &
      J, TOTAL (J), SUM_POS (J), SUM_NEG (J)
  END DO
END SUBROUTINE SUM_RHS
                      
SUBROUTINE CONVECTION_MATRICES (NF, NOD_PER_EL, C, S, & 
           COORD, N_SPACE, IOPT, FLUX, NQ, PT, WT)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!   FACE FLUX MATRICES ON PLANAR FACES (OR ELEMENTS)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 IMPLICIT NONE
 INTEGER, PARAMETER :: N_PARM = 2
!    ALWAYS USED
 INTEGER,  INTENT(IN)  :: NF, NOD_PER_EL, N_SPACE, NQ
 INTEGER,  INTENT(OUT) :: IOPT
 REAL(DP), INTENT(IN)  :: COORD (NOD_PER_EL, N_SPACE), FLUX (NF)
 REAL(DP), INTENT(IN)  :: PT (N_PARM, NQ), WT (NQ)
 REAL(DP), INTENT(OUT) :: C (NOD_PER_EL)
 REAL(DP), INTENT(OUT) :: S (NOD_PER_EL, NOD_PER_EL)
                                                            
 INTEGER  :: I, J, IQ
 REAL(DP) :: F1, F2, F12, DET, DET_WT
!    OPTIONAL FOR NUMERICAL INTEGRATION
 REAL(DP) :: H (NOD_PER_EL ),                       &
             DLH (N_SPACE, NOD_PER_EL ),            &
             AJ (N_SPACE, N_SPACE), AJ_INV (N_SPACE, N_SPACE)
 
! AJ         = JACOBIAN
! AJ_INV     = JACOBIAN INVERSE
! C          = BOUNDARY FLUX COLUMN MATRIX CONTRIBUTIONS
! COORD      = SPATIAL COORDINATES OF SEGMENT NODES
! DLG        = LOCAL DERIVATIVES GEOMETRIC INTERPOLATION
! DLH        = LOCAL DERIVATIVES INTERPOLATION FUNCTIONS
! FLUX       = SPECIFIED BOUNDARY FLUX COMPONENTS
! H          = SOLUTION INTERPOLATION FUNCTIONS
! IOPT       = PROBLEM MATRIX REQUIREMENTS, MUST BE SET.
!            = 0, NO ACTION RETURN
!            = 1, CALCULATE C ONLY
!            = 2, CALCULATE S ONLY
!            = 3, CALCULATE BOTH C AND S
! NOD_PER_EL = NO. OF NODES ON AN ELEMENT BOUNDARY SEGMENT
! NF         = NUMBER OF FLUX COMPONENTS PER NODE
! N_D_FLUX   = N*NF = MAXIMUM NUMBER OF FLUX CONTRIBUTIONS
! N_G_DOF    = NUMBER OF PARAMETERS PER NODE POINT = !
! N_PARM     = PARAMETRIC GEOMETRY SPACE, = N_SPACE USUALLY
! NQ         = NUMBER OF QUADRATURE POINTS
! N_SPACE    = DIMENSION OF SOLUTION SPACE
! PT         = QUADRATURE COORDINATES
! S          = BOUNDARY FLUX SQUARE MATRIX
! WT         = QUADRATURE WEIGHTS
 
  IOPT = 3        ! S & C computed
  F1 = FLUX (1)   ! convection coefficient
  F2 = FLUX (2)   ! reference temperature
  F12 = F1 * F2

!  COMPUTE MATRICES BY NUMERICAL INTEGRATION (ISOPARAMETRIC)
  C = 0.d0 ; S = 0.d0 !  ZERO WORKSPACES
  DO IQ = 1, NQ
    CALL SCALAR_SHAPES (PT (:, IQ),   H)
    CALL SCALAR_DERIVS (PT (:, IQ), DLH)
    AJ = MATMUL ( DLH, COORD ) 
    CALL INVERT_JACOBIAN (AJ, AJ_INV, DET, N_SPACE)
    DET_WT = ABS (DET) * WT (IQ)
    
!      PLANAR FACE (THUS ELEMENT) CONVECTION
    DO I = 1, NOD_PER_EL
      C (I) = C (I) + H (I) * F12 * DET_WT
      DO J = 1, NOD_PER_EL
        S (I, J) = S (I, J) + H (I) * H (J) * DET_WT * F1
      END DO
    END DO
  END DO
END SUBROUTINE CONVECTION_MATRICES
                                 
FUNCTION INVERSION_OF (A_Given)  RESULT (A_inv)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! INVERSION OF NONSYMMETRIC MATRIX A_Given (N_size, N_size)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN) :: A_Given (:, :)
!          Automatic Arrays
 REAL(DP) :: A_inv ( Size (A_Given, 1), Size (A_Given, 1) )
 REAL(DP) :: B ( Size (A_Given, 1) ), C ( Size (A_Given, 1) ), D
 INTEGER  :: I, J, K, M, N_size
! A   IS DESTROYED AND REPLACED BY ITS INVERSE
! B,C = ARE WORKING SPACE VECTORS
! N_size   = SIZE OF GIVEN MATRIX
                                 
  N_size = Size (A_Given, 1)
  A_inv = A_Given         ! Copy into inverse
  IF ( A_inv (1, 1) == 0.0 ) STOP 'ZERO PIVOT IN INVERSION_OF'
  A_inv (1, 1) = 1.d0 / A_inv (1, 1)
                                    
!  NOTE: several loops could be replaced with matmul
  DO M = 1, (N_size - 1)
    K = M + 1
    DO I = 1, M
      B (I) = 0.d0
      DO J = 1, M
        B (I) = B (I) + A_inv (I, J) * A_inv (J, K)
      END DO
    END DO
    D = 0.d0
    DO I = 1, M
      D = D + A_inv (K, I) * B (I)
    END DO
    D = A_inv (K, K) - D
    IF ( D == 0.d0 ) STOP 'ZERO PIVOT IN INVERSION_OF'
    A_inv (K, K) = 1.d0 / D
    DO I = 1, M
      A_inv (I, K) = - B (I) * A_inv (K, K)
    END DO
    DO J = 1, M
      C (J) = 0.d0
      DO I = 1, M
        C (J) = C (J) + A_inv (K, I) * A_inv (I, J)
      END DO
    END DO
    DO J = 1, M
      A_inv (K, J) = - C (J) * A_inv (K, K)
    END DO
    DO I = 1, M
      DO J = 1, M
        A_inv (I, J) = A_inv (I, J) - B (I) * A_inv (K, J)
      END DO
    END DO
  END DO
END FUNCTION INVERSION_OF
                         
!SUBROUTINE REACTEL_ (IE, NOD_PER_EL, N_EL_FRE, N_G_DOF, N_FILE, &
!                     ELEM_NODES, IOPT, D)
!! * * * * * * * * * * * * * * * * * * * * * * * * * * *
!!   GET REACTIONS (FLUXES) AT AN ELEMENTS NODES
!! * * * * * * * * * * * * * * * * * * * * * * * * * * *
!Use Precision_Module
! IMPLICIT NONE
! INTEGER,  INTENT(IN) :: IE, NOD_PER_EL, N_EL_FRE, N_G_DOF
! INTEGER,  INTENT(IN) :: N_FILE, IOPT
! INTEGER,  INTENT(IN) :: ELEM_NODES (NOD_PER_EL)
! REAL(DP), INTENT(IN) :: D (N_EL_FRE)
!                                     
! INTEGER  :: IG, IN, IROW, L
! REAL(DP) :: S (N_EL_FRE, N_EL_FRE), C (N_EL_FRE)
! REAL(DP) :: SUM_PT (N_G_DOF), REACT, ROW
! REAL(DP) :: SUM_EL (N_G_DOF, NOD_PER_EL)
!           
!! C          = ELEMENT COLUMN MATRIX
!! D          = KNOWN SOLUTION CAUSING THE REACTIONS
!! IE         = CURRENT ELEMENT NUMBER
!! IOPT       = WRITE OPTION CODE 1-ONCE, 2-TWICE
!! ELEM_NODES = ELEMENT TOPOLOGY LIST
!! NOD_PER_EL = NUMBER OF NODES PER ELEMENT
!! N_EL_FRE   = NUMBER OF ELEMENT DOF, N_G_DOF*N
!! N_G_DOF    = NUMBER OF GENERALIZED UNKNOWNS PER NODE
!! N_FILE     = UNIT TO HOLDING S & C FROM ELEMENT IE
!! S          = ELEMENT SQUARE MATRIX
!! SUM_EL     = ELEMENT WORKSPACE RESULT
!! SUM_PT     = NODAL WORKSPACE  ESULT
!                                 
!  IF ( N_FILE < 1)  STOP 'INVALID UNIT IN REACTEL'
!!   GIVE REACTIONS, R = S*D - C
!  SUM_PT = 0. ; SUM_EL = 0.
!!  READ STIFFNESS, AND SOURCE FOR ELEM REACTIONS
!!  SINGLE READ
!  IF ( IOPT == 1 ) THEN
!    READ (N_FILE) S, C
!!  DOUBLE READ
!  ELSEIF (IOPT.EQ.2) THEN
!    READ (N_FILE) S ; READ (N_FILE) C
!  END IF
!  DO IN = 1, NOD_PER_EL
!    IF (IN == 1) THEN
!  WRITE (6,  * ) '        ELEMENT', IE, ' REACTIONS'
!  WRITE (6,  * ) ' NODE IDOF    REACTION       SOURCE'
!    END IF
!    DO IG = 1, N_G_DOF
!      IROW           = N_G_DOF * (IN - 1) + IG
!      SUM_EL (IG, 1) = SUM_EL (IG, 1) + C (IROW)
!      ROW            = 0.D0
!      DO L = 1, N_EL_FRE
!        ROW = ROW + S (IROW, L) * D (L) !9 use matmul
!      END DO
!      REACT       = ROW - C (IROW)
!      SUM_PT (IG) = SUM_PT (IG) + REACT
!      PRINT "( I5, I5, 1PE15.5, 1PE15.5 )", &
!               ELEM_NODES (IN), IG, REACT, C (IROW)
!    END DO
!  END DO
!  DO IG = 1, N_G_DOF
!    PRINT "(' SUM:', I5, 1PE15.5, 1PE15.5 )", &
!              IG, SUM_PT (IG), SUM_EL (IG, 1)
!  END DO
!END SUBROUTINE REACTEL_
                       
  SUBROUTINE CAL_PRT (H, PRT_L_PT, VALUES) 
! * * * * * * * * * * * * * * * * * * * * * * * * * *               
!   CALCULATE N_NP_FLO PROPERTIES AT A LOCAL PT USING                 
! ELEMENT'S NODAL PROPERTIES, PRT_L_PT, AND THE NOD_PER_EL     
!      INTERPOLATION FUNCTIONS, H, AT THE POINT                     
! * * * * * * * * * * * * * * * * * * * * * * * * * *               
Use System_Constants
IMPLICIT NONE
  REAL(DP), INTENT(IN)  :: H (NOD_PER_EL)
  REAL(DP), INTENT(IN)  :: PRT_L_PT (NOD_PER_EL, N_NP_FLO)
  REAL(DP), INTENT(OUT) :: VALUES (N_NP_FLO) 
  INTEGER  :: I, J
  REAL(DP) :: SUM   

! NOD_PER_EL = NUMBER OF NODES PER ELEMENT                              
! N_NP_FLO   = NO. OF FLOATING POINT NODAL PROPERTIES                   
! H          = C^0 INTERPOLATION FUNCTIONS FOR AN ELEMENT               
! PRT_L_PT   = FLOATING PT PROPS OF ELEMENT'S NODES                     
! VALUES     = LOCAL VALUES OF PROPERTIES                               
  IF ( N_NP_FLO < 1) STOP 'N_NP_FLO = 0, CAL_PRT' 
  DO I = 1, N_NP_FLO 
    SUM = 0.0 
    DO J = 1, NOD_PER_EL 
      SUM = SUM + H (J) * PRT_L_PT (J, I) 
    END DO 
    VALUES (I) = SUM 
  END DO  !9 use VALUES (I) = SUM ( H(:)*PRT_L_PT (:, I) )
  END SUBROUTINE CAL_PRT                         

!  SUBROUTINE GET_MATERIAL_PROPS (NUM, N_PR_MAT, FLT_MISC, PRT_MAT) 
!    out of date material recovery, see FLT_MAT
!! * * * * * * * * * * * * * * * * * * * * * * * * * *               
!!    EXTRACT REAL MATERIAL PROPERTIES OF MATERIAL                   
!!   NUM FROM MISCELLANEOUS REAL SYSTEM PROPERTIES                   
!! * * * * * * * * * * * * * * * * * * * * * * * * * *               
!Use System_Constants
!IMPLICIT NONE
!  INTEGER,  INTENT(IN)  :: NUM, N_PR_MAT
!  REAL(DP), INTENT(IN)  :: FLT_MISC (MISC_FL)
!  REAL(DP), INTENT(OUT) :: PRT_MAT (N_PR_MAT) 
!  INTEGER :: I, ISTART, N_MAX
!
!! NUM      = MATERIAL NUMBER                                          
!! N_MAX    = MAXIMUM ALLOWABLE MATERIAL NUMBER                        
!! N_PR_MAT = NUMBER OF REAL PROPERTIES PER MATERIAL           
!! MISC_FL  = NUMBER OF MISC REAL SYSTEM PROPERTIES
!! FLT_MISC = SYSTEM STORAGE FOR MISC REAL PROP                        
!! PRT_MAT  = REAL PROPERTY ARRAY FOR MATERIAL NUM, 
!
!!  PROPERTIES STORED IN FLT_MISC IN ORDER OF MATERIAL NUMBER
!  IF ( N_PR_MAT < 1 .OR. N_PR_MAT > MISC_FL ) &
!     STOP 'BAD N_PR_MAT, IN GET_MATERIAL_PROPS' 
!  N_MAX = MISC_FL / N_PR_MAT  ! Max material number
!  IF ( NUM > N_MAX ) STOP 'DIMENSIONS EXCEEDED, GET_MATERIAL_PROPS' 
!  IF ( N_MAX * N_PR_MAT /= MISC_FL ) PRINT *, &
!    'WARNING: MISC_FL > MATERIAL PROPERTY STORAGE'
!  N_WARN = N_WARN + 1 ! INCREMENT WARNING
!  ISTART = N_PR_MAT * (NUM - 1) 
!  DO I = 1, N_PR_MAT 
!    PRT_MAT (I) = FLT_MISC (ISTART + I) 
!  END DO 
!  END SUBROUTINE GET_MATERIAL_PROPS 

!  SUBROUTINE LIST_ALL_MATERIALS (N_PR_MAT, FLT_MISC) 
! out of date method
!  !b SUBROUTINE LIST_ALL_MATERIALS (N_PR_MAT, FLT_MISC, PRT_MAT) 
!! * * * * * * * * * * * * * * * * * * * * * * * * *                 
!!      LIST REAL PROPERTIES BY MATERIAL NUMBER                      
!! * * * * * * * * * * * * * * * * * * * * * * * * *                 
!Use System_Constants
!IMPLICIT NONE
!  INTEGER,  INTENT(IN)  :: N_PR_MAT   
!  REAL(DP), INTENT(IN)  :: FLT_MISC (MISC_FL)
!  !b REAL(DP), INTENT(IN)  :: PRT_MAT (MISC_FL)  
!  REAL(DP) :: PRT_MAT (N_PR_MAT) 
!  INTEGER :: I, J, N_MAX
!  
!! N_PR_MAT = NUMBER OF REAL PROPERTIES PER MATERIAL
!!            <= MISC_FL
!! MISC_FL  = NUMBER OF MISC REAL SYSTEM PROP                          
!! FLT_MISC = SYSTEM STORAGE OF MISC REAL PROP                         
!! PRT_MAT  = REAL PROP ARRAY FOR MATERIAL NUM                         
!! N_MAX    = MAXIMUM ALLOWABLE MATERIAL NUMBER                        
!
!  IF ( N_PR_MAT < 1 .OR. N_PR_MAT > MISC_FL) THEN 
!    WRITE (N_PRT, * ) 'N_PR_MAT, MISC_FL ', N_PR_MAT, MISC_FL 
!    STOP 'BAD N_PR_MAT, IN LIST_ALL_MATERIALS' 
!  ENDIF 
!  N_MAX = MISC_FL / N_PR_MAT 
!  WRITE (N_PRT, 50) N_PR_MAT 
!  50 FORMAT (/, '** LIST OF REAL MATERIAL PROPERTIES **',/, &
!  'MATERIAL,',I4,' REAL PROPERTIES')                                
!  DO I = 1, N_MAX 
!    CALL GET_MATERIAL_PROPS (I, N_PR_MAT, FLT_MISC, PRT_MAT) 
!    WRITE (N_PRT, 5010) I, (PRT_MAT (J), J = 1, N_PR_MAT) 
!    5010 FORMAT ( I5, (10(1PE12.4)) ) 
!  END DO 
!  END SUBROUTINE LIST_ALL_MATERIALS                         

SUBROUTINE PT_CODE (JPT, KODE, KODES)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! EXTRACT B.C. INDICATORS AT NODE NUMBER JPT
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 INTEGER, INTENT(IN)  :: JPT, KODE
 INTEGER, INTENT(OUT) :: KODES (N_G_DOF)
 INTEGER :: I, II, IK, INEW, IOLD, ISUM

! JPT     = NODE NO.
! N_G_DOF = NO. PARAMETERS PER NODE
! KODE    = (N_G_DOF) DIGIT INTEGER CONTAINING BC INDICATORS
! KODES   = VECTOR CONTAINING N_G_DOF INTEGER CODES (0 OR I)
!             0  IMPLIES NO B. C.
!             J  IMPLIES A B. C.  OF TYPE J
  IOLD = KODE
  ISUM = 0
  DO I = 1, N_G_DOF
    II   = N_G_DOF + 1 - I
    INEW = IOLD / 10
    IK   = IOLD - INEW * 10
    ISUM = ISUM + IK * 10** (I - 1)
    IOLD = INEW
    KODES (II) = IK
    IF ( IK > MAX_TYP ) THEN
      PRINT *, 'INVALID BC CODE, NODE = ', JPT, ', CODE = ', KODE
      PRINT *, 'INVALID BC TYPE, NODE = ', JPT, ', PARM = ', I
      PRINT *, 'Set keyword MAX_TYP to ', IK 
      STOP     'INVALID BC TYPE, PT_CODE'
    END IF
  END DO
!  WAS DATA RIGHT JUSTIFIED?
  IF ( KODE > ISUM ) THEN
    WRITE (N_BUG, * ) 'WARNING: BC NOT RIGHT JUSTIFIED AT NODE ', JPT
    N_WARN = N_WARN + 1 ! INCREMENT WARNING
  END IF
END SUBROUTINE PT_CODE
                    
SUBROUTINE COUNT_L_AT_NODE (N_ELEMS, NOD_PER_EL, MAX_NP, &
                             NODES, L_TO_N_SUM)
! * * * * * * * * * * * * * * * * * * * * * * * * * *              
! COUNT NUMBER OF ELEMENTS ADJACENT TO EACH NODE
!         (TO SIZE ELEM ADJACENT TO ELEM LIST)
! * * * * * * * * * * * * * * * * * * * * * * * * * *             
 IMPLICIT NONE
 INTEGER, INTENT(IN)  :: N_ELEMS, NOD_PER_EL, MAX_NP
 INTEGER, INTENT(IN)  :: NODES (N_ELEMS, NOD_PER_EL)
 INTEGER, INTENT(OUT) :: L_TO_N_SUM (MAX_NP)

 INTEGER :: ELEM_NODES (NOD_PER_EL)
 INTEGER :: IE, IN, N_TEMP

! ELEM_NODES     = INCIDENCES ARRAY FOR A SINGLE ELEMENT     
! L_TO_N_SUM (I) = NUMBER OF ELEM NEIGHBORS OF NODE I
! MAX_NP         = TOTAL NUMBER OF NODES                          
! N_ELEMS        = TOTAL NUMBER OF ELEMENTS
! NOD_PER_EL     = NUMBER OF NODES PER ELEMENT               
! NODES          = SYSTEM ARRAY OF INCIDENCES OF ALL ELEMENTS      

  L_TO_N_SUM = 0               ! INITIALIZE
  DO IE = 1, N_ELEMS           ! LOOP OVER ALL ELEMENTS 
!        EXTRACT INCIDENCES LIST FOR ELEMENT IE   
    CALL ELEMENT_NODES (IE, NOD_PER_EL, NODES, ELEM_NODES) 
    DO IN = 1, NOD_PER_EL      ! loop over each node
      N_TEMP = ELEM_NODES (IN)
      IF ( N_TEMP < 1 )  CYCLE ! to a real node
        L_TO_N_SUM (N_TEMP) = L_TO_N_SUM (N_TEMP) + 1
    END DO ! over IN
  END DO ! over elements
END SUBROUTINE COUNT_L_AT_NODE 

SUBROUTINE FORM_L_ADJACENT_NODES (N_ELEMS, NOD_PER_EL, MAX_NP, &
                                NODES, NEIGH_N, L_TO_N_NEIGH)
! * * * * * * * * * * * * * * * * * * * * * * * * * *              
!       TABULATE ELEMENTS ADJACENT TO EACH NODE
!         (TO SIZE ELEM ADJACENT TO ELEM LIST)
! * * * * * * * * * * * * * * * * * * * * * * * * * *             
 IMPLICIT NONE
 INTEGER, INTENT(IN)  :: N_ELEMS, NOD_PER_EL, MAX_NP, NEIGH_N
 INTEGER, INTENT(IN)  :: NODES (N_ELEMS, NOD_PER_EL)
 INTEGER, INTENT(OUT) :: L_TO_N_NEIGH (NEIGH_N, MAX_NP)

 INTEGER :: ELEM_NODES (NOD_PER_EL)     ! scratch area
 INTEGER :: COUNT (MAX_NP)              ! scratch area
 INTEGER :: IE, IN, N_TEMP

! ELEM_NODES     = INCIDENCES ARRAY FOR A SINGLE ELEMENT     
! L_TO_N_SUM (I) = NUMBER OF ELEM NEIGHBORS OF NODE I
! MAX_NP         = TOTAL NUMBER OF NODES                          
! N_ELEMS        = TOTAL NUMBER OF ELEMENTS
! NOD_PER_EL     = NUMBER OF NODES PER ELEMENT               
! NODES          = SYSTEM ARRAY OF INCIDENCES OF ALL ELEMENTS      

  L_TO_N_NEIGH = 0             ! INITIALIZE
  COUNT = 0                    ! INITIALIZE
  DO IE = 1, N_ELEMS           ! LOOP OVER ALL ELEMENTS 

!        EXTRACT INCIDENCES LIST FOR ELEMENT IE   
    CALL ELEMENT_NODES (IE, NOD_PER_EL, NODES, ELEM_NODES) 

    DO IN = 1, NOD_PER_EL      ! loop over each node
      N_TEMP = ELEM_NODES (IN)
      IF ( N_TEMP < 1 )  CYCLE ! to a real node
        COUNT (N_TEMP) = COUNT (N_TEMP) + 1
        IF ( COUNT (N_TEMP) > NEIGH_N ) &
          STOP 'INCONSISTENT NEIGH_N; FORM_L_ADJACENT_NODES'
          L_TO_N_NEIGH (COUNT (N_TEMP), N_TEMP) = IE
    END DO ! over IN
  END DO ! over elements
END SUBROUTINE FORM_L_ADJACENT_NODES

!SUBROUTINE COUNT_L_AT_ELEM (N_ELEMS, NOD_PER_EL, MAX_NP, &
!                            NODES, L_TO_N_SUM, L_TO_L_SUM)
!! * * * * * * * * * * * * * * * * * * * * * * * * * *              
!! -------This is incorrect except for quadratic elements-------
!! COUNT NUMBER OF ELEMENTS ADJACENT TO OTHER ELEMENTS              
!!     (SEE ALTERNATE: COUNT_ELEMS_AT_ELEM)
!! * * * * * * * * * * * * * * * * * * * * * * * * * *             
! IMPLICIT NONE
! INTEGER, INTENT(IN)  :: N_ELEMS, NOD_PER_EL, MAX_NP
! INTEGER, INTENT(IN)  :: NODES (N_ELEMS, NOD_PER_EL)
! INTEGER, INTENT(OUT) :: L_TO_N_SUM (MAX_NP)
! INTEGER, INTENT(OUT) :: L_TO_L_SUM (N_ELEMS)
!
! INTEGER :: ELEM_NODES (NOD_PER_EL)
! INTEGER :: FOUND, IE, IN, N_TEST
!
!! ELEM_NODES     = INCIDENCES ARRAY FOR A SINGLE ELEMENT     
!! L_TO_L_SUM (I) = NUMBER OF ELEM NEIGHBORS OF ELEMENT I
!! L_TO_N_SUM (I) = NUMBER OF ELEM NEIGHBORS OF NODE I
!! MAX_NP         = TOTAL NUMBER OF NODES                          
!! N_ELEMS        = TOTAL NUMBER OF ELEMENTS
!! NOD_PER_EL     = NUMBER OF NODES PER ELEMENT               
!! NODES          = SYSTEM ARRAY OF INCIDENCES OF ALL ELEMENTS      
!
!  DO IE = 1, N_ELEMS      ! LOOP OVER ALL ELEMENTS 
!if (ie ==1 ) stop'COUNT_L_AT_ELEM is incorrect'
!    FOUND = 0
!
!!        EXTRACT INCIDENCES LIST FOR ELEMENT IE   
!    CALL ELEMENT_NODES (IE, NOD_PER_EL, NODES, ELEM_NODES) 
!
!    DO IN = 1, NOD_PER_EL  ! LOOP OVER NODES OF ELEMENT
!      N_TEST = ELEM_NODES (IN)
!      IF ( N_TEST > 0 ) FOUND = FOUND + L_TO_N_SUM (N_TEST)
!    END DO ! over in
!
!!       REMOVE ELEMENTS COUNTED TWICE
!    FOUND = FOUND - 2 * COUNT ( ELEM_NODES > 0 )
!    L_TO_L_SUM (IE) = FOUND
!  END DO ! over all elements
!END SUBROUTINE COUNT_L_AT_ELEM 

SUBROUTINE COUNT_ELEMS_AT_ELEM (N_ELEMS, NOD_PER_EL, MAX_NP,   &
                                L_FIRST, L_LAST, NODES, NEEDS, &
                                L_TO_L_SUM, N_WARN)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *    
! COUNT NUMBER OF ELEMENTS ADJACENT TO OTHER ELEMENTS              
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *     
 IMPLICIT NONE
 INTEGER, INTENT(IN)    :: N_ELEMS, NOD_PER_EL, MAX_NP, NEEDS
 INTEGER, INTENT(IN)    :: L_FIRST (MAX_NP), L_LAST (MAX_NP)
 INTEGER, INTENT(IN)    :: NODES (N_ELEMS, NOD_PER_EL)
 INTEGER, INTENT(OUT)   :: L_TO_L_SUM (N_ELEMS)
 INTEGER, INTENT(INOUT) :: N_WARN

 INTEGER :: ELEM_NODES (NOD_PER_EL)
 INTEGER :: NEIG_NODES (NOD_PER_EL)
 INTEGER :: FOUND, IE, IN, L_TEST, L_START, L_STOP, N_TEST, NULLS
 INTEGER :: NEED, KOUNT

! ELEM_NODES     = INCIDENCES ARRAY FOR A SINGLE ELEMENT     
! KOUNT          = NUMBER OF COMMON NODES
! L_FIRST (I)    = ELEMENT WHERE NODE I FIRST APPEARS              
! L_LAST (I)     = ELEMENT WHERE NODE I LAST APPEARS             
! L_TO_L_SUM (I) = NUMBER OF ELEM NEIGHBORS OF ELEMENT I
! MAX_NP         = TOTAL NUMBER OF NODES                          
! NEEDS          = NUMBER OF COMMON NODES TO BE A NEIGHBOR
! N_ELEMS        = TOTAL NUMBER OF ELEMENTS
! NOD_PER_EL     = NUMBER OF NODES PER ELEMENT               
! NODES          = SYSTEM ARRAY OF INCIDENCES OF ALL ELEMENTS      

  NEED       = MAX (1, NEEDS)   ! INITIALIZE
  L_TO_L_SUM = 0

  MAIN : DO IE = 1, N_ELEMS      ! LOOP OVER ALL ELEMENTS 
    FOUND = 0 

!        EXTRACT INCIDENCES LIST FOR ELEMENT IE   
    CALL ELEMENT_NODES (IE, NOD_PER_EL, NODES, ELEM_NODES) 

!        ESTABLISH RANGE OF POSSIBLE ELEMENT NEIGHBORS    
    L_START = N_ELEMS ; L_STOP = 0 
    DO IN = 1, NOD_PER_EL 
      N_TEST  = ELEM_NODES (IN) !b
      IF ( N_TEST < 1 ) CYCLE! to a real node !b
        L_START = MIN (L_START, L_FIRST (N_TEST) )          !b
        L_STOP  = MAX (L_STOP,  L_LAST  (N_TEST) )          !b
    END DO 

!        LOOP OVER POSSIBLE ELEMENT NEIGHBORS            
    IF ( L_START <= L_STOP) THEN 
      RANGE : DO L_TEST = L_START, L_STOP 
        IF ( L_TEST /= IE) THEN 
          KOUNT = 0  ! NO COMMON NODES

!      LOOP OVER INCIDENCES OF POSSIBLE ELEMENT NEIGHBOR        
          CALL ELEMENT_NODES (L_TEST, NOD_PER_EL, NODES, NEIG_NODES) 
          LOCAL : DO IN = 1, NOD_PER_EL 
            N_TEST = NEIG_NODES (IN) 
            IF ( N_TEST < 1 .OR. N_TEST > MAX_NP )  THEN
              PRINT *, 'COUNT_ELEMS_AT_ELEM: WARNING, ', &
                       'SKIPPING INVALID NODE ', N_TEST, &
                       ' OF ELEMENT ', L_TEST
              PRINT *, NEIG_NODES
              N_WARN = N_WARN + 1 ! INCREMENT WARNING
              CYCLE LOCAL ! to a real node
            END IF ! IMPOSSIBLE NODE
            IF ( L_FIRST (N_TEST)  > IE ) CYCLE LOCAL ! to next node
            IF ( L_LAST  (N_TEST)  < IE ) CYCLE LOCAL ! to next node

!              COMPARE WITH INCIDENCES OF ELEMENT IE             
              IF ( ANY ( ELEM_NODES == N_TEST ) ) THEN
                KOUNT = KOUNT + 1
                IF ( KOUNT == NEED ) THEN ! IS A NEIGHBOR 
                  FOUND = FOUND + 1
                  EXIT LOCAL ! this L_TEST element search loop
                END IF ! NUMBER NEEDED
              END IF
          END DO LOCAL ! over in
        END IF 
      END DO RANGE ! over candidate element L_TEST
    END IF ! a possible candidate
    L_TO_L_SUM (IE) = FOUND
  END DO MAIN ! over all elements

  PRINT *, 'NOTE: MAXIMUM NUMBER OF ELEMENT NEIGHBORS = ', &
            MAXVAL ( L_TO_L_SUM )
  NULLS = COUNT ( L_TO_L_SUM == 0 ) !  CHECK DATA
  IF ( NULLS > 0 ) THEN
    PRINT *, 'WARNING, ', NULLS, ' ELEMENTS HAVE NO ELEMENT NEIGHBORS'
    N_WARN = N_WARN + 1 ! INCREMENT WARNING
  END IF
END SUBROUTINE COUNT_ELEMS_AT_ELEM                         

SUBROUTINE FIRST_LAST_L_AT_PT (NODES, L_FIRST, L_LAST, I_WARN)
!* * * * * * * * * * * * * * * * * * * * * * * * * * * *
! FIND ELEMENT OF FIRST & LAST APPEARANCE OF EACH NODE
!* * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_type_data
Use Interface_header
 IMPLICIT NONE
 INTEGER, INTENT(IN)  :: I_WARN
 INTEGER, INTENT(IN)  :: NODES (L_S_TOT, NOD_PER_EL)
 INTEGER, INTENT(OUT) :: L_FIRST (MAX_NP), L_LAST (MAX_NP)

 INTEGER, ALLOCATABLE :: EL_NODES (:) 
 INTEGER :: I, IE, L, LT

! EL_NODES = INCIDENCES ARRAY OF A SINGLE ELEMENT   
! I_WARN     > 0  WARN OF NODES WITH NO ELEMENTS
! L_FIRST(I) = ELEMENT OF FIRST APPEARANCE OF NODE I
! L_LAST(I)  = ELEMENT OF LAST APPEARANCE OF NODE I
! LT_N       = MAXIMUM NUMBER OF NODES FOR ELEMENT TYPE
! L_S_TOT    = TOTAL GENERALIZED ELEMENTS, N_ELEMS + N_MIXED + N_F_SEG
! MAX_NP     = TOTAL NUMBER OF NODES IN SYSTEM
! NODES      = SYSTEN ARRAY OF ALL ELEMENT INCIDENCES
! NOD_PER_EL = NUMBER OF NODES PER ELEMENT
! N_ELEMS    = NUMBER OF ELEMENTS IN THE SYSTEM
! N_L_TYPE   = NUMBER OF DIFFERENT ELEMENT TYPES USED

  L_FIRST = 0 ; L_LAST = 0 ! initialize
  LT = 1 ! INITIALIZE ELEMENT TYPE
  DO IE = 1, N_ELEMS
!-->     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
!      GET CONTROLS FOR THIS TYPE
      CALL GET_ELEM_TYPE_DATA (LT) !b SET_ELEM_TYPE_INFO (LT)
      LAST_LT = LT

!-->        Allocate Arrays for this type
      IF ( IE > 1 ) DEALLOCATE ( EL_NODES )    !bu
      ALLOCATE ( EL_NODES (LT_N) )   !bu
    END IF ! a new element type

!    EXTRACT ELEMENT'S NODES
    EL_NODES = GET_ELEM_NODES (IE, LT_N, NODES)   

!    SCAN THE NODES
    DO I = 1, LT_N      
      L = EL_NODES (I)   
      IF ( L < 1 ) CYCLE ! to next node
      IF ( L_FIRST (L) == 0 ) L_FIRST (L) = IE
      L_LAST (L) = IE
    END DO
  END DO
  DEALLOCATE ( EL_NODES )   

!  WARN OF NODES WITH NO ELEMENT CONNECTIONS
  IF ( I_WARN == 0 .OR. COUNT( L_FIRST == 0 ) == 0 ) RETURN
  DO I = 1, MAX_NP
    IF ( L_FIRST (I) > 0 ) CYCLE ! to next node
      PRINT *, ' '
      WRITE (6, 5) I ; 5 FORMAT ('WARNING, NODE ', I6, &
      & ' DOES NOT OCCUR IN THE ELEMENT CONNECTION LIST')
      N_WARN = N_WARN + 1 ! INCREMENT WARNING
  END DO
END SUBROUTINE FIRST_LAST_L_AT_PT

SUBROUTINE LIST_FIRST_LAST (L_FIRST, L_LAST)
!* * * * * * * * * * * * * * * * * * * * * * * * * * * *
! LIST ELEMENT OF FIRST & LAST APPEARANCE OF EACH NODE
!* * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 INTEGER, INTENT(IN) :: L_FIRST (MAX_NP), L_LAST (MAX_NP)
 INTEGER :: I

! L_FIRST(I) = ELEMENT OF FIRST APPEARANCE OF NODE I
! L_LAST(I)  = ELEMENT OF LAST APPEARANCE OF NODE I
! MAX_NP     = TOTAL NUMBER OF NODES IN SYSTEM

  WRITE (N_PRT, *) ' '
  WRITE (N_PRT, *) '*** ELEMENTS AT NODE LIST ***'
  WRITE (N_PRT, *) 'NODE  FIRST_ELEM  LAST_ELEM'  
  DO I = 1, MAX_NP
    WRITE (N_PRT, '(I6, 5X, I6, 5X, I6)') &
           I, L_FIRST (I), L_LAST (I)
  END DO
  WRITE (N_PRT, *) ' '
END SUBROUTINE LIST_FIRST_LAST

SUBROUTINE LIST_L_AT_NODE_INFO (L_TO_N_SUM, L_FIRST, L_LAST)
!* * * * * * * * * * * * * * * * * * * * * * * * * * * *
! LIST ELEMENT OF FIRST & LAST APPEARANCE OF EACH NODE
!        AND THE NUMBER OF ELEMENTS ATTACHED
!* * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 INTEGER, INTENT(IN) :: L_FIRST (MAX_NP), L_LAST (MAX_NP)
 INTEGER, INTENT(IN) :: L_TO_N_SUM (MAX_NP)
 INTEGER :: I

! L_FIRST(I)     = ELEMENT OF FIRST APPEARANCE OF NODE I
! L_LAST(I)      = ELEMENT OF LAST APPEARANCE OF NODE I
! L_TO_N_SUM (I) = NUMBER OF ELEM NEIGHBORS OF NODE I
! MAX_NP         = TOTAL NUMBER OF NODES IN SYSTEM

  WRITE (N_PRT, *) ' '
  WRITE (N_PRT, *) '*** ELEMENTS AT NODE SUMMARY ***'
  WRITE (N_PRT, *) '  NODE  TOTAL_ELEM FIRST_ELEM  LAST_ELEM'  
  DO I = 1, MAX_NP
    WRITE (N_PRT, '(I6, 1X, I6, 5X, I6, 5X, I6)') &
           I, L_TO_N_SUM (I), L_FIRST (I), L_LAST (I)
  END DO
  WRITE (N_PRT, *) ' '
END SUBROUTINE LIST_L_AT_NODE_INFO

SUBROUTINE LIST_ELEMENTS_AT_EL (L_TO_L_NEIGH)
!* * * * * * * * * * * * * * * * * * * * * * * * * * * *
! LIST ELEMENT NEIGHBORS OF EACH ELEMENT
!* * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants  ! for NEIGH_L, N_ELEMS
 IMPLICIT NONE
 INTEGER, INTENT(IN) :: L_TO_L_NEIGH (NEIGH_L, N_ELEMS)
 INTEGER :: IE

! L_TO_L_NEIGH (J, I) = ELEMENT NEIGHBOR J OF ELEMENT I
! NEIGH_L             = MAX NUMBER OF ELEM NEIGHBORS AT ANY ELEM
! N_ELEMS             = MAXIMUM NUMBER OF ELEMENTS IN SYSTEM

  WRITE (N_PRT, *) ' '
  WRITE (N_PRT, *) '*** NEIGHBORING ELEMENTS AT EACH ELEMENT ***'
  WRITE (N_PRT, *) 'ELEMENT,    ADJACENT ELEMENT LIST'      
  DO IE = 1, N_ELEMS
    WRITE (N_PRT, 5) IE, L_TO_L_NEIGH (1:NEIGH_L, IE)
    5 FORMAT (I6, ',', 10I6, (/, 7X, 10I6) )
  END DO
  WRITE (N_PRT, *) ' '
END SUBROUTINE LIST_ELEMENTS_AT_EL

SUBROUTINE LIST_ELEMENTS_AT_PT (L_TO_N_NEIGH)
!* * * * * * * * * * * * * * * * * * * * * * * * * * * *
!        LIST ELEMENT NEIGHBORS OF EACH NODE
!* * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants   ! for NEIGH_N, MAX_NP
 IMPLICIT NONE
 INTEGER, INTENT(IN) :: L_TO_N_NEIGH (NEIGH_N, MAX_NP)
 INTEGER :: IP

! L_TO_N_NEIGH (J, I) = ELEMENT NEIGHBOR J OF NODE PT I
! NEIGH_N             = MAXIMUM NUMBER OF ELEM NEIGHBORS AT ANY ELEM
! MAX_NP              = MAXIMUM NUMBER OF NODE PTS IN SYSTEM

  WRITE (N_PRT, *) ' '
  WRITE (N_PRT, *) '*** NEIGHBORING ELEMENTS AT EACH NODE ***'
  WRITE (N_PRT, *) ' NODE,    ADJACENT ELEMENT LIST'      
  DO IP = 1, MAX_NP 
    WRITE (N_PRT, 5) IP, L_TO_N_NEIGH (1:NEIGH_N, IP)
    5 FORMAT (I8, ',   ', 10I6, (/, 12X, 10I6) )
  END DO
  WRITE (N_PRT, *) ' '
END SUBROUTINE LIST_ELEMENTS_AT_PT

SUBROUTINE FORM_ELEMS_AT_EL (N_ELEMS, NOD_PER_EL, MAX_NP, &
                            L_FIRST, L_LAST, NODES, N_SPACE, &
                            L_TO_L_SUM, L_TO_L_NEIGH,     &
                            NEIGH_L, NEEDS, ON_BOUNDARY, U_ON_B)
! * * * * * * * * * * * * * * * * * * * * * * * * * *         
! FORM LIST OF ELEMENTS ADJACENT TO OTHER ELEMENTS              
! * * * * * * * * * * * * * * * * * * * * * * * * * *        
 IMPLICIT NONE
 INTEGER, INTENT(IN)    :: N_ELEMS, NOD_PER_EL, MAX_NP, NEIGH_L
 INTEGER, INTENT(IN)    :: L_FIRST (MAX_NP), L_LAST (MAX_NP)
 INTEGER, INTENT(IN)    :: NODES (N_ELEMS, NOD_PER_EL)
 INTEGER, INTENT(IN)    :: L_TO_L_SUM   (N_ELEMS)
 INTEGER, INTENT(IN)    :: N_SPACE, NEEDS, U_ON_B
 INTEGER, INTENT(OUT)   :: L_TO_L_NEIGH (NEIGH_L, N_ELEMS)
 LOGICAL, INTENT(INOUT) :: ON_BOUNDARY  (N_ELEMS)

 INTEGER :: ELEM_NODES (NOD_PER_EL)
 INTEGER :: NEIG_NODES (NOD_PER_EL)
 INTEGER :: IE, IN, L_TEST, L_START, L_STOP, N_TEST
 INTEGER :: FOUND, NEXT ! add NEED for face or edge sub-sets
 INTEGER :: SUM_L_TO_L
 INTEGER :: IO_1, KOUNT, NEED, N_FACES, WHERE !b

! ON_BOUNDARY         = TRUE IF ELEMENT HAS A FACE ON THE BOUNDARY
! ELEM_NODES          = INCIDENCES ARRAY FOR A SINGLE ELEMENT     
! FOUND               = CURRENT NUMBER OF LOCAL NEIGHBORS
! KOUNT               = CURRENT NUMBER OF COMMON NODES
! L_FIRST (I)         = ELEMENT WHERE NODE I FIRST APPEARS      
! L_LAST (I)          = ELEMENT WHERE NODE I LAST APPEARS
! L_TO_L_NEIGH (J, I) = ELEM NEIGHBOR J OF ELEMENT I
! L_TO_L_SUM (I)      = NUMBER OF ELEM NEIGHBORS OF ELEMENT I
! NEEDS               = NUMBER OF COMMON NODES TO BE A NEIGHBOR
! NEIGH_L             = MAXIMUM NUMBER OF NEIGHBORS AT AN ELEMENT 
! MAX_NP              = TOTAL NUMBER OF NODES                  
! N_ELEMS             = TOTAL NUMBER OF ELEMENTS
! NOD_PER_EL          = NUMBER OF NODES PER ELEMENT               
! NODES               = SYSTEM ARRAY OF INCIDENCES OF ALL ELEMENTS
! WHERE               = LOCATION TO INSERT NEIGHBOR, <= MAX_FACES

  NEED         = MAX (1, NEEDS)   ! INITIALIZE
  L_TO_L_NEIGH = 0                ! INITIALIZE
  MAIN : DO IE = 1, N_ELEMS      ! LOOP OVER ALL ELEMENTS 

    SUM_L_TO_L = L_TO_L_SUM (IE)               ! MAX POSSIBLE NEIGHBORS
    FOUND = COUNT ( L_TO_L_NEIGH (:, IE) > 0 ) ! PREVIOUSLY FOUND
    IF ( FOUND == SUM_L_TO_L ) CYCLE MAIN      ! ALL HAVE BEEN FOUND

!        EXTRACT INCIDENCES LIST FOR ELEMENT IE   
    CALL ELEMENT_NODES (IE, NOD_PER_EL, NODES, ELEM_NODES) 

!        ESTABLISH RANGE OF POSSIBLE ELEMENT NEIGHBORS
    L_START = N_ELEMS + 1 ; L_STOP = 0 
    DO IN = 1, NOD_PER_EL 
      L_TEST  = ELEM_NODES (IN)    !b
      IF ( L_TEST < 1 ) CYCLE ! TO NEXT NODE    !b
      L_START = MIN (L_START, L_FIRST (L_TEST))              !b
      L_STOP  = MAX (L_STOP,  L_LAST  (L_TEST))              !b
      !b L_START = MIN (L_START, L_FIRST (ELEM_NODES (IN)) )
      !b L_STOP  = MAX (L_STOP,  L_LAST  (ELEM_NODES (IN)) )
    END DO 
    L_START = MAX (L_START, IE+1) ! SEARCH ABOVE IE ONLY

!        LOOP OVER POSSIBLE ELEMENT NEIGHBORS
    IF ( L_START <= L_STOP) THEN 
      RANGE : DO L_TEST = L_START, L_STOP 
        KOUNT = 0  ! NO COMMON NODES

!         EXTRACT NODES OF L_TEST
        CALL ELEMENT_NODES (L_TEST, NOD_PER_EL, NODES, NEIG_NODES) 

!         LOOP OVER INCIDENCES OF POSSIBLE ELEMENT NEIGHBOR        
        LOCAL : DO IN = 1, NOD_PER_EL 
          N_TEST = NEIG_NODES (IN) 
          IF ( N_TEST <  1 ) CYCLE ! to a real node
          IF ( L_FIRST (N_TEST)  > IE ) CYCLE LOCAL ! to next node
          IF ( L_LAST  (N_TEST)  < IE ) CYCLE LOCAL ! to next node

!            COMPARE WITH INCIDENCES OF ELEMENT IE             
            IF ( ANY ( ELEM_NODES == N_TEST ) ) THEN
              KOUNT = KOUNT + 1         ! SHARED NODE COUNT
              IF ( KOUNT == NEED ) THEN ! NEIGHBOR PAIR FOUND
                FOUND = FOUND + 1       ! INSERT THE PAIR

!               NOTE: THIS INSERT IS NOT ORDERED. IF A SPECIFIC FACE
!                     ORDER IS REQUIRED DETERMINE WHERE HERE

                WHERE = FOUND ! OR ORDER THE CURRENT FACE
                L_TO_L_NEIGH (WHERE, IE) = L_TEST     ! ONE OF TWO

                NEXT = COUNT ( L_TO_L_NEIGH(:, L_TEST) > 0 ) !b /= ??
                WHERE = NEXT+1 ! OR ORDER THE NEIGHBOR FACE

                IF ( L_TO_L_SUM (L_TEST) > NEXT )  &
                  L_TO_L_NEIGH (NEXT+1, L_TEST) = IE  ! TWO OF TWO
                IF ( SUM_L_TO_L == FOUND ) CYCLE MAIN ! ALL FOUND
                CYCLE RANGE ! this L_TEST element search loop
              END IF ! NUMBER NEEDED
            END IF ! SHARE AT LEAST ONE COMMON NODE
        END DO LOCAL ! over N_TEST
      END DO RANGE ! over candidate element L_TEST
    END IF ! a possible candidate
  END DO MAIN! over all elements

  IF ( NEED >= N_SPACE ) THEN ! DATA ARE EDGE OR FACE NEIGHBORS
!b call at(2698)
    ON_BOUNDARY = .FALSE.
!      SAVE THE ELEMENT NUMBERS THAT FACE THE BOUNDARY
    DO IE = 1, N_ELEMS
      CALL GET_LT_FACES (IE, N_FACES) 
      IF ( N_FACES > 0 ) THEN ! MIGHT BE ON THE BOUNDARY
        IF ( NEIGH_L < N_FACES ) THEN
          PRINT *,'NOTE: REDUCING FACES FROM ', N_FACES, &
                & ' TO ', NEIGH_L
          N_FACES = NEIGH_L ! Possibly valid mesh
        END IF ! validation
        IF ( ANY (L_TO_L_NEIGH (1:N_FACES, IE) == 0) ) THEN ! YES
          ON_BOUNDARY (IE) = .TRUE.
        END IF  ! ON BOUNDARY
      END IF ! POSSIBLE ELEMENT
    END DO ! OVER ELEMENTS
  END IF ! SEARCH OF FACING NEIGHBORS

  IF ( ANY ( ON_BOUNDARY ) ) THEN ! SAVE ELEMENTS ON BOUNDARY
    OPEN (U_ON_B, FILE='boundary_face_elems.tmp',        &
          ACTION='WRITE', STATUS='REPLACE', IOSTAT = IO_1) 
    IF ( IO_1 > 0 ) PRINT *, 'OPEN FAILED FOR boundary_face_elems.tmp'
    DO IE = 1, N_ELEMS
      IF ( ON_BOUNDARY (IE) ) WRITE (U_ON_B, *) IE
    END DO
    PRINT *,'NOTE: Saved ', COUNT ( ON_BOUNDARY ), &
            ' elements on the boundary.'
  END IF ! SAVE ELEMENTS ON BOUNDARY
END SUBROUTINE FORM_ELEMS_AT_EL                         

!!SUBROUTINE FORM_ELEMS_AT_EL_OLD (N_ELEMS, NOD_PER_EL, MAX_NP,    &
!                                L_FIRST, L_LAST, NODES,         &
!                                L_TO_L_SUM, L_TO_L_NEIGH, NEIGH_L)
!! * * * * * * * * * * * * * * * * * * * * * * * * * *         
!! FORM LIST OF ELEMENTS ADJACENT TO OTHER ELEMENTS              
!! * * * * * * * * * * * * * * * * * * * * * * * * * *        
! IMPLICIT NONE
! INTEGER, INTENT(IN)  :: N_ELEMS, NOD_PER_EL, MAX_NP, NEIGH_L
! INTEGER, INTENT(IN)  :: L_FIRST (MAX_NP), L_LAST (MAX_NP)
! INTEGER, INTENT(IN)  :: NODES (N_ELEMS, NOD_PER_EL)
! INTEGER, INTENT(IN)  :: L_TO_L_SUM (N_ELEMS)
! INTEGER, INTENT(OUT) :: L_TO_L_NEIGH (NEIGH_L, N_ELEMS)
!
! INTEGER :: ELEM_NODES (NOD_PER_EL)
! INTEGER :: NEIG_NODES (NOD_PER_EL)
! INTEGER :: IE, IN, L_TEST, L_START, L_STOP, N_TEST
! INTEGER :: FOUND
!
!! ELEM_NODES          = INCIDENCES ARRAY FOR A SINGLE ELEMENT     
!! FOUND               = CURRENT NUMBER OF LOCAL NEIGHBORS
!! L_FIRST (I)         = ELEMENT WHERE NODE I FIRST APPEARS      
!! L_LAST (I)          = ELEMENT WHERE NODE I LAST APPEARS
!! L_TO_L_NEIGH (J, I) = ELEM NEIGHBOR J OF ELEMENT I
!! L_TO_L_SUM (I)      = NUMBER OF ELEM NEIGHBORS OF ELEMENT I
!! NEIGH_L           = MAXIMUM NUMBER OF NEIGHBORS AT AN ELEMENT 
!! MAX_NP              = TOTAL NUMBER OF NODES                  
!! N_ELEMS             = TOTAL NUMBER OF ELEMENTS
!! NOD_PER_EL          = NUMBER OF NODES PER ELEMENT               
!! NODES               = SYSTEM ARRAY OF INCIDENCES OF ALL ELEMENTS
!
!  L_TO_L_NEIGH = 0        ! INITIALIZE
!  DO IE = 1, N_ELEMS      ! LOOP OVER ALL ELEMENTS 
!    FOUND = 0
!
!!        EXTRACT INCIDENCES LIST FOR ELEMENT IE   
!    CALL ELEMENT_NODES (IE, NOD_PER_EL, NODES, ELEM_NODES) 
!
!!        ESTABLISH RANGE OF POSSIBLE ELEMENT NEIGHBORS
!    L_START = HUGE (L_START) ; L_STOP = 0 
!    DO IN = 1, NOD_PER_EL 
!! XXX will get zero for mixed element types XXXX
!      IF ( L_START > L_FIRST (ELEM_NODES (IN) ) ) &
!           L_START = L_FIRST (ELEM_NODES (IN) )
!      IF ( L_STOP  < L_LAST  (ELEM_NODES (IN) ) ) &
!           L_STOP  = L_LAST  (ELEM_NODES (IN) )
!    END DO 
!
!!        LOOP OVER POSSIBLE ELEMENT NEIGHBORS
!    IF ( L_START <= L_STOP) THEN 
!      DO L_TEST = L_START, L_STOP 
!        IF ( L_TEST /= IE) THEN 
!
!!      LOOP OVER INCIDENCES OF POSSIBLE ELEMENT NEIGHBOR        
!          CALL ELEMENT_NODES (L_TEST, NOD_PER_EL, NODES, NEIG_NODES) 
!          DO IN = 1, NOD_PER_EL 
!            N_TEST = NEIG_NODES (IN) 
!            IF ( N_TEST <  1 ) CYCLE ! to a real node
!            IF ( L_FIRST (N_TEST)  > IE ) CYCLE ! to next node
!            IF ( L_LAST  (N_TEST)  < IE ) CYCLE ! to next node
!
!!              COMPARE WITH INCIDENCES OF ELEMENT IE             
!              IF ( ANY ( ELEM_NODES == N_TEST ) ) THEN
!                FOUND = FOUND + 1
!                L_TO_L_NEIGH (FOUND, IE) = L_TEST
!                EXIT ! this L_TEST element search loop
!              END IF
!          END DO ! over N_TEST
!        END IF ! not itself
!        IF ( L_TO_L_SUM (IE) == FOUND ) EXIT ! have last candidate
!      END DO ! over candidate element L_TEST
!    END IF ! a possible candidate
!  END DO ! over all elements
!END SUBROUTINE FORM_ELEMS_AT_EL_OLD                         

SUBROUTINE SAVE_QP_LOOP_COUNT (ELEM_QP)
! ..............................................................
!  SAVE NUMBER OF QUADRATURE POINTS THAT AN ELEMENT IS
!  USING (FOR SCP RECOVERY AND/OR POST-PROCESSING)
! ..............................................................
Use System_Constants ! For U_FLUX unit number
 IMPLICIT NONE
 INTEGER,  INTENT(IN) :: ELEM_QP

! ELEM_QP = NUMBER OF QUADRATURE POINTS USED IN THIS LOOP

  IF ( U_FLUX > 0 ) WRITE (U_FLUX) ELEM_QP
END SUBROUTINE SAVE_QP_LOOP_COUNT 

SUBROUTINE SAVE_QP_LOOP_ARRAYS (XYZ, E, B)
! ..............................................................
!  SAVE ELEM POSITION, CONSTITUTIVE, & OPERATOR MATRICES AT QP
!  USING (FOR SCP RECOVERY AND/OR POST-PROCESSING)
! ..............................................................
Use System_Constants ! For U_FLUX, N_SPACE, N_R_B
Use Elem_Type_Data  ! for: LT_FREE
 IMPLICIT NONE
 REAL(DP), INTENT(IN) :: XYZ (N_SPACE)
 REAL(DP), INTENT(IN) :: E (N_R_B, N_R_B)
 REAL(DP), INTENT(IN) :: B (N_R_B, LT_FREE)

  IF ( U_FLUX > 0 ) WRITE (U_FLUX) XYZ, E, B
END SUBROUTINE SAVE_QP_LOOP_ARRAYS 

SUBROUTINE GET_QP_LOOP_COUNT (ELEM_QP)
! ..............................................................
!  GET NUMBER OF QUADRATURE POINTS THAT AN ELEMENT IS
!  USING (FOR SCP RECOVERY AND/OR POST-PROCESSING)
! ..............................................................
Use System_Constants ! For U_FLUX unit number
 IMPLICIT NONE
 INTEGER,  INTENT(OUT) :: ELEM_QP

! ELEM_QP = NUMBER OF QUADRATURE POINTS USED IN THIS LOOP

  IF ( U_FLUX > 0 ) THEN
    READ (U_FLUX) ELEM_QP
  ELSE
    ELEM_QP = 0
    PRINT *, 'GET_QP_LOOP_COUNT: WARNING, U_FLUX IS ZERO'
    N_WARN = N_WARN + 1
  END IF
END SUBROUTINE GET_QP_LOOP_COUNT 

SUBROUTINE GET_QP_LOOP_ARRAYS (XYZ, E, B)
! ..............................................................
!  GET ELEM POSITION, CONSTITUTIVE, & OPERATOR MATRICES AT QP
!  USING (FOR SCP RECOVERY AND/OR POST-PROCESSING)
! ..............................................................
Use System_Constants ! For U_FLUX, N_SPACE, N_R_B
Use Elem_Type_Data  ! for: LT_FREE
 IMPLICIT NONE
 REAL(DP), INTENT(OUT) :: XYZ (N_SPACE)
 REAL(DP), INTENT(OUT) :: E (N_R_B, N_R_B)
 REAL(DP), INTENT(OUT) :: B (N_R_B, LT_FREE)

  IF ( U_FLUX > 0 ) THEN
    READ (U_FLUX) XYZ, E, B
  ELSE
    XYZ = 0.d0 ; B = 0.d0
    CALL REAL_IDENTITY (N_R_B, E)
    PRINT *, 'GET_QP_LOOP_ARRAYS: WARNING, U_FLUX IS ZERO'
    N_WARN = N_WARN + 1
  END IF
END SUBROUTINE GET_QP_LOOP_ARRAYS 

SUBROUTINE GET_NODE_AT_DOF (N_DOF, N_G_FREE, N_PT)
! ..............................................................
! IDENTIFY NODE POINT, N_PT, ASSOCIATED WITH DEGREE OF FREEDOM
! ..............................................................
 IMPLICIT NONE
 INTEGER, INTENT(IN)  :: N_DOF, N_G_FREE
 INTEGER, INTENT(OUT) :: N_PT
! N_DOF    = EQUATION NUMBER
! N_G_FREE = NUMBER OF UNKNOWNS PER NODE
! N_PT     = GLOBAL NODE NUMBER WHERE EQUATION OCCURS
  N_PT = 1 + (N_DOF -1)/N_G_FREE
END SUBROUTINE GET_NODE_AT_DOF 

SUBROUTINE STORE_FLUX_POINT_COUNT ! Save LT_QP
! ..............................................................
! Store element flux point count
! ..............................................................
Use System_Constants ! For U_FLUX, N_WARN
Use Elem_Type_Data   ! For LT_QP
 LOGICAL, SAVE :: WARN = .FALSE.
 
  IF ( U_FLUX > 0 ) THEN
    WRITE (U_FLUX) LT_QP  ! Declare el flux save loop size
  ELSE
    IF ( NO_SCP_AVE ) RETURN
    IF ( .NOT. WARN ) THEN
      PRINT *,'WARNING: STORE_FLUX_POINT_COUNT, FILE NOT ACTIVE'
      PRINT *,'IS no_scp_ave IN KEYWORD CONTROLS ?'
      N_WARN = N_WARN + 1
      WARN = .TRUE.
    END IF
  END IF
END SUBROUTINE STORE_FLUX_POINT_COUNT ! Save LT_QP

SUBROUTINE STORE_FLUX_POINT_DATA (XYZ, E, B) ! Save LT_QP
! ..............................................................
! Store element flux point data for averaging or post processing
! ..............................................................
Use System_Constants ! For U_FLUX, N_WARN, N_R_B, N_SPACE
Use Elem_Type_Data   ! For LT_FREE
 REAL(DP), INTENT(IN) :: XYZ (N_SPACE), E (N_R_B, N_R_B), &
                         B (N_R_B, LT_FREE)
 LOGICAL, SAVE :: WARN = .FALSE.
 
  IF ( U_FLUX > 0 ) THEN
!-->   WRITE COORDS, E AND DERIVATIVE MATRIX, FOR POST PROCESSING
    WRITE (U_FLUX) XYZ, E, B
  ELSE
    IF ( NO_SCP_AVE ) RETURN
    IF ( .NOT. WARN ) THEN
      PRINT *,'WARNING: STORE_FLUX_POINT_DATA,  FILE NOT ACTIVE'
      PRINT *,'IS no_scp_ave IN KEYWORD CONTROLS ?'
      N_WARN = N_WARN + 1
      WARN = .TRUE.
    END IF
  END IF
END SUBROUTINE STORE_FLUX_POINT_DATA  ! Save LT_QP

SUBROUTINE READ_FLUX_POINT_COUNT (N_IP) ! Save LT_QP
! ..............................................................
! Read element flux point count
! ..............................................................
Use System_Constants ! For U_FLUX, N_WARN
Use Elem_Type_Data   ! For LT_QP
 INTEGER, INTENT(OUT) :: N_IP
 LOGICAL, SAVE :: WARN = .FALSE.
 
  IF ( U_FLUX > 0 ) THEN
    READ  (U_FLUX) N_IP   ! Declare el flux save loop size
  ELSE
    N_IP = 0
    IF ( NO_SCP_AVE ) RETURN
    IF ( .NOT. WARN ) THEN
      PRINT *,'WARNING: READ_FLUX_POINT_COUNT, FILE NOT ACTIVE'
      PRINT *,'IS no_scp_ave IN KEYWORD CONTROLS ?'
      N_WARN = N_WARN + 1
      WARN = .TRUE.
    END IF
  END IF
END SUBROUTINE READ_FLUX_POINT_COUNT ! Save LT_QP

SUBROUTINE READ_FLUX_POINT_DATA (XYZ, E, B) ! Save LT_QP
! ..............................................................
! Read element flux point data for averaging or post processing
! ..............................................................
Use System_Constants ! For U_FLUX, N_WARN, N_R_B, N_SPACE
Use Elem_Type_Data   ! For LT_FREE
 REAL(DP), INTENT(OUT) :: XYZ (N_SPACE), E (N_R_B, N_R_B), &
                          B (N_R_B, LT_FREE)
 LOGICAL, SAVE :: WARN = .FALSE.
 
  IF ( U_FLUX > 0 ) THEN
!-->   READ  COORDS, E AND DERIVATIVE MATRIX, FOR POST PROCESSING
    READ  (U_FLUX) XYZ, E, B
  ELSE
    XYZ = 0.d0 ; E = 0.d0 ; B = 0.d0
    IF ( NO_SCP_AVE ) RETURN
    IF ( .NOT. WARN ) THEN
      PRINT *,'WARNING: READ_FLUX_POINT_DATA,  FILE NOT ACTIVE'
      PRINT *,'IS no_scp_ave IN KEYWORD CONTROLS ?'
      N_WARN = N_WARN + 1
      WARN = .TRUE.
    END IF
  END IF
END SUBROUTINE READ_FLUX_POINT_DATA  ! Save LT_QP

! End File: model_lib.f
