/[MITgcm]/MITgcm_contrib/torge/itd/code/seaice_growth.F

Diff of /MITgcm_contrib/torge/itd/code/seaice_growth.F

Parent Directory | Revision Log | View Revision Graph Revision Graph | View Patch Patch

-revision 1.2 by dimitri,
Fri Apr 27 22:25:23 2012 UTC
+revision 1.14 by torge,
Mon Dec 10 22:19:49 2012 UTC
 Line 91 
 C     i,j,bi,bj :: Loop counters
        INTEGER i, j, bi, bj
  C     number of surface interface layer
        INTEGER kSurface
+ C     IT :: ice thickness category index (MULTICATEGORIES and ITD code)
+       INTEGER IT
+ C     msgBuf      :: Informational/error message buffer
+ #ifdef ALLOW_BALANCE_FLUXES
+       CHARACTER*(MAX_LEN_MBUF) msgBuf
+ #elif (defined (SEAICE_DEBUG))
+       CHARACTER*(MAX_LEN_MBUF) msgBuf
+       CHARACTER*12 msgBufForm
+ #endif
  C     constants
+       _RL pFac
        _RL tempFrz, ICE2SNOW, SNOW2ICE
        _RL QI, QS, recip_QI
+       _RL lhSublim
+ C conversion factors to go from Q (W/m2) to HEFF (ice meters)
+       _RL convertQ2HI, convertHI2Q
+ C conversion factors to go from precip (m/s) unit to HEFF (ice meters)
+       _RL convertPRECIP2HI, convertHI2PRECIP
+ C     Factor by which we increase the upper ocean friction velocity (u*) when
+ C     ice is absent in a grid cell  (dimensionless)
+       _RL MixedLayerTurbulenceFactor
+ C     wind speed square
+       _RL SPEED_SQ
+ C     Regularization values squared
+       _RL area_reg_sq, hice_reg_sq
+ C     pathological cases thresholds
+       _RL heffTooHeavy
+ C     Helper variables: reciprocal of some constants
+       _RL recip_multDim
+       _RL recip_deltaTtherm
+       _RL recip_rhoIce
+ C     local value (=1/HO or 1/HO_south)
+       _RL recip_HO
+ C     local value (=1/ice thickness)
+       _RL recip_HH
+ C     facilitate multi-category snow implementation
+       _RL pFacSnow
+ C     temporary variables available for the various computations
+       _RL tmpscal0, tmpscal1, tmpscal2, tmpscal3, tmpscal4
+ #ifdef SEAICE_ITD
+       _RL tmpscal1itd(1:sNx,1:sNy), tmpscal2itd(1:sNx,1:sNy)
+       _RL tmpscal3itd(1:sNx,1:sNy)
+ #endif
+ #ifdef ALLOW_SITRACER
+       INTEGER iTr
+ #ifdef ALLOW_DIAGNOSTICS
+       CHARACTER*8   diagName
+ #endif
+ #endif /* ALLOW_SITRACER */
+ #ifdef ALLOW_AUTODIFF_TAMC
+       INTEGER ilockey
+ #endif
+ C==   local arrays ==
  C--   TmixLoc        :: ocean surface/mixed-layer temperature (in K)
        _RL TmixLoc       (1:sNx,1:sNy)
+ C     actual ice thickness (with upper and lower limit)
+       _RL heffActual          (1:sNx,1:sNy)
+ C     actual snow thickness
+       _RL hsnowActual         (1:sNx,1:sNy)
+ C     actual ice thickness (with lower limit only) Reciprocal
+       _RL recip_heffActual    (1:sNx,1:sNy)
+ C     AREA_PRE :: hold sea-ice fraction field before any seaice-thermo update
+       _RL AREApreTH           (1:sNx,1:sNy)
+       _RL HEFFpreTH           (1:sNx,1:sNy)
+       _RL HSNWpreTH           (1:sNx,1:sNy)
+ #ifdef SEAICE_ITD
+       _RL AREAITDpreTH        (1:sNx,1:sNy,1:nITD)
+       _RL HEFFITDpreTH        (1:sNx,1:sNy,1:nITD)
+       _RL HSNWITDpreTH        (1:sNx,1:sNy,1:nITD)
+       _RL areaFracFactor      (1:sNx,1:sNy,1:nITD)
+       _RL leadIceThickMin
+ #endif
+ C     wind speed
+       _RL UG                  (1:sNx,1:sNy)
+ C     temporary variables available for the various computations
+       _RL tmparr1             (1:sNx,1:sNy)
+ #ifdef SEAICE_VARIABLE_SALINITY
+       _RL saltFluxAdjust      (1:sNx,1:sNy)
+ #endif
+       _RL ticeInMult          (1:sNx,1:sNy,MULTDIM)
+       _RL ticeOutMult         (1:sNx,1:sNy,MULTDIM)
+       _RL heffActualMult      (1:sNx,1:sNy,MULTDIM)
+       _RL hsnowActualMult     (1:sNx,1:sNy,MULTDIM)
+ #ifdef SEAICE_ITD
+       _RL recip_heffActualMult(1:sNx,1:sNy,MULTDIM)
+ #endif
+       _RL a_QbyATMmult_cover  (1:sNx,1:sNy,MULTDIM)
+       _RL a_QSWbyATMmult_cover(1:sNx,1:sNy,MULTDIM)
+       _RL a_FWbySublimMult    (1:sNx,1:sNy,MULTDIM)
+ #ifdef SEAICE_ITD
+       _RL r_QbyATMmult_cover  (1:sNx,1:sNy,MULTDIM)
+       _RL r_FWbySublimMult    (1:sNx,1:sNy,MULTDIM)
+ #endif
  C     a_QbyATM_cover :: available heat (in W/m^2) due to the interaction of
  C             the atmosphere and the ocean surface - for ice covered water
  C     a_QbyATM_open  :: same but for open water
-Line 111 
 C     a_QSWbyATM_open   - short wave hea
+Line 210 
 C     a_QSWbyATM_open   - short wave hea
  C     a_QSWbyATM_cover  - short wave heat flux under ice in W/m^2
        _RL a_QSWbyATM_open     (1:sNx,1:sNy)
        _RL a_QSWbyATM_cover    (1:sNx,1:sNy)
- C     a_QbyOCN :: available heat (in in W/m^2) due to the
+ C     a_QbyOCN :: available heat (in W/m^2) due to the
  C             interaction of the ice pack and the ocean surface
  C     r_QbyOCN :: residual of a_QbyOCN after freezing/melting
  C             processes have been accounted for
        _RL a_QbyOCN            (1:sNx,1:sNy)
        _RL r_QbyOCN            (1:sNx,1:sNy)
- C conversion factors to go from Q (W/m2) to HEFF (ice meters)
+ C     The change of mean ice thickness due to out-of-bounds values following
-       _RL convertQ2HI, convertHI2Q
+ C     sea ice dyhnamics
- C conversion factors to go from precip (m/s) unit to HEFF (ice meters)
-       _RL convertPRECIP2HI, convertHI2PRECIP
- #ifdef ALLOW_DIAGNOSTICS
- C ICE/SNOW stocks tendencies associated with the various melt/freeze processes
-       _RL d_AREAbyATM         (1:sNx,1:sNy)
-       _RL d_AREAbyOCN         (1:sNx,1:sNy)
-       _RL d_AREAbyICE         (1:sNx,1:sNy)
- #endif
- #ifdef SEAICE_ALLOW_AREA_RELAXATION
- C     ICE/SNOW stocks tendency associated with relaxation towards observation
-       _RL d_AREAbyRLX         (1:sNx,1:sNy)
- c     The change of mean ice thickness due to relaxation
-       _RL d_HEFFbyRLX         (1:sNx,1:sNy)
- #endif
- #ifdef SEAICE_ITD
- c     The change of mean ice area due to out-of-bounds values following
- c     sea ice dynamics
-       _RL d_AREAbyNEG         (1:sNx,1:sNy)
- #endif
- c     The change of mean ice thickness due to out-of-bounds values following
- c     sea ice dynamics
        _RL d_HEFFbyNEG         (1:sNx,1:sNy)
- c     The change of mean ice thickness due to turbulent ocean-sea ice heat fluxes
+ C     The change of mean ice thickness due to turbulent ocean-sea ice heat fluxes
        _RL d_HEFFbyOCNonICE    (1:sNx,1:sNy)
- c     The sum of mean ice thickness increments due to atmospheric fluxes over the open water
+ C     The sum of mean ice thickness increments due to atmospheric fluxes over
- c     fraction and ice-covered fractions of the grid cell
+ C     the open water fraction and ice-covered fractions of the grid cell
        _RL d_HEFFbyATMonOCN    (1:sNx,1:sNy)
- c     The change of mean ice thickness due to flooding by snow
+ C     The change of mean ice thickness due to flooding by snow
        _RL d_HEFFbyFLOODING    (1:sNx,1:sNy)
- c     The mean ice thickness increments due to atmospheric fluxes over the open water
+ C     The mean ice thickness increments due to atmospheric fluxes over the open
- c     fraction and ice-covered fractions of the grid cell, respectively
+ C     water fraction and ice-covered fractions of the grid cell, respectively
        _RL d_HEFFbyATMonOCN_open(1:sNx,1:sNy)
        _RL d_HEFFbyATMonOCN_cover(1:sNx,1:sNy)
-Line 166 
 c     fraction and ice-covered fractions
+Line 241 
 c     fraction and ice-covered fractions
        _RL d_HSNWbyRAIN        (1:sNx,1:sNy)
        _RL d_HFRWbyRAIN        (1:sNx,1:sNy)
- C
  C     a_FWbySublim :: fresh water flux implied by latent heat of
  C                     sublimation to atmosphere, same sign convention
  C                     as EVAP (positive upward)
