! copyright 2005,  J. E. Akin, all rights reserved.
! file: printing_lib.f

SUBROUTINE AT (LINE)
! * * * * * * * * * * * * * * * * * * * * *
! DEBUGGING AID, PRINT AN INTEGER
! * * * * * * * * * * * * * * * * * * * * *
 INTEGER, INTENT(IN) :: LINE
    WRITE (6, *) ' -------->>> HERE AT ', LINE
END SUBROUTINE AT
                 
SUBROUTINE NODAL_BAR_CHART (IBAR, I_PARM, X, D)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! PRINT-PLOT BAR CHARTS OF NODAL PARAMETER I_PARM AT
! THE NODES IN ARRAY NODE AND SCALE THE RELATIVE
! DISTANCE BETWEEN THE POINTS
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 INTEGER,  PARAMETER  :: MID = 25, LINE = 2 * MID
 REAL(DP), INTENT(IN) :: X (MAX_NP, N_SPACE), D (N_D_FRE)
 INTEGER,  INTENT(IN) :: IBAR, I_PARM
!           Automatic Arrays
 INTEGER :: NODE_BAR (IBAR)
!            Locals
 INTEGER :: I, IS, INDEX, IPMNG, ISPACE, J, JY, JZ, &
             K, KOUNT, L, LIMIT, NLAST, NOPLOT
 REAL(DP)    :: DIST, DMIN, DMINSQ, DTEST, RANGE, YMAX, YMIN
 CHARACTER ALINE (LINE+1), SKIP (LINE+1)
 CHARACTER(1), SAVE :: BLANK = ' ', DOT = '+', DASH = '-', &
                       A_X = 'X', PLUS = '+', A_0 = '0'
  DATA NOPLOT / 0 /
! MAX_NP   = TOTAL NUMBER OF NODES IN SYSTEM
! I_PARM   = NODAL PARAMETER TO BE GRAPHED,1<=I_PARM<=N_G_DOF
! NODE_BAR = LIST OF NODES TO BE USED (WHEN IBAR > 1)
! NO_DIST  = 1  OMIT DISTANCES BETWEEN NODAL BAR LINES
! X        = ARRAY OF GLOBAL COORDINATES OF ALL NODES
! N_D_FRE  = TOTAL NO. OF DEGREES OF FREEDOM IN SYSTEM
! D        = ARRAY OF ALL NODAL PARAMETERS IN THE SYSTEM
! N_G_DOF  = NUMBER OF PARAMETERS PER NODE
! N_SPACE  = DIMENSION OF SOLUTION SPACE
! IBAR     = NUMBER OF NODES TO BE INCLUDED IN BAR CHART
!            IF IBAR=1  USE ALL NODES, NODE_BAR NOT USED
  NOPLOT = NOPLOT + 1
  LIMIT = IBAR
  IF ( IBAR == 1) LIMIT = MAX_NP
  WRITE (N_PRT, 5000) NOPLOT, DOF_NAME (I_PARM), LIMIT
  5000 FORMAT ('*** PRINT PLOT NUMBER',I3,' ***',/,   &
               A12, ' EVALUATED AT ',I5,' NODE POINTS')
  !b IPMNG = I_PARM - N_G_DOF
  IF ( NO_DIST == 0) THEN
! FIND MINIMUM DISTANCE BETWEEN POINTS
    DMINSQ = 0.d0
    DO IS = 1, N_SPACE
      DMINSQ = DMINSQ + (X (2, IS) - X (1, IS) ) **2
    END DO
    DO J = 2, LIMIT
      I = J
      IF ( IBAR > 1) I = NODE_BAR (J)
      DTEST = 0.0
      DO IS = 1, N_SPACE
        DTEST = DTEST + (X (I, IS) - X (I - 1, IS) ) **2
      END DO
      IF ( DTEST < DMINSQ) DMINSQ = DTEST
    END DO
    DMIN = SQRT (DMINSQ)
  END IF
!-->  ESTABLISH GRAPH LIMITING VALUES
!   INITALIZE MAX MIN VALUES OF PARAMETER
  L = 1
  IF ( IBAR /= 1) L = NODE_BAR (1)
  !b INDEX = N_G_DOF * L + IPMNG
  INDEX = N_G_DOF * (L - 1) + I_PARM
  YMAX = D (INDEX)
  YMIN = YMAX
  DO I = 1, LIMIT
    L = I
    IF ( IBAR /= 1) L = NODE_BAR (I)
    !b INDEX = N_G_DOF * L + IPMNG
    INDEX = N_G_DOF * (L - 1) + I_PARM
    DTEST = D (INDEX)
    IF ( DTEST > YMAX) YMAX = DTEST
    IF ( DTEST < YMIN) YMIN = DTEST
  END DO
  KOUNT = 1
  WRITE (N_PRT, 5010) YMIN, YMAX
  5010 FORMAT (' RANGE ON GRAPH IS ',1PE12.5,' TO ',1PE12.5/)
  IF ( YMIN == YMAX ) THEN
    WRITE (N_PRT, 5011)
    5011 FORMAT ('CONSTANT VALUE, PLOT SKIPPED')
    RETURN
  END IF
  RANGE = YMAX - YMIN
!  FLAG ZERO VALUE
  IF ( YMIN < 0.0 .AND. YMAX > 0.0 ) THEN
    JZ = 1*(YMAX - 0.0) / RANGE                               &
       + (LINE+1)*(0.0 - YMIN) / RANGE + 0.1
  ELSE
    JZ = 0
  END IF
! SCALING COMPLETE
  SKIP  (2:LINE) = BLANK
  ALINE (2:LINE) = DASH
  SKIP (1) = PLUS
  SKIP (MID+1) = DOT
  SKIP (LINE+1) = PLUS
  ALINE (1:LINE+1:5) = PLUS
  ALINE (MID+1) = DOT
  WRITE (N_PRT, 5020) ALINE
  5020 FORMAT('    NODE    VALUE   ', (101A1) )
  L = 1
  IF ( IBAR > 1) L = NODE_BAR (1)
  DO K = 1, LIMIT
    NLAST = L
    L = K
    IF ( IBAR > 1) L = NODE_BAR (K)
    !b INDEX = N_G_DOF * L + IPMNG
    INDEX = N_G_DOF * (L - 1) + I_PARM
    DTEST = D (INDEX)
    ALINE (2:LINE) = BLANK
    ALINE (1) = PLUS
    ALINE (LINE+1) = PLUS
    IF ( JZ > 0) ALINE (JZ) = A_0
    ISPACE = 1
    
    IF ( NO_DIST == 0) THEN
!-->      FIND DISTANCE BETWEEN TWO POINTS
      DIST = SUM ((X (L, :) - X (NLAST, :))**2) !9
!     MINIMUM DISTANCE IS ONE SPACE
      DIST = SQRT (DIST)
      IF ( K /= 1 ) THEN
        IF ( DMIN > 0.0 ) THEN
          ISPACE = DIST / DMIN + 0.5
        ELSE
          ISPACE = 1
        END IF
      END IF
      IF ( ISPACE > 10) ISPACE = 5
    END IF
    !b IF ( NO_DIST > 0 .AND. K > 1 ) THEN
    IF ( NO_DIST == 0 .AND. K > 1 ) THEN
      DO I = 1, ISPACE
        WRITE (N_PRT, 5030) SKIP
        5030 FORMAT (20X, (101A1) )
      END DO
    END IF
! LINEAR INTERPOLATION FOR COLUMN NUMBER
    JY = 1 * (YMAX - DTEST) / RANGE                    &
       + (LINE+1) * (DTEST - YMIN) / RANGE + 0.1
    ALINE (1:JY) = A_X
    WRITE (N_PRT, 5040) L, D (INDEX), ALINE
    5040 FORMAT (I7,2X,1PE10.3,1X,101A1)
  END DO
  ALINE (2:LINE) = DASH
  ALINE (1:LINE+1:5) = PLUS
  ALINE (MID+1) = DOT
  WRITE (N_PRT, 5020) ALINE
END SUBROUTINE NODAL_BAR_CHART

SUBROUTINE  G_PRINT (A, NR, NC)
! ------------------------------------------------------
!   F90 PRINTING OF REAL MATRIX A(NR,NC) WITH G FORMAT
! ------------------------------------------------------
Use Precision_Module
 IMPLICIT NONE
 INTEGER, INTENT(IN) :: NC, NR
 REAL(DP), INTENT(IN) :: A(NR*NC)
 INTEGER, PARAMETER  :: MAX = 7
 INTEGER             :: NCOL(MAX), I, IC, J, K, L, MAXCOL
 INTEGER             :: ROW, NH, NL
!
! A    = REAL ARRAY, PASSED BY COLUMNS
! MAX  = MAX NUMBER OF COLUMNS PER SCREEN
! NR   = NUMBER OF ROWS IN A
! NC   = NUMBER OF COLUMNS IN A
!
!   LOOP OVER A SCREEN FULL OF COLUMNS
  DO J = 1, NC, MAX
    MAXCOL = 1
    K      = NC - J - 1
    MAXCOL = MIN0 (K,MAX)
    IF ( NC > MAX )  THEN
!       FILL THE TEMPORARY COLUMN NUMBERS TO PRINT
      DO L = 1, MAXCOL
        NCOL(L) = L + J - 1
      END DO ! OF L COLUMNS
      PRINT '("ROW/COL", I7, 6I12)', ( NCOL(IC),IC=1,MAXCOL )
    END IF ! OF NC VS SCREEN WIDTH
!     PRINT ROWS IN CURRENT SCREEN COLUMNS
    DO ROW = 1, NR
      NL = ROW + (J-1)*NR
      NH = NL + (MAXCOL - 1)*NR
      IF ( NC > MAX )  THEN
!         WITH ROW NUMBER WIDER THAN SCREEN
        PRINT '(I3, 1X, 7(1PG12.5))', ROW,( A(I),I=NL,NH,NR )
      ELSE ! NC
!         WITHOUT ROW NUMBER FITS SCREEN
        PRINT '(7(1PG12.5))', ( A(I),I=NL,NH,NR )
      END IF ! OF NC VS SCREEN WIDTH
    END DO ! OF ROW OVER ROWS
  END DO ! OF J
  PRINT *, ' '
  END SUBROUTINE G_PRINT
                  
SUBROUTINE IPRINT (m, NR, NC)
! * * * * * * * * * * * * * * * * * * * * * * *
!    PRINTING OF AN INTEGER ARRAY M(NR,NC)
! * * * * * * * * * * * * * * * * * * * * * * *
 IMPLICIT NONE
 INTEGER, PARAMETER  :: N_PRT = 6, MAX = 10
 INTEGER, INTENT(IN) :: NR, NC
 INTEGER, INTENT(IN) :: m (nr*nc)
!          Automatic Arrays
 INTEGER :: NCOL (MAX)
 INTEGER :: I, IC, IR, J, L, MAXCOL, NH, NL
! m    = INTEGER ARRAY
! NR   = NUMBER OF ROWS IN m
! NC   = NUMBER OF COLUMNS IN m
  DO J = 1, NC, MAX
    MAXCOL = MIN0 ((NC - (J - 1)), MAX)
    DO L = 1, MAXCOL
      NCOL (L) = L + (J - 1)
    END DO
    WRITE (N_PRT, 5000) (NCOL (IC), IC = 1, MAXCOL)
    5000 FORMAT ('ROW/COL',I7, 9I10 )
    DO IR = 1, NR
      NL = IR + (J - 1) * NR
      NH = NL + (MAXCOL - 1) * NR
      WRITE (N_PRT, 5010) IR, (m (I), I = NL, NH, NR)
      5010 FORMAT ( I4, 10I10 )
    END DO
  END DO
END SUBROUTINE IPRINT
                        
SUBROUTINE MAX_AND_MIN_RANGE (DD)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!  LIST EXTREME RANGE OF VALUES OF THE N_G_DOF NODAL DOF
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Interface_Header
 IMPLICIT NONE
 REAL(DP), INTENT(IN) :: DD (N_D_FRE)
!           Automatic Arrays
 REAL(DP) :: RANGE   (N_G_DOF, 2), DD_TEST
 INTEGER  :: N_RANGE (N_G_DOF, 2), INDEX (N_G_DOF)
 INTEGER  :: J, I

! MAX_NP  = NUMBER OF NODES IN SYSTEM
! N_G_DOF = NUMBER OF PARAMETERS (DOF) PER NODE
! N_D_FRE = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM
! RANGE   : 1-MAXIMUM VALUE, 2-MINIMUM VALUE
! DD      = ARRAY OF SYSTEM DEGREES OF FREEDOM
! INDEX   = LIST OF SYSTEM DOF NOS FOR DOF AT NODE
! N_RANGE = ARRAY OF NODE NOS OF EXTREME VALUE POINTS

  RANGE   (1:N_G_DOF, 1) = DD (1:N_G_DOF)  ! Initialize
  RANGE   (1:N_G_DOF, 2) = DD (1:N_G_DOF)  ! Initialize
  N_RANGE (1:N_G_DOF, 1) = 1               ! Initialize
  N_RANGE (1:N_G_DOF, 2) = 1               ! Initialize

  DO I = 1, MAX_NP
    INDEX = GET_INDEX_AT_PT (I) 
    DO J = 1, N_G_DOF
      DD_TEST = DD (INDEX (J) )
      IF ( DD_TEST > RANGE (J, 1) ) THEN ! CURRENT MAXIMUM
        RANGE (J, 1)  = DD_TEST ; N_RANGE (J, 1) = I
        CYCLE ! TO NEXT J
      ELSE IF ( DD_TEST < RANGE (J, 2) ) THEN ! CURRENT MINIMUM
        RANGE (J, 2)  = DD_TEST ; N_RANGE (J, 2) = I
      END IF
    END DO
  END DO

