C $Header: /u/gcmpack/MITgcm/pkg/timeave/timeave_cumul_2v.F,v 1.7 2005/08/19 22:50:26 heimbach Exp $
C $Name:  $
#include "TIMEAVE_OPTIONS.h"

CStartofinterface
      SUBROUTINE TIMEAVE_CUMUL_2V( 
     O   fldtave,
     I   fld1, fld2, Ksize, dir, deltaTloc,
     I   bi, bj, myThid )
C     /==========================================================*
C     | SUBROUTINE TIMEAVE_CUMUL_2V
C     | o Sum over time a product of two arrays depending on the
C     |   relative position of the 2 fields.
C     |      (tracer point, u, v, w ...)
C     \==========================================================*
      IMPLICIT NONE

C     == Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "GRID.h"

C     == Routine arguments ==
C     myThid - Thread number for this instance of the routine.
C     fldtave - time averaged Field
C     fld1,fld2  - Input Field
C     dir - type of grid for 2nd array relatively to the 1rst array
C     0: same grid ; 1: dX/2 shift ; 2: dY/2 shift ; 3: dr/2 shift
C        (2 digits => also shift the 1rst array)
C     Ksize - 3rd dimension of local arrays (Input and Output fields)
      INTEGER Ksize, dir
      _RL fld1(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Ksize,nSx,nSy)
      _RL fld2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Ksize,nSx,nSy)
      _RL fldtave(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Ksize,nSx,nSy)
      _RL deltaTloc
      INTEGER bi, bj, myThid

CEndofinterface

#ifdef ALLOW_TIMEAVE
C     == Local variables ==
C     i,j,k,bi,bj  - Loop counters
      INTEGER i, j, k
      INTEGER km1
 
      IF ( dir.eq.0 ) THEN 
c-    both fields at the same location :

C     DO bj = myByLo(myThid), myByHi(myThid)
C      DO bi = myBxLo(myThid), myBxHi(myThid)
        DO k=1,Ksize
         DO j=1,sNy
          DO i=1,sNx
             fldtave(i,j,k,bi,bj)= fldtave(i,j,k,bi,bj)
     &       +  fld1(i,j,k,bi,bj)*fld2(i,j,k,bi,bj)*deltaTloc
          ENDDO
         ENDDO
        ENDDO
C      ENDDO
C     ENDDO

      ELSEIF ( dir.eq.1 )  THEN
c-    2nd field shifted by -dX/2 (e.g.: 1=T, 2=U) :
 
        DO k=1,Ksize
         DO j=1,sNy
          DO i=1,sNx
             fldtave(i,j,k,bi,bj)= fldtave(i,j,k,bi,bj)
     &       + .5 * ( fld1(i-1,j,k,bi,bj) + fld1(i,j,k,bi,bj) )
     &                *   fld2(i,j,k,bi,bj)
     &                *   deltaTloc
          ENDDO
         ENDDO
        ENDDO

      ELSEIF ( dir.eq.2 ) THEN
c-    2nd field shifted by -dY/2 (e.g.: 1=T, 2=V) :
 
        DO k=1,Ksize
         DO j=1,sNy
          DO i=1,sNx
             fldtave(i,j,k,bi,bj)= fldtave(i,j,k,bi,bj)
     &       + .5 * ( fld1(i,j-1,k,bi,bj) + fld1(i,j,k,bi,bj) )
     &                *   fld2(i,j,k,bi,bj)
     &                *   deltaTloc
          ENDDO
         ENDDO
        ENDDO

      ELSEIF ( dir.eq.3 ) THEN
c-    2nd field shifted by -dR/2 (e.g.: 1=T, 2=W) :
 
        DO k=1,Ksize
         km1 = MAX(k-1,1)
         DO j=1,sNy
          DO i=1,sNx
             fldtave(i,j,k,bi,bj)= fldtave(i,j,k,bi,bj)
     &       + .5 * ( fld1(i,j,km1,bi,bj) + fld1(i,j,k,bi,bj) )
     &                *   fld2(i,j,k,bi,bj)
     &                *   deltaTloc
          ENDDO
         ENDDO
        ENDDO

      ELSEIF ( dir.eq.12 ) THEN