-Line 179 
 C                     as EVAP (positive
+Line 254 
 C                     as EVAP (positive
  C     The latent heat flux which will sublimate all snow and ice
  C     over one time step
        _RL latentHeatFluxMax   (1:sNx,1:sNy)
-       _RL latentHeatFluxMaxMult  (1:sNx,1:sNy,MULTDIM)
+       _RL latentHeatFluxMaxMult(1:sNx,1:sNy,MULTDIM)
- #endif
- C     actual ice thickness (with upper and lower limit)
-       _RL heffActual          (1:sNx,1:sNy)
- C     actual snow thickness
-       _RL hsnowActual         (1:sNx,1:sNy)
- C     actual ice thickness (with lower limit only) Reciprocal
-       _RL recip_heffActual    (1:sNx,1:sNy)
- C     local value (=1/HO or 1/HO_south)
-       _RL recip_HO
- C     local value (=1/ice thickness)
-       _RL recip_HH
- C     AREA_PRE :: hold sea-ice fraction field before any seaice-thermo update
-       _RL AREApreTH           (1:sNx,1:sNy)
-       _RL HEFFpreTH           (1:sNx,1:sNy)
-       _RL HSNWpreTH           (1:sNx,1:sNy)
- #ifdef SEAICE_ITD
-       _RL AREAITDpreTH        (1:sNx,1:sNy,1:nITD)
-       _RL HEFFITDpreTH        (1:sNx,1:sNy,1:nITD)
-       _RL HSNWITDpreTH        (1:sNx,1:sNy,1:nITD)
-       _RL areaFracFactor      (1:sNx,1:sNy,1:nITD)
-       _RL heffFracFactor      (1:sNx,1:sNy,1:nITD)
  #endif
- C     wind speed
+ #ifdef EXF_ALLOW_SEAICE_RELAX
-       _RL UG                  (1:sNx,1:sNy)
+ C     ICE/SNOW stocks tendency associated with relaxation towards observation
- #ifdef ALLOW_ATM_WIND
+       _RL d_AREAbyRLX         (1:sNx,1:sNy)
-       _RL SPEED_SQ
+ C     The change of mean ice thickness due to relaxation
- #endif
+       _RL d_HEFFbyRLX         (1:sNx,1:sNy)
- C     Regularization values squared
-       _RL area_reg_sq, hice_reg_sq
- C     pathological cases thresholds
-       _RL heffTooHeavy
-       _RL lhSublim
- C temporary variables available for the various computations
-       _RL tmpscal0, tmpscal1, tmpscal2, tmpscal3, tmpscal4
-       _RL tmparr1             (1:sNx,1:sNy)
- #ifdef SEAICE_VARIABLE_SALINITY
-       _RL saltFluxAdjust      (1:sNx,1:sNy)
  #endif
-       INTEGER ilockey
- CToM<<<
- C      INTEGER it
-       INTEGER IT, K
- C>>>ToM
-       _RL pFac
-       _RL ticeInMult          (1:sNx,1:sNy,MULTDIM)
-       _RL ticeOutMult         (1:sNx,1:sNy,MULTDIM)
-       _RL heffActualMult      (1:sNx,1:sNy,MULTDIM)
- #ifdef SEAICE_ITD
-       _RL hsnowActualMult     (1:sNx,1:sNy,MULTDIM)
-       _RL recip_heffActualMult(1:sNx,1:sNy,MULTDIM)
- #endif
-       _RL a_QbyATMmult_cover  (1:sNx,1:sNy,MULTDIM)
-       _RL a_QSWbyATMmult_cover(1:sNx,1:sNy,MULTDIM)
-       _RL a_FWbySublimMult    (1:sNx,1:sNy,MULTDIM)
  #ifdef SEAICE_ITD
-       _RL r_QbyATMmult_cover  (1:sNx,1:sNy,MULTDIM)
+       _RL d_HEFFbySublim_ITD         (1:sNx,1:sNy,1:nITD)
-       _RL r_FWbySublimMult    (1:sNx,1:sNy,MULTDIM)
+       _RL d_HSNWbySublim_ITD         (1:sNx,1:sNy,1:nITD)
+       _RL d_HEFFbyOCNonICE_ITD       (1:sNx,1:sNy,1:nITD)
+       _RL d_HSNWbyATMonSNW_ITD       (1:sNx,1:sNy,1:nITD)
+       _RL d_HEFFbyATMonOCN_ITD       (1:sNx,1:sNy,1:nITD)
+       _RL d_HEFFbyATMonOCN_cover_ITD (1:sNx,1:sNy,1:nITD)
+       _RL d_HEFFbyATMonOCN_open_ITD  (1:sNx,1:sNy,1:nITD)
+       _RL d_HSNWbyRAIN_ITD           (1:sNx,1:sNy,1:nITD)
+       _RL d_HSNWbyOCNonSNW_ITD       (1:sNx,1:sNy,1:nITD)
+       _RL d_HEFFbyFLOODING_ITD       (1:sNx,1:sNy,1:nITD)
  #endif
- C     Helper variables: reciprocal of some constants
-       _RL recip_multDim
-       _RL recip_deltaTtherm
-       _RL recip_rhoIce
- C     Factor by which we increase the upper ocean friction velocity (u*) when
- C     ice is absent in a grid cell  (dimensionless)
-       _RL MixedLayerTurbulenceFactor
- #ifdef ALLOW_SITRACER
-       INTEGER iTr
-       CHARACTER*8   diagName
- #endif
  #ifdef ALLOW_DIAGNOSTICS
- c     Helper variables for diagnostics
+ C ICE/SNOW stocks tendencies associated with the various melt/freeze processes
+       _RL d_AREAbyATM         (1:sNx,1:sNy)
+       _RL d_AREAbyOCN         (1:sNx,1:sNy)
+       _RL d_AREAbyICE         (1:sNx,1:sNy)
+ C     Helper variables for diagnostics
        _RL DIAGarrayA    (1:sNx,1:sNy)
        _RL DIAGarrayB    (1:sNx,1:sNy)
        _RL DIAGarrayC    (1:sNx,1:sNy)
        _RL DIAGarrayD    (1:sNx,1:sNy)
- #endif
+ #endif /* ALLOW_DIAGNOSTICS */
+       _RL SItflux     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
+       _RL SIatmQnt    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
+       _RL SIatmFW     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
+ #ifdef ALLOW_BALANCE_FLUXES
+       _RL FWFsiTile(nSx,nSy)
+       _RL FWFsiGlob
+       _RL HFsiTile(nSx,nSy)
+       _RL HFsiGlob
+       _RL FWF2HFsiTile(nSx,nSy)
+       _RL FWF2HFsiGlob
+ #endif
  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
-Line 282 
 C ======================================
+Line 314 
 C ======================================
        ENDIF
  C     avoid unnecessary divisions in loops
+ c#ifdef SEAICE_ITD
+ CToM this is now set by MULTDIM = nITD in SEAICE_SIZE.h
+ C    (see SEAICE_SIZE.h and seaice_readparms.F)
+ c     SEAICE_multDim = nITD
+ c#endif
        recip_multDim        = SEAICE_multDim
        recip_multDim        = ONE / recip_multDim
  C     above/below: double/single precision calculation of recip_multDim
-Line 332 
 C conversion factors to go from precip (
+Line 369 
 C conversion factors to go from precip (
       &                       + act4*max1*max2*max3
  #endif /* ALLOW_AUTODIFF_TAMC */
  C array initializations
  C =====================
-Line 355 
 C =====================
+Line 391 
 C =====================
            d_AREAbyOCN(I,J)           = 0.0 _d 0
  #endif
- #ifdef SEAICE_ALLOW_AREA_RELAXATION
+ #ifdef EXF_ALLOW_SEAICE_RELAX
            d_AREAbyRLX(I,J)           = 0.0 _d 0
            d_HEFFbyRLX(I,J)           = 0.0 _d 0
  #endif
- #ifdef SEAICE_ITD
-           d_AREAbyNEG(I,J)           = 0.0 _d 0
- #endif
            d_HEFFbyNEG(I,J)           = 0.0 _d 0
            d_HEFFbyOCNonICE(I,J)      = 0.0 _d 0
            d_HEFFbyATMonOCN(I,J)      = 0.0 _d 0
-Line 382 
 C =====================
+Line 415 
 C =====================
  #ifdef SEAICE_CAP_SUBLIM
            latentHeatFluxMax(I,J)     = 0.0 _d 0
  #endif
- c
            d_HFRWbyRAIN(I,J)          = 0.0 _d 0
            tmparr1(I,J)               = 0.0 _d 0
  #ifdef SEAICE_VARIABLE_SALINITY
            saltFluxAdjust(I,J)        = 0.0 _d 0
  #endif
-Line 396 
 c
+Line 426 
 c
              a_QbyATMmult_cover(I,J,IT)    = 0.0 _d 0
              a_QSWbyATMmult_cover(I,J,IT)  = 0.0 _d 0
              a_FWbySublimMult(I,J,IT)      = 0.0 _d 0
- #ifdef SEAICE_ITD
-             r_QbyATMmult_cover (I,J,IT)   = 0.0 _d 0
-             r_FWbySublimMult(I,J,IT)      = 0.0 _d 0
- #endif
  #ifdef SEAICE_CAP_SUBLIM
              latentHeatFluxMaxMult(I,J,IT) = 0.0 _d 0
  #endif
+ #ifdef SEAICE_ITD
+             d_HEFFbySublim_ITD(I,J,IT)         = 0.0 _d 0
+             d_HSNWbySublim_ITD(I,J,IT)         = 0.0 _d 0
+             d_HEFFbyOCNonICE_ITD(I,J,IT)       = 0.0 _d 0
+             d_HSNWbyATMonSNW_ITD(I,J,IT)       = 0.0 _d 0
+             d_HEFFbyATMonOCN_ITD(I,J,IT)       = 0.0 _d 0
+             d_HEFFbyATMonOCN_cover_ITD(I,J,IT) = 0.0 _d 0
+             d_HEFFbyATMonOCN_open_ITD(I,J,IT)  = 0.0 _d 0
+             d_HSNWbyRAIN_ITD(I,J,IT)           = 0.0 _d 0
+             d_HSNWbyOCNonSNW_ITD(I,J,IT)       = 0.0 _d 0
+             d_HEFFbyFLOODING_ITD(I,J,IT)       = 0.0 _d 0
+             r_QbyATMmult_cover(I,J,IT) = 0.0 _d 0
+             r_FWbySublimMult(I,J,IT)   = 0.0 _d 0
+ #endif
            ENDDO
           ENDDO
          ENDDO
-Line 414 
 c
+Line 454 
 c
          ENDDO
  #endif
  C =====================================================================
  C ===========PART 1: treat pathological cases (post advdiff)===========
  C =====================================================================
-Line 437 
 C ======================================
+Line 476 
 C ======================================
           ENDDO
          ENDDO
  #ifdef SEAICE_ITD
-         DO K=1,nITD
+         DO IT=1,nITD
           DO J=1,sNy
            DO I=1,sNx
-            HEFFITDpreTH(I,J,K)=HEFFITD(I,J,K,bi,bj)
+            HEFFITDpreTH(I,J,IT)=HEFFITD(I,J,IT,bi,bj)
-            HSNWITDpreTH(I,J,K)=HSNOWITD(I,J,K,bi,bj)
+            HSNWITDpreTH(I,J,IT)=HSNOWITD(I,J,IT,bi,bj)
-            AREAITDpreTH(I,J,K)=AREAITD(I,J,K,bi,bj)
+            AREAITDpreTH(I,J,IT)=AREAITD(I,J,IT,bi,bj)
-            IF (AREA(I,J,bi,bj).GT.0) THEN
-             areaFracFactor(I,J,K)=AREAITD(I,J,K,bi,bj)/AREA(I,J,bi,bj)
-            ELSE
-             areaFracFactor(I,J,K)=ZERO
-            ENDIF
-            IF (HEFF(I,J,bi,bj).GT.0) THEN
-             heffFracFactor(I,J,K)=HEFFITD(I,J,K,bi,bj)/HEFF(I,J,bi,bj)
-            ELSE
-             heffFracFactor(I,J,K)=ZERO
-            ENDIF
            ENDDO
           ENDDO
          ENDDO
-Line 465 
 Cgf no dependency through pathological c
+Line 494 
 Cgf no dependency through pathological c
          IF ( SEAICEadjMODE.EQ.0 ) THEN
  #endif
- #ifdef SEAICE_ALLOW_AREA_RELAXATION
+ #ifdef EXF_ALLOW_SEAICE_RELAX
  CADJ STORE heff(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
  CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
  C 0) relax sea ice concentration towards observation
-Line 485 
 C           d_HEFFbyRLX(i,j) = 1. _d 1 *
+Line 514 
 C           d_HEFFbyRLX(i,j) = 1. _d 1 *
              d_HEFFbyRLX(i,j) = 1. _d 1 * siEps
             ENDIF
  #ifdef SEAICE_ITD
              AREAITD(I,J,1,bi,bj) = AREAITD(I,J,1,bi,bj)
       &                           +  d_AREAbyRLX(i,j)
              HEFFITD(I,J,1,bi,bj) = HEFFITD(I,J,1,bi,bj)
       &                           +  d_HEFFbyRLX(i,j)
-Line 495 
 C           d_HEFFbyRLX(i,j) = 1. _d 1 *
+Line 524 
 C           d_HEFFbyRLX(i,j) = 1. _d 1 *
           ENDDO
          ENDDO
          ENDIF
- #endif /* SEAICE_ALLOW_AREA_RELAXATION */
+ #endif /* EXF_ALLOW_SEAICE_RELAX */
  C 1) treat the case of negative values:
-Line 504 
 CADJ STORE heff(:,:,bi,bj)  = comlev1_bi
+Line 533 
 CADJ STORE heff(:,:,bi,bj)  = comlev1_bi
  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
  CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
  #endif /* ALLOW_AUTODIFF_TAMC */
+ #ifdef SEAICE_ITD
+         DO IT=1,nITD
+ #endif
          DO J=1,sNy
           DO I=1,sNx
  #ifdef SEAICE_ITD
-           DO K=1,nITD
+           tmpscal2=0. _d 0
-            tmpscal2=0. _d 0
+           tmpscal3=0. _d 0
-            tmpscal3=0. _d 0
+           tmpscal2=MAX(-HEFFITD(I,J,IT,bi,bj),0. _d 0)
-            tmpscal4=0. _d 0
+           HEFFITD(I,J,IT,bi,bj)=HEFFITD(I,J,IT,bi,bj)+tmpscal2
-            tmpscal2=MAX(-HEFFITD(I,J,K,bi,bj),0. _d 0)
+           d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2
-            HEFFITD(I,J,K,bi,bj)=HEFFITD(I,J,K,bi,bj)+tmpscal2
+           tmpscal3=MAX(-HSNOWITD(I,J,IT,bi,bj),0. _d 0)
-            d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2
+           HSNOWITD(I,J,IT,bi,bj)=HSNOWITD(I,J,IT,bi,bj)+tmpscal3
-            tmpscal3=MAX(-HSNOWITD(I,J,K,bi,bj),0. _d 0)
+           d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3
-            HSNOWITD(I,J,K,bi,bj)=HSNOWITD(I,J,K,bi,bj)+tmpscal3
+           AREAITD(I,J,IT,bi,bj)=MAX(AREAITD(I,J,IT,bi,bj),0. _d 0)
-            d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3
+ CToM      AREA, HEFF, and HSNOW will be updated at end of PART 1
-            tmpscal4=MAX(-AREAITD(I,J,K,bi,bj),0. _d 0)
+ C         by calling SEAICE_ITD_SUM
-            AREAITD(I,J,K,bi,bj)=AREAITD(I,J,K,bi,bj)+tmpscal4
-            d_AREAbyNEG(I,J)=d_AREAbyNEG(I,J)+tmpscal4
-           ENDDO
-           AREA(I,J,bi,bj)=AREA(I,J,bi,bj)+d_AREAbyNEG(I,J)
  #else
            d_HEFFbyNEG(I,J)=MAX(-HEFF(I,J,bi,bj),0. _d 0)
+           HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+d_HEFFbyNEG(I,J)
            d_HSNWbyNEG(I,J)=MAX(-HSNOW(I,J,bi,bj),0. _d 0)
+           HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+d_HSNWbyNEG(I,J)
            AREA(I,J,bi,bj)=MAX(AREA(I,J,bi,bj),0. _d 0)
  #endif
-           HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+d_HEFFbyNEG(I,J)
-           HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+d_HSNWbyNEG(I,J)
           ENDDO
          ENDDO
+ #ifdef SEAICE_ITD
+         ENDDO
+ #endif
  C 1.25) treat the case of very thin ice:
  #ifdef ALLOW_AUTODIFF_TAMC
  CADJ STORE heff(:,:,bi,bj)  = comlev1_bibj, key = iicekey,byte=isbyte
  #endif /* ALLOW_AUTODIFF_TAMC */
+ #ifdef SEAICE_ITD
+         DO IT=1,nITD
+ #endif
          DO J=1,sNy
           DO I=1,sNx
- #ifdef SEAICE_ITD
-           DO K=1,nITD
-            tmpscal2=0. _d 0
-            tmpscal3=0. _d 0
-            IF (HEFFITD(I,J,K,bi,bj).LE.siEps) THEN
-             tmpscal2=-HEFFITD(I,J,K,bi,bj)
-             tmpscal3=-HSNOWITD(I,J,K,bi,bj)
-             TICES(I,J,K,bi,bj)=celsius2K
-             HEFFITD(I,J,K,bi,bj) =HEFFITD(I,J,K,bi,bj) +tmpscal2
-             HSNOWITD(I,J,K,bi,bj)=HSNOWITD(I,J,K,bi,bj)+tmpscal3
- c
-             TICE(I,J,bi,bj)=celsius2K
- c
-             HEFF(I,J,bi,bj) =HEFF(I,J,bi,bj) +tmpscal2
-             d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2
-             HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+tmpscal3
-             d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3
-            ENDIF
-           ENDDO
- #else
            tmpscal2=0. _d 0
            tmpscal3=0. _d 0
+ #ifdef SEAICE_ITD
+           IF (HEFFITD(I,J,IT,bi,bj).LE.siEps) THEN
+            tmpscal2=-HEFFITD(I,J,IT,bi,bj)
+            tmpscal3=-HSNOWITD(I,J,IT,bi,bj)
+            TICES(I,J,IT,bi,bj)=celsius2K
+ CToM       TICE will be updated at end of Part 1 together with AREA and HEFF
+           ENDIF
+           HEFFITD(I,J,IT,bi,bj) =HEFFITD(I,J,IT,bi,bj) +tmpscal2
+           HSNOWITD(I,J,IT,bi,bj)=HSNOWITD(I,J,IT,bi,bj)+tmpscal3
+ #else
            IF (HEFF(I,J,bi,bj).LE.siEps) THEN
-             tmpscal2=-HEFF(I,J,bi,bj)
+            tmpscal2=-HEFF(I,J,bi,bj)
-             tmpscal3=-HSNOW(I,J,bi,bj)
+            tmpscal3=-HSNOW(I,J,bi,bj)
-             TICE(I,J,bi,bj)=celsius2K
+            TICE(I,J,bi,bj)=celsius2K
-             DO IT=1,SEAICE_multDim
+            DO IT=1,SEAICE_multDim
-               TICES(I,J,IT,bi,bj)=celsius2K
+             TICES(I,J,IT,bi,bj)=celsius2K
-             ENDDO
+            ENDDO
            ENDIF
            HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+tmpscal2
-           d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2
            HSNOW(I,J,bi,bj)=HSNOW(I,J,bi,bj)+tmpscal3
-           d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3
  #endif
+           d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2
+           d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3
           ENDDO
          ENDDO
+ #ifdef SEAICE_ITD
+         ENDDO
+ #endif
  C 1.5) treat the case of area but no ice/snow:
-Line 583 
 C 1.5) treat the case of area but no ice
+Line 610 
 C 1.5) treat the case of area but no ice
  CADJ STORE heff(:,:,bi,bj)  = comlev1_bibj, key = iicekey,byte=isbyte
  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
  #endif /* ALLOW_AUTODIFF_TAMC */
+ #ifdef SEAICE_ITD
+         DO IT=1,nITD
+ #endif
          DO J=1,sNy
           DO I=1,sNx
- CToM<<<
-           IF ((HEFF(i,j,bi,bj).EQ.0. _d 0).AND.
- C     &        (HSNOW(i,j,bi,bj).EQ.0. _d 0)) AREA(I,J,bi,bj)=0. _d 0
-      &        (HSNOW(i,j,bi,bj).EQ.0. _d 0)) THEN
-            AREA(I,J,bi,bj)=0. _d 0
  #ifdef SEAICE_ITD
-            DO K=1,nITD
+           IF ((HEFFITD(I,J,IT,bi,bj).EQ.0. _d 0).AND.
-             AREAITD(I,J,K,bi,bj)=0. _d 0
+      &        (HSNOWITD(I,J,IT,bi,bj).EQ.0. _d 0))
-             HEFFITD(I,J,K,bi,bj)=0. _d 0
+      &     AREAITD(I,J,IT,bi,bj)=0. _d 0
-             HSNOWITD(I,J,K,bi,bj)=0. _d 0
+ #else
-            ENDDO
+           IF ((HEFF(i,j,bi,bj).EQ.0. _d 0).AND.
+      &        (HSNOW(i,j,bi,bj).EQ.0. _d 0)) AREA(I,J,bi,bj)=0. _d 0
  #endif
-           ENDIF
- C>>>ToM
           ENDDO
          ENDDO
+ #ifdef SEAICE_ITD
+         ENDDO
+ #endif
  C 2) treat the case of very small area:
-Line 608 
 C 2) treat the case of very small area:
+Line 635 
 C 2) treat the case of very small area:
  #ifdef ALLOW_AUTODIFF_TAMC
  CADJ STORE area(:,:,bi,bj)  = comlev1_bibj, key = iicekey,byte=isbyte
  #endif /* ALLOW_AUTODIFF_TAMC */
+ #ifdef SEAICE_ITD
+         DO IT=1,nITD
+ #endif
          DO J=1,sNy
           DO I=1,sNx
-           IF ((HEFF(i,j,bi,bj).GT.0).OR.(HSNOW(i,j,bi,bj).GT.0)) THEN
  #ifdef SEAICE_ITD
-            tmpscal2=AREA(I,J,bi,bj)
+           IF ((HEFFITD(I,J,IT,bi,bj).GT.0).OR.
- #endif
+      &        (HSNOWITD(I,J,IT,bi,bj).GT.0)) THEN
+ CToM       SEAICE_area_floor*nITD cannot be allowed to exceed 1
+ C          hence use SEAICE_area_floor devided by nITD
+ C          (or install a warning in e.g. seaice_readparms.F)
+            AREAITD(I,J,IT,bi,bj)=
+      &        MAX(AREAITD(I,J,IT,bi,bj),SEAICE_area_floor/float(nITD))
+           ENDIF
+ #else
+           IF ((HEFF(i,j,bi,bj).GT.0).OR.(HSNOW(i,j,bi,bj).GT.0)) THEN
             AREA(I,J,bi,bj)=MAX(AREA(I,J,bi,bj),SEAICE_area_floor)
- #ifdef SEAICE_ITD
- c          ice area added (tmpscal3 is .ge.0):
-            tmpscal3=AREA(I,J,bi,bj)-tmpscal2
- c          distribute this gain proportionally over categories
-            DO K=1,nITD
-             tmpscal4=AREAITD(I,J,K,bi,bj)/tmpscal2*tmpscal3
-             AREAITD(I,J,K,bi,bj)=AREAITD(I,J,K,bi,bj)+tmpscal4
-            ENDDO
- #endif
            ENDIF
+ #endif
           ENDDO
          ENDDO
+ #ifdef SEAICE_ITD
+         ENDDO
+ #endif
  #endif /* DISABLE_AREA_FLOOR */
  C 2.5) treat case of excessive ice cover, e.g., due to ridging:
+ CToM   for SEAICE_ITD this case is treated in SEAICE_ITD_REDIST,
+ C      which is called at end of PART 1 below
+ #ifndef SEAICE_ITD
  #ifdef ALLOW_AUTODIFF_TAMC
  CADJ STORE area(:,:,bi,bj)  = comlev1_bibj, key = iicekey,byte=isbyte
  #endif /* ALLOW_AUTODIFF_TAMC */
          DO J=1,sNy
           DO I=1,sNx
- #ifdef SEAICE_ITD
-           tmpscal2=AREA(I,J,bi,bj)
- #endif
  #ifdef ALLOW_DIAGNOSTICS
            DIAGarrayA(I,J) = AREA(I,J,bi,bj)
  #endif
