function xza_FillNegative_xza( xza_Var )
implicit none
real(8), intent(inout) :: xza_Var(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
real(8) :: xza_FillNegative_xza(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
real(8) :: xza_DQFILL(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
real(8) :: xza_QSUMPN(DimXMin:DimXMax, DimZMin:DimZMax, SpcNum)
real(8), parameter :: EPS = 1.0d-60 !零での割り算を防ぐ
integer :: i, k, s
!初期化
xza_QSUMPN = 0.0d0
xza_DQFILL = 0.0d0
! do s = 1, SpcNum
! write(*,*) s, xza_Min(2:20,48,s)
! end do
do s = 1, SpcNum
do k = DimZMin+2, DimZMax-2
do i = DimXMin+2, DimXMax-2
xza_QSUMPN(i,k,s) = 1.0d0 / ( ( MAX( 0.0d0, xza_Basic(i-1,k,s) + xza_Var(i-1,k,s) ) * xz_Dens(i-1,k) + MAX( 0.0d0, xza_Basic(i+1,k,s) + xza_Var(i+1,k,s) ) * xz_Dens(i+1,k) + MAX( 0.0d0, xza_Basic(i,k-1,s) + xza_Var(i,k-1,s) ) * xz_Dens(i,k-1) + MAX( 0.0d0, xza_Basic(i,k+1,s) + xza_Var(i,k+1,s) ) * xz_Dens(i,k+1) ) * 0.75d0 + ( MAX( 0.0d0, xza_Basic(i-2,k,s) + xza_Var(i-2,k,s) ) * xz_Dens(i-2,k) + MAX( 0.0d0, xza_Basic(i+2,k,s) + xza_Var(i+2,k,s) ) * xz_Dens(i+2,k) + MAX( 0.0d0, xza_Basic(i,k-2,s) + xza_Var(i,k-2,s) ) * xz_Dens(i,k-2) + MAX( 0.0d0, xza_Basic(i,k+2,s) + xza_Var(i,k+2,s) ) * xz_Dens(i,k+2) ) * 0.25d0 + EPS )
end do
end do
end do
do s = 1, SpcNum
do k = DimZMin+2, DimZMax-2
do i = DimXMin+2, DimXMax-2
xza_DQFILL(i,k,s) = - MIN( 0.0d0, xza_Basic(i,k,s) + xza_Var(i,k,s) ) + MAX( 0.0d0, xza_Basic(i,k,s) + xza_Var(i,k,s) ) * ( ( MIN( 0.0d0, xza_Basic(i-1,k,s) + xza_Var(i-1,k,s) ) * xza_QSUMPN(i-1,k,s) * xz_Dens(i-1,k) + MIN( 0.0d0, xza_Basic(i+1,k,s) + xza_Var(i+1,k,s) ) * xza_QSUMPN(i+1,k,s) * xz_Dens(i+1,k) + MIN( 0.0d0, xza_Basic(i,k-1,s) + xza_Var(i,k-1,s) ) * xza_QSUMPN(i,k-1,s) * xz_Dens(i,k-1) + MIN( 0.0d0, xza_Basic(i,k+1,s) + xza_Var(i,k+1,s) ) * xza_QSUMPN(i,k+1,s) * xz_Dens(i,k+1) ) * 0.75d0 + ( MIN( 0.0d0, xza_Basic(i-2,k,s) + xza_Var(i-2,k,s) ) * xza_QSUMPN(i-2,k,s) * xz_Dens(i-2,k) + MIN( 0.0d0, xza_Basic(i+2,k,s) + xza_Var(i+2,k,s) ) * xza_QSUMPN(i+2,k,s) * xz_Dens(i+2,k) + MIN( 0.0d0, xza_Basic(i,k-2,s) + xza_Var(i,k-2,s) ) * xza_QSUMPN(i,k-2,s) * xz_Dens(i,k-2) + MIN( 0.0d0, xza_Basic(i,k+2,s) + xza_Var(i,k+2,s) ) * xza_QSUMPN(i,k+2,s) * xz_Dens(i,k+2) ) * 0.25d0 )
end do
end do
! write(*,*) sum( xza_DQFILL(DimXMin:DimXMax,DimZMin:DimZMax,s) )
end do
!出力
xza_FillNegative_xza = xza_Var + xza_DQFILL
end function xza_FillNegative_xza