!======================================================================!
!                                                                      !
!    Software Name : FrontCOMP_cure   Ver. 3.1                         !
!                                                                      !
!      Module Name : m_static_LIB_3d2TL                                !
!      Category    : Static Analysis                                   !
!                                                                      !
!      Developed based on "FrontSTR" of RSS21 project                  !
!                                                                      !
!                     Written by Toshio Nagashima,   2006/06/01        !
!                                Yasuji Fukahori,    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_LIB_3d2TL
!C================================================================C

!***********************************************************************
!  3D quadratic Element:
!
!  STF_C3D20TL(ILOG,XX,YY,ZZ,DISP,EE,PP,STR,SS,QF)
!  STF_C3D15TL(ILOG,XX,YY,ZZ,DISP,EE,PP,STR,SS,QF)
!  STF_C3D10TL(ILOG,XX,YY,ZZ,DISP,EE,PP,STR,SS,QF    &
!             ,idx_elpl,h_dash,FLG_ELPL,STR_YLD      &
!             ,TT,CURE3D,ccc_rl,ccc_rl1,ccc_rl2      &
!             ,yeta,yeta1,yeta2,idx_visco,DT0        &
!             ,STRAIN_V1,STRAIN_V2,count)
!  UPDATE_C3D20TL(ILOG,XX,YY,ZZ,DISP,DU,EE,PP,STRESS)
!  UPDATE_C3D15TL(ILOG,XX,YY,ZZ,DISP,DU,EE,PP,STRESS)
!  UPDATE_C3D10TL(ILOG,XX,YY,ZZ,DISP,DU,EE,PP,STRESS                     &
!                ,idx_elpl,h_dash,FLG_ELPL,STR_YLD ,STRAIN,EP_BAR ,iter  &
!                ,TT,TTI,CURE3D,CURE3D_I,ccc_rl,ramd_chm                 &
!                ,ccc_rl1,ccc_rl2                                        &
!                ,yeta,yeta1,yeta2,idx_visco,DT0                         &
!                ,STRAIN_V1,STRAIN_V2,count)
!  NRCV_C3D20TL(ILOG,XX,YY,ZZ,TT,T0,EE,PP,ALP,DISP,STR,ARRAY)
!  NRCV_C3D15TL(ILOG,XX,YY,ZZ,TT,T0,EE,PP,ALP,DISP,STR,ARRAY)
!  NRCV_C3D10TL(ILOG,XX,YY,ZZ,TT,T0,EE,PP,ALP,DISP,STRESS,ARRAY,  &
!               STRAIN)
!  TLOAD_C3D10TL(ILOG,XX,YY,ZZ,DISP,DU,EE,PP,STRESS                     &
!               ,idx_elpl,h_dash,FLG_ELPL,STR_YLD ,STRAIN,EP_BAR ,iter  &
!               ,TT,TTI,CURE3D,CURE3D_I,ccc_rl,ramd_chm,VECT)
!***********************************************************************

module m_static_LIB_3d2TL

use m_static_NLGEOM_elpl
use m_utilities



contains


!----------------------------------------------------------------------*
      SUBROUTINE STF_C3D20TL(ILOG,XX,YY,ZZ,DISP,EE,PP,STR,SS,QF   &
                            ,idx_elpl,h_dash,FLG_ELPL,STR_YLD)
!----------------------------------------------------------------------*
!
! CALCULATION 3D 20 NODE PLANE ELEMENT
!
      use hecmw
      use gauss_integration 
      IMPLICIT NONE
! GAUSS POINT
      INTEGER(kind=kint) NG
      PARAMETER(NG=3)
! I/F VARIABLES
      INTEGER(kind=kint) ILOG
      INTEGER(kind=kint) idx_elpl,FLG_ELPL(NG,NG,NG)
      REAL(kind=kreal) h_dash,STR_YLD(NG,NG,NG)
      REAL(kind=kreal) XX(*),YY(*),ZZ(*),DISP(*),EE,PP,STR(6,NG,NG,NG),SS(*),QF(*)
! LOCAL VARIABLES
      INTEGER(kind=kint) NN,NDOF,ISIZE
      PARAMETER(NN=20,NDOF=3,ISIZE=(NN*NDOF)*(NN*NDOF+1)/2)
      INTEGER(kind=kint) I,J,K,NUM
      INTEGER(kind=kint) LX,LY,LZ
      INTEGER(kind=kint) HEAD
      REAL(kind=kreal) RI,SI,TI,RP,SP,TP,RM,SM,TM
      REAL(kind=kreal) H(NN),HR(NN),HS(NN),HT(NN)
      REAL(kind=kreal) XJ11,XJ21,XJ31,XJ12,XJ22,XJ32,XJ13,XJ23,XJ33,DET,WG,DUM
      REAL(kind=kreal) XJI11,XJI21,XJI31,XJI12,XJI22,XJI32,XJI13,XJI23,XJI33
      REAL(kind=kreal) DMAT(6,6),SMAT(9,9)
      REAL(kind=kreal) DNDX,DNDY,DNDZ
      REAL(kind=kreal) DUDX,DUDY,DUDZ
      REAL(kind=kreal) DVDX,DVDY,DVDZ
      REAL(kind=kreal) DWDX,DWDY,DWDZ
      REAL(kind=kreal) BL0(6,NN*NDOF),BL1(6,NN*NDOF),BNL(9,NN*NDOF)
      REAL(kind=kreal) DB(9,NN*NDOF)
      REAL(kind=kreal) UX(NN),UY(NN),UZ(NN)
      REAL(kind=kreal) S11,S22,S33,S13,S23,S12

      INTEGER(kind=kint) MAT_ISO 
      REAL(kind=kreal) BIJ(3,3),MAT_OTHO1(9) 
!* zero clear stiffness
      DO I=1,ISIZE
        SS(I)=0.0
      ENDDO
!* zero clear internal force
      DO I=1,NN*NDOF
        QF(I)=0.0
      ENDDO
!* Total displacement 
      DO I=1,NN
        UX(I)=DISP(NDOF*I-2)
        UY(I)=DISP(NDOF*I-1)
        UZ(I)=DISP(NDOF*I  )
      ENDDO
! SET DMAT
      DO J=1,6
        DO I=1,6
          DMAT(I,J)=0.0
        ENDDO
      ENDDO 
      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
! LOOP FOR INTEGRATION POINTS      
      DO LX=1,NG
        RI=XG(NG,LX)
        DO LY=1,NG
          SI=XG(NG,LY)
          DO LZ=1,NG
            TI=XG(NG,LZ)
            RP=1.0+RI
            SP=1.0+SI
            TP=1.0+TI
            RM=1.0-RI
            SM=1.0-SI
            TM=1.0-TI
!  INTERPOLATION FUNCTION
            H(1)=-0.125*RM*SM*TM*(2.0+RI+SI+TI)
            H(2)=-0.125*RP*SM*TM*(2.0-RI+SI+TI)
            H(3)=-0.125*RP*SP*TM*(2.0-RI-SI+TI)
            H(4)=-0.125*RM*SP*TM*(2.0+RI-SI+TI)
            H(5)=-0.125*RM*SM*TP*(2.0+RI+SI-TI)
            H(6)=-0.125*RP*SM*TP*(2.0-RI+SI-TI)
            H(7)=-0.125*RP*SP*TP*(2.0-RI-SI-TI)
            H(8)=-0.125*RM*SP*TP*(2.0+RI-SI-TI)
            H(9)=0.25*(1.0-RI**2)*SM*TM
            H(10)=0.25*RP*(1.0-SI**2)*TM
            H(11)=0.25*(1.0-RI**2)*SP*TM
            H(12)=0.25*RM*(1.0-SI**2)*TM
            H(13)=0.25*(1.0-RI**2)*SM*TP
            H(14)=0.25*RP*(1.0-SI**2)*TP
            H(15)=0.25*(1.0-RI**2)*SP*TP
            H(16)=0.25*RM*(1.0-SI**2)*TP
            H(17)=0.25*RM*SM*(1.0-TI**2)
            H(18)=0.25*RP*SM*(1.0-TI**2)
            H(19)=0.25*RP*SP*(1.0-TI**2)
            H(20)=0.25*RM*SP*(1.0-TI**2)
!  DERIVATIVE OF INTERPOLATION FUNCTION
!  FOR R-COORDINATE
            HR(1)=-0.125*RM*SM*TM+0.125*SM*TM*(2.0+RI+SI+TI)
            HR(2)=+0.125*RP*SM*TM-0.125*SM*TM*(2.0-RI+SI+TI)
            HR(3)=+0.125*RP*SP*TM-0.125*SP*TM*(2.0-RI-SI+TI)
            HR(4)=-0.125*RM*SP*TM+0.125*SP*TM*(2.0+RI-SI+TI)
            HR(5)=-0.125*RM*SM*TP+0.125*SM*TP*(2.0+RI+SI-TI)
            HR(6)=+0.125*RP*SM*TP-0.125*SM*TP*(2.0-RI+SI-TI)
            HR(7)=+0.125*RP*SP*TP-0.125*SP*TP*(2.0-RI-SI-TI)
            HR(8)=-0.125*RM*SP*TP+0.125*SP*TP*(2.0+RI-SI-TI)
            HR(9 )=-0.50*RI*SM*TM
            HR(10)=+0.25*(1.0-SI**2)*TM
            HR(11)=-0.50*RI*SP*TM
            HR(12)=-0.25*(1.0-SI**2)*TM
            HR(13)=-0.50*RI*SM*TP
            HR(14)=+0.25*(1.0-SI**2)*TP
            HR(15)=-0.50*RI*SP*TP
            HR(16)=-0.25*(1.0-SI**2)*TP
            HR(17)=-0.25*SM*(1.0-TI**2)
            HR(18)=+0.25*SM*(1.0-TI**2)
            HR(19)=+0.25*SP*(1.0-TI**2)
            HR(20)=-0.25*SP*(1.0-TI**2)
!  FOR S-COORDINATE
            HS(1)=-0.125*RM*SM*TM+0.125*RM*TM*(2.0+RI+SI+TI)
            HS(2)=-0.125*RP*SM*TM+0.125*RP*TM*(2.0-RI+SI+TI)
            HS(3)=+0.125*RP*SP*TM-0.125*RP*TM*(2.0-RI-SI+TI)
            HS(4)=+0.125*RM*SP*TM-0.125*RM*TM*(2.0+RI-SI+TI)
            HS(5)=-0.125*RM*SM*TP+0.125*RM*TP*(2.0+RI+SI-TI)
            HS(6)=-0.125*RP*SM*TP+0.125*RP*TP*(2.0-RI+SI-TI)
            HS(7)=+0.125*RP*SP*TP-0.125*RP*TP*(2.0-RI-SI-TI)
            HS(8)=+0.125*RM*SP*TP-0.125*RM*TP*(2.0+RI-SI-TI)
            HS(9)=-0.25*(1.0-RI**2)*TM
            HS(10)=-0.50*RP*SI*TM
            HS(11)=+0.25*(1.0-RI**2)*TM
            HS(12)=-0.50*RM*SI*TM
            HS(13)=-0.25*(1.0-RI**2)*TP
            HS(14)=-0.50*RP*SI*TP
            HS(15)=+0.25*(1.0-RI**2)*TP
            HS(16)=-0.50*RM*SI*TP
            HS(17)=-0.25*RM*(1.0-TI**2)
            HS(18)=-0.25*RP*(1.0-TI**2)
            HS(19)=+0.25*RP*(1.0-TI**2)
            HS(20)=+0.25*RM*(1.0-TI**2)
!  FOR T-COORDINATE
            HT(1)=-0.125*RM*SM*TM+0.125*RM*SM*(2.0+RI+SI+TI)
            HT(2)=-0.125*RP*SM*TM+0.125*RP*SM*(2.0-RI+SI+TI)
            HT(3)=-0.125*RP*SP*TM+0.125*RP*SP*(2.0-RI-SI+TI)
            HT(4)=-0.125*RM*SP*TM+0.125*RM*SP*(2.0+RI-SI+TI)
            HT(5)=+0.125*RM*SM*TP-0.125*RM*SM*(2.0+RI+SI-TI)
            HT(6)=+0.125*RP*SM*TP-0.125*RP*SM*(2.0-RI+SI-TI)
            HT(7)=+0.125*RP*SP*TP-0.125*RP*SP*(2.0-RI-SI-TI)
            HT(8)=+0.125*RM*SP*TP-0.125*RM*SP*(2.0+RI-SI-TI)
            HT(9)=-0.25*(1.0-RI**2)*SM
            HT(10)=-0.25*RP*(1.0-SI**2)
            HT(11)=-0.25*(1.0-RI**2)*SP
            HT(12)=-0.25*RM*(1.0-SI**2)
            HT(13)=0.25*(1.0-RI**2)*SM
            HT(14)=0.25*RP*(1.0-SI**2)
            HT(15)=0.25*(1.0-RI**2)*SP
            HT(16)=0.25*RM*(1.0-SI**2)
            HT(17)=-0.5*RM*SM*TI
            HT(18)=-0.5*RP*SM*TI
            HT(19)=-0.5*RP*SP*TI
            HT(20)=-0.5*RM*SP*TI
!  JACOBI MATRIX 
            XJ11=0.0
            XJ21=0.0
            XJ31=0.0
            XJ12=0.0
            XJ22=0.0
            XJ32=0.0
            XJ13=0.0
            XJ23=0.0
            XJ33=0.0
            DO I=1,NN
              XJ11=XJ11+HR(I)*XX(I)
              XJ21=XJ21+HS(I)*XX(I)
              XJ31=XJ31+HT(I)*XX(I)
              XJ12=XJ12+HR(I)*YY(I)
              XJ22=XJ22+HS(I)*YY(I)
              XJ32=XJ32+HT(I)*YY(I)
              XJ13=XJ13+HR(I)*ZZ(I)
              XJ23=XJ23+HS(I)*ZZ(I)
              XJ33=XJ33+HT(I)*ZZ(I)
            ENDDO
!DETERMINANT OF JACOBIAN
            DET=XJ11*XJ22*XJ33   &
               +XJ12*XJ23*XJ31   &
               +XJ13*XJ21*XJ32   &
               -XJ13*XJ22*XJ31   &
               -XJ12*XJ21*XJ33   &
               -XJ11*XJ23*XJ32
! INVERSION OF JACOBIAN
            DUM=1.0/DET
            XJI11=DUM*( XJ22*XJ33-XJ23*XJ32)
            XJI21=DUM*(-XJ21*XJ33+XJ23*XJ31)
            XJI31=DUM*( XJ21*XJ32-XJ22*XJ31)
            XJI12=DUM*(-XJ12*XJ33+XJ13*XJ32)
            XJI22=DUM*( XJ11*XJ33-XJ13*XJ31)
            XJI32=DUM*(-XJ11*XJ32+XJ12*XJ31)
            XJI13=DUM*( XJ12*XJ23-XJ13*XJ22)
            XJI23=DUM*(-XJ11*XJ23+XJ13*XJ21)
            XJI33=DUM*( XJ11*XJ22-XJ12*XJ21)
!  WEIGHT VALUE AT GAUSSIAN POINT
            WG=WGT(NG,LX)*WGT(NG,LY)*WGT(NG,LZ)*DET
!* DUDX DUDY DUDZ ....
            DUDX=0.0
            DUDY=0.0
            DUDZ=0.0
            DVDX=0.0
            DVDY=0.0
            DVDZ=0.0
            DWDX=0.0
            DWDY=0.0
            DWDZ=0.0
            DO J=1,NN
              DNDX=XJI11*HR(J)+XJI12*HS(J)+XJI13*HT(J)
              DNDY=XJI21*HR(J)+XJI22*HS(J)+XJI23*HT(J)
              DNDZ=XJI31*HR(J)+XJI32*HS(J)+XJI33*HT(J)
              DUDX=DUDX+DNDX*UX(J)
              DUDY=DUDY+DNDY*UX(J)
              DUDZ=DUDZ+DNDZ*UX(J)
              DVDX=DVDX+DNDX*UY(J)
              DVDY=DVDY+DNDY*UY(J)
              DVDZ=DVDZ+DNDZ*UY(J)
              DWDX=DWDX+DNDX*UZ(J)
              DWDY=DWDY+DNDY*UZ(J)
              DWDZ=DWDZ+DNDZ*UZ(J)
            ENDDO
!** Fill BMAT
            HEAD=1
            DO J=1,NN
!* DNDX DNDY DNDZ
              DNDX=XJI11*HR(J)+XJI12*HS(J)+XJI13*HT(J)
              DNDY=XJI21*HR(J)+XJI22*HS(J)+XJI23*HT(J)
              DNDZ=XJI31*HR(J)+XJI32*HS(J)+XJI33*HT(J)
!** BL0
              BL0(1,HEAD )=DNDX
              BL0(2,HEAD )=0.0
              BL0(3,HEAD )=0.0
              BL0(4,HEAD )=DNDY
              BL0(5,HEAD )=0.0
              BL0(6,HEAD )=DNDZ
              BL0(1,HEAD+1)=0.0
              BL0(2,HEAD+1)=DNDY
              BL0(3,HEAD+1)=0.0
              BL0(4,HEAD+1)=DNDX
              BL0(5,HEAD+1)=DNDZ
              BL0(6,HEAD+1)=0.0
              BL0(1,HEAD+2)=0.0
              BL0(2,HEAD+2)=0.0
              BL0(3,HEAD+2)=DNDZ
              BL0(4,HEAD+2)=0.0
              BL0(5,HEAD+2)=DNDY
              BL0(6,HEAD+2)=DNDX
!** BL1
              BL1(1,HEAD )=DUDX*DNDX
              BL1(2,HEAD )=DUDY*DNDY
              BL1(3,HEAD )=DUDZ*DNDZ
              BL1(4,HEAD )=DUDY*DNDX+DUDX*DNDY
              BL1(5,HEAD )=DUDY*DNDZ+DUDZ*DNDY
              BL1(6,HEAD )=DUDZ*DNDX+DUDX*DNDZ
              BL1(1,HEAD+1)=DVDX*DNDX
              BL1(2,HEAD+1)=DVDY*DNDY
              BL1(3,HEAD+1)=DVDZ*DNDZ
              BL1(4,HEAD+1)=DVDY*DNDX+DVDX*DNDY
              BL1(5,HEAD+1)=DVDY*DNDZ+DVDZ*DNDY
              BL1(6,HEAD+1)=DVDZ*DNDX+DVDX*DNDZ
              BL1(1,HEAD+2)=DWDX*DNDX
              BL1(2,HEAD+2)=DWDY*DNDY
              BL1(3,HEAD+2)=DWDZ*DNDZ
              BL1(4,HEAD+2)=DWDY*DNDX+DWDX*DNDY
              BL1(5,HEAD+2)=DWDY*DNDZ+DWDZ*DNDY
              BL1(6,HEAD+2)=DWDZ*DNDX+DWDX*DNDZ
!** BNL
              BNL(1,HEAD )=DNDX
              BNL(2,HEAD )=0.0
              BNL(3,HEAD )=0.0
              BNL(4,HEAD )=DNDY
              BNL(5,HEAD )=0.0
              BNL(6,HEAD )=0.0
              BNL(7,HEAD )=DNDZ
              BNL(8,HEAD )=0.0
              BNL(9,HEAD )=0.0
              BNL(1,HEAD+1)=0.0
              BNL(2,HEAD+1)=DNDX
              BNL(3,HEAD+1)=0.0
              BNL(4,HEAD+1)=0.0
              BNL(5,HEAD+1)=DNDY
              BNL(6,HEAD+1)=0.0
              BNL(7,HEAD+1)=0.0
              BNL(8,HEAD+1)=DNDZ
              BNL(9,HEAD+1)=0.0
              BNL(1,HEAD+2)=0.0
              BNL(2,HEAD+2)=0.0
              BNL(3,HEAD+2)=DNDX
              BNL(4,HEAD+2)=0.0
              BNL(5,HEAD+2)=0.0
              BNL(6,HEAD+2)=DNDY
              BNL(7,HEAD+2)=0.0
              BNL(8,HEAD+2)=0.0
              BNL(9,HEAD+2)=DNDZ
              HEAD=HEAD+NDOF
            ENDDO
!** Linear stiffness
!ogasawara till next ogasawara
!**   elasto-plastic matrix
      if(idx_elpl .ge. 1) then
        call fstr_elpl_DMAT(DMAT,STR(:,LX,LY,LZ),FLG_ELPL(LX,LY,LZ), h_dash, EE, PP   &
                         ,BIJ,MAT_OTHO1,MAT_ISO)
      end if
!ogasawara till next ogasawara

            DO J=1,6
              DO K=1,NN*NDOF
                DB(J,K)=0.0
              ENDDO
            ENDDO
            DO J=1,6
              DO K=1,NN*NDOF
                DB(J,K)=DB(J,K)                  &
                +DMAT(J,1)*(BL0(1,K)+BL1(1,K))   &
                +DMAT(J,2)*(BL0(2,K)+BL1(2,K))   &
                +DMAT(J,3)*(BL0(3,K)+BL1(3,K))   &
                +DMAT(J,4)*(BL0(4,K)+BL1(4,K))   &
                +DMAT(J,5)*(BL0(5,K)+BL1(5,K))   &
                +DMAT(J,6)*(BL0(6,K)+BL1(6,K))
              ENDDO
            ENDDO
            DO J=1,NN*NDOF
              DO K=1,J
                NUM=(J-1)*J/2+K
                SS(NUM)=SS(NUM)                       &
                +( (BL0(1,K)+BL1(1,K))*DB(1,J)        &
                +  (BL0(2,K)+BL1(2,K))*DB(2,J)        &
                +  (BL0(3,K)+BL1(3,K))*DB(3,J)        &
                +  (BL0(4,K)+BL1(4,K))*DB(4,J)        &
                +  (BL0(5,K)+BL1(5,K))*DB(5,J)        &
                +  (BL0(6,K)+BL1(6,K))*DB(6,J) )*WG
              ENDDO
            ENDDO
!** NonLinear stiffness
!* Current stress
         if(idx_elpl .lt. 2) then
!
            S11=STR(1,LX,LY,LZ)
            S22=STR(2,LX,LY,LZ)
            S33=STR(3,LX,LY,LZ)
            S12=STR(4,LX,LY,LZ)
            S23=STR(5,LX,LY,LZ)
            S13=STR(6,LX,LY,LZ)
            DO J=1,9
              DO K=1,9
                SMAT(J,K)=0.0
              ENDDO
            ENDDO
            DO J=1,3
              SMAT(J  ,J  )=S11
              SMAT(J+3,J+3)=S22
              SMAT(J+6,J+6)=S33
              SMAT(J+3,J+6)=S23
              SMAT(J  ,J+3)=S12
              SMAT(J  ,J+6)=S13
              SMAT(J+6,J+3)=S23
              SMAT(J+3,J  )=S12
              SMAT(J+6,J  )=S13
            ENDDO
            DO J=1,9
              DO K=1,NN*NDOF
                DB(J,K)=0.0
              ENDDO
            ENDDO
            DO J=1,9
              DO K=1,NN*NDOF
                DB(J,K)=DB(J,K)         &
                +SMAT(J,1)*BNL(1,K)     &
                +SMAT(J,2)*BNL(2,K)     &
                +SMAT(J,3)*BNL(3,K)     &
                +SMAT(J,4)*BNL(4,K)     &
                +SMAT(J,5)*BNL(5,K)     &
                +SMAT(J,6)*BNL(6,K)     &
                +SMAT(J,7)*BNL(7,K)     &
                +SMAT(J,8)*BNL(8,K)     &
                +SMAT(J,9)*BNL(9,K)
              ENDDO
            ENDDO
            DO J=1,NN*NDOF
              DO K=1,J
                NUM=(J-1)*J/2+K
                SS(NUM)=SS(NUM)+        &
                (BNL(1,K)*DB(1,J)       &
                +BNL(2,K)*DB(2,J)       &
                +BNL(3,K)*DB(3,J)       &
                +BNL(4,K)*DB(4,J)       &
                +BNL(5,K)*DB(5,J)       &
                +BNL(6,K)*DB(6,J)       &
                +BNL(7,K)*DB(7,J)       &
                +BNL(8,K)*DB(8,J)       &
                +BNL(9,K)*DB(9,J) )*WG
              ENDDO
            ENDDO
!
         end if
!* Add Q force
            DO K=1,NN*NDOF
              QF(K)=QF(K)                      &
              +(S11*(BL0(1,K)+BL1(1,K))        &
              + S22*(BL0(2,K)+BL1(2,K))        &
              + S33*(BL0(3,K)+BL1(3,K))        &
              + S12*(BL0(4,K)+BL1(4,K))        &
              + S23*(BL0(5,K)+BL1(5,K))        &
              + S13*(BL0(6,K)+BL1(6,K)) )*WG
            ENDDO
          ENDDO
        ENDDO
      ENDDO
      RETURN
      END SUBROUTINE STF_C3D20TL
!----------------------------------------------------------------------*
      SUBROUTINE STF_C3D15TL(ILOG,XX,YY,ZZ,DISP,EE,PP,STR,SS,QF   &
                            ,idx_elpl,h_dash,FLG_ELPL,STR_YLD)
!----------------------------------------------------------------------*
!
! CALCULATION 3D 15 NODE PLANE ELEMENT
!
      use hecmw
      use gauss_integration 
      IMPLICIT NONE
! GAUSS POINT
      INTEGER(kind=kint) NG
      PARAMETER(NG=3)
! I/F VARIABLES
      INTEGER(kind=kint) ILOG
      INTEGER(kind=kint) idx_elpl,FLG_ELPL(NG,NG,NG)
      REAL(kind=kreal) h_dash,STR_YLD(NG,NG,NG)
      REAL(kind=kreal) XX(*),YY(*),ZZ(*),DISP(*),EE,PP,STR(6,NG,NG,NG),SS(*),QF(*)
! LOCAL VARIABLES
      INTEGER(kind=kint) NN,NDOF,ISIZE
      PARAMETER(NN=15,NDOF=3,ISIZE=(NN*NDOF)*(NN*NDOF+1)/2)
      INTEGER(kind=kint) I,J,K,NUM
      INTEGER(kind=kint) L1,L2,LZ
      INTEGER(kind=kint) HEAD
      REAL(kind=kreal) RI,SI,TI,RP,SP,TP,RM,SM,TM
      REAL(kind=kreal) H(NN),HL1(NN),HL2(NN),HL3(NN),HZ(NN)
      REAL(kind=kreal) X1,X2,X3,XL1,XL2,ZI
      REAL(kind=kreal) XJ11,XJ21,XJ31,XJ12,XJ22,XJ32,XJ13,XJ23,XJ33,DET,WG,DUM
      REAL(kind=kreal) XJI11,XJI21,XJI31,XJI12,XJI22,XJI32,XJI13,XJI23,XJI33
      REAL(kind=kreal) DMAT(6,6),SMAT(9,9)
      REAL(kind=kreal) DNDX,DNDY,DNDZ
      REAL(kind=kreal) DUDX,DUDY,DUDZ
      REAL(kind=kreal) DVDX,DVDY,DVDZ
      REAL(kind=kreal) DWDX,DWDY,DWDZ
      REAL(kind=kreal) BL0(6,NN*NDOF),BL1(6,NN*NDOF),BNL(9,NN*NDOF)
      REAL(kind=kreal) DB(9,NN*NDOF)
      REAL(kind=kreal) UX(NN),UY(NN),UZ(NN)
      REAL(kind=kreal) S11,S22,S33,S13,S23,S12

      INTEGER(kind=kint) MAT_ISO 
      REAL(kind=kreal) BIJ(3,3),MAT_OTHO1(9) 
!* zero clear stiffness
      DO I=1,ISIZE
        SS(I)=0.0
      ENDDO
!* zero clear internal force
      DO I=1,NN*NDOF
        QF(I)=0.0
      ENDDO
!* Total displacement 
      DO I=1,NN
        UX(I)=DISP(NDOF*I-2)
        UY(I)=DISP(NDOF*I-1)
        UZ(I)=DISP(NDOF*I  )
      ENDDO
! SET DMAT
      DO J=1,6
        DO I=1,6
          DMAT(I,J)=0.0
        ENDDO
      ENDDO 
      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
! LOOP FOR INTEGRATION POINTS      
      DO LZ=1,NG
        ZI=XG(NG,LZ)
        DO L2=1,NG
          XL2=XG(NG,L2)
          X2 =(XL2+1.0)*0.5
          DO L1=1,NG
            XL1=XG(NG,L1)
            X1=0.5*(1.0-X2)*(XL1+1.0)
!  INTERPOLATION FUNCTION
            X3=1.0-X1-X2
            H(1)= 0.5*X1*(2.0*X1-2.-ZI)*(1.0-ZI)
            H(2)= 0.5*X2*(2.0*X2-2.-ZI)*(1.0-ZI)
            H(3)= 0.5*X3*(2.0*X3-2.-ZI)*(1.0-ZI)
            H(4)= 0.5*X1*(2.0*X1-2.+ZI)*(1.0+ZI)
            H(5)= 0.5*X2*(2.0*X2-2.+ZI)*(1.0+ZI)
            H(6)= 0.5*X3*(2.0*X3-2.+ZI)*(1.0+ZI)
            H(7)= 2.0*X1*X2*(1.0-ZI)
            H(8)= 2.0*X2*X3*(1.0-ZI)
            H(9)= 2.0*X1*X3*(1.0-ZI)
            H(10)=2.0*X1*X2*(1.0+ZI)
            H(11)=2.0*X2*X3*(1.0+ZI)
            H(12)=2.0*X1*X3*(1.0+ZI)
            H(13)=X1*(1.0-ZI**2)
            H(14)=X2*(1.0-ZI**2)
            H(15)=X3*(1.0-ZI**2)
