!======================================================================!
!                                                                      !
! Software Name : FrontISTR Ver. 3.6                                   !
!                                                                      !
!      Module Name : lib                                               !
!                                                                      !
!                    Written by Gaku Hashimoto (Univ. of Tokyo)        !
!                                                                      !
!      Contact address :  IIS,The University of Tokyo, CISS            !
!                                                                      !
!      "Structural Analysis for Large Scale Assembly"                  !
!                                                                      !
!======================================================================!

      MODULE m_static_LIB_shell
!####################################################################

      USE hecmw, ONLY : kint, kreal
      USE elementInfo

!--------------------------------------------------------------------

      IMPLICIT NONE

!--------------------------------------------------------------------

      !--------------------------------------------------------------
      !
      ! (Programmer)
      ! Gaku Hashimoto
      ! Department of Human and Engineered Environmental Studies
      ! Graduate School of Frontier Sciences, The University of Tokyo
      ! 5-1-5 Kashiwanoha, Kashiwa, Chiba 277-8563 JAPAN
      !
      ! (Ref.)
      ! [1] Noguchi, H. and Hisada, T.,
      !     "Sensitivity analysis in post-buckling problems of shell
      !     structures,"
      !     Computers & Structures, Vol.47, No.4, pp.699-710, (1993).
      ! [2] Dvorkin, E.N. and Bathe, K.J.,
      !     "A Continuum Mechanics Based Four-node Shell Element for
      !     General Non-linear Analysis,"
      !     Engineering Computations, Vol.1, pp.77-88, (1984).
      ! [3] Bucalem, M.L. and Bathe, K.J.,
      !     "Higher-order MITC general shell element,"
      !     International Journal for Numerical Methods in
      !     Engineering, Vol.36, pp.3729-3754, (1993).
      ! [4] Lee, P.S. and Bathe, K.J.,
      !     "Development of MITC Isotropic Triangular Shell Finite
      !     Elements,"
      !     Computers & Structures, Vol.82, pp.945-962, (2004).
      !
      !--------------------------------------------------------------

!--------------------------------------------------------------------

      CONTAINS


!####################################################################
      SUBROUTINE STF_Shell_MITC                                           &
                 (etype, nn, ndof, ecoord, gausses, stiff, thick, mixflag, nddisp)
!####################################################################

      USE mMechGauss
      USE gauss_integration
      USE m_MatMatrix

!--------------------------------------------------------------------

      INTEGER(KIND = kint), INTENT(IN) :: etype
      INTEGER(KIND = kint), INTENT(IN) :: nn, mixflag
      INTEGER(KIND = kint), INTENT(IN) :: ndof
      REAL(KIND = kreal), INTENT(IN)   :: ecoord(3, nn)
      TYPE(tGaussStatus), INTENT(IN)   :: gausses(:)
      REAL(KIND = kreal), INTENT(OUT)  :: stiff(:, :)
      REAL(KIND = kreal), INTENT(IN)   :: thick

      REAL(KIND = kreal), INTENT(IN), OPTIONAL :: nddisp(3, nn)

!--------------------------------------------------------------------

      INTEGER :: flag, flag_dof, shellmatl
      INTEGER :: i, j, m
      INTEGER :: lx, ly
      INTEGER :: fetype
      INTEGER :: ny
      INTEGER :: ntying
      INTEGER :: npoints_tying(3)
      INTEGER :: it, ip
      INTEGER :: na, nb
      INTEGER :: isize, jsize
      INTEGER :: jsize1, jsize2, jsize3, &
                 jsize4, jsize5, jsize6
      INTEGER :: n_layer,n_total_layer, sstable(24)

      REAL(KIND = kreal) :: D(5, 5), B(5, ndof*nn), DB(5, ndof*nn)
      REAL(KIND = kreal) :: tmpstiff(ndof*nn, ndof*nn)
      REAL(KIND = kreal) :: elem(3, nn)
      REAL(KIND = kreal) :: unode(3, nn)
      REAL(KIND = kreal) :: xi_lx, eta_lx, zeta_ly
      REAL(KIND = kreal) :: w_w_lx, w_ly
      REAL(KIND = kreal) :: B_di(5, ndof*nn, 6, 3, 7)
      REAL(KIND = kreal) :: B1(3, ndof*nn), B2(3, ndof*nn), &
                            B3(3, ndof*nn)
      REAL(KIND = kreal) :: naturalcoord(2)
      REAL(KIND = kreal) :: tpcoord(6, 2, 3)
      REAL(KIND = kreal) :: nncoord(nn, 2)
      REAL(KIND = kreal) :: shapefunc(nn)
      REAL(KIND = kreal) :: shapederiv(nn, 2)
      REAL(KIND = kreal) :: aa1(3), aa2(3), aa3(3)
      REAL(KIND = kreal) :: bb1(3), bb2(3), bb3(3)
      REAL(KIND = kreal) :: cc1(3), cc2(3)
      REAL(KIND = kreal) :: alpha
      REAL(KIND = kreal) :: xxi_lx, eeta_lx
      REAL(KIND = kreal) :: xxi_di(6, 3), eeta_di(6, 3)
      REAL(KIND = kreal) :: h(nn, 3)
      REAL(KIND = kreal) :: v1(3, nn), v2(3, nn), v3(3, nn)
      REAL(KIND = kreal) :: v1_i(3), v2_i(3), v3_i(3)
      REAL(KIND = kreal) :: v1_abs, v2_abs, v3_abs
      REAL(KIND = kreal) :: a_over_2_v3(3, nn)
      REAL(KIND = kreal) :: u_rot(3, nn)
      REAL(KIND = kreal) :: dudxi_rot(3, nn), dudeta_rot(3, nn), &
                            dudzeta_rot(3, nn)
      REAL(KIND = kreal) :: g1(3), g2(3), g3(3)
      REAL(KIND = kreal) :: g1_abs, g2_abs, g3_abs
      REAL(KIND = kreal) :: e_0(3)
      REAL(KIND = kreal) :: cg1(3), cg2(3), cg3(3)
      REAL(KIND = kreal) :: det
      REAL(KIND = kreal) :: det_cg3(3)
      REAL(KIND = kreal) :: det_inv
      REAL(KIND = kreal) :: det_cg3_abs
      REAL(KIND = kreal) :: w_w_w_det
      REAL(KIND = kreal) :: e1_hat(3), e2_hat(3), e3_hat(3)
      REAL(KIND = kreal) :: e1_hat_abs, e2_hat_abs
      REAL(KIND = kreal) :: Cv12(ndof*nn), Cv13(ndof*nn), &
                            Cv21(ndof*nn), Cv23(ndof*nn), &
                            Cv31(ndof*nn), Cv32(ndof*nn)
      REAL(KIND = kreal) :: Cv_theta(ndof*nn), Cv_w(ndof*nn)
      REAL(KIND = kreal) :: Cv(ndof*nn)
      REAL(KIND = kreal) :: sigma_layer, thick_layer

		sstable = 0
		flag_dof = 0
!--------------------------------------------------------------------

      ! MITC4
      IF( etype .EQ. fe_mitc4_shell ) THEN

       fetype = fe_mitc4_shell

       ny = 2

       ntying = 1
       npoints_tying(1)= 4

      ! MITC9
      ELSE IF( etype .EQ. fe_mitc9_shell ) THEN

       fetype = fe_mitc9_shell

       ny = 3

       ntying = 3
       npoints_tying(1)= 6
       npoints_tying(2)= 6
       npoints_tying(3)= 4

      ! MITC3
      ELSE IF( etype .EQ. fe_mitc3_shell ) THEN

       fetype = fe_mitc3_shell

       ny = 2

       ntying = 1
       npoints_tying(1)= 3

      END IF

!--------------------------------------------------------------------

      IF( PRESENT( nddisp ) ) THEN

       unode(:, 1:nn) = nddisp(:, :)

      END IF

!--------------------------------------------------------------------

      flag = gausses(1)%pMaterial%nlgeom_flag

      IF( .NOT. PRESENT( nddisp ) ) flag = INFINITE

!--------------------------------------------------------------------

      elem(:, :) = ecoord(:, :)

!--------------------------------------------------------------------

      tmpstiff(:, :) = 0.0D0

!-------------------------------------------------------------------

      ! xi-coordinate at a node in a local element
      ! eta-coordinate at a node in a local element
      CALL getNodalNaturalCoord(fetype, nncoord)

!-------------------------------------------------------------------

      ! MITC4
      IF( etype .EQ. fe_mitc4_shell ) THEN

       !--------------------------------------------------------

       ! xi-coordinate at a tying pont in a local element
       tpcoord(1, 1, 1) =  0.0D0
       tpcoord(2, 1, 1) =  1.0D0
       tpcoord(3, 1, 1) =  0.0D0
       tpcoord(4, 1, 1) = -1.0D0
       ! eta-coordinate at a tying point in a local element
       tpcoord(1, 2, 1) = -1.0D0
       tpcoord(2, 2, 1) =  0.0D0
       tpcoord(3, 2, 1) =  1.0D0
       tpcoord(4, 2, 1) =  0.0D0

       !--------------------------------------------------------

      ! MITC9
      ELSE IF( etype .EQ. fe_mitc9_shell ) THEN

       !--------------------------------------------------------

       ! xi-coordinate at a tying point in a local element
       tpcoord(1, 1, 1) = -1.0D0*DSQRT( 1.0D0/3.0D0 )
       tpcoord(2, 1, 1) =  1.0D0*DSQRT( 1.0D0/3.0D0 )
       tpcoord(3, 1, 1) =  1.0D0*DSQRT( 1.0D0/3.0D0 )
       tpcoord(4, 1, 1) = -1.0D0*DSQRT( 1.0D0/3.0D0 )
       tpcoord(5, 1, 1) =  1.0D0*DSQRT( 1.0D0/3.0D0 )
       tpcoord(6, 1, 1) = -1.0D0*DSQRT( 1.0D0/3.0D0 )
       ! eta-coordinate at a tying point in a local element
       tpcoord(1, 2, 1) = -1.0D0*DSQRT( 3.0D0/5.0D0 )
       tpcoord(2, 2, 1) = -1.0D0*DSQRT( 3.0D0/5.0D0 )
       tpcoord(3, 2, 1) =  1.0D0*DSQRT( 3.0D0/5.0D0 )
       tpcoord(4, 2, 1) =  1.0D0*DSQRT( 3.0D0/5.0D0 )
       tpcoord(5, 2, 1) =  0.0D0*DSQRT( 3.0D0/5.0D0 )
       tpcoord(6, 2, 1) =  0.0D0*DSQRT( 3.0D0/5.0D0 )

       ! xi-coordinate at a tying point in a local element
       tpcoord(1, 1, 2) = -1.0D0*DSQRT( 3.0D0/5.0D0 )
       tpcoord(2, 1, 2) =  0.0D0*DSQRT( 3.0D0/5.0D0 )
       tpcoord(3, 1, 2) =  1.0D0*DSQRT( 3.0D0/5.0D0 )
       tpcoord(4, 1, 2) =  1.0D0*DSQRT( 3.0D0/5.0D0 )
       tpcoord(5, 1, 2) =  0.0D0*DSQRT( 3.0D0/5.0D0 )
       tpcoord(6, 1, 2) = -1.0D0*DSQRT( 3.0D0/5.0D0 )
       ! eta-coordinate at a tying point in a local element
       tpcoord(1, 2, 2) = -1.0D0*DSQRT( 1.0D0/3.0D0 )
       tpcoord(2, 2, 2) = -1.0D0*DSQRT( 1.0D0/3.0D0 )
       tpcoord(3, 2, 2) = -1.0D0*DSQRT( 1.0D0/3.0D0 )
       tpcoord(4, 2, 2) =  1.0D0*DSQRT( 1.0D0/3.0D0 )
       tpcoord(5, 2, 2) =  1.0D0*DSQRT( 1.0D0/3.0D0 )
       tpcoord(6, 2, 2) =  1.0D0*DSQRT( 1.0D0/3.0D0 )

       ! xi-coordinate at a tying point in a local element
       tpcoord(1, 1, 3) = -1.0D0*DSQRT( 1.0D0/3.0D0 )
       tpcoord(2, 1, 3) =  1.0D0*DSQRT( 1.0D0/3.0D0 )
       tpcoord(3, 1, 3) =  1.0D0*DSQRT( 1.0D0/3.0D0 )
       tpcoord(4, 1, 3) = -1.0D0*DSQRT( 1.0D0/3.0D0 )
       ! eta-coordinate at a tying point in a local element
       tpcoord(1, 2, 3) = -1.0D0*DSQRT( 1.0D0/3.0D0 )
       tpcoord(2, 2, 3) = -1.0D0*DSQRT( 1.0D0/3.0D0 )
       tpcoord(3, 2, 3) =  1.0D0*DSQRT( 1.0D0/3.0D0 )
       tpcoord(4, 2, 3) =  1.0D0*DSQRT( 1.0D0/3.0D0 )

       !--------------------------------------------------------

       ! Xi-coordinate at a tying point in a local element
       xxi_di(1, 1)  = -1.0D0
       xxi_di(2, 1)  =  1.0D0
       xxi_di(3, 1)  =  1.0D0
       xxi_di(4, 1)  = -1.0D0
       xxi_di(5, 1)  =  1.0D0
       xxi_di(6, 1)  = -1.0D0
       ! Eta-coordinate at a tying point in a local element
       eeta_di(1, 1) = -1.0D0
       eeta_di(2, 1) = -1.0D0
       eeta_di(3, 1) =  1.0D0
       eeta_di(4, 1) =  1.0D0
       eeta_di(5, 1) =  0.0D0
       eeta_di(6, 1) =  0.0D0

       ! Xi-coordinate at a tying point in a local element
       xxi_di(1, 2)  = -1.0D0
       xxi_di(2, 2)  =  0.0D0
       xxi_di(3, 2)  =  1.0D0
       xxi_di(4, 2)  =  1.0D0
       xxi_di(5, 2)  =  0.0D0
       xxi_di(6, 2)  = -1.0D0
       ! Eta-coordinate at a tying point in a local element
       eeta_di(1, 2) = -1.0D0
       eeta_di(2, 2) = -1.0D0
       eeta_di(3, 2) = -1.0D0
       eeta_di(4, 2) =  1.0D0
       eeta_di(5, 2) =  1.0D0
       eeta_di(6, 2) =  1.0D0

       !--------------------------------------------------------

      ! MITC3
      ELSE IF( etype .EQ. fe_mitc3_shell ) THEN

       !--------------------------------------------------------

       ! xi-coordinate at a tying point in a local element
       tpcoord(1, 1, 1)  = 0.5D0
       tpcoord(2, 1, 1)  = 0.0D0
       tpcoord(3, 1, 1)  = 0.5D0
       ! eta-coordinate at a tying point in a local element
       tpcoord(1, 2, 1) = 0.0D0
       tpcoord(2, 2, 1) = 0.5D0
       tpcoord(3, 2, 1) = 0.5D0

       !--------------------------------------------------------

      END IF

