!======================================================================!
!                                                                      !
!    Software Name : FrontCOMP_cure   Ver. 3.1                         !
!                                                                      !
!      Module Name : m_fstr_solve_NLGEOM                               !
!      Category    : Static Analysis                                   !
!                                                                      !
!      Developed based on "FrontSTR" of RSS21 project                  !
!                                                                      !
!                     Written by Toshio Nagashima,   2007/07/04        !
!                                Yasuji Fukahori,    2007/07/04        !
!                                Tomotaka Ogasawara, 2013/03/26        !
!                                                                      !
!     Contact address :  IIS,The University of Tokyo, CISS             !
!                                                                      !
!    "Composite Material Strength & Reliability Evaluation Simulator"  !
!                                                                      !
!======================================================================!


!C================================================================C
!C-- module m_fstr_solve_NLGEOM
!C================================================================C

module m_fstr_solve_NLGEOM

use m_static_lib
use m_static_make_result
use m_static_get_prop

use m_hecmw2fstr_mesh_conv

contains

!C***
!C*** SOLVE NONLINEAR SOLID MECHANICS (NLGEOM)
!C***
      subroutine FSTR_SOLVE_NLGEOM(hecMESH,hecMAT,fstrSOLID,fstrRESULT,fstrPARAM)

      use m_fstr
      use hecmw

      implicit none

      type (hecmwST_local_mesh) :: hecMESH
      type (hecmwST_matrix    ) :: hecMAT
      type (fstr_param       )  :: fstrPARAM
      type (fstr_solid       )  :: fstrSOLID
      type (hecmwST_result_data):: fstrRESULT
!C
!C Local Variable
!C
!C=====================================================  !!ogasawra til next !ogasawara
!C
!C----  fstrSOLID%STRESS_YLD(n_gauss,n_gauss,n_gauss,neall)   : yield stress
!C----  fstrSOLID%STRAIN3D_P(6,n_gauss,n_gauss,n_gauss,neall) : total strain (not plastic)
!C----  fstrSOLID%EP_BAR(n_gauss,n_gauss,n_gauss,neall) : equivalent plastic strain
!C----  fstrSOLID%FLG_ELPL  (n_gauss,n_gauss,n_gauss,neall) : flag of plasticity
!C                                                            0:elastic  1: plastic
!C----  fstrSOLID%FACTOR(2)                              : factor of incrementation
!C                                                         1:time t  2: time t+dt
!C
!C                               restart_nout : output interval of restart file
!C                                              (if  .gt.0) restart file write  
!C                                              (if  .lt.0) restart file read and write 
!C                               idx_elpl     : type of analysis
!C                                             -1: elastic (geomet. linear   )
!C                                              0: elastic (geomet. nonlinear)
!C                                              1: plastic (geomet. nonlinear)
!C                                              2: plastic (geomet. linear   )
!C                               t_start      : start time of analysis (not used)
!C                               t_end = ETIME: end time of analysis
!C                               n_step =INCMAX: total step number of analysis
!C                               sig_y0       : initial yield stress
!C                               h_dash       : hardening coefficient
!C                               nout         : output interval of result
!C                               nout_monit   : output interval of result monitoring
!C                               node_monit_1 : node of monitoring result
!C                               elem_monit_1 : element of monitoring result
!C                               intg_monit_1 : integration point of monitoring result
!C                               iout_list    : 0: not output  1: output
!!
      integer(kind=kint) :: ndof, nn, ntotal
      integer(kind=kint) :: stat, ios2 
      integer(kind=kint),parameter :: elpl_IR1        =   82
      integer(kind=kint),parameter :: elpl_IW1        =   83
      integer(kind=kint),parameter :: elpl_IW4        =   84
      integer(kind=kint),parameter :: elpl_IW5        =   85

      integer(kind=kint) :: i, j, ii2, jj2, ierror ,istep, substep, intvl1, iter
      integer(kind=kint) :: n_gauss, g1, g2, g3 ,c1 ,iS ,iE ,icount2, ic_type 
      integer(kind=kint) :: itype2, nn2, icel, nodLocal(20)  
      integer(kind=kint) :: i_step, ig0, flag_u, i_flg_1,i_flg_2
      integer(kind=kint) :: my_rank_monit_1, my_rank_monit_2
      real(kind=kreal) :: tt, factor, res, res0, DT0
      real(kind=kreal) :: time_1, time_2
      real(kind=kreal) :: st_1(10), ep_1(7), xx(12)

      integer(kind=kint) :: restrt_step_num, nodes_all_area 
      integer(kind=kint),pointer :: restrt_step(:)
      character(len=30)  :: cdummy

      hecMAT%NDOF = hecMESH%n_dof

      ndof = hecMAT%NDOF
      nn = ndof*ndof

!----------------------------------------------------------------------
!!
!!--- read static parameter temp
!!
!    if( hecMESH%my_rank .eq. 0) then
!--- resin
!--- fiber
!    end if
!!
!!--- read static parameter temp
!!

!!
!!-- initial value
!!
      my_rank_monit_1 = 999999999
      my_rank_monit_2 = 999999999
      i_flg_1 = 0
      i_flg_2 = 0

      if( i_flg_1 .ne. 1 ) then
        do i= 1, hecMESH%nn_internal
        if( fstrPARAM%global_local_id(1,i) .eq. fstrSOLID%node_monit_1 ) then
          fstrSOLID%node_monit_1 = i
          i_flg_1 = 1
          do j= 1, hecMESH%nn_internal
!!          if( hecMESH%node_ID(2*j-1) .eq. i ) then
            if( j .eq. i ) then
              my_rank_monit_1 = hecMESH%node_ID(2*j)
              exit
            end if
          end do
          exit
        end if
        enddo
      end if

      if( i_flg_2 .ne. 1 ) then
        do i= 1, hecMESH%n_elem
!        do i= 1, hecMESH%ne_internal
        if( hecMESH%global_elem_ID(i) .eq. fstrSOLID%elem_monit_1 ) then
          fstrSOLID%elem_monit_1 = i
          i_flg_2 = 1
              my_rank_monit_2 = hecMESH%elem_ID(2*i)
          exit
        end if
        end do
      end if
!
      if( fstrSOLID%restart_nout .eq. 0 ) then
        fstrSOLID%restart_nout = 999999999
      end if
!
      if( fstrSOLID%sub_intvl .gt. 0 ) then
        intvl1 = fstrSOLID%sub_intvl 
      else if( fstrSOLID%sub_intvl .eq. 0 ) then
        intvl1 = 1 
      else if( fstrSOLID%sub_intvl .lt. 0 ) then
        intvl1 = 1 
      end if
!

!!
!C
!C-- file open for local use
!C
      if( hecMESH%my_rank .eq. my_rank_monit_1) then
        OPEN(elpl_IW4,FILE='elpl_disp_p1.out', status = 'replace', iostat=stat)
        if( stat /= 0 ) then
          write(*,*) 'stop due to file opening error <elpl_disp_p1.>'
          call hecmw_abort( hecmw_comm_get_comm())
        end if
      endif
      if( hecMESH%my_rank .eq. my_rank_monit_2) then
        OPEN(elpl_IW5,FILE='elpl_strs_p1.out', status = 'replace', iostat=stat)
        if( stat /= 0 ) then
          write(*,*) 'stop due to file opening error <elpl_strs_p1.>'
          call hecmw_abort( hecmw_comm_get_comm())
        end if
      endif

!!
!C
!C Allocate Matrix and Vector
!C
      allocate (hecMAT%AL(nn*hecMAT%NPL)        ,STAT=ierror )
          if( ierror /= 0 ) then
            write(idbg,*) 'stop due to allocation error <FSTR_SOLVE_NLGEOM, AL>'
            write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
            call flush(idbg)
            call hecmw_abort( hecmw_comm_get_comm())
          end if
      allocate (hecMAT%AU(nn*hecMAT%NPU)        ,STAT=ierror )
          if( ierror /= 0 ) then
            write(idbg,*) 'stop due to allocation error <FSTR_SOLVE_NLGEOM, AU>'
            write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
            call flush(idbg)
            call hecmw_abort( hecmw_comm_get_comm())
          end if
      allocate (hecMAT%B(ndof*hecMAT%NP)          ,STAT=ierror )
          if( ierror /= 0 ) then
            write(idbg,*) 'stop due to allocation error <FSTR_SOLVE_NLGEOM, B>'
            write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
            call flush(idbg)
            call hecmw_abort( hecmw_comm_get_comm())
          end if
      allocate (hecMAT%D(nn*hecMAT%NP)          ,STAT=ierror )
          if( ierror /= 0 ) then
            write(idbg,*) 'stop due to allocation error <FSTR_SOLVE_NLGEOM, D>'
            write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
            call flush(idbg)
            call hecmw_abort( hecmw_comm_get_comm())
          end if
      allocate (hecMAT%X(ndof*hecMAT%NP)          ,STAT=ierror )
          if( ierror /= 0 ) then
            write(idbg,*) 'stop due to allocation error <FSTR_SOLVE_NLGEOM, X>'
            write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
            call flush(idbg)
            call hecmw_abort( hecmw_comm_get_comm())
          end if
      allocate (hecMAT%ALU(nn*hecMAT%N)         ,STAT=ierror )
          if( ierror /= 0 ) then
            write(idbg,*) 'stop due to allocation error <FSTR_SOLVE_NLGEOM, ALU>'
            write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
            call flush(idbg)
            call hecmw_abort( hecmw_comm_get_comm())
          end if
      ntotal=ndof*hecMAT%NP
!C
!C Allocate for NLGEOM analysis
!C
      allocate ( fstrSOLID%GL( ntotal )          ,STAT=ierror ) 
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <FSTR_SOLVE_NLGEOM, GL>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      allocate ( fstrSOLID%TOTAL_DISP( ntotal )  ,STAT=ierror )
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <FSTR_SOLVE_NLGEOM, TOTAL_DISP>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      allocate ( fstrSOLID%QFORCE( ntotal )      ,STAT=ierror ) 
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <FSTR_SOLVE_NLGEOM, QFORCE>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      allocate ( fstrSOLID%QFORCE_RES( ntotal )      ,STAT=ierror ) 
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <FSTR_SOLVE_NLGEOM, QFORCE_RES>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      allocate ( fstrSOLID%FACTOR( 2 )           ,STAT=ierror ) 
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <FSTR_SOLVE_NLGEOM, FACTOR>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
!!! restart  !!!
      allocate( restrt_step(1)                   ,STAT=ierror )
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <FSTR_SOLVE_NLGEOM, restrt_step>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
!
!!!   20111021  ----
      fstrSOLID%QFORCE    =0.0d0
      fstrSOLID%QFORCE_RES=0.0d0

