c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE RDMPUVW                    ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c
      SUBROUTINE RDMPUVW(nx,ny,nz,
     :           u,v,w, ubar,vbar,rhobar, zp,
     :           umix,vmix,wmix, rdmp)
c
c#######################################################################
c
c     PURPOSE:
c
c     Calculate upper level Rayleigh damping on w, and the perturbations 
c     of u and v in momentum equations. The Rayleigh damping terms are 
c     then added to arrays umix, vmix and wmix.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     4/21/1992.  
c
c     MODIFICATION HISTORY:
c
c     5/05/92 (M. Xue)
c     Added full documentation.
c
c     6/2/92 (M. Xue and H. Jin)
c     Further facelift.
C
c#######################################################################
c
c     INPUT :
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical direction.
c
c       u        x component of velocity (m/s)
c       v        y component of velocity (m/s)
c       w        Vertical component of velocity in Cartesian
c                coordinates (m/s).
c
c       ubar     Base state zonal velocity component (m/s)
c       vbar     Base state meridional velocity component (m/s)
c       rhobar   Base state density (kg/m**3)
c
c       zp       Vertical coordinate of grid points in physical space (m)
c
c       umix     Array containing turbulent and computational mixing on u 
c       vmix     Array containing turbulent and computational mixing on v 
c       wmix     Array containing turbulent and computational mixing on w 
c
c     OUTPUT:
c
c       umix     Mixing including rayleigh damping on u
c       vmix     Mixing including rayleigh damping on v
c       wmix     Mixing including rayleigh damping on w
c
c     WORK ARRAYS:
c
c       rdmp     Rayleigh damping coefficient.
c                A temporary array defined locally.
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none 
c
      integer nx,ny,nz          ! Number of grid points in 3 directions
c 
      real u     (nx,ny,nz)     ! Total u-velocity (m/s)
      real v     (nx,ny,nz)     ! Total v-velocity (m/s)
      real w     (nx,ny,nz)     ! Total w-velocity (m/s)

      real ubar  (nx,ny,nz)     ! Base state u-velocity (m/s)
      real vbar  (nx,ny,nz)     ! Base state v-velocity (m/s)
      real rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
c
      real zp    (nx,ny,nz)     ! The physical height coordinate defined at
                                ! w-point of the staggered grid.
c
      real umix  (nx,ny,nz)     ! Total mixing in u-eq.
      real vmix  (nx,ny,nz)     ! Total mixing in v-eq.
      real wmix  (nx,ny,nz)     ! Total mixing in w-eq.
c
      real rdmp  (nx,ny,nz)     ! Rayleigh damping coefficient.
                                ! Temporary array defined locally.
c
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i, j, k
      integer klowest
      real    pi, zpmax
c
c#######################################################################
c
c     Include files: 
c
c#######################################################################
c

      include 'globcst.inc'

c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      IF( raydmp.eq.0 .or. cfrdmp.eq.0.0 ) RETURN
c
c#######################################################################
c
c     Find the lowest model layer (index klowest) that is entirely or 
c     partially contained in the Rayleigh damping (sponge) layers.
c
c     The Rayleigh damping is then applied only to layers with k>=klowest.
c
c#######################################################################
c
      klowest = nz-1

      DO 100 k=nz-1,2,-1

        zpmax = zp(1,1,k)

        DO 110 j=1,ny-1 
        DO 110 i=1,nx-1 
          zpmax = max( zp(i,j,k), zpmax )
110     CONTINUE

        IF( zpmax.lt.zbrdmp ) THEN 

          klowest = max(2, k+1)
          GOTO 120

        ENDIF

100   CONTINUE

120   CONTINUE

c
c#######################################################################
c
c     Calculate the Rayleigh damping coefficient defined at w-point.
c
c#######################################################################
c
      pi = 4.0*atan( 1.0 )

      DO 200 k= klowest,nz
      DO 200 j=1,ny-1
      DO 200 i=1,nx-1

        rdmp(i,j,k) = cfrdmp* 
     :     (1.0-cos(pi*(zp(i,j,k)-zbrdmp)/(ztrdmp-zbrdmp)))*0.5

200   CONTINUE
c
c#######################################################################
c
c     Calculate and add to umix the Rayleigh damping term in u-equation.
c
c#######################################################################
c
      DO 300 k=klowest,nz-2
      DO 300 j=2,ny-2
      DO 300 i=2,nx-1

        umix(i,j,k)=umix(i,j,k) - 
     :  (rdmp(i,j,k)+rdmp(i-1,j,k)+rdmp(i,j,k+1)+rdmp(i-1,j,k+1))*0.25*
     :  (rhobar(i,j,k)+rhobar(i-1,j,k))*0.5*
     :  (u(i,j,k)-ubar(i,j,k))

300   CONTINUE

c
c#######################################################################
c
c     Calculate and add to vmix the Rayleigh damping term in v-equation.
c
c#######################################################################
c
      DO 400 k=klowest,nz-2
      DO 400 j=2,ny-1
      DO 400 i=2,nx-2

        vmix(i,j,k)=vmix(i,j,k) - 
     :  (rdmp(i,j,k)+rdmp(i,j-1,k)+rdmp(i,j,k+1)+rdmp(i,j-1,k+1))*0.25*
     :  (rhobar(i,j,k)+rhobar(i,j-1,k))*0.5*
     :  (v(i,j,k)-vbar(i,j,k))