!--------------------------------------------------------------------

      ! xi-coordinate at the center point in a local element
      ! eta-coordinate at the center point in a local element
      naturalcoord(1) = 0.0D0
      naturalcoord(2) = 0.0D0

      CALL getShapeDeriv(fetype, naturalcoord, shapederiv)

      !--------------------------------------------------------------

      ! Covariant basis vector
      DO i = 1, 3

       g1(i) = 0.0D0

       DO na = 1, nn

        g1(i) = g1(i)+shapederiv(na, 1) &
                      *elem(i, na)

       END DO

      END DO

      e_0(1) = g1(1)
      e_0(2) = g1(2)
      e_0(3) = g1(3)

      !--------------------------------------------------------------

      DO nb = 1, nn

       !--------------------------------------------------------

       naturalcoord(1) = nncoord(nb, 1)
       naturalcoord(2) = nncoord(nb, 2)

       CALL getShapeDeriv(fetype, naturalcoord, shapederiv)

       !--------------------------------------------------------

       ! Covariant basis vector
       DO i = 1, 3

        g1(i) = 0.0D0
        g2(i) = 0.0D0

        DO na = 1, nn

         g1(i) = g1(i)+shapederiv(na, 1) &
                       *elem(i, na)
         g2(i) = g2(i)+shapederiv(na, 2) &
                       *elem(i, na)

        END DO

       END DO

       !------------------------------------------

       det_cg3(1) = g1(2)*g2(3)-g1(3)*g2(2)
       det_cg3(2) = g1(3)*g2(1)-g1(1)*g2(3)
       det_cg3(3) = g1(1)*g2(2)-g1(2)*g2(1)

       det_cg3_abs = DSQRT( det_cg3(1)*det_cg3(1)   &
                           +det_cg3(2)*det_cg3(2)   &
                           +det_cg3(3)*det_cg3(3) )

       v3(1, nb) = det_cg3(1)/det_cg3_abs
       v3(2, nb) = det_cg3(2)/det_cg3_abs
       v3(3, nb) = det_cg3(3)/det_cg3_abs

       !--------------------------------------------------------

       v2(1, nb) = v3(2, nb)*e_0(3)-v3(3, nb)*e_0(2)
       v2(2, nb) = v3(3, nb)*e_0(1)-v3(1, nb)*e_0(3)
       v2(3, nb) = v3(1, nb)*e_0(2)-v3(2, nb)*e_0(1)

       v2_abs = DSQRT( v2(1, nb)*v2(1, nb)   &
                      +v2(2, nb)*v2(2, nb)   &
                      +v2(3, nb)*v2(3, nb) )

       IF( v2_abs .GT. 1.0D-15 ) THEN

        v2(1, nb) = v2(1, nb)/v2_abs
        v2(2, nb) = v2(2, nb)/v2_abs
        v2(3, nb) = v2(3, nb)/v2_abs

        v1(1, nb) = v2(2, nb)*v3(3, nb) &
                   -v2(3, nb)*v3(2, nb)
        v1(2, nb) = v2(3, nb)*v3(1, nb) &
                   -v2(1, nb)*v3(3, nb)
        v1(3, nb) = v2(1, nb)*v3(2, nb) &
                   -v2(2, nb)*v3(1, nb)

        v1_abs = DSQRT( v1(1, nb)*v1(1, nb)   &
                       +v1(2, nb)*v1(2, nb)   &
                       +v1(3, nb)*v1(3, nb) )

        v1(1, nb) = v1(1, nb)/v1_abs
        v1(2, nb) = v1(2, nb)/v1_abs
        v1(3, nb) = v1(3, nb)/v1_abs

       ELSE    ! YX: impossible

        v1(1, nb) =  0.0D0
        v1(2, nb) =  0.0D0
        v1(3, nb) = -1.0D0

        v2(1, nb) = 0.0D0
        v2(2, nb) = 1.0D0
        v2(3, nb) = 0.0D0

       END IF

       !---------------------------------------------------

       v3(1, nb) = v1(2, nb)*v2(3, nb) &
                  -v1(3, nb)*v2(2, nb)
       v3(2, nb) = v1(3, nb)*v2(1, nb) &
                  -v1(1, nb)*v2(3, nb)
       v3(3, nb) = v1(1, nb)*v2(2, nb) &
                  -v1(2, nb)*v2(1, nb)

       v3_abs = DSQRT( v3(1, nb)*v3(1, nb)   &
                      +v3(2, nb)*v3(2, nb)   &
                      +v3(3, nb)*v3(3, nb) )

       v3(1, nb) = v3(1, nb)/v3_abs
       v3(2, nb) = v3(2, nb)/v3_abs
       v3(3, nb) = v3(3, nb)/v3_abs

       !--------------------------------------------------------

       a_over_2_v3(1, nb) = 0.5D0*thick*v3(1, nb)
       a_over_2_v3(2, nb) = 0.5D0*thick*v3(2, nb)
       a_over_2_v3(3, nb) = 0.5D0*thick*v3(3, nb)

       !--------------------------------------------------------

      END DO

