copyright (C) 2001  MSC-RPN COMM  %%%MC2%%%
***s/r divcal2

      subroutine divcal2 ( div,uu,vv,ww,sbxy,gg1,gg2,gg0r,g0ur,g0vr,q3, 3
     $                 odx,ody,lminx,lmaxx,lminy,lmaxy,lnk,id,jd,iff,jf)
      implicit none
*
      integer   lminx,lmaxx,lminy,lmaxy,lnk,id,jd,iff,jf
      real*8 div (lminx:lmaxx,lminy:lmaxy,lnk),
     $        q3 (lminx:lmaxx,lminy:lmaxy,lnk),odx(*),ody(*)
      real uu  (lminx:lmaxx,lminy:lmaxy,lnk),
     $     vv  (lminx:lmaxx,lminy:lmaxy,lnk),
     $     ww  (lminx:lmaxx,lminy:lmaxy,lnk),
     $     sbxy(lminx:lmaxx,lminy:lmaxy),
     $     gg1 (lminx:lmaxx,lminy:lmaxy,lnk),
     $     gg2 (lminx:lmaxx,lminy:lmaxy,lnk),
     $     gg0r(lminx:lmaxx,lminy:lmaxy,lnk),
     $     g0ur(lminx:lmaxx,lminy:lmaxy,lnk),
     $     g0vr(lminx:lmaxx,lminy:lmaxy,lnk)
*
*AUTHORs    C. Girard & M. Desgagne
*
*OBJECT 
*
*    *******************************************************************
*    *                                                                 *
*    *                                                                 *
*    *             CALCULATE GENERALIZED DIVERGENCE                    *
*    *                                                                 *
*    *                                                                 *
*    *    A) in terrain-following coordinate                           *
*    *                                                                 *
*    *                           ~       ~       ~                     *
*    *           div(u,v,w)   =  du/dX + dv/dY + dw/dZ                 *
*    *                                                                 *
*    *                                                                 *
*    *         *******************************************             *
*    *         *    ~                        ________XZ  *             *
*    *         *    du/dX =      du/dX + G0r*G1*du/dZ    *             *
*    *         *                                         *             *
*    *         *    ~                        ________YZ  *             *
*    *         *    dv/dY =      dv/dY + G0r*G2*dv/dZ    *             *
*    *         *                                         *             *
*    *         *    ~                                    *             *
*    *         *    dw/dZ = G0r* dw/dZ                   *             *
*    *         *                                         *             *
*    *         *******************************************             *
*    *                                                                 *
*    *                                                                 *
*    *    B) in fixed coordinate                                       *
*    *                                                                 *
*    *                             u            v            w         *
*    *     div(u,v,w) = G0r * [ d{G *u}/dX + d{G *v}/dY + d{G *w}/dZ ] *
*    *                             0            0            0         *
*    *                                                                 *
*    *                                                                 *
*    *******************************************************************
*
#include "lcldim.cdk"
#include "nbcpu.cdk"
#include "yomdyn.cdk"
*
      integer i,j,k,km1
      real*8 pt5
      parameter (pt5=0.5d0)
*
*-----------------------------------------------------------------------
*
      if(topo_folwing) then
!$omp do
      do k = gnk, 1, -1
         if (k.eq.gnk) then
            do j = jd, jf
            do i = id, iff
                q3(i,j,gnk) = 0.0
               div(i,j,gnk) = 0.0
            end do
            end do
         else
            km1 = max( k-1 , 1 )
            do j = jd, jf
            do i = id, iff
            q3(i,j,k) = sbxy(i,j) * pt5 *
     $                (-gg1(i+1,j,k) * ( uu(i+1,j,k) - uu(i+1,j,km1) )
     $                - gg1(i  ,j,k) * ( uu(i  ,j,k) - uu(i  ,j,km1) )
     $                - gg2(i,j+1,k) * ( vv(i,j+1,k) - vv(i,j+1,km1) )
     $                - gg2(i,j  ,k) * ( vv(i,j  ,k) - vv(i,j  ,km1) ) )
            div(i,j,k)= sbxy(i,j) * (( uu(i+1,j,k) - uu(i,j,k))*odx(1) +
     $                               ( vv(i,j+1,k) - vv(i,j,k))*ody(j) )
     $                + gg0r(i,j,k)* (ww(i,j,k+1) -ww(i,j,k) +
     $                               (q3(i,j,k+1) +q3(i,j,k) ) * pt5 )
            end do
            end do
         endif
      end do
!$omp enddo
*
      else
c solid body
!$omp do
      do k = 1, gnk
         if (k.eq.gnk) then
            do j = jd, jf
            do i = id, iff
               div(i,j,gnk) = 0.0
            end do
            end do
         else
            do j = jd, jf
            do i = id, iff
            div(i,j,k)= sbxy(i,j) * gg0r(i,j,k) * (
     $                ( g0ur(i+1,j,k)*uu(i+1,j,k)
     $                - g0ur(i  ,j,k)*uu(i  ,j,k))*odx(1) +
     $                ( g0vr(i,j+1,k)*vv(i,j+1,k)
     $                - g0vr(i,j  ,k)*vv(i,j  ,k))*ody(j) )
     $                + gg0r(i,j,k)* ( ww(i,j,k+1) -ww(i,j,k) ) 
            end do
            end do
         endif
      end do
!$omp enddo
      endif
*
*-----------------------------------------------------------------------
*
      return
      end