!    PRINT RANGE OF VALUES
  !b WRITE (U_HIST, 5000) 
  !b WRITE (U_LOG , 5000) 
  WRITE (N_PRT,  5000) ; 5000 FORMAT   ( /,               &
    '*** EXTREME VALUES OF THE NODAL PARAMETERS ***',/,   &
    'PARAMETER         MAXIMUM,   NODE      MINIMUM,   NODE')
  DO J = 1, N_G_DOF
    !b WRITE (U_HIST, 10) DOF_NAME(J), RANGE (J, 1), N_RANGE (J, 1),&
                                    !b RANGE (J, 2), N_RANGE (J, 2)
    !b WRITE (U_LOG , 10) DOF_NAME(J), RANGE (J, 1), N_RANGE (J, 1),&
                                    !b RANGE (J, 2), N_RANGE (J, 2)
    WRITE (N_PRT,  10) DOF_NAME(J), RANGE (J, 1), N_RANGE (J, 1),&
                                    RANGE (J, 2), N_RANGE (J, 2)
    10 FORMAT (A12, 2X, 1PE11.4, ',', I7, 2X, 1PE11.4, ',', I7)
  END DO
   
  IF ( .NOT. IS_NULL ) THEN ! Expect non-zero solutions
    DD_TEST = 2 * TINY (DD_TEST)
    IF ( MAXVAL (ABS (RANGE)) < DD_TEST )  THEN ! Abort
         N_WARN = N_WARN + 1
         PRINT *, 'WARNING: MAX_AND_MIN_RANGE NULL SOLUTION, MUST' 
         PRINT *, 'DELETE el_no_col, OR ADD is_null IN KEYWORDS'
         PRINT *, 'ERROR, NULL SOLUTION, CURRENT CONTROLS ARE:'
         CALL LIST_CONTROLS
         STOP 'ERROR, NULL SOLUTION. ADD is_null IN KEYWORDS'
    END IF ! Abort
  END IF ! not expecting null solution
END SUBROUTINE MAX_AND_MIN_RANGE
                                
SUBROUTINE MAX_AND_MIN_SCP_AVE_F90 (SCP_AVERAGES)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!  LIST EXTREME RANGE OF VALUES OF THE SCP AVERAGES     
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Interface_Header
 IMPLICIT NONE
 REAL(DP), INTENT(IN) :: SCP_AVERAGES (MAX_NP, SCP_FIT)
!           Automatic Arrays
 REAL(DP) :: RANGE   (SCP_FIT, 2), AVE_TEST
 INTEGER  :: N_RANGE (SCP_FIT, 2)
 INTEGER  :: J, I

! MAX_NP       = NUMBER OF NODES IN SYSTEM
! SCP_FIT      = NUMBER OF AVERAGED FLUX ITEMS PER NODE
! N_D_FRE      = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM
! RANGE        : 1-MAXIMUM VALUE, 2-MINIMUM VALUE
! SCP_AVERAGES = ARRAY OF SYSTEM DEGREES OF FREEDOM
! N_RANGE      = ARRAY OF NODE NUMBERS OF EXTREME VALUE POINTS

  RANGE   (1:SCP_FIT, 1) = SCP_AVERAGES (1, 1:SCP_FIT) ! Initialize
  RANGE   (1:SCP_FIT, 2) = SCP_AVERAGES (1, 1:SCP_FIT) ! Initialize
  N_RANGE (1:SCP_FIT, 1) = 1                           ! Initialize
  N_RANGE (1:SCP_FIT, 2) = 1                           ! Initialize

  DO J = 1, SCP_FIT
    DO I = 1, MAX_NP
      AVE_TEST = SCP_AVERAGES (I, J)      
      IF ( AVE_TEST > RANGE (J, 1) ) THEN ! CURRENT MAXIMUM
        RANGE (J, 1)  = AVE_TEST ; N_RANGE (J, 1) = I
      ELSE IF ( AVE_TEST < RANGE (J, 2) ) THEN ! CURRENT MINIMUM
        RANGE (J, 2)  = AVE_TEST ; N_RANGE (J, 2) = I
      END IF
    END DO ! on nodes
  END DO ! on items

!    PRINT RANGE OF VALUES
  !b WRITE (U_HIST, 5000) 
  !b WRITE (U_LOG , 5000) 
  WRITE (N_PRT,  5000) ; 5000 FORMAT   ( /,                       &
    '*** EXTREME VALUES OF THE NODAL SCP_AVERAGE  ITEMS ***',/,   &
    'ITEM            MAXIMUM,   NODE      MINIMUM,   NODE')
  DO J = 1, SCP_FIT
    !b WRITE (U_HIST, 10) AVE_NAME (J), RANGE (J, 1), N_RANGE (J, 1),&
                                     !b RANGE (J, 2), N_RANGE (J, 2)
    !b WRITE (U_LOG , 10) AVE_NAME (J), RANGE (J, 1), N_RANGE (J, 1),&
                                     !b RANGE (J, 2), N_RANGE (J, 2)
    WRITE (N_PRT,  10) AVE_NAME (J), RANGE (J, 1), N_RANGE (J, 1),&
                                     RANGE (J, 2), N_RANGE (J, 2)
    10 FORMAT (A12, 1PE11.4, ',', I7, 2X, 1PE11.4, ',', I7)
  END DO

  AVE_TEST = 2 * TINY (AVE_TEST)               ! Small value
  IF ( MAXVAL (ABS (RANGE)) < AVE_TEST )  THEN ! Warn 
       PRINT *, 'WARNING: MAX_AND_MIN_SCP_AVE, ALL AVERAGES ARE ZERO' 
       PRINT *, '         IN F90 VERSION'
       PRINT *, '         ERROR ESTIMATOR WILL BE SKIPPED'
       NULL_SCP_AVE = .TRUE.
!b print *, 'NULL_SCP_AVE = ', NULL_SCP_AVE !b
       SKIP_ERROR   = .TRUE.
       N_WARN = N_WARN + 1                 
  END IF ! Warn 
END SUBROUTINE MAX_AND_MIN_SCP_AVE_F90
                                
SUBROUTINE MAX_AND_MIN_SCP_AVE_F95 (SCP_AVERAGES)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!  LIST EXTREME RANGE OF VALUES OF THE SCP AVERAGES     
!  USING F95 EXTENSION OF MAXLOC INTRINSIC FUNCTION
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Interface_Header
 IMPLICIT NONE
 REAL(DP), INTENT(IN) :: SCP_AVERAGES (MAX_NP, SCP_FIT)

 REAL(DP) :: MAX_AVE_VALUE (SCP_FIT), MIN_AVE_VALUE (SCP_FIT)
 INTEGER  :: MAX_AVE_NODE  (SCP_FIT), MIN_AVE_NODE  (SCP_FIT)
 REAL(DP) :: AVE_TEST
 INTEGER  :: J

! MAX_NP       = NUMBER OF NODES IN SYSTEM
! SCP_FIT      = NUMBER OF AVERAGED FLUX ITEMS PER NODE
! SCP_AVERAGES = ARRAY OF SYSTEM DEGREES OF FREEDOM

   MAX_AVE_VALUE = MAXVAL (SCP_AVERAGES, DIM = 1)
   MAX_AVE_NODE  = MAXLOC (SCP_AVERAGES, DIM = 1)
   MIN_AVE_VALUE = MINVAL (SCP_AVERAGES, DIM = 1)
   MIN_AVE_NODE  = MINLOC (SCP_AVERAGES, DIM = 1)

!    PRINT RANGE OF VALUES
  !b WRITE (U_HIST, 5000) 
  !b WRITE (U_LOG , 5000) 
  WRITE (N_PRT,  5000) ; 5000 FORMAT   ( /,                      &
    '*** EXTREME VALUES OF THE NODAL SCP_AVERAGE  ITEMS ***',/,  &
    'ITEM            MAXIMUM,   NODE      MINIMUM,   NODE')
  DO J = 1, SCP_FIT
    !b WRITE (U_HIST, 6) AVE_NAME (J),                        &
                      !b MAX_AVE_VALUE (J), MAX_AVE_NODE (J), &
                      !b MIN_AVE_VALUE (J), MIN_AVE_NODE (J)
    !b WRITE (U_LOG , 6) AVE_NAME (J),                        &
                      !b MAX_AVE_VALUE (J), MAX_AVE_NODE (J), &
                      !b MIN_AVE_VALUE (J), MIN_AVE_NODE (J)
    WRITE (N_PRT,  6) AVE_NAME (J),                        &
                      MAX_AVE_VALUE (J), MAX_AVE_NODE (J), &
                      MIN_AVE_VALUE (J), MIN_AVE_NODE (J)
    6 FORMAT (A12, 1PE11.4, ',', I7, 2X, 1PE11.4, ',', I7)
  END DO

  AVE_TEST = 2 * TINY (AVE_TEST)               ! Small value
  IF ( (ABS (MAXVAL (MAX_AVE_VALUE)) < AVE_TEST)    .AND.     &
       (ABS (MINVAL (MIN_AVE_VALUE)) < AVE_TEST) )  THEN ! Warn 
    PRINT *, 'WARNING: MAX_AND_MIN_SCP_AVE, ALL AVERAGES ARE ZERO' 
    PRINT *, '         IN F95 VERSION'
    PRINT *, '         ERROR ESTIMATOR WILL BE SKIPPED'
    NULL_SCP_AVE = .TRUE.
    SKIP_ERROR   = .TRUE.
    N_WARN = N_WARN + 1                 
  END IF ! Warn 
END SUBROUTINE MAX_AND_MIN_SCP_AVE_F95
                                
SUBROUTINE RPRINT (A, NR, NC, I_OPT)
! * * * * * * * * * * * * * * * * * * * * * * *
!      PRINTING OF REAL MATRIX A(NR,NC)
! * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 IMPLICIT NONE
 INTEGER, PARAMETER   :: N_PRT = 6, MAX = 10
 INTEGER, INTENT(IN)  :: NR, NC, I_OPT
 REAL(DP), INTENT(IN) :: A (nr*nc)
!          Automatic Arrays
 INTEGER :: NCOL (MAX)
 INTEGER :: I, IC, IR, J, L, MAXCOL, NH, NL
! A    = REAL ARRAY
! NR   = NUMBER OF ROWS IN A
! NC   = NUMBER OF COLUMNS IN A
!b! I_OPT = 0 USE F FORMAT, OTHERWISE USE E FORMAT
! I_OPT = 0 F FORMAT, 1 E FORMAT, 2 F LONG, 3 E LONG
  DO J = 1, NC, MAX
    MAXCOL = MIN0 ((NC - (J - 1)), MAX)
    DO L = 1, MAXCOL
      NCOL (L) = L + (J - 1)
    END DO
    IF ( I_OPT < 2 ) THEN !b
      WRITE (N_PRT, 5000) (NCOL (IC), IC = 1, MAXCOL)
      5000 FORMAT ('ROW/COL', I7, 9I10 )
    ELSE
      WRITE (N_PRT, 5001) (NCOL (IC), IC = 1, MAXCOL)
      5001 FORMAT ('ROW/COL', I11, 9I14 )
    END IF
    DO IR = 1, NR
      NL = IR + (J - 1) * NR
      NH = NL + (MAXCOL - 1)* NR
      IF ( I_OPT == 0 ) THEN
        WRITE (N_PRT, 5010) IR, (A (I), I = NL, NH, NR)
        5010 FORMAT (I4, 10F10.4)
      !b ELSE
      ELSEIF ( I_OPT == 1 ) THEN !b
        WRITE (N_PRT, 5020) IR, (A (I), I = NL, NH, NR)
        5020 FORMAT (I4, 10(1PE10.2) )
      ELSEIF ( I_OPT == 1 ) THEN !b
        WRITE (N_PRT, 5011) IR, (A (I), I = NL, NH, NR) !b
        5011 FORMAT (I4, 10F14.8) !b
      ELSEIF ( I_OPT == 3 ) THEN !b
        WRITE (N_PRT, 5021) IR, (A (I), I = NL, NH, NR) !b
        5021 FORMAT (I4, 10(1PE14.6) ) !b
      END IF
    END DO
  END DO
END SUBROUTINE RPRINT
                      
SUBROUTINE WRITE_BY_ELEMENTS (DD, NODES)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! OUTPUT, BY ELEMS, OF CALCULATED DEGREES OF FREEDOM
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data
Use Interface_Header
 IMPLICIT NONE
 REAL(DP), INTENT(IN) :: DD (N_D_FRE)
 INTEGER, INTENT(IN) :: NODES (L_S_TOT, NOD_PER_EL)
!          Automatic Arrays
 INTEGER :: INDX (N_G_DOF), EL_NODES (NOD_PER_EL) 
 INTEGER :: IE, K, L, NODE

! N_ELEMS    = NUMBER OF ELEMENTS IN SYSTEM
! NOD_PER_EL = NUMBER NODES PER ELEMENT
! N_G_DOF    = NUMBER OF PARAMETERS (DOF) PER NODE
! N_D_FRE    = NUMBER DEGREES OF FREEDOM IN SYSTEM
! DD         = CALCULATED NODAL PARAMETERS (DOF)
! INDX       = SYSTEM DOF NUMBERS FOR ELEMENT PARAMETERS
! NODES      = NODAL INCIDENCES OF ALL ELEMENTS
! EL_NODES   = NODAL INCIDENCES OF AN ELEMENT

  !b WRITE (N_PRT, 5) N_G_DOF
  !b 5 FORMAT (/, '***  OUTPUT OF RESULTS BY ELEMENT ORDER  ***',/, &
               !b 'ELEMENT, NODE, ', I3,' PARAMETERS' )
  IF ( TRANSIENT ) PRINT *,'At TIME = ', TIME
  WRITE (N_PRT, 5) DOF_NAME (1:N_G_DOF)
  5 FORMAT (/, '***  OUTPUT OF RESULTS BY ELEMENT ORDER  ***',/, &
               'ELEMENT, NODE, ', (6A12) )
  DO IE = 1, N_ELEMS
    LT_N = LT_DATA (1, L_TYPE(IE))  !9 use function to get
    EL_NODES = GET_ELEM_NODES (IE, LT_N, NODES)    
    DO K = 1, LT_N
      NODE = EL_NODES (K)