!  DERIVATIVE OF INTERPOLATION FUNCTION
!  FOR L1-COORDINATE
            HL1(1)= 0.5*(4.0*X1-2.-ZI)*(1.0-ZI)
            HL1(2)= 0.0
            HL1(3)= 0.0
            HL1(4)= 0.5*(4.0*X1-2.+ZI)*(1.0+ZI)
            HL1(5)= 0.0
            HL1(6)= 0.0
            HL1(7)= 2.0*X2*(1.0-ZI)
            HL1(8)= 0.0
            HL1(9)= 2.0*X3*(1.0-ZI)
            HL1(10)=2.0*X2*(1.0+ZI)
            HL1(11)=0.0
            HL1(12)=2.0*X3*(1.0+ZI)
            HL1(13)=1.0-ZI**2
            HL1(14)=0.0
            HL1(15)=0.0
!  FOR L2-COORDINATE
            HL2(1)= 0.0
            HL2(2)= 0.5*(4.0*X2-2.-ZI)*(1.0-ZI)
            HL2(3)= 0.0
            HL2(4)= 0.0
            HL2(5)= 0.5*(4.0*X2-2.+ZI)*(1.0+ZI)
            HL2(6)= 0.0
            HL2(7)= 2.0*X1*(1.0-ZI)
            HL2(8)= 2.0*X3*(1.0-ZI)
            HL2(9)= 0.0
            HL2(10)=2.0*X1*(1.0+ZI)
            HL2(11)=2.0*X3*(1.0+ZI)
            HL2(12)=0.0
            HL2(13)=0.0
            HL2(14)=1.0-ZI**2
            HL2(15)=0.0
!  FOR L3-COORDINATE
            HL3(1)= 0.0
            HL3(2)= 0.0
            HL3(3)= 0.5*(4.0*X3-2.-ZI)*(1.0-ZI)
            HL3(4)= 0.0
            HL3(5)= 0.0
            HL3(6)= 0.5*(4.0*X3-2.+ZI)*(1.0+ZI)
            HL3(7)= 0.0
            HL3(8)= 2.0*X2*(1.0-ZI)
            HL3(9)= 2.0*X1*(1.0-ZI)
            HL3(10)=0.0
            HL3(11)=2.0*X2*(1.0+ZI)
            HL3(12)=2.0*X1*(1.0+ZI)
            HL3(13)=0.0
            HL3(14)=0.0
            HL3(15)=1.0-ZI**2
!  FOR Z-COORDINATE
            HZ(1)= 0.5*X1*(-2.0*X1+1.0+2.0*ZI)
            HZ(2)= 0.5*X2*(-2.0*X2+1.0+2.0*ZI)
            HZ(3)= 0.5*X3*(-2.0*X3+1.0+2.0*ZI)
            HZ(4)= 0.5*X1*( 2.0*X1-1.0+2.0*ZI)
            HZ(5)= 0.5*X2*( 2.0*X2-1.0+2.0*ZI)
            HZ(6)= 0.5*X3*( 2.0*X3-1.0+2.0*ZI)
            HZ(7)= -2.0*X1*X2
            HZ(8)= -2.0*X2*X3
            HZ(9)= -2.0*X1*X3
            HZ(10)= 2.0*X1*X2
            HZ(11)= 2.0*X2*X3
            HZ(12)= 2.0*X1*X3
            HZ(13)=-2.0*X1*ZI
            HZ(14)=-2.0*X2*ZI
            HZ(15)=-2.0*X3*ZI
!  JACOBI MATRIX 
            XJ11=0.0
            XJ21=0.0
            XJ31=0.0
            XJ12=0.0
            XJ22=0.0
            XJ32=0.0
            XJ13=0.0
            XJ23=0.0
            XJ33=0.0
            DO I=1,NN
              XJ11=XJ11+(HL1(I)-HL3(I))*XX(I)
              XJ21=XJ21+(HL2(I)-HL3(I))*XX(I)
              XJ31=XJ31+HZ(I)*XX(I)
              XJ12=XJ12+(HL1(I)-HL3(I))*YY(I)
              XJ22=XJ22+(HL2(I)-HL3(I))*YY(I)
              XJ32=XJ32+HZ(I)*YY(I)
              XJ13=XJ13+(HL1(I)-HL3(I))*ZZ(I)
              XJ23=XJ23+(HL2(I)-HL3(I))*ZZ(I)
              XJ33=XJ33+HZ(I)*ZZ(I)
            ENDDO
!DETERMINANT OF JACOBIAN
            DET=XJ11*XJ22*XJ33   &
              +XJ12*XJ23*XJ31    &
              +XJ13*XJ21*XJ32    &
              -XJ13*XJ22*XJ31    &
              -XJ12*XJ21*XJ33    &
              -XJ11*XJ23*XJ32
! INVERSION OF JACOBIAN
            DUM=1.0/DET
            XJI11=DUM*( XJ22*XJ33-XJ23*XJ32)
            XJI21=DUM*(-XJ21*XJ33+XJ23*XJ31)
            XJI31=DUM*( XJ21*XJ32-XJ22*XJ31)
            XJI12=DUM*(-XJ12*XJ33+XJ13*XJ32)
            XJI22=DUM*( XJ11*XJ33-XJ13*XJ31)
            XJI32=DUM*(-XJ11*XJ32+XJ12*XJ31)
            XJI13=DUM*( XJ12*XJ23-XJ13*XJ22)
            XJI23=DUM*(-XJ11*XJ23+XJ13*XJ21)
            XJI33=DUM*( XJ11*XJ22-XJ12*XJ21)
!  WEIGHT VALUE AT GAUSSIAN POINT
            WG=WGT(NG,L1)*WGT(NG,L2)*WGT(NG,LZ)*DET*(1.0-X2)*0.25
!* DUDX DUDY DUDZ ....
            DUDX=0.0
            DUDY=0.0
            DUDZ=0.0
            DVDX=0.0
            DVDY=0.0
            DVDZ=0.0
            DWDX=0.0
            DWDY=0.0
            DWDZ=0.0
            DO J=1,NN
              DNDX=XJI11*(HL1(J)-HL3(J)) +XJI12*(HL2(J)-HL3(J)) +XJI13*HZ(J)
              DNDY=XJI21*(HL1(J)-HL3(J)) +XJI22*(HL2(J)-HL3(J)) +XJI23*HZ(J)
              DNDZ=XJI31*(HL1(J)-HL3(J)) +XJI32*(HL2(J)-HL3(J)) +XJI33*HZ(J)
              DUDX=DUDX+DNDX*UX(J)
              DUDY=DUDY+DNDY*UX(J)
              DUDZ=DUDZ+DNDZ*UX(J)
              DVDX=DVDX+DNDX*UY(J)
              DVDY=DVDY+DNDY*UY(J)
              DVDZ=DVDZ+DNDZ*UY(J)
              DWDX=DWDX+DNDX*UZ(J)
              DWDY=DWDY+DNDY*UZ(J)
              DWDZ=DWDZ+DNDZ*UZ(J)
            ENDDO
!** Fill BMAT
            HEAD=1
            DO J=1,NN
!* DNDX DNDY DNDZ
              DNDX=XJI11*(HL1(J)-HL3(J)) +XJI12*(HL2(J)-HL3(J)) +XJI13*HZ(J)
              DNDY=XJI21*(HL1(J)-HL3(J)) +XJI22*(HL2(J)-HL3(J)) +XJI23*HZ(J)
              DNDZ=XJI31*(HL1(J)-HL3(J)) +XJI32*(HL2(J)-HL3(J)) +XJI33*HZ(J)
!** BL0
              BL0(1,HEAD )=DNDX
              BL0(2,HEAD )=0.0
              BL0(3,HEAD )=0.0
              BL0(4,HEAD )=DNDY
              BL0(5,HEAD )=0.0
              BL0(6,HEAD )=DNDZ
              BL0(1,HEAD+1)=0.0
              BL0(2,HEAD+1)=DNDY
              BL0(3,HEAD+1)=0.0
              BL0(4,HEAD+1)=DNDX
              BL0(5,HEAD+1)=DNDZ
              BL0(6,HEAD+1)=0.0
              BL0(1,HEAD+2)=0.0
              BL0(2,HEAD+2)=0.0
              BL0(3,HEAD+2)=DNDZ
              BL0(4,HEAD+2)=0.0
              BL0(5,HEAD+2)=DNDY
              BL0(6,HEAD+2)=DNDX
!** BL1
              BL1(1,HEAD )=DUDX*DNDX
              BL1(2,HEAD )=DUDY*DNDY
              BL1(3,HEAD )=DUDZ*DNDZ
              BL1(4,HEAD )=DUDY*DNDX+DUDX*DNDY
              BL1(5,HEAD )=DUDY*DNDZ+DUDZ*DNDY
              BL1(6,HEAD )=DUDZ*DNDX+DUDX*DNDZ
              BL1(1,HEAD+1)=DVDX*DNDX
              BL1(2,HEAD+1)=DVDY*DNDY
              BL1(3,HEAD+1)=DVDZ*DNDZ
              BL1(4,HEAD+1)=DVDY*DNDX+DVDX*DNDY
              BL1(5,HEAD+1)=DVDY*DNDZ+DVDZ*DNDY
              BL1(6,HEAD+1)=DVDZ*DNDX+DVDX*DNDZ
              BL1(1,HEAD+2)=DWDX*DNDX
              BL1(2,HEAD+2)=DWDY*DNDY
              BL1(3,HEAD+2)=DWDZ*DNDZ
              BL1(4,HEAD+2)=DWDY*DNDX+DWDX*DNDY
              BL1(5,HEAD+2)=DWDY*DNDZ+DWDZ*DNDY
              BL1(6,HEAD+2)=DWDZ*DNDX+DWDX*DNDZ
!** BNL
              BNL(1,HEAD )=DNDX
              BNL(2,HEAD )=0.0
              BNL(3,HEAD )=0.0
              BNL(4,HEAD )=DNDY
              BNL(5,HEAD )=0.0
              BNL(6,HEAD )=0.0
              BNL(7,HEAD )=DNDZ
              BNL(8,HEAD )=0.0
              BNL(9,HEAD )=0.0
              BNL(1,HEAD+1)=0.0
              BNL(2,HEAD+1)=DNDX
              BNL(3,HEAD+1)=0.0
              BNL(4,HEAD+1)=0.0
              BNL(5,HEAD+1)=DNDY
              BNL(6,HEAD+1)=0.0
              BNL(7,HEAD+1)=0.0
              BNL(8,HEAD+1)=DNDZ
              BNL(9,HEAD+1)=0.0
              BNL(1,HEAD+2)=0.0
              BNL(2,HEAD+2)=0.0
              BNL(3,HEAD+2)=DNDX
              BNL(4,HEAD+2)=0.0
              BNL(5,HEAD+2)=0.0
              BNL(6,HEAD+2)=DNDY
              BNL(7,HEAD+2)=0.0
              BNL(8,HEAD+2)=0.0
              BNL(9,HEAD+2)=DNDZ
                HEAD=HEAD+NDOF
            ENDDO
!** Linear stiffness
!ogasawara till next ogasawara
!**   elasto-plastic matrix
      if(idx_elpl .ge. 1) then
        call fstr_elpl_DMAT(DMAT,STR(:,L1,L2,LZ),FLG_ELPL(L1,L2,LZ), h_dash, EE, PP   &
                         ,BIJ,MAT_OTHO1,MAT_ISO)
      end if
!ogasawara till next ogasawara

            DO J=1,6
              DO K=1,NN*NDOF
                DB(J,K)=0.0
              ENDDO
            ENDDO
            DO J=1,6
              DO K=1,NN*NDOF
                DB(J,K)=DB(J,K)                     &
                +DMAT(J,1)*(BL0(1,K)+BL1(1,K))      &
                +DMAT(J,2)*(BL0(2,K)+BL1(2,K))      &
                +DMAT(J,3)*(BL0(3,K)+BL1(3,K))      &
                +DMAT(J,4)*(BL0(4,K)+BL1(4,K))      &
                +DMAT(J,5)*(BL0(5,K)+BL1(5,K))      &
                +DMAT(J,6)*(BL0(6,K)+BL1(6,K))
              ENDDO
            ENDDO
            DO J=1,NN*NDOF
              DO K=1,J
                NUM=(J-1)*J/2+K
                SS(NUM)=SS(NUM)                        &
                +( (BL0(1,K)+BL1(1,K))*DB(1,J)         &
                +  (BL0(2,K)+BL1(2,K))*DB(2,J)         &
                +  (BL0(3,K)+BL1(3,K))*DB(3,J)         &
                +  (BL0(4,K)+BL1(4,K))*DB(4,J)         &
                +  (BL0(5,K)+BL1(5,K))*DB(5,J)         &
                +  (BL0(6,K)+BL1(6,K))*DB(6,J) )*WG
              ENDDO
            ENDDO
!** NonLinear stiffness
!* Current stress
         if(idx_elpl .lt. 2) then
!
            S11=STR(1,L1,L2,LZ)
            S22=STR(2,L1,L2,LZ)
            S33=STR(3,L1,L2,LZ)
            S12=STR(4,L1,L2,LZ)
            S23=STR(5,L1,L2,LZ)
            S13=STR(6,L1,L2,LZ)      
            DO J=1,9
              DO K=1,9
                SMAT(J,K)=0.0
              ENDDO
            ENDDO
            DO J=1,3
              SMAT(J  ,J  )=S11
              SMAT(J+3,J+3)=S22
              SMAT(J+6,J+6)=S33
              SMAT(J+3,J+6)=S23
              SMAT(J  ,J+3)=S12
              SMAT(J  ,J+6)=S13
              SMAT(J+6,J+3)=S23
              SMAT(J+3,J  )=S12
              SMAT(J+6,J  )=S13
            ENDDO
            DO J=1,9
              DO K=1,NN*NDOF
                DB(J,K)=0.0
              ENDDO
            ENDDO
            DO J=1,9
              DO K=1,NN*NDOF
                DB(J,K)=DB(J,K)         &
                +SMAT(J,1)*BNL(1,K)     &
                +SMAT(J,2)*BNL(2,K)     &
                +SMAT(J,3)*BNL(3,K)     &
                +SMAT(J,4)*BNL(4,K)     &
                +SMAT(J,5)*BNL(5,K)     &
                +SMAT(J,6)*BNL(6,K)     &
                +SMAT(J,7)*BNL(7,K)     &
                +SMAT(J,8)*BNL(8,K)     &
                +SMAT(J,9)*BNL(9,K)
              ENDDO
            ENDDO
            DO J=1,NN*NDOF
              DO K=1,J
                NUM=(J-1)*J/2+K
                SS(NUM)=SS(NUM)+         &
                (BNL(1,K)*DB(1,J)        &
                +BNL(2,K)*DB(2,J)        &
                +BNL(3,K)*DB(3,J)        &
                +BNL(4,K)*DB(4,J)        &
                +BNL(5,K)*DB(5,J)        &
                +BNL(6,K)*DB(6,J)        &
                +BNL(7,K)*DB(7,J)        &
                +BNL(8,K)*DB(8,J)        &
                +BNL(9,K)*DB(9,J) )*WG
              ENDDO
            ENDDO
!
         end if
!* Add Q force
            DO K=1,NN*NDOF
              QF(K)=QF(K)                     &
              +(S11*(BL0(1,K)+BL1(1,K))       &
              + S22*(BL0(2,K)+BL1(2,K))       &
              + S33*(BL0(3,K)+BL1(3,K))       &
              + S12*(BL0(4,K)+BL1(4,K))       &
              + S23*(BL0(5,K)+BL1(5,K))       &
              + S13*(BL0(6,K)+BL1(6,K)) )*WG     
            ENDDO
          ENDDO
        ENDDO
      ENDDO
      RETURN
      END SUBROUTINE STF_C3D15TL
!----------------------------------------------------------------------*
      SUBROUTINE STF_C3D10TL(ILOG,XX,YY,ZZ,DISP,EE,PP,STR,SS,QF,QF_RES    &
                            ,idx_elpl,h_dash,FLG_ELPL,STR_YLD             &
                            ,idx_temp                                     &
                            ,TT,TTI,CURE3D_nod,CURE3D_nod_I,ccc_rl,ramd_chm  &
                            ,ccc_rl1,ccc_rl2                              &
                            ,yeta,yeta1,yeta2,idx_visco,DT0               &
                            ,STRAIN_V1,STRAIN_V2,iter,i_step,count,icel   &
                            ,MID,ALPHA,MAT_ISO,BIJ,ALPHA_B,MAT_OTHO1)
!----------------------------------------------------------------------*
!
! CALCULATION 3D 10 NODE PLANE ELEMENT
!
      use hecmw
      use gauss_integration 
      IMPLICIT NONE
! GAUSS POINT
      INTEGER(kind=kint) NG
      PARAMETER(NG=3)
! I/F VARIABLES
      INTEGER(kind=kint) ILOG,count,icel,MID 
      INTEGER(kind=kint) idx_elpl,FLG_ELPL(NG,NG,NG)
      INTEGER(kind=kint) idx_temp,idx_visco
      INTEGER(kind=kint) iter, i_step
      REAL(kind=kreal) h_dash,STR_YLD(NG,NG,NG)
      REAL(kind=kreal) XX(*),YY(*),ZZ(*),TT(*),TTI(*),DISP(*),EE,PP,STR(6,NG,NG,NG),SS(*),QF(*)
      REAL(kind=kreal) QF_RES(*)
      REAL(kind=kreal) CURE3D_nod(*),CURE3D_nod_I(*),CURE3D_itg,CURE3D_itg_I
      REAL(kind=kreal) ccc_rl(*),ccc_rl1(*),ccc_rl2(*)
      REAL(kind=kreal) yeta(*),yeta1(*),yeta2(*),DT0
      REAL(kind=kreal) STRAIN_V1(6,NG,NG,NG),STRAIN_V2(6,NG,NG,NG) 
! LOCAL VARIABLES
      INTEGER(kind=kint) NN,NDOF,ISIZE
      PARAMETER(NN=10,NDOF=3,ISIZE=(NN*NDOF)*(NN*NDOF+1)/2)
      INTEGER(kind=kint) I,J,K,L,NUM
      INTEGER(kind=kint) L1,L2,L3
      INTEGER(kind=kint) HEAD
      REAL(kind=kreal) RI,SI,TI,RP,SP,TP,RM,SM,TM
      REAL(kind=kreal) H(NN),HL1(NN),HL2(NN),HL3(NN),HL4(NN)
      REAL(kind=kreal) X1,X2,X3,X4,XL1,XL2,XL3
      REAL(kind=kreal) XJ11,XJ21,XJ31,XJ12,XJ22,XJ32,XJ13,XJ23,XJ33,DET,WG,DUM
      REAL(kind=kreal) XJI11,XJI21,XJI31,XJI12,XJI22,XJI32,XJI13,XJI23,XJI33
      REAL(kind=kreal) DMAT (6,6),SMAT(9,9)
      REAL(kind=kreal) DMAT1(6,6),DMAT2(6,6)
        REAL(kind=kreal) DMAT_G(6,6),DEPS_AD(6),DEPS_AM  
        REAL(kind=kreal) C_G,C_K,C_G1,C_G2,C_K1,C_K2,C_GT,C_KT 
      REAL(kind=kreal) YEINV(6,6),YEINV1(6,6),YEINV2(6,6)
      REAL(kind=kreal) DNDX,DNDY,DNDZ
      REAL(kind=kreal) DUDX,DUDY,DUDZ
      REAL(kind=kreal) DVDX,DVDY,DVDZ
      REAL(kind=kreal) DWDX,DWDY,DWDZ
      REAL(kind=kreal) BL0(6,NN*NDOF),BL1(6,NN*NDOF),BNL(9,NN*NDOF)
      REAL(kind=kreal) DB(9,NN*NDOF)
      REAL(kind=kreal) UX(NN),UY(NN),UZ(NN)
      REAL(kind=kreal) S11,S22,S33,S13,S23,S12
        REAL(kind=kreal) SD(6),SMN 
        REAL(kind=kreal) E_V1(6),E_V2(6),ED_V1(6),EMN_V1,ED_V2(6),EMN_V2 
      REAL(kind=kreal) TEMPC,TEMPI,DSTR_A(6),EE1,EE2,PP1,PP2,YETA_G,YETA_K
        REAL(kind=kreal) YETA_G1,YETA_G2,YETA_K1,YETA_K2,T_G,T_G1,T_G2,T_K,T_K1,T_K2
      REAL(kind=kreal) AAA(6,6),BBB(6),BBB1(6),BBB2(6)
      REAL(kind=kreal) DS(6)
      REAL(kind=kreal) ramd_chm,ALPHA,ALP,ALP1,THERMAL_EPS,EPSTH(6)


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

      REAL(kind=kreal) THERMAL_EPS1,THERMAL_EPS2,THERMAL_EPS3 
      REAL(kind=kreal) EPS_G(3,3),EPS_GB(3,3) 

!* zero clear stiffness
      DO I=1,ISIZE
        SS(I)=0.0
      ENDDO
!* zero clear internal force
      DO I=1,NN*NDOF
        QF(I)=0.0
        QF_RES(I)=0.0
      ENDDO

      THERMAL_EPS = 0.0

      DO I=1,6
        EPSTH(I)  = 0.0
        DS(I)     = 0.0
        DSTR_A(I) = 0.0
      END DO
!
!* Total displacement 
      DO I=1,NN
        UX(I)=DISP(NDOF*I-2)
        UY(I)=DISP(NDOF*I-1)
        UZ(I)=DISP(NDOF*I  )
      ENDDO
! LOOP FOR INTEGRATION POINTS      
      DO L3=1,NG
        XL3=XG(NG,L3)
        X3 =(XL3+1.0)*0.5
        DO L2=1,NG
          XL2=XG(NG,L2)
          X2 =(1.0-X3)*(XL2+1.0)*0.5
          DO L1=1,NG
            XL1=XG(NG,L1)
            X1=(1.0-X2-X3)*(XL1+1.0)*0.5
!  INTERPOLATION FUNCTION
            X4=1.0-X1-X2-X3
            H(1)=X1*(2.0*X1-1.0)
            H(2)=X2*(2.0*X2-1.0)
            H(3)=X3*(2.0*X3-1.0)
            H(4)=X4*(2.0*X4-1.0)
            H(5)=4.0*X1*X2
            H(6)=4.0*X2*X3
            H(7)=4.0*X1*X3
            H(8)=4.0*X1*X4
            H(9)=4.0*X2*X4
            H(10)=4.0*X3*X4
! DERIVATIVE OF INTERPOLATION FUNCTION
! FOR L1-COORDINATE
            HL1(1)= 4.0*X1-1.0
            HL1(2)= 0.0
            HL1(3)= 0.0
            HL1(4)= 0.0
            HL1(5)= 4.0*X2
            HL1(6)= 0.0
            HL1(7)= 4.0*X3
            HL1(8)= 4.0*X4
            HL1(9)= 0.0
            HL1(10)= 0.0
!  FOR L2-COORDINATE
            HL2(1)= 0.0 
            HL2(2)= 4.0*X2-1.0
            HL2(3)= 0.0
            HL2(4)= 0.0
            HL2(5)= 4.0*X1
            HL2(6)= 4.0*X3
            HL2(7)= 0.0
            HL2(8)= 0.0
            HL2(9)= 4.0*X4
            HL2(10)= 0.0
!  FOR L3-COORDINATE
            HL3(1)= 0.0
            HL3(2)= 0.0
            HL3(3)= 4.0*X3-1.0
            HL3(4)= 0.0
            HL3(5)= 0.0
            HL3(6)= 4.0*X2
            HL3(7)= 4.0*X1
            HL3(8)= 0.0
            HL3(9)= 0.0
            HL3(10)= 4.0*X4
!  FOR L4-COORDINATE
            HL4(1)= 0.0
            HL4(2)= 0.0
            HL4(3)= 0.0 
            HL4(4)= 4.0*X4-1.0
            HL4(5)= 0.0
            HL4(6)= 0.0
            HL4(7)= 0.0
            HL4(8)= 4.0*X1
            HL4(9)= 4.0*X2
            HL4(10)= 4.0*X3
!  JACOBI MATRIX 
            XJ11=0.0
            XJ21=0.0
            XJ31=0.0
            XJ12=0.0
            XJ22=0.0
            XJ32=0.0
            XJ13=0.0
            XJ23=0.0
            XJ33=0.0
            DO I=1,NN
              XJ11=XJ11+(HL1(I)-HL4(I))*XX(I)
              XJ21=XJ21+(HL2(I)-HL4(I))*XX(I)
              XJ31=XJ31+(HL3(I)-HL4(I))*XX(I)
              XJ12=XJ12+(HL1(I)-HL4(I))*YY(I)
              XJ22=XJ22+(HL2(I)-HL4(I))*YY(I)
              XJ32=XJ32+(HL3(I)-HL4(I))*YY(I)
              XJ13=XJ13+(HL1(I)-HL4(I))*ZZ(I)
              XJ23=XJ23+(HL2(I)-HL4(I))*ZZ(I)
              XJ33=XJ33+(HL3(I)-HL4(I))*ZZ(I)
            ENDDO
!DETERMINANT OF JACOBIAN
            DET=XJ11*XJ22*XJ33     &
                 +XJ12*XJ23*XJ31   &
                 +XJ13*XJ21*XJ32   &
                 -XJ13*XJ22*XJ31   &
                 -XJ12*XJ21*XJ33   &
                 -XJ11*XJ23*XJ32
! INVERSION OF JACOBIAN
            DUM=1.0/DET
            XJI11=DUM*( XJ22*XJ33-XJ23*XJ32)
            XJI21=DUM*(-XJ21*XJ33+XJ23*XJ31)
            XJI31=DUM*( XJ21*XJ32-XJ22*XJ31)
            XJI12=DUM*(-XJ12*XJ33+XJ13*XJ32)
            XJI22=DUM*( XJ11*XJ33-XJ13*XJ31)
            XJI32=DUM*(-XJ11*XJ32+XJ12*XJ31)
            XJI13=DUM*( XJ12*XJ23-XJ13*XJ22)
            XJI23=DUM*(-XJ11*XJ23+XJ13*XJ21)
            XJI33=DUM*( XJ11*XJ22-XJ12*XJ21)
!  
            DET=-DET
!  WEIGT VALUE AT GAUSSIAN POINT
            WG=WGT(NG,L1)*WGT(NG,L2)*WGT(NG,L3)*DET      &
               *(1.0-X3)*(1.0-X2-X3)*0.125
!* DUDX DUDY DUDZ ....
            DUDX=0.0
            DUDY=0.0
            DUDZ=0.0
            DVDX=0.0
            DVDY=0.0
            DVDZ=0.0
            DWDX=0.0
            DWDY=0.0
            DWDZ=0.0
            DO J=1,NN
              DNDX=XJI11*(HL1(J)-HL4(J)) +XJI12*(HL2(J)-HL4(J)) +XJI13*(HL3(J)-HL4(J))
              DNDY=XJI21*(HL1(J)-HL4(J)) +XJI22*(HL2(J)-HL4(J)) +XJI23*(HL3(J)-HL4(J))
              DNDZ=XJI31*(HL1(J)-HL4(J)) +XJI32*(HL2(J)-HL4(J)) +XJI33*(HL3(J)-HL4(J))
              DUDX=DUDX+DNDX*UX(J)
              DUDY=DUDY+DNDY*UX(J)
              DUDZ=DUDZ+DNDZ*UX(J)
              DVDX=DVDX+DNDX*UY(J)
              DVDY=DVDY+DNDY*UY(J)
              DVDZ=DVDZ+DNDZ*UY(J)
              DWDX=DWDX+DNDX*UZ(J)
              DWDY=DWDY+DNDY*UZ(J)
              DWDZ=DWDZ+DNDZ*UZ(J)
            ENDDO
!** Fill BMAT
            HEAD=1
            DO J=1,NN
!* DNDX DNDY DNDZ
              DNDX=XJI11*(HL1(J)-HL4(J)) +XJI12*(HL2(J)-HL4(J)) +XJI13*(HL3(J)-HL4(J))
              DNDY=XJI21*(HL1(J)-HL4(J)) +XJI22*(HL2(J)-HL4(J)) +XJI23*(HL3(J)-HL4(J))
              DNDZ=XJI31*(HL1(J)-HL4(J)) +XJI32*(HL2(J)-HL4(J)) +XJI33*(HL3(J)-HL4(J))