-Line 646 
 CADJ STORE area(:,:,bi,bj)  = comlev1_bi
+Line 678 
 CADJ STORE area(:,:,bi,bj)  = comlev1_bi
            SItrAREA(I,J,bi,bj,1)=AREA(I,J,bi,bj)
  #endif
            AREA(I,J,bi,bj)=MIN(AREA(I,J,bi,bj),SEAICE_area_max)
+          ENDDO
+         ENDDO
+ #endif /* notSEAICE_ITD */
  #ifdef SEAICE_ITD
- c         ice area subtracted (tmpscal3 is .ge.0):
+ CToM catch up with items 1.25 and 2.5 involving category sums AREA and HEFF
-           tmpscal3=tmpscal2-AREA(I,J,bi,bj)
+         DO IT=1,nITD
- c         distribute this loss proportionally over categories
+          DO J=1,sNy
-           DO K=1,nITD
+           DO I=1,sNx
-            AREAITD(I,J,K,bi,bj)=AREAITD(I,J,K,bi,bj)
+ C    TICES was changed above (item 1.25), now update TICE as ice volume
-      &                         -tmpscal3*areaFracFactor(I,J,K)
+ C     weighted average of TICES
-           ENDDO
+ C    also compute total of AREAITD (needed for finishing item 2.5, see below)
+            IF (IT .eq. 1) THEN
+             tmpscal1itd(i,j) = 0. _d 0
+             tmpscal2itd(i,j) = 0. _d 0
+             tmpscal3itd(i,j) = 0. _d 0
+            ENDIF
+            tmpscal1itd(i,j)=tmpscal1itd(i,j) + TICES(I,J,IT,bi,bj)
+      &                                       * HEFFITD(I,J,IT,bi,bj)
+            tmpscal2itd(i,j)=tmpscal2itd(i,j) + HEFFITD(I,J,IT,bi,bj)
+            tmpscal3itd(i,j)=tmpscal3itd(i,j) + AREAITD(I,J,IT,bi,bj)
+            IF (IT .eq. nITD) THEN
+             TICE(I,J,bi,bj)=tmpscal1itd(i,j)/tmpscal2itd(i,j)
+ C    lines of item 2.5 that were omitted:
+ C    in 2.5 these lines are executed before "ridging" is applied to AREA
+ C    hence we execute them here before SEAICE_ITD_REDIST is called
+ C    although this means that AREA has not been completely regularized
+ #ifdef ALLOW_DIAGNOSTICS
+             DIAGarrayA(I,J) = tmpscal3itd(i,j)
+ #endif
+ #ifdef ALLOW_SITRACER
+             SItrAREA(I,J,bi,bj,1)=tmpscal3itd(i,j)
  #endif
+            ENDIF
+           ENDDO
           ENDDO
          ENDDO
+ CToM finally make sure that all categories meet their thickness limits
+ C    which includes ridging as in item 2.5
+ C    and update AREA, HEFF, and HSNOW
+         CALL SEAICE_ITD_REDIST(bi, bj, myTime, myIter, myThid)
+         CALL SEAICE_ITD_SUM(bi, bj, myTime, myIter, myThid)
+ #endif /* SEAICE_ITD */
+ #ifdef SEAICE_DEBUG
+ #ifdef SEAICE_ITD
+         WRITE(msgBufForm,'(A,I2,A)') '(A,',nITD,'F14.10)'
+ #else
+         WRITE(msgBufForm,'(A,A)') '(A,  F14.10)'
+ #endif
+         WRITE(msgBuf,msgBufForm)
+      &    ' SEAICE_GROWTH: Heff increments 0, HEFF = ',
+ #ifdef SEAICE_ITD
+      &     HEFFITD(1,1,:,bi,bj)
+ #else
+      &     HEFF(1,1,bi,bj)
+ #endif
+         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
+      &    SQUEEZE_RIGHT , myThid)
+         WRITE(msgBuf,msgBufForm)
+      &    ' SEAICE_GROWTH: Area increments 0, AREA = ',
  #ifdef SEAICE_ITD
- C If AREAITD is changed due to regularization (but HEFFITD not) then the
+      &     AREAITD(1,1,:,bi,bj)
- C actual ice thickness (HEFFITD/AREAITD) in a category can be changed so
+ #else
- C that it does not fit its category limits anymore and redistribution is necessary
+      &     AREA(1,1,bi,bj)
-         CALL SEAICE_ITD_REDIST(myTime, myIter, myThid)
- C this should not affect the respective sums (AREA, HEFF, ...)
- C ... except a non-conserving redistribution scheme is used; then call:
- c        CALL SEAICE_ITD_SUM(myTime, myIter, myThid)
  #endif
+         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
+      &    SQUEEZE_RIGHT , myThid)
+ #endif /* SEAICE_DEBUG */
  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
+ C        end SEAICEadjMODE.EQ.0 statement:
          ENDIF
  #endif
-Line 689 
 C 3) store regularized values of heff, h
+Line 772 
 C 3) store regularized values of heff, h
           ENDDO
          ENDDO
  #ifdef SEAICE_ITD
-         DO K=1,nITD
+         DO IT=1,nITD
+          DO J=1,sNy
+           DO I=1,sNx
+            HEFFITDpreTH(I,J,IT)=HEFFITD(I,J,IT,bi,bj)
+            HSNWITDpreTH(I,J,IT)=HSNOWITD(I,J,IT,bi,bj)
+            AREAITDpreTH(I,J,IT)=AREAITD(I,J,IT,bi,bj)
+ C memorize areal and volume fraction of each ITD category
+            IF (AREA(I,J,bi,bj) .GT. ZERO) THEN
+             areaFracFactor(I,J,IT)=AREAITD(I,J,IT,bi,bj)/AREA(I,J,bi,bj)
+            ELSE
+ C           if there is no ice, potential growth starts in 1st category
+             IF (IT .EQ. 1) THEN
+              areaFracFactor(I,J,IT)=ONE
+             ELSE
+              areaFracFactor(I,J,IT)=ZERO
+             ENDIF
+            ENDIF
+           ENDDO
+          ENDDO
+         ENDDO
+ #ifdef ALLOW_SITRACER
+ C prepare SItrHEFF to be computed as cumulative sum
+         DO iTr=2,5
           DO J=1,sNy
            DO I=1,sNx
-            HEFFITDpreTH(I,J,K)=HEFFITD(I,J,K,bi,bj)
+            SItrHEFF(I,J,bi,bj,iTr)=ZERO
-            HSNWITDpreTH(I,J,K)=HSNOWITD(I,J,K,bi,bj)
-            AREAITDpreTH(I,J,K)=AREAITD(I,J,K,bi,bj)
            ENDDO
           ENDDO
          ENDDO
+ C prepare SItrAREA to be computed as cumulative sum
+         DO J=1,sNy
+          DO I=1,sNx
+           SItrAREA(I,J,bi,bj,3)=ZERO
+          ENDDO
+         ENDDO
  #endif
+ #endif /* SEAICE_ITD */
  C 4) treat sea ice salinity pathological cases
  #ifdef SEAICE_VARIABLE_SALINITY
-Line 754 
 Cgf no additional dependency of air-sea
+Line 865 
 Cgf no additional dependency of air-sea
            ENDDO
           ENDDO
  #ifdef SEAICE_ITD
-          DO K=1,nITD
+          DO IT=1,nITD
            DO J=1,sNy
             DO I=1,sNx
-             HEFFITDpreTH(I,J,K) = 0. _d 0
+             HEFFITDpreTH(I,J,IT) = 0. _d 0
-             HSNWITDpreTH(I,J,K) = 0. _d 0
+             HSNWITDpreTH(I,J,IT) = 0. _d 0
-             AREAITDpreTH(I,J,K) = 0. _d 0
+             AREAITDpreTH(I,J,IT) = 0. _d 0
             ENDDO
            ENDDO
           ENDDO
-Line 784 
 CADJ STORE HEFFpreTH = comlev1_bibj, key
+Line 895 
 CADJ STORE HEFFpreTH = comlev1_bibj, key
  CADJ STORE HSNWpreTH = comlev1_bibj, key = iicekey, byte = isbyte
  #endif /* ALLOW_AUTODIFF_TAMC */
  #ifdef SEAICE_ITD
-         DO K=1,nITD
+         DO IT=1,nITD
  #endif
          DO J=1,sNy
           DO I=1,sNx
  #ifdef SEAICE_ITD
-           IF (HEFFITDpreTH(I,J,K) .GT. ZERO) THEN
+           IF (HEFFITDpreTH(I,J,IT) .GT. ZERO) THEN
  #ifdef SEAICE_GROWTH_LEGACY
-             tmpscal1          = MAX(SEAICE_area_reg,AREAITDpreTH(I,J,K))
+             tmpscal1          = MAX(SEAICE_area_reg/float(nITD),
-             hsnowActualMult(I,J,K) = HSNWITDpreTH(I,J,K)/tmpscal1
+      &                              AREAITDpreTH(I,J,IT))
-             tmpscal2               = HEFFITDpreTH(I,J,K)/tmpscal1
+             hsnowActualMult(I,J,IT) = HSNWITDpreTH(I,J,IT)/tmpscal1
-             heffActualMult(I,J,K)  = MAX(tmpscal2,SEAICE_hice_reg)
+             tmpscal2               = HEFFITDpreTH(I,J,IT)/tmpscal1
+             heffActualMult(I,J,IT)  = MAX(tmpscal2,SEAICE_hice_reg)
  #else /* SEAICE_GROWTH_LEGACY */
  cif        regularize AREA with SEAICE_area_reg
