C $Header: /u/gcmpack/MITgcm/pkg/thsice/thsice_advection.F,v 1.16 2011/06/07 22:26:37 jmc Exp $
C $Name:  $

#include "THSICE_OPTIONS.h"
#ifdef ALLOW_GENERIC_ADVDIFF
# include "GAD_OPTIONS.h"
#endif

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C !ROUTINE: THSICE_ADVECTION

C !INTERFACE: ==========================================================
      SUBROUTINE THSICE_ADVECTION(
     I     tracerIdentity,
     I     advectionScheme,
     I     useGridArea,
     I     uTrans, vTrans, maskOc, deltaTadvect, iceEps,
     U     iceVol, iceFld,
     O     afx, afy,
     I     bi, bj, myTime, myIter, myThid)

C !DESCRIPTION:
C Calculates the tendency of a sea-ice field due to advection.
C It uses the multi-dimensional method given in \ref{sect:multiDimAdvection}
C and can only be used for the non-linear advection schemes such as the
C direct-space-time method and flux-limiters.
C
C This routine is an adaption of the GAD_ADVECTION for 2D-fields.
C for Area, effective thickness or other sea-ice field,
C  the contribution iceFld*div(u) (that is present in gad_advection)
C  is not included here.
C
C The algorithm is as follows:
C \begin{itemize}
C \item{$\theta^{(n+1/2)} = \theta^{(n)}
C      - \Delta t \partial_x (u\theta^{(n)}) + \theta^{(n)} \partial_x u$}
C \item{$\theta^{(n+2/2)} = \theta^{(n+1/2)}
C      - \Delta t \partial_y (v\theta^{(n+1/2)}) + \theta^{(n)} \partial_y v$}
C \item{$G_\theta = ( \theta^{(n+2/2)} - \theta^{(n)} )/\Delta t$}
C \end{itemize}
C
C The tendency (output) is over-written by this routine.

C !USES: ===============================================================
      IMPLICIT NONE
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
#if ( defined ALLOW_DBUG_THSICE || defined ALLOW_AUTODIFF_TAMC )
# include "THSICE_SIZE.h"
#endif
#ifdef ALLOW_GENERIC_ADVDIFF
# include "GAD.h"
#endif
#ifdef ALLOW_EXCH2
#include "W2_EXCH2_SIZE.h"
#include "W2_EXCH2_TOPOLOGY.h"
#endif /* ALLOW_EXCH2 */
#ifdef ALLOW_AUTODIFF_TAMC
# include "tamc.h"
# include "tamc_keys.h"
#endif