!** BL0
              BL0(1,HEAD )=DNDX
              BL0(2,HEAD )=0.0
              BL0(3,HEAD )=0.0
              BL0(4,HEAD )=DNDY
              BL0(5,HEAD )=0.0
              BL0(6,HEAD )=DNDZ
              BL0(1,HEAD+1)=0.0
              BL0(2,HEAD+1)=DNDY
              BL0(3,HEAD+1)=0.0
              BL0(4,HEAD+1)=DNDX
              BL0(5,HEAD+1)=DNDZ
              BL0(6,HEAD+1)=0.0
              BL0(1,HEAD+2)=0.0
              BL0(2,HEAD+2)=0.0
              BL0(3,HEAD+2)=DNDZ
              BL0(4,HEAD+2)=0.0
              BL0(5,HEAD+2)=DNDY
              BL0(6,HEAD+2)=DNDX
!** BL1
              BL1(1,HEAD )=DUDX*DNDX
              BL1(2,HEAD )=DUDY*DNDY
              BL1(3,HEAD )=DUDZ*DNDZ
              BL1(4,HEAD )=DUDY*DNDX+DUDX*DNDY
              BL1(5,HEAD )=DUDY*DNDZ+DUDZ*DNDY
              BL1(6,HEAD )=DUDZ*DNDX+DUDX*DNDZ
              BL1(1,HEAD+1)=DVDX*DNDX
              BL1(2,HEAD+1)=DVDY*DNDY
              BL1(3,HEAD+1)=DVDZ*DNDZ
              BL1(4,HEAD+1)=DVDY*DNDX+DVDX*DNDY
              BL1(5,HEAD+1)=DVDY*DNDZ+DVDZ*DNDY
              BL1(6,HEAD+1)=DVDZ*DNDX+DVDX*DNDZ
              BL1(1,HEAD+2)=DWDX*DNDX
              BL1(2,HEAD+2)=DWDY*DNDY
              BL1(3,HEAD+2)=DWDZ*DNDZ
              BL1(4,HEAD+2)=DWDY*DNDX+DWDX*DNDY
              BL1(5,HEAD+2)=DWDY*DNDZ+DWDZ*DNDY
              BL1(6,HEAD+2)=DWDZ*DNDX+DWDX*DNDZ
!** BNL
              BNL(1,HEAD )=DNDX
              BNL(2,HEAD )=0.0
              BNL(3,HEAD )=0.0
              BNL(4,HEAD )=DNDY
              BNL(5,HEAD )=0.0
              BNL(6,HEAD )=0.0
              BNL(7,HEAD )=DNDZ
              BNL(8,HEAD )=0.0
              BNL(9,HEAD )=0.0
              BNL(1,HEAD+1)=0.0
              BNL(2,HEAD+1)=DNDX
              BNL(3,HEAD+1)=0.0
              BNL(4,HEAD+1)=0.0
              BNL(5,HEAD+1)=DNDY
              BNL(6,HEAD+1)=0.0
              BNL(7,HEAD+1)=0.0
              BNL(8,HEAD+1)=DNDZ
              BNL(9,HEAD+1)=0.0
              BNL(1,HEAD+2)=0.0
              BNL(2,HEAD+2)=0.0
              BNL(3,HEAD+2)=DNDX
              BNL(4,HEAD+2)=0.0
              BNL(5,HEAD+2)=0.0
              BNL(6,HEAD+2)=DNDY
              BNL(7,HEAD+2)=0.0
              BNL(8,HEAD+2)=0.0
              BNL(9,HEAD+2)=DNDZ
              HEAD=HEAD+NDOF
            ENDDO
!** Linear stiffness
!ogasawara till next ogasawara
!**   elasto-plastic matrix
!!
    if(idx_temp .eq. 1) then
      TEMPC=0.0
      TEMPI=0.0
      CURE3D_itg  =0.0
      CURE3D_itg_I=0.0
      DO J=1,NN
        TEMPC=TEMPC+TT (J)*H(J)
        TEMPI=TEMPI+TTI(J)*H(J)
        CURE3D_itg  =CURE3D_itg  +CURE3D_nod  (J)*H(J)
        CURE3D_itg_I=CURE3D_itg_I+CURE3D_nod_I(J)*H(J)
      END DO  
!
              if(CURE3D_itg_I .lt. 0.0) then
                if(iter .eq. -999) then
                  write(ILOG,*) 'stop : CURE3D_itg_I .lt. 0.0 '
                  write(ILOG,*) 'elem = ',icel,'   CURE3D_itg_I= ',CURE3D_itg_I
                  call flush(ILOG)
!                 call hecmw_abort( hecmw_comm_get_comm())
                end if
               CURE3D_itg_I = 0.0
              else if(CURE3D_itg_I .gt. 1.0) then
                if(iter .eq. -999) then
                   write(ILOG,*) 'stop : CURE3D_itg_I .gt. 0.0 '
                   write(ILOG,*) 'elem = ',icel,'   CURE3D_itg_I= ',CURE3D_itg_I
                   call flush(ILOG)
!                  call hecmw_abort( hecmw_comm_get_comm())
                end if
               CURE3D_itg_I = 1.0
              end if
!
              if(CURE3D_itg .lt. 0.0) then
               CURE3D_itg = 0.0
              else if(CURE3D_itg .gt. 1.0) then
               CURE3D_itg = 1.0
              end if
!
!
    if(MID .ne. 2) then
      if(MAT_ISO .eq. 0) then
        ALP = ALPHA*(TEMPC-TEMPI) 
        THERMAL_EPS1 = ALP
        THERMAL_EPS2 = ALP
        THERMAL_EPS3 = ALP
      else if(MAT_ISO .eq. 1) then
        THERMAL_EPS1 = ALPHA_B(1)*(TEMPC-TEMPI)
        THERMAL_EPS2 = ALPHA_B(2)*(TEMPC-TEMPI)
        THERMAL_EPS3 = ALPHA_B(3)*(TEMPC-TEMPI)
      end if
    else if(MID .eq. 2) then
      EE = ccc_rl(1)*dexp(-ccc_rl(2)*TEMPI)*ccc_rl(3)*dexp(ccc_rl(4)*CURE3D_itg_I)
!
      ALP =   (ccc_rl(5)*TEMPI+ccc_rl(6))*(1.0-CURE3D_itg_I)     &
            + (ccc_rl(7)*TEMPI+ccc_rl(8))*CURE3D_itg_I 
      ALP1 = ALP
      ALP = ALP*(TEMPC-TEMPI) - ramd_chm*(CURE3D_itg-CURE3D_itg_I)
      THERMAL_EPS1 = ALP/3.0
      THERMAL_EPS2 = ALP/3.0
      THERMAL_EPS3 = ALP/3.0
     end if
    end if
!!
      call fstr_elpl_DMAT(DMAT,STR(:,L1,L2,L3),FLG_ELPL(L1,L2,L3), h_dash, EE, PP   &
                         ,BIJ,MAT_OTHO1,MAT_ISO)
!!
!!!  20110805  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Add Q force
!!   DMAT_G, dSa by visco 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
            if(idx_visco .eq. 0 .or. idx_visco .eq. 1 ) then
              DO J=1,6
                DSTR_A(J) = 0.0
              END DO

            else if(idx_visco .eq. 2) then
              if(idx_temp .eq. 1 .and. MID .ne. 2) then
                DO J=1,6
                 DSTR_A(J) = 0.0
                END DO
              else
                DO J=1,6
                 DSTR_A(J) = 0.0
                END DO

                 T_G    = ccc_rl(9)
                 C_G    = 1.0 / (EE/2.0/(1.0+PP))  
                 YETA_G = 1.0 / C_G * T_G 
                 T_K    = ccc_rl(10)
                 C_K    = 1.0 / (EE/3.0/(1.0-2.0*PP)) 
                 YETA_K = 1.0 / C_K * T_K 

                 S11=STR(1,L1,L2,L3)
                 S22=STR(2,L1,L2,L3)
                 S33=STR(3,L1,L2,L3)
                 S12=STR(4,L1,L2,L3)
                 S23=STR(5,L1,L2,L3)
                 S13=STR(6,L1,L2,L3)      

                 SMN = (S11+S22+S33)/3.0 
                 SD(1) = S11 - SMN 
                 SD(2) = S22 - SMN 
                 SD(3) = S33 - SMN 
                 SD(4) = S12 
                 SD(5) = S23 
                 SD(6) = S13 

                 C_GT = C_G + 1.0/2.0/YETA_G * DT0 
                 C_KT = C_K + 1.0/2.0/YETA_K * DT0 

                 DO J=1,6
                  DEPS_AD(J) = (1.0/YETA_G)*SD(J) 
                 END DO

                 DEPS_AM = (1.0/YETA_K)*SMN 
!
               if(C_GT .eq. 0.0D+00) then
                 write(ILOG,*) 'stop : C_GT = 0.0'
                 write(ILOG,*) 'elem = ',icel,'   C_GT ',C_GT 
                 call flush(ILOG)
                 call hecmw_abort( hecmw_comm_get_comm())
               end if
               if(C_KT .eq. 0.0D+00) then
                 write(ILOG,*) 'stop : C_KT = 0.0'
                 write(ILOG,*) 'elem = ',icel,'   C_KT ',C_KT 
                 call flush(ILOG)
                 call hecmw_abort( hecmw_comm_get_comm())
               end if
!!
                 DO J=1,6
                  DO K=1,6
                   DMAT(J,K) = 0.0
                  END DO
                 END DO

                 DMAT(1,1) = 1.0/C_KT + 4.0/3.0/C_GT 
                 DMAT(1,2) = 1.0/C_KT - 2.0/3.0/C_GT 
                 DMAT(1,3) = 1.0/C_KT - 2.0/3.0/C_GT 
                 DMAT(2,1) = DMAT(1,2)
                 DMAT(2,2) = 1.0/C_KT + 4.0/3.0/C_GT 
                 DMAT(2,3) = 1.0/C_KT - 2.0/3.0/C_GT 
                 DMAT(3,1) = DMAT(1,3)
                 DMAT(3,2) = DMAT(2,3)
                 DMAT(3,3) = 1.0/C_KT + 4.0/3.0/C_GT 
                 DMAT(4,4) =          1.0/  C_GT 
                 DMAT(5,5) =          1.0/  C_GT 
                 DMAT(6,6) =          1.0/  C_GT 

!
!!!! 20111009  !!!!!!!!!
                  DO J=1,6
                   DSTR_A(J) = DEPS_AD(J)/C_GT 
                  END DO
                  DO J=1,3
                   DSTR_A(J) = DSTR_A(J) + DEPS_AM/C_KT 
                  END DO
!
                if(i_step .eq. 1 .and. iter .le. 1) then
                  DO J=1,6
                    DSTR_A(J) = DSTR_A(J)*DT0
                  END DO
                else if(i_step .ge. 2 .and. iter .le. 1) then
                  DO J=1,6
                    DSTR_A(J) = DSTR_A(J)*DT0
                  END DO
                else
                  DO J=1,6
                    DSTR_A(J) = 0.0
                  END DO
                end if
!
              end if
            else if(idx_visco .eq. 3) then
            else if(idx_visco .eq. 4) then
            else if(idx_visco .eq. 5) then
            else if(idx_visco .eq. 6) then
              if(idx_temp .eq. 1 .and. MID .ne. 2) then
                DO J=1,6
                 DSTR_A(J) = 0.0
                END DO

              else
                DO J=1,6
                 DSTR_A(J) = 0.0
                END DO

                 T_G    = ccc_rl(9)
                 C_G    = 1.0 / (EE/2.0/(1.0+PP))  
                 YETA_G = 1.0 / C_G * T_G 
                 T_K    = ccc_rl(10)
                 C_K    = 1.0 / (EE/3.0/(1.0-2.0*PP)) 
                 YETA_K = 1.0 / C_K * T_K 

                 T_G1    = T_G
                 C_G1    = C_G 
                 YETA_G1 = YETA_G 
                 T_K1    = T_K
                 C_K1    = C_K 
                 YETA_K1 = YETA_K 

                 T_G2    = T_G
                 C_G2    = C_G 
                 YETA_G2 = YETA_G 
                 T_K2    = T_K
                 C_K2    = C_K 
                 YETA_K2 = YETA_K 

                 S11=STR(1,L1,L2,L3)
                 S22=STR(2,L1,L2,L3)
                 S33=STR(3,L1,L2,L3)
                 S12=STR(4,L1,L2,L3)
                 S23=STR(5,L1,L2,L3)
                 S13=STR(6,L1,L2,L3)      

                 SMN = (S11+S22+S33)/3.0 
                 SD(1) = S11 - SMN 
                 SD(2) = S22 - SMN 
                 SD(3) = S33 - SMN 
                 SD(4) = S12 
                 SD(5) = S23 
                 SD(6) = S13 

                 DO J=1,6
                   E_V1(J) = STRAIN_V1(J,L1,L2,L3) 
                   E_V2(J) = STRAIN_V2(J,L1,L2,L3) 
                 END DO
                   EMN_V1 = (E_V1(1)+E_V1(2)+E_V1(3))/3.0 
                   EMN_V2 = (E_V2(1)+E_V2(2)+E_V2(3))/3.0 
                 DO J=1,3
                   ED_V1(  J) = E_V1(  J) - EMN_V1 
                   ED_V1(3+J) = E_V1(3+J) 
                   ED_V2(  J) = E_V2(  J) - EMN_V2 
                   ED_V2(3+J) = E_V2(3+J) 
                 END DO
                 
                 C_GT = C_G + 1.0/2.0/YETA_G * DT0 + C_G1*(1.0-(1.0-dexp(-DT0/T_G1))*T_G1/DT0)   & 
                                                   + C_G2*(1.0-(1.0-dexp(-DT0/T_G2))*T_G2/DT0) 
                 C_KT = C_K + 1.0/2.0/YETA_K * DT0 + C_K1*(1.0-(1.0-dexp(-DT0/T_K1))*T_K1/DT0)   & 
                                                   + C_K2*(1.0-(1.0-dexp(-DT0/T_K2))*T_K2/DT0) 
                 DO J=1,6
                  DEPS_AD(J) = (1.0/YETA_G + 1.0/YETA_G1*(1.0-dexp(-DT0/T_G1))*T_G1/DT0            & 
                                           + 1.0/YETA_G2*(1.0-dexp(-DT0/T_G2))*T_G2/DT0)*SD(J)     & 
                                           -(2.0/T_G1   *(1.0-dexp(-DT0/T_G1))*T_G1/DT0)*ED_V1(J)  & 
                                           -(2.0/T_G2   *(1.0-dexp(-DT0/T_G2))*T_G2/DT0)*ED_V2(J)   
                 END DO

                 DEPS_AM = (1.0/YETA_K + 1.0/YETA_K1*(1.0-dexp(-DT0/T_K1))*T_K1/DT0          & 
                                       + 1.0/YETA_K2*(1.0-dexp(-DT0/T_K2))*T_K2/DT0)*SMN     & 
                                       -(3.0/T_K1   *(1.0-dexp(-DT0/T_K1))*T_K1/DT0)*EMN_V1  & 
                                       -(3.0/T_K2   *(1.0-dexp(-DT0/T_K2))*T_K2/DT0)*EMN_V2   
!
               if(C_GT .eq. 0.0D+00) then
                 write(ILOG,*) 'stop : C_GT = 0.0'
                 write(ILOG,*) 'elem = ',icel,'   C_GT ',C_GT 
                 call flush(ILOG)
                 call hecmw_abort( hecmw_comm_get_comm())
               end if
               if(C_KT .eq. 0.0D+00) then
                 write(ILOG,*) 'stop : C_KT = 0.0'
                 write(ILOG,*) 'elem = ',icel,'   C_KT ',C_KT 
                 call flush(ILOG)
                 call hecmw_abort( hecmw_comm_get_comm())
               end if
!!
                 DO J=1,6
                  DO K=1,6
                   DMAT(J,K) = 0.0
                  END DO
                 END DO

                 DMAT(1,1) = 1.0/C_KT + 4.0/3.0/C_GT 
                 DMAT(1,2) = 1.0/C_KT - 2.0/3.0/C_GT 
                 DMAT(1,3) = 1.0/C_KT - 2.0/3.0/C_GT 
                 DMAT(2,1) = DMAT(1,2)
                 DMAT(2,2) = 1.0/C_KT + 4.0/3.0/C_GT 
                 DMAT(2,3) = 1.0/C_KT - 2.0/3.0/C_GT 
                 DMAT(3,1) = DMAT(1,3)
                 DMAT(3,2) = DMAT(2,3)
                 DMAT(3,3) = 1.0/C_KT + 4.0/3.0/C_GT 
                 DMAT(4,4) =          1.0/  C_GT 
                 DMAT(5,5) =          1.0/  C_GT 
                 DMAT(6,6) =          1.0/  C_GT 
!
!!!! 20111009  !!!!!!!!!
                  DO J=1,6
                   DSTR_A(J) = DEPS_AD(J)/C_GT 
                  END DO
                  DO J=1,3
                   DSTR_A(J) = DSTR_A(J) + DEPS_AM/C_KT 
                  END DO
!
                if(i_step .eq. 1 .and. iter .le. 1) then
                  DO J=1,6
                    DSTR_A(J) = DSTR_A(J)*DT0
                  END DO
                else if(i_step .ge. 2 .and. iter .le. 1) then
                  DO J=1,6
                    DSTR_A(J) = DSTR_A(J)*DT0
                  END DO
                else
                  DO J=1,6
                    DSTR_A(J) = 0.0
                  END DO
                end if
!
              end if
            end if
!!!  20110805  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!
!** element stiffness
!!
      if(idx_elpl .eq. -1 .or. idx_elpl .eq. 2) then
        DO J=1,6
          DO K=1,NN*NDOF
            BL1(J,K)=0.0
          END DO
        END DO
        DO J=1,9
          DO K=1,NN*NDOF
            BNL(J,K)=0.0
          END DO
        END DO
      end if
!!
            DO J=1,6
              DO K=1,NN*NDOF
                DB(J,K)=0.0
              ENDDO
            ENDDO
            DO J=1,6
              DO K=1,NN*NDOF
                DB(J,K)=DB(J,K)                   &
                +DMAT(J,1)*(BL0(1,K)+BL1(1,K))    &
                +DMAT(J,2)*(BL0(2,K)+BL1(2,K))    &
                +DMAT(J,3)*(BL0(3,K)+BL1(3,K))    &
                +DMAT(J,4)*(BL0(4,K)+BL1(4,K))    &
                +DMAT(J,5)*(BL0(5,K)+BL1(5,K))    &
                +DMAT(J,6)*(BL0(6,K)+BL1(6,K))
              ENDDO
            ENDDO
            DO J=1,NN*NDOF
              DO K=1,J
                NUM=(J-1)*J/2+K
                SS(NUM)=SS(NUM)                       &
                +( (BL0(1,K)+BL1(1,K))*DB(1,J)        &
                +  (BL0(2,K)+BL1(2,K))*DB(2,J)        &
                +  (BL0(3,K)+BL1(3,K))*DB(3,J)        &
                +  (BL0(4,K)+BL1(4,K))*DB(4,J)        &
                +  (BL0(5,K)+BL1(5,K))*DB(5,J)        &
                +  (BL0(6,K)+BL1(6,K))*DB(6,J) )*WG
              ENDDO
            ENDDO
!** NonLinear stiffness
!* Current stress
!
            S11=STR(1,L1,L2,L3)
            S22=STR(2,L1,L2,L3)
            S33=STR(3,L1,L2,L3)
            S12=STR(4,L1,L2,L3)
            S23=STR(5,L1,L2,L3)
            S13=STR(6,L1,L2,L3)      
            DO J=1,9
              DO K=1,9
                SMAT(J,K)=0.0
              ENDDO
            ENDDO
            DO J=1,3
              SMAT(J  ,J  )=S11
              SMAT(J+3,J+3)=S22
              SMAT(J+6,J+6)=S33
              SMAT(J+3,J+6)=S23
              SMAT(J  ,J+3)=S12
              SMAT(J  ,J+6)=S13
              SMAT(J+6,J+3)=S23
              SMAT(J+3,J  )=S12
              SMAT(J+6,J  )=S13
            ENDDO
            DO J=1,9
              DO K=1,NN*NDOF
                DB(J,K)=0.0
              ENDDO
            ENDDO
            DO J=1,9
              DO K=1,NN*NDOF
                DB(J,K)=DB(J,K)         &
                +SMAT(J,1)*BNL(1,K)     &
                +SMAT(J,2)*BNL(2,K)     &
                +SMAT(J,3)*BNL(3,K)     &
                +SMAT(J,4)*BNL(4,K)     &
                +SMAT(J,5)*BNL(5,K)     &
                +SMAT(J,6)*BNL(6,K)     &
                +SMAT(J,7)*BNL(7,K)     &
                +SMAT(J,8)*BNL(8,K)     &
                +SMAT(J,9)*BNL(9,K)
              ENDDO
            ENDDO
            DO J=1,NN*NDOF
              DO K=1,J
                NUM=(J-1)*J/2+K
                SS(NUM)=SS(NUM)+           &
                (BNL(1,K)*DB(1,J)          &
                +BNL(2,K)*DB(2,J)          &
                +BNL(3,K)*DB(3,J)          &
                +BNL(4,K)*DB(4,J)          &
                +BNL(5,K)*DB(5,J)          &
                +BNL(6,K)*DB(6,J)          &
                +BNL(7,K)*DB(7,J)          &
                +BNL(8,K)*DB(8,J)          &
                +BNL(9,K)*DB(9,J) )*WG
              ENDDO
            ENDDO
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!* Add Q force
!!-----  dSa by de_therm ----
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     if(idx_temp .eq. 1) then
       if(iter .le. 1) then
            DO J=1,6
              DS(J)=0.0
            ENDDO
!
            DO I=1,3
             DO J=1,3
              EPS_G (I,J) = 0.0
              EPS_GB(I,J) = 0.0
             END DO
            END DO

              EPS_GB(1,1) = THERMAL_EPS1 
              EPS_GB(2,2) = THERMAL_EPS2 
              EPS_GB(3,3) = THERMAL_EPS3 

           if(MAT_ISO .eq. 0) then
              DO I=1,3
               DO J=1,3
                  EPS_G(I,J) = EPS_GB(I,J) 
               END DO
              END DO
           else if(MAT_ISO .eq. 1) then
              DO I=1,3
               DO J=1,3
                DO K=1,3
                 DO L=1,3
                  EPS_G(I,J) = EPS_G(I,J) + BIJ(I,K)*BIJ(J,L)*EPS_GB(K,L)
                 END DO
                END DO
               END DO
              END DO
           end if

              EPSTH(1)=EPS_G(1,1)
              EPSTH(2)=EPS_G(2,2)
              EPSTH(3)=EPS_G(3,3)
              EPSTH(4)=EPS_G(1,2)
              EPSTH(5)=EPS_G(2,3)
              EPSTH(6)=EPS_G(3,1)
!**
!** SET DSGM  {ds}=[D]{de_therm}
!**
          if(idx_visco .eq. 0 .or. idx_visco .eq. 1 ) then
            DO J=1,6
             DO K=1,6
               DS(K)=DS(K)+DMAT(K,J)*EPSTH(J)
             ENDDO
            ENDDO
          else if(idx_visco .ge. 2) then
           if(MID .ne. 2) then
             DO J=1,6
              DO K=1,6
                DS(K)=DS(K)+DMAT(K,J)*EPSTH(J)
              ENDDO
             ENDDO
           else if(MID .eq. 2) then
             DO J=1,6
              DO K=1,6
                DS(K)=DS(K)+DMAT(K,J)*EPSTH(J)
              ENDDO
             ENDDO
           end if
          end if
!
            DO J=1,6
              DSTR_A(J) = DSTR_A(J) + DS(J)
            ENDDO
       end if
     end if
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
            DO K=1,NN*NDOF
              QF(K)=QF(K)                                 &
              +((S11-DSTR_A(1))*(BL0(1,K)+BL1(1,K))       &
              + (S22-DSTR_A(2))*(BL0(2,K)+BL1(2,K))       &
              + (S33-DSTR_A(3))*(BL0(3,K)+BL1(3,K))       &
              + (S12-DSTR_A(4))*(BL0(4,K)+BL1(4,K))       &
              + (S23-DSTR_A(5))*(BL0(5,K)+BL1(5,K))       &
              + (S13-DSTR_A(6))*(BL0(6,K)+BL1(6,K)) )*WG     
!!!!  20111021  -----
              QF_RES(K)=QF_RES(K)             &
              +(S11*(BL0(1,K)+BL1(1,K))       &
              + S22*(BL0(2,K)+BL1(2,K))       &
              + S33*(BL0(3,K)+BL1(3,K))       &
              + S12*(BL0(4,K)+BL1(4,K))       &
              + S23*(BL0(5,K)+BL1(5,K))       &
              + S13*(BL0(6,K)+BL1(6,K)) )*WG     
!
              if(i_step .eq. 1 .and. iter .ge. 1) then
                QF_RES(K)=QF(K)            
              else if(i_step .ge. 2 .and. iter .ge. 1) then
                QF_RES(K)=QF(K)            
              end if
            ENDDO
!
          ENDDO
        ENDDO
      ENDDO
      RETURN
      END SUBROUTINE STF_C3D10TL
!----------------------------------------------------------------------*
       SUBROUTINE UPDATE_C3D20TL(ILOG,XX,YY,ZZ,DISP,DU,EE,PP,STRESS   &
                           ,idx_elpl,h_dash,FLG_ELPL,STR_YLD ,STRAIN,EP_BAR ,iter)
!----------------------------------------------------------------------*
!
! UPDATE 3D 20 NODE PLANE ELEMENT
!
      use hecmw
      use gauss_integration 
      IMPLICIT NONE
! GAUSS POINT
      INTEGER(kind=kint) NG
      PARAMETER(NG=3)
! I/F VARIABLES
      INTEGER(kind=kint) ILOG
      REAL(kind=kreal) XX(*),YY(*),ZZ(*),DISP(*),DU(*),EE,PP
      REAL(kind=kreal) STRESS(6,NG,NG,NG)
      INTEGER(kind=kint) iter
      INTEGER(kind=kint) idx_elpl,FLG_ELPL(NG,NG,NG)
      REAL(kind=kreal) h_dash,STR_YLD(NG,NG,NG),EP_BAR(NG,NG,NG),STRAIN(6,NG,NG,NG)
! LOCAL VARIABLES
      INTEGER(kind=kint) NN,NDOF
      PARAMETER(NN=20,NDOF=3)
      INTEGER(kind=kint) I,J,K,NUM
      INTEGER(kind=kint) LX,LY,LZ
      INTEGER(kind=kint) HEAD
      REAL(kind=kreal) RI,SI,TI,RP,SP,TP,RM,SM,TM
      REAL(kind=kreal) H(NN),HR(NN),HS(NN),HT(NN)
      REAL(kind=kreal) XJ11,XJ21,XJ31,XJ12,XJ22,XJ32,XJ13,XJ23,XJ33,DET,WG,DUM
      REAL(kind=kreal) XJI11,XJI21,XJI31,XJI12,XJI22,XJI32,XJI13,XJI23,XJI33
      REAL(kind=kreal) DMAT(6,6),SMAT(9,9)
      REAL(kind=kreal) DNDX,DNDY,DNDZ
      REAL(kind=kreal) DUDX,DUDY,DUDZ
      REAL(kind=kreal) DVDX,DVDY,DVDZ
      REAL(kind=kreal) DWDX,DWDY,DWDZ
      REAL(kind=kreal) BL0(6,NN*NDOF),BL1(6,NN*NDOF)
      REAL(kind=kreal) UX(NN),UY(NN),UZ(NN)
      REAL(kind=kreal) S11,S22,S33,S23,S13,S12
      REAL(kind=kreal) DE(6),DS(6)

      INTEGER(kind=kint) MAT_ISO 
      REAL(kind=kreal) BIJ(3,3),MAT_OTHO1(9) 
!* Total displacement
      DO I=1,NN
        UX(I)=DISP(NDOF*I-2)
        UY(I)=DISP(NDOF*I-1)
        UZ(I)=DISP(NDOF*I  )
      ENDDO
! SET DMAT
      DO J=1,6
        DO I=1,6
          DMAT(I,J)=0.0
        ENDDO
      ENDDO 
      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
! LOOP FOR INTEGRATION POINTS      
      DO LX=1,NG
        RI=XG(NG,LX)
        DO LY=1,NG
          SI=XG(NG,LY)
          DO LZ=1,NG
            TI=XG(NG,LZ)
            RP=1.0+RI
            SP=1.0+SI
            TP=1.0+TI
            RM=1.0-RI
            SM=1.0-SI
            TM=1.0-TI
