!= Module FillNegative
!
! Authors::   SUGIYAMA Ko-ichiro, ODAKA Masatsugu
! Version::   $Id: fillnegative.f90,v 1.1 2011-04-26 04:08:11 sugiyama Exp $
! Tag Name::  $Name:  $
! Copyright:: Copyright (C) GFD Dennou Club, 2008. All rights reserved.
! License::   See COPYRIGHT[link:../../COPYRIGHT]
!
!== Overview
!
!負の雲水量などを穴埋めするためのルーチン
!  * 負となった場合には, 周囲の 8 格子点の値を削って穴埋めする.
! 
!== Error Handling
!
!== Bugs
!
!== Note
!
!  * 中島版 deepconv の QFILL.f を改変
!
!== Future Plans
!
!

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

  !モジュール読み込み
  use dc_types,   only: DP
  use dc_message, only: MessageNotify

  use gridset, only: imin,  &! 配列の X 方向の下限
    &                imax,  &! 配列の X 方向の上限
    &                jmin,  &! 配列の Y 方向の下限
    &                jmax,  &! 配列の Y 方向の上限
    &                kmin,  &! 配列の Z 方向の下限
    &                kmax,  &! 配列の Z 方向の上限
    &                nf      ! 化学種の数
 
  implicit none

  public FillNegative_Init
  public xyza_FillNegative_xyza

  real(DP), save, allocatable  :: xyza_Basic(:,:,:,:)
  real(DP), save, allocatable  :: xyz_Dens(:,:,:)
  real(DP), save, allocatable  :: xza_Basic(:,:,:)
  real(DP), save, allocatable  :: xz_Dens(:,:)

