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

SUBROUTINE DUMMY
! dummy  file
END SUBROUTINE DUMMY

SUBROUTINE FORM_FACING_ELEMENTS (N_ELEMS, NOD_PER_EL, MAX_NP,  &
                                 L_FIRST, L_LAST, NODES, NEEDS,&
                                 MAX_FACES, L_TO_L_NEIGH)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *      
! BUILD LIST OF FACING ELEMENTS ADJACENT TO OTHER ELEMENTS        
! (THE LIST IS UNORDERED. ZEROS DENOTE SURFACE LOCATIONS.)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *     
! Draft 2, revised and tested June 11, 1999, JEA
 IMPLICIT NONE
 INTEGER, INTENT(IN)  :: N_ELEMS, NOD_PER_EL, MAX_NP
 INTEGER, INTENT(IN)  :: NEEDS, MAX_FACES
 INTEGER, INTENT(IN)  :: L_FIRST (MAX_NP), L_LAST (MAX_NP)
 INTEGER, INTENT(IN)  :: NODES (N_ELEMS, NOD_PER_EL)
 INTEGER, INTENT(OUT) :: L_TO_L_NEIGH (MAX_FACES, N_ELEMS)

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

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

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

  MAIN : DO IE = 1, N_ELEMS      ! LOOP OVER ALL ELEMENTS 
    FOUND = COUNT ( L_TO_L_NEIGH (:, IE) /= 0 )

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

!        ESTABLISH RANGE OF POSSIBLE ELEMENT NEIGHBORS    
    L_START = N_ELEMS+1 ; L_STOP = 0 
    DO IN = 1, NOD_PER_EL 
      L_START = MIN (L_START, L_FIRST (ELEM_NODES (IN)) )
      L_STOP  = MAX (L_STOP,  L_LAST  (ELEM_NODES (IN)) )
    END DO 
    L_START = MAX (L_START, (IE+1))

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

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

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

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

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

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

                NEXT  = COUNT ( L_TO_L_NEIGH(:, L_TEST) /= 0 )
                WHERE = NEXT+1 ! OR ORDER THE NEIGHBOR FACE
                L_TO_L_NEIGH (WHERE, L_TEST) = IE    ! 2 OF 2
                IF ( MAX_FACES == FOUND ) CYCLE MAIN  ! IE DONE
                CYCLE RANGE ! this L_TEST element search loop
              END IF ! NUMBER NEEDED
            END IF
        END DO LOCAL ! over in
      END DO RANGE ! over candidate element L_TEST
    END IF ! a possible candidate

!    LOCATE MISSING SURFACE NUMBERS, IF ANY, FOR IE
    IF ( ANY ( L_TO_L_NEIGH (:, IE) == 0 ) ) THEN
!     DO SURFACE SEARCH HERE
    END IF ! SURFACE SEARCH

!   NOW ANY ( L_TO_L_NEIGH (:, IE) == 0 ) NEEDS A NO NEIGHBORS OR 
!       MISSING SURFACE WARNING MESSAGE FOR ELEMENT IE
!       (IFF ALL ELEMENTS HAVE THE SAME MAX_FACES.)
  END DO MAIN ! over all elements
END SUBROUTINE FORM_FACING_ELEMENTS                        

SUBROUTINE FORM_FACING_OR_SURFS (N_ELEMS, NOD_PER_EL, MAX_NP,     &
                                 L_FIRST, L_LAST, NODES, NEEDS,   &
                                 MAX_FACES, L_TO_L_NEIGH, N_SURFS,&
                                 NOD_PER_SURF, NODES_SURF)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *      
!       COUNT NUMBER OF FACING ELEMENTS OR SURFACES
!              ADJACENT TO OTHER ELEMENTS        
!  (INSERT NEGATIVE SURFACE CODES FIRST INTO L_TO_L_NEIGH)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *     
 IMPLICIT NONE
 INTEGER, INTENT(IN)  :: N_ELEMS, NOD_PER_EL, MAX_NP
 INTEGER, INTENT(IN)  :: N_SURFS, NOD_PER_SURF 
 INTEGER, INTENT(IN)  :: NEEDS, MAX_FACES
 INTEGER, INTENT(IN)  :: L_FIRST (MAX_NP), L_LAST (MAX_NP)
 INTEGER, INTENT(IN)  :: NODES      (N_ELEMS, NOD_PER_EL)
 INTEGER, INTENT(IN)  :: NODES_SURF (N_SURFS, NOD_PER_SURF+1)
 INTEGER, INTENT(OUT) :: L_TO_L_NEIGH (MAX_FACES, N_ELEMS)

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

! ELEM_NODES     = INCIDENCES ARRAY FOR A SINGLE ELEMENT     
! KOUNT          = NUMBER OF COMMON NODES
! L_FIRST (I)    = ELEMENT WHERE NODE I FIRST APPEARS              
! L_LAST (I)     = ELEMENT WHERE NODE I LAST APPEARS             
! L_TO_L_NEIGH (I) = NUMBER OF ELEM NEIGHBORS OF ELEMENT I
! MAX_FACES      = MAX NUMBER OF FACING ELEMENTS
! MAX_NP         = TOTAL NUMBER OF NODES                          
! NEEDS          = NUMBER OF COMMON NODES TO BE A NEIGHBOR
! NEIG_NODES     = INCIDENCES ARRAY FOR A POSSIBLE NEIGHBOR
! N_ELEMS        = TOTAL NUMBER OF ELEMENTS
! N_SURFS        = NUMBER OF CODED SURFACES IN SYSTEM
! NOD_PER_EL     = NUMBER OF NODES PER ELEMENT               
! NOD_PER_SURF   = NUMBER OF NODES PER SURFACE
! NODES          = SYSTEM ARRAY OF INCIDENCES OF ALL ELEMENTS      
! NODES_SURF     = CODE NUMBER AND INCIDENCES OF ALL SURFACES
! SURF_NODES     = INCIDENCES ARRAY FOR A SINGLE SURFACE     

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

  SURF : DO IE = 1, N_SURFS      ! LOOP OVER ALL SURFENTS 

!        EXTRACT INCIDENCES LIST FOR SURFENT IE   
    SURF_NODES = NODES_SURF (IE, 2:(NOD_PER_SURF+1))

!        ESTABLISH RANGE OF POSSIBLE ELEMENT NEIGHBORS    
    L_START = N_ELEMS+1 ; L_STOP = 0 
    DO IN = 1, NOD_PER_EL 
      L_START = MIN (L_START, L_FIRST (SURF_NODES (IN)) )
      L_STOP  = MAX (L_STOP,  L_LAST  (SURF_NODES (IN)) )
    END DO 

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

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

!            COMPARE WITH INCIDENCES OF SURFENT IE             
            IF ( ANY ( SURF_NODES == N_TEST ) ) THEN
              KOUNT = KOUNT + 1         ! SHARED NODE COUNT
              IF ( KOUNT == NEED ) THEN ! SURF ON ELEM FOUND
                FOUND = FOUND + 1       ! INSERT THE SURF
                !b L_TO_L_NEIGH (FOUND, IE) = L_TEST     
                L_TO_L_NEIGH (FOUND, L_TEST) = NODES_SURF (IE, 1)
                IF ( MAX_FACES == FOUND ) CYCLE SURF  ! IE DONE
                CYCLE CHECK ! this L_TEST element search loop
              END IF ! NUMBER NEEDED
            END IF
        END DO L_2_S ! over in
        !b END IF not IE
      END DO CHECK ! over candidate element L_TEST
    END IF ! a possible candidate
  END DO SURF ! over all elements

  MAIN : DO IE = 1, N_ELEMS      ! LOOP OVER ALL ELEMENTS 
    FOUND = COUNT ( L_TO_L_NEIGH (:, IE) /= 0 )
!b print *, 'FOUND ', FOUND

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

!        ESTABLISH RANGE OF POSSIBLE ELEMENT NEIGHBORS    
    L_START = N_ELEMS+1 ; L_STOP = 0 
    DO IN = 1, NOD_PER_EL 
!     IF ( L_START > L_FIRST (ELEM_NODES (IN) ) ) &
!          L_START = L_FIRST (ELEM_NODES (IN) )
!     IF ( L_STOP  < L_LAST  (ELEM_NODES (IN) ) ) &
!          L_STOP  = L_LAST  (ELEM_NODES (IN) )
      L_START = MIN (L_START, L_FIRST (ELEM_NODES (IN)) )
      L_STOP  = MAX (L_STOP,  L_LAST  (ELEM_NODES (IN)) )
    END DO 
    L_START = MAX (L_START, (IE+1))

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

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