!  INTERPOLATION FUNCTION
            H(1)=-0.125*RM*SM*TM*(2.0+RI+SI+TI)
            H(2)=-0.125*RP*SM*TM*(2.0-RI+SI+TI)
            H(3)=-0.125*RP*SP*TM*(2.0-RI-SI+TI)
            H(4)=-0.125*RM*SP*TM*(2.0+RI-SI+TI)
            H(5)=-0.125*RM*SM*TP*(2.0+RI+SI-TI)
            H(6)=-0.125*RP*SM*TP*(2.0-RI+SI-TI)
            H(7)=-0.125*RP*SP*TP*(2.0-RI-SI-TI)
            H(8)=-0.125*RM*SP*TP*(2.0+RI-SI-TI)
            H(9)=0.25*(1.0-RI**2)*SM*TM
            H(10)=0.25*RP*(1.0-SI**2)*TM
            H(11)=0.25*(1.0-RI**2)*SP*TM
            H(12)=0.25*RM*(1.0-SI**2)*TM
            H(13)=0.25*(1.0-RI**2)*SM*TP
            H(14)=0.25*RP*(1.0-SI**2)*TP
            H(15)=0.25*(1.0-RI**2)*SP*TP
            H(16)=0.25*RM*(1.0-SI**2)*TP
            H(17)=0.25*RM*SM*(1.0-TI**2)
            H(18)=0.25*RP*SM*(1.0-TI**2)
            H(19)=0.25*RP*SP*(1.0-TI**2)
            H(20)=0.25*RM*SP*(1.0-TI**2)
!  DERIVATIVE OF INTERPOLATION FUNCTION
!  FOR R-COORDINATE
            HR(1)=-0.125*RM*SM*TM+0.125*SM*TM*(2.0+RI+SI+TI)
            HR(2)=+0.125*RP*SM*TM-0.125*SM*TM*(2.0-RI+SI+TI)
            HR(3)=+0.125*RP*SP*TM-0.125*SP*TM*(2.0-RI-SI+TI)
            HR(4)=-0.125*RM*SP*TM+0.125*SP*TM*(2.0+RI-SI+TI)
            HR(5)=-0.125*RM*SM*TP+0.125*SM*TP*(2.0+RI+SI-TI)
            HR(6)=+0.125*RP*SM*TP-0.125*SM*TP*(2.0-RI+SI-TI)
            HR(7)=+0.125*RP*SP*TP-0.125*SP*TP*(2.0-RI-SI-TI)
            HR(8)=-0.125*RM*SP*TP+0.125*SP*TP*(2.0+RI-SI-TI)
            HR(9 )=-0.50*RI*SM*TM
            HR(10)=+0.25*(1.0-SI**2)*TM
            HR(11)=-0.50*RI*SP*TM
            HR(12)=-0.25*(1.0-SI**2)*TM
            HR(13)=-0.50*RI*SM*TP
            HR(14)=+0.25*(1.0-SI**2)*TP
            HR(15)=-0.50*RI*SP*TP
            HR(16)=-0.25*(1.0-SI**2)*TP
            HR(17)=-0.25*SM*(1.0-TI**2)
            HR(18)=+0.25*SM*(1.0-TI**2)
            HR(19)=+0.25*SP*(1.0-TI**2)
            HR(20)=-0.25*SP*(1.0-TI**2)
!  FOR S-COORDINATE
            HS(1)=-0.125*RM*SM*TM+0.125*RM*TM*(2.0+RI+SI+TI)
            HS(2)=-0.125*RP*SM*TM+0.125*RP*TM*(2.0-RI+SI+TI)
            HS(3)=+0.125*RP*SP*TM-0.125*RP*TM*(2.0-RI-SI+TI)
            HS(4)=+0.125*RM*SP*TM-0.125*RM*TM*(2.0+RI-SI+TI)
            HS(5)=-0.125*RM*SM*TP+0.125*RM*TP*(2.0+RI+SI-TI)
            HS(6)=-0.125*RP*SM*TP+0.125*RP*TP*(2.0-RI+SI-TI)
            HS(7)=+0.125*RP*SP*TP-0.125*RP*TP*(2.0-RI-SI-TI)
            HS(8)=+0.125*RM*SP*TP-0.125*RM*TP*(2.0+RI-SI-TI)
            HS(9)=-0.25*(1.0-RI**2)*TM
            HS(10)=-0.50*RP*SI*TM
            HS(11)=+0.25*(1.0-RI**2)*TM
            HS(12)=-0.50*RM*SI*TM
            HS(13)=-0.25*(1.0-RI**2)*TP
            HS(14)=-0.50*RP*SI*TP
            HS(15)=+0.25*(1.0-RI**2)*TP
            HS(16)=-0.50*RM*SI*TP
            HS(17)=-0.25*RM*(1.0-TI**2)
            HS(18)=-0.25*RP*(1.0-TI**2)
            HS(19)=+0.25*RP*(1.0-TI**2)
            HS(20)=+0.25*RM*(1.0-TI**2)
!  FOR T-COORDINATE
            HT(1)=-0.125*RM*SM*TM+0.125*RM*SM*(2.0+RI+SI+TI)
            HT(2)=-0.125*RP*SM*TM+0.125*RP*SM*(2.0-RI+SI+TI)
            HT(3)=-0.125*RP*SP*TM+0.125*RP*SP*(2.0-RI-SI+TI)
            HT(4)=-0.125*RM*SP*TM+0.125*RM*SP*(2.0+RI-SI+TI)
            HT(5)=+0.125*RM*SM*TP-0.125*RM*SM*(2.0+RI+SI-TI)
            HT(6)=+0.125*RP*SM*TP-0.125*RP*SM*(2.0-RI+SI-TI)
            HT(7)=+0.125*RP*SP*TP-0.125*RP*SP*(2.0-RI-SI-TI)
            HT(8)=+0.125*RM*SP*TP-0.125*RM*SP*(2.0+RI-SI-TI)
            HT(9)=-0.25*(1.0-RI**2)*SM
            HT(10)=-0.25*RP*(1.0-SI**2)
            HT(11)=-0.25*(1.0-RI**2)*SP
            HT(12)=-0.25*RM*(1.0-SI**2)
            HT(13)=0.25*(1.0-RI**2)*SM
            HT(14)=0.25*RP*(1.0-SI**2)
            HT(15)=0.25*(1.0-RI**2)*SP
            HT(16)=0.25*RM*(1.0-SI**2)
            HT(17)=-0.5*RM*SM*TI
            HT(18)=-0.5*RP*SM*TI
            HT(19)=-0.5*RP*SP*TI
            HT(20)=-0.5*RM*SP*TI
!  JACOBI MATRIX 
            XJ11=0.0
            XJ21=0.0
            XJ31=0.0
            XJ12=0.0
            XJ22=0.0
            XJ32=0.0
            XJ13=0.0
            XJ23=0.0
            XJ33=0.0
            DO I=1,NN
              XJ11=XJ11+HR(I)*XX(I)
              XJ21=XJ21+HS(I)*XX(I)
              XJ31=XJ31+HT(I)*XX(I)
              XJ12=XJ12+HR(I)*YY(I)
              XJ22=XJ22+HS(I)*YY(I)
              XJ32=XJ32+HT(I)*YY(I)
              XJ13=XJ13+HR(I)*ZZ(I)
              XJ23=XJ23+HS(I)*ZZ(I)
              XJ33=XJ33+HT(I)*ZZ(I)
            ENDDO
!DETERMINANT OF JACOBIAN
            DET=XJ11*XJ22*XJ33   &
            +XJ12*XJ23*XJ31      &
            +XJ13*XJ21*XJ32      &
            -XJ13*XJ22*XJ31      &
            -XJ12*XJ21*XJ33      &
            -XJ11*XJ23*XJ32
! INVERSION OF JACOBIAN
            DUM=1.0/DET
            XJI11=DUM*( XJ22*XJ33-XJ23*XJ32)
            XJI21=DUM*(-XJ21*XJ33+XJ23*XJ31)
            XJI31=DUM*( XJ21*XJ32-XJ22*XJ31)
            XJI12=DUM*(-XJ12*XJ33+XJ13*XJ32)
            XJI22=DUM*( XJ11*XJ33-XJ13*XJ31)
            XJI32=DUM*(-XJ11*XJ32+XJ12*XJ31)
            XJI13=DUM*( XJ12*XJ23-XJ13*XJ22)
            XJI23=DUM*(-XJ11*XJ23+XJ13*XJ21)
            XJI33=DUM*( XJ11*XJ22-XJ12*XJ21)
!* DUDX DUDY DUDZ ....
            DUDX=0.0
            DUDY=0.0
            DUDZ=0.0
            DVDX=0.0
            DVDY=0.0
            DVDZ=0.0
            DWDX=0.0
            DWDY=0.0
            DWDZ=0.0
            DO J=1,NN
              DNDX=XJI11*HR(J)+XJI12*HS(J)+XJI13*HT(J)
              DNDY=XJI21*HR(J)+XJI22*HS(J)+XJI23*HT(J)
              DNDZ=XJI31*HR(J)+XJI32*HS(J)+XJI33*HT(J)
              DUDX=DUDX+DNDX*UX(J)
              DUDY=DUDY+DNDY*UX(J)
              DUDZ=DUDZ+DNDZ*UX(J)
              DVDX=DVDX+DNDX*UY(J)
              DVDY=DVDY+DNDY*UY(J)
              DVDZ=DVDZ+DNDZ*UY(J)
              DWDX=DWDX+DNDX*UZ(J)
              DWDY=DWDY+DNDY*UZ(J)
              DWDZ=DWDZ+DNDZ*UZ(J)
            ENDDO
!** Fill BMAT
            HEAD=1
            DO J=1,NN
!* DNDX DNDY DNDZ
              DNDX=XJI11*HR(J)+XJI12*HS(J)+XJI13*HT(J)
              DNDY=XJI21*HR(J)+XJI22*HS(J)+XJI23*HT(J)
              DNDZ=XJI31*HR(J)+XJI32*HS(J)+XJI33*HT(J)
!** BL0
              BL0(1,HEAD )=DNDX
              BL0(2,HEAD )=0.0
              BL0(3,HEAD )=0.0
              BL0(4,HEAD )=DNDY
              BL0(5,HEAD )=0.0
              BL0(6,HEAD )=DNDZ
              BL0(1,HEAD+1)=0.0
              BL0(2,HEAD+1)=DNDY
              BL0(3,HEAD+1)=0.0
              BL0(4,HEAD+1)=DNDX
              BL0(5,HEAD+1)=DNDZ
              BL0(6,HEAD+1)=0.0
              BL0(1,HEAD+2)=0.0
              BL0(2,HEAD+2)=0.0
              BL0(3,HEAD+2)=DNDZ
              BL0(4,HEAD+2)=0.0
              BL0(5,HEAD+2)=DNDY
              BL0(6,HEAD+2)=DNDX
!** BL1
              BL1(1,HEAD )=DUDX*DNDX
              BL1(2,HEAD )=DUDY*DNDY
              BL1(3,HEAD )=DUDZ*DNDZ
              BL1(4,HEAD )=DUDY*DNDX+DUDX*DNDY
              BL1(5,HEAD )=DUDY*DNDZ+DUDZ*DNDY
              BL1(6,HEAD )=DUDZ*DNDX+DUDX*DNDZ
              BL1(1,HEAD+1)=DVDX*DNDX
              BL1(2,HEAD+1)=DVDY*DNDY
              BL1(3,HEAD+1)=DVDZ*DNDZ
              BL1(4,HEAD+1)=DVDY*DNDX+DVDX*DNDY
              BL1(5,HEAD+1)=DVDY*DNDZ+DVDZ*DNDY
              BL1(6,HEAD+1)=DVDZ*DNDX+DVDX*DNDZ
              BL1(1,HEAD+2)=DWDX*DNDX
              BL1(2,HEAD+2)=DWDY*DNDY
              BL1(3,HEAD+2)=DWDZ*DNDZ
              BL1(4,HEAD+2)=DWDY*DNDX+DWDX*DNDY
              BL1(5,HEAD+2)=DWDY*DNDZ+DWDZ*DNDY
              BL1(6,HEAD+2)=DWDZ*DNDX+DWDX*DNDZ
              HEAD=HEAD+NDOF
            ENDDO
!** DE & DS ***
            DO J=1,6
              DE(J)=0.0
              DS(J)=0.0
            ENDDO
            DO J=1,NN*NDOF
              DO K=1,6
                DE(K)=DE(K)+(BL0(K,J)+BL1(K,J))*DU(J)
              ENDDO
            ENDDO

            DO J=1,6
              STRAIN(J,LX,LY,LZ)=STRAIN(J,LX,LY,LZ)+DE(J)
            ENDDO
!** Update stress
!ogasawara till next ogasawara
!---  elastic analysis
      if(idx_elpl .eq. 0) then
            DO J=1,6
              DO K=1,6
                DS(K)=DS(K)+DMAT(K,J)*DE(J)
              ENDDO
            ENDDO
            DO J=1,6
              STRESS(J,LX,LY,LZ)=STRESS(J,LX,LY,LZ)+DS(J)
            ENDDO
!---  elasto-plastic analysis
!      yield judgement
!      stress integration
      else if(idx_elpl .ge. 1) then
        call fstr_elpl_INTG(DMAT,STRESS(:,LX,LY,LZ),STR_YLD(LX,LY,LZ)     &
                         ,FLG_ELPL(LX,LY,LZ), h_dash, EP_BAR(LX,LY,LZ)    &
                         ,EE, PP, DE, iter,BIJ,MAT_OTHO1,MAT_ISO)
      end if
!ogasawara till next ogasawara
!
          ENDDO
        ENDDO
      ENDDO
      RETURN
      END SUBROUTINE UPDATE_C3D20TL
!----------------------------------------------------------------------*
       SUBROUTINE UPDATE_C3D15TL(ILOG,XX,YY,ZZ,DISP,DU,EE,PP,STRESS   &
                           ,idx_elpl,h_dash,FLG_ELPL,STR_YLD ,STRAIN,EP_BAR ,iter)
!----------------------------------------------------------------------*
!
! UPDATE 3D 6 NODE PLANE ELEMENT
!
      use hecmw
      use gauss_integration 
      IMPLICIT NONE
! GAUSS POINT
      INTEGER(kind=kint) NG
      PARAMETER(NG=3)
! I/F VARIABLES
      INTEGER(kind=kint) ILOG
      REAL(kind=kreal) XX(*),YY(*),ZZ(*),DISP(*),DU(*),EE,PP
      REAL(kind=kreal) STRESS(6,NG,NG,NG)
      INTEGER(kind=kint) iter
      INTEGER(kind=kint) idx_elpl,FLG_ELPL(NG,NG,NG)
      REAL(kind=kreal) h_dash,STR_YLD(NG,NG,NG),EP_BAR(NG,NG,NG),STRAIN(6,NG,NG,NG)
! LOCAL VARIABLES
      INTEGER(kind=kint) NN,NDOF,ISIZE
      PARAMETER(NN=15,NDOF=3,ISIZE=(NN*NDOF)*(NN*NDOF+1)/2)
      INTEGER(kind=kint) I,J,K,NUM
      INTEGER(kind=kint) L1,L2,LZ
      INTEGER(kind=kint) HEAD
      REAL(kind=kreal) RI,SI,TI,RP,SP,TP,RM,SM,TM
      REAL(kind=kreal) H(NN),HL1(NN),HL2(NN),HL3(NN),HZ(NN)
      REAL(kind=kreal) X1,X2,X3,XL1,XL2,ZI
      REAL(kind=kreal) XJ11,XJ21,XJ31,XJ12,XJ22,XJ32,XJ13,XJ23,XJ33,DET,WG,DUM
      REAL(kind=kreal) XJI11,XJI21,XJI31,XJI12,XJI22,XJI32,XJI13,XJI23,XJI33
      REAL(kind=kreal) DMAT(6,6),SMAT(9,9)
      REAL(kind=kreal) DNDX,DNDY,DNDZ
      REAL(kind=kreal) DUDX,DUDY,DUDZ
      REAL(kind=kreal) DVDX,DVDY,DVDZ
      REAL(kind=kreal) DWDX,DWDY,DWDZ
      REAL(kind=kreal) BL0(6,NN*NDOF),BL1(6,NN*NDOF)
      REAL(kind=kreal) UX(NN),UY(NN),UZ(NN)
      REAL(kind=kreal) S11,S22,S33,S23,S13,S12
      REAL(kind=kreal) DE(6),DS(6)

      INTEGER(kind=kint) MAT_ISO 
      REAL(kind=kreal) BIJ(3,3),MAT_OTHO1(9) 
!* Total displacement
      DO I=1,NN
        UX(I)=DISP(NDOF*I-2)
        UY(I)=DISP(NDOF*I-1)
        UZ(I)=DISP(NDOF*I  )
      ENDDO
! SET DMAT
      DO J=1,6
        DO I=1,6
          DMAT(I,J)=0.0
        ENDDO
      ENDDO 
      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
! LOOP FOR INTEGRATION POINTS      
      DO LZ=1,NG
        ZI=XG(NG,LZ)
        DO L2=1,NG
          XL2=XG(NG,L2)
          X2 =(XL2+1.0)*0.5
          DO L1=1,NG
            XL1=XG(NG,L1)
            X1=0.5*(1.0-X2)*(XL1+1.0)
!  INTERPOLATION FUNCTION
            X3=1.0-X1-X2
            H(1)= 0.5*X1*(2.0*X1-2.-ZI)*(1.0-ZI)
            H(2)= 0.5*X2*(2.0*X2-2.-ZI)*(1.0-ZI)
            H(3)= 0.5*X3*(2.0*X3-2.-ZI)*(1.0-ZI)
            H(4)= 0.5*X1*(2.0*X1-2.+ZI)*(1.0+ZI)
            H(5)= 0.5*X2*(2.0*X2-2.+ZI)*(1.0+ZI)
            H(6)= 0.5*X3*(2.0*X3-2.+ZI)*(1.0+ZI)
            H(7)= 2.0*X1*X2*(1.0-ZI)
            H(8)= 2.0*X2*X3*(1.0-ZI)
            H(9)= 2.0*X1*X3*(1.0-ZI)
            H(10)=2.0*X1*X2*(1.0+ZI)
            H(11)=2.0*X2*X3*(1.0+ZI)
            H(12)=2.0*X1*X3*(1.0+ZI)
            H(13)=X1*(1.0-ZI**2)
            H(14)=X2*(1.0-ZI**2)
            H(15)=X3*(1.0-ZI**2)
!  DERIVATIVE OF INTERPOLATION FUNCTION
!  FOR L1-COORDINATE
            HL1(1)= 0.5*(4.0*X1-2.-ZI)*(1.0-ZI)
            HL1(2)= 0.0
            HL1(3)= 0.0
            HL1(4)= 0.5*(4.0*X1-2.+ZI)*(1.0+ZI)
            HL1(5)= 0.0
            HL1(6)= 0.0
            HL1(7)= 2.0*X2*(1.0-ZI)
            HL1(8)= 0.0
            HL1(9)= 2.0*X3*(1.0-ZI)
            HL1(10)=2.0*X2*(1.0+ZI)
            HL1(11)=0.0
            HL1(12)=2.0*X3*(1.0+ZI)
            HL1(13)=1.0-ZI**2
            HL1(14)=0.0
            HL1(15)=0.0
!  FOR L2-COORDINATE
            HL2(1)= 0.0
            HL2(2)= 0.5*(4.0*X2-2.-ZI)*(1.0-ZI)
            HL2(3)= 0.0
            HL2(4)= 0.0
            HL2(5)= 0.5*(4.0*X2-2.+ZI)*(1.0+ZI)
            HL2(6)= 0.0
            HL2(7)= 2.0*X1*(1.0-ZI)
            HL2(8)= 2.0*X3*(1.0-ZI)
            HL2(9)= 0.0
            HL2(10)=2.0*X1*(1.0+ZI)
            HL2(11)=2.0*X3*(1.0+ZI)
            HL2(12)=0.0
            HL2(13)=0.0
            HL2(14)=1.0-ZI**2
            HL2(15)=0.0
!  FOR L3-COORDINATE
            HL3(1)= 0.0
            HL3(2)= 0.0
            HL3(3)= 0.5*(4.0*X3-2.-ZI)*(1.0-ZI)
            HL3(4)= 0.0
            HL3(5)= 0.0
            HL3(6)= 0.5*(4.0*X3-2.+ZI)*(1.0+ZI)
            HL3(7)= 0.0
            HL3(8)= 2.0*X2*(1.0-ZI)
            HL3(9)= 2.0*X1*(1.0-ZI)
            HL3(10)=0.0
            HL3(11)=2.0*X2*(1.0+ZI)
            HL3(12)=2.0*X1*(1.0+ZI)
            HL3(13)=0.0
            HL3(14)=0.0
            HL3(15)=1.0-ZI**2
!  FOR Z-COORDINATE
            HZ(1)= 0.5*X1*(-2.0*X1+1.0+2.0*ZI)
            HZ(2)= 0.5*X2*(-2.0*X2+1.0+2.0*ZI)
            HZ(3)= 0.5*X3*(-2.0*X3+1.0+2.0*ZI)
            HZ(4)= 0.5*X1*( 2.0*X1-1.0+2.0*ZI)
            HZ(5)= 0.5*X2*( 2.0*X2-1.0+2.0*ZI)
            HZ(6)= 0.5*X3*( 2.0*X3-1.0+2.0*ZI)
            HZ(7)= -2.0*X1*X2
            HZ(8)= -2.0*X2*X3
            HZ(9)= -2.0*X1*X3
            HZ(10)= 2.0*X1*X2
            HZ(11)= 2.0*X2*X3
            HZ(12)= 2.0*X1*X3
            HZ(13)=-2.0*X1*ZI
            HZ(14)=-2.0*X2*ZI
            HZ(15)=-2.0*X3*ZI
!  JACOBI MATRIX 
            XJ11=0.0
            XJ21=0.0
            XJ31=0.0
            XJ12=0.0
            XJ22=0.0
            XJ32=0.0
            XJ13=0.0
            XJ23=0.0
            XJ33=0.0
            DO I=1,NN
              XJ11=XJ11+(HL1(I)-HL3(I))*XX(I)
              XJ21=XJ21+(HL2(I)-HL3(I))*XX(I)
              XJ31=XJ31+HZ(I)*XX(I)
              XJ12=XJ12+(HL1(I)-HL3(I))*YY(I)
              XJ22=XJ22+(HL2(I)-HL3(I))*YY(I)
              XJ32=XJ32+HZ(I)*YY(I)
              XJ13=XJ13+(HL1(I)-HL3(I))*ZZ(I)
              XJ23=XJ23+(HL2(I)-HL3(I))*ZZ(I)
              XJ33=XJ33+HZ(I)*ZZ(I)
            ENDDO
!DETERMINANT OF JACOBIAN
            DET=XJ11*XJ22*XJ33    &
              +XJ12*XJ23*XJ31     &
              +XJ13*XJ21*XJ32     &
              -XJ13*XJ22*XJ31     &
              -XJ12*XJ21*XJ33     &
              -XJ11*XJ23*XJ32
! INVERSION OF JACOBIAN
            DUM=1.0/DET
            XJI11=DUM*( XJ22*XJ33-XJ23*XJ32)
            XJI21=DUM*(-XJ21*XJ33+XJ23*XJ31)
            XJI31=DUM*( XJ21*XJ32-XJ22*XJ31)
            XJI12=DUM*(-XJ12*XJ33+XJ13*XJ32)
            XJI22=DUM*( XJ11*XJ33-XJ13*XJ31)
            XJI32=DUM*(-XJ11*XJ32+XJ12*XJ31)
            XJI13=DUM*( XJ12*XJ23-XJ13*XJ22)
            XJI23=DUM*(-XJ11*XJ23+XJ13*XJ21)
            XJI33=DUM*( XJ11*XJ22-XJ12*XJ21)
!* DUDX DUDY DUDZ ....
            DUDX=0.0
            DUDY=0.0
            DUDZ=0.0
            DVDX=0.0
            DVDY=0.0
            DVDZ=0.0
            DWDX=0.0
            DWDY=0.0
            DWDZ=0.0
            DO J=1,NN
              DNDX=XJI11*(HL1(J)-HL3(J)) +XJI12*(HL2(J)-HL3(J)) +XJI13*HZ(J)
              DNDY=XJI21*(HL1(J)-HL3(J)) +XJI22*(HL2(J)-HL3(J)) +XJI23*HZ(J)
              DNDZ=XJI31*(HL1(J)-HL3(J)) +XJI32*(HL2(J)-HL3(J)) +XJI33*HZ(J)
              DUDX=DUDX+DNDX*UX(J)
              DUDY=DUDY+DNDY*UX(J)
              DUDZ=DUDZ+DNDZ*UX(J)
              DVDX=DVDX+DNDX*UY(J)
              DVDY=DVDY+DNDY*UY(J)
              DVDZ=DVDZ+DNDZ*UY(J)
              DWDX=DWDX+DNDX*UZ(J)
              DWDY=DWDY+DNDY*UZ(J)
              DWDZ=DWDZ+DNDZ*UZ(J)
            ENDDO
!** Fill BMAT
            HEAD=1
            DO J=1,NN
!* DNDX DNDY DNDZ
              DNDX=XJI11*(HL1(J)-HL3(J)) +XJI12*(HL2(J)-HL3(J)) +XJI13*HZ(J)
              DNDY=XJI21*(HL1(J)-HL3(J)) +XJI22*(HL2(J)-HL3(J)) +XJI23*HZ(J)
              DNDZ=XJI31*(HL1(J)-HL3(J)) +XJI32*(HL2(J)-HL3(J)) +XJI33*HZ(J)
!** BL0
              BL0(1,HEAD )=DNDX
              BL0(2,HEAD )=0.0
              BL0(3,HEAD )=0.0
              BL0(4,HEAD )=DNDY
              BL0(5,HEAD )=0.0
              BL0(6,HEAD )=DNDZ
              BL0(1,HEAD+1)=0.0
              BL0(2,HEAD+1)=DNDY
              BL0(3,HEAD+1)=0.0
              BL0(4,HEAD+1)=DNDX
              BL0(5,HEAD+1)=DNDZ
              BL0(6,HEAD+1)=0.0
              BL0(1,HEAD+2)=0.0
              BL0(2,HEAD+2)=0.0
              BL0(3,HEAD+2)=DNDZ
              BL0(4,HEAD+2)=0.0
              BL0(5,HEAD+2)=DNDY
              BL0(6,HEAD+2)=DNDX
!** BL1
              BL1(1,HEAD )=DUDX*DNDX
              BL1(2,HEAD )=DUDY*DNDY
              BL1(3,HEAD )=DUDZ*DNDZ
              BL1(4,HEAD )=DUDY*DNDX+DUDX*DNDY
              BL1(5,HEAD )=DUDY*DNDZ+DUDZ*DNDY
              BL1(6,HEAD )=DUDZ*DNDX+DUDX*DNDZ
              BL1(1,HEAD+1)=DVDX*DNDX
              BL1(2,HEAD+1)=DVDY*DNDY
              BL1(3,HEAD+1)=DVDZ*DNDZ
              BL1(4,HEAD+1)=DVDY*DNDX+DVDX*DNDY
              BL1(5,HEAD+1)=DVDY*DNDZ+DVDZ*DNDY
              BL1(6,HEAD+1)=DVDZ*DNDX+DVDX*DNDZ
              BL1(1,HEAD+2)=DWDX*DNDX
              BL1(2,HEAD+2)=DWDY*DNDY
              BL1(3,HEAD+2)=DWDZ*DNDZ
              BL1(4,HEAD+2)=DWDY*DNDX+DWDX*DNDY
              BL1(5,HEAD+2)=DWDY*DNDZ+DWDZ*DNDY
              BL1(6,HEAD+2)=DWDZ*DNDX+DWDX*DNDZ
              HEAD=HEAD+NDOF
            ENDDO
!** DE & DS ***
            DO J=1,6
              DE(J)=0.0
              DS(J)=0.0
            ENDDO
            DO J=1,NN*NDOF
              DO K=1,6
                DE(K)=DE(K)+(BL0(K,J)+BL1(K,J))*DU(J)
              ENDDO
            ENDDO

            DO J=1,6
              STRAIN(J,L1,L2,LZ)=STRAIN(J,L1,L2,LZ)+DE(J)
            ENDDO
!** Update stress
!ogasawara till next ogasawara
!---  elastic analysis
      if(idx_elpl .eq. 0) then
            DO J=1,6
              DO K=1,6
                DS(K)=DS(K)+DMAT(K,J)*DE(J)
              ENDDO
            ENDDO
            DO J=1,6
              STRESS(J,L1,L2,LZ)=STRESS(J,L1,L2,LZ)+DS(J)
            ENDDO
!---  elasto-plastic analysis
!      yield judgement
!      stress integration
      else if(idx_elpl .ge. 1) then
        call fstr_elpl_INTG(DMAT,STRESS(:,L1,L2,LZ),STR_YLD(L1,L2,LZ)     &
                         ,FLG_ELPL(L1,L2,LZ), h_dash, EP_BAR(L1,L2,LZ)    &
                         ,EE, PP, DE, iter,BIJ,MAT_OTHO1,MAT_ISO)
      end if
!ogasawara till next ogasawara
!
          ENDDO
        ENDDO
      ENDDO
      RETURN
      END SUBROUTINE UPDATE_C3D15TL
!----------------------------------------------------------------------*
       SUBROUTINE UPDATE_C3D10TL(ILOG,XX,YY,ZZ,DISP,DU,EE,PP,STRESS                &
                           ,idx_elpl,h_dash,FLG_ELPL,STR_YLD ,STRAIN,EP_BAR        &
                           ,iter,i_step,idx_temp                                   &
                           ,TT,TTI,CURE3D_nod,CURE3D_nod_I,ccc_rl,ramd_chm         &
                           ,ccc_rl1,ccc_rl2                                        &
                           ,yeta,yeta1,yeta2,idx_visco,DT0                         &
                           ,STRAIN_V1,STRAIN_V2,count,icel,MID,ALPHA               &
                           ,STRESS_nd,STRAIN_nd,STRESS_ndP1,SURFF_nd,idx_stres     & 
                           ,MAT_ISO,BIJ,ALPHA_B,MAT_OTHO1)