-            tmpscal1 = SQRT(AREAITDpreTH(I,J,K) * AREAITDpreTH(I,J,K)
+            tmpscal1 = SQRT(AREAITDpreTH(I,J,IT) * AREAITDpreTH(I,J,IT)
       &                     + area_reg_sq)
  cif        heffActual calculated with the regularized AREA
-            tmpscal2 = HEFFITDpreTH(I,J,K) / tmpscal1
+            tmpscal2 = HEFFITDpreTH(I,J,IT) / tmpscal1
  cif        regularize heffActual with SEAICE_hice_reg (add lower bound)
-            heffActualMult(I,J,K) = SQRT(tmpscal2 * tmpscal2
+            heffActualMult(I,J,IT) = SQRT(tmpscal2 * tmpscal2
       &                                  + hice_reg_sq)
  cif        hsnowActual calculated with the regularized AREA
-            hsnowActualMult(I,J,K) = HSNWITDpreTH(I,J,K) / tmpscal1
+            hsnowActualMult(I,J,IT) = HSNWITDpreTH(I,J,IT) / tmpscal1
  #endif /* SEAICE_GROWTH_LEGACY */
  cif        regularize the inverse of heffActual by hice_reg
-            recip_heffActualMult(I,J,K)  = AREAITDpreTH(I,J,K) /
+            recip_heffActualMult(I,J,IT)  = AREAITDpreTH(I,J,IT) /
-      &                 sqrt(HEFFITDpreTH(I,J,K) * HEFFITDpreTH(I,J,K)
+      &                 sqrt(HEFFITDpreTH(I,J,IT) * HEFFITDpreTH(I,J,IT)
       &                      + hice_reg_sq)
  cif       Do not regularize when HEFFpreTH = 0
            ELSE
-             heffActualMult(I,J,K) = ZERO
+             heffActualMult(I,J,IT) = ZERO
-             hsnowActualMult(I,J,K) = ZERO
+             hsnowActualMult(I,J,IT) = ZERO
-             recip_heffActualMult(I,J,K)  = ZERO
+             recip_heffActualMult(I,J,IT)  = ZERO
            ENDIF
- #else
+ #else /* SEAICE_ITD */
            IF (HEFFpreTH(I,J) .GT. ZERO) THEN
  #ifdef SEAICE_GROWTH_LEGACY
              tmpscal1         = MAX(SEAICE_area_reg,AREApreTH(I,J))
-Line 826 
 cif       Do not regularize when HEFFpre
+Line 938 
 cif       Do not regularize when HEFFpre
              tmpscal2         = HEFFpreTH(I,J)/tmpscal1
              heffActual(I,J)  = MAX(tmpscal2,SEAICE_hice_reg)
  #else /* SEAICE_GROWTH_LEGACY */
- cif        regularize AREA with SEAICE_area_reg
+ Cif        regularize AREA with SEAICE_area_reg
             tmpscal1 = SQRT(AREApreTH(I,J)* AREApreTH(I,J) + area_reg_sq)
- cif        heffActual calculated with the regularized AREA
+ Cif        heffActual calculated with the regularized AREA
             tmpscal2 = HEFFpreTH(I,J) / tmpscal1
- cif        regularize heffActual with SEAICE_hice_reg (add lower bound)
+ Cif        regularize heffActual with SEAICE_hice_reg (add lower bound)
             heffActual(I,J) = SQRT(tmpscal2 * tmpscal2 + hice_reg_sq)
- cif        hsnowActual calculated with the regularized AREA
+ Cif        hsnowActual calculated with the regularized AREA
             hsnowActual(I,J) = HSNWpreTH(I,J) / tmpscal1
  #endif /* SEAICE_GROWTH_LEGACY */
- cif        regularize the inverse of heffActual by hice_reg
+ Cif        regularize the inverse of heffActual by hice_reg
             recip_heffActual(I,J)  = AREApreTH(I,J) /
       &                 sqrt(HEFFpreTH(I,J)*HEFFpreTH(I,J) + hice_reg_sq)
- cif       Do not regularize when HEFFpreTH = 0
+ Cif       Do not regularize when HEFFpreTH = 0
            ELSE
              heffActual(I,J) = ZERO
              hsnowActual(I,J) = ZERO
              recip_heffActual(I,J)  = ZERO
            ENDIF
- #endif
+ #endif /* SEAICE_ITD */
           ENDDO
          ENDDO
-Line 862 
 cif       Do not regularize when HEFFpre
+Line 974 
 cif       Do not regularize when HEFFpre
  C 5) COMPUTE MAXIMUM LATENT HEAT FLUXES FOR THE CURRENT ICE
  C    AND SNOW THICKNESS
  #ifdef SEAICE_ITD
-         DO K=1,nITD
+         DO IT=1,nITD
  #endif
          DO J=1,sNy
           DO I=1,sNx
- c        The latent heat flux over the sea ice which
+ C        The latent heat flux over the sea ice which
- c        will sublimate all of the snow and ice over one time
+ C        will sublimate all of the snow and ice over one time
- c        step (W/m^2)
+ C        step (W/m^2)
  #ifdef SEAICE_ITD
-           IF (HEFFITDpreTH(I,J,K) .GT. ZERO) THEN
+           IF (HEFFITDpreTH(I,J,IT) .GT. ZERO) THEN
-             latentHeatFluxMaxMult(I,J,K) = lhSublim*recip_deltaTtherm *
+             latentHeatFluxMaxMult(I,J,IT) = lhSublim*recip_deltaTtherm *
-      &         (HEFFITDpreTH(I,J,K)*SEAICE_rhoIce +
+      &         (HEFFITDpreTH(I,J,IT)*SEAICE_rhoIce +
-      &          HSNWITDpreTH(I,J,K)*SEAICE_rhoSnow)/AREAITDpreTH(I,J,K)
+      &          HSNWITDpreTH(I,J,IT)*SEAICE_rhoSnow)
+      &         /AREAITDpreTH(I,J,IT)
            ELSE
-             latentHeatFluxMaxMult(I,J,K) = ZERO
+             latentHeatFluxMaxMult(I,J,IT) = ZERO
            ENDIF
  #else
            IF (HEFFpreTH(I,J) .GT. ZERO) THEN
-Line 925 
 cCADJ STORE TmixLoc = comlev1_bibj, key
+Line 1038 
 cCADJ STORE TmixLoc = comlev1_bibj, key
  C determine available heat due to the atmosphere -- for ice covered water
  C =======================================================================
- #ifdef ALLOW_ATM_WIND
+         IF (useRelativeWind.AND.useAtmWind) THEN
-         IF (useRelativeWind) THEN
  C     Compute relative wind speed over sea ice.
           DO J=1,sNy
            DO I=1,sNx
-Line 947 
 C     Compute relative wind speed over s
+Line 1059 
 C     Compute relative wind speed over s
            ENDDO
           ENDDO
          ENDIF
- #endif /* ALLOW_ATM_WIND */
  #ifdef ALLOW_AUTODIFF_TAMC
  CADJ STORE tice(:,:,bi,bj)
-Line 962 
 CADJ &                       key = iicek
+Line 1073 
 CADJ &                       key = iicek
  #endif /* ALLOW_AUTODIFF_TAMC */
  C--   Start loop over multi-categories
-         DO IT=1,SEAICE_multDim
+ #ifdef SEAICE_ITD
- c        homogeneous distribution between 0 and 2 x heffActual
+         DO IT=1,nITD
- #ifndef SEAICE_ITD
-          pFac = (2.0 _d 0*real(IT)-1.0 _d 0)*recip_multDim
- #endif
           DO J=1,sNy
            DO I=1,sNx
- #ifndef SEAICE_ITD
+ CToM for SEAICE_ITD heffActualMult and latentHeatFluxMaxMult are calculated above
- CToM for SEAICE_ITD heffActualMult and latentHeatFluxMaxMult are calculated above
  C    (instead of heffActual and latentHeatFluxMax)
+            ticeInMult(I,J,IT)=TICES(I,J,IT,bi,bj)
+            ticeOutMult(I,J,IT)=TICES(I,J,IT,bi,bj)
+            TICE(I,J,bi,bj) = ZERO
+            TICES(I,J,IT,bi,bj) = ZERO
+           ENDDO
+          ENDDO
+         ENDDO
+ #else
+         DO IT=1,SEAICE_multDim
+ C        homogeneous distribution between 0 and 2 x heffActual
+          pFac = (2.0 _d 0*IT - 1.0 _d 0)*recip_multDim
+          pFacSnow = 1. _d 0
+          IF ( SEAICE_useMultDimSnow ) pFacSnow=pFac
+          DO J=1,sNy
+           DO I=1,sNx
             heffActualMult(I,J,IT)= heffActual(I,J)*pFac
+            hsnowActualMult(I,J,IT)=hsnowActual(I,J)*pFacSnow
  #ifdef SEAICE_CAP_SUBLIM
             latentHeatFluxMaxMult(I,J,IT) = latentHeatFluxMax(I,J)*pFac
  #endif
- #endif
             ticeInMult(I,J,IT)=TICES(I,J,IT,bi,bj)
             ticeOutMult(I,J,IT)=TICES(I,J,IT,bi,bj)
             TICE(I,J,bi,bj) = ZERO
-Line 984 
 C    (instead of heffActual and latentHe
+Line 1106 
 C    (instead of heffActual and latentHe
            ENDDO
           ENDDO
          ENDDO
+ #endif
  #ifdef ALLOW_AUTODIFF_TAMC
  CADJ STORE heffActualMult = comlev1_bibj, key = iicekey, byte = isbyte
+ CADJ STORE hsnowActualMult= comlev1_bibj, key = iicekey, byte = isbyte
  CADJ STORE ticeInMult     = comlev1_bibj, key = iicekey, byte = isbyte
  # ifdef SEAICE_CAP_SUBLIM
  CADJ STORE latentHeatFluxMaxMult
-Line 1002 
 CADJ &     comlev1_bibj, key = iicekey,
+Line 1126 
 CADJ &     comlev1_bibj, key = iicekey,
          DO IT=1,SEAICE_multDim
           CALL SEAICE_SOLVE4TEMP(
- #ifdef SEAICE_ITD
       I        UG, heffActualMult(1,1,IT), hsnowActualMult(1,1,IT),
- #else
-      I        UG, heffActualMult(1,1,IT), hsnowActual,
- #endif
  #ifdef SEAICE_CAP_SUBLIM
       I        latentHeatFluxMaxMult(1,1,IT),
  #endif
       U        ticeInMult(1,1,IT), ticeOutMult(1,1,IT),
-      O        a_QbyATMmult_cover(1,1,IT), a_QSWbyATMmult_cover(1,1,IT),
+      O        a_QbyATMmult_cover(1,1,IT),
+      O        a_QSWbyATMmult_cover(1,1,IT),
       O        a_FWbySublimMult(1,1,IT),
       I        bi, bj, myTime, myIter, myThid )
          ENDDO
  #ifdef ALLOW_AUTODIFF_TAMC
  CADJ STORE heffActualMult = comlev1_bibj, key = iicekey, byte = isbyte
+ CADJ STORE hsnowActualMult= comlev1_bibj, key = iicekey, byte = isbyte
  CADJ STORE ticeOutMult    = comlev1_bibj, key = iicekey, byte = isbyte
  # ifdef SEAICE_CAP_SUBLIM
  CADJ STORE latentHeatFluxMaxMult
-Line 1037 
 CADJ &     comlev1_bibj, key = iicekey,
+Line 1159 
 CADJ &     comlev1_bibj, key = iicekey,
  C     update TICE & TICES
  #ifdef SEAICE_ITD
  C     calculate area weighted mean
+ C     (although the ice temperature relates to its energy content
+ C      and hence should be averaged weighted by ice volume,
+ C      the temperature here is a result of the fluxes through the ice surface
+ C      computed individually for each single category in SEAICE_SOLVE4TEMP
+ C      and hence is averaged area weighted [areaFracFactor])
             TICE(I,J,bi,bj) = TICE(I,J,bi,bj)
-      &          +  ticeOutMult(I,J,IT)*areaFracFactor(I,J,K)
+      &          +  ticeOutMult(I,J,IT)*areaFracFactor(I,J,IT)
  #else
             TICE(I,J,bi,bj) = TICE(I,J,bi,bj)
       &          +  ticeOutMult(I,J,IT)*recip_multDim
-Line 1047 
 C     calculate area weighted mean
+Line 1174 
 C     calculate area weighted mean
  C     average over categories
  #ifdef SEAICE_ITD
  C     calculate area weighted mean
+ C     (fluxes are per unit (ice surface) area and are thus area weighted)
             a_QbyATM_cover   (I,J) = a_QbyATM_cover(I,J)
-      &          + a_QbyATMmult_cover(I,J,IT)*areaFracFactor(I,J,K)
+      &          + a_QbyATMmult_cover(I,J,IT)*areaFracFactor(I,J,IT)
             a_QSWbyATM_cover (I,J) = a_QSWbyATM_cover(I,J)
-      &          + a_QSWbyATMmult_cover(I,J,IT)*areaFracFactor(I,J,K)
+      &          + a_QSWbyATMmult_cover(I,J,IT)*areaFracFactor(I,J,IT)
             a_FWbySublim     (I,J) = a_FWbySublim(I,J)
-      &          + a_FWbySublimMult(I,J,IT)*areaFracFactor(I,J,K)
+      &          + a_FWbySublimMult(I,J,IT)*areaFracFactor(I,J,IT)
  #else
             a_QbyATM_cover   (I,J) = a_QbyATM_cover(I,J)
       &          + a_QbyATMmult_cover(I,J,IT)*recip_multDim
-Line 1069 
 C     calculate area weighted mean
+Line 1197 
 C     calculate area weighted mean
  # ifdef ALLOW_DIAGNOSTICS
          DO J=1,sNy
           DO I=1,sNx
- c          The actual latent heat flux realized by SOLVE4TEMP
+ C          The actual latent heat flux realized by SOLVE4TEMP
             DIAGarrayA(I,J) = a_FWbySublim(I,J) * lhSublim
           ENDDO
          ENDDO
- cif     The actual vs. maximum latent heat flux
+ Cif     The actual vs. maximum latent heat flux
          IF ( useDiagnostics ) THEN
            CALL DIAGNOSTICS_FILL(DIAGarrayA,
       &     'SIactLHF',0,1,3,bi,bj,myThid)
-Line 1093 
 CADJ STORE a_FWbySublim    = comlev1_bib
+Line 1221 
 CADJ STORE a_FWbySublim    = comlev1_bib
  #endif /* ALLOW_AUTODIFF_TAMC */
  C switch heat fluxes from W/m2 to 'effective' ice meters
+ #ifdef SEAICE_ITD
+         DO IT=1,nITD
+          DO J=1,sNy
+           DO I=1,sNx
+            a_QbyATMmult_cover(I,J,IT)   = a_QbyATMmult_cover(I,J,IT)
+      &          * convertQ2HI * AREAITDpreTH(I,J,IT)
+            a_QSWbyATMmult_cover(I,J,IT) = a_QSWbyATMmult_cover(I,J,IT)
+      &          * convertQ2HI * AREAITDpreTH(I,J,IT)
+ C and initialize r_QbyATMmult_cover
+            r_QbyATMmult_cover(I,J,IT)=a_QbyATMmult_cover(I,J,IT)
+ C     Convert fresh water flux by sublimation to 'effective' ice meters.
+ C     Negative sublimation is resublimation and will be added as snow.
+ #ifdef SEAICE_DISABLE_SUBLIM
+            a_FWbySublimMult(I,J,IT) = ZERO
+ #endif
+            a_FWbySublimMult(I,J,IT) = SEAICE_deltaTtherm*recip_rhoIce
+      &            * a_FWbySublimMult(I,J,IT)*AREAITDpreTH(I,J,IT)
+            r_FWbySublimMult(I,J,IT)=a_FWbySublimMult(I,J,IT)
+           ENDDO
+          ENDDO
+         ENDDO
+         DO J=1,sNy
+          DO I=1,sNx
+           a_QbyATM_open(I,J)    = a_QbyATM_open(I,J)
+      &         * convertQ2HI * ( ONE - AREApreTH(I,J) )
+           a_QSWbyATM_open(I,J)  = a_QSWbyATM_open(I,J)
+      &         * convertQ2HI * ( ONE - AREApreTH(I,J) )
+ C and initialize r_QbyATM_open
+           r_QbyATM_open(I,J)=a_QbyATM_open(I,J)
+          ENDDO
+         ENDDO
+ #else /* SEAICE_ITD */
          DO J=1,sNy
           DO I=1,sNx
            a_QbyATM_cover(I,J)   = a_QbyATM_cover(I,J)
-Line 1109 
 C and initialize r_QbyATM_cover/r_QbyATM
+Line 1269 
 C and initialize r_QbyATM_cover/r_QbyATM
  C     Convert fresh water flux by sublimation to 'effective' ice meters.
  C     Negative sublimation is resublimation and will be added as snow.
  #ifdef SEAICE_DISABLE_SUBLIM
- cgf just for those who may need to omit this term to reproduce old results
+ Cgf just for those who may need to omit this term to reproduce old results
            a_FWbySublim(I,J) = ZERO
- #endif
+ #endif /* SEAICE_DISABLE_SUBLIM */
            a_FWbySublim(I,J) = SEAICE_deltaTtherm*recip_rhoIce
       &           * a_FWbySublim(I,J)*AREApreTH(I,J)
            r_FWbySublim(I,J)=a_FWbySublim(I,J)
           ENDDO
          ENDDO
- #ifdef SEAICE_ITD
+ #endif /* SEAICE_ITD */
-         DO K=1,nITD
-          DO J=1,sNy
-           DO I=1,sNx
-            a_QbyATMmult_cover(I,J,K)   = a_QbyATMmult_cover(I,J,K)
-      &          * convertQ2HI * AREAITDpreTH(I,J,K)
-            a_QSWbyATMmult_cover(I,J,K) = a_QSWbyATMmult_cover(I,J,K)
-      &          * convertQ2HI * AREAITDpreTH(I,J,K)
-            r_QbyATMmult_cover(I,J,K)=a_QbyATMmult_cover(I,J,K)
- #ifdef SEAICE_DISABLE_SUBLIM
-            a_FWbySublimMult(I,J,K) = ZERO
- #endif
-            a_FWbySublimMult(I,J,K) = SEAICE_deltaTtherm*recip_rhoIce
-      &            * a_FWbySublimMult(I,J,K)*AREAITDpreTH(I,J,K)
-            r_FWbySublimMult(I,J,K)=a_FWbySublimMult(I,J,K)
-           ENDDO
-          ENDDO
-         ENDDO
- #endif
  #ifdef ALLOW_AUTODIFF_TAMC
  CADJ STORE AREApreTH       = comlev1_bibj, key = iicekey, byte = isbyte
-Line 1152 
 CADJ STORE r_FWbySublim    = comlev1_bib
+Line 1294 
 CADJ STORE r_FWbySublim    = comlev1_bib
  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
  Cgf no additional dependency through ice cover
          IF ( SEAICEadjMODE.GE.3 ) THEN
+ #ifdef SEAICE_ITD
+          DO IT=1,nITD
+           DO J=1,sNy
+            DO I=1,sNx
+             a_QbyATMmult_cover(I,J,IT)   = 0. _d 0
+             r_QbyATMmult_cover(I,J,IT)   = 0. _d 0
+             a_QSWbyATMmult_cover(I,J,IT) = 0. _d 0
+            ENDDO
+           ENDDO
+          ENDDO
+ #else
           DO J=1,sNy
            DO I=1,sNx
             a_QbyATM_cover(I,J)   = 0. _d 0
-Line 1159 
 Cgf no additional dependency through ice
+Line 1312 
 Cgf no additional dependency through ice
             a_QSWbyATM_cover(I,J) = 0. _d 0
            ENDDO
           ENDDO
- #ifdef SEAICE_ITD
-          DO K=1,nITD
-           DO J=1,sNy
-            DO I=1,sNx
-             a_QbyATMmult_cover(I,J,K)   = 0. _d 0
-             r_QbyATMmult_cover(I,J,K)   = 0. _d 0
-             a_QSWbyATMmult_cover(I,J,K) = 0. _d 0
-            ENDDO
-           ENDDO
-          ENDDO
  #endif
          ENDIF
  #endif
-Line 1186 
 CADJ &                       key = iicek
+Line 1329 
 CADJ &                       key = iicek
          DO J=1,sNy
           DO I=1,sNx
- c         FREEZING TEMP. OF SEA WATER (deg C)
+ C         FREEZING TEMP. OF SEA WATER (deg C)
            tempFrz = SEAICE_tempFrz0 +
       &              SEAICE_dTempFrz_dS *salt(I,J,kSurface,bi,bj)
- c efficiency of turbulent fluxes : dependency to sign of THETA-TBC
+ C efficiency of turbulent fluxes : dependency to sign of THETA-TBC
            IF ( theta(I,J,kSurface,bi,bj) .GE. tempFrz ) THEN
             tmpscal1 = SEAICE_mcPheePiston
            ELSE
             tmpscal1 =SEAICE_frazilFrac*drF(kSurface)/SEAICE_deltaTtherm
            ENDIF
- c efficiency of turbulent fluxes : dependency to AREA (McPhee cases)
+ C efficiency of turbulent fluxes : dependency to AREA (McPhee cases)
            IF ( (AREApreTH(I,J) .GT. 0. _d 0).AND.
       &         (.NOT.SEAICE_mcPheeStepFunc) ) THEN
             MixedLayerTurbulenceFactor = ONE -
-Line 1206 
 c efficiency of turbulent fluxes : depen
+Line 1349 
 c efficiency of turbulent fluxes : depen
            ELSE
             MixedLayerTurbulenceFactor = ONE
            ENDIF
- c maximum turbulent flux, in ice meters
+ C maximum turbulent flux, in ice meters
            tmpscal2= - (HeatCapacity_Cp*rhoConst * recip_QI)
       &         * (theta(I,J,kSurface,bi,bj)-tempFrz)
       &         * SEAICE_deltaTtherm * maskC(i,j,kSurface,bi,bj)
- c available turbulent flux
+ C available turbulent flux
            a_QbyOCN(i,j) =
       &         tmpscal1 * tmpscal2 * MixedLayerTurbulenceFactor
            r_QbyOCN(i,j) = a_QbyOCN(i,j)
-Line 1221 
 c available turbulent flux
+Line 1364 
 c available turbulent flux
          CALL ZERO_ADJ_1D( sNx*sNy, r_QbyOCN, myThid)
  #endif
  C ===================================================================
  C =========PART 3: determine effective thicknesses increments========
  C ===================================================================
-Line 1234 
 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
+Line 1376 
 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
  CADJ STORE r_FWbySublim     = comlev1_bibj,key=iicekey,byte=isbyte
  #endif /* ALLOW_AUTODIFF_TAMC */
  #ifdef SEAICE_ITD
-         DO K=1,nITD
+         DO IT=1,nITD
  #endif
          DO J=1,sNy
           DO I=1,sNx
  C     First sublimate/deposite snow
            tmpscal2 =
  #ifdef SEAICE_ITD
-      &     MAX(MIN(r_FWbySublimMult(I,J,K),HSNOWITD(I,J,K,bi,bj)
+      &     MAX(MIN(r_FWbySublimMult(I,J,IT),HSNOWITD(I,J,IT,bi,bj)
       &             *SNOW2ICE),ZERO)
+           d_HSNWbySublim_ITD(I,J,IT) = - tmpscal2 * ICE2SNOW
  C         accumulate change over ITD categories
-           d_HSNWbySublim(I,J)     = d_HSNWbySublim(I,J)     - tmpscal2
+           d_HSNWbySublim(I,J)     = d_HSNWbySublim(I,J)      - tmpscal2
       &                                                       *ICE2SNOW
-           HSNOWITD(I,J,K,bi,bj)   = HSNOWITD(I,J,K,bi,bj)   - tmpscal2
+           r_FWbySublimMult(I,J,IT)= r_FWbySublimMult(I,J,IT) - tmpscal2
-      &                                                       *ICE2SNOW
-           r_FWbySublimMult(I,J,K) = r_FWbySublimMult(I,J,K) - tmpscal2
- C         keep total up to date, too
-           r_FWbySublim(I,J)       = r_FWbySublim(I,J)       - tmpscal2
  #else
       &     MAX(MIN(r_FWbySublim(I,J),HSNOW(I,J,bi,bj)*SNOW2ICE),ZERO)
            d_HSNWbySublim(I,J) = - tmpscal2 * ICE2SNOW
-Line 1268 
 CADJ STORE r_FWbySublim    = comlev1_bib
+Line 1407 
 CADJ STORE r_FWbySublim    = comlev1_bib
  C     If anything is left, sublimate ice
            tmpscal2 =
  #ifdef SEAICE_ITD
-      &     MAX(MIN(r_FWbySublimMult(I,J,K),HEFFITD(I,J,K,bi,bj)),ZERO)
+      &     MAX(MIN(r_FWbySublimMult(I,J,IT),HEFFITD(I,J,IT,bi,bj)),ZERO)
-           d_HSNWbySublim(I,J)     = d_HSNWbySublim(I,J)     - tmpscal2
+           d_HEFFbySublim_ITD(I,J,IT) = - tmpscal2
-           HEFFITD(I,J,K,bi,bj)    = HEFFITD(I,J,K,bi,bj)    - tmpscal2
+ C         accumulate change over ITD categories
-           r_FWbySublimMult(I,J,K) = r_FWbySublimMult(I,J,K) - tmpscal2
+           d_HEFFbySublim(I,J)      = d_HEFFbySublim(I,J)      - tmpscal2
- C         keep total up to date, too
+           r_FWbySublimMult(I,J,IT) = r_FWbySublimMult(I,J,IT) - tmpscal2
-           r_FWbySublim(I,J)       = r_FWbySublim(I,J)       - tmpscal2
  #else
       &     MAX(MIN(r_FWbySublim(I,J),HEFF(I,J,bi,bj)),ZERO)
            d_HEFFbySublim(I,J) = - tmpscal2
-Line 1288 
 C     If anything is left, it will be ev
+Line 1426 
 C     If anything is left, it will be ev
  C     Since a_QbyATM_cover was computed for sublimation, not simple evaporation, we need to
  C     remove the fusion part for the residual (that happens to be precisely r_FWbySublim).
  #ifdef SEAICE_ITD
-           a_QbyATMmult_cover(I,J,K) = a_QbyATMmult_cover(I,J,K)
+           a_QbyATMmult_cover(I,J,IT) = a_QbyATMmult_cover(I,J,IT)
-      &                              - r_FWbySublimMult(I,J,K)
+      &                               - r_FWbySublimMult(I,J,IT)
-           r_QbyATMmult_cover(I,J,K) = r_QbyATMmult_cover(I,J,K)
+           r_QbyATMmult_cover(I,J,IT) = r_QbyATMmult_cover(I,J,IT)
-      &                              - r_FWbySublimMult(I,J,K)
+      &                               - r_FWbySublimMult(I,J,IT)
-          ENDDO
+ #else
-         ENDDO
- C       end K loop
-         ENDDO
- C       then update totals
-         DO J=1,sNy
-          DO I=1,sNx
- #endif
            a_QbyATM_cover(I,J) = a_QbyATM_cover(I,J)-r_FWbySublim(I,J)
            r_QbyATM_cover(I,J) = r_QbyATM_cover(I,J)-r_FWbySublim(I,J)
+ #endif
           ENDDO
          ENDDO
+ #ifdef SEAICE_ITD
+ C       end IT loop
+         ENDDO
+ #endif
+ #ifdef SEAICE_DEBUG
+ c ToM<<< debug seaice_growth
+         WRITE(msgBuf,msgBufForm)
+      &    ' SEAICE_GROWTH: Hsnow increments 1, d_HSNWySublim = ',
+ #ifdef SEAICE_ITD
+      &     d_HSNWbySublim_ITD(1,1,:)
+ #else
+      &     d_HSNWbySublim(1,1)
+ #endif
+         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
+      &    SQUEEZE_RIGHT , myThid)
+         WRITE(msgBuf,msgBufForm)
+      &    ' SEAICE_GROWTH: Heff increments 1, d_HEFFbySublim = ',
+ #ifdef SEAICE_ITD
+      &     d_HEFFbySublim_ITD(1,1,:)
+ #else
+      &     d_HEFFbySublim(1,1)
+ #endif
+         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
+      &    SQUEEZE_RIGHT , myThid)
+ c ToM>>>
+ #endif /* SEAICE_DEBUG */
  C compute ice thickness tendency due to ice-ocean interaction
  C ===========================================================
-Line 1314 
 CADJ STORE r_QbyOCN = comlev1_bibj,key=i
+Line 1472 
 CADJ STORE r_QbyOCN = comlev1_bibj,key=i
  #endif /* ALLOW_AUTODIFF_TAMC */
  #ifdef SEAICE_ITD
-         DO K=1,nITD
+         DO IT=1,nITD
           DO J=1,sNy
            DO I=1,sNx
- C          ice growth/melt due to ocean heat is equally distributed under the ice
+ C          ice growth/melt due to ocean heat r_QbyOCN (W/m^2) is
- C          and hence weighted by fractional area of each thickness category
+ C          equally distributed under the ice and hence weighted by
-            tmpscal1=MAX(r_QbyOCN(i,j)*areaFracFactor(I,J,K),
+ C          fractional area of each thickness category
-      &                               -HEFFITD(I,J,K,bi,bj))
+            tmpscal1=MAX(r_QbyOCN(i,j)*areaFracFactor(I,J,IT),
-            HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,bi,bj) + tmpscal1
+      &                               -HEFFITD(I,J,IT,bi,bj))
-            d_HEFFbyOCNonICE(I,J)=d_HEFFbyOCNonICE(I,J) + tmpscal1
+            d_HEFFbyOCNonICE_ITD(I,J,IT)=tmpscal1
+            d_HEFFbyOCNonICE(I,J) = d_HEFFbyOCNonICE(I,J) + tmpscal1
            ENDDO
           ENDDO
          ENDDO
+ #ifdef ALLOW_SITRACER
+         DO J=1,sNy
+          DO I=1,sNx
+           SItrHEFF(I,J,bi,bj,2) = HEFFpreTH(I,J)
+      &                          + d_HEFFbySublim(I,J)
+      &                          + d_HEFFbyOCNonICE(I,J)
+          ENDDO
+         ENDDO
  #endif
          DO J=1,sNy
           DO I=1,sNx
- #ifndef SEAICE_ITD
+           r_QbyOCN(I,J)=r_QbyOCN(I,J)-d_HEFFbyOCNonICE(I,J)
+          ENDDO
+         ENDDO
+ #else /* SEAICE_ITD */
+         DO J=1,sNy
+          DO I=1,sNx
            d_HEFFbyOCNonICE(I,J)=MAX(r_QbyOCN(i,j), -HEFF(I,J,bi,bj))
- #endif
            r_QbyOCN(I,J)=r_QbyOCN(I,J)-d_HEFFbyOCNonICE(I,J)
            HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj) + d_HEFFbyOCNonICE(I,J)
  #ifdef ALLOW_SITRACER
-Line 1339 
 C          and hence weighted by fractio
+Line 1510 
 C          and hence weighted by fractio
  #endif
           ENDDO
          ENDDO
+ #endif /* SEAICE_ITD */
+ #ifdef SEAICE_DEBUG
+ c ToM<<< debug seaice_growth
+         WRITE(msgBuf,msgBufForm)
+      &    ' SEAICE_GROWTH: Heff increments 2, d_HEFFbyOCNonICE = ',
+ #ifdef SEAICE_ITD
+      &     d_HEFFbyOCNonICE_ITD(1,1,:)
+ #else
+      &     d_HEFFbyOCNonICE(1,1)
+ #endif
+         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
+      &    SQUEEZE_RIGHT , myThid)
+ c ToM>>>
+ #endif /* SEAICE_DEBUG */
  C compute snow melt tendency due to snow-atmosphere interaction
  C ==================================================================
-Line 1349 
 CADJ STORE r_QbyATM_cover = comlev1_bibj
+Line 1534 
 CADJ STORE r_QbyATM_cover = comlev1_bibj
  #endif /* ALLOW_AUTODIFF_TAMC */
  #ifdef SEAICE_ITD
-         DO K=1,nITD
+         DO IT=1,nITD
           DO J=1,sNy
            DO I=1,sNx
  C     Convert to standard units (meters of ice) rather than to meters
  C     of snow. This appears to be more robust.
-           tmpscal1=MAX(r_QbyATMmult_cover(I,J,K),-HSNOWITD(I,J,K,bi,bj)
+            tmpscal1=MAX(r_QbyATMmult_cover(I,J,IT),
-      &                                           *SNOW2ICE)
+      &                  -HSNOWITD(I,J,IT,bi,bj)*SNOW2ICE)
-           tmpscal2=MIN(tmpscal1,0. _d 0)
+            tmpscal2=MIN(tmpscal1,0. _d 0)
  #ifdef SEAICE_MODIFY_GROWTH_ADJ
  Cgf no additional dependency through snow
-           IF ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
+            IF ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
  #endif
-           d_HSNWbyATMonSNW(I,J)=d_HSNWbyATMonSNW(I,J)+tmpscal2*ICE2SNOW
+            d_HSNWbyATMonSNW_ITD(I,J,IT) = tmpscal2*ICE2SNOW
-           r_QbyATMmult_cover(I,J,K)=r_QbyATMmult_cover(I,J,K) - tmpscal2
+            d_HSNWbyATMonSNW(I,J) = d_HSNWbyATMonSNW(I,J)
- C         keep the total up to date, too
+      &                           + tmpscal2*ICE2SNOW
-           r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J) - tmpscal2
+            r_QbyATMmult_cover(I,J,IT)=r_QbyATMmult_cover(I,J,IT)
-           ENDDO
+      &                           - tmpscal2
-          ENDDO
+           ENDDO
-         ENDDO
+          ENDDO
- #else
+         ENDDO
+ #else /* SEAICE_ITD */
          DO J=1,sNy
           DO I=1,sNx
  C     Convert to standard units (meters of ice) rather than to meters