!--------------------------------------------------------------------
!		MODIFIED to LAMINATED SHELL ANALYSIS
!--------------------------------------------------------------------

   n_total_layer =  int(gausses(1)%pMaterial%variables(M_TOTAL_LAYER))
   DO n_layer=1,n_total_layer
      DO ly = 1, ny

       !--------------------------------------------------------

       ! MITC4
       IF( etype .EQ. fe_mitc4_shell ) THEN

        zeta_ly = 0.0D0

       ! MITC9
       ELSE IF( etype .EQ. fe_mitc9_shell ) THEN

        zeta_ly = xg(ny, ly)

       ! MITC3
       ELSE IF( etype .EQ. fe_mitc3_shell )THEN

        zeta_ly = 0.0D0

       END IF

       !--------------------------------------------------------

       DO it = 1, ntying

        DO ip = 1, npoints_tying(it)

         !-------------------------------------------------

         naturalcoord(1) = tpcoord(ip, 1, it)
         naturalcoord(2) = tpcoord(ip, 2, it)

         CALL getShapeFunc(fetype, naturalcoord, shapefunc)

         CALL getShapeDeriv(fetype, naturalcoord, shapederiv)

         !-------------------------------------------------

         DO na = 1, nn

          DO i = 1, 3

           u_rot(i, na)                      &
           = shapefunc(na)                   &
             *( zeta_ly*a_over_2_v3(i, na) )

           dudxi_rot(i, na)                  &
           = shapederiv(na, 1)               &
             *( zeta_ly*a_over_2_v3(i, na) )
           dudeta_rot(i, na)                 &
           = shapederiv(na, 2)               &
             *( zeta_ly*a_over_2_v3(i, na) )
           dudzeta_rot(i, na)                &
           = shapefunc(na)                   &
             *( a_over_2_v3(i, na) )

          END DO

         END DO

         !-------------------------------------------------

         ! Covariant basis vector
         DO i = 1, 3

          g1(i) = 0.0D0
          g2(i) = 0.0D0
          g3(i) = 0.0D0

          DO na = 1, nn

           g1(i) = g1(i)+shapederiv(na, 1)  &
                         *elem(i, na)       &
                        +dudxi_rot(i, na)
           g2(i) = g2(i)+shapederiv(na, 2)  &
                         *elem(i, na)       &
                        +dudeta_rot(i, na)
           g3(i) = g3(i)+dudzeta_rot(i, na)

          END DO

         END DO

         !-------------------------------------------------

         ! [ B L ] matrix
         DO nb = 1, nn

          jsize1 = ndof*(nb-1)+1
          jsize2 = ndof*(nb-1)+2
          jsize3 = ndof*(nb-1)+3
          jsize4 = ndof*(nb-1)+4
          jsize5 = ndof*(nb-1)+5
          jsize6 = ndof*(nb-1)+6

          aa1(1) = dudxi_rot(2, nb)  *g1(3)-dudxi_rot(3, nb)  *g1(2)
          aa1(2) = dudxi_rot(3, nb)  *g1(1)-dudxi_rot(1, nb)  *g1(3)
          aa1(3) = dudxi_rot(1, nb)  *g1(2)-dudxi_rot(2, nb)  *g1(1)

          aa2(1) = dudxi_rot(2, nb)  *g2(3)-dudxi_rot(3, nb)  *g2(2)
          aa2(2) = dudxi_rot(3, nb)  *g2(1)-dudxi_rot(1, nb)  *g2(3)
          aa2(3) = dudxi_rot(1, nb)  *g2(2)-dudxi_rot(2, nb)  *g2(1)

          aa3(1) = dudxi_rot(2, nb)  *g3(3)-dudxi_rot(3, nb)  *g3(2)
          aa3(2) = dudxi_rot(3, nb)  *g3(1)-dudxi_rot(1, nb)  *g3(3)
          aa3(3) = dudxi_rot(1, nb)  *g3(2)-dudxi_rot(2, nb)  *g3(1)

          bb1(1) = dudeta_rot(2, nb) *g1(3)-dudeta_rot(3, nb) *g1(2)
          bb1(2) = dudeta_rot(3, nb) *g1(1)-dudeta_rot(1, nb) *g1(3)
          bb1(3) = dudeta_rot(1, nb) *g1(2)-dudeta_rot(2, nb) *g1(1)

          bb2(1) = dudeta_rot(2, nb) *g2(3)-dudeta_rot(3, nb) *g2(2)
          bb2(2) = dudeta_rot(3, nb) *g2(1)-dudeta_rot(1, nb) *g2(3)
          bb2(3) = dudeta_rot(1, nb) *g2(2)-dudeta_rot(2, nb) *g2(1)

          bb3(1) = dudeta_rot(2, nb) *g3(3)-dudeta_rot(3, nb) *g3(2)
          bb3(2) = dudeta_rot(3, nb) *g3(1)-dudeta_rot(1, nb) *g3(3)
          bb3(3) = dudeta_rot(1, nb) *g3(2)-dudeta_rot(2, nb) *g3(1)

          cc1(1) = dudzeta_rot(2, nb)*g1(3)-dudzeta_rot(3, nb)*g1(2)
          cc1(2) = dudzeta_rot(3, nb)*g1(1)-dudzeta_rot(1, nb)*g1(3)
          cc1(3) = dudzeta_rot(1, nb)*g1(2)-dudzeta_rot(2, nb)*g1(1)

          cc2(1) = dudzeta_rot(2, nb)*g2(3)-dudzeta_rot(3, nb)*g2(2)
          cc2(2) = dudzeta_rot(3, nb)*g2(1)-dudzeta_rot(1, nb)*g2(3)
          cc2(3) = dudzeta_rot(1, nb)*g2(2)-dudzeta_rot(2, nb)*g2(1)

          B_di(1, jsize1, ip, it, ly) = shapederiv(nb, 1)*g1(1)
          B_di(2, jsize1, ip, it, ly) = shapederiv(nb, 2)*g2(1)
          B_di(3, jsize1, ip, it, ly) = shapederiv(nb, 1)*g2(1) &
                                       +shapederiv(nb, 2)*g1(1)
          B_di(4, jsize1, ip, it, ly) = shapederiv(nb, 2)*g3(1)
          B_di(5, jsize1, ip, it, ly) = shapederiv(nb, 1)*g3(1)

          B_di(1, jsize2, ip, it, ly) = shapederiv(nb, 1)*g1(2)
          B_di(2, jsize2, ip, it, ly) = shapederiv(nb, 2)*g2(2)
          B_di(3, jsize2, ip, it, ly) = shapederiv(nb, 1)*g2(2) &
                                       +shapederiv(nb, 2)*g1(2)
          B_di(4, jsize2, ip, it, ly) = shapederiv(nb, 2)*g3(2)
          B_di(5, jsize2, ip, it, ly) = shapederiv(nb, 1)*g3(2)

          B_di(1, jsize3, ip, it, ly) = shapederiv(nb, 1)*g1(3)
          B_di(2, jsize3, ip, it, ly) = shapederiv(nb, 2)*g2(3)
          B_di(3, jsize3, ip, it, ly) = shapederiv(nb, 1)*g2(3) &
                                       +shapederiv(nb, 2)*g1(3)
          B_di(4, jsize3, ip, it, ly) = shapederiv(nb, 2)*g3(3)
          B_di(5, jsize3, ip, it, ly) = shapederiv(nb, 1)*g3(3)

          B_di(1, jsize4, ip, it, ly) = aa1(1)
          B_di(2, jsize4, ip, it, ly) = bb2(1)
          B_di(3, jsize4, ip, it, ly) = aa2(1)+bb1(1)
          B_di(4, jsize4, ip, it, ly) = bb3(1)+cc2(1)
          B_di(5, jsize4, ip, it, ly) = aa3(1)+cc1(1)

          B_di(1, jsize5, ip, it, ly) = aa1(2)
          B_di(2, jsize5, ip, it, ly) = bb2(2)
          B_di(3, jsize5, ip, it, ly) = aa2(2)+bb1(2)
          B_di(4, jsize5, ip, it, ly) = bb3(2)+cc2(2)
          B_di(5, jsize5, ip, it, ly) = aa3(2)+cc1(2)

          B_di(1, jsize6, ip, it, ly) = aa1(3)
          B_di(2, jsize6, ip, it, ly) = bb2(3)
          B_di(3, jsize6, ip, it, ly) = aa2(3)+bb1(3)
          B_di(4, jsize6, ip, it, ly) = bb3(3)+cc2(3)
          B_di(5, jsize6, ip, it, ly) = aa3(3)+cc1(3)

         END DO

         !-------------------------------------------------

        END DO

       END DO

       !--------------------------------------------------------

        IF ( n_total_layer .ge. 2 ) THEN
          shellmatl = int(gausses(1)%pMaterial%variables(M_SHELL_MATLTYPE))
          if (shellmatl == 0)then
            sigma_layer = 0.0D0
            DO m = 1, n_layer
              sigma_layer = sigma_layer + 2 * gausses(1)%pMaterial%variables(100+3*m) / thick
            END DO
            zeta_ly = -1 + sigma_layer - gausses(1)%pMaterial%variables(100+3*n_layer) / thick * (1-xg(ny, ly))
          elseif (shellmatl == 1)then
            sigma_layer = 0.0D0
            DO m = 1, n_layer
              sigma_layer = sigma_layer + 2 * gausses(1)%pMaterial%variables(100+8*m-5) / thick
            END DO
            zeta_ly = -1 + sigma_layer - gausses(1)%pMaterial%variables(100+8*n_layer-5) / thick * (1-xg(ny, ly))
          else
            write(*,*)"ERROR : shellmatl isnot correct"; stop
          endif
        ELSE
          zeta_ly = xg(ny, ly)
        ENDIF

        w_ly    = wgt(ny, ly)

       !--------------------------------------------------------

       DO lx = 1, NumOfQuadPoints(fetype)

        !--------------------------------------------------

        CALL getQuadPoint(fetype, lx, naturalcoord)

        xi_lx  = naturalcoord(1)
        eta_lx = naturalcoord(2)

        w_w_lx = getWeight(fetype, lx)

        CALL getShapeFunc(fetype, naturalcoord, shapefunc)

        CALL getShapeDeriv(fetype, naturalcoord, shapederiv)

        !--------------------------------------------------

        DO i = 1, 3

         v1_i(i) = 0.0D0
         v2_i(i) = 0.0D0
         v3_i(i) = 0.0D0

         DO na = 1, nn

          v1_i(i) = v1_i(i)+shapefunc(na)*v1(i, na)
          v2_i(i) = v2_i(i)+shapefunc(na)*v2(i, na)
          v3_i(i) = v3_i(i)+shapefunc(na)*v3(i, na)

         END DO

        END DO

        !--------------------------------------------------

        DO na = 1, nn

         DO i = 1, 3

          u_rot(i, na)                      &
          = shapefunc(na)                   &
            *( zeta_ly*a_over_2_v3(i, na) )

          dudxi_rot(i, na)                  &
          = shapederiv(na, 1)               &
            *( zeta_ly*a_over_2_v3(i, na) )
          dudeta_rot(i, na)                 &
          = shapederiv(na, 2)               &
            *( zeta_ly*a_over_2_v3(i, na) )
          dudzeta_rot(i, na)                &
          = shapefunc(na)                   &
            *( a_over_2_v3(i, na) )

         END DO

        END DO

        !--------------------------------------------------

        ! Covariant basis vector
        DO i = 1, 3

         g1(i) = 0.0D0
         g2(i) = 0.0D0
         g3(i) = 0.0D0

         DO na = 1, nn

          g1(i) = g1(i)+shapederiv(na, 1)  &
                        *elem(i, na)       &
                       +dudxi_rot(i, na)
          g2(i) = g2(i)+shapederiv(na, 2)  &
                        *elem(i, na)       &
                       +dudeta_rot(i, na)
          g3(i) = g3(i)+dudzeta_rot(i, na)

         END DO

        END DO

        !--------------------------------------------------

        ! Jacobian
        det = g1(1)*( g2(2)*g3(3)-g2(3)*g3(2) ) &
             +g1(2)*( g2(3)*g3(1)-g2(1)*g3(3) ) &
             +g1(3)*( g2(1)*g3(2)-g2(2)*g3(1) )

        det_inv = 1.0D0/det

        !--------------------------------------------------

        ! Contravariant basis vector
        cg1(1) = det_inv                      &
                 *( g2(2)*g3(3)-g2(3)*g3(2) )
        cg1(2) = det_inv                      &
                 *( g2(3)*g3(1)-g2(1)*g3(3) )
        cg1(3) = det_inv                      &
                 *( g2(1)*g3(2)-g2(2)*g3(1) )
        cg2(1) = det_inv                      &
                 *( g3(2)*g1(3)-g3(3)*g1(2) )
        cg2(2) = det_inv                      &
                 *( g3(3)*g1(1)-g3(1)*g1(3) )
        cg2(3) = det_inv                      &
                 *( g3(1)*g1(2)-g3(2)*g1(1) )
        cg3(1) = det_inv                      &
                 *( g1(2)*g2(3)-g1(3)*g2(2) )
        cg3(2) = det_inv                      &
                 *( g1(3)*g2(1)-g1(1)*g2(3) )
        cg3(3) = det_inv                      &
                 *( g1(1)*g2(2)-g1(2)*g2(1) )

        !--------------------------------------------------

        g3_abs = DSQRT( g3(1)*g3(1)   &
                       +g3(2)*g3(2)   &
                       +g3(3)*g3(3) )

        !--------------------------------------------------

        ! Orthonormal vectors

        e3_hat(1) = g3(1)/g3_abs
        e3_hat(2) = g3(2)/g3_abs
        e3_hat(3) = g3(3)/g3_abs

        e1_hat(1) = g2(2)*e3_hat(3) &
                   -g2(3)*e3_hat(2)
        e1_hat(2) = g2(3)*e3_hat(1) &
                   -g2(1)*e3_hat(3)
        e1_hat(3) = g2(1)*e3_hat(2) &
                   -g2(2)*e3_hat(1)
        e1_hat_abs = DSQRT( e1_hat(1)*e1_hat(1)   &
                           +e1_hat(2)*e1_hat(2)   &
                           +e1_hat(3)*e1_hat(3) )
        e1_hat(1) = e1_hat(1)/e1_hat_abs
        e1_hat(2) = e1_hat(2)/e1_hat_abs
        e1_hat(3) = e1_hat(3)/e1_hat_abs

        e2_hat(1) = e3_hat(2)*e1_hat(3) &
                   -e3_hat(3)*e1_hat(2)
        e2_hat(2) = e3_hat(3)*e1_hat(1) &
                   -e3_hat(1)*e1_hat(3)
        e2_hat(3) = e3_hat(1)*e1_hat(2) &
                   -e3_hat(2)*e1_hat(1)
        e2_hat_abs = DSQRT( e2_hat(1)*e2_hat(1)   &
                           +e2_hat(2)*e2_hat(2)   &
                           +e2_hat(3)*e2_hat(3) )
        e2_hat(1) = e2_hat(1)/e2_hat_abs
        e2_hat(2) = e2_hat(2)/e2_hat_abs
        e2_hat(3) = e2_hat(3)/e2_hat_abs

        !--------------------------------------------------

        CALL MatlMatrix_Shell                        &
             (gausses(lx), Shell, D,                 &
              e1_hat, e2_hat, e3_hat, cg1, cg2, cg3, &
              alpha, n_layer)

        !--------------------------------------------------

        ! [ B L ] matrix
        DO nb = 1, nn

         jsize1 = ndof*(nb-1)+1
         jsize2 = ndof*(nb-1)+2
         jsize3 = ndof*(nb-1)+3
         jsize4 = ndof*(nb-1)+4
         jsize5 = ndof*(nb-1)+5
         jsize6 = ndof*(nb-1)+6

         aa1(1) = dudxi_rot(2, nb)  *g1(3)-dudxi_rot(3, nb)  *g1(2)
         aa1(2) = dudxi_rot(3, nb)  *g1(1)-dudxi_rot(1, nb)  *g1(3)
         aa1(3) = dudxi_rot(1, nb)  *g1(2)-dudxi_rot(2, nb)  *g1(1)

         aa2(1) = dudxi_rot(2, nb)  *g2(3)-dudxi_rot(3, nb)  *g2(2)
         aa2(2) = dudxi_rot(3, nb)  *g2(1)-dudxi_rot(1, nb)  *g2(3)
         aa2(3) = dudxi_rot(1, nb)  *g2(2)-dudxi_rot(2, nb)  *g2(1)

         aa3(1) = dudxi_rot(2, nb)  *g3(3)-dudxi_rot(3, nb)  *g3(2)
         aa3(2) = dudxi_rot(3, nb)  *g3(1)-dudxi_rot(1, nb)  *g3(3)
         aa3(3) = dudxi_rot(1, nb)  *g3(2)-dudxi_rot(2, nb)  *g3(1)

         bb1(1) = dudeta_rot(2, nb) *g1(3)-dudeta_rot(3, nb) *g1(2)
         bb1(2) = dudeta_rot(3, nb) *g1(1)-dudeta_rot(1, nb) *g1(3)
         bb1(3) = dudeta_rot(1, nb) *g1(2)-dudeta_rot(2, nb) *g1(1)

         bb2(1) = dudeta_rot(2, nb) *g2(3)-dudeta_rot(3, nb) *g2(2)
         bb2(2) = dudeta_rot(3, nb) *g2(1)-dudeta_rot(1, nb) *g2(3)
         bb2(3) = dudeta_rot(1, nb) *g2(2)-dudeta_rot(2, nb) *g2(1)

         bb3(1) = dudeta_rot(2, nb) *g3(3)-dudeta_rot(3, nb) *g3(2)
         bb3(2) = dudeta_rot(3, nb) *g3(1)-dudeta_rot(1, nb) *g3(3)
         bb3(3) = dudeta_rot(1, nb) *g3(2)-dudeta_rot(2, nb) *g3(1)

         cc1(1) = dudzeta_rot(2, nb)*g1(3)-dudzeta_rot(3, nb)*g1(2)
         cc1(2) = dudzeta_rot(3, nb)*g1(1)-dudzeta_rot(1, nb)*g1(3)
         cc1(3) = dudzeta_rot(1, nb)*g1(2)-dudzeta_rot(2, nb)*g1(1)

         cc2(1) = dudzeta_rot(2, nb)*g2(3)-dudzeta_rot(3, nb)*g2(2)
         cc2(2) = dudzeta_rot(3, nb)*g2(1)-dudzeta_rot(1, nb)*g2(3)
         cc2(3) = dudzeta_rot(1, nb)*g2(2)-dudzeta_rot(2, nb)*g2(1)

         B(1, jsize1) = shapederiv(nb, 1)*g1(1)
         B(2, jsize1) = shapederiv(nb, 2)*g2(1)
         B(3, jsize1) = shapederiv(nb, 1)*g2(1) &
                       +shapederiv(nb, 2)*g1(1)
         B(4, jsize1) = shapederiv(nb, 2)*g3(1)
         B(5, jsize1) = shapederiv(nb, 1)*g3(1)

         B(1, jsize2) = shapederiv(nb, 1)*g1(2)
         B(2, jsize2) = shapederiv(nb, 2)*g2(2)
         B(3, jsize2) = shapederiv(nb, 1)*g2(2) &
                       +shapederiv(nb, 2)*g1(2)
         B(4, jsize2) = shapederiv(nb, 2)*g3(2)
         B(5, jsize2) = shapederiv(nb, 1)*g3(2)

         B(1, jsize3) = shapederiv(nb, 1)*g1(3)
         B(2, jsize3) = shapederiv(nb, 2)*g2(3)
         B(3, jsize3) = shapederiv(nb, 1)*g2(3) &
                       +shapederiv(nb, 2)*g1(3)
         B(4, jsize3) = shapederiv(nb, 2)*g3(3)
         B(5, jsize3) = shapederiv(nb, 1)*g3(3)

         B(1, jsize4) = aa1(1)
         B(2, jsize4) = bb2(1)
         B(3, jsize4) = aa2(1)+bb1(1)
         B(4, jsize4) = bb3(1)+cc2(1)
         B(5, jsize4) = aa3(1)+cc1(1)

         B(1, jsize5) = aa1(2)
         B(2, jsize5) = bb2(2)
         B(3, jsize5) = aa2(2)+bb1(2)
         B(4, jsize5) = bb3(2)+cc2(2)
         B(5, jsize5) = aa3(2)+cc1(2)

         B(1, jsize6) = aa1(3)
         B(2, jsize6) = bb2(3)
         B(3, jsize6) = aa2(3)+bb1(3)
         B(4, jsize6) = bb3(3)+cc2(3)
         B(5, jsize6) = aa3(3)+cc1(3)

        END DO

        !--------------------------------------------------

        ! MITC4
        IF( etype .EQ. fe_mitc4_shell ) THEN

         DO jsize = 1, ndof*nn

          B(4, jsize) = 0.0D0
          B(5, jsize) = 0.0D0

          ! B_as(4, jsize)
          B(4, jsize)                                       &
          = 0.5D0*( 1.0D0-xi_lx  )*B_di(4, jsize, 4, 1, ly) &
           +0.5D0*( 1.0D0+xi_lx  )*B_di(4, jsize, 2, 1, ly)
          ! B_as(5, jsize)
          B(5, jsize)                                       &
          = 0.5D0*( 1.0D0-eta_lx )*B_di(5, jsize, 1, 1, ly) &
           +0.5D0*( 1.0D0+eta_lx )*B_di(5, jsize, 3, 1, ly)

         END DO

        ! MITC9
        ELSE IF( etype .EQ. fe_mitc9_shell ) THEN

         xxi_lx  = xi_lx /DSQRT( 1.0D0/3.0D0 )
         eeta_lx = eta_lx/DSQRT( 3.0D0/5.0D0 )

         DO ip = 1, npoints_tying(1)

          h(ip, 1)                                     &
          = ( 0.5D0*( 1.0D0+xxi_di(ip, 1)*xxi_lx ) )   &
            *( ( 0.5D0*eeta_di(ip, 1)*eeta_lx )        &
               *( 1.0D0+eeta_di(ip, 1)*eeta_lx )       &
              +( 1.0D0-eeta_di(ip, 1)*eeta_di(ip, 1) ) &
               *( 1.0D0-eeta_lx*eeta_lx ) )

         END DO

         xxi_lx  = xi_lx /DSQRT( 3.0D0/5.0D0 )
         eeta_lx = eta_lx/DSQRT( 1.0D0/3.0D0 )

         DO ip = 1, npoints_tying(2)

          h(ip, 2)                                      &
          = ( ( 0.5D0*xxi_di(ip, 2) *xxi_lx  )          &
              *( 1.0D0+xxi_di(ip, 2) *xxi_lx  )         &
             +( 1.0D0-xxi_di(ip, 2) *xxi_di(ip, 2)  )   &
              *( 1.0D0-xxi_lx*xxi_lx ) )                &
            *( 0.5D0*( 1.0D0+eeta_di(ip, 2)*eeta_lx ) )

         END DO

         xxi_lx  = xi_lx /DSQRT( 1.0D0/3.0D0 )
         eeta_lx = eta_lx/DSQRT( 1.0D0/3.0D0 )

         DO ip = 1, npoints_tying(3)

          h(ip, 3)                                      &
          = ( 0.5D0*( 1.0D0+xxi_di(ip, 1)*xxi_lx ) )    &
            *( 0.5D0*( 1.0D0+eeta_di(ip, 1)*eeta_lx ) )

         END DO

         DO jsize = 1, ndof*nn

          B(1, jsize) = 0.0D0
          B(2, jsize) = 0.0D0
          B(3, jsize) = 0.0D0
          B(4, jsize) = 0.0D0
          B(5, jsize) = 0.0D0

          DO ip = 1, npoints_tying(1)

           ! B_as(1, jsize)
           B(1, jsize)                                      &
           = B(1, jsize)+h(ip, 1)*B_di(1, jsize, ip, 1, ly)
           ! B_as(5, jsize)
           B(5, jsize)                                      &
           = B(5, jsize)+h(ip, 1)*B_di(5, jsize, ip, 1, ly)

          END DO

          DO ip = 1, npoints_tying(2)

           ! B_as(2, jsize)
           B(2, jsize)                                      &
           = B(2, jsize)+h(ip, 2)*B_di(2, jsize, ip, 2, ly)
           ! B_as(4, jsize)
           B(4, jsize)                                      &
           = B(4, jsize)+h(ip, 2)*B_di(4, jsize, ip, 2, ly)

          END DO

          DO ip = 1, npoints_tying(3)

           ! B_as(3, jsize)
           B(3, jsize)                                      &
           = B(3, jsize)+h(ip, 3)*B_di(3, jsize, ip, 3, ly)

          END DO

         END DO

        ! MITC3
        ELSE IF( etype .EQ. fe_mitc3_shell ) THEN

         DO jsize = 1, ndof*nn

          B(4, jsize) = 0.0D0
          B(5, jsize) = 0.0D0

          ! B_as(4, jsize)
          B(4, jsize)                                 &
          = ( 1.0D0-xi_lx  )*B_di(4, jsize, 2, 1, ly) &
           +xi_lx *B_di(5, jsize, 1, 1, ly)           &
           +xi_lx *( B_di(4, jsize, 3, 1, ly)         &
                    -B_di(5, jsize, 3, 1, ly)  )

          ! B_as(5, jsize)
          B(5, jsize)                                 &
          = eta_lx*B_di(4, jsize, 2, 1, ly)           &
           +( 1.0D0-eta_lx )*B_di(5, jsize, 1, 1, ly) &
           -eta_lx*( B_di(4, jsize, 3, 1, ly)         &
                    -B_di(5, jsize, 3, 1, ly)  )

         END DO

        END IF

        !--------------------------------------------------

        w_w_w_det = w_w_lx*w_ly*det

        !--------------------------------------------------

        DB(1:5, 1:ndof*nn) = MATMUL( D, B(1:5, 1:ndof*nn ) )

        !--------------------------------------------------

        shellmatl = int(gausses(1)%pMaterial%variables(M_SHELL_MATLTYPE))
        if (shellmatl == 0)then
          if ( n_total_layer .ge. 2 ) THEN
            FORALL( isize=1:ndof*nn, jsize=1:ndof*nn )
              tmpstiff(isize, jsize)                               &
              = tmpstiff(isize, jsize)                             &
              +w_w_w_det*gausses(1)%pMaterial%variables(100+3*n_layer)/thick*DOT_PRODUCT( B(:, isize), DB(:, jsize) )
            END FORALL
          else
            FORALL( isize=1:ndof*nn, jsize=1:ndof*nn )
              tmpstiff(isize, jsize)                               &
              = tmpstiff(isize, jsize) +w_w_w_det *DOT_PRODUCT( B(:, isize), DB(:, jsize) )
            END FORALL
          endif
        elseif (shellmatl == 1)then
          FORALL( isize=1:ndof*nn, jsize=1:ndof*nn )
            tmpstiff(isize, jsize)                               &
            = tmpstiff(isize, jsize)                             &
            +w_w_w_det*gausses(1)%pMaterial%variables(100+8*n_layer-5)/thick*DOT_PRODUCT( B(:, isize), DB(:, jsize) )
          END FORALL
        else
          write(*,*)"ERROR : shellmatl isnot correct"; stop
        endif

        !--------------------------------------------------

        ! [ B_{i} ] matrix
        DO nb = 1, nn

         jsize1 = ndof*(nb-1)+1
         jsize2 = ndof*(nb-1)+2
         jsize3 = ndof*(nb-1)+3
         jsize4 = ndof*(nb-1)+4
         jsize5 = ndof*(nb-1)+5
         jsize6 = ndof*(nb-1)+6

         B1(1, jsize1) =  shapederiv(nb, 1)
         B1(2, jsize1) =  0.0D0
         B1(3, jsize1) =  0.0D0
         B1(1, jsize2) =  0.0D0
         B1(2, jsize2) =  shapederiv(nb, 1)
         B1(3, jsize2) =  0.0D0
         B1(1, jsize3) =  0.0D0
         B1(2, jsize3) =  0.0D0
         B1(3, jsize3) =  shapederiv(nb, 1)
         B1(1, jsize4) =  0.0D0
         B1(2, jsize4) = -dudxi_rot(3, nb)
         B1(3, jsize4) =  dudxi_rot(2, nb)
         B1(1, jsize5) =  dudxi_rot(3, nb)
         B1(2, jsize5) =  0.0D0
         B1(3, jsize5) = -dudxi_rot(1, nb)
         B1(1, jsize6) = -dudxi_rot(2, nb)
         B1(2, jsize6) =  dudxi_rot(1, nb)
         B1(3, jsize6) =  0.0D0

         B2(1, jsize1) =  shapederiv(nb, 2)
         B2(2, jsize1) =  0.0D0
         B2(3, jsize1) =  0.0D0
         B2(1, jsize2) =  0.0D0
         B2(2, jsize2) =  shapederiv(nb, 2)
         B2(3, jsize2) =  0.0D0
         B2(1, jsize3) =  0.0D0
         B2(2, jsize3) =  0.0D0
         B2(3, jsize3) =  shapederiv(nb, 2)
         B2(1, jsize4) =  0.0D0
         B2(2, jsize4) = -dudeta_rot(3, nb)
         B2(3, jsize4) =  dudeta_rot(2, nb)
         B2(1, jsize5) =  dudeta_rot(3, nb)
         B2(2, jsize5) =  0.0D0
         B2(3, jsize5) = -dudeta_rot(1, nb)
         B2(1, jsize6) = -dudeta_rot(2, nb)
         B2(2, jsize6) =  dudeta_rot(1, nb)
         B2(3, jsize6) =  0.0D0

         B3(1, jsize1) =  0.0D0
         B3(2, jsize1) =  0.0D0
         B3(3, jsize1) =  0.0D0
         B3(1, jsize2) =  0.0D0
         B3(2, jsize2) =  0.0D0
         B3(3, jsize2) =  0.0D0
         B3(1, jsize3) =  0.0D0
         B3(2, jsize3) =  0.0D0
         B3(3, jsize3) =  0.0D0
         B3(1, jsize4) =  0.0D0
         B3(2, jsize4) = -dudzeta_rot(3, nb)
         B3(3, jsize4) =  dudzeta_rot(2, nb)
         B3(1, jsize5) =  dudzeta_rot(3, nb)
         B3(2, jsize5) =  0.0D0
         B3(3, jsize5) = -dudzeta_rot(1, nb)
         B3(1, jsize6) = -dudzeta_rot(2, nb)
         B3(2, jsize6) =  dudzeta_rot(1, nb)
         B3(3, jsize6) =  0.0D0

        END DO

        !--------------------------------------------------

        ! { C_{ij} } vector
        DO jsize = 1, ndof*nn

         Cv12(jsize) = ( cg1(1)*B1(2, jsize)   &
                        +cg2(1)*B2(2, jsize)   &
                        +cg3(1)*B3(2, jsize) ) &
                      -( cg1(2)*B1(1, jsize)   &
                        +cg2(2)*B2(1, jsize)   &
                        +cg3(2)*B3(1, jsize) )
         Cv13(jsize) = ( cg1(1)*B1(3, jsize)   &
                        +cg2(1)*B2(3, jsize)   &
                        +cg3(1)*B3(3, jsize) ) &
                      -( cg1(3)*B1(1, jsize)   &
                        +cg2(3)*B2(1, jsize)   &
                        +cg3(3)*B3(1, jsize) )
         Cv21(jsize) = ( cg1(2)*B1(1, jsize)   &
                        +cg2(2)*B2(1, jsize)   &
                        +cg3(2)*B3(1, jsize) ) &
                      -( cg1(1)*B1(2, jsize)   &
                        +cg2(1)*B2(2, jsize)   &
                        +cg3(1)*B3(2, jsize) )
         Cv23(jsize) = ( cg1(2)*B1(3, jsize)   &
                        +cg2(2)*B2(3, jsize)   &
                        +cg3(2)*B3(3, jsize) ) &
                      -( cg1(3)*B1(2, jsize)   &
                        +cg2(3)*B2(2, jsize)   &
                        +cg3(3)*B3(2, jsize) )
         Cv31(jsize) = ( cg1(3)*B1(1, jsize)   &
                        +cg2(3)*B2(1, jsize)   &
                        +cg3(3)*B3(1, jsize) ) &
                      -( cg1(1)*B1(3, jsize)   &
                        +cg2(1)*B2(3, jsize)   &
                        +cg3(1)*B3(3, jsize) )
         Cv32(jsize) = ( cg1(3)*B1(2, jsize)   &
                        +cg2(3)*B2(2, jsize)   &
                        +cg3(3)*B3(2, jsize) ) &
                      -( cg1(2)*B1(3, jsize)   &
                        +cg2(2)*B2(3, jsize)   &
                        +cg3(2)*B3(3, jsize) )

        END DO

        !--------------------------------------------------

        ! { Cw } vector
        DO nb = 1, nn

         DO j = 1, ndof

          jsize = ndof*(nb-1)+j

          Cv_w(jsize)                       &
          = v1_i(1)*Cv12(jsize)*v2_i(2) &
           +v1_i(1)*Cv13(jsize)*v2_i(3) &
           +v1_i(2)*Cv21(jsize)*v2_i(1) &
           +v1_i(2)*Cv23(jsize)*v2_i(3) &
           +v1_i(3)*Cv31(jsize)*v2_i(1) &
           +v1_i(3)*Cv32(jsize)*v2_i(2)

         END DO

        END DO

        ! { Ctheta } vector
        DO nb = 1, nn

         jsize1 = ndof*(nb-1)+1
         jsize2 = ndof*(nb-1)+2
         jsize3 = ndof*(nb-1)+3
         jsize4 = ndof*(nb-1)+4
         jsize5 = ndof*(nb-1)+5
         jsize6 = ndof*(nb-1)+6

         Cv_theta(jsize1) = 0.0D0
         Cv_theta(jsize2) = 0.0D0
         Cv_theta(jsize3) = 0.0D0
         Cv_theta(jsize4) = v3_i(1)*shapefunc(nb)
         Cv_theta(jsize5) = v3_i(2)*shapefunc(nb)
         Cv_theta(jsize6) = v3_i(3)*shapefunc(nb)

        END DO

        ! { C } vector
        DO jsize = 1, ndof*nn

         Cv(jsize) = Cv_theta(jsize)-0.5D0*Cv_w(jsize)

        END DO

        !--------------------------------------------------

        ! [ K L ] matrix
        shellmatl = int(gausses(1)%pMaterial%variables(M_SHELL_MATLTYPE))
        if (shellmatl == 0)then
          if ( n_total_layer .ge. 2 ) then
            DO jsize = 1, ndof*nn
              DO isize = 1, ndof*nn

          tmpstiff(isize, jsize)                &
          = tmpstiff(isize, jsize)              &
           +w_w_w_det*gausses(1)%pMaterial%variables(100+3*n_layer)/thick*alpha*Cv(isize)*Cv(jsize)

              END DO
            END DO
          else

            DO jsize = 1, ndof*nn
              DO isize = 1, ndof*nn

          tmpstiff(isize, jsize)                &
          = tmpstiff(isize, jsize)              &
           +w_w_w_det*alpha*Cv(isize)*Cv(jsize)

              END DO
            END DO

          endif
        elseif (shellmatl == 1)then
          DO jsize = 1, ndof*nn
          DO isize = 1, ndof*nn

        tmpstiff(isize, jsize)                &
        = tmpstiff(isize, jsize)              &
         +w_w_w_det*gausses(1)%pMaterial%variables(100+8*n_layer-5)/thick*alpha*Cv(isize)*Cv(jsize)

          END DO
          END DO
        else
          write(*,*)"ERROR : shellmatl isnot correct"; stop
        endif

        !--------------------------------------------------

       END DO

       !--------------------------------------------------------

      END DO

      !--------------------------------------------------------------

      stiff(1:nn*ndof, 1:nn*ndof) = tmpstiff(1:nn*ndof, 1:nn*ndof)