!      ALLOW FOR OMITTED NODES
      IF ( NODE > 0 ) THEN
        INDX = GET_INDEX_AT_PT (NODE) 
        WRITE (N_PRT, 6) IE, NODE, (DD (INDX (L) ), L = 1, N_G_DOF)
        !b 6 FORMAT (I5, I8, 2X, (6(1X, 1PE12.5)) )
        6 FORMAT (I5, I8, 2X, (6(1X, ES11.4)) ) !b
      END IF
    END DO
  END DO
END SUBROUTINE WRITE_BY_ELEMENTS
                                
SUBROUTINE WRITE_BY_NODES (X, DD)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! OUTPUT, BY NODES, OF CALCULATED DEGREES OF FREEDOM
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Interface_Header
 IMPLICIT NONE
 REAL(DP), INTENT(IN) :: X (MAX_NP, N_SPACE), DD (N_D_FRE)
 INTEGER  :: I, K, L
!          Automatic Arrays
 INTEGER  :: INDEX (N_G_DOF)

! MAX_NP  = NUMBER OF NODES IN SYSTEM
! N_G_DOF = NUMBER OF PARAMETERS (DOF) PER NODE
! N_D_FRE = NUMBER OF DOF IN THE SYSTEM
! N_SPACE = DIMENSION OF SPACE
! X       = SYSTEM COORDINATES OF ALL NODES
! DD      = CALCULATED NODAL PARAMETERS
! INDEX   = SYSTEM DOF NOS OF PARAMETERS ON A NODE

  IF ( TRANSIENT ) PRINT *,'At TIME = ', TIME
IF ( N_SPACE==3 .AND. N_G_DOF == 6 ) THEN ! space frame or shell
  WRITE (N_PRT, 30) XYZ_NAME(1:3), DOF_NAME(1:3)
  WRITE (N_PRT, 31)                DOF_NAME(4:6)
  30 FORMAT( /, '***  OUTPUT OF RESULTS IN NODAL ORDER  ***', &
             /, ' NODE, ', (6A12) )
  31 FORMAT(    ' NODE, ', 36X, (3A12) )
  DO I = 1, MAX_NP
    INDEX = GET_INDEX_AT_PT (I) 
    WRITE (N_PRT, 41) I, (X (I, L), L = 1, 3),    &
                         (DD (INDEX (K)), K = 1, 3)
    WRITE (N_PRT, 42) I, (DD (INDEX (K)), K = 4, 6)
    41 FORMAT ( I5, (9(1X, ES11.4)) ) !b
    42 FORMAT ( I5, 36X, (3(1X, ES11.4)) ) !b
  END DO
ELSE
  WRITE (N_PRT, 50) XYZ_NAME(1:N_SPACE), DOF_NAME(1:N_G_DOF)
  50 FORMAT( /, '***  OUTPUT OF RESULTS IN NODAL ORDER  ***', &
       /, '    NODE, ', (9A12) )
  DO I = 1, MAX_NP
    INDEX = GET_INDEX_AT_PT (I) 
    WRITE (N_PRT, 51) I, (X (I, L), L = 1, N_SPACE), &
                        (DD (INDEX (K)), K = 1, N_G_DOF)
    51 FORMAT ( I8, (9(1X, ES11.4)) ) !b
  END DO
END IF ! space frame or shell
END SUBROUTINE WRITE_BY_NODES
                             
SUBROUTINE WRITE_BY_NODES_WITH_EXACT (X, DD)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! OUTPUT, BY NODES, OF CALCULATED DEGREES OF FREEDOM
!            AND THE EXACT ANSWERS
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Select_Source
Use Interface_Header
 IMPLICIT NONE
 REAL(DP), INTENT(IN) :: X (MAX_NP, N_SPACE), DD (N_D_FRE)
 INTEGER  :: I, K, L
!          Automatic Arrays
 INTEGER  :: INDEX     (N_G_DOF)
 REAL(DP) :: EXACT_SOL (N_G_DOF)    ! PT EXACT VALUE

! MAX_NP  = NUMBER OF NODES IN SYSTEM
! N_G_DOF = NUMBER OF PARAMETERS (DOF) PER NODE
! N_D_FRE = NUMBER OF DOF IN THE SYSTEM
! N_SPACE = DIMENSION OF SPACE
! X       = SYSTEM COORDINATES OF ALL NODES
! DD      = CALCULATED NODAL PARAMETERS
! INDEX   = SYSTEM DOF NOS OF PARAMETERS ON A NODE

  WRITE (N_PRT, 50) XYZ_NAME   (1:N_SPACE), DOF_NAME (1:N_G_DOF),   &
                    EXACT_NAME (1:N_G_DOF) ; 50 FORMAT( /,          &
  '***  OUTPUT OF RESULTS AND EXACT VALUES IN NODAL ORDER  ***', /, &
  '  NODE, ', (9A12) )
  DO I = 1, MAX_NP
    CALL SELECT_EXACT_SOLUTION (X (I, :), EXACT_SOL)
    INDEX = GET_INDEX_AT_PT (I) 
    IF ( N_SPACE < 3 ) THEN
      WRITE (N_PRT, 51) I, (X (I, L), L = 1, N_SPACE), &
        (DD (INDEX (K)), K = 1, N_G_DOF), EXACT_SOL
    ELSE
      WRITE (N_PRT, 51) I, (X (I, L), L = 1, N_SPACE), &
        (DD (INDEX (K)), K = 1, N_G_DOF)
      WRITE (N_PRT, 52) I,  EXACT_SOL
    END IF
    51 FORMAT ( I6, (9(1X, ES11.4)) ) !b
    52 FORMAT ( I6, 36X, (6(1X, ES11.4)) ) !b
  END DO
END SUBROUTINE WRITE_BY_NODES_WITH_EXACT
                             
SUBROUTINE DISPLAY_D_MATRIX ( A )
! ------------------------------------------------------
! PRETTY PRINTING OF DP MATRIX A(NR,NC) WITH E FORMAT
! ------------------------------------------------------
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)   :: A(:,:)        ! THE ARRAY
 
 INTEGER, PARAMETER :: MAX_COL = 7    ! COLS PER PAGE
 INTEGER :: NCOL(MAX_COL)               ! COLUMN NUMBERS
 INTEGER :: NC, NR, COL, N, J, K, FULL
 INTEGER :: PART, COL_SIZE, PAGES

! A    = REAL ARRAY, PASSED BY COLUMNS
! COL  = MAX NUMBER OF COLUMNS PER SCREEN
! NC   = NUMBER OF COLUMNS IN A
! NCOL = VECTOR SUBSCRIPT FOR COLUMNS TO BE PRINTED
! NR   = NUMBER OF ROWS IN A

   NR  = SIZE (A, 1)         ! FIND NUMBER OF ROWS
   NC  = SIZE (A, 2)         ! FIND NUMBER OF COLUMNS
nr = SIZE (A, DIM=1)
NC  = SIZE (A, DIM=2)
if ( NR < 1 .OR. NC < 1 ) STOP 'DISPLAY_D_MATRIX: zero size'
   COL = MIN0(NC,MAX_COL)
!       LOOP OVER A SCREEN FULL OF COLUMNS
   PAGES = CEILING(REAL(NC)/COL)
   FULL = NC/COL
   PART = NC - FULL*COL
   DO J = 1, PAGES
!        FILL THE TEMPORARY COLUMN NUMBERS TO PRINT
     IF ( J > FULL )  THEN ! LAST PARTIAL PAGE
       COL_SIZE = PART
     ELSE ! FULL PAGE
       COL_SIZE = COL
     END IF ! STATUS CHECK
     !b NCOL = (/ (((J-1)*COL+K), K = 1,COL_SIZE) /)
     DO K = 1, COL_SIZE
       NCOL (K) = (J-1)*COL + K
     END DO
     WRITE (*,5)  ( NCOL(K), K = 1,COL_SIZE)
     5 FORMAT (' R\C', I6, 11I12 )
!         PRINT ROWS IN CURRENT SCREEN COLUMNS
     DO N = 1, NR
       WRITE (*,10)  N, ( A(N,NCOL(K)), K = 1,COL_SIZE )
       !b 10 FORMAT ( I3, 11(1PE12.4) )
       10 FORMAT ( I3, 11(ES12.4) ) !b
     END DO ! OF N OVER ROWS
   END DO ! OF J COLUMNS
   WRITE (*,*) ' '
END SUBROUTINE DISPLAY_D_MATRIX

SUBROUTINE DISPLAY_D_VECTOR ( A )
! ------------------------------------------------------
! PRETTY PRINTING OF DP VECTOR A(NC) WITH E FORMAT
! ------------------------------------------------------
Use Precision_Module
 IMPLICIT NONE
 REAL(DP), INTENT(IN)   :: A(:)        ! THE ARRAY
 
 INTEGER, PARAMETER :: MAX_COL = 7    ! COLS PER PAGE
 INTEGER :: NCOL(MAX_COL)               ! COLUMN NUMBERS
 INTEGER :: NC, MAX, N=1, J, K, FULL
 INTEGER :: PART, COL_SIZE, PAGES

! A    = REAL ARRAY, PASSED BY COLUMNS
! MAX  = MAX NUMBER OF COLUMNS PER SCREEN
! NC   = NUMBER OF COLUMNS IN A
! NCOL = VECTOR SUBSCRIPT FOR COLUMNS TO BE PRINTED

   NC  = SIZE (A)         ! FIND NUMBER OF ROWS
if ( NC < 1 ) STOP 'DISPLAY_D_VECTOR: zero size'
   MAX = MIN0(NC,MAX_COL)
!       LOOP OVER A SCREEN FULL OF COLUMNS
   PAGES = CEILING(REAL(NC)/MAX)
   FULL = NC/MAX
   PART = NC - FULL*MAX
   DO J = 1, PAGES
!        FILL THE TEMPORARY COLUMN NUMBERS TO PRINT
     IF ( J > FULL )  THEN ! LAST PARTIAL PAGE
       COL_SIZE = PART
     ELSE ! FULL PAGE
       COL_SIZE = MAX
     END IF ! STATUS CHECK
     !b NCOL = (/ ((J-1)*MAX+K, K = 1,COL_SIZE) /)
     DO K = 1, COL_SIZE
       NCOL (K) = (J-1)*MAX + K
     END DO
     WRITE (*,5)  ( NCOL(K), K = 1,COL_SIZE)
     5 FORMAT (' R\C', I6, 11I12 )
!         PRINT ROW IN CURRENT SCREEN COLUMNS
     WRITE (*,10)  N, ( A(NCOL(K)), K = 1,COL_SIZE )
     10 FORMAT ( I3, 11(1PE12.4) )
   END DO ! OF J COLUMNS
   WRITE (*,*) ' '
END SUBROUTINE DISPLAY_D_VECTOR

SUBROUTINE DISPLAY_R_MATRIX ( A )
! ------------------------------------------------------
! PRETTY PRINTING OF REAL MATRIX A(NR,NC) WITH E FORMAT
! ------------------------------------------------------
 IMPLICIT NONE
 REAL, INTENT(IN)   :: A(:,:)        ! THE ARRAY
 
 INTEGER, PARAMETER :: MAX_COL = 7    ! COLS PER PAGE
 INTEGER :: NCOL(MAX_COL)               ! COLUMN NUMBERS
 INTEGER :: NC, NR, MAX, N, J, K, FULL
 INTEGER :: PART, COL_SIZE, PAGES

! A    = REAL ARRAY, PASSED BY COLUMNS
! MAX  = MAX NUMBER OF COLUMNS PER SCREEN
! NC   = NUMBER OF COLUMNS IN A
! NCOL = VECTOR SUBSCRIPT FOR COLUMNS TO BE PRINTED
! NR   = NUMBER OF ROWS IN A

   NR  = SIZE (A, 1)         ! FIND NUMBER OF ROWS
   NC  = SIZE (A, 2)         ! FIND NUMBER OF COLUMNS
   MAX = MIN0(NC,MAX_COL)
!       LOOP OVER A SCREEN FULL OF COLUMNS
   PAGES = CEILING(REAL(NC)/MAX)
   FULL = NC/MAX
   PART = NC - FULL*MAX
   DO J = 1, PAGES
!        FILL THE TEMPORARY COLUMN NUMBERS TO PRINT
     IF ( J > FULL )  THEN ! LAST PARTIAL PAGE
       COL_SIZE = PART
     ELSE ! FULL PAGE
       COL_SIZE = MAX
     END IF ! STATUS CHECK
     !b NCOL = (/ ((J-1)*MAX+K, K = 1,COL_SIZE) /)
     DO K = 1, COL_SIZE
       NCOL (K) = (J-1)*MAX + K
     END DO
     WRITE (*,5)  ( NCOL(K), K = 1,COL_SIZE)
     5 FORMAT (' R\C', I6, 11I12 )
