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

module m_heat_solve_TRAN


!!    private :: heat_output_monit

   contains
!C
!C** TRANSIENT
!C
   subroutine heat_solve_TRAN ( hecMESH,hecMAT,fstrRESULT,fstrPARAM,fstrHEAT,ISTEP,CTIME )

      use m_fstr
      use m_heat_mat_ass_conductivity
      use m_heat_mat_ass_capacity
      use m_heat_mat_ass_boundary
      use m_heat_init
      use m_heat_make_result
      use m_heat_output
      use m_hecmw2fstr_mesh_conv

      implicit none
      integer(kind=kint) ISTEP,ITM,incr,iend,iterALL,iterCUR,i,inod,mnod,nd,id,ndof,bup_n_dof
      integer(kind=kint) g1,g2,g3 
      real(kind=kreal)   CTIME,ST,DTIME,EETIME,DELMAX,DELMIN,TT,BETA,VAL,CHK,tmpmax,dltmp,tmpmax_myrank
!C file name
      character(len=HECMW_HEADER_LEN) :: header
      character(len=HECMW_NAME_LEN) :: label
      character(len=HECMW_NAME_LEN) :: nameID

      type (hecmwST_local_mesh ) :: hecMESH
      type (hecmwST_matrix     ) :: hecMAT
      type (hecmwST_result_data) :: fstrRESULT
      type (fstr_param         ) :: fstrPARAM
      type (fstr_heat          ) :: fstrHEAT

!!----- oga2009 -------------------------------------------------

      integer(kind=kint),parameter :: heat_IR1        =   86
      integer(kind=kint),parameter :: heat_IW1        =   87
      integer(kind=kint),parameter :: heat_IW4        =   88
      integer(kind=kint),parameter :: heat_IW5        =   89

      character(len=30) :: cdummy

      integer(kind=kint) :: j, ierror, stat, ios2 
      integer(kind=kint) :: i_flg_1,i_flg_2,icont2
      integer(kind=kint) :: my_rank_monit_1, my_rank_monit_2
      integer(kind=kint) :: n_gauss, ic_type
      integer(kind=kint) :: restrt_step_num, nodes_all_area 
      integer(kind=kint),pointer :: restrt_step(:)
      real(kind=kreal)  ,pointer :: restrt_time(:)
      real(kind=kreal)   :: CURE3D_0, xx(12)

!!----- oga2009 -------------------------------------------------

      if( ISTEP .eq. 1 ) then
        ST = 0.0d0
      else
        ST = 0.0d0
        do i = 1, ISTEP - 1
          ST = ST + fstrHEAT%STEP_EETIME(i)
        enddo
      endif   

      DTIME  = fstrHEAT%STEP_DLTIME(ISTEP)
      EETIME = fstrHEAT%STEP_EETIME(ISTEP)
      DELMIN = fstrHEAT%STEP_DELMIN(ISTEP)
      DELMAX = fstrHEAT%STEP_DELMAX(ISTEP)
      ITM    = fstrPARAM%itmax(ISTEP)
      EPS    = fstrPARAM%eps(ISTEP)
      TT     = CTIME
      if( NPRINT.eq.0 ) NPRINT = 1
      write(*,*) ' DTIME=',DTIME
      write(*,*) 'EETIME=',EETIME
      write(*,*) 'DELMIN=',DELMIN
      write(*,*) 'DELMAX=',DELMAX
      write(*,*) '   ITM=',ITM
      write(*,*) '   EPS=',EPS
      write(*,*) '    TT=',TT
      write(*,*) 'NPRINT=',NPRINT

      BETA = 0.5
      incr = 0
      iend = 0
      hecMAT%NDOF = 1
      hecMAT%Iarray(98) = 1   !Assmebly complete

!!----- oga2009 -------------------------------------------------

!!--- read heat 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. fstrHEAT%node_monit_1 ) then
          fstrHEAT%node_monit_1 = i
          i_flg_1 = 1
          do j= 1, hecMESH%nn_internal
            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
        if( hecMESH%global_elem_ID(i) .eq. fstrHEAT%elem_monit_1 ) then
          fstrHEAT%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( fstrHEAT%restart_nout .eq. 0 ) then
        fstrHEAT%restart_nout = 999999999
      end if
