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

SUBROUTINE BT_D_B_OPTIONS (D, B, S, ROWS, COLS, IOPT, CONST)
! * * * * * * * * * * * * * * * * * * * * * *
!  SPECIAL MATRIX MULTIPLICATION OPERATION
! IF IOPT=0,  S = (B)T*D*B*COEFF       PRODUCT
! IF IOPT=1,  S = (B)T*D*B*COEFF + S   NUM. INTG.
! * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 IMPLICIT NONE
 INTEGER,  INTENT(IN)  :: ROWS, COLS, IOPT
 REAL(DP), OPTIONAL    :: CONST        ! default = 1
 REAL(DP), INTENT(IN)  :: D (ROWS, ROWS), B (ROWS, COLS)
 REAL(DP), INTENT(OUT) :: S (COLS, COLS)
 ! INTEGER  :: I, J, K, L
 ! REAL(DP) :: SUM, DBIK, COEFF
 REAL(DP) :: COEFF

! D(ROWS,M)    = SYMMETRIC SQUARE MATRIX
! B(ROWS,COLS) = RECTANGLUAR ARRAY
! S(COLS,COLS) = RETURNED SYMMETRIC SQ MATRIX
! CONST        = SCALAR COEFFICIENT
  IF ( PRESENT (CONST) ) THEN
    COEFF = CONST
  ELSE 
    COEFF = 1.0_DP
  END IF
  SELECT CASE (IOPT)
    CASE (0)
      S = COEFF * MATMUL ( (MATMUL( TRANSPOSE(B),D)), B)
    CASE (1)
      S = S + COEFF * MATMUL ( (MATMUL( TRANSPOSE(B),D)), B)
    CASE DEFAULT ; STOP 'INVALID OPTION, BT_D_B_OPTIONS'
  END SELECT
  !9 use matmul and transpose here
  !DO L = 1, COLS
  !  DO K = 1, COLS
  !    SUM = 0.0_DP
  !    DO I = 1, ROWS
  !      DBIK = 0.0_DP
  !      DO J = 1, ROWS
! !         USE SYMMETRY OF D
  !        DBIK = DBIK + D (J, I) * B (J, K)
  !      END DO
  !      SUM = SUM + B (I, L) * DBIK
  !    END DO
  !    IF ( IOPT == 0 ) THEN
  !      S (L, K) = SUM * COEFF
  !    ELSE
  !      S (L, K) = S (L, K) + SUM * COEFF
  !    END IF
  !  END DO
  !END DO
END SUBROUTINE BT_D_B_OPTIONS
                     
SUBROUTINE BT_DIAG_B_OPTIONS (DIAG, B, S, ROWS, COLS, IOPT, CONST)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!  SPECIAL DIAGONIAL MATRIX MULTIPLICATION OPERATION
!  IF IOPT=0,  S = (B)T*DIAG*B*COEFF       PRODUCT ONLY
!  IF IOPT=1,  S = (B)T*DIAG*B*COEFF + S   NUM. INTEG.
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 IMPLICIT NONE
 INTEGER,  INTENT(IN)  :: ROWS, COLS, IOPT
 REAL(DP), OPTIONAL    :: CONST        ! default = 1
 REAL(DP), INTENT(IN)  :: DIAG (ROWS), B (ROWS, COLS)
 REAL(DP), INTENT(OUT) :: S (COLS, COLS)
 INTEGER  :: I, K, L
 REAL(DP) :: SUM, COEFF

! DIAG(ROWS)   = DIAGONAL MATRIX
! B(ROWS,COLS) = RECTANGLUAR ARRAY
! S(COLS,COLS) = RETURNED SYMMETRIC SQ MATRIX
! CONST        = SCALAR COEFFICIENT
  IF ( PRESENT (CONST) ) THEN
    COEFF = CONST
  ELSE 
    COEFF = 1.0_DP
  END IF
  !9 use matmul and transpose here
  DO L = 1, COLS
    DO K = 1, COLS
      SUM = 0.0_DP
      DO I = 1, ROWS
        SUM = SUM + B (I, L) * DIAG (I) * B (I, K)
      END DO
      IF ( IOPT == 0 ) THEN
        S (L, K) = SUM * COEFF
      ELSE
        S (L, K) = S (L, K) + SUM * COEFF
      END IF
    END DO
  END DO
END SUBROUTINE BT_DIAG_B_OPTIONS
                        
SUBROUTINE CONDENSE_MATRICES (N_TOTAL, N_EL_FRE, S, C)
! * * * * * * * * * * * * * * * * * * * * * * *
! CONDENSATION OF ELEMENT MATRICES TO REMOVE
!            INTERNAL DEGREES OF FREEDOM
! * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 IMPLICIT NONE
 INTEGER,  INTENT(IN)    :: N_TOTAL, N_EL_FRE
 REAL(DP), INTENT(INOUT) :: S (N_TOTAL, N_TOTAL), C (N_TOTAL)
 REAL(DP), PARAMETER     :: ZERO = 0.d0
 INTEGER  :: I, J, K, L, M, NELIM
 REAL(DP) :: CK, SKK, SLKSKK
                                           
! INTERNAL DEGREES OF FREEDOM *MUST* COME LAST
!    :  SAA : SAB  : : DA :     : CA :
!    :......:......: :....:  =  :....:
!    :  SBA : SBB  : : DB :     : CB :
! ENTER FULL S ; RETURN CONDENSED IN SAA AND CA
! DIMENSION  SAA(N_EL_FRE,N_EL_FRE), CA(N_EL_FRE)
! SAA* = (SAA) - (SAB)*(SBB)I*(SAB)T
! CA*  = (CA) - (SAB)*(SBB)I*(CB)
                                 
! N_TOTAL  = ORIG. NO. OF D.O.F. OF ELEMENT
! N_EL_FRE = FINAL NO. OF D.O.F. OF ELEMENT
! NELIM    = NUMBER OF DOF TO ELIMINATE
! S        = SQUARE ELEMENT MATRIX
! C        = ELEMENT COLUMN MATRIX
  NELIM = N_TOTAL - N_EL_FRE
  DO I = 1, NELIM
    J = N_TOTAL - I
    K = J + 1
    SKK = S (K, K)
    CK  = C (K)
    IF ( SKK /= ZERO ) THEN
      C (K) = CK / SKK
      DO L = 1, J
        SLKSKK = S (L, K) / SKK
        S (L, K) = SLKSKK
        DO M = L, J
          S (L, M) = S (L, M) - S (K, M) * SLKSKK
          S (M, L) = S (L, M)
        END DO
        C (L) = C (L) - CK * SLKSKK
      END DO
    END IF
  END DO
END SUBROUTINE CONDENSE_MATRICES

SUBROUTINE CONTOUR_ISOPAR_SURF (L_TYPE, MAX_PTS, PT_INOUT, V, COORD, &
                                XYZ_PTS, KOUNT)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! GET COORDINATES OF POINTS ON A CONTOUR ON AN ISOPARAMETRIC SURFACE
!        Int. J. Num. Meth. Engr., v. 14, n. 3, 1979
! * * * * * * ** * * * * * *  * * * * * * * * * * * * * * * * * * *
Use System_Constants ! for DP, NOD_PER_EL, N_SPACE
 IMPLICIT NONE

  INTEGER,  INTENT(IN)    :: L_TYPE, MAX_PTS
  REAL(DP), INTENT(IN)    :: COORD (NOD_PER_EL, N_SPACE), V (NOD_PER_EL)
  REAL(DP), INTENT(INOUT) :: PT_INOUT (2)      
  INTEGER,  INTENT(OUT)   :: KOUNT
  REAL(DP), INTENT(OUT)   :: XYZ_PTS (N_SPACE, MAX_PTS)

  REAL(DP) :: H_L (NOD_PER_EL), DH_L (2, NOD_PER_EL), DER (2)
  REAL(DP) :: PT_L (2, MAX_PTS), R_NEW, S_NEW
  REAL(DP) :: TRY_DL, DL, DV_DR, DV_DS, GRAD

  INTEGER  :: I, J
  LOGICAL  :: OUTSIDE 
  EQUIVALENCE (DER (1), DV_DR), (DER (2), DV_DS) 