!!--- 20120927  --------
!! Read covariant base vector in material coord. at node (all) 
!!---  FrontSTR_FRP_29_A --------

    if ( fstrSOLID%mat_iso_num(1) .ne. -1) then
        allocate(fstrSOLID%GBARU_ALL(3,3,1)  ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_InitialState, GBARU_ALL>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
    else if ( fstrSOLID%mat_iso_num(1) .eq. -1) then
      if( intvl1 .gt. -99999 ) then
         read (matl_vec_IW1,*)  cdummy 
         write(idbg        ,*)  cdummy
         read (matl_vec_IW1,*)  nodes_all_area 
         write(idbg        ,*)  nodes_all_area

        allocate(fstrSOLID%GBARU_ALL(3,3,nodes_all_area)  ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_InitialState, GBARU_ALL>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        fstrSOLID%GBARU_ALL=0.0d0
        xx                 =0.0d0

        do icount2= 1, nodes_all_area 
         read (matl_vec_IW1,*)  ((fstrSOLID%GBARU_ALL(i,j,icount2),j=1,3),i=1,3),(xx(i),i=1,3) 
         write(idbg        ,*)  icount2
         write(idbg        ,*)  ((fstrSOLID%GBARU_ALL(i,j,icount2),j=1,3),i=1,3),(xx(i),i=1,3) 
        end do
      end if
!!---  FrontSTR_FRP_29_B --------

      if( intvl1 .eq. -99999 ) then
        xx                 = 0.0D+00
        read (matl_vec_IW1,*)  cdummy

        nodes_all_area = 0
        do
         read (matl_vec_IW1,*,iostat=ios2) (xx(i),i=1,12)
         if(ios2 .lt. 0) exit
         nodes_all_area = nodes_all_area + 1
        end do
        rewind(matl_vec_IW1)

        allocate(fstrSOLID%GBARU_ALL(3,3,nodes_all_area)  ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_InitialState, GBARU_ALL>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        fstrSOLID%GBARU_ALL=0.0d0
        xx                 =0.0d0

         read (matl_vec_IW1,*)  cdummy
         write(idbg        ,*)  cdummy
         write(idbg        ,*)  nodes_all_area
        do icount2= 1, nodes_all_area 
         read (matl_vec_IW1,*)  ((fstrSOLID%GBARU_ALL(i,j,icount2),j=1,3),i=1,3),(xx(i),i=1,3) 
         write(idbg        ,*)  icount2
         write(idbg        ,*)  ((fstrSOLID%GBARU_ALL(i,j,icount2),j=1,3),i=1,3),(xx(i),i=1,3) 
        end do
      end if
    end if
!     return
!!--------------------------------------------------------------
!C Initialize 

      call fstrNLGEOM_InitialState(hecMESH,hecMAT,fstrSOLID)

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

!! Read Temperature data

    if(fstrSOLID%idx_temp .eq. 1 ) then
!!   if( hecMESH%my_rank .eq. 0) then
     call hecmw_update_1_R ( hecMESH,fstrSOLID%TEMPI        ,hecMESH%n_node )
     call hecmw_update_1_R ( hecMESH,fstrSOLID%CURE3D_nod_I ,hecMESH%n_node )
     call hecmw_update_1_R ( hecMESH,fstrSOLID%TEMPC      ,hecMESH%n_node )
     call hecmw_update_1_R ( hecMESH,fstrSOLID%CURE3D_nod ,hecMESH%n_node )

      call fstrNLGEOM_ReadTemp(hecMESH,hecMAT,fstrSOLID,fstrPARAM    &
                              ,heat_sta_IW1,heat_sta_IW2,i_step, DT0)

     call hecmw_update_1_R ( hecMESH,fstrSOLID%TEMPI        ,hecMESH%n_node )
     call hecmw_update_1_R ( hecMESH,fstrSOLID%CURE3D_nod_I ,hecMESH%n_node )
     call hecmw_update_1_R ( hecMESH,fstrSOLID%TEMPC      ,hecMESH%n_node )
     call hecmw_update_1_R ( hecMESH,fstrSOLID%CURE3D_nod ,hecMESH%n_node )
!!   end if
    end if

!---------------------------------------------
!!! restart  !!!
      restrt_step(1) = 1
      restrt_step_num = 1

    if(fstrSOLID%restart_nout .lt. 0 ) then
      n_gauss=2
          ic_type= hecMESH%elem_type_item(1)
      if( hecMESH%n_dof .eq. 3 ) then
          if (ic_type.eq.342.or.ic_type.eq.352.or.ic_type.eq.362) then 
            n_gauss=3
          end if
      else if( hecMESH%n_dof .eq. 2 ) then
          if (ic_type.eq.242 .or. ic_type.eq.232 ) then
            n_gauss=3
          end if
      end if
      call hecmw_restart_open()
      call hecmw_restart_read_int(restrt_step)
      call hecmw_restart_read_real(fstrSOLID%TOTAL_DISP)
!     call hecmw_restart_read_real(fstrSOLID%QFORCE)
      call hecmw_restart_read_real(fstrSOLID%GL)
      if ( hecMESH%n_dof .eq. 3 ) then
        do g1= 1, n_gauss
         do g2= 1, n_gauss
          do g3= 1, n_gauss
            call hecmw_restart_read_real(fstrSOLID%EP_BAR(g1,g2,g3,:))
            call hecmw_restart_read_real(fstrSOLID%STRESS_YLD(g1,g2,g3,:))
            call hecmw_restart_read_int(fstrSOLID%FLG_ELPL(g1,g2,g3,:))
            do c1= 1, 6
              call hecmw_restart_read_real(fstrSOLID%STRESS3D(c1,g1,g2,g3,:))
              call hecmw_restart_read_real(fstrSOLID%STRAIN3D_P(c1,g1,g2,g3,:))
             if(fstrSOLID%idx_visco .ge. 4) then
              call hecmw_restart_read_real(fstrSOLID%STRAIN3D_V1(c1,g1,g2,g3,:))
             end if
             if(fstrSOLID%idx_visco .eq. 6) then
              call hecmw_restart_read_real(fstrSOLID%STRAIN3D_V2(c1,g1,g2,g3,:))
             end if
            end do
          end do
         end do
        end do
      else if( hecMESH%n_dof .eq. 2 ) then
        do g1= 1, n_gauss
         do g2= 1, n_gauss
            call hecmw_restart_read_real(fstrSOLID%EP_BAR(g1,g2,1,:))
            call hecmw_restart_read_real(fstrSOLID%STRESS_YLD(g1,g2,1,:))
            call hecmw_restart_read_int(fstrSOLID%FLG_ELPL(g1,g2,1,:))
            do c1= 1, 4
              call hecmw_restart_read_real(fstrSOLID%STRESS2D(c1,g1,g2,:))
              call hecmw_restart_read_real(fstrSOLID%STRAIN3D_P(c1,g1,g2,1,:))
            end do
         end do
        end do
      end if

      call hecmw_restart_close()
      restrt_step_num = restrt_step(1) + 1
    end if

!C
!C-- time step loop
!C
      i_step  = 0
      substep = 0
      iter    = 0
      DT0 = DT
      tt = (restrt_step_num -1)*DT0
!!
!!    step = 0
!!

!C-- output initial condition
      if( restrt_step_num .eq. 1 ) then
        write(ILOG,'(a,i10,5x,a,1pe12.4)') 'STEP =',0,'TIME =',0.
        call fstrNLGEOM_Output(hecMESH,hecMAT,fstrSOLID,tt,fstrPARAM,DT0,iter, i_step)
!C-- output result of monitoring node
        if( hecMESH%n_dof .eq. 3 ) then
         call fstrNLGEOM_output_monit   (hecMESH, fstrPARAM, fstrSOLID, elpl_IW4, elpl_IW5, i_step,  &
                                         DT0, my_rank_monit_1, my_rank_monit_2)
        else if( hecMESH%n_dof .eq. 2 ) then
         call fstrNLGEOM_output_monit_2d(hecMESH, fstrPARAM, fstrSOLID, elpl_IW4, elpl_IW5, i_step,  &
                                         DT0, my_rank_monit_1, my_rank_monit_2)
        end if
      end if

!!--------------------------------------------------------------
!!
!!    step = 1,2,....,n_step
!!

      do i_step = restrt_step_num, INCMAX
!
        time_1 =  hecmw_Wtime()
!
!     do while ( tt. lt. ETIME-0.5*DT )

!!    time history of factor
      ig0 = 1
      flag_u = 1
      tt=tt+DT

        call table_nlsta(hecMESH,fstrSOLID,i_step-1,ig0,factor,flag_u)
          fstrSOLID%FACTOR(1) = factor
        call table_nlsta(hecMESH,fstrSOLID,i_step,ig0,factor,flag_u)
          fstrSOLID%FACTOR(2) = factor
!!
!! Read Temperature data
!!
   if(fstrSOLID%idx_temp .eq. 1 ) then
     call hecmw_update_1_R ( hecMESH,fstrSOLID%TEMPI        ,hecMESH%n_node )
     call hecmw_update_1_R ( hecMESH,fstrSOLID%CURE3D_nod_I ,hecMESH%n_node )
     call hecmw_update_1_R ( hecMESH,fstrSOLID%TEMPC      ,hecMESH%n_node )
     call hecmw_update_1_R ( hecMESH,fstrSOLID%CURE3D_nod ,hecMESH%n_node )
     call hecmw_update_1_R ( hecMESH,fstrSOLID%DELTEMP       ,hecMESH%n_node )
     call hecmw_update_1_R ( hecMESH,fstrSOLID%DELCURE3D_nod ,hecMESH%n_node )
!
     do i= 1, hecMESH%nn_internal
       fstrSOLID%TEMPI(i) = fstrSOLID%TEMPC(i)
     enddo
!
     do i= 1, hecMESH%nn_internal
       fstrSOLID%CURE3D_nod_I(i) = fstrSOLID%CURE3D_nod(i)
     enddo
!
!!  if( hecMESH%my_rank .eq. 0) then
     call hecmw_update_1_R ( hecMESH,fstrSOLID%TEMPI        ,hecMESH%n_node )
     call hecmw_update_1_R ( hecMESH,fstrSOLID%CURE3D_nod_I ,hecMESH%n_node )
     call hecmw_update_1_R ( hecMESH,fstrSOLID%TEMPC      ,hecMESH%n_node )
     call hecmw_update_1_R ( hecMESH,fstrSOLID%CURE3D_nod ,hecMESH%n_node )

    if( mod((i_step-1),intvl1) == 0 ) then

      if( fstrSOLID%sub_intvl .ge. 0 ) then
       call fstrNLGEOM_ReadTemp(hecMESH,hecMAT,fstrSOLID,fstrPARAM    &
                               ,heat_sta_IW1,heat_sta_IW2,i_step, DT0)
      else if( fstrSOLID%sub_intvl .lt. 0 ) then
       do i= 1, - fstrSOLID%sub_intvl 
         call fstrNLGEOM_ReadTemp(hecMESH,hecMAT,fstrSOLID,fstrPARAM    &
                                 ,heat_sta_IW1,heat_sta_IW2,i_step, DT0)
       end do
      end if

      call hecmw_update_1_R ( hecMESH,fstrSOLID%TEMPI        ,hecMESH%n_node )
      call hecmw_update_1_R ( hecMESH,fstrSOLID%CURE3D_nod_I ,hecMESH%n_node )
      call hecmw_update_1_R ( hecMESH,fstrSOLID%TEMPC      ,hecMESH%n_node )
      call hecmw_update_1_R ( hecMESH,fstrSOLID%CURE3D_nod ,hecMESH%n_node )
!
      do i= 1, hecMESH%nn_internal
       fstrSOLID%DELCURE3D_nod(i) = (fstrSOLID%CURE3D_nod(i) - fstrSOLID%CURE3D_nod_I(i))/intvl1 
       fstrSOLID%DELTEMP(i)       = (fstrSOLID%TEMPC(i)      - fstrSOLID%TEMPI(i))       /intvl1 
      enddo

      call hecmw_update_1_R ( hecMESH,fstrSOLID%TEMPI        ,hecMESH%n_node )
      call hecmw_update_1_R ( hecMESH,fstrSOLID%CURE3D_nod_I ,hecMESH%n_node )
      call hecmw_update_1_R ( hecMESH,fstrSOLID%TEMPC        ,hecMESH%n_node )
      call hecmw_update_1_R ( hecMESH,fstrSOLID%CURE3D_nod   ,hecMESH%n_node )
      call hecmw_update_1_R ( hecMESH,fstrSOLID%DELTEMP       ,hecMESH%n_node )
      call hecmw_update_1_R ( hecMESH,fstrSOLID%DELCURE3D_nod ,hecMESH%n_node )
    end if
!
      do i= 1, hecMESH%nn_internal
        fstrSOLID%CURE3D_nod(i) = fstrSOLID%CURE3D_nod_I(i) + fstrSOLID%DELCURE3D_nod(i)
        fstrSOLID%TEMPC(i)      = fstrSOLID%TEMPI(i)        + fstrSOLID%DELTEMP(i) 
      enddo
!
      call hecmw_update_1_R ( hecMESH,fstrSOLID%TEMPI        ,hecMESH%n_node )
      call hecmw_update_1_R ( hecMESH,fstrSOLID%CURE3D_nod_I ,hecMESH%n_node )
      call hecmw_update_1_R ( hecMESH,fstrSOLID%TEMPC        ,hecMESH%n_node )
      call hecmw_update_1_R ( hecMESH,fstrSOLID%CURE3D_nod   ,hecMESH%n_node )
      call hecmw_update_1_R ( hecMESH,fstrSOLID%DELTEMP       ,hecMESH%n_node )
      call hecmw_update_1_R ( hecMESH,fstrSOLID%DELCURE3D_nod ,hecMESH%n_node )
!!  end if
   end if

!C
!C Set Load Vector 
!C
      call fstrNLGEOM_SetLoad(hecMESH,hecMAT,fstrSOLID,iter,i_step)

!!--------------------------------------------------------------
!C
!C Newton Raphson Iteration
!C
        do iter=1,ITMAX
          write(ILOG,*) '###ITERATION:',iter
!C
!C Set Internal force & stiffness
!C
          call fstrNLGEOM_SetKQ(hecMESH,hecMAT,fstrSOLID,DT0,iter,i_step)
!C
!C Set residual
!C
          call fstrNLGEOM_SetResidual(hecMESH,hecMAT,fstrSOLID,factor,res )
!C
!C Gather global residual
!C
          call hecmw_allREDUCE_R1(hecMESH,res,hecmw_sum)

        if(fstrSOLID%idx_temp .eq. 0 ) then
         if(fstrSOLID%idx_visco .lt. 2 ) then
           res=sqrt(res)
           if( iter.eq.1 ) then
             res0=res
           endif
         else if(fstrSOLID%idx_visco .ge. 2 ) then
           res=sqrt(res)
           if( iter.le.1 ) then
             res0=res
           endif
         end if
        else if(fstrSOLID%idx_temp .eq. 1 ) then
           res=sqrt(res)
           if( iter.le.1 ) then
             res0=res
           endif
        end if
!C
!C check convergence
!C
          if( mod(i_step,fstrSOLID%nout_monit) == 0 ) then
           if( hecMESH%my_rank.eq.0) then
            write(ISTA,'(''iter='',I5,''  res/res0='',E15.7,''  res='',E15.7)')iter,res/res0, res
           endif
          end if
!!!!  20111024  -------
         if(fstrSOLID%idx_visco .lt. 2 ) then
          if(i_step .ge. 1 .and. iter .eq. 1) then
            if( res .lt. res0*EPS ) exit
          else if(i_step .ge. 1 .and. iter .ge. 2) then
            if( res .lt. res0*EPS ) exit
          end if
         else if(fstrSOLID%idx_visco .ge. 2 ) then
          if(i_step .ge. 1 .and. iter .eq. 1) then
            if( res .lt. res0*EPS ) exit
          else if(i_step .ge. 1 .and. iter .ge. 2) then
            if( res .lt. res0*EPS ) exit
            if( res .lt. EPS ) exit
          end if
         end if

!C
!C Solve System Eqs. using ILU method
!C
          call fstrNLGEOM_AddBC(hecMESH,hecMAT,fstrSOLID,iter)
!C
!C-- SOLVE [Kt]{du}={R} 
!C
!!
      if( i_step .eq. restrt_step_num .and. iter .eq. 1 ) then
        hecMAT%Iarray(98) = 1   !Assmebly complete
      end if

      if( iter .eq. 1 ) then
        hecMAT%Iarray(97) = 1   !Need numerical factorization
      else if( iter .gt. 1 ) then
        hecMAT%Iarray(97) = 1   !Need numerical factorization
      end if
!
        CALL solve_LINEQ(hecMESH,hecMAT,imsg)

          if( hecMESH%n_dof .eq. 3 ) then
             call hecmw_update_3_R (hecMESH, hecMAT%X, hecMAT%NP)
             if( hecMESH%my_rank.eq.0) then
             end if
           else if( hecMESH%n_dof.eq. 2 ) then
             call hecmw_update_2_R (hecMESH, hecMAT%X, hecMAT%NP)
             if( hecMESH%my_rank.eq.0) then
             end if
          endif

!C
!C update U,E,S
!C
          call fstrNLGEOM_Update(hecMESH,hecMAT,fstrSOLID,factor,iter,i_step,DT0)
!
        enddo
!C-- end of iteration loop

!    not convergence
!!!!  20111024  -------
       if(fstrSOLID%idx_visco .lt. 2 ) then
         if( res .ge. res0*EPS ) then
           do i = 1,10
             write(ILOG,'("### Not Convergence")')
           end do
           if( hecMESH%my_rank.eq.0) then
             write(ISTA,'("### Not Convergence")')
             write(*   ,'("### Not Convergence")')
           end if
           exit
         end if
       else if(fstrSOLID%idx_visco .ge. 2 ) then
         if( res .ge. res0*EPS ) then
         if( res .ge. EPS ) then
           do i = 1,10
             write(ILOG,'("### Not Convergence")')
           end do
           if( hecMESH%my_rank.eq.0) then
             write(ISTA,'("### Not Convergence")')
             write(*   ,'("### Not Convergence")')
           end if
           exit
         end if
         end if
       end if

!    convergence
      if( mod(i_step,fstrSOLID%nout_monit) == 0 ) then
        if( hecMESH%my_rank.eq.0) then
         write(ISTA,'("### Convergence Achived : tt=",E12.4,"  iter=",I5)') tt,iter
        endif
      end if
!!
!!! restart  !!!
    if(fstrSOLID%restart_nout .lt. 0) then
      fstrSOLID%restart_nout = - fstrSOLID%restart_nout
    end if
    if( mod(i_step,fstrSOLID%restart_nout) == 0 ) then
      n_gauss=2
          ic_type= hecMESH%elem_type_item(1)
      if( hecMESH%n_dof .eq. 3 ) then
          if (ic_type.eq.342.or.ic_type.eq.352.or.ic_type.eq.362) then 
            n_gauss=3
          end if
      else if( hecMESH%n_dof .eq. 2 ) then
          if (ic_type.eq.242 .or. ic_type.eq.232 ) then
            n_gauss=3
          end if
      end if
      restrt_step(1) = i_step
      call hecmw_restart_add_int(restrt_step,size(restrt_step))
      call hecmw_restart_add_real(fstrSOLID%TOTAL_DISP,size(fstrSOLID%TOTAL_DISP))
      call hecmw_restart_add_real(fstrSOLID%GL,size(fstrSOLID%GL))

      if( hecMESH%n_dof .eq. 3 ) then
        do g1= 1, n_gauss
         do g2= 1, n_gauss
          do g3= 1, n_gauss
            call hecmw_restart_add_real(fstrSOLID%EP_BAR(g1,g2,g3,:),size(fstrSOLID%EP_BAR(g1,g2,g3,:)))
            call hecmw_restart_add_real(fstrSOLID%STRESS_YLD(g1,g2,g3,:),size(fstrSOLID%STRESS_YLD(g1,g2,g3,:)))
            call hecmw_restart_add_int(fstrSOLID%FLG_ELPL(g1,g2,g3,:),size(fstrSOLID%FLG_ELPL(g1,g2,g3,:)))
            do c1= 1, 6
              call hecmw_restart_add_real(fstrSOLID%STRESS3D(c1,g1,g2,g3,:),size(fstrSOLID%STRESS3D(c1,g1,g2,g3,:)))
              call hecmw_restart_add_real(fstrSOLID%STRAIN3D_P(c1,g1,g2,g3,:),size(fstrSOLID%STRAIN3D_P(c1,g1,g2,g3,:)))
             if(fstrSOLID%idx_visco .ge. 4) then
              call hecmw_restart_add_real(fstrSOLID%STRAIN3D_V1(c1,g1,g2,g3,:),size(fstrSOLID%STRAIN3D_V1(c1,g1,g2,g3,:)))
             end if
             if(fstrSOLID%idx_visco .eq. 6) then
              call hecmw_restart_add_real(fstrSOLID%STRAIN3D_V2(c1,g1,g2,g3,:),size(fstrSOLID%STRAIN3D_V2(c1,g1,g2,g3,:)))
             end if
            end do
          end do
         end do
        end do
      else if( hecMESH%n_dof .eq. 2 ) then
        do g1= 1, n_gauss
         do g2= 1, n_gauss
            call hecmw_restart_add_real(fstrSOLID%EP_BAR(g1,g2,1,:),size(fstrSOLID%EP_BAR(g1,g2,1,:)))
            call hecmw_restart_add_real(fstrSOLID%STRESS_YLD(g1,g2,1,:),size(fstrSOLID%STRESS_YLD(g1,g2,1,:)))
            call hecmw_restart_add_int(fstrSOLID%FLG_ELPL(g1,g2,1,:),size(fstrSOLID%FLG_ELPL(g1,g2,1,:)))
            do c1= 1, 4
              call hecmw_restart_add_real(fstrSOLID%STRESS2D(c1,g1,g2,:),size(fstrSOLID%STRESS2D(c1,g1,g2,:)))
              call hecmw_restart_add_real(fstrSOLID%STRAIN3D_P(c1,g1,g2,1,:),size(fstrSOLID%STRAIN3D_P(c1,g1,g2,1,:)))
            end do
         end do
        end do
      end if

      call hecmw_restart_write()
    end if
!C
!C-- output new displacement,strain and stress
!C
      if( mod(i_step,fstrSOLID%nout) == 0 ) then
        write(ILOG,'(a,i10,5x,a,1pe12.4)') 'STEP =',i_step,'TIME =',DT0*i_step
!C
!C Output results
!C
        call fstrNLGEOM_Output(hecMESH,hecMAT,fstrSOLID,tt,fstrPARAM,DT0,iter, i_step)
!C
!C Summary & Result
!C
        call fstrNLGEOM_Post(hecMESH,hecMAT,fstrSOLID,tt,i_step)
!
!C-- POST PROCESSING VIA MEMORY
       if( IVISUAL.eq. 1 ) then
!!---  20100203   ----------------
         call fstrNLGEOM_make_result(hecMESH,fstrSOLID,fstrRESULT)
           call fstr2hecmw_mesh_conv(hecMESH)
         call hecmw_visualize_init
!
         istep = i_step
         idummy=1
         call hecmw_visualize(hecMESH,fstrRESULT,istep,INCMAX,idummy)
         call hecmw_visualize_finalize
           call hecmw2fstr_mesh_conv(hecMESH)
         call hecmw_result_free(fstrRESULT)
       endif

      end if
!C
!C-- output result of monitoring node
      if( mod(i_step,fstrSOLID%nout_monit) == 0 ) then
        if( hecMESH%n_dof .eq. 3 ) then
         call fstrNLGEOM_output_monit   (hecMESH, fstrPARAM, fstrSOLID, elpl_IW4, elpl_IW5, i_step,  &
                                         DT0, my_rank_monit_1, my_rank_monit_2)
        else if( hecMESH%n_dof .eq. 2 ) then
         call fstrNLGEOM_output_monit_2d(hecMESH, fstrPARAM, fstrSOLID, elpl_IW4, elpl_IW5, i_step,  &
                                         DT0, my_rank_monit_1, my_rank_monit_2)
        end if
      end if
!!
       time_2 =  hecmw_Wtime()
!
      enddo
!C-- end of time step loop
!C
!C freeing allocated memory
!C
      deallocate(hecMAT%AL                  ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, AL>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      deallocate(hecMAT%AU                  ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, AU>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      deallocate(hecMAT%B                   ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, B>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      deallocate(hecMAT%D                   ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, D>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      deallocate(hecMAT%X                   ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, X>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      deallocate(hecMAT%ALU                 ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, ALU>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
!
      deallocate(fstrSOLID%GL               ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, GL>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      deallocate(fstrSOLID%TOTAL_DISP       ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, TOTAL_DISP>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      deallocate(fstrSOLID%QFORCE           ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, QFORCE>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      deallocate(fstrSOLID%QFORCE_RES           ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, QFORCE_RES>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      deallocate(fstrSOLID%FACTOR           ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, FACTOR>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
!!! restart  !!!
      deallocate( restrt_step               ,STAT=ierror )
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, restrt_step>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
!
      if( hecMESH%n_dof .eq. 3 ) then
        deallocate(fstrSOLID%STRESS3D       ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, STRESS3D>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        deallocate(fstrSOLID%SURFF       ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, SURFF>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        deallocate(fstrSOLID%STRESS3D_nd       ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, STRESS3D_nd>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        deallocate(fstrSOLID%STRESS3D_ndP1       ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, STRESS3D_ndP1>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        deallocate(fstrSOLID%STRESS_nd       ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, STRESS_nd>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        deallocate(fstrSOLID%STRESS_PRIS       ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, STRESS_PRIS>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        deallocate(fstrSOLID%SURFF_nd       ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, SURFF_nd>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        deallocate(fstrSOLID%STRSURF_nd       ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, STRSURF_nd>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        deallocate(fstrSOLID%STRAIN3D_P     ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, STRAIN3D_P>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        deallocate(fstrSOLID%STRAIN3D_P_nd     ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, STRAIN3D_P_nd>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        deallocate(fstrSOLID%STRAIN_nd     ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, STRAIN_nd>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        deallocate(fstrSOLID%STRAIN_PRIS     ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, STRAIN_PRIS>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        deallocate(fstrSOLID%EP_BAR         ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, EP_BAR>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        deallocate(fstrSOLID%STRESS_YLD     ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, STRESS_YLD>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        deallocate(fstrSOLID%FLG_ELPL       ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, FLG_ELPL>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        deallocate(fstrSOLID%mat_iso       ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, mat_iso>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        deallocate(fstrSOLID%CURE3D_nod         ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, CURE3D_nod>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        deallocate(fstrSOLID%CURE3D_nod_I       ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, CURE3D_nod_I>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        deallocate(fstrSOLID%DELCURE3D_nod         ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, DELCURE3D_nod>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        deallocate(fstrSOLID%TEMPI          ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, TEMPI>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        deallocate(fstrSOLID%TEMPC          ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, TEMPC>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        deallocate(fstrSOLID%DELTEMP          ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, DELTEMP>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        deallocate(fstrSOLID%GBARU          ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, GBARU>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        deallocate(fstrSOLID%GBARU_ALL      ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, GBARU_ALL>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
       if(fstrSOLID%idx_visco .ge. 4) then
        deallocate(fstrSOLID%STRAIN3D_V1     ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, STRAIN3D_V1>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
       end if
       if(fstrSOLID%idx_visco .eq. 6) then
        deallocate(fstrSOLID%STRAIN3D_V2     ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, STRAIN3D_V2>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
       end if

      else if( hecMESH%n_dof .eq. 2 ) then
        deallocate(fstrSOLID%STRESS2D       ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, STRESS2D>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        deallocate(fstrSOLID%STRAIN3D_P     ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, STRAIN3D_P>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        deallocate(fstrSOLID%EP_BAR         ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, EP_BAR>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        deallocate(fstrSOLID%STRESS_YLD     ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, STRESS_YLD>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        deallocate(fstrSOLID%FLG_ELPL       ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <FSTR_SOLVE_NLGEOM, FLG_ELPL>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      endif
!
!C-- finalization
!C-- file close for local use
      if( hecMESH%my_rank .eq. my_rank_monit_1) then
        CLOSE(elpl_IW4)
      endif
      if( hecMESH%my_rank .eq. my_rank_monit_2) then
        CLOSE(elpl_IW5)
      endif
!C-- end of finalization

!C
!C message
!C 
      IF(myrank .EQ. 0)THEN
       WRITE(IMSG,'("### FSTR_SOLVE_NLGEOM FINISHED!")')
      ENDIF
      return

      end subroutine FSTR_SOLVE_NLGEOM
!C---------------------------------------------------------------------*
      subroutine fstrNLGEOM_SetLoad(hecMESH, hecMAT, fstrSOLID, iter, i_step)
!C---------------------------------------------------------------------*
!C*********************
!C  Set Load  
!C*********************
      use m_fstr

      implicit none

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

!C Local Variable
      real(kind=kreal) xx(20),yy(20),zz(20)
      real(kind=kreal) params(0:6)
      real(kind=kreal) vect(60)
      integer(kind=kint) iwk(60)
      integer(kind=kint) nodLocal(20)
      integer(kind=kint) iter, count, ndof 
      integer(kind=kint) i_step 
      real(kind=kreal) tt(20),tt0(20),tti(20)
      real(kind=kreal),pointer:: temp(:)
      real(kind=kreal) disp(60), ddisp(60)

      integer(kind=kint) ig0,ig,ityp,iS0,iE0,ik,in,ltype,i,icel,ic_type 
      integer(kind=kint) nn,iS,iE,j,isect,nsize,itype,id,iset,ierror 
      real(kind=kreal) val,ee,pp,rho,alpha,thick,pa1 
      logical :: fg_surf

!C
!C CLEAR LOAD
!C 
!     fstrSOLID%GL=0.0d0

      if( hecMESH%n_dof .eq. 3 ) then
!C
!C CLOAD
!C 
        do ig0= 1, fstrSOLID%CLOAD_ngrp_tot
          ig= fstrSOLID%CLOAD_ngrp_ID(ig0)
          ityp= fstrSOLID%CLOAD_ngrp_DOF(ig0)
          val= fstrSOLID%CLOAD_ngrp_val(ig0)
!C START & END
          iS0= hecMESH%node_group%grp_index(ig-1) + 1
          iE0= hecMESH%node_group%grp_index(ig  )
          do ik= iS0, iE0
            in   = hecMESH%node_group%grp_item(ik)
            fstrSOLID%GL(3*(in-1)+ityp)=fstrSOLID%GL(3*(in-1)+ityp)        &
                                        +val*(fstrSOLID%FACTOR(2)-fstrSOLID%FACTOR(1))
          enddo
        enddo
!C
!C DLOAD
!C
        do ig0= 1, fstrSOLID%DLOAD_ngrp_tot
          ig= fstrSOLID%DLOAD_ngrp_ID(ig0)
          ltype= fstrSOLID%DLOAD_ngrp_LID(ig0)
          do i=0,6
            params(i)= fstrSOLID%DLOAD_ngrp_params(i,ig0)
          enddo
!C START & END
          fg_surf = (ltype == 100)
          if( fg_surf ) then                            ! surface group
            iS0= hecMESH%surf_group%grp_index(ig-1) + 1
            iE0= hecMESH%surf_group%grp_index(ig  )
          else                                          ! element group
            iS0= hecMESH%elem_group%grp_index(ig-1) + 1
            iE0= hecMESH%elem_group%grp_index(ig  )
          endif
!
          do ik= iS0, iE0

            if( fg_surf ) then                          ! surface group
              ltype  = hecMESH%surf_group%grp_item(2*ik)*10
              icel   = hecMESH%surf_group%grp_item(2*ik-1)
              ic_type= hecMESH%elem_type(icel)
            else                                        ! element group
              icel   = hecMESH%elem_group%grp_item(ik)
              ic_type= hecMESH%elem_type(icel)
            endif
!C** Create local stiffness		
            nn = hecmw_get_max_node(ic_type)
!C** node ID
            iS= hecMESH%elem_node_index(icel-1)
            do j=1,nn
              nodLOCAL(j)= hecMESH%elem_node_item (iS+j)
!C** nodal coordinate
              xx(j)=hecMESH%node(3*nodLOCAL(j)-2)
              yy(j)=hecMESH%node(3*nodLOCAL(j)-1)
              zz(j)=hecMESH%node(3*nodLOCAL(j)  )
!C** create iwk array ***
              iwk(3*(j-1)+1)=3*(nodLOCAL(j)-1)+1
              iwk(3*(j-1)+2)=3*(nodLOCAL(j)-1)+2
              iwk(3*(j-1)+3)=3*(nodLOCAL(j)-1)+3
            enddo  
!C** section  ID
            isect= hecMESH%section_ID(icel)
!C** Get Properties
            call fstrNLGEOM_get_prop(hecMESH,isect,ee,pp,rho,alpha,thick)
!C** Create local stiffness		
            if (ic_type.eq.341) then
              call DL_C3D4(xx,yy,zz,rho,ltype,params,vect,nsize)
            else if ( ic_type.eq.351) then
              call DL_C3D6(xx,yy,zz,rho,ltype,params,vect,nsize)
            else if ( ic_type.eq.361) then
              call DL_C3D8(xx,yy,zz,rho,ltype,params,vect,nsize)
            else if ( ic_type.eq.342) then
              call DL_C3D10(xx,yy,zz,rho,ltype,params,vect,nsize)
            else if ( ic_type.eq.352) then
              call DL_C3D15(xx,yy,zz,rho,ltype,params,vect,nsize)
            else if ( ic_type.eq.362) then
              call DL_C3D20(xx,yy,zz,rho,ltype,params,vect,nsize)
            endif
!C** Add vector
            do j=1,nsize
              fstrSOLID%GL( iwk(j) )=fstrSOLID%GL( iwk(j) )        &
                                    +vect(j)*(fstrSOLID%FACTOR(2)-fstrSOLID%FACTOR(1))
            enddo
          enddo
        enddo
!C
!C THERMAL LOAD USING TEMPERATURE  ( for one step analysis)
!C 
!C Set Temperature
!C 
        if( fstrSOLID%TEMP_ngrp_tot > 0 ) then
          allocate ( temp(hecMESH%n_node)       ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_SetLoad, temp>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
          temp=0
          do ig0= 1, fstrSOLID%TEMP_ngrp_tot
            ig= fstrSOLID%TEMP_ngrp_ID(ig0)
            val=fstrSOLID%TEMP_ngrp_val(ig0)
!C START & END
            iS0= hecMESH%node_group%grp_index(ig-1) + 1
            iE0= hecMESH%node_group%grp_index(ig  )
            do ik= iS0, iE0
              in   = hecMESH%node_group%grp_item(ik)
              temp( in ) = val
            enddo
          enddo
!C
!C +-------------------------------+
!C | ELEMENT-by-ELEMENT ASSEMBLING | 
!C | according to ELEMENT TYPE     |
!C +-------------------------------+
!C===
          do itype= 1, hecMESH%n_elem_type
            iS= hecMESH%elem_type_index(itype-1) + 1
            iE= hecMESH%elem_type_index(itype  )
            ic_type= hecMESH%elem_type_item(itype)
            if (.not. hecmw_is_etype_solid(ic_type)) cycle
!C** Set number of nodes		
            nn = hecmw_get_max_node(ic_type)
!C element loop
            do icel= iS, iE
!C** node ID
              iS= hecMESH%elem_node_index(icel-1)
              do j=1,nn
                nodLOCAL(j)=hecMESH%elem_node_item(iS+j)
!C** nodal coordinate
                xx(j)=hecMESH%node(3*nodLOCAL(j)-2)
                yy(j)=hecMESH%node(3*nodLOCAL(j)-1)
                zz(j)=hecMESH%node(3*nodLOCAL(j)  )
                tt(j)=temp( nodLOCAL(j) ) 
                tt0(j)=ref_temp
!C** create iwk array ***
                iwk(3*(j-1)+1)=3*(nodLOCAL(j)-1)+1
                iwk(3*(j-1)+2)=3*(nodLOCAL(j)-1)+2
                iwk(3*(j-1)+3)=3*(nodLOCAL(j)-1)+3
              enddo  
!C** section  ID
              isect= hecMESH%section_ID(icel)
!C** Get Properties
              call fstrNLGEOM_get_prop(hecMESH,isect,ee,pp,rho,alpha,thick)
!C** Create local stiffness		
              if (ic_type.eq.341) then
                call TLOAD_C3D4(xx,yy,zz,tt,tt0,ee,pp,alpha,vect)       
              else if ( ic_type.eq.351) then
                call TLOAD_C3D6(xx,yy,zz,tt,tt0,ee,pp,alpha,vect)
              else if ( ic_type.eq.361) then
                call TLOAD_C3D8(xx,yy,zz,tt,tt0,ee,pp,alpha,vect)
              else if ( ic_type.eq.342) then
                call TLOAD_C3D10(xx,yy,zz,tt,tt0,ee,pp,alpha,vect)  
              else if ( ic_type.eq.352) then
                call TLOAD_C3D15(xx,yy,zz,tt, tt0,ee,pp,alpha,vect)     
              else if ( ic_type.eq.362) then
                call TLOAD_C3D20(xx,yy,zz,tt,tt0,ee,pp,alpha,vect)
              endif
!C** Add vector
             do j=1,3*nn
               fstrSOLID%GL( iwk(j) )=fstrSOLID%GL( iwk(j) )+vect(j)
             enddo
            enddo
          enddo
          deallocate ( temp            ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <fstrNLGEOM_SetLoad, temp>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        endif
!C---------------------------------------------------------------------*
      else if( hecMESH%n_dof .eq. 2 ) then
!C
!C CLOAD
!C 
        do ig0= 1, fstrSOLID%CLOAD_ngrp_tot
          ig= fstrSOLID%CLOAD_ngrp_ID(ig0)
          ityp= fstrSOLID%CLOAD_ngrp_DOF(ig0)
          val= fstrSOLID%CLOAD_ngrp_val(ig0)
!C START & END
          iS0= hecMESH%node_group%grp_index(ig-1) + 1
          iE0= hecMESH%node_group%grp_index(ig  )
          do ik= iS0, iE0
            in   = hecMESH%node_group%grp_item(ik)
            fstrSOLID%GL(2*(in-1)+ityp)=fstrSOLID%GL(2*(in-1)+ityp)+val
          enddo
       enddo
!C
!C DLOAD
!C
       do ig0= 1, fstrSOLID%DLOAD_ngrp_tot
         ig= fstrSOLID%DLOAD_ngrp_ID(ig0)
         ltype= fstrSOLID%DLOAD_ngrp_LID(ig0)
         do i=0,6
           params(i)= fstrSOLID%DLOAD_ngrp_params(i,ig0)
         enddo
!C START & END
          fg_surf = (ltype == 100)
          if( fg_surf ) then                            ! surface group
            iS0= hecMESH%surf_group%grp_index(ig-1) + 1
            iE0= hecMESH%surf_group%grp_index(ig  )
          else                                          ! element group
            iS0= hecMESH%elem_group%grp_index(ig-1) + 1
            iE0= hecMESH%elem_group%grp_index(ig  )
          endif
!
          do ik= iS0, iE0
            if( fg_surf ) then                          ! surface group
              ltype  = hecMESH%surf_group%grp_item(2*ik)*10
              icel   = hecMESH%surf_group%grp_item(2*ik-1)
              ic_type= hecMESH%elem_type(icel)
            else                                        ! element group
              icel   = hecMESH%elem_group%grp_item(ik)
              ic_type= hecMESH%elem_type(icel)
            endif

!C** Create local stiffness		
           nn = hecmw_get_max_node(ic_type)
!C** node ID
           iS= hecMESH%elem_node_index(icel-1)
           do j=1,nn
             nodLOCAL(j)= hecMESH%elem_node_item (iS+j)
!C** nodal coordinate
             xx(j)=hecMESH%node(3*nodLOCAL(j)-2)
             yy(j)=hecMESH%node(3*nodLOCAL(j)-1)
!C** create iwk array ***
             iwk(2*(j-1)+1)=2*(nodLOCAL(j)-1)+1
             iwk(2*(j-1)+2)=2*(nodLOCAL(j)-1)+2
           enddo
!C** section  ID
           isect= hecMESH%section_ID(icel)
!C** Get Properties
           call fstrNLGEOM_get_prop(hecMESH,isect,ee,pp,rho,alpha,thick)
!C** Section Data
           id=hecMESH%section%sect_opt(isect)
           if( id.eq.0 ) then
             iset=1
           else if( id.eq.1) then
             iset=0
           else if( id.eq.2) then
             iset=2
           endif
           pa1=1.0
!C** Create local stiffness
           if (ic_type.eq.241) then
             call DL_C2D4(xx,yy,rho,pa1,ltype,params,vect,nsize,iset)
           else if ( ic_type.eq.242) then
             call DL_C2D8(xx,yy,rho,pa1,ltype,params,vect,nsize,iset)
           else if ( ic_type.eq.231) then
             call DL_C2D3(xx,yy,rho,pa1,ltype,params,vect,nsize,iset)
           else if ( ic_type.eq.232) then
             call DL_C2D6(xx,yy,rho,pa1,ltype,params,vect,nsize,iset)
           endif
!C** Add vector
           do j=1,nsize
             fstrSOLID%GL( iwk(j) )=fstrSOLID%GL( iwk(j) )+vect(j)
           enddo
         enddo
       enddo
!C
!C THERMAL LOAD USING TEMPERATURE
!C 
!C Set Temperature
!C 
       if( fstrSOLID%TEMP_ngrp_tot > 0 ) then
         allocate ( temp(hecMESH%n_node)             ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_SetLoad, temp>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
         temp=0
         do ig0= 1, fstrSOLID%TEMP_ngrp_tot
           ig= fstrSOLID%TEMP_ngrp_ID(ig0)
           val=fstrSOLID%TEMP_ngrp_val(ig0)
!C START & END
           iS0= hecMESH%node_group%grp_index(ig-1) + 1
           iE0= hecMESH%node_group%grp_index(ig  )
           do ik= iS0, iE0
             in   = hecMESH%node_group%grp_item(ik)
             temp( in ) = val
           enddo
         enddo
!C
!C +-------------------------------+
!C | ELEMENT-by-ELEMENT ASSEMBLING | 
!C | according to ELEMENT TYPE     |
!C +-------------------------------+
!C===
         do itype= 1, hecMESH%n_elem_type
           iS= hecMESH%elem_type_index(itype-1) + 1
           iE= hecMESH%elem_type_index(itype  )
           ic_type= hecMESH%elem_type_item(itype)
            if (.not. hecmw_is_etype_surface(ic_type)) cycle
!C** Set number of nodes
           nn = hecmw_get_max_node(ic_type)
!C element loop
           do icel= iS, iE
!C** node ID
             iS= hecMESH%elem_node_index(icel-1)
             do j=1,nn
               nodLOCAL(j)=hecMESH%elem_node_item (iS+j) 
!C** nodal coordinate
               xx(j)=hecMESH%node(3*nodLOCAL(j)-2)
               yy(j)=hecMESH%node(3*nodLOCAL(j)-1)
               tt(j)=temp( nodLOCAL(j) ) 
               tt0(j)=ref_temp
!C** create iwk array ***
               iwk(2*(j-1)+1)=2*(nodLOCAL(j)-1)+1
               iwk(2*(j-1)+2)=2*(nodLOCAL(j)-1)+2
             enddo  
!C** section  ID
             isect= hecMESH%section_ID(icel)
!C** Get Properties
             call fstrNLGEOM_get_prop(hecMESH,isect,ee,pp,rho,alpha,thick)
!C** Section Data
             id=hecMESH%section%sect_opt(isect)
             if( id.eq.0 ) then
               iset=1
             else if( id.eq.1) then
               iset=0
             else if( id.eq.2) then
               iset=2
             endif
             pa1=1.0
!C** Create local stiffness
             if (ic_type.eq.241) then
               call TLOAD_C2D4(xx,yy,tt,tt0,ee,pp,alpha,pa1,iset,vect)
             else if ( ic_type.eq.242) then
               call TLOAD_C2D8(xx,yy,tt,tt0,ee,pp,alpha,pa1,iset,vect)
             else if ( ic_type.eq.231) then
               call TLOAD_C2D3(xx,yy,tt,tt0,ee,pp,alpha,pa1,iset,vect)
             else if ( ic_type.eq.232) then
               call TLOAD_C2D6(xx,yy,tt,tt0,ee,pp,alpha,pa1,iset,vect)
             endif
!C** Add vector
             do j=1,2*nn
               fstrSOLID%GL( iwk(j) )=fstrSOLID%GL( iwk(j) )+vect(j)
             enddo
            enddo
          enddo
          deallocate ( temp              ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <fstrNLGEOM_SetLoad, temp>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        endif
      endif
!C
!C Update for fstrSOLID%GL
!C 
      if( hecMESH%n_dof .eq. 3 ) then
        call hecmw_update_3_R (hecMESH,fstrSOLID%GL,hecMESH%n_node)
      else if( hecMESH%n_dof .eq. 2 ) then
        call hecmw_update_2_R (hecMESH,fstrSOLID%GL,hecMESH%n_node)
      endif
!C
!C message
!C
      if( hecMESH%my_rank.eq.0) then
!!       write(IMSG,'("### fstrNLGEOM_SetLoad FINISHED!")')
      end if
      return
      end subroutine fstrNLGEOM_SetLoad
!C---------------------------------------------------------------------*
      subroutine fstrNLGEOM_InitialState(hecMESH,hecMAT,fstrSOLID)
!C---------------------------------------------------------------------*
!C*********************
!C  Initial State   
!C*********************
      use m_fstr
      implicit none
      type (hecmwST_local_mesh) :: hecMESH
      type (hecmwST_matrix)     :: hecMAT
      type (fstr_solid       )  :: fstrSOLID
!C* Local variables
      integer(kind=kint) :: neall,n_gauss,itype,iS,iE,ic_type,icel,ii 
      integer(kind=kint) :: ierror 
!C* Clear total_disp
      fstrSOLID%FACTOR    =0.0
      fstrSOLID%TOTAL_DISP=0.0
      fstrSOLID%GL=0.0d0
!C* 
      neall=0
      n_gauss=2
!C Allocate element stress
      if( hecMESH%n_dof .eq. 3 ) then
!C for 3D Problem
        do itype= 1, hecMESH%n_elem_type
          iS= hecMESH%elem_type_index(itype-1) + 1
          iE= hecMESH%elem_type_index(itype  )
          ic_type= hecMESH%elem_type_item(itype)
!C** Set number of nodes	
          if (ic_type.eq.342.or.ic_type.eq.352.or.ic_type.eq.362) then 
            n_gauss=3
          endif
!C element loop
          do icel= iS, iE
            neall=neall+1
          enddo
        enddo
        allocate(fstrSOLID%STRESS3D(6,n_gauss,n_gauss,n_gauss,neall)    ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_InitialState, STRESS3D>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        allocate(fstrSOLID%SURFF(8,n_gauss,n_gauss,n_gauss,neall)    ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_InitialState, SURFF>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
!       allocate(fstrSOLID%STRESS3D_nd(6,10,    1)    ,STAT=ierror)
        allocate(fstrSOLID%STRESS3D_nd(6,10,neall)    ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_InitialState, STRESS3D_nd>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
!       allocate(fstrSOLID%STRESS3D_ndP1(6,10,    1)    ,STAT=ierror)
        allocate(fstrSOLID%STRESS3D_ndP1(6,10,neall)    ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_InitialState, STRESS3D_ndP1>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        allocate(fstrSOLID%STRESS_nd(7*hecMESH%n_node)    ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_InitialState, STRESS_nd>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        allocate(fstrSOLID%STRESS_PRIS(12*hecMESH%n_node)    ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_InitialState, STRESS_PRIS>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
!       allocate(fstrSOLID%SURFF_nd(8,10,    1)    ,STAT=ierror)
        allocate(fstrSOLID%SURFF_nd(8,10,neall)    ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_InitialState, SURFF_nd>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        allocate(fstrSOLID%STRSURF_nd(8*hecMESH%n_node)    ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_InitialState, STRSURF_nd>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        allocate(fstrSOLID%STRAIN3D_P(6,n_gauss,n_gauss,n_gauss,neall)  ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_InitialState, STRAIN3D_P>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
!       allocate(fstrSOLID%STRAIN3D_P_nd(6,10,    1)  ,STAT=ierror)
        allocate(fstrSOLID%STRAIN3D_P_nd(6,10,neall)  ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_InitialState, STRAIN3D_P_nd>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        allocate(fstrSOLID%STRAIN_nd(6*hecMESH%n_node)  ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_InitialState, STRAIN_nd>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        allocate(fstrSOLID%STRAIN_PRIS(12*hecMESH%n_node)  ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_InitialState, STRAIN_PRIS>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        allocate(fstrSOLID%EP_BAR(n_gauss,n_gauss,n_gauss,neall)        ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_InitialState, EP_BAR>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        allocate(fstrSOLID%STRESS_YLD(n_gauss,n_gauss,n_gauss,neall)    ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_InitialState, STRESS_YLD>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        allocate(fstrSOLID%FLG_ELPL  (n_gauss,n_gauss,n_gauss,neall)    ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_InitialState, FLG_ELPL>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        allocate(fstrSOLID%mat_iso  (neall)    ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_InitialState, mat_iso>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        allocate(fstrSOLID%CURE3D_nod  (hecMESH%n_node)                     ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_InitialState, CURE3D_nod>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        allocate(fstrSOLID%CURE3D_nod_I  (hecMESH%n_node)                     ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_InitialState, CURE3D_nod_I>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        allocate(fstrSOLID%DELCURE3D_nod  (hecMESH%n_node)                     ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_InitialState, DELCURE3D_nod>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        allocate(fstrSOLID%TEMPI  (hecMESH%n_node)                     ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_InitialState, TEMPI>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        allocate(fstrSOLID%TEMPC  (hecMESH%n_node)                     ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_InitialState, TEMPC>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        allocate(fstrSOLID%DELTEMP  (hecMESH%n_node)                     ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_InitialState, DELTEMP>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        allocate(fstrSOLID%GBARU(3,3,hecMESH%n_node)  ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_InitialState, GBARU>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if

        if(fstrSOLID%idx_visco .ge. 4) then
          allocate(fstrSOLID%STRAIN3D_V1(6,n_gauss,n_gauss,n_gauss,neall)  ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_InitialState, STRAIN3D_V1>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
          fstrSOLID%STRAIN3D_V1=0.0
        end if
        if(fstrSOLID%idx_visco .eq. 6) then
          allocate(fstrSOLID%STRAIN3D_V2(6,n_gauss,n_gauss,n_gauss,neall)  ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_InitialState, STRAIN3D_V2>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
          fstrSOLID%STRAIN3D_V2=0.0
        end if

        fstrSOLID%STRESS3D       =0.0
        fstrSOLID%SURFF          =0.0
        fstrSOLID%STRESS3D_nd    =0.0
        fstrSOLID%STRESS3D_ndP1  =0.0
        fstrSOLID%STRESS_nd      =0.0
        fstrSOLID%STRESS_PRIS    =0.0
        fstrSOLID%SURFF_nd       =0.0
        fstrSOLID%STRSURF_nd     =0.0
        fstrSOLID%STRAIN3D_P   =0.0
        fstrSOLID%STRAIN3D_P_nd=0.0
        fstrSOLID%STRAIN_nd    =0.0
        fstrSOLID%STRAIN_PRIS  =0.0
        fstrSOLID%EP_BAR=0.0
        fstrSOLID%STRESS_YLD=fstrSOLID%sig_y0
        fstrSOLID%FLG_ELPL  =0
        fstrSOLID%mat_iso  =0
        fstrSOLID%CURE3D_nod   =0.0
        fstrSOLID%CURE3D_nod_I =0.0
        fstrSOLID%DELCURE3D_nod=0.0
        fstrSOLID%TEMPI  =0.0
        fstrSOLID%TEMPC  =0.0
        fstrSOLID%DELTEMP=0.0
        fstrSOLID%GBARU    =0.0
!!-- 20100330 ----!
        call fstrNLGEOM_init_gbar( hecMESH,hecMAT,fstrSOLID)
!!-- 20100330 ----!

      else if( hecMESH%n_dof .eq. 2) then
!C for 2D Problem
        do itype= 1, hecMESH%n_elem_type
          iS= hecMESH%elem_type_index(itype-1) + 1
          iE= hecMESH%elem_type_index(itype  )
          ic_type= hecMESH%elem_type_item(itype)
!C** Set number of nodes	
          if ( ic_type.eq.242 .or. ic_type.eq.232 ) then
            n_gauss=3
          endif
!C element loop
          do icel= iS, iE
            neall=neall+1
          enddo
        enddo
        allocate(fstrSOLID%STRESS2D(4,n_gauss,n_gauss,neall)            ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_InitialState, STRESS2D>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        allocate(fstrSOLID%STRAIN3D_P(4,n_gauss,n_gauss,1,neall)  ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_InitialState, STRAIN3D_P>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        allocate(fstrSOLID%EP_BAR(n_gauss,n_gauss,1,neall)        ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_InitialState, EP_BAR>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        allocate(fstrSOLID%STRESS_YLD(n_gauss,n_gauss,1,neall)    ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_InitialState, STRESS_YLD>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        allocate(fstrSOLID%FLG_ELPL  (n_gauss,n_gauss,1,neall)    ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_InitialState, FLG_ELPL>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        fstrSOLID%STRESS2D  =0.0
        fstrSOLID%STRAIN3D_P=0.0
        fstrSOLID%EP_BAR=0.0
        fstrSOLID%STRESS_YLD=fstrSOLID%sig_y0
        fstrSOLID%FLG_ELPL  =0
      endif
!
      if( hecMESH%my_rank.eq.0) then
       write(IMSG,*) '####fstrNLGEOM_InitialState finished'
      end if
      return
      end subroutine fstrNLGEOM_InitialState
!----------------------------------------------------------------------*
      subroutine fstrNLGEOM_SetKQ(hecMESH,hecMAT,fstrSOLID,DT0,iter,i_step)
!----------------------------------------------------------------------*
!C*******************************************
!C  Set internal force & tangential stiffness
!C*******************************************
      use m_fstr
      use m_static_mat_ass_main
      implicit none
      type (hecmwST_local_mesh) :: hecMESH
      type (hecmwST_matrix    ) :: hecMAT
      type (fstr_solid       )  :: fstrSOLID
!C
!C Local variables for 2D and 3D elements
!C
      integer(kind=kint) :: ndof, count, itype, iS, iE, ic_type, nn, icel
      integer(kind=kint) :: iiS, i, j, k, l, idx, idx_e, isect, id, iset
      integer(kind=kint) :: inod, inod_e
      integer(kind=kint) :: iter, i_step 
      real(kind=kreal) :: xx(20), yy(20), zz(20), ee, pp, local_stf(1830), stiffness(20*6, 20*6)
      real(kind=kreal) :: disp(60), qf(60), qf_res(60), tt(20), tti(20), CUR_nod(20), CUR_nod_I(20)
      integer(kind=kint) :: nodLOCAL(20)
      integer(kind=kint) :: g1, g2, g3, mid 

      real(kind=kreal) :: rho, alpha, thick, param1 
      real(kind=kreal) :: DT0, STRAIN_V1(6,3,3,3),STRAIN_V2(6,3,3,3)

      real(kind=kreal) :: BIJ(3,3),GU(3,3),GUB(3,3),ALPHA_B(3),MAT_OTHO1(9) 
!C
!C +-------+
!C | INIT. |
!C +-------+
!C====
      hecMAT%AL= 0.d0
      hecMAT%AU= 0.d0
      hecMAT%X = 0.d0
      hecMAT%D = 0.d0
!C*
!C*ZERO CLEAR ARRAY QFORCE
!C*
!!!!  20111026  -----
      fstrSOLID%QFORCE_RES=0.0d0
      do i=1, hecMESH%n_node*  hecMESH%n_dof
        fstrSOLID%QFORCE_RES(i)=0.0d0
      enddo
!
      fstrSOLID%QFORCE    =0.0d0
!
      STRAIN_V1=0.0d0
      STRAIN_V2=0.0d0
!C
!C +-------------------------------+
!C | ELEMENT-by-ELEMENT ASSEMBLING | 
!C | according to ELEMENT TYPE     |
!C +-------------------------------+
      ndof = hecMAT%NDOF

      count=0
      do itype= 1, hecMESH%n_elem_type
        iS= hecMESH%elem_type_index(itype-1) + 1
        iE= hecMESH%elem_type_index(itype  )
        ic_type= hecMESH%elem_type_item(itype)
!C** Ignore link elements
        if (hecmw_is_etype_link(ic_type)) cycle
!C** Set number of nodes	
        nn = hecmw_get_max_node(ic_type)
!C element loop
        do icel= iS, iE
!C** node ID
          iiS= hecMESH%elem_node_index(icel-1)
          do j=1,nn
            nodLOCAL(j)= hecMESH%elem_node_item (iiS+j)
!C** nodal coordinate
            xx(j)=hecMESH%node(3*nodLOCAL(j)-2)
            yy(j)=hecMESH%node(3*nodLOCAL(j)-1)
            zz(j)=hecMESH%node(3*nodLOCAL(j)  )
            tt (j)=fstrSOLID%TEMPC(nodLOCAL(j) )
            tti(j)=fstrSOLID%TEMPI(nodLOCAL(j) )
            CUR_nod  (j)=fstrSOLID%CURE3D_nod  (nodLOCAL(j) )
            CUR_nod_I(j)=fstrSOLID%CURE3D_nod_I(nodLOCAL(j) )
            idx_e = ndof*(j-1)
            idx = ndof*(nodLOCAL(j)-1)
            do i = 1, ndof
              disp(idx_e+i)=fstrSOLID%TOTAL_DISP(idx+i)
            enddo
          enddo  
!C** section  ID
          isect= hecMESH%section_ID(icel)
!
          mid=hecMESH%section%sect_mat_ID_item(isect) 
!C** Get Properties
          call fstrNLGEOM_get_prop(hecMESH,isect,ee,pp,rho,alpha,thick)
!C** Section Data
          if (ndof .eq. 2) then
            id=hecMESH%section%sect_opt(isect)
            if( id.eq.0 ) then
              iset=1
            else if( id.eq.1) then
              iset=0
            else if( id.eq.2) then
              iset=2
            endif
            param1=1.0
          endif
!C** Create local stiffness
          count=count+1
!!-- 20100330 ----!
         do i = 1, 3
           ALPHA_B(i) = 0.d0
           do j = 1, 3
             GU (i,j) = 0.d0
             GUB(i,j) = 0.d0
             BIJ(i,j) = 0.d0
           end do
         end do
         do i = 1, 3
           GU (i,i) = 1.d0
         end do
         do i = 1, 9
           MAT_OTHO1(i) = 7.d0
         end do
!
!!----- 20120927 ----!
       if ( fstrSOLID%mat_iso_num(1) .ne. -1) then
         do i = 1, 4
          if ( mid .eq. i) then
            do l=1,nn
             do j = 1, 3
              do k = 1, 3
               fstrSOLID%GBARU(j,k,nodLOCAL(l)) =  fstrSOLID%GBARU_matrl(j,k,i)
              end do
             end do
            end do
          end if
         end do
       else if ( fstrSOLID%mat_iso_num(1) .eq. -1) then
            do l=1,nn
             do j = 1, 3
              do k = 1, 3
               fstrSOLID%GBARU(j,k,nodLOCAL(l)) =  fstrSOLID%GBARU_ALL(j,k,hecMESH%global_node_ID(nodLOCAL(l)))
              end do
             end do
            end do
       end if

         do l=1,nn
          do j = 1, 3
           do k = 1, 3
            GUB(j,k) = GUB(j,k) + fstrSOLID%GBARU(j,k,nodLOCAL(l))/nn 
           end do
          end do
         end do
!
          do i = 1, 3
           do j = 1, 3
            do k = 1, 3
             BIJ(i,j) = BIJ(i,j) + GU (i,k)*GUB(j,k)
            end do
           end do
          end do
!
         if( fstrSOLID%mat_iso(count) .eq. 1 ) then
           do i = 1, 3
             ALPHA_B(i)   = fstrSOLID%mat_beta(i,1)
           end do
           do i = 1, 9
             MAT_OTHO1(i) = fstrSOLID%mat_otho (i,1)
           end do
         end if
!!-- 20100330 ----!

          if ( ic_type.eq.241) then
            call STF_C2D4TL(ILOG,xx,yy,disp,ee,pp,param1,                     &
                          iset,fstrSOLID%STRESS2D(:,:,:,count),local_stf,qf,  &
                          fstrSOLID%idx_elpl,fstrSOLID%h_dash,                &
                          fstrSOLID%FLG_ELPL(:,:,:,count),                    &
                          fstrSOLID%STRESS_YLD(:,:,:,count))
          else if ( ic_type.eq.242) then
            call STF_C2D8TL(ILOG,xx,yy,disp,ee,pp,param1,                     &
                          iset,fstrSOLID%STRESS2D(:,:,:,count), local_stf,qf, &
                          fstrSOLID%idx_elpl,fstrSOLID%h_dash,                &
                          fstrSOLID%FLG_ELPL(:,:,:,count),                    &
                          fstrSOLID%STRESS_YLD(:,:,:,count))
          else if ( ic_type.eq.231) then
            call STF_C2D3TL(ILOG,xx,yy,disp,ee,pp,param1,                     &
                          iset,fstrSOLID%STRESS2D(:,:,:,count),local_stf,qf,  &
                          fstrSOLID%idx_elpl,fstrSOLID%h_dash,                &
                          fstrSOLID%FLG_ELPL(:,:,:,count),                    &
                          fstrSOLID%STRESS_YLD(:,:,:,count))
          else if ( ic_type.eq.232) then
            call STF_C2D6TL(ILOG,xx,yy,disp,ee,pp,param1,                     &
                          iset,fstrSOLID%STRESS2D(:,:,:,count),local_stf,qf,  &
                          fstrSOLID%idx_elpl,fstrSOLID%h_dash,                &
                          fstrSOLID%FLG_ELPL(:,:,:,count),                    &
                          fstrSOLID%STRESS_YLD(:,:,:,count))
!
          else if (ic_type.eq.341) then  
            call STF_C3D4TL(ILOG,xx,yy,zz,disp,ee,pp,                      &
                          fstrSOLID%STRESS3D(:,:,:,:,count),local_stf,qf,  &
                          fstrSOLID%idx_elpl,fstrSOLID%h_dash,             &
                          fstrSOLID%FLG_ELPL(:,:,:,count),fstrSOLID%STRESS_YLD(:,:,:,count))
          else if ( ic_type.eq.351) then
            call STF_C3D6TL(ILOG,xx,yy,zz,disp,ee,pp,                      &
                          fstrSOLID%STRESS3D(:,:,:,:,count),local_stf,qf,  &
                          fstrSOLID%idx_elpl,fstrSOLID%h_dash,             &
                          fstrSOLID%FLG_ELPL(:,:,:,count),fstrSOLID%STRESS_YLD(:,:,:,count))
          else if ( ic_type.eq.361) then
            call STF_C3D8TL(ILOG,xx,yy,zz,disp,ee,pp,                      &
                          fstrSOLID%STRESS3D(:,:,:,:,count),local_stf,qf,  &
                          fstrSOLID%idx_elpl,fstrSOLID%h_dash,             &
                          fstrSOLID%FLG_ELPL(:,:,:,count),fstrSOLID%STRESS_YLD(:,:,:,count))
          else if ( ic_type.eq.342) then
!!!!!  visco  ---------
            if(fstrSOLID%idx_visco .ge. 4) then
              do g1= 1, 3
               do g2= 1, 3
                do g3= 1, 3
                 do  i= 1, 6
                   STRAIN_V1(i,g1,g2,g3) = fstrSOLID%STRAIN3D_V1(i,g1,g2,g3,count)
                 end do
                end do
               end do
              end do
            end if
            if(fstrSOLID%idx_visco .eq. 6) then
              do g1= 1, 3
               do g2= 1, 3
                do g3= 1, 3
                 do  i= 1, 6
                   STRAIN_V2(i,g1,g2,g3) = fstrSOLID%STRAIN3D_V2(i,g1,g2,g3,count)
                 end do
                end do
               end do
              end do
            end if
!!!!!
            call STF_C3D10TL(ILOG,xx,yy,zz,disp,ee,pp,                                          &
                          fstrSOLID%STRESS3D(:,:,:,:,count),local_stf,qf,qf_res,                &
                          fstrSOLID%idx_elpl,fstrSOLID%h_dash,                                  &
                          fstrSOLID%FLG_ELPL(:,:,:,count),fstrSOLID%STRESS_YLD(:,:,:,count),    &
                          fstrSOLID%idx_temp,                                                   &
                          tt,tti,CUR_nod, CUR_nod_I,                                            &
                          fstrSOLID%ccc_rl,fstrSOLID%ramd_chm,fstrSOLID%ccc_rl1,                &
                          fstrSOLID%ccc_rl2,fstrSOLID%yeta,fstrSOLID%yeta1,fstrSOLID%yeta2,     &
                          fstrSOLID%idx_visco,DT0,STRAIN_V1,STRAIN_V2,iter,i_step,count,icel,   &
                          mid,alpha,fstrSOLID%mat_iso(count),BIJ,ALPHA_B,MAT_OTHO1)
!!!!!
          else if ( ic_type.eq.352) then
            call STF_C3D15TL(ILOG,xx,yy,zz,disp,ee,pp,                     &
                          fstrSOLID%STRESS3D(:,:,:,:,count),local_stf,qf,  &
                          fstrSOLID%idx_elpl,fstrSOLID%h_dash,             &
                          fstrSOLID%FLG_ELPL(:,:,:,count),fstrSOLID%STRESS_YLD(:,:,:,count))
          else if ( ic_type.eq.362) then
            call STF_C3D20TL(ILOG,xx,yy,zz,disp,ee,pp,                     &
                          fstrSOLID%STRESS3D(:,:,:,:,count),local_stf,qf,  &
                          fstrSOLID%idx_elpl,fstrSOLID%h_dash,             &
                          fstrSOLID%FLG_ELPL(:,:,:,count),fstrSOLID%STRESS_YLD(:,:,:,count))
          else
            write(*,*) '###ERROR### : Element type not supported for nonlinear static analysis'
            write(*,*) ' ic_type = ', ic_type
            call hecmw_abort(hecmw_comm_get_comm())
          endif
!C** restore local stiffness matrix ***
          call fstr_local_stf_restore(local_stf, nn*ndof, stiffness)
!C== CONSTRUCT the GLOBAL MATRIX STARTED
          call hecmw_mat_ass_elem(hecMAT, nn, nodLOCAL, stiffness)
!C* add local to global ***
          do inod_e= 1, nn
            inod = nodLOCAL(inod_e)
            idx = ndof*(inod-1)
            idx_e = ndof*(inod_e-1)
            do i = 1, ndof
!!!   20111021  ----
              fstrSOLID%QFORCE    (idx+i) = fstrSOLID%QFORCE    (idx+i)+qf    (idx_e+i)
              fstrSOLID%QFORCE_RES(idx+i) = fstrSOLID%QFORCE_RES(idx+i)+qf_res(idx_e+i)
            enddo
          enddo
        enddo
      enddo

!* for EQUATION
      call hecmw_mat_ass_equation ( hecMESH, hecMAT )

!C
!C Update for fstrSOLID%QFORCE
!C 
      if( ndof .eq. 3 ) then
        call hecmw_update_3_R(hecMESH,fstrSOLID%QFORCE    ,hecMESH%n_node)
        call hecmw_update_3_R(hecMESH,fstrSOLID%QFORCE_RES,hecMESH%n_node)
      else if( ndof .eq. 2 ) then
        call hecmw_update_2_R(hecMESH,fstrSOLID%QFORCE    ,hecMESH%n_node)
        call hecmw_update_2_R(hecMESH,fstrSOLID%QFORCE_RES,hecMESH%n_node)
      endif
!C
!C Message
!C
      if( hecMESH%my_rank.eq.0) then
!!       write(IMSG,*) '####fstrNLGEOM_SetKQ finished'
      end if
      return
      end subroutine fstrNLGEOM_SetKQ
!----------------------------------------------------------------------*
      subroutine fstrNLGEOM_AddBC(hecMESH,hecMAT,fstrSOLID,iter)
!----------------------------------------------------------------------*
!C*******************************************
!C  Add Essential Boundary Conditions
!C*******************************************
      use m_fstr
      implicit none

      type (hecmwST_local_mesh) :: hecMESH
      type (hecmwST_matrix)     :: hecMAT
      type (fstr_solid        ) :: fstrSOLID
      integer(kind=kint) :: iter

!C Local variables
      integer(kind=kint) :: ig0, ig, ityp, idofS, idofE, idof, iS0, iE0, ik, in
      real(kind=kreal) :: RHS 

!!
      do ig0 = 1, fstrSOLID%BOUNDARY_ngrp_tot
        ig   = fstrSOLID%BOUNDARY_ngrp_ID(ig0)
        RHS  = fstrSOLID%BOUNDARY_ngrp_val(ig0)

!ogasawara till next ogasawara
        RHS= RHS*(fstrSOLID%FACTOR(2)-fstrSOLID%FACTOR(1))

        if( iter .gt. 1 ) then
           RHS= 0.0
        end if
!ogasawara

        ityp = fstrSOLID%BOUNDARY_ngrp_type(ig0)

        idofS = ityp/10
        idofE = ityp - idofS*10

        iS0 = hecMESH%node_group%grp_index(ig-1) + 1
        iE0 = hecMESH%node_group%grp_index(ig  )

        do ik = iS0, iE0
          in = hecMESH%node_group%grp_item(ik)

          do idof = idofS, idofE
            call hecmw_mat_ass_bc(hecMAT, in, idof, RHS)
          enddo
        enddo
      enddo
!C
!C Message
!C
      if( hecMESH%my_rank.eq.0) then
!!       write(IMSG,*) '####fstrNLGEOM_AddBC finished'
      end if
      return
      end subroutine fstrNLGEOM_AddBC
!C---------------------------------------------------------------------*
      subroutine fstrNLGEOM_SetResidual(hecMESH,hecMAT,fstrSOLID,factor,res )
!C---------------------------------------------------------------------*
!C*******************************************
!C  Set resudual 
!C*******************************************
      use m_fstr
      implicit none

      type (hecmwST_local_mesh) :: hecMESH
      type (hecmwST_matrix)     :: hecMAT
      type (hecmwST_matrix)     :: MAT2
      type (fstr_solid        ) :: fstrSOLID

      real(kind=kreal) factor,res
!C Local variables
      integer(kind=kint) i,ndof,ig0,ig,ityp,iS0,iE0,ik,in,idof1,idof2,idof 
      real(kind=kreal) RHS 
      integer(kind=kint) ierror 

!C**
!C** Set residual load
!C**
      hecMAT%NDOF = hecMESH%n_dof
      ndof = hecMAT%NDOF
!
      allocate (MAT2%B(ndof*hecMAT%NP)          ,STAT=ierror )
          if( ierror /= 0 ) then
            write(idbg,*) 'stop due to allocation error <fstrNLGEOM_SetResidual, B>'
            write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
            call flush(idbg)
            call hecmw_abort( hecmw_comm_get_comm())
          end if
!!
!!!   20111021  ----
      do i=1, hecMESH%n_node*  hecMESH%n_dof
        hecMAT%B(i)=fstrSOLID%GL(i)-fstrSOLID%QFORCE(i)
        MAT2  %B(i)=fstrSOLID%GL(i)-fstrSOLID%QFORCE_RES(i)
      enddo
!C**
!C** Consider SPC condition
!C**
      ndof = hecMAT%NDOF
      do ig0= 1, fstrSOLID%BOUNDARY_ngrp_tot
        ig= fstrSOLID%BOUNDARY_ngrp_ID(ig0)
        RHS= fstrSOLID%BOUNDARY_ngrp_val(ig0)
        ityp= fstrSOLID%BOUNDARY_ngrp_type(ig0)
        iS0= hecMESH%node_group%grp_index(ig-1) + 1
        iE0= hecMESH%node_group%grp_index(ig  )
        do ik= iS0, iE0
          in   = hecMESH%node_group%grp_item(ik)
          idof1 = ityp/10
          idof2 = ityp - idof1*10
          do idof=idof1,idof2
            hecMAT%B( ndof*(in-1) + idof ) = factor*RHS-fstrSOLID%TOTAL_DISP(ndof*(in-1) + idof)
            MAT2  %B( ndof*(in-1) + idof ) = factor*RHS-fstrSOLID%TOTAL_DISP(ndof*(in-1) + idof)
          enddo
        enddo
      enddo
!C**
!C** Calculate residual
!C**
      res=0.0
!
      if(fstrSOLID%idx_temp .eq. 0 ) then
       if(fstrSOLID%idx_visco .lt. 2 ) then
         do i=1, hecMESH%nn_internal *  hecMESH%n_dof
           res=res+ ( MAT2  %B(i) )**2
         enddo
       else if(fstrSOLID%idx_visco .ge. 2 ) then
         do i=1, hecMESH%nn_internal *  hecMESH%n_dof
           res=res+ ( MAT2  %B(i) )**2
         enddo
       end if
      else if(fstrSOLID%idx_temp .eq. 1 ) then
         do i=1, hecMESH%nn_internal *  hecMESH%n_dof
           res=res+ ( MAT2  %B(i) )**2
         enddo
      end if
!!!!!
      if( hecMESH%my_rank.eq.0) then
!!       write(IMSG,*) '####fstrNLGEOM_SetResidual finished'
      end if
!
      deallocate(MAT2%B                   ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <fstrNLGEOM_SetResidual, B>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      return
      end subroutine fstrNLGEOM_SetResidual
!C---------------------------------------------------------------------*
      subroutine fstrNLGEOM_Update(hecMESH,hecMAT,fstrSOLID,factor,iter,i_step,DT0)
!C---------------------------------------------------------------------*
!C*******************************************
!C  Update Stress Strain Displacement
!C*******************************************
      use m_fstr
      implicit none

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

      real(kind=kreal) :: factor
      integer(kind=kint) :: iter, i_step
!C* Local variable
      integer(kind=kint) :: ndof, count, itype, iS, iE, ic_type, nn, icel
      integer(kind=kint) :: i, j, k, l, isect, id, iset
      integer(kind=kint) :: nodLocal(20), ID_area 
      integer(kind=kint) :: g1, g2, g3, mid 
      real(kind=kreal) :: xx(20), yy(20), zz(20), tt(20), tti(20), ee, pp
      real(kind=kreal) :: disp(60), ddisp(60), CUR_nod(20), CUR_nod_I(20)
      real(kind=kreal) :: rho, alpha, thick, param1 
      real(kind=kreal) :: DT0, STRAIN_V1(6,3,3,3),STRAIN_V2(6,3,3,3)

      real(kind=kreal) :: BIJ(3,3),GU(3,3),GUB(3,3),ALPHA_B(3),MAT_OTHO1(9) 

!C +-------------------------------+
!C | ELEMENT-by-ELEMENT ASSEMBLING |
!C | according to ELEMENT TYPE     |
!C +-------------------------------+
      ndof = hecMAT%NDOF
!
      STRAIN_V1=0.0d0
      STRAIN_V2=0.0d0
!
      count=0
      do itype= 1, hecMESH%n_elem_type
        iS= hecMESH%elem_type_index(itype-1) + 1
        iE= hecMESH%elem_type_index(itype  )
        ic_type= hecMESH%elem_type_item(itype)
        if (hecmw_is_etype_link(ic_type)) cycle
!C** Set number of nodes
        nn = hecmw_get_max_node(ic_type)
!C element loop
        do icel= iS, iE
!C** node ID
          iS= hecMESH%elem_node_index(icel-1)
          do j=1,nn
            nodLOCAL(j)= hecMESH%elem_node_item (iS+j)
!C** nodal coordinate
            xx (j)=hecMESH%node(3*nodLOCAL(j)-2)
            yy (j)=hecMESH%node(3*nodLOCAL(j)-1)
            zz (j)=hecMESH%node(3*nodLOCAL(j)  )
            tt (j)=fstrSOLID%TEMPC(nodLOCAL(j) )
            tti(j)=fstrSOLID%TEMPI(nodLOCAL(j) )
            CUR_nod  (j)=fstrSOLID%CURE3D_nod  (nodLOCAL(j) )
            CUR_nod_I(j)=fstrSOLID%CURE3D_nod_I(nodLOCAL(j) )
!C** set disp & ddisp
            do i = 1, ndof
              disp(ndof*(j-1)+i)=fstrSOLID%TOTAL_DISP(ndof*(nodLOCAL(j)-1)+i)
              ddisp(ndof*(j-1)+i)=hecMAT%X(ndof*(nodLOCAL(j)-1)+i)
            enddo
          enddo  
!C** section  ID
          isect= hecMESH%section_ID(icel)
          mid  = hecMESH%section%sect_mat_ID_item(isect) 
!C** Get Properties
          call fstrNLGEOM_get_prop(hecMESH,isect,ee,pp,rho,alpha,thick)
!C** Section Data
          if (ndof .eq. 2) then
            id=hecMESH%section%sect_opt(isect)
            if( id.eq.0 ) then
              iset=1
            else if( id.eq.1) then
              iset=0
            else if( id.eq.2) then
              iset=2
            endif
            param1=1.0
          endif
!C** Update Stress
          count=count+1
!!-- 20100330 ----!
         do i = 1, 3
           ALPHA_B(i) = 0.d0
           do j = 1, 3
             GU (i,j) = 0.d0
             GUB(i,j) = 0.d0
             BIJ(i,j) = 0.d0
           end do
         end do
         do i = 1, 3
           GU (i,i) = 1.d0
         end do
         do i = 1, 9
           MAT_OTHO1(i) = 7.d0
         end do
!
!!----- 20120927 ----!
       if ( fstrSOLID%mat_iso_num(1) .ne. -1) then
         do i = 1, 4
          if ( mid .eq. i) then
            do l=1,nn
             do j = 1, 3
              do k = 1, 3
               fstrSOLID%GBARU(j,k,nodLOCAL(l)) =  fstrSOLID%GBARU_matrl(j,k,i)
              end do
             end do
            end do
          end if
         end do
       else if ( fstrSOLID%mat_iso_num(1) .eq. -1) then
            do l=1,nn
             do j = 1, 3
              do k = 1, 3
               fstrSOLID%GBARU(j,k,nodLOCAL(l)) =  fstrSOLID%GBARU_ALL(j,k,hecMESH%global_node_ID(nodLOCAL(l)))
              end do
             end do
            end do
       end if

         do l=1,nn
          do j = 1, 3
           do k = 1, 3
            GUB(j,k) = GUB(j,k) + fstrSOLID%GBARU(j,k,nodLOCAL(l))/nn 
           end do
          end do
         end do
!
          do i = 1, 3
           do j = 1, 3
            do k = 1, 3
             BIJ(i,j) = BIJ(i,j) + GU (i,k)*GUB(j,k)
            end do
           end do
          end do
!
         if( fstrSOLID%mat_iso(count) .eq. 1 ) then
           do i = 1, 3
             ALPHA_B(i)   = fstrSOLID%mat_beta(i,1)
           end do
           do i = 1, 9
             MAT_OTHO1(i) = fstrSOLID%mat_otho (i,1)
           end do
         end if

!!-- 20100330 ----!
          if ( ic_type.eq.241) then
            call UPDATE_C2D4TL(ILOG,xx,yy,disp,ddisp,                             &
                               ee,pp,param1,iset,fstrSOLID%STRESS2D(:,:,:,count), &
                               fstrSOLID%idx_elpl,fstrSOLID%h_dash,               &
                               fstrSOLID%FLG_ELPL(:,:,:,count),                   &
                               fstrSOLID%STRESS_YLD(:,:,:,count),                 &
                               fstrSOLID%STRAIN3D_P(:,:,:,:,count),               &
                               fstrSOLID%EP_BAR(:,:,:,count),iter)
          else if ( ic_type.eq.242) then
            call UPDATE_C2D8TL(ILOG,xx,yy,disp,ddisp,                             &
                               ee,pp,param1,iset,fstrSOLID%STRESS2D(:,:,:,count), &
                               fstrSOLID%idx_elpl,fstrSOLID%h_dash,               &
                               fstrSOLID%FLG_ELPL(:,:,:,count),                   &
                               fstrSOLID%STRESS_YLD(:,:,:,count),                 &
                               fstrSOLID%STRAIN3D_P(:,:,:,:,count),               &
                               fstrSOLID%EP_BAR(:,:,:,count),iter)
          else if ( ic_type.eq.231) then
            call UPDATE_C2D3TL(ILOG,xx,yy,disp,ddisp,                             &
                               ee,pp,param1,iset,fstrSOLID%STRESS2D(:,:,:,count), &
                               fstrSOLID%idx_elpl,fstrSOLID%h_dash,               &
                               fstrSOLID%FLG_ELPL(:,:,:,count),                   &
                               fstrSOLID%STRESS_YLD(:,:,:,count),                 &
                               fstrSOLID%STRAIN3D_P(:,:,:,:,count),               &
                               fstrSOLID%EP_BAR(:,:,:,count),iter)
          else if ( ic_type.eq.232) then
            call UPDATE_C2D6TL(ILOG,xx,yy,disp,ddisp,                             &
                               ee,pp,param1,iset,fstrSOLID%STRESS2D(:,:,:,count), &
                               fstrSOLID%idx_elpl,fstrSOLID%h_dash,               &
                               fstrSOLID%FLG_ELPL(:,:,:,count),                   &
                               fstrSOLID%STRESS_YLD(:,:,:,count),                 &
                               fstrSOLID%STRAIN3D_P(:,:,:,:,count),               &
                               fstrSOLID%EP_BAR(:,:,:,count),iter)
!
          else if (ic_type.eq.341) then
            call UPDATE_C3D4TL(ILOG,xx,yy,zz,disp,ddisp,                &
                               ee,pp,fstrSOLID%STRESS3D(:,:,:,:,count), &
                               fstrSOLID%idx_elpl,fstrSOLID%h_dash,     &
                               fstrSOLID%FLG_ELPL(:,:,:,count),         &
                               fstrSOLID%STRESS_YLD(:,:,:,count),       &
                               fstrSOLID%STRAIN3D_P(:,:,:,:,count),     &
                               fstrSOLID%EP_BAR(:,:,:,count),iter)
          else if ( ic_type.eq.351) then
            call UPDATE_C3D6TL(ILOG,xx,yy,zz,disp,ddisp,                &
                               ee,pp,fstrSOLID%STRESS3D(:,:,:,:,count), &
                               fstrSOLID%idx_elpl,fstrSOLID%h_dash,     &
                               fstrSOLID%FLG_ELPL(:,:,:,count),         &
                               fstrSOLID%STRESS_YLD(:,:,:,count),       &
                               fstrSOLID%STRAIN3D_P(:,:,:,:,count),     &
                               fstrSOLID%EP_BAR(:,:,:,count),iter)
          else if ( ic_type.eq.361) then
            call UPDATE_C3D8TL(ILOG,xx,yy,zz,disp,ddisp,                &
                               ee,pp,fstrSOLID%STRESS3D(:,:,:,:,count), &
                               fstrSOLID%idx_elpl,fstrSOLID%h_dash,     &
                               fstrSOLID%FLG_ELPL(:,:,:,count),         &
                               fstrSOLID%STRESS_YLD(:,:,:,count),       &
                               fstrSOLID%STRAIN3D_P(:,:,:,:,count),     &
                               fstrSOLID%EP_BAR(:,:,:,count),iter)
          else if ( ic_type.eq.342) then
!!!!!  visco  ---------
            if(fstrSOLID%idx_visco .ge. 4) then
              do g1= 1, 3
               do g2= 1, 3
                do g3= 1, 3
                 do  i= 1, 6
                   STRAIN_V1(i,g1,g2,g3) = fstrSOLID%STRAIN3D_V1(i,g1,g2,g3,count)
                 end do
                end do
               end do
              end do
            end if
            if(fstrSOLID%idx_visco .eq. 6) then
              do g1= 1, 3
               do g2= 1, 3
                do g3= 1, 3
                 do  i= 1, 6
                   STRAIN_V2(i,g1,g2,g3) = fstrSOLID%STRAIN3D_V2(i,g1,g2,g3,count)
                 end do
                end do
               end do
              end do
            end if
!!!!!
            call UPDATE_C3D10TL(ILOG,xx,yy,zz,disp,ddisp,                           &
                               ee,pp,fstrSOLID%STRESS3D(:,:,:,:,count),             &
                               fstrSOLID%idx_elpl,fstrSOLID%h_dash,                 &
                               fstrSOLID%FLG_ELPL(:,:,:,count),                     &
                               fstrSOLID%STRESS_YLD(:,:,:,count),                   &
                               fstrSOLID%STRAIN3D_P(:,:,:,:,count),                 &
                               fstrSOLID%EP_BAR(:,:,:,count),iter,i_step,           &
                               fstrSOLID%idx_temp,                                  &
                               tt,tti,CUR_nod, CUR_nod_I,                           &
                               fstrSOLID%ccc_rl,fstrSOLID%ramd_chm,                 &
                               fstrSOLID%ccc_rl1,fstrSOLID%ccc_rl2,fstrSOLID%yeta,  &
                               fstrSOLID%yeta1,fstrSOLID%yeta2,fstrSOLID%idx_visco, &
                               DT0,STRAIN_V1,STRAIN_V2,count,icel,mid,alpha,        &
                               fstrSOLID%STRESS3D_nd(:,:,count),                    &
                               fstrSOLID%STRAIN3D_P_nd(:,:,count),                  & 
                               fstrSOLID%STRESS3D_ndP1(:,:,count),                  &  
                               fstrSOLID%SURFF_nd(:,:,count),                       &  
                               fstrSOLID%idx_stres,                                 & 
                               fstrSOLID%mat_iso(count),BIJ,ALPHA_B,MAT_OTHO1)
!                              fstrSOLID%STRESS3D_nd(:,:,1),                        &
!                              fstrSOLID%STRAIN3D_P_nd(:,:,1),                      & 
!                              fstrSOLID%STRESS3D_ndP1(:,:,1),                      &  
!                              fstrSOLID%SURFF_nd(:,:,1) ) 
!!!!!  visco  ---------
            if(fstrSOLID%idx_visco .ge. 4) then
              do g1= 1, 3
               do g2= 1, 3
                do g3= 1, 3
                 do  i= 1, 6
                   fstrSOLID%STRAIN3D_V1(i,g1,g2,g3,count) = STRAIN_V1(i,g1,g2,g3)
                 end do
                end do
               end do
              end do
            end if
            if(fstrSOLID%idx_visco .eq. 6) then
              do g1= 1, 3
               do g2= 1, 3
                do g3= 1, 3
                 do  i= 1, 6
                   fstrSOLID%STRAIN3D_V2(i,g1,g2,g3,count) = STRAIN_V2(i,g1,g2,g3)
                 end do
                end do
               end do
              end do
            end if
!!!!!
          else if ( ic_type.eq.352) then
            call UPDATE_C3D15TL(ILOG,xx,yy,zz,disp,ddisp,                &
                               ee,pp,fstrSOLID%STRESS3D(:,:,:,:,count), &
                               fstrSOLID%idx_elpl,fstrSOLID%h_dash,     &
                               fstrSOLID%FLG_ELPL(:,:,:,count),         &
                               fstrSOLID%STRESS_YLD(:,:,:,count),       &
                               fstrSOLID%STRAIN3D_P(:,:,:,:,count),     &
                               fstrSOLID%EP_BAR(:,:,:,count),iter)
          else if ( ic_type.eq.362) then
            call UPDATE_C3D20TL(ILOG,xx,yy,zz,disp,ddisp,                &
                               ee,pp,fstrSOLID%STRESS3D(:,:,:,:,count), &
                               fstrSOLID%idx_elpl,fstrSOLID%h_dash,     &
                               fstrSOLID%FLG_ELPL(:,:,:,count),         &
                               fstrSOLID%STRESS_YLD(:,:,:,count),       &
                               fstrSOLID%STRAIN3D_P(:,:,:,:,count),     &
                               fstrSOLID%EP_BAR(:,:,:,count),iter)
          endif
        enddo
      enddo
!C*
!C* Update Displacement
!C*
      do i=1,hecMESH%n_node * ndof
        fstrSOLID%TOTAL_DISP(i) = fstrSOLID%TOTAL_DISP(i)+hecMAT%X(i)
        fstrSOLID%DISP(i)       =                         hecMAT%X(i)
      enddo
!C
!C Message
!C
      if( hecMESH%my_rank.eq.0) then
!!       write(IMSG,*) '####fstrNLGEOM_Update finished'
      end if
      return
      end subroutine fstrNLGEOM_Update
!C----------------------------------------------------------------------*
      subroutine fstrNLGEOM_Output(hecMESH,hecMAT,fstrSOLID,tt,fstrPARAM,DT0,iter, i_step)
!C----------------------------------------------------------------------*
!C*******************************************
!C   Output displacement
!C*******************************************
      use m_fstr

      implicit none

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

      real(kind=kreal) tt
      real(kind=kreal) :: DT0
      integer(kind=kint) :: iter, i_step 
! Local variable
      integer(kind=kint) i,jj,ii,k,j 
      write(ILOG,'(''#### Resul time='',E12.4)') tt
!C
!C-- DISPLACEMENT
!C
      if( hecMESH%n_dof .eq. 3 ) then
        do i= 1, hecMESH%n_node
          fstrSOLID%DISP(3*i-2)=fstrSOLID%TOTAL_DISP(3*i-2)
          fstrSOLID%DISP(3*i-1)=fstrSOLID%TOTAL_DISP(3*i-1)
          fstrSOLID%DISP(3*i  )=fstrSOLID%TOTAL_DISP(3*i  )
        enddo
      else if( hecMESH%n_dof .eq. 2 ) then
        do i= 1, hecMESH%n_node
          fstrSOLID%DISP(3*i-2)=fstrSOLID%TOTAL_DISP(2*i-1)
          fstrSOLID%DISP(3*i-1)=fstrSOLID%TOTAL_DISP(2*i  )
          fstrSOLID%DISP(3*i  )=0.0
        enddo
      endif
!
     if( fstrSOLID%iout_list(1) .eq. 1 ) then
      if( hecMESH%n_dof .eq. 3 ) then
        write(ILOG,*) '#### TOTAL DISPLACEMENT 3D'
        do i= 1, hecMESH%nn_internal
          jj=fstrPARAM%global_local_id(1,i)
          ii=fstrPARAM%global_local_id(2,i)
          write(ILOG,'(i10,3e15.7)') jj,(fstrSOLID%TOTAL_DISP(3*ii-3+k),k=1,3)
        enddo
      else if( hecMESH%n_dof .eq. 2 ) then
        write(ILOG,*) '#### TOTAL DISPLACEMENT 2D'
        do i= 1, hecMESH%nn_internal
          jj=fstrPARAM%global_local_id(1,i)
          ii=fstrPARAM%global_local_id(2,i)
          write(ILOG,'(i10,2e15.7)') jj,(fstrSOLID%TOTAL_DISP(2*ii-2+k),k=1,2)
        enddo
      endif
     end if
!C
!C-- REACTION FORCE
!C
      if( hecMESH%n_dof .eq. 3 ) then
        write(ILOG,*) '#### REACTION FORCE 3D'
        call fstrNLGEOM_reaction_force3d( hecMESH, fstrSOLID,DT0,iter, i_step)
      else if( hecMESH%n_dof .eq. 2 ) then
        write(ILOG,*) '#### REACTION FORCE 2D'
        call fstrNLGEOM_reaction_force2d( hecMESH, fstrSOLID)
      endif
!C
!C-- NODAL STRESS
!C
      if( hecMESH%n_dof .eq. 3 ) then
        write(ILOG,*) '#### NODAL STRESS 3D'
!!       call fstr2hecmw_mesh_conv(hecMESH)
        call fstrNLGEOM_nodal_stress3d( hecMESH, fstrSOLID, i_step, tt)
!!       call hecmw2fstr_mesh_conv(hecMESH)
      else if( hecMESH%n_dof .eq. 2 ) then
        write(ILOG,*) '#### NODAL STRESS 2D'
        call fstrNLGEOM_nodal_stress2d( hecMESH, fstrSOLID)
      endif
!C
!C*** Update Strain and Stress
!C
      call hecmw_update_m_R(hecMESH,fstrSOLID%STRAIN,hecMESH%n_node,6)
      call hecmw_update_m_R(hecMESH,fstrSOLID%STRESS,hecMESH%n_node,7)
      call hecmw_update_m_R(hecMESH,fstrSOLID%STRAIN_PRIS,hecMESH%n_node,12)
      call hecmw_update_m_R(hecMESH,fstrSOLID%STRESS_PRIS,hecMESH%n_node,12)
!!!
      call hecmw_update_m_R(hecMESH,fstrSOLID%STRAIN_nd,hecMESH%n_node,6)
      call hecmw_update_m_R(hecMESH,fstrSOLID%STRESS_nd,hecMESH%n_node,7)
      call hecmw_update_m_R(hecMESH,fstrSOLID%STRSURF_nd,hecMESH%n_node,8)
!C
!C*** Show Results
!C
     if( fstrSOLID%iout_list(5) .eq. 1 ) then
      write(ILOG,*) '#### STRAIN@POST'
      do i= 1, hecMESH%nn_internal
        j=hecMESH%global_node_ID(i)
        write(ILOG,'(i10,6e10.2)') j,(fstrSOLID%STRAIN(6*(i-1)+k),k=1,6)
      enddo
      write(ILOG,*) '#### STRAIN_PRINCIPL@POST'
      write(ILOG,*) '#### r1  r2  r3  vect(r1)  vect(r2)  vect(r3)'
      do i= 1, hecMESH%nn_internal
        j=hecMESH%global_node_ID(i)
        write(ILOG,'(i10,12e10.2)') j,(fstrSOLID%STRAIN_PRIS(12*(i-1)+k),k=1,12)
      enddo
     end if
!
     if( fstrSOLID%iout_list(6) .eq. 1 ) then
      write(ILOG,*) '#### STRESS@POST'
      do i= 1, hecMESH%nn_internal
        j=hecMESH%global_node_ID(i)
        write(ILOG,'(i10,7e10.2)') j,(fstrSOLID%STRESS(7*(i-1)+k),k=1,7)
      enddo
      write(ILOG,*) '#### STRESS_PRINCIPL@POST'
      write(ILOG,*) '#### r1  r2  r3  vect(r1)  vect(r2)  vect(r3)'
      do i= 1, hecMESH%nn_internal
        j=hecMESH%global_node_ID(i)
        write(ILOG,'(i10,12e10.2)') j,(fstrSOLID%STRESS_PRIS(12*(i-1)+k),k=1,12)
      enddo
     end if
!!---  20100203   ----------------
      call fstrNLGEOM_output_elem(hecMESH,fstrSOLID,fstrPARAM)
!!
      return
      end subroutine fstrNLGEOM_Output
!C----------------------------------------------------------------------*
      subroutine fstrNLGEOM_reaction_force3d(hecMESH,fstrSOLID,DT0,iter, i_step)
!C----------------------------------------------------------------------*
!C
!C REACTION FORCE 3D
!C
      use m_fstr

      implicit none

      type (hecmwST_local_mesh) :: hecMESH
      type (fstr_solid)      :: fstrSOLID
!C** Local variables
      real(kind=kreal) xx(20),yy(20),zz(20),tt(20),tti(20),local_stf(1830)
      real(kind=kreal) disp(60),qf(60),qf_res(60) , CUR_nod(20), CUR_nod_I(20)
      integer(kind=kint) nodLOCAL(20)
      real(kind=kreal), dimension(:), allocatable :: spcForce
      integer (kind=KINT), dimension(:), allocatable :: id_spc
!!
      integer(kind=kint) NDOF,i,j,icount,ig0,ig,iS0,iE0,ik,in,nspc,itype
      integer(kind=kint) iS,iE,ic_type,nn,icel,iflag,isect,nid,k,l 
      integer(kind=kint) ierror 
      integer(kind=kint) iter, i_step 
      integer(kind=kint) g1, g2, g3, mid 
      real(kind=kreal) ee,pp,rho,alpha,thick 
      real(kind=kreal) DT0, STRAIN_V1(6,3,3,3),STRAIN_V2(6,3,3,3)


      real(kind=kreal) :: BIJ(3,3),GU(3,3),GUB(3,3),ALPHA_B(3),MAT_OTHO1(9) 

!C** Freedom
      data NDOF/3/
!C Allocate array
      allocate ( id_spc( hecMESH%n_node )             ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_reaction_force3d, id_spc>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      do i=1,hecMESH%n_node
        id_spc(i)=0
      enddo
!C
!C  SPC Nodal Check Flag
!C
!
      STRAIN_V1=0.0d0
      STRAIN_V2=0.0d0
!
      icount=0
      do ig0= 1, fstrSOLID%BOUNDARY_ngrp_tot
        ig= fstrSOLID%BOUNDARY_ngrp_ID(ig0)
        iS0= hecMESH%node_group%grp_index(ig-1) + 1
        iE0= hecMESH%node_group%grp_index(ig  )
        do ik= iS0, iE0
          in   = hecMESH%node_group%grp_item(ik)
          if( id_spc(in) .eq. 0 ) then
            icount=icount+1
            id_spc(in)=icount
          endif
        enddo
      enddo
      icount=0
      do i=1,hecMESH%n_node
        if( id_spc(i).ge.1 ) icount=icount+1
      enddo
      nspc=icount
!C Allocate array
      if( nspc .gt. 0 ) then
         allocate ( spcForce( nspc * NDOF )             ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_reaction_force3d, spcForce>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
         do i=1,nspc*NDOF
           spcForce(i)=0.0
         enddo
      else
        deallocate ( id_spc              ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <fstrNLGEOM_reaction_force3d, id_spc>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        return
      endif
!C +-------------------------------+
!C | according to ELEMENT TYPE     |
!C +-------------------------------+
      icount=0
      do itype= 1, hecMESH%n_elem_type
        iS= hecMESH%elem_type_index(itype-1) + 1
        iE= hecMESH%elem_type_index(itype  )
        ic_type= hecMESH%elem_type_item(itype)
        if (.not. hecmw_is_etype_solid(ic_type)) cycle
!C** Set number of nodes	
        nn = hecmw_get_max_node(ic_type)
!C element loop
        do icel= iS, iE
!C** countUP 
            icount=icount+1
!C** node ID
          iS= hecMESH%elem_node_index(icel-1)
          do j=1,nn
            nodLOCAL(j)= hecMESH%elem_node_item (iS+j)
!C** nodal coordinate
            xx(j)=hecMESH%node(3*nodLOCAL(j)-2)
            yy(j)=hecMESH%node(3*nodLOCAL(j)-1)
            zz(j)=hecMESH%node(3*nodLOCAL(j)  )
            tt (j)=fstrSOLID%TEMPC(nodLOCAL(j) )
            tti(j)=fstrSOLID%TEMPI(nodLOCAL(j) )
            CUR_nod  (j)=fstrSOLID%CURE3D_nod  (nodLOCAL(j) )
            CUR_nod_I(j)=fstrSOLID%CURE3D_nod_I(nodLOCAL(j) )
            disp(3*j-2)=fstrSOLID%TOTAL_DISP(3*nodLOCAL(j)-2)
            disp(3*j-1)=fstrSOLID%TOTAL_DISP(3*nodLOCAL(j)-1)
            disp(3*j  )=fstrSOLID%TOTAL_DISP(3*nodLOCAL(j)  )
          enddo
!C** check node id
          iflag=0
          do j=1,nn
            if( id_spc(  nodLOCAL(j) ) .ge. 1 ) iflag=1
          enddo
          if( iflag .eq. 1 ) then
!C** section  ID
            isect= hecMESH%section_ID(icel)
            mid  = hecMESH%section%sect_mat_ID_item(isect) 
!C** Get Properties
            call fstrNLGEOM_get_prop(hecMESH,isect,ee,pp,rho,alpha,thick)
!C** Create local stiffness

!!-- 20100330 ----!
         do i = 1, 3
           ALPHA_B(i) = 0.d0
           do j = 1, 3
             GU (i,j) = 0.d0
             GUB(i,j) = 0.d0
             BIJ(i,j) = 0.d0
           end do
         end do
         do i = 1, 3
           GU (i,i) = 1.d0
         end do
         do i = 1, 9
           MAT_OTHO1(i) = 7.d0
         end do
!
!!----- 20120927 ----!
       if ( fstrSOLID%mat_iso_num(1) .ne. -1) then
         do i = 1, 4
          if ( mid .eq. i) then
            do l=1,nn
             do j = 1, 3
              do k = 1, 3
               fstrSOLID%GBARU(j,k,nodLOCAL(l)) =  fstrSOLID%GBARU_matrl(j,k,i)
              end do
             end do
            end do
          end if
         end do
       else if ( fstrSOLID%mat_iso_num(1) .eq. -1) then
            do l=1,nn
             do j = 1, 3
              do k = 1, 3
               fstrSOLID%GBARU(j,k,nodLOCAL(l)) =  fstrSOLID%GBARU_ALL(j,k,hecMESH%global_node_ID(nodLOCAL(l)))
              end do
             end do
            end do
       end if

         do l=1,nn
          do j = 1, 3
           do k = 1, 3
            GUB(j,k) = GUB(j,k) + fstrSOLID%GBARU(j,k,nodLOCAL(l))/nn 
           end do
          end do
         end do
!
          do i = 1, 3
           do j = 1, 3
            do k = 1, 3
             BIJ(i,j) = BIJ(i,j) + GU (i,k)*GUB(j,k)
            end do
           end do
          end do
!
         if( fstrSOLID%mat_iso(icount) .eq. 1 ) then
           do i = 1, 3
             ALPHA_B(i)   = fstrSOLID%mat_beta(i,1)
           end do
           do i = 1, 9
             MAT_OTHO1(i) = fstrSOLID%mat_otho (i,1)
           end do
         end if
!!-- 20100330 ----!

            if (ic_type.eq.341) then
              call STF_C3D4TL(ILOG,xx,yy,zz,disp,ee,pp,                      &
                            fstrSOLID%STRESS3D(:,:,:,:,icount),local_stf,qf, &
                            fstrSOLID%idx_elpl,fstrSOLID%h_dash,             &
                            fstrSOLID%FLG_ELPL(:,:,:,icount),                &
                            fstrSOLID%STRESS_YLD(:,:,:,icount))
            else if ( ic_type.eq.351) then
              call STF_C3D6TL(ILOG,xx,yy,zz,disp,ee,pp,                      &
                            fstrSOLID%STRESS3D(:,:,:,:,icount),local_stf,qf, &
                            fstrSOLID%idx_elpl,fstrSOLID%h_dash,             &
                            fstrSOLID%FLG_ELPL(:,:,:,icount),                &
                            fstrSOLID%STRESS_YLD(:,:,:,icount))
            else if ( ic_type.eq.361) then
              call STF_C3D8TL(ILOG,xx,yy,zz,disp,ee,pp,                      &
                            fstrSOLID%STRESS3D(:,:,:,:,icount),local_stf,qf, &
                            fstrSOLID%idx_elpl,fstrSOLID%h_dash,             &
                            fstrSOLID%FLG_ELPL(:,:,:,icount),                &
                            fstrSOLID%STRESS_YLD(:,:,:,icount))
            else if ( ic_type.eq.342) then
!!!!!  visco  ---------
            if(fstrSOLID%idx_visco .ge. 4) then
              do g1= 1, 3
               do g2= 1, 3
                do g3= 1, 3
                 do  j= 1, 6
                   STRAIN_V1(j,g1,g2,g3) = fstrSOLID%STRAIN3D_V1(j,g1,g2,g3,icount)
                 end do
                end do
               end do
              end do
            end if
            if(fstrSOLID%idx_visco .eq. 6) then
              do g1= 1, 3
               do g2= 1, 3
                do g3= 1, 3
                 do  j= 1, 6
                   STRAIN_V2(j,g1,g2,g3) = fstrSOLID%STRAIN3D_V2(j,g1,g2,g3,icount)
                 end do
                end do
               end do
              end do
            end if
!!!!!
              call STF_C3D10TL(ILOG,xx,yy,zz,disp,ee,pp,                                            &
                            fstrSOLID%STRESS3D(:,:,:,:,icount),local_stf,qf,qf_res,                 &
                            fstrSOLID%idx_elpl,fstrSOLID%h_dash,                                    &
                            fstrSOLID%FLG_ELPL(:,:,:,icount),fstrSOLID%STRESS_YLD(:,:,:,icount),    &
                            fstrSOLID%idx_temp,                                                     &
                            tt,tti,CUR_nod, CUR_nod_I,                                              &
                            fstrSOLID%ccc_rl,fstrSOLID%ramd_chm,fstrSOLID%ccc_rl1,                  &
                            fstrSOLID%ccc_rl2,fstrSOLID%yeta,fstrSOLID%yeta1,fstrSOLID%yeta2,       &
                            fstrSOLID%idx_visco,DT0,STRAIN_V1,STRAIN_V2,iter,i_step,icount,icel,    &
                            mid,alpha,fstrSOLID%mat_iso(icount),BIJ,ALPHA_B,MAT_OTHO1)
!!!!!
            else if ( ic_type.eq.352) then
              call STF_C3D15TL(ILOG,xx,yy,zz,disp,ee,pp,                      &
                            fstrSOLID%STRESS3D(:,:,:,:,icount),local_stf,qf, &
                            fstrSOLID%idx_elpl,fstrSOLID%h_dash,             &
                            fstrSOLID%FLG_ELPL(:,:,:,icount),                &
                            fstrSOLID%STRESS_YLD(:,:,:,icount))
            else if ( ic_type.eq.362) then
              call STF_C3D20TL(ILOG,xx,yy,zz,disp,ee,pp,                      &
                            fstrSOLID%STRESS3D(:,:,:,:,icount),local_stf,qf, &
                            fstrSOLID%idx_elpl,fstrSOLID%h_dash,             &
                            fstrSOLID%FLG_ELPL(:,:,:,icount),                &
                            fstrSOLID%STRESS_YLD(:,:,:,icount))
            endif
!C*** Add SPC FORCE ****
            do j=1,nn
              nid=id_spc(  nodLOCAL(j) )
              if( nid .ge. 1 ) then
                do k=1,NDOF
                  spcForce(NDOF*(nid-1)+k)=spcForce(NDOF*(nid-1)+k) +qf(NDOF*(j-1)+k)
                enddo
              endif
            enddo
          endif
        enddo
      enddo
!C*** Show
     if( fstrSOLID%iout_list(4) .eq. 1 ) then
      do i=1,hecMESH%nn_internal
        j=hecMESH%global_node_ID(i)
        nid=id_spc(i)
        if( nid.ge.1 ) then
          write(ILOG,'(i8,3e11.3)') j,                        &
                                    spcForce(NDOF*(nid-1)+1), &
                                    spcForce(NDOF*(nid-1)+2), &
                                    spcForce(NDOF*(nid-1)+3)
        endif
      enddo
     end if
!C Deallocate array
      deallocate ( spcForce             ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <fstrNLGEOM_reaction_force3d, spcForce>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      deallocate ( id_spc               ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <fstrNLGEOM_reaction_force3d, id_spc>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      end subroutine fstrNLGEOM_reaction_force3d
!C----------------------------------------------------------------------*
      subroutine fstrNLGEOM_reaction_force2d(hecMESH,fstrSOLID)
!C----------------------------------------------------------------------*
!C
!C REACTION FORCE 2D
!C
      use m_fstr

      implicit none

      type (hecmwST_local_mesh) :: hecMESH
      type (fstr_solid)         :: fstrSOLID
!** Local variables
      real(kind=kreal) xx(8),yy(8),local_stf(136)
      real(kind=kreal) disp(16),qf(16)
      integer(kind=kint) nodLOCAL(8)
      real    (kind=KREAL),dimension(:), allocatable :: spcForce
      integer (kind=KINT), dimension(:), allocatable :: id_spc
!!
      integer(kind=kint) NDOF,i,icount,ig0,ig,iS0,iE0,ik,in,nspc,itype,iS,iE 
      integer(kind=kint) ic_type,nn,icel,j,iflag,isect,id,iset,nid,k 
      integer(kind=kint) ierror 
      real(kind=kreal) ee,pp,rho,alpha,thick,pa1 
!** Freedom
      data NDOF/2/
!C Allocate array
      allocate ( id_spc( hecMESH%n_node )                ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_reaction_force2d, id_spc>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      do i=1,hecMESH%n_node
        id_spc(i)=0
      enddo
!C
!C  SPC Nodal Check Flag
!C 
      icount=0
      do ig0= 1, fstrSOLID%BOUNDARY_ngrp_tot
        ig= fstrSOLID%BOUNDARY_ngrp_ID(ig0)
        iS0= hecMESH%node_group%grp_index(ig-1) + 1
        iE0= hecMESH%node_group%grp_index(ig  )
        do ik= iS0, iE0
          in   = hecMESH%node_group%grp_item(ik)
          if( id_spc(in) .eq. 0 ) then
            icount=icount+1
            id_spc(in)=icount
          endif
        enddo
      enddo
      icount=0
      do i=1,hecMESH%n_node
        if( id_spc(i).ge.1 ) icount=icount+1
      enddo
      nspc=icount
!C Allocate array
      if( nspc .gt. 0 ) then
         allocate ( spcForce( nspc * NDOF )                 ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_reaction_force2d, spcForce>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
         do i=1,nspc*NDOF
           spcForce(i)=0.0
         enddo
      else
        deallocate ( id_spc                  ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <fstrNLGEOM_reaction_force2d, id_spc>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
        return
      endif
!C +-------------------------------+
!C | according to ELEMENT TYPE     |
!C +-------------------------------+
      icount=0
      do itype= 1, hecMESH%n_elem_type
        iS= hecMESH%elem_type_index(itype-1) + 1
        iE= hecMESH%elem_type_index(itype  )
        ic_type= hecMESH%elem_type_item(itype)
        if (.not. hecmw_is_etype_surface(ic_type)) cycle
!C** Set number of nodes	
        nn = hecmw_get_max_node(ic_type)
!C element loop
        do icel= iS, iE
!C** countUP 
            icount=icount+1
!C** node ID
          iS= hecMESH%elem_node_index(icel-1)
          do j=1,nn
            nodLOCAL(j)= hecMESH%elem_node_item (iS+j)
!C** nodal coordinate
            xx(j)=hecMESH%node(3*nodLOCAL(j)-2)
            yy(j)=hecMESH%node(3*nodLOCAL(j)-1)
            disp(2*j-1)=fstrSOLID%TOTAL_DISP(2*nodLOCAL(j)-1)
            disp(2*j  )=fstrSOLID%TOTAL_DISP(2*nodLOCAL(j)  )
          enddo
!C** check node id
          iflag=0
          do j=1,nn
            if( id_spc(  nodLOCAL(j) ) .ge. 1 ) iflag=1
          enddo
          if( iflag .eq. 1 ) then
!C** section  ID
            isect= hecMESH%section_ID(icel)
!C** Get Properties
            call fstrNLGEOM_get_prop(hecMESH,isect,ee,pp,rho,alpha,thick)
!C** Section Data
            id=hecMESH%section%sect_opt(isect)
            if( id.eq.0 ) then
              iset=1
            else if( id.eq.1) then
              iset=0
            else if( id.eq.2) then
              iset=2
            endif
            pa1=1.0
!C** Create local stiffness
            if (ic_type.eq.241) then
              call STF_C2D4TL(ILOG,xx,yy,disp,ee,pp,pa1,                        &
                            iset,fstrSOLID%STRESS2D(:,:,:,icount),local_stf,qf, &
                            fstrSOLID%idx_elpl,fstrSOLID%h_dash,                &
                            fstrSOLID%FLG_ELPL(:,:,:,icount),                   &
                            fstrSOLID%STRESS_YLD(:,:,:,icount))
            else if ( ic_type.eq.242) then
              call STF_C2D8TL(ILOG,xx,yy,disp,ee,pp,pa1,                        &
                            iset,fstrSOLID%STRESS2D(:,:,:,icount),local_stf,qf, &
                            fstrSOLID%idx_elpl,fstrSOLID%h_dash,                &
                            fstrSOLID%FLG_ELPL(:,:,:,icount),                   &
                            fstrSOLID%STRESS_YLD(:,:,:,icount))
            else if ( ic_type.eq.231) then
              call STF_C2D3TL(ILOG,xx,yy,disp,ee,pp,pa1,                        &
                            iset,fstrSOLID%STRESS2D(:,:,:,icount),local_stf,qf, &
                            fstrSOLID%idx_elpl,fstrSOLID%h_dash,                &
                            fstrSOLID%FLG_ELPL(:,:,:,icount),                   &
                            fstrSOLID%STRESS_YLD(:,:,:,icount))
            else if ( ic_type.eq.232) then
              call STF_C2D6TL(ILOG,xx,yy,disp,ee,pp,pa1,                        &
                            iset,fstrSOLID%STRESS2D(:,:,:,icount),local_stf,qf, &
                            fstrSOLID%idx_elpl,fstrSOLID%h_dash,                &
                            fstrSOLID%FLG_ELPL(:,:,:,icount),                   &
                            fstrSOLID%STRESS_YLD(:,:,:,icount))
            endif
!C*** Add SPC FORCE ****
            do j=1,nn
              nid=id_spc( nodLOCAL(j) )
              if( nid .ge. 1 ) then
                do k=1,NDOF
                  spcForce(NDOF*(nid-1)+k)=spcForce(NDOF*(nid-1)+k) +qf(NDOF*(j-1)+k)
                enddo
              endif
            enddo
          endif
        enddo
      enddo
!C*** Show
     if( fstrSOLID%iout_list(4) .eq. 1 ) then
      do i=1,hecMESH%nn_internal
        j=hecMESH%global_node_ID(i)
        nid=id_spc(i)
        if( nid.ge.1 ) then
          write(ILOG,'(i8,2e11.3)') j,spcForce(NDOF*(nid-1)+1), &
                                      spcForce(NDOF*(nid-1)+2) 
        endif
      enddo
     end if
!C Deallocate array
      deallocate ( spcForce                  ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <fstrNLGEOM_reaction_force2d, spcForce>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      deallocate ( id_spc                    ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <fstrNLGEOM_reaction_force2d, id_spc>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      end subroutine fstrNLGEOM_reaction_force2d
!C----------------------------------------------------------------------*
      subroutine fstrNLGEOM_nodal_stress3d(hecMESH,fstrSOLID,i_step,ctime)
!C----------------------------------------------------------------------*
!C
!C NODAL STRESS 3D
!C
      use m_fstr

      implicit none
      type (hecmwST_local_mesh) :: hecMESH
      type (fstr_solid)      :: fstrSOLID
!** Local variables
      real(kind=kreal) xx(20),yy(20),zz(20)
      real(kind=kreal) disp(60)
      real(kind=kreal) tt(20),tt0(20)
      integer(kind=kint) nodLOCAL(20)
      real(kind=kreal) array(13,20),array_pris(24,20)
      real(kind=kreal) array_tmp1(13,27),ctime

      integer(kind=kint) ig0,ig,iS0,iE0,ik,in,ic,itype,iS,iE,ic_type,nn,icel,j 
      integer(kind=kint) isect,k,l,i,ielem,ID_area,mid  
      integer(kind=kint) i_step 
      integer(kind=kint) ierror 
      real(kind=kreal) val,ee,pp,rho,alpha,thick,s11,s22,s33,s12,s23,s13,ps,smises 
!** Array for nodal recovery
      real    (kind=KREAL), dimension(:,:), allocatable :: arrayTotal
      real    (kind=KREAL), dimension(:,:), allocatable :: arrayTotal_pris
      real    (kind=KREAL), dimension(:,:), allocatable :: arrayTotal_nd
      real    (kind=KREAL), dimension(:,:), allocatable :: arrayTotal_SUR_nd
      integer (kind=KINT), dimension(:), allocatable :: number
      real    (kind=KREAL), dimension(:), allocatable :: temp


      real(kind=kreal) :: BIJ(3,3),GU(3,3),GUB(3,3) 

!*Allocate array
      allocate ( arrayTotal( hecMESH%n_node, 13  )         ,STAT=ierror) 
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_nodal_stress3d, arrayTotal>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      allocate ( arrayTotal_pris( hecMESH%n_node, 24  )         ,STAT=ierror) 
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_nodal_stress3d, arrayTotal_pris>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      allocate ( arrayTotal_nd( hecMESH%n_node, 13  )         ,STAT=ierror) 
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_nodal_stress3d, arrayTotal_nd>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      allocate ( arrayTotal_SUR_nd( hecMESH%n_node, 8  )         ,STAT=ierror) 
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_nodal_stress3d, arrayTotal_SUR_nd>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      allocate ( number(hecMESH%n_node)                    ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_nodal_stress3d, number>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
!*ZERO CLEAR 
      arrayTotal=0
      arrayTotal_pris=0
      arrayTotal_nd=0
      arrayTotal_SUR_nd=0
      number=0
!C
!C Set Temperature
!C
      allocate ( temp(hecMESH%n_node)                      ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_nodal_stress3d, temp>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      temp=0
      if( fstrSOLID%TEMP_ngrp_tot > 0 ) then
        do ig0= 1, fstrSOLID%TEMP_ngrp_tot
          ig= fstrSOLID%TEMP_ngrp_ID(ig0)
          val=fstrSOLID%TEMP_ngrp_val(ig0)
!C START & END
          iS0= hecMESH%node_group%grp_index(ig-1) + 1
          iE0= hecMESH%node_group%grp_index(ig  )
           do ik= iS0, iE0
             in   = hecMESH%node_group%grp_item(ik)
             temp( in ) = val
          enddo
        enddo
      endif
!C
!C +-------------------------------+
!C | according to ELEMENT TYPE     |
!C +-------------------------------+
      ic=0
      do itype= 1, hecMESH%n_elem_type
        iS= hecMESH%elem_type_index(itype-1) + 1
        iE= hecMESH%elem_type_index(itype  )
        ic_type= hecMESH%elem_type_item(itype)
        if (.not. hecmw_is_etype_solid(ic_type)) cycle
!C** Set number of nodes	
        nn = hecmw_get_max_node(ic_type)
!C element loop
        do icel= iS, iE
!C** node ID
          iS= hecMESH%elem_node_index(icel-1)
          do j=1,nn
            nodLOCAL(j)= hecMESH%elem_node_item (iS+j)
!C** nodal coordinate
            xx(j)=hecMESH%node(3*nodLOCAL(j)-2)
            yy(j)=hecMESH%node(3*nodLOCAL(j)-1)
            zz(j)=hecMESH%node(3*nodLOCAL(j)  )
!!            tt(j)=temp( nodLOCAL(j) )
!!            tt0(j)=ref_temp
            disp(3*j-2)=fstrSOLID%TOTAL_DISP(3*nodLOCAL(j)-2)
            disp(3*j-1)=fstrSOLID%TOTAL_DISP(3*nodLOCAL(j)-1)
            disp(3*j  )=fstrSOLID%TOTAL_DISP(3*nodLOCAL(j)  )
          enddo  
!C** section  ID
          isect= hecMESH%section_ID(icel)
          mid  = hecMESH%section%sect_mat_ID_item(isect) 
!C** Get Properties
          call fstrNLGEOM_get_prop(hecMESH,isect,ee,pp,rho,alpha,thick)
!C** Create local stiffness
          ic=ic+1
!!-- 20100330 ----!
         do i = 1, 3
           do j = 1, 3
             GU (i,j) = 0.d0
             GUB(i,j) = 0.d0
             BIJ(i,j) = 0.d0
           end do
         end do
         do i = 1, 3
           GU (i,i) = 1.d0
         end do
!
!!----- 20120927 ----!
       if ( fstrSOLID%mat_iso_num(1) .ne. -1) then
         do i = 1, 4
          if ( mid .eq. i) then
            do l=1,nn
             do j = 1, 3
              do k = 1, 3
               fstrSOLID%GBARU(j,k,nodLOCAL(l)) =  fstrSOLID%GBARU_matrl(j,k,i)
              end do
             end do
            end do
          end if
         end do
       else if ( fstrSOLID%mat_iso_num(1) .eq. -1) then
            do l=1,nn
             do j = 1, 3
              do k = 1, 3
               fstrSOLID%GBARU(j,k,nodLOCAL(l)) =  fstrSOLID%GBARU_ALL(j,k,hecMESH%global_node_ID(nodLOCAL(l)))
              end do
             end do
            end do
       end if

         do l=1,nn
          do j = 1, 3
           do k = 1, 3
            GUB(j,k) = GUB(j,k) + fstrSOLID%GBARU(j,k,nodLOCAL(l))/nn 
           end do
          end do
         end do
!
          do i = 1, 3
           do j = 1, 3
            do k = 1, 3
             BIJ(i,j) = BIJ(i,j) + GU (i,k)*GUB(j,k)
            end do
           end do
          end do
!!-- 20100330 ----!

          if (ic_type.eq.341) then
            call NRCV_C3D4TL(ILOG,xx,yy,zz,tt,tt0,ee,pp,alpha,       &
                          disp,fstrSOLID%STRESS3D(:,:,:,:,ic),array)
          else if ( ic_type.eq.351) then
            call NRCV_C3D6TL(ILOG,xx,yy,zz,tt,tt0,ee,pp,alpha,       &
                          disp,fstrSOLID%STRESS3D(:,:,:,:,ic),array)
          else if ( ic_type.eq.361) then
            call NRCV_C3D8TL(ILOG,xx,yy,zz,tt,tt0,ee,pp,alpha,       &
                         disp,fstrSOLID%STRESS3D(:,:,:,:,ic),array)
          else if ( ic_type.eq.342) then
            call NRCV_C3D10TL(ILOG,xx,yy,zz,tt,tt0,ee,pp,alpha,      &
                         disp,fstrSOLID%STRESS3D(:,:,:,:,ic),array,  &
                         fstrSOLID%SURFF(:,:,:,:,ic),                &
                         fstrSOLID%STRAIN3D_P(:,:,:,:,ic),ic,i_step, &
                         fstrSOLID%idx_stres,array_tmp1,array_pris,  &
                         fstrSOLID%idx_out_cood_st,fstrSOLID%mat_iso(ic),BIJ, & 
                         hecMESH%global_elem_ID(ic))
          else if ( ic_type.eq.352) then
            call NRCV_C3D15TL(ILOG,xx,yy,zz,tt,tt0,ee,pp,alpha,      &
                         disp,fstrSOLID%STRESS3D(:,:,:,:,ic),array)
          else if ( ic_type.eq.362) then
            call NRCV_C3D20TL(ILOG,xx,yy,zz,tt,tt0,ee,pp,alpha,      &
                         disp,fstrSOLID%STRESS3D(:,:,:,:,ic),array)
          endif

!!---  20100203   ----------------
          ielem = ic
          ID_area = hecMESH%elem_ID(ic*2)
!!------  output of stress and strain of int. point at last step  -------!!
          if(i_step .eq. INCMAX ) then
           if( ID_area.eq.hecMESH%my_rank ) then
             do j=1,27
              if( j .eq. 1 ) then
                write(IDBG,'(i10,1pe13.4e3,'' ic ='',i10,'' global_elem_ID(ic) ='',i10)')  &
                      i_step, ctime, ic, hecMESH%global_elem_ID(ic)
                write(IDBG,'(''  nodLOCAL(k) = '',10i10)') (hecMESH%global_node_ID(nodLOCAL(k)),k=1,10)
                write(IDBG,'(33x,i10,1p13e13.4e3)'           )       &
                       j,(array_tmp1(k,j),k=1,13)
              else if( j .ge. 2 ) then
                write(IDBG,'(33x,i10,1p13e13.4e3)'           )       &
                       j,(array_tmp1(k,j),k=1,13)
              end if
             end do
           end if
          end if
!!
!!!       if( ID_area.eq.hecMESH%my_rank ) then
            do j=1,6
             fstrSOLID%ESTRAIN(6*ielem-(j-1)) = 0.0D0
            end do
            do j=1,7
             fstrSOLID%ESTRESS(7*ielem-(j-1)) = 0.0D0
            end do
!!!       end if
!!
          do j=1,nn
            do k=1,13                                             
              arrayTotal(nodLOCAL(j),k)= arrayTotal(nodLOCAL(j),k)+array(k,j)

              if(k .le. 6) then
                arrayTotal_nd(nodLOCAL(j),k)= arrayTotal_nd(nodLOCAL(j),k)+fstrSOLID%STRESS3D_ndP1(k,j,ic)
              else if(k .ge. 8 .and. k .le. 13) then
                arrayTotal_nd(nodLOCAL(j),k)= arrayTotal_nd(nodLOCAL(j),k)+fstrSOLID%STRAIN3D_P_nd(k-7,j,ic)
              end if
              if(k .le. 8) then
                arrayTotal_SUR_nd(nodLOCAL(j),k)= arrayTotal_SUR_nd(nodLOCAL(j),k)+fstrSOLID%SURFF_nd(k,j,ic)
              end if

            enddo

            do k=1,24                                             
              arrayTotal_pris(nodLOCAL(j),k)= arrayTotal_pris(nodLOCAL(j),k)+array_pris(k,j)
            enddo

            number( nodLOCAL(j) )=number( nodLOCAL(j) )+1

!!!!!
!!!         if( ID_area.eq.hecMESH%my_rank ) then
              do k=1,6                                             
               fstrSOLID%ESTRAIN(6*(ielem-1)+k) = fstrSOLID%ESTRAIN(6*(ielem-1)+k) + array(k,j)/nn
              end do
              do k=1,6                                             
               fstrSOLID%ESTRESS(7*(ielem-1)+k) = fstrSOLID%ESTRESS(7*(ielem-1)+k) + array(6+k,j)/nn
              end do
!!!         end if
!!!!!
          enddo
!!!!!
!!!       if( ID_area.eq.hecMESH%my_rank ) then
            s11=fstrSOLID%ESTRESS(7*ielem-6)
            s22=fstrSOLID%ESTRESS(7*ielem-5)
            s33=fstrSOLID%ESTRESS(7*ielem-4)
            s12=fstrSOLID%ESTRESS(7*ielem-3)
            s23=fstrSOLID%ESTRESS(7*ielem-2)
            s13=fstrSOLID%ESTRESS(7*ielem-1)
            ps=(s11+s22+s33)/3.0
            smises=0.5*( (s11-ps)**2+(s22-ps)**2+(s33-ps)**2 )    &
                  +s12**2+s23**2+s13**2
            fstrSOLID%ESTRESS(7*ielem)=sqrt(3.0*smises)
!!!       end if
!!!!!
        enddo
      enddo
!** Average over nodes
      do i=1,hecMESH%n_node
        do j=1,13
          arrayTotal(i,j)=arrayTotal(i,j)/FLOAT(number(i))
          arrayTotal_nd(i,j)=arrayTotal_nd(i,j)/FLOAT(number(i))
          if(j .le. 8) then
            arrayTotal_SUR_nd(i,j)=arrayTotal_SUR_nd(i,j)/FLOAT(number(i))
          end if
        enddo
        do j=1,24
          arrayTotal_pris(i,j)=arrayTotal_pris(i,j)/FLOAT(number(i))
        enddo
      enddo
!** RE-CALCULATE von MISES stress
      do i=1,hecMESH%n_node
        s11=arrayTotal(i,7)
        s22=arrayTotal(i,8)
        s33=arrayTotal(i,9)
        s12=arrayTotal(i,10)
        s23=arrayTotal(i,11)
        s13=arrayTotal(i,12)
        ps=(s11+s22+s33)/3.0
        smises=0.5*( (s11-ps)**2+(s22-ps)**2+(s33-ps)**2 ) +s12**2+s23**2+s13**2
        arrayTotal(i,13)=sqrt(3.0*smises)
      enddo
      do i=1,hecMESH%n_node
        s11=arrayTotal_nd(i,1)
        s22=arrayTotal_nd(i,2)
        s33=arrayTotal_nd(i,3)
        s12=arrayTotal_nd(i,4)
        s23=arrayTotal_nd(i,5)
        s13=arrayTotal_nd(i,6)
        ps=(s11+s22+s33)/3.0
        smises=0.5*( (s11-ps)**2+(s22-ps)**2+(s33-ps)**2 ) +s12**2+s23**2+s13**2
        arrayTotal_nd(i,7)=sqrt(3.0*smises)
      enddo
!** Set Array 
      do i=1,hecMESH%n_node
        do j=1,6
          fstrSOLID%STRAIN(6*(i-1)+j)=arrayTotal(i,j)
        enddo
        do j=1,7
          fstrSOLID%STRESS(7*(i-1)+j)=arrayTotal(i,j+6)
        enddo
        do j=1,12
          fstrSOLID%STRAIN_PRIS(12*(i-1)+j)=arrayTotal_pris(i,j)
          fstrSOLID%STRESS_PRIS(12*(i-1)+j)=arrayTotal_pris(i,j+12)
        enddo
      enddo
!
      do i=1,hecMESH%n_node
        do j=1,6
          fstrSOLID%STRAIN_nd(6*(i-1)+j)=arrayTotal_nd(i,j+7)
        enddo
        do j=1,7
          fstrSOLID%STRESS_nd(7*(i-1)+j)=arrayTotal_nd(i,j  )
        enddo
      enddo
      do i=1,hecMESH%n_node
        do j=1,8
          fstrSOLID%STRSURF_nd(8*(i-1)+j)=arrayTotal_SUR_nd(i,j)
        enddo
      enddo
!*Deallocate array
      deallocate( arrayTotal                 ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <fstrNLGEOM_nodal_stress3d, arrayTotal>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      deallocate( arrayTotal_pris                 ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <fstrNLGEOM_nodal_stress3d, arrayTotal_pris>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      deallocate( arrayTotal_nd                 ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <fstrNLGEOM_nodal_stress3d, arrayTotal_nd>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      deallocate( arrayTotal_SUR_nd                 ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <fstrNLGEOM_nodal_stress3d, arrayTotal_SUR_nd>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      deallocate( number                     ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <fstrNLGEOM_nodal_stress3d, number>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      deallocate( temp                       ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <fstrNLGEOM_nodal_stress3d, temp>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      return
      end subroutine fstrNLGEOM_nodal_stress3d
!C----------------------------------------------------------------------*
      subroutine fstrNLGEOM_nodal_stress2d(hecMESH,fstrSOLID)
!C----------------------------------------------------------------------*
!C
!C NODAL STRESS 2D
!C
      use m_fstr

      implicit none

      type (hecmwST_local_mesh) :: hecMESH
      type (fstr_solid)      :: fstrSOLID
!** Local variables
      real(kind=kreal) xx(8),yy(8)
      real(kind=kreal) disp(16)
      real(kind=kreal) tt(8),tt0(8)
      integer(kind=kint) nodLOCAL(8)
      real(kind=kreal) array(9,8)

      integer(kind=kint) ig0,ig,iS0,iE0,ik,in,ic,itype,iS,iE,ic_type,nn,icel,j,isect 
      integer(kind=kint) id,iset,k,i 
      integer(kind=kint) ierror 
      real(kind=kreal) val,ee,pp,rho,alpha,thick,pa1,s11,s22,s33,s12,s23,s13,ps,smises 
!** Array for nodal recovery
      real    (kind=KREAL), dimension(:,:), allocatable :: arrayTotal
      integer (kind=KINT), dimension(:), allocatable :: number
      real (kind=KREAL), dimension(:), allocatable :: temp
!*Allocate array
      allocate ( arrayTotal( hecMESH%n_node, 9  )          ,STAT=ierror) 
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_nodal_stress2d, arrayTotal>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      allocate ( number(hecMESH%n_node)                    ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_nodal_stress2d, number>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
!*ZERO CLEAR 
      arrayTotal=0
      number=0
!C
!C Set Temperature
!C
      allocate ( temp(hecMESH%n_node)                      ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to allocation error <fstrNLGEOM_nodal_stress2d, temp>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      temp=0
      if( fstrSOLID%TEMP_ngrp_tot > 0 ) then
        do ig0= 1, fstrSOLID%TEMP_ngrp_tot
          ig= fstrSOLID%TEMP_ngrp_ID(ig0)
          val=fstrSOLID%TEMP_ngrp_val(ig0)
!C START & END
          iS0= hecMESH%node_group%grp_index(ig-1) + 1
          iE0= hecMESH%node_group%grp_index(ig  )
           do ik= iS0, iE0
             in   = hecMESH%node_group%grp_item(ik)
             temp( in ) = val
          enddo
        enddo
      endif
!C
!C +-------------------------------+
!C | according to ELEMENT TYPE     |
!C +-------------------------------+
      ic=0
      do itype= 1, hecMESH%n_elem_type
        iS= hecMESH%elem_type_index(itype-1) + 1
        iE= hecMESH%elem_type_index(itype  )
        ic_type= hecMESH%elem_type_item(itype)
        if (.not. hecmw_is_etype_surface(ic_type)) cycle
!C** Set number of nodes	
        nn = hecmw_get_max_node(ic_type)
!C element loop
        do icel= iS, iE
!C** node ID
          iS= hecMESH%elem_node_index(icel-1)
          do j=1,nn
            nodLOCAL(j)= hecMESH%elem_node_item (iS+j)
!C** nodal coordinate
            xx(j)=hecMESH%node(3*nodLOCAL(j)-2)
            yy(j)=hecMESH%node(3*nodLOCAL(j)-1)
            tt(j)=temp( nodLOCAL(j) )
            tt0(j)=ref_temp
            disp(2*j-1)=fstrSOLID%TOTAL_DISP(2*nodLOCAL(j)-1)
            disp(2*j  )=fstrSOLID%TOTAL_DISP(2*nodLOCAL(j)  )
          enddo  
!C** section  ID
          isect= hecMESH%section_ID(icel)
!C** Get Properties
          call fstrNLGEOM_get_prop(hecMESH,isect,ee,pp,rho,alpha,thick)
!C** Section Data
          id=hecMESH%section%sect_opt(isect)
          if( id.eq.0 ) then
            iset=1
          else if( id.eq.1) then
            iset=0
          else if( id.eq.2) then
            iset=2
          endif
          pa1=1.0     
!C** Create local stiffness
          ic=ic+1
          if (ic_type.eq.241) then
            call NRCV_C2D4TL(ILOG,xx,yy,tt,tt0,ee,pp,alpha,                 &
                          pa1,iset,disp,fstrSOLID%STRESS2D(:,:,:,ic),array)
          else if ( ic_type.eq.242) then
            call NRCV_C2D8TL(ILOG,xx,yy,tt,tt0,ee,pp,alpha,                 &
                          pa1,iset,disp,fstrSOLID%STRESS2D(:,:,:,ic),array)
          else if ( ic_type.eq.231) then
            call NRCV_C2D3TL(ILOG,xx,yy,tt,tt0,ee,pp,alpha,                 &
                          pa1,iset,disp,fstrSOLID%STRESS2D(:,:,:,ic),array)
          else if ( ic_type.eq.232) then
            call NRCV_C2D6TL(ILOG,xx,yy,tt,tt0,ee,pp,alpha,                 &
                          pa1,iset,disp,fstrSOLID%STRESS2D(:,:,:,ic),array)
          endif
          do j=1,nn
            do k=1,9
              arrayTotal(nodLOCAL(j),k)= arrayTotal(nodLOCAL(j),k)+array(k,j)
            enddo
            number( nodLOCAL(j) )=number( nodLOCAL(j) )+1    
          enddo 
        enddo
      enddo
!** Average over nodes
      do i=1,hecMESH%n_node
        do j=1,9
          arrayTotal(i,j)=arrayTotal(i,j)/FLOAT(number(i))
        enddo
      enddo
!** Show
!C      do i=1,hecMESH%n_node
!C       j=hecMESH%global_node_ID(i)
!C       write(*,'(i8)') j
!C       write(*,'(''E '',4e11.3)') (arrayTotal(i,k),k=1,4)
!C       write(*,'(''S '',5e11.3)') (arrayTotal(i,k),k=5,9)
!C      enddo
!** RE-CALCULATE von MISES stress
      do i=1,hecMESH%n_node
        s11=arrayTotal(i,5)
        s22=arrayTotal(i,6)
        s33=arrayTotal(i,8)
        s12=arrayTotal(i,7)
        s23=0.0
        s13=0.0
        ps=(s11+s22+s33)/3.0
        smises=0.5*( (s11-ps)**2+(s22-ps)**2+(s33-ps)**2 ) +s12**2+s23**2+s13**2
        arrayTotal(i,9)=sqrt(3.0*smises)
      enddo
!** Set Array 
      do i=1,hecMESH%n_node
        fstrSOLID%STRAIN(6*(i-1)+1)=arrayTotal(i,1)
        fstrSOLID%STRAIN(6*(i-1)+2)=arrayTotal(i,2)
        fstrSOLID%STRAIN(6*(i-1)+3)=arrayTotal(i,4)
        fstrSOLID%STRAIN(6*(i-1)+4)=arrayTotal(i,3)
        fstrSOLID%STRAIN(6*(i-1)+5)=0.0
        fstrSOLID%STRAIN(6*(i-1)+6)=0.0
        fstrSOLID%STRESS(7*(i-1)+1)=arrayTotal(i,5)
        fstrSOLID%STRESS(7*(i-1)+2)=arrayTotal(i,6)
        fstrSOLID%STRESS(7*(i-1)+3)=arrayTotal(i,8)
        fstrSOLID%STRESS(7*(i-1)+4)=arrayTotal(i,7)
        fstrSOLID%STRESS(7*(i-1)+5)=0.0
        fstrSOLID%STRESS(7*(i-1)+6)=0.0
        fstrSOLID%STRESS(7*(i-1)+7)=arrayTotal(i,9)
      enddo
!C**Deallocate array
      deallocate( arrayTotal                            ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <fstrNLGEOM_nodal_stress2d, arrayTotal>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      deallocate( number                                ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <fstrNLGEOM_nodal_stress2d, number>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      deallocate( temp                                  ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <fstrNLGEOM_nodal_stress2d, temp>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      return
      end subroutine fstrNLGEOM_nodal_stress2d
!C----------------------------------------------------------------------*
      subroutine fstrNLGEOM_output_elem(hecMESH,fstrSOLID,fstrPARAM)
!C----------------------------------------------------------------------*
!     output of strain & stress in element
!!---  20100203  ---------------------

      use m_fstr
      implicit none
      type (hecmwST_local_mesh) :: hecMESH
      type (fstr_solid        ) :: fstrSOLID
      type (fstr_param        ) :: fstrPARAM

!** Local variables
      integer(kind=kint) i,k,ii,jj,ID_area   
!!

        write(ILOG,*) '#### 3D STRAIN @Element'
        write(ILOG,'(a,a)')                                            &
                    '  ELEMENT     E11         E22         E33     '   &
                             ,'    E12         E23         E31     '
        write(ILOG,'(a,a)')                                            &
                    '  -------+-----------+-----------+-----------+'   &
                             ,'-----------+-----------+-------------'
        do i=1,hecMESH%n_elem
          jj = hecMESH%global_elem_ID(i)
!!!       ii = hecMESH%elem_ID(i*2-1)
          ii = i
          ID_area = hecMESH%elem_ID(i*2)
          if( ID_area.eq.hecMESH%my_rank ) then
            write(ILOG,'(i10,1p6e12.4)') jj,                           &
                                 (fstrSOLID%ESTRAIN(6*(ii-1)+k),k=1,6)
          endif
        enddo
        write(ILOG,*) '#### 3D STRESS @Element'
        write(ILOG,'(a,a)')                                            &
                    '  ELEMENT     S11         S22         S33     '   &
                 ,'    S12         S23         S31        MISES'
        write(ILOG,'(a,a)')                                            &
                    '  -------+-----------+-----------+-----------+'   &
                 ,'-----------+-----------+-----------+-------------'
        do i=1,hecMESH%n_elem
          jj = hecMESH%global_elem_ID(i)
!!!       ii = hecMESH%elem_ID(i*2-1)
          ii = i
          ID_area = hecMESH%elem_ID(i*2)
          if( ID_area.eq.hecMESH%my_rank ) then
            write(ILOG,'(i10,1p7e12.4)') jj,                            &
                                 (fstrSOLID%ESTRESS(7*(ii-1)+k),k=1,7)
          endif
        enddo

      return
      end subroutine fstrNLGEOM_output_elem 
!C----------------------------------------------------------------------*
      subroutine fstrNLGEOM_Post(hecMESH,hecMAT,fstrSOLID,tt,i_step)
!C----------------------------------------------------------------------*
!C*******************************************
!C  Post process
!C*******************************************
      use m_fstr

      implicit none

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

      integer(kind=kint) i_step 
      real(kind=kreal) tt
!C file name
      character(len=HECMW_HEADER_LEN) :: header
      character(len=HECMW_NAME_LEN) :: label
      character(len=HECMW_NAME_LEN) :: nameID 
!** Local variables
!C Max Min
      real(kind=kreal) Umax(3),Umin(3),Emax(6),Emin(6),Smax(7),Smin(7)
      real(kind=kreal) E_PRmax(3),E_PRmin(3),S_PRmax(3),S_PRmin(3)
      integer(kind=kint) IUmax(3),IUmin(3),IEmax(6),IEmin(6),ISmax(7),ISmin(7)
      integer(kind=kint) IE_PRmax(3),IE_PRmin(3),IS_PRmax(3),IS_PRmin(3)

      integer(kind=kint) j,i,k,ii,MDOF,ID_area,nd,id,ndof 

      real(kind=kreal)    EEmax(6), EEmin(6) ,ESmax(7), ESmin(7)
      integer(kind=kint) IEEmax(6),IEEmin(6),IESmax(7),IESmin(7)
!C*-------- solver control -----------*
      logical :: ds = .false. !using Direct Solver or not

! in case of direct solver
      if (hecMAT%Iarray(99) .eq. 2) then
        ds = .true.
      end if

!C Time
      write(ILOG,'(''#### Resul time='',E12.4)') tt
!C Clear
      Umax=0.0d0
      Umin=0.0d0
      Emax=0.0d0
      Emin=0.0d0
      Smax=0.0d0
      Smin=0.0d0
       E_PRmax=0.0d0
       E_PRmin=0.0d0
       S_PRmax=0.0d0
       S_PRmin=0.0d0
      IUmax=0
      IUmin=0
      IEmax=0
      IEmin=0
      ISmax=0
      ISmin=0
       IE_PRmax=0
       IE_PRmin=0
       IS_PRmax=0
       IS_PRmin=0

      EEmax=0.0d0
      EEmin=0.0d0
      ESmax=0.0d0
      ESmin=0.0d0
      IEEmax=0
      IEEmin=0
      IESmax=0
      IESmin=0
!C*** Displacement 
      write(ILOG,*) '#### DISPLACEMENT@POST'
      do i= 1, hecMESH%nn_internal
        j=hecMESH%global_node_ID(i)
!C Max & Min
        if( i.eq.1) then
          do k=1,3 
            Umax(k)=fstrSOLID%DISP(3*i-3+k)
            Umin(k)=fstrSOLID%DISP(3*i-3+k)
            IUmax(k)=j
            IUmin(k)=j
          enddo
        else
          do k=1,3 
            if( fstrSOLID%DISP(3*i-3+k) .gt. Umax(k) ) then
              Umax(k)=fstrSOLID%DISP(3*i-3+k)
              IUmax(k)=j
            endif
            if( fstrSOLID%DISP(3*i-3+k) .lt. Umin(k) ) then
              Umin(k)=fstrSOLID%DISP(3*i-3+k)
              IUmin(k)=j
            endif
          enddo
        endif
      enddo
!C***Strain
      write(ILOG,*) '#### STRAIN@POST'
!! @node
      do i= 1, hecMESH%nn_internal
        j=hecMESH%global_node_ID(i)
!C Max & Min
        if( i.eq.1) then
          do k=1,6 
            Emax(k)=fstrSOLID%STRAIN(6*i-6+k)
            Emin(k)=fstrSOLID%STRAIN(6*i-6+k)
            IEmax(k)=j
            IEmin(k)=j
          enddo
          do k=1,3 
            E_PRmax(k)=fstrSOLID%STRAIN_PRIS(12*i-12+k)
            E_PRmin(k)=fstrSOLID%STRAIN_PRIS(12*i-12+k)
            IE_PRmax(k)=j
            IE_PRmin(k)=j
          enddo
        else
          do k=1,6
            if( fstrSOLID%STRAIN(6*i-6+k) .gt. Emax(k) ) then
              Emax(k)=fstrSOLID%STRAIN(6*i-6+k)
              IEmax(k)=j
            endif
            if( fstrSOLID%STRAIN(6*i-6+k) .lt. Emin(k) ) then
              Emin(k)=fstrSOLID%STRAIN(6*i-6+k)
              IEmin(k)=j
            endif
          enddo
          do k=1,3
            if( fstrSOLID%STRAIN_PRIS(12*i-12+k) .gt. E_PRmax(k) ) then
              E_PRmax(k)=fstrSOLID%STRAIN_PRIS(12*i-12+k)
              IE_PRmax(k)=j
            endif
            if( fstrSOLID%STRAIN_PRIS(12*i-12+k) .lt. E_PRmin(k) ) then
              E_PRmin(k)=fstrSOLID%STRAIN_PRIS(12*i-12+k)
              IE_PRmin(k)=j
            endif
          enddo
        endif
      enddo
!! @element
      MDOF = 6
      do i= 1, hecMESH%n_elem
!!!     ii = hecMESH%elem_ID(i*2-1)
        ii = i
        ID_area = hecMESH%elem_ID(i*2)
!       if( ID_area.eq.hecMESH%my_rank ) then
          j=hecMESH%global_elem_ID(ii)
          if( i.eq.1) then
            do k=1,MDOF
              EEmax(k)=fstrSOLID%ESTRAIN(MDOF*(ii-1)+k)
              EEmin(k)=fstrSOLID%ESTRAIN(MDOF*(ii-1)+k)
              IEEmax(k)=j
              IEEmin(k)=j
            enddo
          else
            do k=1,MDOF
              if( fstrSOLID%ESTRAIN(MDOF*(ii-1)+k) .gt. EEmax(k) ) then
                EEmax(k)=fstrSOLID%ESTRAIN(MDOF*(ii-1)+k)
                IEEmax(k)=j
              endif
              if( fstrSOLID%ESTRAIN(MDOF*(ii-1)+k) .lt. EEmin(k) ) then
                EEmin(k)=fstrSOLID%ESTRAIN(MDOF*(ii-1)+k)
                IEEmin(k)=j
              endif
            enddo
          endif
!       endif
      enddo
!C*** Stress
      write(ILOG,*) '#### STRESS@POST'
!! @node
      do i= 1, hecMESH%nn_internal
        j=hecMESH%global_node_ID(i)
!C Max & Min
        if( i.eq.1) then
          do k=1,7
            Smax(k)=fstrSOLID%STRESS(7*i-7+k)
            Smin(k)=fstrSOLID%STRESS(7*i-7+k)
            ISmax(k)=j
            ISmin(k)=j
          enddo
          do k=1,3 
            S_PRmax(k)=fstrSOLID%STRESS_PRIS(12*i-12+k)
            S_PRmin(k)=fstrSOLID%STRESS_PRIS(12*i-12+k)
            IS_PRmax(k)=j
            IS_PRmin(k)=j
          enddo
        else
          do k=1,7
            if( fstrSOLID%STRESS(7*i-7+k) .gt. Smax(k) ) then
              Smax(k)=fstrSOLID%STRESS(7*i-7+k)
              ISmax(k)=j
            endif
            if( fstrSOLID%STRESS(7*i-7+k) .lt. Smin(k) ) then
              Smin(k)=fstrSOLID%STRESS(7*i-7+k)
              ISmin(k)=j
            endif
          enddo
          do k=1,3
            if( fstrSOLID%STRESS_PRIS(12*i-12+k) .gt. S_PRmax(k) ) then
              S_PRmax(k)=fstrSOLID%STRESS_PRIS(12*i-12+k)
              IS_PRmax(k)=j
            endif
            if( fstrSOLID%STRESS_PRIS(12*i-12+k) .lt. S_PRmin(k) ) then
              S_PRmin(k)=fstrSOLID%STRESS_PRIS(12*i-12+k)
              IS_PRmin(k)=j
            endif
          enddo
        endif
      enddo
!! @element
      MDOF = 7
      do i= 1, hecMESH%n_elem
!!!     ii = hecMESH%elem_ID(i*2-1)
        ii = i
        ID_area = hecMESH%elem_ID(i*2)
!       if( ID_area.eq.hecMESH%my_rank ) then
          j=hecMESH%global_elem_ID(ii)
          if( i.eq.1) then
            do k=1,MDOF
              ESmax(k)=fstrSOLID%ESTRESS(MDOF*(ii-1)+k)
              ESmin(k)=fstrSOLID%ESTRESS(MDOF*(ii-1)+k)
              IESmax(k)=j
              IESmin(k)=j
            enddo
          else
            do k=1,MDOF
              if( fstrSOLID%ESTRESS(MDOF*(ii-1)+k) .gt. ESmax(k) ) then
                ESmax(k)=fstrSOLID%ESTRESS(MDOF*(ii-1)+k)
                IESmax(k)=j
              endif
              if( fstrSOLID%ESTRESS(MDOF*(ii-1)+k) .lt. ESmin(k) ) then
                ESmin(k)=fstrSOLID%ESTRESS(MDOF*(ii-1)+k)
                IESmin(k)=j
              endif
            enddo
          endif
!       endif
      enddo
!!
      write(ILOG,*) '##### Local Summary @Node :Max/IdMax/Min/IdMin####'
      write(ILOG,*) '//U1 ',Umax(1),IUmax(1),Umin(1),IUmin(1)
      write(ILOG,*) '//U2 ',Umax(2),IUmax(2),Umin(2),IUmin(2)
      write(ILOG,*) '//U3 ',Umax(3),IUmax(3),Umin(3),IUmin(3)
      write(ILOG,*) '//E11',Emax(1),IEmax(1),Emin(1),IEmin(1)
      write(ILOG,*) '//E22',Emax(2),IEmax(2),Emin(2),IEmin(2)
      write(ILOG,*) '//E33',Emax(3),IEmax(3),Emin(3),IEmin(3)
      write(ILOG,*) '//E12',Emax(4),IEmax(4),Emin(4),IEmin(4)
      write(ILOG,*) '//E23',Emax(5),IEmax(5),Emin(5),IEmin(5)
      write(ILOG,*) '//E13',Emax(6),IEmax(6),Emin(6),IEmin(6)
      write(ILOG,*) '//E1 ',E_PRmax(1),IE_PRmax(1),E_PRmin(1),IE_PRmin(1)
      write(ILOG,*) '//E2 ',E_PRmax(2),IE_PRmax(2),E_PRmin(2),IE_PRmin(2)
      write(ILOG,*) '//E3 ',E_PRmax(3),IE_PRmax(3),E_PRmin(3),IE_PRmin(3)
      write(ILOG,*) '//S11',Smax(1),ISmax(1),Smin(1),ISmin(1)
      write(ILOG,*) '//S22',Smax(2),ISmax(2),Smin(2),ISmin(2)
      write(ILOG,*) '//S33',Smax(3),ISmax(3),Smin(3),ISmin(3)
      write(ILOG,*) '//S12',Smax(4),ISmax(4),Smin(4),ISmin(4)
      write(ILOG,*) '//S23',Smax(5),ISmax(5),Smin(5),ISmin(5)
      write(ILOG,*) '//S13',Smax(6),ISmax(6),Smin(6),ISmin(6)
      write(ILOG,*) '//SMS',Smax(7),ISmax(7),Smin(7),ISmin(7)
      write(ILOG,*) '//S1 ',S_PRmax(1),IS_PRmax(1),S_PRmin(1),IS_PRmin(1)
      write(ILOG,*) '//S2 ',S_PRmax(2),IS_PRmax(2),S_PRmin(2),IS_PRmin(2)
      write(ILOG,*) '//S3 ',S_PRmax(3),IS_PRmax(3),S_PRmin(3),IS_PRmin(3)
      write(ILOG,*) '##### Local Summary @Element :Max/IdMax/Min/IdMin####'
      write(ILOG,*) '//E11',EEmax(1),IEEmax(1),EEmin(1),IEEmin(1)
      write(ILOG,*) '//E22',EEmax(2),IEEmax(2),EEmin(2),IEEmin(2)
      write(ILOG,*) '//E33',EEmax(3),IEEmax(3),EEmin(3),IEEmin(3)
      write(ILOG,*) '//E12',EEmax(4),IEEmax(4),EEmin(4),IEEmin(4)
      write(ILOG,*) '//E23',EEmax(5),IEEmax(5),EEmin(5),IEEmin(5)
      write(ILOG,*) '//E13',EEmax(6),IEEmax(6),EEmin(6),IEEmin(6)
      write(ILOG,*) '//S11',ESmax(1),IESmax(1),ESmin(1),IESmin(1)
      write(ILOG,*) '//S22',ESmax(2),IESmax(2),ESmin(2),IESmin(2)
      write(ILOG,*) '//S33',ESmax(3),IESmax(3),ESmin(3),IESmin(3)
      write(ILOG,*) '//S12',ESmax(4),IESmax(4),ESmin(4),IESmin(4)
      write(ILOG,*) '//S23',ESmax(5),IESmax(5),ESmin(5),IESmin(5)
      write(ILOG,*) '//S13',ESmax(6),IESmax(6),ESmin(6),IESmin(6)
      write(ILOG,*) '//SMS',ESmax(7),IESmax(7),ESmin(7),IESmin(7)

!C*** Show Summary
!!---- oga2009 --    if (ds .ne. .true.) then !In case of Direct Solver prevent MPI
     if (ds .neqv. .true.) then !In case of Direct Solver prevent MPI
       call hecmw_allREDUCE_R(hecMESH,Umax,3,hecmw_max)
       call hecmw_allREDUCE_R(hecMESH,Umin,3,hecmw_min)
       call hecmw_allREDUCE_R(hecMESH,Emax,6,hecmw_max)
       call hecmw_allREDUCE_R(hecMESH,Emin,6,hecmw_min)
       call hecmw_allREDUCE_R(hecMESH,Smax,7,hecmw_max)
       call hecmw_allREDUCE_R(hecMESH,Smin,7,hecmw_min)
       call hecmw_allREDUCE_R(hecMESH,E_PRmax,3,hecmw_max)
       call hecmw_allREDUCE_R(hecMESH,E_PRmin,3,hecmw_min)
       call hecmw_allREDUCE_R(hecMESH,S_PRmax,3,hecmw_max)
       call hecmw_allREDUCE_R(hecMESH,S_PRmin,3,hecmw_min)
       call hecmw_allREDUCE_R(hecMESH,EEmax,6,hecmw_max)
       call hecmw_allREDUCE_R(hecMESH,EEmin,6,hecmw_min)
       call hecmw_allREDUCE_R(hecMESH,ESmax,7,hecmw_max)
       call hecmw_allREDUCE_R(hecMESH,ESmin,7,hecmw_min)
     end if 

     if( hecMESH%my_rank .eq. 0 ) then
        write(ILOG,*) '##### Global Summary @Node :Max/Min####'
        write(ILOG,*) '//U1 ',Umax(1),Umin(1)
        write(ILOG,*) '//U2 ',Umax(2),Umin(2)
        write(ILOG,*) '//U3 ',Umax(3),Umin(3)
        write(ILOG,*) '//E11',Emax(1),Emin(1)
        write(ILOG,*) '//E22',Emax(2),Emin(2)
        write(ILOG,*) '//E33',Emax(3),Emin(3)
        write(ILOG,*) '//E12',Emax(4),Emin(4)
        write(ILOG,*) '//E23',Emax(5),Emin(5)
        write(ILOG,*) '//E13',Emax(6),Emin(6)
        write(ILOG,*) '//E1 ',E_PRmax(1),E_PRmin(1)
        write(ILOG,*) '//E2 ',E_PRmax(2),E_PRmin(2)
        write(ILOG,*) '//E3 ',E_PRmax(3),E_PRmin(3)
        write(ILOG,*) '//S11',Smax(1),Smin(1)
        write(ILOG,*) '//S22',Smax(2),Smin(2)
        write(ILOG,*) '//S33',Smax(3),Smin(3)
        write(ILOG,*) '//S12',Smax(4),Smin(4)
        write(ILOG,*) '//S23',Smax(5),Smin(5)
        write(ILOG,*) '//S13',Smax(6),Smin(6)
        write(ILOG,*) '//SMS',Smax(7),Smin(7)
        write(ILOG,*) '//S1 ',S_PRmax(1),S_PRmin(1)
        write(ILOG,*) '//S2 ',S_PRmax(2),S_PRmin(2)
        write(ILOG,*) '//S3 ',S_PRmax(3),S_PRmin(3)
        write(ILOG,*) '##### Global Summary @Element :Max/Min####'
        write(ILOG,*) '//E11',EEmax(1),EEmin(1)
        write(ILOG,*) '//E22',EEmax(2),EEmin(2)
        write(ILOG,*) '//E33',EEmax(3),EEmin(3)
        write(ILOG,*) '//E12',EEmax(4),EEmin(4)
        write(ILOG,*) '//E23',EEmax(5),EEmin(5)
        write(ILOG,*) '//E13',EEmax(6),EEmin(6)
        write(ILOG,*) '//S11',ESmax(1),ESmin(1)
        write(ILOG,*) '//S22',ESmax(2),ESmin(2)
        write(ILOG,*) '//S33',ESmax(3),ESmin(3)
        write(ILOG,*) '//S12',ESmax(4),ESmin(4)
        write(ILOG,*) '//S23',ESmax(5),ESmin(5)
        write(ILOG,*) '//S13',ESmax(6),ESmin(6)
        write(ILOG,*) '//SMS',ESmax(7),ESmin(7)
     endif
!C
        if( IRESULT.eq.0 ) return
!C*** Write Result File
!C*** INITIALIZE
        header='*fstrresult'
        nd = i_step
        call hecmw_result_init(hecMESH,nd,header)
!C*** ADD
!C*** DISPLACEMENT
          if( fstrSOLID%iout_list(1) .eq. 1 ) then
            id = 1
            ndof=3
            label='DISPLACEMENT'
            call hecmw_result_add(id,ndof,label,fstrSOLID%DISP)
          end if
!C
        if( i_step .gt. 0) then
          if( fstrSOLID%iout_list(5) .eq. 1 ) then
!C*** STRAIN @node
            id = 1
            ndof=6
            label='STRAIN'
            call hecmw_result_add(id,ndof,label,fstrSOLID%STRAIN)

            id = 1
            ndof=12
            label='PRINCP_STRAIN'
            call hecmw_result_add(id,ndof,label,fstrSOLID%STRAIN_PRIS)
!C*** STRAIN @element
!!---  20100203   ----------------
            id = 2
            ndof=6
            label='ESTRAIN'
            call hecmw_result_add(id,ndof,label,fstrSOLID%ESTRAIN)
          end if
!
          if( fstrSOLID%iout_list(6) .eq. 1 ) then
!C*** STRESS @node
            id = 1
            ndof=7
            label='STRESS'
            call hecmw_result_add(id,ndof,label,fstrSOLID%STRESS)

            id = 1
            ndof=12
            label='PRINCP_STRESS'
            call hecmw_result_add(id,ndof,label,fstrSOLID%STRESS_PRIS)
!C*** STRESS @element
            id = 2
            ndof=7
            label='ESTRESS'
            call hecmw_result_add(id,ndof,label,fstrSOLID%ESTRESS)
          end if
        end if
!C*** WRITE NOW
        nameID='fstrRES'
        call hecmw_result_write_by_name(nameID)
!C*** FINALIZE
        call hecmw_result_finalize
!C
      return
      end subroutine fstrNLGEOM_Post

!C----------------------------------------------------------------------*
      subroutine fstrNLGEOM_make_result(hecMESH,fstrSOLID,fstrRESULT)
!C----------------------------------------------------------------------*
!     output of strain & stress in element
!!---  20100203  ---------------------

      use m_fstr
      implicit none
      type (hecmwST_local_mesh) :: hecMESH
      type (fstr_solid)         :: fstrSOLID
      type (hecmwST_result_data):: fstrRESULT
!** Local variables
      integer(kind=kint) NDOF,i,k  
      integer(kind=kint) ielem,ID_area   

!C***
      NDOF = hecMESH%n_dof
      if( NDOF.eq.3  .or. NDOF.eq.2 ) then
        fstrRESULT%nn_component=5
        fstrRESULT%ne_component=2
!C @node
        allocate( fstrRESULT%nn_dof(5) )
        allocate( fstrRESULT%node_label(5) )
        allocate( fstrRESULT%node_val_item(40*hecMESH%n_node))
        fstrRESULT%nn_dof(1)=3
        fstrRESULT%nn_dof(2)=6
        fstrRESULT%nn_dof(3)=7
        fstrRESULT%nn_dof(4)=12
        fstrRESULT%nn_dof(5)=12
        fstrRESULT%node_label(1)='DISPLACEMENT'
        fstrRESULT%node_label(2)='STRAIN'
        fstrRESULT%node_label(3)='STRESS'
        fstrRESULT%node_label(4)='PRINCP_STRAIN'
        fstrRESULT%node_label(5)='PRINCP_STRESS'
!C @element
        allocate( fstrRESULT%ne_dof(2) )
        allocate( fstrRESULT%elem_label(2) )
        allocate( fstrRESULT%elem_val_item(13*hecMESH%n_elem))
        fstrRESULT%ne_dof(1)=6
        fstrRESULT%ne_dof(2)=7
        fstrRESULT%elem_label(1)='ESTRAIN'
        fstrRESULT%elem_label(2)='ESTRESS'
        fstrRESULT%elem_val_item=0.0
!C*** Set Displacement/Strain/Stress
!C @node
        do i= 1, hecMESH%n_node
          do k=1,3
            fstrRESULT%node_val_item(40*(i-1)+k)    &
            =fstrSOLID%DISP(3*i-3+k)
          enddo
          do k=1,6
            fstrRESULT%node_val_item(40*(i-1)+k+3)  &
            =fstrSOLID%STRAIN(6*i-6+k)
          enddo
          do k=1,7
            fstrRESULT%node_val_item(40*(i-1)+k+9)  &
            =fstrSOLID%STRESS(7*i-7+k)
          enddo
          do k=1,12
            fstrRESULT%node_val_item(40*(i-1)+k+16)  &
            =fstrSOLID%STRAIN_PRIS(12*i-12+k)
          enddo
          do k=1,12
            fstrRESULT%node_val_item(40*(i-1)+k+28)  &
            =fstrSOLID%STRESS_PRIS(12*i-12+k)
          enddo
        enddo
!C @element
        do i= 1, hecMESH%n_elem
!!!       ielem   = hecMESH%elem_ID( i*2-1 )
          ielem   = i
          ID_area = hecMESH%elem_ID( i*2 )
!!!          if( ID_area /= hecMESH%my_rank ) cycle
          do k=1,6
            fstrRESULT%elem_val_item(13*(i-1)+k)     &
            =fstrSOLID%ESTRAIN(6*ielem-6+k)
          enddo
          do k=1,7
            fstrRESULT%elem_val_item(13*(i-1)+k+6)   &
            =fstrSOLID%ESTRESS(7*ielem-7+k)
          enddo
        enddo
      end if
!
      return
      end subroutine fstrNLGEOM_make_result 
!C----------------------------------------------------------------------*
!----------------------------------------------------------------------*
      subroutine fstrNLGEOM_init_gbar( hecMESH,hecMAT,fstrSOLID) 
!----------------------------------------------------------------------*
!
!C  Set covariant base vector & mechanical properties in element of material coord. 
!
      use m_fstr
      implicit none
      type (hecmwST_local_mesh) :: hecMESH
      type (hecmwST_matrix    ) :: hecMAT
      type (fstr_solid       )  :: fstrSOLID
!C
!C Local variables for 2D and 3D elements
!C
      integer(kind=kint) :: ndof, count, itype, iS, iE, ic_type, nn, icel
      integer(kind=kint) :: iiS, i, j, k, idx, idx_e, isect, id, iset
      integer(kind=kint) :: inod, inod_e
      integer(kind=kint) :: iter, i_step 
      real(kind=kreal) :: xx(20), yy(20), zz(20), ee, pp, local_stf(1830), stiffness(20*6, 20*6)
      real(kind=kreal) :: disp(60), qf(60), tt(20), tti(20), CUR_nod(20), CUR_nod_I(20)
      integer(kind=kint) :: nodLOCAL(20)
      integer(kind=kint) :: g1, g2, g3, mid 

      real(kind=kreal) :: rho, alpha, thick, param1 
      real(kind=kreal) :: DT0, STRAIN_V1(6,3,3,3),STRAIN_V2(6,3,3,3)
!C
!C +-------+
!C | INIT. |
!C +-------+
!C====
!C
!C +-------------------------------+
!C | ELEMENT-by-ELEMENT ASSEMBLING | 
!C | according to ELEMENT TYPE     |
!C +-------------------------------+
      ndof = hecMAT%NDOF

      count=0
      do itype= 1, hecMESH%n_elem_type
        iS= hecMESH%elem_type_index(itype-1) + 1
        iE= hecMESH%elem_type_index(itype  )
        ic_type= hecMESH%elem_type_item(itype)
!C** Ignore link elements
        if (hecmw_is_etype_link(ic_type)) cycle
!C** Set number of nodes	
        nn = hecmw_get_max_node(ic_type)
!C element loop
        do icel= iS, iE
          count=count+1
!C** node ID
          iiS= hecMESH%elem_node_index(icel-1)
!C** section  ID
          isect= hecMESH%section_ID(icel)
          mid=hecMESH%section%sect_mat_ID_item(isect) 
!C** Get Properties
!         call fstrNLGEOM_get_prop(hecMESH,isect,ee,pp,rho,alpha,thick)
!!!!!   part material  !!!!!
     if(count .eq. -200) then
      if(count .eq. 0) then
        write(ILOG,*) 'check data '
      end if
        write(ILOG,*) 'elem = ',count,  ' mid = ',mid
        call flush(ILOG)
     end if
!!!!!   part material  !!!!!

          do i = 1, 4
           if ( mid .eq. i) then
             fstrSOLID%mat_iso  (count) =  fstrSOLID%mat_iso_num( i)
              if ( fstrSOLID%mat_iso(count) .eq. -1) then
                fstrSOLID%mat_iso  (count) =  1 
              end if
           end if
          end do
!
        enddo
      enddo
!
      return
      end subroutine fstrNLGEOM_init_gbar 
!----------------------------------------------------------------------*
end module m_fstr_solve_NLGEOM
