copyright (C) 2001  MSC-RPN COMM  %%%MC2%%%
*

      subroutine rhs3_qp (sol, rhs, s, q3, ww, div, invs1, nx, ny, 1,1
     $                      lminx, lmaxx, lminy, lmaxy, lnk, niter)
      implicit none
*
      integer nx, ny, lminx, lmaxx, lminy, lmaxy, lnk, niter
      real*8 sol (nx, ny, *), rhs (nx, ny, *),invs1(nx, ny, *),  
     $	     div (lminx:lmaxx,lminy:lmaxy,lnk),
     $        q3 (lminx:lmaxx,lminy:lmaxy,lnk)
      real    ww (lminx:lmaxx,lminy:lmaxy,lnk)
      real   s(1-(niter-1):nx+(niter-1),1-(niter-1):ny+(niter-1))
*
*AUTHORs    C. Girard & M. Desgagne
*
*OBJECT
*
*    *******************************************************************
*    *                                                                 *
*    *     CALCULATES                                                  *
*    *                modified right-hand-side terms                   *
*    *                Qw*, Qp*                                         *
*    *       FROM                                                      *
*    *                right-hand-side terms:                           *
*    *                Qu, Qv, Qw, Qb, Qp                               *
*    *                stored in                                        *
*    *                up, vp, wp, bp, pp                               *
*    *      STORES                                                     *
*    *                results in                                       *
*    *                 wp=Qw*, rhs=Qp*                                 *
*    *                                                                 *
*    *     PERFORMS   essentially the following calculations:          *
*    *                                                                 *
*    *                                                                 *
*    *      wp =   [  wp + dtp * bp /(1+alfa) ] * nu0b                 *
*    *                                                     _Z          *
*    *     rhs = - pp + dtp * div(up,vp,wp) - dtp * gc02 * w           *
*    *                                                                 *
*    *                                                                 *
*    *     SEE  - divcal for calculation of DIV                        *
*    *                                                  -1             *
*    *          - nacbar for calculation of ap1r=(1+alfa) , nu0b, gc02 *
*    *                                                                 *
*    *                                                                 *
*    *******************************************************************
*
#include "dynmem.cdk"
#include "physcom.cdk"
#include "yomdyn.cdk"
#include "yomdyn1.cdk"
#include "dynpar.cdk"
#include "dtmdtp.cdk"
#include "nbcpu.cdk"
#include "alfmem.cdk"
*
      integer i, j, k, km1, kp1, nkt
      real*8  pt5, one, coe, hdiv, con, t1, t2
      parameter (pt5 = 0.5d0, one = 1.0d0)
*
*-----------------------------------------------------------------------
*
      con  =  dtp * dtp
      nkt = gnk-1
      if(flextop) nkt = gnk
      ww = 0.
*
!$omp do
      do k = 1, gnk
*
         do j = 1, ny
         do i = 1, nx
            wwp(i,j,k) = ( wwp(i,j,k) + dtp*bbp(i,j,k)*ap1r(i,j,k) )
     $                                               * nu0b(i,j,k)
            ww (i,j,k) = wwp(i,j,k)
         end do
         end do
         if (k.eq.1) then
         if (topo_folwing) then
            do j = 1, ny
            do i = 1, nx
               t1 = sbxy(i,j) * pt5 *
     $              (-gg1(i+1,j,1)*uup(i+1,j,1) - gg1(i,j,1)*uup(i,j,1)
     $              - gg2(i,j+1,1)*vvp(i,j+1,1) - gg2(i,j,1)*vvp(i,j,1) )
               ww(i,j,1) = dhdt(i,j,1) - t1
            end do
            end do
            else
            do j = 1, ny
            do i = 1, nx
               ww(i,j,1) = 0.0
            end do
            end do
            endif
            do j = 1-(niter-1), ny+(niter-1)
            do i = 1-(niter-1), nx+(niter-1)
               s(i,j) = con * sbxy(i,j)
            end do
            end do
         endif
      end do
!$omp enddo
*
      if (.not.flextop) then
!$omp do
         do j = 1, ny
         do i = 1, nx
            ww(i,j,gnk) = 0.
         end do
         end do
!$omp enddo
      endif
*
      call divcal2 (div,uup,vvp,ww,sbxy,gg1,gg2,gg0r,g0ur,g0vr,q3,
     $                  odx,ody,minx,maxx,miny,maxy,gnk,1,1,nx,ny)
*
!$omp do
      do k = 1, nkt
         kp1 = min( k+1 , gnk )
         do j = 1, ny
         do i = 1, nx
            rhs(i,j,k) = - ppp(i,j,k) + dtp * div(i,j,k)
     $                   - dtp*pt5*gc02(i,j,k)*(ww(i,j,kp1)+ww(i,j,k))
            sol(i,j,k) = pp0(i,j,k)
         end do
         end do
      end do
!$omp enddo
*
*-----------------------------------------------------------------------
      return
      end