C !INPUT PARAMETERS: ===================================================
C  tracerIdentity  :: tracer identifier (required only for OBCS)
C  advectionScheme :: advection scheme to use (Horizontal plane)
C  useGridArea     :: use grid-cell Area (instead of "iceVol" field)
C  uTrans,vTrans   :: volume transports at U,V points
C  maskOc          :: oceanic mask
C  iceVol          :: sea-ice volume
C  iceFld          :: sea-ice field
C  deltaTadvect    :: time-step used for advection [s]
C  iceEps          :: small volume (to avoid division by zero if iceVol==0)
C  bi,bj           :: tile indices
C  myTime          :: current time in simulation [s]
C  myIter          :: current iteration number
C  myThid          :: my thread Id. number
      INTEGER tracerIdentity
      INTEGER advectionScheme
      LOGICAL useGridArea
      _RL uTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL vTrans(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RS maskOc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL iceFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL iceVol(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL deltaTadvect, iceEps
      INTEGER bi,bj
      _RL myTime
      INTEGER myIter
      INTEGER myThid

C !OUTPUT PARAMETERS: ==================================================
C  iceVol (Updated):: sea-ice volume
C  iceFld (Updated):: sea-ice field
C  afx             :: horizontal advective flux, x direction
C  afy             :: horizontal advective flux, y direction
      _RL afx   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL afy   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)

#ifdef ALLOW_GENERIC_ADVDIFF
C !LOCAL VARIABLES: ====================================================
C  maskLocC      :: 2-D array mask at grid-cell center
C  maskLocW      :: 2-D array for mask at West points
C  maskLocS      :: 2-D array for mask at South points
C  iMin,iMax,    :: loop range for called routines
C  jMin,jMax     :: loop range for called routines
C [iMin,iMax]Upd :: loop range to update sea-ice field
C [jMin,jMax]Upd :: loop range to update sea-ice field
C  i,j           :: loop indices
C  uCfl          :: CFL number, zonal direction
C  vCfl          :: CFL number, meridional direction
C  af            :: 2-D array for horizontal advective flux
C  calc_fluxes_X :: logical to indicate to calculate fluxes in X dir
C  calc_fluxes_Y :: logical to indicate to calculate fluxes in Y dir
C  interiorOnly  :: only update the interior of myTile, but not the edges
C  overlapOnly   :: only update the edges of myTile, but not the interior
C  nipass        :: number of passes in multi-dimensional method
C  ipass         :: number of the current pass being made
C  myTile        :: variables used to determine which cube face
C  nCFace        :: owns a tile for cube grid runs using
C                :: multi-dim advection.
C [N,S,E,W]_edge :: true if N,S,E,W edge of myTile is an Edge of the cube
      _RS maskLocC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RS maskLocW(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RS maskLocS(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      INTEGER iMin,iMax,jMin,jMax
      INTEGER iMinUpd,iMaxUpd,jMinUpd,jMaxUpd
      INTEGER i,j,k
      _RL uCfl    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL vCfl    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL af      (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
      _RL tmpVol
      LOGICAL calc_fluxes_X, calc_fluxes_Y, withSigns
      LOGICAL interiorOnly, overlapOnly
      INTEGER nipass,ipass
      INTEGER nCFace
      LOGICAL N_edge, S_edge, E_edge, W_edge
#ifdef ALLOW_EXCH2
      INTEGER myTile
#endif
#ifdef ALLOW_DBUG_THSICE
      LOGICAL dBugFlag
      INTEGER idb,jdb,biDb
      _RL tmpFac
#endif /* ALLOW_DBUG_THSICE */
CEOP

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

      k = 1
#ifdef ALLOW_AUTODIFF_TAMC
      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
      ticekey = (act1 + 1) + act2*max1
     &                     + act3*max1*max2
     &                     + act4*max1*max2*max3
#endif /* ALLOW_AUTODIFF_TAMC */

#ifdef ALLOW_DBUG_THSICE
      dBugFlag = debugLevel.GE.debLevC
     &     .AND. myIter.EQ.nIter0
     &     .AND. ( tracerIdentity.EQ.GAD_SI_HICE .OR.
     &             tracerIdentity.EQ.GAD_SI_QICE2 )
c    &     .AND. tracerIdentity.EQ.GAD_SI_FRAC
      idb  = MIN(13,sNx)
      jdb  = MOD(17,sNy)
      biDb = nSx/2
#endif /* ALLOW_DBUG_THSICE */

C--   Set up work arrays with valid (i.e. not NaN) values
C     These inital values do not alter the numerical results. They
C     just ensure that all memory references are to valid floating
C     point numbers. This prevents spurious hardware signals due to
C     uninitialised but inert locations.

C--   Set tile-specific parameters for horizontal fluxes
      IF (useCubedSphereExchange) THEN
       nipass=3
#ifdef ALLOW_EXCH2
       myTile = W2_myTileList(bi,bj)
       nCFace = exch2_myFace(myTile)
       N_edge = exch2_isNedge(myTile).EQ.1
       S_edge = exch2_isSedge(myTile).EQ.1
       E_edge = exch2_isEedge(myTile).EQ.1
       W_edge = exch2_isWedge(myTile).EQ.1
#else
       nCFace = bi
       N_edge = .TRUE.
       S_edge = .TRUE.
       E_edge = .TRUE.
       W_edge = .TRUE.
#endif
      ELSE
       nipass=2
       nCFace = bi
       N_edge = .FALSE.
       S_edge = .FALSE.
       E_edge = .FALSE.
       W_edge = .FALSE.
      ENDIF

      iMin = 1-OLx
      iMax = sNx+OLx
      jMin = 1-OLy
      jMax = sNy+OLy

C--   Start horizontal fluxes

C--   set mask West & South (and local Centered mask)
      DO j=1-OLy,sNy+OLy
       maskLocW(1-Olx,j) = 0.
       DO i=2-OLx,sNx+OLx
        maskLocW(i,j) = MIN( maskOc(i-1,j), maskOc(i,j) )
#ifdef ALLOW_OBCS
     &                * maskInW(i,j,bi,bj)
#endif /* ALLOW_OBCS */
       ENDDO
      ENDDO
      DO i=1-OLx,sNx+OLx
       maskLocS(i,1-Oly) = 0.
      ENDDO
      DO j=2-OLy,sNy+OLy
       DO i=1-OLx,sNx+OLx
        maskLocS(i,j) = MIN( maskOc(i,j-1), maskOc(i,j) )
#ifdef ALLOW_OBCS
     &                * maskInS(i,j,bi,bj)
#endif /* ALLOW_OBCS */
       ENDDO
      ENDDO
C     maskLocC is just a local copy of Ocean mask (maskOc) except if using OBCS:
C     use "maksInC" to prevent updating tracer field in OB regions
      DO j=1-OLy,sNy+OLy
       DO i=1-OLx,sNx+OLx
        maskLocC(i,j) = maskOc(i,j)
#ifdef ALLOW_OBCS
     &                * maskInC(i,j,bi,bj)
#endif /* ALLOW_OBCS */
       ENDDO
      ENDDO

cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
      IF (useCubedSphereExchange) THEN
       withSigns = .FALSE.
       CALL FILL_CS_CORNER_UV_RS(
     &      withSigns, maskLocW,maskLocS, bi,bj, myThid )
      ENDIF
cph-exch2#endif

C--   Multiple passes for different directions on different tiles
C--   For cube need one pass for each of red, green and blue axes.
      DO ipass=1,nipass
#ifdef ALLOW_AUTODIFF_TAMC
         ikey_4 = ipass
     &         + nipass*act1
     &         + nipass*max1*act2
     &         + nipass*max1*max2*act3
     &         + nipass*max1*max2*max3*act4
#endif

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE iceFld(:,:) = comlev1_thsice_4, key=ikey_4, byte=isbyte
CADJ STORE iceVol(:,:) = comlev1_thsice_4, key=ikey_4, byte=isbyte
CADJ STORE af(:,:)     = comlev1_thsice_4, key=ikey_4, byte=isbyte
#endif

       interiorOnly = .FALSE.
       overlapOnly  = .FALSE.
       IF (useCubedSphereExchange) THEN
C--   CubedSphere : pass 3 times, with partial update of local seaice field
        IF (ipass.EQ.1) THEN
         overlapOnly  = MOD(nCFace,3).EQ.0
         interiorOnly = MOD(nCFace,3).NE.0
         calc_fluxes_X = nCFace.EQ.6 .OR. nCFace.EQ.1 .OR. nCFace.EQ.2
         calc_fluxes_Y = nCFace.EQ.3 .OR. nCFace.EQ.4 .OR. nCFace.EQ.5
        ELSEIF (ipass.EQ.2) THEN
         overlapOnly  = MOD(nCFace,3).EQ.2
         interiorOnly = MOD(nCFace,3).EQ.1
         calc_fluxes_X = nCFace.EQ.2 .OR. nCFace.EQ.3 .OR. nCFace.EQ.4
         calc_fluxes_Y = nCFace.EQ.5 .OR. nCFace.EQ.6 .OR. nCFace.EQ.1
        ELSE
         interiorOnly = .TRUE.
         calc_fluxes_X = nCFace.EQ.5 .OR. nCFace.EQ.6
         calc_fluxes_Y = nCFace.EQ.2 .OR. nCFace.EQ.3
        ENDIF
       ELSE
C--   not CubedSphere
        calc_fluxes_X = MOD(ipass,2).EQ.1
        calc_fluxes_Y = .NOT.calc_fluxes_X
       ENDIF
#ifdef ALLOW_DBUG_THSICE
       IF (dBugFlag.AND.bi.EQ.biDb ) WRITE(6,'(A,3I4,2I5,4L5)')
     &   'ICE_adv:', tracerIdentity, ipass, bi, idb, jdb,
     &   calc_fluxes_X, calc_fluxes_Y, overlapOnly, interiorOnly
#endif /* ALLOW_DBUG_THSICE */

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

       IF (calc_fluxes_X) THEN

C-     Do not compute fluxes if
C       a) needed in overlap only
C   and b) the overlap of myTile are not cube-face Edges
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE iceFld(:,:) = comlev1_thsice_4, key=ikey_4, byte=isbyte
CADJ STORE iceVol(:,:) = comlev1_thsice_4, key=ikey_4, byte=isbyte
CADJ STORE af(:,:)     = comlev1_thsice_4, key=ikey_4, byte=isbyte
#endif
        IF ( .NOT.overlapOnly .OR. N_edge .OR. S_edge ) THEN

C-     Advective flux in X
         DO j=1-Oly,sNy+Oly
          DO i=1-Olx,sNx+Olx
            af(i,j) = 0.
          ENDDO
         ENDDO

cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
C-     Internal exchange for calculations in X
         IF ( overlapOnly ) THEN
          CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
     &                               iceFld, bi,bj, myThid )
          IF ( .NOT.useGridArea )
     &    CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
     &                               iceVol, bi,bj, myThid )
         ENDIF