! COORD      = GLOBAL COORD. OF NODES OF ELEMENT  
! DH_L       = LOCAL COORD. ELEM_DERIVATIVES OF H_L
! DL         = LOCAL COORD. LENGTH OF EACH SEGMENT     
! H_L        = ELEMENT INTERPOLATION FUNCTIONS                         
! KOUNT      = NO. OF PTS. ON CONTOUR CURVE           
! L_TYPE     = ELEM. TYPE, 0=QUADRILATERAL, 1=TRIANGLE                
! MAX_PTS    = MAX. NO. CONTOUR POINTS                  
! NOD_PER_EL = NUMBER OF NODES PER ELEMENT                    
! N_SPACE    = DIMENSION OF GLOBAL SPACE, 2 OR 3
! PT_INOUT   = LOCAL COORD OF 1ST (IN) OR LAST (OUT) POINT.
!              USUALLY THE INPUT POINT IS WHERE THE CONTOUR CROSSES
!              A POINT ON THE ELEMENT BOUNDARY.
! PT_L       = LOCAL COORDINATES OF THE CONTOUR PATH 
! V          = NODAL VALUES OF QUANTITY TO BE CONTOURED                
! XYZ_PTS    = GLOBAL COORD. ARRAY FOR CONTOUR PTS.

! INITIALIZE
   KOUNT = 1   ; TRY_DL = 2.5d0 / MAX_PTS
   DL = TRY_DL ; OUTSIDE = .FALSE.
   IF ( L_TYPE == 0 ) DL = DL * 2.d0 
   PT_L (1:2, 1) = PT_INOUT (1:2)

!      MARCH ALONG CONTOUR (ALLOWING FOR ONE DIRECTION REVERSAL)
   DO J = 1, MAX_PTS + 2
!      FORM SHAPE FUNCTION LOCAL ELEM_DERIVATIVES                      
     CALL SCALAR_DERIVS (PT_L (:, KOUNT), DH_L)

!      FORM LOCAL ELEM_DERIVATIVE OF VARIABLE, DER = DH*V              
     DER = MATMUL ( DH_L, V )
     GRAD = SQRT ( DOT_PRODUCT (DER, DER) )

!      FIND LOCAL COORD. OF SEGMENT END                                 
     R_NEW = PT_L (1, KOUNT) - DL * DV_DS / GRAD 
     S_NEW = PT_L (2, KOUNT) + DL * DV_DR / GRAD 

!      IS NEXT POINT IN THE ELEMENT                                     
     IF (L_TYPE == 0) THEN 
!         QUADRILATERAL (-1 TO +1): IS POINT OUTSIDE ?                  
       IF ( ABS (R_NEW) > 1.d0 .OR. &
            ABS (S_NEW) > 1.d0 ) OUTSIDE = .TRUE.
     ELSE 
!         TRIANGLE (UNIT COORDINATES): OUTSIDE ?                        
       IF ( R_NEW < 0.d0 .OR. S_NEW < 0.d0 .OR. &
           (R_NEW + S_NEW) > 1.d0) OUTSIDE = .TRUE.
     ENDIF 

     IF ( .NOT. OUTSIDE ) THEN ! ADD POINT TO CONTOUR LIST
       KOUNT = KOUNT + 1 
       PT_L (1, KOUNT) = R_NEW  ; PT_L (2, KOUNT) = S_NEW 
       IF ( KOUNT == MAX_PTS ) EXIT ! THE LOOP
     ELSE ! DID IT GO OUT ON FIRST STEP ?
       IF ( KOUNT == 1 ) THEN
         DL = -1.1d0 * DL  ! REVERSE DIRECTION
!           IS THE CONTOUR A SINGLE POINT (TWO REVERSALS)
         IF ( J > 1 ) EXIT ! THE LOOP
       ELSE ! LINE COMPLETED IN LOCAL COORDINATES
         EXIT ! THE LOOP
       END IF
     END IF
   END DO ! OVER J TRIES

!      RETURN LOCAL COORD. OF LAST POINT                                
   PT_INOUT = PT_L (1:2, KOUNT) 

!-->    CONVERT TO GLOBAL COORDINATES                                   
   DO I = 1, KOUNT 
!         EVALUATE SHAPE FUNCTIONS                                      
     CALL SCALAR_SHAPES (PT_L (:, I), H_L)
!         EVALUATE GLOBAL COORDINATES, XYZ_PTS = H*COORD         
     XYZ_PTS (1:N_SPACE, I) = MATMUL ( H_L (1:NOD_PER_EL), COORD )
   END DO 
!-->    CONTOUR FINISHED, RETURN FOR PLOTTING                           
END SUBROUTINE CONTOUR_ISOPAR_SURF                         

SUBROUTINE GET_VALUE_FROM_TABLE (T, VALUE, T_TABLE, V_TABLE, NT, NV) 
! * * * * * * * * * * * * * * * * * * * * * * * * * *               
! INTERPOLATION FOR VALUES AT T FROM TABULATED DATA                 
! * * * * * * * * * * * * * * * * * * * * * * * * * *               
Use System_Constants ! for DP
 IMPLICIT NONE
 INTEGER,  INTENT(IN)  :: NT, NV
 REAL(DP), INTENT(IN)  :: T, T_TABLE (NT), V_TABLE (NT, NV)
 REAL(DP), INTENT(OUT) :: VALUE (NV)
 INTEGER  :: FIRST, I, LAST
 REAL(DP) :: DIFF, RATIO
      
! NT      = NO INDEPENDENT VARIABLES IN T_TABLE                   
! NV      = NO DEPENDENT VARS FOR EACH INDEPENDENT VARIABLE        
! T       = GIVEN INDEPENDENT VARIABLE                               
! T_TABLE = INDEPENDENT VARIABLE TABLE                               
! V_TABLE =   DEPENDENT VARIABLE TABLE                               
! VALUE   = INTERPOLATED DEPENDENT VAR. AT T                          

  IF ( T <= T_TABLE (1) ) THEN   
    VALUE (:) = V_TABLE (1, :) ; RETURN
  END IF
  IF ( T >= T_TABLE (NT) ) THEN   
    VALUE (:) = V_TABLE (NT, :) ; RETURN
  END IF
  DO I = 2, NT 
    LAST = I 
    IF ( T <= T_TABLE (I) ) EXIT ! THIS LOOP
  END DO 
  FIRST = LAST - 1 
  DIFF   = T_TABLE (LAST) - T_TABLE (FIRST) 
  IF ( DIFF == 0.d0 ) THEN
    RATIO = 0.d0
  ELSE
    RATIO = (T - T_TABLE (FIRST) ) / DIFF  
  END IF 
  VALUE (:) = (1.d0 - RATIO) * V_TABLE (FIRST, :) &
            +          RATIO * V_TABLE (LAST, :) 
END SUBROUTINE GET_VALUE_FROM_TABLE                         
                                
SUBROUTINE G_DERIV (AJ_INV, DELTA, GLOBAL)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! N_SPACE GLOBAL DERIVATIVES OF LT_FREE INTERPOLATION
!           FUNCTIONS AT A LOCAL POINT.
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data
 REAL(DP), INTENT(IN)  :: AJ_INV (N_SPACE, N_SPACE)
 REAL(DP), INTENT(IN)  :: DELTA  (N_SPACE, LT_FREE)
 REAL(DP), INTENT(OUT) :: GLOBAL (N_SPACE, LT_FREE)

! N_SPACE    = DIMENSION OF SPACE
! LT_FREE    = NUMBER OF DEGREES OF FREEDOM PER ELEMENT
! AJ_INV     = INVERSE JACOBIAN MATRIX AT LOCAL POINT
! DELTA      = LOCAL COORD DERIV AT POINT OF INTEREST
! GLOBAL     = GLOBAL DERIVATIVES MATRIX AT LOCAL POINT