!----------------------------------------------------------------------*
!
! UPDATE 3D 10 NODE PLANE ELEMENT
!
      use hecmw
      use gauss_integration 
      IMPLICIT NONE
! GAUSS POINT
      INTEGER(kind=kint) NG
      PARAMETER(NG=3)
! I/F VARIABLES
      INTEGER(kind=kint) ILOG,count,icel,MID
      INTEGER(kind=kint) idx_temp,idx_visco
      INTEGER(kind=kint) idx_stres
      REAL(kind=kreal) XX(*),YY(*),ZZ(*),TT(*),TTI(*),DISP(*),DU(*),EE,PP
      REAL(kind=kreal) STRESS(6,NG,NG,NG)
      INTEGER(kind=kint) iter, i_step
      INTEGER(kind=kint) idx_elpl,FLG_ELPL(NG,NG,NG)
      INTEGER(kind=kint) FLG_ELPL_nd(10)
      REAL(kind=kreal) h_dash,STR_YLD(NG,NG,NG),EP_BAR(NG,NG,NG),STRAIN(6,NG,NG,NG)
      REAL(kind=kreal) STRESS_nd(6,10),STRAIN_nd(6,10),STRESS_ndP1(6,10)
      REAL(kind=kreal) SURFF_nd(8,10) 
      REAL(kind=kreal) X_NODE(10),Y_NODE(10),Z_NODE(10)
      REAL(kind=kreal) STR_YLD_nd(10),EP_BAR_nd(10)
! LOCAL VARIABLES
      INTEGER(kind=kint) NN,NDOF
      PARAMETER(NN=10,NDOF=3)
      INTEGER(kind=kint) I,J,K,L,NUM
      INTEGER(kind=kint) L1,L2,L3
      INTEGER(kind=kint) HEAD
      REAL(kind=kreal) RI,SI,TI,RP,SP,TP,RM,SM,TM
      REAL(kind=kreal) H(NN),HL1(NN),HL2(NN),HL3(NN),HL4(NN)
      REAL(kind=kreal) X1,X2,X3,X4,XL1,XL2,XL3
      REAL(kind=kreal) XJ11,XJ21,XJ31,XJ12,XJ22,XJ32,XJ13,XJ23,XJ33,DET,WG,DUM
      REAL(kind=kreal) XJI11,XJI21,XJI31,XJI12,XJI22,XJI32,XJI13,XJI23,XJI33
      REAL(kind=kreal) DMAT(6,6),SMAT(9,9)
      REAL(kind=kreal) DNDX,DNDY,DNDZ
      REAL(kind=kreal) DUDX,DUDY,DUDZ
      REAL(kind=kreal) DVDX,DVDY,DVDZ
      REAL(kind=kreal) DWDX,DWDY,DWDZ
      REAL(kind=kreal) BL0(6,NN*NDOF),BL1(6,NN*NDOF)
      REAL(kind=kreal) UX(NN),UY(NN),UZ(NN)
      REAL(kind=kreal) S11,S22,S33,S23,S13,S12
        REAL(kind=kreal) SD(6),SMN 
        REAL(kind=kreal) DSD(6),DSMN 
        REAL(kind=kreal) E_V1(6),E_V2(6),ED_V1(6),EMN_V1,ED_V2(6),EMN_V2 
        REAL(kind=kreal) DE_V1(6),DE_V2(6),DED_V1(6),DEMN_V1,DED_V2(6),DEMN_V2 
      REAL(kind=kreal) DE(6),DS(6)
      REAL(kind=kreal) DE0(6)
      REAL(kind=kreal) TEMPC,TEMPI,ccc_rl(*)
      REAL(kind=kreal) CURE3D_nod(*),CURE3D_nod_I(*),CURE3D_itg,CURE3D_itg_I
      REAL(kind=kreal) ccc_rl1(*),ccc_rl2(*)
      REAL(kind=kreal) ramd_chm,ALPHA,ALP,ALP1,THERMAL_EPS,EPSTH(6)
      REAL(kind=kreal) yeta(*),yeta1(*),yeta2(*),DT0
      REAL(kind=kreal) STRAIN_V1(6,NG,NG,NG),STRAIN_V2(6,NG,NG,NG)
      REAL(kind=kreal) DSTR_A(6),EE1,EE2,PP1,PP2,YETA_G,YETA_K
        REAL(kind=kreal) YETA_G1,YETA_G2,YETA_K1,YETA_K2,T_G,T_G1,T_G2,T_K,T_K1,T_K2
      REAL(kind=kreal) AAA(6,6),BBB(6),BBB1(6),BBB2(6)
      REAL(kind=kreal) DMAT1(6,6),DMAT2(6,6)
        REAL(kind=kreal) DMAT_G(6,6),DEPS_AD(6),DEPS_AM  
        REAL(kind=kreal) C_G,C_K,C_G1,C_G2,C_K1,C_K2,C_GT,C_KT 
      REAL(kind=kreal) YEINV(6,6),YEINV1(6,6),YEINV2(6,6)
      REAL(kind=kreal) DETF,FMAT(3,3),FMATINV(3,3),PK2MAT(3,3),WK(3,3),TMAT(3,3)
      REAL(kind=kreal) PIMAT(3,3),NORML0(3),NORML1(3) 

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

      REAL(kind=kreal) THERMAL_EPS1,THERMAL_EPS2,THERMAL_EPS3 
      REAL(kind=kreal) EPS_G(3,3),EPS_GB(3,3) 
!* Total displacement
      DO I=1,NN
        UX(I)=DISP(NDOF*I-2)
        UY(I)=DISP(NDOF*I-1)
        UZ(I)=DISP(NDOF*I  )
      ENDDO
!!
      THERMAL_EPS = 0.0

      DO I=1,6
        DSTR_A(I) = 0.0
      END DO
!
      DO I=1,10
        FLG_ELPL_nd(I) = 0
        STR_YLD_nd (I) = 0.0
        EP_BAR_nd  (I) = 0.0
      END DO


! LOOP FOR INTEGRATION POINTS
      DO L3=1,NG
        XL3=XG(NG,L3)
        X3 =(XL3+1.0)*0.5
        DO L2=1,NG
          XL2=XG(NG,L2)
          X2 =(1.0-X3)*(XL2+1.0)*0.5
          DO L1=1,NG
            XL1=XG(NG,L1)
            X1=(1.0-X2-X3)*(XL1+1.0)*0.5
!  INTERPOLATION FUNCTION
            X4=1.0-X1-X2-X3
            H(1)=X1*(2.0*X1-1.0)
            H(2)=X2*(2.0*X2-1.0)
            H(3)=X3*(2.0*X3-1.0)
            H(4)=X4*(2.0*X4-1.0)
            H(5)=4.0*X1*X2
            H(6)=4.0*X2*X3
            H(7)=4.0*X1*X3
            H(8)=4.0*X1*X4
            H(9)=4.0*X2*X4
            H(10)=4.0*X3*X4
! DERIVATIVE OF INTERPOLATION FUNCTION
! FOR L1-COORDINATE
            HL1(1)= 4.0*X1-1.0
            HL1(2)= 0.0
            HL1(3)= 0.0
            HL1(4)= 0.0
            HL1(5)= 4.0*X2
            HL1(6)= 0.0
            HL1(7)= 4.0*X3
            HL1(8)= 4.0*X4
            HL1(9)= 0.0
            HL1(10)= 0.0
!  FOR L2-COORDINATE
            HL2(1)= 0.0 
            HL2(2)= 4.0*X2-1.0
            HL2(3)= 0.0
            HL2(4)= 0.0
            HL2(5)= 4.0*X1
            HL2(6)= 4.0*X3
            HL2(7)= 0.0
            HL2(8)= 0.0
            HL2(9)= 4.0*X4
            HL2(10)= 0.0
!  FOR L3-COORDINATE
            HL3(1)= 0.0
            HL3(2)= 0.0
            HL3(3)= 4.0*X3-1.0
            HL3(4)= 0.0
            HL3(5)= 0.0
            HL3(6)= 4.0*X2
            HL3(7)= 4.0*X1
            HL3(8)= 0.0
            HL3(9)= 0.0
            HL3(10)= 4.0*X4
!  FOR L4-COORDINATE
            HL4(1)= 0.0
            HL4(2)= 0.0
            HL4(3)= 0.0 
            HL4(4)= 4.0*X4-1.0
            HL4(5)= 0.0
            HL4(6)= 0.0
            HL4(7)= 0.0
            HL4(8)= 4.0*X1
            HL4(9)= 4.0*X2
            HL4(10)= 4.0*X3
!  JACOBI MATRIX 
            XJ11=0.0
            XJ21=0.0
            XJ31=0.0
            XJ12=0.0
            XJ22=0.0
            XJ32=0.0
            XJ13=0.0
            XJ23=0.0
            XJ33=0.0
            DO I=1,NN
              XJ11=XJ11+(HL1(I)-HL4(I))*XX(I)
              XJ21=XJ21+(HL2(I)-HL4(I))*XX(I)
              XJ31=XJ31+(HL3(I)-HL4(I))*XX(I)
              XJ12=XJ12+(HL1(I)-HL4(I))*YY(I)
              XJ22=XJ22+(HL2(I)-HL4(I))*YY(I)
              XJ32=XJ32+(HL3(I)-HL4(I))*YY(I)
              XJ13=XJ13+(HL1(I)-HL4(I))*ZZ(I)
              XJ23=XJ23+(HL2(I)-HL4(I))*ZZ(I)
              XJ33=XJ33+(HL3(I)-HL4(I))*ZZ(I)
            ENDDO
!DETERMINANT OF JACOBIAN
            DET=XJ11*XJ22*XJ33        &
                 +XJ12*XJ23*XJ31      &
                 +XJ13*XJ21*XJ32      &
                 -XJ13*XJ22*XJ31      &
                 -XJ12*XJ21*XJ33      &
                 -XJ11*XJ23*XJ32
! INVERSION OF JACOBIAN
            DUM=1.0/DET
            XJI11=DUM*( XJ22*XJ33-XJ23*XJ32)
            XJI21=DUM*(-XJ21*XJ33+XJ23*XJ31)
            XJI31=DUM*( XJ21*XJ32-XJ22*XJ31)
            XJI12=DUM*(-XJ12*XJ33+XJ13*XJ32)
            XJI22=DUM*( XJ11*XJ33-XJ13*XJ31)
            XJI32=DUM*(-XJ11*XJ32+XJ12*XJ31)
            XJI13=DUM*( XJ12*XJ23-XJ13*XJ22)
            XJI23=DUM*(-XJ11*XJ23+XJ13*XJ21)
            XJI33=DUM*( XJ11*XJ22-XJ12*XJ21)
!* DUDX DUDY DUDZ ....
            DUDX=0.0
            DUDY=0.0
            DUDZ=0.0
            DVDX=0.0
            DVDY=0.0
            DVDZ=0.0
            DWDX=0.0
            DWDY=0.0
            DWDZ=0.0
            DO J=1,NN
              DNDX=XJI11*(HL1(J)-HL4(J)) +XJI12*(HL2(J)-HL4(J)) +XJI13*(HL3(J)-HL4(J))
              DNDY=XJI21*(HL1(J)-HL4(J)) +XJI22*(HL2(J)-HL4(J)) +XJI23*(HL3(J)-HL4(J))
              DNDZ=XJI31*(HL1(J)-HL4(J)) +XJI32*(HL2(J)-HL4(J)) +XJI33*(HL3(J)-HL4(J))
              DUDX=DUDX+DNDX*UX(J)
              DUDY=DUDY+DNDY*UX(J)
              DUDZ=DUDZ+DNDZ*UX(J)
              DVDX=DVDX+DNDX*UY(J)
              DVDY=DVDY+DNDY*UY(J)
              DVDZ=DVDZ+DNDZ*UY(J)
              DWDX=DWDX+DNDX*UZ(J)
              DWDY=DWDY+DNDY*UZ(J)
              DWDZ=DWDZ+DNDZ*UZ(J)
            ENDDO
!** Fill BMAT
            HEAD=1
            DO J=1,NN
!* DNDX DNDY DNDZ
              DNDX=XJI11*(HL1(J)-HL4(J)) +XJI12*(HL2(J)-HL4(J)) +XJI13*(HL3(J)-HL4(J))
              DNDY=XJI21*(HL1(J)-HL4(J)) +XJI22*(HL2(J)-HL4(J)) +XJI23*(HL3(J)-HL4(J))
              DNDZ=XJI31*(HL1(J)-HL4(J)) +XJI32*(HL2(J)-HL4(J)) +XJI33*(HL3(J)-HL4(J))
!** BL0
              BL0(1,HEAD )=DNDX
              BL0(2,HEAD )=0.0
              BL0(3,HEAD )=0.0
              BL0(4,HEAD )=DNDY
              BL0(5,HEAD )=0.0
              BL0(6,HEAD )=DNDZ
              BL0(1,HEAD+1)=0.0
              BL0(2,HEAD+1)=DNDY
              BL0(3,HEAD+1)=0.0
              BL0(4,HEAD+1)=DNDX
              BL0(5,HEAD+1)=DNDZ
              BL0(6,HEAD+1)=0.0
              BL0(1,HEAD+2)=0.0
              BL0(2,HEAD+2)=0.0
              BL0(3,HEAD+2)=DNDZ
              BL0(4,HEAD+2)=0.0
              BL0(5,HEAD+2)=DNDY
              BL0(6,HEAD+2)=DNDX
!** BL1
              BL1(1,HEAD )=DUDX*DNDX
              BL1(2,HEAD )=DUDY*DNDY
              BL1(3,HEAD )=DUDZ*DNDZ
              BL1(4,HEAD )=DUDY*DNDX+DUDX*DNDY
              BL1(5,HEAD )=DUDY*DNDZ+DUDZ*DNDY
              BL1(6,HEAD )=DUDZ*DNDX+DUDX*DNDZ
              BL1(1,HEAD+1)=DVDX*DNDX
              BL1(2,HEAD+1)=DVDY*DNDY
              BL1(3,HEAD+1)=DVDZ*DNDZ
              BL1(4,HEAD+1)=DVDY*DNDX+DVDX*DNDY
              BL1(5,HEAD+1)=DVDY*DNDZ+DVDZ*DNDY
              BL1(6,HEAD+1)=DVDZ*DNDX+DVDX*DNDZ
              BL1(1,HEAD+2)=DWDX*DNDX
              BL1(2,HEAD+2)=DWDY*DNDY
              BL1(3,HEAD+2)=DWDZ*DNDZ
              BL1(4,HEAD+2)=DWDY*DNDX+DWDX*DNDY
              BL1(5,HEAD+2)=DWDY*DNDZ+DWDZ*DNDY
              BL1(6,HEAD+2)=DWDZ*DNDX+DWDX*DNDZ
              HEAD=HEAD+NDOF
            ENDDO

!** DE & DS ***
!!
   if(idx_temp .eq. 1) then
!
      TEMPC=0.0
      TEMPI=0.0
      CURE3D_itg  =0.0
      CURE3D_itg_I=0.0
      DO J=1,NN
        TEMPC=TEMPC+TT (J)*H(J)
        TEMPI=TEMPI+TTI(J)*H(J)
        CURE3D_itg  =CURE3D_itg  +CURE3D_nod  (J)*H(J)
        CURE3D_itg_I=CURE3D_itg_I+CURE3D_nod_I(J)*H(J)
      END DO  
!
              if(CURE3D_itg_I .lt. 0.0) then
                if(iter .eq. -200) then
                  write(ILOG,*) 'stop : CURE3D_itg_I .lt. 0.0 '
                  write(ILOG,*) 'elem = ',icel,'   CURE3D_itg_I= ',CURE3D_itg_I
                  call flush(ILOG)
!                 call hecmw_abort( hecmw_comm_get_comm())
                end if
               CURE3D_itg_I = 0.0
              else if(CURE3D_itg_I .gt. 1.0) then
                if(iter .eq. -200) then
                  write(ILOG,*) 'stop : CURE3D_itg_I .gt. 0.0 '
                  write(ILOG,*) 'elem = ',icel,'   CURE3D_itg_I= ',CURE3D_itg_I
                  call flush(ILOG)
!                 call hecmw_abort( hecmw_comm_get_comm())
                end if
               CURE3D_itg_I = 1.0
              end if
!
              if(CURE3D_itg .lt. 0.0) then
               CURE3D_itg = 0.0
              else if(CURE3D_itg .gt. 1.0) then
               CURE3D_itg = 1.0
              end if
!
!
    if(MID .ne. 2) then
      if(MAT_ISO .eq. 0) then
        ALP = ALPHA*(TEMPC-TEMPI) 
        THERMAL_EPS1 = ALP
        THERMAL_EPS2 = ALP
        THERMAL_EPS3 = ALP
      else if(MAT_ISO .eq. 1) then
        THERMAL_EPS1 = ALPHA_B(1)*(TEMPC-TEMPI)
        THERMAL_EPS2 = ALPHA_B(2)*(TEMPC-TEMPI)
        THERMAL_EPS3 = ALPHA_B(3)*(TEMPC-TEMPI)
      end if
    else if(MID .eq. 2) then
      EE  = ccc_rl(1)*dexp(-ccc_rl(2)*TEMPI)*ccc_rl(3)*dexp(ccc_rl(4)*CURE3D_itg_I)
!
      ALP =   (ccc_rl(5)*TEMPI+ccc_rl(6))*(1.0-CURE3D_itg_I)     &
            + (ccc_rl(7)*TEMPI+ccc_rl(8))*CURE3D_itg_I 
      ALP1 = ALP
      ALP = ALP*(TEMPC-TEMPI) - ramd_chm*(CURE3D_itg-CURE3D_itg_I)
      THERMAL_EPS1 = ALP/3.0
      THERMAL_EPS2 = ALP/3.0
      THERMAL_EPS3 = ALP/3.0
!!
     end if
    end if
!!
      call fstr_elpl_DMAT(DMAT,STRESS(:,L1,L2,L3),FLG_ELPL(L1,L2,L3), h_dash, EE, PP  &
                         ,BIJ,MAT_OTHO1,MAT_ISO)
!!
!!!  20110805  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!! Add Q force
!!   DMAT_G, dSa by visco 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
            if(idx_visco .eq. 0 .or. idx_visco .eq. 1 ) then
              DO J=1,6
                DSTR_A(J) = 0.0
              END DO

            else if(idx_visco .eq. 2) then
              if(idx_temp .eq. 1 .and. MID .ne. 2) then
                DO J=1,6
                 DSTR_A(J) = 0.0
                END DO
              else
                DO J=1,6
                 DSTR_A(J) = 0.0
                END DO

                 T_G    = ccc_rl(9)
                 C_G    = 1.0 / (EE/2.0/(1.0+PP))  
                 YETA_G = 1.0 / C_G * T_G 
                 T_K    = ccc_rl(10)
                 C_K    = 1.0 / (EE/3.0/(1.0-2.0*PP)) 
                 YETA_K = 1.0 / C_K * T_K 

                 S11=STRESS(1,L1,L2,L3)
                 S22=STRESS(2,L1,L2,L3)
                 S33=STRESS(3,L1,L2,L3)
                 S12=STRESS(4,L1,L2,L3)
                 S23=STRESS(5,L1,L2,L3)
                 S13=STRESS(6,L1,L2,L3)      

                 SMN = (S11+S22+S33)/3.0 
                 SD(1) = S11 - SMN 
                 SD(2) = S22 - SMN 
                 SD(3) = S33 - SMN 
                 SD(4) = S12 
                 SD(5) = S23 
                 SD(6) = S13 

                 C_GT = C_G + 1.0/2.0/YETA_G * DT0 
                 C_KT = C_K + 1.0/2.0/YETA_K * DT0 

                 DO J=1,6
                  DEPS_AD(J) = (1.0/YETA_G)*SD(J) 
                 END DO

                 DEPS_AM = (1.0/YETA_K)*SMN 
!
               if(C_GT .eq. 0.0D+00) then
                 write(ILOG,*) 'stop : C_GT = 0.0'
                 write(ILOG,*) 'elem = ',icel,'   C_GT ',C_GT 
                 call flush(ILOG)
                 call hecmw_abort( hecmw_comm_get_comm())
               end if
               if(C_KT .eq. 0.0D+00) then
                 write(ILOG,*) 'stop : C_KT = 0.0'
                 write(ILOG,*) 'elem = ',icel,'   C_KT ',C_KT 
                 call flush(ILOG)
                 call hecmw_abort( hecmw_comm_get_comm())
               end if
!!
                 DO J=1,6
                  DO K=1,6
                   DMAT(J,K) = 0.0
                  END DO
                 END DO

                 DMAT(1,1) = 1.0/C_KT + 4.0/3.0/C_GT 
                 DMAT(1,2) = 1.0/C_KT - 2.0/3.0/C_GT 
                 DMAT(1,3) = 1.0/C_KT - 2.0/3.0/C_GT 
                 DMAT(2,1) = DMAT(1,2)
                 DMAT(2,2) = 1.0/C_KT + 4.0/3.0/C_GT 
                 DMAT(2,3) = 1.0/C_KT - 2.0/3.0/C_GT 
                 DMAT(3,1) = DMAT(1,3)
                 DMAT(3,2) = DMAT(2,3)
                 DMAT(3,3) = 1.0/C_KT + 4.0/3.0/C_GT 
                 DMAT(4,4) =          1.0/  C_GT 
                 DMAT(5,5) =          1.0/  C_GT 
                 DMAT(6,6) =          1.0/  C_GT 
!
!!!! 20111009  !!!!!!!!!
                  DO J=1,6
                   DSTR_A(J) = DEPS_AD(J)/C_GT 
                  END DO
                  DO J=1,3
                   DSTR_A(J) = DSTR_A(J) + DEPS_AM/C_KT 
                  END DO

                if(i_step .eq. 1 .and. iter .le. 1) then
                  DO J=1,6
                    DSTR_A(J) = DSTR_A(J)*DT0
                  END DO
                else if(i_step .ge. 2 .and. iter .le. 1) then
                  DO J=1,6
                    DSTR_A(J) = DSTR_A(J)*DT0
                  END DO
                else
                    DSTR_A(J) = DSTR_A(J)*DT0
                end if
              end if
            else if(idx_visco .eq. 3) then
            else if(idx_visco .eq. 4) then
            else if(idx_visco .eq. 5) then
            else if(idx_visco .eq. 6) then
              if(idx_temp .eq. 1 .and. MID .ne. 2) then
                DO J=1,6
                 DSTR_A(J) = 0.0
                END DO

              else
                DO J=1,6
                 DSTR_A(J) = 0.0
                END DO

                 T_G    = ccc_rl(9)
                 C_G    = 1.0 / (EE/2.0/(1.0+PP))  
                 YETA_G = 1.0 / C_G * T_G 
                 T_K    = ccc_rl(10)
                 C_K    = 1.0 / (EE/3.0/(1.0-2.0*PP)) 
                 YETA_K = 1.0 / C_K * T_K 

                 T_G1    = T_G
                 C_G1    = C_G 
                 YETA_G1 = YETA_G 
                 T_K1    = T_K
                 C_K1    = C_K 
                 YETA_K1 = YETA_K 

                 T_G2    = T_G
                 C_G2    = C_G 
                 YETA_G2 = YETA_G 
                 T_K2    = T_K
                 C_K2    = C_K 
                 YETA_K2 = YETA_K 

                 S11=STRESS(1,L1,L2,L3)
                 S22=STRESS(2,L1,L2,L3)
                 S33=STRESS(3,L1,L2,L3)
                 S12=STRESS(4,L1,L2,L3)
                 S23=STRESS(5,L1,L2,L3)
                 S13=STRESS(6,L1,L2,L3)      

                 SMN = (S11+S22+S33)/3.0 
                 SD(1) = S11 - SMN 
                 SD(2) = S22 - SMN 
                 SD(3) = S33 - SMN 
                 SD(4) = S12 
                 SD(5) = S23 
                 SD(6) = S13 

                 DO J=1,6
                   E_V1(J) = STRAIN_V1(J,L1,L2,L3) 
                   E_V2(J) = STRAIN_V2(J,L1,L2,L3) 
                 END DO
                   EMN_V1 = (E_V1(1)+E_V1(2)+E_V1(3))/3.0 
                   EMN_V2 = (E_V2(1)+E_V2(2)+E_V2(3))/3.0 
                 DO J=1,3
                   ED_V1(  J) = E_V1(  J) - EMN_V1 
                   ED_V1(3+J) = E_V1(3+J) 
                   ED_V2(  J) = E_V2(  J) - EMN_V2 
                   ED_V2(3+J) = E_V2(3+J) 
                 END DO
                 
                 C_GT = C_G + 1.0/2.0/YETA_G * DT0 + C_G1*(1.0-(1.0-dexp(-DT0/T_G1))*T_G1/DT0)   & 
                                                   + C_G2*(1.0-(1.0-dexp(-DT0/T_G2))*T_G2/DT0) 
                 C_KT = C_K + 1.0/2.0/YETA_K * DT0 + C_K1*(1.0-(1.0-dexp(-DT0/T_K1))*T_K1/DT0)   & 
                                                   + C_K2*(1.0-(1.0-dexp(-DT0/T_K2))*T_K2/DT0) 
                 DO J=1,6
                  DEPS_AD(J) = (1.0/YETA_G + 1.0/YETA_G1*(1.0-dexp(-DT0/T_G1))*T_G1/DT0            & 
                                           + 1.0/YETA_G2*(1.0-dexp(-DT0/T_G2))*T_G2/DT0)*SD(J)     & 
                                           -(2.0/T_G1   *(1.0-dexp(-DT0/T_G1))*T_G1/DT0)*ED_V1(J)  & 
                                           -(2.0/T_G2   *(1.0-dexp(-DT0/T_G2))*T_G2/DT0)*ED_V2(J)   
                 END DO

                 DEPS_AM = (1.0/YETA_K + 1.0/YETA_K1*(1.0-dexp(-DT0/T_K1))*T_K1/DT0          & 
                                       + 1.0/YETA_K2*(1.0-dexp(-DT0/T_K2))*T_K2/DT0)*SMN     & 
                                       -(3.0/T_K1   *(1.0-dexp(-DT0/T_K1))*T_K1/DT0)*EMN_V1  & 
                                       -(3.0/T_K2   *(1.0-dexp(-DT0/T_K2))*T_K2/DT0)*EMN_V2   
!
               if(C_GT .eq. 0.0D+00) then
                 write(ILOG,*) 'stop : C_GT = 0.0'
                 write(ILOG,*) 'elem = ',icel,'   C_GT ',C_GT 
                 call flush(ILOG)
                 call hecmw_abort( hecmw_comm_get_comm())
               end if
               if(C_KT .eq. 0.0D+00) then
                 write(ILOG,*) 'stop : C_KT = 0.0'
                 write(ILOG,*) 'elem = ',icel,'   C_KT ',C_KT 
                 call flush(ILOG)
                 call hecmw_abort( hecmw_comm_get_comm())
               end if
!!
                 DO J=1,6
                  DO K=1,6
                   DMAT(J,K) = 0.0
                  END DO
                 END DO

                 DMAT(1,1) = 1.0/C_KT + 4.0/3.0/C_GT 
                 DMAT(1,2) = 1.0/C_KT - 2.0/3.0/C_GT 
                 DMAT(1,3) = 1.0/C_KT - 2.0/3.0/C_GT 
                 DMAT(2,1) = DMAT(1,2)
                 DMAT(2,2) = 1.0/C_KT + 4.0/3.0/C_GT 
                 DMAT(2,3) = 1.0/C_KT - 2.0/3.0/C_GT 
                 DMAT(3,1) = DMAT(1,3)
                 DMAT(3,2) = DMAT(2,3)
                 DMAT(3,3) = 1.0/C_KT + 4.0/3.0/C_GT 
                 DMAT(4,4) =          1.0/  C_GT 
                 DMAT(5,5) =          1.0/  C_GT 
                 DMAT(6,6) =          1.0/  C_GT 
!
!!!! 20111009  !!!!!!!!!
                  DO J=1,6
                   DSTR_A(J) = DEPS_AD(J)/C_GT 
                  END DO
                  DO J=1,3
                   DSTR_A(J) = DSTR_A(J) + DEPS_AM/C_KT 
                  END DO

                if(i_step .eq. 1 .and. iter .le. 1) then
                  DO J=1,6
                    DSTR_A(J) = DSTR_A(J)*DT0
                  END DO
                else if(i_step .ge. 2 .and. iter .le. 1) then
                  DO J=1,6
                    DSTR_A(J) = DSTR_A(J)*DT0
                  END DO
                else
                    DSTR_A(J) = DSTR_A(J)*DT0
                end if
              end if
!
            end if
!!!  20110805  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!
      if(idx_elpl .eq. -1 .or. idx_elpl .eq. 2) then
        DO J=1,6
          DO K=1,NN*NDOF
            BL1(J,K)=0.0
          END DO
        END DO
      end if
