pkg/thsice/thsice_advdiff.F

a4eca6e929 Jean*0001 #include "THSICE_OPTIONS.h"
e7f517b398 Jean*0002 #ifdef ALLOW_GENERIC_ADVDIFF
                0003 # include "GAD_OPTIONS.h"
                0004 #endif
6b47d550f4 Mart*0005 #ifdef ALLOW_AUTODIFF
                0006 # include "AUTODIFF_OPTIONS.h"
                0007 #endif
a4eca6e929 Jean*0008
                0009 CBOP
                0010 C !ROUTINE: THSICE_ADVDIFF
                0011
                0012 C !INTERFACE: ==========================================================
                0013       SUBROUTINE THSICE_ADVDIFF(
                0014      U                  uIce, vIce,
                0015      I                  bi, bj, myTime, myIter, myThid )
                0016
                0017 C !DESCRIPTION: \bv
                0018 C     *===========================================================*
                0019 C     | SUBROUTINE THSICE_ADVDIFF
                0020 C     | o driver for different advection routines
                0021 C     |   calls an adaption of gad_advection to call different
                0022 C     |   advection routines of pkg/generic_advdiff
                0023 C     *===========================================================*
                0024 C \ev
                0025
                0026 C !USES: ===============================================================
                0027       IMPLICIT NONE
                0028
                0029 C     === Global variables ===
                0030 C   oceFWfx   :: fresh water flux to the ocean  [kg/m^2/s]
ba0b047096 Mart*0031 C   oceSflx   :: salt flux to the ocean         [g/m^2/s]
a4eca6e929 Jean*0032 C   oceQnet   :: heat flux to the ocean         [W/m^2]
                0033
                0034 #include "SIZE.h"
                0035 #include "EEPARAMS.h"
                0036 #include "PARAMS.h"
                0037 #include "GRID.h"
                0038 #include "THSICE_SIZE.h"
                0039 #include "THSICE_PARAMS.h"
                0040 #include "THSICE_VARS.h"
e7f517b398 Jean*0041 #ifdef ALLOW_GENERIC_ADVDIFF
                0042 # include "GAD.h"
                0043 #endif
d6f06800ae Patr*0044 #ifdef ALLOW_AUTODIFF_TAMC
                0045 # include "tamc.h"
                0046 #endif