!              GLOBAL = AJ_INV*DELTA
    GLOBAL = MATMUL (AJ_INV, DELTA)
END SUBROUTINE G_DERIV
                     
SUBROUTINE GLOBAL_DERIV_OF_H (GLOBAL)
!b SUBROUTINE GLOBAL_DERIV_OF_H (COORD, DLH, GLOBAL)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! N_SPACE GLOBAL DERIVATIVES OF LT_FREE INTERPOLATION
!           FUNCTIONS AT A LOCAL POINT.
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
Use Elem_Type_Data ! for COORD (LT_N, N_SPACE), DLH (N_SPACE, LT_N)
 !b REAL(DP), INTENT(IN)  :: COORD  (LT_N, N_SPACE)
 !b REAL(DP), INTENT(IN)  :: DLH  (N_SPACE, LT_N)
 REAL(DP), INTENT(OUT) :: GLOBAL (N_SPACE, LT_N)

 REAL(DP) :: AJ (N_SPACE, N_SPACE), AJ_INV (N_SPACE, N_SPACE), DET

! AJ         = JACOBIAN MATRIX AT LOCAL POINT
! AJ_INV     = INVERSE JACOBIAN MATRIX AT LOCAL POINT
! COORD      = COORDINATES OF ALL NODES ON ELEMENT
! DLH      = LOCAL COORD DERIV AT POINT OF INTEREST
! GLOBAL     = GLOBAL DERIVATIVES MATRIX AT LOCAL POINT
! LT_N       = MAXIMUM NUMBER OF NODES FOR ELEMENT TYPE
! N_SPACE    = DIMENSION OF SPACE

! AJ = DLH*COORD,  GLOBAL = AJ_INV*DLH
    AJ     = MATMUL (DLH, COORD)
    CALL INVERT_JACOBIAN (AJ, AJ_INV, DET, N_SPACE)
    GLOBAL = MATMUL (AJ_INV, DLH)
END SUBROUTINE GLOBAL_DERIV_OF_H
                     
!wSUBROUTINE GLOBAL_2ND_DERIV_OF_H (GLOBAL, GLOBAL_2)
!w!b SUBROUTINE GLOBAL_2ND_DERIV_OF_H (COORD, DLH, D2LH, &
                                  !w!b GLOBAL, GLOBAL_2)
!w! * * * * * * * * * * * * * * * * * * * * * * * * * *
! N_SPACE GLOBAL FIRST AND SECOND DERIVATIVES OF LT_N
!w!  INTERPOLATION FUNCTIONS AT A LOCAL POINT.
!w! * * * * * * * * * * * * * * * * * * * * * * * * * *
!wUse System_Constants ! for  N_SPACE, N_2_DER
!wUse Elem_Type_Data ! for LT_N,COORD (LT_N, N_SPACE),
                   !w! DLH (N_SPACE, LT_N), D2LH (N_2_DER, LT_N)
 !w!b REAL(DP), INTENT(IN)  :: COORD  (LT_N, N_SPACE)
 !w!b REAL(DP), INTENT(IN)  :: DLH    (N_SPACE, LT_N)
 !w!b REAL(DP), INTENT(IN)  :: D2LH  (N_2_DER, LT_N) !w
 !wREAL(DP), INTENT(OUT) :: GLOBAL   (N_SPACE, LT_N)
 !wREAL(DP), INTENT(OUT) :: GLOBAL_2 (N_2_DER, LT_N)
!w
 !wREAL(DP) :: AJ (N_SPACE, N_SPACE), AJ_INV (N_SPACE, N_SPACE), DET
 !wREAL(DP) :: P_AJ (N_2_DER, N_2_DER), P_AJ_INV (N_2_DER, N_2_DER)
 !wINTEGER  :: I_ERROR = 0
!w
!w! AJ         = JACOBIAN MATRIX AT LOCAL POINT
!w! AJ_INV     = INVERSE JACOBIAN MATRIX AT LOCAL POINT
!w! DLH       = LOCAL COORD DERIV AT POINT OF INTEREST
!w! D1LH      = LOCAL COORD SECOND DERIV AT POINT OF INTEREST
!w! GLOBAL     = GLOBAL DERIVATIVES MATRIX AT LOCAL POINT
!w! GLOBAL_2   = GLOBAL SECOND DERIVATIVES MATRIX AT LOCAL POINT
!w! LT_N       = MAXIMUM NUMBER OF NODES FOR ELEMENT TYPE
!w! N_SPACE    = DIMENSION OF SPACE
!w! N_2_DER    = NUMBER OF SECOND DERIVATIVES, 1, 3, OR 6
!w! P_AJ       = PRODUCT OF JACOBIAN TERMS 
!w! P_AJ_INV   = INVERSE OF PRODUCT OF JACOBIAN TERMS
!w
!w! AJ = DLH*COORD,  GLOBAL = AJ_INV*DLH
    !wAJ     = MATMUL (DLH, COORD)
    !wCALL INVERT_JACOBIAN (AJ, AJ_INV, DET, N_SPACE)
    !wGLOBAL = MATMUL (AJ_INV, DLH)
!w
!w! NOW GET SECOND DERIVATIVE: 
!w! GLOBAL_2 = P_AJ_INV (D2HL - HESSIAN * GLOBAL)
!w! HESSIAN  = D2HL * COORD
    !wCALL JACOBIAN_PRODUCTS (AJ, P_AJ)
    !wCALL INV_SMALL_MAT (N_2_DER, P_AJ, P_AJ_INV, I_ERROR)
    !wIF ( I_ERROR > 0 ) THEN
      !wPRINT *, 'WARNING: GLOBAL_2ND_DERIV_OF_H, BAD INVERSE'
      !wN_WARN = N_WARN + 1
    !wEND IF
    !wGLOBAL_2 = MATMUL (P_AJ_INV,                                &
                 !w(D2LH - MATMUL ( MATMUL (D2LH, COORD), GLOBAL)))
!wEND SUBROUTINE GLOBAL_2ND_DERIV_OF_H
                     
FUNCTION GET_REAL_IDENTITY (N)  RESULT (EYE)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!          COMPUTE A REAL N BY N IDENTITY MATRIX
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 INTEGER,  INTENT(IN)  :: N            ! SIZE OF MATRIX
 REAL(DP)              :: EYE (N, N)   ! IDENTITY MATRIX
 
 INTEGER :: J
 REAL(DP), PARAMETER :: I_2 (2, 2)                       &
                      = RESHAPE ((/ 1.d0, 0.d0,          &
                                    0.d0, 1.d0 /),       &
                                    (/ 2, 2/) )
 REAL(DP), PARAMETER :: I_3 (3, 3)                       &
                      = RESHAPE ((/ 1.d0, 0.d0, 0.d0,    &
                                    0.d0, 1.d0, 0.d0,    &
                                    0.d0, 0.d0, 1.d0 /), &
                                    (/ 3, 3/) )
  SELECT CASE (N)
    CASE (:0) ; STOP 'INVALID SIZE, GET_REAL_IDENTITY'
    CASE (1)  ; EYE = 1.D0
    CASE (2)  ; EYE = I_2
    CASE (3)  ; EYE = I_3
    CASE DEFAULT
      EYE = 0.D0
      DO J = 1, N ; EYE (J, J) = 1.D0 ; END DO
  END SELECT
END FUNCTION GET_REAL_IDENTITY 

