! copyright 2005,  J. E. Akin, all rights reserved.

SUBROUTINE INPUT_CONSTRAINT_EQS (N_RES_EQ, COEF_C_EQ, NDEX_C_EQ)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!    READ NODAL PARAMETER CONSTRAINT EQUATION DATA
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 INTEGER,  INTENT(IN)  :: N_RES_EQ  (MAX_ACT)
 REAL(DP), INTENT(OUT) :: COEF_C_EQ (MAX_ACT, N_CEQ)
 INTEGER,  INTENT(OUT) :: NDEX_C_EQ (MAX_ACT, N_CEQ)
 INTEGER   :: IEQ, IN, NTEST, NEQ, OK
 INTEGER   :: NODE_1, IPAR_1, NODE_2, IPAR_2, NODE_3, IPAR_3
 REAL(DP)  :: A_1, A_2, A_3, A_4

! COEF_C_EQ = (I,J) CONSTRAINT COEFF I OF EQUATION J
! MAX_ACT   = NUMBER OF ACTIVE ESSENTIAL BC'S & CONSTRAINTS, >= 1
! N_G_DOF   = NO. PARAMETERS PER NODE
! NDEX_C_EQ = (I,J) CONSTRAINT DOF NO I OF EQUATION J
! N_RES_EQ  = (I) NUMBER OF CONSTRAINTS OF TYPE I

!==>> DATA TYPE 7,  Essential Boundary Conditions
!
!    Input data for a Type 1 nodal constraint
!
!    NODE_1 = system nodal point number.
!    IPAR_1 = Constrainted parameter number, 1 <= IPAR_1 <= N_G_DOF
!    A_1    = Value assigned to the parameter.
!
!  NOTES *** These constraints are of the form D(J) = A_1, where D(J) 
!  is the system degree of freedom corresponding to parameter IPAR_1 at
!  node NODE_1. That is, J = N_G_DOF*(NODE_1-1)+IPAR_1.  There will be 
!  one TYPE 7 card for each Type 1 constraint flagged in I_BC(NODE_1).
!
!==>> DATA TYPE 8, Type 2 linear constraint, if any. 
!
!    NODE_1 = First system node number.
!    IPAR_1 = Constrained parameter number at NODE_1.
!    NODE_2 = System node number of second node.
!    IPAR_2 = Constrained parameter number at NODE_2.
!    A_1    = Coefficient of first parameter.
!    A_2    = Coefficient of second parameter.
!    A_3    = Coefficient of right hand side.
!
!  NOTES *** These constraints are of the form A_1*D(J1)+A_2*D(J2)=A_3.
!  There will be one TYPE 8 card for each Type 2 constraint equation.
!
!==>> DATA TYPE 9, Type 3 linear constraint, if any. 
!
!    NODE_1 = First node point number.
!    IPAR_1 = First constrained parameter number.
!    NODE_2 = Second nodal point.
!    IPAR_2 = Second constrained parameter number.
!    NODE_3 = Third nodal number.
!    IPAR_3 = Third constrained parameter.
!    A_1    = Coefficient of first parameter.
!    A_2    = Coefficient of second parameter.
!    A_3    = Coefficient of third parameter.
!    A_4    = Constant on right hand side.
!
!  NOTES: Constraints are of the form 
!   A_1 * D(J1) + A_2 * D(J2) + A_3*D(J3) = A_4   There will 
!   be one TYPE 9 card for each Type 3 constraint equation.

  IF (ECHO_BC) WRITE (N_PRT, "(/, '*** CONSTRAINT EQUATION DATA ***')")
  !b IF ( I_BUG > 0 ) PRINT *, 'N_RES_EQ = ', N_RES_EQ
  
!  WARN IF NO CONSTRAINTS
  IF ( MAX_ACT < 1 ) THEN
    WRITE (N_PRT,  * ) 'WARNING, NO CONSTRAINTS  !!!!!!!'
    N_WARN = N_WARN + 1 ! INCREMENT WARNING
    RETURN
  END IF
  IEQ = 0

  DO IN = 1, MAX_ACT  ! LOOP ON ACTIVE CONSTRAINT TYPES
    NTEST = N_RES_EQ (IN)
    IF ( NTEST == 0 ) CYCLE ! to next IN

    SELECT CASE ( IN) ! BRANCH ON CONSTRAINT TYPE

     CASE ( 1 ) !-->         TYPE 1     D(L1) = A_1

       IF (ECHO_BC) WRITE (N_PRT, 10)
       10 FORMAT ('CONSTRAINT TYPE_ONE: (PAR_1 @ NODE_1) =  A_1.' &
              ,/, 'EQ. NO.   NODE_1   PAR_1          A_1')
!   NOTE: Can compute value as a function of coordinates here
       DO NEQ = 1, NTEST
         IEQ = IEQ + 1
         READ (N_CRD, *, IOSTAT = OK ) NODE_1, IPAR_1, A_1
         CALL INCREMENT_INPUT_LINE ! Update position count
         CALL READ_STATUS (OK, NODE_1, IPAR_1)
         IF (ECHO_BC) WRITE (N_PRT, 30) IEQ, NODE_1, IPAR_1, A_1
         30 FORMAT ( 3I7 ,2X, 1PE12.5 )
         !b IF ( NODE_1 > MAX_NP .OR. IPAR_1 > N_G_DOF) THEN
         IF ( NODE_1 > MAX_NP .OR. IPAR_1 > N_G_DOF .OR. &
              NODE_1 < 1      .OR. IPAR_1 < 1      ) THEN
           PRINT *, 'WARNING: IMPOSSIBLE EQUATION NUMBER COMPUTED'
           N_WARN = N_WARN + 1
           PRINT *, 'INPUT_CONSTRAINT_EQS: ERROR, CHECK elems KEYWORD,'
           PRINT *, 'OR THE NODE AND DOF NUMBERS AT THIS LINE'
           CALL CITE_INPUT_LINE ! List current input position
           STOP 'ERROR IN SUBROUTINE INPUT_CONSTRAINT_EQS'
         END IF
         NDEX_C_EQ (1, IEQ) = N_G_DOF * (NODE_1 - 1) + IPAR_1
         COEF_C_EQ (1, IEQ) = A_1
       END DO
       CYCLE ! to next IN

     CASE ( 2 ) !-->         TYPE 2   A_1*D(L1)+A_2*D(L2)=A_3

       IF (ECHO_BC) WRITE (N_PRT, 50)
       50 FORMAT ('CONSTRAINT TYPE_TWO: ',                        &
       'A_1*(PAR_1 @ NODE_1) + A_2*(PAR_2 @ NODE_2) = A_3.' &
       ,/,     'EQ. NO.  NODE_1  PAR_1  NODE_2  PAR_2',       &
       '        A_1        A_2          A_3')
       DO NEQ = 1, NTEST
         IEQ = IEQ + 1
         READ (N_CRD, *, IOSTAT=OK) NODE_1, IPAR_1, NODE_2, IPAR_2, & 
                                    A_1, A_2, A_3
         CALL INCREMENT_INPUT_LINE ! Update position count
         CALL READ_STATUS (OK, NODE_1, IPAR_1)
         IF (ECHO_BC) WRITE (N_PRT, 70) IEQ, NODE_1, IPAR_1, &
                                        NODE_2, IPAR_2,      &
                                        A_1, A_2, A_3
         70 FORMAT (5I5, 3(2X, 1PE10.3))
         IF ( NODE_1 > MAX_NP  .OR. NODE_2 > MAX_NP .OR.        &
              IPAR_1 > N_G_DOF .OR. IPAR_2 > N_G_DOF ) THEN
           CALL CITE_INPUT_LINE ! List current input position
           STOP 'DATA ERROR IN SUBROUTINE INPUT_CONSTRAINT_EQS'
         END IF
         NDEX_C_EQ (1, IEQ) = N_G_DOF * (NODE_1 - 1) + IPAR_1
         NDEX_C_EQ (2, IEQ) = N_G_DOF * (NODE_2 - 1) + IPAR_2
         COEF_C_EQ (1, IEQ) = A_2 / A_1
         COEF_C_EQ (2, IEQ) = A_3 / A_1
       END DO
       CYCLE ! to next IN

     CASE ( 3 ) !        TYPE 3   A_1*D(L1)+A_2*D(L2)+A_3*D(L3)=A_4

       IF (ECHO_BC) WRITE (N_PRT, 80)
       80 FORMAT ('CONSTRAINT TYPE_THREE:', /,                       &
       'A_1*(PAR_1 @ NODE_1) + A_2*(PAR_2 @ NODE_2) + ',             &
       'A_3*DOF(PAR_3 @ NODE_3) = A_4.', /,                          &
       'EQ. NO.   NODE_1   PAR_1   NODE_2   PAR_2   NODE_3   PAR_3', &
       '          A_1          A_2          A_3          A_4')
       DO NEQ = 1, NTEST
         IEQ = IEQ + 1
         READ (N_CRD, *, IOSTAT=OK ) NODE_1, IPAR_1, NODE_2, IPAR_2, &
                          NODE_3, IPAR_3, A_1, A_2, A_3, A_4
         CALL INCREMENT_INPUT_LINE ! Update position count
         CALL READ_STATUS (OK, NODE_1, IPAR_1)
         IF (ECHO_BC) WRITE (N_PRT, 90) IEQ, NODE_1, IPAR_1, &
                                        NODE_2, IPAR_2,      &
                                        NODE_3, IPAR_3,      &
                                        A_1, A_2, A_3, A_4
         90 FORMAT ( 7I5, 4(2X,1PE10.3))
         IF ( NODE_1 > MAX_NP  .OR. NODE_2 > MAX_NP  .OR.            &
              NODE_3 > MAX_NP  .OR. IPAR_3 > N_G_DOF .OR.            &
              IPAR_1 > N_G_DOF .OR. IPAR_2 > N_G_DOF) THEN
           CALL CITE_INPUT_LINE ! List current input position
           STOP 'DATA ERROR IN SUBROUTINE INPUT_CONSTRAINT_EQS'
         END IF
         NDEX_C_EQ (1, IEQ) = N_G_DOF * (NODE_1 - 1) + IPAR_1
         NDEX_C_EQ (2, IEQ) = N_G_DOF * (NODE_2 - 1) + IPAR_2
         NDEX_C_EQ (3, IEQ) = N_G_DOF * (NODE_3 - 1) + IPAR_3
         COEF_C_EQ (1, IEQ) = A_2 / A_1
         COEF_C_EQ (2, IEQ) = A_3 / A_1
         COEF_C_EQ (3, IEQ) = A_4 / A_1
       END DO
       CYCLE ! to next IN
  
     CASE DEFAULT !        OTHER TYPES NOT TREATED
        STOP 'UNSUPPORTED OPTION, INPUT_CONSTRAINT_EQS'
    END SELECT ! BC TYPE
  END DO
END SUBROUTINE INPUT_CONSTRAINT_EQS
                                   
SUBROUTINE INPUT_EXACT_BC_AND_MPC (X, N_RES_EQ, COEF_C_EQ, NDEX_C_EQ)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!    READ NODAL PARAMETER CONSTRAINT EQUATIONS (IF ANY) AND USE
!    THE EXACT SOLUTION TO ENFORCE ESSENTIAL BOUNDARY CONDITIONS
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Select_Source
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: X         (MAX_NP, N_SPACE)
 INTEGER,  INTENT(IN)  :: N_RES_EQ  (MAX_ACT)
 REAL(DP), INTENT(OUT) :: COEF_C_EQ (MAX_ACT, N_CEQ)
 INTEGER,  INTENT(OUT) :: NDEX_C_EQ (MAX_ACT, N_CEQ)

 REAL(DP)  :: XYZ       (N_SPACE)    ! POINT IN SPACE
 REAL(DP)  :: EXACT_SOL (N_G_DOF)    ! PT EXACT VALUE
 INTEGER   :: IEQ, IN, NTEST, NEQ, OK
 INTEGER   :: NODE_1, IPAR_1, NODE_2, IPAR_2, NODE_3, IPAR_3
 REAL(DP)  :: A_1, A_2, A_3, A_4

! COEF_C_EQ = (I,J) CONSTRAINT COEFF I OF EQUATION J
! MAX_ACT   = NUMBER OF ACTIVE ESSENTIAL BC'S & CONSTRAINTS, >= 1
! N_G_DOF   = NO. PARAMETERS PER NODE
! NDEX_C_EQ = (I,J) CONSTRAINT DOF NO I OF EQUATION J
! N_RES_EQ  = (I) NUMBER OF CONSTRAINTS OF TYPE I