a4eca6e929 Jean*0047
                0048 C !INPUT PARAMETERS: ===================================================
                0049 C     === Routine arguments ===
                0050 C     uIce/vIce :: ice velocity on C-grid [m/s]
                0051 C     bi,bj     :: Tile indices
                0052 C     myTime    :: Current time in simulation (s)
                0053 C     myIter    :: Current iteration number
                0054 C     myThid    :: My Thread Id. number
                0055       _RL     uIce  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0056       _RL     vIce  (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0057       INTEGER bi,bj
                0058       _RL     myTime
                0059       INTEGER myIter
                0060       INTEGER myThid
                0061 CEndOfInterface
                0062
                0063 #ifdef ALLOW_THSICE
                0064 C !LOCAL VARIABLES: ====================================================
                0065 C     === Local variables ===
                0066 C     i,j,      :: Loop counters
                0067 C     uTrans    :: sea-ice area transport, x direction
                0068 C     vTrans    :: sea-ice area transport, y direction
                0069 C     uTrIce    :: sea-ice volume transport, x direction
                0070 C     vTrIce    :: sea-ice volume transport, y direction
                0071 C     afx       :: horizontal advective flux, x direction
                0072 C     afy       :: horizontal advective flux, y direction
                0073 C     iceFrc    :: (new) sea-ice fraction
                0074 C     iceVol    :: temporary array used in advection S/R
                0075 C     oldVol    :: (old) sea-ice volume
7c0d02fed1 Jean*0076 C     msgBuf    :: Informational/error message buffer
a4eca6e929 Jean*0077       INTEGER i, j
                0078       LOGICAL thSIce_multiDimAdv
                0079       CHARACTER*(MAX_LEN_MBUF) msgBuf
                0080
                0081       _RL uTrans    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0082       _RL vTrans    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0083       _RL uTrIce    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0084       _RL vTrIce    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0085       _RL afx       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0086       _RL afy       (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0087       _RS maskOce   (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0088       _RL iceFrc    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
7c0d02fed1 Jean*0089       _RL iceVol    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0090       _RL oldVol    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
68ae87c30f Patr*0091       _RL iceTmp    (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
8db81c2271 Jean*0092       _RL r_minArea
a4eca6e929 Jean*0093       _RL meanCellArea, areaEpsil, vol_Epsil
                0094 #ifdef ALLOW_DIAGNOSTICS
                0095       CHARACTER*8 diagName
                0096       CHARACTER*4 THSICE_DIAG_SUFX, diagSufx
                0097       EXTERNAL    THSICE_DIAG_SUFX
0e322d6e56 Jean*0098       LOGICAL  DIAGNOSTICS_IS_ON
                0099       EXTERNAL DIAGNOSTICS_IS_ON
                0100       _RL tmpFac
a4eca6e929 Jean*0101 #endif
7c50f07931 Mart*0102 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0103 C     tkey :: tape key (depends on tiles)
                0104       INTEGER tkey
7c50f07931 Mart*0105 #endif
03026c2540 Jean*0106 #ifdef ALLOW_DBUG_THSICE
                0107       _RL tmpVar, sumVar1, sumVar2
                0108 #endif
                0109       LOGICAL dBugFlag
                0110 #include "THSICE_DEBUG.h"
a4eca6e929 Jean*0111 CEOP
                0112
                0113 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
d6f06800ae Patr*0114
                0115 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0116       tkey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
d6f06800ae Patr*0117 #endif /* ALLOW_AUTODIFF_TAMC */
                0118
a4eca6e929 Jean*0119 C     areaEpsil, vol_Epsil are 2 small numbers for ice area & ice volume:
                0120 C     if ice area (=ice fraction * grid-cell area) or ice volume (= effective
                0121 C     thickness * grid-cell area) are too small (i.e.: < areaEpsil,vol_Epsil)
                0122 C     will assume that ice is gone, and will loose mass or energy.
                0123 C     However, if areaEpsil,vol_Epsil are much smaller than minimun ice area
c2c29a499f Jean*0124 C     (iceMaskMin*rAc) and minimum ice volume (iceMaskMin*hIceMin*rAc),
a4eca6e929 Jean*0125 C     good chance that this will never happen within 1 time step.
                0126
ae605e558b Jean*0127       dBugFlag = debugLevel.GE.debLevC
a4eca6e929 Jean*0128
                0129 C-    definitively not an accurate computation of mean grid-cell area;
                0130 C     but what matter here is just to have the right order of magnitude.
                0131       meanCellArea = Nx*Ny
                0132       meanCellArea = globalArea / meanCellArea
                0133       areaEpsil = 1. _d -10 * meanCellArea
                0134       vol_Epsil = 1. _d -15 * meanCellArea
                0135
                0136       r_minArea = 0. _d 0
8db81c2271 Jean*0137       IF ( iceMaskMin.GT.0. _d 0 ) r_minArea = 1. _d 0 / iceMaskMin
a4eca6e929 Jean*0138
                0139       thSIce_multiDimAdv = .TRUE.
e7f517b398 Jean*0140 #ifdef ALLOW_GENERIC_ADVDIFF
a4eca6e929 Jean*0141       IF ( thSIceAdvScheme.EQ.ENUM_CENTERED_2ND
                0142      & .OR.thSIceAdvScheme.EQ.ENUM_UPWIND_3RD
                0143      & .OR.thSIceAdvScheme.EQ.ENUM_CENTERED_4TH ) THEN
                0144        thSIce_multiDimAdv = .FALSE.
                0145       ENDIF
e7f517b398 Jean*0146 #endif /* ALLOW_GENERIC_ADVDIFF */
a4eca6e929 Jean*0147
7c0d02fed1 Jean*0148 #ifdef ALLOW_DIAGNOSTICS
                0149       IF ( useDiagnostics ) THEN
                0150         CALL DIAGNOSTICS_FILL(iceMask,'SI_AdvFr',0,1,1,bi,bj,myThid)
                0151 C-     Ice-fraction weighted quantities:
                0152         tmpFac = 1. _d 0
                0153         CALL DIAGNOSTICS_FRACT_FILL(
                0154      I                   iceHeight, iceMask,tmpFac,1,'SI_AdvHi',
                0155      I                   0,1,1,bi,bj,myThid)
                0156         CALL DIAGNOSTICS_FRACT_FILL(
                0157      I                   snowHeight,iceMask,tmpFac,1,'SI_AdvHs',
                0158      I                   0,1,1,bi,bj,myThid)
                0159 C-     Ice-Volume weighted quantities:
                0160         IF ( DIAGNOSTICS_IS_ON('SI_AdvQ1',myThid) .OR.
                0161      &       DIAGNOSTICS_IS_ON('SI_AdvQ2',myThid) ) THEN
                0162          DO j=1,sNy
                0163           DO i=1,sNx
68ae87c30f Patr*0164            iceTmp(i,j) = iceMask(i,j,bi,bj)*iceHeight(i,j,bi,bj)
7c0d02fed1 Jean*0165           ENDDO
                0166          ENDDO
                0167          CALL DIAGNOSTICS_FRACT_FILL(
                0168      I                   Qice1(1-OLx,1-OLy,bi,bj),
68ae87c30f Patr*0169      I                   iceTmp,tmpFac,1,'SI_AdvQ1',
7c0d02fed1 Jean*0170      I                   0,1,2,bi,bj,myThid)
                0171          CALL DIAGNOSTICS_FRACT_FILL(
                0172      I                   Qice2(1-OLx,1-OLy,bi,bj),
68ae87c30f Patr*0173      I                   iceTmp,tmpFac,1,'SI_AdvQ2',
7c0d02fed1 Jean*0174      I                   0,1,2,bi,bj,myThid)
                0175         ENDIF
                0176       ENDIF
                0177 #endif /* ALLOW_DIAGNOSTICS */
                0178
a4eca6e929 Jean*0179 C--   Initialisation (+ build oceanic mask)
                0180       DO j=1-OLy,sNy+OLy
                0181        DO i=1-OLx,sNx+OLx
                0182          maskOce(i,j) = 0. _d 0
                0183          IF ( hOceMxL(i,j,bi,bj).GT.0. ) maskOce(i,j) = 1.
                0184          iceVol(i,j) = 0. _d 0
                0185          uTrans(i,j) = 0. _d 0
                0186          vTrans(i,j) = 0. _d 0
                0187          uTrIce(i,j) = 0. _d 0
                0188          vTrIce(i,j) = 0. _d 0
                0189          oceFWfx(i,j,bi,bj) = 0. _d 0
                0190          oceSflx(i,j,bi,bj) = 0. _d 0
                0191          oceQnet(i,j,bi,bj) = 0. _d 0
                0192        ENDDO
                0193       ENDDO
                0194
d6f06800ae Patr*0195 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0196 CADJ STORE iceHeight (:,:,bi,bj) = comlev1_bibj,key=tkey,byte=isbyte
                0197 CADJ STORE snowHeight(:,:,bi,bj) = comlev1_bibj,key=tkey,byte=isbyte
d6f06800ae Patr*0198 #endif
a4eca6e929 Jean*0199       IF ( thSIce_diffK .GT. 0. ) THEN
                0200         CALL THSICE_DIFFUSION(
                0201      I              maskOce,
                0202      U              uIce, vIce,
                0203      I              bi, bj, myTime, myIter, myThid )
                0204       ENDIF
                0205
                0206       IF ( thSIce_multiDimAdv ) THEN
                0207
                0208 C-    Calculate ice transports through tracer cell faces.
7c0d02fed1 Jean*0209         DO j=1-OLy,sNy+OLy
                0210          DO i=1-OLx+1,sNx+OLx
a4eca6e929 Jean*0211           uTrIce(i,j) = uIce(i,j)*_dyG(i,j,bi,bj)
                0212      &                *maskOce(i-1,j)*maskOce(i,j)
                0213          ENDDO
                0214         ENDDO
7c0d02fed1 Jean*0215         DO j=1-OLy+1,sNy+OLy
                0216          DO i=1-OLx,sNx+OLx
a4eca6e929 Jean*0217           vTrIce(i,j) = vIce(i,j)*_dxG(i,j,bi,bj)
                0218      &                *maskOce(i,j-1)*maskOce(i,j)
                0219          ENDDO
                0220         ENDDO
                0221
                0222 C--   Fractional area
7c0d02fed1 Jean*0223         DO j=1-OLy,sNy+OLy
                0224          DO i=1-OLx,sNx+OLx
a4eca6e929 Jean*0225           iceFrc(i,j) = iceMask(i,j,bi,bj)
                0226          ENDDO
                0227         ENDDO
d6f06800ae Patr*0228 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0229 CADJ STORE utrice = comlev1_bibj, key=tkey, byte=isbyte
                0230 CADJ STORE vtrice = comlev1_bibj, key=tkey, byte=isbyte
d6f06800ae Patr*0231 #endif
a4eca6e929 Jean*0232         CALL THSICE_ADVECTION(
                0233      I       GAD_SI_FRAC,  thSIceAdvScheme, .TRUE.,
                0234      I       uTrIce, vTrIce, maskOce, thSIce_deltaT, areaEpsil,
                0235      U       iceVol, iceFrc,
                0236      O       uTrans, vTrans,
                0237      I       bi, bj, myTime, myIter, myThid )
                0238
                0239 C--   Snow thickness
7c0d02fed1 Jean*0240         DO j=1-OLy,sNy+OLy
                0241          DO i=1-OLx,sNx+OLx
a4eca6e929 Jean*0242           iceVol(i,j) = iceMask(i,j,bi,bj)*rA(i,j,bi,bj)
                0243          ENDDO
                0244         ENDDO
a85293d087 Mart*0245 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0246 CADJ STORE uTrans                = comlev1_bibj,key=tkey,byte=isbyte
                0247 CADJ STORE vTrans                = comlev1_bibj,key=tkey,byte=isbyte
a85293d087 Mart*0248 #endif
a4eca6e929 Jean*0249         CALL THSICE_ADVECTION(
                0250      I       GAD_SI_HSNOW, thSIceAdvScheme, .FALSE.,
                0251      I       uTrans, vTrans, maskOce, thSIce_deltaT, areaEpsil,
7c0d02fed1 Jean*0252      U       iceVol, snowHeight(1-OLx,1-OLy,bi,bj),
a4eca6e929 Jean*0253      O       afx, afy,
                0254      I       bi, bj, myTime, myIter, myThid )
                0255
d6f06800ae Patr*0256 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0257 CADJ STORE iceHeight (:,:,bi,bj) = comlev1_bibj,key=tkey,byte=isbyte
                0258 CADJ STORE iceMask   (:,:,bi,bj) = comlev1_bibj,key=tkey,byte=isbyte
d6f06800ae Patr*0259 #endif
a4eca6e929 Jean*0260 C--   sea-ice Thickness
7c0d02fed1 Jean*0261         DO j=1-OLy,sNy+OLy
                0262          DO i=1-OLx,sNx+OLx
a4eca6e929 Jean*0263           iceVol(i,j) = iceMask(i,j,bi,bj)*rA(i,j,bi,bj)
                0264           oldVol(i,j) = iceVol(i,j)*iceHeight(i,j,bi,bj)
                0265          ENDDO
                0266         ENDDO
                0267         CALL THSICE_ADVECTION(
                0268      I       GAD_SI_HICE,  thSIceAdvScheme, .FALSE.,
                0269      I       uTrans, vTrans, maskOce, thSIce_deltaT, areaEpsil,
7c0d02fed1 Jean*0270      U       iceVol, iceHeight(1-OLx,1-OLy,bi,bj),
a4eca6e929 Jean*0271      O       uTrIce, vTrIce,
                0272      I       bi, bj, myTime, myIter, myThid )
d6f06800ae Patr*0273
                0274 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0275 CADJ STORE utrice      = comlev1_bibj, key=tkey, byte=isbyte
                0276 CADJ STORE vtrice      = comlev1_bibj, key=tkey, byte=isbyte
d6f06800ae Patr*0277 #endif
                0278
a4eca6e929 Jean*0279 #ifdef ALLOW_DBUG_THSICE
                0280         IF ( dBugFlag ) THEN
                0281          sumVar1 = 0.
                0282          sumVar2 = 0.
                0283          DO j=1,sNy
                0284           DO i=1,sNx
                0285 C-      Check that updated iceVol = iceFrc*rA
                0286            tmpVar = ABS(iceVol(i,j)-iceFrc(i,j)*rA(i,j,bi,bj))
                0287            IF ( tmpVar.GT.0. ) THEN
                0288              sumVar1 = sumVar1 + 1.
                0289              sumVar2 = sumVar2 + tmpVar
                0290            ENDIF
                0291            IF ( tmpVar.GT.vol_Epsil ) THEN
                0292             WRITE(6,'(A,2I4,2I2,I12)') 'ARE_ADV: ij,bij,it=',
                0293      &                                  i,j,bi,bj,myIter
                0294             WRITE(6,'(2(A,1P2E14.6))') 'ARE_ADV: iceVol,iceFrc*rA=',
                0295      &                        iceVol(i,j),iceFrc(i,j)*rA(i,j,bi,bj),
                0296      &            ' , diff=', tmpVar
                0297            ENDIF
                0298            IF ( dBug(i,j,bi,bj) ) THEN
                0299             WRITE(6,'(A,2I4,2I2,I12)') 'ICE_ADV: ij,bij,it=',
                0300      &                                  i,j,bi,bj,myIter
                0301             WRITE(6,'(2(A,1P2E14.6))')
                0302      &       'ICE_ADV: uIce=', uIce(i,j), uIce(i+1,j),
                0303      &             ' , vIce=', vIce(i,j), vIce(i,j+1)
                0304             WRITE(6,'(2(A,1P2E14.6))')
0e322d6e56 Jean*0305      &       'ICE_ADV: area_b,a=', iceMask(i,j,bi,bj), iceFrc(i,j),
                0306      &       ' , Heff_b,a=', oldVol(i,j)*recip_rA(i,j,bi,bj),
                0307      &                       iceHeight(i,j,bi,bj)*iceFrc(i,j)
a4eca6e929 Jean*0308            ENDIF
                0309           ENDDO
                0310          ENDDO
                0311          IF ( sumVar2.GT.vol_Epsil )
                0312      &   WRITE(6,'(A,2I2,I10,A,I4,1P2E14.6)') 'ARE_ADV: bij,it:',
                0313      &                    bi,bj,myIter, ' ; Npts,aveDiff,Epsil=',
                0314      &                    INT(sumVar1),sumVar2/sumVar1,vol_Epsil
                0315         ENDIF
                0316 #endif
                0317 #ifdef ALLOW_DIAGNOSTICS
                0318 C--     Diagnosse advective fluxes (ice-fraction, snow & ice thickness):
                0319         IF ( useDiagnostics ) THEN
                0320           diagSufx = THSICE_DIAG_SUFX( GAD_SI_FRAC, myThid )
                0321           diagName = 'ADVx'//diagSufx
                0322           CALL DIAGNOSTICS_FILL( uTrans, diagName, 1,1,2,bi,bj, myThid )
                0323           diagName = 'ADVy'//diagSufx
                0324           CALL DIAGNOSTICS_FILL( vTrans, diagName, 1,1,2,bi,bj, myThid )
                0325
                0326           diagSufx = THSICE_DIAG_SUFX( GAD_SI_HSNOW, myThid )
                0327           diagName = 'ADVx'//diagSufx
                0328           CALL DIAGNOSTICS_FILL( afx, diagName, 1,1,2,bi,bj, myThid )
                0329           diagName = 'ADVy'//diagSufx
                0330           CALL DIAGNOSTICS_FILL( afy, diagName, 1,1,2,bi,bj, myThid )
                0331
                0332           diagSufx = THSICE_DIAG_SUFX( GAD_SI_HICE, myThid )
                0333           diagName = 'ADVx'//diagSufx
                0334           CALL DIAGNOSTICS_FILL( uTrIce, diagName, 1,1,2,bi,bj, myThid )
                0335           diagName = 'ADVy'//diagSufx
                0336           CALL DIAGNOSTICS_FILL( vTrIce, diagName, 1,1,2,bi,bj, myThid )
                0337         ENDIF
                0338 #endif
                0339
                0340 C--   Enthalpy in layer 1
7c0d02fed1 Jean*0341         DO j=1-OLy,sNy+OLy
                0342          DO i=1-OLx,sNx+OLx
a4eca6e929 Jean*0343           iceVol(i,j) = oldVol(i,j)
                0344          ENDDO
                0345         ENDDO
                0346         CALL THSICE_ADVECTION(
                0347      I       GAD_SI_QICE1, thSIceAdvScheme, .FALSE.,
                0348      I       uTrIce, vTrIce, maskOce, thSIce_deltaT, vol_Epsil,
7c0d02fed1 Jean*0349      U       iceVol, Qice1(1-OLx,1-OLy,bi,bj),
a4eca6e929 Jean*0350      O       afx, afy,
                0351      I       bi, bj, myTime, myIter, myThid )
                0352 #ifdef ALLOW_DBUG_THSICE
                0353         IF ( dBugFlag ) THEN
                0354          DO j=1,sNy
                0355           DO i=1,sNx
                0356            IF ( dBug(i,j,bi,bj) ) THEN
                0357 c           WRITE(6,'(A,1P4E14.6)') 'ICE_ADV: Qice1_b,a=',
                0358 c    &             Qice1(i,j,bi,bj),
                0359 c    &        ( iceFld(i,j) + thSIce_deltaT * gFld(i,j)
                0360 c    &          ) * recip_heff(i,j)
                0361 c           WRITE(6,'(A,1P4E14.6)') 'ICE_ADV: q1Fx=', gFld(i,j)
                0362            ENDIF
                0363           ENDDO
                0364          ENDDO
                0365         ENDIF
                0366 #endif
                0367 #ifdef ALLOW_DIAGNOSTICS
                0368         IF ( useDiagnostics ) THEN
                0369           diagSufx = THSICE_DIAG_SUFX( GAD_SI_QICE1, myThid )
                0370           diagName = 'ADVx'//diagSufx
                0371           CALL DIAGNOSTICS_FILL( afx, diagName, 1,1,2,bi,bj, myThid )
                0372           diagName = 'ADVy'//diagSufx
                0373           CALL DIAGNOSTICS_FILL( afy, diagName, 1,1,2,bi,bj, myThid )
                0374         ENDIF
                0375 #endif
                0376
                0377 C--   Enthalpy in layer 2
7c0d02fed1 Jean*0378         DO j=1-OLy,sNy+OLy
                0379          DO i=1-OLx,sNx+OLx
a4eca6e929 Jean*0380           iceVol(i,j) = oldVol(i,j)
                0381          ENDDO
                0382         ENDDO
                0383         CALL THSICE_ADVECTION(
                0384      I       GAD_SI_QICE2, thSIceAdvScheme, .FALSE.,
                0385      I       uTrIce, vTrIce, maskOce, thSIce_deltaT, vol_Epsil,
7c0d02fed1 Jean*0386      U       iceVol, Qice2(1-OLx,1-OLy,bi,bj),
a4eca6e929 Jean*0387      O       afx, afy,
                0388      I       bi, bj, myTime, myIter, myThid )
                0389 #ifdef ALLOW_DBUG_THSICE
                0390         IF ( dBugFlag ) THEN
                0391          sumVar1 = 0.
                0392          sumVar2 = 0.
                0393          DO j=1,sNy
                0394           DO i=1,sNx
                0395 C-      Check that updated iceVol = Hic*Frc*rA
                0396            tmpVar = ABS(iceVol(i,j)
                0397      &             -iceHeight(i,j,bi,bj)*iceFrc(i,j)*rA(i,j,bi,bj))
                0398            IF ( tmpVar.GT.0. ) THEN
                0399              sumVar1 = sumVar1 + 1.
                0400              sumVar2 = sumVar2 + tmpVar
                0401            ENDIF
                0402            IF ( tmpVar.GT.vol_Epsil ) THEN
                0403             WRITE(6,'(A,2I4,2I2,I12)') 'VOL_ADV: ij,bij,it=',
                0404      &                                  i,j,bi,bj,myIter
                0405             WRITE(6,'(2(A,1P2E14.6))') 'VOL_ADV: iceVol,Hic*Frc*rA=',
                0406      &      iceVol(i,j),iceHeight(i,j,bi,bj)*iceFrc(i,j)*rA(i,j,bi,bj),
                0407      &             ' , diff=', tmpVar
                0408            ENDIF
                0409            IF ( dBug(i,j,bi,bj) ) THEN
                0410 c           WRITE(6,'(A,1P4E14.6)') 'ICE_ADV: Qice2_b,a=',
                0411 c    &             Qice2(i,j,bi,bj),
                0412 c    &        ( iceFld(i,j) + thSIce_deltaT * gFld(i,j)
                0413 c    &          ) * recip_heff(i,j)
                0414 c           WRITE(6,'(A,1P4E14.6)') 'ICE_ADV: q2Fx=', gFld(i,j)
                0415            ENDIF
                0416           ENDDO
                0417          ENDDO
                0418          IF ( sumVar2.GT.vol_Epsil )
                0419      &   WRITE(6,'(A,2I2,I10,A,I4,1P2E14.6)') 'VOL_ADV: bij,it:',
                0420      &                    bi,bj,myIter, ' ; Npts,aveDiff,Epsil=',
                0421      &                    INT(sumVar1),sumVar2/sumVar1,vol_Epsil
                0422         ENDIF
                0423 #endif
                0424 #ifdef ALLOW_DIAGNOSTICS
                0425         IF ( useDiagnostics ) THEN
                0426           diagSufx = THSICE_DIAG_SUFX( GAD_SI_QICE2, myThid )
                0427           diagName = 'ADVx'//diagSufx
                0428           CALL DIAGNOSTICS_FILL( afx, diagName, 1,1,2,bi,bj, myThid )
                0429           diagName = 'ADVy'//diagSufx
                0430           CALL DIAGNOSTICS_FILL( afy, diagName, 1,1,2,bi,bj, myThid )
                0431         ENDIF
                0432 #endif
                0433
d6f06800ae Patr*0434 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0435 CADJ STORE iceHeight (:,:,bi,bj) = comlev1_bibj,key=tkey,byte=isbyte
                0436 CADJ STORE snowHeight(:,:,bi,bj) = comlev1_bibj,key=tkey,byte=isbyte
                0437 CADJ STORE iceFrc                = comlev1_bibj,key=tkey,byte=isbyte
d6f06800ae Patr*0438 #endif
                0439
8db81c2271 Jean*0440 C--   Update Ice Fraction: ensure that fraction is > iceMaskMin & < 1
                0441 C      and adjust Ice thickness and snow thickness accordingly
a4eca6e929 Jean*0442         DO j=1,sNy
                0443          DO i=1,sNx
8db81c2271 Jean*0444           IF ( iceFrc(i,j) .GT. 1. _d 0 ) THEN
0e322d6e56 Jean*0445 c           IF ( dBug(i,j,bi,bj) )
a4eca6e929 Jean*0446             iceMask(i,j,bi,bj)    = 1. _d 0
8db81c2271 Jean*0447             iceHeight(i,j,bi,bj)  = iceHeight(i,j,bi,bj) *iceFrc(i,j)
a4eca6e929 Jean*0448             snowHeight(i,j,bi,bj) = snowHeight(i,j,bi,bj)*iceFrc(i,j)
8db81c2271 Jean*0449           ELSEIF ( iceFrc(i,j) .LT. iceMaskMin ) THEN
0e322d6e56 Jean*0450 c           IF ( dBug(i,j,bi,bj) )
8db81c2271 Jean*0451             iceMask(i,j,bi,bj)    = iceMaskMin
                0452             iceHeight(i,j,bi,bj)  = iceHeight(i,j,bi,bj)
                0453      &                             *iceFrc(i,j)*r_minArea
a4eca6e929 Jean*0454             snowHeight(i,j,bi,bj) = snowHeight(i,j,bi,bj)
                0455      &                             *iceFrc(i,j)*r_minArea
                0456           ELSE
8db81c2271 Jean*0457             iceMask(i,j,bi,bj)    = iceFrc(i,j)
                0458           ENDIF
                0459          ENDDO
                0460         ENDDO
a85293d087 Mart*0461 #ifdef ALLOW_AUTODIFF_TAMC
                0462 C     Store directives to avoid recomputing of this routine in ad-code
edb6656069 Mart*0463 CADJ STORE iceHeight (:,:,bi,bj) = comlev1_bibj,key=tkey,byte=isbyte
                0464 CADJ STORE snowHeight(:,:,bi,bj) = comlev1_bibj,key=tkey,byte=isbyte
                0465 CADJ STORE iceMask   (:,:,bi,bj) = comlev1_bibj,key=tkey,byte=isbyte
                0466 CADJ STORE qIce1     (:,:,bi,bj) = comlev1_bibj,key=tkey,byte=isbyte
                0467 CADJ STORE qIce2     (:,:,bi,bj) = comlev1_bibj,key=tkey,byte=isbyte
a85293d087 Mart*0468 #endif
d680aaaecd Jean*0469 C-     adjust sea-ice state if ice is too thin.
8db81c2271 Jean*0470         DO j=1,sNy
                0471          DO i=1,sNx
c2c29a499f Jean*0472           IF ( iceHeight(i,j,bi,bj).LT.hIceMin ) THEN
d680aaaecd Jean*0473            iceVol(i,j) = iceMask(i,j,bi,bj)*iceHeight(i,j,bi,bj)
                0474 c          IF ( dBug(i,j,bi,bj) )
                0475            IF ( iceVol(i,j).GE.hIceMin*iceMaskMin ) THEN
                0476             iceMask(i,j,bi,bj)    = iceVol(i,j)/hIceMin
                0477             snowHeight(i,j,bi,bj) = snowHeight(i,j,bi,bj)
                0478      &                             *hIceMin/iceHeight(i,j,bi,bj)
                0479             iceHeight(i,j,bi,bj)  = hIceMin
                0480            ELSE
03026c2540 Jean*0481 C-    Not enough ice, melt the tiny amount of snow & ice:
8db81c2271 Jean*0482 C     and return fresh-water, salt & energy to the ocean (flx > 0 = into ocean)
03026c2540 Jean*0483 C- -  Note: using 1rst.Order Upwind, I can get the same results as when
                0484 C     using seaice_advdiff (with SEAICEadvScheme=1) providing I comment
                0485 C     out the following lines (and then loose conservation).
                0486 C- -
8db81c2271 Jean*0487             oceFWfx(i,j,bi,bj) =  ( rhos*snowHeight(i,j,bi,bj)
                0488      &                             +rhoi*iceHeight(i,j,bi,bj)
                0489      &                            )*iceMask(i,j,bi,bj)/thSIce_deltaT
c2c29a499f Jean*0490             oceSflx(i,j,bi,bj) =    rhoi*iceHeight(i,j,bi,bj)*saltIce
8db81c2271 Jean*0491      &                             *iceMask(i,j,bi,bj)/thSIce_deltaT
                0492             oceQnet(i,j,bi,bj) = -( rhos*snowHeight(i,j,bi,bj)*qsnow
                0493      &                             +rhoi*iceHeight(i,j,bi,bj)
                0494      &                                  *( Qice1(i,j,bi,bj)
                0495      &                                    +Qice2(i,j,bi,bj) )*0.5 _d 0
                0496      &                            )*iceMask(i,j,bi,bj)/thSIce_deltaT
03026c2540 Jean*0497 C- -
a4eca6e929 Jean*0498 c           flx2oc (i,j) = flx2oc (i,j) +
                0499 c           frw2oc (i,j) = frw2oc (i,j) +
                0500 c           fsalt  (i,j) = fsalt  (i,j) +
                0501             iceMask   (i,j,bi,bj) = 0. _d 0
                0502             iceHeight (i,j,bi,bj) = 0. _d 0
                0503             snowHeight(i,j,bi,bj) = 0. _d 0
                0504             Qice1     (i,j,bi,bj) = 0. _d 0
                0505             Qice2     (i,j,bi,bj) = 0. _d 0
                0506             snowAge   (i,j,bi,bj) = 0. _d 0
d680aaaecd Jean*0507            ENDIF
a4eca6e929 Jean*0508           ENDIF
                0509          ENDDO
                0510         ENDDO
                0511
                0512 #ifdef ALLOW_DBUG_THSICE
                0513         IF ( dBugFlag ) THEN
                0514          DO j=1,sNy
                0515           DO i=1,sNx
                0516            IF ( dBug(i,j,bi,bj) ) THEN
                0517             WRITE(6,'(2(A,1P2E14.6))')
                0518 c    &       'ICE_ADV: area_b,a=', AREA(i,j,2,bi,bj),AREA(i,j,1,bi,bj)
                0519 c           WRITE(6,'(A,1P4E14.6)') 'ICE_ADV: mFx=', gFld(i,j)
                0520            ENDIF
                0521           ENDDO
                0522          ENDDO
                0523         ENDIF
                0524 #endif
                0525
                0526       ELSE
                0527 C---  if not multiDimAdvection
                0528
                0529         WRITE(msgBuf,'(2A)') 'S/R THSICE_ADVDIFF: ',
                0530      &       'traditional advection/diffusion not yet implemented'
                0531         CALL PRINT_ERROR( msgBuf , myThid)
                0532         WRITE(msgBuf,'(2A)') '                    ',
                0533      &       'for ThSice variable Qice1, Qice2, SnowHeight. Sorry!'
                0534         CALL PRINT_ERROR( msgBuf , myThid)
                0535           STOP 'ABNORMAL: END: S/R THSICE_ADVDIFF'
                0536
                0537 C---  end if multiDimAdvection
                0538       ENDIF
                0539
                0540 #endif /* ALLOW_THSICE */
                0541
                0542       RETURN
                0543       END