!            COMPARE WITH INCIDENCES OF ELEMENT IE             
            IF ( ANY ( ELEM_NODES == N_TEST ) ) THEN
              KOUNT = KOUNT + 1         ! SHARED NODE COUNT
!b print *, 'IE, L_TEST, KOUNT       ', IE, L_TEST, KOUNT      
              IF ( KOUNT == NEED ) THEN ! NEIGHBOR PAIR FOUND
                FOUND = FOUND + 1       ! INSERT THE PAIR
                L_TO_L_NEIGH (FOUND, IE) = L_TEST     ! 1 OF 2
                NEXT = COUNT ( L_TO_L_NEIGH(:, L_TEST) /= 0 )
!b print *, 'IE, L_TEST, FOUND, NEXT ', IE, L_TEST, FOUND, NEXT
                L_TO_L_NEIGH (NEXT+1, L_TEST) = IE    ! 2 OF 2
                IF ( MAX_FACES == FOUND ) CYCLE MAIN  ! IE DONE
                CYCLE RANGE ! this L_TEST element search loop
              END IF ! NUMBER NEEDED
            END IF
        END DO LOCAL ! over in
        !b END IF not IE
      END DO RANGE ! over candidate element L_TEST
    END IF ! a possible candidate
  END DO MAIN ! over all elements
END SUBROUTINE FORM_FACING_OR_SURFS                        

SUBROUTINE TRAPEZOIDAL_TIME_RULE (I_DIAG, S_SKY, M_SKY, S_SKY_LOWER)
!                          not used
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!   APPLY GENERALIZED TRAPEZOIDAL RULE TO REFORM TWO SQUARE MATRICES
!                     M_SKY * D_DOT + S_SKY * D = C
!                 INTO THE TWO USED IN TIME INTEGRATION
! (M_SKY/T_STEP + TS_B*S_SKY) D_1 = (M_SKY/T_STEP - (1-TS_B)*S_SKY) D_0
!                                 + TS_B*C_1 + (1-TS_B)*C_0
!                                 OR
!         M_SKY D_1 = S_SKY D_0 + TS_B*C_1 + (1-TS_B)*C_0
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants  ! for N_COEFF, N_D_FRE, NOT_SYM, TS_B, TIME_STEP
 IMPLICIT NONE
 INTEGER,  INTENT(IN)    :: I_DIAG (N_D_FRE)
 REAL(DP), INTENT(INOUT) :: M_SKY (N_COEFF)
 REAL(DP), INTENT(INOUT) :: S_SKY (N_COEFF)
 REAL(DP), INTENT(INOUT) :: S_SKY_LOWER (NOT_SYM)

 INTEGER :: I
! N_D_FRE         = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM
! N_COEFF         = NUMBER OF NON-ZERO COEFFICIENTS IN S_SKY
! NOT_SYM         = NUMBER OF NON-ZERO COEFF IN S_SKY_LOWER,
!                   = 0 OR N_COEFF
! M_SKY           = SYS. EQ. SQ. MASS MATRIX (SYMMETRIC)
! S_SKY           = SYS. EQ. SQ. MATRIX
! S_SKY_LOWER     = SYS. EQ. SQ. MATRIX

  IF ( I_BUG == 1 ) PRINT *, 'Entering TRAPEZOIDAL_TIME_RULE'

  M_SKY = M_SKY / TIME_STEP + TS_B * S_SKY
  S_SKY = M_SKY - S_SKY
END SUBROUTINE TRAPEZOIDAL_TIME_RULE

SUBROUTINE SKY_BC_NULL (N_RES_EQ, NDEX_C_EQ, M_SKY, I_DIAG)
!b SUBROUTINE SKY_BC_NULL (N_RES_EQ, NDEX_C_EQ, M_SKY, CC, I_DIAG)
!                          not used
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!   FOR TRANSIENT OR DYNAMIC ANALYSIS, ZERO ALL ROWS AND COLUMNS 
!     IN THE SYSTEM MASS MATRIX WHERE CONSTRAINTS ARE APPLIED
!        (THIS ASSUMES VELOCITES OF CONSTRAINTS ARE ZERO)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants
 IMPLICIT NONE
 INTEGER,  INTENT(IN)    :: N_RES_EQ  (MAX_ACT),             &
                            NDEX_C_EQ (MAX_ACT, N_CEQ),      &
                            I_DIAG (N_D_FRE)
 !b REAL(DP), INTENT(INOUT) :: M_SKY (N_COEFF), CC (N_D_FRE)
 REAL(DP), INTENT(INOUT) :: M_SKY (N_COEFF)
 
 INTEGER  :: IC, IEQ, NEQ, NTEST, K
 !b REAL(DP) :: ZERO (MAX_ACT)               ! NULL BC COEFF
 REAL(DP) :: ZERO (MAX_ACT), CC (N_D_FRE) ! BOTH ZERO
                                      
! N_D_FRE         = TOTAL NUMBER OF SYSTEM DEGREES OF FREEDOM
! N_COEFF         = NUMBER OF NON-ZERO COEFFICIENTS IN M_SKY
! M_SKY           = SYS. EQ. SQ. MASS MATRIX (SYMMETRIC)
! CC              = SYS. EQ. COL. MATRIX
! N_G_DOF         = NO. PARAMETERS PER NODE
! N_RES_EQ        = NO OF CONSTRAINT EQS OF EACH TYPE
! NDEX_C_EQ (I,J) = CONSTRAINT EQS DOF NO I FOR EQ J
! N_CEQ           = NUMBER OF CONSTRAINT EQUATIONS

  ZERO (MAX_ACT) = 0.d0 ; CC (N_D_FRE) = 0.d0
  IF ( I_BUG == 1 ) THEN
    PRINT *, 'Entering SKY_BC_NULL'
  END IF
  IF ( STATIC ) THEN
    PRINT *, 'WARNING: SKY_BC_NULL CALLED IN STATIC RUN'
    N_WARN = N_WARN + 1
  END IF
 
!  WORK BACKWARD FROM LARGEST CONSTRAINT TYPE (DO TYPE 1 LAST)
  IF ( SUM(N_RES_EQ) < 1) THEN
    PRINT *, 'SKY_BC_NULL: WARNING, NO RESTRAINTS APPLIED'
    N_WARN = N_WARN + 1 ! INCREMENT WARNING
  END IF
  IEQ = SUM ( N_RES_EQ ) + 1  
  DO IC = MAX_ACT, 1, -1
    NTEST = N_RES_EQ (IC)
    IF ( NTEST > 0 ) THEN ! Apply constraint of type IC
      IF ( IC == 1 ) THEN
!-->         TYPE 1     D(L1) = C1
        DO NEQ = 1, NTEST
          IEQ = IEQ - 1 
          CALL SKY_TYPE_1_SYM (NDEX_C_EQ (1, IEQ), ZERO,     &
                               M_SKY, CC, I_DIAG)
!           ALSO ZERO DIAGONAL TERM
          M_SKY ( I_DIAG (NDEX_C_EQ (1, IEQ)) ) = 0.d0
        END DO ! FOR EACH SINGLE DOF

      ELSE ! IC >= 2
!-->         TYPE 2   D(L1) + C1*D(L2) = C2
!-->         TYPE 3   D(L1) + C1*D(L2) + C2*D(L3) = C3, ETC.
        DO NEQ = 1, NTEST
          IEQ = IEQ - 1 
          DO K = 1, IC
            CALL SKY_TYPE_1_SYM (NDEX_C_EQ (K, IEQ), ZERO,     &
                                 M_SKY, CC, I_DIAG)
!           ALSO ZERO DIAGONAL TERM
            M_SKY ( I_DIAG (NDEX_C_EQ (K, IEQ)) ) = 0.d0
          END DO ! FOR EACH COUPLED DOF
        END DO ! FOR COUPLED DOF SETS
      END IF
    END IF
  END DO
END SUBROUTINE SKY_BC_NULL

SUBROUTINE SET_L_SHAPE_FACE_NEEDS
Use System_Constants
Use Elem_Type_Data
   SELECT CASE (L_SHAPE) ! SHAPE OF ELEMENT
     CASE (1) ; NEEDS = 1 ; MAX_FACES = 2  ! LINE ELEMENT
     CASE (2) ; NEEDS = 2 ; MAX_FACES = 3  ! TRIANGULAR ELEMENT
     CASE (3) ; NEEDS = 2 ; MAX_FACES = 4  ! QUADRILATERAL ELEMENT
     CASE (4) ; NEEDS = 4 ; MAX_FACES = 6  ! HEXAHEDRA ELEMENT
     CASE (5) ; NEEDS = 3 ; MAX_FACES = 4  ! TETRAHEDRA ELEMENT
     CASE (6) ; NEEDS = 3 ; MAX_FACES = 5  ! WEDGE ELEMENT
     CASE DEFAULT ; NEEDS = 1 ; MAX_FACES = 6
   END SELECT
