!======================================================================!
!                                                                      !
!    Software Name : FrontCOMP_cure   Ver. 3.1                         !
!                                                                      !
!      Module Name : m_static_NLGEOM_elpl                              !
!      Category    : Static Analysis                                   !
!                                                                      !
!      Developed based on "FrontSTR" of RSS21 project                  !
!                                                                      !
!                     Written by Tomotaka Ogasawara, 2006/06/01        !
!                                Tomotaka Ogasawara, 2013/03/26        !
!                                                                      !
!     Contact address :  IIS,The University of Tokyo, CISS             !
!                                                                      !
!    "Composite Material Strength & Reliability Evaluation Simulator"  !
!                                                                      !
!======================================================================!

!C================================================================C
!C-- module m_static_NLGEOM_elpl
!C================================================================C

module m_static_NLGEOM_elpl

contains



!C***
!C*** subroutines for NONLINEAR SOLID MECHANICS
!C***

!C================================================================C
!C-- subroutine fstrNLGEOM_output_monit
!C================================================================C
      subroutine fstrNLGEOM_output_monit(hecMESH, fstrPARAM, fstrSOLID,    &
                                         elpl_IW4, elpl_IW5, i_step, DT0,  &
                                         my_rank_monit_1, my_rank_monit_2)
!
      use m_fstr
!
      implicit none
!C
!C-- global variable
!C
      type ( hecmwST_local_mesh  ) :: hecMESH
      type ( fstr_solid          ) :: fstrSOLID
      type ( fstr_param          ) :: fstrPARAM
!C
!C-- local variable
!C
      integer(kind=kint) :: elpl_IW4, elpl_IW5
      integer(kind=kint) :: ii,jj
      integer(kind=kint) :: nb_intg, nb_node, nb_elem
      integer(kind=kint) :: i_step
      integer(kind=kint) :: my_rank_monit_1, my_rank_monit_2
      real(kind=kreal) :: DT0
      real(kind=kreal) :: s11, s22, s33, s12, s23, s31, sbar
      real(kind=kreal) :: sd11, sd22, sd33, sig_y
      real(kind=kreal) :: e11, e22, e33, e12, e23, e31, epbar

!!
      nb_intg = fstrSOLID%intg_monit_1
      nb_node = fstrSOLID%node_monit_1
      nb_elem = fstrSOLID%elem_monit_1

      if( hecMESH%my_rank .eq. my_rank_monit_1) then
        jj=fstrPARAM%global_local_id(1,nb_node)
        ii=fstrPARAM%global_local_id(2,nb_node)
!!
!C-- displacement
!!
        if( fstrSOLID%iout_list(1) .eq. 1 ) then
           write(elpl_IW4,'(i10,1pe13.4e3,i10,1p6e13.4e3)') i_step &
                  ,DT0*i_step &
                  ,jj ,fstrSOLID%TOTAL_DISP( hecMESH%n_dof*(ii-1)+1 : hecMESH%n_dof*ii )
        end if
      end if
!
      if( hecMESH%my_rank .eq. my_rank_monit_2) then

        s11 = fstrSOLID%STRESS3D(1,nb_intg,nb_intg,nb_intg,nb_elem)
        s22 = fstrSOLID%STRESS3D(2,nb_intg,nb_intg,nb_intg,nb_elem)
        s33 = fstrSOLID%STRESS3D(3,nb_intg,nb_intg,nb_intg,nb_elem)
        s12 = fstrSOLID%STRESS3D(4,nb_intg,nb_intg,nb_intg,nb_elem)
        s23 = fstrSOLID%STRESS3D(5,nb_intg,nb_intg,nb_intg,nb_elem)
        s31 = fstrSOLID%STRESS3D(6,nb_intg,nb_intg,nb_intg,nb_elem)

        sd11 = s11 -(s11+s22+s33)/3.0
        sd22 = s22 -(s11+s22+s33)/3.0
        sd33 = s33 -(s11+s22+s33)/3.0

        sbar = dsqrt(3.0/2.0*(sd11**2 +sd22**2 +sd33**2 +2.0*s12**2 +2.0*s23**2 +2.0*s31**2))

        e11 = fstrSOLID%STRAIN3D_P(1,nb_intg,nb_intg,nb_intg,nb_elem)
        e22 = fstrSOLID%STRAIN3D_P(2,nb_intg,nb_intg,nb_intg,nb_elem)
        e33 = fstrSOLID%STRAIN3D_P(3,nb_intg,nb_intg,nb_intg,nb_elem)
        e12 = fstrSOLID%STRAIN3D_P(4,nb_intg,nb_intg,nb_intg,nb_elem)
        e23 = fstrSOLID%STRAIN3D_P(5,nb_intg,nb_intg,nb_intg,nb_elem)
        e31 = fstrSOLID%STRAIN3D_P(6,nb_intg,nb_intg,nb_intg,nb_elem)

        epbar = fstrSOLID%EP_BAR(nb_intg,nb_intg,nb_intg,nb_elem)
!!
!C-- stress , strain ,epbar
!!
        if( fstrSOLID%iout_list(6) .eq. 1 ) then
           write(elpl_IW5,'(i10,1pe13.4e3,i10,1p14e13.4e3)') i_step &
                  ,DT0*i_step ,hecMESH%global_elem_ID(nb_elem)      &
!                 ,s11 ,sbar, epbar
                  ,s11, s22, s33, s12, s23, s31, sbar ,e11, e22, e33, e12, e23, e31, epbar
        end if
      end if

      end subroutine fstrNLGEOM_output_monit


!C================================================================C
!C-- subroutine fstrNLGEOM_output_monit_2d
!C================================================================C
      subroutine fstrNLGEOM_output_monit_2d(hecMESH, fstrPARAM, fstrSOLID,    &
                                            elpl_IW4, elpl_IW5, i_step, DT0,  &
                                            my_rank_monit_1, my_rank_monit_2)
!
      use m_fstr
!
      implicit none
!C
!C-- global variable
!C
      type ( hecmwST_local_mesh  ) :: hecMESH
      type ( fstr_solid          ) :: fstrSOLID
      type ( fstr_param          ) :: fstrPARAM
!C
!C-- local variable
!C
      integer(kind=kint) :: elpl_IW4, elpl_IW5
      integer(kind=kint) :: ii,jj
      integer(kind=kint) :: nb_intg, nb_node, nb_elem
      integer(kind=kint) :: i_step
      integer(kind=kint) :: my_rank_monit_1, my_rank_monit_2
      real(kind=kreal) :: DT0
      real(kind=kreal) :: s11, s22, s33, s12, s23, s31, sbar
      real(kind=kreal) :: sd11, sd22, sd33, sig_y
      real(kind=kreal) :: epbar

!!
      nb_intg = fstrSOLID%intg_monit_1
      nb_node = fstrSOLID%node_monit_1
      nb_elem = fstrSOLID%elem_monit_1

      if( hecMESH%my_rank .eq. my_rank_monit_1) then
        jj=fstrPARAM%global_local_id(1,nb_node)
        ii=fstrPARAM%global_local_id(2,nb_node)
!!
!C-- displacement
!!
        if( fstrSOLID%iout_list(1) .eq. 1 ) then
           write(elpl_IW4,'(i10,1pe13.4e3,i10,1p6e13.4e3)') i_step &
                  ,DT0*i_step &
                  ,jj ,fstrSOLID%TOTAL_DISP( hecMESH%n_dof*(ii-1)+1 : hecMESH%n_dof*ii )
        end if
      end if
!
      if( hecMESH%my_rank .eq. my_rank_monit_2) then

        s11 = fstrSOLID%STRESS2D(1,nb_intg,nb_intg,nb_elem)
        s22 = fstrSOLID%STRESS2D(2,nb_intg,nb_intg,nb_elem)
        s33 = fstrSOLID%STRESS2D(4,nb_intg,nb_intg,nb_elem)
        s12 = fstrSOLID%STRESS2D(3,nb_intg,nb_intg,nb_elem)
        s23 = 0.0
        s31 = 0.0

        sd11 = s11 -(s11+s22+s33)/3.0
        sd22 = s22 -(s11+s22+s33)/3.0
        sd33 = s33 -(s11+s22+s33)/3.0

        sbar = dsqrt(3.0/2.0*(sd11**2 +sd22**2 +sd33**2 +2.0*s12**2 +2.0*s23**2 +2.0*s31**2))

        epbar = fstrSOLID%EP_BAR(nb_intg,nb_intg,1,nb_elem)
!!
!C-- stress , strain ,epbar
!!
        if( fstrSOLID%iout_list(6) .eq. 1 ) then
           write(elpl_IW5,'(i10,1pe13.4e3,i10,1p3e13.4e3)') i_step &
                  ,DT0*i_step &
                  ,hecMESH%global_elem_ID(nb_elem) ,s11 ,sbar, epbar
        end if
      end if

      end subroutine fstrNLGEOM_output_monit_2d


!C================================================================C
!C-- subroutine fstr_elpl_INTG
!C
!C-- yield judgement
!C-- stress integration
!C================================================================C
      subroutine fstr_elpl_INTG(DMAT,STR,STR_YLD, FLG_ELPL, h_dash, EP_BAR, EE, PP, DE &
                              , iter,BIJ,MAT_OTHO1,MAT_ISO)