!==>> DATA TYPE 7,  Essential Boundary Conditions
!
!    Input data for a Type 1 nodal constraint
!
!    NODE_1 = system nodal point number.
!    IPAR_1 = Constrainted parameter number, 1 <= IPAR_1 <= N_G_DOF
!    A_1    = Value assigned to the parameter (ignored, replaced by
!             the exact value.
!
!  NOTES *** These constraints are of the form D(J) = A_1, where D(J) 
!  is the system degree of freedom corresponding to parameter IPAR_1 at
!  node NODE_1. That is, J = N_G_DOF*(NODE_1-1)+IPAR_1.  There will be 
!  one TYPE 7 card for each Type 1 constraint flagged in I_BC(NODE_1).
!  See  DATA TYPE 5.
!
!==>> DATA TYPE 8, Type 2 linear constraint, if any. 
!
!    NODE_1 = First system node number.
!    IPAR_1 = Constrained parameter number at NODE_1.
!    NODE_2 = System node number of second node.
!    IPAR_2 = Constrained parameter number at NODE_2.
!    A_1    = Coefficient of first parameter.
!    A_2    = Coefficient of second parameter.
!    A_3    = Coefficient of right hand side.
!
!  NOTES *** These constraints are of the form A_1*D(J1)+A_2*D(J2)=A_3.
!  There will be one TYPE 8 card for each Type 2 constraint equation.
!
!==>> DATA TYPE 9, Type 3 linear constraint, if any. 
!
!    NODE_1 = First node point number.
!    IPAR_1 = First constrained parameter number.
!    NODE_2 = Second nodal point.
!    IPAR_2 = Second constrained parameter number.
!    NODE_3 = Third nodal number.
!    IPAR_3 = Third constrained parameter.
!    A_1    = Coefficient of first parameter.
!    A_2    = Coefficient of second parameter.
!    A_3    = Coefficient of third parameter.
!    A_4    = Constant on right hand side.
!
!  NOTES: Constraints are of the form A_1*D(J1)+A_2*D(J2)+A_3*D(J3)=A_4
!  There will be one TYPE 9 card for each Type 3 constraint equation.

  IF (ECHO_BC) WRITE (N_PRT, "(/, '*** CONSTRAINT EQUATION DATA ***')")
  IF ( I_BUG > 0 ) THEN
    PRINT *, 'N_RES_EQ = ', N_RES_EQ
  END IF
!  WARN IF NO CONSTRAINTS
  IF ( MAX_ACT < 1 ) THEN
    WRITE (N_PRT,  * ) 'WARNING, NO CONSTRAINTS  !!!!!!!'
    N_WARN = N_WARN + 1 ! INCREMENT WARNING
    RETURN
  END IF
  IF ( USE_EXACT_BC .AND. (EXACT_REALS == 0)) WRITE (N_PRT,  * ) &
     & 'NOTE: REAL PROPERTIES NOT VAILABLE YET TO EXACT SOLUTIONS'
  IEQ = 0

  DO IN = 1, MAX_ACT  ! LOOP ON ACTIVE CONSTRAINT TYPES
    NTEST = N_RES_EQ (IN)
    IF ( NTEST == 0 ) CYCLE ! to next IN

    SELECT CASE (IN) ! BRANCH ON CONSTRAINT TYPE

     CASE ( 1 ) !-->         TYPE 1     D(L1) = A_1

       IF ( ECHO_BC ) WRITE (N_PRT, 10) ; 10 FORMAT (             &
         'EXACT CONSTRAINT TYPE_ONE: DOF(PAR_1 @ NODE_1) =  A_1.' &
         ,/, 'EQ. NO.   NODE_1   PAR_1    EXACT_A_1')

       DO NEQ = 1, NTEST
         IEQ = IEQ + 1
         READ (N_CRD, *, IOSTAT=OK) NODE_1, IPAR_1, A_1
         CALL INCREMENT_INPUT_LINE ! Update position count
         CALL READ_STATUS (OK, NODE_1, IPAR_1)
         IF ( NODE_1 > MAX_NP .OR. IPAR_1 > N_G_DOF) THEN   
           CALL CITE_INPUT_LINE ! List current input position
           PRINT *,'ERROR: IMPOSSIBLE NODE OR DOF NUMBER'
           PRINT *, NODE_1, ' OR ', IPAR_1
           STOP 'ERROR IN SUBROUTINE INPUT_EXACT_BC_AND_MPC'
         END IF

!           GET THE COORDINATES AND EXACT VALUE
         CALL PT_COORD (NODE_1, N_SPACE, X, XYZ)
         CALL SELECT_EXACT_SOLUTION (XYZ, EXACT_SOL) !b

!           ASSIGN VALUE AND LIST
         A_1 = EXACT_SOL (IPAR_1)
         IF ( ECHO_BC ) WRITE (N_PRT, 30) IEQ, NODE_1, IPAR_1, A_1
         30 FORMAT ( 3I7 ,2X, 1PE12.5 )

!           COPY INTO SYSTEM ARRAYS
         NDEX_C_EQ (1, IEQ) = N_G_DOF * (NODE_1 - 1) + IPAR_1
         COEF_C_EQ (1, IEQ) = A_1
       END DO
       CYCLE ! to next IN

     CASE ( 2 ) !-->         TYPE 2   A_1*D(L1)+A_2*D(L2)=A_3

       IF ( ECHO_BC ) WRITE (N_PRT, 50)
       50 FORMAT ('CONSTRAINT TYPE_TWO: ',                        &
       'A_1*DOF(PAR_1 @ NODE_1) + A_2*DOF(PAR_2 @ NODE_2) = A_3.' &
       ,/,     'EQ. NO.   NODE_1   PAR_1   NODE_2   PAR_2',       &
       '          A_1          A_2            A_3')
       DO NEQ = 1, NTEST
         IEQ = IEQ + 1
         READ (N_CRD, *, IOSTAT=OK) NODE_1, IPAR_1, NODE_2, IPAR_2, &
                                    A_1, A_2, A_3
         CALL INCREMENT_INPUT_LINE ! Update position count
         CALL READ_STATUS (OK, NODE_1, IPAR_1)
         IF ( ECHO_BC ) WRITE (N_PRT, 70) IEQ, NODE_1, IPAR_1, &
                               NODE_2, IPAR_2, A_1, A_2, A_3
         70 FORMAT (5I7, 3(2X, 1PE12.5))
         IF ( NODE_1 > MAX_NP  .OR. NODE_2 > MAX_NP .OR.        &
              IPAR_1 > N_G_DOF .OR. IPAR_2 > N_G_DOF ) THEN
           CALL CITE_INPUT_LINE ! List current input position
           PRINT *,'ERROR: IMPOSSIBLE NODE OR DOF NUMBER'
           PRINT *, NODE_1, ' OR ', IPAR_1
           PRINT *, NODE_2, ' OR ', IPAR_2
           STOP 'DATA ERROR IN SUBROUTINE INPUT_EXACT_BC_AND_MPC'
         END IF 
         NDEX_C_EQ (1, IEQ) = N_G_DOF * (NODE_1 - 1) + IPAR_1
         NDEX_C_EQ (2, IEQ) = N_G_DOF * (NODE_2 - 1) + IPAR_2
         COEF_C_EQ (1, IEQ) = A_2 / A_1
         COEF_C_EQ (2, IEQ) = A_3 / A_1
       END DO
       CYCLE ! to next IN

     CASE ( 3 ) !        TYPE 3   A_1*D(L1)+A_2*D(L2)+A_3*D(L3)=A_4

       IF ( ECHO_BC ) WRITE (N_PRT, 80)
       80 FORMAT ('CONSTRAINT TYPE_THREE:', /,                       &
       'A_1*DOF(PAR_1 @ NODE_1) + A_2*DOF(PAR_2 @ NODE_2) + ',       &
       'A_3*DOF(PAR_3 @ NODE_3) = A_4.', /,                          &
       'EQ. NO.   NODE_1   PAR_1   NODE_2   PAR_2   NODE_3   PAR_3', &
       '          A_1          A_2          A_3          A_4')
       DO NEQ = 1, NTEST
         IEQ = IEQ + 1
         READ (N_CRD, *, IOSTAT=OK) NODE_1, IPAR_1, NODE_2, IPAR_2,  &
                          NODE_3, IPAR_3, A_1, A_2, A_3, A_4
         CALL INCREMENT_INPUT_LINE ! Update position count
         CALL READ_STATUS (OK, NODE_1, IPAR_1)
         IF ( ECHO_BC ) WRITE (N_PRT, 90) IEQ, NODE_1, IPAR_1, &
                               NODE_2, IPAR_2, NODE_3, IPAR_3, &
                               A_1, A_2, A_3, A_4
         90 FORMAT ( 7I7, 4(2X,1PE12.5))
         IF ( NODE_1 > MAX_NP  .OR. NODE_2 > MAX_NP  .OR.            &
              NODE_3 > MAX_NP  .OR. IPAR_3 > N_G_DOF .OR.            &
              IPAR_1 > N_G_DOF .OR. IPAR_2 > N_G_DOF) THEN
           CALL CITE_INPUT_LINE ! List current input position
           PRINT *,'ERROR: IMPOSSIBLE NODE OR DOF NUMBER'
           PRINT *, NODE_1, ' OR ', IPAR_1
           PRINT *, NODE_2, ' OR ', IPAR_2
           PRINT *, NODE_3, ' OR ', IPAR_3
           STOP 'DATA ERROR IN SUBROUTINE INPUT_EXACT_BC_AND_MPC'
         END IF
         NDEX_C_EQ (1, IEQ) = N_G_DOF * (NODE_1 - 1) + IPAR_1
         NDEX_C_EQ (2, IEQ) = N_G_DOF * (NODE_2 - 1) + IPAR_2
         NDEX_C_EQ (3, IEQ) = N_G_DOF * (NODE_3 - 1) + IPAR_3
         COEF_C_EQ (1, IEQ) = A_2 / A_1
         COEF_C_EQ (2, IEQ) = A_3 / A_1
         COEF_C_EQ (3, IEQ) = A_4 / A_1
       END DO
       CYCLE ! to next IN
  
     CASE DEFAULT !        OTHER TYPES NOT TREATED
        CALL CITE_INPUT_LINE ! List current input position
        PRINT *,'ERROR: UNKNOWN CONSTRAINT TYPE = ', IN
        STOP 'UNSUPPORTED OPTION, INPUT_EXACT_BC_AND_MPC'
    END SELECT ! BC TYPE
  END DO
END SUBROUTINE INPUT_EXACT_BC_AND_MPC
                                   
SUBROUTINE INPUT_ELEM_TYPES ()
! * * * * * * * * * * * * * * * * * * * * * * * * * * *
!            SET OR READ ELEMENT TYPE DATA ITEMS
! * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants  ! for MAX_DAT, N_L_TYPE
Use Elem_Type_Data    ! for LT_DATA (MAX_DAT, N_L_TYPE)
 IMPLICIT NONE
 INTEGER :: BIG_QP, EXPECT, IT, LT, OK ! READ STATUS
                  
! MAX_DAT  = Number of items in an element type (8)
! N_L_TYPE = NUMBER OF ELEMENT TYPES
! 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_SHAP  = ELEMENT SHAPE FLAG NUMBER
!   6, LT_FACES = NUMBER OF ELEMENT (FACING) NEIGHBORS 
!   7, LT_SONS  = NUMBER OF SONS FOR H_ADAPTIVITY 
!   8, LT_USER  = APPLICATION DEPENDENT OPTIONAL USER ITEM 

   IF ( TYPES_ALLOCATED ) THEN ! get from multiple keywords  !b
     IF ( ALL (TYPE_NODES == 0 ) ) THEN
       EXPECT = 2 ; IF ( NOD_PER_EL > 0 ) EXPECT = NOD_PER_EL
       PRINT *,'WARNING: NODES PER ELEMENT TYPE DEFAULTED TO ', &
                EXPECT
       PRINT *,'OVERRIDE WITH KEYWORDS type_nodes IN CONTROL,'
       PRINT *,'OR OVERRIDE WITH KEYWORD el_nodes IN CONTROL'
       N_WARN = N_WARN + 1 ; TYPE_NODES = 2
       IF ( NOD_PER_EL < EXPECT ) NOD_PER_EL = EXPECT
     END IF
     IF ( ALL (TYPE_GAUSS == 0 ) ) THEN
       EXPECT = 0 ;  IF ( N_QP > 0 ) EXPECT = N_QP
       PRINT *,'WARNING: GAUSS POINTS PER TYPE DEFAULTED TO ', EXPECT
       PRINT *,'OVERRIDE WITH KEYWORDS type_gauss IN CONTROL,'
       PRINT *,'OR OVERRIDE WITH KEYWORD gauss IN CONTROL'
       N_WARN = N_WARN + 1 ; TYPE_GAUSS = EXPECT 
     END IF
     IF ( ALL (TYPE_GEOM  == 0 ) ) THEN
       PRINT *,'NOTE: GEOMETRY NODES DEFAULTED TO ANALYSIS NODES'
       PRINT *,'OVERRIDE WITH KEYWORDS type_geom IN CONTROL'
       TYPE_GEOM  = TYPE_NODES ; N_GEOM = MAXVAL (TYPE_GEOM)
     END IF
     IF ( ALL (TYPE_SHAPE == 0 ) ) THEN
       EXPECT = 1 ; IF ( L_SHAPE > 0 ) EXPECT = L_SHAPE
       PRINT *,'WARNING: ALL ELEMENT TYPE SHAPES DEFAULTED TO ', EXPECT
       PRINT *,'OVERRIDE WITH KEYWORDS type_shape IN CONTROL,'
       PRINT *,'OR OVERRIDE WITH KEYWORD shape IN CONTROL'
       N_WARN = N_WARN + 1 ; TYPE_SHAPE = EXPECT ; L_SHAPE = EXPECT
     END IF
     IF ( ALL (TYPE_PARM  == 0 ) ) THEN
       PRINT *,'NOTE: ALL TYPE PARAMETERIC SHAPES DEFAULTED'
       PRINT *,'OVERRIDE WITH KEYWORDS type_parm IN CONTROL'
     END IF
     IF ( ALL (TYPE_FACES == 0 ) ) THEN
       PRINT *,'NOTE: ALL ELEMENT TYPE FACE COUNTS DEFAULTED'
     END IF
     IF ( ALL (TYPE_SONS  == 0 ) ) THEN
       PRINT *,'NOTE: ALL ELEMENT TYPE SONS COUNTS DEFAULTED'
     END IF

     DO IT = 1, N_L_TYPE

!        CHECK FACES NEEDED FOR SCP
       SELECT CASE (TYPE_SHAPE (IT))
         CASE (0) ; EXPECT = 1  ! POINT ELEMENT
         CASE (1) ; EXPECT = 2  ! LINE ELEMENT
         CASE (2) ; EXPECT = 3  ! TRIANGULAR ELEMENT 
         CASE (3) ; EXPECT = 4  ! QUADRILATERAL ELEMENT 
         CASE (4) ; EXPECT = 6  ! HEXAHEDRONAL ELEMENT
         CASE (5) ; EXPECT = 4  ! TETRAHEDRONAL ELEMENT
         CASE (6) ; EXPECT = 5  ! WEDGE ELEMENT
         CASE DEFAULT ! USER SUPPLIED
           EXPECT = 2
       !b  IF ( TYPE_FACES (IT) < 2 ) THEN
       !b    PRINT *,'WARNING: INPUT_ELEM_TYPES CHECK NUMBER OF FACES'
       !b    N_WARN = N_WARN + 1
       !b  END IF
       END SELECT
       IF ( TYPE_FACES (IT) /= EXPECT ) THEN
         PRINT *,'NOTE: INPUT_ELEM_TYPES ', IT,' GIVEN ', &
               TYPE_FACES (IT), ', BUT EXPECTED ', EXPECT, ' FACES'
       !b N_WARN = N_WARN + 1
         IF ( TYPE_FACES (IT) == 0 ) THEN
           PRINT *,'WILL USE ', EXPECT, ' FACES'
           TYPE_FACES (IT) = EXPECT
         END IF 
       END IF ! Change faces
       MAX_FACES = MAXVAL (TYPE_FACES)

!        CHECK PARAMETRIC SPACES
       SELECT CASE (TYPE_SHAPE (IT))
         CASE (0) ; EXPECT = 0  ! POINT ELEMENT
         CASE (1) ; EXPECT = 1  ! LINE ELEMENT
         CASE (2) ; EXPECT = 2  ! TRIANGULAR ELEMENT 
         CASE (3) ; EXPECT = 2  ! QUADRILATERAL ELEMENT 
         CASE (4) ; EXPECT = 3  ! HEXAHEDRONAL ELEMENT
         CASE (5) ; EXPECT = 3  ! TETRAHEDRONAL ELEMENT
         CASE (6) ; EXPECT = 3  ! WEDGE ELEMENT
         CASE DEFAULT ! USER SUPPLIED
           EXPECT = 1
           IF ( TYPE_PARM  (IT) < 1 ) THEN
             PRINT *,'WARNING: INPUT_ELEM_TYPES CHECK NUMBER OF PARMS'
             N_WARN = N_WARN + 1
           END IF
       END SELECT
       IF ( TYPE_PARM  (IT) /= EXPECT ) THEN
         PRINT *,'WARNING: INPUT_ELEM_TYPES ', IT,' GIVEN ', &
               TYPE_PARM  (IT), ', BUT EXPECTED ', EXPECT, ' PARMS'
         N_WARN = N_WARN + 1
         IF ( TYPE_PARM  (IT) == 0 ) THEN
           PRINT *,'WILL USE ', EXPECT, ' PARMS'
           TYPE_PARM  (IT) = EXPECT
         END IF 
       END IF ! Change parametric space
       N_PARM = MAXVAL (TYPE_PARM)

!        CHECK CHILDREN COUNTS  
       SELECT CASE (TYPE_SHAPE (IT))
         CASE (0) ; EXPECT = 0  ! POINT ELEMENT
         CASE (1) ; EXPECT = 2  ! LINE ELEMENT
         CASE (2) ; EXPECT = 4  ! TRIANGULAR ELEMENT 
         CASE (3) ; EXPECT = 4  ! QUADRILATERAL ELEMENT 
         CASE (4) ; EXPECT = 16 ! HEXAHEDRONAL ELEMENT
         CASE (5) ; EXPECT = 16 ! TETRAHEDRONAL ELEMENT ! XXX ?
         CASE (6) ; EXPECT = 12 ! WEDGE ELEMENT
         CASE DEFAULT ! USER SUPPLIED
           EXPECT = 1
   !b      IF ( TYPE_SONS  (IT) < 1 ) THEN
   !b        PRINT *,'WARNING: INPUT_ELEM_TYPES CHECK NUMBER OF SONS'
   !b        N_WARN = N_WARN + 1
   !b      END IF
       END SELECT
       IF ( TYPE_SONS  (IT) /= EXPECT ) THEN
   !b    PRINT *,'WARNING: INPUT_ELEM_TYPES ', IT,' GIVEN ', &
   !b          TYPE_SONS  (IT), ', BUT EXPECTED ', EXPECT, ' SONS'
   !b    N_WARN = N_WARN + 1
         IF ( TYPE_SONS  (IT) == 0 ) THEN
   !b      PRINT *,'WILL USE ', EXPECT, ' SONS'
           TYPE_SONS  (IT) = EXPECT
         END IF 
       END IF ! Change children count   

       LT_DATA (1, IT) = TYPE_NODES (IT) ! NODES PER ELEM TYPE
       LT_DATA (2, IT) = TYPE_GAUSS (IT) ! QUADRATURE POINTS
       LT_DATA (3, IT) = TYPE_GEOM  (IT) ! GEOMETRIC CONTROL NODES
       LT_DATA (4, IT) = TYPE_PARM  (IT) ! PARAMETRIC SPACES
       LT_DATA (5, IT) = TYPE_SHAPE (IT) ! SHAPE TYPE
       LT_DATA (6, IT) = TYPE_FACES (IT) ! NEIGHBOR COUNT
       LT_DATA (7, IT) = TYPE_SONS  (IT) ! SONS WHEN DIVIDED
       LT_DATA (8, IT) = TYPE_USER  (IT) ! FOR USER APPLICATION
     END DO
   END IF  !b need an else here??

  IF ( N_L_TYPE == 1 ) THEN

    LT      = 1  !b
   IF ( .NOT. TYPES_ALLOCATED ) THEN ! get from keywords & reads  !b

!    DEFAULT IS SINGLE TYPE ONLY
    LT_USER = 0                    ! FOR USER APPLICATION
    LT_DATA (1, 1) = NOD_PER_EL    ! NODES PER ELEM TYPE
    LT_DATA (2, 1) = N_QP          ! QUADRATURE POINTS
    LT_DATA (3, 1) = NOD_PER_EL    ! GEOMETRIC CONTROL NODES
    SELECT CASE (L_SHAPE)          
      CASE (1)  ! LINE ELEMENT
        LT_DATA (4, 1) = 1         ! PARAMETRIC SPACES
        LT_DATA (6, 1) = 2         ! NEIGHBOR COUNT
        LT_DATA (7, 1) = 2         ! SONS WHEN DIVIDED
      CASE (2)  ! TRIANGULAR ELEMENT
        LT_DATA (4, 1) = 2         ! PARAMETRIC SPACES
        LT_DATA (6, 1) = 3         ! NEIGHBOR COUNT
        LT_DATA (7, 1) = 4         ! SONS WHEN DIVIDED
      CASE (3)  ! QUADRILATERAL ELEMENT
        LT_DATA (4, 1) = 2         ! PARAMETRIC SPACES
        LT_DATA (6, 1) = 4         ! NEIGHBOR COUNT
        LT_DATA (7, 1) = 4         ! SONS WHEN DIVIDED
      CASE (4)  ! HEXAHEDRONAL ELEMENT
        LT_DATA (4, 1) = 3         ! PARAMETRIC SPACES
        LT_DATA (6, 1) = 6         ! NEIGHBOR COUNT
        LT_DATA (7, 1) = 8         ! SONS WHEN DIVIDED
      CASE (5)  ! TETRAHEDRONAL ELEMENT
        LT_DATA (4, 1) = 3         ! PARAMETRIC SPACES
        LT_DATA (6, 1) = 4         ! NEIGHBOR COUNT
        LT_DATA (7, 1) = 8         ! SONS WHEN DIVIDED
      CASE (6)  ! WEDGE ELEMENT
        LT_DATA (4, 1) = 3         ! PARAMETRIC SPACES
        LT_DATA (6, 1) = 5         ! NEIGHBOR COUNT
        LT_DATA (7, 1) = 8         ! SONS WHEN DIVIDED
      CASE DEFAULT ! USER SUPPLIED
        PRINT *, 'WARNING: PARAMETRIC SPACE DEFAULTED TO ', &
                  N_SPACE, ' ; INPUT_ELEM_TYPES'
        PRINT *, 'AND FACING NEIGHBOR COUNT DEFAULTS TO 6'   
        N_WARN = N_WARN + 1 ! INCREMENT WARNING
        LT_DATA (4, 1) = N_SPACE       
        LT_DATA (6, 1) = 6
        LT_DATA (7, 1) = 0
    END SELECT
    LT_DATA (5, 1) = L_SHAPE       ! SHAPE TYPE
    LT_DATA (8, 1) = LT_USER       ! FOR USER APPLICATION

   END IF ! types keywords  !b
   IF (ECHO_EL) WRITE (N_PRT, 20) LT, LT_DATA (1:MAX_DAT, LT)      

  ELSE  ! multiple element types

    BIG_QP = 0
    DO IT = 1, N_L_TYPE
   IF ( TYPES_ALLOCATED ) THEN ! get from multiple keywords  !b

    LT = IT
    LT_DATA (1, IT) = TYPE_NODES (IT) ! NODES PER ELEM TYPE
    LT_DATA (2, IT) = TYPE_GAUSS (IT) ! QUADRATURE POINTS
    LT_DATA (3, IT) = TYPE_GEOM  (IT) ! GEOMETRIC CONTROL NODES
    LT_DATA (4, IT) = TYPE_PARM  (IT) ! PARAMETRIC SPACES
    LT_DATA (5, IT) = TYPE_SHAPE (IT) ! SHAPE TYPE
    LT_DATA (6, IT) = TYPE_FACES (IT) ! NEIGHBOR COUNT
    LT_DATA (7, IT) = TYPE_SONS  (IT) ! SONS WHEN DIVIDED
    LT_DATA (8, IT) = TYPE_USER  (IT) ! FOR USER APPLICATION

!     VALIDATE SOME QUADRATURE DATA
    IF ( TYPE_SHAPE (IT) == 2 ) THEN  ! Tri
      SELECT CASE ( TYPE_GAUSS (IT) )
        CASE (1, 3, 4, 6, 7, 12, 13 ) ! ok
        CASE DEFAULT
          PRINT *,'WARNING: FOR ELEMENT TYPE ', IT, ' (A TRIANGLE),'
          PRINT *,'EXPECTED A GAUSS RULE OF 1, 3, 4, 6, 7, 12, OR 13'
          N_WARN = N_WARN + 1
      END SELECT
    END IF

    IF ( TYPE_SHAPE (IT) == 3 ) THEN  ! Quad
      SELECT CASE ( TYPE_GAUSS (IT) )
        CASE (1, 4, 9, 16, 25, 36 ) ! ok
        CASE DEFAULT
          PRINT *,'WARNING: FOR ELEMENT TYPE ', IT, ' (A QUAD),'
          PRINT *,'EXPECTED A GAUSS RULE OF 1, 4, 9, 16, 25, OR 36'
          N_WARN = N_WARN + 1
      END SELECT
    END IF

    IF ( TYPE_SHAPE (IT) == 4 ) THEN  ! Hex
      SELECT CASE ( TYPE_GAUSS (IT) )
        CASE (1, 8, 27, 64, 125) ! ok
        CASE DEFAULT
          PRINT *,'WARNING: FOR ELEMENT TYPE ', IT, ' (A HEX),'
          PRINT *,'EXPECTED A GAUSS RULE OF 1, 8, 27, 64, OR 125'
          N_WARN = N_WARN + 1
      END SELECT
    END IF

    IF ( TYPE_SHAPE (IT) == 5 ) THEN  ! tet
      SELECT CASE ( TYPE_GAUSS (IT) )
        CASE (1, 4, 5, 10, 11, 15, 24, 31, 45 ) ! ok
        CASE DEFAULT
          PRINT *,'WARNING: FOR ELEMENT TYPE ', IT, ' (A TET), EXPECTED'
          PRINT *,'A GAUSS RULE OF 1, 4, 5, 10, 11, 15, 24, 31, OR 45'
          N_WARN = N_WARN + 1
      END SELECT
    END IF

   ELSE ! read user input from OLD STYLE data  !b

!    READ SEVERAL DIFFERENT TYPES
      READ (N_CRD, *, IOSTAT=OK) LT, LT_N, LT_QP, LT_GEOM, LT_PARM, &
                       LT_SHAP, LT_FACES, LT_SONS, LT_USER
      CALL INCREMENT_INPUT_LINE ! Update position count
      CALL READ_STATUS (OK, LT, 0)
      IF ( LT_QP > BIG_QP ) BIG_QP = LT_QP !b

!       ERROR CHECK EACH INPUT TYPE
      IF ( LT < 1 .OR. LT > N_L_TYPE ) THEN
        WRITE (N_PRT, *) 'INPUT_ELEM_TYPES: ERROR, ELEMENT TYPE OUT', &
        ' OF RANGE OF 1 TO ', N_L_TYPE
        WRITE (N_PRT, *) 'RESET el_types IN KEYWORDS'
        CALL CITE_INPUT_LINE ! List current input position
        STOP  'ELEMENT TYPE OUT OF BOUNTS, INPUT_ELEM_TYPES'
      END IF

      IF ( LT_N > NOD_PER_EL .OR. LT_QP > N_QP .OR. LT_GEOM > N_GEOM &
                             .OR. LT_PARM > N_PARM ) THEN 
        WRITE (N_PRT, *) 'WARNING: ELEMENT TYPE DATA EXCEEDS LIMITS',&
        & ' IN INPUT_ELEM_TYPES TYPE = ', LT
        N_WARN = N_WARN + 1 ! INCREMENT WARNING
      END IF ! OUT OF RANGE

      IF ( LT_SHAP == 1 .AND. (LT_GEOM /= LT_N) ) THEN 
        WRITE (N_PRT, *) 'INPUT_ELEM_TYPES: WARNING, FOR LINE ', &
        'ELEMENTS USUALLY ALL NODES NEEDED FOR GEOMETRY.'
        WRITE (N_PRT, *) 'SET LT_GEOM = ', LT_N, ' FOR TYPE ', IT
        N_WARN = N_WARN + 1 ! INCREMENT WARNING
      END IF ! PROBABLE GEOMETRY ERROR

      IF ( LT_N < 1 ) LT_N = NOD_PER_EL

      IF ( LT_PARM > N_PARM ) THEN
        WRITE (N_PRT, *) 'INPUT_ELEM_TYPES: WARNING, PARAMETRIC ', &
        'SPACE = ', LT_PARM, ' TOO LARGE IN TYPE ', IT
        WRITE (N_PRT, *) 'CHECK shape IN KEYWORDS. USING ', N_PARM
        LT_PARM = N_PARM ; N_WARN = N_WARN + 1 ! INCREMENT WARNING
      END IF ! OUT OF RANGE

      IF ( LT_N > NOD_PER_EL ) THEN
        WRITE (N_PRT, *) 'INPUT_ELEM_TYPES: WARNING, NODES PER ', &
        'ELEMENT = ', LT_N, ' IS TOO LARGE IN ELEMENT TYPE ', IT
        WRITE (N_PRT, *) 'CHECK el_nodes IN KEYWORDS'
        N_WARN = N_WARN + 1 ! INCREMENT WARNING
      END IF ! OUT OF RANGE

      IF ( LT_QP > N_QP ) THEN
        WRITE (N_PRT, *) 'WARNING: INPUT_ELEM_TYPES, QUADRATURE ', &
        'RULE = ', LT_QP, ' IS TOO LARGE IN TYPE ', IT
        WRITE (N_PRT, *) 'RESET gauss IN KEYWORDS, WILL USE ', LT_QP
        N_QP = LT_QP ; N_WARN = N_WARN + 1 ! INCREMENT WARNING
      !b ELSE 
        !b IF ( LT_QP > BIG_QP ) BIG_QP = LT_QP
      END IF ! OUT OF RANGE
      
      IF ( LT_GEOM < 1 ) LT_GEOM = LT_N   !b

!        SET ELEMENT TYPE DATA
      LT_DATA (1, LT) = LT_N
      LT_DATA (2, LT) = LT_QP
      LT_DATA (3, LT) = LT_GEOM
      LT_DATA (4, LT) = LT_PARM
      LT_DATA (5, LT) = LT_SHAP
      LT_DATA (6, LT) = LT_FACES
      LT_DATA (7, LT) = LT_SONS 
      LT_DATA (8, LT) = LT_USER

!        CHECK FACES NEEDED FOR SCP
      SELECT CASE (LT_SHAP)
        CASE (0) ; EXPECT = 1  ! POINT ELEMENT
        CASE (1) ; EXPECT = 2  ! LINE ELEMENT
        CASE (2) ; EXPECT = 3  ! TRIANGULAR ELEMENT 
        CASE (3) ; EXPECT = 4  ! QUADRILATERAL ELEMENT 
        CASE (4) ; EXPECT = 6  ! HEXAHEDRONAL ELEMENT
        CASE (5) ; EXPECT = 4  ! TETRAHEDRONAL ELEMENT
        CASE (6) ; EXPECT = 5  ! WEDGE ELEMENT
        CASE DEFAULT ! USER SUPPLIED
          EXPECT = 2
          IF ( LT_FACES < 2 ) THEN
            PRINT *,'NOTE: INPUT_ELEM_TYPES, CHECK NUMBER OF FACES'
          END IF
      END SELECT
      IF ( LT_FACES /= EXPECT ) THEN
        PRINT *,'NOTE: INPUT_ELEM_TYPES ', LT,' GIVEN ', &
                 LT_FACES, ', BUT EXPECTED ', EXPECT, ' FACES'
        N_WARN = N_WARN + 1
        IF ( LT_FACES == 0 ) THEN
          PRINT *,'WILL USE ', EXPECT, ' FACES'
          LT_DATA (6, LT) = EXPECT
        END IF 
      END IF

    END IF ! keyword inputs !b

!b print *, 'LT, MAX_DAT, IT ', LT, MAX_DAT, IT
      IF (ECHO_EL) WRITE (N_PRT, 20) LT, LT_DATA (1:MAX_DAT, LT)   
      20 FORMAT ( /, 'ELEMENT TYPE NUMBER =', I3, /,       &
      'NUMBER OF NODES PER ELEMENT ..........', I3, /,     &
      'NUMBER OF QUADRATURE POINTS ..........', I3, /,     &
      'NUMBER OF GEOMETRIC CONTROL NODES ....', I3, /,     &
      'NUMBER OF PARAMETRIC DIMENSIONS ......', I3, /,     &
      'ELEMENT SHAPE NUMBER .................', I3, /,     &
      'NUMBER OF FACES ......................', I3, /,     &
      'NUMBER OF CHILDREN ELEMENTS ..........', I3, /,     &
      'APPLICATION DEPENDENT USER FLAG ......', I3)
    END DO ! OVER ELEMENT TYPES
    
    IF ( .NOT. TYPES_ALLOCATED .AND.   &
         (BIG_QP /=  N_QP) ) THEN ! UNDER INTEGRATION ?
    !b IF ( BIG_QP /=  N_QP ) THEN ! UNDER INTEGRATION ?
      PRINT *,'WARNING: INPUT_ELEM_TYPES, GIVEN ', BIG_QP, &
              ', BUT EXPECTED ', N_QP, ' QUADRATURE POINTS'
      N_WARN = N_WARN + 1
    END IF
  END IF ! MULTIPLE TYPES USED
END SUBROUTINE INPUT_ELEM_TYPES
                            
SUBROUTINE INPUT_MESH (X, I_BC, NODES, N_RES_EQ)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!    READ BASIC PROBLEM DATA
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data
 IMPLICIT NONE
 REAL(DP), INTENT(OUT)   :: X    (MAX_NP, N_SPACE)
 INTEGER,  INTENT(OUT)   :: I_BC (MAX_NP), NODES (L_S_TOT, NOD_PER_EL)
 INTEGER,  INTENT(INOUT) :: N_RES_EQ (MAX_TYP) !b

 REAL(DP)          :: X_J     (N_SPACE)         ! Automatic array
 INTEGER           :: NODES_J (NOD_PER_EL)      ! Automatic array
 CHARACTER(LEN=26) :: BC_F                      ! BC FLAG FORMAT
 INTEGER           :: I, J, K, I_BC_J, L_TYPE_J, OK, KOUNT

! MAX_NP     = NUMBER OF NODES IN SYSTEM
! NOD_PER_EL = NUMBER OF NODES PER ELEMENT
! N_ELEMS    = NUMBER OF ELEMENTS IN SYSTEM
! N_G_DOF    = NUMBER OF PARAMETERS (DOF) PER NODE
! N_SPACE    = DIMENSION OF SOLUTION SPACE
! X          = SYSTEM COORDINATES OF ALL NODES
! I_BC       = PACKED NODAL CONSTRAINT INDICATOR
! NODES      = SYSTEM ARRAY OF ELEMENT INCIDENCES

!==>> DATA TYPE 5, Nodal point data.
!                                     
!    J       = System node number.
!    I_BC(J) = Parameters constraint indicators at node J. (Right
! 	       justified) Digit K equals constraint type of
!              parameter K.
!    X(J,K)  = Spatial coordinates (K=1,N_SPACE) of node J
!
!  NOTE *** There will be one TYPE 5 card for each nodal point in the
!  system. Example, I_BC = 12213 and I_BC = 10330 implies 3 Type 1,
!  1 Type 2, 1 Type 3.  (Mesh generation programs are available)
!  This type terminates when J = MAX_NP.
!
!==>> DATA TYPE 6,  Element Connectivity Data
!                                
!    If N_L_TYPE = 1:
!      J          = Element number.
!      NODES(J,K) = Nodal incidences (K=1,NOD_PER_EL) of elem. J
!    If N_L_TYPE > 1:
!      J          = Element number.
!      L_TYPE(J)  = Element type number
!      NODES(J,K) = Nodal incidences (K=1,NOD_PER_EL) of elem. J
!
!  NOTE *** There will be one TYPE 6 card for each elem.

  X = 0 ; I_BC = 0 ; NODES = 0 ; L_TYPE = 0 ; LT_DATA = 0

!      CREATE FORMAT TO DISPLAY ALL ZEROS IN I_BC FLAG
!      ( I7, I(N_G_DOF+2).N_G_DOF, (6F10.4) )
!      1234567         89.12     3456789.123456 = 26 Characters
  WRITE (BC_F, '("( I7, I", I2, ".", I2, ", (6F10.4) )")' ) &
         (N_G_DOF+2), N_G_DOF

!-->   READ NODAL POINT DATA
  IF (ECHO_PTS) WRITE (N_PRT, 500) XYZ_NAME (1:N_SPACE) (1:10) ! 10 of 12 char
  500 FORMAT ( /, '*** NODAL POINT DATA ***',/, &
                  '   NODE, BC_FLAG,', (3A10) )
  KOUNT = 0
  DO I = 1, MAX_NP
    READ  (N_CRD, *, IOSTAT=OK) J, I_BC_J, (X_J (K), K = 1, N_SPACE)
    CALL INCREMENT_INPUT_LINE ! Update position count
    CALL READ_STATUS (OK, J, I_BC_J)
    IF ( J < 1 .OR. J > MAX_NP ) THEN
      WRITE (N_BUG, * ) 'INPUT_MESH: INVALID NODE NUMBER', J
      CALL CITE_INPUT_LINE ! List current input position
      STOP 'INPUT_MESH: INVALID NODE NUMBER'
    ELSE 
      I_BC (J) = I_BC_J
      X (J, :) = X_J (:)
    END IF ! INVALID DATA
    IF (ECHO_PTS) WRITE (N_PRT, BC_F ) &
                  J, I_BC (J), (X (J, K), K = 1, N_SPACE)
    KOUNT = KOUNT + 1
    IF ( J == MAX_NP ) THEN
      IF ( KOUNT /= MAX_NP ) THEN
        PRINT *,'WARNING: READ ONLY ', KOUNT,' VALUES INSTEAD OF ', MAX_NP
        PRINT *,'MOVE THE LAST VALUE TO END OF THE INPUT LIST'
        N_WARN = N_WARN + 1
      END IF
      EXIT ! THE NODE INPUT LOOP
    END IF
  END DO

!      SUMMARIZE NODAL CONSTRAINT DATA !b
  CALL CONSTRAINT_COUNT (N_RES_EQ, I_BC) !b

!-->   READ ELEMENT TYPE DATA, OR DEFAULT THEM IF ONLY ONE
  CALL INPUT_ELEM_TYPES ()

!-->   READ ELEMENT CONNECTIVITY DATA
  KOUNT = 0
  IF ( N_L_TYPE == 1  ) THEN
    L_TYPE = 1 ! initialize array for one type
    IF (ECHO_EL) WRITE (N_PRT, 30) NOD_PER_EL
    30   FORMAT ( /, '*** ELEMENT CONNECTIVITY DATA ***',/,     &
                  'ELEMENT, ', I2, ' NODAL INCIDENCES.')
    DO I = 1, N_ELEMS
      KOUNT = KOUNT + 1
      READ  (N_CRD, *, IOSTAT=OK) J, (NODES_J (K), K = 1, NOD_PER_EL)
      CALL INCREMENT_INPUT_LINE ! Update position count
      CALL READ_STATUS (OK, J, 0)
      IF ( J < 1 .OR. J > N_ELEMS ) THEN
        WRITE (N_BUG, * ) 'INPUT_MESH: INVALID ELEMENT NUMBER ', J
        CALL CITE_INPUT_LINE ! List current input position
        STOP 'INPUT_MESH: INVALID ELEMENT NUMBER'
      ELSE
        NODES (J, :) = NODES_J (:)
      END IF ! INVALID DATA
      IF (ECHO_EL) WRITE (N_PRT, 51) &
                   J, (NODES (J, K), K = 1, NOD_PER_EL)
      51 FORMAT ( I7, (15I6) )
    END DO
  ELSE ! MULTIPLE ELEMENT TYPES
    IF (ECHO_EL) WRITE (N_PRT, 35)
    35 FORMAT ( /, '*** ELEMENT CONNECTIVITY DATA ***',/,     &
                     'ELEMENT, TYPE, NODAL INCIDENCES.')
    DO I = 1, N_ELEMS
      KOUNT = KOUNT + 1
      READ  (N_CRD, *, IOSTAT=OK ) J, L_TYPE_J
      CALL INCREMENT_INPUT_LINE ! Update position count
      CALL READ_STATUS (OK, J, L_TYPE_J)
      IF ( J < 1 .OR. J > N_ELEMS ) THEN
        WRITE (N_BUG, * ) 'INPUT_MESH: INVALID ELEMENT NUMBER ', J
        CALL CITE_INPUT_LINE ! List current input position
        STOP 'INPUT_MESH: INVALID ELEMENT NUMBER'
      ELSE
        IF ( L_TYPE_J > N_L_TYPE ) THEN
          WRITE (N_BUG, * ) 'INPUT_MESH:  ELEMENT ', J,     &
          ' INVALID TYPE ', L_TYPE_J
          CALL CITE_INPUT_LINE ! List current input position
          STOP 'INPUT_MESH: INVALID ELEMENT TYPE'
        END IF ! TYPE CHECK
        IF ( L_TYPE_J == 0 ) L_TYPE_J = 1
        L_TYPE (J) = L_TYPE_J
      END IF ! INVALID DATA
      LT_N = LT_DATA (1, L_TYPE_J)
!b print *,'LT_N, TYPE ', LT_N, L_TYPE_J
!b print *, LT_DATA (1:5,  L_TYPE_J)
      BACKSPACE (N_CRD) !b
      READ  (N_CRD, *, IOSTAT=OK) J, L_TYPE_J,          &
                                  (NODES (J, K), K = 1, LT_N) !b
      CALL INCREMENT_INPUT_LINE ! Update position count
      CALL READ_STATUS (OK, J, L_TYPE_J)
      IF (ECHO_EL) WRITE (N_PRT, 50) &
                   J, L_TYPE_J, (NODES (J, K), K = 1, LT_N) !b
      50 FORMAT ( (15I6) )
    END DO
  END IF
END SUBROUTINE INPUT_MESH

SUBROUTINE SAVE_MESH (X, I_BC, L_TYPE, NODES)
!b SUBROUTINE SAVE_MESH (X, I_BC, NODES)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!          SAVE BASIC PROBLEM DATA FOR PLOTS 
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 REAL(DP), INTENT(IN) :: X    (MAX_NP, N_SPACE)
 INTEGER,  INTENT(IN) :: I_BC (MAX_NP), L_TYPE (L_S_TOT)
 INTEGER,  INTENT(IN) :: NODES (L_S_TOT, NOD_PER_EL)
 INTEGER :: I, IO_1, IO_2
 CHARACTER(LEN=21) :: BC_F               ! BC FLAG FORMAT

! I_BC       = PACKED NODAL CONSTRAINT INDICATOR
! L_TYPE     = ELEMENT TYPE NUMBER
! MAX_NP     = NUMBER OF NODES IN SYSTEM
! NODES      = SYSTEM ARRAY OF ELEMENT INCIDENCES
! NOD_PER_EL = NUMBER OF NODES PER ELEMENT
! N_ELEMS    = NUMBER OF ELEMENTS IN SYSTEM
! N_SPACE    = DIMENSION OF SOLUTION SPACE
! U_PLT1     = UNIT TO OUTPUT NODE COORDINATES FOR PLOTS
! U_PLT2     = UNIT TO OUTPUT ELEMENT CONNECTIVITY FOR PLOTS
! X          = SYSTEM COORDINATES OF ALL NODES

  IF ( U_PLT1 < 1 .OR. U_PLT2 < 1 ) RETURN ! output turned off
  OPEN (U_PLT1, FILE='msh_bc_xyz.tmp', ACTION='WRITE', &
        STATUS='REPLACE', IOSTAT = IO_1)
  OPEN (U_PLT2, FILE='msh_typ_nodes.tmp', ACTION='WRITE', &
        STATUS='REPLACE', IOSTAT = IO_2)
  IF ( IO_1 > 0 .OR. IO_2 > 0 ) THEN
    PRINT *, 'WARNING: FILE OPEN FAILED. SKIPPING SAVE_MESH'
    N_WARN = N_WARN + 1 ; RETURN
  END IF ! OPEN FAILURE

!      CREATE FORMAT TO DISPLAY ALL ZEROS IN I_BC FLAG
!      ( I(N_G_DOF+2).N_G_DOF, (6F10.4) ) ! 21 characters
!b print *,'Draft BC_F '
!b print '("( I", I2, ".", I2, ", 6(1PE13.5))")', (N_G_DOF+2), N_G_DOF
  WRITE (BC_F, '("( I", I2, ".", I2, ", 6(1PE13.5))")' ) &
         (N_G_DOF+2), N_G_DOF
!b print *,'Draft BC_F ', BC_F
  DO I = 1, MAX_NP
    !b WRITE (U_PLT1, 5) I_BC (I), X (I, :)
    !b 5 FORMAT ( I6, 6(1PE13.5) )
    WRITE (U_PLT1, BC_F) I_BC (I), X (I, :) !b
  END DO
  IF ( ANY ( L_TYPE < 1 ) ) THEN
    PRINT *,'WARNING: SAVE_MESH HAS INVALID ELEMENT TYPES'
    IF ( N_MIXED > 0 .OR. N_F_SEG > 0 )            &
      PRINT *,'PROBABLY CALLED BEFORE SEGMENT INPUT'
    PRINT *,'EXAMINE OUTPUT SAVED IN msh_typ_nodes.tmp'
  END IF ! DATA ERROR
  !b DO I = 1, N_ELEMS
  DO I = 1, L_S_TOT
    !B WRITE (U_PLT2, 6) NODES (I, :) 
    WRITE (U_PLT2, 6) L_TYPE (I), NODES (I, :)  !b
    6 FORMAT ( (17I5) ) 
  END DO
  CLOSE (U_PLT1) 
  !b  CLOSE (U_PLT2)
  PRINT *, 'NOTE: NODE FLAG & COORDINATES SAVED TO msh_bc_xyz.tmp'
  PRINT *, 'NOTE: ELEMENT CONNECTIVITY SAVED TO msh_typ_nodes.tmp'
END SUBROUTINE SAVE_MESH 

SUBROUTINE SAVE_RESULTS (DD) 
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!          SAVE NODAL RESULTS DATA FOR PLOTS 
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Interface_Header
 IMPLICIT NONE
 REAL(DP), INTENT(IN) :: DD (N_D_FRE)
 INTEGER :: I, K, IO_1
!          Automatic Arrays
 INTEGER  :: INDEX (N_G_DOF)

! DD      = CALCULATED NODAL PARAMETERS
! MAX_NP  = NUMBER OF NODES IN SYSTEM
! N_D_FRE = NUMBER OF DOF IN THE SYSTEM
! N_G_DOF = NUMBER OF PARAMETERS (DOF) PER NODE
! U_PLT3  = UNIT TO OUTPUT NODE COORDINATES FOR PLOTS

  IF ( U_PLT3 < 1 ) RETURN ! output turned off
  OPEN (U_PLT3, FILE='node_results.tmp', ACTION='WRITE', &
        STATUS='REPLACE', IOSTAT = IO_1)
  IF ( IO_1 > 0 ) THEN
    PRINT *, 'WARNING: FILE OPEN FAILED. SKIPPING SAVE_RESULTS'
    N_WARN = N_WARN + 1 ; RETURN
  END IF ! OPEN FAILURE
  DO I = 1, MAX_NP
    INDEX = GET_INDEX_AT_PT (I)
    WRITE (U_PLT3, 5) (DD (INDEX (K)), K = 1, N_G_DOF)
    5 FORMAT ( (10(1X, 1PE12.5)) )
  END DO
  CLOSE (U_PLT3) ; PRINT *, ' '
  PRINT *, 'NOTE: SOLUTION RESULTS SAVED TO node_results.tmp'
END SUBROUTINE SAVE_RESULTS
                         
SUBROUTINE SAVE_RESULTS_WITH_EXACT (DD, X) 
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!      SAVE NODAL AND EXACT RESULTS DATA FOR PLOTS 
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Select_Source
Use Interface_Header
 IMPLICIT NONE
 REAL(DP), INTENT(IN) :: DD (N_D_FRE)
 REAL(DP), INTENT(IN) :: X  (MAX_NP, N_SPACE)
 INTEGER :: I, K, IO_1
!          Automatic Arrays
 INTEGER  :: INDEX     (N_G_DOF)
 REAL(DP) :: EXACT_SOL (N_G_DOF)

! DD        = CALCULATED NODAL PARAMETERS
! EXACT_SOL = USER SUPPLIED SOLUTION AT A POINT 
! MAX_NP    = NUMBER OF NODES IN SYSTEM
! N_D_FRE   = NUMBER OF DOF IN THE SYSTEM
! N_G_DOF   = NUMBER OF PARAMETERS (DOF) PER NODE
! U_PLT3    = UNIT TO OUTPUT NODE COORDINATES FOR PLOTS
! X         = COORDINATES OF ALL SYSTEM NODES

  IF ( U_PLT3 < 1 ) RETURN ! output turned off
  OPEN (U_PLT3, FILE='node_results.tmp', ACTION='WRITE', &
        STATUS='REPLACE', IOSTAT = IO_1)
  IF ( IO_1 > 0 ) THEN
    PRINT *, 'WARNING: OPEN FAILED. SKIPPING SAVE_RESULTS_WITH_EXACT'
    N_WARN = N_WARN + 1
    RETURN
  END IF ! OPEN FAILURE
  DO I = 1, MAX_NP
    !b CALL EXACT_SOLUTION (X(I, :), EXACT_SOL)
    CALL SELECT_EXACT_SOLUTION (X(I, :), EXACT_SOL) !b
    INDEX = GET_INDEX_AT_PT (I)
    WRITE (U_PLT3, 5) (DD (INDEX (K)), K = 1, N_G_DOF),  &
                      EXACT_SOL
    5 FORMAT ( (10(1X, 1PE12.5)) )
  END DO
  CLOSE (U_PLT3) 
  PRINT *, 'NOTE: FE & EXACT RESULTS SAVED TO node_results.tmp'
END SUBROUTINE SAVE_RESULTS_WITH_EXACT
                         
SUBROUTINE SAVE_EXACT_FLUXES_AT_NODES (X)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!       SAVE EXACT FLUXES TO A FILE
! * * * * * * * * * * * * * * * * * * * * * * * * * *
 Use System_Constants
Use Select_Source
 IMPLICIT NONE
 REAL(DP), INTENT(IN) :: X (MAX_NP, N_SPACE)
 REAL(DP) :: FLUX (N_R_B)
 INTEGER  :: IP, IO_1

! FLUX         = EXACT FLUX COMPONENTS AT A POINT
! MAX_NP       = NUMBER OF SYSTEM NODES
! X            = COORDINATES OF ALL SYSTEM NODES
! U_PLT9       = UNIT TO OUTPUT NODAL GRADIENTS, IF U_SCPR>0

  IF ( .NOT. USE_EXACT_FLUX ) THEN
    PRINT *, 'WARNING, SAVE_EXACT_FLUXES_AT_NODES: ', &
             'NO EXACT SOLUTION USED'
    N_WARN = N_WARN + 1 ; RETURN
  END IF
  IF ( U_PLT9 < 1 ) RETURN  ! output turned off
  OPEN (U_PLT9, FILE='exact_node_flux.tmp', ACTION='WRITE', &
        STATUS='REPLACE', IOSTAT = IO_1)
  IF ( IO_1 > 0 ) THEN
    PRINT *, 'OPEN FAILED FOR exact_node_flux.tmp'
    RETURN
  END IF

  DO IP = 1, MAX_NP
    !b CALL EXACT_FLUX (X (IP, :), FLUX)
    CALL SELECT_EXACT_FLUX (X (IP, :), FLUX) !b
    WRITE (U_PLT9, 5)  FLUX
    5 FORMAT ( (8(1PE13.5)) )
  END DO ! FOR ALL NODES
  CLOSE (U_PLT9)
  PRINT *, 'NOTE: EXACT NODAL FLUXES SAVED TO exact_node_flux.tmp'
END SUBROUTINE SAVE_EXACT_FLUXES_AT_NODES

SUBROUTINE SAVE_EXACT_SOLUTION_AT_NODES (X) 
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!      SAVE NODAL AND EXACT RESULTS DATA FOR PLOTS 
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Select_Source
Use Interface_Header
 IMPLICIT NONE
 REAL(DP), INTENT(IN) :: X  (MAX_NP, N_SPACE)
 INTEGER :: I, IO_1
!          Automatic Arrays
 REAL(DP) :: EXACT_SOL (N_G_DOF)

! EXACT_SOL = USER SUPPLIED SOLUTION AT A POINT 
! MAX_NP    = NUMBER OF NODES IN SYSTEM
! N_D_FRE   = NUMBER OF DOF IN THE SYSTEM
! N_G_DOF   = NUMBER OF PARAMETERS (DOF) PER NODE
! U_PLT8    = UNIT TO OUTPUT NODE SOLUTION FOR PLOTS
! X         = COORDINATES OF ALL SYSTEM NODES

  IF ( U_PLT8 < 1 ) RETURN ! output turned off
  OPEN (U_PLT8, FILE='exact_node_solution.tmp', ACTION='WRITE', &
        STATUS='REPLACE', IOSTAT = IO_1)
  IF ( IO_1 > 0 ) THEN
    PRINT *,'WARNING: OPEN FAILED IN SAVE_EXACT_SOLUTION_AT_NODES'
    N_WARN = N_WARN + 1 ; RETURN
  END IF ! OPEN FAILURE
  DO I = 1, MAX_NP
    !b CALL EXACT_SOLUTION (X (I, :), EXACT_SOL)
    CALL SELECT_EXACT_SOLUTION (X (I, :), EXACT_SOL) !b
    WRITE (U_PLT8, 5) EXACT_SOL
    5 FORMAT ( (10(1X, 1PE12.5)) )
  END DO
  CLOSE (U_PLT8) 
  PRINT *, 'NOTE: EXACT NODAL VALUES SAVED TO exact_node_solution.tmp'
END SUBROUTINE SAVE_EXACT_SOLUTION_AT_NODES
                         
SUBROUTINE INPUT_PROPERTIES ()
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!  INPUT NODAL POINT, ELEMENT, BOUNDARY SEGMENT,
!     AND MISCELLANEOUS SYSTEM PROPERTIES
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Sys_Properties_Data
 IMPLICIT NONE

! MAX_NP   = NUMBER OF SYSTEM NODES
! N_F_SEG  = NUMBER OF BOUNDARY SEGMENTS WITH FLUX
! N_ELEMS  = NUMBER OF ELEMENTS IN SYSTEM
! FLT_NP, FLT_EL, FLT_MISC, FLT_SP = REAL PROPERTIES OF
!            SYSTEM NODES, ELEMENTS, SEGMENTS, AND MISC.
! NP_FIX, LP_FIX, MISC_FIX, SP_FIX = INTEGER PROPERTIES OF
!            SYSTEM NODES, ELEMENTS, SEGMENTS, AND MISC.
! L_HOMO   = 1, IF ALL ELEMENTS & SEGMENTS SAME AS FIRST
! N_HOMO   = 1, IF ALL NODES SAME AS FIRST ONE
! N_BS_FLO = NUMBER OF REAL SEGMENT PROPERTIES
! N_NP_FLO = NUMBER OF REAL NODAL PROPERTIES
! N_BS_FIX = NUMBER OF INTEGER SEGMENT PROPERTIES
! N_NP_FIX = NUMBER OF INTEGER NODAL PROPERTIES
! N_LP_FIX = NUMBER OF INTEGER ELEMENT PROPERTIES
! N_LP_FLO = NUMBER OF REAL ELEMENT PROPERTIES
! MISC_FL  = NUMBER OF MISC. REAL SYSTEM PROPERTIES
! MISC_FX  = NUMBER OF MISC. INTEGER SYSTEM PROPERTIES
!
!  --------------------------------------------------------------
!  GENERAL NOTE: If a property set is flagged as non-homogeneous
!  (with ?_HOMO = 1) the array is initially filled with copies of 
!  the properties of the first entry. Then only different entries
!  and the last one need be input.
!  --------------------------------------------------------------
!
!==>> DATA TYPE 10, Integer nodal properties.
!  Utilized only when N_NP_FIX > 0 on DATA TYPE 3. 
!    J           = Nodal point number.
!    NP_FIX(J,K) = Integer properties. K=1,N_NP_FIX
!  NOTE *** Reading terminated when J=M, unless N_HOMO = 1. Then 
!   only J=1 is read.  
!
!==>> DATA TYPE 11, Real nodal properties.
!  Utilized only when N_NP_FLO > 0 on DATA TYPE 3. 
!    J           = Nodal point number.
!    FLT_NP(J,K) = real properties. K=1,N_NP_FLO
!  NOTE *** Reading terminated when J=M, unless N_HOMO = 1. Then only
!  J=1 is read.  Omitted values default to zero when N_HOMO=0.

   IF ( N_NP_FIX > 0 .OR. N_NP_FLO > 0 ) CALL   &
     READ_NODAL_POINT_PROPERTIES
!
!==>> DATA TYPE 12,   Integer element properties.
!  Utilizied only when N_LP_FIX > 0 on DATA TYPE 3. 
!    J           = Element number.
!    LP_FIX(J,K) = Integer properties. K=1,N_LP_FIX
!  NOTE *** Reading terminated when J=N_ELEMS, unless L_HOMO = 1. Then
!  only J=1 is read.  Omitted values default to zero when L_HOMO=0.
!
!==>> DATA TYPE 13, Real element properties.
!  Utilized only when N_LP_FLO > 0 on DATA TYPE 3. 
!    J           = Element number.
!    FLT_EL(J,K) = real properties. K=1,N_LP_FLO
!  NOTE *** Reading terminated when J=N_ELEMS, unless L_HOMO = 1. Then
!  only J=1 is read.  
!
   IF ( N_LP_FIX > 0 .OR. N_LP_FLO > 0 ) CALL READ_ELEMENT_PROPERTIES

!==>> DATA TYPE 14,   Integer Mixed_BC  segment properties.
!  Utilizied only when N_MX_FIX > 0 in keywords
!    J           = Segment number.
!    MX_FIX(J,K) = Integer properties. K=1,N_MX_FIX
! NOTE *** Reading terminated when J=N_MIXED, unless L_HOMO = 1.
!   Then only J=1 is read.
!
!==>> DATA TYPE 15,   Real Mixed_BC segment properties.
!  Utilizied only when N_MX_FLO > 0 in keywords
!    J           = Segment number.
!    FLT_MX(J,K) = real properties. K=1,N_MX_FLO
! NOTE *** Reading terminated when J=N_MIXED, unless L_HOMO = 1.
!   Then only J=1 is read.

    IF ( N_MX_FIX > 0 .OR. N_MX_FLO > 0 ) CALL READ_MIXED_BC_PROPERTIES

!==>> DATA TYPE 16,   Integer flux segment properties.
!  Utilizied only when N_BS_FIX > 0 on DATA TYPE 3. 
!    J           = Segment number.
!    SP_FIX(J,K) = Integer properties. K=1,N_BS_FIX
!  NOTE *** Reading terminated when J=N_F_SEG, unless L_HOMO = 1. Then 
!  only J=1 is read.  
!
!==>> DATA TYPE 17,   Real flux segment properties.
!  Utilized only when N_BS_FLO > 0 on DATA TYPE 3. 
!    J           = Segment number.
!    FLT_SP(J,K) = real properties. K=1,N_BS_FLO
!  NOTE *** Reading terminated when J=N_F_SEG, unless L_HOMO = 1. Then 
!  only J=1 is read.  
!

     IF ( N_BS_FIX > 0 .OR. N_BS_FLO > 0 ) CALL    &
       READ_FLUX_SEGMENT_PROPERTIES

!==>> DATA TYPE 18, Integer miscellaneous properties.
!  Utilized only if MISC_FX > 0 on DATA TYPE 3. 
!    MISC_FIX(K) = Misc. fixed point properties. K=1,MISC_FX
!
!==>> DATA TYPE 19,   Real miscellaneous properties.
!  Utilized only when MISC_FL > 0 on DATA TYPE 3. 
!    FLT_MISC(K) = Misc. floating point properties. K=1,MISC_FL
!

   IF ( MISC_FX > 0 .OR. MISC_FL > 0 ) CALL  &
     READ_MISCELLANEOUS_PROPERTIES

!==>> DATA TYPE 20, Real material properties.
!  Utilized only when N_MAT and  MAT_FLO > 0 in keywords.
!    J             = Material number.
!    FLT_MAT(J,K) = real properties. K=1,MAT_FLO
!  NOTE *** Reading terminated when J=N_MAT.

  IF ( MAT_FLO > 0 ) CALL READ_MATERIAL_PROPERTIES

! List any system wide properties set by control wordss XXX
  IF ( GLOBAL_PROPERTY ) THEN 
    IF ( AREA_THICK /= 1.d0 ) PRINT *, 'NOTE: 2-D DOMAIN ', &
      & 'THICKNESS SET TO ', AREA_THICK 
    IF ( ANY ( BODY_FORCE (1:N_SPACE) /= 0.d0 ) ) PRINT *,'NOTE: ',&
      & 'BODY_FORCE VECTOR = ', BODY_FORCE (1:N_SPACE)
  END IF

END SUBROUTINE INPUT_PROPERTIES

SUBROUTINE INPUT_REMARKS (N_SAY)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! READ AND PRINT N_SAY DATAS FOR USER REMARKS
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 INTEGER, INTENT(IN) :: N_SAY
 CHARACTER(4) :: REMARK (16)
 INTEGER      :: I ! loop

!==>> DATA TYPE 4,    User comments.
!    Utilized only when N_SAY > 0 on DATA TYPE 2. (20A4)
!    REMARK = Additional information to be printed with problem
!    output. Used to describe meaning of prameters, properties.
!  NOTE *** There will be N_SAY TYPE 4 cards.
!
  WRITE (6, * ) ' '
  WRITE (6, * ) 'THE NEXT ', N_SAY, ' LINES ARE USER REMARKS'
  DO I = 1, N_SAY
    READ  (5, 10) REMARK    ; 10 FORMAT (16A4)
    WRITE (6, 11) I, REMARK ; 11 FORMAT (I2, 1X, 16A4)
    CALL INCREMENT_INPUT_LINE ! Update position count
  END DO
END SUBROUTINE INPUT_REMARKS
                               
SUBROUTINE INPUT_RHS (CC)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! INPUT SPECIFIED VALUES IN FORCING VECTOR, CC
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Interface_Header
 IMPLICIT NONE
 REAL(DP), INTENT(OUT) :: CC (N_D_FRE)
!           Automatic Arrays
 REAL(DP) :: TOTAL (N_G_DOF), SUM_POS (N_G_DOF), SUM_NEG (N_G_DOF)
 REAL(DP) :: VALUE
 INTEGER  :: I, NODE, I_PARM, J, OK

! 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

!==>> DATA TYPE 18, Initial system column matrix
!  Utilized only if IN_RHS > 0 on DATA TYPE 2.  
!    NODE   = System node number of the point.
!    I_PARM = Parameter number associated with the value.
!    VALUE  = Specified value to be added to column matrix.
!
!  NOTES *** Addition is of the form of C(J) = C(J) + VALUE where
!  J=N_G_DOF*(NODE-1)+I_PARM.  Reading is terminated when NODE=M and
!  I_PARM=N_G_DOF.  Omitted values default to zero.

  CC = 0.d0 ; TOTAL = 0.d0 ; SUM_POS = 0.d0 ; SUM_NEG = 0.d0 ! Initial
  IF (ECHO_RHS) WRITE (N_PRT, &
     & "(/, '*** INITIAL FORCING VECTOR DATA ***',/,       &
     &      '     NODE   PARAMETER     VALUE    EQUATION')")
  DO I = 1, N_D_FRE
    READ (N_CRD, *, IOSTAT=OK ) NODE, I_PARM, VALUE
    CALL INCREMENT_INPUT_LINE ! Update position count
    CALL READ_STATUS (OK, NODE, I_PARM)
    IF ( OK > 0 ) PRINT *,'AT LOAD FOR EQUATION ', I

!    FIND CORRESPONDING DEGREE OF FREEDOM NUMBER
    J      = GET_PT_DOF_INDEX (NODE, I_PARM) 
    CC (J) = VALUE
    !b TOTAL (I_PARM) = TOTAL (I_PARM) + VALUE
    IF ( VALUE >= 0.d0 ) THEN
      SUM_POS (I_PARM) = SUM_POS (I_PARM) + VALUE
    ELSE
      SUM_NEG (I_PARM) = SUM_NEG (I_PARM) + VALUE
    END IF

!   LIST INPUT DATA
    IF (ECHO_RHS) WRITE (N_PRT, 20) NODE, I_PARM, VALUE, J
    20 FORMAT ( 2I10, 1X, 1PE13.5, I8)
    IF ( NODE > MAX_NP .OR. I_PARM > N_G_DOF) THEN
      CALL CITE_INPUT_LINE ! List current input position
      STOP 'NODE OR I_PARM DATA ERROR, INPUT_RHS'
    END IF
    IF ( J == N_D_FRE) EXIT ! THIS LOOP
  END DO
  IF (ECHO_RHS) WRITE (6, 8) 
  8 FORMAT ('*** RESULTANTS ***' ,/, &
            'COMPONENT        SUM        POSITIVE      NEGATIVE')
  TOTAL = SUM_POS + SUM_NEG
  DO J = 1, N_G_DOF
    IF (ECHO_RHS) WRITE (6, 7) INPUT_NAME (J), TOTAL (J), &
                               SUM_POS (J), SUM_NEG (J)
    !b 7 FORMAT ( I3, 2X, 1PE12.4, 2X, 1PE12.4, 2X, 1PE12.4 )
    7 FORMAT (A12, 1PE12.4, 2X, 1PE12.4, 2X, 1PE12.4 )
  END DO
END SUBROUTINE INPUT_RHS
                        
SUBROUTINE INPUT_DAMPENER (DAMPENER)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! INPUT SPECIFIED VALUES FOR POINT DAMPENER
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Interface_Header
 IMPLICIT NONE
 REAL(DP), INTENT(OUT) :: DAMPENER (N_D_FRE)
!           Automatic Arrays
 REAL(DP) :: VALUE, TOTAL (N_G_DOF)
 INTEGER  :: I, NODE, I_PARM, J, OK

! N_D_FRE  = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM
! N_G_DOF  = NUMBER OF PARAMETERS PER NODE
! DAMPENER = POINT DAMPENER ADDED TO ASSEMBLED DAMPING MATRIX
! MAX_NP   = TOTAL NUMBER OF NODES IN THE SYSTEM

!==>> DATA TYPE 18a, Initial system diiagonal mass matrix
!  Utilized only if PT_DAMPENER IS TRUE
!    NODE   = System node number of the point.
!    I_PARM = Parameter number associated with the value.
!    VALUE  = Specified value to be added to column matrix.
!
!  NOTES *** Addition is DAMPENER(J) = DAMPENER(J) + VALUE where
!  J=N_G_DOF*(NODE-1)+I_PARM.  Reading is terminated when NODE=M and
!  I_PARM=N_G_DOF.  Omitted values default to zero.

  DAMPENER = 0.d0 ; TOTAL = 0.d0
  IF (ECHO_RHS) WRITE (N_PRT, &
     & "(/, '*** INITIAL POINT DAMPENER DATA ***',/,       &
     &      '     NODE   PARAMETER   DAMPENER   EQUATION')")
  DO I = 1, N_D_FRE
    READ (N_CRD, *, IOSTAT=OK ) NODE, I_PARM, VALUE
    CALL INCREMENT_INPUT_LINE ! Update position count
    CALL READ_STATUS (OK, NODE, I_PARM)
    IF ( OK > 0 ) PRINT *,'AT LOAD FOR EQUATION ', I

!    FIND CORRESPONDING DEGREE OF FREEDOM NUMBER
    J      = GET_PT_DOF_INDEX (NODE, I_PARM) 
    DAMPENER (J) = VALUE
    IF ( VALUE < 0.d0 ) THEN
      PRINT *,'WARNING: NEGATIVE DAMPENER ARE USUALLY A FATAL ERROR'
      N_WARN = N_WARN + 1
    END IF

!   LIST INPUT DATA
    IF (ECHO_RHS) WRITE (N_PRT, 20) NODE, I_PARM, VALUE, J
    20 FORMAT ( 2I10, 1X, 1PE13.5, I8)
    IF ( NODE > MAX_NP .OR. I_PARM > N_G_DOF) THEN
      CALL CITE_INPUT_LINE ! List current input position
      STOP 'NODE OR I_PARM DATA ERROR, INPUT_DAMPENER'
    END IF
    TOTAL (I_PARM) = TOTAL (I_PARM) + VALUE
    IF ( J == N_D_FRE) EXIT ! THIS LOOP
  END DO ! over all dof
  IF (ECHO_RHS) WRITE (6, 8) 
  8 FORMAT ('*** INPUT POINT DAMPENER TOTALS ***' ,/, &
            'COMPONENT        SUM        ')
  DO J = 1, N_G_DOF
    IF (ECHO_RHS) WRITE (6, 7) INPUT_NAME (J), TOTAL (J)
    7 FORMAT (A12, 1PE12.4)
  END DO ! over components
END SUBROUTINE INPUT_DAMPENER
                        
SUBROUTINE INPUT_MASSES (MASSES)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! INPUT SPECIFIED VALUES FOR POINT MASSES
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Interface_Header
 IMPLICIT NONE
 REAL(DP), INTENT(OUT) :: MASSES (N_D_FRE)
!           Automatic Arrays
 REAL(DP) :: VALUE, TOTAL (N_G_DOF)
 INTEGER  :: I, NODE, I_PARM, J, OK

! N_D_FRE = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM
! N_G_DOF = NUMBER OF PARAMETERS PER NODE
! MASSES  = POINT MASSES TO BE ADDED TO ASSEMBLED MASS MATRIX
! MAX_NP  = TOTAL NUMBER OF NODES IN THE SYSTEM

!==>> DATA TYPE 18a, Initial system diiagonal mass matrix
!  Utilized only if PT_MASSES IS TRUE
!    NODE   = System node number of the point.
!    I_PARM = Parameter number associated with the value.
!    VALUE  = Specified value to be added to column matrix.
!
!  NOTES *** Addition is MASSES(J) = MASSES(J) + VALUE where
!  J=N_G_DOF*(NODE-1)+I_PARM.  Reading is terminated when NODE=M and
!  I_PARM=N_G_DOF.  Omitted values default to zero.

  MASSES = 0.d0 ; TOTAL = 0.d0
  IF (ECHO_RHS) WRITE (N_PRT, &
     & "(/, '*** INITIAL POINT MASSES DATA ***',/,       &
     &      '     NODE   PARAMETER     MASS     EQUATION')")
  DO I = 1, N_D_FRE
    READ (N_CRD, *, IOSTAT=OK ) NODE, I_PARM, VALUE
    CALL INCREMENT_INPUT_LINE ! Update position count
    CALL READ_STATUS (OK, NODE, I_PARM)
    IF ( OK > 0 ) PRINT *,'AT LOAD FOR EQUATION ', I

!    FIND CORRESPONDING DEGREE OF FREEDOM NUMBER
    J      = GET_PT_DOF_INDEX (NODE, I_PARM) 
    MASSES (J) = VALUE
    IF ( VALUE < 0.d0 ) THEN
      PRINT *,'WARNING: NEGATIVE MASSES ARE USUALLY A FATAL ERROR'
      N_WARN = N_WARN + 1
    END IF

!   LIST INPUT DATA
    IF (ECHO_RHS) WRITE (N_PRT, 20) NODE, I_PARM, VALUE, J
    20 FORMAT ( 2I10, 1X, 1PE13.5, I8)
    IF ( NODE > MAX_NP .OR. I_PARM > N_G_DOF) THEN
      CALL CITE_INPUT_LINE ! List current input position
      STOP 'NODE OR I_PARM DATA ERROR, INPUT_MASSES'
    END IF
    TOTAL (I_PARM) = TOTAL (I_PARM) + VALUE
    IF ( J == N_D_FRE) EXIT ! THIS LOOP
  END DO ! over all dof
  IF (ECHO_RHS) WRITE (6, 8) 
  8 FORMAT ('*** INPUT POINT MASS TOTALS ***' ,/, &
            'COMPONENT        SUM        ')
  DO J = 1, N_G_DOF
    IF (ECHO_RHS) WRITE (6, 7) INPUT_NAME (J), TOTAL (J)
    7 FORMAT (A12, 1PE12.4)
  END DO ! over components
END SUBROUTINE INPUT_MASSES
                        
SUBROUTINE INPUT_FLUX_SEG  (L_TYPE, NODES) 
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!    READ FLUX BOUNDARY DATA & STORE THEM
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Sys_Properties_Data   ! for S_FLUX (L_B_N, N_G_FLUX, N_F_SEG)
 IMPLICIT NONE
 INTEGER,  INTENT(INOUT) :: L_TYPE (L_S_TOT)
 INTEGER,  INTENT(INOUT) :: NODES  (L_S_TOT, NOD_PER_EL)
!       Automatic Arrays
 REAL(DP) :: FLUX (L_B_N, N_G_FLUX)    !b (N_R_B) ?
 INTEGER  :: ELEM_NODES (L_B_N)
 INTEGER  :: J_SEG, J_TYPE, L, ISEG, K, IS, LAST, L_B_T, OK

! FLUX     = COMPONENTS OF GIVEN BOUNDARY FLUX ON THE SEGMENT
! LAST     = LAST SQ MATRIX ELEMENT, N_ELEMS + N_MIXED
! L_S_TOT  = TOTAL NUMBER OF ELEMENTS & THEIR SEGMENTS
! L_TYPE   = ELEMENT TYPE NUMBER
! L_B_N    = MAXIMUM NODES ON AN ELEMENT BOUNDARY SEGMENT
! L_B_T    = NO. OF NODES ON AN ELEMENT BOUNDARY SEGMENT TYPE
! MAX_NP   = NO. OF SYSTEM NODES,
! N_D_FLUX = L_B_N * N_G_FLUX NUMBER OF SEGMENT FLUXES
! N_ELEMS  = NUMBER OF ELEMENTS
! N_F_SEG  = NO. OF ELEMENT BOUNDARY SEGMENTS IN SYSTEM
! N_G_FLUX = NO. OF FLUX COMPONENTS PER NODE
! N_G_DOF  = NO. OF DOF PER NODE
! N_MIXED  = NUMBER OF SEGMENTS WITH MIXED BOUNDARY CONDITIONS
! N_SPACE  = DIMENSION OF SPACE
! S_FLUX   = SPECIFIED COMPONENTS OF FLUX AT NODES ON SEGMENTS

!=>> DATA TYPE 19, Specified boundary flux.
!  Utilized only when N_F_SEG > 0 and L_B_N > 0 
!  TYPE 19A  
!    If NL_TYPE = 1, (single element type, usually wrong here):
!      J             = Segment number, 1 <= J <= N_F_SEG
!      ELEM_NODES(J) = List of nodes on the segment. J=1,L_B_N.
!    
!     Else if NL_TYPE > 1, (allows a mixture of element types):
!       J             = Segment number,      1 <= J      <= N_F_SEG.
!       L_TYPE        = Segment type number, 1 <= L_TYPE <= NL_TYPE.
!       ELEM_NODES(J) = List of nodes on the segment. J=1,L_B_N.
!  TYPE 19B 
!    (FLUX(K,L),L=1,N_G_FLUX = The N_G_FLUX components of flux at each
!                              of the above nodes. K=1,L_B_N.
!
!  NOTES *** These data are related to problem dependent calculations.
!  Thus the user assigns meanings to the nodal order and component
!  order.  There will be N_F_SEG TYPE 19 DATA SETs.

  IF ( N_L_TYPE == 1 ) THEN 
   WRITE (N_PRT, "('*******', /,                                     &
   &'WARNING: SEGMENTS ALMOST ALWAYS REQUIRE MULTIPLE ELEMENT TYPES',&
   &/, '*******')")
    N_WARN = N_WARN + 1 ! INCREMENT WARNING
  END IF
  IF (ECHO_FLUX) WRITE (N_PRT, "(/, '*** ELEMENT BOUNDARY FLUXES ***',/)")
  IF ( FLUX_NORMAL ) THEN ! CONSTANT NORMAL FLUX ON ALL SEGMENTS
    WRITE (N_PRT, "('NOTE: CONSTANT NORMAL FLUX ON ALL SEGMENTS')")
    WRITE (N_PRT, "('USING A VALUE OF ', 1PE12.4 )") Q_NORMAL_SEG  
    IF (ECHO_FLUX) WRITE (N_PRT, "('SEGMENT, TYPE, ',I2,       &
                          & ' MAX NODES ON THE SEGMENT')") L_B_N
  ELSE ! NEED VARIABLE FLUX DATA
    IF (ECHO_FLUX) THEN
      WRITE (N_PRT, "('SEGMENT, TYPE, ',I2,      &
             & ' MAX NODES ON THE SEGMENT')")  L_B_N
      WRITE (N_PRT, "('SEGMENT, NODE, ',I2,                 &
             &      ' FLUX COMPONENTS PER NODE')") N_G_FLUX
    END IF ! ECHO
  END IF ! CONSTANT FLUX

  LAST = N_ELEMS + N_MIXED ; J_TYPE = 1 ; L_B_T = L_B_N   ! INITIALIZE
  ELEM_NODES = 0
  DO ISEG = 1, N_F_SEG

!-->     READ BOUNDARY NODES
    IF ( N_L_TYPE == 1  ) THEN
      READ (N_CRD, *, IOSTAT=OK) J_SEG, ELEM_NODES (1:L_B_N) !b
      CALL INCREMENT_INPUT_LINE ! Update position count
      CALL READ_STATUS (OK, J_SEG, 0)
    ELSE
      !b READ (N_CRD, *, IOSTAT=OK) J_SEG, J_TYPE, ELEM_NODES (1:L_B_N)
      READ (N_CRD, *, IOSTAT=OK) J_SEG, J_TYPE !b
      CALL READ_STATUS (OK, J_SEG, J_TYPE) !b
      L_B_T = TYPE_NODES (J_TYPE) !b
      READ (N_CRD, *, IOSTAT=OK) J_SEG, J_TYPE, ELEM_NODES (1:L_B_T) !b
      CALL INCREMENT_INPUT_LINE ! Update position count
      CALL READ_STATUS (OK, J_SEG, J_TYPE)
    END IF
    IF ( J_SEG < 1 .OR. J_SEG > N_F_SEG ) THEN         
      PRINT *,'INPUT_FLUX_SEG: INVALID SEGMENT NUMBER ', J_SEG
      CALL CITE_INPUT_LINE ! List current input position
      STOP 'INVALID SEGMENT NUMBER IN INPUT_FLUX_SEG'
    END IF
    IF ( J_TYPE < 1 .OR. J_TYPE > N_L_TYPE ) THEN         
      PRINT *,'INPUT_FLUX_SEG: INVALID ELEMENT TYPE ', J_TYPE
      CALL CITE_INPUT_LINE ! List current input position
      STOP 'INVALID ELEMENT TYPE NUMBER IN INPUT_FLUX_SEG'
    END IF
    IF ( MINVAL(ELEM_NODES(1:L_B_T)) < 1 .OR.       &
         MAXVAL(ELEM_NODES(1:L_B_T)) > MAX_NP ) THEN         
      PRINT *,'INPUT_FLUX_SEG: INVALID NODES AT SEGMENT ', J_SEG
      CALL CITE_INPUT_LINE ! List current input position
      STOP 'INVALID NODE NUMBER IN INPUT_FLUX_SEG'
    END IF
    IF (ECHO_FLUX) WRITE (N_PRT, 30) J_SEG, J_TYPE, ELEM_NODES (1:L_B_T)
    30 FORMAT ( I7, I6, 3X, (16I5) )

!-->     READ BOUNDARY FLUX
    IF ( .NOT. FLUX_NORMAL ) THEN ! NEED VARIABLE DATA
      READ (N_CRD, *, IOSTAT=OK) ((FLUX (K, L), L = 1, N_G_FLUX), &
                                 K = 1, L_B_T)
      CALL INCREMENT_INPUT_LINE ! Update position count
      IF ( OK /= 0 ) WRITE (N_PRT, *) 'ERROR READING BOUNDARY FLUX'
      IF (ECHO_FLUX)  THEN
        DO L = 1, L_B_T
          WRITE (N_PRT, 40) ISEG, ELEM_NODES (L),          &
                            (FLUX (L, IS), IS = 1, N_G_FLUX)
          40 FORMAT ( I7, I6, 4X, (6(1PE11.3)) )
        END DO
      END IF ! ECHO
    END IF ! CONSTANT FLUX

!-->   TRANSFER INTO SYSTEM STORAGE

    L_TYPE (LAST + J_SEG)               = J_TYPE
    NODES  (LAST + J_SEG, 1:L_B_T)      = ELEM_NODES (1:L_B_T)
!bb S_FLUX (1:L_B_T, 1:N_G_FLUX, J_SEG) = FLUX (1:L_B_T, 1:N_G_FLUX)
  END DO
END SUBROUTINE INPUT_FLUX_SEG
                             
SUBROUTINE INPUT_EXACT_FLUX_SEG  (L_TYPE, NODES, X) 
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! READ FLUX BOUNDARY DATA, REPLACE WITH EXACT & STORE THEM
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Select_Source
Use Sys_Properties_Data   ! for S_FLUX (L_B_N, N_G_FLUX, N_F_SEG)
 IMPLICIT NONE
 INTEGER,  INTENT(INOUT) :: L_TYPE (L_S_TOT)
 INTEGER,  INTENT(INOUT) :: NODES  (L_S_TOT, NOD_PER_EL)
 REAL(DP), INTENT(IN) :: X (MAX_NP, N_SPACE)
!       Automatic Arrays
 REAL(DP) :: FLUX (L_B_N, N_G_FLUX)    !b (N_R_B) ?
 REAL(DP) :: FLUX_E (N_G_FLUX)    !b (N_R_B) ?
 REAL(DP) :: XYZ       (N_SPACE)    ! POINT IN SPACE
 INTEGER  :: ELEM_NODES (L_B_N)
 INTEGER  :: J_SEG, J_TYPE, L, ISEG, K, IS, LAST, L_B_T, OK

! FLUX     = COMPONENTS OF GIVEN BOUNDARY FLUX ON THE SEGMENT
! LAST     = LAST SQ MATRIX ELEMENT, N_ELEMS + N_MIXED
! L_S_TOT  = TOTAL NUMBER OF ELEMENTS & THEIR SEGMENTS
! L_TYPE   = ELEMENT TYPE NUMBER
! L_B_N    = MAXIMUM NODES ON AN ELEMENT BOUNDARY SEGMENT
! L_B_T    = NO. OF NODES ON AN ELEMENT BOUNDARY SEGMENT TYPE
! MAX_NP   = NO. OF SYSTEM NODES,
! N_D_FLUX = L_B_N * N_G_FLUX NUMBER OF SEGMENT FLUXES
! N_ELEMS  = NUMBER OF ELEMENTS
! N_F_SEG  = NO. OF ELEMENT BOUNDARY SEGMENTS IN SYSTEM
! N_G_FLUX = NO. OF FLUX COMPONENTS PER NODE
! N_G_DOF  = NO. OF DOF PER NODE
! N_MIXED  = NUMBER OF SEGMENTS WITH MIXED BOUNDARY CONDITIONS
! N_SPACE  = DIMENSION OF SPACE
! S_FLUX   = SPECIFIED COMPONENTS OF FLUX AT NODES ON SEGMENTS

!=>> DATA TYPE 19, Specified boundary flux.
!  Utilized only when N_F_SEG > 0 and L_B_N > 0 
!  TYPE 19A  
!    If NL_TYPE = 1, (single element type, usually wrong here):
!      J             = Segment number, 1 <= J <= N_F_SEG
!      ELEM_NODES(J) = List of nodes on the segment. J=1,L_B_N.
!    
!     Else if NL_TYPE > 1, (allows a mixture of element types):
!       J             = Segment number,      1 <= J      <= N_F_SEG.
!       L_TYPE        = Segment type number, 1 <= L_TYPE <= NL_TYPE.
!       ELEM_NODES(J) = List of nodes on the segment. J=1,L_B_N.
!  TYPE 19B 
!    (FLUX(K,L),L=1,N_G_FLUX = The N_G_FLUX components of flux at each
!                              of the above nodes. K=1,L_B_N.
!
!  NOTES *** These data are related to problem dependent calculations.
!  Thus the user assigns meanings to the nodal order and component
!  order.  There will be N_F_SEG TYPE 19 DATA SETs.

  IF ( N_L_TYPE == 1 ) THEN 
   WRITE (N_PRT, "('*******', /,                                     &
   &'WARNING: SEGMENTS ALMOST ALWAYS REQUIRE MULTIPLE ELEMENT TYPES',&
   &/, '*******')")
    N_WARN = N_WARN + 1 ! INCREMENT WARNING
  END IF
  IF (ECHO_FLUX) WRITE (N_PRT, &
    & "(/, '*** ELEMENT BOUNDARY FLUXES ***',/,                &
    &'SEGMENT, TYPE, ',I2,' MAX NODES ON THE SEGMENT', /,      &
    &'SEGMENT, NODE, ',I2,' EXACT FLUX COMPONENTS PER NODE')") &
    L_B_N, N_G_FLUX

  LAST = N_ELEMS + N_MIXED ; J_TYPE = 1      ! INITIALIZE
  DO ISEG = 1, N_F_SEG

!-->     READ BOUNDARY NODES
!b print *,'ISEG ', ISEG, N_L_TYPE !b
    IF ( N_L_TYPE == 1  ) THEN
      READ (N_CRD, *, IOSTAT=OK) J_SEG, ELEM_NODES (1:L_B_N) !b
      CALL INCREMENT_INPUT_LINE ! Update position count
      CALL READ_STATUS (OK, J_SEG, 0)
    ELSE
      !b READ (N_CRD, *, IOSTAT=OK) J_SEG, J_TYPE, ELEM_NODES (1:L_B_N)
      READ (N_CRD, *, IOSTAT=OK) J_SEG, J_TYPE !b
      CALL READ_STATUS (OK, J_SEG, J_TYPE) !b
      L_B_T = TYPE_NODES (J_TYPE) !b
      READ (N_CRD, *, IOSTAT=OK) J_SEG, J_TYPE, ELEM_NODES (1:L_B_T) !b
      CALL INCREMENT_INPUT_LINE ! Update position count
      CALL READ_STATUS (OK, J_SEG, J_TYPE)
    END IF
    IF ( J_SEG < 1 .OR. J_SEG > N_F_SEG ) THEN         
      PRINT *,'INPUT_EXACT_FLUX_SEG: INVALID SEGMENT NUMBER ', J_SEG
      CALL CITE_INPUT_LINE ! List current input position
      STOP 'INVALID SEGMENT NUMBER IN INPUT_EXACT_FLUX_SEG'
    END IF
    IF ( J_TYPE < 1 .OR. J_TYPE > N_L_TYPE ) THEN         
      PRINT *,'INPUT_EXACT_FLUX_SEG: INVALID ELEMENT TYPE ', J_TYPE
      CALL CITE_INPUT_LINE ! List current input position
      STOP 'INVALID ELEMENT TYPE NUMBER IN INPUT_EXACT_FLUX_SEG'
    END IF
    IF ( MINVAL(ELEM_NODES(1:L_B_T)) < 1 .OR.       &
         MAXVAL(ELEM_NODES(1:L_B_T)) > MAX_NP ) THEN         
      PRINT *,'INPUT_EXACT_FLUX_SEG: INVALID NODES AT SEGMENT ', J_SEG
      CALL CITE_INPUT_LINE ! List current input position
      STOP 'INVALID NODE NUMBER IN INPUT_EXACT_FLUX_SEG'
    END IF
    IF (ECHO_FLUX) WRITE (N_PRT, 30)      &
        J_SEG, J_TYPE, ELEM_NODES (1:L_B_T)
    30 FORMAT ( I7, I6, 3X, (16I5) )

!-->     READ BOUNDARY FLUX

    READ (N_CRD, *, IOSTAT=OK) ((FLUX (K, L), L = 1, N_G_FLUX), &
                               K = 1, L_B_T)
    CALL INCREMENT_INPUT_LINE ! Update position count
    IF ( OK /= 0 ) WRITE (N_PRT, *) 'ERROR READING BOUNDARY FLUX'
!      REPLACE INPUT VALUES WITH EXACT VALUES
    DO L = 1, L_B_T
      CALL PT_COORD (ELEM_NODES (L),  N_SPACE, X, XYZ)
      IF ( N_G_FLUX > 1 ) THEN ! USE ALL COMPONENTS
        !b CALL EXACT_FLUX (XYZ, FLUX_E)
        CALL SELECT_EXACT_FLUX (XYZ, FLUX_E) !b
      ELSE ! USE NORMAL FLUX ONLY
        !b CALL EXACT_NORMAL_FLUX (XYZ, FLUX_E)
        CALL SELECT_EXACT_NORMAL_FLUX (XYZ, FLUX_E (1)) !b
      END IF
      FLUX (L, 1:N_G_FLUX) = FLUX_E
      IF (ECHO_FLUX) WRITE (N_PRT, 40) ISEG, ELEM_NODES (L),          &
                        (FLUX (L, IS), IS = 1, N_G_FLUX)
      40 FORMAT ( I7, I6, 4X, (6(1PE11.3)) )
    END DO

!-->   TRANSFER INTO SYSTEM STORAGE

    L_TYPE (LAST + J_SEG)               = J_TYPE
    NODES  (LAST + J_SEG, 1:L_B_T)      = ELEM_NODES (1:L_B_T)
!bb S_FLUX (1:L_B_T, 1:N_G_FLUX, J_SEG) = FLUX (1:L_B_T, 1:N_G_FLUX)
  END DO
END SUBROUTINE INPUT_EXACT_FLUX_SEG
                             
SUBROUTINE INPUT_MIXED_BC_SEG  (L_TYPE, NODES) 
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!  READ & STORE MIXED BOUNDARY CONDITION SEGMENT NODES 
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants      ! for TYPE_NODES (N_L_TYPE)
Use Sys_Properties_Data   ! for S_FLUX (M_B_N, N_G_FLUX, N_MIXED)
 IMPLICIT NONE
 INTEGER,  INTENT(INOUT) :: L_TYPE (L_S_TOT)
 INTEGER,  INTENT(INOUT) :: NODES  (L_S_TOT, NOD_PER_EL)
!       Automatic Arrays
 INTEGER  :: ELEM_NODES (M_B_N)
 INTEGER  :: J_SEG, J_TYPE, ISEG, LAST, M_B_T, OK

! LAST     = LAST SQ MATRIX ELEMENT, N_ELEMS 
! L_HOMO   = 1, IF ALL ELEMENTS & SEGMENTS SAME AS FIRST
! L_S_TOT  = TOTAL NUMBER OF ELEMENTS & THEIR SEGMENTS
! L_TYPE   = ELEMENT TYPE NUMBER
! M_B_N    = MAXIMUM NUMBER OF NODES ON A MIXED BC SEGMENT     
! M_B_T    = NUMBER OF NODES ON A MIXED BC SEGMENT TYPE
! MAX_NP   = NO. OF SYSTEM NODES,
! N_ELEMS  = NUMBER OF ELEMENTS
! N_MIXED  = NUMBER OF SEGMENTS WITH MIXED BOUNDARY CONDITIONS

!=>> DATA TYPE XX, Mixed boundary condition segment nodes.
!  Utilized only when N_MIXED > 0 and M_B_N > 0 on DATA TYPE 2.
!  TYPE XXA  
!    ELEM_NODES(K) = List of nodes on the segment. K=1,M_B_N.
! NOTE *** Reading terminated when J=N_MIXED, unless L_HOMO = 1.
!   Then only J=1 is read.

  IF ( N_L_TYPE == 1 ) THEN
    WRITE (N_PRT, "('*******', /,                 &
   &'WARNING: SEGMENTS ALMOST ALWAYS REQUIRE MULTIPLE ELEMENT TYPES',&
   &/, '*******')")
    WRITE (N_PRT, *) 'PUT KEYWORD el_types 2 IN CONTROL SEGMENT'
    N_WARN = N_WARN + 1 ! INCREMENT WARNING
  END IF
  IF (ECHO_MIXED) WRITE                                       &
    (N_PRT, "(/, '*** MIXED BOUNDARY CONDITION SEGMENTS ***')") 
  IF ( M_B_N <= 0 ) THEN
    CALL CITE_INPUT_LINE ! List current input position
    PRINT *, 'ERROR: INPUT_MIXED_BC_SEG, NO NODES ON A MIXED BC SEGMENT'
    PRINT *, 'PUT KEYWORD el_segment AND VALUE IN CONTROL'
    STOP     'ERROR: INPUT_MIXED_BC_SEG, NO NODES ON A MIXED BC SEGMENT'
  END IF

  IF ( CONVECTION ) THEN ! CONSTANT CONVECTION ON ALL MIXED SEGMENTS
    PRINT *,'NOTE: CONSTANT CONVECTION ON ALL MIXED SEGMENTS USING'
    PRINT *,'CONVECTION COEFFICIENT = ', CONVECT_COEF
    PRINT *,'CONVECTION TEMPERATURE = ', CONVECT_TEMP 
    IF ( CONVECT_THICK /= 1.d0 )                               &
         PRINT *,'LINE CONVECTION THICKNESS   = ', CONVECT_THICK
  END IF ! CONSTANT CONVECTION

  IF ( ROBIN_1_SEG /= 0.d0 ) THEN ! CONSTANT ROBIN ON ALL MIXED SEGMENTS
    PRINT *,'NOTE: CONSTANT ROBIN BC ON ALL MIXED SEGMENTS USING'
    PRINT *,'K_n * U,n + ROBIN_1 * U + ROBIN_2 = 0, WITH'
    PRINT *,'ROBIN 1, 2 TERMS = ', ROBIN_1_SEG, ROBIN_2_SEG
    IF ( ROBIN_THICK /= 1.d0 )                        &
         PRINT *,'ROBIN BC THICKNESS   = ', ROBIN_THICK
  END IF ! CONSTANT ROBIN

  IF (ECHO_MIXED) WRITE                                           &
    (N_PRT, "('SEGMENT, TYPE, ',I2,' MAX NODES ON THE SEGMENT')") M_B_N

  LAST = N_ELEMS ; J_TYPE = 1 ; M_B_T = M_B_N      ! INITIALIZE
  ELEM_NODES = 0
  DO ISEG = 1, N_MIXED

!-->     READ BOUNDARY NODES
    IF ( N_L_TYPE == 1  ) THEN
      READ (N_CRD, *, IOSTAT=OK) J_SEG, ELEM_NODES (1:M_B_N) !b
      CALL INCREMENT_INPUT_LINE ! Update position count
      CALL READ_STATUS (OK, J_SEG, 0)
    ELSE
      !b READ (N_CRD, *, IOSTAT=OK) J_SEG, J_TYPE, ELEM_NODES (1:M_B_N)
      READ (N_CRD, *, IOSTAT=OK) J_SEG, J_TYPE !b
      CALL READ_STATUS (OK, J_SEG, J_TYPE) !b
      M_B_T = TYPE_NODES (J_TYPE)        !b
!b print *,'J_SEG, J_TYPE, M_B_T ',J_SEG, J_TYPE, M_B_T
      BACKSPACE (N_CRD) !b
      READ (N_CRD, *, IOSTAT=OK) J_SEG, J_TYPE, ELEM_NODES (1:M_B_T) !b
!b print *, 'ELEM_NODES (1:M_B_T) ', ELEM_NODES (1:M_B_T)
      CALL INCREMENT_INPUT_LINE ! Update position count
      CALL READ_STATUS (OK, J_SEG, J_TYPE)
    END IF
    IF ( ECHO_MIXED ) WRITE (N_PRT, 30)    &
         J_SEG, J_TYPE, ELEM_NODES (1:M_B_T)
    30 FORMAT ( I7, I6, 3X, (16I5) )
    IF( J_SEG < 1 .OR. J_SEG > N_MIXED) THEN
      CALL CITE_INPUT_LINE ! List current input position
      STOP 'INVALID SEGMENT NUMBER IN INPUT_MIXED_BC_SEG'
    END IF
    IF( J_TYPE < 1 .OR. J_TYPE > N_L_TYPE) THEN
      CALL CITE_INPUT_LINE ! List current input position
      STOP 'INVALID ELEMENT TYPE NUMBER IN INPUT_MIXED_BC_SEG'
    END IF
    IF( MINVAL(ELEM_NODES(1:M_B_T)) < 1 .OR.       &
        MAXVAL(ELEM_NODES(1:M_B_T)) > MAX_NP ) THEN 
        CALL CITE_INPUT_LINE ! List current input position
        PRINT *, 'J_SEG, J_TYPE, ELEM_NODES (1:M_B_T)',   &
                  J_SEG, J_TYPE, ELEM_NODES (1:M_B_T)
        STOP 'INVALID NODE NUMBER IN INPUT_MIXED_BC_SEG'
    END IF

!-->   TRANSFER INTO SYSTEM STORAGE
    L_TYPE (LAST + J_SEG)               = J_TYPE
    NODES  (LAST + J_SEG, 1:M_B_T)      = ELEM_NODES (1:M_B_T)
  END DO ! OVER MIXED BC
END SUBROUTINE INPUT_MIXED_BC_SEG
                             
SUBROUTINE LIST_CONTROLS
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!             LIST THE PROBLEM CONTROL DATA
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
IMPLICIT NONE
INTEGER, SAVE :: COUNT, EX_BC = 0, EX_F = 0, EX_S = 0, AXI_S = 0

!         PROBLEM CLASS
!b  WRITE (U_HIST, 21) I_SAY, N_SPACE, N_R_B, N_ITER,              &
!b                     MAX_NP, N_ELEMS, N_G_DOF,                   &
!b                     NOD_PER_EL, N_MIXED, M_B_N, N_F_SEG, L_B_N, &
!b                     N_G_FLUX, N_L_TYPE, L_SHAPE, N_QP
!b  WRITE (U_LOG , 21) I_SAY, N_SPACE, N_R_B, N_ITER,              &
!b                     MAX_NP, N_ELEMS, N_G_DOF,                   &
!b                     NOD_PER_EL, N_MIXED, M_B_N, N_F_SEG, L_B_N, &
!b                     N_G_FLUX, N_L_TYPE, L_SHAPE, N_QP
  WRITE (6, 21) I_SAY, N_SPACE, N_R_B, N_ITER,              &
                MAX_NP, N_ELEMS, N_G_DOF,                   &
                NOD_PER_EL, N_MIXED, M_B_N, N_F_SEG, L_B_N, &
                N_G_FLUX, N_L_TYPE, L_SHAPE, N_QP
 21 FORMAT (/,'****  PROBLEM CLASS: (DEFAULT) VALUE ****',/,  &
   'NUMBER OF USER REMARKS LINES ..................(0)',I5,/, &
   'DIMENSION OF SPACE ............................(1)',I5,/, &
   'NUMBER OF ROWS IN B MATRIX ....................(1)',I5,/, &
   'NUMBER OF ITERATIONS TO BE RUN ................(1)',I5,/, &
   'NUMBER OF NODAL POINTS IN SYSTEM ..............(2)',I5,/, &
   'NUMBER OF ELEMENTS IN SYSTEM ..................(1)',I5,/, &
   'NUMBER OF PARAMETERS PER NODE .................(1)',I5,/, &
   'MAXIMUM NUMBER OF NODES PER ELEMENT ...........(2)',I5,/, &
   'NUMBER OF MIXED BOUNDARY CONDITION SEGMENTS ...(0)',I5,/, &
   'NUMBER OF NODES PER MIXED BC SEGMENT ..........(0)',I5,/, &
   'NUMBER OF BOUNDARIES WITH GIVEN FLUX ..........(0)',I5,/, &
   'NUMBER OF NODES ON FLUX BOUNDARY SEGMENT ......(0)',I5,/, &
   'NUMBER OF FLUX COMPONENTS PER FLUX SEG NODE ...(0)',I5,/, &
   'NUMBER OF DIFFERENT ELEMENT TYPES .............(1)',I5,/, &
   'SHAPE 1=LINE 2=TRI 3=QUAD 4=HEX 5=TET 6=WEDG ..(1)',I5,/, &
   'MAXIMUM NUMBER OF QUADRATURE POINTS ...........(0)',I5) 

!       PROBLEM PROPERTIES
  WRITE (6, 22) N_LP_FIX, N_LP_FLO, N_BS_FIX, N_BS_FLO, &
                MISC_FX, MISC_FL, N_MX_FIX, N_MX_FLO, &
                N_NP_FIX, N_NP_FLO, N_MAT, MAT_FLO
 22 FORMAT (/,'****  PROBLEM PROPERTIES: (DEFAULT) VALUE ****',/,&
   'NUMBER OF INTEGER PROPERTIES PER ELEMENT ......(0)',I5,/, &
   'NUMBER OF REAL    PROPERTIES PER ELEMENT ......(0)',I5,/, &
   'NUMBER OF INTEGER PROPERTIES PER FLUX SEGMENT..(0)',I5,/, &
   'NUMBER OF REAL    PROPERTIES PER FLUX SEGMENT..(0)',I5,/, &
   'NUMBER OF INTEGER MISCELLANEOUS PROPERTIES ....(0)',I5,/, &
   'NUMBER OF REAL    MISCELLANEOUS PROPERTIES ....(0)',I5,/, &
   'NUMBER OF INTEGER PROPERTIES PER MIXED BC .....(0)',I5,/, &
   'NUMBER OF REAL    PROPERTIES PER MIXED BC .....(0)',I5,/, &
   'NUMBER OF INTEGER PROPERTIES PER NODE .........(0)',I5,/, &
   'NUMBER OF REAL    PROPERTIES PER NODE .........(0)',I5,/, &
   'NUMBER OF MATERIAL PROPERTY SETS ..............(0)',I5,/, &
   'NUMBER OF REAL    PROPERTIES PER MATERIAL SET..(0)',I5)

!        PROBLEM OPTIONS
  IF ( USE_EXACT_BC   ) EX_BC = 1
  IF ( USE_EXACT_FLUX ) EX_F  = 1
  IF ( USE_EXACT      ) EX_S  = 1
  IF ( AXISYMMETRIC   ) AXI_S = 1 !b
  WRITE (6, 23) AXI_S, IN_RHS, NUL_COL, MODE, L_HOMO, N_HOMO,  &
                LEM_WRT, NPT_WRT, LIST_LR, EX_BC, EX_F, &
                EX_S, N_FILE1, N_FILE2
 23 FORMAT (/,'****  OPTIONS: (DEFAULT) VALUE ****',/, &
   'AXISYMMETRIC DOMAIN:         0=FALSE, 1=TRUE ..(0)',I5,/, &
   'INITIAL FORCING VECTOR:      0=OMIT,  1=READ ..(0)',I5,/, &
   'NULL ELEM SOURCE VECTORS:    0=FALSE, 1=TRUE ..(0)',I5,/, &
   'UNSYMMETRIC MATRIX MODE:     0=FALSE, 1=TRUE ..(0)',I5,/, &
   'HOMOGENEOUS ELEM PROPERTIES: 0=FALSE, 1=TRUE ..(0)',I5,/, &
   'HOMOGENEOUS NODE PROPERTIES: 0=FALSE, 1=TRUE ..(0)',I5,/, &
   'LIST RESULTS BY ELEMS:       0=FALSE, 1=TRUE ..(0)',I5,/, &
   'LIST RESULTS BY NODES:       0=FALSE, 1=TRUE ..(0)',I5,/, &
   'LIST ELEMENT REACTIONS:      0=OMIT,  1=LIST ..(0)',I5,/, &
   'USE EXACT BC IN FE SOLUTION: 0=FALSE, 1=TRUE ..(0)',I5,/, &
   'USE & LIST EXACT FLUXES:     0=FALSE, 1=TRUE ..(0)',I5,/, &
   'USE & LIST EXACT SOLUTION:   0=FALSE, 1=TRUE ..(0)',I5,/, &
   'USER POST_PROCESSING UNIT FLAG_1, ACTIVE > 0 ..(0)',I5,/, &
   'USER POST_PROCESSING UNIT FLAG_2, ACTIVE > 0 ..(0)',I5)

!        PROBLEM SUMMARY
  WRITE (6, 24) N_EL_FRE, N_D_FLUX, (M_B_N*N_G_DOF), N_D_FRE
 24 FORMAT (/,'****  SUMMARY:  (DEFAULT) VALUE ****',/, &
  'NUMBER OF D.O.F. FOR ELEMENT ..................(2)',I5,/,  &
  'NUMBER OF D.O.F. ON MIXED BC SEGMENT ..........(0)',I5,/,  &
  'NUMBER OF D.O.F. ON FLUX SEGMENT ..............(0)',I5,/,  &
  'NUMBER OF D.O.F. IN TOTAL SYSTEM ..............(2)',I5)


 !bWRITE (6, 30) MAX_NP, N_ELEMS, N_G_DOF, NOD_PER_EL, N_SPACE, &
                 !bN_F_SEG, L_B_N, N_ITER, LIST_LR, IN_RHS, I_SAY,&
                 !bN_R_B, N_QP, L_SHAPE, N_L_TYPE, MODE,          &
                  !bN_MIXED, M_B_N, N_MX_FIX, N_MX_FLO
                  !b!b N_R_B, N_QP, L_SHAPE, N_L_TYPE, MODE
 !30 FORMAT (/,'****  PROBLEM PARAMETERS (DEFAULT IF 0) ****',/, &
 ! 'NUMBER OF NODAL POINTS IN SYSTEM ..............(2)',I5,/, &
 ! 'NUMBER OF ELEMENTS IN SYSTEM ..................(1)',I5,/, &
 ! 'NUMBER OF PARAMETERS PER NODE .................(1)',I5,/, &
 ! 'MAXIMUM NUMBER OF NODES PER ELEMENT ...........(2)',I5,/, &
 ! 'DIMENSION OF SPACE ............................(1)',I5,/, &
 ! 'NUMBER OF BOUNDARIES WITH GIVEN FLUX ..........(0)',I5,/, &
 ! 'NUMBER OF NODES ON FLUX BOUNDARY SEGMENT ......(0)',I5,/, &
 ! 'NUMBER OF ITERATIONS TO BE RUN ................(1)',I5,/, &
 ! 'LIST ELEMENT REACTIONS: 0=OMIT, 1=LIST ........(0)',I5,/, &
 ! 'INITIAL FORCING VECTOR: 0=OMIT, 1=READ ........(0)',I5,/, &
 ! 'NUMBER OF USER REMARKS LINES ..................(0)',I5,/, &
 ! 'NUMBER OF ROWS IN B MATRIX ....................(1)',I5,/, &
 ! 'MAXIMUM NUMBER OF QUADRATURE POINTS ...........(0)',I5,/, &
 ! 'SHAPE 1=LINE 2=TRI 3=QUAD 4=HEX 5=TET 6=WEDG ..(1)',I5,/, &
 ! 'NUMBER OF DIFFERENT ELEMENT TYPES .............(1)',I5,/, &
 ! 'UNSYMMETRIC MATRIX MODE: 0=FALSE, 1=TRUE ......(0)',I5,/, &
 ! 'NUMBER OF MIXED BOUNDARY CONDITION SEGMENTS ...(0)',I5,/, &
 ! 'NUMBER OF NODES PER MIXED BC SEGMENT ..........(0)',I5,/, &
 ! 'NUMBER OF INTEGER PROPERTIES PER MIXED BC .....(0)',I5,/, &
 ! 'NUMBER OF REAL    PROPERTIES PER MIXED BC .....(0)',I5)
   !b 'UNSYMMETRIC MATRIX MODE: 0=FALSE, 1=TRUE ......(0)',I5)

  !bWRITE (6, 80) N_NP_FIX, N_NP_FLO, N_LP_FIX, N_LP_FLO,     &
                  !bMISC_FX, MISC_FL, N_HOMO, L_HOMO, NPT_WRT,  &
                  !bLEM_WRT, N_FILE1, N_FILE2, N_G_FLUX,        &
                  !bNUL_COL, N_BS_FIX, N_BS_FLO
 !80 FORMAT ( /,                                            &
 ! 'NUMBER OF INTEGER PROPERTIES PER NODE .........(0)',I5,/, &
 ! 'NUMBER OF REAL    PROPERTIES PER NODE .........(0)',I5,/, &
 ! 'NUMBER OF INTEGER PROPERTIES PER ELEMENT ......(0)',I5,/, &
 ! 'NUMBER OF REAL    PROPERTIES PER ELEMENT ......(0)',I5,/, &
 ! 'NUMBER OF INTEGER MISCELLANEOUS PROPERTIES ....(0)',I5,/, &
 ! 'NUMBER OF REAL    MISCELLANEOUS PROPERTIES ....(0)',I5,/, &
 ! 'HOMOGENEOUS NODE PROPERTIES: 0=FALSE, 1=TRUE ..(0)',I5,/, &
 ! 'HOMOGENEOUS ELEM PROPERTIES: 0=FALSE, 1=TRUE ..(0)',I5,/, &
 ! 'LIST RESULTS BY NODES: 0=FALSE, 1=TRUE ........(0)',I5,/, &
 ! 'LIST RESULTS BY ELEMS: 0=FALSE, 1=TRUE ........(0)',I5,/, &
 ! 'POST_PROCESSING UNIT FLAG_1, ACTIVE > 0 .......(0)',I5,/, &
 ! 'POST_PROCESSING UNIT FLAG_2, ACTIVE > 0 .......(0)',I5,/, &
 ! 'NUMBER OF FLUX COMPONENTS PER FLUX SEG NODE ...(0)',I5,/, &
 ! 'NULL ELEM SOURCE VECTORS: 0=FALSE, 1=TRUE .....(0)',I5,/, & 
 ! 'NUMBER OF INTEGER PROPERTIES PER FLUX SEGMENT..(0)',I5,/, &
 ! 'NUMBER OF REAL    PROPERTIES PER FLUX SEGMENT..(0)',I5)

  !bWRITE (6, 81) N_EL_FRE, N_D_FLUX, (M_B_N*N_G_DOF), N_D_FRE
 !81 FORMAT ( /,                                            &
 !'NUMBER OF D.O.F. FOR ELEMENT ..................(2)',I5,/,  &
 !'NUMBER OF D.O.F. ON MIXED BC SEGMENT ..........(0)',I5,/,  &
 !'NUMBER OF D.O.F. ON FLUX SEGMENT ..............(0)',I5,/,  &
 !'NUMBER OF D.O.F. IN TOTAL SYSTEM ..............(2)',I5)
!
  IF ( I_BUG > 0 )  WRITE (6, * ) 'DEBUG PRINTS TURNED ON'

  COUNT = N_FILE1 + N_FILE2 
  IF ( COUNT > 0 ) THEN
    PRINT *, ' '
    PRINT *, 'APPLICATION DEPENDENT POST-PROCESSING IS ACTIVE'
  END IF
  IF ( N_QP > 0 .AND. U_SCPR > 0 ) PRINT *,  &
    'SUPER_CONVERGENT PATCH IS ACTIVE'
  COUNT = N_FILE3 + N_FILE4 + N_FILE5
  IF ( COUNT > 0 ) WRITE (6, 88) N_FILE3, N_FILE4, N_FILE5
    88 FORMAT (/, 'OPTIONAL UNIT NUMBERS (UTILIZED IF > 0)',/,&
    'N_FILE3 = ',I2,', N_FILE4 = ',I2,/,'N_FILE5 = ',I2)

END SUBROUTINE LIST_CONTROLS

SUBROUTINE FIXED_CONTROL_INPUT
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!  Old style fixed input (replaces OLD_CONTROL_INPUT) for
!  old dusty data files. See newer PARSE_INPUT
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data
 IMPLICIT NONE

!  ==>> DATA TYPE 1, Title
!    TITLE, A general description of the problem.  Enclose in
!           double quotes or use_string_with_no_blanks.

!  ==>> DATA TYPE 2, Control parameters.

!    MAX_NP     = Number of nodal points in the system.
!    N_ELEMS    = Number of elements in the system.
!    N_G_DOF    = Number of parameters per node.
!    NOD_PER_EL = Number of nodes per element.
!
!    N_SPACE  = Dimension of solution space.
!    N_F_SEG  = No. of elem boundary segments with given flux.
!    L_B_N    = Number of nodes on an element boundary segment.
!    N_ITER   = Number of iterations to be run. (If any.)
!
!    LIST_LR  = List reactions on elements, 0=no, 1=yes
!    IN_RHS   = Initial system column matrix flag.
!             = 0, Not to be input. (Default)
!             = 1, Will be input by DATA TYPE 18.
!    I_SAY    = Number of user remark cards to be read and
!               printed. See DATA TYPE 4.
!    N_R_B    = Number of rows in element B matrix. (Default=1)
!
!    N_QP     = Number of quadrature points. (Default=0)
!    L_SHAPE  = Shape 1-line 2-tri 3-quad 4-hex 5-tet 6-wedge
!    N_L_TYPE = Number of element types. (Default=1)
!    MODE     = Unsymmetric storage flag 0=false, 1=true (Default=0)
!
!    NOTE: A negative value of N_ITER turns on debug printing.
!    The absolute value of N_ITER is used.

!  ==>> DATA TYPE 3, Control parameters.

!    N_NP_FIX = Number of integer properties per node.
!    N_NP_FLO = Number of real properties per node.
!    N_LP_FIX = Number of integer properties per elem.
!    N_LP_FLO = Number of real properties per elem.
!
!    MISC_FX  = Number of miscellaneous integer properties.
!    MISC_FL  = Number of miscellaneous real properties
!    N_HOMO   = Nodal properties flag.
!              = 1, Properties are homogeneous. Only values for 
!                the first node will be input.
!              = 0, Non-homogeneous properties. (Default)
!                The property array is initially filled with copies 
!                of the properties of the first entry. Then only 
!                different entries and the last node set are input.
!    L_HOMO   = Element properties flag.
!              = 1, Properties are homogeneous. Only values for 
!                the first element will be input.
!              = 0, Non-homogeneous properties. (Default)
!                The property array is initially filled with copies 
!                of the properties of the first entry. Then only 
!                different entries and the last element set are input.
!
!    NPT_WRT  = Nodal parameter print flag.
!              = 1, List calculated values by nodes. 
!                Otherwise omit nodal list.                (Default=0)
!    LEM_WRT  = Nodal parameter print flag.
!              = 1, List node and calculated nodal values for 
!                each element.  Otherwise omit list.        (Default=0)
!    N_FILE1  = I/O Unit flag for post solution calculation data.
!              = 0, Not used in post solution calculations. (Default=0)
!              = 1, Unit to be utilized.
!    N_FILE2  = I/O Unit flag for post solution calculation data.
!              = 0, Not used in post solution calculations. (Default=0)
!              = 1, Unit to be utilized.
!
!    N_G_FLUX = Number of generalized flux components (N_F_SEG>0)
!    NUL_COL  = Element column matrix flag.
!              = 0, Problem requires the matrix.            (Default=0)
!              = 1, The matrix is always a null vector.
!    N_BS_FIX = Number of integer properties per segment.
!    N_BS_FLO = Number of real properties per segment.
!
!    NOTE: User provided post-processing is activated if either
!          N_FILE1 or N_FILE2 is > 0.  That combination allows the
!          user at least two options within their post-processing.

!      PRINT AUTHOR CREDITS
! CALL  TOOT

  CALL SET_CONTROL_DEFAULTS  ! Initialize

!-->    ** READ AND PRINT TITLE AND CONTROL DATA **
  READ (5, * ) TITLE 
  CALL INCREMENT_INPUT_LINE ! Update position count
  READ (5, * ) MAX_NP,  N_ELEMS, N_G_DOF,  NOD_PER_EL,  &
               N_SPACE, N_F_SEG,   L_B_N,    N_ITER,    &
               LIST_LR, IN_RHS,  I_SAY,    N_R_B,       &
               N_QP,    L_SHAPE, N_L_TYPE, MODE
  CALL INCREMENT_INPUT_LINE ! Update position count
  READ (5, * ) N_NP_FIX, N_NP_FLO, N_LP_FIX, N_LP_FLO,  &
               MISC_FX,  MISC_FL,  N_HOMO,   L_HOMO,    &
               NPT_WRT,  LEM_WRT,  N_FILE1,  N_FILE2,   &
               N_G_FLUX, NUL_COL,  N_BS_FIX, N_BS_FLO
  CALL INCREMENT_INPUT_LINE ! Update position count

! TRANSLATE TO KEYWORD FORM FOR FUTURE USE
  CALL OLD_F77_CONTROL_TO_KEYWORDS

END SUBROUTINE FIXED_CONTROL_INPUT

SUBROUTINE GET_LT_FACES (IE, N_FACES)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!   GET THE NUMBER OF FACES ON ELEMENT NUMBER IE
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Elem_Type_Data     ! for LT_FACES
 IMPLICIT NONE
  INTEGER, INTENT(IN)  :: IE          ! ELEMENT NUMBER
  INTEGER, INTENT(OUT) :: N_FACES     ! FACE COUNT
    N_FACES = LT_FACES
if ( N_FACES < 1 ) then !b
 print *,'ERROR, GET_LT_FACES, LT_FACES UNDEFINED ELEMENT ', IE
 stop 'GET_LT_FACES: LT_FACES UNDEFINED'
end if !b
END SUBROUTINE GET_LT_FACES 

SUBROUTINE VALIDATE_AND_LIST_CONTROLS
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!      VALIDATE_AND_LIST_CONTROL DATA
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data
 IMPLICIT NONE

  L_S_TOT = N_ELEMS + N_MIXED + N_F_SEG ! Total elements, and segments

!   CHECK DEFAULTS
  !b IF ( N_ANAL      < 1 )  N_ANAL = 1  !b
  !b IF ( N_GEOM > N_ANAL )                           &  !b
    !b PRINT *, 'NOTE: VALIDATE_AND_LIST_CONTROLS, ', &  !b
    !b 'MORE GEOMETRY NODES THAN ANALYSIS NODES'  !b
  !b NOD_PER_EL = MAX (NOD_PER_EL, N_GEOM, N_ANAL)  !b
  IF ( L_SHAPE    < 1)  L_SHAPE    = 1
  IF ( NOD_PER_EL < 1)  NOD_PER_EL = 2
  IF ( N_G_DOF    < 1)  N_G_DOF    = 1
  IF ( N_L_TYPE   < 1)  N_L_TYPE   = 1
  IF ( N_R_B      < 1)  N_R_B      = 1
  IF ( N_SPACE    < 1)  N_SPACE    = 1
  
SELECT CASE ( N_SPACE )  ! Number of second derivatives
  CASE ( 1 ); N_2_DER = 1 ; N_3_DER = 1
  CASE ( 2 ); N_2_DER = 3 ; N_3_DER = 4
  CASE ( 3 ); N_2_DER = 6 ; N_3_DER = 9
  CASE DEFAULT
    PRINT *, 'WARNING: N_2_DER and N_3_DER = 0 DEFAULT SET'
    PRINT *, 'Can set N_2_DER and N_3_DER in Keywords'
    N_2_DER = 0 ; N_3_DER = 0 ; N_WARN = N_WARN + 1
END SELECT

  CALL CHECK_EL_SHAPE       !  CHECK ELEMENT SHAPE DATA

  !b IF ( LIST_LR == 0 ) L_REACT = 0 ! Element reaction data storage unit
  IF ( LIST_LR == 1 ) L_REACT = 21 ! Element reaction data storage unit !b

  IF ( N_QP < 1 .AND. U_SCPR > 0 ) THEN
    PRINT *,'NOTE: SCP Recovery turned off since N_QP = 0'
    U_SCPR = 0              ! No active SCP direct file unit
    U_FLUX = 0              ! No fluxes to save to average
U_GRAD = U_FLUX !bb
    U_PLT5 = 0              ! No node gradient output for plotters
    U_PLT6 = 0              ! No element error output for plotters
    U_PLT7 = 0              ! No pt ave error est output for plotters
    !b U_PLT8 = 0              ! No exact solution output for plots
    !b U_PLT9 = 0              ! No exact flux output for plots
  END IF ! No SCP Recovery

  SCP_FIT  = N_R_B          ! default terms for nodal averaging
  IF ( SCP_2ND_DERIV ) SCP_FIT = N_R_B * (N_SPACE + 1)
  IF ( SCP_3RD_DERIV ) SCP_FIT = N_R_B * (2 * N_SPACE + 1)
 
!  SET ADVANCED USER DEFAULTS
  IF ( N_CORNER == 0 ) N_CORNER = NOD_PER_EL  ! Number of corner nodes
  IF ( N_GEOM   == 0 ) N_GEOM   = NOD_PER_EL  ! Number of geometry nodes
  IF ( N_PARM   == 0 ) N_PARM   = N_SPACE     ! Parametric dimensions 

!  CHECK FOR EXTRA OUTPUT REQUESTS
  !b I_BUG = 0
  IF ( N_ITER < 0 ) THEN
    N_ITER = IABS (N_ITER)
    I_BUG  = N_ITER           ! allow levels of debug
  END IF
  IF ( N_ITER < 1 ) N_ITER = 1

  IF ( L_B_N > NOD_PER_EL ) THEN  
    WRITE (6, *) 'WARNING, NUMBER OF FLUX SEGMENT NODES > ELEMENT',&
                 ' NODES, RESET TO ', NOD_PER_EL
    WRITE (6, *) 'REVISE el_segment IN KEYWORDS'
    N_WARN = N_WARN + 1 ! INCREMENT WARNING
    L_B_N = NOD_PER_EL
  END IF
  N_EL_FRE = NOD_PER_EL * N_G_DOF
  !b N_EL_FRE = N_ANAL * N_G_DOF !b
  N_D_FRE  = MAX_NP * N_G_DOF
  N_D_FLUX = L_B_N  * N_G_DOF
  IF ( N_F_SEG > 0 .AND. N_G_FLUX < 1 ) THEN
    WRITE (6, *) 'WARNING, NUMBER OF FLUX COMPONENTS SET = ', N_G_DOF
    WRITE (6, *) 'INCLUDE seg_pt_flux IN KEYWORDS'
    N_WARN   = N_WARN + 1 ! INCREMENT WARNING
    N_G_FLUX = N_G_DOF
  END IF

! IP_TEST  > 0, SOME PROPERTIES ARE DEFINED
  IP_TEST = N_NP_FIX + N_NP_FLO + N_LP_FIX + N_LP_FLO &
          + N_BS_FIX + N_BS_FLO + MISC_FX  + MISC_FL  &
          + N_MX_FIX + N_MX_FLO + MAT_FLO + EXACT_FL + EXACT_FX !b
          !b + N_MX_FIX + N_MX_FLO + MAT_FLO !b

! LP_TEST  > 0, ELEMENT PROPERTIES ARE DEFINED
  LP_TEST = N_LP_FIX + N_LP_FLO

!  INITIALIZE ELEMENT TYPE DATA
  LAST_LT = 0
  LT_FREE = N_EL_FRE
  LT_GEOM = N_GEOM
  LT_N    = NOD_PER_EL
  LT_PARM = N_PARM
  LT_QP   = N_QP
  LT_SHAP = L_SHAPE
  LT_FACES = 0 ; LT_SONS = 0 ; LT_USER = 0

! EIGENPROBLEMS
  IF ( EIGEN_SHOW > N_D_FRE ) EIGEN_SHOW = N_D_FRE
  IF ( EIGEN_SCP  > N_D_FRE ) EIGEN_SCP  = N_D_FRE

  !b WRITE (U_HIST, *) 'TITLE: ', TITLE 
  !b WRITE (U_LOG , *) 'TITLE: ', TITLE 
  WRITE (6,      *) 'TITLE: ', TITLE 
  IF ( ECHO_INPUT ) CALL LIST_CONTROLS
  IF ( I_SAY > 0) CALL INPUT_REMARKS (I_SAY)
  IF ( DEBUG_UNITS ) THEN 
    CALL LIST_SYSTEM_UNIT_NUMBERS
    CALL LIST_OPEN_IO_UNITS
  END IF

END SUBROUTINE VALIDATE_AND_LIST_CONTROLS

subroutine parse_input
!---------------------------------------------------------
! Parse a general input stream of words and data
!---------------------------------------------------------
 use keyword_buffer  ! MAX_KEY, ECHO_KEY, ECHO_OPEN, ECHO_UNIT
 Use System_Constants ! for LINE_IN and INCREMENT_INPUT_LINE 
 implicit none
 character(len=MAX_KEY)   :: key    
 logical, save            :: enough = .false.   
 integer                  :: ok       ! file status

  !b ECHO_KEY = .true.             ! to echo keyword inputs 
  !b call open_echo_unit           ! file keyword_input.echo
  if ( .not. ECHO_OPEN ) then
    open (ECHO_UNIT, file = 'keyword_input.echo', &
          action ='write',  status = 'unknown',   &
          iostat = ok) 
    if ( ok == 0 ) ECHO_OPEN = .true.
  end if

!...INPUT LOOP       
  print *, 'Beginning user keyword inputs'
  !b if ( .NOT. ECHO_KEY ) print *,'NOTE: See file keyword_input.echo'
  do while (.not. enough)
    CALL INCREMENT_INPUT_LINE ! Update position count
    call get_key (key)  ! first keyword on a line

    select case ( key ) ! Check allowed keys
      Case ('abort') ; stop 'user exit with keyword abort'
      Case ('done')  ; enough = .true.
      Case ('EOF')   ; enough = .true.
      Case ('end')   ; enough = .true.
      Case ('exit')  ; enough = .true.
      Case ('quit')  ; enough = .true.
      Case Default   ! general or application control
        call apply_key ( key ) ! *** application dependent ***
        !b if ( ECHO_KEY ) write (ECHO_UNIT, *) ' ' ! must end line
        write (ECHO_UNIT, *) ' ' ! must end line !b
    end select
    if ( enough ) print *, 'Keyword input terminated by user.'
  end do

! Validate multiple element types if used
      if ( (N_F_SEG > 0 .OR. N_MIXED > 0 ) .AND. L_B_N == 0 ) then
         PRINT *,'WARNING: NODES ON BOUNDARY SEGMENT DEFAULTED TO ', &
                  NOD_PER_EL ; N_WARN = N_WARN + 1
         L_B_N = NOD_PER_EL
       if ( N_GEOM == 0 ) N_GEOM = NOD_PER_EL
      end if
end subroutine parse_input

!SUBROUTINE OLD_CONTROL_TO_KEYWORDS
!! now in keyword_lib.f
!! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!!      CONVERT OLD F90 FIXED CONTROL DATA TO KEYWORD FORM
!!               STORED IN FILE old_control.key
!! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
!Use System_Constants
!  implicit none
!  integer :: unit, last
!  !b integer :: get_next_io_unit
!
!  interface 
!    function get_next_io_unit () result (next)
!      integer :: next
!    end function get_next_io_unit
!  end interface
!
!!    get a free unit and open the file to save (over write)
!  unit = get_next_io_unit ()
!  !b print *, 'unit = ', unit !b
!  open (unit, file = "old_control.key", status = "unknown")
!
!!    save the old controls and title in keyword style
!   last = len_trim (title)
!   write (unit, *) 'title "', title (1:last), '" ! begin keywords'
!
!!   first old control line
!   write (unit, *) 'nodes      ', MAX_NP
!   write (unit, *) 'elems      ', N_ELEMS
!   write (unit, *) 'dof        ', N_G_DOF
!   write (unit, *) 'el_nodes   ', NOD_PER_EL
!   write (unit, *) 'space      ', N_SPACE
!   write (unit, *) 'b_rows     ', N_R_B
!   write (unit, *) 'shape      ', L_SHAPE
!   if ( N_F_SEG  > 0 ) write (unit, *) 'segments   ', N_F_SEG
!   if ( L_B_N    > 0 ) write (unit, *) 'el_segment ', L_B_N
!   if ( N_ITER   > 1 ) write (unit, *) 'iters      ', N_ITER
!   if ( N_ITER   < 0 ) then
!                       write (unit, *) 'iters      ', IABS(N_ITER)
!                       write (unit, *) 'debug      '
!   end if
!   if ( LIST_LR  > 0 ) write (unit, *) 'el_react   '! LIST_LR
!   if ( IN_RHS   > 0 ) write (unit, *) 'loads      '! IN_RHS
!   if ( I_SAY    > 0 ) write (unit, *) 'remarks    ', I_SAY
!   if ( N_QP     > 0 ) write (unit, *) 'gauss      ', N_QP
!   if ( N_L_TYPE > 0 ) write (unit, *) 'el_types   ', N_L_TYPE
!   if ( MODE == 0 ) then
!     write (unit, *) 'symmetric '
!   else
!     write (unit, *) 'unsymmetric'
!   end if
!
!!   second old control line
!   if ( N_NP_FIX > 0 ) write (unit, *) 'pt_int     ', N_NP_FIX
!   if ( N_NP_FLO > 0 ) write (unit, *) 'pt_real    ', N_NP_FLO
!   if ( N_LP_FIX > 0 ) write (unit, *) 'el_int     ', N_LP_FIX
!   if ( N_LP_FLO > 0 ) write (unit, *) 'el_real    ', N_LP_FLO
!   if ( MISC_FX  > 0 ) write (unit, *) 'integers   ', MISC_FX
!   if ( MISC_FL  > 0 ) write (unit, *) 'reals      ', MISC_FL
!   if ( N_HOMO   > 0 ) write (unit, *) 'pt_homo    '! N_HOMO
!   if ( L_HOMO   > 0 ) write (unit, *) 'el_homo    '! L_HOMO
!   if ( NPT_WRT  > 0 ) write (unit, *) 'pt_list    '! NPT_WRT
!   if ( LEM_WRT  > 0 ) write (unit, *) 'el_list    '! LEM_WRT
!   if ( N_FILE1  > 0 ) write (unit, *) 'post_1     '! N_FILE1
!   if ( N_FILE2  > 0 ) write (unit, *) 'post_2     '! N_FILE2  
!   if ( N_G_FLUX > 0 ) write (unit, *) 'seg_pt_flux', N_G_FLUX
!   if ( NUL_COL  > 0 ) write (unit, *) 'el_no_col  '! NUL_COL
!   if ( N_BS_FIX > 0 ) write (unit, *) 'seg_int    ', N_BS_FIX
!   if ( N_BS_FLO > 0 ) write (unit, *) 'seg_real   ', N_BS_FLO
!
!!    terminate keyword control input
!   write (unit, *) 'quit ! keyword input'
!   close (unit)
!END SUBROUTINE OLD_CONTROL_TO_KEYWORDS
 
SUBROUTINE READ_NODAL_POINT_PROPERTIES
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!  INPUT NODAL POINT INTEGER AND/OR REAL PROPERTIES
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Sys_Properties_Data
 IMPLICIT NONE
 INTEGER :: MAX, I, J, K, I_TMP, OK

! MAX_NP   = NUMBER OF SYSTEM NODES
! FLT_NP   = REAL PROPERTIES OF SYSTEM NODES
! NP_FIX   = INTEGER PROPERTIES OF SYSTEM NODES
! N_HOMO   = 1, IF ALL NODES SAME AS FIRST ONE
! N_NP_FLO = NUMBER OF REAL NODAL PROPERTIES
! N_NP_FIX = NUMBER OF INTEGER NODAL PROPERTIES
 
!==>> DATA TYPE 10, Integer nodal properties.
!  Utilized only when N_NP_FIX > 0 on DATA TYPE 3. 
!    J           = Nodal point number.
!    NP_FIX(J,K) = Integer properties. K=1,N_NP_FIX
!  NOTE *** Reading terminated when J=M, unless N_HOMO = 1. Then 
!   only J=1 is read.  
!
!==>> DATA TYPE 11, Real nodal properties.
!  Utilized only when N_NP_FLO > 0 on DATA TYPE 3. 
!    J           = Nodal point number.
!    FLT_NP(J,K) = real properties. K=1,N_NP_FLO
!  NOTE *** Reading terminated when J=M, unless N_HOMO = 1. Then only
!  J=1 is read.  Omitted values default to zero when N_HOMO=0.

!-->     READ NODAL POINT PROPERTIES

  IF ( N_NP_FIX > 0 .OR. N_NP_FLO > 0 ) THEN ! data supplied
    IF (ECHO_PROP) WRITE (N_PRT, "(/, '***   NODAL  POINT  PROPERTIES   ***')")
    MAX = MAX_NP
    IF ( N_HOMO == 1) MAX = 1
    IF ( N_NP_FIX > 0 ) THEN
      IF (ECHO_PROP) WRITE (N_PRT, 25) N_NP_FIX
      25 FORMAT ('   NODE, ', I2,' PROPERTY & INTEGER_VALUE PAIRS')
      DO I = 1, MAX
        READ (N_CRD, *, IOSTAT=OK) J, (NP_FIX (J, K), K = 1, N_NP_FIX)
!b call at(2189)
        CALL INCREMENT_INPUT_LINE ! Update position count
        CALL READ_STATUS (OK, J, 0)
        IF ( J < 1 .OR. J > MAX ) THEN
          PRINT *,'READ_NODAL_POINT_PROPERTIES: INPUT ERROR'  
          PRINT *,'ENTRY FOR ', J,' IS AN INVALID NODE'
          CALL CITE_INPUT_LINE ! List current input position
          STOP 'READ_NODAL_POINT_PROPERTIES: INTEGER PROPERTY ERROR'
        END IF ! NODE OUT OF BOUNDS
        IF ( I == 1 .AND. MAX > 1 ) THEN ! FILL ARRAY
          DO I_TMP = 1, MAX
            NP_FIX (I_TMP, 1:N_NP_FIX) = NP_FIX (J, 1:N_NP_FIX)
          END DO
        END IF
        IF ( J == MAX_NP) EXIT ! THIS LOOP
      END DO
      DO J = 1, MAX
        IF (ECHO_PROP) WRITE (N_PRT, 30) J, (K, NP_FIX (J, K), K = 1, N_NP_FIX)
        30 FORMAT (I6, (10(1X, I2, 1X, I6)))
      END DO
      IF (ECHO_PROP) WRITE (N_PRT, "('END INTEGER PROPERTIES OF NODES', /)")
    END IF

    IF ( N_NP_FLO > 0 ) THEN
      IF (ECHO_PROP) WRITE (N_PRT, 35) N_NP_FLO
      35 FORMAT ('   NODE, ', I2,' PROPERTY & REAL_VALUE PAIRS')
      DO I = 1, MAX
        READ (N_CRD, *, IOSTAT=OK) J, (FLT_NP (J, K), K = 1, N_NP_FLO)
        CALL INCREMENT_INPUT_LINE ! Update position count
        CALL READ_STATUS (OK, J, 0)
        IF ( J < 1 .OR. J > MAX ) THEN
          PRINT *,'READ_NODAL_POINT_PROPERTIES: INPUT ERROR'  
          PRINT *,'ENTRY FOR ', J,' IS AN INVALID NODE'
          CALL CITE_INPUT_LINE ! List current input position
          STOP 'READ_NODAL_POINT_PROPERTIES: REAL PROPERTY ERROR'
        END IF ! NODE OUT OF BOUNDS
        IF ( I == 1 .AND. MAX > 1 ) THEN ! FILL ARRAY
          DO I_TMP = 1, MAX
            FLT_NP (I_TMP, 1:N_NP_FLO) = FLT_NP (J, 1:N_NP_FLO)
          END DO
        END IF
        IF ( J == MAX_NP) EXIT ! THIS LOOP
      END DO
      DO J = 1, MAX
        IF (ECHO_PROP) WRITE (N_PRT, 60) J, (K, FLT_NP (J, K), K = 1, N_NP_FLO)
        60 FORMAT (I6, 4(1X, I2, 1X, 1PE12.5),      &
           (/, (6X, 4(1X, I2, 1X, 1PE12.5))),       & !b 
           (/, (6X, 4(1X, I2, 1X, 1PE12.5))),       & !b 
           (/, (6X, 4(1X, I2, 1X, 1PE12.5))))
           !b (/, (6X, 4(1X, I2, 1X, 1PE12.5))))
      END DO
      IF (ECHO_PROP) WRITE (N_PRT, "('END OF REAL PROPERTIES OF NODES', /)")
    END IF
  END IF ! nodal data supplied
END SUBROUTINE READ_NODAL_POINT_PROPERTIES

SUBROUTINE READ_ELEMENT_PROPERTIES
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!  INPUT ELEMENT POINT INTEGER AND/OR REAL PROPERTIES
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Sys_Properties_Data
 IMPLICIT NONE
 INTEGER :: MAX, I, J, K, I_TMP, OK

! N_ELEMS  = NUMBER OF ELEMENTS IN SYSTEM
! FLT_EL   = REAL PROPERTIES OF SYSTEM ELEMENTS
! LP_FIX   = INTEGER PROPERTIES OF SYSTEM ELEMENTS
! L_HOMO   = 1, IF ALL ELEMENTS & SEGMENTS SAME AS FIRST
! N_LP_FIX = NUMBER OF INTEGER ELEMENT PROPERTIES
! N_LP_FLO = NUMBER OF REAL ELEMENT PROPERTIES
!
!==>> DATA TYPE 12,   Integer element properties.
!  Utilizied only when N_LP_FIX > 0 on DATA TYPE 3.
!    J           = Element number.
!    LP_FIX(J,K) = Integer properties. K=1,N_LP_FIX
!  NOTE *** Reading terminated when J=N_ELEMS, unless L_HOMO = 1. Then
!  only J=1 is read.  Omitted values default to zero when L_HOMO=0.
!
!==>> DATA TYPE 13, Real element properties.
!  Utilized only when N_LP_FLO > 0 on DATA TYPE 3.
!    J           = Element number.
!    FLT_EL(J,K) = real properties. K=1,N_LP_FLO
!  NOTE *** Reading terminated when J=N_ELEMS, unless L_HOMO = 1. Then
!  only J=1 is read.
!

!-->     READ ELEMENT PROPERTIES

  IF ( N_LP_FIX > 0 .OR. N_LP_FLO > 0 ) THEN ! data supplied
    IF (ECHO_PROP) WRITE (N_PRT, "(/, '***  ELEMENT  PROPERTIES   ***')")
    MAX = N_ELEMS
    IF ( L_HOMO == 1) MAX = 1
    IF ( N_LP_FIX > 0 ) THEN
      IF (ECHO_PROP) WRITE (N_PRT, 45) N_LP_FIX
      45 FORMAT ('ELEMENT, ', I2,' PROPERTY & INTEGER_VALUE PAIRS')
      DO I = 1, MAX
        READ (N_CRD, *, IOSTAT=OK) J, (LP_FIX (J, K), K = 1, N_LP_FIX)
!b call at(2284)
        CALL INCREMENT_INPUT_LINE ! Update position count
        CALL READ_STATUS (OK, J, 0)
        IF ( J < 1 .OR. J > MAX ) THEN
          PRINT *,'READ_ELEMENT_PROPERTIES: INPUT ERROR'  
          PRINT *,'ENTRY FOR ', J,' IS AN INVALID NODE'
          CALL CITE_INPUT_LINE ! List current input position
          STOP 'READ_ELEMENT_PROPERTIES: INTEGER PROPERTY ERROR'
        END IF ! ELEMENT OUT OF BOUNDS
        IF ( I == 1 .AND. MAX > 1 ) THEN ! FILL ARRAY
          DO I_TMP = 1, MAX
            LP_FIX (I_TMP, 1:N_LP_FIX) = LP_FIX (1,1:N_LP_FIX)
          END DO
        END IF
        IF ( J == N_ELEMS ) EXIT ! THIS LOOP
      END DO
      DO J = 1, MAX
        IF (ECHO_PROP) WRITE (N_PRT, 30) J, (K, LP_FIX (J, K), K = 1, N_LP_FIX)
        30 FORMAT (I6, (10(1X, I2, 1X, I6)))
      END DO
      IF (ECHO_PROP) WRITE (N_PRT, "('END INTEGER PROPERTIES OF ELEMENTS', /)")
    END IF

    IF ( N_LP_FLO > 0 ) THEN
      IF (ECHO_PROP) WRITE (N_PRT, 55) N_LP_FLO
      55 FORMAT ('ELEMENT, ', I2,' PROPERTY & REAL_VALUE PAIRS')
      DO I = 1, MAX
        READ (N_CRD, *, IOSTAT=OK) J, (FLT_EL (J, K), K = 1, N_LP_FLO)
        CALL INCREMENT_INPUT_LINE ! Update position count
        CALL READ_STATUS (OK, J, 0)
        IF ( J < 1 .OR. J > MAX ) THEN
          PRINT *,'READ_ELEMENT_PROPERTIES: INPUT ERROR'  
          PRINT *,'ENTRY FOR ', J,' IS AN INVALID NODE'
          CALL CITE_INPUT_LINE ! List current input position
          STOP 'READ_ELEMENT_PROPERTIES: REAL PROPERTY ERROR'
        END IF ! ELEMENT OUT OF BOUNDS
        IF ( I == 1 .AND. MAX > 1 ) THEN ! FILL ARRAY
          DO I_TMP = 1, MAX
            FLT_EL (I_TMP, 1: N_LP_FLO) = FLT_EL (J, 1: N_LP_FLO)
          END DO
        END IF
        IF ( J == N_ELEMS ) EXIT ! THIS LOOP
      END DO
      DO J = 1, MAX
        IF (ECHO_PROP) WRITE (N_PRT, 60) J, (K, FLT_EL (J, K), K = 1, N_LP_FLO)
        60 FORMAT (I6, 4(1X, I2, 1X, 1PE12.5),      &
           (/, (6X, 4(1X, I2, 1X, 1PE12.5))),       & !b
           (/, (6X, 4(1X, I2, 1X, 1PE12.5))),       & !b
           (/, (6X, 4(1X, I2, 1X, 1PE12.5))))
           !b (/, (6X, 4(1X, I2, 1X, 1PE12.5))))
      END DO
      IF (ECHO_PROP) WRITE (N_PRT, "('END REAL PROPERTIES OF ELEMENTS', /)")
    END IF
  END IF ! element data supplied