END SUBROUTINE SET_L_SHAPE_FACE_NEEDS

SUBROUTINE SET_EL_TYPE_SHAPE_FACE_NEEDS (LT_SHAP, LT_NEEDS, LT_FACES)
 IMPLICIT NONE
 INTEGER, INTENT (IN)  :: LT_SHAP
 INTEGER, INTENT (OUT) :: LT_NEEDS, LT_FACES
   SELECT CASE (LT_SHAP) ! SHAPE OF THIS ELEMENT TYPE
     CASE (1) ; LT_NEEDS = 1 ; LT_FACES = 2  ! LILT_NE ELEMENT
     CASE (2) ; LT_NEEDS = 2 ; LT_FACES = 3  ! TRIANGULAR ELEMENT
     CASE (3) ; LT_NEEDS = 2 ; LT_FACES = 4  ! QUADRILATERAL ELEMENT
     CASE (4) ; LT_NEEDS = 4 ; LT_FACES = 6  ! HEXAHEDRA ELEMENT
     CASE (5) ; LT_NEEDS = 3 ; LT_FACES = 4  ! TETRAHEDRA ELEMENT
     CASE (6) ; LT_NEEDS = 3 ; LT_FACES = 5  ! WEDGE ELEMENT
     CASE DEFAULT ; LT_NEEDS = 1 ; LT_FACES = 6
   END SELECT
END SUBROUTINE SET_EL_TYPE_SHAPE_FACE_NEEDS

SUBROUTINE SAVE_QP_SOURCE (N_FILE, XYZ, SOURCE)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
! SAVE QUADRATURE POINT COORDINATES AND SOURCE TERM(S)
! * * * * * * * * * * * * * * * * * * * * * * * * * * * *
Use System_Constants  ! for N_SPACE
 IMPLICIT NONE
  INTEGER,   INTENT(IN) :: N_FILE
  REAL(DP),  INTENT(IN) :: XYZ (N_SPACE), SOURCE

   WRITE (N_FILE, 40)  XYZ, SOURCE
   40 FORMAT ( (10(1PE13.5)) )
END SUBROUTINE SAVE_QP_SOURCE

SUBROUTINE OPEN_QP_SOURCE (N_FILE)
! OPEN TO SAVE QUADRATURE POINT COORDINATES AND SOURCE TERM
Use System_Constants  ! for N_WARN
 IMPLICIT NONE
  INTEGER, INTENT(IN) :: N_FILE ; INTEGER :: IO_1

   if ( N_FILE > 0 )  then !b
     OPEN (N_FILE, FILE='el_qp_xyz_sources.tmp', ACTION='WRITE', &
       STATUS='REPLACE', IOSTAT = IO_1) !b
     IF ( IO_1 > 0 ) THEN !b
       PRINT *, 'NOTE: FILE OPEN', N_FILE, ' FAILED' !b
     ELSE !b
       PRINT *, 'NOTE: el_qp_xyz_sources.tmp OPEN AS ', N_FILE !b
     END IF ! open failed !b
   else
    N_WARN = N_WARN + 1
    PRINT *, 'WARNING: OPEN_QP_SOURCE, el_qp_xyz_sources.tmp FAILED'
   end if ! N_FILE 
END SUBROUTINE OPEN_QP_SOURCE 

! ============================================================

SUBROUTINE ELEM_CENTER_GEOMETRY (CENTER, RADIAL)
 Use System_Constants ! for N_SPACE
 Use Elem_Type_Data   ! for LT_N, COORD (LT_N, N_SPACE)
 IMPLICIT NONE
  REAL(DP), INTENT (OUT) :: CENTER (N_SPACE)       ! average of nodes
  REAL(DP), INTENT (OUT) :: RADIAL (LT_N, N_SPACE) ! relative positions
  INTEGER :: J ! loops

  CENTER = SUM ( COORD, DIM=1 ) / LT_N   ! average nodal positions
  DO J = 1, N_SPACE
    RADIAL (:, J) = COORD (:, J) - CENTER (J) ! From center
  END DO
END SUBROUTINE ELEM_CENTER_GEOMETRY 

SUBROUTINE GET_RADIALS_FROM_CENTER (CENTER, RADIAL)
 Use System_Constants ! for N_SPACE
 Use Elem_Type_Data   ! for LT_N, COORD (LT_N, N_SPACE)
 IMPLICIT NONE
  REAL(DP), INTENT (IN ) :: CENTER (N_SPACE)       ! average of nodes
  REAL(DP), INTENT (OUT) :: RADIAL (LT_N, N_SPACE) ! relative positions
  INTEGER :: J ! loops

  DO J = 1, N_SPACE
    RADIAL (:, J) = COORD (:, J) - CENTER (J) ! From center
  END DO
END SUBROUTINE GET_RADIALS_FROM_CENTER 

SUBROUTINE UPWIND_DIST_AND_LOGIC (RADIAL, UNIT_V, DOWN, &
                                   LENGTH, IS_DOWNWIND)
 Use System_Constants ! for N_SPACE
 Use Elem_Type_Data   ! for LT_N
 IMPLICIT NONE
  REAL(DP), INTENT (IN ) :: RADIAL (LT_N, N_SPACE) ! relative positions
  REAL(DP), INTENT (IN ) :: UNIT_V (N_SPACE)       ! vel unit vector
  REAL(DP), INTENT (OUT) :: DOWN   (LT_N), LENGTH  ! downwind if > 0
  LOGICAL,  INTENT (OUT) :: IS_DOWNWIND (LT_N)     ! true if downwind 
  REAL(DP) :: L_POS, L_NEG 
  INTEGER :: K
  LOGICAL :: LOGIC (LT_N)   ! debug
  
!       L_POS  = MAXVAL (DOWN, IS_DOWNWIND)
!       L_NEG  = ABS (MINVAL (DOWN, .NOT. IS_DOWNWIND))

   IF ( ALL (UNIT_V == 0.d0) ) STOP 'Zero unit vector. u_d_a_l'
   IF ( ALL (RADIAL == 0.d0) ) STOP 'Zero radial array. u_d_a_l'

!          IDENTIFY DOWNWIND DISTANCES, AND LOGIC
   DOWN  = MATMUL ( RADIAL, UNIT_V )    ! POSITIVE & NEGATIVE
   IF ( ALL (DOWN == 0.d0) ) STOP 'Zero DW vector. '

   L_POS = DOWN (1) ; L_NEG = DOWN (1) 
   DO K = 1, LT_N
     IF ( DOWN (K) >= 0.d0 ) THEN
       IS_DOWNWIND (K) = .TRUE.
     ELSE
       IS_DOWNWIND (K) = .FALSE.
     END IF
     IF ( L_POS < DOWN (K) ) L_POS = DOWN (K)
     IF ( L_NEG > DOWN (K) ) L_NEG = DOWN (K)
   END DO

!        GET MAXIMUM STREAMLINE DISTANCES
   L_NEG = ABS (L_NEG) ; LENGTH = L_NEG + L_POS

!       WHERE ( DOWN >= 0.d0 ) ; IS_DOWNWIND = .TRUE.
!       ELSEWHERE              ; IS_DOWNWIND = .FALSE.  ; END WHERE
   WHERE ( DOWN >= 0.d0 ) ; LOGIC = .TRUE.
   ELSEWHERE              ; LOGIC = .FALSE.  ; END WHERE

if ( any (IS_DOWNWIND .NEQV. LOGIC) ) then
  print *,'ERROR inconsistent logic for element ', THIS_EL
  print *,'IS_DW ', IS_DOWNWIND
  print *,'LOGIC ', LOGIC 
  STOP 'UPWIND_DIST_AND_LOGIC: ERROR in logic'
end if

if ( all(IS_DOWNWIND) ) then
  print *,'ERROR all nodes downwind for element ', THIS_EL
  print *,'RADIAL '; call rprint (RADIAL, LT_N, N_SPACE, 1)
  print *,'UNIT_V ', UNIT_V
  STOP 'UPWIND_DIST_AND_LOGIC: ERROR all nodes downwind'
end if

!if ( THIS_EL == 1 ) then
!  print *,'============== UPWIND_DIST_AND_LOGIC ================'
!  print *,'RADIAL '; call rprint (RADIAL, LT_N, N_SPACE, 1)
!  print *,'UNIT_V ', UNIT_V
!  print *,'DOWN ', DOWN
!  print *,'L_NEG, L_POS, LENGTH ', L_NEG, L_POS, LENGTH
!  print *,'IS_DW ', IS_DOWNWIND
!  print *,'LOGIC ', LOGIC 
!end if

    WHERE ( IS_DOWNWIND ) ; DOWN = DOWN + L_NEG
    ELSEWHERE             ; DOWN = ABS (DOWN) + L_POS ; END WHERE