SUBROUTINE REAL_IDENTITY (N, EYE)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!          COMPUTE A REAL N BY N IDENTITY MATRIX
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 INTEGER,   INTENT(IN)  :: N            ! SIZE OF MATRIX
 REAL(DP),  INTENT(OUT) :: EYE (N, N)   ! IDENTITY MATRIX
 
 INTEGER :: J
 REAL(DP), PARAMETER :: I_2 (2, 2)                       &
                      = RESHAPE ((/ 1.d0, 0.d0,          &
                                    0.d0, 1.d0 /),       &
                                    (/ 2, 2/) )
 REAL(DP), PARAMETER :: I_3 (3, 3)                       &
                      = RESHAPE ((/ 1.d0, 0.d0, 0.d0,    &
                                    0.d0, 1.d0, 0.d0,    &
                                    0.d0, 0.d0, 1.d0 /), &
                                    (/ 3, 3/) )
  SELECT CASE (N)
    CASE (:0) ; STOP 'INVALID SIZE, GET_REAL_IDENTITY'
    CASE (1)  ; EYE = 1.D0
    CASE (2)  ; EYE = I_2
    CASE (3)  ; EYE = I_3
    CASE DEFAULT
      EYE = 0.D0
      DO J = 1, N ; EYE (J, J) = 1.D0 ; END DO
  END SELECT
END SUBROUTINE REAL_IDENTITY 

SUBROUTINE INVERT_JACOBIAN (AJ, AJ_INV, DET, P_SPACE)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! FIND INVERSE AND DETERMINANT OF GEOMETRY JACOBIAN
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 INTEGER,  INTENT(IN ) :: P_SPACE
 REAL(DP), INTENT(IN ) :: AJ (P_SPACE, P_SPACE)
 REAL(DP), INTENT(OUT) :: AJ_INV (P_SPACE, P_SPACE), DET
 INTEGER :: I_ERROR = 0
                                                       
! P_SPACE = NUMBER OF SPATIAL DIMENSIONS
! AJ      = JACOBIAN MATRIX AT A POINT
! AJ_INV   = INVERSE OF AJ
! DET     = DETERMINANT OF AJ
                             
  SELECT CASE (P_SPACE)
!-->  1-D
    CASE (1)
      DET = AJ(1,1)
      IF ( ABS(DET) <= SPACING(DET) ) THEN
        WRITE (N_PRT, *) 'WARNING: INVERT 1x1, NEARLY SINGULAR'
        WRITE (N_PRT, *) 'ELEMENT, DET = ', THIS_EL, DET
        N_WARN = N_WARN + 1
        I_ERROR = 1
        !b IF ( DET == 0.d0 ) STOP 'BAD DET, INVERT_JACOBIAN'
      ELSE
        AJ_INV(1,1) = 1.d0/DET
      END IF
!-->  2-D
    CASE (2)
      CALL INVERT_2BY2 (AJ, AJ_INV, DET, I_ERROR)
!-->  3-D
    CASE (3)
      CALL  INVERT_3BY3 (AJ, AJ_INV, DET, I_ERROR)
!-->  Error
    CASE DEFAULT
      PRINT *,'FATAL ERROR IN ELEMENT ', THIS_EL
      STOP 'BAD PARAMETRIC SPACE, INVERT_JACOBIAN'
  END SELECT
  IF ( DET <= 0.d0 ) THEN
    WRITE (N_PRT, *) 'WARNING: CORRECTED NEGATIVE JACOBIAN'
    WRITE (N_PRT, *) 'ELEMENT, DET = ', THIS_EL, DET
    N_WARN = N_WARN + 1
  END IF
END SUBROUTINE INVERT_JACOBIAN
                              
SUBROUTINE INVERT_3BY3 (A, AINV, DET, I_ERROR)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! FIND INVERSE AND DETERMINANT OF MATRIX A(3,3)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: A(3,3)
 REAL(DP), INTENT(OUT) :: AINV(3,3), DET
 INTEGER,  INTENT(INOUT) :: I_ERROR
                                        
! A    = ORIGINAL MATRIX
! AINV = INVERSE OF MATRIX A
! DET  = DETERMINANT OF A
                         
  AINV(1,1) =  A(2,2)*A(3,3) - A(3,2)*A(2,3)
  AINV(2,1) = -A(2,1)*A(3,3) + A(3,1)*A(2,3)
  AINV(3,1) =  A(2,1)*A(3,2) - A(3,1)*A(2,2)
  AINV(1,2) = -A(1,2)*A(3,3) + A(3,2)*A(1,3)
  AINV(2,2) =  A(1,1)*A(3,3) - A(3,1)*A(1,3)
  AINV(3,2) = -A(1,1)*A(3,2) + A(3,1)*A(1,2)
  AINV(1,3) =  A(1,2)*A(2,3) - A(2,2)*A(1,3)
  AINV(2,3) = -A(1,1)*A(2,3) + A(2,1)*A(1,3)
  AINV(3,3) =  A(1,1)*A(2,2) - A(2,1)*A(1,2)
  DET = A(1,1)*AINV(1,1) + A(1,2)*AINV(2,1) &
      + A(1,3)*AINV(3,1)
  IF ( ABS(DET) <= SPACING(DET) ) THEN
    WRITE (N_PRT, *) 'WARNING: INVERT_3BY3, NEARLY SINGULAR MATRIX'
    WRITE (N_PRT, *) 'ELEMENT, DET = ', THIS_EL, DET
    N_WARN = N_WARN + 1
    I_ERROR = 1
    !b IF ( DET == 0.d0 ) STOP 'SINGULAR 3X3 MATRIX IN INVERT_3BY3'
  ELSE
    AINV = AINV/DET
  END IF
END SUBROUTINE INVERT_3BY3
                          
SUBROUTINE INVERT_2BY2 (A, AINV, DET, I_ERROR)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! CALCULATE THE DETERMINANT AND INVERSE OF A(2,2)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 implicit none
 REAL(DP), INTENT(IN)  :: A (2, 2)
 REAL(DP), INTENT(OUT) :: AINV (2, 2), DET
 INTEGER,  INTENT(INOUT) :: I_ERROR

! A    = ORIGINAL MATRIX
! AINV = INVERSE OF MATRIX A
! DET  = DETERMINANT OF A

  DET = A (1, 1) * A (2, 2) - A (1, 2) * A (2, 1)
  IF ( ABS(DET) <= SPACING(DET) ) THEN
    WRITE (N_PRT, *) 'WARNING: INVERT_2BY2, NEARLY SINGULAR MATRIX'
    WRITE (N_PRT, *) 'ELEMENT, DET = ', THIS_EL, DET
    N_WARN = N_WARN + 1
    I_ERROR = 1
  END IF
  IF ( DET == 0.d0 ) THEN
    PRINT *, "SINGULAR MATRIX:" ; CALL RPRINT (A, 2, 2, 1)
    STOP 'SINGULAR 2X2 MATRIX, INVERT_2BY2'
  ELSE
    AINV (1, 1) =   A (2, 2) / DET
    AINV (1, 2) = - A (1, 2) / DET
    AINV (2, 1) = - A (2, 1) / DET
    AINV (2, 2) =   A (1, 1) / DET
  END IF ! singular
END SUBROUTINE INVERT_2BY2
                          
SUBROUTINE JACOBIAN_AT_PT (DELTA, COORD, AJ)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
!  CALCULATE THE JACOBIAN MATRIX AT A LOCAL POINT
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants  ! for N_SPACE, NOD_PER_EL
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: DELTA (N_SPACE,    NOD_PER_EL)
 REAL(DP), INTENT(IN)  :: COORD (NOD_PER_EL, N_SPACE)
 REAL(DP), INTENT(OUT) :: AJ    (N_SPACE,    N_SPACE)
 
! NOD_PER_EL = NUMBER OF NODES PER ELEMENT
! N_SPACE    = DIMENSION OF SPACE
! DELTA      = LOCAL  DERIVATIVES OF NOD_PER_EL INTERPOLATION
!              FUNCTIONS AT POINT OF INTEREST.
! COORD      = SPATIAL COORDINATES OF ELEMENT'S NODES
! AJ         = JACOBIAN MATRIX = DELTA*COORD
   AJ = MATMUL (DELTA, COORD)
END SUBROUTINE JACOBIAN_AT_PT
                    