!--------------------------------------------------------------------
    END DO     !< LAMINATED SHELL ANALYSIS

!   write(*,"(24E25.16)") stiff

!******************** Shell-Solid mixed analysis ********************
!		mixglaf = 0 < natural shell (6 dof)
!		mixglaf = 1 < mixed 361 (3*2 dof)*8 nod
!		mixglaf = 2 < mixed 351 (3*2 dof)*6 nod
       if( mixflag == 1 )then

!       write(*,*) 'convert for shell-solid mixed analysis'
        sstable(1) = 1
        sstable(2) = 2
        sstable(3) = 3
        sstable(4) = 7
        sstable(5) = 8
        sstable(6) = 9
        sstable(7) = 13
        sstable(8) = 14
        sstable(9) = 15
        sstable(10)= 19
        sstable(11)= 20
        sstable(12)= 21
        sstable(13)= 4
        sstable(14)= 5
        sstable(15)= 6
        sstable(16)= 10
        sstable(17)= 11
        sstable(18)= 12
        sstable(19)= 16
        sstable(20)= 17
        sstable(21)= 18
        sstable(22)= 22
        sstable(23)= 23
        sstable(24)= 24

        tmpstiff(1:nn*ndof, 1:nn*ndof) = stiff(1:nn*ndof, 1:nn*ndof)

       DO i = 1, nn*ndof
       DO j = 1, nn*ndof
        stiff(i,j) = tmpstiff(sstable(i),sstable(j))
       ENDDO
       ENDDO

    elseif( mixflag == 2 )then
    !		write(*,*) 'convert for shell-solid mixed analysis 351'
        sstable(1) = 1
        sstable(2) = 2
        sstable(3) = 3
        sstable(4) = 7
        sstable(5) = 8
        sstable(6) = 9
        sstable(7) = 13
        sstable(8) = 14
        sstable(9) = 15
        sstable(10)= 4
        sstable(11)= 5
        sstable(12)= 6
        sstable(13)= 10
        sstable(14)= 11
        sstable(15)= 12
        sstable(16)= 16
        sstable(17)= 17
        sstable(18)= 18

        tmpstiff(1:nn*ndof, 1:nn*ndof) = stiff(1:nn*ndof, 1:nn*ndof)

       DO i = 1, nn*ndof
       DO j = 1, nn*ndof
        stiff(i,j) = tmpstiff(sstable(i),sstable(j))
       ENDDO
       ENDDO

     endif

      RETURN

!####################################################################
      END SUBROUTINE STF_Shell_MITC
!####################################################################


!####################################################################
      SUBROUTINE ElementStress_Shell_MITC                  &
                 (etype, nn, ndof, ecoord, gausses, edisp, &
                  strain, stress, thick, zeta, n_layer, n_total_layer)
!####################################################################

      USE mMechGauss
      USE gauss_integration
      USE m_MatMatrix

!--------------------------------------------------------------------

      INTEGER(KIND = kint), INTENT(IN) :: etype
      INTEGER(KIND = kint), INTENT(IN) :: nn
      INTEGER(KIND = kint), INTENT(IN) :: ndof
      REAL(KIND = kreal), INTENT(IN)   :: ecoord(3, nn)
      TYPE(tGaussStatus), INTENT(IN)   :: gausses(:)
      REAL(KIND = kreal), INTENT(IN)   :: edisp(6, nn)
      REAL(KIND = kreal), INTENT(OUT)  :: strain(:,:)
      REAL(KIND = kreal), INTENT(OUT)  :: stress(:,:)
      REAL(KIND = kreal), INTENT(IN)   :: thick
      REAL(KIND = kreal), INTENT(IN)   :: zeta

!--------------------------------------------------------------------

      INTEGER :: i, j, m
      INTEGER :: lx
      INTEGER :: fetype
      INTEGER :: ntying
      INTEGER :: npoints_tying(3)
      INTEGER :: it, ip
      INTEGER :: na, nb
      INTEGER :: isize, jsize
      INTEGER :: jsize1, jsize2, jsize3, &
                 jsize4, jsize5, jsize6
      INTEGER :: n_layer, n_total_layer, shellmatl

      REAL(KIND = kreal) :: D(5, 5)
      REAL(KIND = kreal) :: elem(3, nn)
      REAL(KIND = kreal) :: xi_lx, eta_lx, zeta_ly
      REAL(KIND = kreal) :: B_di(5, ndof*nn, 6, 3)
      REAL(KIND = kreal) :: B1(3, ndof*nn), B2(3, ndof*nn), &
                            B3(3, ndof*nn)
      REAL(KIND = kreal) :: naturalcoord(2)
      REAL(KIND = kreal) :: tpcoord(6, 2, 3)
      REAL(KIND = kreal) :: nncoord(nn, 2)
      REAL(KIND = kreal) :: shapefunc(nn)
      REAL(KIND = kreal) :: shapederiv(nn, 2)
      REAL(KIND = kreal) :: alpha
      REAL(KIND = kreal) :: xxi_lx, eeta_lx
      REAL(KIND = kreal) :: xxi_di(6, 3), eeta_di(6, 3)
      REAL(KIND = kreal) :: h(nn, 3)
      REAL(KIND = kreal) :: v1(3, nn), v2(3, nn), v3(3, nn)
      REAL(KIND = kreal) :: v1_abs, v2_abs, v3_abs
      REAL(KIND = kreal) :: a_over_2_v3(3, nn)
      REAL(KIND = kreal) :: a_over_2_theta_cross_v3(3, nn)
      REAL(KIND = kreal) :: u_rot(3, nn)
      REAL(KIND = kreal) :: theta(3, nn)
      REAL(KIND = kreal) :: dudxi(3), dudeta(3), dudzeta(3)
      REAL(KIND = kreal) :: dudxi_rot(3, nn), dudeta_rot(3, nn), &
                            dudzeta_rot(3, nn)
      REAL(KIND = kreal) :: g1(3), g2(3), g3(3)
      REAL(KIND = kreal) :: g1_abs, g2_abs, g3_abs
      REAL(KIND = kreal) :: e_0(3)
      REAL(KIND = kreal) :: cg1(3), cg2(3), cg3(3)
      REAL(KIND = kreal) :: det
      REAL(KIND = kreal) :: det_cg3(3)
      REAL(KIND = kreal) :: det_inv
      REAL(KIND = kreal) :: det_cg3_abs
      REAL(KIND = kreal) :: e1_hat(3), e2_hat(3), e3_hat(3)
      REAL(KIND = kreal) :: e1_hat_abs, e2_hat_abs
      REAL(KIND = kreal) :: e11, e22, e12_2, e23_2, e31_2
      REAL(KIND = kreal) :: e11_di(6, 3), e22_di(6, 3),     &
                            e12_di_2(6, 3), e23_di_2(6, 3), &
                            e31_di_2(6, 3)
      REAL(KIND = kreal) :: E(3, 3), Ev(5)
      REAL(KIND = kreal) :: S(3, 3), Sv(5)
      REAL(KIND = kreal) :: sigma_layer

!--------------------------------------------------------------------

      ! for lamina stress

      ! MITC4
      IF( etype .EQ. fe_mitc4_shell ) THEN

       fetype = fe_mitc4_shell

       ntying = 1
       npoints_tying(1)= 4

      ! MITC9
      ELSE IF( etype .EQ. fe_mitc9_shell ) THEN

       fetype = fe_mitc9_shell

       ntying = 3
       npoints_tying(1)= 6
       npoints_tying(2)= 6
       npoints_tying(3)= 4

      ! MITC3
      ELSE IF( etype .EQ. fe_mitc3_shell ) THEN

       fetype = fe_mitc3_shell

       ntying = 1
       npoints_tying(1)= 3

      END IF

