      MODULE shape_hex27n
!####################################################################
      
      IMPLICIT NONE
      
!--------------------------------------------------------------------
      
      INTEGER, PARAMETER, PRIVATE :: kreal = KIND( 0.0D0 )
      
!--------------------------------------------------------------------
      
      CONTAINS
      
!####################################################################
      SUBROUTINE ShapeFunc_hex27n(lcoord, func)
!####################################################################
      
      REAL(KIND = kreal), INTENT(IN)  :: lcoord(3)
      REAL(KIND = kreal), INTENT(OUT) :: func(27)
      
!--------------------------------------------------------------------
      
      REAL(KIND = kreal) :: xi(27), eta(27), zeta(27)
      REAL(KIND = kreal) :: n_xi(27), n_eta(27), n_zeta(27)
      
      INTEGER :: na
      
!--------------------------------------------------------------------
      
      ! xi-coordinate at a node in a local element
      xi(1)    = -1.0D0
      xi(2)    =  1.0D0
      xi(3)    =  1.0D0
      xi(4)    = -1.0D0
      xi(5)    = -1.0D0
      xi(6)    =  1.0D0
      xi(7)    =  1.0D0
      xi(8)    = -1.0D0
      xi(9)    =  0.0D0
      xi(10)   =  1.0D0
      xi(11)   =  0.0D0
      xi(12)   = -1.0D0
      xi(13)   =  0.0D0
      xi(14)   =  1.0D0
      xi(15)   =  0.0D0
      xi(16)   = -1.0D0
      xi(17)   = -1.0D0
      xi(18)   =  1.0D0
      xi(19)   =  1.0D0
      xi(20)   = -1.0D0
      xi(21)   =  0.0D0
      xi(22)   =  0.0D0
      xi(23)   =  0.0D0
      xi(24)   =  1.0D0
      xi(25)   =  0.0D0
      xi(26)   = -1.0D0
      xi(27)   =  0.0D0
      ! eta-coordinate at a node in a local element
      eta(1)   = -1.0D0
      eta(2)   = -1.0D0
      eta(3)   =  1.0D0
      eta(4)   =  1.0D0
      eta(5)   = -1.0D0
      eta(6)   = -1.0D0
      eta(7)   =  1.0D0
      eta(8)   =  1.0D0
      eta(9)   = -1.0D0
      eta(10)  =  0.0D0
      eta(11)  =  1.0D0
      eta(12)  =  0.0D0
      eta(13)  = -1.0D0
      eta(14)  =  0.0D0
      eta(15)  =  1.0D0
      eta(16)  =  0.0D0
      eta(17)  = -1.0D0
      eta(18)  = -1.0D0
      eta(19)  =  1.0D0
      eta(20)  =  1.0D0
      eta(21)  =  0.0D0
      eta(22)  =  0.0D0
      eta(23)  = -1.0D0
      eta(24)  =  0.0D0
      eta(25)  =  1.0D0
      eta(26)  =  0.0D0
      eta(27)  =  0.0D0
      ! zeta-coordinate at a node in a local element
      zeta(1)  = -1.0D0
      zeta(2)  = -1.0D0
      zeta(3)  = -1.0D0
      zeta(4)  = -1.0D0
      zeta(5)  =  1.0D0
      zeta(6)  =  1.0D0
      zeta(7)  =  1.0D0
      zeta(8)  =  1.0D0
      zeta(9)  = -1.0D0
      zeta(10) = -1.0D0
      zeta(11) = -1.0D0
      zeta(12) = -1.0D0
      zeta(13) =  1.0D0
      zeta(14) =  1.0D0
      zeta(15) =  1.0D0
      zeta(16) =  1.0D0
      zeta(17) =  0.0D0
      zeta(18) =  0.0D0
      zeta(19) =  0.0D0
      zeta(20) =  0.0D0
      zeta(21) = -1.0D0
      zeta(22) =  1.0D0
      zeta(23) =  0.0D0
      zeta(24) =  0.0D0
      zeta(25) =  0.0D0
      zeta(26) =  0.0D0
      zeta(27) =  0.0D0
      
