C*********************************************************************
C*********************************************************************
C
      SUBROUTINE CLWADV
     I   ( UM1, VM1, WM1,
     I     UM2, VM2, WM2,
     O     DUADV, DVADV, DWADV,
     W     DIVMU, DIVMW,
     B     DENS0, FDZ, FDZM,
     G     DX,
     D     NX, NZ, NXB, NZB      )
C
C*********************************************************************
      IMPLICIT REAL*8 ( A-H, O-Z )
C*********************************************************************
C
      DIMENSION
     \   UM1 ( -2:NXB, -2:NZB ), UM2 ( -2:NXB, -2:NZB ),
     \   VM1 ( -2:NXB, -2:NZB ), VM2 ( -2:NXB, -2:NZB ),
     \   WM1 ( -2:NXB, -2:NZB ), WM2 ( -2:NXB, -2:NZB ),
     \   DUADV( -2:NXB, -2:NZB ),
     \   DVADV( -2:NXB, -2:NZB ),
     \   DWADV( -2:NXB, -2:NZB ),
     \   DIVMU( -2:NXB, -2:NZB ), DIVMW( -2:NXB, -2:NZB )
C
      DIMENSION
     \   DENS0(-2:NZB), FDZ(-2:NZB), FDZM(-2:NZB)
C
      DIMENSION
     \   DENSM(-2:1000)
C
C=====================================================================
C
      DO 10 IZ = 0, NZ
         DENSM(IZ) = 0.5D0 * ( DENS0(IZ) + DENS0(IZ-1) )
   10 CONTINUE
C
C
      DO 11 IZ = 0, NZ
      DO 11 IX = 1, NX
         DIVMW(IX,IZ) =
     \      ( UM1(IX+1,IZ) - UM1(IX,IZ) )
     \    + (   DENS0(IZ  ) * WM2(IX,IZ  )
     \        - DENS0(IZ-1) * WM2(IX,IZ-1) )
     \               / FDZM(IZ) / DENSM(IZ)
C        DIVMW(IX,IZ) =
C    \      ( UM1(IX+1,IZ) - UM1(IX,IZ) )
C    \               * FDZM(IZ) * DENSM(IZ)
C    \    + (   DENS0(IZ  ) * WM2(IX,IZ  )
C    \        - DENS0(IZ-1) * WM2(IX,IZ-1) )
   11 CONTINUE
C
      DO 12 IZ = 0, NZ-1
      DO 12 IX = 1, NX
         DIVMU(IX,IZ) =
     \      ( UM2(IX,IZ) - UM2(IX-1,IZ) )
     \    + (   DENSM(IZ+1) * WM1(IX,IZ+1)
     \        - DENSM(IZ  ) * WM1(IX,IZ  ) )
     \               / FDZ (IZ) / DENS0(IZ)
C        DIVMU(IX,IZ) =
C    \      ( UM2(IX,IZ) - UM2(IX-1,IZ) )
C    \               * FDZ (IZ) * DENS0(IZ)
C    \    + (   DENSM(IZ+1) * WM1(IX,IZ+1)
C    \        - DENSM(IZ  ) * WM1(IX,IZ  ) )
   12 CONTINUE
C
      CALL BOUND ( DIVMW, NX, NZ, NXB, NZB, 1, 1, 0, 1 )
      CALL BOUND ( DIVMU, NX, NZ, NXB, NZB, 1, 1, 1, 1 )
C
C--------------------------------------------------------------------
C
      DO 21 IZ = 0, NZ-1
      DO 21 IX = 1, NX
         DUADV(IX,IZ) =
     \  - ( ( UM2(IX  ,IZ) * (  UM2(IX  ,IZ) - DIVMU(IX,IZ)/2.0D0 )
     \      - UM2(IX-1,IZ) * (  UM2(IX-1,IZ) + DIVMU(IX,IZ)/2.0D0 ) )
     \     + (  UM1(IX,IZ+1) * WM1(IX,IZ+1) * DENSM(IZ+1)
     \        - UM1(IX,IZ  ) * WM1(IX,IZ  ) * DENSM(IZ  ) )
     \        / DENS0(IZ) / FDZ(IZ)                          ) / DX