!--------------------------------------------------------------------

      elem(:, :) = ecoord(:, :)

!--------------------------------------------------------------------

      DO na = 1, nn

       theta(1, na) = edisp(4, na)
       theta(2, na) = edisp(5, na)
       theta(3, na) = edisp(6, na)

      END DO

!-------------------------------------------------------------------

      ! xi-coordinate at a node in a local element
      ! eta-coordinate at a node in a local element
      CALL getNodalNaturalCoord(fetype, nncoord)

!-------------------------------------------------------------------

      ! MITC4
      IF( etype .EQ. fe_mitc4_shell ) THEN

       !--------------------------------------------------------

       ! xi-coordinate at a tying pont in a local element
       tpcoord(1, 1, 1) =  0.0D0
       tpcoord(2, 1, 1) =  1.0D0
       tpcoord(3, 1, 1) =  0.0D0
       tpcoord(4, 1, 1) = -1.0D0
       ! eta-coordinate at a tying point in a local element
       tpcoord(1, 2, 1) = -1.0D0
       tpcoord(2, 2, 1) =  0.0D0
       tpcoord(3, 2, 1) =  1.0D0
       tpcoord(4, 2, 1) =  0.0D0

       !--------------------------------------------------------

      ! MITC9
      ELSE IF( etype .EQ. fe_mitc9_shell ) THEN

       !--------------------------------------------------------

       ! xi-coordinate at a tying point in a local element
       tpcoord(1, 1, 1) = -1.0D0*DSQRT( 1.0D0/3.0D0 )
       tpcoord(2, 1, 1) =  1.0D0*DSQRT( 1.0D0/3.0D0 )
       tpcoord(3, 1, 1) =  1.0D0*DSQRT( 1.0D0/3.0D0 )
       tpcoord(4, 1, 1) = -1.0D0*DSQRT( 1.0D0/3.0D0 )
       tpcoord(5, 1, 1) =  1.0D0*DSQRT( 1.0D0/3.0D0 )
       tpcoord(6, 1, 1) = -1.0D0*DSQRT( 1.0D0/3.0D0 )
       ! eta-coordinate at a tying point in a local element
       tpcoord(1, 2, 1) = -1.0D0*DSQRT( 3.0D0/5.0D0 )
       tpcoord(2, 2, 1) = -1.0D0*DSQRT( 3.0D0/5.0D0 )
       tpcoord(3, 2, 1) =  1.0D0*DSQRT( 3.0D0/5.0D0 )
       tpcoord(4, 2, 1) =  1.0D0*DSQRT( 3.0D0/5.0D0 )
       tpcoord(5, 2, 1) =  0.0D0*DSQRT( 3.0D0/5.0D0 )
       tpcoord(6, 2, 1) =  0.0D0*DSQRT( 3.0D0/5.0D0 )

       ! xi-coordinate at a tying point in a local element
       tpcoord(1, 1, 2) = -1.0D0*DSQRT( 3.0D0/5.0D0 )
       tpcoord(2, 1, 2) =  0.0D0*DSQRT( 3.0D0/5.0D0 )
       tpcoord(3, 1, 2) =  1.0D0*DSQRT( 3.0D0/5.0D0 )
       tpcoord(4, 1, 2) =  1.0D0*DSQRT( 3.0D0/5.0D0 )
       tpcoord(5, 1, 2) =  0.0D0*DSQRT( 3.0D0/5.0D0 )
       tpcoord(6, 1, 2) = -1.0D0*DSQRT( 3.0D0/5.0D0 )
       ! eta-coordinate at a tying point in a local element
       tpcoord(1, 2, 2) = -1.0D0*DSQRT( 1.0D0/3.0D0 )
       tpcoord(2, 2, 2) = -1.0D0*DSQRT( 1.0D0/3.0D0 )
       tpcoord(3, 2, 2) = -1.0D0*DSQRT( 1.0D0/3.0D0 )
       tpcoord(4, 2, 2) =  1.0D0*DSQRT( 1.0D0/3.0D0 )
       tpcoord(5, 2, 2) =  1.0D0*DSQRT( 1.0D0/3.0D0 )
       tpcoord(6, 2, 2) =  1.0D0*DSQRT( 1.0D0/3.0D0 )

       ! xi-coordinate at a tying point in a local element
       tpcoord(1, 1, 3) = -1.0D0*DSQRT( 1.0D0/3.0D0 )
       tpcoord(2, 1, 3) =  1.0D0*DSQRT( 1.0D0/3.0D0 )
       tpcoord(3, 1, 3) =  1.0D0*DSQRT( 1.0D0/3.0D0 )
       tpcoord(4, 1, 3) = -1.0D0*DSQRT( 1.0D0/3.0D0 )
       ! eta-coordinate at a tying point in a local element
       tpcoord(1, 2, 3) = -1.0D0*DSQRT( 1.0D0/3.0D0 )
       tpcoord(2, 2, 3) = -1.0D0*DSQRT( 1.0D0/3.0D0 )
       tpcoord(3, 2, 3) =  1.0D0*DSQRT( 1.0D0/3.0D0 )
       tpcoord(4, 2, 3) =  1.0D0*DSQRT( 1.0D0/3.0D0 )

       !--------------------------------------------------------

       ! Xi-coordinate at a tying point in a local element
       xxi_di(1, 1)  = -1.0D0
       xxi_di(2, 1)  =  1.0D0
       xxi_di(3, 1)  =  1.0D0
       xxi_di(4, 1)  = -1.0D0
       xxi_di(5, 1)  =  1.0D0
       xxi_di(6, 1)  = -1.0D0
       ! Eta-coordinate at a tying point in a local element
       eeta_di(1, 1) = -1.0D0
       eeta_di(2, 1) = -1.0D0
       eeta_di(3, 1) =  1.0D0
       eeta_di(4, 1) =  1.0D0
       eeta_di(5, 1) =  0.0D0
       eeta_di(6, 1) =  0.0D0

       ! Xi-coordinate at a tying point in a local element
       xxi_di(1, 2)  = -1.0D0
       xxi_di(2, 2)  =  0.0D0
       xxi_di(3, 2)  =  1.0D0
       xxi_di(4, 2)  =  1.0D0
       xxi_di(5, 2)  =  0.0D0
       xxi_di(6, 2)  = -1.0D0
       ! Eta-coordinate at a tying point in a local element
       eeta_di(1, 2) = -1.0D0
       eeta_di(2, 2) = -1.0D0
       eeta_di(3, 2) = -1.0D0
       eeta_di(4, 2) =  1.0D0
       eeta_di(5, 2) =  1.0D0
       eeta_di(6, 2) =  1.0D0

       !--------------------------------------------------------

      ! MITC3
      ELSE IF( etype .EQ. fe_mitc3_shell ) THEN

       !--------------------------------------------------------

       ! xi-coordinate at a tying point in a local element
       tpcoord(1, 1, 1)  = 0.5D0
       tpcoord(2, 1, 1)  = 0.0D0
       tpcoord(3, 1, 1)  = 0.5D0
       ! eta-coordinate at a tying point in a local element
       tpcoord(1, 2, 1) = 0.0D0
       tpcoord(2, 2, 1) = 0.5D0
       tpcoord(3, 2, 1) = 0.5D0

       !--------------------------------------------------------

      END IF

!--------------------------------------------------------------------

      ! xi-coordinate at the center point in a local element
      ! eta-coordinate at the center point in a local element
      naturalcoord(1) = 0.0D0
      naturalcoord(2) = 0.0D0

      CALL getShapeDeriv(fetype, naturalcoord, shapederiv)

      !--------------------------------------------------------------

      ! Covariant basis vector
      DO i = 1, 3

       g1(i) = 0.0D0

       DO na = 1, nn

        g1(i) = g1(i)+shapederiv(na, 1) &
                      *elem(i, na)

       END DO

      END DO

      e_0(1) = g1(1)
      e_0(2) = g1(2)
      e_0(3) = g1(3)

      !--------------------------------------------------------------

      DO nb = 1, nn

       !--------------------------------------------------------

       naturalcoord(1) = nncoord(nb, 1)
       naturalcoord(2) = nncoord(nb, 2)

       CALL getShapeDeriv(fetype, naturalcoord, shapederiv)

       !--------------------------------------------------------

       ! Covariant basis vector
       DO i = 1, 3

        g1(i) = 0.0D0
        g2(i) = 0.0D0

        DO na = 1, nn

         g1(i) = g1(i)+shapederiv(na, 1) &
                       *elem(i, na)
         g2(i) = g2(i)+shapederiv(na, 2) &
                       *elem(i, na)

        END DO

       END DO

       !--------------------------------------------------------

       det_cg3(1) = g1(2)*g2(3)-g1(3)*g2(2)
       det_cg3(2) = g1(3)*g2(1)-g1(1)*g2(3)
       det_cg3(3) = g1(1)*g2(2)-g1(2)*g2(1)

       det_cg3_abs = DSQRT( det_cg3(1)*det_cg3(1)   &
                           +det_cg3(2)*det_cg3(2)   &
                           +det_cg3(3)*det_cg3(3) )

       v3(1, nb) = det_cg3(1)/det_cg3_abs
       v3(2, nb) = det_cg3(2)/det_cg3_abs
       v3(3, nb) = det_cg3(3)/det_cg3_abs

       !--------------------------------------------------------

       v2(1, nb) = v3(2, nb)*e_0(3)-v3(3, nb)*e_0(2)
       v2(2, nb) = v3(3, nb)*e_0(1)-v3(1, nb)*e_0(3)
       v2(3, nb) = v3(1, nb)*e_0(2)-v3(2, nb)*e_0(1)

       v2_abs = DSQRT( v2(1, nb)*v2(1, nb)   &
                      +v2(2, nb)*v2(2, nb)   &
                      +v2(3, nb)*v2(3, nb) )

       IF( v2_abs .GT. 1.0D-15 ) THEN

        v2(1, nb) = v2(1, nb)/v2_abs
        v2(2, nb) = v2(2, nb)/v2_abs
        v2(3, nb) = v2(3, nb)/v2_abs

        v1(1, nb) = v2(2, nb)*v3(3, nb) &
                   -v2(3, nb)*v3(2, nb)
        v1(2, nb) = v2(3, nb)*v3(1, nb) &
                   -v2(1, nb)*v3(3, nb)
        v1(3, nb) = v2(1, nb)*v3(2, nb) &
                   -v2(2, nb)*v3(1, nb)

        v1_abs = DSQRT( v1(1, nb)*v1(1, nb)   &
                       +v1(2, nb)*v1(2, nb)   &
                       +v1(3, nb)*v1(3, nb) )

        v1(1, nb) = v1(1, nb)/v1_abs
        v1(2, nb) = v1(2, nb)/v1_abs
        v1(3, nb) = v1(3, nb)/v1_abs

       ELSE

        v1(1, nb) =  0.0D0
        v1(2, nb) =  0.0D0
        v1(3, nb) = -1.0D0

        v2(1, nb) = 0.0D0
        v2(2, nb) = 1.0D0
        v2(3, nb) = 0.0D0

       END IF

       !--------------------------------------------------------

       v3(1, nb) = v1(2, nb)*v2(3, nb) &
                  -v1(3, nb)*v2(2, nb)
       v3(2, nb) = v1(3, nb)*v2(1, nb) &
                  -v1(1, nb)*v2(3, nb)
       v3(3, nb) = v1(1, nb)*v2(2, nb) &
                  -v1(2, nb)*v2(1, nb)

       v3_abs = DSQRT( v3(1, nb)*v3(1, nb)   &
                      +v3(2, nb)*v3(2, nb)   &
                      +v3(3, nb)*v3(3, nb) )

       v3(1, nb) = v3(1, nb)/v3_abs
       v3(2, nb) = v3(2, nb)/v3_abs
       v3(3, nb) = v3(3, nb)/v3_abs

       !--------------------------------------------------------

       a_over_2_v3(1, nb) = 0.5D0*thick*v3(1, nb)
       a_over_2_v3(2, nb) = 0.5D0*thick*v3(2, nb)
       a_over_2_v3(3, nb) = 0.5D0*thick*v3(3, nb)

       !--------------------------------------------------------

       a_over_2_theta_cross_v3(1, nb)    &
       = theta(2, nb)*a_over_2_v3(3, nb) &
        -theta(3, nb)*a_over_2_v3(2, nb)
       a_over_2_theta_cross_v3(2, nb)    &
       = theta(3, nb)*a_over_2_v3(1, nb) &
        -theta(1, nb)*a_over_2_v3(3, nb)
       a_over_2_theta_cross_v3(3, nb)    &
       = theta(1, nb)*a_over_2_v3(2, nb) &
        -theta(2, nb)*a_over_2_v3(1, nb)

       !--------------------------------------------------------

      END DO