FUNCTION GEOMETRIC_JACOBIAN () RESULT (AJ)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!  GET JACOBIAN MATRIX AT A POINT, LT_GEOM <= LT_N
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants   ! For N_SPACE
Use Elem_Type_Data  ! for: LT_GEOM, LT_PARM,
    ! GEOMETRY (LT_GEOM, N_SPACE),         
    ! DLG (LT_PARM, LT_GEOM)
 IMPLICIT NONE
 REAL(DP) :: AJ    (LT_PARM, LT_PARM)  ! Out
 
! LT_GEOM = NUMBER OF NODES PER ELEMENT
! LT_PARM = DIMENSION OF PARAMETRIC SPACE
! N_SPACE = DIMENSION OF SPACE
! DLG   = LOCAL  DERIVATIVES OF LT_GEOM INTERPOLATION
!              FUNCTIONS AT POINT OF INTEREST.
! GEOMETRY   = SPATIAL COORDINATES OF GEOMETRIC NODES
! AJ      = JACOBIAN MATRIX = DLG*GEOMETRY

   AJ = MATMUL (DLG, GEOMETRY (:, 1:LT_PARM))
END FUNCTION GEOMETRIC_JACOBIAN
                    
SUBROUTINE DIAGONALIZE_SQ_MATRIX (NL, SQ, DSQ)
! * * * * * * * * * * * * * * * * * * * * *
! Diagonialize a consistent square matrix
! * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 IMPLICIT NONE
 INTEGER,  INTENT(IN)  :: NL
 REAL(DP), INTENT(IN)  :: SQ (NL, NL)
 REAL(DP), INTENT(OUT) :: DSQ (NL)
 INTEGER  :: I
 REAL(DP) :: DIAG, TOTAL, RATIO
                                             
! NL  = size of array
! SQ  = Original square matrix
! DSQ = Diagonal form

  DIAG  = 0.0_DP
  TOTAL = SUM ( SQ )              ! sum of all terms
  DO I = 1, NL                    !  Count diagonal sum
    DIAG = DIAG + SQ (I, I)
  END DO
!  Scale diagonal
  IF (DIAG <= 0.0_DP) STOP 'ERROR IN DIAGONALIZE_SQ_MATRIX'
    RATIO = TOTAL / DIAG
    DO I = 1, NL
      DSQ (I) = RATIO * SQ (I, I)
    END DO
END SUBROUTINE DIAGONALIZE_SQ_MATRIX
                      
SUBROUTINE LINE_2D_GEOM (I_OPT, COORD, DX, DY, DS,    &
                         X_NORM, Y_NORM, X_TAN, Y_TAN)
! * * * * * * * * * * * * * * * * * * * * * * * * * *               
! GEOMETRIC PROPERTIES OF A STRAIGHT LINE IN X-Y SPACE              
! * * * * * * * * * * * * * * * * * * * * * * * * * *               
Use System_Constants   ! for DP
 IMPLICIT NONE 
  INTEGER,  INTENT(IN)  :: I_OPT
  REAL(DP), INTENT(IN)  :: COORD (2, 2) 
  REAL(DP), INTENT(OUT) :: DX, DY, DS, X_NORM, Y_NORM, X_TAN, Y_TAN

!  COORD         = GIVEN COORDINATES OF TWO DEFINING POINTS
!  DX,DY,DS      = X,Y COMPONENTS AND TOTAL LENGTH     N
!  X_NORM,Y_NORM = COMPONENTS OF UNIT NORMAL           1>>>2  V
!  X_TAN,Y_TAN   = COMPONENTS OF UNIT TANGENT FROM 1 TO 2         
!  I_OPT         = 0, CALCULATE DISTANCES ONLY              
!  I_OPT         = 1, CALCULATE DISTANCES AND NORMAL VECTOR  
!  I_OPT         = 2, CALCULATE DIST, NORMAL, TANGENT

   DX = COORD (2, 1) - COORD (1, 1) 
   DY = COORD (2, 2) - COORD (1, 2) 
   DS = SQRT (DX * DX + DY * DY) 
   IF (I_OPT <= 0) RETURN 
   IF (DS == 0.d0) THEN 
     WRITE (6, *) 'ZERO LENGTH IN LINE_2D_GEOM'
     STOP         'ZERO LENGTH IN LINE_2D_GEOM'
   ELSE
     X_NORM = DY / DS 
     Y_NORM = - DX / DS 
     IF (I_OPT == 1) RETURN 
       X_TAN = DX / DS 
       Y_TAN = DY / DS 
   END IF
END SUBROUTINE LINE_2D_GEOM                         

SUBROUTINE SQ_TO_DIAGONAL_SQ_MATRIX (NL, SQ, DSQ)
! * * * * * * * * * * * * * * * * * * * * *
! Diagonialize a consistent square matrix
! * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 IMPLICIT NONE
 INTEGER,  INTENT(IN)  :: NL
 REAL(DP), INTENT(IN)  :: SQ  (NL, NL)
 REAL(DP), INTENT(OUT) :: DSQ (NL, NL)
 INTEGER  :: I
 REAL(DP) :: DIAG, TOTAL, RATIO
                                             
! NL  = size of array
! SQ  = Original square matrix
! DSQ = Diagonal form

  DIAG  = 0.0_DP
  TOTAL = SUM ( SQ )              ! sum of all terms
  DO I = 1, NL                    !  Count diagonal sum
    DIAG = DIAG + SQ (I, I)
  END DO
!  Scale diagonal
  IF (DIAG <= 0.0_DP) STOP 'ERROR IN SQ_TO_DIAGONAL_SQ_MATRIX'
    RATIO = TOTAL / DIAG
    DO I = 1, NL
      DSQ (I, I) = RATIO * SQ (I, I)
    END DO
END SUBROUTINE SQ_TO_DIAGONAL_SQ_MATRIX
                      
FUNCTION OUTER_PRODUCT (A, B) RESULT (C)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! OUTER PRODUCT OF TWO VECTORS, C_IJ = A_I * B_J
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 REAL(DP), INTENT(IN) :: A (:), B (:)
 REAL(DP)             :: C (SIZE(A), SIZE(B))

   C = SPREAD (A, DIM=2, NCOPIES=SIZE(B)) &
     * SPREAD (B, DIM=1, NCOPIES=SIZE(A))
END FUNCTION OUTER_PRODUCT 

FUNCTION OUTER_DIVIDE  (A, B) RESULT (C)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! OUTER QUOTIENT OF TWO VECTORS, C_IJ = A_I / B_J
! * * * * * * * * * * * * * * * * * * * * * * * * * *
Use Precision_Module
 REAL(DP), INTENT(IN) :: A (:), B (:)
 REAL(DP)             :: C (SIZE(A), SIZE(B))

   C = SPREAD (A, DIM=2, NCOPIES=SIZE(B)) &
     / SPREAD (B, DIM=1, NCOPIES=SIZE(A))
END FUNCTION OUTER_DIVIDE  

FUNCTION OUTER_AND (A, B) RESULT (C)
! * * * * * * * * * * * * * * * * * * * * * * * * * *
! OUTER LOGICAL 'AND' OF TWO VECTORS, C_IJ = A_I & B_J
! * * * * * * * * * * * * * * * * * * * * * * * * * *
 LOGICAL, INTENT(IN) :: A (:), B (:)
 LOGICAL             :: C (SIZE(A), SIZE(B))

   C =     SPREAD (A, DIM=2, NCOPIES=SIZE(B)) &
     .AND. SPREAD (B, DIM=1, NCOPIES=SIZE(A))
END FUNCTION OUTER_AND