!
      use hecmw
      use m_fstr
!
      implicit none
!C
!C-- I/F variable
!C
      INTEGER(kind=kint) :: iter
      INTEGER(kind=kint) :: FLG_ELPL
      REAL(kind=kreal) :: DMAT(6,6),STR(6),STR_YLD,h_dash,EE,PP
      REAL(kind=kreal) :: DE(6),DS(6)
      REAL(kind=kreal) :: STR_trial(6)
      REAL(kind=kreal) :: EP_BAR
!C
!C-- local variable
!C
      INTEGER(kind=kint) :: i,j,k
      real(kind=kreal) :: GG, GG2
      real(kind=kreal) ::  s11,  s22,  s33,  s12,  s23,  s31, sbar, sbar2
      real(kind=kreal) :: ds11, ds22, ds33, ds12, ds23, ds31
      real(kind=kreal) ::  sd11,  sd22,  sd33, sig_y
      real(kind=kreal) :: dsd11, dsd22, dsd33
      real(kind=kreal) :: x_ganma, x_A, x_B, x_E
      real(kind=kreal) :: depbar




      INTEGER(kind=kint) :: MAT_ISO
      REAL(kind=kreal) :: BIJ(3,3),MAT_OTHO1(9)


!!
!!-- check parameters
!!
      if(h_dash .lt. 0.0) then
         if(myrank .eq. 0) then
            write(imsg,*) 'stop due to fstrSOLID%h_dash < 0.0'
         end if
         call hecmw_abort( hecmw_comm_get_comm())
      end if
!!
      GG = EE/(1.0+PP)/2.0
!!
!!----   elastic
!!
      if( FLG_ELPL .eq. 0) then
         call fstr_elpl_DMAT(DMAT,STR,FLG_ELPL, h_dash, EE, PP,BIJ,MAT_OTHO1,MAT_ISO)

         DS = 0.0
         do i = 1,6
           do j = 1,6
             DS(i)=DS(i)+DMAT(i,j)*DE(j)
           end do
         end do

      STR_trial = STR + DS

      s11 = STR_trial(1)
      s22 = STR_trial(2)
      s33 = STR_trial(3)
      s12 = STR_trial(4)
      s23 = STR_trial(5)
      s31 = STR_trial(6)

      sd11 = s11 -(s11+s22+s33)/3.0
      sd22 = s22 -(s11+s22+s33)/3.0
      sd33 = s33 -(s11+s22+s33)/3.0

      sbar2 = dsqrt(3.0/2.0*(sd11**2 +sd22**2 +sd33**2 +2.0*s12**2 +2.0*s23**2 +2.0*s31**2))

        if( sbar2 .lt. STR_YLD .or. iter .gt. 1) then

           STR = STR_trial

        else if( sbar2 .ge. STR_YLD ) then

           FLG_ELPL = 1

           s11 = STR(1)
           s22 = STR(2)
           s33 = STR(3)
           s12 = STR(4)
           s23 = STR(5)
           s31 = STR(6)
           sd11 = s11 -(s11+s22+s33)/3.0
           sd22 = s22 -(s11+s22+s33)/3.0
           sd33 = s33 -(s11+s22+s33)/3.0

           ds11 = DS(1)
           ds22 = DS(2)
           ds33 = DS(3)
           ds12 = DS(4)
           ds23 = DS(5)
           ds31 = DS(6)
           dsd11 = ds11 -(ds11+ds22+ds33)/3.0
           dsd22 = ds22 -(ds11+ds22+ds33)/3.0
           dsd33 = ds33 -(ds11+ds22+ds33)/3.0

           x_A = 0.5*(dsd11**2 +dsd22**2 +dsd33**2 +2.0*ds12**2 +2.0*ds23**2 &
                     +2.0*ds31**2)
           x_E = 0.5*(sd11**2 +sd22**2 +sd33**2 +2.0*s12**2 +2.0*s23**2 +2.0*s31**2)
           x_E = x_E -1.0/3.0*STR_YLD**2
           x_B = 0.5*(sd11*dsd11 +sd22*dsd22 +sd33*dsd33 +2.0*s12*ds12 +2.0*s23*ds23 &
                     +2.0*s31*ds31)
           x_ganma = (dsqrt(x_B**2-x_A*x_E)-x_B)/x_A

           STR = STR + x_ganma*DS

           call fstr_elpl_DMAT(DMAT,STR,FLG_ELPL, h_dash, EE, PP,BIJ,MAT_OTHO1,MAT_ISO)

           DS = 0.0

           do i = 1,6
             do j = 1,6
               DS(i)=DS(i)+DMAT(i,j)*(1.0-x_ganma)*DE(j)
             end do
           end do

           STR_trial = STR + DS

           s11 = STR_trial(1)
           s22 = STR_trial(2)
           s33 = STR_trial(3)
           s12 = STR_trial(4)
           s23 = STR_trial(5)
           s31 = STR_trial(6)

           sd11 = s11 -(s11+s22+s33)/3.0
           sd22 = s22 -(s11+s22+s33)/3.0
           sd33 = s33 -(s11+s22+s33)/3.0

           sbar2 = dsqrt(3.0/2.0*(sd11**2 +sd22**2 +sd33**2 +2.0*s12**2 +2.0*s23**2 +2.0*s31**2))
!
           DMAT = 0.0
!--  elastic matrix
           DMAT(1,1)=EE*(1.0-PP)/(1.0-2.0*PP)/(1.0+PP)
           DMAT(1,2)=EE*PP/(1.0-2.0*PP)/(1.0+PP)
           DMAT(1,3)=DMAT(1,2)
           DMAT(2,1)=DMAT(1,2)
           DMAT(2,2)=DMAT(1,1)
           DMAT(2,3)=DMAT(1,2)
           DMAT(3,1)=DMAT(1,3)
           DMAT(3,2)=DMAT(2,3)
           DMAT(3,3)=DMAT(1,1)
           DMAT(4,4)=EE/(1.0+PP)*0.5
           DMAT(5,5)=EE/(1.0+PP)*0.5
           DMAT(6,6)=EE/(1.0+PP)*0.5
!
           DS = 0.0

           do i = 1,6
             do j = 1,6
               DS(i)=DS(i)+DMAT(i,j)*(1.0-x_ganma)*DE(j)
             end do
           end do

           ds11 = DS(1)
           ds22 = DS(2)
           ds33 = DS(3)
           ds12 = DS(4)
           ds23 = DS(5)
           ds31 = DS(6)

           depbar = sd11*ds11 +sd22*ds22 +sd33*ds33 +2.0*s12*ds12 +2.0*s23*ds23 &
                   +2.0*s31*ds31
           depbar = depbar/(2.0*GG+2.0/3.0*h_dash)/sbar2

           EP_BAR = EP_BAR + depbar

           sig_y = STR_YLD + h_dash*depbar

           STR = sig_y/sbar2 * STR_trial

           STR_YLD = sig_y

        end if
!!
!!----   plastic
!!
      else if( FLG_ELPL .eq. 1) then
         call fstr_elpl_DMAT(DMAT,STR,FLG_ELPL, h_dash, EE, PP,BIJ,MAT_OTHO1,MAT_ISO)

           DS = 0.0
           do i = 1,6
             do j = 1,6
               DS(i)=DS(i)+DMAT(i,j)*DE(j)
             end do
           end do

           STR_trial = STR + DS

           s11 = STR_trial(1)
           s22 = STR_trial(2)
           s33 = STR_trial(3)
           s12 = STR_trial(4)
           s23 = STR_trial(5)
           s31 = STR_trial(6)

           sd11 = s11 -(s11+s22+s33)/3.0
           sd22 = s22 -(s11+s22+s33)/3.0
           sd33 = s33 -(s11+s22+s33)/3.0

           sbar2 = dsqrt(3.0/2.0*(sd11**2 +sd22**2 +sd33**2 +2.0*s12**2 +2.0*s23**2 +2.0*s31**2))
!
        if( sbar2 .ge. STR_YLD .or. iter .gt. 1) then
!
           DMAT = 0.0
!--  elastic matrix
           DMAT(1,1)=EE*(1.0-PP)/(1.0-2.0*PP)/(1.0+PP)
           DMAT(1,2)=EE*PP/(1.0-2.0*PP)/(1.0+PP)
           DMAT(1,3)=DMAT(1,2)
           DMAT(2,1)=DMAT(1,2)
           DMAT(2,2)=DMAT(1,1)
           DMAT(2,3)=DMAT(1,2)
           DMAT(3,1)=DMAT(1,3)
           DMAT(3,2)=DMAT(2,3)
           DMAT(3,3)=DMAT(1,1)
           DMAT(4,4)=EE/(1.0+PP)*0.5
           DMAT(5,5)=EE/(1.0+PP)*0.5
           DMAT(6,6)=EE/(1.0+PP)*0.5