contains

  subroutine FillNegative_Init( xyza_VarBasic, xyz_DensBasic )
    !
    ! 基準値を取得. 
    ! 基本場と擾乱成分の和がゼロよりも小さくなることは許容しない
    !

    implicit none

    !入力変数
    real(DP), intent(in) :: xyza_VarBasic(imin:imax,jmin:jmax,kmin:kmax, nf)
    real(DP), intent(in) :: xyz_DensBasic(imin:imax,jmin:jmax,kmin:kmax)

    !初期化
    allocate( xyza_Basic(imin:imax,jmin:jmax,kmin:kmax, nf) )
    allocate( xyz_Dens  (imin:imax,jmin:jmax,kmin:kmax) )
    allocate( xza_Basic(imin:imax,kmin:kmax, nf) )
    allocate( xz_Dens(imin:imax,kmin:kmax ) )
    
    !値の代入
    xyza_Basic = xyza_VarBasic
    xyz_Dens   = xyz_DensBasic

    xza_Basic = xyza_VarBasic(:,1,:,:)
    xz_Dens   = xyz_DensBasic(:,1,:)

  end subroutine FillNegative_Init
  

  function xyza_FillNegative_xyza( xyza_Var )
    
    implicit none

    real(DP), intent(in) :: xyza_Var(imin:imax,jmin:jmax,kmin:kmax, nf)
    real(DP) :: xyza_FillNegative_xyza(imin:imax,jmin:jmax,kmin:kmax, nf)
    real(DP) :: xyza_DQFILL(imin:imax,jmin:jmax,kmin:kmax, nf)
    real(DP) :: xyza_QSUMPN(imin:imax,jmin:jmax,kmin:kmax, nf)
    real(DP) :: xza_Var(imin:imax, kmin:kmax, nf)
    real(DP) :: xza_DQFILL(imin:imax, kmin:kmax, nf)
    real(DP) :: xza_QSUMPN(imin:imax, kmin:kmax, nf)
    real(DP), parameter     :: EPS = 1.0d-60  !零での割り算を防ぐ
    integer                :: i, j, k, s

    
    if (jmin==jmax) then 

      !初期化
      xza_QSUMPN = 0.0d0
      xza_DQFILL = 0.0d0
      xza_Var(:,:,:) = xyza_Var(:,1,:,:)
      
      do s = 1, nf
        do k = kmin+2, kmax-2
          do i = imin+2, imax-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, nf
        do k = kmin+2, kmax-2
          do i = imin+2, imax-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
      end do
      
      !出力
      xyza_FillNegative_xyza(:,1,:,:) = xza_Var + xza_DQFILL
      
    else
      
      !初期化
      xyza_QSUMPN = 0.0d0
      xyza_DQFILL = 0.0d0
      
      do s = 1, nf
        do k = kmin+2, kmax-2
          do j = jmin+2, jmax-2
            do i = imin+2, imax-2
              xyza_QSUMPN(i,j,k,s) = 1.0d0 /                                       &
                & ( (    MAX( 0.0d0, xyza_Basic(i-1,j,k,s) + xyza_Var(i-1,j,k,s) ) &
                &         * xyz_Dens(i-1,j,k)                                      &
                &      + MAX( 0.0d0, xyza_Basic(i+1,j,k,s) + xyza_Var(i+1,j,k,s) ) &
                &         * xyz_Dens(i+1,j,k)                                      &
                &      + MAX( 0.0d0, xyza_Basic(i,j-1,k,s) + xyza_Var(i,j-1,k,s) ) &
                &         * xyz_Dens(i,j-1,k)                                      &
                &      + MAX( 0.0d0, xyza_Basic(i,j+1,k,s) + xyza_Var(i,j+1,k,s) ) &
                &         * xyz_Dens(i,j+1,k)                                      &
                &      + MAX( 0.0d0, xyza_Basic(i,j,k-1,s) + xyza_Var(i,j,k-1,s) ) &
                &         * xyz_Dens(i,j,k-1)                                      &
                &      + MAX( 0.0d0, xyza_Basic(i,j,k+1,s) + xyza_Var(i,j,k+1,s) ) &
                &         * xyz_Dens(i,j,k+1)                                      &
                &     ) * 0.75d0                                                   &
                &  + (   MAX( 0.0d0, xyza_Basic(i-2,j,k,s) + xyza_Var(i-2,j,k,s) ) &
                &         * xyz_Dens(i-2,j,k)                                      &
                &      + MAX( 0.0d0, xyza_Basic(i+2,j,k,s) + xyza_Var(i+2,j,k,s) ) &
                &         * xyz_Dens(i+2,j,k)                                      &
                &      + MAX( 0.0d0, xyza_Basic(i,j-2,k,s) + xyza_Var(i,j-2,k,s) ) &
                &         * xyz_Dens(i,j-2,k)                                      &
                &      + MAX( 0.0d0, xyza_Basic(i,j+2,k,s) + xyza_Var(i,j+2,k,s) ) &
                &         * xyz_Dens(i,j+2,k)                                      &
                &      + MAX( 0.0d0, xyza_Basic(i,j,k-2,s) + xyza_Var(i,j,k-2,s) ) &
                &         * xyz_Dens(i,j,k-2)                                      &
                &      + MAX( 0.0d0, xyza_Basic(i,j,k+2,s) + xyza_Var(i,j,k+2,s) ) &
                &         * xyz_Dens(i,j,k+2)                                      &
                &     ) * 0.25d0                                              &
                &  + EPS ) 
            end do
          end do
        end do
      end do
      
      do s = 1, nf
        do k = kmin+2, kmax-2
          do j = jmin+2, jmax-2
            do i = imin+2, imax-2
              xyza_DQFILL(i,j,k,s) =                                              &
                &  - MIN( 0.0d0, xyza_Basic(i,j,k,s) + xyza_Var(i,j,k,s) )        &
                &  + MAX( 0.0d0, xyza_Basic(i,j,k,s) + xyza_Var(i,j,k,s) )        &
                &    * ( ( MIN( 0.0d0, xyza_Basic(i-1,j,k,s) + xyza_Var(i-1,j,k,s) ) &
                &           * xyza_QSUMPN(i-1,j,k,s)                               &
                &           * xyz_Dens(i-1,j,k)                                    &
                &        + MIN( 0.0d0, xyza_Basic(i+1,j,k,s) + xyza_Var(i+1,j,k,s) ) &
                &           * xyza_QSUMPN(i+1,j,k,s)                               &
                &           * xyz_Dens(i+1,j,k)                                    &
                &        + MIN( 0.0d0, xyza_Basic(i,j-1,k,s) + xyza_Var(i,j-1,k,s) ) &
                &           * xyza_QSUMPN(i,j-1,k,s)                               &
                &           * xyz_Dens(i,j-1,k)                                    &
                &        + MIN( 0.0d0, xyza_Basic(i,j+1,k,s) + xyza_Var(i,j+1,k,s) ) &
                &           * xyza_QSUMPN(i,j+1,k,s)                               &
                &           * xyz_Dens(i,j+1,k)                                    &  
                &        + MIN( 0.0d0, xyza_Basic(i,j,k-1,s) + xyza_Var(i,j,k-1,s) ) &
                &           * xyza_QSUMPN(i,j,k-1,s)                               &
                &           * xyz_Dens(i,j,k-1)                                    &
                &        + MIN( 0.0d0, xyza_Basic(i,j,k+1,s) + xyza_Var(i,j,k+1,s) ) &
                &           * xyza_QSUMPN(i,j,k+1,s)                               &
                &           * xyz_Dens(i,j,k+1)                                    &  
                &       ) * 0.75d0                                                 &
                &      + ( MIN( 0.0d0, xyza_Basic(i-2,j,k,s) + xyza_Var(i-2,j,k,s) ) &
                &           * xyza_QSUMPN(i-2,j,k,s)                               &
                &           * xyz_Dens(i-2,j,k)                                    &
                &        + MIN( 0.0d0, xyza_Basic(i+2,j,k,s) + xyza_Var(i+2,j,k,s) ) &
                &           * xyza_QSUMPN(i+2,j,k,s)                               &
                &           * xyz_Dens(i+2,j,k)                                    &
                &        + MIN( 0.0d0, xyza_Basic(i,j-2,k,s) + xyza_Var(i,j-2,k,s) ) &
                &           * xyza_QSUMPN(i,j-2,k,s)                               &
                &           * xyz_Dens(i,j-2,k)                                    &
                &        + MIN( 0.0d0, xyza_Basic(i,j+2,k,s) + xyza_Var(i,j+2,k,s) ) &
                &           * xyza_QSUMPN(i,j+2,k,s)                               &
                &           * xyz_Dens(i,j+2,k)                                    &
                &        + MIN( 0.0d0, xyza_Basic(i,j,k-2,s) + xyza_Var(i,j,k-2,s) ) &
                &           * xyza_QSUMPN(i,j,k-2,s)                               &
                &           * xyz_Dens(i,j,k-2)                                    &
                &        + MIN( 0.0d0, xyza_Basic(i,j,k+2,s) + xyza_Var(i,j,k+2,s) ) &
                &           * xyza_QSUMPN(i,j,k+2,s)                                &
                &           * xyz_Dens(i,j,k+2)                                    &
                &       ) * 0.25d0                                                 &
                &     )
            end do
          end do
        end do
      end do
      
      !出力
      xyza_FillNegative_xyza = xyza_Var + xyza_DQFILL
    
    end if

  end function xyza_FillNegative_xyza

end module FillNegative