-Line 1380 
 Cgf no additional dependency through sno
+Line 1566 
 Cgf no additional dependency through sno
            IF ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
  #endif
            d_HSNWbyATMonSNW(I,J)= tmpscal2*ICE2SNOW
+           HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + tmpscal2*ICE2SNOW
            r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J) - tmpscal2
           ENDDO
          ENDDO
  #endif /* SEAICE_ITD */
-         DO J=1,sNy
+ #ifdef SEAICE_DEBUG
-          DO I=1,sNx
+ c ToM<<< debug seaice_growth
-           HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + d_HSNWbyATMonSNW(I,J)
+         WRITE(msgBuf,msgBufForm)
-          ENDDO
+      &    ' SEAICE_GROWTH: Hsnow increments 3, d_HSNWbyATMonSNW = ',
-         ENDDO
+ #ifdef SEAICE_ITD
+      &    d_HSNWbyATMonSNW_ITD(1,1,:)
+ #else
+      &    d_HSNWbyATMonSNW(1,1)
+ #endif
+         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
+      &    SQUEEZE_RIGHT , myThid)
+ c ToM>>>
+ #endif /* SEAICE_DEBUG */
  C compute ice thickness tendency due to the atmosphere
  C ====================================================
-Line 1404 
 Cgf the v1.81=>v1.82 revision would chan
+Line 1599 
 Cgf the v1.81=>v1.82 revision would chan
  Cgf warming conditions, the lab_sea results were not changed.
  #ifdef SEAICE_ITD
-         DO K=1,nITD
+         DO IT=1,nITD
           DO J=1,sNy
            DO I=1,sNx
+            tmpscal1 = HEFFITDpreTH(I,J,IT)
+      &              + d_HEFFbySublim_ITD(I,J,IT)
+      &              + d_HEFFbyOCNonICE_ITD(I,J,IT)
  #ifdef SEAICE_GROWTH_LEGACY
-            tmpscal2 =MAX(-HEFFITD(I,J,K,bi,bj),r_QbyATMmult_cover(I,J,K))
+            tmpscal2 = MAX(-tmpscal1,
+      &                     r_QbyATMmult_cover(I,J,IT))
  #else
-            tmpscal2 =MAX(-HEFFITD(I,J,K,bi,bj),r_QbyATMmult_cover(I,J,K)
+            tmpscal2 = MAX(-tmpscal1,
+      &                     r_QbyATMmult_cover(I,J,IT)
  c         Limit ice growth by potential melt by ocean
-      &     + AREAITDpreTH(I,J,K) * r_QbyOCN(I,J)*areaFracFactor(I,J,K))
+      &              + AREAITDpreTH(I,J,IT) * r_QbyOCN(I,J))
  #endif /* SEAICE_GROWTH_LEGACY */
+            d_HEFFbyATMonOCN_cover_ITD(I,J,IT) = tmpscal2
             d_HEFFbyATMonOCN_cover(I,J) = d_HEFFbyATMonOCN_cover(I,J)
       &                                 + tmpscal2
+            d_HEFFbyATMonOCN_ITD(I,J,IT) = d_HEFFbyATMonOCN_ITD(I,J,IT)
+      &                                 + tmpscal2
             d_HEFFbyATMonOCN(I,J)       = d_HEFFbyATMonOCN(I,J)
       &                                 + tmpscal2
-            r_QbyATM_cover(I,J)         = r_QbyATM_cover(I,J)
+            r_QbyATMmult_cover(I,J,IT)  = r_QbyATMmult_cover(I,J,IT)
       &                                 - tmpscal2
-            HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,bi,bj) + tmpscal2
+           ENDDO
-           ENDDO
+          ENDDO
-          ENDDO
+         ENDDO
-         ENDDO
+ #ifdef ALLOW_SITRACER
- #else
+         DO J=1,sNy
+          DO I=1,sNx
+           SItrHEFF(I,J,bi,bj,3) = SItrHEFF(I,J,bi,bj,2)
+      &                          + d_HEFFbyATMonOCN_cover(I,J)
+          ENDDO
+         ENDDO
+ #endif
+ #else /* SEAICE_ITD */
          DO J=1,sNy
           DO I=1,sNx
-Line 1432 
 c         Limit ice growth by potential
+Line 1642 
 c         Limit ice growth by potential
            tmpscal2 = MAX(-HEFF(I,J,bi,bj),r_QbyATM_cover(I,J))
  #else
            tmpscal2 = MAX(-HEFF(I,J,bi,bj),r_QbyATM_cover(I,J)+
- c         Limit ice growth by potential melt by ocean
+ C         Limit ice growth by potential melt by ocean
       &        AREApreTH(I,J) * r_QbyOCN(I,J))
  #endif /* SEAICE_GROWTH_LEGACY */
            d_HEFFbyATMonOCN_cover(I,J)=tmpscal2
            d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal2
            r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J)-tmpscal2
-          ENDDO
+           HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal2
-         ENDDO
- #endif /* SEAICE_ITD */
-         DO J=1,sNy
-          DO I=1,sNx
-           HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) +d_HEFFbyATMonOCN_cover(I,J)
  #ifdef ALLOW_SITRACER
            SItrHEFF(I,J,bi,bj,3)=HEFF(I,J,bi,bj)
  #endif
           ENDDO
          ENDDO
+ #endif /* SEAICE_ITD */
+ #ifdef SEAICE_DEBUG
+ c ToM<<< debug seaice_growth
+         WRITE(msgBuf,msgBufForm)
+      &   ' SEAICE_GROWTH: Heff increments 4, d_HEFFbyATMonOCN_cover = ',
+ #ifdef SEAICE_ITD
+      &     d_HEFFbyATMonOCN_cover_ITD(1,1,:)
+ #else
+      &     d_HEFFbyATMonOCN_cover(1,1)
+ #endif
+         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
+      &    SQUEEZE_RIGHT , myThid)
+         WRITE(msgBuf,msgBufForm)
+      &    ' SEAICE_GROWTH: Heff increments 4, d_HEFFbyATMonOCN = ',
+ #ifdef SEAICE_ITD
+      &     d_HEFFbyATMonOCN_ITD(1,1,:)
+ #else
+      &     d_HEFFbyATMonOCN(1,1)
+ #endif
+         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
+      &    SQUEEZE_RIGHT , myThid)
+ c ToM>>>
+ #endif /* SEAICE_DEBUG */
- C attribute precip to fresh water or snow stock,
+ C add snow precipitation to HSNOW.
- C depending on atmospheric conditions.
  C =================================================
  #ifdef ALLOW_ATM_TEMP