!
           DS = 0.0
           do i = 1,6
             do j = 1,6
               DS(i)=DS(i)+DMAT(i,j)*DE(j)
             end do
           end do

           ds11 = DS(1)
           ds22 = DS(2)
           ds33 = DS(3)
           ds12 = DS(4)
           ds23 = DS(5)
           ds31 = DS(6)

           depbar = sd11*ds11 +sd22*ds22 +sd33*ds33 +2.0*s12*ds12 +2.0*s23*ds23 &
                   +2.0*s31*ds31
           depbar = depbar/(2.0*GG+2.0/3.0*h_dash)/sbar2

           EP_BAR = EP_BAR + depbar

           sig_y = STR_YLD + h_dash*depbar

           STR = sig_y/sbar2 * STR_trial

           STR_YLD = sig_y

        else if( sbar2 .lt. STR_YLD ) then
           FLG_ELPL = 0

           call fstr_elpl_DMAT(DMAT,STR,FLG_ELPL, h_dash, EE, PP,BIJ,MAT_OTHO1,MAT_ISO)

           DS = 0.0
           do i = 1,6
             do j = 1,6
               DS(i)=DS(i)+DMAT(i,j)*DE(j)
             end do
           end do

           STR = STR + DS

        end if

      end if
!
  100 continue
!
      end subroutine fstr_elpl_INTG


!C================================================================C
!C-- subroutine fstr_elpl_INTG_2d
!C
!C-- yield judgement
!C-- stress integration
!C================================================================C
      subroutine fstr_elpl_INTG_2d(DMAT,STR,STR_YLD, FLG_ELPL, h_dash, EP_BAR, EE, PP, DE &
                                  ,iter, ISET)
!
      use hecmw
      use m_fstr
!
      implicit none
!C
!C-- I/F variable
!C
      INTEGER(kind=kint) :: iter 
      INTEGER(kind=kint) :: FLG_ELPL, ISET
      REAL(kind=kreal) :: DMAT(4,4),STR(4),STR_YLD,h_dash,EE,PP
      REAL(kind=kreal) :: DE(4),DS(4)
      REAL(kind=kreal) :: STR_trial(4)
      REAL(kind=kreal) :: EP_BAR
!C
!C-- local variable
!C
      INTEGER(kind=kint) :: i,j,k
      real(kind=kreal) :: GG, GG2, COEF1, COEF2
      real(kind=kreal) ::  s11,  s22,  s33,  s12,  s23,  s31, sbar, sbar2
      real(kind=kreal) :: ds11, ds22, ds33, ds12, ds23, ds31
      real(kind=kreal) ::  sd11,  sd22,  sd33, sig_y
      real(kind=kreal) :: dsd11, dsd22, dsd33
      real(kind=kreal) :: x_ganma, x_A, x_B, x_E
      real(kind=kreal) :: depbar
!!
!!-- check parameters
!!
      if(h_dash .lt. 0.0) then
         if(myrank .eq. 0) then
            write(imsg,*) 'stop due to fstrSOLID%h_dash < 0.0'
         end if
         call hecmw_abort( hecmw_comm_get_comm())
      end if
!!
      GG = EE/(1.0+PP)/2.0
!!
!!----   elastic
!!
      if( FLG_ELPL .eq. 0) then
         call fstr_elpl_DMAT_2d(DMAT,STR,FLG_ELPL, h_dash, EE, PP,ISET)

         DS = 0.0
         do i = 1,4
           do j = 1,4
             DS(i)=DS(i)+DMAT(i,j)*DE(j)
           end do
         end do

      STR_trial = STR + DS

      s11 = STR_trial(1)
      s22 = STR_trial(2)
      s33 = STR_trial(4)
      s12 = STR_trial(3)
      s23 = 0.0
      s31 = 0.0

      sd11 = s11 -(s11+s22+s33)/3.0
      sd22 = s22 -(s11+s22+s33)/3.0
      sd33 = s33 -(s11+s22+s33)/3.0

      sbar2 = dsqrt(3.0/2.0*(sd11**2 +sd22**2 +sd33**2 +2.0*s12**2 +2.0*s23**2 +2.0*s31**2))

        if( sbar2 .lt. STR_YLD .or. iter .gt. 1) then

           STR = STR_trial

        else if( sbar2 .ge. STR_YLD ) then

           FLG_ELPL = 1

           s11 = STR(1)
           s22 = STR(2)
           s33 = STR(4)
           s12 = STR(3)
           s23 = 0.0
           s31 = 0.0
           sd11 = s11 -(s11+s22+s33)/3.0
           sd22 = s22 -(s11+s22+s33)/3.0
           sd33 = s33 -(s11+s22+s33)/3.0

           ds11 = DS(1)
           ds22 = DS(2)
           ds33 = DS(4)
           ds12 = DS(3)
           ds23 = 0.0
           ds31 = 0.0
           dsd11 = ds11 -(ds11+ds22+ds33)/3.0
           dsd22 = ds22 -(ds11+ds22+ds33)/3.0
           dsd33 = ds33 -(ds11+ds22+ds33)/3.0

           x_A = 0.5*(dsd11**2 +dsd22**2 +dsd33**2 +2.0*ds12**2 +2.0*ds23**2 &
                     +2.0*ds31**2)
           x_E = 0.5*(sd11**2 +sd22**2 +sd33**2 +2.0*s12**2 +2.0*s23**2 +2.0*s31**2)
           x_E = x_E -1.0/3.0*STR_YLD**2
           x_B = 0.5*(sd11*dsd11 +sd22*dsd22 +sd33*dsd33 +2.0*s12*ds12 +2.0*s23*ds23 &
                     +2.0*s31*ds31)
           x_ganma = (dsqrt(x_B**2-x_A*x_E)-x_B)/x_A

           STR = STR + x_ganma*DS

           call fstr_elpl_DMAT_2d(DMAT,STR,FLG_ELPL, h_dash, EE, PP,ISET)

           DS = 0.0

           do i = 1,4
             do j = 1,4
               DS(i)=DS(i)+DMAT(i,j)*(1.0-x_ganma)*DE(j)
             end do
           end do

           STR_trial = STR + DS

           s11 = STR_trial(1)
           s22 = STR_trial(2)
           s33 = STR_trial(4)
           s12 = STR_trial(3)
           s23 = 0.0
           s31 = 0.0

           sd11 = s11 -(s11+s22+s33)/3.0
           sd22 = s22 -(s11+s22+s33)/3.0
           sd33 = s33 -(s11+s22+s33)/3.0

           sbar2 = dsqrt(3.0/2.0*(sd11**2 +sd22**2 +sd33**2 +2.0*s12**2 +2.0*s23**2 +2.0*s31**2))
!
           DMAT = 0.0
!--  elastic matrix
      IF(ISET.EQ.0) THEN
!*PLANE STRAIN
        COEF1=EE/((1.0+PP)*(1.0-2.0*PP))
        COEF2=EE/(2.0*(1.0+PP))
        DMAT(1,1)=COEF1*(1.0-PP)
        DMAT(1,2)=COEF1*PP
        DMAT(1,3)=0.0
        DMAT(2,1)=DMAT(1,2)
        DMAT(2,2)=DMAT(1,1)
        DMAT(2,3)=0.0
        DMAT(3,1)=0.0
        DMAT(3,2)=0.0
        DMAT(3,3)=COEF2
      ELSE IF(ISET.EQ.1) THEN
!*PLANE STRESS
        COEF1=EE/(1.0-PP*PP)
        COEF2=0.5*(1.0-PP)
        DMAT(1,1)=COEF1
        DMAT(1,2)=COEF1*PP
        DMAT(1,3)=0.0
        DMAT(2,1)=DMAT(1,2)
        DMAT(2,2)=DMAT(1,1)
        DMAT(2,3)=0.0
        DMAT(3,1)=0.0
        DMAT(3,2)=0.0
        DMAT(3,3)=COEF1*COEF2
      ELSE IF(ISET.EQ.2) THEN
!*AXISYMMETRIC
        COEF1=EE*(1.0-PP)/((1.0+PP)*(1.0-2.0*PP))
        COEF2=(1.0-2.0*PP)/(2.0*(1.0-PP))
        DMAT(1,1)=COEF1
        DMAT(1,2)=COEF1*PP/(1.0-PP)
        DMAT(1,3)=0.0
        DMAT(1,4)=DMAT(1,2)
        DMAT(2,1)=DMAT(1,2)
        DMAT(2,2)=DMAT(1,1)
        DMAT(2,3)=0.0
        DMAT(2,4)=DMAT(1,2)
        DMAT(3,1)=0.0
        DMAT(3,2)=0.0
        DMAT(3,3)=COEF1*COEF2
        DMAT(3,4)=0.0
        DMAT(4,1)=DMAT(1,4)
        DMAT(4,2)=DMAT(2,4)
        DMAT(4,3)=0.0
        DMAT(4,4)=DMAT(1,1)
      END IF
!
           DS = 0.0

           do i = 1,4
             do j = 1,4
               DS(i)=DS(i)+DMAT(i,j)*(1.0-x_ganma)*DE(j)
             end do
           end do

           ds11 = DS(1)
           ds22 = DS(2)
           ds33 = DS(4)
           ds12 = DS(3)
           ds23 = 0.0
           ds31 = 0.0

           depbar = sd11*ds11 +sd22*ds22 +sd33*ds33 +2.0*s12*ds12 +2.0*s23*ds23 &
                   +2.0*s31*ds31
           depbar = depbar/(2.0*GG+2.0/3.0*h_dash)/sbar2

           EP_BAR = EP_BAR + depbar

           sig_y = STR_YLD + h_dash*depbar

           STR = sig_y/sbar2 * STR_trial

           STR_YLD = sig_y

        end if
