C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/model/src/impldiff.F,v 1.11 2000/09/11 20:50:39 heimbach Exp $

#include "CPP_OPTIONS.h"

C     /==========================================================\
C     | S/R IMPLDIFF                                             |
C     | o Solve implicit diffusion equation for vertical         |
C     |   diffusivity.                                           |
C     \==========================================================/
      SUBROUTINE IMPLDIFF( bi, bj, iMin, iMax, jMin, jMax,
     I                     deltaTX,KappaRX,recip_hFac,
     U                     gXNm1,
     I                     myThid )
      IMPLICIT NONE
C     == Global data ==
#include "SIZE.h"
#include "DYNVARS.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"

#ifdef ALLOW_AUTODIFF_TAMC
#include "tamc_keys.h"
#endif

C     == Routine Arguments ==
      INTEGER bi,bj,iMin,iMax,jMin,jMax
      _RL deltaTX
      _RL KappaRX(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)
      _RS recip_hFac(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
      _RL gXnm1(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr,nSx,nSy)
      INTEGER myThid

C     == Local variables ==
      INTEGER i,j,k
      _RL a(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL b(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL c(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL ckm1(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL bet(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
      _RL gam(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nr)

#ifdef ALLOW_AUTODIFF_TAMC
      INTEGER kkey
#endif

C--   Only need do anything if Nr>1
      IF (Nr.GT.1) THEN

C--    Beginning of forward sweep (top level)
       DO j=jMin,jMax
        DO i=iMin,iMax
         c(i,j)=-deltaTX*recip_hFac(i,j,1,bi,bj)*recip_drF(1)
     &         *KappaRX(i,j,2)*recip_drC(2)
         b(i,j)=1.-c(i,j)
         bet(i,j)=0.
         IF (b(i,j).NE.0.) bet(i,j)=1. / b(i,j)
        ENDDO
       ENDDO

      ENDIF

#ifdef ALLOW_AUTODIFF_TAMC
CADJ store bet                 = comlev1_impl, key = idkey
CADJ store gXNm1(:,:,:,bi,bj)  = comlev1_impl, key = idkey
#endif

       DO j=jMin,jMax
        DO i=iMin,iMax
         gXNm1(i,j,1,bi,bj) = gXNm1(i,j,1,bi,bj)*bet(i,j)
        ENDDO
       ENDDO

C--   Middle of forward sweep
      IF (Nr.GT.2) THEN

       DO k=2,Nr-1

#ifdef ALLOW_AUTODIFF_TAMC
        kkey = (idkey-1)*(Nr-2) + k-1
#endif

        DO j=jMin,jMax
         DO i=iMin,iMax
          ckm1(i,j)=c(i,j)
          a(i,j)=-deltaTX*recip_hFac(i,j,k,bi,bj)*recip_drF(k)
     &         *KappaRX(i,j, k )*recip_drC( k )
          c(i,j)=-deltaTX*recip_hFac(i,j,k,bi,bj)*recip_drF(k)
     &         *KappaRX(i,j,k+1)*recip_drC(k+1)
          b(i,j)=1.-c(i,j)-a(i,j)
         ENDDO
        ENDDO

#ifdef ALLOW_AUTODIFF_TAMC
CADJ store ckm1, bet = comlev1_impl_k, key = kkey
#endif

        DO j=jMin,jMax
         DO i=iMin,iMax
          gam(i,j,k)=ckm1(i,j)*bet(i,j)
         ENDDO
        ENDDO

        DO j=jMin,jMax
         DO i=iMin,iMax
          bet(i,j)=b(i,j)-a(i,j)*gam(i,j,k)
          IF (bet(i,j).NE.0.) bet(i,j)=1. / bet(i,j)
         ENDDO
        ENDDO

#ifdef ALLOW_AUTODIFF_TAMC
CADJ store bet                    = comlev1_impl_k, key = kkey
CADJ store gXNm1(:,:,k-1:k,bi,bj) = comlev1_impl_k, key = kkey
#endif

        DO j=jMin,jMax
         DO i=iMin,iMax
          gXnm1(i,j,k,bi,bj)=(gXnm1(i,j,k,bi,bj)
     &       -a(i,j)*gXnm1(i,j,k-1,bi,bj))*bet(i,j)
         ENDDO
        ENDDO

       ENDDO

      ENDIF


      IF (Nr.GT.1) THEN

C--    End of forward sweep (bottom level)
       DO j=jMin,jMax
        DO i=iMin,iMax
         ckm1(i,j)=c(i,j)
         a(i,j)=-deltaTX*recip_hFac(i,j,Nr,bi,bj)*recip_drF(Nr)
     &        *KappaRX(i,j, Nr )*recip_drC( Nr )
         b(i,j)=1.-a(i,j)
        ENDDO
       ENDDO

#ifdef ALLOW_AUTODIFF_TAMC
CADJ store ckm1 = comlev1_impl, key = idkey
CADJ store a,b  = comlev1_impl, key = idkey
CADJ store bet  = comlev1_impl, key = idkey
#endif

       DO j=jMin,jMax
        DO i=iMin,iMax
         gam(i,j,Nr)=ckm1(i,j)*bet(i,j)
        ENDDO
       ENDDO
       DO j=jMin,jMax
        DO i=iMin,iMax
         bet(i,j)=b(i,j)-a(i,j)*gam(i,j,Nr)
         IF (bet(i,j).NE.0.) bet(i,j)=1. / bet(i,j)
        ENDDO
       ENDDO

#ifdef ALLOW_AUTODIFF_TAMC
CADJ store a,bet              = comlev1_impl, key = idkey
CADJ store gXnm1(:,:,:,bi,bj) = comlev1_impl, key = idkey
#endif

       DO j=jMin,jMax
        DO i=iMin,iMax
         gXnm1(i,j,Nr,bi,bj)=(gXnm1(i,j,Nr,bi,bj)
     &      -a(i,j)*gXnm1(i,j,Nr-1,bi,bj))*bet(i,j)
        ENDDO
       ENDDO

#ifdef ALLOW_AUTODIFF_TAMC
CADJ store gam   = comlev1_impl, key = idkey
#endif

C--    Backward sweep
       DO k=Nr-1,1,-1
        DO j=jMin,jMax
         DO i=iMin,iMax
          gXnm1(i,j,k,bi,bj)=gXnm1(i,j,k,bi,bj)
     &              -gam(i,j,k+1)*gXnm1(i,j,k+1,bi,bj)
         ENDDO
        ENDDO
       ENDDO

      ENDIF

      RETURN
      END