C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_advdiff.F,v 1.74 2014/10/20 03:20:57 gforget Exp $
C $Name:  $

#include "SEAICE_OPTIONS.h"
#ifdef ALLOW_AUTODIFF
# include "AUTODIFF_OPTIONS.h"
#endif

CBOP
C !ROUTINE: SEAICE_ADVDIFF

C !INTERFACE: ==========================================================
      SUBROUTINE SEAICE_ADVDIFF(
     I     uc, vc, myTime, myIter, myThid )

C !DESCRIPTION: \bv
C     *===========================================================*
C     | SUBROUTINE SEAICE_ADVDIFF
C     | o driver for different advection routines
C     |   calls an adaption of gad_advection to call different
C     |   advection routines of pkg/generic_advdiff
C     *===========================================================*
C \ev

C !USES: ===============================================================
      IMPLICIT NONE

C     === Global variables ===
C     UICE/VICE :: ice velocity
C     HEFF      :: scalar field to be advected
C     HEFFM     :: mask for scalar field
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#include "GAD.h"
#include "SEAICE_SIZE.h"
#include "SEAICE_PARAMS.h"
#include "SEAICE.h"
#include "SEAICE_TRACER.h"

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

C !INPUT PARAMETERS: ===================================================
C     === Routine arguments ===
C     myTime    :: current time
C     myIter    :: iteration number
C     myThid    :: Thread no. that called this routine.
      _RL myTime
      INTEGER myIter
      INTEGER myThid
CEndOfInterface

C !LOCAL VARIABLES: ====================================================
C     === Local variables ===
C     i,j,bi,bj :: Loop counters
#ifdef SEAICE_ITD
C     it        :: Loop counter for ice thickness categories
#endif /* SEAICE_ITD */
C     ks        :: surface level index
C     uc/vc     :: current ice velocity on C-grid
C     uTrans    :: volume transport, x direction
C     vTrans    :: volume transport, y direction
C     afx       :: horizontal advective flux, x direction
C     afy       :: horizontal advective flux, y direction
C     gFld      :: tendency of seaice field
C     xA,yA     :: "areas" of X and Y face of tracer cells
      INTEGER i, j, bi, bj
#ifdef SEAICE_ITD
      INTEGER it
#endif /* SEAICE_ITD */
      INTEGER ks
      LOGICAL SEAICEmultiDimAdvection
#ifdef ALLOW_AUTODIFF_TAMC
      INTEGER itmpkey
#endif /* ALLOW_AUTODIFF_TAMC */
#ifdef ALLOW_SITRACER
      _RL hEffNm1    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL areaNm1    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      INTEGER iTr, SEAICEadvSchSItr
      _RL SEAICEdiffKhSItr
      _RL SItrExt   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL tmpscal1, tmpscal2