!** DE & DS ***
!! Update strain & stress 
!! Update strain
            DO J=1,6
              DE(J)=0.0
              DS(J)=0.0
              EPSTH(J)=0.0
            ENDDO

            DO J=1,NN*NDOF
              DO K=1,6
                DE(K)=DE(K)+(BL0(K,J)+BL1(K,J))*DU(J)
              ENDDO
            ENDDO

            DO K=1,6
              DE0(K)=DE(K)
            ENDDO
!
            if(iter .eq. 1) then
!
              DO I=1,3
               DO J=1,3
                EPS_G (I,J) = 0.0
                EPS_GB(I,J) = 0.0
               END DO
              END DO

              EPS_GB(1,1) = THERMAL_EPS1 
              EPS_GB(2,2) = THERMAL_EPS2 
              EPS_GB(3,3) = THERMAL_EPS3 

            if(MAT_ISO .eq. 0) then
              DO I=1,3
               DO J=1,3
                  EPS_G(I,J) = EPS_GB(I,J) 
               END DO
              END DO
            else if(MAT_ISO .eq. 1) then
              DO I=1,3
               DO J=1,3
                DO K=1,3
                 DO L=1,3
                  EPS_G(I,J) = EPS_G(I,J) + BIJ(I,K)*BIJ(J,L)*EPS_GB(K,L)
                 END DO
                END DO
               END DO
              END DO
            end if

              EPSTH(1)=EPS_G(1,1)
              EPSTH(2)=EPS_G(2,2)
              EPSTH(3)=EPS_G(3,3)
              EPSTH(4)=EPS_G(1,2)
              EPSTH(5)=EPS_G(2,3)
              EPSTH(6)=EPS_G(3,1)
!
              DO J=1,6
               DE(J)=DE(J)-EPSTH(J)
              ENDDO
            end if
!!!!!!!!!!!
            DO J=1,6
              STRAIN(J,L1,L2,L3)=STRAIN(J,L1,L2,L3)+DE0(J)
            ENDDO
!
!! Update stress
!---  elastic analysis
      if(idx_elpl .le. 0) then
            DO J=1,6
              DO K=1,6
                DS(K)=DS(K)+DMAT(K,J)*DE(J)
              ENDDO
            ENDDO
!
!!!! 20111009  !!!!!!!!!
         if(i_step .eq. 1 .and. iter .le. 1) then
            DO J=1,6
              DS(J) = DS(J) - DSTR_A(J)
            END DO
         else if(i_step .ge. 2 .and. iter .le. 1) then
            DO J=1,6
              DS(J) = DS(J) - DSTR_A(J)
            END DO
         end if
!!
         if(iter .gt. 0) then
            DO J=1,6
              STRESS(J,L1,L2,L3)=STRESS(J,L1,L2,L3)+DS(J)
            ENDDO
         end if
!---  elasto-plastic analysis
!      yield judgement
!      stress integration
      else if(idx_elpl .ge. 1) then
        call fstr_elpl_INTG(DMAT,STRESS(:,L1,L2,L3),STR_YLD(L1,L2,L3)     &
                         ,FLG_ELPL(L1,L2,L3), h_dash, EP_BAR(L1,L2,L3)    &
                         ,EE, PP, DE, iter,BIJ,MAT_OTHO1,MAT_ISO)
      end if
!! Update Voigt element strain
            if(idx_visco .eq. 0 .or. idx_visco .eq. 1 ) then
            else if(idx_visco .eq. 2) then
            else if(idx_visco .eq. 3) then
            else if(idx_visco .eq. 4) then
            else if(idx_visco .eq. 5) then
            else if(idx_visco .eq. 6) then
              if(idx_temp .eq. 1 .and. MID .ne. 2) then

              else
                 DSMN = (DS(1)+DS(2)+DS(3))/3.0 
                 DSD(1) = DS(1) - DSMN 
                 DSD(2) = DS(2) - DSMN 
                 DSD(3) = DS(3) - DSMN 
                 DSD(4) = DS(4) 
                 DSD(5) = DS(5) 
                 DSD(6) = DS(6) 

                 DEMN_V1 = C_K1/3.0*(1.0-dexp(-DT0/T_K1))*SMN    & 
                         + C_K1/3.0*(1.0-(1.0-dexp(-DT0/T_K1))*T_K1/DT0)*DSMN - (1.0-dexp(-DT0/T_K1))*EMN_V1 
                 DEMN_V2 = C_K2/3.0*(1.0-dexp(-DT0/T_K2))*SMN    & 
                         + C_K2/3.0*(1.0-(1.0-dexp(-DT0/T_K2))*T_K2/DT0)*DSMN - (1.0-dexp(-DT0/T_K2))*EMN_V2 
                 DO J=1,6
                   DED_V1(J) = C_G1/2.0*(1.0-dexp(-DT0/T_G1))*SD(J)    & 
                             + C_G1/2.0*(1.0-(1.0-dexp(-DT0/T_G1))*T_G1/DT0)*DSD(J) - (1.0-dexp(-DT0/T_G1))*ED_V1(J) 
                   DED_V2(J) = C_G2/2.0*(1.0-dexp(-DT0/T_G2))*SD(J)    & 
                             + C_G2/2.0*(1.0-(1.0-dexp(-DT0/T_G2))*T_G2/DT0)*DSD(J) - (1.0-dexp(-DT0/T_G2))*ED_V2(J) 
                 END DO

                 DO J=1,3
                   STRAIN_V1(  J,L1,L2,L3) = STRAIN_V1(  J,L1,L2,L3) +    (DED_V1(  J) + DEMN_V1) 
                   STRAIN_V1(3+J,L1,L2,L3) = STRAIN_V1(3+J,L1,L2,L3) + 2.0*DED_V1(3+J) 
                   STRAIN_V2(  J,L1,L2,L3) = STRAIN_V2(  J,L1,L2,L3) +    (DED_V2(  J) + DEMN_V2) 
                   STRAIN_V2(3+J,L1,L2,L3) = STRAIN_V2(3+J,L1,L2,L3) + 2.0*DED_V2(3+J) 
                 END DO
              end if
            end if
!
          ENDDO
        ENDDO
      ENDDO

      RETURN
      END SUBROUTINE UPDATE_C3D10TL
!----------------------------------------------------------------------*
      SUBROUTINE NRCV_C3D20TL(ILOG,XX,YY,ZZ,TT,T0,EE,PP,ALP,DISP,STRESS,ARRAY)
!----------------------------------------------------------------------*
!
! NODAL RECOVERY FOR 3D 20 NODE ELEMENT /TOTAL LAGRANGE
!
      use hecmw
      use gauss_integration 
      IMPLICIT NONE
! GAUSS POINT
      INTEGER(kind=kint) NG
      PARAMETER(NG=3)
! I/F VARIABLES
      INTEGER(kind=kint) ILOG
      REAL(kind=kreal) XX(*),YY(*),ZZ(*),TT(*),T0(*),DISP(*),ARRAY(13,*)
      REAL(kind=kreal) STRESS(6,NG,NG,NG)
      REAL(kind=kreal) EE,PP,ALP
! LOCAL VARIABLES
      INTEGER(kind=kint) NN,NDOF
      PARAMETER(NN=20,NDOF=3)
      INTEGER(kind=kint) I,J,K,L
      INTEGER(kind=kint) LX,LY,LZ
      REAL(kind=kreal) RI,SI,TI,RP,SP,TP,RM,SM,TM
      REAL(kind=kreal) H(NN),HR(NN),HS(NN),HT(NN)
      REAL(kind=kreal) XJ11,XJ21,XJ31,XJ12,XJ22,XJ32,XJ13,XJ23,XJ33,DET,WG,DUM
      REAL(kind=kreal) XJI11,XJI21,XJI31,XJI12,XJI22,XJI32,XJI13,XJI23,XJI33
      REAL(kind=kreal) DNDX,DNDY,DNDZ
      REAL(kind=kreal) DUDX,DUDY,DUDZ
      REAL(kind=kreal) DVDX,DVDY,DVDZ
      REAL(kind=kreal) DWDX,DWDY,DWDZ
      REAL(kind=kreal) BL0(6,NN*NDOF),BL1(6,NN*NDOF)
      REAL(kind=kreal) UX(NN),UY(NN),UZ(NN)
      REAL(kind=kreal) S11,S22,S33,S23,S13,S12
      REAL(kind=kreal) DETF,FMAT(3,3),SMAT(3,3),WK(3,3),TMAT(3,3),STRAIN_L(6)
      REAL(kind=kreal) SGM(6),EPS(6),PS,SMISES
      REAL(kind=kreal) ARRAY_TEMP(13,NG*NG*NG)
      INTEGER(kind=kint) IC
!*************************
!  EVALUATION POINTS
!*************************
      REAL(kind=kreal) HG(NG*NG*NG),RX(NN),SX(NN),TX(NN)
      DATA RX/-1., 1., 1.,-1.,-1., 1., 1.,-1.,  &
               0., 1., 0.,-1., 0., 1., 0.,-1.,  &
              -1., 1., 1.,-1./
      DATA SX/-1.,-1., 1., 1.,-1.,-1., 1., 1.,  &
              -1., 0., 1., 0.,-1., 0., 1., 0.,  &
              -1.,-1., 1., 1./
      DATA TX/-1.,-1.,-1.,-1., 1., 1., 1., 1.,  &
              -1.,-1.,-1.,-1., 1.,1.,1.,1.,     &
               0., 0., 0., 0./
!* Total displacement 
      DO I=1,NN
        UX(I)=DISP(NDOF*I-2)
        UY(I)=DISP(NDOF*I-1)
        UZ(I)=DISP(NDOF*I  )
      ENDDO
      IC=0
! LOOP FOR INTEGRATION POINTS      
      DO LX=1,NG
        RI=XG(NG,LX)
        DO LY=1,NG
          SI=XG(NG,LY)
          DO LZ=1,NG
            TI=XG(NG,LZ)
            RP=1.0+RI
            SP=1.0+SI
            TP=1.0+TI
            RM=1.0-RI
            SM=1.0-SI
            TM=1.0-TI
!  INTERPOLATION FUNCTION
            H(1)=-0.125*RM*SM*TM*(2.0+RI+SI+TI)
            H(2)=-0.125*RP*SM*TM*(2.0-RI+SI+TI)
            H(3)=-0.125*RP*SP*TM*(2.0-RI-SI+TI)
            H(4)=-0.125*RM*SP*TM*(2.0+RI-SI+TI)
            H(5)=-0.125*RM*SM*TP*(2.0+RI+SI-TI)
            H(6)=-0.125*RP*SM*TP*(2.0-RI+SI-TI)
            H(7)=-0.125*RP*SP*TP*(2.0-RI-SI-TI)
            H(8)=-0.125*RM*SP*TP*(2.0+RI-SI-TI)
            H(9)=0.25*(1.0-RI**2)*SM*TM
            H(10)=0.25*RP*(1.0-SI**2)*TM
            H(11)=0.25*(1.0-RI**2)*SP*TM
            H(12)=0.25*RM*(1.0-SI**2)*TM
            H(13)=0.25*(1.0-RI**2)*SM*TP
            H(14)=0.25*RP*(1.0-SI**2)*TP
            H(15)=0.25*(1.0-RI**2)*SP*TP
            H(16)=0.25*RM*(1.0-SI**2)*TP
            H(17)=0.25*RM*SM*(1.0-TI**2)
            H(18)=0.25*RP*SM*(1.0-TI**2)
            H(19)=0.25*RP*SP*(1.0-TI**2)
            H(20)=0.25*RM*SP*(1.0-TI**2)
!  DERIVATIVE OF INTERPOLATION FUNCTION
!  FOR R-COORDINATE
            HR(1)=-0.125*RM*SM*TM+0.125*SM*TM*(2.0+RI+SI+TI)
            HR(2)=+0.125*RP*SM*TM-0.125*SM*TM*(2.0-RI+SI+TI)
            HR(3)=+0.125*RP*SP*TM-0.125*SP*TM*(2.0-RI-SI+TI)
            HR(4)=-0.125*RM*SP*TM+0.125*SP*TM*(2.0+RI-SI+TI)
            HR(5)=-0.125*RM*SM*TP+0.125*SM*TP*(2.0+RI+SI-TI)
            HR(6)=+0.125*RP*SM*TP-0.125*SM*TP*(2.0-RI+SI-TI)
            HR(7)=+0.125*RP*SP*TP-0.125*SP*TP*(2.0-RI-SI-TI)
            HR(8)=-0.125*RM*SP*TP+0.125*SP*TP*(2.0+RI-SI-TI)
            HR(9 )=-0.50*RI*SM*TM
            HR(10)=+0.25*(1.0-SI**2)*TM
            HR(11)=-0.50*RI*SP*TM
            HR(12)=-0.25*(1.0-SI**2)*TM
            HR(13)=-0.50*RI*SM*TP
            HR(14)=+0.25*(1.0-SI**2)*TP
            HR(15)=-0.50*RI*SP*TP
            HR(16)=-0.25*(1.0-SI**2)*TP
            HR(17)=-0.25*SM*(1.0-TI**2)
            HR(18)=+0.25*SM*(1.0-TI**2)
            HR(19)=+0.25*SP*(1.0-TI**2)
            HR(20)=-0.25*SP*(1.0-TI**2)
!  FOR S-COORDINATE
            HS(1)=-0.125*RM*SM*TM+0.125*RM*TM*(2.0+RI+SI+TI)
            HS(2)=-0.125*RP*SM*TM+0.125*RP*TM*(2.0-RI+SI+TI)
            HS(3)=+0.125*RP*SP*TM-0.125*RP*TM*(2.0-RI-SI+TI)
            HS(4)=+0.125*RM*SP*TM-0.125*RM*TM*(2.0+RI-SI+TI)
            HS(5)=-0.125*RM*SM*TP+0.125*RM*TP*(2.0+RI+SI-TI)
            HS(6)=-0.125*RP*SM*TP+0.125*RP*TP*(2.0-RI+SI-TI)
            HS(7)=+0.125*RP*SP*TP-0.125*RP*TP*(2.0-RI-SI-TI)
            HS(8)=+0.125*RM*SP*TP-0.125*RM*TP*(2.0+RI-SI-TI)
            HS(9)=-0.25*(1.0-RI**2)*TM
            HS(10)=-0.50*RP*SI*TM
            HS(11)=+0.25*(1.0-RI**2)*TM
            HS(12)=-0.50*RM*SI*TM
            HS(13)=-0.25*(1.0-RI**2)*TP
            HS(14)=-0.50*RP*SI*TP
            HS(15)=+0.25*(1.0-RI**2)*TP
            HS(16)=-0.50*RM*SI*TP
            HS(17)=-0.25*RM*(1.0-TI**2)
            HS(18)=-0.25*RP*(1.0-TI**2)
            HS(19)=+0.25*RP*(1.0-TI**2)
            HS(20)=+0.25*RM*(1.0-TI**2)
!  FOR T-COORDINATE
            HT(1)=-0.125*RM*SM*TM+0.125*RM*SM*(2.0+RI+SI+TI)
            HT(2)=-0.125*RP*SM*TM+0.125*RP*SM*(2.0-RI+SI+TI)
            HT(3)=-0.125*RP*SP*TM+0.125*RP*SP*(2.0-RI-SI+TI)
            HT(4)=-0.125*RM*SP*TM+0.125*RM*SP*(2.0+RI-SI+TI)
            HT(5)=+0.125*RM*SM*TP-0.125*RM*SM*(2.0+RI+SI-TI)
            HT(6)=+0.125*RP*SM*TP-0.125*RP*SM*(2.0-RI+SI-TI)
            HT(7)=+0.125*RP*SP*TP-0.125*RP*SP*(2.0-RI-SI-TI)
            HT(8)=+0.125*RM*SP*TP-0.125*RM*SP*(2.0+RI-SI-TI)
            HT(9)=-0.25*(1.0-RI**2)*SM
            HT(10)=-0.25*RP*(1.0-SI**2)
            HT(11)=-0.25*(1.0-RI**2)*SP
            HT(12)=-0.25*RM*(1.0-SI**2)
            HT(13)=0.25*(1.0-RI**2)*SM
            HT(14)=0.25*RP*(1.0-SI**2)
            HT(15)=0.25*(1.0-RI**2)*SP
            HT(16)=0.25*RM*(1.0-SI**2)
            HT(17)=-0.5*RM*SM*TI
            HT(18)=-0.5*RP*SM*TI
            HT(19)=-0.5*RP*SP*TI
            HT(20)=-0.5*RM*SP*TI
!  JACOBI MATRIX 
            XJ11=0.0
            XJ21=0.0
            XJ31=0.0
            XJ12=0.0
            XJ22=0.0
            XJ32=0.0
            XJ13=0.0
            XJ23=0.0
            XJ33=0.0
            DO I=1,NN
              XJ11=XJ11+HR(I)*XX(I)
              XJ21=XJ21+HS(I)*XX(I)
              XJ31=XJ31+HT(I)*XX(I)
              XJ12=XJ12+HR(I)*YY(I)
              XJ22=XJ22+HS(I)*YY(I)
              XJ32=XJ32+HT(I)*YY(I)
              XJ13=XJ13+HR(I)*ZZ(I)
              XJ23=XJ23+HS(I)*ZZ(I)
              XJ33=XJ33+HT(I)*ZZ(I)
            ENDDO
!DETERMINANT OF JACOBIAN
            DET=XJ11*XJ22*XJ33     &
               +XJ12*XJ23*XJ31     &
               +XJ13*XJ21*XJ32     &
               -XJ13*XJ22*XJ31     &
               -XJ12*XJ21*XJ33     &
               -XJ11*XJ23*XJ32
! INVERSION OF JACOBIAN
            DUM=1.0/DET
              XJI11=DUM*( XJ22*XJ33-XJ23*XJ32)
              XJI21=DUM*(-XJ21*XJ33+XJ23*XJ31)
              XJI31=DUM*( XJ21*XJ32-XJ22*XJ31)
              XJI12=DUM*(-XJ12*XJ33+XJ13*XJ32)
              XJI22=DUM*( XJ11*XJ33-XJ13*XJ31)
              XJI32=DUM*(-XJ11*XJ32+XJ12*XJ31)
              XJI13=DUM*( XJ12*XJ23-XJ13*XJ22)
              XJI23=DUM*(-XJ11*XJ23+XJ13*XJ21)
              XJI33=DUM*( XJ11*XJ22-XJ12*XJ21)
!* DUDX DUDY DUDZ ....
            DUDX=0.0
            DUDY=0.0
            DUDZ=0.0
            DVDX=0.0
            DVDY=0.0
            DVDZ=0.0
            DWDX=0.0
            DWDY=0.0
            DWDZ=0.0
            DO J=1,NN
              DNDX=XJI11*HR(J)+XJI12*HS(J)+XJI13*HT(J)
              DNDY=XJI21*HR(J)+XJI22*HS(J)+XJI23*HT(J)
              DNDZ=XJI31*HR(J)+XJI32*HS(J)+XJI33*HT(J)
              DUDX=DUDX+DNDX*UX(J)
              DUDY=DUDY+DNDY*UX(J)
              DUDZ=DUDZ+DNDZ*UX(J)
              DVDX=DVDX+DNDX*UY(J)
              DVDY=DVDY+DNDY*UY(J)
              DVDZ=DVDZ+DNDZ*UY(J)
              DWDX=DWDX+DNDX*UZ(J)
              DWDY=DWDY+DNDY*UZ(J)
              DWDZ=DWDZ+DNDZ*UZ(J)
            ENDDO
!** Fmat
            FMAT(1,1)=1.0+DUDX
            FMAT(1,2)=DUDY
            FMAT(1,3)=DUDZ
            FMAT(2,1)=DVDX
            FMAT(2,2)=1.0+DVDY
            FMAT(2,3)=DVDZ
            FMAT(3,1)=DWDX
            FMAT(3,2)=DWDY
            FMAT(3,3)=1.0+DWDZ
!* DETERMINANT
            DETF=FMAT(1,1)*FMAT(2,2)*FMAT(3,3)   &
                +FMAT(1,2)*FMAT(2,3)*FMAT(3,1)   &
                +FMAT(1,3)*FMAT(2,1)*FMAT(3,2)   &
                -FMAT(1,3)*FMAT(2,2)*FMAT(3,1)   &
                -FMAT(1,2)*FMAT(2,1)*FMAT(3,3)   &
                -FMAT(1,1)*FMAT(2,3)*FMAT(3,2)
!** Logarithmic strain 
            CALL GET_LE(ILOG,FMAT,STRAIN_L)
            EPS(1)=STRAIN_L(1)
            EPS(2)=STRAIN_L(2)
            EPS(3)=STRAIN_L(3)
            EPS(4)=2.0*STRAIN_L(6)
            EPS(5)=2.0*STRAIN_L(5)
            EPS(6)=2.0*STRAIN_L(4)
!** PK2 to Cauchy
            SMAT(1,1)=STRESS(1,LX,LY,LZ)
            SMAT(1,2)=STRESS(4,LX,LY,LZ)
            SMAT(1,3)=STRESS(6,LX,LY,LZ)
            SMAT(2,1)=STRESS(4,LX,LY,LZ)
            SMAT(2,2)=STRESS(2,LX,LY,LZ)
            SMAT(2,3)=STRESS(5,LX,LY,LZ)
            SMAT(3,1)=STRESS(6,LX,LY,LZ)
            SMAT(3,2)=STRESS(5,LX,LY,LZ)
            SMAT(3,3)=STRESS(3,LX,LY,LZ)
            DO J=1,3
              DO K=1,3
                WK(J,K)=0.0
                DO L=1,3
                  WK(J,K)=WK(J,K)+SMAT(J,L)*FMAT(K,L)
                ENDDO
              ENDDO
            ENDDO
            DO J=1,3
              DO K=1,3
                TMAT(J,K)=0.0
                DO L=1,3
                  TMAT(J,K)=TMAT(J,K)+FMAT(J,L)*WK(L,K)
                ENDDO
              ENDDO
            ENDDO
            DO J=1,3
              DO K=1,3
                TMAT(J,K)=TMAT(J,K)/DETF
              ENDDO
            ENDDO
!**
!** CAUCHY STRESS
!**
!ogasawara till next ogasawara
!           SGM(1)=TMAT(1,1)
!           SGM(2)=TMAT(2,2)
!           SGM(3)=TMAT(3,3)
!           SGM(4)=TMAT(2,3)
!           SGM(5)=TMAT(1,3)
!           SGM(6)=TMAT(1,2)
!
            SGM(1)=SMAT(1,1)
            SGM(2)=SMAT(2,2)
            SGM(3)=SMAT(3,3)
            SGM(4)=SMAT(2,3)
            SGM(5)=SMAT(1,3)
            SGM(6)=SMAT(1,2)
!ogasawara
!**
!** MISES
!**
            PS=(SGM(1)+SGM(2)+SGM(3))/3.0
            SMISES= DSQRT( 3.0D0*( 0.5D0                         &
            *( (SGM(1)-PS)**2 +(SGM(2)-PS)**2+(SGM(3)-PS)**2) )  &
            +3.0*SGM(4)**2 +3.0*SGM(5)**2+3.0*SGM(6)**2 )
            IC=IC+1
            ARRAY_TEMP(1,IC)=EPS(1)
            ARRAY_TEMP(2,IC)=EPS(2)
            ARRAY_TEMP(3,IC)=EPS(3)
            ARRAY_TEMP(4,IC)=EPS(4)
            ARRAY_TEMP(5,IC)=EPS(5)
            ARRAY_TEMP(6,IC)=EPS(6)
            ARRAY_TEMP(7,IC)=SGM(1)
            ARRAY_TEMP(8,IC)=SGM(2)
            ARRAY_TEMP(9,IC)=SGM(3)
            ARRAY_TEMP(10,IC)=SGM(4)
            ARRAY_TEMP(11,IC)=SGM(5)
            ARRAY_TEMP(12,IC)=SGM(6)
            ARRAY_TEMP(13,IC)=SMISES
          ENDDO
        ENDDO
      ENDDO
!**
!** NODAL VALUE
!**
      DO J=1,NN
        RI=RX(J)/XG(NG,NG)
        SI=SX(J)/XG(NG,NG)
        TI=TX(J)/XG(NG,NG)
        RP=1.0+RI
        SP=1.0+SI
        TP=1.0+TI
        RM=1.0-RI
        SM=1.0-SI
        TM=1.0-TI
!  INTERPOLATION FUNCTION (GAUSS POINT)
!  INTERPOLATION FUNCTION (GAUSS POINT)
        HG( 1)=-0.125*RM*SM*TM*(2.0+RI+SI+TI)
        HG(19)=-0.125*RP*SM*TM*(2.0-RI+SI+TI)
        HG(25)=-0.125*RP*SP*TM*(2.0-RI-SI+TI)
        HG( 7)=-0.125*RM*SP*TM*(2.0+RI-SI+TI)
        HG( 3)=-0.125*RM*SM*TP*(2.0+RI+SI-TI)
        HG(21)=-0.125*RP*SM*TP*(2.0-RI+SI-TI)
        HG(27)=-0.125*RP*SP*TP*(2.0-RI-SI-TI)
        HG( 9)=-0.125*RM*SP*TP*(2.0+RI-SI-TI)
        HG(10)=0.25*(1.0-RI**2)*SM*TM
        HG(22)=0.25*RP*(1.0-SI**2)*TM
        HG(16)=0.25*(1.0-RI**2)*SP*TM
!!        HG( 9)=0.25*RM*(1.0-SI**2)*TM
        HG( 4)=0.25*RM*(1.0-SI**2)*TM       !oga
        HG(12)=0.25*(1.0-RI**2)*SM*TP
        HG(24)=0.25*RP*(1.0-SI**2)*TP
        HG(18)=0.25*(1.0-RI**2)*SP*TP
        HG( 6)=0.25*RM*(1.0-SI**2)*TP
        HG( 2)=0.25*RM*SM*(1.0-TI**2)
        HG(20)=0.25*RP*SM*(1.0-TI**2)
        HG(26)=0.25*RP*SP*(1.0-TI**2)
        HG( 8)=0.25*RM*SP*(1.0-TI**2)
        HG( 5)=0.0
        HG(11)=0.0
        HG(13)=0.0
        HG(14)=0.0
        HG(15)=0.0
        HG(17)=0.0
        HG(23)=0.0
        DO I=1,13
          ARRAY(I,J)=0.0
          DO K=1,NG*NG*NG
            ARRAY(I,J)=ARRAY(I,J)+HG(K)*ARRAY_TEMP(I,K)
          ENDDO
        ENDDO
      ENDDO      
      RETURN
      END SUBROUTINE NRCV_C3D20TL
!----------------------------------------------------------------------*
      SUBROUTINE NRCV_C3D15TL(ILOG,XX,YY,ZZ,TT,T0,EE,PP,ALP,DISP,STRESS,ARRAY)
!----------------------------------------------------------------------*
!
! NODAL RECOVERY FOR 3D 15 NODE PLANE ELEMENT /TOTAL LAGRANGE
!
      use hecmw
      use gauss_integration 
      IMPLICIT NONE
! GAUSS POINT
      INTEGER(kind=kint) NG
      PARAMETER(NG=3)
! I/F VARIABLES
      INTEGER(kind=kint) ILOG
      REAL(kind=kreal) XX(*),YY(*),ZZ(*),TT(*),T0(*),DISP(*),ARRAY(13,*)
      REAL(kind=kreal) STRESS(6,NG,NG,NG)
      REAL(kind=kreal) EE,PP,ALP
! LOCAL VARIABLES
      INTEGER(kind=kint) NN,NDOF
      PARAMETER(NN=15,NDOF=3)
      INTEGER(kind=kint) I,J,K,L
      INTEGER(kind=kint) L1,L2,LZ
      REAL(kind=kreal) X1,X2,X3,ZI
      REAL(kind=kreal) XL1,XL2
      REAL(kind=kreal) RI,SI,TI,RP,SP,TP,RM,SM,TM
      REAL(kind=kreal) H(NN),HL1(NN),HL2(NN),HL3(NN),HZ(NN)
      REAL(kind=kreal) XJ11,XJ21,XJ31,XJ12,XJ22,XJ32,XJ13,XJ23,XJ33,DET,WG,DUM
      REAL(kind=kreal) XJI11,XJI21,XJI31,XJI12,XJI22,XJI32,XJI13,XJI23,XJI33
      REAL(kind=kreal) DNDX,DNDY,DNDZ
      REAL(kind=kreal) DUDX,DUDY,DUDZ
      REAL(kind=kreal) DVDX,DVDY,DVDZ
      REAL(kind=kreal) DWDX,DWDY,DWDZ
      REAL(kind=kreal) BL0(6,NN*NDOF),BL1(6,NN*NDOF)
      REAL(kind=kreal) UX(NN),UY(NN),UZ(NN)
      REAL(kind=kreal) S11,S22,S33,S23,S13,S12
      REAL(kind=kreal) DETF,FMAT(3,3),SMAT(3,3),WK(3,3),TMAT(3,3),STRAIN_L(6)
      REAL(kind=kreal) SGM(6),EPS(6),PS,SMISES
      REAL(kind=kreal) ARRAY_TEMP(13,NG*NG*NG)
      INTEGER(kind=kint) IC