!--------------------------------------------------------------------
      
      DO na = 1, 27
       
       n_xi(na)  =  ( 0.5D0*xi(na) *lcoord(1) )    &
                    *( 1.0D0+xi(na) *lcoord(1) )   &
                   +( 1.0D0-xi(na) *xi(na)  )      &
                    *( 1.0D0-lcoord(1)*lcoord(1) ) 
       n_eta(na) =  ( 0.5D0*eta(na)*lcoord(2) )    &
                    *( 1.0D0+eta(na)*lcoord(2) )   &
                   +( 1.0D0-eta(na)*eta(na) )      &
                    *( 1.0D0-lcoord(2)*lcoord(2) ) 
       n_zeta(na) = ( 0.5D0*zeta(na)*lcoord(3) )   &
                    *( 1.0D0+zeta(na)*lcoord(3) )  &
                   +( 1.0D0-zeta(na)*zeta(na) )    &
                    *( 1.0D0-lcoord(3)*lcoord(3) ) 
       
      END DO
      
!--------------------------------------------------------------------
      
      DO na = 1, 27
       
       func(na) = n_xi(na)*n_eta(na)*n_zeta(na)
       
      END DO
      
!--------------------------------------------------------------------
      
      RETURN
      
!####################################################################
      END SUBROUTINE ShapeFunc_hex27n
!####################################################################
      
      
!####################################################################
      SUBROUTINE ShapeDeriv_hex27n(lcoord, func)
!####################################################################
      
      REAL(KIND = kreal), INTENT(IN)  :: lcoord(3)
      REAL(KIND = kreal), INTENT(OUT) :: func(27, 3)
      
!--------------------------------------------------------------------
      
      REAL(KIND = kreal) :: xi(27), eta(27), zeta(27)
      REAL(KIND = kreal) :: n_xi(27), n_eta(27), n_zeta(27)
      REAL(KIND = kreal) :: dn_xi(27), dn_eta(27), dn_zeta(27)
      
      INTEGER :: na
      
!--------------------------------------------------------------------
      
      ! xi-coordinate at a node in a local element
      xi(1)    = -1.0D0
      xi(2)    =  1.0D0
      xi(3)    =  1.0D0
      xi(4)    = -1.0D0
      xi(5)    = -1.0D0
      xi(6)    =  1.0D0
      xi(7)    =  1.0D0
      xi(8)    = -1.0D0
      xi(9)    =  0.0D0
      xi(10)   =  1.0D0
      xi(11)   =  0.0D0
      xi(12)   = -1.0D0
      xi(13)   =  0.0D0
      xi(14)   =  1.0D0
      xi(15)   =  0.0D0
      xi(16)   = -1.0D0
      xi(17)   = -1.0D0
      xi(18)   =  1.0D0
      xi(19)   =  1.0D0
      xi(20)   = -1.0D0
      xi(21)   =  0.0D0
      xi(22)   =  0.0D0
      xi(23)   =  0.0D0
      xi(24)   =  1.0D0
      xi(25)   =  0.0D0
      xi(26)   = -1.0D0
      xi(27)   =  0.0D0
      ! eta-coordinate at a node in a local element
      eta(1)   = -1.0D0
      eta(2)   = -1.0D0
      eta(3)   =  1.0D0
      eta(4)   =  1.0D0
      eta(5)   = -1.0D0
      eta(6)   = -1.0D0
      eta(7)   =  1.0D0
      eta(8)   =  1.0D0
      eta(9)   = -1.0D0
      eta(10)  =  0.0D0
      eta(11)  =  1.0D0
      eta(12)  =  0.0D0
      eta(13)  = -1.0D0
      eta(14)  =  0.0D0
      eta(15)  =  1.0D0
      eta(16)  =  0.0D0
      eta(17)  = -1.0D0
      eta(18)  = -1.0D0
      eta(19)  =  1.0D0
      eta(20)  =  1.0D0
      eta(21)  =  0.0D0
      eta(22)  =  0.0D0
      eta(23)  = -1.0D0
      eta(24)  =  0.0D0
      eta(25)  =  1.0D0
      eta(26)  =  0.0D0
      eta(27)  =  0.0D0
      ! zeta-coordinate at a node in a local element
      zeta(1)  = -1.0D0
      zeta(2)  = -1.0D0
      zeta(3)  = -1.0D0
      zeta(4)  = -1.0D0
      zeta(5)  =  1.0D0
      zeta(6)  =  1.0D0
      zeta(7)  =  1.0D0
      zeta(8)  =  1.0D0
      zeta(9)  = -1.0D0
      zeta(10) = -1.0D0
      zeta(11) = -1.0D0
      zeta(12) = -1.0D0
      zeta(13) =  1.0D0
      zeta(14) =  1.0D0
      zeta(15) =  1.0D0
      zeta(16) =  1.0D0
      zeta(17) =  0.0D0
      zeta(18) =  0.0D0
      zeta(19) =  0.0D0
      zeta(20) =  0.0D0
      zeta(21) = -1.0D0
      zeta(22) =  1.0D0
      zeta(23) =  0.0D0
      zeta(24) =  0.0D0
      zeta(25) =  0.0D0
      zeta(26) =  0.0D0
      zeta(27) =  0.0D0
      