!         PRINT ROWS IN CURRENT SCREEN COLUMNS
     DO N = 1, NR
       WRITE (*,10)  N, ( A(N,NCOL(K)), K = 1,COL_SIZE )
       10 FORMAT ( I3, 11(1PE12.4) )
     END DO ! OF N OVER ROWS
   END DO ! OF J COLUMNS
   WRITE (*,*) ' '
END SUBROUTINE DISPLAY_R_MATRIX

SUBROUTINE  DISPLAY_R_VECTOR ( A )
! ------------------------------------------------------
! PRETTY PRINTING OF REAL MATRIX A(NR,NC) WITH E FORMAT
! ------------------------------------------------------
 IMPLICIT NONE
 REAL, INTENT(IN)   :: A(:)        ! THE ARRAY
 
 INTEGER, PARAMETER :: MAX_COL = 7    ! COLS PER PAGE
 INTEGER :: NCOL(MAX_COL)               ! COLUMN NUMBERS
 INTEGER :: NC, MAX, N=1, J, K, FULL
 INTEGER :: PART, COL_SIZE, PAGES

! A    = REAL ARRAY, PASSED BY COLUMNS
! MAX  = MAX NUMBER OF COLUMNS PER SCREEN
! NC   = NUMBER OF COLUMNS IN A
! NCOL = VECTOR SUBSCRIPT FOR COLUMNS TO BE PRINTED

   NC  = SIZE (A)         ! FIND NUMBER OF ROWS
   MAX = MIN0(NC,MAX_COL)
!       LOOP OVER A SCREEN FULL OF COLUMNS
   PAGES = CEILING(REAL(NC)/MAX)
   FULL = NC/MAX
   PART = NC - FULL*MAX
   DO J = 1, PAGES
!        FILL THE TEMPORARY COLUMN NUMBERS TO PRINT
     IF ( J > FULL )  THEN ! LAST PARTIAL PAGE
       COL_SIZE = PART
     ELSE ! FULL PAGE
       COL_SIZE = MAX
     END IF ! STATUS CHECK
     !b NCOL = (/ ((J-1)*MAX+K, K = 1,COL_SIZE) /)
     DO K = 1, COL_SIZE
       NCOL (K) = (J-1)*MAX + K
     END DO
     WRITE (*,5)  ( NCOL(K), K = 1,COL_SIZE)
     5 FORMAT (' R\C', I6, 11I12 )
!         PRINT ROW IN CURRENT SCREEN COLUMNS
     WRITE (*,10)  N, ( A(NCOL(K)), K = 1,COL_SIZE )
     10 FORMAT ( I3, 11(1PE12.4) )
   END DO ! OF J COLUMNS
   WRITE (*,*) ' '
END SUBROUTINE DISPLAY_R_VECTOR

SUBROUTINE DISPLAY_I_MATRIX ( A )
! ------------------------------------------------------
! PRETTY PRINTING OF INTEGER MATRIX A(NR,NC) WITH I FORMAT
! ------------------------------------------------------
 IMPLICIT NONE
 INTEGER, INTENT(IN)  :: A(:,:)        ! THE ARRAY

 INTEGER, PARAMETER   :: MAX_COL = 7    ! COLS PER PAGE
 INTEGER, ALLOCATABLE :: NCOL(:)       ! COLUMN NUMBERS
 INTEGER  NC, NR, MAX, N, J, K, FULL, PART, COL_SIZE, PAGES

! A    = INTEGER ARRAY
! MAX  = MAX NUMBER OF COLUMNS PER SCREEN
! NC   = NUMBER OF COLUMNS IN A
! NCOL = VECTOR SUBSCRIPT FOR COLUMNS TO BE PRINTED
! NR   = NUMBER OF ROWS IN A

   NR  = SIZE (A, 1)         ! FIND NUMBER OF ROWS
   NC  = SIZE (A, 2)         ! FIND NUMBER OF COLUMNS
   MAX = MIN0(NC,MAX_COL)
   ALLOCATE( NCOL(MAX) ) 
!       LOOP OVER A SCREEN FULL OF COLUMNS
   PAGES = CEILING(REAL(NC)/MAX)
   FULL = NC/MAX
   PART = NC - FULL*MAX
   DO J = 1, PAGES
!        FILL THE TEMPORARY COLUMN NUMBERS TO PRINT
     IF ( J > FULL )  THEN ! LAST PARTIAL PAGE
       COL_SIZE = PART
     ELSE ! FULL PAGE
       COL_SIZE = MAX
     END IF ! STATUS CHECK
     !b NCOL = (/ ((J-1)*MAX+K, K = 1,COL_SIZE) /)
     DO K = 1, COL_SIZE
       NCOL (K) = (J-1)*MAX + K
     END DO
     WRITE (*,5)  ( NCOL(K), K = 1,COL_SIZE)
     5 FORMAT (' R\C', I6, 11I11 )
!         PRINT ROWS IN CURRENT SCREEN COLUMNS
     DO N = 1, NR
       WRITE (*,10)  N, ( A(N,NCOL(K)), K = 1,COL_SIZE )
       10 FORMAT ( I3, 11( I11 ) )
     END DO ! OF N OVER ROWS
   END DO ! OF J COLUMNS
   WRITE (*,*) ' '
END SUBROUTINE DISPLAY_I_MATRIX

SUBROUTINE DISPLAY_I_VECTOR ( A )
! ------------------------------------------------------
! PRETTY PRINTING OF INTEGER VECTOR A(NC) WITH I FORMAT
! ------------------------------------------------------
 IMPLICIT NONE
 INTEGER, INTENT(IN)  :: A(:)        ! THE ARRAY

 INTEGER, PARAMETER   :: MAX_COL = 7    ! COLS PER PAGE
 INTEGER, ALLOCATABLE :: NCOL(:)       ! COLUMN NUMBERS
 INTEGER ::  NC, MAX, N=1, J, K, FULL, PART, COL_SIZE, PAGES

! A    = INTEGER ARRAY
! MAX  = MAX NUMBER OF COLUMNS PER SCREEN
! NC   = NUMBER OF COLUMNS IN A
! NCOL = VECTOR SUBSCRIPT FOR COLUMNS TO BE PRINTED

   NC  = SIZE (A)         ! FIND NUMBER OF ROWS
   MAX = MIN0(NC,MAX_COL)
   ALLOCATE( NCOL(MAX) ) 
!       LOOP OVER A SCREEN FULL OF COLUMNS
   PAGES = CEILING(REAL(NC)/MAX)
   FULL = NC/MAX
   PART = NC - FULL*MAX
   DO J = 1, PAGES
!        FILL THE TEMPORARY COLUMN NUMBERS TO PRINT
     IF ( J > FULL )  THEN ! LAST PARTIAL PAGE
       COL_SIZE = PART
     ELSE ! FULL PAGE
       COL_SIZE = MAX
     END IF ! STATUS CHECK
     !b NCOL = (/ ((J-1)*MAX+K, K = 1,COL_SIZE) /)
     DO K = 1, COL_SIZE
       NCOL (K) = (J-1)*MAX + K
     END DO
     WRITE (*,5)  ( NCOL(K), K = 1,COL_SIZE)
     5 FORMAT (' R\C', I6, 11I11 )
!         PRINT ROWS IN CURRENT SCREEN COLUMNS
     WRITE (*,10)  N, ( A(NCOL(K)), K = 1,COL_SIZE )
     10 FORMAT ( I3, 11( I11 ) )
   END DO ! OF J COLUMNS
   WRITE (*,*) ' '
END SUBROUTINE DISPLAY_I_VECTOR

SUBROUTINE  ELEM_BAR_CHART (I_BAR, I_PARM, G_PARM, PARM)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!   PRINT-PLOT BAR CHARTS OF ELEM PARAMETER I_PARM 
!   FROM A TOTAL OF G_PARM PER ELEMENT
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 INTEGER,  PARAMETER  :: MID = 25, LINE = 2 * MID
 INTEGER,  INTENT(IN) :: I_BAR, I_PARM, G_PARM
 REAL(DP), INTENT(IN) :: PARM (N_ELEMS * G_PARM)
!           Automatic Arrays
 INTEGER :: ELEM_BAR (I_BAR)
!            Locals
 INTEGER :: I, INDEX, IPMNG, ISPACE, JY, JZ, &
             K, KOUNT, L, LIMIT, NLAST, NOPLOT
 REAL(DP)    :: DTEST, RANGE, YMAX, YMIN
 CHARACTER ALINE (LINE+1), SKIP (LINE+1)
 CHARACTER(1), SAVE :: BLANK = ' ', DOT = '+', DASH = '-', &
                       A_X = 'X', PLUS = '+', A_0 = '0'
  DATA NOPLOT / 0 /
! I_PARM   = NODAL PARAMETER TO BE GRAPHED,1<=I_PARM<=G_PARM
! ELEM_BAR = LIST OF ELEMENTS TO BE USED (WHEN I_BAR > 1)
! PARM     = ARRAY OF ALL NODAL PARAMETERS IN THE SYSTEM
! G_PARM   = NUMBER OF PARAMETERS PER ELEMENT
! N_SPACE  = DIMENSION OF SOLUTION SPACE
! I_BAR    = NUMBER OF ELEMENTS TO BE INCLUDED IN BAR CHART
!          IF I_BAR=1  USE ALL, ELEM_BAR NOT USED
  NOPLOT = NOPLOT + 1
  LIMIT = I_BAR
  IF ( I_BAR == 1) LIMIT = N_ELEMS
  WRITE (N_PRT, 5000) NOPLOT, I_PARM, LIMIT
  5000 FORMAT ('*** PRINT PLOT NUMBER',I3,' ***',/,           &
               'ELEMENT PARAMETER', I3,                       &
               ', EVALUATED AT ',I5,' ELEMENTS',/)
  IPMNG = I_PARM - G_PARM
!-->  ESTABLISH GRAPH LIMITING VALUES
!   INITALIZE MAX MIN VALUES OF PARAMETER
  L = 1
  IF ( I_BAR /= 1) L = ELEM_BAR (1)
  INDEX = G_PARM * L + IPMNG
  YMAX = PARM (INDEX)
  YMIN = YMAX
  DO I = 1, LIMIT
    L = I
    IF ( I_BAR /= 1) L = ELEM_BAR (I)
    INDEX = G_PARM * L + IPMNG
    DTEST = PARM (INDEX)
    IF ( DTEST > YMAX) YMAX = DTEST
    IF ( DTEST < YMIN) YMIN = DTEST
  END DO
  KOUNT = 1
  WRITE (N_PRT, 5010) YMIN, YMAX
  5010 FORMAT (' RANGE ON GRAPH IS ',1PE12.5,' TO ',1PE12.5/)
  IF ( YMIN == YMAX ) THEN
    WRITE (N_PRT, 5011)
    5011 FORMAT ('CONSTANT VALUE, PLOT SKIPPED')
    RETURN
  END IF
  RANGE = YMAX - YMIN
!  FLAG ZERO VALUE
  IF ( YMIN < 0.0 .AND. YMAX > 0.0 ) THEN
    JZ = 1*(YMAX - 0.0) / RANGE                               &
       + (LINE+1)*(0.0 - YMIN) / RANGE + 0.1
  ELSE
    JZ = 0
  END IF
! SCALING COMPLETE
  SKIP  (2:LINE) = BLANK
  ALINE (2:LINE) = DASH
  SKIP (1) = PLUS
  SKIP (MID+1) = DOT
  SKIP (LINE+1) = PLUS
  ALINE (1:LINE+1:5) = PLUS
  ALINE (MID+1) = DOT
  WRITE (N_PRT, 5020) ALINE
  5020 FORMAT(' ELEMENT    VALUE   ', (101A1) )
  L = 1
  IF ( I_BAR > 1) L = ELEM_BAR (1)
  DO K = 1, LIMIT
    NLAST = L
    L = K
    IF ( I_BAR > 1) L = ELEM_BAR (K)
    INDEX = G_PARM * L + IPMNG
    DTEST = PARM (INDEX)
    ALINE (2:LINE) = BLANK
    ALINE (1) = PLUS
    ALINE (LINE+1) = PLUS
    IF ( JZ > 0) ALINE (JZ) = A_0
    ISPACE = 1
! LINEAR INTERPOLATION FOR COLUMN NUMBER
    JY = 1 * (YMAX - DTEST) / RANGE                    &
       + (LINE+1) * (DTEST - YMIN) / RANGE + 0.1
    ALINE (1:JY) = A_X
    WRITE (N_PRT, 5040) L, PARM (INDEX), ALINE
    5040 FORMAT (I7,2X,1PE10.3,1X,101A1)
  END DO
  ALINE (2:LINE) = DASH
  ALINE (1:LINE+1:5) = PLUS
  ALINE (MID+1) = DOT
  WRITE (N_PRT, 5020) ALINE
END SUBROUTINE  ELEM_BAR_CHART

! now in scp_lib SAVE_QP_FLUX_FOR_SCP (IE, IQ, XYZ, E, STRAIN)

SUBROUTINE LIST_ELEM_GRADIENTS (N_FILE, IE)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! LIST ELEMENT GRADIENTS AT QUADRATURE POINTS, ON N_FILE       
! GENERALLY USED ONLY FOR ELEM POST_PROCESSING SINCE
! FLUXES ARE NOW STANDARD OUTPUTS
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data   ! for COORD (LT_N, N_SPACE), D (LT_FREE)
 IMPLICIT NONE
!    ALWAYS USED
 INTEGER,  INTENT(IN)    :: N_FILE, IE