!SUBROUTINE NGRAND (WT, DET, V, DGV, XPT, NOD_PER_EL, N_SPACE, &
!                   N_EL_FRE, COL, SQ, N_FILE )
!! * * * * * * * * * * * * * * * * * * * * * * * * * *
!!     PROBLEM DEPENDENT INTEGRAND EVALUATION IN
!!     AN ISOPARAMETRIC OR SUBPARAMETRIC ELEMENT
!! * * * * * * * * * * * * * * * * * * * * * * * * * *
!Use Precision_Module
!IMPLICIT NONE
!INTEGER,  INTENT(IN)    :: NOD_PER_EL, N_SPACE, N_EL_FRE, N_FILE
!REAL(DP), INTENT(IN)    :: V (N_EL_FRE), DGV (N_SPACE, N_EL_FRE)
!REAL(DP), INTENT(IN)    :: XPT (N_SPACE), WT, DET
!REAL(DP), INTENT(INOUT) :: COL (N_EL_FRE), SQ (N_EL_FRE, N_EL_FRE)
! 
!! NOD_PER_EL = NUMBER OF NODES PER ELEMENT
!! N_SPACE    = NUMBER OF SPATIAL DIMENSIONS
!! N_EL_FRE   = NUMBER OF ELEMENT DEGREES OF FREEDOM
!! V          = ELEMENT INTERPOLATION FUNCTIONS
!! DGV        = GLOBAL DERIVATIVES OF V
!! XPT        = GLOBAL COORDS OF THE POINT
!! WT         = QUADRATURE WEIGHT AT POINT
!! DET        = JACOBIAN DETERMINANT AT POINT
!! COL        = PROB DEP COLUMN MATRIX INTEGRAND
!! SQ         = PROB DEP SQUARE MATRIX INTEGRAND
!! N_FILE      = STORAGE UNIT FOR POST SOLUTION DATA
!! ....................................................
!! *** NGRAND PROBLEM DEPENDENT STATEMENTS FOLLOW ***
!! ....................................................
!  !inc include 'my_ngrand_inc'
!END SUBROUTINE NGRAND

!FUNCTION COSH (X)  RESULT (VALUE) ! Hyperbolic cosine
!Use Precision_module
!IMPLICIT NONE
! REAL(DP), INTENT(IN) :: X
! REAL(DP)             :: VALUE
!   VALUE = (EXP (X) + EXP(-X)) * 0.5d0
!END FUNCTION COSH 
!
!FUNCTION SINH (X)  RESULT (VALUE) ! Hyperbolic sine
!Use Precision_module
!IMPLICIT NONE
! REAL(DP), INTENT(IN) :: X
! REAL(DP)             :: VALUE
!   VALUE = (EXP (X) - EXP(-X)) * 0.5d0
!END FUNCTION SINH 

SUBROUTINE SYMMETRIC_INVERSE (A, A_INV, I_ERROR)
! ------------------------------------------------------
! INVERT SYMMETRIC MATRIX A(N,N) 
! ------------------------------------------------------
Use Precision_Module
 IMPLICIT NONE
 REAL (DP), INTENT (IN)  :: A (:, :)
 REAL (DP), INTENT (OUT) :: A_INV (SIZE(A, 1), SIZE(A, 1))
 INTEGER,   INTENT (OUT) :: I_ERROR  ! FATAL RESULT IF /= 0

 INTEGER   :: I, J, K, N
 REAL (DP) :: D

! A  = ORIGINAL MATRIX
! A_INV = INVERSE OF MATRIX A
! D  = A DIAGONAL ELEMENT OF A
! N  = SIZE OF A

  I_ERROR = 0 ; N = SIZE(A, 1) ; A_INV = A       ! INITIALIZE

  DO K = 1,N                       ! LOOP OVER ROWS
    D = A_INV (K,K)                ! DIAGONAL IS PIVOT TERM
    !b IF ( D == 0.d0 ) THEN          ! FATAL CONDITION
    IF ( ABS(D) <= SPACING(D) ) THEN          ! FATAL CONDITION
      PRINT *, 'ERROR: ZERO PIVOT IN SYMMETRIC_INVERSE'
      I_ERROR = 1 ; RETURN         ! ABORT
    END IF ! FATAL ERROR
    A_INV (K, :) = -A_INV (K, :)/D ! SCALE ROW
    DO I = 1, N                    ! LOOP OVER ROW AGAIN
      IF ( I /= K )  THEN          ! OFF-DIAGONAL TERMS
        DO J = 1, N                ! LOOP OVER COLUMNS
          IF ( J /= K ) A_INV (I,J) =                       &
                        A_INV (I,J) + A_INV (I,K)*A_INV (K,J) 
        END DO ! OVER J
      END IF ! OFF DIAGONAL I
      A_INV (I,K) = A_INV (I,K)/D
    END DO ! OVER ROW I   
    A_INV (K,K) = 1.d0/D
  END DO ! OVER K-TH COLUMN
END SUBROUTINE SYMMETRIC_INVERSE

SUBROUTINE  INVERT_SMALL_MATRIX (A, A_INV, I_ERROR)
! ------------------------------------------------------
! INVERT A SMALL MATRIX A(N,N) 
! ------------------------------------------------------
Use Precision_Module
 IMPLICIT NONE
 REAL (DP), INTENT (IN)  :: A (:, :)
 REAL (DP), INTENT (OUT) :: A_INV (SIZE(A, 1), SIZE(A, 1))
 INTEGER,   INTENT (OUT) :: I_ERROR  ! FATAL RESULT IF /= 0

 INTEGER   :: I, K, N
 REAL (DP) :: D

! A  = ORIGINAL MATRIX
! A_INV = INVERSE OF MATRIX A
! D  = A DIAGONAL ELEMENT OF A
! N  = SIZE OF A

  I_ERROR = 0 ;N = SIZE(A, 1) ; A_INV = A       ! INITIALIZE

  DO  K = 1,N                      ! LOOP OVER DIAGONALS
    D = A_INV (K, K)               ! PIVOT TERM
    !b IF ( D == 0.d0 ) THEN          ! FATAL CONDITION
    IF ( ABS(D) <= SPACING(D) ) THEN          ! FATAL CONDITION
      PRINT *, 'ERROR: ZERO PIVOT IN INVERT_SMALL_MATRIX'
      I_ERROR = 1 ; RETURN         ! ABORT
    END IF ! FATAL ERROR
    A_INV (K, K) = 1.d0
    A_INV (K, :) = A_INV (K, :)/D  ! SCALE ROW
    DO I = 1, N                    ! LOOP OVER ROW AGAIN
      IF ( I == K ) CYCLE ! TO NEXT I
        D            = A_INV (I, K)                
        A_INV (I, K) = 0.d0                
        A_INV (I, :) = A_INV (I, :) - D*A_INV (K, :) 
    END DO ! OVER ROW I   
  END DO ! OVER K-TH COLUMN
END SUBROUTINE INVERT_SMALL_MATRIX

SUBROUTINE INV_SMALL_MAT (N, A, A_INV, I_ERROR)
! ------------------------------------------------------
! INVERT A SMALL MATRIX A(N,N) 
! ------------------------------------------------------
Use System_Constants
 IMPLICIT NONE
 INTEGER,   INTENT (IN)  :: N
 REAL (DP), INTENT (IN)  :: A     (N, N)
 REAL (DP), INTENT (OUT) :: A_INV (N, N)
 INTEGER,   INTENT (OUT) :: I_ERROR  ! FATAL RESULT IF /= 0

 INTEGER   :: I, K
 REAL (DP) :: D, DET