END SUBROUTINE READ_ELEMENT_PROPERTIES

SUBROUTINE READ_MIXED_BC_PROPERTIES
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!  INPUT MIXED BC SEGMENTS INTEGER AND/OR REAL PROPERTIES
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Sys_Properties_Data
 IMPLICIT NONE
 INTEGER :: MAX, I, J, K, I_TMP, OK

! FLT_MX    = ARRAY OF REAL MIXED_SEGMENT PROPERTIES (N_MIXED, N_MX_FLO)
! L_HOMO    = 1, IF ALL ELEMENTS & SEGMENTS SAME AS FIRST
! MX_FIX    = ARRAY OF INTEGER MIXED_SEG PROPERTIES (N_MIXED, N_MX_FIX)
! M_B_N     = NUMBER OF NODES ON A MIXED BC SEGMENT
! N_MIXED   = NUMBER OF SEGMENTS WITH MIXED BOUNDARY CONDITIONS
! N_MX_FIX  = NUMBER OF FIXED POINT MIXED_SEGMENT PROPERTIES
! N_MX_FLO  = NUMBER OF FLOATING POINT MIXED_SEGMENT PROPERTIES

!==>> DATA TYPE 14,   Integer Mixed_BC  segment properties.
!  Utilizied only when N_MX_FIX > 0 in keywords
!    J           = Segment number.
!    MX_FIX(J,K) = Integer properties. K=1,N_MX_FIX
! NOTE *** Reading terminated when J=N_MIXED, unless L_HOMO = 1.
!   Then only J=1 is read.
!
!==>> DATA TYPE 15,   Real Mixed_BC segment properties.
!  Utilizied only when N_MX_FLO > 0 in keywords
!    J           = Segment number.
!    FLT_MX(J,K) = real properties. K=1,N_MX_FLO
! NOTE *** Reading terminated when J=N_MIXED, unless L_HOMO = 1.
!   Then only J=1 is read.