!if ( THIS_EL == 1 ) then
!  print *,'DOWN ', DOWN
!end if
END SUBROUTINE UPWIND_DIST_AND_LOGIC 

SUBROUTINE GET_MAX_DOWNWIND_DIST (DOWN, DOWN_MAX)
!b SUBROUTINE GET_NODAL_DOWNWIND_DIST (DOWN, DOWN_NODAL) -> MAX
 Use System_Constants ! for N_SPACE
 Use Elem_Type_Data   ! for LT_N
 IMPLICIT NONE
  REAL(DP), INTENT (IN ) :: DOWN     (LT_N) ! downwind wrt center
  REAL(DP), INTENT (OUT) :: DOWN_MAX (LT_N) ! downwind total 
  REAL(DP) :: L_POS, L_NEG 
  
!        GET AVERAGE STREAMLINE DISTANCES
        L_POS  = MAXVAL (DOWN)
        L_NEG  = ABS (MINVAL (DOWN))
if ( L_POS <= 0.d0 ) then
  print *,'ERROR, no DW dist for element ', THIS_EL
  print *,'DOWN dist ', DOWN
  stop    'GET_MAX_DOWNWIND_DIST: ERROR, no DW dist for element '
end if
        
!        GET TOTAL DOWNWIND MAX DISTANCES
        WHERE ( DOWN >= 0.d0 ) 
          DOWN_MAX = DOWN + L_NEG
        ELSEWHERE             
          DOWN_MAX = ABS (DOWN) + L_POS ; END WHERE
!if ( THIS_EL == 1 ) then
!  print *,'DOWN     ', DOWN    
!  print *,'L_NEG, L_POS ',  L_NEG, L_POS
!  print *,'DOWN_MAX ', DOWN_MAX
!end if
END SUBROUTINE GET_MAX_DOWNWIND_DIST 

SUBROUTINE GET_AVE_DOWNWIND_DIST (DOWN, DOWN_AVE)
!b SUBROUTINE GET_AVE_DOWNWIND_DIST (DOWN, IS_DOWNWIND, DOWN_AVE)
 Use System_Constants ! for N_SPACE
 Use Elem_Type_Data   ! for LT_N
 IMPLICIT NONE
  REAL(DP), INTENT (IN ) :: DOWN        (LT_N)  ! downwind if > 0
  !b LOGICAL,  INTENT (IN ) :: IS_DOWNWIND (LT_N)  ! true if downwind 
  REAL(DP), INTENT (OUT) :: DOWN_AVE    (LT_N)  ! downwind total      
  LOGICAL  :: IS_DOWNWIND (LT_N)  ! true if downwind 
  LOGICAL  :: LOGIC (LT_N)
  REAL(DP) :: AVE_POS, AVE_NEG 
  INTEGER  :: C_DW, C_UW, K  ! down or up wind
  
!if ( THIS_EL == 1 ) then
!  print *,'GET_AVE_DOWNWIND_DIST =======', THIS_EL !b , IS_DOWNWIND
!  print *,'DOWN ', DOWN
!end if

AVE_POS = 0 ; AVE_NEG = 0 ; C_DW = 0 ; C_UW = 0
DO K = 1, LT_N
  IF ( DOWN (K) >= 0.d0 ) THEN
    AVE_POS = AVE_POS + DOWN (K)
    C_DW    = C_DW + 1
    LOGIC (K) = .TRUE.
  ELSE
    AVE_NEG = AVE_NEG + DOWN (K)
    C_UW    = C_UW + 1
    LOGIC (K) = .FALSE.
  END IF
  !b IF ( LOGIC (K) .NEQV. IS_DOWNWIND (K) ) PRINT *, &
    !b 'ERROR, Logic mismatch at node, el: ',         &
     !b K, THIS_EL, IS_DOWNWIND (K), LOGIC (K)
END DO
AVE_POS = AVE_POS / C_DW
AVE_NEG = ABS (AVE_NEG) / C_UW

!        GET AVERAGE STREAMLINE DISTANCES
!b !  C_UW = COUNT ((.NOT. IS_DOWNWIND))
!b !  C_DW = LT_N - C_UW ! by pass compiler bug COUNT (IS_DOWNWIND)

   IF ( C_DW < 1 .OR. C_DW >= LT_N .OR. &
        C_UW < 1 .OR. C_UW >= LT_N ) THEN
     PRINT *,'IMPOSSIBLE DOWNWIND COUNT, LT, UP', C_DW, LT_N, C_UW
     PRINT *,'THIS_EL     ', THIS_EL    
     !b PRINT *,'IS_DOWNWIND ', IS_DOWNWIND
     PRINT *,'LOGIC       ', LOGIC
     PRINT *,'DOWN dist', DOWN
     STOP 'GET_AVE_DOWNWIND_DIST: ERROR, IMPOSSIBLE DOWNWIND COUNT'
   END IF

!b !  AVE_POS = SUM (DOWN, IS_DOWNWIND) / C_DW
!b !  AVE_NEG = ABS (SUM (DOWN, .NOT. IS_DOWNWIND)) / C_UW

if ( AVE_POS <= 0.d0 ) then
  print *,'ERROR, no DW dist for element ', THIS_EL, AVE_POS
  print *,'DOWN dist ', DOWN
  !b PRINT *,'IS_DOWNWIND ', IS_DOWNWIND
  stop    'GET_AVE_DOWNWIND_DIST: ERROR, no DW dist for element '
end if

!b call at(1488)
!        GET TOTAL DOWNWIND AVE DISTANCES
        WHERE ( DOWN >= 0.d0 ) 
          DOWN_AVE = DOWN + AVE_NEG
        ELSEWHERE             
          DOWN_AVE = ABS (DOWN) + AVE_POS 
        END WHERE
!b call at(1495)

!if ( THIS_EL == 1 ) then
!  print *,'GET_AVE_DOWNWIND_DIST ', THIS_EL
!  print *,'POS_AVE, NEG_AVE ', AVE_POS, AVE_NEG
!  PRINT *,'DOWN dist', DOWN
!  print *,'DOWN_AVE ', DOWN_AVE
!  !b PRINT *,'IS_DOWNWIND ', IS_DOWNWIND
!  PRINT *,'LOGIC       ', LOGIC
!end if
END SUBROUTINE GET_AVE_DOWNWIND_DIST 

SUBROUTINE GET_DOWNWIND_LOGIC (RADIAL, UNIT_V, DOWN, IS_DOWNWIND)
 Use System_Constants ! for N_SPACE
 Use Elem_Type_Data   ! for LT_N
 IMPLICIT NONE
  REAL(DP), INTENT (IN ) :: RADIAL (LT_N, N_SPACE) ! relative positions
  REAL(DP), INTENT (IN ) :: UNIT_V (N_SPACE)       ! vel unit vector
  REAL(DP), INTENT (OUT) :: DOWN   (LT_N)          ! downwind if > 0
  LOGICAL,  INTENT (OUT) :: IS_DOWNWIND (LT_N)     ! true if downwind 
  INTEGER :: K
  !b LOGICAL :: LOGIC (LT_N)   ! debug

   IF ( ALL (UNIT_V == 0.d0) ) STOP 'Zero unit vector. g-dw-l' 
   IF ( ALL (RADIAL == 0.d0) ) STOP 'Zero radial array. g-dw-l' 
  
!          IDENTIFY DOWNWIND DISTANCES, AND LOGIC
    DOWN  = MATMUL ( RADIAL, UNIT_V )    ! POSITIVE & NEGATIVE
    IF ( ALL (DOWN == 0.d0) ) STOP 'Zero DW vector. g-dw-l'

    WHERE ( DOWN >= 0.d0 ) ; IS_DOWNWIND = .TRUE.
    ELSEWHERE              ; IS_DOWNWIND = .FALSE.  ; END WHERE

   !b DO K = 1, LT_N
     !b IF ( DOWN (K) >= 0.d0 ) THEN
       !b LOGIC (K) = .TRUE.
     !b ELSE
       !b LOGIC (K) = .FALSE.
     !b END IF
   !b END DO

if ( all(IS_DOWNWIND) ) then
  print *,'ERROR all nodes downwind for element ', THIS_EL
  print *,'RADIAL '; call rprint (RADIAL, LT_N, N_SPACE, 1)
  print *,'UNIT_V ', UNIT_V
  STOP 'GET_DOWNWIND_LOGIC: ERROR all nodes downwind'
