!---------------------------------------------------------------------
!     Copyright (C) GFD Dennou Club, 2005. All rights reserved.
!---------------------------------------------------------------------
!= Subroutine FillNegative
!
!   * Developer: SUGIYAMA Ko-ichiro
!   * Version: $Id: $
!   * Tag Name: $Name:  $
!   * Change History: 
!
!== Overview
!
!負の雲水量などを穴埋めするためのルーチン
!  * 負となった場合には, 周囲の 8 格子点の値を削って穴埋めする.
! 
!== Error Handling
!
!== Bugs
!
!== Note
!
!  * 中島版 deepconv の QFILL.f を改変
!
!== Future Plans
!
!

subroutine FillNegative( xz_DensBasicZ, xz_Var )
  !
  !負の雲水量などを穴埋めするためのルーチン
  !  * 負となった場合には, 周囲の 8 格子点の値を削って穴埋めする.
  ! 

  !モジュール読み込み
  use gridset,    only: DimXMin,  &! 配列の X 方向の下限
    &                   DimXMax,  &! 配列の X 方向の上限
    &                   DimZMin,  &! 配列の Z 方向の下限
    &                   DimZMax    ! 配列の Z 方向の上限
!  use basicset, only: xz_DensBasicZ
  
  implicit none

  real(8), intent(inout) :: xz_Var(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), intent(in)    :: xz_DensBasicZ(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8)                :: xz_VarFill(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8)                :: xz_DQFILL(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8)                :: xz_QSUMPN(DimXMin:DimXMax, DimZMin:DimZMax)
  real(8), parameter     :: EPS = 1.0d-60  !零での割り算を防ぐ
  integer                :: i, k

  !初期化
  xz_QSUMPN = 0.0d0
  xz_DQFILL = 0.0d0


  do k = DimZMin+2, DimZMax-2
    do i = DimXMin+2, DimXMax-2
      xz_QSUMPN(i,k) = 1.0d0 / &
        &  ( (   MAX( 0.0d0, xz_Var(i-1,k) )* xz_DensBasicZ(i-1,k) &
        &      + MAX( 0.0d0, xz_Var(i+1,k) )* xz_DensBasicZ(i+1,k) &
        &      + MAX( 0.0d0, xz_Var(i,k-1) )* xz_DensBasicZ(i,k-1) &
        &      + MAX( 0.0d0, xz_Var(i,k+1) )* xz_DensBasicZ(i,k+1) &
        &     ) * 0.75d0   &
        &  + (   MAX( 0.0d0, xz_Var(i-2,k) )* xz_DensBasicZ(i-2,k) &
        &      + MAX( 0.0d0, xz_Var(i+2,k) )* xz_DensBasicZ(i+2,k) &
        &      + MAX( 0.0d0, xz_Var(i,k-2) )* xz_DensBasicZ(i,k-2) &
        &      + MAX( 0.0d0, xz_Var(i,k+2) )* xz_DensBasicZ(i,k+2) &
        &     ) * 0.25d0 + EPS ) 
    end do
  end do

  do k = DimZMin+2, DimZMax-2
    do i = DimXMin+2, DimXMax-2
      xz_DQFILL(i,k) =         &
        &  - MIN( 0.0d0, xz_Var(i,k) ) &
        &  + MAX( 0.0d0, xz_Var(i,k) ) &
        &    * ( ( MIN( 0.0d0, xz_Var(i-1,k) ) * xz_QSUMPN(i-1,k) &
        &          * xz_DensBasicZ(i-1,k)     &
        &        + MIN( 0.0d0, xz_Var(i+1,k) ) * xz_QSUMPN(i+1,k) &
        &          * xz_DensBasicZ(i+1,k)     &
        &        + MIN( 0.0d0, xz_Var(i,k-1) ) * xz_QSUMPN(i,k-1) &
        &          * xz_DensBasicZ(i,k-1)   &
        &        + MIN( 0.0d0, xz_Var(i,k+1) ) * xz_QSUMPN(i,k+1) &
        &          * xz_DensBasicZ(i,k+1)   & 
        &         ) * 0.75d0                    &
        &      + ( MIN( 0.0d0, xz_Var(i-2,k) ) * xz_QSUMPN(i-2,k) &
        &           * xz_DensBasicZ(i-2,k)    &
        &        + MIN( 0.0d0, xz_Var(i+2,k) ) * xz_QSUMPN(i+2,k) &
        &           * xz_DensBasicZ(i+2,k)    &
        &        + MIN( 0.0d0, xz_Var(i,k-2) ) * xz_QSUMPN(i,k-2) &
        &           * xz_DensBasicZ(i,k-2)  &
        &        + MIN( 0.0d0, xz_Var(i,k+2) ) * xz_QSUMPN(i,k+2) &
        &           * xz_DensBasicZ(i,k+2)  &
        &        ) * 0.25d0 &
        &      )
    end do
  end do

  xz_VarFill = xz_Var + xz_DQFILL

  !出力
  xz_Var = xz_VarFill


end subroutine FillNegative