!-->     READ MIXED_BC SEGMENT PROPERTIES

  IF ( N_MX_FIX > 0 .OR. N_MX_FLO > 0 ) THEN ! data supplied
    IF (ECHO_PROP) WRITE (N_PRT, "(/, '***  MIXED_BC SEGMENT PROPERTIES  ***')")
    IF ( ROBIN_THICK /= 1.d0 ) PRINT *, 'NOTE: DEFAULT THICKNESS ',&
       & 'OF ROBIN SEGMENT IS ', ROBIN_THICK
    MAX = N_MIXED
    IF ( L_HOMO == 1) MAX = 1
    IF ( N_MX_FIX > 0 ) THEN
      IF (ECHO_PROP) WRITE (N_PRT, 65) N_MX_FIX
      65 FORMAT ('SEGMENT, ', I2,' PROPERTY & INTEGER_VALUE PAIRS')
      DO I = 1, MAX
        READ (N_CRD, *, IOSTAT=OK) J, (MX_FIX (J, K), K = 1, N_MX_FIX)
        CALL INCREMENT_INPUT_LINE ! Update position count
        CALL READ_STATUS (OK, J, 0)
        IF ( J < 1 .OR. J > MAX ) THEN
          PRINT *,'READ_MIXED_BC_PROPERTIES: INPUT ERROR'  
          PRINT *,'ENTRY FOR ', J,' IS AN INVALID SEGMENT NUMBER'
          CALL CITE_INPUT_LINE ! List current input position
          STOP 'READ_MIXED_BC_PROPERTIES: INTEGER PROPERTY ERROR'
        END IF ! INVALID SEGMENT
        IF ( I == 1 .AND. MAX > 1 ) THEN ! FILL ARRAY
          DO I_TMP = 1, MAX
            MX_FIX (I_TMP, 1:N_MX_FIX) = MX_FIX (J, 1:N_MX_FIX)
          END DO
        END IF
        IF ( J == N_MIXED ) EXIT ! THIS DO
      END DO
      DO J = 1, MAX
        IF (ECHO_PROP) WRITE (N_PRT, 30) J, (K, MX_FIX (J, K), K = 1, N_MX_FIX)
        30 FORMAT (I6, (10(1X, I2, 1X, I6)))
      END DO
      IF (ECHO_PROP) WRITE (N_PRT,"('END INTEGER PROPERTIES OF MIXED_BC SEGMENTS', /)")
    END IF

    IF ( N_MX_FLO > 0 ) THEN
      IF (ECHO_PROP) WRITE (N_PRT, 85) N_MX_FLO
      85 FORMAT ('SEGMENT, ', I2,' PROPERTY & REAL_VALUE PAIRS')
      DO I = 1, MAX
        READ (N_CRD, *, IOSTAT=OK) J, (FLT_MX (J, K), K = 1, N_MX_FLO)
        CALL INCREMENT_INPUT_LINE ! Update position count
        CALL READ_STATUS (OK, J, 0)
        IF ( J < 1 .OR. J > MAX ) THEN
          PRINT *,'READ_MIXED_BC_PROPERTIES: INPUT ERROR'  
          PRINT *,'ENTRY FOR ', J,' IS AN INVALID SEGMENT NUMBER'
          CALL CITE_INPUT_LINE ! List current input position
          STOP 'READ_MIXED_BC_PROPERTIES: REAL PROPERTY ERROR'
        END IF ! INVALID SEGMENT
        IF ( I == 1 .AND. MAX > 1 ) THEN ! FILL ARRAY
          DO I_TMP = 1, MAX
            FLT_MX (I_TMP, 1:N_MX_FLO) = FLT_MX (J, 1:N_MX_FLO)
          END DO
        END IF
        IF ( J == N_MIXED ) EXIT ! THIS DO
      END DO
      DO J = 1, MAX
        IF (ECHO_PROP) WRITE (N_PRT, 60) J, (K, FLT_MX (J, K), K = 1, N_MX_FLO)
        60 FORMAT (I6, 4(1X, I2, 1X, 1PE12.5),      &
           (/, (6X, 4(1X, I2, 1X, 1PE12.5))),       & !b
           (/, (6X, 4(1X, I2, 1X, 1PE12.5))),       & !b
           (/, (6X, 4(1X, I2, 1X, 1PE12.5))))
           !b (/, (6X, 4(1X, I2, 1X, 1PE12.5))))
      END DO
      IF (ECHO_PROP) WRITE (N_PRT,"('END REAL PROPERTIES OF MIXED_BC SEGMENTS', /)")
    END IF
  END IF ! segment data supplied