cph-exch2#endif

C-     Compute CFL number
         IF ( useGridArea ) THEN
          DO j=1-Oly,sNy+Oly
           DO i=2-Olx,sNx+Olx
            uCfl(i,j) = deltaTadvect*(
     &          MAX( uTrans(i,j), 0. _d 0 )*recip_rA(i-1,j,bi,bj)
     &         +MAX(-uTrans(i,j), 0. _d 0 )*recip_rA( i ,j,bi,bj)
     &                               )
           ENDDO
          ENDDO
         ELSE
          DO j=1-Oly,sNy+Oly
           DO i=2-Olx,sNx+Olx
            uCfl(i,j) = deltaTadvect*(
     &        MAX( uTrans(i,j), 0. _d 0 )/MAX( iceVol(i-1,j), iceEps )
     &       +MAX(-uTrans(i,j), 0. _d 0 )/MAX( iceVol( i ,j), iceEps )
     &                               )
           ENDDO
          ENDDO
         ENDIF

         IF ( advectionScheme.EQ.ENUM_UPWIND_1RST
     &        .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
          CALL GAD_DST2U1_ADV_X( bi,bj,k, advectionScheme, .FALSE.,
     I             deltaTadvect, uTrans, uCfl, iceFld,
     O             af, myThid )
#ifdef ALLOW_DBUG_THSICE
          IF (dBugFlag.AND.bi.EQ.biDb ) WRITE(6,'(A,1P4E14.6)')
     &      'ICE_adv: xFx=', af(idb,jdb), iceFld(idb,jdb),
     &       uTrans(idb,jdb), af(idb,jdb)/uTrans(idb,jdb)
#endif /* ALLOW_DBUG_THSICE */
         ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
          CALL GAD_FLUXLIMIT_ADV_X( bi,bj,k, .FALSE., deltaTadvect,
     I             uTrans, uCfl, maskLocW, iceFld,
     O             af, myThid )
         ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
          CALL GAD_DST3_ADV_X(      bi,bj,k, .FALSE., deltaTadvect,
     I             uTrans, uCfl, maskLocW, iceFld,
     O             af, myThid )
         ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
          CALL GAD_DST3FL_ADV_X(    bi,bj,k, .FALSE., deltaTadvect,
     I             uTrans, uCfl, maskLocW, iceFld,
     O             af, myThid )
         ELSE
          STOP
     & 'THSICE_ADVECTION: adv. scheme incompatibale with multi-dim'
         ENDIF

cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
C--   Internal exchange for next calculations in Y
         IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
          CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
     &                               iceFld, bi,bj, myThid )
          IF ( .NOT.useGridArea )
     &    CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
     &                               iceVol, bi,bj, myThid )
         ENDIF
