C***********************************************************
C***********************************************************
C
      SUBROUTINE  NLDIFV
     I   (   QADV, SCALE, zscale, dt, dx,
     I       IDF,
     O       DTQDNL,
     W       QLAP, CDQNLX, CDQNLZ,
     B       zgrid, 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 ),
     \     DTQDNL( -2:NXB, -2:NZB ),
     \     QLAP  ( -2:NXB, -2:NZB ),
     \     CDQNLX( -2:NXB, -2:NZB ),
     \     CDQNLZ( -2:NXB, -2:NZB ),
     \     DENS0 (-2:NZB), FDZ(-2:NZB),
     \     ZGRID (-2:NZB)
C
C=====================================================================
C     COMPUTATIONAL DIFFUSION TERM ( ARTIFICIAL SMOOTHING )
C --------------------------------------------------------------------
C
C *** Calculation of pseudo Laplacian ***
C
      DO 11 IZ = 0, NZ-1
         sclap = scale * exp ( zgrid(iz) * zscale )
      DO 11 IX = 1, NX

         QLAP(IX,IZ) =
     \ (3.*ABS(QADV(IX+1,IZ,IDF)+QADV(IX-1,IZ,IDF)-2.*QADV(IX,IZ,IDF))
     \ +ABS(QADV(IX,IZ+1,IDF)+QADV(IX,IZ-1,IDF)-2.*QADV(IX,IZ,IDF)) )
     \ / sclap

   11 CONTINUE
C
      CALL BOUND ( QLAP, NX, NZ, NXB, NZB, 1, 1, 1, 1 )
C
C *** Calculation of Diffusivity ***
C
      CD00 =    1.0D0 
      CDMAX =   ( dx**2 / dt ) / 10.0d0
C
      DO 21 IZ = 0, NZ
      DO 21 IX = 1, NX
         CDQNLX(IX,IZ) = MIN ( CDMAX, (
     \      CD00 * ( ( QLAP(IX,IZ) + QLAP(IX-1,IZ) ) / 2. )  ) )
         CDQNLZ(IX,IZ) = MIN ( CDMAX, (
     \      CD00 * ( ( QLAP(IX,IZ) + QLAP(IX,IZ-1) ) / 2. )  ) )
   21 CONTINUE
C
      CALL BOUND( CDQNLX, NX, NZ, NXB, NZB, 1, 1, 1, 1 )
      CALL BOUND( CDQNLZ, NX, NZ, NXB, NZB, 1, 1, 1, 1 )
C
C *** Time derivatives due to Non-linear Diffusion ***
C
      DO 31 IZ = 0, NZ-1
      DO 31 IX = 1, NX

         DQDINX = 
     \   - CDQNLX(IX  ,IZ) * ( QADV(IX  ,IZ,IDF) - QADV(IX-1,IZ,IDF) )
     \   + CDQNLX(IX+1,IZ) * ( QADV(IX+1,IZ,IDF) - QADV(IX  ,IZ,IDF) )

C
         DQDINZ =
     \    (  ( DENS0(IZ+1) + DENS0(IZ) ) / 2.0D0 
     \         * CDQNLZ(IX,IZ+1) 
     \         * ( QADV(IX,IZ+1,IDF) - QADV(IX,IZ,IDF) )
     \     - ( DENS0(IZ-1) + DENS0(IZ) ) / 2.0D0 
     \         * CDQNLZ(IX,IZ  ) 
     \         * ( QADV(IX,IZ,IDF) - QADV(IX,IZ-1,IDF) ) )
     \    / DENS0(IZ) / FDZ(IZ)

         DTQDNL(IX,IZ) = ( dqdinx + dqdinz ) / dx ** 2

   31 CONTINUE
C
      CALL BOUND( DTQDNL, NX, NZ, NXB, NZB, 1, 1, 1, 1 )
C
C================================================================
C
      RETURN
      E N D