!--------------------------------------------------------------------
!				Modified stress in laminated shell
!--------------------------------------------------------------------
      ! MITC4
      IF( etype .EQ. fe_mitc4_shell ) THEN

       zeta_ly = 0.0D0

      ! MITC9
      ELSE IF( etype .EQ. fe_mitc9_shell ) THEN

       zeta_ly = zeta

      ! MITC3
      ELSE IF( etype .EQ. fe_mitc3_shell ) THEN

       zeta_ly = 0.0D0

      END IF

      !---------------------------------------------------------

      DO it = 1, ntying

       DO ip = 1, npoints_tying(it)

        !-------------------------------------------------

        naturalcoord(1) = tpcoord(ip, 1, it)
        naturalcoord(2) = tpcoord(ip, 2, it)

        CALL getShapeFunc(fetype, naturalcoord, shapefunc)

        CALL getShapeDeriv(fetype, naturalcoord, shapederiv)

        !-------------------------------------------------

        DO na = 1, nn

         DO i = 1, 3

          u_rot(i, na)                      &
          = shapefunc(na)                   &
            *( zeta_ly*a_over_2_v3(i, na) )

          dudxi_rot(i, na)                  &
          = shapederiv(na, 1)               &
            *( zeta_ly*a_over_2_v3(i, na) )
          dudeta_rot(i, na)                 &
          = shapederiv(na, 2)               &
            *( zeta_ly*a_over_2_v3(i, na) )
          dudzeta_rot(i, na)                &
          = shapefunc(na)                   &
            *( a_over_2_v3(i, na) )

         END DO

        END DO

        !-------------------------------------------------

        ! Covariant basis vector
        DO i = 1, 3

         g1(i) = 0.0D0
         g2(i) = 0.0D0
         g3(i) = 0.0D0

         DO na = 1, nn

          g1(i) = g1(i)+shapederiv(na, 1)  &
                        *elem(i, na)       &
                      +dudxi_rot(i, na)
          g2(i) = g2(i)+shapederiv(na, 2)  &
                        *elem(i, na)       &
                       +dudeta_rot(i, na)
          g3(i) = g3(i)+dudzeta_rot(i, na)

         END DO

        END DO

        !---------------------------------------------

        DO i = 1, 3

         dudxi(i)   = 0.0D0
         dudeta(i)  = 0.0D0
         dudzeta(i) = 0.0D0

         DO na = 1, nn

          dudxi(i)                                      &
          = dudxi(i)                                    &
           +shapederiv(na, 1)                           &
            *( edisp(i, na)                             &
              +zeta_ly*a_over_2_theta_cross_v3(i, na) )
          dudeta(i)                                     &
          = dudeta(i)                                   &
           +shapederiv(na, 2)                           &
            *( edisp(i, na)                             &
              +zeta_ly*a_over_2_theta_cross_v3(i, na) )
          dudzeta(i)                                    &
          = dudzeta(i)                                  &
           +shapefunc(na)                               &
            *( a_over_2_theta_cross_v3(i, na) )

         END DO

        END DO

        !---------------------------------------------

        ! Infinitesimal strain tensor
        e11_di(ip, it)                             &
        = 0.5D0                                    &
          *( ( g1(1)*dudxi(1) +dudxi(1) *g1(1) )   &
            +( g1(2)*dudxi(2) +dudxi(2) *g1(2) )   &
            +( g1(3)*dudxi(3) +dudxi(3) *g1(3) ) )
        e22_di(ip, it)                             &
        = 0.5D0                                    &
          *( ( g2(1)*dudeta(1)+dudeta(1)*g2(1) )   &
            +( g2(2)*dudeta(2)+dudeta(2)*g2(2) )   &
            +( g2(3)*dudeta(3)+dudeta(3)*g2(3) ) )
        e12_2                                   &
        = ( g1(1)*dudeta(1) +dudxi(1)  *g2(1) ) &
         +( g1(2)*dudeta(2) +dudxi(2)  *g2(2) ) &
         +( g1(3)*dudeta(3) +dudxi(3)  *g2(3) )
        e23_di_2(ip, it)                        &
        = ( g2(1)*dudzeta(1)+dudeta(1) *g3(1) ) &
         +( g2(2)*dudzeta(2)+dudeta(2) *g3(2) ) &
         +( g2(3)*dudzeta(3)+dudeta(3) *g3(3) )
        e31_di_2(ip, it)                        &
        = ( g3(1)*dudxi(1)  +dudzeta(1)*g1(1) ) &
         +( g3(2)*dudxi(2)  +dudzeta(2)*g1(2) ) &
         +( g3(3)*dudxi(3)  +dudzeta(3)*g1(3) )

        !-------------------------------------------------

       END DO

      END DO

      !--------------------------------------------------------

      IF ( n_total_layer.ge. 2 ) THEN
        shellmatl = int(gausses(1)%pMaterial%variables(M_SHELL_MATLTYPE))
        if (shellmatl == 0)then
          sigma_layer = 0.0D0
          DO m = 1, n_layer
            sigma_layer = sigma_layer + 2 * gausses(1)%pMaterial%variables(100+3*m) / thick
          END DO
          zeta_ly = -1 + sigma_layer - gausses(1)%pMaterial%variables(100+3*n_layer) / thick * (1-zeta)
        elseif (shellmatl == 1)then
          sigma_layer = 0.0D0
          DO m = 1, n_layer
            sigma_layer = sigma_layer + 2 * gausses(1)%pMaterial%variables(100+8*m-5) / thick
          END DO
          zeta_ly = -1 + sigma_layer - gausses(1)%pMaterial%variables(100+8*n_layer-5) / thick * (1-zeta)
        else
        write(*,*)"ERROR : shellmatl isnot correct"; stop
        endif
      ELSE
        zeta_ly = zeta
      ENDIF

      !--------------------------------------------------------

      DO lx = 1, nn

       !--------------------------------------------------

       naturalcoord(1) = nncoord(lx, 1)
       naturalcoord(2) = nncoord(lx, 2)

       xi_lx  = naturalcoord(1)
       eta_lx = naturalcoord(2)

       CALL getShapeFunc(fetype, naturalcoord, shapefunc)

       CALL getShapeDeriv(fetype, naturalcoord, shapederiv)

       !--------------------------------------------------

       DO na = 1, nn

        DO i = 1, 3

         u_rot(i, na)                      &
         = shapefunc(na)                   &
           *( zeta_ly*a_over_2_v3(i, na) )

         dudxi_rot(i, na)                  &
         = shapederiv(na, 1)               &
           *( zeta_ly*a_over_2_v3(i, na) )
         dudeta_rot(i, na)                 &
        = shapederiv(na, 2)                &
           *( zeta_ly*a_over_2_v3(i, na) )
         dudzeta_rot(i, na)                &
         = shapefunc(na)                   &
           *( a_over_2_v3(i, na) )

        END DO

       END DO

       !--------------------------------------------------

       ! Covariant basis vector
       DO i = 1, 3

        g1(i) = 0.0D0
        g2(i) = 0.0D0
        g3(i) = 0.0D0

        DO na = 1, nn

         g1(i) = g1(i)+shapederiv(na, 1)  &
                       *elem(i, na)       &
                      +dudxi_rot(i, na)
         g2(i) = g2(i)+shapederiv(na, 2)  &
                       *elem(i, na)       &
                      +dudeta_rot(i, na)
         g3(i) = g3(i)+dudzeta_rot(i, na)

        END DO
       END DO

       !--------------------------------------------------

       ! Jacobian
       det = g1(1)*( g2(2)*g3(3)-g2(3)*g3(2) ) &
            +g1(2)*( g2(3)*g3(1)-g2(1)*g3(3) ) &
            +g1(3)*( g2(1)*g3(2)-g2(2)*g3(1) )

      if(det == 0.0d0) then
        write(*,*)"ERROR:LIB Shell in l2009 Not Jacobian"
        stop
      endif

       det_inv = 1.0D0/det

       !--------------------------------------------------

       ! Contravariant basis vector
       cg1(1) = det_inv                      &
                *( g2(2)*g3(3)-g2(3)*g3(2) )
       cg1(2) = det_inv                      &
                *( g2(3)*g3(1)-g2(1)*g3(3) )
       cg1(3) = det_inv                      &
                *( g2(1)*g3(2)-g2(2)*g3(1) )
       cg2(1) = det_inv                      &
                *( g3(2)*g1(3)-g3(3)*g1(2) )
       cg2(2) = det_inv                      &
                *( g3(3)*g1(1)-g3(1)*g1(3) )
       cg2(3) = det_inv                      &
                *( g3(1)*g1(2)-g3(2)*g1(1) )
       cg3(1) = det_inv                      &
                *( g1(2)*g2(3)-g1(3)*g2(2) )
       cg3(2) = det_inv                      &
                *( g1(3)*g2(1)-g1(1)*g2(3) )
       cg3(3) = det_inv                      &
                *( g1(1)*g2(2)-g1(2)*g2(1) )

       !--------------------------------------------------

       g3_abs = DSQRT( g3(1)*g3(1)   &
                      +g3(2)*g3(2)   &
                      +g3(3)*g3(3) )

       !--------------------------------------------------

       ! Orthonormal vectors

       e3_hat(1) = g3(1)/g3_abs
       e3_hat(2) = g3(2)/g3_abs
       e3_hat(3) = g3(3)/g3_abs

       e1_hat(1) = g2(2)*e3_hat(3) &
                  -g2(3)*e3_hat(2)
       e1_hat(2) = g2(3)*e3_hat(1) &
                  -g2(1)*e3_hat(3)
       e1_hat(3) = g2(1)*e3_hat(2) &
                  -g2(2)*e3_hat(1)
       e1_hat_abs = DSQRT( e1_hat(1)*e1_hat(1)   &
                          +e1_hat(2)*e1_hat(2)   &
                          +e1_hat(3)*e1_hat(3) )
       e1_hat(1) = e1_hat(1)/e1_hat_abs
       e1_hat(2) = e1_hat(2)/e1_hat_abs
       e1_hat(3) = e1_hat(3)/e1_hat_abs

       e2_hat(1) = e3_hat(2)*e1_hat(3) &
                  -e3_hat(3)*e1_hat(2)
       e2_hat(2) = e3_hat(3)*e1_hat(1) &
                  -e3_hat(1)*e1_hat(3)
       e2_hat(3) = e3_hat(1)*e1_hat(2) &
                  -e3_hat(2)*e1_hat(1)
       e2_hat_abs = DSQRT( e2_hat(1)*e2_hat(1)   &
                          +e2_hat(2)*e2_hat(2)   &
                          +e2_hat(3)*e2_hat(3) )
       e2_hat(1) = e2_hat(1)/e2_hat_abs
       e2_hat(2) = e2_hat(2)/e2_hat_abs
       e2_hat(3) = e2_hat(3)/e2_hat_abs

       !--------------------------------------------------

       DO i = 1, 3

        dudxi(i)   = 0.0D0
        dudeta(i)  = 0.0D0
        dudzeta(i) = 0.0D0

        DO na = 1, nn

         dudxi(i)                                      &
         = dudxi(i)                                    &
          +shapederiv(na, 1)                           &
           *( edisp(i, na)                             &
             +zeta_ly*a_over_2_theta_cross_v3(i, na) )
         dudeta(i)                                     &
         = dudeta(i)                                   &
          +shapederiv(na, 2)                           &
           *( edisp(i, na)                             &
             +zeta_ly*a_over_2_theta_cross_v3(i, na) )
         dudzeta(i)                                    &
         = dudzeta(i)                                  &
          +shapefunc(na)                               &
           *( a_over_2_theta_cross_v3(i, na) )

        END DO

       END DO

       !--------------------------------------------------

       ! Infinitesimal strain tensor
       e11                                        &
       = 0.5D0                                    &
         *( ( g1(1)*dudxi(1) +dudxi(1) *g1(1) )   &
           +( g1(2)*dudxi(2) +dudxi(2) *g1(2) )   &
           +( g1(3)*dudxi(3) +dudxi(3) *g1(3) ) )
       e22                                        &
       = 0.5D0                                    &
         *( ( g2(1)*dudeta(1)+dudeta(1)*g2(1) )   &
           +( g2(2)*dudeta(2)+dudeta(2)*g2(2) )   &
           +( g2(3)*dudeta(3)+dudeta(3)*g2(3) ) )
       e12_2                                   &
       = ( g1(1)*dudeta(1) +dudxi(1)  *g2(1) ) &
        +( g1(2)*dudeta(2) +dudxi(2)  *g2(2) ) &
        +( g1(3)*dudeta(3) +dudxi(3)  *g2(3) )
       e23_2                                   &
       = ( g2(1)*dudzeta(1)+dudeta(1) *g3(1) ) &
        +( g2(2)*dudzeta(2)+dudeta(2) *g3(2) ) &
        +( g2(3)*dudzeta(3)+dudeta(3) *g3(3) )
       e31_2                                   &
       = ( g3(1)*dudxi(1)  +dudzeta(1)*g1(1) ) &
        +( g3(2)*dudxi(2)  +dudzeta(2)*g1(2) ) &
        +( g3(3)*dudxi(3)  +dudzeta(3)*g1(3) )

       ! MITC4
       IF( etype .EQ. fe_mitc4_shell ) THEN

        e23_2 = 0.0D0
        e31_2 = 0.0D0

        ! e23_as_2
        e23_2                                   &
        = 0.5D0*( 1.0D0-xi_lx  )*e23_di_2(4, 1) &
         +0.5D0*( 1.0D0+xi_lx  )*e23_di_2(2, 1)
        ! e31_as_2
        e31_2                                   &
        = 0.5D0*( 1.0D0-eta_lx )*e31_di_2(1, 1) &
         +0.5D0*( 1.0D0+eta_lx )*e31_di_2(3, 1)

       ! MITC9
       ELSE IF( etype .EQ. fe_mitc9_shell ) THEN

        xxi_lx  = xi_lx /DSQRT( 1.0D0/3.0D0 )
        eeta_lx = eta_lx/DSQRT( 3.0D0/5.0D0 )

        DO ip = 1, npoints_tying(1)

         h(ip, 1)                                     &
         = ( 0.5D0*( 1.0D0+xxi_di(ip, 1)*xxi_lx ) )   &
           *( ( 0.5D0*eeta_di(ip, 1)*eeta_lx )        &
              *( 1.0D0+eeta_di(ip, 1)*eeta_lx )       &
             +( 1.0D0-eeta_di(ip, 1)*eeta_di(ip, 1) ) &
              *( 1.0D0-eeta_lx*eeta_lx ) )

        END DO

        xxi_lx  = xi_lx /DSQRT( 3.0D0/5.0D0 )
        eeta_lx = eta_lx/DSQRT( 1.0D0/3.0D0 )

        DO ip = 1, npoints_tying(2)

         h(ip, 2)                                      &
         = ( ( 0.5D0*xxi_di(ip, 2) *xxi_lx  )          &
             *( 1.0D0+xxi_di(ip, 2) *xxi_lx  )         &
            +( 1.0D0-xxi_di(ip, 2) *xxi_di(ip, 2)  )   &
             *( 1.0D0-xxi_lx*xxi_lx ) )                &
           *( 0.5D0*( 1.0D0+eeta_di(ip, 2)*eeta_lx ) )

        END DO

        xxi_lx  = xi_lx /DSQRT( 1.0D0/3.0D0 )
        eeta_lx = eta_lx/DSQRT( 1.0D0/3.0D0 )

        DO ip = 1, npoints_tying(3)

         h(ip, 3)                                      &
         = ( 0.5D0*( 1.0D0+xxi_di(ip, 1)*xxi_lx ) )    &
           *( 0.5D0*( 1.0D0+eeta_di(ip, 1)*eeta_lx ) )

        END DO

        e11   = 0.0D0
        e31_2 = 0.0D0

        ! e11_as, e31_as_2
        DO ip = 1, npoints_tying(1)

         e11   = e11  +h(ip, 1)*e11_di(ip, 1)
         e31_2 = e31_2+h(ip, 1)*e31_di_2(ip, 1)

        END DO

        e22   = 0.0D0
        e23_2 = 0.0D0

        ! e22_as, e23_as_2
        DO ip = 1, npoints_tying(2)

         e22   = e22  +h(ip, 2)*e22_di(ip, 2)
         e23_2 = e23_2+h(ip, 2)*e23_di_2(ip, 2)

        END DO

        e12_2 = 0.0D0

        ! e12_as_2
        DO ip = 1, npoints_tying(3)

         e12_2 = e12_2+h(ip, 3)*e12_di_2(ip, 3)

        END DO

       ! MITC3
       ELSE IF( etype .EQ. fe_mitc3_shell ) THEN

        e23_2 = 0.0D0
        e31_2 = 0.0D0

        ! e23_as_2
        e23_2                                       &
         = ( 1.0D0-xi_lx )*e23_di_2(2, 1)           &
          +xi_lx *e31_di_2(1, 1)                    &
          +xi_lx *( e23_di_2(3, 1)-e31_di_2(3, 1) )
        ! e31_as_2
        e31_2                                       &
         = eta_lx*e23_di_2(2, 1)                    &
          +( 1.0D0-eta_lx )*e31_di_2(1, 1)          &
          -eta_lx*( e23_di_2(3, 1)-e31_di_2(3, 1) )

       END IF

       !--------------------------------------------------

       ! { E } vector
       Ev(1) = e11
       Ev(2) = e22
       Ev(3) = e12_2
       Ev(4) = e23_2
       Ev(5) = e31_2

       ! Infinitesimal strain tensor
       ! [ E ] matrix
       E(1, 1) = Ev(1)
       E(2, 2) = Ev(2)
       E(3, 3) = 0.0D0
       E(1, 2) = 0.5D0*Ev(3)
       E(2, 1) = 0.5D0*Ev(3)
       E(2, 3) = 0.5D0*Ev(4)
       E(3, 2) = 0.5D0*Ev(4)
       E(3, 1) = 0.5D0*Ev(5)
       E(1, 3) = 0.5D0*Ev(5)

       !--------------------------------------------------
    !   write(*,*) 'Stress_n_layer', n_layer
       CALL MatlMatrix_Shell                        &
            (gausses(lx), Shell, D,                 &
             e1_hat, e2_hat, e3_hat, cg1, cg2, cg3, &
             alpha, n_layer)

       !--------------------------------------------------

       Sv = MATMUL( D, Ev )

       ! Infinitesimal stress tensor
       ! [ S ] matrix
       S(1, 1) = Sv(1)
       S(2, 2) = Sv(2)
       S(3, 3) = 0.0D0
       S(1, 2) = Sv(3)
       S(2, 1) = Sv(3)
       S(2, 3) = Sv(4)
       S(3, 2) = Sv(4)
       S(3, 1) = Sv(5)
       S(1, 3) = Sv(5)

       !------------------------------------------------

       stress(lx,1)                  &
        = ( S(1, 1)*g1(1)*g1(1)   &
           +S(1, 2)*g1(1)*g2(1)   &
           +S(1, 3)*g1(1)*g3(1)   &
           +S(2, 1)*g2(1)*g1(1)   &
           +S(2, 2)*g2(1)*g2(1)   &
           +S(2, 3)*g2(1)*g3(1)   &
           +S(3, 1)*g3(1)*g1(1)   &
           +S(3, 2)*g3(1)*g2(1) )
        stress(lx,2)                 &
        = ( S(1, 1)*g1(2)*g1(2)   &
           +S(1, 2)*g1(2)*g2(2)   &
           +S(1, 3)*g1(2)*g3(2)   &
           +S(2, 1)*g2(2)*g1(2)   &
           +S(2, 2)*g2(2)*g2(2)   &
           +S(2, 3)*g2(2)*g3(2)   &
           +S(3, 1)*g3(2)*g1(2)   &
           +S(3, 2)*g3(2)*g2(2) )
        stress(lx,3)                 &
        = ( S(1, 1)*g1(3)*g1(3)   &
           +S(1, 2)*g1(3)*g2(3)   &
           +S(1, 3)*g1(3)*g3(3)   &
           +S(2, 1)*g2(3)*g1(3)   &
           +S(2, 2)*g2(3)*g2(3)   &
           +S(2, 3)*g2(3)*g3(3)   &
           +S(3, 1)*g3(3)*g1(3)   &
           +S(3, 2)*g3(3)*g2(3) )
        stress(lx,4)                 &
        = ( S(1, 1)*g1(1)*g1(2)   &
           +S(1, 2)*g1(1)*g2(2)   &
           +S(1, 3)*g1(1)*g3(2)   &
           +S(2, 1)*g2(1)*g1(2)   &
           +S(2, 2)*g2(1)*g2(2)   &
           +S(2, 3)*g2(1)*g3(2)   &
           +S(3, 1)*g3(1)*g1(2)   &
           +S(3, 2)*g3(1)*g2(2) )
        stress(lx,5)                 &
        = ( S(1, 1)*g1(2)*g1(3)   &
           +S(1, 2)*g1(2)*g2(3)   &
           +S(1, 3)*g1(2)*g3(3)   &
           +S(2, 1)*g2(2)*g1(3)   &
           +S(2, 2)*g2(2)*g2(3)   &
           +S(2, 3)*g2(2)*g3(3)   &
           +S(3, 1)*g3(2)*g1(3)   &
           +S(3, 2)*g3(2)*g2(3) )
        stress(lx,6)                 &
        = ( S(1, 1)*g1(3)*g1(1)   &
           +S(1, 2)*g1(3)*g2(1)   &
           +S(1, 3)*g1(3)*g3(1)   &
           +S(2, 1)*g2(3)*g1(1)   &
           +S(2, 2)*g2(3)*g2(1)   &
           +S(2, 3)*g2(3)*g3(1)   &
           +S(3, 1)*g3(3)*g1(1)   &
           +S(3, 2)*g3(3)*g2(1) )

        strain(lx,1)                   &
        = ( E(1, 1)*cg1(1)*cg1(1)   &
           +E(1, 2)*cg1(1)*cg2(1)   &
           +E(1, 3)*cg1(1)*cg3(1)   &
           +E(2, 1)*cg2(1)*cg1(1)   &
           +E(2, 2)*cg2(1)*cg2(1)   &
           +E(2, 3)*cg2(1)*cg3(1)   &
           +E(3, 1)*cg3(1)*cg1(1)   &
           +E(3, 2)*cg3(1)*cg2(1) )
        strain(lx,2)                   &
        = ( E(1, 1)*cg1(2)*cg1(2)   &
           +E(1, 2)*cg1(2)*cg2(2)   &
           +E(1, 3)*cg1(2)*cg3(2)   &
           +E(2, 1)*cg2(2)*cg1(2)   &
           +E(2, 2)*cg2(2)*cg2(2)   &
           +E(2, 3)*cg2(2)*cg3(2)   &
           +E(3, 1)*cg3(2)*cg1(2)   &
           +E(3, 2)*cg3(2)*cg2(2) )
        strain(lx,3)                   &
        = ( E(1, 1)*cg1(3)*cg1(3)   &
           +E(1, 2)*cg1(3)*cg2(3)   &
           +E(1, 3)*cg1(3)*cg3(3)   &
           +E(2, 1)*cg2(3)*cg1(3)   &
           +E(2, 2)*cg2(3)*cg2(3)   &
           +E(2, 3)*cg2(3)*cg3(3)   &
           +E(3, 1)*cg3(3)*cg1(3)   &
           +E(3, 2)*cg3(3)*cg2(3) )
        strain(lx,4)                   &
        = ( E(1, 1)*cg1(1)*cg1(2)   &
           +E(1, 2)*cg1(1)*cg2(2)   &
           +E(1, 3)*cg1(1)*cg3(2)   &
           +E(2, 1)*cg2(1)*cg1(2)   &
           +E(2, 2)*cg2(1)*cg2(2)   &
           +E(2, 3)*cg2(1)*cg3(2)   &
           +E(3, 1)*cg3(1)*cg1(2)   &
           +E(3, 2)*cg3(1)*cg2(2) )
        strain(lx,5)                   &
        = ( E(1, 1)*cg1(2)*cg1(3)   &
           +E(1, 2)*cg1(2)*cg2(3)   &
           +E(1, 3)*cg1(2)*cg3(3)   &
           +E(2, 1)*cg2(2)*cg1(3)   &
           +E(2, 2)*cg2(2)*cg2(3)   &
           +E(2, 3)*cg2(2)*cg3(3)   &
           +E(3, 1)*cg3(2)*cg1(3)   &
           +E(3, 2)*cg3(2)*cg2(3) )
        strain(lx,6)                   &
        = ( E(1, 1)*cg1(3)*cg1(1)   &
           +E(1, 2)*cg1(3)*cg2(1)   &
           +E(1, 3)*cg1(3)*cg3(1)   &
           +E(2, 1)*cg2(3)*cg1(1)   &
           +E(2, 2)*cg2(3)*cg2(1)   &
           +E(2, 3)*cg2(3)*cg3(1)   &
           +E(3, 1)*cg3(3)*cg1(1)   &
           +E(3, 2)*cg3(3)*cg2(1) )

       !--------------------------------------------------

      END DO