end if

!b if ( any (IS_DOWNWIND .NEQV. LOGIC) ) then
  !b print *,'ERROR inconsistent logic for element ', THIS_EL
  !b print *,'IS_DW ', IS_DOWNWIND
  !b print *,'LOGIC ', LOGIC
  !b STOP 'GET_DOWNWIND_LOGIC: ERROR in logic'
!b end if

!if ( THIS_EL == 1 ) then
!  print *,'============== GET_DOWNWIND_LOGIC ================'
!  print *,'RADIAL '; call rprint (RADIAL, LT_N, N_SPACE, 1)
!  print *,'UNIT_V ', UNIT_V
!  print *,'DOWN ', DOWN
!  print *,'IS_DW ', IS_DOWNWIND
!  !b print *,'LOGIC ', LOGIC
!end if
END SUBROUTINE GET_DOWNWIND_LOGIC 

FUNCTION VOLUME_UPWIND_LENGTH (H_INTG, DOWN_NODAL)  RESULT (LENGTH)
! Get DOWN from UPWIND_DIST_AND_LOGIC
 Use System_Constants ! for N_SPACE
 Use Elem_Type_Data   ! for LT_N
 IMPLICIT NONE
  REAL(DP), INTENT (IN ) :: H_INTG     (LT_N)  ! H integral    
  REAL(DP), INTENT (IN ) :: DOWN_NODAL (LT_N)  ! nodal downwind lengths
  REAL(DP)               :: LENGTH         ! volume average downwind
  REAL(DP) :: VOLUME                       ! element volume          
  
   IF ( ANY (DOWN_NODAL < 0.d0) ) THEN
     PRINT *,'WARNING: VOLUME_UPWIND_LENGTH, NEGATIVE DATA.'
     PRINT *,'DOWN_NODAL ', DOWN_NODAL
     N_WARN = N_WARN + 1
   END IF ! negative lengths
   VOLUME = SUM (H_INTG)
   IF ( VOLUME <= 0.d0 ) STOP 'VOLUME_UPWIND_LENGTH: VOLUME <= 0'

   LENGTH = DOT_PRODUCT (H_INTG, DOWN_NODAL) / VOLUME ! Integral average
!b if ( THIS_EL == 1 ) print *,'VOLUME_UPWIND_LENGTH ', LENGTH !b
END FUNCTION VOLUME_UPWIND_LENGTH 

SUBROUTINE STOCASTIC_VOLATILITY_VECTOR (XYZ, U)
USE SYSTEM_CONSTANTS   ! FOR N_SPACE, THIS_EL
Use Sys_Properties_Data
 IMPLICIT NONE
 REAL(DP), INTENT (IN)  :: XYZ    (N_SPACE)        ! point
 REAL(DP), INTENT (OUT) :: U      (N_SPACE)        ! velocity

 !  Variable E, u items
  REAL(DP), SAVE :: Sigma, Kappa, Rho, rate     ! Global data
  REAL(DP), SAVE :: Theta, Lamda, source        ! Global data

  IF ( THIS_EL == 1 ) THEN ! Get global data
! Misc reals are: 1-Sigma 2-Kappa 3-Rho 4-r 5-Theta 6-Lamda 7-Q
     Sigma = 0.9d0 ; Kappa = 5 ; Rho = 0.1d0    ! Clarke-Parrott
     Rate = 0.1d0 ; Theta = 0.16d0              ! Clarke-Parrott
     Lamda = 0 ; Source = 0                     ! Clarke-Parrott
     IF ( REALS > 0 ) Sigma  = GET_REAL_MISC (1)
     IF ( REALS > 1 ) Kappa  = GET_REAL_MISC (2)
     IF ( REALS > 2 ) Rho    = GET_REAL_MISC (3)
     IF ( REALS > 3 ) Rate   = GET_REAL_MISC (4)
     IF ( REALS > 4 ) Theta  = GET_REAL_MISC (5)
     IF ( REALS > 5 ) Lamda  = GET_REAL_MISC (6)
     IF ( REALS > 6 ) Source = GET_REAL_MISC (7)
  END IF ! first el

  U (1) = -(Rate - XYZ (2) - Rho * Sigma / 2) * XYZ (1)
  U (2) = -Kappa * (Theta - XYZ (2)) + Sigma * Sigma / 2  &
        + (Lamda + Rho * Sigma / 2) * XYZ (2)

END SUBROUTINE STOCASTIC_VOLATILITY_VECTOR 

SUBROUTINE VELOCITY_AT_POINT (XYZ, U, UNIT_V, SPEED)
USE SYSTEM_CONSTANTS   ! FOR N_SPACE
Use Sys_Properties_Data
 IMPLICIT NONE
 REAL(DP), INTENT (IN)  :: XYZ    (N_SPACE)        ! point
 REAL(DP), INTENT (OUT) :: U      (N_SPACE)        ! velocity
 REAL(DP), INTENT (OUT) :: UNIT_V (N_SPACE), SPEED ! unit vec, speed

 INTEGER  :: FLOW_CASE !b , SCALE
 REAL(DP) :: RATE, ANGLE, RATIO, Y_ZERO, SCALE
 REAL(DP), PARAMETER :: NIL = EPSILON (1.d0)
 REAL(DP), PARAMETER :: y_1 = 0.25d0, y_2 = 0.33333333333333d0, &
                        y_12 = 0.0833333333333d03

  FLOW_CASE = 0 ; IF ( INTEGERS > 0 ) FLOW_CASE = GET_INTEGER_MISC (1)
  SCALE     = 1 ; IF ( INTEGERS > 1 ) SCALE     = GET_INTEGER_MISC (2)
  if ( SCALE ==  0.d0 ) print *,'NOTE SCALE = 0 gives zero velocity'
  !b IF ( SCALE < 0.d0 ) SCALE = SQRT ( ABS (SCALE) )

  SELECT CASE (FLOW_CASE)
    CASE (1) ! Smith-Hutton
      U(1) =  2.d0 * xyz(2) * (1.d0-xyz(1)**2) ! Smith-Hutton
      U(2) = -2.d0 * xyz(1) * (1.d0-xyz(2)**2) ! Smith-Hutton
    CASE (2) ! Donald's money
      RATE = 1.d0 ; IF ( REALS > 2 ) RATE = GET_REAL_MISC (4)
      U    = -RATE * XYZ
    CASE (3) ! ! Hughes cosine hill
      U(1) = - xyz(2) ; U(2) = xyz(1)          ! Brooks-Hughes
    CASE (4) ! Hughes advection skew
      ANGLE = 0.d0 ; IF ( INTEGERS > 2 ) ANGLE =  GET_INTEGER_MISC (3)
      ANGLE = ANGLE / 57.3d0     ! radians
      y_zero = Y_1 + XYZ(1) * TAN (ANGLE)
      IF ( XYZ (2) <= y_zero ) THEN
        U(1) = cos (angle) ; U(2) = sin (angle)
      ELSE IF ( XYZ (2) >= (y_zero + y_12) ) THEN
        U = 0.d0
      ELSE
        RATIO = (y_12 - (XYZ (2) - y_zero))/y_12
        U(1) = cos (angle) * ratio
        U(2) = sin (angle) * ratio
      END IF
    CASE (5) ! Constant advection skew
      ANGLE = 0.d0 ; IF ( INTEGERS > 2 ) ANGLE =  GET_INTEGER_MISC (3)
      ANGLE = ANGLE / 57.3d0     ! radians
      U(1) = cos (angle) ; U(2) = sin (angle)
    CASE DEFAULT ! Unit flow in x
      U = 0.d0 ; U(1) = 1.d0
  END SELECT
  U = U * SCALE

!       GET SPEED, AND UNIT VECTOR
   SPEED = SQRT ( SUM (U **2) )        ! SPEED
   !b IF ( N_SPACE == 1 ) SPEED = SIGN (SPEED, U(1))
   IF ( SPEED > NIL ) THEN             ! no divide by 0
     UNIT_V = U / SPEED                ! UPWIND UNIT VECTOR
   ELSE ! Null vectors, place at equal angles to axes
     SPEED = NIL   
     UNIT_V = 1.d0 / SQRT (DBLE(N_SPACE)) ! for N_SPACE components
     U = NIL * UNIT_V                     ! for N_SPACE components
   END IF
END SUBROUTINE VELOCITY_AT_POINT