# ifdef ALLOW_SITRACER_ADVCAP
      _RL SItrPrev  (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
# endif
# ifdef ALLOW_SITRACER_DEBUG_DIAG
      _RL DIAGarray (1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)
# endif
#endif /* ALLOW_SITRACER */

      _RL uc        (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL vc        (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL fldNm1    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL uTrans    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL vTrans    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL afx       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL afy       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL gFld      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RS xA        (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RS yA        (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
      _RL recip_heff(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
CEOP

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

      ks = 1

C--   make a local copy of the velocities for compatibility with B-grid
C--   alternatively interpolate to C-points if necessary
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
#ifndef SEAICE_CGRID /* not SEAICE_CGRID = BGRID */
C--   hack to ensure backward compatibility: 
C     average B-grid seaice velocity to C-grid
        DO j=1-OLy,sNy+OLy-1
         DO i=1-OLx,sNx+OLx-1
          uc(i,j,bi,bj)=.5 _d 0*(UICE(i,j,bi,bj)+UICE(i,j+1,bi,bj))
          vc(i,j,bi,bj)=.5 _d 0*(VICE(i,j,bi,bj)+VICE(i+1,j,bi,bj))
         ENDDO
        ENDDO
#endif /* SEAICE_CGRID */
C-    compute cell areas used by all tracers
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          xA(i,j,bi,bj) = _dyG(i,j,bi,bj)*_maskW(i,j,ks,bi,bj)
          yA(i,j,bi,bj) = _dxG(i,j,bi,bj)*_maskS(i,j,ks,bi,bj)
         ENDDO
        ENDDO
       ENDDO
      ENDDO

#ifndef SEAICE_CGRID
C     Do we need this? I am afraid so.
      CALL EXCH_UV_XY_RL(uc,vc,.TRUE.,myThid)
#endif /* not SEAICE_CGRID */

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE uc   = comlev1, key = ikey_dynamics, kind=isbyte
CADJ STORE vc   = comlev1, key = ikey_dynamics, kind=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */

      SEAICEmultidimadvection = .TRUE.
      IF ( SEAICEadvScheme.EQ.ENUM_CENTERED_2ND
     & .OR.SEAICEadvScheme.EQ.ENUM_UPWIND_3RD
     & .OR.SEAICEadvScheme.EQ.ENUM_CENTERED_4TH ) THEN
       SEAICEmultiDimAdvection = .FALSE.
      ENDIF

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE heffm  = comlev1, key = ikey_dynamics, kind=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
      IF ( SEAICEmultiDimAdvection ) THEN

      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)
C---   loops on tile indices bi,bj

#ifdef ALLOW_AUTODIFF_TAMC
C     Initialise for TAF
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          gFld(i,j)       = 0. _d 0
         ENDDO
        ENDDO
C
        act1 = bi - myBxLo(myThid)
        max1 = myBxHi(myThid) - myBxLo(myThid) + 1
        act2 = bj - myByLo(myThid)
        max2 = myByHi(myThid) - myByLo(myThid) + 1
        act3 = myThid - 1
        max3 = nTx*nTy
        act4 = ikey_dynamics - 1
        itmpkey = (act1 + 1) + act2*max1
     &                       + act3*max1*max2
     &                       + act4*max1*max2*max3
C
CADJ STORE area(:,:,bi,bj)  = comlev1_bibj,
CADJ &     key = itmpkey, kind=isbyte
CADJ STORE heff(:,:,bi,bj)  = comlev1_bibj,
CADJ &     key = itmpkey, kind=isbyte
CADJ STORE heffm(:,:,bi,bj)  = comlev1_bibj,
CADJ &     key = itmpkey, kind=isbyte
CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
CADJ &     key = itmpkey, kind=isbyte
# ifdef SEAICE_VARIABLE_SALINITY
CADJ STORE hsalt(:,:,bi,bj)  = comlev1_bibj,
CADJ &     key = itmpkey, kind=isbyte
# endif
#endif /* ALLOW_AUTODIFF_TAMC */

        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
#ifdef ALLOW_SITRACER
          hEffNm1(i,j,bi,bj) = HEFF(i,j,bi,bj)
          areaNm1(i,j,bi,bj) = AREA(i,j,bi,bj)
#endif
          recip_heff(i,j)    = 1. _d 0
         ENDDO
        ENDDO

C-    Calculate "volume transports" through tracer cell faces.
        DO j=1-OLy,sNy+OLy
         DO i=1-OLx,sNx+OLx
          uTrans(i,j) = uc(i,j,bi,bj)*xA(i,j,bi,bj)
          vTrans(i,j) = vc(i,j,bi,bj)*yA(i,j,bi,bj)
         ENDDO
        ENDDO

C--   Effective Thickness (Volume)
        IF ( SEAICEadvHeff ) THEN
#ifdef SEAICE_ITD
        DO it=1,SEAICE_multDim
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           HEFF(i,j,bi,bj)=HEFFITD(i,j,it,bi,bj)
          ENDDO
         ENDDO
#endif /* SEAICE_ITD */
          CALL SEAICE_ADVECTION(
     I         GAD_HEFF, SEAICEadvSchHeff,
     I         uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj),
     I         uTrans, vTrans, HEFF(1-OLx,1-OLy,bi,bj), recip_heff,
     O         gFld, afx, afy,
     I         bi, bj, myTime, myIter, myThid )
         IF ( SEAICEdiffKhHeff .GT. 0. _d 0 ) THEN
C-    Add tendency due to diffusion
          CALL SEAICE_DIFFUSION(
     I         GAD_HEFF, SEAICEdiffKhHeff, ONE,
     I         HEFF(1-OLx,1-OLy,bi,bj), HEFFM,
     I         xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
     U         gFld,
     I         bi, bj, myTime, myIter, myThid )
         ENDIF
C     now do the "explicit" time step
         DO j=1,sNy
          DO i=1,sNx
           HEFF(i,j,bi,bj) = HEFFM(i,j,bi,bj) * (
     &          HEFF(i,j,bi,bj) + SEAICE_deltaTtherm * gFld(i,j)
     &          )
          ENDDO
         ENDDO
#ifdef SEAICE_ITD
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           HEFFITD(i,j,it,bi,bj)=HEFF(i,j,bi,bj)
          ENDDO
         ENDDO
        ENDDO
#endif /* SEAICE_ITD */
        ENDIF

C--   Fractional area
        IF ( SEAICEadvArea ) THEN
#ifdef SEAICE_ITD
        DO it=1,SEAICE_multDim
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           AREA(i,j,bi,bj)=AREAITD(i,j,it,bi,bj)
          ENDDO
         ENDDO
#endif
          CALL SEAICE_ADVECTION(
     I         GAD_AREA, SEAICEadvSchArea,
     I         uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj),
     I         uTrans, vTrans, AREA(1-OLx,1-OLy,bi,bj), recip_heff,
     O         gFld, afx, afy,
     I         bi, bj, myTime, myIter, myThid )
         IF ( SEAICEdiffKhArea .GT. 0. _d 0 ) THEN
C-    Add tendency due to diffusion
          CALL SEAICE_DIFFUSION(
     I         GAD_AREA, SEAICEdiffKhArea, ONE,
     I         AREA(1-OLx,1-OLy,bi,bj), HEFFM,
     I         xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
     U         gFld,
     I         bi, bj, myTime, myIter, myThid )
         ENDIF
C     now do the "explicit" time step
         DO j=1,sNy
          DO i=1,sNx
           AREA(i,j,bi,bj) = HEFFM(i,j,bi,bj) * (
     &          AREA(i,j,bi,bj) + SEAICE_deltaTtherm * gFld(i,j)
     &          )
          ENDDO
         ENDDO
#ifdef SEAICE_ITD
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           AREAITD(i,j,it,bi,bj)=AREA(i,j,bi,bj)
          ENDDO
         ENDDO
        ENDDO
C     open water fraction needs to be advected for the ridging scheme
        CALL SEAICE_ADVECTION(
     I       GAD_AREA, SEAICEadvSchArea,
     I       uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj),
     I       uTrans, vTrans, opnWtrFrac(1-OLx,1-OLy,bi,bj), recip_heff,
     O       gFld, afx, afy,
     I       bi, bj, myTime, myIter, myThid )
        IF ( SEAICEdiffKhArea .GT. 0. _d 0 ) THEN
C-    Add tendency due to diffusion
         CALL SEAICE_DIFFUSION(
     I        GAD_AREA, SEAICEdiffKhArea, ONE,
     I        opnWtrFrac(1-OLx,1-OLy,bi,bj), HEFFM,
     I        xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
     U        gFld,
     I        bi, bj, myTime, myIter, myThid )
        ENDIF
C     now do the "explicit" time step
        DO j=1,sNy
         DO i=1,sNx
          opnWtrFrac(i,j,bi,bj) = HEFFM(i,j,bi,bj) * (
     &         opnWtrFrac(i,j,bi,bj) + SEAICE_deltaTtherm * gFld(i,j)
     &         )
         ENDDO
        ENDDO
#endif /* SEAICE_ITD */
        ENDIF

C--   Effective Snow Thickness (Volume)
        IF ( SEAICEadvSnow ) THEN
#ifdef SEAICE_ITD
        DO it=1,SEAICE_multDim
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           HSNOW(i,j,bi,bj)=HSNOWITD(i,j,it,bi,bj)
          ENDDO
         ENDDO
#endif /* SEAICE_ITD */
          CALL SEAICE_ADVECTION(
     I         GAD_SNOW, SEAICEadvSchSnow,
     I         uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj),
     I         uTrans, vTrans, HSNOW(1-OLx,1-OLy,bi,bj), recip_heff,
     O         gFld, afx, afy,
     I         bi, bj, myTime, myIter, myThid )
         IF ( SEAICEdiffKhSnow .GT. 0. _d 0 ) THEN
C--   Add tendency due to diffusion
          CALL SEAICE_DIFFUSION(
     I         GAD_SNOW, SEAICEdiffKhSnow, ONE,
     I         HSNOW(1-OLx,1-OLy,bi,bj), HEFFM,
     I         xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
     U         gFld,
     I         bi, bj, myTime, myIter, myThid )
         ENDIF
C     now do the "explicit" time step
         DO j=1,sNy
          DO i=1,sNx
           HSNOW(i,j,bi,bj) = HEFFM(i,j,bi,bj) * (
     &          HSNOW(i,j,bi,bj) + SEAICE_deltaTtherm * gFld(i,j)
     &          )
          ENDDO
         ENDDO
#ifdef SEAICE_ITD
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
           HSNOWITD(i,j,it,bi,bj)=HSNOW(i,j,bi,bj)
          ENDDO
         ENDDO
        ENDDO
#endif /* SEAICE_ITD */
        ENDIF

#ifdef SEAICE_ITD
C     update mean ice thickness HEFF and total ice concentration AREA
C     to match single category values
C     (necessary here because updated HEFF is used below for SItracer)
        CALL SEAICE_ITD_SUM(bi, bj, myTime, myIter, myThid)
#endif /* SEAICE_ITD */

#ifdef SEAICE_VARIABLE_SALINITY
C--   Effective Sea Ice Salinity (Mass of salt)
        IF ( SEAICEadvSalt ) THEN
          CALL SEAICE_ADVECTION(
     I         GAD_SALT, SEAICEadvSchSalt,
     I         uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj),
     I         uTrans, vTrans, HSALT(1-OLx,1-OLy,bi,bj), recip_heff,
     O         gFld, afx, afy,
     I         bi, bj, myTime, myIter, myThid )
         IF ( SEAICEdiffKhSalt .GT. 0. _d 0 ) THEN