!         Automatic Arrays 
 REAL(DP) :: XYZ (N_SPACE), DGH (N_SPACE, LT_FREE),          &
             STRAIN (N_R_B + 2), E (N_R_B, N_R_B), RMS
 REAL(DP) :: B (N_R_B, LT_FREE)        !b

 INTEGER,  SAVE :: TEST_IE, TEST_IP, J, N_IP, EOF, IO_1
 REAL(DP), SAVE :: DERIV_MAX = 0.d0

!                 VARIABLES:
! B         = GRADIENT VERSUS DOF MATRIX (N_R_B, LT_FREE)
! COORD     = SPATIAL COORDINATES OF ELEMENT'S NODES
! D         = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT
! DGH       = GLOBAL DERIVATIVES INTERPOLATION FUNCTIONS
! E         = ELEMENT CONSTITUTIVE MATRIX AT GAUSS POINT
! G         = GEOMETRIC INTERPOLATION FUNCTIONS
! H         = SOLUTION INTERPOLATION FUNCTIONS
! LT_N      = NUMBER OF NODES PER ELEMENT
! LT_FREE   = NUMBER OF DEGREES OF FREEDOM PER ELEMENT
! N_ELEMS   = TOTAL NUMBER OF ELEMENTS                        
! N_R_B     = NUMBER OF ROWS IN B AND E MATRICES
! N_SPACE   = DIMENSION OF SPACE
! N_FILE    = UNIT FOR POST SOLUTION MATRICES STORAGE
! STRAIN    = STRAIN OR GRADIENT VECTOR
! XYZ       = SPACE COORDINATES AT A POINT
! U_PLT4    = UNIT TO STORE PLOT DATA, IF > 0           

!-->     PRINT TITLES ON THE FIRST CALL                        
  IF ( IE == 1) THEN 
    REWIND (N_FILE)
    IF ( LIST_QP_FLUX ) THEN
      WRITE (N_PRT, 5) N_SPACE, N_R_B  
      5 FORMAT (/,'*** GRADIENTS AT ELEMENT INTEGRATION POINTS ***',/, &
      ' ELEMENT, POINT, ', I2, ' COORDINATES', I2, ' GRADIENTS', /)
    END IF

!   OPEN PLOT FILE IF ACTIVE
    IF ( U_PLT4 > 0 ) OPEN (U_PLT4, FILE='el_qp_xyz_grads.tmp', &
       ACTION='WRITE', STATUS='REPLACE', IOSTAT = IO_1)
    IF ( IO_1 > 0 ) PRINT *, 'OPEN FAILED FOR el_qp_xyz_grads.tmp'
  ENDIF 

!   IS THIS AN ELEMENT, OR A BOUNDARY SEGMENT ?
  IF ( IE <= N_ELEMS ) THEN ! ELEMENT RESULTS

!      READ NUMBER OF POINTS TO BE CALCULATED                    
    READ (N_FILE, IOSTAT = EOF) N_IP 
    IF ( EOF /= 0 ) THEN
      PRINT *, 'LIST_ELEM_GRADIENTS EOF AT ELEMENT ', IE, EOF
      PRINT *, 'WAS EXPECTING: N_IP. WAS UNIT ', N_FILE, ' REWOUND ?'
      STOP     'ERROR, EOF IN LIST_ELEM_GRADIENTS'
    END IF

!-->     READ COORDS. AND DERIVATIVE MATRIX                        
    DO J = 1, N_IP 
      READ (N_FILE, IOSTAT = EOF) XYZ, E, B   !b  DGH 
      IF ( EOF /= 0 ) THEN
        PRINT *, 'LIST_ELEM_GRADIENTS EOF AT ELEMENT ', IE, EOF
        PRINT *, 'WAS EXPECTING: XYZ, E, B AT POINT ', J, ' OF ', N_IP
        !b PRINT *, 'WAS EXPECTING: XYZ, E, DGH'
        CALL INQUIRE_ABOUT_UNIT (N_FILE)
        STOP     'ERROR, EOF IN LIST_ELEM_GRADIENTS'
      END IF

!        CALCULATE DERIVATIVES,  STRAIN = B*D                   
      STRAIN (1:N_R_B) = MATMUL (B, D) 

!        SCALAR MEASURE (VALID FOR VECTORS)
      RMS = SQRT ( SUM (STRAIN (1:N_R_B)**2) )

!-->     PRINT COORDINATES AND GRADIENT AT THE POINT      
      IF ( LIST_QP_FLUX ) THEN
        WRITE (N_PRT, 20) IE, J, XYZ, STRAIN (1:N_R_B)
        20 FORMAT ( I5, I3, 10(1PE11.3), (8X, 10(1PE11.3)))
      END IF
      IF ( RMS > DERIV_MAX ) THEN
        DERIV_MAX = RMS
        TEST_IE   = IE ; TEST_IP = J
      END IF

!-->     STORE RESULTS TO BE PLOTTED LATER
      IF ( U_PLT4 > 0 .AND. IO_1 == 0 )        &
        WRITE (U_PLT4, 40) XYZ, STRAIN (1:N_R_B)
        40 FORMAT ( (10(1PE13.5)) )

    END DO 
    PRINT *, ' '

!-->   ARE GRADIENT CALCULATIONS COMPLETE FOR ALL ELEMENTS 
    IF (IE == N_ELEMS) THEN !      LIST LARGEST GRADIENT
      WRITE (N_PRT, "('LARGEST GRADIENT = ', 1PE13.5)") DERIV_MAX
      WRITE (N_PRT, "('AT ELEM =', I6, ', POINT = ', I2)") &
             TEST_IE, TEST_IP
      MAX_DERIVATIVE = DERIV_MAX ! plot scale
      IF ( U_PLT4 > 0 ) CLOSE (U_PLT4)
    END IF ! LAST ELEMENT
  END IF ! ELEMENT OR BOUNDARY SEGMENT
END SUBROUTINE LIST_ELEM_GRADIENTS

SUBROUTINE LIST_ELEM_SOLUTION_INTEGRAL (N_FILE, IE, EACH)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! LIST INTEGRAL OF SOLUTION FROM DATA POINTS, ON N_FILE       
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data    ! for  D (LT_FREE)
 IMPLICIT NONE
!    ALWAYS USED
 INTEGER,  INTENT(IN) :: N_FILE, IE
 LOGICAL,  INTENT(IN) :: EACH          ! list each if true

 REAL(DP), SAVE :: VALUE, TOTAL 
 REAL(DP)       :: H_INTG (LT_FREE)
 INTEGER        :: EOF
 
!                 VARIABLES:
! D       = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT
! H       = SOLUTION INTERPOLATION FUNCTIONS
! H_INTG  = INTEGRAL OF INTERPOLATION FUNCTIONS
! LT_FREE = NUMBER OF DEGREES OF FREEDOM PER ELEMENT
! N_FILE  = UNIT FOR POST SOLUTION MATRICES STORAGE

!-->     PRINT TITLES ON THE FIRST CALL AND INITIALIZE         
  IF ( IE == 1) THEN 
    REWIND (N_FILE)
    TOTAL = 0.d0
    IF ( EACH ) WRITE (N_PRT, 5) 
    5 FORMAT (/,'** VOLUME INTEGRAL OF DEPENDENT VARIABLE **',/, &
                'ELEMENT      INTEGRAL')
  ENDIF 
  IF ( IE <= N_ELEMS ) THEN ! ELEMENT RESULTS

!-->   CALCULATE ELEMENT CONTRIBUTION
    READ (N_FILE, IOSTAT = EOF) H_INTG
    IF ( EOF /= 0 ) THEN
      PRINT *, 'LIST_ELEM_SOLUTION_INTEGRAL EOF AT ELEMENT ', IE
      PRINT *, 'WAS EXPECTING: H_INTG'
      STOP     'ERROR, EOF IN LIST_ELEM_SOLUTION_INTEGRAL'
    END IF
    VALUE = DOT_PRODUCT (H_INTG, D)
    TOTAL = TOTAL + VALUE
    IF ( EACH ) WRITE (N_PRT, '(I7, 1PE18.6)') IE, VALUE
    IF (IE == N_ELEMS) WRITE (N_PRT, &
      "('TOTAL SOLUTION VOLUME INTEGRAL = ', 1PE16.6, /)")  TOTAL
  END IF
END SUBROUTINE LIST_ELEM_SOLUTION_INTEGRAL 

SUBROUTINE LIST_ELEM_FLUXES (N_FILE, IE)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! LIST ELEMENT FLUXES AT QUADRATURE POINTS, ON N_FILE       
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data   ! for LT_FREE, D (LT_FREE)
 IMPLICIT NONE
!    ALWAYS USED
 INTEGER,  INTENT(IN)    :: N_FILE, IE

!         Automatic Arrays 
 !b REAL(DP) :: DGH (N_SPACE, LT_FREE), STRAIN (N_R_B + 2),       &
 REAL(DP) :: STRAIN (N_R_B + 2),       &
             XYZ (N_SPACE), E (N_R_B, N_R_B), STRESS (N_R_B + 2)
 REAL(DP) :: B (N_R_B, LT_FREE), RMS   !b

 INTEGER,  SAVE :: TEST_IE, TEST_IP, J, N_IP, EOF, IO_1
 REAL(DP), SAVE :: DERIV_MAX = 0.d0

!                 VARIABLES:
! D         = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT
! E         = ELEMENT CONSTITUTIVE MATRIX AT GAUSS POINT
! G         = GEOMETRIC INTERPOLATION FUNCTIONS
! H         = SOLUTION INTERPOLATION FUNCTIONS
! LT_N      = NUMBER OF NODES PER ELEMENT
! LT_FREE   = NUMBER OF DEGREES OF FREEDOM PER ELEMENT
! N_ELEMS   = TOTAL NUMBER OF ELEMENTS                        
! N_R_B     = NUMBER OF ROWS IN B AND E MATRICES
! N_SPACE   = DIMENSION OF SPACE
! N_FILE    = UNIT FOR POST SOLUTION MATRICES STORAGE
! RMS       = ROOT MEAN SQUARE VALUE
! STRAIN    = STRAIN OR FLUX VECTOR
! XYZ       = SPACE COORDINATES AT A POINT
! U_PLT4    = UNIT TO STORE PLOT DATA, IF > 0           

!-->     PRINT TITLES ON THE FIRST CALL                        
  IF ( IE == 1) THEN 
    REWIND (N_FILE)
    IF ( LIST_QP_FLUX ) THEN
      WRITE (N_PRT, 5) XYZ_NAME (1:N_SPACE), FLUX_NAME (1:N_R_B)  !b
      5 FORMAT (/,                                               &
      '*** FLUX COMPONENTS AT ELEMENT INTEGRATION POINTS ***',/, &
      '(SAVED FOR AVERAGING OR APPLICATION POST PROCESSING)',/,  &
      'ELEMENT, PT, ', (6A12) )
    END IF
    RECORD_NUMBER = 0

!   OPEN FLUX PLOT FILE IF ACTIVE
    IF ( N_FILE5 > 0 ) OPEN (N_FILE5, FILE='el_qp_xyz_grads.tmp', &
       ACTION='WRITE', STATUS='REPLACE', IOSTAT = IO_1)
    IF ( U_PLT4 > 0 ) OPEN (U_PLT4, FILE='el_qp_xyz_fluxes.tmp', &
       ACTION='WRITE', STATUS='REPLACE', IOSTAT = IO_1)
    IF ( IO_1 > 0 ) PRINT *, 'OPEN FAILED FOR el_qp_xyz_fluxes.tmp'
  ENDIF 
!b print *,'N_FILE, IO_1 ', N_FILE, IO_1 !b
!b call LIST_SYSTEM_UNIT_NUMBERS !b

!   IS THIS AN ELEMENT, OR A BOUNDARY SEGMENT ?
  IF ( IE <= N_ELEMS ) THEN ! ELEMENT RESULTS, or IF (IS_ELEMENT)

!      READ NUMBER OF INTEGRATION POINTS TO BE EMPLOYED
    READ (N_FILE, IOSTAT = EOF) N_IP  ! expect N_FILE == U_FLUX
    IF ( EOF /= 0 ) THEN ! HIT END OF FILE
      PRINT *, 'LIST_ELEM_FLUXES EOF AT ELEMENT ', IE, EOF
      PRINT *, 'WAS EXPECTING: N_IP. WAS UNIT ', N_FILE, ' REWOUND ?'
      CALL INQUIRE_ABOUT_UNIT (N_FILE)
      STOP     'ERROR, EOF IN LIST_ELEM_FLUXES'
    END IF
    IF ( DEBUG_SCP ) THEN
     PRINT *,'ELEM, DOF ', IE, D
    END IF

!-->     READ COORDS, CONSTITUTIVE,  AND DERIVATIVE MATRIX             
    DO J = 1, N_IP  ! OVER ALL INTEGRATION POINTS
      READ (N_FILE, IOSTAT = EOF) XYZ, E, B  ! expect N_FILE == U_FLUX
      IF ( EOF /= 0 ) THEN
        PRINT *, 'LIST_ELEM_FLUXES EOF AT ELEMENT ', IE, EOF
        PRINT *, 'WAS EXPECTING: XYZ, E, B AT POINT ', J, ' OF ', N_IP
        PRINT *, 'AT UNIT ', N_FILE
        CALL INQUIRE_ABOUT_UNIT (N_FILE)
        STOP     'ERROR, EOF IN LIST_ELEM_FLUXES'
      END IF

!        CALCULATE DERIVATIVES,  STRAIN = B * D                   
      STRAIN (1:N_R_B) = MATMUL (B, D)  !b 
      IF ( DEBUG_SCP ) THEN
       PRINT *,'J, XYZ  ', J, XYZ
       PRINT *,'B MATRIX'; CALL RPRINT (B, N_R_B, LT_FREE, 1)
       PRINT *,'E MATRIX'; CALL RPRINT (E, N_R_B, N_R_B,   1)
       PRINT *,'STRAINS '; CALL RPRINT (STRAIN(1:N_R_B), 1, N_R_B, 1)
      END IF