END SUBROUTINE READ_MIXED_BC_PROPERTIES

SUBROUTINE READ_FLUX_SEGMENT_PROPERTIES
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!  INPUT FLUX SEGMENT INTEGER AND/OR REAL PROPERTIES
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Sys_Properties_Data
 IMPLICIT NONE
 INTEGER :: MAX, I, J, K, I_TMP, OK

! N_F_SEG  = NUMBER OF BOUNDARY SEGMENTS WITH FLUX
! FLT_SP = REAL PROPERTIES OF FLUX SEGMENTS
! SP_FIX = INTEGER PROPERTIES OF FLUX SEGMENTS
! N_BS_FLO = NUMBER OF REAL SEGMENT PROPERTIES
! N_BS_FIX = NUMBER OF INTEGER SEGMENT PROPERTIES
! L_HOMO   = 1, IF ALL ELEMENTS & SEGMENTS SAME AS FIRST

!==>> DATA TYPE 16,   Integer flux segment properties.
!  Utilizied only when N_BS_FIX > 0 on DATA TYPE 3. 
!    J           = Segment number.
!    SP_FIX(J,K) = Integer properties. K=1,N_BS_FIX
!  NOTE *** Reading terminated when J=N_F_SEG, unless L_HOMO = 1. Then 
!  only J=1 is read.  
!
!==>> DATA TYPE 17,   Real flux segment properties.
!  Utilized only when N_BS_FLO > 0 on DATA TYPE 3. 
!    J           = Segment number.
!    FLT_SP(J,K) = real properties. K=1,N_BS_FLO
!  NOTE *** Reading terminated when J=N_F_SEG, unless L_HOMO = 1. Then 
!  only J=1 is read.