C--   Add tendency due to diffusion
          CALL SEAICE_DIFFUSION(
     I         GAD_SALT, SEAICEdiffKhSalt, ONE,
     I         HSALT(1-OLx,1-OLy,bi,bj), HEFFM,
     I         xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
     U         gFld,
     I         bi, bj, myTime, myIter, myThid )
         ENDIF
C     now do the "explicit" time step
         DO j=1,sNy
          DO i=1,sNx
           HSALT(i,j,bi,bj) = HEFFM(i,j,bi,bj) * (
     &          HSALT(i,j,bi,bj) + SEAICE_deltaTtherm * gFld(i,j)
     &          )
          ENDDO
         ENDDO
        ENDIF
#endif /* SEAICE_VARIABLE_SALINITY */

#ifdef ALLOW_SITRACER
C--   Sea Ice Tracers
        DO iTr = 1, SItrNumInUse
        IF ( (SEAICEadvHEFF.AND.(SItrMate(iTr).EQ.'HEFF')).OR.
     &       (SEAICEadvAREA.AND.(SItrMate(iTr).EQ.'AREA')) ) THEN
C--   scale to effective value
         IF (SItrMate(iTr).EQ.'HEFF') THEN
          SEAICEadvSchSItr=SEAICEadvSchHEFF
          SEAICEdiffKhSItr=SEAICEdiffKhHEFF
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
            SItrExt(i,j,bi,bj) = HEFFM(i,j,bi,bj) *
     &          SItracer(i,j,bi,bj,iTr) * hEffNm1(i,j,bi,bj)
           ENDDO
          ENDDO