! A     = ORIGINAL MATRIX
! A_INV = INVERSE OF MATRIX A
! D     = A DIAGONAL ELEMENT OF A
! N     = SIZE OF A

  I_ERROR = 0 ; A_INV = A       ! INITIALIZE

  SELECT CASE (N) !-->  1-D
    CASE (1)
      DET = A (1,1)
      IF ( ABS(DET) <= SPACING(DET) ) THEN
        WRITE (N_PRT, *) 'WARNING: INVERT 1x1, NEARLY SINGULAR'
        WRITE (N_PRT, *) 'DET = ', DET
        N_WARN = N_WARN + 1
        I_ERROR = 1
      ELSE
        A_INV(1,1) = 1.d0/DET
      END IF
    CASE (2) !-->  2-D
      CALL INVERT_2BY2 (A, A_INV, DET, I_ERROR)
    CASE (3) !-->  3-D
      CALL INVERT_3BY3 (A, A_INV, DET, I_ERROR)
    CASE DEFAULT
      DO  K = 1,N                      ! LOOP OVER DIAGONALS
        D = A_INV (K, K)               ! PIVOT TERM
        IF ( ABS(D) <= SPACING(D) ) THEN          ! FATAL CONDITION
          PRINT *, 'WARNING: ZERO PIVOT IN INV_SMALL_MAT'
          N_WARN = N_WARN + 1
          I_ERROR = 1 
          EXIT ! THIS DO ON K  !b
        END IF ! FATAL ERROR
        A_INV (K, K) = 1.d0
        A_INV (K, :) = A_INV (K, :)/D  ! SCALE ROW
        DO I = 1, N                    ! LOOP OVER ROW AGAIN
          IF ( I == K ) CYCLE ! TO NEXT I
            D            = A_INV (I, K)        
            A_INV (I, K) = 0.d0        
            A_INV (I, :) = A_INV (I, :) - D*A_INV (K, :) 
        END DO ! OVER ROW I   
      END DO ! OVER K-TH COLUMN
  END SELECT
  IF ( I_ERROR > 0 ) THEN
    PRINT *, 'ERROR: INV_SMALL_MAT FAILED, IDENTITY MATRIX RETURNED'
    CALL REAL_IDENTITY (N, A_INV)
  END IF ! INVERSION FAILED
END SUBROUTINE INV_SMALL_MAT

SUBROUTINE JACOBIAN_PRODUCTS (AJ, P_AJ)
!* * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! GET JACOBIAN PRODUCTS FOR SECOND DERIVATIVE CALCULATION
!* * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
  REAL (DP), INTENT (IN)  :: AJ   (N_SPACE, N_SPACE)
  REAL (DP), INTENT (OUT) :: P_AJ (N_2_DER, N_2_DER)

! AJ    = JACOBIAN MATRIX AT A POINT
! P_AJ  = JACOBIAN PRODUCTS AT POINT IN TERMS OF THE 2ND GLOBAL 
!         DERIVATIVES (GENERALLY WANT P_AJ INVERSE)
! V_,AB = P_AJ_IAJB * V_,IJ  (INDEX NOTATION: _ PRECEEDS SUBSCRIPTS)
!         WHERE P_AJ_IAJB = X_I,A * X_J,B

   SELECT CASE ( N_2_DER ) ! NUMBER OF SECOND DERIVATIVES
     CASE ( 1 ) ! One-dimensional space
!         | V,AA | = (P_AJ) | V,II |    A,B IMPLY LOCAL COORD
!                                       I,J IMPLY GLOBAL COORD
       P_AJ (1, 1) = AJ (1, 1) **2

     CASE ( 3 ) ! Two-dimensional space
!         | V,AA |          | V,II |    A,B IMPLY LOCAL COORD
!         | V,AB | = (P_AJ) | V,IJ |    I,J IMPLY GLOBAL COORD
!         | V,BB |          | V,JJ |

       P_AJ (1, 1) = AJ (1, 1) **2
       P_AJ (2, 1) = AJ (1, 1) * AJ (2, 1)
       P_AJ (3, 1) = AJ (2, 1) **2
       P_AJ (1, 2) = 2.d0 * AJ (1, 1) * AJ (1, 2)
       P_AJ (2, 2) = AJ (1, 1) * AJ (2, 2) + AJ (2, 1) * AJ (1, 2)
       P_AJ (3, 2) = 2.d0 * AJ (2, 1) * AJ (2, 2)
       P_AJ (1, 3) = AJ (1, 2) **2
       P_AJ (2, 3) = AJ (1, 2) * AJ (2, 2)
       P_AJ (3, 3) = AJ (2, 2) **2

     CASE ( 6 ) ! Three-dimensional space
!         | V,AA |          | V,II |    A,B,C IMPLY LOCAL        
!         | V,AB |          | V,IJ |    I,J,K IMPLY GLOBAL        
!         | V,AC | = (P_AJ) | V,IK |        
!         | V,BB |          | V,JJ |        
!         | V,BC |          | V,JK |        
!         | V,CC |          | V,KK |

       P_AJ (1, 1) = AJ (1, 1) **2 
       P_AJ (2, 1) = AJ (1, 2) * AJ (2, 2) 
       P_AJ (3, 1) = AJ (1, 1) * AJ (3, 1) 
       P_AJ (4, 1) = AJ (2, 1) **2 
       P_AJ (5, 1) = AJ (2, 1) * AJ (3, 1) 
       P_AJ (6, 1) = AJ (3, 1) **2 
       P_AJ (1, 2) = AJ (1, 1) * AJ (1, 2) * 2.d0 
       P_AJ (2, 2) = AJ (1, 1) * AJ (2, 2) + AJ (1, 2) * AJ (2, 1) 
       P_AJ (3, 2) = AJ (1, 1) * AJ (3, 2) + AJ (1, 2) * AJ (3, 1) 
       P_AJ (4, 2) = AJ (2, 1) * AJ (2, 2) * 2.d0 
       P_AJ (5, 2) = AJ (2, 1) * AJ (3, 2) + AJ (2, 2) * AJ (3, 1) 
       P_AJ (6, 2) = AJ (3, 1) * AJ (3, 2) * 2.d0 
       P_AJ (1, 3) = AJ (1, 1) * AJ (1, 3) * 2.d0 
       P_AJ (2, 3) = AJ (1, 1) * AJ (2, 3) + AJ (1, 3) * AJ (2, 1) 
       P_AJ (3, 3) = AJ (1, 1) * AJ (3, 3) + AJ (1, 3) * AJ (3, 1) 
       P_AJ (4, 3) = AJ (2, 1) * AJ (2, 3) * 2.d0 
       P_AJ (5, 3) = AJ (2, 1) * AJ (3, 3) + AJ (2, 3) * AJ (3, 1) 
       P_AJ (6, 3) = AJ (3, 1) * AJ (3, 3) * 2.d0 
       P_AJ (1, 4) = AJ (1, 2) **2 
       P_AJ (2, 4) = AJ (1, 2) * AJ (2, 2) 
       P_AJ (3, 4) = AJ (1, 2) * AJ (3, 2) 
       P_AJ (4, 4) = AJ (2, 2) **2 
       P_AJ (5, 4) = AJ (2, 2) * AJ (3, 2) 
       P_AJ (6, 4) = AJ (3, 2) **2 
       P_AJ (1, 5) = AJ (1, 2) * AJ (1, 3) * 2.d0 
       P_AJ (2, 5) = AJ (1, 2) * AJ (2, 3) + AJ (1, 3) * AJ (2, 2) 
       P_AJ (3, 5) = AJ (1, 2) * AJ (3, 3) + AJ (1, 3) * AJ (3, 2) 
       P_AJ (4, 5) = AJ (2, 2) * AJ (2, 3) * 2.d0 
       P_AJ (5, 5) = AJ (2, 2) * AJ (3, 3) + AJ (2, 3) * AJ (3, 2) 
       P_AJ (6, 5) = AJ (3, 2) * AJ (3, 3) * 2.d0 
       P_AJ (1, 6) = AJ (1, 3) **2 
       P_AJ (2, 6) = AJ (1, 3) * AJ (2, 3) 
       P_AJ (3, 6) = AJ (1, 3) * AJ (3, 3) 
       P_AJ (4, 6) = AJ (2, 3) **2 
       P_AJ (5, 6) = AJ (2, 3) * AJ (3, 3) 
       P_AJ (6, 6) = AJ (3, 3) **2
     CASE DEFAULT
       PRINT *, 'JACOBIAN_PRODUCTS: UNKNOWN N_2_DER = ', N_2_DER
       STOP     'JACOBIAN_PRODUCTS: UNKNOWN N_2_DER'
   END SELECT