!-->     READ SEGMENT PROPERTIES, AND FLUX DATA

  IF ( N_BS_FIX > 0 .OR. N_BS_FLO > 0 ) THEN ! data supplied
    IF (ECHO_PROP) WRITE (N_PRT, "(/, '***  FLUX SEGMENT PROPERTIES   ***')")
    IF ( FLUX_THICK /= 1.d0 ) PRINT *, 'NOTE: DEFAULT THICKNESS ',&
       & 'OF FLUX SEGMENT IS ', FLUX_THICK
    MAX = N_F_SEG
    IF ( L_HOMO == 1) MAX = 1
    IF ( N_BS_FIX > 0 ) THEN
      IF (ECHO_PROP) WRITE (N_PRT, 65) N_BS_FIX
      65 FORMAT ('SEGMENT, ', I2,' PROPERTY & INTEGER_VALUE PAIRS')
      DO I = 1, MAX
        READ (N_CRD, *, IOSTAT=OK) J, (SP_FIX (J, K), K = 1, N_BS_FIX)
        CALL INCREMENT_INPUT_LINE ! Update position count
        CALL READ_STATUS (OK, J, 0)
        IF ( J < 1 .OR. J > MAX ) THEN
          PRINT *,'READ_FLUX_SEGMENT_PROPERTIES: INPUT ERROR'  
          PRINT *,'ENTRY FOR ', J,' IS AN INVALID SEGMENT NUMBER'
          CALL CITE_INPUT_LINE ! List current input position
          STOP 'READ_FLUX_SEGMENT_PROPERTIES: INTEGER PROPERTY ERROR'
        END IF ! INVALID SEGMENT
        IF ( I == 1 .AND. MAX > 1 ) THEN ! FILL ARRAY
          DO I_TMP = 1, MAX
            SP_FIX (I_TMP, 1:N_BS_FIX) = SP_FIX (J, 1:N_BS_FIX)
          END DO
        END IF
        IF ( J == N_F_SEG ) EXIT ! THIS DO
      END DO
      DO J = 1, MAX
        IF (ECHO_PROP) WRITE (N_PRT, 30) J, (K, SP_FIX (J, K), K = 1, N_BS_FIX)
        30 FORMAT (I6, (10(1X, I2, 1X, I6)))
      END DO
      IF (ECHO_PROP) WRITE (N_PRT, "('END INTEGER PROPERTIES OF SEGMENTS', /)")
    END IF

    IF ( N_BS_FLO > 0 ) THEN
      IF (ECHO_PROP) WRITE (N_PRT, 85) N_BS_FLO
      85 FORMAT ('SEGMENT, ', I2,' PROPERTY & REAL_VALUE PAIRS')
      DO I = 1, MAX
        READ (N_CRD, *, IOSTAT=OK) J, (FLT_SP (J, K), K = 1, N_BS_FLO)
        CALL INCREMENT_INPUT_LINE ! Update position count
        CALL READ_STATUS (OK, J, 0)
        IF ( J < 1 .OR. J > MAX ) THEN
          PRINT *,'READ_FLUX_SEGMENT_PROPERTIES: INPUT ERROR'  
          PRINT *,'ENTRY FOR ', J,' IS AN INVALID SEGMENT NUMBER'
          CALL CITE_INPUT_LINE ! List current input position
          STOP 'READ_FLUX_SEGMENT_PROPERTIES: REAL PROPERTY ERROR'
        END IF ! INVALID SEGMENT
        IF ( I == 1 .AND. MAX > 1 ) THEN ! FILL ARRAY
          DO I_TMP = 1, MAX
            FLT_SP (I_TMP, 1:N_BS_FLO) = FLT_SP (J, 1:N_BS_FLO)
          END DO
        END IF
        IF ( J == N_F_SEG ) EXIT ! THIS DO
      END DO
      DO J = 1, MAX
        IF (ECHO_PROP) WRITE (N_PRT, 60) J, (K, FLT_SP (J, K), K = 1, N_BS_FLO)
        60 FORMAT (I6, 4(1X, I2, 1X, 1PE12.5),      &
           (/, (6X, 4(1X, I2, 1X, 1PE12.5))),       & !b
           (/, (6X, 4(1X, I2, 1X, 1PE12.5))),       & !b
           (/, (6X, 4(1X, I2, 1X, 1PE12.5))))
           !b (/, (6X, 4(1X, I2, 1X, 1PE12.5))))
      END DO
      IF (ECHO_PROP) WRITE (N_PRT, "('END REAL PROPERTIES OF SEGMENTS', /)")
    END IF
  END IF ! segment data supplied