!        FLUX FROM CONSTITUTIVE DATA
      STRESS (1:N_R_B) = MATMUL (E, STRAIN (1:N_R_B))
      !bb IF ( GRAD_BASE_ERROR ) STRESS (1:N_R_B) = STRAIN (1:N_R_B)
      !b IF ( GRAD_BASE_ERROR ) THEN
        !b STRESS (1:N_R_B) = STRAIN (1:N_R_B)
      !b ELSE
        !b STRESS (1:N_R_B) = MATMUL (E, STRAIN (1:N_R_B))
      !b END IF ! GRAD VS FLUX

!        SCALAR MEASURE (VALID FOR VECTORS)
      RMS = SQRT ( SUM (STRESS (1:N_R_B)**2) )

!-->     PRINT COORDINATES AND FLUX AT THE POINT      
      IF ( LIST_QP_FLUX ) WRITE (N_PRT, '(I7, I3, 10(ES12.4))') &
           IE, J, XYZ, STRESS (1:N_R_B)
      IF ( RMS > DERIV_MAX ) THEN
        DERIV_MAX = RMS                 ! LARGEST VALUE
        TEST_IE   = IE ; TEST_IP = J    ! AT ELEM AND POINT
      END IF

!-->     STORE FLUX RESULTS TO BE PLOTTED LATER, IF USED
!b print *,'U_PLT4, IO_1 ', U_PLT4, IO_1 !b
      IF ( U_PLT4 > 0 .AND. IO_1 == 0 )        &
        WRITE (U_PLT4, '( (10(1PE13.5)) )') XYZ, STRESS (1:N_R_B)
      IF ( N_FILE5 > 0 .AND. IO_1 == 0 )        &
        WRITE (N_FILE5, '( (10(1PE13.5)) )') XYZ, STRAIN (1:N_R_B)

!        SAVE COORDINATES & FLUX COMPONENTS FOR SCP FLUX AVERAGING
      IF ( U_SCPR > 0 ) THEN
        RECORD_NUMBER             = RECORD_NUMBER + 1
        SCP_RECORD_NUMBER (IE, J) = RECORD_NUMBER
        !b WRITE (U_SCPR, REC = RECORD_NUMBER) XYZ, STRESS (1:SCP_FIT)
        IF ( GRAD_BASE_ERROR ) THEN !b JEA 4/6/02 XXX
          WRITE (U_SCPR, REC = RECORD_NUMBER) XYZ, STRAIN (1:SCP_FIT)
        ELSE
          WRITE (U_SCPR, REC = RECORD_NUMBER) XYZ, STRESS (1:SCP_FIT)
        END IF ! GRAD VS FLUX
!       NOTE: THE ERROR ESTIMATOR WILL NEED E ALSO
     !b WRITE (U_SCPR, REC = RECORD_NUMBER) XYZ, E, STRESS (1:SCP_FIT)
     !b change SCP to read E, and change record length !b
      END IF ! SCP RECOVERY
    END DO  ! OVER INTEGRATION POINTS

!b print *,'rec numb ', RECORD_NUMBER !b
    PRINT *, ' '
    !b IF ( I_BUG > 0 ) CALL INQUIRE_ABOUT_SCP_UNIT !b

!-->   ARE FLUX CALCULATIONS COMPLETE FOR ALL ELEMENTS 
    IF ( IE == N_ELEMS ) THEN !      LIST LARGEST FLUX
      WRITE (N_PRT, "('LARGEST FLUX RMS_VALUE = ', 1PE13.5)") &
             DERIV_MAX
      WRITE (N_PRT, "('AT ELEM =', I6, ', POINT = ', I2)")   &
             TEST_IE, TEST_IP
      MAX_DERIVATIVE = DERIV_MAX ! FOR POSSIBLE PLOT SCALE
      IF ( U_PLT4 > 0 ) CLOSE (U_PLT4)
      !b IF ( U_PLT4 > 0 ) THEN
        !b CALL INQUIRE_ABOUT_UNIT (U_PLT4 )
        !b PRINT *,'el_qp_xyz_fluxes CLOSED AT ELEMENT ', IE
        !b CLOSE (U_PLT4)
      !b END IF
    END IF ! LAST ELEMENT

  END IF ! ELEMENT OR BOUNDARY SEGMENT

!  FLAG IF ANY SCP DATA WERE SAVED FOR AVERAGING   
  CALL UPDATE_SCP_STATUS
  IF ( IE == N_ELEMS .AND. IO_1 == 0 ) PRINT *,         &
     'NOTE: SAVED ELEMENT FLUXES IN el_qp_xyz_fluxes.tmp'
END SUBROUTINE LIST_ELEM_FLUXES

SUBROUTINE LIST_ELEM_AND_EXACT_FLUXES (N_FILE, IE)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! LIST ELEMENT FLUXES AT QUADRATURE POINTS, ON N_FILE       
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Select_Source
Use Elem_Type_Data   ! for COORD (LT_N, N_SPACE), D (LT_FREE)
 IMPLICIT NONE
!    ALWAYS USED
 INTEGER,  INTENT(IN)    :: N_FILE, IE

!         Automatic Arrays 
 !b REAL(DP) :: DGH (N_SPACE, LT_N), STRAIN (N_R_B + 2),       & !b
 REAL(DP) :: STRAIN (N_R_B + 2),       & !b
             XYZ (N_SPACE), E (N_R_B, N_R_B), STRESS (N_R_B + 2)
 REAL(DP) :: B (N_R_B, LT_FREE), FLUX (N_R_B), RMS   !b

 INTEGER,  SAVE :: TEST_IE, TEST_IP, J, N_IP, EOF, IO_1
 REAL(DP), SAVE :: DERIV_MAX = 0.d0

!                 VARIABLES:
! COORD     = SPATIAL COORDINATES OF ELEMENT'S NODES
! D         = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT
! E         = ELEMENT CONSTITUTIVE MATRIX AT GAUSS POINT
! FLUX      = EXACT FLUX COMPONENTS AT THE POINT
! G         = GEOMETRIC INTERPOLATION FUNCTIONS
! H         = SOLUTION INTERPOLATION FUNCTIONS
! LT_N      = NUMBER OF NODES PER ELEMENT
! LT_FREE   = NUMBER OF DEGREES OF FREEDOM PER ELEMENT
! N_ELEMS   = TOTAL NUMBER OF ELEMENTS                        
! N_R_B     = NUMBER OF ROWS IN B AND E MATRICES
! N_SPACE   = DIMENSION OF SPACE
! N_FILE    = UNIT FOR POST SOLUTION MATRICES STORAGE
! RMS       = ROOT MEAN SQUARE VALUE
! STRAIN    = STRAIN OR FLUX VECTOR
! XYZ       = SPACE COORDINATES AT A POINT
! U_PLT4    = UNIT TO STORE PLOT DATA, IF > 0           

!-->     PRINT TITLES ON THE FIRST CALL                        
  IF ( IE == 1) THEN 
    FLUX = 0.d0
    REWIND (N_FILE)
    WRITE (N_PRT, 4) ; 4 FORMAT (/,                                 &
    '*** FE AND EXACT FLUX COMPONENTS AT INTEGRATION POINTS ***',/, &
    '(SAVED FOR AVERAGING OR APPLICATION POST PROCESSING)' )
    RECORD_NUMBER = 0

!   OPEN FLUX PLOT FILE IF ACTIVE
    IF ( U_PLT4 > 0 ) OPEN (U_PLT4, FILE='el_qp_xyz_fluxes.tmp', &
       ACTION='WRITE', STATUS='REPLACE', IOSTAT = IO_1)
    IF ( IO_1 > 0 ) THEN
      PRINT *, 'WARNING: OPEN FAILED FOR el_qp_xyz_fluxes.tmp'
      N_WARN = N_WARN + 1
    END IF 
  ENDIF 

!   IS THIS AN ELEMENT, OR A BOUNDARY SEGMENT ?
  IF ( IE <= N_ELEMS ) THEN ! ELEMENT RESULTS

    IF ( LIST_QP_FLUX ) THEN
      WRITE (N_PRT, 5) XYZ_SHORT (1:N_SPACE), FLUX_SHORT (1:N_R_B)  !b
      WRITE (N_PRT, 6) XYZ_SHORT (1:N_SPACE), EXACT_SHORT (1:N_R_B) !b
      5 FORMAT ('ELEMENT, PT, ', (9A10) )
      6 FORMAT ('         PT, ', (9A10) )
    END IF

!      READ NUMBER OF INTEGRATION POINTS TO BE EMPLOYED
    READ (N_FILE, IOSTAT = EOF) N_IP 
    IF ( EOF /= 0 ) THEN ! HIT END OF FILE
      PRINT *, 'LIST_ELEM_AND_EXACT_FLUXES EOF AT ELEMENT ', IE, EOF
      PRINT *, 'WAS EXPECTING: N_IP. WAS UNIT ', N_FILE, ' REWOUND ?'
      PRINT *, 'MAY NEED KEYWORD no_scp_ave IN CONTROL, OR MAY NEED'
      PRINT *, 'TO DELETE KEYWORD list_exact_flux, OR MAY NEED TO HAVE'
      PRINT *, 'STORE_FLUX_POINT_COUNT CALL IN my_el_sq_inc'
      CALL INQUIRE_ABOUT_UNIT (N_FILE)
      STOP     'ERROR, EOF IN LIST_ELEM_AND_EXACT_FLUXES'
    END IF

!-->     READ COORDS, CONSTITUTIVE, AND DERIVATIVE MATRIX             
    DO J = 1, N_IP  ! OVER ALL INTEGRATION POINTS
      READ (N_FILE, IOSTAT = EOF) XYZ, E, B   
      IF ( EOF /= 0 ) THEN
        PRINT *, 'LIST_ELEM_AND_EXACT_FLUXES EOF AT ELEMENT ', IE, EOF
        PRINT *, 'WAS EXPECTING: XYZ, E, B AT POINT ', J, ' OF ', N_IP
        PRINT *, 'AT UNIT ', N_FILE
        PRINT *, 'MAY NEED KEYWORD no_scp_ave IN CONTROL, OR MAY NEED'
        PRINT *, 'TO DELETE KEYWORD list_exact_flux, OR MAY NEED'
        PRINT *, 'STORE_FLUX_POINT_DATA CALL IN my_el_sq_inc'
        CALL INQUIRE_ABOUT_UNIT (N_FILE)
        STOP 'ERROR, EOF IN LIST_ELEM_AND_EXACT_FLUXES'
      END IF

!        CALCULATE DERIVATIVES,  STRAIN = B * D                   
      STRAIN (1:N_R_B) = MATMUL (B, D)  !b 

!        FLUX FROM CONSTITUTIVE DATA
      STRESS (1:N_R_B) = MATMUL (E, STRAIN (1:N_R_B))
if ( DEBUG_POST_EL .AND. EXAMPLE == 202 .AND. THIS_EL < 3 ) then
 print *, 'E ', THIS_EL; CALL RPRINT (E, N_R_B, N_R_B, 1)
 print *, 'STRAIN '  ; CALL RPRINT (STRAIN, 1, N_R_B, 1)
 print *, 'STRESS '  ; CALL RPRINT (STRESS, 1, N_R_B, 1)
end if

!        SCALAR MEASURE (VALID FOR VECTORS)
      RMS = SQRT ( SUM (STRESS (1:N_R_B)**2) )

!        GET EXACT FLUX VALUES
      CALL SELECT_EXACT_FLUX (XYZ, FLUX)
      
      IF ( DEBUG_B .AND. IE == 1 ) THEN
        PRINT *,'J ', J
        PRINT *,'E ' ; CALL RPRINT (E, N_R_B, N_R_B,   1)
        PRINT *,'B ' ; CALL RPRINT (B, N_R_B, LT_FREE, 1)
      END IF ! debug

!-->     PRINT COORDINATES AND FLUX AT THE POINT      
      IF ( LIST_QP_FLUX ) THEN 
        WRITE (N_PRT, '(I8, I3, 10(1PE11.3), (8X, 10(1PE11.3)))')    &
               IE, J, XYZ,  STRESS (1:N_R_B)
        WRITE (N_PRT, '(8X, I3, 10(1PE11.3), (8X, 10(1PE11.3)))')    &
               J, XYZ,  FLUX
      END IF
      IF ( RMS > DERIV_MAX ) THEN
        DERIV_MAX = RMS                 ! LARGEST VALUE
        TEST_IE   = IE ; TEST_IP = J    ! AT ELEM AND POINT
      END IF

!-->     STORE FLUX RESULTS TO BE PLOTTED LATER, IF USED
      IF ( U_PLT4 > 0 .AND. IO_1 == 0 )        &
        WRITE (U_PLT4, '( (10(1PE13.5)) )') XYZ, STRESS (1:N_R_B)

!        SAVE COORDINATES & FLUX COMPONENTS FOR SCP FLUX AVERAGING
      IF ( U_SCPR > 0 ) THEN
        RECORD_NUMBER             = RECORD_NUMBER + 1
        SCP_RECORD_NUMBER (IE, J) = RECORD_NUMBER
        IF (RECORD_NUMBER == 1 .AND. DEBUG_SCP)         &
          PRINT *,'Debug N_R_B, SCP_FIT ', N_R_B, SCP_FIT
        IF ( GRAD_BASE_ERROR ) THEN !b JEA 4/6/02 XXX
          WRITE (U_SCPR, REC = RECORD_NUMBER) XYZ, STRAIN (1:SCP_FIT)
        ELSE
          WRITE (U_SCPR, REC = RECORD_NUMBER) XYZ, STRESS (1:N_R_B) !b
          !b WRITE (U_SCPR, REC = RECORD_NUMBER) XYZ, STRESS (1:SCP_FIT)
        END IF ! GRAD VS FLUX
      END IF ! SCP RECOVERY

    END DO  ! OVER INTEGRATION POINTS