SUBROUTINE UNIT_VECTOR_FORM (U, UNIT_V, SPEED)
USE SYSTEM_CONSTANTS   ! FOR N_SPACE
Use Sys_Properties_Data
 IMPLICIT NONE
 REAL(DP), INTENT (OUT) :: U      (N_SPACE)        ! velocity
 REAL(DP), INTENT (OUT) :: UNIT_V (N_SPACE), SPEED ! unit vec, speed

 REAL(DP), PARAMETER :: NIL = EPSILON (1.d0)
 REAL(DP), PARAMETER :: ROOT_2 = 0.70710678d0, ROOT_3 = 0.57735021d0

!       GET SPEED, AND UNIT VECTOR
   SPEED = SQRT ( SUM (U **2) )        ! SPEED
   IF ( SPEED > NIL ) THEN             ! no divide by 0
     UNIT_V = U / SPEED                ! UPWIND UNIT VECTOR
   ELSE ! Null vectors, place at equal angles to axes
     SELECT CASE (N_SPACE)
       CASE (1) ; UNIT_V = 1.d0
       CASE (2) ; UNIT_V = ROOT_2
       CASE (3) ; UNIT_V = ROOT_3
       CASE DEFAULT
         UNIT_V = 1.d0 / SQRT (DBLE(N_SPACE)) ! for N_SPACE components
     END SELECT
     SPEED = NIL   
     U = NIL * UNIT_V                     ! for N_SPACE components
   END IF
END SUBROUTINE UNIT_VECTOR_FORM

SUBROUTINE DIFFUSION_UPWIND (E, UNIT_V, E_UPWIND, E_CROSS)
 Use System_Constants ! for N_SPACE, N_R_B
 IMPLICIT NONE
  REAL(DP), INTENT (IN)  :: E (N_R_B, N_R_B) ! constitutive law
  REAL(DP), INTENT (IN)  :: UNIT_V (N_SPACE) ! unit vector
  REAL(DP), INTENT (OUT) :: E_UPWIND, E_CROSS
  REAL(DP) :: E_LOCAL (N_R_B, N_R_B), LAMDA (N_SPACE, N_SPACE)

   IF ( N_SPACE == 1 ) THEN ! SPECIAL CASE
    E_UPWIND = E (1, 1) ; E_CROSS = 0.d0
    RETURN
   END IF ! 1-D

   IF ( N_SPACE /= 2 .OR. N_R_B /= 2 )                         &
     STOP 'DIFFUSION_UPWIND invalid if N_SPACE /= 2, N_R_B /= 2'

!   SET DIRECTION COSINES
   LAMDA (1, 1:2) = UNIT_V (1:2)
   LAMDA (2, 1) = -UNIT_V (2) ; LAMDA (2, 2) = UNIT_V (1)

!   TRANSFORM
   E_LOCAL = MATMUL (LAMDA, MATMUL (E, TRANSPOSE (LAMDA)))
   E_UPWIND = E_LOCAL (1, 1)
   E_CROSS  = E_LOCAL (2, 2)
END SUBROUTINE DIFFUSION_UPWIND 

SUBROUTINE DIFFUSION_LOCAL (E, UNIT_V, E_LOCAL)
! Get diffusion value in the direction of a unit vector
! Vector is usually // to streamline, but may be otherwise
 Use System_Constants ! for N_SPACE, N_R_B
 IMPLICIT NONE
  REAL(DP), INTENT (IN)  :: E (N_R_B, N_R_B) ! constitutive law
  REAL(DP), INTENT (IN)  :: UNIT_V (N_SPACE) ! unit vector
  REAL(DP), INTENT (OUT) :: E_LOCAL          ! principle value

   IF ( N_SPACE /= N_R_B )                            &
     STOP 'DIFFUSION_LOCAL invalid if N_SPACE /= N_R_B'

   IF ( N_SPACE == 1 ) THEN ! SPECIAL CASE
    E_LOCAL = E (1, 1) ; RETURN
   END IF ! 1-D

!   TRANSFORM (FOR ONE PRINCIPLE VALUE)
   E_LOCAL = DOT_PRODUCT (UNIT_V, MATMUL (E, UNIT_V))
END SUBROUTINE DIFFUSION_LOCAL 

SUBROUTINE CENTER_DGH_D2GH (DGH, D2GH, DET_WT)
! GET PHYSICAL DERIVATIVES AT THE ELEMENT CENTER
 Use System_Constants ! for N_SPACE, N_2_DER, CONSTANT_J
 Use Elem_Type_Data   ! for LT_N, LT_PARM, DLH (LT_PARM, LT_N),
                      ! H (LT_N), COORD (LT_N, N_SPACE),
                      ! D2LH (N_2_DER, LT_N)
 IMPLICIT NONE
  REAL(DP), INTENT (OUT) :: DGH  (N_SPACE, LT_N)    
  REAL(DP), INTENT (OUT) :: D2GH (N_2_DER, LT_N), DET_WT
!         Automatic Arrays
  REAL(DP) :: AJ (N_SPACE, N_SPACE), AJ_INV (N_SPACE, N_SPACE), DET
  REAL(DP) :: ONE_PT (LT_PARM), ONE_WT         ! 1 pt rule     
!!        Second derivatives
!  REAL(DP) :: P_AJ      (N_2_DER, N_2_DER)    ! AJ PRODUCTS
!  REAL(DP) :: P_AJ_INV  (N_2_DER, N_2_DER)    ! PRODUCTS INVERSE
!  INTEGER  :: I_ERROR
!  LOGICAL  :: CONSTANT_J = .FALSE. !b move to System_Constants XXX

!        LOCAL AND GLOBAL FIRST DERIVATIVES AT CENTROID
    CALL GET_ONE_PT_RULE (ONE_PT, ONE_WT)       ! local point
    CALL SCALAR_DERIVS (ONE_PT, DLH)     ! SCALAR DERIVATIVES
    AJ  = MATMUL (DLH, COORD)                      ! JACOBIAN
    CALL INVERT_JACOBIAN (AJ, AJ_INV, DET, N_SPACE) ! INVERSE
    DGH = MATMUL (AJ_INV, DLH)         ! PHYSICAL DERIVATIVES
    DET_WT = DET * ONE_WT
!b if ( THIS_EL == 1 ) print *,'DET, ONE_WT ', DET, ONE_WT

D2GH = 0.d0  ! default here, enhance later, uncomment below
!!       SECOND DERIVATIVES
!    CALL SCALAR_2ND_DERIVS (ONE_PT, D2LH) 
!
!    CALL JACOBIAN_PRODUCTS (AJ, P_AJ) ! ALWAYS NEEDED
!    ! NOTE: D2GH     = P_AJ_INV (D2HL - HESSIAN * DGH)
!    !       HESSIAN  = D2HL * COORD
!    CALL INV_SMALL_MAT (N_2_DER, P_AJ, P_AJ_INV, I_ERROR)
!    IF ( I_ERROR > 0 ) THEN
!      PRINT *, 'ERROR: CENTER_DGH_D2GH, NO PRODUCT INVERSE'
!      STOP     'ERROR: CENTER_DGH_D2GH, NO PRODUCT INVERSE'
!    END IF ! INVERSION ERROR
!    IF ( CONSTANT_J ) THEN ! HESSIAN = 0
!      D2GH = MATMUL (P_AJ_INV, D2LH)
!    ELSE
!      D2GH = MATMUL (P_AJ_INV, (D2LH            &
!           - MATMUL (MATMUL (D2LH, COORD), DGH)))
!    END IF ! ZERO HESSIAN
END SUBROUTINE CENTER_DGH_D2GH

SUBROUTINE SECOND_GLOBAL_DERIV_H (ONE_PT, DGH, D2GH)
! ---------------------------------------------------------
!    GET PHYSICAL DERIVATIVES AT A POINT, RE_BUILD JACOBIAN
! ---------------------------------------------------------
 Use System_Constants ! for N_SPACE, N_2_DER, CONSTANT_J
 Use Elem_Type_Data   ! for LT_N, LT_PARM, DLH (LT_PARM, LT_N),
                      ! H (LT_N), COORD (LT_N, N_SPACE),
                      ! D2LH (N_2_DER, LT_N)
  REAL(DP), INTENT (IN)  :: ONE_PT (LT_PARM)
  REAL(DP), INTENT (OUT) :: DGH  (N_SPACE, LT_N)    
  REAL(DP), INTENT (OUT) :: D2GH (N_2_DER, LT_N)
!         Automatic Arrays
  REAL(DP) :: AJ (N_SPACE, N_SPACE), AJ_INV (N_SPACE, N_SPACE), DET
!        Second derivatives
  REAL(DP) :: P_AJ      (N_2_DER, N_2_DER)    ! AJ PRODUCTS
  REAL(DP) :: P_AJ_INV  (N_2_DER, N_2_DER)    ! PRODUCTS INVERSE
  INTEGER, SAVE :: I_ERROR = 0

