!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2005. All rights reserved.
!---------------------------------------------------------------------
                                                                 !=begin
!= Subroutine FillNegative
!
!   * Developer: SUGIYAMA Ko-ichiro
!   * Version: $Id: $
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview
!
!負の雲水量などを穴埋めするためのルーチン
!  * 負となった場合には, 周囲の 8 格子点の値を削って穴埋めする.
! 
!== Error Handling
!
!== Bugs
!
!== Note
!
!  * 中島版 deepconv の QFILL.f を改変
!  * 元々のプログラムの ΔZ は不等間隔なので, ΔZ を乗じた形で書いてある
!  * 中島さん曰く, 「適当に作ったよ, 保存性はチェックしたよ」
!
!== Future Plans
!
!
subroutine FillNegative( ss_Var )

  use gridset,  only: DimXMin, DimXMax, DimZMin, DimZMax, DelZ
  use basicset, only: ss_DensBasicZ
  use bcset,    only: ss_BC

  implicit none

  real(8), intent(inout) :: ss_Var(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8)                :: ss_VarFill(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8)                :: ss_DQFILL(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8)                :: ss_QSUMPN(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), parameter     :: EPS = 1.0d-60  !零での割り算を防ぐ
  integer                :: i, k

  do k = DimZMin+2, DimZMax-2
     do i = DimXMin+2, DimXMax-2
        ss_QSUMPN(i,k) = 1.0d0 / &
             &  ( (   MAX( 0.0d0, ss_Var(i-1,k) )* ss_DensBasicZ(i-1,k)* DelZ &
             &      + MAX( 0.0d0, ss_Var(i+1,k) )* ss_DensBasicZ(i+1,k)* DelZ &
             &      + MAX( 0.0d0, ss_Var(i,k-1) )* ss_DensBasicZ(i,k-1)* DelZ &
             &      + MAX( 0.0d0, ss_Var(i,k+1) )* ss_DensBasicZ(i,k+1)* DelZ &
             &     ) * 0.75d0   &
             &  + (   MAX( 0.0d0, ss_Var(i-2,k) )* ss_DensBasicZ(i-2,k)* DelZ &
             &      + MAX( 0.0d0, ss_Var(i+2,k) )* ss_DensBasicZ(i+2,k)* DelZ &
             &      + MAX( 0.0d0, ss_Var(i,k-2) )* ss_DensBasicZ(i,k-2)* DelZ &
             &      + MAX( 0.0d0, ss_Var(i,k+2) )* ss_DensBasicZ(i,k+2)* DelZ &
             &     ) * 0.25d0 + EPS ) 
     end do
  end do
!  do k = DimZMin, DimZMax 
!     write(*,*) k, "QSUMPN: ", ss_QSUMPN(50,k)
!     write(*,*) k, "Var:    ", ss_Var(50,k), "( ", ss_Var(48,k), " )" 
!     write(*,*) k, "MAX:    ", MAX( 0.0d0, ss_Var(49,k) ),   &
!          &                    MAX( 0.0d0, ss_Var(51,k) ),   &
!          &                    MAX( 0.0d0, ss_Var(50,k-1) ), &
!          &                    MAX( 0.0d0, ss_Var(50,k+1) ), &
!          &                    MAX( 0.0d0, ss_Var(48,k) ),   &
!          &                    MAX( 0.0d0, ss_Var(52,k) ),   &
!          &                    MAX( 0.0d0, ss_Var(50,k-2) ), &
!          &                    MAX( 0.0d0, ss_Var(50,k+2) )
!  end do
  call boundary( ss_Bc, ss_QSUMPN )

  do k = DimZMin+2, DimZMax-2
     do i = DimXMin+2, DimXMax-2
        ss_DQFILL(i,k) =         &
             &  - MIN( 0.0d0, ss_Var(i,k) ) &
             &  + MAX( 0.0d0, ss_Var(i,k) ) &
             &    * ( ( MIN( 0.0d0, ss_Var(i-1,k) ) * ss_QSUMPN(i-1,k) &
             &          * ss_DensBasicZ(i-1,k) * DelZ    &
             &        + MIN( 0.0d0, ss_Var(i+1,k) ) * ss_QSUMPN(i+1,k) &
             &          * ss_DensBasicZ(i+1,k) * DelZ    &
             &        + MIN( 0.0d0, ss_Var(i,k-1) ) * ss_QSUMPN(i,k-1) &
             &          * ss_DensBasicZ(i,k-1) * DelZ  &
             &        + MIN( 0.0d0, ss_Var(i,k+1) ) * ss_QSUMPN(i,k+1) &
             &          * ss_DensBasicZ(i,k+1) * DelZ  & 
             &         ) * 0.75d0                    &
             &      + ( MIN( 0.0d0, ss_Var(i-2,k) ) * ss_QSUMPN(i-2,k) &
             &           * ss_DensBasicZ(i-2,k) * DelZ   &
             &        + MIN( 0.0d0, ss_Var(i+2,k) ) * ss_QSUMPN(i+2,k) &
             &           * ss_DensBasicZ(i+2,k) * DelZ   &
             &        + MIN( 0.0d0, ss_Var(i,k-2) ) * ss_QSUMPN(i,k-2) &
             &           * ss_DensBasicZ(i,k-2) * DelZ &
             &        + MIN( 0.0d0, ss_Var(i,k+2) ) * ss_QSUMPN(i,k+2) &
             &           * ss_DensBasicZ(i,k+2) * DelZ &
             &        ) * 0.25d0 &
             &      )
     end do
  end do
  call boundary( ss_Bc, ss_DQFILL )


  ss_VarFill = ss_Var + ss_DQFILL
  call boundary( ss_Bc, ss_VarFill )

  !出力
  ss_Var = ss_VarFill


end subroutine FillNegative
