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