!*************************
!  EVALUATION POINTS
!*************************
      REAL(kind=kreal) HG(NG*NG*NG),RX(NN),SX(NN),TX(NN)
      DATA RX/ 1., 0.,-1., 1., 0.,-1.,   &
               1.,-1., 0., 1.,-1., 0.,   &
               1., 0.,-1./
      DATA SX/-1., 1.,-1.,-1., 1.,-1.,   &
               0., 0.,-1., 0., 0.,-1.,   &
              -1., 1.,-1./
      DATA TX/-1.,-1.,-1., 1., 1., 1.,   &
              -1.,-1.,-1., 1.,1.,1.,     &
               0., 0., 0./
!* Total displacement 
      DO I=1,NN
        UX(I)=DISP(NDOF*I-2)
        UY(I)=DISP(NDOF*I-1)
        UZ(I)=DISP(NDOF*I  )
      ENDDO
      IC=0
! LOOP FOR INTEGRATION POINTS      
      DO LZ=1,NG
        ZI=XG(NG,LZ)
        DO L2=1,NG
          XL2=XG(NG,L2)
          X2 =(XL2+1.0)*0.5
          DO L1=1,NG
            XL1=XG(NG,L1)
            X1=0.5*(1.0-X2)*(XL1+1.0)
!  INTERPOLATION FUNCTION
            X3=1.0-X1-X2
            H(1)= 0.5*X1*(2.0*X1-2.-ZI)*(1.0-ZI)
            H(2)= 0.5*X2*(2.0*X2-2.-ZI)*(1.0-ZI)
            H(3)= 0.5*X3*(2.0*X3-2.-ZI)*(1.0-ZI)
            H(4)= 0.5*X1*(2.0*X1-2.+ZI)*(1.0+ZI)
            H(5)= 0.5*X2*(2.0*X2-2.+ZI)*(1.0+ZI)
            H(6)= 0.5*X3*(2.0*X3-2.+ZI)*(1.0+ZI)
            H(7)= 2.0*X1*X2*(1.0-ZI)
            H(8)= 2.0*X2*X3*(1.0-ZI)
            H(9)= 2.0*X1*X3*(1.0-ZI)
            H(10)=2.0*X1*X2*(1.0+ZI)
            H(11)=2.0*X2*X3*(1.0+ZI)
            H(12)=2.0*X1*X3*(1.0+ZI)
            H(13)=X1*(1.0-ZI**2)
            H(14)=X2*(1.0-ZI**2)
            H(15)=X3*(1.0-ZI**2)
!  DERIVATIVE OF INTERPOLATION FUNCTION
!  FOR L1-COORDINATE
            HL1(1)= 0.5*(4.0*X1-2.-ZI)*(1.0-ZI)
            HL1(2)= 0.0
            HL1(3)= 0.0
            HL1(4)= 0.5*(4.0*X1-2.+ZI)*(1.0+ZI)
            HL1(5)= 0.0
            HL1(6)= 0.0
            HL1(7)= 2.0*X2*(1.0-ZI)
            HL1(8)= 0.0
            HL1(9)= 2.0*X3*(1.0-ZI)
            HL1(10)=2.0*X2*(1.0+ZI)
            HL1(11)=0.0
            HL1(12)=2.0*X3*(1.0+ZI)
            HL1(13)=1.0-ZI**2
            HL1(14)=0.0
            HL1(15)=0.0
!  FOR L2-COORDINATE
            HL2(1)= 0.0
            HL2(2)= 0.5*(4.0*X2-2.-ZI)*(1.0-ZI)
            HL2(3)= 0.0
            HL2(4)= 0.0
            HL2(5)= 0.5*(4.0*X2-2.+ZI)*(1.0+ZI)
            HL2(6)= 0.0
            HL2(7)= 2.0*X1*(1.0-ZI)
            HL2(8)= 2.0*X3*(1.0-ZI)
            HL2(9)= 0.0
            HL2(10)=2.0*X1*(1.0+ZI)
            HL2(11)=2.0*X3*(1.0+ZI)
            HL2(12)=0.0
            HL2(13)=0.0
            HL2(14)=1.0-ZI**2
            HL2(15)=0.0
!  FOR L3-COORDINATE
            HL3(1)= 0.0
            HL3(2)= 0.0
            HL3(3)= 0.5*(4.0*X3-2.-ZI)*(1.0-ZI)
            HL3(4)= 0.0
            HL3(5)= 0.0
            HL3(6)= 0.5*(4.0*X3-2.+ZI)*(1.0+ZI)
            HL3(7)= 0.0
            HL3(8)= 2.0*X2*(1.0-ZI)
            HL3(9)= 2.0*X1*(1.0-ZI)
            HL3(10)=0.0
            HL3(11)=2.0*X2*(1.0+ZI)
            HL3(12)=2.0*X1*(1.0+ZI)
            HL3(13)=0.0
            HL3(14)=0.0
            HL3(15)=1.0-ZI**2
!  FOR Z-COORDINATE
            HZ(1)= 0.5*X1*(-2.0*X1+1.0+2.0*ZI)
            HZ(2)= 0.5*X2*(-2.0*X2+1.0+2.0*ZI)
            HZ(3)= 0.5*X3*(-2.0*X3+1.0+2.0*ZI)
            HZ(4)= 0.5*X1*( 2.0*X1-1.0+2.0*ZI)
            HZ(5)= 0.5*X2*( 2.0*X2-1.0+2.0*ZI)
            HZ(6)= 0.5*X3*( 2.0*X3-1.0+2.0*ZI)
            HZ(7)= -2.0*X1*X2
            HZ(8)= -2.0*X2*X3
            HZ(9)= -2.0*X1*X3
            HZ(10)= 2.0*X1*X2
            HZ(11)= 2.0*X2*X3
            HZ(12)= 2.0*X1*X3
            HZ(13)=-2.0*X1*ZI
            HZ(14)=-2.0*X2*ZI
            HZ(15)=-2.0*X3*ZI
!  JACOBI MATRIX 
            XJ11=0.0
            XJ21=0.0
            XJ31=0.0
            XJ12=0.0
            XJ22=0.0
            XJ32=0.0
            XJ13=0.0
            XJ23=0.0
            XJ33=0.0
            DO I=1,NN
              XJ11=XJ11+(HL1(I)-HL3(I))*XX(I)
              XJ21=XJ21+(HL2(I)-HL3(I))*XX(I)
              XJ31=XJ31+HZ(I)*XX(I)
              XJ12=XJ12+(HL1(I)-HL3(I))*YY(I)
              XJ22=XJ22+(HL2(I)-HL3(I))*YY(I)
              XJ32=XJ32+HZ(I)*YY(I)
              XJ13=XJ13+(HL1(I)-HL3(I))*ZZ(I)
              XJ23=XJ23+(HL2(I)-HL3(I))*ZZ(I)
              XJ33=XJ33+HZ(I)*ZZ(I)
            ENDDO
!DETERMINANT OF JACOBIAN
            DET=XJ11*XJ22*XJ33   &
              +XJ12*XJ23*XJ31    &
              +XJ13*XJ21*XJ32    &
              -XJ13*XJ22*XJ31    &
              -XJ12*XJ21*XJ33    &
              -XJ11*XJ23*XJ32
! INVERSION OF JACOBIAN
            DUM=1.0/DET
            XJI11=DUM*( XJ22*XJ33-XJ23*XJ32)
            XJI21=DUM*(-XJ21*XJ33+XJ23*XJ31)
            XJI31=DUM*( XJ21*XJ32-XJ22*XJ31)
            XJI12=DUM*(-XJ12*XJ33+XJ13*XJ32)
            XJI22=DUM*( XJ11*XJ33-XJ13*XJ31)
            XJI32=DUM*(-XJ11*XJ32+XJ12*XJ31)
            XJI13=DUM*( XJ12*XJ23-XJ13*XJ22)
            XJI23=DUM*(-XJ11*XJ23+XJ13*XJ21)
            XJI33=DUM*( XJ11*XJ22-XJ12*XJ21)
!* DUDX DUDY DUDZ ....
            DUDX=0.0
            DUDY=0.0
            DUDZ=0.0
            DVDX=0.0
            DVDY=0.0
            DVDZ=0.0
            DWDX=0.0
            DWDY=0.0
            DWDZ=0.0
            DO J=1,NN
              DNDX=XJI11*(HL1(J)-HL3(J)) +XJI12*(HL2(J)-HL3(J)) +XJI13*HZ(J)
              DNDY=XJI21*(HL1(J)-HL3(J)) +XJI22*(HL2(J)-HL3(J)) +XJI23*HZ(J)
              DNDZ=XJI31*(HL1(J)-HL3(J)) +XJI32*(HL2(J)-HL3(J)) +XJI33*HZ(J)
              DUDX=DUDX+DNDX*UX(J)
              DUDY=DUDY+DNDY*UX(J)
              DUDZ=DUDZ+DNDZ*UX(J)
              DVDX=DVDX+DNDX*UY(J)
              DVDY=DVDY+DNDY*UY(J)
              DVDZ=DVDZ+DNDZ*UY(J)
              DWDX=DWDX+DNDX*UZ(J)
              DWDY=DWDY+DNDY*UZ(J)
              DWDZ=DWDZ+DNDZ*UZ(J)
            ENDDO
!** Fmat
            FMAT(1,1)=1.0+DUDX
            FMAT(1,2)=DUDY
            FMAT(1,3)=DUDZ
            FMAT(2,1)=DVDX
            FMAT(2,2)=1.0+DVDY
            FMAT(2,3)=DVDZ
            FMAT(3,1)=DWDX
            FMAT(3,2)=DWDY
            FMAT(3,3)=1.0+DWDZ
!* DETERMINANT
            DETF=FMAT(1,1)*FMAT(2,2)*FMAT(3,3)      &
                +FMAT(1,2)*FMAT(2,3)*FMAT(3,1)      &
                +FMAT(1,3)*FMAT(2,1)*FMAT(3,2)      &
                -FMAT(1,3)*FMAT(2,2)*FMAT(3,1)      &
                -FMAT(1,2)*FMAT(2,1)*FMAT(3,3)      &
                -FMAT(1,1)*FMAT(2,3)*FMAT(3,2)
!** Logarithmic strain 
            CALL GET_LE(ILOG,FMAT,STRAIN_L)
            EPS(1)=STRAIN_L(1)
            EPS(2)=STRAIN_L(2)
            EPS(3)=STRAIN_L(3)
            EPS(4)=2.0*STRAIN_L(6)
            EPS(5)=2.0*STRAIN_L(5)
            EPS(6)=2.0*STRAIN_L(4)
!** PK2 to Cauchy
            SMAT(1,1)=STRESS(1,L1,L2,LZ)
            SMAT(1,2)=STRESS(4,L1,L2,LZ)
            SMAT(1,3)=STRESS(6,L1,L2,LZ)
            SMAT(2,1)=STRESS(4,L1,L2,LZ)
            SMAT(2,2)=STRESS(2,L1,L2,LZ)
            SMAT(2,3)=STRESS(5,L1,L2,LZ)
            SMAT(3,1)=STRESS(6,L1,L2,LZ)
            SMAT(3,2)=STRESS(5,L1,L2,LZ)
            SMAT(3,3)=STRESS(3,L1,L2,LZ)
            DO J=1,3
              DO K=1,3
                WK(J,K)=0.0
                DO L=1,3
                  WK(J,K)=WK(J,K)+SMAT(J,L)*FMAT(K,L)
                ENDDO
              ENDDO
            ENDDO
            DO J=1,3
              DO K=1,3
                TMAT(J,K)=0.0
                DO L=1,3
                  TMAT(J,K)=TMAT(J,K)+FMAT(J,L)*WK(L,K)
                ENDDO
              ENDDO
            ENDDO
            DO J=1,3
              DO K=1,3
                TMAT(J,K)=TMAT(J,K)/DETF
              ENDDO
            ENDDO
!**
!** CAUCHY STRESS
!**
!ogasawara till next ogasawara
!           SGM(1)=TMAT(1,1)
!           SGM(2)=TMAT(2,2)
!           SGM(3)=TMAT(3,3)
!           SGM(4)=TMAT(2,3)
!           SGM(5)=TMAT(1,3)
!           SGM(6)=TMAT(1,2)
!
            SGM(1)=SMAT(1,1)
            SGM(2)=SMAT(2,2)
            SGM(3)=SMAT(3,3)
            SGM(4)=SMAT(2,3)
            SGM(5)=SMAT(1,3)
            SGM(6)=SMAT(1,2)
!ogasawara
!**
!** MISES
!**
            PS=(SGM(1)+SGM(2)+SGM(3))/3.0
            SMISES= DSQRT( 3.0D0*( 0.5D0                          &
            *( (SGM(1)-PS)**2 +(SGM(2)-PS)**2+(SGM(3)-PS)**2) )   &
            +3.0*SGM(4)**2 +3.0*SGM(5)**2+3.0*SGM(6)**2 )
            IC=IC+1
            ARRAY_TEMP(1,IC)=EPS(1)
            ARRAY_TEMP(2,IC)=EPS(2)
            ARRAY_TEMP(3,IC)=EPS(3)
            ARRAY_TEMP(4,IC)=EPS(4)
            ARRAY_TEMP(5,IC)=EPS(5)
            ARRAY_TEMP(6,IC)=EPS(6)
            ARRAY_TEMP(7,IC)=SGM(1)
            ARRAY_TEMP(8,IC)=SGM(2)
            ARRAY_TEMP(9,IC)=SGM(3)
            ARRAY_TEMP(10,IC)=SGM(4)
            ARRAY_TEMP(11,IC)=SGM(5)
            ARRAY_TEMP(12,IC)=SGM(6)
            ARRAY_TEMP(13,IC)=SMISES
          ENDDO
        ENDDO
      ENDDO
!**
!** NODAL VALUE
!**
      DO J=1,NN
        RI=RX(J)/XG(NG,NG)
        SI=SX(J)/XG(NG,NG)
        TI=TX(J)/XG(NG,NG)
        RP=1.0+RI
        SP=1.0+SI
        TP=1.0+TI
        RM=1.0-RI
        SM=1.0-SI
        TM=1.0-TI
!  INTERPOLATION FUNCTION (GAUSS POINT)
!oga
!!        HG(1)=0.125*RM*SM*TM
!!        HG(5)=0.125*RP*SM*TM
!!        HG(7)=0.125*RP*SP*TM
!!        HG(3)=0.125*RM*SP*TM
!!        HG(2)=0.125*RM*SM*TP
!!        HG(6)=0.125*RP*SM*TP
!!        HG(8)=0.125*RP*SP*TP
!!        HG(4)=0.125*RM*SP*TP
        HG( 1)=-0.125*RM*SM*TM*(2.0+RI+SI+TI)
        HG(19)=-0.125*RP*SM*TM*(2.0-RI+SI+TI)
        HG(25)=-0.125*RP*SP*TM*(2.0-RI-SI+TI)
        HG( 7)=-0.125*RM*SP*TM*(2.0+RI-SI+TI)
        HG( 3)=-0.125*RM*SM*TP*(2.0+RI+SI-TI)
        HG(21)=-0.125*RP*SM*TP*(2.0-RI+SI-TI)
        HG(27)=-0.125*RP*SP*TP*(2.0-RI-SI-TI)
        HG( 9)=-0.125*RM*SP*TP*(2.0+RI-SI-TI)
        HG(10)=0.25*(1.0-RI**2)*SM*TM
        HG(22)=0.25*RP*(1.0-SI**2)*TM
        HG(16)=0.25*(1.0-RI**2)*SP*TM
!!        HG( 9)=0.25*RM*(1.0-SI**2)*TM
        HG( 4)=0.25*RM*(1.0-SI**2)*TM
        HG(12)=0.25*(1.0-RI**2)*SM*TP
        HG(24)=0.25*RP*(1.0-SI**2)*TP
        HG(18)=0.25*(1.0-RI**2)*SP*TP
        HG( 6)=0.25*RM*(1.0-SI**2)*TP
        HG( 2)=0.25*RM*SM*(1.0-TI**2)
        HG(20)=0.25*RP*SM*(1.0-TI**2)
        HG(26)=0.25*RP*SP*(1.0-TI**2)
        HG( 8)=0.25*RM*SP*(1.0-TI**2)
        HG( 5)=0.0
        HG(11)=0.0
        HG(13)=0.0
        HG(14)=0.0
        HG(15)=0.0
        HG(17)=0.0
        HG(23)=0.0
!oga
        DO I=1,13
          ARRAY(I,J)=0.0
          DO K=1,NG*NG*NG
            ARRAY(I,J)=ARRAY(I,J)+HG(K)*ARRAY_TEMP(I,K)
          ENDDO
        ENDDO
      ENDDO      
      RETURN
      END SUBROUTINE NRCV_C3D15TL
!----------------------------------------------------------------------*
      SUBROUTINE NRCV_C3D10TL(ILOG,XX,YY,ZZ,TT,T0,EE,PP,ALP,DISP,STRESS,ARRAY,  &
                              SURFF,STRAIN,elm_ic,i_step,idx_stres,ARRAY_TEMP,  &
                              ARRAY_PRIS,IDX_OUT_COOD_ST,MAT_ISO,BIJ,G_ELM_ID)
!----------------------------------------------------------------------*
!
! NODAL RECOVERY FOR 3D 10 NODE PLANE ELEMENT /TOTAL LAGRANGE
!
!     use m_fstr
      use hecmw
      use gauss_integration 
      IMPLICIT NONE

! GAUSS POINT
      INTEGER(kind=kint) NG
      PARAMETER(NG=3)
! I/F VARIABLES
      INTEGER(kind=kint) ILOG
      INTEGER(kind=kint) idx_stres
      REAL(kind=kreal) XX(*),YY(*),ZZ(*),TT(*),T0(*),DISP(*),ARRAY(13,*),ARRAY_PRIS(24,*)
      REAL(kind=kreal) ARRAY_2(6),ARRAY_PRIS_2(12) 
      REAL(kind=kreal) STRESS(6,NG,NG,NG),STRAIN(6,NG,NG,NG)
      REAL(kind=kreal) SURFF(8,NG,NG,NG)
      REAL(kind=kreal) EE,PP,ALP
! LOCAL VARIABLES
      INTEGER(kind=kint) NN,NDOF
      PARAMETER(NN=10,NDOF=3)
      INTEGER(kind=kint) I,J,K,L,M 
      INTEGER(kind=kint) LX,LY,LZ
      INTEGER(kind=kint) L1,L2,L3
      REAL(kind=kreal) X1,X2,X3,X4
      REAL(kind=kreal) XL1,XL2,XL3
      REAL(kind=kreal) RI,SI,TI,RP,SP,TP,RM,SM,TM
      REAL(kind=kreal) H(NN),HL1(NN),HL2(NN),HL3(NN),HL4(NN)
      REAL(kind=kreal) XJ11,XJ21,XJ31,XJ12,XJ22,XJ32,XJ13,XJ23,XJ33,DET,WG,DUM
      REAL(kind=kreal) XJI11,XJI21,XJI31,XJI12,XJI22,XJI32,XJI13,XJI23,XJI33
      REAL(kind=kreal) DNDX,DNDY,DNDZ
      REAL(kind=kreal) DUDX,DUDY,DUDZ
      REAL(kind=kreal) DVDX,DVDY,DVDZ
      REAL(kind=kreal) DWDX,DWDY,DWDZ
      REAL(kind=kreal) BL0(6,NN*NDOF),BL1(6,NN*NDOF)
      REAL(kind=kreal) UX(NN),UY(NN),UZ(NN)
      REAL(kind=kreal) S11,S22,S33,S23,S13,S12
      REAL(kind=kreal) DETF,FMAT(3,3),FMATINV(3,3),SMAT(3,3),WK(3,3),TMAT(3,3),STRAIN_L(6)
      REAL(kind=kreal) PIMAT(3,3),NORML0(3),NORML1(3) 
      REAL(kind=kreal) SGM(6),EPS(6),PS,SMISES
      REAL(kind=kreal) ARRAY_TEMP(13,NG*NG*NG)
      INTEGER(kind=kint) elm_ic,G_ELM_ID,IC
      INTEGER(kind=kint) i_step
!*************************
!  EVALUATION POINTS
!*************************
      REAL(kind=kreal) HG(NG*NG*NG),RX(NN),SX(NN),TX(NN)
      DATA RX/ 1., 0., 0.,-1.,          &
               1., 0., 1.,              &
               0.,-1.,-1./
      DATA SX/-1., 1., 0.,-1.,          &
               0., 1.,-1.,              &
              -1., 0.,-1./
      DATA TX/-1.,-1., 1.,-1.,          &
              -1., 0., 0.,              &
              -1.,-1., 0./
!
      INTEGER(kind=kint) IDX_OUT_COOD_ST,MAT_ISO 
      REAL(kind=kreal) BIJ(3,3) 

      REAL(kind=kreal) EPS_G(3,3),EPS_GB(3,3),SIG_G(3,3),SIG_GB(3,3) 

!* Total displacement 
      DO I=1,NN
        UX(I)=DISP(NDOF*I-2)
        UY(I)=DISP(NDOF*I-1)
        UZ(I)=DISP(NDOF*I  )
      ENDDO
!
      IC=0
! LOOP FOR INTEGRATION POINTS      
       DO L3=1,NG
        XL3=XG(NG,L3)
        X3 =(XL3+1.0)*0.5
        DO L2=1,NG
          XL2=XG(NG,L2)
          X2 =(1.0-X3)*(XL2+1.0)*0.5
          DO L1=1,NG
            XL1=XG(NG,L1)
            X1=(1.0-X2-X3)*(XL1+1.0)*0.5

!  INTERPOLATION FUNCTION
            X4=1.0-X1-X2-X3
            H(1)=X1*(2.0*X1-1.0)
            H(2)=X2*(2.0*X2-1.0)
            H(3)=X3*(2.0*X3-1.0)
            H(4)=X4*(2.0*X4-1.0)
            H(5)=4.0*X1*X2
            H(6)=4.0*X2*X3
            H(7)=4.0*X1*X3
            H(8)=4.0*X1*X4
            H(9)=4.0*X2*X4
            H(10)=4.0*X3*X4
! DERIVATIVE OF INTERPOLATION FUNCTION
! FOR L1-COORDINATE
            HL1(1)= 4.0*X1-1.0
            HL1(2)= 0.0
            HL1(3)= 0.0
            HL1(4)= 0.0
            HL1(5)= 4.0*X2
            HL1(6)= 0.0
            HL1(7)= 4.0*X3
            HL1(8)= 4.0*X4
            HL1(9)= 0.0
            HL1(10)= 0.0
!  FOR L2-COORDINATE
            HL2(1)= 0.0 
            HL2(2)= 4.0*X2-1.0
            HL2(3)= 0.0
            HL2(4)= 0.0
            HL2(5)= 4.0*X1
            HL2(6)= 4.0*X3
            HL2(7)= 0.0
            HL2(8)= 0.0
            HL2(9)= 4.0*X4
            HL2(10)= 0.0
!  FOR L3-COORDINATE
            HL3(1)= 0.0
            HL3(2)= 0.0
            HL3(3)= 4.0*X3-1.0
            HL3(4)= 0.0
            HL3(5)= 0.0
            HL3(6)= 4.0*X2
            HL3(7)= 4.0*X1
            HL3(8)= 0.0
            HL3(9)= 0.0
            HL3(10)= 4.0*X4
!  FOR L4-COORDINATE
            HL4(1)= 0.0
            HL4(2)= 0.0
            HL4(3)= 0.0 
            HL4(4)= 4.0*X4-1.0
            HL4(5)= 0.0
            HL4(6)= 0.0
            HL4(7)= 0.0
            HL4(8)= 4.0*X1
            HL4(9)= 4.0*X2
            HL4(10)= 4.0*X3
!  JACOBI MATRIX 
            XJ11=0.0
            XJ21=0.0
            XJ31=0.0
            XJ12=0.0
            XJ22=0.0
            XJ32=0.0
            XJ13=0.0
            XJ23=0.0
            XJ33=0.0
            DO I=1,NN
              XJ11=XJ11+(HL1(I)-HL4(I))*XX(I)
              XJ21=XJ21+(HL2(I)-HL4(I))*XX(I)
              XJ31=XJ31+(HL3(I)-HL4(I))*XX(I)
              XJ12=XJ12+(HL1(I)-HL4(I))*YY(I)
              XJ22=XJ22+(HL2(I)-HL4(I))*YY(I)
              XJ32=XJ32+(HL3(I)-HL4(I))*YY(I)
              XJ13=XJ13+(HL1(I)-HL4(I))*ZZ(I)
              XJ23=XJ23+(HL2(I)-HL4(I))*ZZ(I)
              XJ33=XJ33+(HL3(I)-HL4(I))*ZZ(I)
            ENDDO
!DETERMINANT OF JACOBIAN
            DET=XJ11*XJ22*XJ33      &
                 +XJ12*XJ23*XJ31    &
                 +XJ13*XJ21*XJ32    &
                 -XJ13*XJ22*XJ31    &
                 -XJ12*XJ21*XJ33    &
                 -XJ11*XJ23*XJ32
! INVERSION OF JACOBIAN
            DUM=1.0/DET
            XJI11=DUM*( XJ22*XJ33-XJ23*XJ32)
            XJI21=DUM*(-XJ21*XJ33+XJ23*XJ31)
            XJI31=DUM*( XJ21*XJ32-XJ22*XJ31)
            XJI12=DUM*(-XJ12*XJ33+XJ13*XJ32)
            XJI22=DUM*( XJ11*XJ33-XJ13*XJ31)
            XJI32=DUM*(-XJ11*XJ32+XJ12*XJ31)
            XJI13=DUM*( XJ12*XJ23-XJ13*XJ22)
            XJI23=DUM*(-XJ11*XJ23+XJ13*XJ21)
            XJI33=DUM*( XJ11*XJ22-XJ12*XJ21)
!* DUDX DUDY DUDZ ....
            DUDX=0.0
            DUDY=0.0
            DUDZ=0.0
            DVDX=0.0
            DVDY=0.0
            DVDZ=0.0
            DWDX=0.0
            DWDY=0.0
            DWDZ=0.0
            DO J=1,NN
              DNDX=XJI11*(HL1(J)-HL4(J)) +XJI12*(HL2(J)-HL4(J)) +XJI13*(HL3(J)-HL4(J))
              DNDY=XJI21*(HL1(J)-HL4(J)) +XJI22*(HL2(J)-HL4(J)) +XJI23*(HL3(J)-HL4(J))
              DNDZ=XJI31*(HL1(J)-HL4(J)) +XJI32*(HL2(J)-HL4(J)) +XJI33*(HL3(J)-HL4(J))
              DUDX=DUDX+DNDX*UX(J)
              DUDY=DUDY+DNDY*UX(J)
              DUDZ=DUDZ+DNDZ*UX(J)
              DVDX=DVDX+DNDX*UY(J)
              DVDY=DVDY+DNDY*UY(J)
              DVDZ=DVDZ+DNDZ*UY(J)
              DWDX=DWDX+DNDX*UZ(J)
              DWDY=DWDY+DNDY*UZ(J)
              DWDZ=DWDZ+DNDZ*UZ(J)
            ENDDO
!** Fmat
            FMAT(1,1)=1.0+DUDX
            FMAT(1,2)=DUDY
            FMAT(1,3)=DUDZ
            FMAT(2,1)=DVDX
            FMAT(2,2)=1.0+DVDY
            FMAT(2,3)=DVDZ
            FMAT(3,1)=DWDX
            FMAT(3,2)=DWDY
            FMAT(3,3)=1.0+DWDZ
!* DETERMINANT
            DETF=FMAT(1,1)*FMAT(2,2)*FMAT(3,3)   &
                +FMAT(1,2)*FMAT(2,3)*FMAT(3,1)   &
                +FMAT(1,3)*FMAT(2,1)*FMAT(3,2)   &
                -FMAT(1,3)*FMAT(2,2)*FMAT(3,1)   &
                -FMAT(1,2)*FMAT(2,1)*FMAT(3,3)   &
                -FMAT(1,1)*FMAT(2,3)*FMAT(3,2)
! INVERSION OF FMAT
            IF(DETF.EQ.0.0) STOP ' STOP DETF '
            DUM=1.0/DETF
            FMATINV(1,1)=DUM*( FMAT(2,2)*FMAT(3,3)-FMAT(2,3)*FMAT(3,2))
            FMATINV(2,1)=DUM*(-FMAT(2,1)*FMAT(3,3)+FMAT(2,3)*FMAT(3,1))
            FMATINV(3,1)=DUM*( FMAT(2,1)*FMAT(3,2)-FMAT(2,2)*FMAT(3,1))
            FMATINV(1,2)=DUM*(-FMAT(1,2)*FMAT(3,3)+FMAT(1,3)*FMAT(3,2))
            FMATINV(2,2)=DUM*( FMAT(1,1)*FMAT(3,3)-FMAT(1,3)*FMAT(3,1))
            FMATINV(3,2)=DUM*(-FMAT(1,1)*FMAT(3,2)+FMAT(1,2)*FMAT(3,1))
            FMATINV(1,3)=DUM*( FMAT(1,2)*FMAT(2,3)-FMAT(1,3)*FMAT(2,2))
            FMATINV(2,3)=DUM*(-FMAT(1,1)*FMAT(2,3)+FMAT(1,3)*FMAT(2,1))
            FMATINV(3,3)=DUM*( FMAT(1,1)*FMAT(2,2)-FMAT(1,2)*FMAT(2,1))