!
!C
!C-- file open for local use
!C
      if( hecMESH%my_rank .eq. my_rank_monit_1) then
        OPEN(heat_IW4,FILE='heat_temp_p1.out', status = 'replace', iostat=stat)
        if( stat /= 0 ) then
          write(*,*) 'stop due to file opening error <heat_temp_p1.>'
          call hecmw_abort( hecmw_comm_get_comm())
        end if
      endif
      if( hecMESH%my_rank .eq. my_rank_monit_1) then
        OPEN(heat_IW5,FILE='heat_cure_p1.out', status = 'replace', iostat=stat)
        if( stat /= 0 ) then
          write(*,*) 'stop due to file opening error <heat_cure_p1.>'
          call hecmw_abort( hecmw_comm_get_comm())
        end if
      endif
!
!C
!C Allocate Matrix and Vector
!C
      allocate (fstrHEAT%mat_iso(hecMESH%n_elem)          ,STAT=ierror )
          if( ierror /= 0 ) then
            write(idbg,*) 'stop due to allocation error <heat_solve_TRAN, mat_iso>'
            write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
            call flush(idbg)
            call hecmw_abort( hecmw_comm_get_comm())
          end if
      allocate (fstrHEAT%idx_mat_node(hecMESH%n_node)          ,STAT=ierror )
          if( ierror /= 0 ) then
            write(idbg,*) 'stop due to allocation error <heat_solve_TRAN, idx_mat_node>'
            write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
            call flush(idbg)
            call hecmw_abort( hecmw_comm_get_comm())
          end if
      allocate (fstrHEAT%CURE3D_nod(hecMESH%n_node)          ,STAT=ierror )
          if( ierror /= 0 ) then
            write(idbg,*) 'stop due to allocation error <heat_solve_TRAN, CURE3D_nod>'
            write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
            call flush(idbg)
            call hecmw_abort( hecmw_comm_get_comm())
          end if
      allocate (fstrHEAT%CURE3D_nod_I(hecMESH%n_node)          ,STAT=ierror )
          if( ierror /= 0 ) then
            write(idbg,*) 'stop due to allocation error <heat_solve_TRAN, CURE3D_nod_I>'
            write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
            call flush(idbg)
            call hecmw_abort( hecmw_comm_get_comm())
          end if
      allocate (fstrHEAT%CURE3D_nod_cfrp(hecMESH%n_node)          ,STAT=ierror )
          if( ierror /= 0 ) then
            write(idbg,*) 'stop due to allocation error <heat_solve_TRAN, CURE3D_nod_cfrp>'
            write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
            call flush(idbg)
            call hecmw_abort( hecmw_comm_get_comm())
          end if
      allocate (fstrHEAT%D_CURE3D_nod(hecMESH%n_node)          ,STAT=ierror )
          if( ierror /= 0 ) then
            write(idbg,*) 'stop due to allocation error <heat_solve_TRAN, D_CURE3D_nod>'
            write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
            call flush(idbg)
            call hecmw_abort( hecmw_comm_get_comm())
          end if
      allocate (fstrHEAT%GBARU(3,3,hecMESH%n_node)          ,STAT=ierror )
          if( ierror /= 0 ) then
            write(idbg,*) 'stop due to allocation error <heat_solve_TRAN, GBARU>'
            write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
            call flush(idbg)
            call hecmw_abort( hecmw_comm_get_comm())
          end if