cph-exch2#endif

C--   Advective flux in X : done
        ENDIF

C-     Update the local seaice field where needed:

C     update in overlap-Only
        IF ( overlapOnly ) THEN
         iMinUpd = 1-OLx+1
         iMaxUpd = sNx+OLx-1
C--   notes: these 2 lines below have no real effect (because recip_hFac=0
C            in corner region) but safer to keep them.
         IF ( W_edge ) iMinUpd = 1
         IF ( E_edge ) iMaxUpd = sNx

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE iceFld(:,:) = comlev1_thsice_4, key=ikey_4, byte=isbyte
CADJ STORE iceVol(:,:) = comlev1_thsice_4, key=ikey_4, byte=isbyte
#endif
         IF ( S_edge ) THEN
          IF ( useGridArea ) THEN
           DO j=1-OLy,0
            DO i=iMinUpd,iMaxUpd
             iceFld(i,j) = iceFld(i,j)
     &                    -deltaTadvect*maskLocC(i,j)
     &                                 *recip_rA(i,j,bi,bj)
     &                                 *( af(i+1,j)-af(i,j) )
            ENDDO
           ENDDO
          ELSE
           DO j=1-OLy,0
            DO i=iMinUpd,iMaxUpd
             tmpVol = iceVol(i,j)
             iceVol(i,j) = iceVol(i,j)
     &                    -deltaTadvect*maskLocC(i,j)
     &                           *( uTrans(i+1,j)-uTrans(i,j) )
             IF ( iceVol(i,j).GT.iceEps )
     &       iceFld(i,j) = ( iceFld(i,j)*tmpVol
     &                      -deltaTadvect*maskLocC(i,j)
     &                                  *( af(i+1,j)-af(i,j) )
     &                     )/iceVol(i,j)
            ENDDO
           ENDDO
          ENDIF