!        GLOBAL FIRST DERIVATIVES 
    AJ  = MATMUL (DLH, COORD)                       ! JACOBIAN
    CALL INVERT_JACOBIAN (AJ, AJ_INV, DET, N_SPACE) ! INVERSE
    DGH = MATMUL (AJ_INV, DLH)         ! PHYSICAL DERIVATIVES

!       SECOND DERIVATIVES
    CALL SCALAR_2ND_DERIVS (ONE_PT, D2LH) 

    CALL JACOBIAN_PRODUCTS (AJ, P_AJ) ! ALWAYS NEEDED
    ! NOTE: D2GH     = P_AJ_INV (D2HL - HESSIAN * DGH)
    !       HESSIAN  = D2HL * COORD
    CALL INV_SMALL_MAT (N_2_DER, P_AJ, P_AJ_INV, I_ERROR)
    IF ( I_ERROR > 0 ) THEN
      PRINT *, 'ERROR: SECOND_GLOBAL_DERIV_H, NO PRODUCT INVERSE'
      STOP     'ERROR: SECOND_GLOBAL_DERIV_H, NO PRODUCT INVERSE'
    END IF ! INVERSION ERROR
    IF ( CONSTANT_J ) THEN ! HESSIAN = 0
      D2GH = MATMUL (P_AJ_INV, D2LH)
    ELSE
      D2GH = MATMUL (P_AJ_INV, (D2LH            &
           - MATMUL (MATMUL (D2LH, COORD), DGH)))
    END IF ! ZERO HESSIAN
END SUBROUTINE SECOND_GLOBAL_DERIV_H

!SUBROUTINE GET_PT_UPWIND_DATA (XYZ, RADIAL, DGH, E_UPWIND, &
!                               PECLET, UP_WIND, GEOM_H  )
! Use System_Constants ! for N_SPACE, N_2_DER, CONSTANT_J
! Use Elem_Type_Data   ! for LT_N, LT_PARM, DLH (LT_PARM, LT_N),
!                      ! H (LT_N), COORD (LT_N, N_SPACE),
!                      ! D2LH (N_2_DER, LT_N)
! IMPLICIT NONE
!  REAL(DP), INTENT (IN) :: XYZ    (N_SPACE)       ! pt in space
!  REAL(DP), INTENT (IN) :: RADIAL (LT_N, N_SPACE) ! relative positions
!  REAL(DP), INTENT (IN) :: DGH    (N_SPACE, LT_N) ! global deriv of H
!  REAL(DP), INTENT (IN) :: E_UPWIND               ! diffusion value
!  REAL(DP), INTENT (OUT) :: PECLET                ! local Peclet number
!  REAL(DP), INTENT (OUT) :: GEOM_H                ! geometric length
!  REAL(DP), INTENT (OUT) :: UP_WIND (LT_N)        ! upwind additions
!
!!          Automatic arrays
!  REAL(DP) :: U           (N_SPACE)          ! VELOCITY
!  REAL(DP) :: UNIT_V      (N_SPACE), SPEED   ! Unit vec, speed
!  REAL(DP) :: U_DGH       (LT_N)             ! streamline gradient
!  REAL(DP) :: DOWN        (LT_N)             ! downwind if > 0
!  LOGICAL  :: IS_DOWNWIND (LT_N)             ! true if downwind node
!  !b REAL(DP) :: D2GH   (N_2_DER, LT_N)     ! SUPG
!
!   CALL VELOCITY_AT_POINT (XYZ, U, UNIT_V, SPEED) 
!!       GEOMETRY BASED TAU
!   CALL UPWIND_DIST_AND_LOGIC (RADIAL, UNIT_V, DOWN, &
!                               GEOM_H  , IS_DOWNWIND)
!   PECLET = ABS (SPEED) * GEOM_H   / E_UPWIND
!if ( THIS_EL == 1 ) print *,'PECLET, GEOM_H   ', PECLET, GEOM_H   
!
!!        DOT VELOCITY WITH SOLUTION GRADIENT
!   U_DGH  = MATMUL (U, DGH) ! v_x dp/dx + v_y dp/dy + v_z dp/dz
!
!   UP_WIND (:) = DOWN (:) * U_DGH (:) / SPEED ! later * ALPHA/2
!END SUBROUTINE GET_PT_UPWIND_DATA

SUBROUTINE PECLET_OPTIMAL_RULE (P, R)
Use System_Constants ! for DP
IMPLICIT NONE
  REAL(DP), INTENT (IN ) :: P         ! local Peclet number
  REAL(DP), INTENT (OUT) :: R         ! optimal rule 0 <= R < 1

! Optimal Rule: R = coth(P) - 1/P, underestimated < 1 %
  IF ( P > 3.d0 ) THEN ! Critical approximation
    R = 1 - 1 / P
  ELSE ! continued fraction fit
    R = P * (315 + 14 * P **2) / (945 + 105 * P **2 + P **4)
  END IF
END SUBROUTINE PECLET_OPTIMAL_RULE 

SUBROUTINE INTEGRAL_H_ON_L3 (LENGTH, C_E)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! INTEGRATE SHAPE FUNCTIONS OF A 3 NODE LINE ELEMENT
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module  ! for DP
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: LENGTH   ! PHYSICAL LENGTH
 REAL(DP), INTENT(OUT) :: C_E (3)  ! INTEGRAL OF H
! LOCAL NODE COORD. ARE -1,0,+1    1-----2-----3

  C_E = LENGTH * (/ 1.d0, 4.d0, 1.d0 /) / 6.d0
END SUBROUTINE INTEGRAL_H_ON_L3

SUBROUTINE INTEGRAL_HT_H_ON_L3 (LENGTH, MASS_E)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! INTEGRATE H-TRANSPOSE H ON A 3 NODE LINE ELEMENT
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module  ! for DP
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: LENGTH         ! PHYSICAL LENGTH
 REAL(DP), INTENT(OUT) :: MASS_E (3, 3)  ! INTEGRAL OF H' H
 REAL(DP), PARAMETER   :: MASS (3, 3) = RESHAPE (       &
   (/ 4, 2, -1, 2, 16, 2, -1, 2, 4 /), (/3, 3/) ) / 30.d0 
! LOCAL NODE COORD. ARE -1,0,+1    1-----2-----3

  MASS_E = LENGTH * MASS
END SUBROUTINE INTEGRAL_HT_H_ON_L3

SUBROUTINE INTEGRAL_DHDXT_DHDX_ON_L3 (LENGTH, STIFF_E)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! INTEGRATE DHDX-TRANSPOSE DHDX ON A 3 NODE LINE ELEMENT
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module  ! for DP
 IMPLICIT NONE
 REAL(DP), INTENT(IN)  :: LENGTH         ! PHYSICAL LENGTH  
 REAL(DP), INTENT(OUT) :: STIFF_E (3, 3) ! INTEGRAL DHDX' DHDX
 REAL(DP), PARAMETER   :: STIFF (3, 3) = RESHAPE (       &
   (/ 7, -8, 1, -8, 16, -8, 1, -8, 7/), (/3, 3/) ) / 3.d0 
! LOCAL NODE COORD. ARE -1,0,+1    1-----2-----3

  STIFF_E = STIFF / LENGTH
END SUBROUTINE INTEGRAL_DHDXT_DHDX_ON_L3

SUBROUTINE INTEGRAL_HT_DHDX_ON_L3 (U_E)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! INTEGRATE H-TRANSPOSE DHDX ON A 3 NODE LINE ELEMENT
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use Precision_Module  ! for DP
 IMPLICIT NONE
 REAL(DP), INTENT(OUT) :: U_E (3, 3) ! INTEGRAL H' DHDX
 REAL(DP), PARAMETER   :: U (3, 3) = RESHAPE (       &
   (/ -3, -4, 1, 4, 0, -4, -1, 4, 3/), (/3, 3/) ) / 6.d0 
! LOCAL NODE COORD. ARE -1,0,+1    1-----2-----3

  U_E = U 
END SUBROUTINE INTEGRAL_HT_DHDX_ON_L3

SUBROUTINE SPACE_TIME_EBC_UPDATE (NDEX_C_EQ, DD, COEF_C_EQ)
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
! PUSH CONTINUOUS SPACE-TIME RESULTS, DD, INTO PREVIOUS
!  ESSENTIAL BOUNDARY CONDITION VALUES, COEF_C_EQ
! *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-* *-*
Use System_Constants  ! for DP, MAX_ACT, N_CEQ, N_D_FRE
 IMPLICIT NONE
 INTEGER,  INTENT(IN)    :: NDEX_C_EQ (MAX_ACT, N_CEQ)
 REAL(DP), INTENT(IN)    :: DD        (N_D_FRE)
 REAL(DP), INTENT(INOUT) :: COEF_C_EQ (MAX_ACT, N_CEQ)
 INTEGER :: ID, IEQ, IT, KOUNT