c TAF?   ELSEIF (SItrMate(iTr).EQ.'AREA') THEN
         ELSE
          SEAICEadvSchSItr=SEAICEadvSchAREA
          SEAICEdiffKhSItr=SEAICEdiffKhAREA
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
            SItrExt(i,j,bi,bj) = HEFFM(i,j,bi,bj) *
     &          SItracer(i,j,bi,bj,iTr) * areaNm1(i,j,bi,bj)
           ENDDO
          ENDDO
         ENDIF
C--   store a couple things
          DO j=1-OLy,sNy+OLy
           DO i=1-OLx,sNx+OLx
#ifdef ALLOW_SITRACER_ADVCAP
C--   store previous value for spurious maxima treament
            SItrPrev(i,j,bi,bj)=SItracer(i,j,bi,bj,iTr)
#endif
#ifdef ALLOW_SITRACER_DEBUG_DIAG
            diagArray(I,J,2+(iTr-1)*5) = SItrExt(i,j,bi,bj)
#endif
           ENDDO
          ENDDO
C--   compute advective tendency
          CALL SEAICE_ADVECTION(
     I         GAD_SITR+iTr-1, SEAICEadvSchSItr,
     I         uc(1-OLx,1-OLy,bi,bj), vc(1-OLx,1-OLy,bi,bj),
     I         uTrans, vTrans, SItrExt(1-OLx,1-OLy,bi,bj),
     I         recip_heff,
     O         gFld, afx, afy,
     I         bi, bj, myTime, myIter, myThid )
          IF ( SEAICEdiffKhHeff .GT. 0. _d 0 ) THEN
