#include "CPP_EEOPTIONS.h"

CStartOfInterFace
      SUBROUTINE CALC_QGPV(
     I          bi, bj, iMin,iMax,jMin,jMax, pH, nsquare,
     I          K13, K23, 
     I          myTime, myIter,
     O          q, qr, qs, KpvU, KpvV, gUpv, gVpv, dqdyU,
     I          myThid )

C     /==========================================================\
C     | SUBROUTINE CALC_QGPV                                     |
C     | o Calculate the quasigeostrophic potential vorticity     |
C     | o and eddy flux terms of qgpv etc                        |
C     |==========================================================|
C     |  Richard Wardle   7/98                                   |
C     |  Daniel Jamous    1/00 N^2 entered as an argument        |
C     |                        instead of being calculated from  |
C     |                        temperature                       |
C     |                        Updated to take into account      |
C     |                        arbitrary topography              |
C     |                   3/00 compute PV fluxes using a local   |
C     |                        N^2 (as opposed to an horizontal  |
C     |                        average)--see key N2local         |
C     |  IMPORTANT NOTE: For now, the routine has been tested    |
C     |  only when there are no interpolation of gsf and dTdz    |
C     |  at horizontal boundaries, when qgpv is only made of     |
C     |  qs the stretching term, and when the Kpv are constants. |
C     |  For a more general case, some more coding might be      |
C     |  necessary.                                              |
C     \==========================================================/
C
C     == Global variables ==
#include "SIZE.h"
#include "DYNVARS.h"
#include "GRID.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "CG2D.h"
#include "FFIELDS.h"
C 
C     ========= Local variables that we may want to export =============
C     q    - quasigeostrophic potential vorticity
C     ==================================================================
C
      _RL pH         (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL q          (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL qf         (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL qr         (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL qs         (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL KpvV       (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL KpvU       (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL gUpv       (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
      _RL gVpv       (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
      _RL nsquare       (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
      _RL K13        (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL K23        (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
C
      INTEGER bi,bj,iMin,iMax,jMin,jMax
      INTEGER myThid
      INTEGER myIter
      _RL myTime
C
CEndOfInterface
C
C     define local variables here:
      _RL ptotal     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL sfn        (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL gsfn       (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL tmpphi     (lShare8,MAX_NO_THREADS)
      _RL tmpphi_area     (lShare8,MAX_NO_THREADS)
c     _RL tbarxy     (Nr)
      _RL sfnbarxy   (Nr)
      _RL dTdz       (Nr)
      _RL dTdz3d     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL dVdx       (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL dUdy       (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL dVdxbarxy  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL dUdybarxy  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL dgsfndz    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL arr        (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
      _RL arrprime   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr+1)
      _RL arr2       (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
C
      _RL dqdyV      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL dqdyT      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL dqdyU      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
C
      _RL dqdxU      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL dqdxT      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL dqdxV      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
C
      _RL pvFacT     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL pvFacU     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL pvFacV     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
C
      _RL d2gsfndz2  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL interp1    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL d2Tdz2     (Nr)
      _RL interp2    (Nr)
C
      _RL pmask      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL umask      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
      _RL vmask      (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
C
      _RL num_i     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL num_b     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL denom     (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL kbotindex (1-OLx:sNx+OLx,1-OLy:sNy+OLy)

      INTEGER i,j,k
      INTEGER iG, jG
      INTEGER ku, kv, kup1, kbot
ctest
      CHARACTER*(MAX_LEN_MBUF) suff
ctest
      _RL tmpphiS
      _RL Kpvref
      _RL tmp, fac
      _RL mldu, mldv
      _RL int1, int2, int3, int4, int5, int6
      _RL int2km1, signtest
  
C     ==================================================================
C--   iG and jG are the global indices
C     ==================================================================      
      Kpvref = 1000.
C     ==================================================================
C--   Compute mask at pressure and velocity points
C     ==================================================================

      DO k=1,Nr
       DO j=iMin,jMax
        DO i=jMin,iMax
          pmask(i,j,k) = 1.
          IF (_hFacC(i,j,k,bi,bj).eq.0.) pmask(i,j,k)=0.
        ENDDO
       ENDDO
       DO j=iMin,jMax
        DO i=jMin,iMax
          umask(i,j,k) = pmask(i-1,j,k)*pmask(i,j,k)
          vmask(i,j,k) = pmask(i,j-1,k)*pmask(i,j,k)
        ENDDO
       ENDDO
      ENDDO

      DO j=iMin,jMax
       DO i=jMin,iMax
         kbotindex(i,j) = 0.
         IF ( pmask(i,j,Nr) .eq. 1. ) THEN
           kbotindex(i,j) = float(Nr)
         ELSE
           DO K = Nr-1,1,-1
             IF ( pmask(i,j,k+1) .EQ. 0. .AND. 
     &            pmask(i,j,k) .EQ. 1. ) THEN
               kbotindex(i,j) = float(k)
             ENDIF
           ENDDO
         ENDIF
       ENDDO
      ENDDO

C
C     ==================================================================
C     Compute the geostrophic streamfunction: gsf ======================
C     ==================================================================
C
C       -------v--------
C          |       |
C          |       |
C          u   x   u       Streamfunction is located at p points
C          |  gsfn |
C          |       |
C       -------v--------
C
C     ==================================================================
C     Compute the total pressure field:
      DO k=1,Nr
       DO j=iMin,jMax
        DO i=jMin,iMax
          ptotal(i,j,k) = pH(i,j,k) +
     &                    cg2d_x(i,j,bi,bj) * (gBaro * rhonil)
        ENDDO
       ENDDO
      ENDDO
C     Compute the streamfunction:
      DO k=1,Nr
       DO j=iMin,jMax
        DO i=iMin,iMax
           sfn(i,j,k) = ptotal(i,j,k) / ( rhonil * f0 )
        ENDDO
       ENDDO
      ENDDO
C
C     Compute the streamfunction:
      DO k=1,Nr
       DO j=iMin,jMax
        DO i=iMin,iMax
           sfn(i,j,k) = ptotal(i,j,k) / ( rhonil * f0 )
        ENDDO
       ENDDO
      ENDDO
C     ==================================================================
C     1. Evaluate the global horizontal mean of sf:
c     write(0,*)'Evaluate the horizontal mean'
      do k=1,Nr
       tmpphi(1,myThid)=0.
       tmpphi_area(1,myThid)=0.
       do j=1,sNy
       do i=1,sNx
          tmpphi(1,myThid) = tmpphi(1,myThid) + sfn(i,j,k)*
     &                       _rA(i,j,bi,bj)*pmask(i,j,k)
          tmpphi_area(1,myThid) = tmpphi_area(1,myThid) + 
     &                       _rA(i,j,bi,bj)*pmask(i,j,k)
       enddo
       enddo
       _GLOBAL_SUM_R8( tmpphi, myThid )
       _GLOBAL_SUM_R8( tmpphi_area, myThid )
        sfnbarxy(k)=tmpphi(1,myThid)/tmpphi_area(1,myThid)
c      write(0,*)'level :',k,'  horizontal mean :',sfnbarxy(k)
      enddo
C
C     2. Subtract sgbarxy from sf to give the geostrophic streamfunction:
       DO k=1,Nr
          DO j=jMin,jMax
            DO i=iMin,iMax
              gsfn(i,j,k) = ( sfn(i,j,k) - sfnbarxy(k) ) * pmask(i,j,k) 
            ENDDO
          ENDDO
       ENDDO 
C     ==================================================================
C--   Interpolate gsf at horizontal boundaries:
C     ==================================================================
C--   Vertical gradient of gsfn:
      DO k=2,Nr
       DO j=jMin,jMax
        DO i=iMin,iMax       
          dgsfndz(i,j,k) = recip_drC(k)*( gsfn(i,j,k-1)-gsfn(i,j,k) )
        ENDDO
       ENDDO
      ENDDO      
C--   Vertical gradient of dgsfndz:
      DO k=2,Nr
       DO j=jMin,jMax
        DO i=iMin,iMax       
          d2gsfndz2(i,j,k) = recip_drC(k)*(dgsfndz(i,j,k)-dgsfndz(i,j,k+1))
        ENDDO
       ENDDO
      ENDDO 
C--
      DO k=2,3
       DO j=jMin,jMax
        DO i=iMin,iMax       
          d2gsfndz2(i,j,k) = d2gsfndz2(i,j,4)
        ENDDO
       ENDDO
      ENDDO       
C
C--   Evaluate the gradients of gsfn
       DO j=jMin,jMax
        DO i=iMin,iMax  
          interp1(i,j,3) = dgsfndz(i,j,4) + drC(3)*d2gsfndz2(i,j,3)
        ENDDO
       ENDDO
C
       DO j=jMin,jMax
        DO i=iMin,iMax  
          interp1(i,j,2) = interp1(i,j,3) + drC(2)*d2gsfndz2(i,j,2)
        ENDDO
       ENDDO
      DO k=4,Nr
       DO j=jMin,jMax
        DO i=iMin,iMax       
          interp1(i,j,k) = dgsfndz(i,j,k)
        ENDDO
       ENDDO
      ENDDO  
C     ==================================================================
C     ==================================================================
C     ===  Global mean temperature profile: dTdz3d =====================
C     ==================================================================
C     Evaluate the global horizontal mean profile of nsquare   
c     write(0,*)'mean N2'
      do k=1,Nr
       tmpphi(1,myThid)=0.
       tmpphi_area(1,myThid)=0.
       do j=1,sNy
        do i=1,sNx
         tmpphi(1,myThid)=tmpphi(1,myThid)+nsquare(i,j,k,bi,bj)*
     &                       _rA(i,j,bi,bj)*pmask(i,j,k)
         tmpphi_area(1,myThid)=tmpphi_area(1,myThid) +
     &                       _rA(i,j,bi,bj)*pmask(i,j,k)
        enddo
       enddo
       _GLOBAL_SUM_R8( tmpphi, myThid )
       _GLOBAL_SUM_R8( tmpphi_area, myThid )
       dTdz(K)=tmpphi(1,myThid)/tmpphi_area(1,myThid)
c      write(0,*)'level :',k,'  horizontal mean :',dTdz(k)
      enddo
C
C     Vertical gradient of tbarxy
c     DO k=2,Nr
c       dTdz(k) = recip_drC(k)*( tbarxy(k-1)-tbarxy(k) )
c     ENDDO
C     Avoid dividing by zero:
      dTdz(1) = dTdz(2)
C     ==================================================================
C     interpolate dTdz at horizontal boundaries: 
C     ==================================================================
C--   Vertical gradient of dTdz:
      DO k=2,Nr
       DO j=jMin,jMax
        DO i=iMin,iMax       
          d2Tdz2(k) = recip_drC(k)*(dTdz(k)-dTdz(k+1))
        ENDDO
       ENDDO
      ENDDO 
C--
      DO k=2,3
       DO j=jMin,jMax
        DO i=iMin,iMax       
          d2Tdz2(k) = d2Tdz2(4)
        ENDDO
       ENDDO
      ENDDO       
C
C--   Evaluate the gradients of gsfn
       DO j=jMin,jMax
        DO i=iMin,iMax  
          interp2(3) = dTdz(4) + drC(3)*d2Tdz2(3)
        ENDDO
       ENDDO
C
       DO j=jMin,jMax
        DO i=iMin,iMax  
          interp2(2) = interp2(3) + drC(2)*d2Tdz2(2)
        ENDDO
       ENDDO
C
      DO k=4,Nr
       DO j=jMin,jMax
        DO i=iMin,iMax       
          interp2(k) = dTdz(k)
        ENDDO
       ENDDO
      ENDDO  
C
C     ==================================================================
C     Create a threedimensional array of dTdz:
      DO k=1,Nr
         DO j=jMin,jMax
           DO i=iMin,iMax
              dTdz3d(i,j,k) = dTdz(k)   
c            dTdz3d(i,j,k) = interp2(k)
           ENDDO
         ENDDO
      ENDDO 
C     ================================================================== 
C     ================================================================== 
C     ======== Calculate the quasigeostrophic potential vorticity ======
C     ================================================================== 
C
C     q = quasigeostrophic potential vorticity
C
C       -------v--------
C          |       |
C          |       |
C          u   x   u         qgpv is located at p points
C          |  qgpv |
C          |       |
C       -------v--------
C
C     q =    f  + (gsf)_xx + (gsf)_yy + fo^2 ( (gsf)_z / N^2 )_z
C
C          = qf +     qr   +    qs
C
C     ==============================
C             1: Coriolis 
C     ==============================
C     although fCori is defined at u-points they are at the same 
C     latitude as p-points and so do not need to be inerpolated:
      DO j=jMin,jMax
       DO i=iMin,iMax
          qf(i,j) = fCori(i,j,bi,bj)
       ENDDO
      ENDDO
C
C     ==============================
C            2: relative 
C     ==============================
C     !!     for now use gradients of velocity field    !!
C     ==============================
C--   Zonal gradient of vVel
C     ==============================
      DO k=1,Nr
        DO j=jMin,jMax
          DO i=iMin,iMax
            dVdx(i,j,k) = _recip_dxV(i,j,bi,bj) *
     &      ( vVel(i,j,k,bi,bj) - vVel(i-1,j,k,bi,bj) )
          ENDDO
        ENDDO
      ENDDO
C
C--   Interpolate term to center point
      DO k=1,Nr
       DO j=jMin,jMax
        DO i=iMin,iMax
         iG = myXGlobalLo-1+(bi-1)*sNx+I
C
         dVdxbarxy(i,j,k) = 0.25 *
     &    (  dVdx(i,j  ,k) + dVdx(i+1,j  ,k)
     &     + dVdx(i,j+1,k) + dVdx(i+1,j+1,k) )
C
C--      free-slip b.c.'s:
         IF ( iG .EQ. 1 ) THEN
          dVdxbarxy(i,j,k) = 0.5 *
     &    (  dVdx(i+1,j  ,k) + dVdx(i+1,j+1,k) )
         ENDIF
         IF ( iG .EQ. Nx-1 ) THEN
          dVdxbarxy(i,j,k) = 0.5 *
     &    (  dVdx(i,j,k) + dVdx(i,j+1,k) )
         ENDIF
C
        ENDDO
       ENDDO
      ENDDO
C--------------------------------------------------------
C     ==============================
C--   Meridional gradient of uVel
C     ==============================
      DO k=1,Nr
        DO j=jMin,jMax
          DO i=iMin,iMax
            dUdy(i,j,k) = _recip_dyU(i,j,bi,bj) *
     &      ( uVel(i,j,k,bi,bj) - uVel(i,j-1,k,bi,bj) )
          ENDDO
        ENDDO
      ENDDO
C
C--   Interpolate term to center point
      DO k=1,Nr
       DO j=jMin,jMax
        DO i=iMin,iMax
         jG = myYGlobalLo-1+(bj-1)*sNy+J
C
         dUdybarxy(i,j,k) = 0.25 *
     &    (  dUdy(i,j  ,k) + dUdy(i+1,j  ,k)
     &     + dUdy(i,j+1,k) + dUdy(i+1,j+1,k) )
C
C--      free-slip b.c.'s
         IF ( jG .EQ. 1 ) THEN
              dUdybarxy(i,j,k) = 0.5 *
     &        (  dUdy(i, j+1,k) + dUdy(i+1,j+1,k) )
         ENDIF
         IF ( jG .EQ. Ny-1 ) THEN
              dUdybarxy(i,j,k) = 0.5 *
     &        (  dUdy(i, j,k) + dUdy(i+1,j,k) )
         ENDIF
C
        ENDDO
       ENDDO
      ENDDO
C--   =================================
      DO  k=1,Nr
       DO j=jMin,jMax
        DO i=iMin,iMax
          qr(i,j,k) = ( dVdxbarxy(i,j,k) - dUdybarxy(i,j,k) )
        ENDDO
       ENDDO
      ENDDO
C--   =================================
C           3. stretching
C--   =================================
C--   ====a. Vertical gradient of gsfn:
      DO k=2,Nr
       DO j=jMin,jMax
        DO i=iMin,iMax       
          dgsfndz(i,j,k) = recip_drC(k)*( gsfn(i,j,k-1)-gsfn(i,j,k) )*pmask(i,j,k)
        ENDDO
       ENDDO
      ENDDO
C
C--   ====b. Evaluate (gsf)_z / N^2 
      DO k=2,Nr
       DO j=jMin,jMax
        DO i=iMin,iMax  
c         arr(i,j,k) = interp1(i,j,k)  
          arr(i,j,k) = dgsfndz(i,j,k)  
     &                / dTdz3d(i,j,k) 
        ENDDO
       ENDDO
      ENDDO
C     set to zero in the upper and lower layers:
C     to include the pv sheets
      DO j=jMin,jMax
        DO i=iMin,iMax  
          arr(i,j,1)    = 0.
          arr(i,j,Nr+1) = 0.
        ENDDO
      ENDDO
C
C--   ====c. Evaluate z-derivative of (gsf)_z / N^2
      DO k=1,Nr
       DO j=jMin,jMax
        DO i=iMin,iMax       
          arr2(i,j,k) = recip_drF(k)*( arr(i,j,k)-arr(i,j,k+1) )
        ENDDO
       ENDDO
      ENDDO
C
C--   ====d. mulitply by f0^2
      DO k=1,Nr
        DO j=jMin,jMax
          DO i=iMin,iMax    
            qs(i,j,k) = ( f0**2. ) * arr2(i,j,k) * pmask(i,j,k)
          ENDDO
        ENDDO
      ENDDO
C     ==================================================================
C     ==================================================================
C     Sum for q: the quasigeostrophic potential vorticity
C     ==================================================================
      DO k=1,Nr
       DO j=jMin,jMax
        DO i=iMin,iMax  
c         q(i,j,k)= qf(i,j) + qr(i,j,k) + qs(i,j,k)
         q(i,j,k) = qs(i,j,k)
        ENDDO
       ENDDO
      ENDDO
c     write(0,*) (q(10,j,2),j=1,sNy)

C     ==================================================================
C     ==================================================================
C--
      IF ( pvForcing ) THEN 
C--
C--    =================================================================
C--    ================Compute the PV gradients=========================
C--    =================================================================
C--
       IF ( N2local ) THEN
C
C      Interpolate K13 and K23 at horizontal u and v positions but vertical
C      w level
         DO k=1,Nr
          DO j=jMin,jMax
           DO i=iMin,iMax
             arr(i,j,k) = - 0.5 *
     &                      ( K23(i-1,j,k)+K23(i,j,k) ) * umask(i,j,k)
             arrprime(i,j,k) = - 0.5 *
     &                           ( K13(i,j-1,k)+K13(i,j,k) ) * vmask(i,j,k)
           ENDDO
          ENDDO
         ENDDO
         DO j=jMin,jMax
          DO i=iMin,iMax
            arr(i,j,Nr+1) = 0.
            arrprime(i,j,Nr+1) = 0.
          ENDDO
         ENDDO
C      Take z-derivative of arr and arrprime
         DO k=1,Nr
          DO j=jMin,jMax
           DO i=iMin,iMax
             dqdyU(i,j,k) = fCori(i,j,bi,bj) * recip_drF(k) *
     &                             (arr(i,j,k)-arr(i,j,k+1))
             dqdxV(i,j,k) = 0.5 * ( fCori(i,j-1,bi,bj)+fCori(i,j,bi,bj) ) * recip_drF(k) *
     &                           (arrprime(i,j,k)-arrprime(i,j,k+1))
             IF ( PVsheetmld ) THEN
c              ku = int( min(mldindex(i-1,j),mldindex(i,j)) )
c              mldu = min( mld(i-1,j),mld(i,j) )
               ku = int( max(mldindex(i-1,j),mldindex(i,j)) )
               mldu = max( mld(i-1,j),mld(i,j) )
               IF ( k .le. ku .AND. mldu .ne. 0. )
     &           dqdyU(i,j,k) = fCori(i,j,bi,bj)*(-arr(i,j,ku+1))/mldu
c              kv = int( min(mldindex(i,j-1),mldindex(i,j)) )
c              mldv = min( mld(i,j-1),mld(i,j) )
               kv = int( max(mldindex(i,j-1),mldindex(i,j)) )
               mldv = max( mld(i,j-1),mld(i,j) )
               IF ( k .le. kv .AND. mldv .ne. 0. )
     &           dqdxV(i,j,k) = 0.5*(fCori(i,j-1,bi,bj)+fCori(i,j,bi,bj))
     &                                  *(-arrprime(i,j,kv+1))/mldv
             ENDIF
             IF ( betaterm ) THEN
               dqdyU(i,j,k) = dqdyU(i,j,k) + dfdy(i,j,bi,bj)*umask(i,j,k)
             ENDIF
           ENDDO
          ENDDO
         ENDDO
C
        ELSE
C
C      ========================
C--    Meridional gradient of q 
C      ========================
c      write(0,*)'Meridional gradient of q'
c--    At V-points:
       DO k=1,Nr
         DO j=jMin,jMax
          DO i=iMin,iMax
            dqdyV(i,j,k) = _recip_dyU(i,j,bi,bj) 
     &                     * (q(i,j,k)-q(i,j-1,k))
          ENDDO
         ENDDO
       ENDDO
C
C--    ----------------------------
C--    Interpolate term to U points
C--    ----------------------------
C--    ::T-points first
        DO k=1,Nr
         DO j=jMin,jMax
          DO i=iMin,iMax
           jG = myYGlobalLo-1+(bj-1)*sNy+J
           dqdyT(i,j,k) = 0.
           tmp = vmask(i,j,k) + vmask(i,j+1,k)
           IF ( tmp .ne. 0. ) THEN
             dqdyT(i,j,k) = ( dqdyV(i,j,k)*vmask(i,j,k) + dqdyV(i,j+1,k)*vmask(i,j+1,k) )/( vmask(i,j,k) + vmask(i,j+1,k) )
           ENDIF
C
C--        bondaries:
c          IF ( jG .EQ. 1 ) THEN
c             dqdyT(i,j,k) = dqdyV(i, j+1,k) 
c          ENDIF
c          IF ( jG .EQ. sNy-1 ) THEN
c             dqdyT(i,j,k) = dqdyV(i,j,k) 
c          ENDIF
C--         
          ENDDO
         ENDDO
        ENDDO
C--    ::U-points
        DO k=1,Nr
         DO j=jMin,jMax
          DO i=iMin,iMax
           iG = myXGlobalLo-1+(bi-1)*sNx+I
           dqdyU(i,j,k) = 0.
           IF ( umask(i,j,k) .ne. 0. ) THEN
             dqdyU(i,j,k) = ( dqdyT(i,j,k)*pmask(i,j,k) + dqdyT(i-1,j,k)*pmask(i-1,j,k) )/( pmask(i,j,k) + pmask(i-1,j,k) )
           ENDIF
C
C--        bondaries:
c          IF ( iG .EQ. 1 ) THEN
c             dqdyU(i,j,k) = dqdyT(i+1,j,k) 
c          ENDIF
c          IF ( iG .EQ. sNx-1 ) THEN
c             dqdyU(i,j,k) = dqdyT(i,j,k) 
c          ENDIF
C--         
          ENDDO
         ENDDO
        ENDDO

C      ===================
C--    Zonal gradient of q 
C      ===================
c--    At U-points:
c      write(0,*)'Zonal gradient of q'
       DO k=1,Nr
         DO j=jMin,jMax
          DO i=iMin,iMax
            dqdxU(i,j,k) = _recip_dxC(i,j,bi,bj) 
     &                     * ( q(i,j,k)-q(i-1,j,k) )
          ENDDO
         ENDDO
       ENDDO
C
C--    ----------------------------
C--    Interpolate term to V points
C--    ----------------------------
C--    ::T-points first
        DO k=1,Nr
         DO j=jMin,jMax
          DO i=iMin,iMax
           iG = myXGlobalLo-1+(bi-1)*sNx+I
           dqdxT(i,j,k) = 0.
           tmp = umask(i,j,k) + umask(i+1,j,k)
           IF ( tmp .ne. 0. ) THEN
             dqdxT(i,j,k) = ( dqdxU(i,j,k)*umask(i,j,k) + dqdxU(i+1,j,k)*umask(i+1,j,k) )/( umask(i,j,k) + umask(i+1,j,k) )
           ENDIF
C
C--        bondaries:
c          IF ( iG .EQ. 1 ) THEN
c             dqdxT(i,j,k) = dqdxU(i+1,j,k) 
c          ENDIF
c          IF ( iG .EQ. sNx-1 ) THEN
c             dqdxT(i,j,k) = dqdxU(i ,j,k) 
c          ENDIF
C--         
          ENDDO
         ENDDO
        ENDDO
C--    ::V-points
        DO k=1,Nr
         DO j=jMin,jMax
          DO i=iMin,iMax
           jG = myYGlobalLo-1+(bj-1)*sNy+J
           dqdxV(i,j,k) = 0.
           IF ( vmask(i,j,k) .ne. 0. ) THEN
             dqdxV(i,j,k) = ( dqdxT(i,j,k)*pmask(i,j,k) + dqdxT(i,j-1,k)*pmask(i,j-1,k) )/( pmask(i,j,k) + pmask(i,j-1,k) )
           ENDIF
C
C--        bondaries:
c          IF ( jG .EQ. 1 ) THEN
c             dqdxV(i,j,k) = dqdxT(i, j+1,k) 
c          ENDIF
c          IF ( jG .EQ. sNy-1 ) THEN
c             dqdxV(i,j,k) = dqdxT(i,j,k) 
c          ENDIF
C--         
          ENDDO
         ENDDO
        ENDDO
C
       ENDIF
C--
C      ===============
C--    Evaluate Kpv :
C      ===============
C-- 
       IF ( N2local .AND. betaterm ) THEN
c      IF ( N2local ) THEN
c
        DO i=iMin,iMax
         DO j=jMin,jMax
           kbot = int(kbotindex(i,j))
c zonal component
           int1 = 0.
           int2 = 0.
           int3 = 0.
           ku = int( max(mldindex(i-1,j),mldindex(i,j)) )
           DO k=1,ku
             int1 = int1 + delz(k)*dqdyU(i,j,k)
           ENDDO
           IF ( ku .LT. kbot ) int3 = delz(kbot)*dqdyU(i,j,kbot)
           IF ( ku .LT. kbot-1 ) THEN
             DO k=ku+1,kbot-1
               int2 = int2 + delz(k)*dqdyU(i,j,k)
             ENDDO
           ENDIF
           qymld(i,j) = int1
           qyint(i,j) = int2
           qybot(i,j) = int3
c meridional component
           int4 = 0.
           int5 = 0.
           int6 = 0.
           kv = int( max(mldindex(i,j-1),mldindex(i,j)) )
           DO k=1,kv
             int4 = int4 + delz(k)*dqdxV(i,j,k)
           ENDDO
           IF ( kv .LT. kbot ) int6 = delz(kbot)*dqdxV(i,j,kbot)
           IF ( kv .LT. kbot-1 ) THEN
             DO k=kv+1,kbot-1
               int5 = int5 + delz(k)*dqdxV(i,j,k)
             ENDDO
           ENDIF
           qxmld(i,j) = int4
           qxint(i,j) = int5
           qxbot(i,j) = int6
c compute the K
           num_i(i,j) = int3*int4 - int1*int6
           num_b(i,j) = int1*int5 - int2*int4
           denom(i,j) = int2*int6 - int3*int5
           IF ( denom(i,j) .ne. 0. ) THEN
             ratioqy(i,j) = num_i(i,j)/denom(i,j)
             ratio_b(i,j) = num_b(i,j)/denom(i,j)
           ENDIF
           IF ( abs(ratioqy(i,j)-25.) .ge. 25. .OR. 
     &          abs(ratio_b(i,j)-25.) .ge. 25. ) THEN
             fac = 0.
           ELSE
             fac = 1.
           ENDIF
           DO k=1,kbot-1
             IF ( k .le. ku ) THEN
               KpvU(i,j,k) = fac * Kpvref * umask(i,j,k)
             ELSE
               KpvU(i,j,k) = fac * ratioqy(i,j) * Kpvref * umask(i,j,k)
             ENDIF
             IF ( k .le. kv ) THEN 
               KpvV(i,j,k) = fac * Kpvref * vmask(i,j,k)
             ELSE
               KpvV(i,j,k) = fac * ratioqy(i,j) * Kpvref * vmask(i,j,k)
             ENDIF
           ENDDO
           KpvU(i,j,kbot) = fac * ratio_b(i,j) * Kpvref * umask(i,j,kbot)
           KpvV(i,j,kbot) = fac * ratio_b(i,j) * Kpvref * vmask(i,j,kbot)
         ENDDO
        ENDDO
c
       ELSE
c
        DO i=iMin,iMax       
         DO j=jMin,jMax
           pvFacT(i,j) = 0. _d 0
         ENDDO
        ENDDO
C      =========================
C--    Lateral variation of Kpv:
        DO i=iMin,iMax       
         DO j=jMin,jMax
crmw          pvFacT(i,j) = qp2(i,j,bi,bj) * 1.0 _d 0
           pvFacT(i,j) = 1.0 _d 0
         ENDDO
        ENDDO
C--
C      ::At U-point
        DO i=iMin,iMax       
         DO j=jMin,jMax
           iG = myXGlobalLo-1+(bi-1)*sNx+I
C--
           pvFacU(i,j) = 0.5 * ( pvFacT(i,j) + pvFacT(i+1,j) )
C--
C--        bondaries:
           IF ( iG .EQ. 1 ) THEN
              pvFacU(i,j) = pvFacT(i+1,j) 
           ENDIF
C--
         ENDDO
        ENDDO             
C--
C      ::At V-point
        DO i=iMin,iMax       
         DO j=jMin,jMax
           jG = myYGlobalLo-1+(bj-1)*sNy+J
C--
           pvFacV(i,j) = 0.5 * ( pvFacT(i,j) + pvFacT(i,j+1) )
C--
C--        bondaries:
           IF ( jG .EQ. 1 ) THEN
              pvFacV(i,j) = pvFacT(i,j+1) 
           ENDIF
C--
         ENDDO
        ENDDO             
C--
C      =========================
       DO K=1,Nr
        DO i=iMin,iMax
         DO j=jMin,jMax
             KpvU(i,j,K) = Kpvref * pvFacU(i,j)
             KpvV(i,j,K) = Kpvref * pvFacV(i,j)
c          IF ( I .EQ. 10 .AND. J .EQ. 10) THEN         
c          write(0,*) 't0', Kpvref, pvFacU(i,j), KpvU(i,j,K)
c          ENDIF
         ENDDO
        ENDDO
       ENDDO
c
      ENDIF
C--
C      =========================================
C--
C      =======================
C--    Evaluate gUpv and gVpv:
C      =======================
C--
c        write(0,*)' Evaluate gUpv and gVpv'
C
         DO k=1,Nr
          DO j=jMin,jMax
           DO i=iMin,iMax  
             gUpv(i,j,k,bi,bj) = - KpvU(i,j,K) * dqdyU(i,j,k)
c            gUpv(i,j,k,bi,bj) = - Kpvref * ( dqdyU(i,j,k) - dfdy(i,j,bi,bj) )
             gVpv(i,j,k,bi,bj) =   KpvV(i,j,k) * dqdxV(i,j,k)
c            gVpv(i,j,k,bi,bj) =   Kpvref * dqdxV(i,j,k)
           ENDDO
          ENDDO
         ENDDO
C
C--       
       ELSE                   !!!!!!!(pvForcing = .FALSE.)
C--
        IF ( GMadv ) THEN
C
C     Compute horizontal eddy-induced velocities and save them in 
C     arrays gUpv, gVpv
         DO k=1,Nr
          DO j=jMin,jMax
           DO i=iMin,iMax
             arr(i,j,k) = - Kpvref * 0.5 *
     &                      ( K13(i-1,j,k)+K13(i,j,k) ) * umask(i,j,k)
             arrprime(i,j,k) = - Kpvref * 0.5 *
     &                           ( K23(i,j-1,k)+K23(i,j,k) ) * vmask(i,j,k)
           ENDDO
          ENDDO
         ENDDO
         DO j=jMin,jMax
          DO i=iMin,iMax
            arr(i,j,Nr+1) = 0.
            arrprime(i,j,Nr+1) = 0.
          ENDDO
         ENDDO
C      Take z-derivative of arr and arrprime
         DO k=1,Nr
          DO j=jMin,jMax
           DO i=iMin,iMax
             KpvU(i,j,k)          =  0. _d 0
             KpvV(i,j,k)          =  0. _d 0
             gUpv(i,j,k,bi,bj) = recip_drF(k) * (arr(i,j,k)-arr(i,j,k+1))
             gVpv(i,j,k,bi,bj) = recip_drF(k) * (arrprime(i,j,k)-arrprime(i,j,k+1))
           ENDDO
          ENDDO
         ENDDO
C
        ELSE
C
            DO k=1,Nr
              DO j=jMin,jMax
                DO i=iMin,iMax  
                  KpvU(i,j,k)          =  0. _d 0
                  KpvV(i,j,k)          =  0. _d 0
                  gUpv(i,j,k,bi,bj)    =  0. _d 0
                  gVpv(i,j,k,bi,bj)    =  0. _d 0
                ENDDO
              ENDDO
            ENDDO
C
        ENDIF
C--
      ENDIF                  

ctest
       _BARRIER
       _BEGIN_MASTER( myThid )
 
        WRITE(suff,'(I10.10)') myIter           

C--     Write model fields
        CALL WRITE_FLD_XY_RL('kbot.',suff,kbotindex,myIter,myThid)
c       CALL WRITE_FLD_XY_RL('mldindex.',suff,mldindex,myIter,myThid)
c       CALL WRITE_FLD_XY_RL('mld.',suff,mld,myIter,myThid)
c       CALL WRITE_FLD_XY_RL('qymld.',suff,qymld,myIter,myThid)
c       CALL WRITE_FLD_XY_RL('qyint.',suff,qyint,myIter,myThid)
c       CALL WRITE_FLD_XY_RL('qybot.',suff,qybot,myIter,myThid)
c       CALL WRITE_FLD_XY_RL('qxmld.',suff,qxmld,myIter,myThid)
c       CALL WRITE_FLD_XY_RL('qxint.',suff,qxint,myIter,myThid)
c       CALL WRITE_FLD_XY_RL('qxbot.',suff,qxbot,myIter,myThid)
c       CALL WRITE_FLD_XY_RL('ratioqy.',suff,ratioqy,myIter,myThid)
c       CALL WRITE_FLD_XY_RL('ratio_b.',suff,ratio_b,myIter,myThid)
c       CALL WRITE_FLD_XYZ_RL( 'dqdyU.',suff,dqdyU,myIter,myThid)
c       CALL WRITE_FLD_XYZ_RL( 'dqdxV.',suff,dqdxV,myIter,myThid)
c       CALL WRITE_FLD_XYZ_RL( 'KpvU.',suff,KpvU,myIter,myThid)
c       CALL WRITE_FLD_XYZ_RL( 'KpvV.',suff,KpvV,myIter,myThid)
c       CALL WRITE_FLD_XYZ_RL( 'gUpv.',suff,gUpv,myIter,myThid)
c       CALL WRITE_FLD_XYZ_RL( 'gVpv.',suff,gVpv,myIter,myThid)
c       CALL WRITE_FLD_XY_RL('num_i.',suff,num_i,myIter,myThid)
c       CALL WRITE_FLD_XY_RL('num_b.',suff,num_b,myIter,myThid)
c       CALL WRITE_FLD_XY_RL('denom.',suff,denom,myIter,myThid)

       _END_MASTER( myThid )
       _BARRIER 
ctest
C--
C     =========================================================
c      DO k=1,Nr
c       DO j=jMin,jMax
c        DO i=iMin,iMax  
c         write(0,*) 'TTTT ',  Kpvref, pvFacU(i,j)
c        ENDDO
c       ENDDO
c      ENDDO
C--
C     =========================================================

      RETURN
      END