C-    keep advective flux (for diagnostics)
          DO j=1-OLy,0
           DO i=1-OLx+1,sNx+OLx
             afx(i,j) = af(i,j)
           ENDDO
          ENDDO
C-    end if South Edge
         ENDIF
         IF ( N_edge ) THEN
          IF ( useGridArea ) THEN
           DO j=sNy+1,sNy+OLy
            DO i=iMinUpd,iMaxUpd
             iceFld(i,j) = iceFld(i,j)
     &                    -deltaTadvect*maskLocC(i,j)
     &                                 *recip_rA(i,j,bi,bj)
     &                                 *( af(i+1,j)-af(i,j) )
            ENDDO
           ENDDO
          ELSE
           DO j=sNy+1,sNy+OLy
            DO i=iMinUpd,iMaxUpd
             tmpVol = iceVol(i,j)
             iceVol(i,j) = iceVol(i,j)
     &                    -deltaTadvect*maskLocC(i,j)
     &                           *( uTrans(i+1,j)-uTrans(i,j) )
             IF ( iceVol(i,j).GT.iceEps )
     &       iceFld(i,j) = ( iceFld(i,j)*tmpVol
     &                      -deltaTadvect*maskLocC(i,j)
     &                                   *( af(i+1,j)-af(i,j) )
     &                     )/iceVol(i,j)
            ENDDO
           ENDDO
          ENDIF
C-    keep advective flux (for diagnostics)
          DO j=sNy+1,sNy+OLy
           DO i=1-OLx+1,sNx+OLx
             afx(i,j) = af(i,j)
           ENDDO
          ENDDO
C-    end if North Edge
         ENDIF

        ELSE