400   CONTINUE

c
c#######################################################################
c
c     Calculate and add to wmix the Rayleigh damping term in w-equation.
c
c#######################################################################
c
      DO 500 k=klowest,nz-1
      DO 500 j=2,ny-2
      DO 500 i=2,nx-2

        wmix(i,j,k)=wmix(i,j,k) - rdmp(i,j,k)*
     :              (rhobar(i,j,k)+rhobar(i,j,k-1))*0.5*
     :              w(i,j,k)

500   CONTINUE

      RETURN
      END

C
c
c     ##################################################################
c     ##################################################################
c     ######                                                      ######
c     ######                SUBROUTINE RDMPPT                     ######
c     ######                                                      ######
c     ######                Copyright (c) 1992                    ######
c     ######  Center for the Analysis and Prediction of Storms    ######
c     ######    University of Oklahoma.  All rights reserved.     ######
c     ######                                                      ######
c     ##################################################################
c     ##################################################################
c
      SUBROUTINE RDMPPT(nx,ny,nz,
     :                  ptprt ,rhobar, zp,
     :                  ptmix, rdmp)
c
c#######################################################################
c
c     PURPOSE:
c
c     Calculate upper level Rayleigh damping on potential temperature
c     perturbations. The term is accumulated in array ptmix.
c
c#######################################################################
c
c     AUTHOR: Ming Xue
c     4/21/1992.  
c
c     MODIFICATION HISTORY:
c
c     5/05/92 (M. Xue)
c     Added full documentation.
c
c     6/2/92 (M. Xue and H. Jin)
c     Further facelift.
C
c#######################################################################
c
c     INPUT :
c
c       nx       Number of grid points in the x-direction (east/west)
c       ny       Number of grid points in the y-direction (north/south)
c       nz       Number of grid points in the vertical direction.
c
c       ptprt    Perturbation potential temperature (K)
c
c       rhobar   Base state air density (kg/m**3)
c
c       zp       Vertical coordinate of grid points in physical space (m)
c
c     OUTPUT:
c
c       ptmix    Total mixing in potential temperature equation.
c
c
c     WORK ARRAYS:
c
c       rdmp     Rayleigh damping coefficient.
c                A temporary array defined locally.
c
c#######################################################################
c

c
c#######################################################################
c
c     Variable Declarations.
c
c#######################################################################
c
      implicit none
c
      integer nx,ny,nz          ! Number of grid points in 3 directions
c 
      real ptprt (nx,ny,nz)     ! Perturbation potential temperature (K)
      real rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
      real zp    (nx,ny,nz)     ! The physical height coordinate defined at
                                ! w-point of the staggered grid.
c
      real ptmix (nx,ny,nz)     ! Total mixing on potential temperature
c
      real rdmp  (nx,ny,nz)     ! Rayleigh damping coefficient.
c
c#######################################################################
c
c     Misc. local variables:
c
c#######################################################################
c
      integer i, j, k
      integer klowest
      real    pi, zpmax
c
c#######################################################################
c
c     Include files: 
c
c#######################################################################
c

      include 'globcst.inc'

c
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
C
C     Beginning of executable code...
C
C@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
c
      IF( raydmp.eq.0 .or. cfrdmp.eq.0.0 ) RETURN
c
c#######################################################################
c
c     Find the lowest model layer (index klowest) that is entirely or
c     partly contained in the Rayleigh damping (sponge) layers.
c
c     The Rayleigh damping is then only applied to layers with k>=klowest.
c
c#######################################################################
c
      klowest = nz-1

      DO 100 k=nz-1,2,-1

        zpmax = zp(1,1,k)

        DO 110 j=1,ny-1
        DO 110 i=1,nx-1
          zpmax = max( zp(i,j,k), zpmax )
110     CONTINUE

        IF( zpmax.lt.zbrdmp ) THEN

          klowest = max(2, k+1)
          GOTO 120

        ENDIF

100   CONTINUE

120   CONTINUE

c
c#######################################################################
c
c     Calculate the Rayleigh damping coefficient defined at w-point.
c
c#######################################################################
c
      pi = 4.0*atan( 1.0 )

      DO 200 k= klowest,nz
      DO 200 j=1,ny-1
      DO 200 i=1,nx-1

        rdmp(i,j,k) = cfrdmp*
     :     (1.0-cos(pi*(zp(i,j,k)-zbrdmp)/(ztrdmp-zbrdmp)))*0.5

200   CONTINUE
c
c#######################################################################
c
c     Calculate the Rayleigh damping and add it to the ptmix term
c     in pt-equation.
c
c#######################################################################
c
      DO 300 k=klowest,nz-2
      DO 300 j=2,ny-2
      DO 300 i=2,nx-2

        ptmix(i,j,k)=ptmix(i,j,k) - 
     :              (rdmp(i,j,k)+rdmp(i,j,k+1))*0.5*
     :               rhobar(i,j,k)*ptprt(i,j,k)

300   CONTINUE


      RETURN
      END