!--------------------------------------------------------------------

      RETURN

!####################################################################
      END SUBROUTINE ElementStress_Shell_MITC
!####################################################################


!####################################################################
      SUBROUTINE DL_Shell                                  &
                 (etype, nn, ndof, xx, yy, zz, rho, thick, &
                  ltype, params, vect, nsize, gausses)
!####################################################################

      USE hecmw
      USE m_utilities
      USE mMechGauss
      USE gauss_integration

      TYPE(tGaussStatus), INTENT(IN)   :: gausses(:)
!--------------------------------------------------------------------

      INTEGER(KIND = kint), INTENT(IN) :: etype
      INTEGER(KIND = kint), INTENT(IN) :: nn
      INTEGER(KIND = kint), INTENT(IN) :: ndof
      REAL(KIND = kreal), INTENT(IN) :: xx(*), yy(*), zz(*)
      REAL(KIND = kreal), INTENT(IN) :: rho
      REAL(KIND = kreal), INTENT(IN) :: thick
      REAL(KIND = kreal), INTENT(IN) :: params(*)
      REAL(KIND = kreal), INTENT(OUT) :: vect(*)
      INTEGER(KIND = kint), INTENT(OUT) :: nsize

!--------------------------------------------------------------------

      INTEGER :: ivol, isurf
      INTEGER :: lx, ly
      INTEGER :: fetype
      INTEGER :: ny
      INTEGER :: i, m
      INTEGER :: na, nb
      INTEGER :: isize
      INTEGER :: jsize1, jsize2, jsize3, &
                 jsize4, jsize5, jsize6
      INTEGER :: ltype
      INTEGER :: n_total_layer, n_layer, shellmatl

      REAL(KIND = kreal) :: elem(3, nn)
      REAL(KIND = kreal) :: val
      REAL(KIND = kreal) :: ax, ay, az
      REAL(KIND = kreal) :: rx, ry, rz
      REAL(KIND = kreal) :: xi_lx, eta_lx, zeta_ly
      REAL(KIND = kreal) :: w_w_lx, w_ly
      REAL(KIND = kreal) :: naturalcoord(2)
      REAL(KIND = kreal) :: nncoord(nn, 2)
      REAL(KIND = kreal) :: shapefunc(nn)
      REAL(KIND = kreal) :: shapederiv(nn, 2)
      REAL(KIND = kreal) :: v1(3, nn), v2(3, nn), v3(3, nn)
      REAL(KIND = kreal) :: v1_abs, v2_abs, v3_abs
      REAL(KIND = kreal) :: a_over_2_v3(3, nn)
      REAL(KIND = kreal) :: u_rot(3, nn)
      REAL(KIND = kreal) :: dudxi_rot(3, nn), dudeta_rot(3, nn), &
                            dudzeta_rot(3, nn)
      REAL(KIND = kreal) :: g1(3), g2(3), g3(3)
      REAL(KIND = kreal) :: g1_abs, g2_abs, g3_abs
      REAL(KIND = kreal) :: g1_cross_g2(3)
      REAL(KIND = kreal) :: e_0(3)
      REAL(KIND = kreal) :: cg1(3), cg2(3), cg3(3)
      REAL(KIND = kreal) :: det
      REAL(KIND = kreal) :: det_cg3(3)
      REAL(KIND = kreal) :: det_inv
      REAL(KIND = kreal) :: det_cg3_abs
      REAL(KIND = kreal) :: w_w_w_det
      REAL(KIND = kreal) :: e3_hat(3)
      REAL(KIND = kreal) :: N(3, ndof*nn)
      REAL(KIND = kreal) :: hx, hy, hz
      REAL(KIND = kreal) :: phx, phy, phz
      REAL(KIND = kreal) :: coefx, coefy, coefz
      REAL(KIND = kreal) :: x, y, z
      REAL(KIND = kreal) :: sigma_layer

!--------------------------------------------------------------------

      !   BX   LTYPE=1  :BODY FORCE IN X-DIRECTION
      !   BY   LTYPE=2  :BODY FORCE IN Y-DIRECTION
      !   BZ   LTYPE=3  :BODY FORCE IN Z-DIRECTION
      !   CRAV LTYPE=4  :GRAVITY FORCE
      !   CENT LTYPE=5  :CENTRIFUGAL LOAD
      !   P    LTYPE=10 :TRACTION IN NORMAL-DIRECTION FOR SHELL SURFACE

!--------------------------------------------------------------------

      nsize = ndof*nn

!--------------------------------------------------------------------

      val = params(1)
      ax  = params(2)
      ay  = params(3)
      az  = params(4)
      rx  = params(5)
      ry  = params(6)
      rz  = params(7)

!--------------------------------------------------------------------

      ! MITC4
      IF( etype .EQ. fe_mitc4_shell ) THEN

        fetype = fe_mitc4_shell

        ny = 2

      ! MITC9
      ELSE IF( etype .EQ. fe_mitc9_shell ) THEN

        fetype = fe_mitc9_shell

        ny = 3

      ! MITC3
      ELSE IF( etype .EQ. fe_mitc3_shell ) THEN

        fetype = fe_mitc3_shell

        ny = 2

      END IF

!--------------------------------------------------------------------

      DO na = 1, nn

        elem(1, na) = xx(na)
        elem(2, na) = yy(na)
        elem(3, na) = zz(na)

      END DO

!-------------------------------------------------------------------

      ! xi-coordinate at a node in a local element
      ! eta-coordinate at a node in a local element
      CALL getNodalNaturalCoord(fetype, nncoord)

!--------------------------------------------------------------------

      ! Local load vector
      DO isize = 1, ndof*nn

        vect(isize) = 0.0D0

      END DO

!--------------------------------------------------------------------

      ! xi-coordinate at the center point in a local element
      ! eta-coordinate at the center point in a local element
      naturalcoord(1) = 0.0D0
      naturalcoord(2) = 0.0D0

      CALL getShapeDeriv(fetype, naturalcoord, shapederiv)

      !--------------------------------------------------------------

      ! Covariant basis vector
      DO i = 1, 3

        g1(i) = 0.0D0

        DO na = 1, nn

          g1(i) = g1(i)+shapederiv(na, 1) &
                      *elem(i, na)

        END DO

      END DO

      e_0(1) = g1(1)
      e_0(2) = g1(2)
      e_0(3) = g1(3)

      !--------------------------------------------------------------

      DO nb = 1, nn

        !--------------------------------------------------------

        naturalcoord(1) = nncoord(nb, 1)
        naturalcoord(2) = nncoord(nb, 2)

        CALL getShapeDeriv(fetype, naturalcoord, shapederiv)

       !--------------------------------------------------------

       ! Covariant basis vector
        DO i = 1, 3

          g1(i) = 0.0D0
          g2(i) = 0.0D0

          DO na = 1, nn

            g1(i) = g1(i)+shapederiv(na, 1) &
                       *elem(i, na)
            g2(i) = g2(i)+shapederiv(na, 2) &
                       *elem(i, na)

          END DO

        END DO

       !--------------------------------------------------------

        det_cg3(1) = g1(2)*g2(3)-g1(3)*g2(2)
        det_cg3(2) = g1(3)*g2(1)-g1(1)*g2(3)
        det_cg3(3) = g1(1)*g2(2)-g1(2)*g2(1)

        det_cg3_abs = DSQRT( det_cg3(1)*det_cg3(1)   &
                           +det_cg3(2)*det_cg3(2)   &
                           +det_cg3(3)*det_cg3(3) )

        v3(1, nb) = det_cg3(1)/det_cg3_abs
        v3(2, nb) = det_cg3(2)/det_cg3_abs
        v3(3, nb) = det_cg3(3)/det_cg3_abs

        !--------------------------------------------------------

        v2(1, nb) = v3(2, nb)*e_0(3)-v3(3, nb)*e_0(2)
        v2(2, nb) = v3(3, nb)*e_0(1)-v3(1, nb)*e_0(3)
        v2(3, nb) = v3(1, nb)*e_0(2)-v3(2, nb)*e_0(1)

        v2_abs = DSQRT( v2(1, nb)*v2(1, nb)   &
                      +v2(2, nb)*v2(2, nb)   &
                      +v2(3, nb)*v2(3, nb) )

        IF( v2_abs .GT. 1.0D-15 ) THEN

          v2(1, nb) = v2(1, nb)/v2_abs
          v2(2, nb) = v2(2, nb)/v2_abs
          v2(3, nb) = v2(3, nb)/v2_abs

          v1(1, nb) = v2(2, nb)*v3(3, nb) &
                   -v2(3, nb)*v3(2, nb)
          v1(2, nb) = v2(3, nb)*v3(1, nb) &
                   -v2(1, nb)*v3(3, nb)
          v1(3, nb) = v2(1, nb)*v3(2, nb) &
                   -v2(2, nb)*v3(1, nb)

          v1_abs = DSQRT( v1(1, nb)*v1(1, nb)   &
                       +v1(2, nb)*v1(2, nb)   &
                       +v1(3, nb)*v1(3, nb) )

          v1(1, nb) = v1(1, nb)/v1_abs
          v1(2, nb) = v1(2, nb)/v1_abs
          v1(3, nb) = v1(3, nb)/v1_abs

        ELSE

          v1(1, nb) =  0.0D0
          v1(2, nb) =  0.0D0
          v1(3, nb) = -1.0D0

          v2(1, nb) = 0.0D0
          v2(2, nb) = 1.0D0
          v2(3, nb) = 0.0D0

        END IF

       !--------------------------------------------------------

        v3(1, nb) = v1(2, nb)*v2(3, nb) &
                  -v1(3, nb)*v2(2, nb)
        v3(2, nb) = v1(3, nb)*v2(1, nb) &
                  -v1(1, nb)*v2(3, nb)
        v3(3, nb) = v1(1, nb)*v2(2, nb) &
                  -v1(2, nb)*v2(1, nb)

        v3_abs = DSQRT( v3(1, nb)*v3(1, nb)   &
                      +v3(2, nb)*v3(2, nb)   &
                      +v3(3, nb)*v3(3, nb) )

        v3(1, nb) = v3(1, nb)/v3_abs
        v3(2, nb) = v3(2, nb)/v3_abs
        v3(3, nb) = v3(3, nb)/v3_abs

       !--------------------------------------------------------

        a_over_2_v3(1, nb) = 0.5D0*thick*v3(1, nb)
        a_over_2_v3(2, nb) = 0.5D0*thick*v3(2, nb)
        a_over_2_v3(3, nb) = 0.5D0*thick*v3(3, nb)

       !--------------------------------------------------------

      END DO