- #ifdef ALLOW_AUTODIFF_TAMC
+ # ifdef ALLOW_AUTODIFF_TAMC
  CADJ STORE a_QbyATM_cover = comlev1_bibj,key=iicekey,byte=isbyte
  CADJ STORE PRECIP(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
  CADJ STORE AREApreTH = comlev1_bibj,key=iicekey,byte=isbyte
- #endif /* ALLOW_AUTODIFF_TAMC */
+ # endif /* ALLOW_AUTODIFF_TAMC */
-         DO J=1,sNy
+         IF ( snowPrecipFile .NE. ' ' ) THEN
-          DO I=1,sNx
+ C add snowPrecip to HSNOW
+          DO J=1,sNy
+           DO I=1,sNx
+            d_HSNWbyRAIN(I,J) = convertPRECIP2HI * ICE2SNOW *
+      &          snowPrecip(i,j,bi,bj) * AREApreTH(I,J)
+            d_HFRWbyRAIN(I,J) = -convertPRECIP2HI *
+      &          ( PRECIP(I,J,bi,bj) - snowPrecip(I,J,bi,bj) ) *
+      &          AREApreTH(I,J)
+            HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + d_HSNWbyRAIN(I,J)
+           ENDDO
+          ENDDO
+         ELSE
+ C attribute precip to fresh water or snow stock,
+ C depending on atmospheric conditions.
+          DO J=1,sNy
+           DO I=1,sNx
  C possible alternatives to the a_QbyATM_cover criterium
  c          IF (TICE(I,J,bi,bj) .LT. TMIX) THEN
  c          IF (atemp(I,J,bi,bj) .LT. celsius2K) THEN
-           IF ( a_QbyATM_cover(I,J).GE. 0. _d 0 ) THEN
+            IF ( a_QbyATM_cover(I,J).GE. 0. _d 0 ) THEN
  C           add precip as snow
              d_HFRWbyRAIN(I,J)=0. _d 0
              d_HSNWbyRAIN(I,J)=convertPRECIP2HI*ICE2SNOW*
       &            PRECIP(I,J,bi,bj)*AREApreTH(I,J)
-           ELSE
+            ELSE
  C           add precip to the fresh water bucket
              d_HFRWbyRAIN(I,J)=-convertPRECIP2HI*
       &            PRECIP(I,J,bi,bj)*AREApreTH(I,J)
              d_HSNWbyRAIN(I,J)=0. _d 0
-           ENDIF
+            ENDIF
-           HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + d_HSNWbyRAIN(I,J)
           ENDDO
          ENDDO
  #ifdef SEAICE_ITD
-         DO K=1,nITD
+         DO IT=1,nITD
           DO J=1,sNy
            DO I=1,sNx
-            HSNOWITD(I,J,K,bi,bj) = HSNOWITD(I,J,K,bi,bj)
+            d_HSNWbyRAIN_ITD(I,J,IT)
-      &     + d_HSNWbyRAIN(I,J)*areaFracFactor(I,J,K)
+      &     = d_HSNWbyRAIN(I,J)*areaFracFactor(I,J,IT)
            ENDDO
           ENDDO
          ENDDO
+ #else
+         DO J=1,sNy
+          DO I=1,sNx
+           HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj) + d_HSNWbyRAIN(I,J)
+          ENDDO
+         ENDDO
  #endif
  Cgf note: this does not affect air-sea heat flux,
  Cgf since the implied air heat gain to turn
  Cgf rain to snow is not a surface process.
+ C end of IF statement snowPrecipFile:
+         ENDIF
  #endif /* ALLOW_ATM_TEMP */
+ #ifdef SEAICE_DEBUG
+ c ToM<<< debug seaice_growth
+         WRITE(msgBuf,msgBufForm)
+      &    ' SEAICE_GROWTH: Hsnow increments 5, d_HSNWbyRAIN = ',
+ #ifdef SEAICE_ITD
+      &     d_HSNWbyRAIN_ITD(1,1,:)
+ #else
+      &     d_HSNWbyRAIN(1,1)
+ #endif
+         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
+      &    SQUEEZE_RIGHT , myThid)
+ c ToM>>>
+ #endif /* SEAICE_DEBUG */
  C compute snow melt due to heat available from ocean.
  C =================================================================
-Line 1507 
 CADJ STORE r_QbyOCN = comlev1_bibj,key=i
+Line 1769 
 CADJ STORE r_QbyOCN = comlev1_bibj,key=i
  #endif /* ALLOW_AUTODIFF_TAMC */
  #ifdef SEAICE_ITD
-         DO K=1,nITD
+         DO IT=1,nITD
           DO J=1,sNy
            DO I=1,sNx
-            tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW*areaFracFactor(I,J,K),
+            tmpscal4 = HSNWITDpreTH(I,J,IT)
-      &                                        -HSNOW(I,J,bi,bj))
+      &              + d_HSNWbySublim_ITD(I,J,IT)
+      &              + d_HSNWbyATMonSNW_ITD(I,J,IT)
+      &              + d_HSNWbyRAIN_ITD(I,J,IT)
+            tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW*areaFracFactor(I,J,IT),
+      &                  -tmpscal4)
             tmpscal2=MIN(tmpscal1,0. _d 0)
  #ifdef SEAICE_MODIFY_GROWTH_ADJ
  Cgf no additional dependency through snow
             if ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
  #endif
-            HSNOWITD(I,J,K,bi,bj) = HSNOWITD(I,J,K,bi,bj) + tmpscal2
+            d_HSNWbyOCNonSNW_ITD(I,J,IT) = tmpscal2
             d_HSNWbyOCNonSNW(I,J) = d_HSNWbyOCNonSNW(I,J) + tmpscal2
+            r_QbyOCN(I,J)=r_QbyOCN(I,J) - tmpscal2*SNOW2ICE
            ENDDO
           ENDDO
          ENDDO
- #endif
+ #else /* SEAICE_ITD */
          DO J=1,sNy
           DO I=1,sNx
- #ifndef SEAICE_ITD
            tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW, -HSNOW(I,J,bi,bj))
            tmpscal2=MIN(tmpscal1,0. _d 0)
  #ifdef SEAICE_MODIFY_GROWTH_ADJ
-Line 1533 
 Cgf no additional dependency through sno
+Line 1799 
 Cgf no additional dependency through sno
            if ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
  #endif
            d_HSNWbyOCNonSNW(I,J) = tmpscal2
- #endif
            r_QbyOCN(I,J)=r_QbyOCN(I,J)
       &                               -d_HSNWbyOCNonSNW(I,J)*SNOW2ICE
            HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)+d_HSNWbyOCNonSNW(I,J)
           ENDDO
          ENDDO
+ #endif /* SEAICE_ITD */
  #endif /* SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING */
  Cph)
+ #ifdef SEAICE_DEBUG
+ c ToM<<< debug seaice_growth
+         WRITE(msgBuf,msgBufForm)
+      &    ' SEAICE_GROWTH: Hsnow increments 6, d_HSNWbyOCNonSNW = ',
+ #ifdef SEAICE_ITD
+      &     d_HSNWbyOCNonSNW_ITD(1,1,:)
+ #else
+      &     d_HSNWbyOCNonSNW(1,1)
+ #endif
+         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
+      &    SQUEEZE_RIGHT , myThid)
+ c ToM>>>
+ #endif /* SEAICE_DEBUG */
  C gain of new ice over open water
  C ===============================
-Line 1554 
 CADJ STORE a_QSWbyATM_open = comlev1_bib
+Line 1833 
 CADJ STORE a_QSWbyATM_open = comlev1_bib
          DO J=1,sNy
           DO I=1,sNx
- c           Initial ice growth is triggered by open water
+ #ifdef SEAICE_ITD
- c           heat flux overcoming potential melt by ocean
+ C         HEFF will be updated at the end of PART 3,
+ C         hence sum of tendencies so far is needed
+           tmpscal4 = HEFFpreTH(I,J)
+      &             + d_HEFFbySublim(I,J)
+      &             + d_HEFFbyOCNonICE(I,J)
+      &             + d_HEFFbyATMonOCN(I,J)
+ #else
+ C         HEFF is updated step by step throughout seaice_growth
+           tmpscal4 = HEFF(I,J,bi,bj)
+ #endif
+ C           Initial ice growth is triggered by open water
+ C           heat flux overcoming potential melt by ocean
              tmpscal1=r_QbyATM_open(I,J)+r_QbyOCN(i,j) *
       &            (1.0 _d 0 - AREApreTH(I,J))
- c           Penetrative shortwave flux beyond first layer
+ C           Penetrative shortwave flux beyond first layer
- c           that is therefore not available to ice growth/melt
+ C           that is therefore not available to ice growth/melt
              tmpscal2=SWFracB * a_QSWbyATM_open(I,J)
  C           impose -HEFF as the maxmum melting if SEAICE_doOpenWaterMelt
  C           or 0. otherwise (no melting if not SEAICE_doOpenWaterMelt)
              tmpscal3=facOpenGrow*MAX(tmpscal1-tmpscal2,
-      &           -HEFF(I,J,bi,bj)*facOpenMelt)*HEFFM(I,J,bi,bj)
+      &           -tmpscal4*facOpenMelt)*HEFFM(I,J,bi,bj)
+ #ifdef SEAICE_ITD
+ C         ice growth in open water adds to first category
+           d_HEFFbyATMonOCN_open_ITD(I,J,1)=tmpscal3
+           d_HEFFbyATMonOCN_ITD(I,J,1)     =d_HEFFbyATMonOCN_ITD(I,J,1)
+      &                                    +tmpscal3
+ #endif
            d_HEFFbyATMonOCN_open(I,J)=tmpscal3
            d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal3
            r_QbyATM_open(I,J)=r_QbyATM_open(I,J)-tmpscal3
- #ifdef SEAICE_ITD
- C         determine thickness of new ice
- C         considering the entire open water area to refreeze
-           tmpscal4 = tmpscal3/(ONE-AREApreTH(I,J))
- C         then add new ice volume to appropriate thickness category
-           DO K=1,nITD
-            IF (tmpscal4.LT.Hlimit(K)) THEN
-             HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,bi,bj) + tmpscal3
-             AREAITD(I,J,K,bi,bj) = AREAITD(I,J,K,bi,bj)
-      &                           + ONE-AREApreTH(I,J)
-            ENDIF
-           ENDDO
- C         in this case no open water is left after this step
-           AREA(I,J,bi,bj) = ONE
- #endif
            HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal3
           ENDDO
          ENDDO
-Line 1590 
 C         in this case no open water is
+Line 1871 
 C         in this case no open water is
  #ifdef ALLOW_SITRACER
          DO J=1,sNy
           DO I=1,sNx
- c needs to be here to allow use also with LEGACY branch
+ C needs to be here to allow use also with LEGACY branch
+ #ifdef SEAICE_ITD
+           SItrHEFF(I,J,bi,bj,4)=SItrHEFF(I,J,bi,bj,3)
+      &                         +d_HEFFbyATMonOCN_open(I,J)
+ #else
            SItrHEFF(I,J,bi,bj,4)=HEFF(I,J,bi,bj)
+ #endif
           ENDDO
          ENDDO
  #endif /* ALLOW_SITRACER */
+ #ifdef SEAICE_DEBUG
+ c ToM<<< debug seaice_growth
+         WRITE(msgBuf,msgBufForm)
+      &    ' SEAICE_GROWTH: Heff increments 7, d_HEFFbyATMonOCN_open = ',
+ #ifdef SEAICE_ITD
+      &     d_HEFFbyATMonOCN_open_ITD(1,1,:)
+ #else
+      &     d_HEFFbyATMonOCN_open(1,1)
+ #endif
+         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
+      &    SQUEEZE_RIGHT , myThid)
+         WRITE(msgBuf,msgBufForm)
+      &    ' SEAICE_GROWTH: Heff increments 7, d_HEFFbyATMonOCN = ',
+ #ifdef SEAICE_ITD
+      &     d_HEFFbyATMonOCN_ITD(1,1,:)
+ #else
+      &     d_HEFFbyATMonOCN(1,1)
+ #endif
+         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
+      &    SQUEEZE_RIGHT , myThid)
+ c ToM>>>
+ #endif /* SEAICE_DEBUG */
  C convert snow to ice if submerged.
  C =================================
-Line 1607 
 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
+Line 1915 
 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
  #endif /* ALLOW_AUTODIFF_TAMC */
          IF ( SEAICEuseFlooding ) THEN
  #ifdef SEAICE_ITD
-          DO K=1,nITD
+          DO IT=1,nITD
            DO J=1,sNy
             DO I=1,sNx
-             tmpscal0 = (HSNOWITD(I,J,K,bi,bj)*SEAICE_rhoSnow
+             tmpscal3 = HEFFITDpreTH(I,J,IT)
-      &               +HEFFITD(I,J,K,bi,bj)*SEAICE_rhoIce)*recip_rhoConst
+      &               + d_HEFFbySublim_ITD(I,J,IT)
-             tmpscal1 = MAX( 0. _d 0, tmpscal0 - HEFFITD(I,J,K,bi,bj))
+      &               + d_HEFFbyOCNonICE_ITD(I,J,IT)
-             d_HEFFbyFLOODING(I,J) = d_HEFFbyFLOODING(I,J) + tmpscal1
+      &               + d_HEFFbyATMonOCN_ITD(I,J,IT)
-             HEFFITD(I,J,K,bi,bj)  = HEFFITD(I,J,K,bi,bj)  + tmpscal1
+             tmpscal4 = HSNWITDpreTH(I,J,IT)
-             HSNOWITD(I,J,K,bi,bj) = HSNOWITD(I,J,K,bi,bj) - tmpscal1
+      &               + d_HSNWbySublim_ITD(I,J,IT)
-      &                            * ICE2SNOW
+      &               + d_HSNWbyATMonSNW_ITD(I,J,IT)
-            ENDDO
+      &               + d_HSNWbyRAIN_ITD(I,J,IT)
-           ENDDO
+             tmpscal0 = (tmpscal4*SEAICE_rhoSnow
-          ENDDO
+      &               +  tmpscal3*SEAICE_rhoIce)
+      &               * recip_rhoConst
+             tmpscal1 = MAX( 0. _d 0, tmpscal0 - tmpscal3)
+             d_HEFFbyFLOODING_ITD(I,J,IT) = tmpscal1
+             d_HEFFbyFLOODING(I,J) = d_HEFFbyFLOODING(I,J)  + tmpscal1
+            ENDDO
+           ENDDO
+          ENDDO
  #else
           DO J=1,sNy
            DO I=1,sNx
-Line 1627 
 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
+Line 1942 
 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
       &              +HEFF(I,J,bi,bj)*SEAICE_rhoIce)*recip_rhoConst
             tmpscal1 = MAX( 0. _d 0, tmpscal0 - HEFF(I,J,bi,bj))
             d_HEFFbyFLOODING(I,J)=tmpscal1
-           ENDDO
-          ENDDO
- #endif
-          DO J=1,sNy
-           DO I=1,sNx
             HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj)+d_HEFFbyFLOODING(I,J)
             HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)-
       &                           d_HEFFbyFLOODING(I,J)*ICE2SNOW
            ENDDO
           ENDDO
+ #endif
          ENDIF
  #endif /* SEAICE_GROWTH_LEGACY */
+ #ifdef SEAICE_DEBUG
+ c ToM<<< debug seaice_growth
+         WRITE(msgBuf,msgBufForm)
+      &    ' SEAICE_GROWTH: Heff increments 8, d_HEFFbyFLOODING = ',
+ #ifdef SEAICE_ITD
+      &     d_HEFFbyFLOODING_ITD(1,1,:)
+ #else
+      &     d_HEFFbyFLOODING(1,1)
+ #endif
+         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
+      &    SQUEEZE_RIGHT , myThid)
+ c ToM>>>
+ #endif /* SEAICE_DEBUG */
+ #ifdef SEAICE_ITD
+ C apply ice and snow thickness changes
+ C =================================================================
+          DO IT=1,nITD
+           DO J=1,sNy
+            DO I=1,sNx
+             HEFFITD(I,J,IT,bi,bj) = HEFFITD(I,J,IT,bi,bj)
+      &                            + d_HEFFbySublim_ITD(I,J,IT)
+      &                            + d_HEFFbyOCNonICE_ITD(I,J,IT)
+      &                            + d_HEFFbyATMonOCN_ITD(I,J,IT)
+      &                            + d_HEFFbyFLOODING_ITD(I,J,IT)
+             HSNOWITD(I,J,IT,bi,bj) = HSNOWITD(I,J,IT,bi,bj)
+      &                            + d_HSNWbySublim_ITD(I,J,IT)
+      &                            + d_HSNWbyATMonSNW_ITD(I,J,IT)
+      &                            + d_HSNWbyRAIN_ITD(I,J,IT)
+      &                            + d_HSNWbyOCNonSNW_ITD(I,J,IT)
+      &                            - d_HEFFbyFLOODING_ITD(I,J,IT)
+      &                            * ICE2SNOW
+            ENDDO
+           ENDDO
+          ENDDO
+ #endif
+ c ToM<<< debug seaice_growth
+         WRITE(msgBuf,msgBufForm)
+      &    ' SEAICE_GROWTH: Heff increments 9, HEFF = ',
+ #ifdef SEAICE_ITD
+      &     HEFFITD(1,1,:,bi,bj)
+ #else
+      &     HEFF(1,1,bi,bj)
+ #endif
+         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
+      &    SQUEEZE_RIGHT , myThid)
+         WRITE(msgBuf,msgBufForm)
+      &    ' SEAICE_GROWTH: Area increments 9, AREA = ',
+ #ifdef SEAICE_ITD
+      &     AREAITD(1,1,:,bi,bj)
+ #else
+      &     AREA(1,1,bi,bj)
+ #endif
+         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
+      &    SQUEEZE_RIGHT , myThid)
+ c ToM>>>
  C ===================================================================
  C ==========PART 4: determine ice cover fraction increments=========-
-Line 1665 
 CADJ STORE AREA(:,:,bi,bj) = comlev1_bib
+Line 2031 
 CADJ STORE AREA(:,:,bi,bj) = comlev1_bib
  #endif /* ALLOW_AUTODIFF_TAMC */
  #ifdef SEAICE_ITD