END SUBROUTINE JACOBIAN_PRODUCTS 

subroutine Sort_Reals (lines, database)
!------------------------------------------------------------
! Bubble Sort of (changed) Real Database
!------------------------------------------------------------
use System_Constants ! for DP
implicit none
 integer, intent(in)    :: lines             ! number of records 
 !b real,    intent(inout) :: database (lines)  ! records in database
 real(DP), intent(inout) :: database (lines)  ! records in database

 integer :: swaps_Made       ! number of swaps made in one pass
 integer :: count            ! loop variable
 !b real    :: temp             ! temporary holder for making swap
 real(DP):: temp             ! temporary holder for making swap

   do ! Repeat this loop forever... (until we break out of it)
     swaps_Made = 0             ! Initially, we've made no swaps
     ! Make one pass of the bubble sort algorithm
     do count = 1, (lines - 1)
       ! If item is greater than the one after it, swap them
       if ( database (count) > database (count + 1) ) then
              temp                 = database (count)
              database (count)     = database (count + 1)
              database (count + 1) = temp
              swaps_Made           = swaps_Made + 1
        end if
      end do
      ! If we made no swaps, break out of the loop.
      if ( swaps_Made == 0 ) exit ! do count swaps
    end do
end subroutine Sort_Reals

subroutine Real_Vector_Order (lines, database, order)
!------------------------------------------------------------
! Ordered Bubble Sort of (Unchanged) Real Vector
!------------------------------------------------------------
use System_Constants ! for DP
implicit none
 integer,  intent(in)  :: lines             ! number of records
 real(dp), intent(in)  :: database (lines)  ! records in database
 integer,  intent(out) :: order    (lines)  ! the order array

 integer :: swaps_Made       ! number of swaps made in one pass
 integer :: count            ! loop variable
 integer :: temp             ! temporary holder for making swap

   order = (/ (count, count = 1, lines) /) ! default order
   do ! Repeat this loop forever... (until we break out of it)
     swaps_Made = 0             ! Initially, we've made no swaps
     ! Make one pass of the bubble sort algorithm
     do count = 1, (lines - 1)
       ! If item is greater than the one after it, swap them
       if ( database (order (count)) >     &
            database (order (count + 1)) ) then
              temp              = order (count)
              order (count)     = order (count + 1)
              order (count + 1) = temp
              swaps_Made        = swaps_Made + 1
        end if
      end do
      ! If we made no swaps, break out of the loop.
      if ( swaps_Made == 0 ) exit ! do count swaps
    end do
end subroutine Real_Vector_Order

subroutine Sort_string (lines, database)
!------------------------------------------------------------
! Bubble Sort of (Changed) String Database
!------------------------------------------------------------
implicit none

 integer,          intent(in)     :: lines            ! input size
 character(len=*), intent(inout)  :: database (lines) ! records

 character (len = len(database (1))) :: temp ! swap holder 
 integer :: swaps_Made      ! number of swaps in a pass
 integer :: count           ! loop variable

   interface ! to_lower
     function  to_lower (string)  result (new_string) 
       character (len = *), intent(in) :: string     
       character (len = len(string))   :: new_string  
     end function to_lower
   end interface ! to_lower

   do ! Repeat this loop forever... (until we break out of it)
     swaps_Made = 0            ! Initially, we've made no swaps
     ! Make one pass of the bubble sort algorithm
     do count = 1, (lines - 1)
       ! If the element is greater than the one after it, swap them
       if ( LGT (to_lower (database (count    )),   &
                 to_lower (database (count + 1))) ) then
           temp                 = database (count    )
           database (count    ) = database (count + 1)
           database (count + 1) = temp
           swaps_Made           = swaps_Made + 1
       end if
     end do
     ! If we made no swaps, berak out of the loop.
     if ( swaps_Made == 0) exit ! do count swaps
   end do
end subroutine Sort_string

function  to_lower (string)  result (new_string) 
!------------------------------------------------------------
!      Convert a string or character to lower case
!       (valid for ASCII or EBCDIC processors)
!------------------------------------------------------------
implicit none
 character (len = *), intent(in) :: string      ! unknown length
 character (len = len(string))   :: new_string  ! same length 
 character (len=26), parameter   ::  &
           UPPER = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ',  &
           lower = 'abcdefghijklmnopqrstuvwxyz' 
 integer :: k    ! loop counter
 integer :: loc  ! position in alphabet

   new_string = string       ! copy everything
   do k = 1, len(string)     ! to change letters
      loc = index ( UPPER, string(k:k))
      if (loc /= 0 ) new_string(k:k) = lower(loc:loc)
   end do ! over string characters
end function to_lower

subroutine Integer_Sort (lines, database, order)
!------------------------------------------------------------
! Ordered Bubble Sort of (Unchanged) Integer Database
!------------------------------------------------------------
implicit none
 integer, intent(in)  :: lines             ! number of records 
 integer, intent(in)  :: database (lines)  ! records in database
 integer, intent(out) :: order    (lines)  ! the order array

 integer :: swaps_Made       ! number of swaps made in one pass
 integer :: count            ! loop variable
 integer :: temp             ! temporary holder for making swap

   order = (/ (count, count = 1, lines) /) ! default order
   do ! Repeat this loop forever... (until we break out of it)
     swaps_Made = 0             ! Initially, we've made no swaps
     ! Make one pass of the bubble sort algorithm
     do count = 1, (lines - 1)
       ! If item is greater than the one after it, swap them
       if ( database (order (count)) >     &
            database (order (count + 1)) ) then
              temp              = order (count)
              order (count)     = order (count + 1)
              order (count + 1) = temp
              swaps_Made        = swaps_Made + 1
        end if
      end do
      ! If we made no swaps, break out of the loop.
      if ( swaps_Made == 0 ) exit ! do count swaps
    end do
end subroutine Integer_Sort

FUNCTION GET_ABS_NORM (A) RESULT (VALUE)
Use System_Constants    ! for DP
 IMPLICIT NONE
  REAL(DP), INTENT (IN) :: A (:, :)
  REAL(DP)              :: VALUE
    VALUE = SUM ( ABS ( A ) )
END FUNCTION GET_ABS_NORM

FUNCTION GET_INF_NORM (A) RESULT (VALUE)
Use System_Constants    ! for DP
 IMPLICIT NONE
  REAL(DP), INTENT (IN) :: A (:, :)
  REAL(DP)              :: VALUE, TEST
  INTEGER :: I
    VALUE = 0.d0
    DO I = 1, SIZE (A, 1)
      TEST = SUM ( ABS ( A (I, :) ) )
      IF ( TEST > VALUE ) VALUE = TEST
    END DO ! over rows
END FUNCTION GET_INF_NORM

FUNCTION GET_MAX_NORM (A) RESULT (VALUE)
Use System_Constants    ! for DP
 IMPLICIT NONE
  REAL(DP), INTENT (IN) :: A (:, :)
  REAL(DP)              :: VALUE
    VALUE =  MAXVAL ( ABS (A) )
END FUNCTION GET_MAX_NORM

FUNCTION GET_SQ_NORM (A) RESULT (VALUE)
Use System_Constants    ! for DP
 IMPLICIT NONE
  REAL(DP), INTENT (IN) :: A (:, :)
  REAL(DP)              :: VALUE
    VALUE = SQRT ( SUM ( A **2) ) / SIZE (A, 1)
END FUNCTION GET_SQ_NORM

FUNCTION GET_1_NORM (A) RESULT (VALUE)
Use System_Constants    ! for DP
 IMPLICIT NONE
  REAL(DP), INTENT (IN) :: A (:, :)
  REAL(DP)              :: VALUE, TEST
  INTEGER :: I
    VALUE = 0.d0
    DO I = 1, SIZE (A, 2)
      TEST = SUM ( ABS ( A (:, I) ) )
      IF ( TEST > VALUE ) VALUE = TEST
    END DO ! over rows
END FUNCTION GET_1_NORM