!b print *,'rec numb ', RECORD_NUMBER !b
    !b PRINT *, ' '
    !b IF ( DEBUG_SCP ) CALL INQUIRE_ABOUT_SCP_UNIT !b

!-->   ARE FLUX CALCULATIONS COMPLETE FOR ALL ELEMENTS 
    IF ( IE == N_ELEMS ) THEN !      LIST LARGEST FLUX
      WRITE (N_PRT, "('LARGEST FLUX RMS_VALUE = ', 1PE12.4)") &
             DERIV_MAX
      WRITE (N_PRT, "('AT ELEM =', I6, ', POINT = ', I2)")   &
             TEST_IE, TEST_IP
      MAX_DERIVATIVE = DERIV_MAX ! FOR POSSIBLE PLOT SCALE
      IF ( U_PLT4 > 0 ) CLOSE (U_PLT4)
    END IF ! LAST ELEMENT

  END IF ! ELEMENT OR BOUNDARY SEGMENT

!  FLAG IF ANY SCP DATA WERE SAVED FOR AVERAGING   
  CALL UPDATE_SCP_STATUS
  IF ( IE == N_ELEMS .AND. IO_1 == 0 ) PRINT *,         &
     'NOTE: SAVED ELEMENT FLUXES (ONLY) IN el_qp_xyz_fluxes.tmp'

END SUBROUTINE LIST_ELEM_AND_EXACT_FLUXES

SUBROUTINE LIST_ELEM_CONVECTION_LOSS (N_FILE, IE, EACH)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! LIST CONVECTION HEAT LOSS IN EACH ELEMENT, ASSUMING
! CONSTANT CONVECTION COEFFICIENT AND EXTERNAL TEMPERATURE
!   LOSS = integral ( h * H * (T - T_e))
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data    ! for  D (LT_FREE)
 IMPLICIT NONE

 INTEGER,  INTENT(IN) :: N_FILE, IE
 LOGICAL,  INTENT(IN) :: EACH          ! list each if true

 REAL(DP), SAVE :: VALUE, TOTAL, T_E
 REAL(DP)       :: H_INTG (LT_FREE)    ! integral h * H
 INTEGER        :: EOF
 
! D       = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT
! EACH    = .TRUE. TO LIST EVERY VALUE, ELSE GIVE ONLY TOTAL
! H       = SOLUTION INTERPOLATION FUNCTIONS
! H_INTG  = INTEGRAL OF INTERPOLATION FUNCTIONS
! IE      = ELEMENT, OR FACE NUMBER
! LT_FREE = NUMBER OF DEGREES OF FREEDOM PER ELEMENT
! N_FILE  = UNIT FOR POST SOLUTION MATRICES STORAGE
! T_E     = EXTERNAL TEMPERATURE FOR THE ELEMENT

!-->     PRINT TITLES ON THE FIRST CALL AND INITIALIZE         
  IF ( IE == 1) THEN 
    REWIND (N_FILE)
    TOTAL = 0.d0
    IF ( EACH ) WRITE (N_PRT, 5) 
    5 FORMAT (/,'** CONVECTIVE HEAT LOSS RESULTS **',/, &
                'ELEMENT      LOSS')
  ENDIF 
!b print *, 'N_FILE ', N_FILE  !b
!b call LIST_SYSTEM_UNIT_NUMBERS !b
  IF ( IE <= N_ELEMS ) THEN ! ELEMENT RESULTS

!-->   GET ELEMENT EXTERNAL TEMPERATURE & CONVECTION INTEGRAL
    READ (N_FILE, IOSTAT = EOF) T_E, H_INTG
    IF ( EOF /= 0 ) THEN
      PRINT *, 'LIST_ELEM_CONVECTION_LOSS EOF AT ELEMENT ', IE
      PRINT *, 'WAS EXPECTING: T_E, H_INTG'
      STOP     'ERROR, EOF IN LIST_ELEM_CONVECTION_LOSS'
    END IF

!-->   CALCULATE ELEMENT CONTRIBUTION
    VALUE = DOT_PRODUCT (H_INTG, (D - T_E))
    TOTAL = TOTAL + VALUE
    IF ( EACH ) WRITE (N_PRT, '(I7, 1PE18.6)') IE, VALUE
    IF (IE == N_ELEMS) WRITE (N_PRT, &
      "(/, 'TOTAL CONVECTIVE HEAT LOSS     = ', 1PE16.6, /)") TOTAL
  END IF
END SUBROUTINE LIST_ELEM_CONVECTION_LOSS 

SUBROUTINE LIST_REACTIONS_AT_ELEMS (NODES, DD)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!       ELEMENT LEVEL REACTION RECOVERY
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data
Use Sys_Properties_Data  ! for DELETED
Use Interface_Header
 IMPLICIT NONE
! Material option for future uses. Add to arguments, etc.
!    ALWAYS USED
 REAL(DP), INTENT(IN) :: DD (N_D_FRE)
 INTEGER,  INTENT(IN) :: NODES (L_S_TOT, NOD_PER_EL)
!    Automatic Arrays
 INTEGER  :: IE, LT

!                 VARIABLES:
! C            = ELEMENT INTERNAL SOURCE VECTOR
! D            = NODAL PARAMETERS ASSOCIATED WITH AN ELEMENT
! DD           = ARRAY OF SYSTEM DEGREES OF FREEDOM
! INDEX        = SYSTEM DOF NOS ASSOCIATED WITH ELEMENT
! ELEM_NODES   = THE NOD_PER_EL INCIDENCES OF THE ELEMENT
! NOD_PER_EL   = NUMBER OF NODES PER ELEMENT
! N_D_FRE      = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM
! N_ELEMS      = NUMBER OF ELEMENTS IN SYSTEM
! N_G_DOF      = NUMBERS PARAMETERS PER NODE
! NODES        = ELEMENT INCIDENCES OF ALL ELEMENTS
! L_REACT      = UNIT FOR POST SOLUTION MATRICES STORAGE
! S            = ELEMENT SQUARE MATRIX
                                        
!  L_REACT MUST BE > 0
  IF ( L_REACT ==  0 ) THEN
    STOP 'NO L_REACT IN LIST_REACTIONS_AT_ELEMS'
  ELSE
    IF ( L_REACT > 0 ) REWIND (L_REACT)
  END IF
  LT = 1 ! INITIALIZE ELEMENT TYPES 
 
!-->   LOOP OVER ELEMENTS
  DO IE = 1, N_ELEMS ! for elements 

!-->     GET ELEMENT TYPE NUMBER
    IF ( N_L_TYPE > 1) LT = L_TYPE (IE)
!    SAME AS LAST TYPE ?
    IF ( LT /= LAST_LT ) THEN ! this is a new type
!      GET CONTROLS FOR THIS TYPE
      CALL SET_ELEM_TYPE_INFO (LT)
    END IF ! a new element type

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

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

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

!-->     EXTRACT NODAL PARAMETERS OF THE ELEMENT
      !b CALL EL_FRE (N_D_FRE, LT_FREE, D, DD, INDEX)
      !b CALL GET_ELEM_DOF (D, DD, INDEX)
      D = GET_ELEM_DOF (DD)

!--> LIST THE REACTIONS, ELEMENT SOURCES, AND THEIR SUMS
      !b CALL GET_EL_REACTIONS (IE, LT_N, LT_FREE, N_G_DOF, &
                               !b L_REACT, ELEM_NODES, D)
      IF ( LIST_LR > 0 .AND. L_REACT > 0 )                   &
          CALL GET_EL_REACTIONS (IE, LT_N, LT_FREE, N_G_DOF, & !b
                                 L_REACT, ELEM_NODES, D) !b
    END DO ! over all elements and boundary segments

!--> COMPLETED, FREE TYPE ARRAYS
  !b CALL DEALLOCATE_TYPE_INTERPOLATIONS !bb

contains

SUBROUTINE GET_EL_REACTIONS (IE, LT_N, LT_FREE, N_G_DOF, &
                               N_FILE, ELEM_NODES, D)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!         GET REACTIONS (FLUXES) AT AN ELEMENT'S NODES,  
!                            R = S*D - C 
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 IMPLICIT NONE
 INTEGER,  INTENT(IN) :: IE, LT_N, LT_FREE, N_G_DOF
 INTEGER,  INTENT(IN) :: N_FILE
 INTEGER,  INTENT(IN) :: ELEM_NODES (LT_N)
 REAL(DP), INTENT(IN) :: D (LT_FREE)
                                     
 INTEGER  :: IG, IN, IROW, EOF
 REAL(DP) :: S (LT_FREE, LT_FREE), C (LT_FREE)
 REAL(DP) :: ELEM_REACTS (LT_FREE)
 REAL(DP) :: SUM_REACT   (N_G_DOF)
 REAL(DP) :: SUM_SOURCE  (N_G_DOF)
           
! C          = ELEMENT COLUMN MATRIX
! D          = KNOWN SOLUTION CAUSING THE REACTIONS
! IE         = CURRENT ELEMENT NUMBER
! ELEM_NODES = ELEMENT TOPOLOGY LIST
! ELEM_REACTS = ELEMENT REACTIONS VECTOR
! LT_N       = NUMBER OF NODES PER ELEMENT
! LT_FREE    = NUMBER OF ELEMENT DOF, N_G_DOF*N
! N_G_DOF    = NUMBER OF GENERALIZED UNKNOWNS PER NODE
! N_FILE     = UNIT TO HOLDING S & C FROM ELEMENT IE
! S          = ELEMENT SQUARE MATRIX
! SUM_SOURCE     = ELEMENT WORKSPACE RESULT
! SUM_REACT     = NODAL WORKSPACE RESULT
                                 
!  READ STIFFNESS, AND SOURCE FOR ELEM REACTIONS
  IF ( N_FILE < 1)  STOP 'INVALID UNIT IN GET_EL_REACTIONS'
  IF ( IE == 1 )  WRITE (N_PRT, &
    "(/, '** ELEMENT REACTION, INTERNAL SOURCES AND SUMMATIONS **')")
  READ (N_FILE, IOSTAT = EOF) S, C
  IF ( EOF /= 0 ) THEN
    PRINT *, 'GET_EL_REACTIONS: EOF AT ELEMENT ', IE
    PRINT *, 'WAS EXPECTING: S, C'
    STOP     'ERROR, EOF IN GET_EL_REACTIONS'
  END IF
  WRITE (N_PRT, * ) '               ELEMENT ', IE
  WRITE (N_PRT, * ) ' NODE  DOF   REACTION       ELEM_SOURCE    SUMS'

!   GIVE REACTIONS, R = S*D - C 
  ELEM_REACTS = MATMUL (S, D) - C 
  SUM_REACT = 0. ; SUM_SOURCE = 0.
  DO IN = 1, LT_N
    DO IG = 1, N_G_DOF
      IROW            = N_G_DOF * (IN - 1) + IG
      SUM_SOURCE (IG) = SUM_SOURCE (IG) + C (IROW)
      SUM_REACT  (IG) = SUM_REACT  (IG) + ELEM_REACTS (IROW) 
      WRITE (N_PRT, "( I5, I5, 1PE15.5, 1PE15.5 )") &
               ELEM_NODES (IN), IG, ELEM_REACTS (IROW), C (IROW)
    END DO
  END DO
  DO IG = 1, N_G_DOF
    PRINT "(' SUM:', I5, 3(1PE15.5))", IG, SUM_REACT (IG), &
    SUM_SOURCE (IG), (SUM_REACT (IG) + SUM_SOURCE (IG))
  END DO
 END SUBROUTINE GET_EL_REACTIONS
END SUBROUTINE LIST_REACTIONS_AT_ELEMS

!b see library file_utility_lib.f for INQUIRE_ABOUT_UNIT
!b see library file_utility_lib.f for INQUIRE_ABOUT_FILE 