C     do not only update the overlap
          jMinUpd = 1-OLy
          jMaxUpd = sNy+OLy
          IF ( interiorOnly .AND. S_edge ) jMinUpd = 1
          IF ( interiorOnly .AND. N_edge ) jMaxUpd = sNy
          IF ( useGridArea ) THEN
           DO j=jMinUpd,jMaxUpd
            DO i=1-OLx+1,sNx+OLx-1
             iceFld(i,j) = iceFld(i,j)
     &                    -deltaTadvect*maskLocC(i,j)
     &                                 *recip_rA(i,j,bi,bj)
     &                                 *( af(i+1,j)-af(i,j) )
            ENDDO
           ENDDO
          ELSE
           DO j=jMinUpd,jMaxUpd
            DO i=1-OLx+1,sNx+OLx-1
             tmpVol = iceVol(i,j)
             iceVol(i,j) = iceVol(i,j)
     &                    -deltaTadvect*maskLocC(i,j)
     &                           *( uTrans(i+1,j)-uTrans(i,j) )
             IF ( iceVol(i,j).GT.iceEps )
     &       iceFld(i,j) = ( iceFld(i,j)*tmpVol
     &                      -deltaTadvect*maskLocC(i,j)
     &                                   *( af(i+1,j)-af(i,j) )
     &                     )/iceVol(i,j)
            ENDDO
           ENDDO
          ENDIF
C-    keep advective flux (for diagnostics)
          DO j=jMinUpd,jMaxUpd
           DO i=1-OLx+1,sNx+OLx
             afx(i,j) = af(i,j)
           ENDDO
          ENDDO

C-     end if/else update overlap-Only
        ENDIF

C--   End of X direction
       ENDIF

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

       IF (calc_fluxes_Y) THEN

C-     Do not compute fluxes if
C       a) needed in overlap only
C   and b) the overlap of myTile are not cube-face edges
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE iceFld(:,:) = comlev1_thsice_4, key=ikey_4, byte=isbyte
CADJ STORE iceVol(:,:) = comlev1_thsice_4, key=ikey_4, byte=isbyte
CADJ STORE af(:,:)     = comlev1_thsice_4, key=ikey_4, byte=isbyte
#endif
        IF ( .NOT.overlapOnly .OR. E_edge .OR. W_edge ) THEN

C-     Advective flux in Y
         DO j=1-OLy,sNy+OLy
          DO i=1-OLx,sNx+OLx
            af(i,j) = 0.
          ENDDO
         ENDDO

cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
C-     Internal exchange for calculations in Y
         IF ( overlapOnly ) THEN
          CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
     &                               iceFld, bi,bj, myThid )
          IF ( .NOT.useGridArea )
     &    CALL FILL_CS_CORNER_TR_RL( 2, .FALSE.,
     &                               iceVol, bi,bj, myThid )
         ENDIF
cph-exch2#endif

C-     Compute CFL number
         IF ( useGridArea ) THEN
          DO j=2-Oly,sNy+Oly
           DO i=1-Olx,sNx+Olx
            vCfl(i,j) = deltaTadvect*(
     &          MAX( vTrans(i,j), 0. _d 0 )*recip_rA(i,j-1,bi,bj)
     &         +MAX(-vTrans(i,j), 0. _d 0 )*recip_rA(i, j ,bi,bj)
     &                               )
           ENDDO
          ENDDO
         ELSE
          DO j=2-Oly,sNy+Oly
           DO i=1-Olx,sNx+Olx
            vCfl(i,j) = deltaTadvect*(
     &        MAX( vTrans(i,j), 0. _d 0 )/MAX( iceVol(i,j-1), iceEps )
     &       +MAX(-vTrans(i,j), 0. _d 0 )/MAX( iceVol(i, j ), iceEps )
     &                               )
           ENDDO
          ENDDO
         ENDIF

         IF ( advectionScheme.EQ.ENUM_UPWIND_1RST
     &        .OR. advectionScheme.EQ.ENUM_DST2 ) THEN
          CALL GAD_DST2U1_ADV_Y( bi,bj,k, advectionScheme, .FALSE.,
     I             deltaTadvect, vTrans, vCfl, iceFld,
     O             af, myThid )
#ifdef ALLOW_DBUG_THSICE
          IF (dBugFlag.AND.bi.EQ.biDb ) WRITE(6,'(A,1P4E14.6)')
     &      'ICE_adv: yFx=', af(idb,jdb), iceFld(idb,jdb),
     &       vTrans(idb,jdb), af(idb,jdb)/vTrans(idb,jdb)