END SUBROUTINE READ_FLUX_SEGMENT_PROPERTIES

SUBROUTINE READ_EXACT_SOLUTION_PROPERTIES
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!  INPUT EXACT_SOLUTION INTEGER AND/OR REAL PROPERTIES
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Sys_Properties_Data
 IMPLICIT NONE
 INTEGER :: K, OK

! FLT_EXACT = REAL EXACT_SOLUTION PROPERTIES
! EXACT_FIX = INTEGER EXACT_SOLUTION PROPERTIES
! EXACT_FL  = NUMBER OF EXACT. REAL SYSTEM PROPERTIES
! EXACT_FX  = NUMBER OF EXACT. INTEGER SYSTEM PROPERTIES

!==>> DATA TYPE NEW, Integer exact solution properties.
!  Utilized only if EXACT_FX > 0 on DATA TYPE 3.
!    EXACT_FIX(K) = Exact fixed point properties. K=1,EXACT_FX
!
!==>> DATA TYPE NEW,   Real exact solution properties.
!  Utilized only when EXACT_FL > 0 on DATA TYPE 3.
!    FLT_EXACT(K) = Exact floating point properties. K=1,EXACT_FL

!-->     READ EXACT. SYSTEM PROPERTIES

  IF ( EXACT_FX > 0 .OR. EXACT_FL > 0 ) THEN ! data supplied
    IF (ECHO_PROP) WRITE (N_PRT, "(/, '***  EXACT_SOLUTION SYSTEM PROPERTIES  ***')")
    IF ( EXACT_FX > 0 ) THEN
      IF (ECHO_PROP) WRITE (N_PRT, "('PROPERTY    INTEGER_VALUE')")
      READ (N_CRD, *, IOSTAT=OK) (EXACT_FIX (K), K = 1, EXACT_FX)
      CALL INCREMENT_INPUT_LINE ! Update position count
      CALL READ_STATUS (OK, K, 0)
      DO K = 1, EXACT_FX
        IF (ECHO_PROP) WRITE (N_PRT, 40) K, EXACT_FIX (K)
        40 FORMAT ( I8, 3X, I10 )
      END DO
      IF (ECHO_PROP) WRITE (N_PRT, "('END OF EXACT INTEGER PROPERTIES OF SYSTEM')")
    END IF
    IF ( EXACT_FL > 0 ) THEN
      IF (ECHO_PROP) WRITE (N_PRT, "('PROPERTY    REAL_VALUE')")
      READ (N_CRD, *, IOSTAT=OK) (FLT_EXACT (K), K = 1, EXACT_FL)
      CALL INCREMENT_INPUT_LINE ! Update position count
      CALL READ_STATUS (OK, K, 0)
      DO K = 1, EXACT_FL
        IF (ECHO_PROP) WRITE (N_PRT, 70) K, FLT_EXACT (K)
        70 FORMAT ( I8, 3X, 1PE13.5 )
      END DO
      IF (ECHO_PROP) WRITE (N_PRT, "('END OF EXACT REAL PROPERTIES OF SYSTEM')")
    END IF
  END IF ! exact solution data supplied