C--   add diffusive tendency
          CALL SEAICE_DIFFUSION(
     I         GAD_SITR+iTr-1, SEAICEdiffKhSItr, ONE,
     I         SItrExt(1-OLx,1-OLy,bi,bj), HEFFM,
     I         xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
     U         gFld,
     I         bi, bj, myTime, myIter, myThid )
          ENDIF
C--   apply tendency
          DO j=1,sNy
           DO i=1,sNx
            SItrExt(i,j,bi,bj) = HEFFM(i,j,bi,bj) * (
     &        SItrExt(i,j,bi,bj) + SEAICE_deltaTtherm * gFld(i,j) )
           ENDDO
          ENDDO
C--   scale back to actual value, or move effective value to ocean bucket
         IF (SItrMate(iTr).EQ.'HEFF') THEN
          DO j=1,sNy
           DO i=1,sNx
            if (HEFF(I,J,bi,bj).GE.siEps) then
            SItracer(i,j,bi,bj,iTr)=SItrExt(i,j,bi,bj)/HEFF(I,J,bi,bj)
            SItrBucket(i,j,bi,bj,iTr)=0. _d 0
            else
            SItracer(i,j,bi,bj,iTr)=0. _d 0
            SItrBucket(i,j,bi,bj,iTr)=SItrExt(i,j,bi,bj)
            endif
#ifdef ALLOW_SITRACER_ADVCAP
C hack to try avoid 'spontaneous generation' of maxima, which supposedly would
C occur less frequently if we advected SItr with uXheff instead SItrXheff with u
             tmpscal1=max(SItrPrev(i,j,bi,bj),
     &       SItrPrev(i+1,j,bi,bj),SItrPrev(i-1,j,bi,bj),
     &       SItrPrev(i,j+1,bi,bj),SItrPrev(i,j-1,bi,bj))
             tmpscal2=MAX(ZERO,SItracer(i,j,bi,bj,iTr)-tmpscal1)
             SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)-tmpscal2
             SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr)
     &           +tmpscal2*HEFF(I,J,bi,bj)
