!=====================================================================
! $B29EYJ,I[7W;;%b%8%e!<%k(B
!  pending
!   - x = 1 and IMAX $B$N%(%M%k%.!<N.F~NL$O%+%&%s%H$7$F$J$$(B
! 2003/12/16 
! 2004/01/01  energy conservation $B$r9MN8$7$?9M$(J}$KJQ99(B
! 2004/01/15  wgt_dammy(KMAX+1)$B$rIU$12C$($k(B
! 2004/01/19  Ms $B$,J#?t$N3J;R$K8Y$k>l9g$ND4@a$,$&$^$/$$$C$F$J$$(B
!=====================================================================
module T_pack
  
  use prm, only : IMAX, JMAX, KMAX,dt, dt_T, dx
  
  use prmphys, only: k_ice, gamma, rho, c_p, Npw, grav, A_con, beta, L_ice
  use com_var, only: zk, wgt, zk_half, Apara
  use Tbound,  only: Ts, Ms


  implicit none
  private
  
  public get_T
  public mk_energy

  public adjust_energy
  public energy2Temp

  public get_Qstrain
  public homo_T


contains

!=====================================================================
! $B3F3J;R$N%(%M%k%.!<$N<0(B
!=====================================================================
subroutine mk_energy(Temp1, energy)
  implicit none
  
  real(8), intent(IN)  :: Temp1(1:IMAX, 1:JMAX, 1:KMAX)
  real(8), intent(OUT) :: energy(1:IMAX, 1:JMAX, 1:KMAX)

  energy =  Temp1 * rho  * c_p

end subroutine mk_energy

!=====================================================================
! $BG.$N<0$N0lIt(B($B3H;69`$N7W;;$rJQ99(B)
! $BDlLL29EY$,M;2rE@0J>e$K$J$k$H2<$+$i$NG.%U%i%C%/%9$r(B 0 $B$H$9$k(B.
!=====================================================================
subroutine get_T(Hthk2, Temp1, Qstrain, Temp2, gamma_flg)
  implicit none

  real(8), parameter   :: one=1.0D0
  real(8), intent(IN)  :: Hthk2(1:IMAX, 1:JMAX)
  real(8), intent(IN)  :: Temp1(1:IMAX, 1:JMAX, 1:KMAX)
  real(8), intent(IN)  :: Qstrain(1:IMAX, 1:JMAX, 1:KMAX)
  real(8), intent(OUT) :: Temp2(1:IMAX, 1:JMAX, 1:KMAX)
  real(8)              :: Temp2L(1:IMAX/2, 1:JMAX, 1:KMAX)
  real(8)              :: Temp2R(1:IMAX/2, 1:JMAX, 1:KMAX)
  real(8)              :: TempH(1:IMAX, 1:JMAX, 1:KMAX)
  real(8)              :: Tpm(1:IMAX, 1:JMAX, KMAX)
!  real(8)              :: Qlat(1:IMAX, 1:JMAX, 1:KMAX)

  !work $BJQ?t(B
  real(8)              :: T_WORK1(1:IMAX, 1:JMAX,1:KMAX)
  real(8)              :: T_WORK2(1:IMAX, 1:JMAX,1:KMAX)
  real(8)              :: T_WORK3(1:IMAX, 1:JMAX,1:KMAX)
  real(8)              :: T_WORK0(1:IMAX, 1:JMAX,1:KMAX)


!  real(8)              :: DQ(1:IMAX, 1:JMAX, 1:KMAX)

  integer              :: i, j, k 
  integer, intent(in)              :: gamma_flg



  call homo_T(Hthk2,  Temp1, Tpm,TempH)
  !  call get_Qlat(Hthk2, Tpm, TempH, Qlat, Temp1)
  

  do j = 1, JMAX
     do i = 1, IMAX
        if(Hthk2(i, j) > epsilon(0.0))then

           ! k=1
           T_WORK1(i, j, 1) = 0.0D0
           
           T_WORK3(i, j, 1) = - k_ice * dt_T/Hthk2(i, j)&
                &  /(zk(2) - zk(1))
           ! $B2<It6-3&$,M;2r29EY0J2<$N>l9g(B
           if((TempH(i, j, 1) < Tpm(i, j, 1)).or.(gamma_flg == 0))then
              T_WORK2(i, j, 1) = Hthk2(i, j)*wgt(1)*rho*c_p &
                   & + k_ice * dt_T/Hthk2(i, j)/wgt(1)        
              T_WORK0(i, j, 1) = Qstrain(i, j, 1)&
                   & + Hthk2(i, j)*wgt(1)*rho*c_p*Temp1(i, j, 1)&
                   & + gamma * dt_T
           ! $B2<It6-3&$,M;2r29EY0J>e$K$J$k>l9g(B 
           else
              T_WORK2(i, j, 1) = Hthk2(i, j)*wgt(1)*rho*c_p &
                   & + k_ice * dt_T/Hthk2(i, j)/wgt(1)        &
                   &   + k_ice *dt_T/Hthk2(i, j)/wgt(1)/2.0D0

              T_WORK0(i, j, 1) = Qstrain(i, j, 1) &
                   & + Hthk2(i, j)*wgt(1)*rho*c_p*Temp1(i, j, 1)&
                   & +k_ice*dt_T*Tpm(i, j, 1)/Hthk2(i, j)/wgt(1)/2.0D0 
           end if


           ! k=KMAX
!           T_WORK1(i, j, KMAX) = 0.0D0
!           T_WORK2(i, j, KMAX) = (Hsfc2(i,j)-Hbnd(i, j))*wgt(KMAX)*rho*c_p
!           T_WORK3(i, j, KMAX) = 0.0D0
!           T_WORK0(i, j, KMAX) = (Hsfc2(i,j)-Hbnd(i, j))*wgt(KMAX)*rho*c_p*Temp1(i, j, KMAX)

           T_WORK1(i, j, KMAX) = -k_ice*dt_T/Hthk2(i, j)&
                & /(zk(KMAX) -zk(KMAX-1))
           T_WORK2(i, j, KMAX) = Hthk2(i, j) * wgt(KMAX) * rho *c_p  &
                & + k_ice*dt_T*(1.0D0/Hthk2(i, j)/wgt(KMAX)/2.0D0 &
                &   + 1.0D0/Hthk2(i, j)/(zk(KMAX)-zk(KMAX-1)))
           T_WORK3(i, j, KMAX) = 0.0D0
           T_WORK0(i, j, KMAX) = Qstrain(i, j, KMAX)&
                & + Hthk2(i, j)*wgt(KMAX)&
                &           * rho*c_p*Temp1(i, j, KMAX)  &
                & + k_ice*dt_T*Ts(i, j)/Hthk2(i, j)/wgt(KMAX)/2.0D0 
           
             
  
           do k = 2, KMAX-1
              T_WORK1(i, j, k) = -k_ice*dt_T/Hthk2(i, j)/(zk(k) -zk(k-1))
              T_WORK2(i, j, k) =  Hthk2(i, j)* wgt(k)*rho *c_p  &
                   & + k_ice*dt_T*(1.0D0/Hthk2(i, j)/(zk(k+1)-zk(k)) &
                   &   + 1.0D0/Hthk2(i, j)/(zk(k)-zk(k-1)))
              T_WORK3(i, j, k) = -k_ice*dt_T/Hthk2(i, j)/(zk(k+1) -zk(k))
              T_WORK0(i, j, k) = Qstrain(i, j, k) &
                   & + Hthk2(i, j)*wgt(k)* rho*c_p*Temp1(i, j, k)
             
          end do


           ! $B9TNs7W;;(B
           call tridag(T_WORK1(i, j, :), T_WORK2(i, j, :), T_WORK3(i, j, :), &
                &   T_WORK0(i, j, :), Temp2(i, j, :))
           
        else
           Temp2(i, j, :) = Ts(i, j)
        end if
        do k =1, KMAX
           ! $BDI5-(B(2004/02/14)
           if(Temp2(i, j, k)  > Tpm(i, j,k))then
              Temp2(i, j, k) = Tpm(i, j, k)
           end if
              
              
           if(Temp2(i, j, k ) < 0.0D0)then
              write(6, *)'T_pack: temperature is error'
              write(6, *)'T_pack: temperature ', 'i=', i, 'j=', j, temp2(i, j, :)
              
              stop
           end if
        end do
        ! energy $BJ]B8$N3NG'(B??
        ! k=1
!        DQ(i ,j, 1) = dx*dt_T*k_ice*(Temp2(i, j, 2)-Temp2(i, j, 1))& 
!                & /(Hsfc2(i, j)-Hbnd(i, j))/(zk(2) - zk(1))   & ! $B>e$X(B
!                + gamma *dt_T *dx                              ! $B2<$+$i(B
       ! k=KMAX
!        DQ(i, j, KMAX) = 0.0D0

!        do k = 2, KMAX-1
!           DQ(i, j, k) = dx*dt_T*k_ice*(Temp2(i, j, k+1)-Temp2(i, j, k))& 
!                & /(Hsfc2(i, j)-Hbnd(i, j))/(zk(k+1) - zk(k))             & ! $B>e$X(B
!                & - dx*dt_T*k_ice*(Temp2(i, j, k)-Temp2(i, j, k-1))& 
!                & /(Hsfc2(i, j)-Hbnd(i, j)/(zk(k) - zk(k-1))               ! $B2<$+$i(B
!        end do

!        if(sum(DQ(i, j, :)) == gamma * dt_T * dx)then
!           if((i==16).and.(j==2))then
!              write(6, *)'DQ=', DQ(i, j, :)
!           end if
!       end if
        
     end do
  end do

  Temp2L(:, :, :)= Temp2(1:15, :, :)
    do i = IMAX/2+2, IMAX
       Temp2R(IMAX-i+1, :, :)= Temp2(i, :, :)
    end do
    
    do j =1, JMAX
       do k=1, KMAX
          do i =1, IMAX/2
             if((Temp2L(i, j, k)/Temp2R(i, j, k)  < 0.99999).or.  &
             (Temp2L(i, j, k)/Temp2R(i, j, k)  > 1.0000001))then
                write(6, *)'Temp2 is not symmetry', i, j, k, Temp2(i, j, k)
             end if
          end do
       end do
    end do

    call monit_add_3('Temp      ', Temp2(:, : , :), one)
    call homo_T(Hthk2, Temp2, Tpm, TempH)
    call monit_add_3('HOMOTEMP  ', TempH(:, : , :), one)

end subroutine get_T


subroutine get_Qstrain(Hsfc2, Hthk2, Qstrain)
  implicit none
  real(8), intent(IN) :: Hsfc2(1:IMAX, 1:JMAX)
  real(8), intent(IN) :: Hthk2(1:IMAX, 1:JMAX)
  real(8), intent(OUT) :: Qstrain(1:IMAX, 1:JMAX, 1:KMAX) ! J/m (Pa m^2)

  integer              :: i, j, k

  Qstrain(1,    :, :) = 0.0D0
  Qstrain(IMAX, :, :) = 0.0D0

  do j = 1, JMAX
     do i = 2,IMAX-1
        do k=1,KMAX
           Qstrain(i, j, k) =  2.0D0 * (-1.0D0)**(Npw+1.0D0) &
                & * (rho*grav)**(Npw+1.0D0)                  &
                & * ((Hsfc2(i+1, j) - Hsfc2(i-1, j))/2.0D0/dx)**(Npw+1.0D0)&
                & * Hthk2(i, j)**(Npw+2.0D0)&
                & * dt_T * wgt(k)/2.0D0             &
                & * ((Apara(i, j, k)*(1.0D0 - zk_half(k))**(Npw+1.0D0)) &
                &      + (Apara(i, j, k)*(1.0D0 - zk_half(k-1))**(Npw+1.0D0)))

        end do
     end do
  end do
  
end subroutine get_Qstrain

!=====================================================================
! $B3J;R$ND4@a(B
!=====================================================================
subroutine adjust_energy(Temp1, Uflgrd, Hthk1, Hthk2, energy_fl, energy2)
  implicit none

  real(8), intent(IN)     :: Uflgrd(   0:IMAX, 1:JMAX, 1:KMAX)   ! m^2/s
  real(8), intent(IN)     :: Temp1(1:IMAX, 1:JMAX, 1:KMAX)       ! K

  real(8), intent(IN)     :: Hthk1(1:IMAX, 1:JMAX)               ! m
  real(8), intent(IN)     :: Hthk2(1:IMAX, 1:JMAX)               ! m
  real(8)                 :: energy1(1:IMAX, 1:JMAX, 1:KMAX+1)! $BD4@aA0$N%(%M%k%.!<(B
  real(8)                 :: energy1_WORK(1:IMAX, 1:JMAX, 1:KMAX+1)! $BD4@aA0$N%(%M%k%.!<(B
  real(8),intent(OUT)      :: energy2(1:IMAX, 1:JMAX, 1:KMAX) ! $BD4@a8e(B

!  real(8)                 :: energy_fl(1:IMAX, 1:JMAX, 1:KMAX)      ! J/m
  real(8),intent(OUT)                 :: energy_fl(1:IMAX, 1:JMAX, 1:KMAX)      ! J/m
  real(8)                 :: energy_org(1:IMAX, 1:JMAX, 1:KMAX)      ! J/m
  real(8)                 :: energy_Ms(1:IMAX, 1:JMAX)      ! J/m
  real(8)                 :: energy_tot(1:IMAX, 1:JMAX)      ! J/m

  real(8)                 :: Uflgrdl_pls(0:IMAX, 1:JMAX, 1:KMAX) !m^2/s
  real(8)                 :: Uflgrdl_mns(0:IMAX, 1:JMAX, 1:KMAX) !m^2/s
  real(8)                 :: Uflgrdr_pls(0:IMAX, 1:JMAX, 1:KMAX) !m^2/s
  real(8)                 :: Uflgrdr_mns(0:IMAX, 1:JMAX, 1:KMAX) !m^2/s
  real(8)                 :: mass_fl(1:IMAX, 1:JMAX, 1:KMAX)
  real(8)                 :: ratio_ab(1:IMAX, 1:JMAX, 1:KMAX)
  real(8)                 :: mass_org(1:IMAX, 1:JMAX, 1:KMAX)
  real(8)                 :: mass(1:IMAX, 1:JMAX, 1:KMAX)
  real(8)                 :: mass_Ms(1:IMAX, 1:JMAX)

  real(8)                 :: wgt_dammy(1:KMAX+1)
  real(8)                 :: wgt_sum


  real(8), parameter      :: non = -9.99E33
  integer                 :: flg 
  integer                 :: i, j, k
  integer                :: kk_str, kk

  
  ! x $BJ}8~6-3&(B
  energy1(1, :, :)    = 0.0D0
  energy1(IMAX, :, :) = 0.0D0
  energy2             = 0.0D0
  energy_fl           = 0.0D0

  !$B$+$sM\$K$h$k<ANLJQ2=(B
  mass_Ms   = Ms*dt*dx*rho
  energy_Ms = Ms*dt*dx*rho*c_p*Ts

  do j =1, JMAX
     do i = 2, IMAX-1
        kk_str = KMAX+1

        do k = 1, KMAX
           ! $B<ANL%U%i%C%/%9$NId9f$rJ,$1$k(B
           ! $B3J;R(B i $B$N:8B&(B
           if(Uflgrd(i-1, j, k) >= 0.0D0)then ! $B3J;R(B i$B$XN.F~(B
              Uflgrdl_pls(i, j, k) = Uflgrd(i-1, j, k) !* Temp1(i-1, j, k)*c_p
              Uflgrdl_mns(i, j, k) = 0.0D0
           else ! $BN.=P(B
              Uflgrdl_pls(i, j, k) = 0.0D0
              Uflgrdl_mns(i, j, k) = Uflgrd(i-1, j, k) !* Temp1(i, j, k)*c_p
          end if
           ! $B3J;R(B i $B$N1&B&(B
           if(Uflgrd(i, j, k) >= 0.0D0)then ! $BN.=P(B
              Uflgrdr_pls(i, j, k) = 0.0D0
              Uflgrdr_mns(i, j, k) = - Uflgrd(i, j, k) ! * Temp1(i, j, k)*c_p
           else ! $BN.F~(B
              Uflgrdr_pls(i, j, k) = - Uflgrd(i, j, k) ! * Temp1(i+1, j, k)c_p
              Uflgrdr_mns(i, j, k) = 0.0D0
           end if

           energy_fl(i, j, k) = &
                &   + (Uflgrdl_pls(i, j, k) * Temp1(i-1, j, k)     &
                &   + (Uflgrdl_mns(i, j, k)+Uflgrdr_mns(i, j, k))*Temp1(i,j,k)&
                &   + Uflgrdr_pls(i, j, k) * Temp1(i+1, j, k))     &
                &   * dt * rho * c_p
           ! $BN.F0$K$h$k<ANLJQ2=(B
           mass_fl(i, j, k) =  &
                &   + (Uflgrdl_pls(i, j, k)+Uflgrdl_mns(i, j, k)   &
                &   + Uflgrdr_pls(i, j, k)+Uflgrdr_mns(i, j, k))   &
                &      * dt * rho 

           ! $B;~9o(B t $B$N;~$N<ANL(B
           mass_org(i, j, k) = Hthk1(i, j)*wgt(k)*dx*rho 

           ! $B;~9o(B t $B$N;~$N%(%M%k%.!<(B
           energy_org(i, j, k)= mass_org(i, j, k)*c_p*Temp1(i, j, k)
           ! $B;~9o(B t+1$B$N;~$N<ANL(B
           mass(i, j, k)  = Hthk2(i, j)*wgt(k)*dx*rho 

           ! $B;~9o(B t+1 $B$N;~$N%(%M%k%.!<(B($B3J;RD4@aA0(B)

        end do
        ! $B>CLW$7$F$$$k>l9g(B
        if((mass_Ms(i,j) < 0.0D0).and.(Hthk2(i, j) > 0.0D0))then

           do k = KMAX, 1, -1
              if(mass_org(i, j, k)+mass_fl(i, j, k)+mass_Ms(i,j) > 0.0D0)then
                 ! $B$=$l$>$l$N<ANLHf$G>CLWNL$r=P$9(B
                 ratio_ab(i, j, k) =  mass_org(i, j, k)           &
                      & /(mass_org(i, j, k)+mass_fl(i, j, k))
                 mass_org(i, j, k) = mass_org(i,j, k)&
                      & + ratio_ab(i, j, k)*mass_Ms(i, j)
                 energy_org(i, j,k)= &
                      & mass_org(i, j, k)*c_p*Temp1(i, j, k)
                 energy_fl(i, j, k) = &
                      & energy_fl(i, j, k)                       &
                      & +(1.0D0-ratio_ab(i,j,k))* mass_Ms(i, j)&
                      & /mass_fl(i, j, k) *energy_fl(i, j, k)
                 mass_fl(i, j, k) = mass_fl(i,j, k)            &
                      &  +(1.0D0-ratio_ab(i,j,k))* mass_Ms(i, j)
                 ! $BCM$r%j%;%C%H(B
                 mass_Ms(i, j) = 0.0D0
                 energy_Ms(i, j) =0.0D0
              else ! 2 $BAX0J>e$K$o$?$k>l9g(B
                 mass_Ms(i, j)=mass_org(i, j, k)+mass_fl(i, j, k)+mass_Ms(i,j)
                 mass_org(i,j ,k) = 0.0D0
                 mass_fl(i, j, k) = 0.0D0
                 energy_org(i, j, k) = 0.0D0
                 energy_fl(i, j,k) = 0.0D0
              end if
           end do
           if((sum(mass_org(i, j, :))+sum(mass_fl(i, j, :)) < 0.0D0))then
              write(6, *)'mass - error', i, j
              stop
           end if
        end if



        ! $BD4@aA0$N%(%M%k%.!<(B
        energy1(i, j, 1:KMAX)=energy_org(i, j, 1:KMAX)+energy_fl(i, j, 1:KMAX)
        energy1(i, j, KMAX+1)=energy_Ms(i, j)



        energy_tot(i, j) = sum(energy1(i, j, :))
!        if((i == 2).and.(j==1))then
!           write(6, *)'energy, before', energy1(i, j, :)
!        end if
        ! $BI9>2$N=PMh$O$8$a(B

        if(Hthk2(i, j) > 0.0D0)then
           wgt_dammy(1:KMAX)=&
                & (mass_fl(i,j,1:KMAX)+mass_org(i, j, 1:KMAX))  &
                &          /sum(mass(i, j, :))
           ! KMAX+1 $BAXL\$O>e$+$i$+$sM\$,2C$o$k(B
           wgt_dammy(KMAX+1) = mass_Ms(i, j)/sum(mass(i, j, :))
        else
           wgt_dammy(:) = non
           energy1(i, j, :) = 0.0D0
           cycle
        end if

! $B%A%'%C%/(B
!           if((i==5).and.(j==3))then
!              write(6, *)i, j,  'energy', sum(energy1(i,j,:)), energy1(i, j, :)
!              write(6, *)i, j,  'wgt_dammy', sum(wgt_dammy), wgt_dammy
!             write(6, *)i, j,  'Hsfc1', Hsfc1(i, j)
!             write(6, *)i, j,  'Hsfc2', Hsfc2(i, j)
!          end if

        energy1_WORK(i, j, :) = energy1(i, j, :)  
        wgt_sum = wgt_dammy(kk_str)     ! $B@c$rBeF~(B

        flg = 0
        do k=KMAX, 1,  -1               ! $BD4@a8e$N3J;R(B

           do kk = kk_str, 1, -1        ! $BD4@aA0$N3J;R(B

              if(wgt(k) >= wgt_sum)then  ! wgt_dammy(kk)$B$NA4It(B
                 flg = 0                 !  $B0lIt$,F~$k%U%i%0$rLa$9(B
                 energy2(i, j, k) = energy2(i, j, k)&
                      & +energy1_WORK(i, j, kk)
                 if(kk >= 2)then
                 ! $BF~$k$J$i2<$NAX$rB-$7$F$^$?%A%'%C%/(B
                    wgt_sum       = wgt_sum + wgt_dammy(kk-1)
                 else
                 !kk=1 $B$rB-$7$F$bM>$k>l9g(B
                    if((wgt_sum/wgt(k) > 1.000000001).or.&
                         (wgt_sum/wgt(k) > 0.9999999999))then
                       write(6, *)' error!! mass is loss', wgt_sum
                       stop
                    end if

                 end if

                 
              else                      ! wgt_dammy(kk)$B$N0lIt$,F~$k(B
                 if(flg == 0)then ! 1 $B2sL\(B
                    energy2(i, j, k) = energy2(i, j, k)&
                         & +(wgt(k)-(wgt_sum-wgt_dammy(kk)) )/wgt_dammy(kk) &
                         &   * energy1(i, j, kk)
                    
                    ! $B;D$j(B ! wgt(k) >= wgt_sum $B$G2C$($i$l$k(Benergy1
                    energy1_WORK(i, j, kk) = &
                         & (wgt_sum-wgt(k))/wgt_dammy(kk) &
                         &   * energy1(i, j, kk)  
                    wgt_sum = wgt_sum-wgt(k) ! $BM>$j(B


                    kk_str = kk
                    flg = 1
                    exit  
                 else             ! 2 $B2sL\0J9_(B
                    energy2(i, j, k) = energy2(i, j, k)&
                         & +wgt(k)/wgt_dammy(kk) * energy1(i, j, kk)
                    ! $B;D$j(B ! wgt(k) >= wgt_sum $B$G2C$($i$l$k(Benergy1
                    energy1_WORK(i, j, kk) = &
                         & (wgt_sum-wgt(k))/wgt_dammy(kk) &
                         &   * energy1(i, j, kk)  
                    wgt_sum = wgt_sum-wgt(k) ! $BM>$j(B
                    kk_str = kk
                    flg = 1
                 
                    exit  

                 end if
              end if
                 
           end do
              
        end do
        ! $B3J;R$N@Z$jJ,$1$NA08e$G$NJ]B8$N%A%'%C%/(B
        if(sum(energy1(i, j, :)) ==0.0D0)cycle
        if(((sum(energy2(i, j, :))/energy_tot(i,j) > 1.0D0+epsilon(0.0)).or. &
             &(sum(energy2(i, j, :))/energy_tot(i,j) < 1.0D0-epsilon(0.0))).and.(energy_tot(i, j) > epsilon(0.0)))then
           write(6, *)'energy is not conserved. postion is', i, j
           write(6, *) 'energy = '  , sum(energy2(i, j, :))/energy_tot(i,j)
           write(6, *)'energy before', energy_tot(i, j)
           write(6, *)'energy after', sum(energy2(i, j, :)), energy2(i, j, :)
           stop
        end if
        
     end do
  end do

end subroutine adjust_energy

!=====================================================================
!$BFbIt%(%M%k%.!<$+$i29EY$r;;=P(B
!=====================================================================
subroutine energy2Temp(energy, Hthk2,  Temp)
  implicit none
  
  real(8), INTENT(IN) :: energy(1:IMAX, 1:JMAX, 1:KMAX)
  real(8), INTENT(IN) :: Hthk2(1:IMAX, 1:JMAX)
  real(8), INTENT(OUT):: Temp(1:IMAX, 1:JMAX, 1:KMAX)
  real(8)             :: mass(1:IMAX, 1:JMAX, 1:KMAX)
  integer             :: i, j, k

  do i =1, IMAX
     do j = 1, JMAX
        if(Hthk2(i, j) > 0.0D0)then
           do k =1, KMAX
              mass(i, j, k) = Hthk2(i, j)*wgt(k)*dx*rho
              Temp(i, j, k) = energy(i, j, k)/mass(i, j,k)/c_p
           end do
        else
           Temp(i, j, :) = Ts(i, j)
        end if
     end do
  end do

  
end subroutine energy2Temp
  
!=====================================================================
!$B05NOM;2rE@$r4p=`$K$7$?29EY(B
!=====================================================================
subroutine homo_T(Hthk2, Temp2, Tpm, TempH)
  implicit none

  real(8) ,intent(IN) :: Hthk2(1:IMAX, 1:JMAX)
  real(8) ,intent(IN) :: Temp2(1:IMAX, 1:JMAX, 1:KMAX)
  real(8) ,intent(OUT) :: TempH(1:IMAX, 1:JMAX, 1:KMAX)
  real(8) ,intent(OUT) :: Tpm(1:IMAX, 1:JMAX, 1:KMAX)

  integer              :: k

  do k =1, KMAX
     TempH(:, :, k) = Temp2(:, :, k) &
          & + beta * Hthk2(:, :)*(1.0D0 - zk(k))
     Tpm(:, :, k) = 273.15 - beta *Hthk2(:, :)*(1.0D0 - zk(k))
  end do

  
end subroutine homo_T

!=====================================================================
! $B@xG.$K$h$k%(%M%k%.!<(B
!=====================================================================

subroutine get_Qlat(Hthk2, Tpm, TempH, Qlat, Temp2)
  implicit none
  
  real(8), intent(IN)  :: Hthk2(1:IMAX, 1:JMAX)
  real(8), intent(IN)  :: Tpm(1:IMAX, 1:JMAX, 1:KMAX)  ! $BM;2r29EY(B
  real(8), intent(IN)  :: TempH(1:IMAX, 1:JMAX, 1:KMAX)
  real(8), intent(OUT)  :: Qlat(1:IMAX, 1:JMAX, 1:KMAX)
  real(8), intent(INOUT)  :: Temp2(1:IMAX, 1:JMAX, 1:KMAX)

  integer              :: i, j, k
  Qlat = 0.0D0
  do i =1, IMAX
     do j = 1, JMAX
        do k = 1, KMAX
           if(TempH(i, j, k)   > Tpm(i, j, k))then
              Qlat(i, j, k)=L_ice * Hthk2(i, j)*wgt(k)*rho
              Temp2(i, j, k) = Tpm(i, j, k)
           end if
        end do
     end do
  end do



end subroutine get_Qlat




!=====================================================================
! $B9TNs7W;;MQ%5%V%k!<%A%s(B
!=====================================================================
subroutine tridag(a1, a2, a3, vectI, vectX)
  implicit none
  real(8), INTENT(IN)   :: a1(1:KMAX)
  real(8), INTENT(IN)   :: a2(1:KMAX)
  real(8), INTENT(IN)   :: a3(1:KMAX)
  real(8), INTENT(IN)   :: vectI(1:KMAX)
  real(8), INTENT(INOUT)   :: vectX(1:KMAX)

  integer               :: k
  real(8)               :: bet, gam(1:KMAX)

  if(a2(1) == 0)stop 'tridag: rewrite equations'

  bet = a2(1)
  vectX(1)=vectI(1)/bet

  do k = 1, KMAX
     gam(k) = a3(k-1)/bet
     bet=a2(k)-a1(k)*gam(k)
     if(bet == 0.0D0) stop 'tridag failed'
     vectX(k) = (vectI(k)-a1(k)*vectX(k-1))/bet
  end do

  do k = KMAX-1, 1, -1
     vectX(k) = vectX(k) -gam(k+1)*vectX(k+1) 
  end do



end subroutine tridag


end module T_pack