!!
!!----   plastic
!!
      else if( FLG_ELPL .eq. 1) then
         call fstr_elpl_DMAT_2d(DMAT,STR,FLG_ELPL, h_dash, EE, PP,ISET)

           DS = 0.0
           do i = 1,4
             do j = 1,4
               DS(i)=DS(i)+DMAT(i,j)*DE(j)
             end do
           end do

           STR_trial = STR + DS

           s11 = STR_trial(1)
           s22 = STR_trial(2)
           s33 = STR_trial(4)
           s12 = STR_trial(3)
           s23 = 0.0
           s31 = 0.0

           sd11 = s11 -(s11+s22+s33)/3.0
           sd22 = s22 -(s11+s22+s33)/3.0
           sd33 = s33 -(s11+s22+s33)/3.0

           sbar2 = dsqrt(3.0/2.0*(sd11**2 +sd22**2 +sd33**2 +2.0*s12**2 +2.0*s23**2 +2.0*s31**2))
!
        if( sbar2 .ge. STR_YLD .or. iter .gt. 1) then
!
           DMAT = 0.0
!--  elastic matrix
      IF(ISET.EQ.0) THEN
!*PLANE STRAIN
        COEF1=EE/((1.0+PP)*(1.0-2.0*PP))
        COEF2=EE/(2.0*(1.0+PP))
        DMAT(1,1)=COEF1*(1.0-PP)
        DMAT(1,2)=COEF1*PP
        DMAT(1,3)=0.0
        DMAT(2,1)=DMAT(1,2)
        DMAT(2,2)=DMAT(1,1)
        DMAT(2,3)=0.0
        DMAT(3,1)=0.0
        DMAT(3,2)=0.0
        DMAT(3,3)=COEF2
      ELSE IF(ISET.EQ.1) THEN
!*PLANE STRESS
        COEF1=EE/(1.0-PP*PP)
        COEF2=0.5*(1.0-PP)
        DMAT(1,1)=COEF1
        DMAT(1,2)=COEF1*PP
        DMAT(1,3)=0.0
        DMAT(2,1)=DMAT(1,2)
        DMAT(2,2)=DMAT(1,1)
        DMAT(2,3)=0.0
        DMAT(3,1)=0.0
        DMAT(3,2)=0.0
        DMAT(3,3)=COEF1*COEF2
      ELSE IF(ISET.EQ.2) THEN
!*AXISYMMETRIC
        COEF1=EE*(1.0-PP)/((1.0+PP)*(1.0-2.0*PP))
        COEF2=(1.0-2.0*PP)/(2.0*(1.0-PP))
        DMAT(1,1)=COEF1
        DMAT(1,2)=COEF1*PP/(1.0-PP)
        DMAT(1,3)=0.0
        DMAT(1,4)=DMAT(1,2)
        DMAT(2,1)=DMAT(1,2)
        DMAT(2,2)=DMAT(1,1)
        DMAT(2,3)=0.0
        DMAT(2,4)=DMAT(1,2)
        DMAT(3,1)=0.0
        DMAT(3,2)=0.0
        DMAT(3,3)=COEF1*COEF2
        DMAT(3,4)=0.0
        DMAT(4,1)=DMAT(1,4)
        DMAT(4,2)=DMAT(2,4)
        DMAT(4,3)=0.0
        DMAT(4,4)=DMAT(1,1)
      END IF
!
           DS = 0.0
           do i = 1,4
             do j = 1,4
               DS(i)=DS(i)+DMAT(i,j)*DE(j)
             end do
           end do

           ds11 = DS(1)
           ds22 = DS(2)
           ds33 = DS(4)
           ds12 = DS(3)
           ds23 = 0.0
           ds31 = 0.0

           depbar = sd11*ds11 +sd22*ds22 +sd33*ds33 +2.0*s12*ds12 +2.0*s23*ds23 &
                   +2.0*s31*ds31
           depbar = depbar/(2.0*GG+2.0/3.0*h_dash)/sbar2

           EP_BAR = EP_BAR + depbar

           sig_y = STR_YLD + h_dash*depbar

           STR = sig_y/sbar2 * STR_trial

           STR_YLD = sig_y

        else if( sbar2 .lt. STR_YLD ) then
           FLG_ELPL = 0

           call fstr_elpl_DMAT_2d(DMAT,STR,FLG_ELPL, h_dash, EE, PP,ISET)

           DS = 0.0
           do i = 1,4
             do j = 1,4
               DS(i)=DS(i)+DMAT(i,j)*DE(j)
             end do
           end do

           STR = STR + DS

        end if

      end if
!
  100 continue
!
      end subroutine fstr_elpl_INTG_2d


!C================================================================C
!C-- subroutine fstr_elpl_DMAT
!C
!C-- set elasto-plastic matrix for stress integration
!C================================================================C
      subroutine fstr_elpl_DMAT(DMAT,STR,FLG_ELPL, h_dash,EE,PP,BIJ,MAT_OTHO1,MAT_ISO)
!
      use hecmw
      use m_fstr
!
      implicit none
!C
!C-- I/F variable
!C
      INTEGER(kind=kint) :: FLG_ELPL
      REAL(kind=kreal) :: DMAT(6,6),STR(6),h_dash,EE,PP
!C
!C-- local variable
!C
      INTEGER(kind=kint) :: i,j,k,l,i2,j2,k2,l2 
      real(kind=kreal) :: GG, GG2
      real(kind=kreal) :: s11, s22, s33, s12, s23, s31, sbar
      real(kind=kreal) :: sd11, sd22, sd33

      INTEGER(kind=kint) :: MAT_ISO
      REAL(kind=kreal) :: BIJ(3,3),MAT_OTHO1(9)
      REAL(kind=kreal) :: CIJKL(3,3,3,3),CIJKLB(3,3,3,3)
      real(kind=kreal) :: E1,E2,E3,G12,G23,G31,nyu12,nyu23,nyu13 
      real(kind=kreal) :: nyu21,nyu32,nyu31,delta1 

!!
      DMAT = 0.0

      do i = 1,3
       do j = 1,3
        do k = 1,3
         do l = 1,3
          CIJKL (i,j,k,l) = 0.0
          CIJKLB(i,j,k,l) = 0.0
         end do
        end do
       end do
      end do

!!
!--  elastic matrix
!
      if( MAT_ISO .eq. 0) then
         DMAT(1,1)=EE*(1.0-PP)/(1.0-2.0*PP)/(1.0+PP)
         DMAT(1,2)=EE*PP/(1.0-2.0*PP)/(1.0+PP)
         DMAT(1,3)=DMAT(1,2)
         DMAT(2,1)=DMAT(1,2)
         DMAT(2,2)=DMAT(1,1)
         DMAT(2,3)=DMAT(1,2)
         DMAT(3,1)=DMAT(1,3)
         DMAT(3,2)=DMAT(2,3)
         DMAT(3,3)=DMAT(1,1)
         DMAT(4,4)=EE/(1.0+PP)*0.5
         DMAT(5,5)=EE/(1.0+PP)*0.5
         DMAT(6,6)=EE/(1.0+PP)*0.5
!     else if( MAT_ISO .eq. 1 .or. MAT_ISO .eq. -1) then
      else if( MAT_ISO .eq. 1) then
         E1    = MAT_OTHO1(1) 
         E2    = MAT_OTHO1(2) 
         E3    = MAT_OTHO1(3) 
         G12   = MAT_OTHO1(4) 
         G23   = MAT_OTHO1(5) 
         G31   = MAT_OTHO1(6) 
         nyu12 = MAT_OTHO1(7) 
         nyu23 = MAT_OTHO1(8) 
         nyu13 = MAT_OTHO1(9) 
         nyu21 = E2/E1*nyu12 
         nyu32 = E3/E2*nyu23 
         nyu31 = E3/E1*nyu13 
         delta1 = (1.0 -nyu12*nyu21 -nyu23*nyu32 -nyu31*nyu13 -2.0*nyu12*nyu23*nyu31)/(E1*E2*E3) 
         CIJKLB(1,1,1,1) = (1.0/E3 -nyu23*nyu23/E2)/(E2*delta1)    
         CIJKLB(2,2,2,2) = (1.0/E3 -nyu13*nyu13/E1)/(E1*delta1)    
         CIJKLB(3,3,3,3) = (1.0/E2 -nyu12*nyu12/E1)/(E1*delta1)    
         CIJKLB(1,1,2,2) = (nyu12/E3 +nyu13*nyu23/E2)/(E1*delta1)    
         CIJKLB(2,2,3,3) = (nyu23/E2 +nyu12*nyu13/E1)/(E1*delta1)    
         CIJKLB(3,3,1,1) = (nyu13 +nyu12*nyu23)/(E1*E2*delta1)    
         CIJKLB(1,2,1,2) = G12 
         CIJKLB(2,3,2,3) = G23 
         CIJKLB(3,1,3,1) = G31 