- C       replaces Hibler '79 scheme and lead closing parameter
+ C--     account for lateral ice growth and melt only in thinnest category
- C       because ITD accounts explicitly for lead openings and
+ C--     use HEFF, ARE, HSNOW, etc. temporarily for 1st category
- C       different melt rates due to varying ice thickness
+ C       (this way we can use same code for ITD and non-ITD case)
- C
+         DO J=1,sNy
- C       only consider ice area loss due to total ice thickness loss
+          DO I=1,sNx
- C       ice area gain due to freezing of open water as handled above
+           HEFF(I,J,bi,bj)=HEFFITD(I,J,1,bi,bj)
- C       under "gain of new ice over open water"
+           AREA(I,J,bi,bj)=AREAITD(I,J,1,bi,bj)
- C
+           HSNOW(I,J,bi,bj)=HSNOWITD(I,J,1,bi,bj)
- C       does not account for lateral melt of ice floes
+           HEFFpreTH(I,J)=HEFFITDpreTH(I,J,1)
- C
+           AREApreTH(I,J)=AREAITDpreTH(I,J,1)
-         DO K=1,nITD
+           recip_heffActual(I,J)=recip_heffActualMult(I,J,1)
-          DO J=1,sNy
+          ENDDO
-           DO I=1,sNx
+         ENDDO
-            IF (HEFFITD(I,J,K,bi,bj).LE.ZERO) THEN
+ C       all other categories only experience basal growth or melt,
-             AREAITD(I,J,K,bi,bj)=ZERO
+ C       i.e. change sin AREA only occur when all ice in a category is melted
-            ENDIF
+         IF (nITD .gt. 1) THEN
-           ENDDO
+          DO IT=2,nITD
-          ENDDO
+           DO J=1,sNy
-         ENDDO
+            DO I=1,sNx
- C       update total AREA, HEFF, HSNOW
+             IF (HEFFITD(I,J,IT,bi,bj).LE.ZERO) THEN
-         CALL SEAICE_ITD_SUM(myTime,myIter,myThid)
+              AREAITD(I,J,IT,bi,bj)=ZERO
- #else
+             ENDIF
+            ENDDO
+           ENDDO
+          ENDDO
+         ENDIF
+ #endif
          DO J=1,sNy
           DO I=1,sNx
-Line 1756 
 C apply tendency
+Line 2127 
 C apply tendency
  #endif /* ALLOW_DIAGNOSTICS */
           ENDDO
          ENDDO
- #endif /* SEAICE_ITD */
+ #ifdef SEAICE_ITD
+ C       transfer 1st category values back into ITD variables
+         DO J=1,sNy
+          DO I=1,sNx
+           HEFFITD(I,J,1,bi,bj)=HEFF(I,J,bi,bj)
+           AREAITD(I,J,1,bi,bj)=AREA(I,J,bi,bj)
+           HSNOWITD(I,J,1,bi,bj)=HSNOW(I,J,bi,bj)
+          ENDDO
+         ENDDO
+ #endif
  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
  Cgf 'bulk' linearization of area=f(HEFF)
          IF ( SEAICEadjMODE.GE.1 ) THEN
  #ifdef SEAICE_ITD
-          DO K=1,nITD
+          DO IT=1,nITD
            DO J=1,sNy
             DO I=1,sNx
-             AREAITD(I,J,K,bi,bj) = AREAITDpreTH(I,J,K) + 0.1 _d 0 *
+             AREAITD(I,J,IT,bi,bj) = AREAITDpreTH(I,J,IT) + 0.1 _d 0 *
-      &               ( HEFFITD(I,J,K,bi,bj) - HEFFITDpreTH(I,J,K) )
+      &               ( HEFFITD(I,J,IT,bi,bj) - HEFFITDpreTH(I,J,IT) )
             ENDDO
            ENDDO
           ENDDO
- C        update total AREA, HEFF, HSNOW
-          CALL SEAICE_ITD_SUM(myTime,myIter,myThid)
  #else
           DO J=1,sNy
            DO I=1,sNx
-Line 1783 
 C            AREA(I,J,bi,bj) = 0.1 _d 0
+Line 2161 
 C            AREA(I,J,bi,bj) = 0.1 _d 0
  #endif
          ENDIF
  #endif
+ #ifdef SEAICE_ITD
+ C check categories for consistency with limits after growth/melt
+         CALL SEAICE_ITD_REDIST(bi, bj, myTime,myIter,myThid)
+ C finally update total AREA, HEFF, HSNOW
+         CALL SEAICE_ITD_SUM(bi, bj, myTime,myIter,myThid)
+ #endif
  C ===================================================================
  C =============PART 5: determine ice salinity increments=============
-Line 1805 
 CADJ &                       key = iicek
+Line 2189 
 CADJ &                       key = iicek
            tmpscal1 = d_HEFFbyNEG(I,J) + d_HEFFbyOCNonICE(I,J) +
       &               d_HEFFbyATMonOCN(I,J) + d_HEFFbyFLOODING(I,J)
       &             + d_HEFFbySublim(I,J)
- #ifdef SEAICE_ALLOW_AREA_RELAXATION
+ #ifdef EXF_ALLOW_SEAICE_RELAX
-                    + d_HEFFbyRLX(I,J)
+      &             + d_HEFFbyRLX(I,J)
  #endif
            tmpscal2 = tmpscal1 * SEAICE_salt0 * HEFFM(I,J,bi,bj)
       &            * recip_deltaTtherm * SEAICE_rhoIce
-Line 1900 
 C set HSALT = 0 if HEFF = 0 and compute
+Line 2284 
 C set HSALT = 0 if HEFF = 0 and compute
          ENDDO
  #endif /* SEAICE_VARIABLE_SALINITY */
  C =======================================================================
  C ==LEGACY PART 6 (LEGACY) treat pathological cases, then do flooding ===
  C =======================================================================
-Line 1983 
 C =================================
+Line 2366 
 C =================================
  #ifdef ALLOW_SITRACER
          DO J=1,sNy
           DO I=1,sNx
- c needs to be here to allow use also with LEGACY branch
+ C needs to be here to allow use also with LEGACY branch
            SItrHEFF(I,J,bi,bj,5)=HEFF(I,J,bi,bj)
           ENDDO
          ENDDO
-Line 1997 
 C compute net heat flux leaving/entering
+Line 2380 
 C compute net heat flux leaving/entering
  C accounting for the part used in melt/freeze processes
  C =====================================================
+ #ifdef SEAICE_ITD
+ C compute total of "mult" fluxes for ocean forcing
+         DO J=1,sNy
+          DO I=1,sNx
+           a_QbyATM_cover(I,J)   = 0.0 _d 0
+           r_QbyATM_cover(I,J)   = 0.0 _d 0
+           a_QSWbyATM_cover(I,J) = 0.0 _d 0
+           r_FWbySublim(I,J)     = 0.0 _d 0
+          ENDDO
+         ENDDO
+         DO IT=1,nITD
+          DO J=1,sNy
+           DO I=1,sNx
+ cToM if fluxes in W/m^2 then
+ c           a_QbyATM_cover(I,J)=a_QbyATM_cover(I,J)
+ c     &      + a_QbyATMmult_cover(I,J,IT) * areaFracFactor(I,J,IT)
+ c           r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J)
+ c     &      + r_QbyATMmult_cover(I,J,IT) * areaFracFactor(I,J,IT)
+ c           a_QSWbyATM_cover(I,J)=a_QSWbyATM_cover(I,J)
+ c     &      + a_QSWbyATMmult_cover(I,J,IT) * areaFracFactor(I,J,IT)
+ c           r_FWbySublim(I,J)=r_FWbySublim(I,J)
+ c     &      + r_FWbySublimMult(I,J,IT) * areaFracFactor(I,J,IT)
+ cToM if fluxes in effective ice meters, i.e. ice volume per area, then
+            a_QbyATM_cover(I,J)=a_QbyATM_cover(I,J)
+      &      + a_QbyATMmult_cover(I,J,IT)
+            r_QbyATM_cover(I,J)=r_QbyATM_cover(I,J)
+      &      + r_QbyATMmult_cover(I,J,IT)
+            a_QSWbyATM_cover(I,J)=a_QSWbyATM_cover(I,J)
+      &      + a_QSWbyATMmult_cover(I,J,IT)
+            r_FWbySublim(I,J)=r_FWbySublim(I,J)
+      &      + r_FWbySublimMult(I,J,IT)
+           ENDDO
+          ENDDO
+         ENDDO
+ #endif
  #ifdef ALLOW_AUTODIFF_TAMC
  CADJ STORE d_hsnwbyneg = comlev1_bibj,key=iicekey,byte=isbyte
  CADJ STORE d_hsnwbyocnonsnw = comlev1_bibj,key=iicekey,byte=isbyte
-Line 2010 
 C in principle a_QSWbyATM_cover should a
+Line 2429 
 C in principle a_QSWbyATM_cover should a
  C for backward compatibility it is left out of the LEGACY branch
       &         +   a_QSWbyATM_cover(I,J)
  #endif /* SEAICE_GROWTH_LEGACY */
-      &         - ( d_HEFFbyOCNonICE(I,J) +
+      &         - ( d_HEFFbyOCNonICE(I,J)
-      &             d_HSNWbyOCNonSNW(I,J)*SNOW2ICE +
+      &           + d_HSNWbyOCNonSNW(I,J)*SNOW2ICE
-      &             d_HEFFbyNEG(I,J) +
+      &           + d_HEFFbyNEG(I,J)
- #ifdef SEAICE_ALLOW_AREA_RELAXATION
+ #ifdef EXF_ALLOW_SEAICE_RELAX
-      &             d_HEFFbyRLX(I,J) +
+      &           + d_HEFFbyRLX(I,J)
  #endif
-      &             d_HSNWbyNEG(I,J)*SNOW2ICE )
+      &           + d_HSNWbyNEG(I,J)*SNOW2ICE
-      &         * maskC(I,J,kSurface,bi,bj)
+      &           - convertPRECIP2HI *
+      &             snowPrecip(i,j,bi,bj) * (ONE-AREApreTH(I,J))
+      &           ) * maskC(I,J,kSurface,bi,bj)
+          ENDDO
+         ENDDO
+         DO J=1,sNy
+          DO I=1,sNx
            QSW(I,J,bi,bj)  = a_QSWbyATM_cover(I,J) + a_QSWbyATM_open(I,J)
           ENDDO
          ENDDO
-Line 2045 
 CADJ STORE d_HSNWbyATMonSNW = comlev1_bi
+Line 2470 
 CADJ STORE d_HSNWbyATMonSNW = comlev1_bi
  CADJ STORE theta(:,:,kSurface,bi,bj) = comlev1_bibj,
  CADJ &                       key = iicekey, byte = isbyte
  # endif /* ALLOW_AUTODIFF_TAMC */
-       IF ( SEAICEheatConsFix ) THEN
+ cgf Unlike for evap and precip, the temperature of gained/lost
- c Unlike for evap and precip, the temperature of gained/lost
+ C ocean liquid water due to melt/freeze of solid water cannot be chosen
- c ocean liquid water due to melt/freeze of solid water cannot be chosen
+ C arbitrarily to be e.g. the ocean SST. Indeed the present seaice model
- c to be e.g. the ocean SST. It must be done at 0degC. The fix below anticipates
+ C implies a constant ice temperature of 0degC. If melt/freeze water is exchanged
- c on external_forcing_surf.F and applies the correction to QNET.
+ C at a different temperature, it leads to a loss of conservation in the
-       IF ((convertFW2Salt.EQ.-1.).OR.(temp_EvPrRn.EQ.UNSET_RL)) THEN
+ C ocean+ice system. While this is mostly a serious issue in the
- c I leave alone the exotic case when onvertFW2Salt.NE.-1 and temp_EvPrRn.NE.UNSET_RL and
+ C real fresh water + non linear free surface framework, a mismatch
- c the small error of the synchronous time stepping case (see external_forcing_surf.F).
+ C between ice and ocean boundary condition can result in all cases.
+ C Below we therefore anticipate on external_forcing_surf.F
+ C to diagnoze and/or apply the correction to QNET.
          DO J=1,sNy
           DO I=1,sNx
- #ifdef ALLOW_DIAGNOSTICS
+ C ocean water going to ice/snow, in precip units
- c store unaltered QNET for diagnostic purposes
+            tmpscal3=rhoConstFresh*maskC(I,J,kSurface,bi,bj)*(
-            DIAGarrayA(I,J)=QNET(I,J,bi,bj)
- #endif
- c compute the ocean water going to ice/snow, in precip units
-            tmpscal3=rhoConstFresh*maskC(I,J,kSurface,bi,bj)*
       &       ( d_HSNWbyATMonSNW(I,J)*SNOW2ICE
       &       + d_HSNWbyOCNonSNW(I,J)*SNOW2ICE
       &       + d_HEFFbyOCNonICE(I,J) + d_HEFFbyATMonOCN(I,J)
       &       + d_HEFFbyNEG(I,J) + d_HSNWbyNEG(I,J)*SNOW2ICE )
       &       * convertHI2PRECIP
- c factor in the heat content that external_forcing_surf.F
+      &       - snowPrecip(i,j,bi,bj) * (ONE-AREApreTH(I,J)) )
- c will associate with EMPMR, and remove it from QNET, so that
+ C factor in the heat content as done in external_forcing_surf.F
- c melt/freez water is in effect consistently gained/lost at 0degC
+       IF ( (temp_EvPrRn.NE.UNSET_RL).AND.useRealFreshWaterFlux
-            IF (temp_EvPrRn.NE.UNSET_RL) THEN
+      &     .AND.(nonlinFreeSurf.NE.0) ) THEN
-              QNET(I,J,bi,bj)=QNET(I,J,bi,bj) - tmpscal3*
+              tmpscal1 = - tmpscal3*
       &         HeatCapacity_Cp * temp_EvPrRn
-            ELSE
+       ELSEIF ( (temp_EvPrRn.EQ.UNSET_RL).AND.useRealFreshWaterFlux
-              QNET(I,J,bi,bj)=QNET(I,J,bi,bj) - tmpscal3*
+      &         .AND.(nonlinFreeSurf.NE.0) ) THEN
+              tmpscal1 = - tmpscal3*
       &         HeatCapacity_Cp * theta(I,J,kSurface,bi,bj)
-            ENDIF
+       ELSEIF ( (temp_EvPrRn.NE.UNSET_RL) ) THEN
+            tmpscal1 = - tmpscal3*HeatCapacity_Cp*
+      &       ( temp_EvPrRn - theta(I,J,kSurface,bi,bj) )
+       ELSEIF ( (temp_EvPrRn.EQ.UNSET_RL) ) THEN
+            tmpscal1 = ZERO
+       ENDIF
  #ifdef ALLOW_DIAGNOSTICS
- c back out the eventual TFLUX adjustement and fill diag
+ C in all cases, diagnoze the boundary condition mismatch to SIaaflux
-            DIAGarrayA(I,J)=QNET(I,J,bi,bj)-DIAGarrayA(I,J)
+            DIAGarrayA(I,J)=tmpscal1
  #endif
+ C remove the mismatch when real fresh water is exchanged (at 0degC here)
+       IF ( useRealFreshWaterFlux.AND.(nonlinFreeSurf.GT.0).AND.
+      &   SEAICEheatConsFix ) QNET(I,J,bi,bj)=QNET(I,J,bi,bj)+tmpscal1
           ENDDO
          ENDDO
-       ENDIF
  #ifdef ALLOW_DIAGNOSTICS
          CALL DIAGNOSTICS_FILL(DIAGarrayA,
       &      'SIaaflux',0,1,3,bi,bj,myThid)
  #endif
-       ENDIF
  #endif /* ndef SEAICE_DISABLE_HEATCONSFIX */
+ C compute the net heat flux, incl. adv. by water, entering ocean+ice
+ C ===================================================================
+          DO J=1,sNy
+           DO I=1,sNx
+ cgf 1) SIatmQnt (analogous to qnet; excl. adv. by water exch.)
+ CML If I consider the atmosphere above the ice, the surface flux
+ CML which is relevant for the air temperature dT/dt Eq
+ CML accounts for sensible and radiation (with different treatment
+ CML according to wave-length) fluxes but not for "latent heat flux",
+ CML since it does not contribute to heating the air.
+ CML So this diagnostic is only good for heat budget calculations within
+ CML the ice-ocean system.
+            SIatmQnt(I,J,bi,bj) =
+      &            maskC(I,J,kSurface,bi,bj)*convertHI2Q*(
+ #ifndef SEAICE_GROWTH_LEGACY
+      &            a_QSWbyATM_cover(I,J) +
+ #endif /* SEAICE_GROWTH_LEGACY */
+      &            a_QbyATM_cover(I,J) + a_QbyATM_open(I,J) )
+ cgf 2) SItflux (analogous to tflux; includes advection by water
+ C             exchanged between atmosphere and ocean+ice)
+ C solid water going to atm, in precip units
+            tmpscal1 = rhoConstFresh*maskC(I,J,kSurface,bi,bj)
+      &       * convertHI2PRECIP * ( - d_HSNWbyRAIN(I,J)*SNOW2ICE
+      &       + a_FWbySublim(I,J) - r_FWbySublim(I,J) )
+ C liquid water going to atm, in precip units
+            tmpscal2=rhoConstFresh*maskC(I,J,kSurface,bi,bj)*
+      &       ( ( EVAP(I,J,bi,bj)-PRECIP(I,J,bi,bj) )
+      &         * ( ONE - AREApreTH(I,J) )
+ #ifdef ALLOW_RUNOFF
+      &         - RUNOFF(I,J,bi,bj)
+ #endif /* ALLOW_RUNOFF */
+      &         + ( d_HFRWbyRAIN(I,J) + r_FWbySublim(I,J) )
+      &         *convertHI2PRECIP )
+ C In real fresh water flux + nonlinFS, we factor in the advected specific
+ C energy (referenced to 0 for 0deC liquid water). In virtual salt flux or
+ C linFS, rain/evap get a special treatment (see external_forcing_surf.F).
+            tmpscal1= - tmpscal1*
+      &       ( -SEAICE_lhFusion + HeatCapacity_Cp * ZERO )
+       IF ( (temp_EvPrRn.NE.UNSET_RL).AND.useRealFreshWaterFlux
+      &     .AND.(nonlinFreeSurf.NE.0) ) THEN
+            tmpscal2= - tmpscal2*
+      &       ( ZERO + HeatCapacity_Cp * temp_EvPrRn )
+       ELSEIF ( (temp_EvPrRn.EQ.UNSET_RL).AND.useRealFreshWaterFlux
+      &        .AND.(nonlinFreeSurf.NE.0) ) THEN
+            tmpscal2= - tmpscal2*
+      &       ( ZERO + HeatCapacity_Cp * theta(I,J,kSurface,bi,bj) )
+       ELSEIF ( (temp_EvPrRn.NE.UNSET_RL) ) THEN
+            tmpscal2= - tmpscal2*HeatCapacity_Cp*
+      &       ( temp_EvPrRn - theta(I,J,kSurface,bi,bj) )
+       ELSEIF ( (temp_EvPrRn.EQ.UNSET_RL) ) THEN
+            tmpscal2= ZERO
+       ENDIF
+            SItflux(I,J,bi,bj)=
+      &            SIatmQnt(I,J,bi,bj)-tmpscal1-tmpscal2
+           ENDDO
+          ENDDO
  C compute net fresh water flux leaving/entering
  C the ocean, accounting for fresh/salt water stocks.
  C ==================================================