! N_CEQ           = NUMBER OF CONSTRAINT EQUATIONS
! COEF_C_EQ (I,J) = COEFF TERM I+1 OF CONSTRAINT EQ J
! NDEX_C_EQ (I,J) = DOF NUMBER OF TERM I OF CONSTRAINT J
! DD              = SYSTEM LIST OF NODAL PARAMETERS
! In subroutine INPUT_CONSTRAINT_EQS: (Type 1)
!    NDEX_C_EQ (1, IEQ) = N_G_DOF * (NODE_1 - 1) + IPAR_1
!    COEF_C_EQ (1, IEQ) = A_1

! Move the last half of the system dof into EBC for
! the first half for the next space-time step
!b call at(1920); PRINT *,'N_CEQ = ', N_CEQ

! Loop over dof
  KOUNT = 0
  DOF: DO ID = (N_D_FRE/2 + 1), N_D_FRE
    ! Search for corresponding EBC
    IT = ID - N_D_FRE/2
!b print *,'ID, IT = ', ID, IT
    DO IEQ = 1, N_CEQ
      !b print *, 'IEQ, NDEX ', IEQ, NDEX_C_EQ (1, IEQ) 
      IF ( NDEX_C_EQ (1, IEQ) == IT ) THEN
        !b print *, 'COEFF, DD ', COEF_C_EQ (1, IEQ), DD (ID)
        COEF_C_EQ (1, IEQ) = DD (ID)
        KOUNT = KOUNT + 1
        CYCLE DOF
      END IF
    END DO ! over constraints
  END DO DOF ! over half system dof
  IF ( KOUNT /= N_D_FRE/2 ) THEN ! error
    PRINT *, 'WARNING, SPACE_TIME_EBC_UPDATE COUNT IS ', KOUNT
    N_WARN = N_WARN + 1
  END IF

END SUBROUTINE SPACE_TIME_EBC_UPDATE 

SUBROUTINE ONE_STEP_TRANSIENT_AT_EL (K, M, F)
! ..............................................................
!          CONVERT S, EL_M, C INTO K, M, F FOR
!   ASSEMBLE OF ELEMENT MATRICES INTO ALGEBRAIC SYSTEM FOR
!          A ONE STEP TIME INTEGRATION METHOD
! ..............................................................
! So element matrix     S * D (t) + EL_M * D_dot (t) = C (t)
! Converts to algebraic K * D (t) = F (t) with
! K a function of S, EL_M, STEP_SIZE
! M a function of S, EL_M, STEP_SIZE, so that
! F is a function of M * D and C
! Accesses globals: TIME_METHOD, TIME_STEP
Use System_Constants
Use Elem_Type_Data     ! for PT (LT_PARM, LT_QP), WT (LT_QP),
                       ! G (LT_GEOM, LT_QP), H (LT_FREE),
                       ! DLG (LT_PARM, LT_GEOM, LT_QP),
                       ! DLH (LT_PARM, LT_FREE, LT_QP),
                       ! S (LT_FREE, LT_FREE), C (LT_FREE),
                       ! EL_M (LT_FREE, LT_FREE), DIAG_M (LT_FREE),
                       ! D (LT_FREE), ELEM_NODES (LT_N)
IMPLICIT NONE
REAL(DP), INTENT(OUT) :: K (LT_FREE, LT_FREE), F (LT_FREE)
REAL(DP), INTENT(OUT) :: M (LT_FREE, LT_FREE)
REAL(DP) :: Del_t
INTEGER  :: J                   

!          SQUARE MATRICES (STIFFNESS & MASS)
!   S = S + MATMUL (TRANSPOSE(DGH), DGH)* WT (IQ) * DX_DR * K_e
!   EL_M = EL_M + OUTER_PRODUCT (H, H)  * WT (IQ) * DX_DR * Rho_e

   IF ( AVERAGE_MASS ) THEN  ! average consistent & diagonal forms
     CALL DIAGONALIZE_SQ_MATRIX (LT_FREE, EL_M, DIAG_M) 
     EL_M = 0.5d0 * EL_M
     DO J = 1, LT_FREE
       EL_M (J, J) = EL_M (J, J) + DIAG_M (J) * 0.5d0
     END DO
   END IF
   IF ( DIAGONAL_MASS ) THEN  ! use the scaled diagonal mass form
     IF ( DEBUG_EL_SQ ) PRINT *,'BEFORE DIAG_M AT ELEMENT ', THIS_EL
     CALL DIAGONALIZE_SQ_MATRIX (LT_FREE, EL_M, DIAG_M) 
     EL_M = 0.d0
     DO J = 1, LT_FREE
       EL_M (J, J) = DIAG_M (J) 
     END DO
   END IF

    IF ( debug_el_sq .AND. THIS_EL == 1 ) THEN
      PRINT *, 'EL_M, EL_S, EL_C for one step method = ', TIME_METHOD
      CALL RPRINT (EL_M, LT_FREE, LT_FREE, 1)
      CALL RPRINT (S, LT_FREE, LT_FREE, 1)
      CALL RPRINT (C, LT_FREE, 1      , 1)
    END IF

   !b print *,'THIS_EL, D ', THIS_EL, D
   Del_t = TIME_STEP
   SELECT CASE ( TIME_METHOD ) ! for integration rule
     CASE ( 2 ) ! Crank-Nicolson
        F = C 
        M = EL_M / Del_t - S / 2 ! multiplies D
        K = EL_M / Del_t + S / 2
     CASE ( 3 ) ! Galerkin continuous in time
        F = C
        M = EL_M / Del_t - S / 3 ! multiplies D
        K = 2 * S / 3.d0 + EL_M / Del_t
     CASE ( 4 ) ! Least Squares in time, assuming C constant
          ! check if F signs are both - not +
        F = -MATMUL (TRANSPOSE(S),C) * Del_t / 2   &
          -  MATMUL (TRANSPOSE(EL_M),C) 
        M = MATMUL (TRANSPOSE(EL_M), EL_M) / Del_t &
          + MATMUL (TRANSPOSE(S),EL_M) / 2         &
          - MATMUL (TRANSPOSE(EL_M),S) / 2         &
          - MATMUL (TRANSPOSE(S),S) * Del_t / 6
        K = MATMUL (TRANSPOSE(S),S) * Del_t / 3    &
          + MATMUL (TRANSPOSE(S),EL_M) / 2         &
          + MATMUL (TRANSPOSE(EL_M),S) / 2         &
          + MATMUL (TRANSPOSE(EL_M),EL_M) / Del_t 
     CASE DEFAULT ! Method 1, Euler forward difference
        F = C 
        M = EL_M / Del_t 
        K = S + EL_M / Del_t
   END SELECT

    IF ( debug_el_sq .AND. THIS_EL == 1 ) THEN
      PRINT *, 'EL_S, EL_M, K_t, M_t, F_t, D_t for method', TIME_METHOD
      CALL RPRINT (S, LT_FREE, LT_FREE, 1)
      CALL RPRINT (EL_M, LT_FREE, LT_FREE, 1)
      CALL RPRINT (K, LT_FREE, LT_FREE, 1)
      CALL RPRINT (M, LT_FREE, LT_FREE, 1)
      CALL RPRINT (F, LT_FREE, 1      , 1)
      CALL RPRINT (D, LT_FREE, 1      , 1)
    END IF

!   APPLY GENERALIZED TRAPEZOIDAL RULE TO REFORM TWO SQUARE MATRICES
!                     M_SKY * D_DOT + S_SKY * D = C
!                 INTO THE TWO USED IN TIME INTEGRATION
! (M_SKY/T_STEP + TS_B*S_SKY)*D_1 = (M_SKY/T_STEP + TS_A*S_SKY)*D_0
!                                 + TS_B*C_1 + TS_A*C_0
!       BUT, TS_A = 1-TS_B, AND C_0 = C_1 =C, ASSUMED CONSTANT
!         S_SKY_new*D_1 = M_SKY_new*D_0 + TS_B*C_1 + TS_A*C_0
!                       = M_SKY_new*D_0 + C 
!!          NOW FORM CURRENT RHS VECTOR               
!    DO J = 1, N_D_FRE
!      TEMP       = CC (J)
!      CC (J)     = CC (J) * TS_B + CC_OLD (J) * TS_A + C_MD (J)
!      CC_OLD (J) = TEMP
!    END DO ! FOR C_1
END SUBROUTINE ONE_STEP_TRANSIENT_AT_EL