C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/CS_ADJOINT_TESTS/diffuse/code_ad/diffuse_theta.F,v 1.1 2004/07/25 22:51:38 heimbach Exp $
C $Name:  $

#include "PACKAGES_CONFIG.h"
#include "CPP_OPTIONS.h"

CBOP
C     !ROUTINE: DIFFUSE_THETA
C     !INTERFACE:
      SUBROUTINE DIFFUSE_THETA(myTime, myIter, myThid)

C     !USES:
      IMPLICIT NONE
C     == Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#include "GRID.h"
#include "GAD.h"

C     !INPUT/OUTPUT PARAMETERS:
C     == Routine arguments ==
C     myTime - Current time in simulation
C     myIter - Current iteration number in simulation
C     myThid - Thread number for this instance of the routine.
      _RL myTime
      INTEGER myIter
      INTEGER myThid

C     !LOCAL VARIABLES:
C     == Local variables
C     xA, yA                 - Per block temporaries holding face areas
C     uTrans, vTrans         - Per block temporaries holding flow
C                              transport
C                              o uTrans: Zonal transport
C                              o vTrans: Meridional transport
C     maskUp                   o maskUp: land/water mask for W points
C     fVer                     o fVer: Vertical flux term - note fVer
C                                      is "pipelined" in the vertical
C                                      so we need an fVer for each
C                                      variable.
C     iMin, iMax     - Ranges and sub-block indices on which calculations
C     jMin, jMax       are applied.
C     bi, bj
C     k, kup,        - Index for layer above and below. kup and kDown
C     kDown, km1       are switched with layer to be the appropriate
C                      index into fVerTerm.
      _RS xA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RS yA      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL uTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL vTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL rTrans  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RS maskUp  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL fVer    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,2)
      INTEGER iMin, iMax
      INTEGER jMin, jMax
      INTEGER bi, bj
      INTEGER i, j
      INTEGER k, km1, kup, kDown

      _RL df    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL localT(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL fZon  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL fMer  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
CEOP

      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          xA(i,j)        = 0. _d 0
          yA(i,j)        = 0. _d 0
          uTrans(i,j)    = 0. _d 0
          vTrans(i,j)    = 0. _d 0
          rTrans (i,j)   = 0. _d 0
          fVer   (i,j,1) = 0. _d 0
          fVer   (i,j,2) = 0. _d 0
         ENDDO
        ENDDO

        DO k=1,Nr
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           gT(i,j,k,bi,bj)   = 0. _d 0
          ENDDO
         ENDDO
        ENDDO

        DO k=Nr,1,-1

         km1  = MAX(1,k-1)
         kup  = 1+MOD(k+1,2)
         kDown= 1+MOD(k,2)

         iMin = 1-OLx
         iMax = sNx+OLx
         jMin = 1-OLy
         jMax = sNy+OLy
C--      Get temporary terms used by tendency routines
         CALL CALC_COMMON_FACTORS (
     I        bi,bj,iMin,iMax,jMin,jMax,k,
     O        xA,yA,uTrans,vTrans,rTrans,maskUp,
     I        myThid)

C--      Make local copy of tracer array
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           localT(i,j)    =theta(i,j,k,bi,bj)
          ENDDO
         ENDDO
C--      Zero out work arrays
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           fZon(i,j)      = 0. _d 0
           fMer(i,j)      = 0. _d 0
           df(i,j)        = 0. _d 0
          ENDDO
         ENDDO

C--      Diffusive fluxes
         CALL GAD_DIFF_X(bi,bj,k,xA,diffKhT,localT,df,myThid)
         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx
           fZon(i,j) = fZon(i,j) + df(i,j)
          ENDDO
         ENDDO
         CALL GAD_DIFF_Y(bi,bj,k,yA,diffKhT,localT,df,myThid)
         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx
           fMer(i,j) = fMer(i,j) + df(i,j)
          ENDDO
         ENDDO

C--      Divergence of fluxes
         DO j=1-Oly,sNy+Oly-1
          DO i=1-Olx,sNx+Olx-1
           gT(i,j,k,bi,bj)=gT(i,j,k,bi,bj)
     &      -_recip_hFacC(i,j,k,bi,bj)*recip_drF(k)*recip_rA(i,j,bi,bj)
     &      *( (fZon(i+1,j)-fZon(i,j))
     &        +(fMer(i,j+1)-fMer(i,j))
     &        +(fVer(i,j,kUp)-fVer(i,j,kDown))*rkFac
     &       )

          ENDDO
         ENDDO

C--      Step field forward
         CALL ADAMS_BASHFORTH2(
     I                        bi, bj, K,
     U                        gT, gTnm1,
     I                        myIter, myThid )
         CALL TIMESTEP_TRACER(
     I        bi,bj,iMin,iMax,jMin,jMax,k,tempAdvScheme,
     I        theta, gT,
     I        myIter, myThid)

        ENDDO
       ENDDO
      ENDDO

C--   Update halo regions
      _EXCH_XYZ_R8( gTNM1 , myThid )
      _EXCH_XYZ_R8( gT    , myThid )

C--   Cycle timestepping arrays
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
        DO k=1,Nr
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           theta(i,j,k,bi,bj)=gT(i,j,k,bi,bj)
          ENDDO
         ENDDO
        ENDDO
       ENDDO
      ENDDO

      CALL PLOT_FIELD_XYZRL( theta, 'Theta after diffuse' ,
     &                        Nr, myIter, myThid )

      RETURN
      END