#endif
C           treat case of potential negative value
            if (HEFF(I,J,bi,bj).GE.siEps) then
              tmpscal1=MIN(0. _d 0,SItracer(i,j,bi,bj,iTr))
              SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)-tmpscal1
              SItrBucket(i,j,bi,bj,iTr)=SItrBucket(i,j,bi,bj,iTr)
     &                                 +HEFF(I,J,bi,bj)*tmpscal1
            endif
#ifdef ALLOW_SITRACER_DEBUG_DIAG
      diagArray(I,J,1+(iTr-1)*5)= - SItrBucket(i,j,bi,bj,iTr)
     &  *HEFFM(I,J,bi,bj)/SEAICE_deltaTtherm*SEAICE_rhoIce
      tmpscal1= ( HEFF(I,J,bi,bj)*SItracer(i,j,bi,bj,iTr)
     &  + SItrBucket(i,j,bi,bj,iTr) )*HEFFM(I,J,bi,bj)
      diagArray(I,J,2+(iTr-1)*5)= tmpscal1-diagArray(I,J,2+(iTr-1)*5)
      diagArray(I,J,3+(iTr-1)*5)=HEFFM(i,j,bi,bj) *
     &  SEAICE_deltaTtherm * gFld(i,j)
#endif
           ENDDO
          ENDDO
c TAF?   ELSEIF (SItrMate(iTr).EQ.'AREA') THEN
         ELSE
          DO j=1,sNy
           DO i=1,sNx
            if (AREA(I,J,bi,bj).GE.SEAICE_area_floor) then
            SItracer(i,j,bi,bj,iTr)=SItrExt(i,j,bi,bj)/AREA(I,J,bi,bj)
            else
            SItracer(i,j,bi,bj,iTr)=0. _d 0
            endif
            SItrBucket(i,j,bi,bj,iTr)=0. _d 0
#ifdef ALLOW_SITRACER_ADVCAP
             tmpscal1=max(SItrPrev(i,j,bi,bj),
     &       SItrPrev(i+1,j,bi,bj),SItrPrev(i-1,j,bi,bj),
     &       SItrPrev(i,j+1,bi,bj),SItrPrev(i,j-1,bi,bj))
             tmpscal2=MAX(ZERO,SItracer(i,j,bi,bj,iTr)-tmpscal1)
             SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)-tmpscal2
#endif
C           treat case of potential negative value
            if (AREA(I,J,bi,bj).GE.SEAICE_area_floor) then
              tmpscal1=MIN(0. _d 0,SItracer(i,j,bi,bj,iTr))
              SItracer(i,j,bi,bj,iTr)=SItracer(i,j,bi,bj,iTr)-tmpscal1
            endif
#ifdef ALLOW_SITRACER_DEBUG_DIAG
      diagArray(I,J,1+(iTr-1)*5)= 0. _d 0
      diagArray(I,J,2+(iTr-1)*5)= - diagArray(I,J,2+(iTr-1)*5)
     & + AREA(I,J,bi,bj)*SItracer(i,j,bi,bj,iTr)*HEFFM(I,J,bi,bj)
      diagArray(I,J,3+(iTr-1)*5)=HEFFM(i,j,bi,bj) *
     &  SEAICE_deltaTtherm * gFld(i,j)
#endif
           ENDDO
          ENDDO
         ENDIF
C--
         ENDIF
        ENDDO
#ifdef ALLOW_SITRACER_DEBUG_DIAG
c     CALL DIAGNOSTICS_FILL(DIAGarray,'UDIAG2  ',0,Nr,2,bi,bj,myThid)
#endif
#endif /* ALLOW_SITRACER */

C---   end bi,bj loops
       ENDDO
      ENDDO

      ELSE