!
         CIJKLB(2,2,1,1) = CIJKLB(1,1,2,2) 
         CIJKLB(3,3,2,2) = CIJKLB(2,2,3,3)    
         CIJKLB(1,1,3,3) = CIJKLB(3,3,1,1)   

         CIJKLB(2,1,1,2) = CIJKLB(1,2,1,2) 
         CIJKLB(2,1,2,1) = CIJKLB(1,2,1,2) 
         CIJKLB(1,2,2,1) = CIJKLB(1,2,1,2) 
         CIJKLB(3,2,2,3) = CIJKLB(2,3,2,3) 
         CIJKLB(3,2,3,2) = CIJKLB(2,3,2,3) 
         CIJKLB(2,3,3,2) = CIJKLB(2,3,2,3) 
         CIJKLB(1,3,3,1) = CIJKLB(3,1,3,1) 
         CIJKLB(1,3,1,3) = CIJKLB(3,1,3,1) 
         CIJKLB(3,1,1,3) = CIJKLB(3,1,3,1) 
!
      do i = 1,3
       do j = 1,3
        do k = 1,3
         do l = 1,3
          do i2 = 1,3
           do j2 = 1,3
            do k2 = 1,3
             do l2 = 1,3
              CIJKL(i,j,k,l) = CIJKL(i,j,k,l) + BIJ(i,i2)*BIJ(j,j2)*BIJ(k,k2)*BIJ(l,l2)*CIJKLB(i2,j2,k2,l2)
             end do
            end do
           end do
          end do
         end do
        end do
       end do
      end do
!
         DMAT(1,1)=CIJKL(1,1,1,1)
         DMAT(1,2)=CIJKL(1,1,2,2)
         DMAT(1,3)=CIJKL(1,1,3,3)
         DMAT(2,1)=CIJKL(2,2,1,1)
         DMAT(2,2)=CIJKL(2,2,2,2)
         DMAT(2,3)=CIJKL(2,2,3,3)
         DMAT(3,1)=CIJKL(3,3,1,1)
         DMAT(3,2)=CIJKL(3,3,2,2)
         DMAT(3,3)=CIJKL(3,3,3,3)
         DMAT(4,4)=CIJKL(1,2,1,2)
         DMAT(5,5)=CIJKL(2,3,2,3)
         DMAT(6,6)=CIJKL(3,1,3,1)

      else if( MAT_ISO .eq. 2) then
      end if
!!
!--  elasto-plastic matrix
!
      if( FLG_ELPL .eq. 1) then

      GG = EE/(1.0+PP)/2.0
!
      s11 = STR(1)
      s22 = STR(2)
      s33 = STR(3)
      s12 = STR(4)
      s23 = STR(5)
      s31 = STR(6)

      sd11 = s11 -(s11+s22+s33)/3.0
      sd22 = s22 -(s11+s22+s33)/3.0
      sd33 = s33 -(s11+s22+s33)/3.0

      sbar = dsqrt(3.0/2.0*(sd11**2 +sd22**2 +sd33**2 +2.0*s12**2 +2.0*s23**2 +2.0*s31**2))

      GG2 = 9.0*GG**2/(3.0*GG+h_dash)/sbar**2
!!
        DMAT(1,1) = DMAT(1,1) - GG2*sd11*sd11
        DMAT(1,2) = DMAT(1,2) - GG2*sd11*sd22
        DMAT(1,3) = DMAT(1,3) - GG2*sd11*sd33
        DMAT(1,4) = DMAT(1,4) - GG2*sd11*s12
        DMAT(1,5) = DMAT(1,5) - GG2*sd11*s23
        DMAT(1,6) = DMAT(1,6) - GG2*sd11*s31

        DMAT(2,1) = DMAT(2,1) - GG2*sd22*sd11
        DMAT(2,2) = DMAT(2,2) - GG2*sd22*sd22
        DMAT(2,3) = DMAT(2,3) - GG2*sd22*sd33
        DMAT(2,4) = DMAT(2,4) - GG2*sd22*s12
        DMAT(2,5) = DMAT(2,5) - GG2*sd22*s23
        DMAT(2,6) = DMAT(2,6) - GG2*sd22*s31

        DMAT(3,1) = DMAT(3,1) - GG2*sd33*sd11
        DMAT(3,2) = DMAT(3,2) - GG2*sd33*sd22
        DMAT(3,3) = DMAT(3,3) - GG2*sd33*sd33
        DMAT(3,4) = DMAT(3,4) - GG2*sd33*s12
        DMAT(3,5) = DMAT(3,5) - GG2*sd33*s23
        DMAT(3,6) = DMAT(3,6) - GG2*sd33*s31

        DMAT(4,1) = DMAT(4,1) - GG2*s12 *sd11
        DMAT(4,2) = DMAT(4,2) - GG2*s12 *sd22
        DMAT(4,3) = DMAT(4,3) - GG2*s12 *sd33
        DMAT(4,4) = DMAT(4,4) - GG2*s12 *s12
        DMAT(4,5) = DMAT(4,5) - GG2*s12 *s23
        DMAT(4,6) = DMAT(4,6) - GG2*s12 *s31

        DMAT(5,1) = DMAT(5,1) - GG2*s23 *sd11
        DMAT(5,2) = DMAT(5,2) - GG2*s23 *sd22
        DMAT(5,3) = DMAT(5,3) - GG2*s23 *sd33
        DMAT(5,4) = DMAT(5,4) - GG2*s23 *s12
        DMAT(5,5) = DMAT(5,5) - GG2*s23 *s23
        DMAT(5,6) = DMAT(5,6) - GG2*s23 *s31

        DMAT(6,1) = DMAT(6,1) - GG2*s31 *sd11
        DMAT(6,2) = DMAT(6,2) - GG2*s31 *sd22
        DMAT(6,3) = DMAT(6,3) - GG2*s31 *sd33
        DMAT(6,4) = DMAT(6,4) - GG2*s31 *s12
        DMAT(6,5) = DMAT(6,5) - GG2*s31 *s23
        DMAT(6,6) = DMAT(6,6) - GG2*s31 *s31
!
      end if
!
      end subroutine fstr_elpl_DMAT


!C================================================================C
!C-- subroutine fstr_elpl_DMAT_2d
!C
!C-- set elasto-plastic matrix for stress integration
!C================================================================C
      subroutine fstr_elpl_DMAT_2d(DMAT,STR,FLG_ELPL, h_dash, EE, PP, ISET)
!
      use hecmw
!
      implicit none
!C
!C-- I/F variable
!C
      INTEGER(kind=kint) :: FLG_ELPL, ISET
      REAL(kind=kreal) :: DMAT(4,4),STR(4),h_dash,EE,PP
!C
!C-- local variable
!C
      real(kind=kreal) :: GG, GG2, COEF1, COEF2
      real(kind=kreal) :: s11, s22, s33, s12, s23, s31, sbar
      real(kind=kreal) :: sd11, sd22, sd33
      real(kind=kreal) :: ss11, ss22, ss12
!!
         DMAT = 0.0
!!
!--  elastic matrix
!
      IF(ISET.EQ.0) THEN
!*PLANE STRAIN
        COEF1=EE/((1.0+PP)*(1.0-2.0*PP))
        COEF2=EE/(2.0*(1.0+PP))
        DMAT(1,1)=COEF1*(1.0-PP)
        DMAT(1,2)=COEF1*PP
        DMAT(1,3)=0.0
        DMAT(2,1)=DMAT(1,2)
        DMAT(2,2)=DMAT(1,1)
        DMAT(2,3)=0.0
        DMAT(3,1)=0.0
        DMAT(3,2)=0.0
        DMAT(3,3)=COEF2
      ELSE IF(ISET.EQ.1) THEN
!*PLANE STRESS
        COEF1=EE/(1.0-PP*PP)
        COEF2=0.5*(1.0-PP)
        DMAT(1,1)=COEF1
        DMAT(1,2)=COEF1*PP
        DMAT(1,3)=0.0
        DMAT(2,1)=DMAT(1,2)
        DMAT(2,2)=DMAT(1,1)
        DMAT(2,3)=0.0
        DMAT(3,1)=0.0
        DMAT(3,2)=0.0
        DMAT(3,3)=COEF1*COEF2
      ELSE IF(ISET.EQ.2) THEN
!*AXISYMMETRIC
        COEF1=EE*(1.0-PP)/((1.0+PP)*(1.0-2.0*PP))
        COEF2=(1.0-2.0*PP)/(2.0*(1.0-PP))
        DMAT(1,1)=COEF1
        DMAT(1,2)=COEF1*PP/(1.0-PP)
        DMAT(1,3)=0.0
        DMAT(1,4)=DMAT(1,2)
        DMAT(2,1)=DMAT(1,2)
        DMAT(2,2)=DMAT(1,1)
        DMAT(2,3)=0.0
        DMAT(2,4)=DMAT(1,2)
        DMAT(3,1)=0.0
        DMAT(3,2)=0.0
        DMAT(3,3)=COEF1*COEF2
        DMAT(3,4)=0.0
        DMAT(4,1)=DMAT(1,4)
        DMAT(4,2)=DMAT(2,4)
        DMAT(4,3)=0.0
        DMAT(4,4)=DMAT(1,1)
      END IF
!!
!--  elasto-plastic matrix
!
      if( FLG_ELPL .eq. 1) then
!
       IF(ISET.EQ.0) THEN
!*PLANE STRAIN

      GG = EE/(1.0+PP)/2.0
