C*********************************************************************
C *** 負の雲水量などを穴埋めするルーチン ***
C*********************************************************************
C
      SUBROUTINE QNFILL
     I   ( QADV, ICL,
     O     DQFILL,
     W     QSUMPN,
     B     DENS0, FDZ,
     D     NX, NZ, NXB, NZB, NROT )
C
C*********************************************************************
      IMPLICIT REAL*8 ( A-H, O-Z )
C*********************************************************************
C
      DIMENSION QADV   ( -2:NXB, -2:NZB, NROT ),
     \          DQFILL ( -2:NXB, -2:NZB ),
     \          QSUMPN ( -2:NXB, -2:NZB )
      DIMENSION DENS0(-2:NZB), FDZ(-2:NZB)
C
C=====================================================================
C
      EPS = 1.0D-60
C
      DO 11 IZ = 0, NZ-1
      DO 11 IX = 1, NX
         QSUMPN(IX,IZ) = 1.0D0 /
     \  ( (   MAX( 0.0D0, QADV(IX-1,IZ  ,ICL) )
     \          * DENS0(IZ  ) * FDZ(IZ  )
     \      + MAX( 0.0D0, QADV(IX+1,IZ  ,ICL) )
     \          * DENS0(IZ  ) * FDZ(IZ  )
     \      + MAX( 0.0D0, QADV(IX  ,IZ-1,ICL) )
     \          * DENS0(IZ-1) * FDZ(IZ-1)
     \      + MAX( 0.0D0, QADV(IX  ,IZ+1,ICL) )
     \          * DENS0(IZ+1) * FDZ(IZ+1)       ) * 0.75D0
     \  + (   MAX( 0.0D0, QADV(IX-2,IZ  ,ICL) )
     \          * DENS0(IZ  ) * FDZ(IZ  )
     \      + MAX( 0.0D0, QADV(IX+2,IZ  ,ICL) )
     \          * DENS0(IZ  ) * FDZ(IZ  )
     \      + MAX( 0.0D0, QADV(IX  ,IZ-2,ICL) )
     \          * DENS0(IZ-2) * FDZ(IZ-2)
     \      + MAX( 0.0D0, QADV(IX  ,IZ+2,ICL) )
     \          * DENS0(IZ+2) * FDZ(IZ+2)       ) * 0.25D0 + EPS )
   11 CONTINUE
C
      CALL BOUND ( QSUMPN, NX, NZ, NXB, NZB, 1, 1, 1, 1 )
C
C
      DO 21 IZ = 0, NZ-1
      DO 21 IX = 1, NX
         DQFILL(IX,IZ) =
     \  - MIN ( 0.0D0, QADV(IX,IZ,ICL) )
     \  + MAX ( 0.0D0, QADV(IX,IZ,ICL) )
     \ * ( 0.75D0*(  MIN( 0.0D0, QADV(IX-1,IZ,ICL) ) * QSUMPN(IX-1,IZ)
     \               * DENS0(IZ  ) * FDZ(IZ  )
     \             + MIN( 0.0D0, QADV(IX+1,IZ,ICL) ) * QSUMPN(IX+1,IZ)
     \               * DENS0(IZ  ) * FDZ(IZ  )
     \             + MIN( 0.0D0, QADV(IX,IZ-1,ICL) ) * QSUMPN(IX,IZ-1)
     \               * DENS0(IZ-1) * FDZ(IZ-1)
     \             + MIN( 0.0D0, QADV(IX,IZ+1,ICL) ) * QSUMPN(IX,IZ+1)
     \               * DENS0(IZ+1) * FDZ(IZ+1)       )
     \   + 0.25D0*(  MIN( 0.0D0, QADV(IX-2,IZ,ICL) ) * QSUMPN(IX-2,IZ)
     \               * DENS0(IZ  ) * FDZ(IZ  )
     \             + MIN( 0.0D0, QADV(IX+2,IZ,ICL) ) * QSUMPN(IX+2,IZ)
     \               * DENS0(IZ  ) * FDZ(IZ  )
     \             + MIN( 0.0D0, QADV(IX,IZ-2,ICL) ) * QSUMPN(IX,IZ-2)
     \               * DENS0(IZ-2) * FDZ(IZ-2)
     \             + MIN( 0.0D0, QADV(IX,IZ+2,ICL) ) * QSUMPN(IX,IZ+2)
     \               * DENS0(IZ+2) * FDZ(IZ+2)       )              )
   21 CONTINUE
C
C
      CALL BOUND ( DQFILL, NX, NZ, NXB, NZB, 1, 1, 1, 1 )
C
      DO 31 IZ = 0, NZ-1
      DO 31 IX = 1, NX
         QADV(IX,IZ,ICL) = QADV(IX,IZ,ICL) + DQFILL(IX,IZ)
   31 CONTINUE
C
C
      CALL BOUND ( QADV, NX, NZ, NXB, NZB, NROT, 1, 1, ICL )
C
C     QFSUM = 0.0D0
C
C     DO 41 IZ = 0, NZ-1
C     DO 41 IX = 1, NX
C        QFSUM = QFSUM + DQFILL(IX,IZ) * DENS0(IZ) * FDZ(IZ)
C  41 CONTINUE
C
C     WRITE(6,*) 'SUM OF FILLING =', QFSUM
C
C======================================================================
C
      RETURN
      E N D