C--   if not multiDimAdvection

       IF ( SEAICEadvHEff ) THEN
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE heff = comlev1, key = ikey_dynamics, kind=isbyte
#endif
#ifdef SEAICE_ITD
        DO it=1,SEAICE_multDim
         DO bj=myByLo(myThid),myByHi(myThid)
          DO bi=myBxLo(myThid),myBxHi(myThid)
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
             HEFF(i,j,bi,bj)=HEFFITD(i,j,it,bi,bj)
            ENDDO
           ENDDO
          ENDDO
         ENDDO
#endif /* SEAICE_ITD */
        CALL ADVECT( uc, vc, hEff, fldNm1, HEFFM, myThid )
        IF ( SEAICEdiffKhHeff .GT. 0. _d 0 ) THEN
C-    Add tendency due to diffusion
         DO bj=myByLo(myThid),myByHi(myThid)
          DO bi=myBxLo(myThid),myBxHi(myThid)
            CALL SEAICE_DIFFUSION(
     I           GAD_HEFF, SEAICEdiffKhHeff, SEAICE_deltaTtherm,
     I           fldNm1(1-OLx,1-OLy,bi,bj), HEFFM,
     I           xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
     U           HEFF(1-OLx,1-OLy,bi,bj),
     I           bi, bj, myTime, myIter, myThid )
          ENDDO
         ENDDO
        ENDIF
#ifdef SEAICE_ITD
         DO bj=myByLo(myThid),myByHi(myThid)
          DO bi=myBxLo(myThid),myBxHi(myThid)
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
             HEFFITD(i,j,it,bi,bj)=HEFF(i,j,bi,bj)
            ENDDO
           ENDDO
          ENDDO
         ENDDO
        ENDDO
#endif /* SEAICE_ITD */
       ENDIF
       IF ( SEAICEadvArea ) THEN
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE area = comlev1, key = ikey_dynamics, kind=isbyte
#endif
#ifdef SEAICE_ITD
        DO it=1,SEAICE_multDim
         DO bj=myByLo(myThid),myByHi(myThid)
          DO bi=myBxLo(myThid),myBxHi(myThid)
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
             AREA(i,j,bi,bj)=AREAITD(i,j,it,bi,bj)
            ENDDO
           ENDDO
          ENDDO
         ENDDO
#endif /* SEAICE_ITD */
        CALL ADVECT( uc, vc, area, fldNm1, HEFFM, myThid )
        IF ( SEAICEdiffKhArea .GT. 0. _d 0 ) THEN
C-    Add tendency due to diffusion
         DO bj=myByLo(myThid),myByHi(myThid)
          DO bi=myBxLo(myThid),myBxHi(myThid)
            CALL SEAICE_DIFFUSION(
     I           GAD_AREA, SEAICEdiffKhArea, SEAICE_deltaTtherm,
     I           fldNm1(1-OLx,1-OLy,bi,bj), HEFFM,
     I           xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
     U           Area(1-OLx,1-OLy,bi,bj),
     I           bi, bj, myTime, myIter, myThid )
          ENDDO
         ENDDO
        ENDIF
#ifdef SEAICE_ITD
         DO bj=myByLo(myThid),myByHi(myThid)
          DO bi=myBxLo(myThid),myBxHi(myThid)
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
             AREAITD(i,j,it,bi,bj)=AREA(i,j,bi,bj)
            ENDDO
           ENDDO
          ENDDO
         ENDDO
        ENDDO
C     open water fraction needs to be advected for the ridging scheme
        CALL ADVECT( uc, vc, opnWtrFrac, fldNm1, HEFFM, myThid )
        IF ( SEAICEdiffKhArea .GT. 0. _d 0 ) THEN