#endif /* ALLOW_DBUG_THSICE */
         ELSEIF (advectionScheme.EQ.ENUM_FLUX_LIMIT) THEN
          CALL GAD_FLUXLIMIT_ADV_Y( bi,bj,k, .FALSE., deltaTadvect,
     I             vTrans, vCfl, maskLocS, iceFld,
     O             af, myThid )
         ELSEIF (advectionScheme.EQ.ENUM_DST3 ) THEN
          CALL GAD_DST3_ADV_Y(      bi,bj,k, .FALSE., deltaTadvect,
     I             vTrans, vCfl, maskLocS, iceFld,
     O             af, myThid )
         ELSEIF (advectionScheme.EQ.ENUM_DST3_FLUX_LIMIT ) THEN
          CALL GAD_DST3FL_ADV_Y(    bi,bj,k, .FALSE., deltaTadvect,
     I             vTrans, vCfl, maskLocS, iceFld,
     O             af, myThid )
         ELSE
          STOP
     &  'THSICE_ADVECTION: adv. scheme incompatibale with mutli-dim'
         ENDIF

cph-exch2#ifndef ALLOW_AUTODIFF_TAMC
         IF ( overlapOnly .AND. ipass.EQ.1 ) THEN
          CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
     &                               iceFld, bi,bj, myThid )
          IF ( .NOT.useGridArea )
     &    CALL FILL_CS_CORNER_TR_RL( 1, .FALSE.,
     &                               iceVol, bi,bj, myThid )
         ENDIF
cph-exch2#endif

C-     Advective flux in Y : done
        ENDIF

C-     Update the local seaice field where needed:

C      update in overlap-Only
        IF ( overlapOnly ) THEN
         jMinUpd = 1-OLy+1
         jMaxUpd = sNy+OLy-1
C- notes: these 2 lines below have no real effect (because recip_hFac=0
C         in corner region) but safer to keep them.
         IF ( S_edge ) jMinUpd = 1
         IF ( N_edge ) jMaxUpd = sNy

#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE iceFld(:,:) = comlev1_thsice_4, key=ikey_4, byte=isbyte
CADJ STORE iceVol(:,:) = comlev1_thsice_4, key=ikey_4, byte=isbyte
#endif
         IF ( W_edge ) THEN
          IF ( useGridArea ) THEN
           DO j=jMinUpd,jMaxUpd
            DO i=1-OLx,0
             iceFld(i,j) = iceFld(i,j)
     &                    -deltaTadvect*maskLocC(i,j)
     &                                 *recip_rA(i,j,bi,bj)
     &                                 *( af(i,j+1)-af(i,j) )
            ENDDO
           ENDDO
          ELSE
           DO j=jMinUpd,jMaxUpd
            DO i=1-OLx,0
             tmpVol = iceVol(i,j)
             iceVol(i,j) = iceVol(i,j)
     &                    -deltaTadvect*maskLocC(i,j)
     &                           *( vTrans(i,j+1)-vTrans(i,j) )
             IF ( iceVol(i,j).GT.iceEps )
     &       iceFld(i,j) = ( iceFld(i,j)*tmpVol
     &                      -deltaTadvect*maskLocC(i,j)
     &                                   *( af(i,j+1)-af(i,j) )
     &                     )/iceVol(i,j)
            ENDDO
           ENDDO
          ENDIF
C-    keep advective flux (for diagnostics)
          DO j=1-OLy+1,sNy+OLy
           DO i=1-OLx,0
             afy(i,j) = af(i,j)
           ENDDO
          ENDDO
C-    end if West Edge
         ENDIF
         IF ( E_edge ) THEN
          IF ( useGridArea ) THEN
           DO j=jMinUpd,jMaxUpd
            DO i=sNx+1,sNx+OLx
             iceFld(i,j) = iceFld(i,j)
     &                    -deltaTadvect*maskLocC(i,j)
     &                                 *recip_rA(i,j,bi,bj)
     &                                 *( af(i,j+1)-af(i,j) )
            ENDDO
           ENDDO
          ELSE
           DO j=jMinUpd,jMaxUpd
            DO i=sNx+1,sNx+OLx
             tmpVol = iceVol(i,j)
             iceVol(i,j) = iceVol(i,j)
     &                    -deltaTadvect*maskLocC(i,j)
     &                           *( vTrans(i,j+1)-vTrans(i,j) )
             IF ( iceVol(i,j).GT.iceEps )
     &       iceFld(i,j) = ( iceFld(i,j)*tmpVol
     &                      -deltaTadvect*maskLocC(i,j)
     &                                   *( af(i,j+1)-af(i,j) )
     &                     )/iceVol(i,j)
            ENDDO
           ENDDO
          ENDIF