!--------------------------------------------------------------------

      ! Selction of load type

      ivol = 0
      isurf = 0

      IF( ltype .LT. 10 ) THEN

        ivol = 1

      ELSE IF( ltype .GE. 10 ) THEN

        isurf = 1

      END IF

!--------------------------------------------------------------------

      !** Surface load
      IF( isurf .EQ. 1 ) THEN

       !--------------------------------------------------------

       DO lx = 1, NumOfQuadPoints(fetype)

        !--------------------------------------------------

        CALL getQuadPoint(fetype, lx, naturalcoord)

        xi_lx  = naturalcoord(1)
        eta_lx = naturalcoord(2)

        w_w_lx = getWeight(fetype, lx)

        CALL getShapeFunc(fetype, naturalcoord, shapefunc)

        CALL getShapeDeriv(fetype, naturalcoord, shapederiv)

        !--------------------------------------------------

        DO na = 1, nn

         DO i = 1, 3

          u_rot(i, na)                    &
          = shapefunc(na)                 &
            *( 0.0D0*a_over_2_v3(i, na) )

          dudxi_rot(i, na)                &
          = shapederiv(na, 1)             &
            *( 0.0D0*a_over_2_v3(i, na) )
          dudeta_rot(i, na)               &
          = shapederiv(na, 2)             &
            *( 0.0D0*a_over_2_v3(i, na) )
          dudzeta_rot(i, na)              &
          = shapefunc(na)                 &
            *( a_over_2_v3(i, na) )

         END DO

        END DO

        !--------------------------------------------------

        ! Covariant basis vector
        DO i = 1, 3

         g1(i) = 0.0D0
         g2(i) = 0.0D0
         !g3(i) = 0.0D0

         DO na = 1, nn

          g1(i) = g1(i)+shapederiv(na, 1)  &
                        *elem(i, na)       &
                       +dudxi_rot(i, na)
          g2(i) = g2(i)+shapederiv(na, 2)  &
                        *elem(i, na)       &
                       +dudeta_rot(i, na)
          !g3(i) = g3(i)+dudzeta_rot(i, na)

         END DO

        END DO

        !--------------------------------------------------

        !g3_abs = DSQRT( g3(1)*g3(1)   &
        !               +g3(2)*g3(2)   &
        !               +g3(3)*g3(3) )

        !--------------------------------------------------

        !e3_hat(1) = g3(1)/g3_abs
        !e3_hat(2) = g3(2)/g3_abs
        !e3_hat(3) = g3(3)/g3_abs

        !--------------------------------------------------

        ! Jacobian
        !det = g1(1)*( g2(2)*g3(3)-g2(3)*g3(2) ) &
        !     +g1(2)*( g2(3)*g3(1)-g2(1)*g3(3) ) &
        !     +g1(3)*( g2(1)*g3(2)-g2(2)*g3(1) )

        !--------------------------------------------------

        g1_cross_g2(1) = g1(2)*g2(3)-g1(3)*g2(2)
        g1_cross_g2(2) = g1(3)*g2(1)-g1(1)*g2(3)
        g1_cross_g2(3) = g1(1)*g2(2)-g1(2)*g2(1)

        !--------------------------------------------------

        DO nb = 1, nn

         jsize1 = ndof*(nb-1)+1
         jsize2 = ndof*(nb-1)+2
         jsize3 = ndof*(nb-1)+3
         jsize4 = ndof*(nb-1)+4
         jsize5 = ndof*(nb-1)+5
         jsize6 = ndof*(nb-1)+6

         N(1, jsize1) = shapefunc(nb)
         N(1, jsize2) = 0.0D0
         N(1, jsize3) = 0.0D0
         N(1, jsize4) = 0.0D0
         N(1, jsize5) = 0.0D0
         N(1, jsize6) = 0.0D0
         N(2, jsize1) = 0.0D0
         N(2, jsize2) = shapefunc(nb)
         N(2, jsize3) = 0.0D0
         N(2, jsize4) = 0.0D0
         N(2, jsize5) = 0.0D0
         N(2, jsize6) = 0.0D0
         N(3, jsize1) = 0.0D0
         N(3, jsize2) = 0.0D0
         N(3, jsize3) = shapefunc(nb)
         N(3, jsize4) = 0.0D0
         N(3, jsize5) = 0.0D0
         N(3, jsize6) = 0.0D0

        END DO

        DO isize = 1, ndof*nn

         vect(isize)                                    &
         = vect(isize)                                  &
          +w_w_lx*( N(1, isize)*g1_cross_g2(1)       &
                      +N(2, isize)*g1_cross_g2(2)       &
                      +N(3, isize)*g1_cross_g2(3) )*val

        END DO

        !--------------------------------------------------

       END DO

       !--------------------------------------------------------

      END IF

!--------------------------------------------------------------------

      !** Volume load
      IF( ivol .EQ. 1 ) THEN

        !--------------------------------------------------------
        n_total_layer =  int(gausses(1)%pMaterial%variables(M_TOTAL_LAYER))
        DO n_layer=1,n_total_layer
          DO ly = 1, ny

            !--------------------------------------------------


            IF ( n_total_layer .ge. 2 ) THEN
              shellmatl = int(gausses(1)%pMaterial%variables(M_SHELL_MATLTYPE))
              if (shellmatl == 0)then
                sigma_layer = 0.0D0
                DO m = 1, n_layer
                  sigma_layer = sigma_layer + 2 * gausses(1)%pMaterial%variables(100+3*m) / thick
                END DO
                zeta_ly = -1 + sigma_layer - gausses(1)%pMaterial%variables(100+3*n_layer) / thick * (1-xg(ny, ly))
              elseif (shellmatl == 1)then
                sigma_layer = 0.0D0
                DO m = 1, n_layer
                  sigma_layer = sigma_layer + 2 * gausses(1)%pMaterial%variables(100+8*m-5) / thick
                END DO
                zeta_ly = -1 + sigma_layer - gausses(1)%pMaterial%variables(100+8*n_layer-5) / thick * (1-xg(ny, ly))
              else
                write(*,*)"ERROR : shellmatl isnot correct"; stop
              endif
            ELSE
              zeta_ly = xg(ny, ly)
            ENDIF

            !zeta_ly = xg(ny, ly)
            w_ly    = wgt(ny, ly)

            !--------------------------------------------------

            DO lx = 1, NumOfQuadPoints(fetype)

              !--------------------------------------------

              CALL getQuadPoint(fetype, lx, naturalcoord)

              xi_lx  = naturalcoord(1)
              eta_lx = naturalcoord(2)

              w_w_lx = getWeight(fetype, lx)

              CALL getShapeFunc(fetype, naturalcoord, shapefunc)

              CALL getShapeDeriv(fetype, naturalcoord, shapederiv)

              !--------------------------------------------

              DO na = 1, nn

                DO i = 1, 3

                  u_rot(i, na)                      &
                    = shapefunc(na)                   &
                    *( zeta_ly*a_over_2_v3(i, na) )

                  dudxi_rot(i, na)                  &
                    = shapederiv(na, 1)               &
                    *( zeta_ly*a_over_2_v3(i, na) )
                  dudeta_rot(i, na)                 &
                    = shapederiv(na, 2)               &
                    *( zeta_ly*a_over_2_v3(i, na) )
                  dudzeta_rot(i, na)                &
                    = shapefunc(na)                   &
                    *( a_over_2_v3(i, na) )

                END DO

              END DO

              !--------------------------------------------

              ! Covariant basis vector
              DO i = 1, 3

                g1(i) = 0.0D0
                g2(i) = 0.0D0
                g3(i) = 0.0D0

                DO na = 1, nn

                  g1(i) = g1(i)+shapederiv(na, 1)  &
                               *elem(i, na)       &
                              +dudxi_rot(i, na)
                  g2(i) = g2(i)+shapederiv(na, 2)  &
                               *elem(i, na)       &
                              +dudeta_rot(i, na)
                  g3(i) = g3(i)+dudzeta_rot(i, na)

                END DO

              END DO

              !--------------------------------------------

              ! Jacobian
              det = g1(1)*( g2(2)*g3(3)-g2(3)*g3(2) ) &
                +g1(2)*( g2(3)*g3(1)-g2(1)*g3(3) ) &
                +g1(3)*( g2(1)*g3(2)-g2(2)*g3(1) )

              !--------------------------------------------

              ! [ N ] matrix
              DO nb = 1, nn

                jsize1 = ndof*(nb-1)+1
                jsize2 = ndof*(nb-1)+2
                jsize3 = ndof*(nb-1)+3
                jsize4 = ndof*(nb-1)+4
                jsize5 = ndof*(nb-1)+5
                jsize6 = ndof*(nb-1)+6

                N(1, jsize1) =  shapefunc(nb)
                N(2, jsize1) =  0.0D0
                N(3, jsize1) =  0.0D0
                N(1, jsize2) =  0.0D0
                N(2, jsize2) =  shapefunc(nb)
                N(3, jsize2) =  0.0D0
                N(1, jsize3) =  0.0D0
                N(2, jsize3) =  0.0D0
                N(3, jsize3) =  shapefunc(nb)
                N(1, jsize4) =  0.0D0
                N(2, jsize4) = -u_rot(3, nb)
                N(3, jsize4) =  u_rot(2, nb)
                N(1, jsize5) =  u_rot(3, nb)
                N(2, jsize5) =  0.0D0
                N(3, jsize5) = -u_rot(1, nb)
                N(1, jsize6) = -u_rot(2, nb)
                N(2, jsize6) =  u_rot(1, nb)
                N(3, jsize6) =  0.0D0

              ENDDO

              !--------------------------------------------

              w_w_w_det = w_w_lx*w_ly*det

              !--------------------------------------------

              IF( ltype .EQ. 1 ) THEN

                DO isize = 1, ndof*nn

                  vect(isize) = vect(isize)+w_w_w_det*N(1, isize)*val

                END DO

              ELSE IF( ltype .EQ. 2 ) THEN

                DO isize = 1, ndof*nn

                  vect(isize) = vect(isize)+w_w_w_det*N(2, isize)*val

                END DO

              ELSE IF( ltype .EQ. 3 ) THEN

                DO isize = 1, ndof*nn

                  vect(isize) = vect(isize)+w_w_w_det*N(3, isize)*val

                END DO

              ELSE IF( ltype .EQ. 4 ) THEN

                DO isize = 1, ndof*nn

                  vect(isize) = vect(isize)+w_w_w_det*rho*ax*N(1, isize)*val
                  vect(isize) = vect(isize)+w_w_w_det*rho*ay*N(2, isize)*val
                  vect(isize) = vect(isize)+w_w_w_det*rho*az*N(3, isize)*val

                END DO

              ELSE IF( ltype .EQ. 5 ) THEN

                x = 0.0D0
                y = 0.0D0
                z = 0.0D0

                DO nb = 1, nn

                  x = x+shapefunc(nb)*elem(1, nb)
                  y = y+shapefunc(nb)*elem(2, nb)
                  z = z+shapefunc(nb)*elem(3, nb)

                END DO

                hx = ax+( (x-ax)*rx+(y-ay)*ry+(z-az)*rz )/( rx**2+ry**2+rz**2 )*rx
                hy = ay+( (x-ax)*rx+(y-ay)*ry+(z-az)*rz )/( rx**2+ry**2+rz**2 )*ry
                hz = az+( (x-ax)*rx+(y-ay)*ry+(z-az)*rz )/( rx**2+ry**2+rz**2 )*rz

                phx = x-hx
                phy = y-hy
                phz = z-hz

                coefx = phx*val*rho*val
                coefy = phy*val*rho*val
                coefz = phz*val*rho*val

                DO isize = 1, ndof*nn

                  vect(isize)                       &
                  = vect(isize)                     &
                  +w_w_w_det*( N(1, isize)*coefx   &
                              +N(2, isize)*coefy   &
                              +N(3, isize)*coefz )

                END DO

              END IF

            !--------------------------------------------

            END DO

          !----------------------------------------------

          END DO

          !----------------------------------------------

        END DO

!--------------------------------------------------------

      END IF
!--------------------------------------------------------------------

      RETURN

!####################################################################
      END SUBROUTINE DL_Shell
!####################################################################


!####################################################################
      SUBROUTINE DL_Shell_33                               &
                 (ic_type, nn, ndof, xx, yy, zz, rho, thick, &
                  ltype, params, vect, nsize, gausses)
!####################################################################

      USE hecmw
      USE m_utilities
      USE mMechGauss
      USE gauss_integration

      TYPE(tGaussStatus) :: gausses(:)
!--------------------------------------------------------------------

      INTEGER(KIND = kint) :: ic_type
      INTEGER(KIND = kint) :: nn
      INTEGER(KIND = kint) :: ndof
      REAL(KIND = kreal) :: xx(*), yy(*), zz(*)
      REAL(KIND = kreal) :: rho
      REAL(KIND = kreal) :: thick
      REAL(KIND = kreal) :: params(*)
      REAL(KIND = kreal) :: vect(*)
      INTEGER(KIND = kint) :: nsize
      INTEGER :: ltype, i, j
      REAL(KIND = kreal) :: tmp(24)
!--------------------------------------------------------------------

      if(ic_type == 761)then
        !ic_type = 731
        !nn = 3
        !ndof = 6
        call DL_Shell(731, 3, 6, xx, yy, zz, rho, thick, ltype, params, vect, nsize, gausses)
        !ic_type = 761
        !nn = 6
        !ndof = 3

        tmp = 0.0
        do i=1,18
          tmp(i) = vect(i)
        enddo

        vect( 1) = tmp(1)
        vect( 2) = tmp(2)
        vect( 3) = tmp(3)
        vect( 4) = tmp(7)
        vect( 5) = tmp(8)
        vect( 6) = tmp(9)
        vect( 7) = tmp(13)
        vect( 8) = tmp(14)
        vect( 9) = tmp(15)
        vect(10) = tmp(4)
        vect(11) = tmp(5)
        vect(12) = tmp(6)
        vect(13) = tmp(10)
        vect(14) = tmp(11)
        vect(15) = tmp(12)
        vect(16) = tmp(16)
        vect(17) = tmp(17)
        vect(18) = tmp(18)

      elseif(ic_type == 781)then
        !ic_type = 741
        !nn = 4
        !ndof = 6
        call DL_Shell(741, 4, 6, xx, yy, zz, rho, thick, ltype, params, vect, nsize, gausses)
        !ic_type = 781
        !nn = 8
        !ndof = 3

        tmp = 0.0
        do i=1,24
          tmp(i) = vect(i)
        enddo

        vect( 1) = tmp(1)
        vect( 2) = tmp(2)
        vect( 3) = tmp(3)
        vect( 4) = tmp(7)
        vect( 5) = tmp(8)
        vect( 6) = tmp(9)
        vect( 7) = tmp(13)
        vect( 8) = tmp(14)
        vect( 9) = tmp(15)
        vect(10) = tmp(19)
        vect(11) = tmp(20)
        vect(12) = tmp(21)
        vect(13) = tmp(4)
        vect(14) = tmp(5)
        vect(15) = tmp(6)
        vect(16) = tmp(10)
        vect(17) = tmp(11)
        vect(18) = tmp(12)
        vect(19) = tmp(16)
        vect(20) = tmp(17)
        vect(21) = tmp(18)
        vect(22) = tmp(22)
        vect(23) = tmp(23)
        vect(24) = tmp(24)

      endif

      END SUBROUTINE DL_Shell_33

!####################################################################
      END MODULE m_static_LIB_shell