!
      s11 = STR(1)
      s22 = STR(2)
      s33 = STR(4)
      s12 = STR(3)
      s23 = 0.0
      s31 = 0.0

      sd11 = s11 -(s11+s22+s33)/3.0
      sd22 = s22 -(s11+s22+s33)/3.0
      sd33 = s33 -(s11+s22+s33)/3.0

      sbar = dsqrt(3.0/2.0*(sd11**2 +sd22**2 +sd33**2 +2.0*s12**2 +2.0*s23**2 +2.0*s31**2))

      GG2 = 9.0*GG**2/(3.0*GG+h_dash)/sbar**2
!!
        DMAT(1,1) = DMAT(1,1) - GG2*sd11*sd11
        DMAT(1,2) = DMAT(1,2) - GG2*sd11*sd22
        DMAT(1,3) = DMAT(1,3) - GG2*sd11*s12

        DMAT(2,1) = DMAT(2,1) - GG2*sd22*sd11
        DMAT(2,2) = DMAT(2,2) - GG2*sd22*sd22
        DMAT(2,3) = DMAT(2,3) - GG2*sd22*s12

        DMAT(3,1) = DMAT(3,1) - GG2*s12 *sd11
        DMAT(3,2) = DMAT(3,2) - GG2*s12 *sd22
        DMAT(3,3) = DMAT(3,3) - GG2*s12 *s12

       ELSE IF(ISET.EQ.1) THEN
!*PLANE STRESS

      GG = EE/(1.0+PP)/2.0
!
      s11 = STR(1)
      s22 = STR(2)
      s33 = STR(4)
      s12 = STR(3)
      s23 = 0.0
      s31 = 0.0

      sd11 = s11 -(s11+s22+s33)/3.0
      sd22 = s22 -(s11+s22+s33)/3.0
      sd33 = s33 -(s11+s22+s33)/3.0

      sbar = dsqrt(3.0/2.0*(sd11**2 +sd22**2 +sd33**2 +2.0*s12**2 +2.0*s23**2 +2.0*s31**2))

      ss11 = DMAT(1,1)*sd11 + DMAT(1,2)*sd22 + 2.0*DMAT(1,3)*s12
      ss22 = DMAT(2,1)*sd11 + DMAT(2,2)*sd22 + 2.0*DMAT(2,3)*s12
      ss12 = DMAT(3,1)*sd11 + DMAT(3,2)*sd22 + 2.0*DMAT(3,3)*s12

      GG2 = 1.0 /((ss11*sd11 + ss22*sd22 + 2.0*ss12*s12) + 4.0/9.0*(sbar**2)*h_dash)
!!
        DMAT(1,1) = DMAT(1,1) - GG2*ss11*ss11
        DMAT(1,2) = DMAT(1,2) - GG2*ss11*ss22
        DMAT(1,3) = DMAT(1,3) - GG2*ss11*ss12

        DMAT(2,1) = DMAT(2,1) - GG2*ss22*ss11
        DMAT(2,2) = DMAT(2,2) - GG2*ss22*ss22
        DMAT(2,3) = DMAT(2,3) - GG2*ss22*ss12

        DMAT(3,1) = DMAT(3,1) - GG2*ss12*ss11
        DMAT(3,2) = DMAT(3,2) - GG2*ss12*ss22
        DMAT(3,3) = DMAT(3,3) - GG2*ss12*ss12

       ELSE IF(ISET.EQ.2) THEN
!*AXISYMMETRIC

      GG = EE/(1.0+PP)/2.0
!
      s11 = STR(1)
      s22 = STR(2)
      s33 = STR(4)
      s12 = STR(3)
      s23 = 0.0
      s31 = 0.0

      sd11 = s11 -(s11+s22+s33)/3.0
      sd22 = s22 -(s11+s22+s33)/3.0
      sd33 = s33 -(s11+s22+s33)/3.0

      sbar = dsqrt(3.0/2.0*(sd11**2 +sd22**2 +sd33**2 +2.0*s12**2 +2.0*s23**2 +2.0*s31**2))

      GG2 = 9.0*GG**2/(3.0*GG+h_dash)/sbar**2
!!
        DMAT(1,1) = DMAT(1,1) - GG2*sd11*sd11
        DMAT(1,2) = DMAT(1,2) - GG2*sd11*sd22
        DMAT(1,3) = DMAT(1,3) - GG2*sd11*s12
        DMAT(1,4) = DMAT(1,4) - GG2*sd11*sd33

        DMAT(2,1) = DMAT(2,1) - GG2*sd22*sd11
        DMAT(2,2) = DMAT(2,2) - GG2*sd22*sd22
        DMAT(2,3) = DMAT(2,3) - GG2*sd22*s12
        DMAT(2,4) = DMAT(2,4) - GG2*sd22*sd33

        DMAT(3,1) = DMAT(3,1) - GG2*s12 *sd11
        DMAT(3,2) = DMAT(3,2) - GG2*s12 *sd22
        DMAT(3,3) = DMAT(3,3) - GG2*s12 *s12
        DMAT(3,4) = DMAT(3,4) - GG2*s12 *sd33

        DMAT(4,1) = DMAT(4,1) - GG2*sd33*sd11
        DMAT(4,2) = DMAT(4,2) - GG2*sd33*sd22
        DMAT(4,3) = DMAT(4,3) - GG2*sd33*s12
        DMAT(4,4) = DMAT(4,4) - GG2*sd33*sd33

       END IF
!
      end if
!
      end subroutine fstr_elpl_DMAT_2d


!C================================================================C
!C-- subroutine table_nlsta
!C================================================================C
    subroutine table_nlsta(hecMESH, fstrSOLID, i_step, ig0, f_t, flag_u)

    use hecmw
    use m_fstr

    implicit none
!C
!C-- global variable
!C
      type ( hecmwST_local_mesh  ) :: hecMESH
      type ( fstr_solid          ) :: fstrSOLID
!C
!C-- local variable
!C
    integer(kind=kint) :: i, ig0, i_step
    integer(kind=kint) :: jj_n_amp, jj1, jj2
    integer(kind=kint) :: s1, s2, flag, flag_u
    real(kind=kreal) :: t_1, t_2, t_t, f_1, f_2, f_t

!C
!C-- zero clear
!C
      i = 0
      s1 = 0
      s2 = 0
      flag = 0
      t_1 = 0.0
      t_2 = 0.0
      t_t = 0.0
      f_1 = 0.0
      f_2 = 0.0
      f_t = 0.0
!!
!      fstrSOLID%NLSTATIC_ngrp_amp = 1
!         jj_n_amp = -1
!!

          jj_n_amp = fstrSOLID%NLSTATIC_ngrp_amp

      if( jj_n_amp .eq. 0 ) then
          t_t = DT*i_step
          f_t = t_t

      else
          jj1 = hecMESH%amp%amp_index(jj_n_amp - 1)
          jj2 = hecMESH%amp%amp_index(jj_n_amp)

      jj1 = jj1 + 2
          t_t = DT*i_step

!      if(jj2 .eq. 0) then
!         f_t = 1.0
      if(t_t .gt. hecMESH%amp%amp_table(jj2)) then
         f_t = hecMESH%amp%amp_val(jj2)
      else if(t_t .le. hecMESH%amp%amp_table(jj2)) then
         do i = jj1, jj2
            if(t_t .le. hecMESH%amp%amp_table(i)) then
              s2 = i
              s1 = i - 1
              flag = 1
            end if
            if( flag .eq. 1 ) exit
         end do

         t_2 = hecMESH%amp%amp_table(s2)
         t_1 = hecMESH%amp%amp_table(s1)
         f_2 = hecMESH%amp%amp_val(s2)
         f_1 = hecMESH%amp%amp_val(s1)
            if( t_2-t_1 .lt. 1.0e-20) then
              if(myrank .eq. 0) then
                 write(imsg,*) 'stop due to t_2-t_1 <= 0'
              end if
              call hecmw_abort( hecmw_comm_get_comm())
            end if
         f_t = ((t_2*f_1 - t_1*f_2) + (f_2 - f_1)*t_t) / (t_2 - t_1)
      end if

      end if

  end subroutine table_nlsta

!C================================================================C
!C-- subroutine fstrNLGEOM_get_prop 
!C================================================================C
!---   get material properties

   SUBROUTINE fstrNLGEOM_get_prop(hecMESH,isect,ee,pp,rho,alpha,thick)

      use m_fstr
      
      IMPLICIT REAL(kind=kreal) (A-H,O-Z)
      TYPE (hecmwST_local_mesh) :: hecMESH

!*EHM June 24 04
!Get thickness
      n_item = hecMESH%section%sect_R_index(isect)-hecMESH%section%sect_R_index(isect-1)
      !Print *,'n_item for thickness: ',n_item
      !PAUSE
      ihead = hecMESH%section%sect_R_index(isect-1)
      !Print *,'ihead for thickness: ',ihead
      !PAUSE
      !do i = 1,n_item
       thick = hecMESH%section%sect_R_item(ihead+1)
       IF(thick.LE.0.0) STOP "Zero thickness <= 0 is illegal"
       !Print *,'cval:',cval
       !PAUSE
      !end do
       
!C** material ID
      mid=hecMESH%section%sect_mat_ID_item(isect)