C-    keep advective flux (for diagnostics)
          DO j=1-OLy+1,sNy+OLy
           DO i=sNx+1,sNx+OLx
             afy(i,j) = af(i,j)
           ENDDO
          ENDDO
C-    end if East Edge
         ENDIF

        ELSE
C     do not only update the overlap
          iMinUpd = 1-OLx
          iMaxUpd = sNx+OLx
          IF ( interiorOnly .AND. W_edge ) iMinUpd = 1
          IF ( interiorOnly .AND. E_edge ) iMaxUpd = sNx
          IF ( useGridArea ) THEN
           DO j=1-OLy+1,sNy+OLy-1
            DO i=iMinUpd,iMaxUpd
             iceFld(i,j) = iceFld(i,j)
     &                    -deltaTadvect*maskLocC(i,j)
     &                                 *recip_rA(i,j,bi,bj)
     &                                 *( af(i,j+1)-af(i,j) )
            ENDDO
           ENDDO
          ELSE
           DO j=1-OLy+1,sNy+OLy-1
            DO i=iMinUpd,iMaxUpd
             tmpVol = iceVol(i,j)
             iceVol(i,j) = iceVol(i,j)
     &                    -deltaTadvect*maskLocC(i,j)
     &                           *( vTrans(i,j+1)-vTrans(i,j) )
             IF ( iceVol(i,j).GT.iceEps )
     &       iceFld(i,j) = ( iceFld(i,j)*tmpVol
     &                      -deltaTadvect*maskLocC(i,j)
     &                                   *( af(i,j+1)-af(i,j) )
     &                     )/iceVol(i,j)
            ENDDO
           ENDDO
          ENDIF
C-    keep advective flux (for diagnostics)
          DO j=1-OLy+1,sNy+OLy
           DO i=iMinUpd,iMaxUpd
             afy(i,j) = af(i,j)
           ENDDO
          ENDDO

C      end if/else update overlap-Only
        ENDIF

C--   End of Y direction
       ENDIF

C--   End of ipass loop
      ENDDO

C-    explicit advection is done ; add some debugging print
#ifdef ALLOW_DBUG_THSICE
      IF ( dBugFlag ) THEN
       DO j=1-OLy,sNy+OLy
        DO i=1-OLx,sNx+OLx
         IF ( i.EQ.idb .AND. j.EQ.jdb .AND. bi.EQ.biDb ) THEN
          tmpFac= deltaTadvect*recip_rA(i,j,bi,bj)
          WRITE(6,'(A,1P4E14.6)') 'ICE_adv:',
     &     afx(i,j)*tmpFac,afx(i+1,j)*tmpFac,
     &     afy(i,j)*tmpFac,afy(i,j+1)*tmpFac
         ENDIF
        ENDDO
       ENDDO
      ENDIF

#ifdef ALLOW_DEBUG
      IF ( debugLevel .GE. debLevC
     &     .AND. tracerIdentity.EQ.GAD_SI_HICE
     &     .AND. k.LE.3 .AND. myIter.EQ.1+nIter0
     &     .AND. nPx.EQ.1 .AND. nPy.EQ.1
     &     .AND. useCubedSphereExchange ) THEN
       CALL DEBUG_CS_CORNER_UV( ' afx,afy from THSICE_ADVECTION',
     &      afx,afy, k, standardMessageUnit,bi,bj,myThid )
      ENDIF
#endif /* ALLOW_DEBUG */
#endif /* ALLOW_DBUG_THSICE */

#endif /* ALLOW_GENERIC_ADVDIFF */

      RETURN
      END