!--------------------------------------------------------------------
      
      DO na = 1, 27
       
       n_xi(na)  =  ( 0.5D0*xi(na) *lcoord(1) )    &
                    *( 1.0D0+xi(na) *lcoord(1) )   &
                   +( 1.0D0-xi(na) *xi(na)  )      &
                    *( 1.0D0-lcoord(1)*lcoord(1) ) 
       n_eta(na)  = ( 0.5D0*eta(na)*lcoord(2) )    &
                    *( 1.0D0+eta(na)*lcoord(2) )   &
                   +( 1.0D0-eta(na)*eta(na) )      &
                    *( 1.0D0-lcoord(2)*lcoord(2) ) 
       n_zeta(na) = ( 0.5D0*zeta(na)*lcoord(3) )   &
                    *( 1.0D0+zeta(na)*lcoord(3) )  &
                   +( 1.0D0-zeta(na)*zeta(na) )    &
                    *( 1.0D0-lcoord(3)*lcoord(3) ) 
       
       dn_xi(na)   = ( 0.5D0*xi(na)   )            &
                     *( 1.0D0+xi(na)  *lcoord(1) ) &
                    +( 0.5D0*xi(na)  *lcoord(1) )  &
                     *xi(na)                       &
                    +( 1.0D0-xi(na)  *xi(na)   )   &
                     *( -2.0D0*lcoord(1) )         
       dn_eta(na)  = ( 0.5D0*eta(na)  )            &
                     *( 1.0D0+eta(na) *lcoord(2) ) &
                    +( 0.5D0*eta(na) *lcoord(2) )  &
                     *eta(na)                      &
                    +( 1.0D0-eta(na) *eta(na)  )   &
                     *( -2.0D0*lcoord(2) )         
       dn_zeta(na) = ( 0.5D0*zeta(na) )            &
                     *( 1.0D0+zeta(na)*lcoord(3) ) &
                    +( 0.5D0*zeta(na)*lcoord(3) )  &
                     *zeta(na)                     &
                    +( 1.0D0-zeta(na)*zeta(na) )   &
                     *( -2.0D0*lcoord(3) )         
       
      END DO
      
!--------------------------------------------------------------------
      
      DO na = 1, 27
       
       func(na, 1) = dn_xi(na)*n_eta(na) *n_zeta(na)
       func(na, 2) = n_xi(na) *dn_eta(na)*n_zeta(na)
       func(na, 3) = n_xi(na) *n_eta(na) *dn_zeta(na)
       
      END DO
      
!--------------------------------------------------------------------
      
      RETURN
      
!####################################################################
      END SUBROUTINE ShapeDeriv_hex27n
!####################################################################
      
!####################################################################
      END MODULE shape_hex27n