!
            NORML0(1) = 0.0
            NORML0(2) = 0.0
            NORML0(3) = 1.0

            DO I=1,3
             NORML1(I) = 0.0
             DO J=1,3
              NORML1(I) = NORML1(I) + FMATINV(J,I)*NORML0(J) 
             END DO
            END DO

            DUM = (NORML1(1)**2 +NORML1(2)**2 +NORML1(3)**2 )**0.5 
            IF(DUM.EQ.0.0) STOP ' STOP NORML1 '
            DO I=1,3
             NORML1(I) = NORML1(I)/DUM 
            END DO

!** Logarithmic strain 
        IF(i_step .EQ. -9999) THEN
            CALL GET_LE(ILOG,FMAT,STRAIN_L)
        END IF

            EPS(1)   =STRAIN(1,L1,L2,L3)
            EPS(2)   =STRAIN(2,L1,L2,L3)
            EPS(3)   =STRAIN(3,L1,L2,L3)
            EPS(4)   =STRAIN(4,L1,L2,L3)
            EPS(5)   =STRAIN(5,L1,L2,L3)
            EPS(6)   =STRAIN(6,L1,L2,L3)
!** PK2 to Cauchy
!! 2nd P-K stress
            SMAT(1,1)=STRESS(1,L1,L2,L3)
            SMAT(1,2)=STRESS(4,L1,L2,L3)
            SMAT(1,3)=STRESS(6,L1,L2,L3)
            SMAT(2,1)=STRESS(4,L1,L2,L3)
            SMAT(2,2)=STRESS(2,L1,L2,L3)
            SMAT(2,3)=STRESS(5,L1,L2,L3)
            SMAT(3,1)=STRESS(6,L1,L2,L3)
            SMAT(3,2)=STRESS(5,L1,L2,L3)
            SMAT(3,3)=STRESS(3,L1,L2,L3)
!!      IF(i_step .EQ. -9999) THEN
            DO J=1,3
              DO K=1,3
                WK(J,K)=0.0
                DO L=1,3
                  WK(J,K)=WK(J,K)+SMAT(J,L)*FMAT(K,L)
                ENDDO
              ENDDO
            ENDDO
!! 1st P-K stress
            DO J=1,3
              DO K=1,3
                PIMAT(J,K)=WK(J,K)
              ENDDO
            ENDDO
!! Cauchy stress
            DO J=1,3
              DO K=1,3
                TMAT(J,K)=0.0
                DO L=1,3
                  TMAT(J,K)=TMAT(J,K)+FMAT(J,L)*WK(L,K)
                ENDDO
              ENDDO
            ENDDO
            DO J=1,3
              DO K=1,3
                TMAT(J,K)=TMAT(J,K)/DETF
              ENDDO
            ENDDO
!!
            SURFF(1,L1,L2,L3) = 0.0
            DO J=1,3
              DO K=1,3
                SURFF(1,L1,L2,L3)=SURFF(1,L1,L2,L3)+TMAT(J,K)*NORML1(J)*NORML1(K)
              ENDDO
            ENDDO
!
            DO J=1,3
              SURFF(1+J,L1,L2,L3) = 0.0
              DO K=1,3
                SURFF(1+J,L1,L2,L3)=SURFF(1+J,L1,L2,L3)+TMAT(K,J)*NORML1(K)
              ENDDO
            ENDDO
            SURFF(5,L1,L2,L3) = (SURFF(2,L1,L2,L3)**2 +SURFF(3,L1,L2,L3)**2 +SURFF(4,L1,L2,L3)**2 )**0.5 
            DO J=1,3
              SURFF(5+J,L1,L2,L3) = NORML1(J) 
            ENDDO

!!      END IF
!**
!** CAUCHY STRESS
!**
!
        IF(idx_stres .EQ. 0) THEN
           SGM(1)=PIMAT(1,1)
           SGM(2)=PIMAT(2,2)
           SGM(3)=PIMAT(3,3)
           SGM(4)=PIMAT(1,2)
           SGM(5)=PIMAT(2,3)
           SGM(6)=PIMAT(3,1)

        ELSE IF(idx_stres .EQ. 1) THEN
           SGM(1)=SMAT(1,1)
           SGM(2)=SMAT(2,2)
           SGM(3)=SMAT(3,3)
           SGM(4)=SMAT(1,2)
           SGM(5)=SMAT(2,3)
           SGM(6)=SMAT(3,1)

        ELSE IF(idx_stres .EQ. 2) THEN
           SGM(1)=TMAT(1,1)
           SGM(2)=TMAT(2,2)
           SGM(3)=TMAT(3,3)
           SGM(4)=TMAT(1,2)
           SGM(5)=TMAT(2,3)
           SGM(6)=TMAT(3,1)

        END IF
!**
!** MISES
!**
            PS=(SGM(1)+SGM(2)+SGM(3))/3.0
            SMISES= DSQRT( 3.0D0*( 0.5D0                         &
            *( (SGM(1)-PS)**2 +(SGM(2)-PS)**2+(SGM(3)-PS)**2) )  &
            +3.0*SGM(4)**2 +3.0*SGM(5)**2+3.0*SGM(6)**2 )
            IC=IC+1
            ARRAY_TEMP(1,IC)=EPS(1)
            ARRAY_TEMP(2,IC)=EPS(2)
            ARRAY_TEMP(3,IC)=EPS(3)
            ARRAY_TEMP(4,IC)=EPS(4)
            ARRAY_TEMP(5,IC)=EPS(5)
            ARRAY_TEMP(6,IC)=EPS(6)
            ARRAY_TEMP(7,IC)=SGM(1)
            ARRAY_TEMP(8,IC)=SGM(2)
            ARRAY_TEMP(9,IC)=SGM(3)
            ARRAY_TEMP(10,IC)=SGM(4)
            ARRAY_TEMP(11,IC)=SGM(5)
            ARRAY_TEMP(12,IC)=SGM(6)
            ARRAY_TEMP(13,IC)=SMISES
          ENDDO
        ENDDO
      ENDDO
!**
!** NODAL VALUE
!**
    IF(i_step .EQ. -9999) THEN
      DO J=1,NN
        RI=RX(J)/XG(NG,NG)
        SI=SX(J)/XG(NG,NG)
        TI=TX(J)/XG(NG,NG)
        RP=1.0+RI
        SP=1.0+SI
        TP=1.0+TI
        RM=1.0-RI
        SM=1.0-SI
        TM=1.0-TI
!  INTERPOLATION FUNCTION (GAUSS POINT)
        HG( 1)=-0.125*RM*SM*TM*(2.0+RI+SI+TI)
        HG( 3)=-0.125*RP*SM*TM*(2.0-RI+SI+TI)
        HG( 9)=-0.125*RP*SP*TM*(2.0-RI-SI+TI)
        HG( 7)=-0.125*RM*SP*TM*(2.0+RI-SI+TI)
        HG(19)=-0.125*RM*SM*TP*(2.0+RI+SI-TI)
        HG(21)=-0.125*RP*SM*TP*(2.0-RI+SI-TI)
        HG(27)=-0.125*RP*SP*TP*(2.0-RI-SI-TI)
        HG(25)=-0.125*RM*SP*TP*(2.0+RI-SI-TI)
        HG( 2)=0.25*(1.0-RI**2)*SM*TM
        HG( 6)=0.25*RP*(1.0-SI**2)*TM
        HG( 8)=0.25*(1.0-RI**2)*SP*TM
        HG( 4)=0.25*RM*(1.0-SI**2)*TM
        HG(20)=0.25*(1.0-RI**2)*SM*TP
        HG(24)=0.25*RP*(1.0-SI**2)*TP
        HG(26)=0.25*(1.0-RI**2)*SP*TP
        HG(22)=0.25*RM*(1.0-SI**2)*TP
        HG(10)=0.25*RM*SM*(1.0-TI**2)
        HG(12)=0.25*RP*SM*(1.0-TI**2)
        HG(18)=0.25*RP*SP*(1.0-TI**2)
        HG(16)=0.25*RM*SP*(1.0-TI**2)
        HG( 5)=0.0
        HG(11)=0.0
        HG(13)=0.0
        HG(14)=0.0
        HG(15)=0.0
        HG(17)=0.0
        HG(23)=0.0
!!
        DO I=1,13
          ARRAY(I,J)=0.0
          DO K=1,NG*NG*NG
            ARRAY(I,J)=ARRAY(I,J)+HG(K)*ARRAY_TEMP(I,K)
          ENDDO
        ENDDO
      ENDDO      
    END IF
!!
        DO I=1,13
          ARRAY(I, 1)=ARRAY_TEMP(I, 3)
          ARRAY(I, 2)=ARRAY_TEMP(I, 8)
          ARRAY(I, 3)=ARRAY_TEMP(I,23)
          ARRAY(I, 4)=ARRAY_TEMP(I, 1)
          ARRAY(I, 5)=ARRAY_TEMP(I, 6)
          ARRAY(I, 6)=ARRAY_TEMP(I,18)
          ARRAY(I, 7)=ARRAY_TEMP(I,12)
          ARRAY(I, 8)=ARRAY_TEMP(I, 2)
          ARRAY(I, 9)=ARRAY_TEMP(I, 4)
          ARRAY(I,10)=ARRAY_TEMP(I,10)
        ENDDO
!!
      DO M=1,NN
       DO I=1,6
        ARRAY_2(I) = ARRAY(I,M) 
       END DO
!!     call fstr_elpl_eig_jcb(ARRAY_2, ARRAY_PRIS_2,i_step,G_ELM_ID) 
       call fstr_elpl_eig_qr (ARRAY_2, ARRAY_PRIS_2,i_step,G_ELM_ID) 
       DO I=1,12
        ARRAY_PRIS(I,M) = ARRAY_PRIS_2(I) 
       END DO
!
       DO I=1,6
        ARRAY_2(I) = ARRAY(6+I,M) 
       END DO
!!     call fstr_elpl_eig_jcb(ARRAY_2, ARRAY_PRIS_2,i_step,G_ELM_ID) 
       call fstr_elpl_eig_qr (ARRAY_2, ARRAY_PRIS_2,i_step,G_ELM_ID) 
       DO I=1,12
        ARRAY_PRIS(12+I,M) = ARRAY_PRIS_2(I) 
       END DO
      END DO
!!
      IF(IDX_OUT_COOD_ST .EQ. 1) THEN
        DO M=1,NN
         DO I=1,3
          DO J=1,3
           EPS_G (I,J) = 0.0
           EPS_GB(I,J) = 0.0
           SIG_G (I,J) = 0.0
           SIG_GB(I,J) = 0.0
          END DO
         END DO

         EPS_G (1,1) = ARRAY( 1, M)
         EPS_G (1,2) = ARRAY( 4, M)
         EPS_G (1,3) = ARRAY( 6, M)
         EPS_G (2,1) = ARRAY( 4, M)
         EPS_G (2,2) = ARRAY( 2, M)
         EPS_G (2,3) = ARRAY( 5, M)
         EPS_G (3,1) = ARRAY( 6, M)
         EPS_G (3,2) = ARRAY( 5, M)
         EPS_G (3,3) = ARRAY( 3, M)
         SIG_G (1,1) = ARRAY( 7, M)
         SIG_G (1,2) = ARRAY(10, M)
         SIG_G (1,3) = ARRAY(12, M)
         SIG_G (2,1) = ARRAY(10, M)
         SIG_G (2,2) = ARRAY( 8, M)
         SIG_G (2,3) = ARRAY(11, M)
         SIG_G (3,1) = ARRAY(12, M)
         SIG_G (3,2) = ARRAY(11, M)
         SIG_G (3,3) = ARRAY( 9, M)

       if(MAT_ISO .eq. 0) then
         DO I=1,3
          DO J=1,3
            EPS_GB(I,J) = EPS_G(I,J) 
            SIG_GB(I,J) = SIG_G(I,J) 
          END DO
         END DO
       else if(MAT_ISO .eq. 1) then
         DO I=1,3
          DO J=1,3
           DO K=1,3
            DO L=1,3
             EPS_GB(I,J) = EPS_GB(I,J) + BIJ(K,I)*BIJ(L,J)*EPS_G(K,L)
             SIG_GB(I,J) = SIG_GB(I,J) + BIJ(K,I)*BIJ(L,J)*SIG_G(K,L)
            END DO
           END DO
          END DO
         END DO
       end if

         PS=(SIG_GB(1,1)+SIG_GB(2,2)+SIG_GB(3,3))/3.0
         SMISES= DSQRT( 3.0D0*( 0.5D0                                           &
            *( (SIG_GB(1,1)-PS)**2 +(SIG_GB(2,2)-PS)**2+(SIG_GB(3,3)-PS)**2) )  &
            +3.0*SIG_GB(1,2)**2 +3.0*SIG_GB(2,3)**2+3.0*SIG_GB(3,1)**2 )
         ARRAY( 1, M) = EPS_GB(1,1) 
         ARRAY( 2, M) = EPS_GB(2,2) 
         ARRAY( 3, M) = EPS_GB(3,3) 
         ARRAY( 4, M) = EPS_GB(1,2) 
         ARRAY( 5, M) = EPS_GB(2,3) 
         ARRAY( 6, M) = EPS_GB(3,1) 
         ARRAY( 7, M) = SIG_GB(1,1) 
         ARRAY( 8, M) = SIG_GB(2,2) 
         ARRAY( 9, M) = SIG_GB(3,3) 
         ARRAY(10, M) = SIG_GB(1,2) 
         ARRAY(11, M) = SIG_GB(2,3) 
         ARRAY(12, M) = SIG_GB(3,1) 
         ARRAY(13, M) = SMISES 
        ENDDO
      END IF
!
      RETURN
      END SUBROUTINE NRCV_C3D10TL
!----------------------------------------------------------------------*
       SUBROUTINE TLOAD_C3D10TL(ILOG,XX,YY,ZZ,DISP,DU,EE,PP,STRESS                 &
                           ,idx_elpl,h_dash,FLG_ELPL,STR_YLD ,STRAIN,EP_BAR ,iter  &
                           ,TT,TTI,CURE3D,CURE3D_I,ccc_rl,ramd_chm,VECT)
!----------------------------------------------------------------------*
      use hecmw
      use gauss_integration 
      IMPLICIT NONE
! GAUSS POINT
      INTEGER(kind=kint) NG
      PARAMETER(NG=3)
! I/F VARIABLES
      INTEGER(kind=kint) ILOG
      REAL(kind=kreal) XX(*),YY(*),ZZ(*),TT(*),TTI(*),DISP(*),DU(*),EE,PP
      REAL(kind=kreal) VECT(*)
      REAL(kind=kreal) STRESS(6,NG,NG,NG)
      INTEGER(kind=kint) iter
      INTEGER(kind=kint) idx_elpl,FLG_ELPL(NG,NG,NG)
      REAL(kind=kreal) h_dash,STR_YLD(NG,NG,NG),EP_BAR(NG,NG,NG),STRAIN(6,NG,NG,NG)
! LOCAL VARIABLES
      INTEGER(kind=kint) NN,NDOF
      PARAMETER(NN=10,NDOF=3)
      INTEGER(kind=kint) I,J,K,NUM
      INTEGER(kind=kint) L1,L2,L3
      INTEGER(kind=kint) HEAD
      REAL(kind=kreal) RI,SI,TI,RP,SP,TP,RM,SM,TM
      REAL(kind=kreal) H(NN),HL1(NN),HL2(NN),HL3(NN),HL4(NN)
      REAL(kind=kreal) X1,X2,X3,X4,XL1,XL2,XL3
      REAL(kind=kreal) XJ11,XJ21,XJ31,XJ12,XJ22,XJ32,XJ13,XJ23,XJ33,DET,WG,DUM
      REAL(kind=kreal) XJI11,XJI21,XJI31,XJI12,XJI22,XJI32,XJI13,XJI23,XJI33
      REAL(kind=kreal) DMAT(6,6),SMAT(9,9)
      REAL(kind=kreal) DNDX,DNDY,DNDZ
      REAL(kind=kreal) DUDX,DUDY,DUDZ
      REAL(kind=kreal) DVDX,DVDY,DVDZ
      REAL(kind=kreal) DWDX,DWDY,DWDZ
      REAL(kind=kreal) BL0(6,NN*NDOF),BL1(6,NN*NDOF)
      REAL(kind=kreal) UX(NN),UY(NN),UZ(NN)
      REAL(kind=kreal) S11,S22,S33,S23,S13,S12
      REAL(kind=kreal) DE(6),DS(6)
      REAL(kind=kreal) TEMPC,TEMPI,CURE3D(NG,NG,NG),CURE3D_I(NG,NG,NG),ccc_rl(*)
      REAL(kind=kreal) ramd_chm,ALP,THERMAL_EPS,EPSTH(6)

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

! CLEAR VECT
      DO I=1,NN*NDOF
        VECT(I)=0.0
      ENDDO
!* Total displacement
      DO I=1,NN
        UX(I)=DISP(NDOF*I-2)
        UY(I)=DISP(NDOF*I-1)
        UZ(I)=DISP(NDOF*I  )
      ENDDO
! LOOP FOR INTEGRATION POINTS
      DO L3=1,NG
        XL3=XG(NG,L3)
        X3 =(XL3+1.0)*0.5
        DO L2=1,NG
          XL2=XG(NG,L2)
          X2 =(1.0-X3)*(XL2+1.0)*0.5
          DO L1=1,NG
            XL1=XG(NG,L1)
            X1=(1.0-X2-X3)*(XL1+1.0)*0.5
!  INTERPOLATION FUNCTION
            X4=1.0-X1-X2-X3
            H(1)=X1*(2.0*X1-1.0)
            H(2)=X2*(2.0*X2-1.0)
            H(3)=X3*(2.0*X3-1.0)
            H(4)=X4*(2.0*X4-1.0)
            H(5)=4.0*X1*X2
            H(6)=4.0*X2*X3
            H(7)=4.0*X1*X3
            H(8)=4.0*X1*X4
            H(9)=4.0*X2*X4
            H(10)=4.0*X3*X4
! DERIVATIVE OF INTERPOLATION FUNCTION
! FOR L1-COORDINATE
            HL1(1)= 4.0*X1-1.0
            HL1(2)= 0.0
            HL1(3)= 0.0
            HL1(4)= 0.0
            HL1(5)= 4.0*X2
            HL1(6)= 0.0
            HL1(7)= 4.0*X3
            HL1(8)= 4.0*X4
            HL1(9)= 0.0
            HL1(10)= 0.0
!  FOR L2-COORDINATE
            HL2(1)= 0.0 
            HL2(2)= 4.0*X2-1.0
            HL2(3)= 0.0
            HL2(4)= 0.0
            HL2(5)= 4.0*X1
            HL2(6)= 4.0*X3
            HL2(7)= 0.0
            HL2(8)= 0.0
            HL2(9)= 4.0*X4
            HL2(10)= 0.0
!  FOR L3-COORDINATE
            HL3(1)= 0.0
            HL3(2)= 0.0
            HL3(3)= 4.0*X3-1.0
            HL3(4)= 0.0
            HL3(5)= 0.0
            HL3(6)= 4.0*X2
            HL3(7)= 4.0*X1
            HL3(8)= 0.0
            HL3(9)= 0.0
            HL3(10)= 4.0*X4
!  FOR L4-COORDINATE
            HL4(1)= 0.0
            HL4(2)= 0.0
            HL4(3)= 0.0 
            HL4(4)= 4.0*X4-1.0
            HL4(5)= 0.0
            HL4(6)= 0.0
            HL4(7)= 0.0
            HL4(8)= 4.0*X1
            HL4(9)= 4.0*X2
            HL4(10)= 4.0*X3
!  JACOBI MATRIX 
            XJ11=0.0
            XJ21=0.0
            XJ31=0.0
            XJ12=0.0
            XJ22=0.0
            XJ32=0.0
            XJ13=0.0
            XJ23=0.0
            XJ33=0.0
            DO I=1,NN
              XJ11=XJ11+(HL1(I)-HL4(I))*XX(I)
              XJ21=XJ21+(HL2(I)-HL4(I))*XX(I)
              XJ31=XJ31+(HL3(I)-HL4(I))*XX(I)
              XJ12=XJ12+(HL1(I)-HL4(I))*YY(I)
              XJ22=XJ22+(HL2(I)-HL4(I))*YY(I)
              XJ32=XJ32+(HL3(I)-HL4(I))*YY(I)
              XJ13=XJ13+(HL1(I)-HL4(I))*ZZ(I)
              XJ23=XJ23+(HL2(I)-HL4(I))*ZZ(I)
              XJ33=XJ33+(HL3(I)-HL4(I))*ZZ(I)
            ENDDO
!DETERMINANT OF JACOBIAN
            DET=XJ11*XJ22*XJ33        &
                 +XJ12*XJ23*XJ31      &
                 +XJ13*XJ21*XJ32      &
                 -XJ13*XJ22*XJ31      &
                 -XJ12*XJ21*XJ33      &
                 -XJ11*XJ23*XJ32
! INVERSION OF JACOBIAN
            DUM=1.0/DET
            XJI11=DUM*( XJ22*XJ33-XJ23*XJ32)
            XJI21=DUM*(-XJ21*XJ33+XJ23*XJ31)
            XJI31=DUM*( XJ21*XJ32-XJ22*XJ31)
            XJI12=DUM*(-XJ12*XJ33+XJ13*XJ32)
            XJI22=DUM*( XJ11*XJ33-XJ13*XJ31)
            XJI32=DUM*(-XJ11*XJ32+XJ12*XJ31)
            XJI13=DUM*( XJ12*XJ23-XJ13*XJ22)
            XJI23=DUM*(-XJ11*XJ23+XJ13*XJ21)
            XJI33=DUM*( XJ11*XJ22-XJ12*XJ21)
!  
            DET=-DET
!  WEIGT VALUE AT GAUSSIAN POINT
            WG=WGT(NG,L1)*WGT(NG,L2)*WGT(NG,L3)*DET      &
               *(1.0-X3)*(1.0-X2-X3)*0.125
!* DUDX DUDY DUDZ ....
            DUDX=0.0
            DUDY=0.0
            DUDZ=0.0
            DVDX=0.0
            DVDY=0.0
            DVDZ=0.0
            DWDX=0.0
            DWDY=0.0
            DWDZ=0.0
            DO J=1,NN
              DNDX=XJI11*(HL1(J)-HL4(J)) +XJI12*(HL2(J)-HL4(J)) +XJI13*(HL3(J)-HL4(J))
              DNDY=XJI21*(HL1(J)-HL4(J)) +XJI22*(HL2(J)-HL4(J)) +XJI23*(HL3(J)-HL4(J))
              DNDZ=XJI31*(HL1(J)-HL4(J)) +XJI32*(HL2(J)-HL4(J)) +XJI33*(HL3(J)-HL4(J))
              DUDX=DUDX+DNDX*UX(J)
              DUDY=DUDY+DNDY*UX(J)
              DUDZ=DUDZ+DNDZ*UX(J)
              DVDX=DVDX+DNDX*UY(J)
              DVDY=DVDY+DNDY*UY(J)
              DVDZ=DVDZ+DNDZ*UY(J)
              DWDX=DWDX+DNDX*UZ(J)
              DWDY=DWDY+DNDY*UZ(J)
              DWDZ=DWDZ+DNDZ*UZ(J)
            ENDDO
!** Fill BMAT
            HEAD=1
            DO J=1,NN
!* DNDX DNDY DNDZ
              DNDX=XJI11*(HL1(J)-HL4(J)) +XJI12*(HL2(J)-HL4(J)) +XJI13*(HL3(J)-HL4(J))
              DNDY=XJI21*(HL1(J)-HL4(J)) +XJI22*(HL2(J)-HL4(J)) +XJI23*(HL3(J)-HL4(J))
              DNDZ=XJI31*(HL1(J)-HL4(J)) +XJI32*(HL2(J)-HL4(J)) +XJI33*(HL3(J)-HL4(J))
!** BL0
              BL0(1,HEAD )=DNDX
              BL0(2,HEAD )=0.0
              BL0(3,HEAD )=0.0
              BL0(4,HEAD )=DNDY
              BL0(5,HEAD )=0.0
              BL0(6,HEAD )=DNDZ
              BL0(1,HEAD+1)=0.0
              BL0(2,HEAD+1)=DNDY
              BL0(3,HEAD+1)=0.0
              BL0(4,HEAD+1)=DNDX
              BL0(5,HEAD+1)=DNDZ
              BL0(6,HEAD+1)=0.0
              BL0(1,HEAD+2)=0.0
              BL0(2,HEAD+2)=0.0
              BL0(3,HEAD+2)=DNDZ
              BL0(4,HEAD+2)=0.0
              BL0(5,HEAD+2)=DNDY
              BL0(6,HEAD+2)=DNDX
!** BL1
              BL1(1,HEAD )=DUDX*DNDX
              BL1(2,HEAD )=DUDY*DNDY
              BL1(3,HEAD )=DUDZ*DNDZ
              BL1(4,HEAD )=DUDY*DNDX+DUDX*DNDY
              BL1(5,HEAD )=DUDY*DNDZ+DUDZ*DNDY
              BL1(6,HEAD )=DUDZ*DNDX+DUDX*DNDZ
              BL1(1,HEAD+1)=DVDX*DNDX
              BL1(2,HEAD+1)=DVDY*DNDY
              BL1(3,HEAD+1)=DVDZ*DNDZ
              BL1(4,HEAD+1)=DVDY*DNDX+DVDX*DNDY
              BL1(5,HEAD+1)=DVDY*DNDZ+DVDZ*DNDY
              BL1(6,HEAD+1)=DVDZ*DNDX+DVDX*DNDZ
              BL1(1,HEAD+2)=DWDX*DNDX
              BL1(2,HEAD+2)=DWDY*DNDY
              BL1(3,HEAD+2)=DWDZ*DNDZ
              BL1(4,HEAD+2)=DWDY*DNDX+DWDX*DNDY
              BL1(5,HEAD+2)=DWDY*DNDZ+DWDZ*DNDY
              BL1(6,HEAD+2)=DWDZ*DNDX+DWDX*DNDZ
              HEAD=HEAD+NDOF
            ENDDO

!** DE & DS ***
!!
      TEMPC=0.0
      TEMPI=0.0
      DO J=1,NN
        TEMPC=TEMPC+TT (J)*H(J)
        TEMPI=TEMPI+TTI(J)*H(J)
      END DO  

!     EE  = ccc_rl(1)*dexp(-ccc_rl(2)*TEMPC)*ccc_rl(3)*dexp(ccc_rl(4)*CURE3D(L1,L2,L3))
      EE  = ccc_rl(1)*dexp(-ccc_rl(2)*TEMPI)*ccc_rl(3)*dexp(ccc_rl(4)*CURE3D_I(L1,L2,L3))
!!      EE = ccc_rl(1)
      ALP =   (ccc_rl(5)*TEMPI+ccc_rl(6))*(1.0-CURE3D_I(L1,L2,L3))     &
            + (ccc_rl(7)*TEMPI+ccc_rl(8))*CURE3D_I(L1,L2,L3) 
      ALP = ALP*(TEMPC-TEMPI) - ramd_chm*(CURE3D(L1,L2,L3)-CURE3D_I(L1,L2,L3))
      THERMAL_EPS = ALP/3.0
!!
      call fstr_elpl_DMAT(DMAT,STRESS(:,L1,L2,L3),FLG_ELPL(L1,L2,L3), h_dash, EE, PP  &
                         ,BIJ,MAT_OTHO1,MAT_ISO)
!!
      if(idx_elpl .eq. -1 .or. idx_elpl .eq. 2) then
        DO J=1,6
          DO K=1,NN*NDOF
            BL1(J,K)=0.0
          END DO
        END DO
      end if
!!
!** DE & DS ***
            DO J=1,6
              DS(J)=0.0
            ENDDO
!
              EPSTH(1)=THERMAL_EPS
              EPSTH(2)=THERMAL_EPS
              EPSTH(3)=THERMAL_EPS
              EPSTH(4)=0.0
              EPSTH(5)=0.0
              EPSTH(6)=0.0
!**
!** SET DSGM  {ds}=[D]{de_therm}
!**
            DO J=1,6
              DO K=1,6
                DS(K)=DS(K)+DMAT(K,J)*EPSTH(J)
              ENDDO
            ENDDO
!**
!** CALCULATE LOAD {dF}=[B]T{de_therm}
!**
            DO J=1,NN*NDOF
              DO K=1,6
                VECT(J)=VECT(J)+(BL0(K,J)+BL1(K,J))*DS(K)*WG
              ENDDO
            ENDDO

!ogasawara till next ogasawara
!
          ENDDO
        ENDDO
      ENDDO
      RETURN
      END SUBROUTINE TLOAD_C3D10TL

end module m_static_LIB_3d2TL