!C** Number of Item
      n_item=hecMESH%material%mat_ITEM_index(mid)-hecMESH%material%mat_ITEM_index(mid-1) 
      
!C** Head possition
      ihead=hecMESH%material%mat_ITEM_index(mid-1)
!C Get ITEM of Meterial (Young's Modulus & Possion's Ratio
      pp=0.0
      if( n_item .lt. 1 ) then
        write(IMSG,*) 'n_item=',n_item
        write(IMSG,*) '###Error 1'
        stop
      endif
      if ( n_item.ge.1 ) then
!C** Number of Sub Item
        n_subitem = hecMESH%material%mat_subITEM_index(ihead+1)                          &
                  - hecMESH%material%mat_subITEM_index(ihead)   
!C** Head possition
        mpos=hecMESH%material%mat_subITEM_index(ihead)
!C Get SUBITEM of Meterial
        if( n_subitem .lt. 1 ) then
          write(IMSG,*) '###Error 2'
          stop
        endif
        if ( n_subitem.ge.1 ) then
          ee=hecMESH%material%mat_val(mpos+1)
        endif
        if ( n_subitem.ge.2 ) then
          pp=hecMESH%material%mat_val(mpos+2)
        endif
      endif
!C Get ITEM of Meterial (Density)
      rho=0.0
      if ( n_item.ge.2 ) then
!C** Number of Sub Item
        n_subitem=hecMESH%material%mat_subITEM_index(ihead+2)                            &
                 -hecMESH%material%mat_subITEM_index(ihead+1)
!C** Head possition
        mpos=hecMESH%material%mat_subITEM_index(ihead+1)
!C Get SUBITEM of Meterial
        if( n_subitem .lt. 1 ) then
          write(IMSG,*) '###Error 3'
          stop
        endif
        if ( n_subitem.ge.1 ) then
          rho=hecMESH%material%mat_val(mpos+1)
        endif
      endif
!C Get ITEM of Meterial (Thermal Expansion)
      alpha=0.0
      if ( n_item.ge.3 ) then
!C** Number of Sub Item
        n_subitem=hecMESH%material%mat_subITEM_index(ihead+3)                            &      
                 -hecMESH%material%mat_subITEM_index(ihead+2)
!C** Head possition
        mpos=hecMESH%material%mat_subITEM_index(ihead+2)
!C Get SUBITEM of Meterial
        if( n_subitem .lt. 1 ) then
          write(IMSG,*) '###Error 4'
          stop
        endif
        if ( n_subitem.ge.1 ) then
          alpha=hecMESH%material%mat_val(mpos+1)
        endif
      endif
   end subroutine fstrNLGEOM_get_prop

!C================================================================C
!C-- subroutine fstr_elpl_eig_jcb
!C
!C-- eigen solver for 3 x 3 matrix
!C--  jacobi method
!C================================================================C
      subroutine fstr_elpl_eig_jcb(ARRAY_2, ARRAY_PRIS_2,i_step,G_ELM_ID)
!
      use hecmw
      use m_fstr
!
      implicit none
!C
!C-- I/F variable
!C
      REAL(kind=kreal) :: ARRAY_2(6),ARRAY_PRIS_2(12)
!C
!C-- local variable
!C
      INTEGER(kind=kint), parameter :: n    = 3 
      INTEGER(kind=kint), parameter :: imax = 100 
      real(kind=kreal),   parameter :: eps_eig  = 1.0e-06 

      INTEGER(kind=kint) :: i,j,k,l,m,i2,j2,k2,l2 
      INTEGER(kind=kint) :: ic, ic2, i_step, G_ELM_ID 
      INTEGER(kind=kint) :: idx_2(3), idx_2_iv(3) 
      real(kind=kreal) :: a(n,n), x(n,n), amax, akm, q, co, si, ww 
      real(kind=kreal) :: amax2 
!!

      do i = 1,12
       ARRAY_PRIS_2(i) = 0.0
      end do
      do i = 1,3
       do j = 1,3
        x(i,j) = 0.0
       end do
      end do
      a(1,1) = ARRAY_2(1) 
      a(1,2) = ARRAY_2(4) 
      a(1,3) = ARRAY_2(6) 
      a(2,1) = ARRAY_2(4) 
      a(2,2) = ARRAY_2(2) 
      a(2,3) = ARRAY_2(5) 
      a(3,1) = ARRAY_2(6) 
      a(3,2) = ARRAY_2(5) 
      a(3,3) = ARRAY_2(3) 
!!!
      do k = 1, n
       do m = 1, n
        if(k .eq. m) then
         x(k,m) = 1.0
        else
         x(k,m) = 0.0
        end if
       end do
      end do
!
      do ic = 1, imax
       ic2 = 0
       amax = 0.0
!
       do k = 1, n-1
        do m = k+1, n
         if(m .le. n) then
           akm =abs(a(k,m)) 
           if(akm .gt. amax) then
             i = k
             j = m
             amax = akm 
           end if
         end if
        end do
       end do
!
       if(amax .lt. eps_eig) then
        ic2 = 1
        exit
       end if
!
       if(a(i,i) .eq. a(j,j)) then
        q = 0.5 * acos(0.0)
       else
        q = 0.5 * atan(2.0*a(j,i)/(a(i,i)-a(j,j)))
       end if
       co = cos(q)
       si = sin(q)
!
       do k = 1, n
         ww = a(i,k) 
         a(i,k) = ww*co + a(j,k)*si
         a(j,k) =-ww*si + a(j,k)*co
         if((k.ne.i) .and. (k.ne.j)) then
           a(k,i) = a(i,k)
           a(k,j) = a(j,k)
         end if
       end do
       a(i,i) = a(i,i)*co + a(i,j)*si
       a(j,j) =-a(j,i)*si + a(j,j)*co
       a(i,j) = 0.0
       a(j,i) = 0.0
!
       do k = 1, n
         ww = x(k,i) 
         x(k,i) = ww*co + x(k,j)*si
         x(k,j) =-ww*si + x(k,j)*co
       end do
      end do 
!
      if(i_step .ge. -100) then
       if(ic2 .eq. 0) then
         write(idbg,'(''not convergence at eigen solver of principal value: ''           &
                     ,''  step ='',i10,'' global_elem_ID(ic) ='',i10,'' imax ='',i10 )') &
                      i_step, G_ELM_ID, ic 
       end if
      end if
!!!
      idx_2(1) = 1
      idx_2_iv(1) = 1
      amax2 = a(1,1)

      do i = 2,3
       if(a(i,i) .gt. amax2) then
        idx_2(i) = 1
        idx_2_iv(1) = i
        amax2 = a(i,i)
       end if
      end do

      if(idx_2_iv(1) .eq. 1) then
       if(a(2,2) .gt. a(3,3)) then
        idx_2(2) = 2
        idx_2(3) = 3
        idx_2_iv(2) = 2
        idx_2_iv(3) = 3
       else if(a(2,2) .le. a(3,3)) then
        idx_2(2) = 3
        idx_2(3) = 2
        idx_2_iv(2) = 3
        idx_2_iv(3) = 2
       end if
      else if(idx_2_iv(1) .eq. 2) then
       if(a(3,3) .gt. a(1,1)) then
        idx_2(3) = 2
        idx_2(1) = 3
        idx_2_iv(2) = 3
        idx_2_iv(3) = 1
       else if(a(3,3) .le. a(1,1)) then
        idx_2(3) = 3
        idx_2(1) = 2
        idx_2_iv(2) = 1
        idx_2_iv(3) = 3
       end if
      else if(idx_2_iv(1) .eq. 3) then
       if(a(1,1) .gt. a(2,2)) then
        idx_2(1) = 2
        idx_2(2) = 3
        idx_2_iv(2) = 1
        idx_2_iv(3) = 2
       else if(a(1,1) .le. a(2,2)) then
        idx_2(1) = 3
        idx_2(2) = 2
        idx_2_iv(2) = 2
        idx_2_iv(3) = 1
       end if
      end if

      do i = 1,3
       ARRAY_PRIS_2( i)  = a(idx_2_iv(i),idx_2_iv(i)) 
      end do
      do i = 1,3
       ARRAY_PRIS_2(3+i) = x(i,idx_2_iv(1)) 
       ARRAY_PRIS_2(6+i) = x(i,idx_2_iv(2)) 
       ARRAY_PRIS_2(9+i) = x(i,idx_2_iv(3)) 
      end do

!!     do i = 1,3
!       ARRAY_PRIS_2( i)  = a(i,i) 
!      end do
!      do i = 1,3
!       ARRAY_PRIS_2(3+i) = x(i,1) 
!       ARRAY_PRIS_2(6+i) = x(i,2) 
!       ARRAY_PRIS_2(9+i) = x(i,3) 
!!     end do

!
      end subroutine fstr_elpl_eig_jcb 

!C================================================================C
!C-- subroutine fstr_elpl_eig_qr 
!C
!C-- eigen solver for 3 x 3 matrix
!C--  jacobi method  (2011.06.01)
!C--  QR method      (not yet included) 
!C================================================================C
      subroutine fstr_elpl_eig_qr(ARRAY_2, ARRAY_PRIS_2,i_step,G_ELM_ID)
!
      use hecmw
      use m_fstr
!
      implicit none
!C
!C-- I/F variable
!C
      REAL(kind=kreal) :: ARRAY_2(6),ARRAY_PRIS_2(12)
!C
!C-- local variable
!C
      INTEGER(kind=kint), parameter :: n    = 3 
      INTEGER(kind=kint), parameter :: imax = 100 
      real(kind=kreal),   parameter :: eps_eig  = 1.0e-06 

      INTEGER(kind=kint) :: i,j,k,l,m,i2,j2,k2,l2 
      INTEGER(kind=kint) :: ic, ic2, i_step, G_ELM_ID 
      INTEGER(kind=kint) :: idx_2(3), idx_2_iv(3) 
      real(kind=kreal) :: a(n,n), x(n,n), amax, akm, q, co, si, ww 
      real(kind=kreal) :: amax2 
!!

      do i = 1,12
       ARRAY_PRIS_2(i) = 0.0
      end do
      do i = 1,3
       do j = 1,3
        x(i,j) = 0.0
       end do
      end do
      a(1,1) = ARRAY_2(1) 
      a(1,2) = ARRAY_2(4) 
      a(1,3) = ARRAY_2(6) 
      a(2,1) = ARRAY_2(4) 
      a(2,2) = ARRAY_2(2) 
      a(2,3) = ARRAY_2(5) 
      a(3,1) = ARRAY_2(6) 
      a(3,2) = ARRAY_2(5) 
      a(3,3) = ARRAY_2(3) 
!!!
      do k = 1, n
       do m = 1, n
        if(k .eq. m) then
         x(k,m) = 1.0
        else
         x(k,m) = 0.0
        end if
       end do
      end do
!
      do ic = 1, imax
       ic2 = 0
       amax = 0.0
!
       do k = 1, n-1
        do m = k+1, n
         if(m .le. n) then
           akm =abs(a(k,m)) 
           if(akm .gt. amax) then
             i = k
             j = m
             amax = akm 
           end if
         end if
        end do
       end do
!
       if(amax .lt. eps_eig) then
        ic2 = 1
        exit
       end if
!
       if(a(i,i) .eq. a(j,j)) then
        q = 0.5 * acos(0.0)
       else
        q = 0.5 * atan(2.0*a(j,i)/(a(i,i)-a(j,j)))
       end if
       co = cos(q)
       si = sin(q)
!
       do k = 1, n
         ww = a(i,k) 
         a(i,k) = ww*co + a(j,k)*si
         a(j,k) =-ww*si + a(j,k)*co
         if((k.ne.i) .and. (k.ne.j)) then
           a(k,i) = a(i,k)
           a(k,j) = a(j,k)
         end if
       end do
       a(i,i) = a(i,i)*co + a(i,j)*si
       a(j,j) =-a(j,i)*si + a(j,j)*co
       a(i,j) = 0.0
       a(j,i) = 0.0
!
       do k = 1, n
         ww = x(k,i) 
         x(k,i) = ww*co + x(k,j)*si
         x(k,j) =-ww*si + x(k,j)*co
       end do
      end do 
!
      if(i_step .ge. -100) then
       if(ic2 .eq. 0) then
         write(idbg,'(''not convergence at eigen solver of principal value: ''           &
                     ,''  step ='',i10,'' global_elem_ID(ic) ='',i10,'' imax ='',i10 )') &
                      i_step, G_ELM_ID, ic 
       end if
      end if
!!!
      idx_2(1) = 1
      idx_2_iv(1) = 1
      amax2 = a(1,1)

      do i = 2,3
       if(a(i,i) .gt. amax2) then
        idx_2(i) = 1
        idx_2_iv(1) = i
        amax2 = a(i,i)
       end if
      end do

      if(idx_2_iv(1) .eq. 1) then
       if(a(2,2) .gt. a(3,3)) then
        idx_2(2) = 2
        idx_2(3) = 3
        idx_2_iv(2) = 2
        idx_2_iv(3) = 3
       else if(a(2,2) .le. a(3,3)) then
        idx_2(2) = 3
        idx_2(3) = 2
        idx_2_iv(2) = 3
        idx_2_iv(3) = 2
       end if
      else if(idx_2_iv(1) .eq. 2) then
       if(a(3,3) .gt. a(1,1)) then
        idx_2(3) = 2
        idx_2(1) = 3
        idx_2_iv(2) = 3
        idx_2_iv(3) = 1
       else if(a(3,3) .le. a(1,1)) then
        idx_2(3) = 3
        idx_2(1) = 2
        idx_2_iv(2) = 1
        idx_2_iv(3) = 3
       end if
      else if(idx_2_iv(1) .eq. 3) then
       if(a(1,1) .gt. a(2,2)) then
        idx_2(1) = 2
        idx_2(2) = 3
        idx_2_iv(2) = 1
        idx_2_iv(3) = 2
       else if(a(1,1) .le. a(2,2)) then
        idx_2(1) = 3
        idx_2(2) = 2
        idx_2_iv(2) = 2
        idx_2_iv(3) = 1
       end if
      end if

      do i = 1,3
       ARRAY_PRIS_2( i)  = a(idx_2_iv(i),idx_2_iv(i)) 
      end do
      do i = 1,3
       ARRAY_PRIS_2(3+i) = x(i,idx_2_iv(1)) 
       ARRAY_PRIS_2(6+i) = x(i,idx_2_iv(2)) 
       ARRAY_PRIS_2(9+i) = x(i,idx_2_iv(3)) 
      end do

!!     do i = 1,3
!       ARRAY_PRIS_2( i)  = a(i,i) 
!      end do
!      do i = 1,3
!       ARRAY_PRIS_2(3+i) = x(i,1) 
!       ARRAY_PRIS_2(6+i) = x(i,2) 
!       ARRAY_PRIS_2(9+i) = x(i,3) 
!!     end do

!
      end subroutine fstr_elpl_eig_qr 

!C================================================================C
!C-- subroutine fstrNLGEOM_ReadTemp 
!C================================================================C
   subroutine fstrNLGEOM_ReadTemp(hecMESH,hecMAT,fstrSOLID,fstrPARAM    &
                                 ,heat_sta_IW1,heat_sta_IW2             &
                                 ,i_step, DT0)
      use m_fstr

      implicit none

      type (hecmwST_local_mesh) :: hecMESH
      type (hecmwST_matrix    ) :: hecMAT
      type (fstr_solid        ) :: fstrSOLID
      type (fstr_param        ) :: fstrPARAM

! Local variable
      integer(kind=kint) :: i,j,inod,ielem,g1,g2,g3
      integer(kind=kint) :: heat_sta_IW1, heat_sta_IW2
      integer(kind=kint) :: i_step, incr 
      real(kind=kreal)   :: bb(27)
      real(kind=kreal)   :: DT0, CTIME 
      character(len= 30)  :: cdummy,cdummy2(6)
      character(len=100)  :: cdummy3
      character(len=4)  :: cdummy4(2)
      character(len=8)  :: cdummy5(2)
      character(len=1)  :: cdummy6(2)
!     character(len=*)  :: cdummy7(6)
      character  :: cdummy7(6)
!C
!C-- TEMPERATURE
!C
!         read (heat_sta_IW1,*)  cdummy2(1)
!!        read (heat_sta_IW1  )  cdummy2(1)
          read (heat_sta_IW1  )  cdummy4(1)
!         write(idbg        ,*)  cdummy4(1)
!         read (heat_sta_IW1,*) (cdummy2(j),j=1,6)
!         read (heat_sta_IW1  ) (cdummy2(j),j=1,6)
          read (heat_sta_IW1  )  cdummy5(1),incr,cdummy5(2),CTIME
!         write(idbg        ,*)  cdummy5(1),incr,cdummy5(2),CTIME
!         write(idbg        ,*) (cdummy2(j),j=1,6)
        do i= 1, hecMESH%nn_internal
          inod=fstrPARAM%global_local_id(1,i)
!         read (heat_sta_IW1,*) inod,fstrSOLID%TEMPC(i)
          read (heat_sta_IW1  ) inod,fstrSOLID%TEMPC(i)
!         write(idbg        ,*) inod,fstrSOLID%TEMPC(i)
        enddo
!C
!C-- CURE
!C
!         read (heat_sta_IW2,*)  cdummy2(1)
!         read (heat_sta_IW2  )  cdummy2(1)
          read (heat_sta_IW2  )  cdummy4(1)
!         write(idbg        ,*)  cdummy4(1)
!         read (heat_sta_IW2,*) (cdummy2(j),j=1,6)
!         read (heat_sta_IW2  ) (cdummy2(j),j=1,6)
          read (heat_sta_IW2  )  cdummy5(1),incr,cdummy5(2),CTIME
!         write(idbg        ,*)  cdummy5(1),incr,cdummy5(2),CTIME
!         write(idbg        ,*) (cdummy2(j),j=1,6)
        do i= 1, hecMESH%nn_internal
          inod=fstrPARAM%global_local_id(1,i)
!         read (heat_sta_IW2,*) inod,fstrSOLID%CURE3D_nod(i)
          read (heat_sta_IW2  ) inod,fstrSOLID%CURE3D_nod(i)
!         write(idbg        ,*) inod,fstrSOLID%CURE3D_nod(i)
        enddo

   end subroutine fstrNLGEOM_ReadTemp 

end module m_static_NLGEOM_elpl
