copyright (C) 2001  MSC-RPN COMM  %%%MC2%%%
***s/r out_qrdd -- Compute vorticity and divergence
*

      subroutine out_qrdd (qr,dd,uv,msf,lminx,lmaxx,lminy,lmaxy,lnk) 1
      implicit none
*
      integer lminx,lmaxx,lminy,lmaxy,lnk
      real qr(lminx:lmaxx,lminy:lmaxy,lnk),
     $     dd(lminx:lmaxx,lminy:lmaxy,lnk),
     $     uv(lminx:lmaxx,lminy:lmaxy,lnk,2),
     $    msf(lminx:lmaxx,lminy:lmaxy)
*
*OBJECT
*     This subroutine computes relative vorticity (qr)
*     and the divergence (dd).
*
**
#include "lcldim.cdk"
#include "yomdyn1.cdk"
*
      integer i,j,k,i0,in,j0,jn
      real c1,c2

*----------------------------------------------------------------------
      c1 = 0.25 / grdx
      c2 = 1.0  / grdx
*
      i0 = 1    - (hx-1)*west 
      j0 = 1    - (hy-1)*south
      in = ldni + (hx-1)*east
      jn = ldnj + (hy-1)*north
*
      do k=1,lnk
*
*   * Computing vorticity
*
         do j=j0,jn
         do i=i0,in
            qr(i,j,k) = msf(i,j) * c1 *
     $             ((uv(i+1,j+1,k,2)+uv(i+1,j,k,2)
     $              -uv(i-1,j+1,k,2)-uv(i-1,j,k,2))
     $             -(uv(i+1,j+1,k,1)+uv(i,j+1,k,1)
     $              -uv(i+1,j-1,k,1)-uv(i,j-1,k,1)))
         end do
         end do
*
*   * Computing divergence
*
         do j=j0-south,jn
         do i=i0-west ,in
            dd(i,j,k) = msf(i,j) * c2 *
     $                  ( (uv(i+1,j,k,1)-uv(i,j,k,1))
     $                  + (uv(i,j+1,k,2)-uv(i,j,k,2)) )
         end do
         end do
*
      end do
*
      if (.not.period_x) then
         if (west.gt.0) then
            do k=1,lnk
            do j=j0,jn
               qr(i0-1,j,k) = 2.*qr(i0,j,k) - qr(i0+1,j,k)
            end do
            end do
         endif
         if (east.gt.0) then
            do k=1,lnk
            do j=j0,jn
               qr(in+1,j,k) = 2.*qr(in,j,k) - qr(in-1,j,k)
            end do
            do j=j0-south,jn
               dd(in+1,j,k) = 2.*dd(in,j,k) - dd(in-1,j,k)
            end do
            end do
         endif
      endif
      if (.not.period_y) then
         if (south.gt.0) then
            do k=1,lnk
            do i=i0-west ,in+east
               qr(i,j0-1,k) = 2.*qr(i,j0,k) - qr(i,j0+1,k)
            end do
            end do
         endif
         if (north.gt.0) then
            do k=1,lnk
            do i=i0-west ,in+east
               qr(i,jn+1,k) = 2.*qr(i,jn,k) - qr(i,jn-1,k)
               dd(i,jn+1,k) = 2.*dd(i,jn,k) - dd(i,jn-1,k)
            end do
            end do
         endif
      endif
*----------------------------------------------------------------------
      return
      end