!
      fstrHEAT%mat_iso         = 0
      fstrHEAT%idx_mat_node    = 0
      fstrHEAT%CURE3D_nod      = fstrHEAT%ccc_rl( 9)
      fstrHEAT%CURE3D_nod_I    = fstrHEAT%ccc_rl( 9)
      fstrHEAT%CURE3D_nod_cfrp = fstrHEAT%ccc_rl( 9)
      fstrHEAT%D_CURE3D_nod    = 0.0D+00
      fstrHEAT%GBARU           = 0.0D+00

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

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

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

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

      if( incr .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 (fstrHEAT%GBARU_ALL(3,3,nodes_all_area)          ,STAT=ierror )
          if( ierror /= 0 ) then
            write(idbg,*) 'stop due to allocation error <heat_solve_TRAN, GBARU_ALL>'
            write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
            call flush(idbg)
            call hecmw_abort( hecmw_comm_get_comm())
          end if
        fstrHEAT%GBARU_ALL = 0.0D+00
        xx                 = 0.0D+00

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

!!-- 20100330 ----!
          call heat_init_gbar ( hecMESH,hecMAT,fstrHEAT )
!!-- 20100330 ----!
        
!!! restart  !!!
      allocate (restrt_step(1)           ,STAT=ierror )
          if( ierror /= 0 ) then
            write(idbg,*) 'stop due to allocation error <heat_solve_TRAN, restrt_step >'
            write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
            call flush(idbg)
            call hecmw_abort( hecmw_comm_get_comm())
          end if
      allocate (restrt_time(1)           ,STAT=ierror )
          if( ierror /= 0 ) then
            write(idbg,*) 'stop due to allocation error <heat_solve_TRAN, restrt_time >'
            write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
            call flush(idbg)
            call hecmw_abort( hecmw_comm_get_comm())
          end if

!!----- oga2009 -------------------------------------------------

!---------------------------------------------
!!! restart  !!!
      restrt_step(1) = 1
      restrt_time(1) = 0.0D+00
      restrt_step_num = 1

    if(fstrHEAT%restart_nout .lt. 0 ) then
       n_gauss=2
          ic_type= hecMESH%elem_type_item(1)
          if (ic_type.eq.342.or.ic_type.eq.352.or.ic_type.eq.362) then
            n_gauss=3
          end if
          if (ic_type.eq.242 .or. ic_type.eq.232 ) then
            n_gauss=3
          end if
      call hecmw_restart_open()
      call hecmw_restart_read_int(restrt_step)
      call hecmw_restart_read_real(restrt_time)
      call hecmw_restart_read_real(fstrHEAT%  TEMP)
      call hecmw_restart_read_real(fstrHEAT%D_TEMP)
      call hecmw_restart_read_real(fstrHEAT%  CURE3D_nod)
      call hecmw_restart_read_real(fstrHEAT%D_CURE3D_nod)
      call hecmw_restart_close()
      restrt_step_num = restrt_step(1) + 1
      incr            = restrt_step(1) 
      TT              = restrt_time(1) 
      CTIME = TT
      do i= 1, hecMESH%n_node
        fstrHEAT%CURE3D_nod_I(i) = fstrHEAT%CURE3D_nod(i)
      enddo
      do i= 1, hecMESH%n_node
        fstrHEAT%TEMP0(i) = fstrHEAT%TEMP(i)
      enddo
    end if

!------ INITIAL OUTPUT
!
    if( fstrHEAT%idx_cure .eq. 1) then
!!   if( hecMESH%my_rank .eq. 0) then
          call hecmw_update_1_R ( hecMESH,fstrHEAT%TEMP       ,hecMESH%n_node )
          call hecmw_update_1_R ( hecMESH,fstrHEAT%TEMP0      ,hecMESH%n_node )
          call hecmw_update_1_R ( hecMESH,fstrHEAT%TEMPC      ,hecMESH%n_node )
          call hecmw_update_1_R ( hecMESH,fstrHEAT%D_TEMP     ,hecMESH%n_node )
          call hecmw_update_1_R ( hecMESH,fstrHEAT%CURE3D_nod  ,hecMESH%n_node )
          call hecmw_update_1_R ( hecMESH,fstrHEAT%CURE3D_nod_I,hecMESH%n_node )
          call hecmw_update_1_R ( hecMESH,fstrHEAT%D_CURE3D_nod,hecMESH%n_node )

          write(heat_sta_IW1) 'step'
          write(heat_sta_IW2) 'step'
          write(heat_sta_IW1) ' INCR  =',incr,' Time  =',CTIME
          write(heat_sta_IW2) ' INCR  =',incr,' Time  =',CTIME
          call heat_output_sta(hecMESH,hecMAT,fstrHEAT,fstrPARAM,heat_sta_IW1   &
                              ,heat_sta_IW2)
!!   end if
    end if

      if( restrt_step_num .eq. 1 ) then
          write(ILOG,*)
          write(ILOG,'(a,i6)')              ' ISTEP =',ISTEP
          write(ILOG,'(a,i10,a,1pe13.4e3)') ' INCR  =',incr ,' Time  =',CTIME
          call heat_output(hecMESH,hecMAT,fstrHEAT,fstrPARAM)
      end if
!!

!C--------------------   START TRANSIET LOOP   ------------------------
      tr_loop: do
!C--------------------   
        incr = incr + 1
        if( TT+DTIME.ge.EETIME ) then
          DTIME = EETIME - TT
          TT = EETIME
          iend = 1
        else
          TT = TT + DTIME
          iend = 0
        endif
        CTIME = TT

       if( mod(incr,fstrHEAT%nout_monit) == 0 ) then
        if( hecMESH%my_rank.eq.0 ) then
            write(IMSG,*) 
            write(IMSG,*) '// INCREMENT NO. =',incr
            write(IMSG,*) '//  CURRENT TIME =',TT
            write(IMSG,*) '//    DELTA TIME =',DTIME
        endif
       end if

        if( DTIME .lt.DELMIN ) then
          if( hecMESH%my_rank.eq.0 ) then
            write(IMSG,*) ' !!! DELTA TIME EXCEEDED TOLERANCE OF TIME INCREMENT'
            call flush(IMSG)
          endif
          call hecmw_abort( hecmw_comm_get_comm() )
        endif

        if( incr .gt. INCMAX ) then
          if( hecMESH%my_rank.eq.0 ) then
            write(IMSG,*) ' !!! NUMBER OF INCREMENTS EXCEEDED MAXIMUM INCREMENTS'
            call flush(IMSG)
          endif
          call hecmw_abort( hecmw_comm_get_comm() )
        endif
!!
!C==============  START OF ITERATION LOOP  ===========
        iterALL= 0
        iterCUR= 0

    if( fstrHEAT%idx_cure .eq. 1) then

!---- CURE BALANCE         at time = t --------- 
        if( mod(incr,fstrHEAT%nout_monit) == 0 ) then
         if( hecMESH%my_rank.eq.0 ) then
           write(IMSG,*) ' -- ITERATION FOR CURE -- '
           call flush(IMSG)
         end if
        end if

        BETA = 0.0
!
        do iterCUR=1,ITM

!C-- MATRIX ASSEMBLING FOR CURE  -----

          hecMAT%AL= 0.d0
          hecMAT%AU= 0.d0
          hecMAT%D = 0.d0
          hecMAT%B = 0.d0
!
          call cure_mat_ass_capacity ( hecMESH,hecMAT,fstrHEAT,DTIME )
 
          call cure_mat_ass_boundary ( hecMESH,hecMAT,fstrHEAT,TT,ST,DTIME,BETA,incr,iterCUR )
!C
!C-- SOLVER FOR CURE    -----
!!!  20110623  !!
          if( incr .eq. restrt_step_num .and. iterCUR .eq. 1 ) then
            hecMAT%Iarray(98) = 1   !Assmebly complete
          end if
          hecMAT%Iarray(97) = 1   !Need numerical factorization

          if( hecMAT%Iarray(99).eq. 1 ) then
            call hecmw_solve_11 ( hecMESH,hecMAT )
          else
            bup_n_dof = hecMESH%n_dof
            hecMESH%n_dof = 1

            if(hecMESH%PETOT.gt.1) then
              call hecmw_solve_direct_parallel(hecMESH,hecMAT,IMSG)
            else
              call hecmw_solve_direct(hecMESH,hecMAT,IMSG)
            end if
!!           hecMAT%X = hecMAT%B -- leading stack overflow (intel9)
             do i=1,hecMAT%NP*hecMESH%n_dof
                 hecMAT%X(i) = hecMAT%B(i)
             end do
            hecMESH%n_dof=bup_n_dof
          endif
!C
!C-- UPDATE FOR CURE   -----

          do i= 1, hecMESH%n_node
            fstrHEAT%TEMPC(i)= fstrHEAT%CURE3D_nod(i)
          enddo
    
          call hecmw_update_1_R ( hecMESH,hecMAT%X,hecMAT%NP )

          icont2 = 0
          do i= 1, hecMESH%n_node
            fstrHEAT%CURE3D_nod (i)= hecMAT%X(i)
          enddo
!C
!C-- GLOBAL RESIDUAL FOR CURE     -----

          VAL= 0.d0 
          do i= 1, hecMESH%nn_internal
            VAL= VAL + (fstrHEAT%CURE3D_nod(i)-fstrHEAT%TEMPC(i))**2
          enddo
          call hecmw_allREDUCE_R1 ( hecMESH,VAL,hecmw_sum )
!C
          CHK= dsqrt(VAL)

          if ( hecMESH%my_rank.eq.0 ) then 
            write(*,'(i8,1pe16.6)') iterCUR, CHK
            if( mod(incr,fstrHEAT%nout_monit) == 0 ) then
              write(IMSG,'(i8,1pe16.6)') iterCUR, CHK
              call flush(IMSG)
            end if
          endif
!C
!!!!  20111014  -------
!!        if ( CHK.lt.EPS ) then
          if ( CHK.lt. 1.0D-06 ) then
            if ( hecMESH%my_rank.eq.0 ) then 
              write(*,*) '  '
              write(*,*) ' !!! CURE CONVERGENCE  '
              if( mod(incr,fstrHEAT%nout_monit) == 0 ) then
               write(IMSG,*) '  '
               write(IMSG,*) ' !!! CURE CONVERGENCE  '
               call flush(IMSG)
              end if
            endif
            exit
          endif
!C
          if ( iterCUR.ge.ITM ) then
            if ( hecMESH%my_rank.eq.0 ) then 
              write(*,*) ' !!! ITERATION COUNT OVER FOR CURE : MAX = ', ITM
              write(IMSG,*) ' !!! ITERATION COUNT OVER FOR CURE : MAX = ', ITM
              call flush(IMSG)
            endif
            call hecmw_abort( hecmw_comm_get_comm() )
          endif
!
        end do
    end if

!---- TEMPERATURE ENERGY BALANCE      at time = t+dt --------- 
!
        if( mod(incr,fstrHEAT%nout_monit) == 0 ) then
         if( hecMESH%my_rank.eq.0 ) then
           write(IMSG,*) ' -- ITERATION FOR TEMPERATURE -- '
           call flush(IMSG)
         end if
        end if
!
        BETA = 1.0
!
        do
!C============== 
          iterALL= iterALL + 1
!C
!C-- MATRIX ASSEMBLING FOR TEMPERATURE   -----

          call heat_mat_ass_conductivity ( hecMESH,hecMAT,fstrHEAT,BETA )
 
          call heat_mat_ass_capacity ( hecMESH,hecMAT,fstrHEAT,DTIME )
 
          call heat_mat_ass_boundary ( hecMESH,hecMAT,fstrHEAT,TT,ST,DTIME,BETA,incr,iterALL )
!C
!C-- SOLVER FOR TEMPERATURE   -----
!!!  20110623  !!
          if( incr .eq. restrt_step_num .and. iterALL .eq. 1 ) then
            hecMAT%Iarray(98) = 1   !Assmebly complete
          end if
          hecMAT%Iarray(97) = 1   !Need numerical factorization

          if( hecMAT%Iarray(99).eq. 1 ) then
            call hecmw_solve_11 ( hecMESH,hecMAT )
          else
            bup_n_dof = hecMESH%n_dof
            hecMESH%n_dof = 1

            if(hecMESH%PETOT.gt.1) then
              call hecmw_solve_direct_parallel(hecMESH,hecMAT,IMSG)
            else
              call hecmw_solve_direct(hecMESH,hecMAT,IMSG)
            end if
!!           hecMAT%X = hecMAT%B -- leading stack overflow (intel9)
             do i=1,hecMAT%NP*hecMESH%n_dof
                 hecMAT%X(i) = hecMAT%B(i)
             end do
            hecMESH%n_dof=bup_n_dof
          endif
!C
!C-- UPDATE  FOR TEMPERATURE   -----

          do i= 1, hecMESH%n_node
            fstrHEAT%TEMPC(i)= fstrHEAT%TEMP(i)
          enddo
    
          call hecmw_update_1_R ( hecMESH,hecMAT%X,hecMAT%NP )

          do i= 1, hecMESH%n_node
            fstrHEAT%TEMP (i)= hecMAT%X(i)
          enddo
!C
!C-- GLOBAL RESIDUAL FOR TEMPERATURE    -----

          VAL= 0.d0 
          do i= 1, hecMESH%nn_internal
            VAL= VAL + (fstrHEAT%TEMP(i)-fstrHEAT%TEMPC(i))**2
          enddo
          call hecmw_allREDUCE_R1 ( hecMESH,VAL,hecmw_sum )
!C
          CHK= dsqrt(VAL)

          if ( hecMESH%my_rank.eq.0 ) then 
            write(*,'(i8,1pe16.6)') iterALL, CHK
            if( mod(incr,fstrHEAT%nout_monit) == 0 ) then
              write(IMSG,'(i8,1pe16.6)') iterALL, CHK
              call flush(IMSG)
            end if
          endif
!C
          if ( CHK.lt.EPS ) then
            if ( hecMESH%my_rank.eq.0 ) then 
              write(*,*) '  '
              write(*,*) ' !!! TEMPERATURE CONVERGENCE '
              if( mod(incr,fstrHEAT%nout_monit) == 0 ) then
               write(IMSG,*) '  '
               write(IMSG,*) ' !!! TEMPERATURE CONVERGENCE '
               call flush(IMSG)
              end if
            endif
            exit
          endif
!C
          if ( iterALL.ge.ITM ) then
            if ( hecMESH%my_rank.eq.0 ) then 
              write(*,*) ' !!! ITERATION COUNT OVER FOR TEMPERATURE : MAX = ', ITM
              write(IMSG,*) ' !!! ITERATION COUNT OVER FOR TEMPERATURE : MAX = ', ITM
              call flush(IMSG)
            endif

            if( DELMIN .gt. 0.d0) then
              TT = TT - DTIME
              DTIME = 0.5*DTIME
              cycle tr_loop
            else
              call hecmw_abort( hecmw_comm_get_comm() )
            endif
          endif
!C============== 
        enddo
!C==============  END OF ITERATION LOOP  ===========
!C

        if( DELMIN .gt. 0.d0 ) then
          tmpmax = 0.d0
          do i= 1, hecMESH%nn_internal
            inod=fstrPARAM%global_local_id(1,i)
            dltmp  = fstrHEAT%TEMP0(i) - fstrHEAT%TEMP(i)
            if( DABS(dltmp).gt.tmpmax ) then
              mnod = inod
              tmpmax = DABS(dltmp)
            endif
          enddo

          tmpmax_myrank = tmpmax
          call hecmw_allREDUCE_R1(hecMESH, tmpmax, hecmw_max)

          if (tmpmax .gt. tmpmax_myrank) mnod = -1
          call hecmw_allREDUCE_I1(hecMESH, mnod, hecmw_max)

          if ( tmpmax .gt. DELMAX ) then
            if( hecMESH%my_rank.eq.0 ) then
              write(IMSG,*) ' !!! EXCEEDED TOLERANCE OF VARIATION IN TEMPERATUTE.'
              write(IMSG,*) ' : NODE NUMBER  = ', mnod
              write(IMSG,*) ' : DELTA TEMP   = ', tmpmax
              call flush(IMSG)
            endif
            TT    = TT - DTIME
            DTIME = 0.5*DTIME
            cycle 
          endif

          if ( iterALL .le. 2 ) DTIME = DTIME * 2.0
        endif

!--- update cure & temperature   ----------------------------
!--- update cure                   at time = t+dt ---
        do i= 1, hecMESH%n_node
          fstrHEAT%D_CURE3D_nod(i) = fstrHEAT%CURE3D_nod(i) - fstrHEAT%CURE3D_nod_I(i) 
          fstrHEAT%CURE3D_nod_I(i) = fstrHEAT%CURE3D_nod(i)
        enddo
!--- update temperature            at time = t+dt  ---
        do i= 1, hecMESH%n_node
          fstrHEAT%D_TEMP(i) = fstrHEAT%TEMP(i) - fstrHEAT%TEMP0(i) 
          fstrHEAT% TEMP0(i) = fstrHEAT%TEMP(i)
        enddo
!------------------------------------------------------------

!!! restart  !!!
    if(fstrHEAT%restart_nout .lt. 0) then
      fstrHEAT%restart_nout = - fstrHEAT%restart_nout
    end if
    if( mod(incr,fstrHEAT%restart_nout) == 0 ) then
      n_gauss=2
          ic_type= hecMESH%elem_type_item(1)
          if (ic_type.eq.342.or.ic_type.eq.352.or.ic_type.eq.362) then
            n_gauss=3
          end if
          if (ic_type.eq.242 .or. ic_type.eq.232 ) then
            n_gauss=3
          end if

      restrt_step(1) = incr 
      restrt_time(1) = CTIME 
      call hecmw_restart_add_int(restrt_step,size(restrt_step))
      call hecmw_restart_add_real(restrt_time,size(restrt_time))
      call hecmw_restart_add_real(fstrHEAT%  TEMP,size(fstrHEAT%  TEMP))
      call hecmw_restart_add_real(fstrHEAT%D_TEMP,size(fstrHEAT%D_TEMP))
      call hecmw_restart_add_real(fstrHEAT%  CURE3D_nod,size(fstrHEAT%  CURE3D_nod))
      call hecmw_restart_add_real(fstrHEAT%D_CURE3D_nod,size(fstrHEAT%D_CURE3D_nod))
      call hecmw_restart_write()
    end if

!C
!C=== OUTPUT
!C
        if( mod(incr,fstrHEAT%nout_sta) == 0 ) then

         if( fstrHEAT%idx_cure .eq. 1) then
!!        if( hecMESH%my_rank .eq. 0) then
          call hecmw_update_1_R ( hecMESH,fstrHEAT%TEMP      ,hecMESH%n_node )
          call hecmw_update_1_R ( hecMESH,fstrHEAT%CURE3D_nod,hecMESH%n_node )

          write(heat_sta_IW1) 'step'
          write(heat_sta_IW2) 'step'
          write(heat_sta_IW1) ' INCR  =',incr,' Time  =',CTIME
          write(heat_sta_IW2) ' INCR  =',incr,' Time  =',CTIME

           call heat_output_sta(hecMESH,hecMAT,fstrHEAT,fstrPARAM,heat_sta_IW1   &
                              ,heat_sta_IW2)
!!        end if
         end if
        end if

        if( mod(incr,fstrHEAT%nout) == 0 .or. iend == 1) then    !--- oga2009 --

          write(ILOG,*)
          write(ILOG,'(a,i6)')              ' ISTEP =',ISTEP
          write(ILOG,'(a,i10,a,1pe13.4e3)') ' INCR  =',incr ,' Time  =',CTIME

          call heat_output(hecMESH,hecMAT,fstrHEAT,fstrPARAM)

!!---  Post by visualizer  

!       for CFRP
         if( fstrHEAT%idx_cure .eq. 1) then
           call cure_output_cfrp(hecMESH,hecMAT,fstrHEAT) 
         end if
!!
          call heat_post(hecMESH,hecMAT,fstrHEAT)

!!---  Write Result File

          if( IRESULT.eq.1 ) then
            nd = incr
            header = '*fstrresult'
            call hecmw_result_init( hecMESH,nd,header )
            if(fstrHEAT%iout_list(1) .eq. 1) then
              id    = 1
              ndof  = 1
              label = 'TEMPERATURE'
              call hecmw_result_add(id,ndof,label,fstrHEAT%TEMP)
            end if
            if(fstrHEAT%iout_list(2) .eq. 1) then
              id    = 1
              ndof  = 1
              label = 'CURE'
              call hecmw_result_add(id,ndof,label,fstrHEAT%CURE3D_nod_cfrp)
            end if

            nameID = 'fstrRES'
            call hecmw_result_write_by_name(nameID)
            call hecmw_result_finalize
            if( hecMESH%my_rank.eq.0 ) then
              write(IMSG,*) '### FSTR output Result_File.'
              call flush(IMSG)
            endif
          endif

!!-- POST PROCESSING VIA MEMORY

          if( IVISUAL.eq.1 ) then
            call heat_init_result ( hecMESH, fstrRESULT )
            call heat_make_result ( hecMESH, fstrHEAT, fstrRESULT )
              call fstr2hecmw_mesh_conv(hecMESH)
            call hecmw_visualize_init
            idummy = 1
            call hecmw_visualize ( hecMESH,fstrRESULT,incr,INCMAX,idummy )
            call hecmw_visualize_finalize
              call hecmw2fstr_mesh_conv(hecMESH)
            call hecmw_result_free(fstrRESULT)
            if( hecMESH%my_rank.eq.0 ) then
              write(IMSG,*) '### FSTR output Visual_File.'
              call flush(IMSG)
            endif
          endif

        endif

!C-- output result of monitoring node    -------  oga2009 --
      if( mod(incr,fstrHEAT%nout_monit) == 0 ) then
         call heat_output_monit   (hecMESH, fstrPARAM, fstrHEAT, heat_IW4, heat_IW5, &
                                   incr, CTIME, my_rank_monit_1, my_rank_monit_2)
      end if
!C===
      if( iend.ne.0 ) exit 
!C--------------------   
      enddo tr_loop
!C--------------------   START TRANSIET LOOP   ------------------------

!!----- oga2009 -------------------------------------------------
      deallocate(fstrHEAT%mat_iso              ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <heat_solve_TRAN, mat_iso>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      deallocate(fstrHEAT%idx_mat_node               ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <heat_solve_TRAN, idx_mat_node>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      deallocate(fstrHEAT%CURE3D_nod              ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <heat_solve_TRAN, CURE3D_nod>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      deallocate(fstrHEAT%CURE3D_nod_I              ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <heat_solve_TRAN, CURE3D_nod_I>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      deallocate(fstrHEAT%CURE3D_nod_cfrp             ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <heat_solve_TRAN, CURE3D_nod_cfrp>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      deallocate(fstrHEAT%D_CURE3D_nod              ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <heat_solve_TRAN, D_CURE3D_nod>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      deallocate(fstrHEAT%GBARU             ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <heat_solve_TRAN, GBARU>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      deallocate(fstrHEAT%GBARU_ALL         ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <heat_solve_TRAN, GBARU_ALL>'
              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 <heat_solve_TRAN, restrt_step>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if
      deallocate(restrt_time              ,STAT=ierror)
            if( ierror /= 0 ) then
              write(idbg,*) 'stop due to deallocation error <heat_solve_TRAN, restrt_time>'
              write(idbg,*) '  rank = ', hecMESH%my_rank,'  ierror = ',ierror
              call flush(idbg)
              call hecmw_abort( hecmw_comm_get_comm())
            end if

!C-- file close for local use
      if( hecMESH%my_rank .eq. my_rank_monit_1) then
        CLOSE(heat_IW4)
      endif
      if( hecMESH%my_rank .eq. my_rank_monit_1) then
        CLOSE(heat_IW5)
      endif

   end subroutine heat_solve_TRAN


!C================================================================C
!C-- subroutine heat_output_monit
!C================================================================C
      subroutine heat_output_monit(hecMESH, fstrPARAM, fstrHEAT,    &
                                         heat_IW4, heat_IW5, i_step, CTIME,  &
                                         my_rank_monit_1, my_rank_monit_2)
!
      use m_fstr
!
      implicit none
!C
!C-- global variable
!C
      type ( hecmwST_local_mesh  ) :: hecMESH
      type ( fstr_param          ) :: fstrPARAM
      type ( fstr_heat           ) :: fstrHEAT
!C
!C-- local variable
!C
      integer(kind=kint) :: heat_IW4, heat_IW5
      integer(kind=kint) :: ii,jj
      integer(kind=kint) :: nb_intg, nb_node, nb_elem
      integer(kind=kint) :: i_step
      integer(kind=kint) :: my_rank_monit_1, my_rank_monit_2
      real(kind=kreal) :: CTIME
      real(kind=kreal) :: cure0

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

      if( hecMESH%my_rank .eq. my_rank_monit_1) then
        jj=fstrPARAM%global_local_id(1,nb_node)
        ii=fstrPARAM%global_local_id(2,nb_node)
!!
!C-- temperature
!!
        if(fstrHEAT%iout_list(1) .eq. 1) then
           write(heat_IW4,'(i10,1pe13.4e3,i10,1p1e13.4e3)') i_step &
                        ,CTIME ,jj ,fstrHEAT%TEMP(ii)
           call flush(heat_IW4)
        end if
      end if
!
      if( hecMESH%my_rank .eq. my_rank_monit_1) then
        jj=fstrPARAM%global_local_id(1,nb_node)
        ii=fstrPARAM%global_local_id(2,nb_node)
!!
!C-- cure
!!
        if(fstrHEAT%iout_list(2) .eq. 1) then
           write(heat_IW5,'(i10,1pe13.4e3,i10,1p1e13.4e3)') i_step &
                        ,CTIME ,jj ,fstrHEAT%CURE3D_nod(ii)
           call flush(heat_IW5)
        end if
      end if

      end subroutine heat_output_monit

end module m_heat_solve_TRAN