c-    1rst & 2nd fields shifted by -dY/2 & -dX/2 
c           (e.g.: 1=U, 2=V, product at the corner) :
 
        DO k=1,Ksize
         DO j=1,sNy
          DO i=1,sNx
           fldtave(i,j,k,bi,bj) = fldtave(i,j,k,bi,bj)
     &      + .25 _d 0*( fld1(i,j-1,k,bi,bj) + fld1(i,j,k,bi,bj) )
     &                *( fld2(i-1,j,k,bi,bj) + fld2(i,j,k,bi,bj) )
     &                * deltaTloc
          ENDDO
         ENDDO
        ENDDO

      ELSEIF ( dir.eq.13 ) THEN
c-    1rst & 2nd fields shifted by -dR/2 & -dX/2 (e.g.: 1=U, 2=W):
 
        DO k=1,Ksize
         km1 = MAX(k-1,1)
         DO j=1,sNy
          DO i=1,sNx
           fldtave(i,j,k,bi,bj) = fldtave(i,j,k,bi,bj)
     &      + .25 _d 0*( fld1(i,j,km1,bi,bj) + fld1(i,j,k,bi,bj) )
     &                *( fld2(i-1,j,k,bi,bj)*rA(i-1,j,bi,bj)
     &                  +fld2( i ,j,k,bi,bj)*rA( i ,j,bi,bj)
     &                 )*recip_rAw(i,j,bi,bj)
     &                * deltaTloc
          ENDDO
         ENDDO
        ENDDO

      ELSEIF ( dir.eq.23 ) THEN
c-    1rst & 2nd fields shifted by -dR/2 & -dY/2 (e.g.: 1=V, 2=W):
 
        DO k=1,Ksize
         km1 = MAX(k-1,1)
         DO j=1,sNy
          DO i=1,sNx
           fldtave(i,j,k,bi,bj) = fldtave(i,j,k,bi,bj)
     &      + .25 _d 0*( fld1(i,j,km1,bi,bj) + fld1(i,j,k,bi,bj) )
     &                *( fld2(i,j-1,k,bi,bj)*rA(i,j-1,bi,bj)
     &                  +fld2(i, j ,k,bi,bj)*rA(i, j ,bi,bj)
     &                 )*recip_rAs(i,j,bi,bj)
     &                * deltaTloc
          ENDDO
         ENDDO
        ENDDO

      ELSEIF ( dir.eq.-13 ) THEN
c-    gradient of the 1rst field * 2nd fields, shifted by -dR/2 & -dX/2 resp.
c-    (e.g.: used for advective form of vertical advection: w.du/dr)
 
        DO k=2,Ksize
         DO j=1,sNy
          DO i=1,sNx
           fldtave(i,j,k,bi,bj) = fldtave(i,j,k,bi,bj)
     &       + .5 _d 0*( fld1(i,j,k-1,bi,bj) - fld1(i,j,k,bi,bj) )
     &                *( fld2(i-1,j,k,bi,bj)*rA(i-1,j,bi,bj)
     &                  +fld2( i ,j,k,bi,bj)*rA( i ,j,bi,bj)
     &                 )*recip_rAw(i,j,bi,bj)
     &                * deltaTloc
          ENDDO
         ENDDO
        ENDDO

      ELSEIF ( dir.eq.-23 ) THEN
c-    gradient of the 1rst field * 2nd fields, shifted by -dR/2 & -dY/2 resp.
c-    (e.g.: used for advective form of vertical advection: w.dv/dr)
 
        DO k=2,Ksize
         DO j=1,sNy
          DO i=1,sNx
           fldtave(i,j,k,bi,bj) = fldtave(i,j,k,bi,bj)
     &       + .5 _d 0*( fld1(i,j,k-1,bi,bj) - fld1(i,j,k,bi,bj) )
     &                *( fld2(i,j-1,k,bi,bj)*rA(i,j-1,bi,bj)
     &                  +fld2(i, j ,k,bi,bj)*rA(i, j ,bi,bj)
     &                 )*recip_rAs(i,j,bi,bj)
     &                * deltaTloc
          ENDDO
         ENDDO
        ENDDO

      ENDIF

#endif /* ALLOW_TIMEAVE */

      RETURN
      END