module height
  use prm, only : IMAX, JMAX, DX, dt
  use prmphys, only :  s_a
  use Tbound, only: Ms
  use com_var, only: Apara
  implicit none

  private 
  public get_H
  public add_Ms
  public mk_thick

contains

  ! $BI98|$r7W;;(B
  subroutine mk_thick(Hsfc1, Hbnd, Hthk1)
    implicit none
    

    real(8), INTENT(IN)  :: Hsfc1(1:IMAX, 1:JMAX)
    real(8), INTENT(IN)  :: Hbnd(1:IMAX, 1:JMAX)
    real(8), INTENT(OUT) :: Hthk1(1:IMAX, 1:JMAX)
    integer              :: i , j 
    
    do i = 1,IMAX
       do j =1, JMAX
          if(Hsfc1(i, j) >= 0.0D0)then
             Hthk1(i, j) = Hsfc1(i, j)-Hbnd(i, j)
          elseif((Hsfc1(i, j) < 0.0D0).and.(Hbnd(i, j) < 0.0D0))then
             Hthk1(i, j) = -(Hsfc1(i, j) - Hbnd(i, j))
          else
             write(6, *)'Hsfc < 0 Hbnd >0 , error !! mk_thick'
             write(6, *)'Hsfc=', i, j, Hsfc1(i, j)
             write(6, *)'Hsfc=', i, j, Hbnd(i, j)
             stop 
          end if
          if(Hthk1(i, j) < 0.0D0)then
             write(6, *)'Height is -minus value. error !!  mk_thick'
             stop
          end if
       end do
    end do


  end subroutine mk_thick

  subroutine add_Ms(Hthk2)
    implicit none


    real(8), parameter   :: one = 1.0D0
    real(8), INTENT(INOUT)  :: Hthk2(1:IMAX, 1:JMAX)
    integer  :: i, j    


    do j = 1, JMAX
       Hthk2(1, j) = 0.D0
       Hthk2(IMAX, j) =0.D0
       do i = 2, IMAX-1
          Hthk2(i, j)= Hthk2(i, j) &
               & + Ms(i, j) * dt
          if(Hthk2(i, j) <= 0.0D0)then
             Hthk2(i, j) = 0.0D0
          end if
       end do
    end do
    call monit_add_2('HIGHT     ', Hthk2, one)

  end subroutine add_Ms

  subroutine get_H(Hthk1, Ufl, Hthk2)
    implicit none
    real(8), INTENT(IN)   :: Hthk1(1:IMAX, 1:JMAX)
    real(8), INTENT(IN)   :: Ufl(0:IMAX, 1:JMAX)
    real(8), INTENT(OUT)  :: Hthk2(1:IMAX, 1:JMAX)

    real(8)               :: Hthk2L(1:IMAX/2, 1:JMAX)
    real(8)               :: Hthk2R(1:IMAX/2, 1:JMAX)

    integer  :: i, j
  
    do j = 1, JMAX
       Hthk2(1, j) = 0.D0
       Hthk2(IMAX, j) =0.D0
       do i = 2, IMAX-1
          Hthk2(i, j)= Hthk1(i, j) &
               &  - dt/dx*(Ufl(i, j) - Ufl(i-1, j))!+ Ms(i, j) * dt
          if(Hthk2(i, j) <= 0.0D0)then
             Hthk2(i, j) = 0.0D0
          end if
       end do
    end do


    call monit_add_2('uflux     ', ufl(1:IMAX, :), dt/s_a)


    Hthk2L(:, :) = Hthk2(1:IMAX/2, :)
    do  i = 1, IMAX/2
       Hthk2R(i, :) = Hthk2(IMAX+1-i, :)
    end do

    do j = 1, JMAX
       do i = 1, IMAX/2
          if((Hthk2L(i, j)/Hthk2R(i, j)  < 0.9999).or.  &
               (Hthk2L(i, j)/Hthk2R(i, j)  > 1.0001))then
             write(6, *)i, j, 'Hthk', Hthk2L(i, j), Hthk2R(i, j)
             write(6, *)'AparaL', Apara(i, j, :)
             write(6, *)'AparaR', Apara(IMAX+1-i, j, :)
          end if
       end do
    end do



!    $BI9>2$NC<N.$l=P$5$J$$LdBj$N%A%'%C%/(B
!    do j = 1, JMAX
!       do i = IMAX, 1, -1
!          if(Hsfc1(i, j)>0D0)then
!             write(6, *)'Hsfc2=', Hsfc2(i-1, j), Hsfc2(i, j), Hsfc2(i+1, j)
!             write(6, *)'Hsfc1=',Hsfc1(i-1, j), Hsfc1(i, j), Hsfc1(i+1, j)
!             write(6, *)'Ufl(i)=', Ufl(i-1, j),  Ufl(i, j), Ufl(i+1, j)
!             write(6, *)'Ufl(i-1)=' , Ufl(i-2, j), Ufl(i-1, j), Ufl(i, j)
!             write(6, *)'-dt*(Ufl(i, j) - Ufl(i-1, j))/dx=', &
!                  & - dt*(Ufl(i-1, j) - Ufl(i-2, j))/dx, &
!                  & - dt*(Ufl(i, j) - Ufl(i-1, j))/dx, &
!                  & - dt*(Ufl(i+1, j) - Ufl(i, j))/dx
!             exit
!          end if
!
!       end do
!    end do


end subroutine get_H


end module height