C        DUADV(IX,IZ) =
C    \  - ( ( UM2(IX,IZ)**2 - UM2(IX-1,IZ)**2
C    \      -   DIVMU(IX,IZ)/4.0D0/ DENS0(IZ) / FDZ(IZ)
C    \        * ( UM2(IX,IZ) + UM2(IX-1,IZ) )             )
C    \     + (  UM1(IX,IZ+1) * WM1(IX,IZ+1) * DENSM(IZ+1)
C    \        - UM1(IX,IZ  ) * WM1(IX,IZ  ) * DENSM(IZ  )
C    \        -   DIVMU(IX,IZ)/4.0D0
C    \          * ( UM1(IX,IZ+1) + UM1(IX,IZ) )           )
C    \        / DENS0(IZ) / FDZ(IZ)                          ) / DX
   21 CONTINUE
C
      DO 22 IZ = 0, NZ
      DO 22 IX = 1, NX
         DVADV(IX,IZ) =
     \  - ( ( VM2(IX  ,IZ) * (  UM2(IX  ,IZ) - DIVMU(IX,IZ)/2.0D0 )
     \      - VM2(IX-1,IZ) * (  UM2(IX-1,IZ) + DIVMU(IX,IZ)/2.0D0 ) )
     \     + (  VM1(IX,IZ+1) * WM1(IX,IZ+1) * DENSM(IZ+1)
     \        - VM1(IX,IZ  ) * WM1(IX,IZ  ) * DENSM(IZ  ) )
     \        / DENS0(IZ) / FDZ(IZ)                          ) / DX
C        DVADV(IX,IZ) =
C    \  - ( ( VM2(IX,IZ)*UM2(IX,IZ) - VM2(IX-1,IZ)*UM2(IX-1,IZ)
C    \      -   DIVMU(IX,IZ)/4.0D0/ DENS0(IZ) / FDZ(IZ)
C    \        * ( VM2(IX,IZ) + VM2(IX-1,IZ) )             )
C    \     + (  VM1(IX,IZ+1) * WM1(IX,IZ+1) * DENSM(IZ+1)
C    \        - VM1(IX,IZ  ) * WM1(IX,IZ  ) * DENSM(IZ  )
C    \        -   DIVMU(IX,IZ)/4.0D0
C    \          * ( VM1(IX,IZ+1) + VM1(IX,IZ) )           )
C    \        / DENS0(IZ) / FDZ(IZ)                          ) / DX
   22 CONTINUE
C
      DO 23 IZ = 0, NZ
      DO 23 IX = 1, NX
         DWADV(IX,IZ) =
     \  - ( (  WM1(IX+1,IZ) * ( UM1(IX+1,IZ) - DIVMW(IX,IZ)/2.0D0 )
     \       - WM1(IX  ,IZ) * ( UM1(IX  ,IZ) + DIVMW(IX,IZ)/2.0D0 ) )
     \     + (  WM2(IX,IZ  ) * WM2(IX,IZ  ) * DENS0(IZ  )
     \        - WM2(IX,IZ-1) * WM2(IX,IZ-1) * DENS0(IZ-1) )
     \        / DENSM(IZ) / FDZM(IZ)                         ) / DX
C        DWADV(IX,IZ) =
C    \  - ( ( WM1(IX+1,IZ)*UM1(IX+1,IZ) - WM1(IX,IZ)*UM1(IX,IZ)
C    \      -   DIVMW(IX,IZ)/4.0D0/ DENSM(IZ) / FDZM(IZ)
C    \        * ( WM1(IX+1,IZ) + WM1(IX,IZ) )             )
C    \     + (  WM2(IX,IZ  )**2 * DENS0(IZ  )
C    \        - WM2(IX,IZ-1)**2 * DENS0(IZ-1)
C    \        -   DIVMW(IX,IZ)/4.0D0
C    \          * ( WM2(IX,IZ) + WM2(IX,IZ-1) )           )
C    \        / DENSM(IZ) / FDZM(IZ)                          ) / DX
   23 CONTINUE
C
C
      CALL BOUND ( DUADV, NX, NZ, NXB, NZB, 1, 1, 1, 1 )
      CALL BOUND ( DVADV, NX, NZ, NXB, NZB, 1, 1, 1, 1 )
      CALL BOUND ( DWADV, NX, NZ, NXB, NZB, 1,-1, 0, 1 )
C
C===================================================================
C
      RETURN
      E N D