C-    Add tendency due to diffusion
         DO bj=myByLo(myThid),myByHi(myThid)
          DO bi=myBxLo(myThid),myBxHi(myThid)
            CALL SEAICE_DIFFUSION(
     I           GAD_AREA, SEAICEdiffKhArea, SEAICE_deltaTtherm,
     I           fldNm1(1-OLx,1-OLy,bi,bj), HEFFM,
     I           xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
     U           opnWtrFrac(1-OLx,1-OLy,bi,bj),
     I           bi, bj, myTime, myIter, myThid )
          ENDDO
         ENDDO
        ENDIF
#endif /* SEAICE_ITD */
       ENDIF
       IF ( SEAICEadvSnow ) THEN
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE hsnow = comlev1, key = ikey_dynamics, kind=isbyte
#endif
#ifdef SEAICE_ITD
        DO it=1,SEAICE_multDim
         DO bj=myByLo(myThid),myByHi(myThid)
          DO bi=myBxLo(myThid),myBxHi(myThid)
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
             HSNOW(i,j,bi,bj)=HSNOWITD(i,j,it,bi,bj)
            ENDDO
           ENDDO
          ENDDO
         ENDDO
#endif /* SEAICE_ITD */
        CALL ADVECT( uc, vc, HSNOW, fldNm1, HEFFM, myThid )
        IF ( SEAICEdiffKhSnow .GT. 0. _d 0 ) THEN
C-    Add tendency due to diffusion
         DO bj=myByLo(myThid),myByHi(myThid)
          DO bi=myBxLo(myThid),myBxHi(myThid)
            CALL SEAICE_DIFFUSION(
     I           GAD_SNOW, SEAICEdiffKhSnow, SEAICE_deltaTtherm,
     I           fldNm1(1-OLx,1-OLy,bi,bj), HEFFM,
     I           xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
     U           HSNOW(1-OLx,1-OLy,bi,bj),
     I           bi, bj, myTime, myIter, myThid )
          ENDDO
         ENDDO
        ENDIF
#ifdef SEAICE_ITD
         DO bj=myByLo(myThid),myByHi(myThid)
          DO bi=myBxLo(myThid),myBxHi(myThid)
           DO j=1-OLy,sNy+OLy
            DO i=1-OLx,sNx+OLx
             HSNOWITD(i,j,it,bi,bj)=HSNOW(i,j,bi,bj)
            ENDDO
           ENDDO
          ENDDO
         ENDDO
        ENDDO
#endif /* SEAICE_ITD */
       ENDIF

#ifdef SEAICE_ITD
       DO bj=myByLo(myThid),myByHi(myThid)
        DO bi=myBxLo(myThid),myBxHi(myThid)
C     update mean ice thickness HEFF and total ice concentration AREA
C     to match single category values
C     (necessary here because updated HEFF is used below for SItracer)
         CALL SEAICE_ITD_SUM(bi, bj, myTime, myIter, myThid)
        ENDDO
       ENDDO
#endif /* SEAICE_ITD */

#ifdef SEAICE_VARIABLE_SALINITY
       IF ( SEAICEadvSalt ) THEN
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE hsalt = comlev1, key = ikey_dynamics, kind=isbyte
#endif
        CALL ADVECT( uc, vc, HSALT, fldNm1, HEFFM, myThid )
        IF ( SEAICEdiffKhSalt .GT. 0. _d 0 ) THEN
C-    Add tendency due to diffusion
         DO bj=myByLo(myThid),myByHi(myThid)
          DO bi=myBxLo(myThid),myBxHi(myThid)
            CALL SEAICE_DIFFUSION(
     I           GAD_SALT, SEAICEdiffKhSalt, SEAICE_deltaTtherm,
     I           fldNm1(1-OLx,1-OLy,bi,bj), HEFFM,
     I           xA(1-OLx,1-OLy,bi,bj), yA(1-OLx,1-OLy,bi,bj),
     U           HSALT(1-OLx,1-OLy,bi,bj),
     I           bi, bj, myTime, myIter, myThid )
          ENDDO
         ENDDO
        ENDIF
       ENDIF
#endif /* SEAICE_VARIABLE_SALINITY */

C--   end if multiDimAdvection
      ENDIF

      RETURN
      END