C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/sciascia/rbcs/rbcs_add_tendency.F,v 1.1 2012/08/08 01:57:14 heimbach Exp $
C $Name:  $

#include "RBCS_OPTIONS.h"

CBOP
C !ROUTINE: RBCS_ADD_TSTENDENCY

C !INTERFACE: ==========================================================
      SUBROUTINE RBCS_ADD_TENDENCY(
     I                    bi, bj, k, tracerNum,
     I                    myTime, myThid )

C !DESCRIPTION:
C     Will update tendencies with terms to relax to 3-D field

C !USES: ===============================================================
      IMPLICIT NONE
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
c#include "GRID.h"
#include "DYNVARS.h"
#include "NH_VARS.h"
#ifdef ALLOW_PTRACERS
#include "PTRACERS_SIZE.h"
#include "PTRACERS_FIELDS.h"
#endif
#include "RBCS_SIZE.h"
#include "RBCS_PARAMS.h"
#include "RBCS_FIELDS.h"

C !INPUT PARAMETERS: ===================================================
C  bi,bj          :: tile indices
C  k              :: vertical level index
C  tracerNum      :: tracer number (1=Temp, 2=Salt, >2 : ptracer)
C  myTime         :: current time
C  myIter         :: current timestep
C  myThid         :: my Thread Id number
      INTEGER bi, bj, k
      INTEGER tracerNum
      _RL myTime
c     INTEGER myIter
      INTEGER myThid

C !LOCAL VARIABLES: ====================================================
C  i,j            :: loop indices
      INTEGER i,j
      _RL rbcsVanishingFac
      INTEGER ium,ivm,iwm 
#ifdef ALLOW_PTRACERS
      INTEGER iTracer
      INTEGER irbc
#endif
CEOP

#ifdef ALLOW_RBCS

      if (rbcsVanishingTime.GT.0. _d 0) then
        rbcsVanishingFac = 
     &  MAX( 0. _d 0 , 1. _d 0 - mytime / rbcsVanishingTime )
      else
        rbcsVanishingFac = 1. _d 0
      endif
 
#ifndef DISABLE_RBCS_MOM
      DO ium=1,UmLEN
       IF ( tracerNum.EQ.-1 .AND. useRBCuVel(ium) ) THEN
        DO j=0,sNy+1
         DO i=0,sNx+1
          gU(i,j,k,bi,bj) = gU(i,j,k,bi,bj)
     &       - RBC_maskU(i,j,k,bi,bj,ium)/tauRelaxU(ium)
     &        *( uVel(i,j,k,bi,bj)- RBCuVel(i,j,k,bi,bj,ium) )
c     &        *maskW(i,j,k,bi,bj)
     &        *rbcsVanishingFac         
         ENDDO
        ENDDO
       ENDIF
      ENDDO
      DO ivm=1,VmLEN
       IF ( tracerNum.EQ.-2 .AND. useRBCvVel(ivm) ) THEN
        DO j=0,sNy+1
         DO i=0,sNx+1
          gV(i,j,k,bi,bj) = gV(i,j,k,bi,bj)
     &       - RBC_maskV(i,j,k,bi,bj,ivm)/tauRelaxV(ivm)
     &        *( vVel(i,j,k,bi,bj)- RBCvVel(i,j,k,bi,bj,ivm) )
c     &        *maskS(i,j,k,bi,bj)
     &        *rbcsVanishingFac         
         ENDDO
        ENDDO
       ENDIF
      ENDDO
#ifdef ALLOW_NONHYDROSTATIC
      DO iwm=1,WmLEN
       IF ( tracerNum.EQ.-3 .AND. useRBCwVel(iwm) ) THEN
        DO j=0,sNy+1
         DO i=0,sNx+1
          gW(i,j,k,bi,bj) = gW(i,j,k,bi,bj)
     &       - RBC_maskW(i,j,k,bi,bj,iwm)/tauRelaxW(iwm)
     &        *( wVel(i,j,k,bi,bj)- RBCwVel(i,j,k,bi,bj,iwm) )
c     &        *maskS(i,j,k,bi,bj)
     &        *rbcsVanishingFac         
         ENDDO
        ENDDO
       ENDIF
      ENDDO
#endif
#endif /* DISABLE_RBCS_MOM */

      IF ( tracerNum.EQ.1 .AND. useRBCtemp ) THEN
        DO j=1,sNy
         DO i=1,sNx
          gT(i,j,k,bi,bj) = gT(i,j,k,bi,bj)
     &       - RBC_mask(i,j,k,bi,bj,1)/tauRelaxT
     &        *( theta(i,j,k,bi,bj)- RBCtemp(i,j,k,bi,bj) )
c    &        *maskC(i,j,k,bi,bj)
     &        *rbcsVanishingFac

         ENDDO
        ENDDO
      ENDIF

      IF ( tracerNum.EQ.2 .AND. useRBCsalt ) THEN
        DO j=1,sNy
         DO i=1,sNx
          gS(i,j,k,bi,bj) = gS(i,j,k,bi,bj)
     &       - RBC_mask(i,j,k,bi,bj,2)/tauRelaxS
     &        *( salt(i,j,k,bi,bj)- RBCsalt(i,j,k,bi,bj) )
c    &        *maskC(i,j,k,bi,bj)
     &        *rbcsVanishingFac
         ENDDO
        ENDDO
      ENDIF

#ifdef ALLOW_PTRACERS
      IF ( usePTRACERS .AND. tracerNum.GT.2 ) THEN
       iTracer = tracerNum-2
       irbc = MIN(maskLEN,tracerNum)
       IF ( useRBCptrnum(iTracer) ) THEN
        DO j=1,sNy
         DO i=1,sNx
          gPtr(i,j,k,bi,bj,iTracer) = gPtr(i,j,k,bi,bj,iTracer)
     &       - RBC_mask(i,j,k,bi,bj,irbc)/tauRelaxPTR(iTracer)
     &        *( pTracer(i,j,k,bi,bj,iTracer)
     &           - RBC_ptracers(i,j,k,bi,bj,iTracer) )
c    &        *maskC(i,j,k,bi,bj)
     &        *rbcsVanishingFac
         ENDDO
        ENDDO
       ENDIF
      ENDIF
#endif /* ALLOW_PTRACERS */

#endif /* ALLOW_RBCS */

      RETURN
      END