-Line 2103 
 C ======================================
+Line 2590 
 C ======================================
       &             +d_HEFFbyOCNonICE(I,J)
       &             +d_HEFFbyATMonOCN(I,J)
       &             +d_HEFFbyNEG(I,J)
- #ifdef SEAICE_ALLOW_AREA_RELAXATION
+ #ifdef EXF_ALLOW_SEAICE_RELAX
       &             +d_HEFFbyRLX(I,J)
  #endif
       &             +d_HSNWbyNEG(I,J)*SNOW2ICE
-Line 2117 
 C     If r_FWbySublim>0, then it is evap
+Line 2604 
 C     If r_FWbySublim>0, then it is evap
  #endif /* ALLOW_RUNOFF */
       &         + tmpscal1*convertHI2PRECIP
       &         )*rhoConstFresh
+ c and the flux leaving/entering the ocean+ice
+            SIatmFW(I,J,bi,bj) = maskC(I,J,kSurface,bi,bj)*(
+      &          EVAP(I,J,bi,bj)*( ONE - AREApreTH(I,J) )
+      &          - PRECIP(I,J,bi,bj)
+ #ifdef ALLOW_RUNOFF
+      &          - RUNOFF(I,J,bi,bj)
+ #endif /* ALLOW_RUNOFF */
+      &           )*rhoConstFresh
+      &     + a_FWbySublim(I,J) * SEAICE_rhoIce * recip_deltaTtherm
           ENDDO
          ENDDO
-Line 2189 
 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
+Line 2686 
 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
           ENDDO
          ENDIF
+ #ifdef ALLOW_BALANCE_FLUXES
+ C Compute tile integrals of heat/fresh water fluxes to/from atm.
+ C ==============================================================
+        FWFsiTile(bi,bj) = 0. _d 0
+        IF ( balanceEmPmR ) THEN
+         DO j=1,sNy
+          DO i=1,sNx
+           FWFsiTile(bi,bj) =
+      &      FWFsiTile(bi,bj) + SIatmFW(i,j,bi,bj)
+      &      * rA(i,j,bi,bj) * maskInC(i,j,bi,bj)
+          ENDDO
+         ENDDO
+        ENDIF
+ c to translate global mean FWF adjustements (see below) we may need :
+        FWF2HFsiTile(bi,bj) = 0. _d 0
+        IF ( balanceEmPmR.AND.(temp_EvPrRn.EQ.UNSET_RL) ) THEN
+         DO j=1,sNy
+          DO i=1,sNx
+           FWF2HFsiTile(bi,bj) = FWF2HFsiTile(bi,bj) +
+      &      HeatCapacity_Cp * theta(I,J,kSurface,bi,bj)
+      &      * rA(i,j,bi,bj) * maskInC(i,j,bi,bj)
+          ENDDO
+         ENDDO
+        ENDIF
+        HFsiTile(bi,bj) = 0. _d 0
+        IF ( balanceQnet ) THEN
+         DO j=1,sNy
+          DO i=1,sNx
+           HFsiTile(bi,bj) =
+      &      HFsiTile(bi,bj) + SItflux(i,j,bi,bj)
+      &      * rA(i,j,bi,bj) * maskInC(i,j,bi,bj)
+          ENDDO
+         ENDDO
+        ENDIF
+ #endif
  C ===================================================================
  C ======================PART 8: diagnostics==========================
  C ===================================================================
-Line 2239 
 C three that actually need intermediate
+Line 2772 
 C three that actually need intermediate
  #ifdef ALLOW_ATM_TEMP
           DO J=1,sNy
            DO I=1,sNx
- CML If I consider the atmosphere above the ice, the surface flux
- CML which is relevant for the air temperature dT/dt Eq
- CML accounts for sensible and radiation (with different treatment
- CML according to wave-length) fluxes but not for "latent heat flux",
- CML since it does not contribute to heating the air.
- CML So this diagnostic is only good for heat budget calculations within
- CML the ice-ocean system.
-            DIAGarrayA(I,J) = maskC(I,J,kSurface,bi,bj)*convertHI2Q*(
- #ifndef SEAICE_GROWTH_LEGACY
-      &            a_QSWbyATM_cover(I,J) +
- #endif /* SEAICE_GROWTH_LEGACY */
-      &            a_QbyATM_cover(I,J) + a_QbyATM_open(I,J) )
- C
             DIAGarrayB(I,J) = maskC(I,J,kSurface,bi,bj) *
       &       a_FWbySublim(I,J) * SEAICE_rhoIce * recip_deltaTtherm
- C
-            DIAGarrayC(I,J) = maskC(I,J,kSurface,bi,bj)*(
-      &          PRECIP(I,J,bi,bj)
-      &          - EVAP(I,J,bi,bj)*( ONE - AREApreTH(I,J) )
- #ifdef ALLOW_RUNOFF
-      &          + RUNOFF(I,J,bi,bj)
- #endif /* ALLOW_RUNOFF */
-      &           )*rhoConstFresh
-      &     - a_FWbySublim(I,J) * SEAICE_rhoIce * recip_deltaTtherm
            ENDDO
           ENDDO
-          CALL DIAGNOSTICS_FILL(DIAGarrayA,
-      &      'SIatmQnt',0,1,3,bi,bj,myThid)
           CALL DIAGNOSTICS_FILL(DIAGarrayB,
       &      'SIfwSubl',0,1,3,bi,bj,myThid)
-          CALL DIAGNOSTICS_FILL(DIAGarrayC,
-      &      'SIatmFW ',0,1,3,bi,bj,myThid)
  C
           DO J=1,sNy
            DO I=1,sNx
-Line 2278 
 C the actual Freshwater flux of sublimat
+Line 2785 
 C the actual Freshwater flux of sublimat
             DIAGarrayA(I,J) = maskC(I,J,kSurface,bi,bj)
       &       * (a_FWbySublim(I,J)-r_FWbySublim(I,J))
       &       * SEAICE_rhoIce * recip_deltaTtherm
- c the residual Freshwater flux of sublimated ice
+ C the residual Freshwater flux of sublimated ice
             DIAGarrayC(I,J) = maskC(I,J,kSurface,bi,bj)
       &       * r_FWbySublim(I,J)
       &       * SEAICE_rhoIce * recip_deltaTtherm
-Line 2295 
 C the latent heat flux
+Line 2802 
 C the latent heat flux
           CALL DIAGNOSTICS_FILL(DIAGarrayA,'SIacSubl',0,1,3,bi,bj,myThid)
           CALL DIAGNOSTICS_FILL(DIAGarrayC,'SIrsSubl',0,1,3,bi,bj,myThid)
           CALL DIAGNOSTICS_FILL(DIAGarrayB,'SIhl    ',0,1,3,bi,bj,myThid)
-          DO J=1,sNy
-           DO I=1,sNx
- c compute ice/snow water going to atm, in precip units
-            tmpscal1 = rhoConstFresh*maskC(I,J,kSurface,bi,bj)
-      &       * convertHI2PRECIP * ( - d_HSNWbyRAIN(I,J)*SNOW2ICE
-      &       + a_FWbySublim(I,J) - r_FWbySublim(I,J) )
- c compute ocean water going to atm, in precip units
-            tmpscal2=rhoConstFresh*maskC(I,J,kSurface,bi,bj)*
-      &       ( ( EVAP(I,J,bi,bj)-PRECIP(I,J,bi,bj) )
-      &         * ( ONE - AREApreTH(I,J) )
- #ifdef ALLOW_RUNOFF
-      &         - RUNOFF(I,J,bi,bj)
- #endif /* ALLOW_RUNOFF */
-      &         + ( d_HFRWbyRAIN(I,J) + r_FWbySublim(I,J) )
-      &         *convertHI2PRECIP )
- c factor in the advected specific energy (referenced to 0 for 0deC liquid water)
-            tmpscal1= - tmpscal1*
-      &       ( -SEAICE_lhFusion + HeatCapacity_Cp * ZERO )
-            IF (temp_EvPrRn.NE.UNSET_RL) THEN
-            tmpscal2= - tmpscal2*
-      &       ( ZERO + HeatCapacity_Cp * temp_EvPrRn )
-            ELSE
-            tmpscal2= - tmpscal2*
-      &       ( ZERO + HeatCapacity_Cp * theta(I,J,kSurface,bi,bj) )
-            ENDIF
- c add to SIatmQnt, leading to SItflux, which is analogous to TFLUX
-            DIAGarrayA(I,J)=maskC(I,J,kSurface,bi,bj)*convertHI2Q*(
- #ifndef SEAICE_GROWTH_LEGACY
-      &            a_QSWbyATM_cover(I,J) +
- #endif
-      &            a_QbyATM_cover(I,J) + a_QbyATM_open(I,J) )
-      &            -tmpscal1-tmpscal2
-           ENDDO
-          ENDDO
-          CALL DIAGNOSTICS_FILL(DIAGarrayA,
-      &      'SItflux ',0,1,3,bi,bj,myThid)
  #endif /* ALLOW_ATM_TEMP */
          ENDIF
-Line 2341 
 C close bi,bj loops
+Line 2811 
 C close bi,bj loops
         ENDDO
        ENDDO
+ C ===================================================================
+ C =========PART 9: HF/FWF global integrals and balancing=============
+ C ===================================================================
+ #ifdef ALLOW_BALANCE_FLUXES
+ c 1)  global sums
+ # ifdef ALLOW_AUTODIFF_TAMC
+ CADJ STORE FWFsiTile = comlev1, key=ikey_dynamics, kind=isbyte
+ CADJ STORE HFsiTile = comlev1, key=ikey_dynamics, kind=isbyte
+ CADJ STORE FWF2HFsiTile = comlev1, key=ikey_dynamics, kind=isbyte
+ # endif /* ALLOW_AUTODIFF_TAMC */
+       FWFsiGlob=0. _d 0
+       IF ( balanceEmPmR )
+      &   CALL GLOBAL_SUM_TILE_RL( FWFsiTile, FWFsiGlob, myThid )
+       FWF2HFsiGlob=0. _d 0
+       IF ( balanceEmPmR.AND.(temp_EvPrRn.EQ.UNSET_RL) ) THEN
+          CALL GLOBAL_SUM_TILE_RL(FWF2HFsiTile, FWF2HFsiGlob, myThid)
+       ELSEIF ( balanceEmPmR ) THEN
+          FWF2HFsiGlob=HeatCapacity_Cp * temp_EvPrRn * globalArea
+       ENDIF
+       HFsiGlob=0. _d 0
+       IF ( balanceQnet )
+      &   CALL GLOBAL_SUM_TILE_RL( HFsiTile, HFsiGlob, myThid )
+ c 2) global means
+ c mean SIatmFW
+       tmpscal0=FWFsiGlob / globalArea
+ c corresponding mean advection by atm to ocean+ice water exchange
+ c        (if mean SIatmFW was removed uniformely from ocean)
+       tmpscal1=FWFsiGlob / globalArea * FWF2HFsiGlob / globalArea
+ c mean SItflux (before potential adjustement due to SIatmFW)
+       tmpscal2=HFsiGlob / globalArea
+ c mean SItflux (after potential adjustement due to SIatmFW)
+       IF ( balanceEmPmR ) tmpscal2=tmpscal2-tmpscal1
+ c 3) balancing adjustments
+       IF ( balanceEmPmR ) THEN
+       DO bj=myByLo(myThid),myByHi(myThid)
+        DO bi=myBxLo(myThid),myBxHi(myThid)
+         DO j=1-OLy,sNy+OLy
+          DO i=1-OLx,sNx+OLx
+             empmr(i,j,bi,bj) = empmr(i,j,bi,bj) - tmpscal0
+             SIatmFW(i,j,bi,bj) = SIatmFW(i,j,bi,bj) - tmpscal0
+ c           adjust SItflux consistently
+             IF ( (temp_EvPrRn.NE.UNSET_RL).AND.
+      &        useRealFreshWaterFlux.AND.(nonlinFreeSurf.NE.0) ) THEN
+             tmpscal1=
+      &       ( ZERO + HeatCapacity_Cp * temp_EvPrRn )
+             ELSEIF ( (temp_EvPrRn.EQ.UNSET_RL).AND.
+      &        useRealFreshWaterFlux.AND.(nonlinFreeSurf.NE.0) ) THEN
+             tmpscal1=
+      &       ( ZERO + HeatCapacity_Cp * theta(I,J,kSurface,bi,bj) )
+             ELSEIF ( (temp_EvPrRn.NE.UNSET_RL) ) THEN
+             tmpscal1=
+      &       HeatCapacity_Cp*(temp_EvPrRn - theta(I,J,kSurface,bi,bj))
+             ELSE
+             tmpscal1=ZERO
+             ENDIF
+             SItflux(i,j,bi,bj) = SItflux(i,j,bi,bj) - tmpscal0*tmpscal1
+ c           no qnet or tflux adjustement is needed
+          ENDDO
+         ENDDO
+        ENDDO
+       ENDDO
+       IF ( balancePrintMean ) THEN
+        _BEGIN_MASTER( myThid )
+        WRITE(msgBuf,'(a,a,e24.17)') 'rm Global mean of ',
+      &      'SIatmFW = ', tmpscal0
+        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
+      &      SQUEEZE_RIGHT , myThid)
+        _END_MASTER( myThid )
+       ENDIF
+       ENDIF
+       IF ( balanceQnet ) THEN
+       DO bj=myByLo(myThid),myByHi(myThid)
+        DO bi=myBxLo(myThid),myBxHi(myThid)
+         DO j=1-OLy,sNy+OLy
+          DO i=1-OLx,sNx+OLx
+             SItflux(i,j,bi,bj) = SItflux(i,j,bi,bj) - tmpscal2
+             qnet(i,j,bi,bj) = qnet(i,j,bi,bj) - tmpscal2
+             SIatmQnt(i,j,bi,bj) = SIatmQnt(i,j,bi,bj) - tmpscal2
+          ENDDO
+         ENDDO
+        ENDDO
+       ENDDO
+       IF ( balancePrintMean ) THEN
+        _BEGIN_MASTER( myThid )
+        WRITE(msgBuf,'(a,a,e24.17)') 'rm Global mean of ',
+      &      'SItflux = ', tmpscal2
+        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
+      &      SQUEEZE_RIGHT , myThid)
+        _END_MASTER( myThid )
+       ENDIF
+       ENDIF
+ #endif /* ALLOW_BALANCE_FLUXES */
+ #ifdef ALLOW_DIAGNOSTICS
+ c these diags need to be done outside of the bi,bj loop so that
+ c we may do potential global mean adjustement to them consistently.
+          CALL DIAGNOSTICS_FILL(SItflux,
+      &      'SItflux ',0,1,0,1,1,myThid)
+          CALL DIAGNOSTICS_FILL(SIatmQnt,
+      &      'SIatmQnt',0,1,0,1,1,myThid)
+ c SIatmFW follows the same convention as empmr -- SIatmFW diag does not
+          tmpscal1= - 1. _d 0
+          CALL DIAGNOSTICS_SCALE_FILL(SIatmFW,
+      &      tmpscal1,1,'SIatmFW ',0,1,0,1,1,myThid)
+ #endif /* ALLOW_DIAGNOSTICS */
        RETURN
        END

 Legend:



Removed from v.1.2
 


changed lines


 
Added in v.1.14
 Legend:



Removed from v.1.2
 


changed lines


 
Added in v.1.14
-Removed from v.1.2
+Added in v.1.14

	ViewVC Help
Powered by ViewVC 1.1.22