SUBROUTINE LIST_SYSTEM_UNIT_NUMBERS
Use System_Constants
  PRINT *, ' '
  PRINT *, '** CURRENT SYSTEM UNIT NUMBER VALUES **  (0 IF UNUSED)'
  PRINT *, ' '
  PRINT *, 'BINARY USER ON/OFF CONTROL:'
  PRINT *, 'USER APPLICATION OPTION ........ N_FILE1 = ', N_FILE1
  PRINT *, 'USER APPLICATION OPTION ........ N_FILE2 = ', N_FILE2
  PRINT *, 'USER APPLICATION OPTION ........ N_FILE3 = ', N_FILE3
  PRINT *, 'USER APPLICATION OPTION ........ N_FILE4 = ', N_FILE4
  PRINT *, 'USER APPLICATION OPTION ........ N_FILE5 = ', N_FILE5
  PRINT *, ' '
  PRINT *, 'FORMATTED I/O DEFAULTS:'
  PRINT *, 'STANDARD INPUT ................. N_CRD   = ', N_CRD 
  PRINT *, 'STANDARD OUTPUT ................ N_PRT   = ', N_PRT 
  PRINT *, 'DEBUG OUTPUT ................... N_BUG   = ', N_BUG 
  PRINT *, ' '
  PRINT *, 'FORMATTED RESULTS (AND MATLAB PLOTS):'
  PRINT *, 'RESULT dof_history_###.tmp ..... U_DOF   = ', U_DOF 
  PRINT *, 'RESULT node_history_###.tmp .... U_NODE  = ', U_NODE
  PRINT *, 'RESULT msh_bc_xyz.tmp .......... U_PLT1  = ', U_PLT1
  PRINT *, 'RESULT msh_typ_nodes.tmp ....... U_PLT2  = ', U_PLT2
  PRINT *, 'RESULT node_results.tmp ........ U_PLT3  = ', U_PLT3
  PRINT *, 'RESULT el_qp_xyz_fluxes.tmp .... U_PLT4  = ', U_PLT4
  PRINT *, 'RESULT scp_node_ave_fluxes.tmp . U_PLT5  = ', U_PLT5
  PRINT *, 'RESULT el_error_est.tmp ........ U_PLT6  = ', U_PLT6
  PRINT *, 'RESULT pt_ave_error_est.tmp .... U_PLT7  = ', U_PLT7
  PRINT *, 'RESULT pt_ave_grad_flux.tmp .... U2_PLT1 = ', U2_PLT1
  PRINT *, 'EXACT  pt_ex_grad_flux.tmp ....  U2_PLT2 = ', U2_PLT2
  PRINT *, 'EXACT exact_node_solution.tmp .. U_PLT8  = ', U_PLT8
  PRINT *, 'EXACT exact_node_flux.tmp ...... U_PLT9  = ', U_PLT9
  PRINT *, 'Note: See CALL TURN_OFF_PLOTS, and keywords sav_*'
  PRINT *, ' '
  PRINT *, 'BINARY INTERNAL DEFAULTS:'
  PRINT *, 'SYSTEM REACTION STORAGE ........ N_REACT = ', N_REACT
  PRINT *, 'ELEMENTS REACTION STORAGE ...... L_REACT = ', L_REACT
  PRINT *, 'SUPER_CONVERGENT_PATCH ......... U_SCPR  = ', U_SCPR
  PRINT *, 'SCP 2ND DERIVATIVES ............ U_2_DER = ', U_2_DER
  PRINT *, 'SCP 3RD DERIVATIVES ............ U_3_DER = ', U_3_DER
  PRINT *, 'Note: See CALL TURN_ON_SCP, CALL TURN_OFF_SCP, no_scp_ave'
  PRINT *, 'ELEMENT FLUX DATA .............. U_FLUX  = ', U_FLUX
  PRINT *, 'ELEMENT INTEGRATIONS ........... U_INTGR = ', U_INTGR
  PRINT *, 'Note: CALL INQUIRE_ABOUT_UNIT (unit) for current status'   
  PRINT *, 'Note: CALL LIST_OPEN_IO_UNITS for current open status'   
  PRINT *, ' '
  PRINT *, 'FORMATED JOB LOG, AND HISTORIES'
  PRINT *, 'CURRENT JOB LOG ................ U_LOG   = ', U_LOG
  PRINT *, 'HISTORY OF RECENT JOB LOGS ..... U_HIST  = ', U_HIST
  PRINT *, 'keyword_input.echo ............ECHO_UNIT = ', ECHO_UNIT
END SUBROUTINE LIST_SYSTEM_UNIT_NUMBERS

SUBROUTINE READ_STATUS (IO_STATUS, ITEM, SUB_ITEM)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! IF ERROR, LIST IO_STATUS, ITEM AND SUB_ITEM NUMBERS
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants ! For N_WARN, LINE_WARN 
IMPLICIT NONE
 INTEGER, INTENT (IN) :: IO_STATUS, ITEM, SUB_ITEM
 !b INTEGER, PARAMETER   :: LIMIT = 5
 INTEGER, SAVE        :: COUNT = 0

   SELECT CASE ( IO_STATUS )  
     CASE (:-1) ! END OF FILE
       CALL CITE_INPUT_LINE ! List current input position
       PRINT *, 'INPUT ERROR: END OF FILE AT ITEM ', ITEM
       IF (SUB_ITEM > 0 ) PRINT *, 'AND SUB_ITEM ', SUB_ITEM
       STOP 'INPUT ERROR: END OF FILE'
     CASE (1:)  ! PROBABLE USER FORMAT ERROR
       CALL CITE_INPUT_LINE ! List current input position
       PRINT *, 'WARNING: USER FORMAT ERROR LIKELY AT ITEM ', ITEM
       IF (SUB_ITEM > 0 ) PRINT *, 'AND SUB_ITEM ', SUB_ITEM
       N_WARN = N_WARN + 1
       COUNT  = COUNT + 1
     CASE DEFAULT ; RETURN  ! NO ERROR
   END SELECT
   !b IF ( COUNT == LIMIT ) STOP 'READ_STATUS: MAX WARNINGS REACHED'
   IF ( COUNT == LINE_WARN ) THEN
     PRINT *,'ERROR: MAX WARNINGS FOR 1 LINE REACHED, ENTER keyword'
     PRINT *,'LINE_WARN WITH A VALUE > ', LINE_WARN
     STOP 'READ_STATUS: MAX WARNINGS REACHED'
   END IF ! avoid infinite loops
END SUBROUTINE READ_STATUS 

SUBROUTINE LIST_2ND_DERIV_ORDER
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! List the storage order used for second derivatives in MODEL
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants

 PRINT *,'NOTE, EXACT AND FE 2ND DERIVATIVES ARE STORED HERE AS:'
 SELECT CASE (N_SPACE)
   CASE (1) ; PRINT *,'EXACT: U,xx    FE: U,xx'
   CASE (2) ; PRINT *,'EXACT 3: U,xx U,yx=U,xy U,yy'
              PRINT *,'FE    4: U,xx U,yx U,xy U,yy'
   CASE (3)
     PRINT *,'EXACT 6: U,xx U,yx=U,xy U,zx=U,xz U,yy U,zy=U,yz U,zz'
     PRINT *,'FE    9: U,xx U,yx U,zx U,xy U,yy U,zy U,xz U,yz U,zz'
   CASE DEFAULT
     PRINT *,'ARE NOT DEFINED IN THIS DIMENSION SPACE'
 END SELECT
END SUBROUTINE LIST_2ND_DERIV_ORDER

SUBROUTINE GET_STEP_RESULTS (DD) 
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!      GET NODAL RESULTS DATA FOR RESTART
!        Requires Global RESTART_STEP >= 0
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Interface_Header
 IMPLICIT NONE
 CHARACTER (LEN=3)  :: T_STRING  ! for file names
 CHARACTER (LEN=20) :: F_NAME    ! file name
 REAL(DP), INTENT(OUT) :: DD (N_D_FRE)
 INTEGER :: I, K, IO_1, U_TMP, GET_NEXT_IO_UNIT
!          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_TMP  = UNIT TO OUTPUT NODE RESULTS FOR PLOTS

  U_TMP = GET_NEXT_IO_UNIT ()
  T_STRING = TO3STR (RESTART_STEP)
  F_NAME = 'node_results_' // T_STRING // '.tmp'
  OPEN (U_TMP, FILE=F_NAME, ACTION='READ', &
        STATUS='OLD', IOSTAT = IO_1)
  IF ( IO_1 > 0 ) THEN
    PRINT *, 'WARNING: D==0, RESTART OPEN FAILED. ', F_NAME   
    N_WARN = N_WARN + 1 ; DD = 0.d0 ; RETURN
  END IF ! OPEN FAILURE
  DO I = 1, MAX_NP
    INDEX = GET_INDEX_AT_PT (I)
    READ (U_TMP, *) (DD (INDEX (K)), K = 1, N_G_DOF)
    !b 5 FORMAT ( (10(1X, 1PE12.5)) )
  END DO
  CLOSE (U_TMP) ; PRINT *, ' '
  PRINT *, 'NOTE: GOT RESTART SOLUTION FROM ', F_NAME
END SUBROUTINE GET_STEP_RESULTS

SUBROUTINE SAVE_STEP_RESULTS (DD) 
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!          SAVE NODAL RESULTS DATA FOR PLOTS 
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Interface_Header
 IMPLICIT NONE
 CHARACTER (LEN=3)  :: T_STRING  ! for file names
 CHARACTER (LEN=20) :: F_NAME    ! file name
 REAL(DP), INTENT(IN) :: DD (N_D_FRE)
 INTEGER :: I, K, IO_1, U_TMP, GET_NEXT_IO_UNIT
!          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_TMP  = UNIT TO OUTPUT NODE RESULTS FOR PLOTS

  U_TMP = GET_NEXT_IO_UNIT ()
  T_STRING = TO3STR (THIS_STEP)
  F_NAME = 'node_results_' // T_STRING // '.tmp'
  OPEN (U_TMP, FILE=F_NAME, ACTION='WRITE', &
        STATUS='REPLACE', IOSTAT = IO_1)
  IF ( IO_1 > 0 ) THEN
    PRINT *, 'WARNING: OPEN FAILED. ', F_NAME   
    N_WARN = N_WARN + 1 ; RETURN
  END IF ! OPEN FAILURE
  DO I = 1, MAX_NP
    INDEX = GET_INDEX_AT_PT (I)
    WRITE (U_TMP, 5) (DD (INDEX (K)), K = 1, N_G_DOF)
    5 FORMAT ( (10(1X, 1PE12.5)) )
  END DO
  CLOSE (U_TMP) 
  PRINT *, 'NOTE: SOLUTION RESULTS SAVED TO ', F_NAME
  PRINT *, ' '
END SUBROUTINE SAVE_STEP_RESULTS
                         
SUBROUTINE SAVE_STEP_RESULTS_WITH_EXACT (DD, X) 
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!      SAVE NODAL AND EXACT RESULTS DATA FOR PLOTS 
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Select_Source
Use Interface_Header
 IMPLICIT NONE
 CHARACTER (LEN=3)  :: T_STRING  ! for file names
 CHARACTER (LEN=18) :: F_NAME    ! file name
 REAL(DP), INTENT(IN) :: DD (N_D_FRE)
 REAL(DP), INTENT(IN) :: X  (MAX_NP, N_SPACE)
 INTEGER :: I, K, IO_1, U_TMP, GET_NEXT_IO_UNIT
!          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

  !b IF ( U_PLT3 < 1 ) RETURN ! output turned off
  U_TMP = GET_NEXT_IO_UNIT ()
  T_STRING = TO3STR (THIS_STEP)
  F_NAME = 'node_exact_' // T_STRING // '.tmp'
  OPEN (U_TMP , FILE=F_NAME, ACTION='WRITE', &
        STATUS='REPLACE', IOSTAT = IO_1)
  IF ( IO_1 > 0 ) THEN
    PRINT *, 'WARNING: OPEN FAILED. ', F_NAME
    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_TMP , 5) (DD (INDEX (K)), K = 1, N_G_DOF),  &
                      EXACT_SOL
    5 FORMAT ( (10(1X, 1PE12.5)) )
  END DO
  CLOSE (U_TMP ) 
  PRINT *, 'NOTE: FE & EXACT RESULTS SAVED TO ', F_NAME
END SUBROUTINE SAVE_STEP_RESULTS_WITH_EXACT
                         
SUBROUTINE SAVE_STEP_EXACT_FLUXES_AT_NODES (X)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!       SAVE EXACT FLUXES TO A FILE
! * * * * * * * * * * * * * * * * * * * * * * * * * *
 Use System_Constants
Use Select_Source
 IMPLICIT NONE
 CHARACTER (LEN=3)  :: T_STRING  ! for file names
 CHARACTER (LEN=23) :: F_NAME    ! file name
 REAL(DP), INTENT(IN) :: X (MAX_NP, N_SPACE)
 REAL(DP) :: FLUX (N_R_B)
 INTEGER  :: IP, IO_1, U_TMP, GET_NEXT_IO_UNIT

! 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_STEP_EXACT_FLUXES_AT_NODES: ', &
             'NO EXACT SOLUTION USED'
    N_WARN = N_WARN + 1 ; RETURN
  END IF
  !b IF ( U_PLT9 < 1 ) RETURN  ! output turned off
  U_TMP = GET_NEXT_IO_UNIT ()
  T_STRING = TO3STR (THIS_STEP)
  F_NAME = 'exact_node_flux_' // T_STRING // '.tmp'
  OPEN (U_TMP, FILE=F_NAME, ACTION='WRITE', &
        STATUS='REPLACE', IOSTAT = IO_1)
  IF ( IO_1 > 0 ) THEN
    PRINT *, 'OPEN FAILED FOR ', F_NAME
    N_WARN = N_WARN + 1 ; 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_TMP, 5)  FLUX
    5 FORMAT ( (8(1PE13.5)) )
  END DO ! FOR ALL NODES
  CLOSE (U_TMP)
  PRINT *, 'NOTE: EXACT NODAL FLUXES SAVED TO ', F_NAME
END SUBROUTINE SAVE_STEP_EXACT_FLUXES_AT_NODES

SUBROUTINE SAVE_STEP_EXACT_NODE_SOL (X) 
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!      SAVE NODAL AND EXACT RESULTS DATA FOR PLOTS 
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Select_Source
Use Interface_Header
 IMPLICIT NONE
 CHARACTER (LEN=3)  :: T_STRING  ! for file names
 CHARACTER (LEN=22) :: F_NAME    ! file name
 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
  T_STRING = TO3STR (THIS_STEP)
  F_NAME = 'exact_node_sol_' // T_STRING // '.tmp'
  OPEN (U_PLT8, FILE=F_NAME, ACTION='WRITE', &
        STATUS='REPLACE', IOSTAT = IO_1)
  IF ( IO_1 > 0 ) THEN
    PRINT *,'WARNING: OPEN FAILED ', F_NAME
    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 ', F_NAME
END SUBROUTINE SAVE_STEP_EXACT_NODE_SOL