END SUBROUTINE READ_EXACT_SOLUTION_PROPERTIES

SUBROUTINE READ_MISCELLANEOUS_PROPERTIES
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!  INPUT MISCELLANEOUS INTEGER AND/OR REAL PROPERTIES
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Sys_Properties_Data
 IMPLICIT NONE
 INTEGER :: K, OK

! FLT_MISC = REAL MISCELLANEOUS PROPERTIES
! MISC_FIX = INTEGER MISCELLANEOUS PROPERTIES
! MISC_FL  = NUMBER OF MISC. REAL SYSTEM PROPERTIES
! MISC_FX  = NUMBER OF MISC. INTEGER SYSTEM PROPERTIES

!==>> DATA TYPE 18, Integer miscellaneous properties.
!  Utilized only if MISC_FX > 0 on DATA TYPE 3.
!    MISC_FIX(K) = Misc. fixed point properties. K=1,MISC_FX
!
!==>> DATA TYPE 19,   Real miscellaneous properties.
!  Utilized only when MISC_FL > 0 on DATA TYPE 3.
!    FLT_MISC(K) = Misc. floating point properties. K=1,MISC_FL

!-->     READ MISC. SYSTEM PROPERTIES

  IF ( MISC_FX > 0 .OR. MISC_FL > 0 ) THEN ! data supplied
    IF (ECHO_PROP) WRITE (N_PRT, "(/, '***  MISCELLANEOUS SYSTEM PROPERTIES  ***')")
    IF ( MISC_FX > 0 ) THEN
      IF (ECHO_PROP) WRITE (N_PRT, "('PROPERTY    INTEGER_VALUE')")
      READ (N_CRD, *, IOSTAT=OK) (MISC_FIX (K), K = 1, MISC_FX)
      CALL INCREMENT_INPUT_LINE ! Update position count
      CALL READ_STATUS (OK, K, 0)
      DO K = 1, MISC_FX
        IF (ECHO_PROP) WRITE (N_PRT, 40) K, MISC_FIX (K)
        40 FORMAT ( I8, 3X, I10 )
      END DO
      IF (ECHO_PROP) WRITE (N_PRT, "('END OF MISC INTEGER PROPERTIES OF SYSTEM')")
    END IF
    IF ( MISC_FL > 0 ) THEN
      IF (ECHO_PROP) WRITE (N_PRT, "('PROPERTY    REAL_VALUE')")
      READ (N_CRD, *, IOSTAT=OK) (FLT_MISC (K), K = 1, MISC_FL)
      CALL INCREMENT_INPUT_LINE ! Update position count
      CALL READ_STATUS (OK, K, 0)
      DO K = 1, MISC_FL
        IF (ECHO_PROP) WRITE (N_PRT, 70) K, FLT_MISC (K)
        70 FORMAT ( I8, 3X, 1PE13.5 )
      END DO
      IF (ECHO_PROP) WRITE (N_PRT, "('END OF MISC. REAL PROPERTIES OF SYSTEM')")
    END IF
  END IF ! miscellaneous data supplied
END SUBROUTINE READ_MISCELLANEOUS_PROPERTIES

SUBROUTINE READ_MATERIAL_PROPERTIES
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!  INPUT MATERIAL REAL PROPERTIES
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Sys_Properties_Data
 IMPLICIT NONE
 INTEGER :: MAX, I, J, K, I_TMP, OK

! FLT_MAT  = REAL PROPERTIES OF SYSTEM MATERIALS, (N_MAT, MAT_FLO)
! N_MAT    = NUMBER OF MATERIAL TYPES
! MAT_FLO  = NUMBER OF REAL MATERIAL PROPERTIES
 
!==>> DATA TYPE 20, Real material properties.
!  Utilized only when N_MAT and  MAT_FLO > 0 in keywords.
!    J             = Material number.
!    FLT_MAT(J,K) = real properties. K=1,MAT_FLO
!  NOTE *** Reading terminated when J=N_MAT.

!-->     READ MATERIAL PROPERTIES

    IF ( MAT_FLO > 0 ) THEN
      IF (ECHO_PROP) WRITE (N_PRT, "(/, '***  MATERIAL PROPERTIES   ***')")
      MAX = N_MAT
      IF (ECHO_PROP) WRITE (N_PRT, 55) MAT_FLO
      55 FORMAT ('MATERIAL, ', I2,' PROPERTY & REAL_VALUE PAIRS')
      DO I = 1, MAX
        READ (N_CRD, *, IOSTAT=OK) J, (FLT_MAT (J, K), K = 1, MAT_FLO)
        CALL INCREMENT_INPUT_LINE ! Update position count
        CALL READ_STATUS (OK, J, K)
        IF ( J < 1 .OR. J > MAX ) THEN
          PRINT *,'READ_MATERIAL_PROPERTIES: INPUT ERROR'  
          PRINT *,'ENTRY FOR ', J,' IS AN INVALID MATERIAL NUMBER'
          CALL CITE_INPUT_LINE ! List current input position
          STOP 'READ_MATERIAL_PROPERTIES: REAL PROPERTY ERROR'
        END IF ! INVALID MATERIAL
        IF ( I == 1 .AND. MAX > 1 ) THEN ! FILL ARRAY
          DO I_TMP = 1, MAX
            FLT_MAT (I_TMP, 1: MAT_FLO) = FLT_MAT (J, 1: MAT_FLO)
          END DO
        END IF
        IF ( J == N_MAT ) EXIT ! THIS LOOP
      END DO
      DO J = 1, MAX
        IF (ECHO_PROP) WRITE (N_PRT, 60) J, (K, FLT_MAT (J, K), K = 1, MAT_FLO)
        60 FORMAT (I6, 4(1X, I2, 1X, 1PE12.5),      &
           (/, (6X, 4(1X, I2, 1X, 1PE12.5))),       & !b
           (/, (6X, 4(1X, I2, 1X, 1PE12.5))),       & !b
           (/, (6X, 4(1X, I2, 1X, 1PE12.5))))
           !b (/, (6X, 4(1X, I2, 1X, 1PE12.5))))
      END DO
      IF (ECHO_PROP) WRITE (N_PRT, "('END REAL PROPERTIES OF MATERIALS', /)")
  END IF ! material data supplied
END SUBROUTINE READ_MATERIAL_PROPERTIES
