/[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.4 by torge,
Fri Sep 28 19:01:17 2012 UTC
+revision 1.10 by torge,
Mon Oct 22 21:12:47 2012 UTC
 Line 58 
 C     !FUNCTIONS:
  C     !LOCAL VARIABLES:
  C     === Local variables ===
+ #ifdef SEAICE_DEBUG
  c ToM<<< debug seaice_growth
  C     msgBuf      :: Informational/error message buffer
        CHARACTER*(MAX_LEN_MBUF) msgBuf
  c ToM>>>
+ #endif
  C
  C unit/sign convention:
  C    Within the thermodynamic computation all stocks, except HSNOW,
-Line 95 
 C     i,j,bi,bj :: Loop counters
+Line 97 
 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
+       _RL pFac
  C     constants
        _RL tempFrz, ICE2SNOW, SNOW2ICE
        _RL QI, QS, recip_QI
+       _RL lhSublim
- C--   TmixLoc        :: ocean surface/mixed-layer temperature (in K)
-       _RL TmixLoc       (1:sNx,1:sNy)
- 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
- C     r_QbyATM_cover :: residual of a_QbyATM_cover after freezing/melting processes
- C     r_QbyATM_open  :: same but for open water
-       _RL a_QbyATM_cover      (1:sNx,1:sNy)
-       _RL a_QbyATM_open       (1:sNx,1:sNy)
-       _RL r_QbyATM_cover      (1:sNx,1:sNy)
-       _RL r_QbyATM_open       (1:sNx,1:sNy)
- C     a_QSWbyATM_open   - short wave heat flux over ocean in W/m^2
- 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             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)
        _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
- #ifdef ALLOW_DIAGNOSTICS
+ C     wind speed square
- C ICE/SNOW stocks tendencies associated with the various melt/freeze processes
+       _RL SPEED_SQ
-       _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
-       _RL d_HEFFbyOCNonICE    (1:sNx,1:sNy)
- c     The sum of mean ice thickness increments due to atmospheric fluxes over the open water
- c     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
-       _RL d_HEFFbyFLOODING    (1:sNx,1:sNy)
- c     The mean ice thickness increments due to atmospheric fluxes over the open water
+ C     Regularization values squared
- c     fraction and ice-covered fractions of the grid cell, respectively
+       _RL area_reg_sq, hice_reg_sq
-       _RL d_HEFFbyATMonOCN_open(1:sNx,1:sNy)
+ C     pathological cases thresholds
-       _RL d_HEFFbyATMonOCN_cover(1:sNx,1:sNy)
+       _RL heffTooHeavy
-       _RL d_HSNWbyNEG         (1:sNx,1:sNy)
+ C     Helper variables: reciprocal of some constants
-       _RL d_HSNWbyATMonSNW    (1:sNx,1:sNy)
+       _RL recip_multDim
-       _RL d_HSNWbyOCNonSNW    (1:sNx,1:sNy)
+       _RL recip_deltaTtherm
-       _RL d_HSNWbyRAIN        (1:sNx,1:sNy)
+       _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
-       _RL d_HFRWbyRAIN        (1:sNx,1:sNy)
+ C     temporary variables available for the various computations
- C
+       _RL tmpscal0, tmpscal1, tmpscal2, tmpscal3, tmpscal4
- C     a_FWbySublim :: fresh water flux implied by latent heat of
- C                     sublimation to atmosphere, same sign convention
- C                     as EVAP (positive upward)
-       _RL a_FWbySublim        (1:sNx,1:sNy)
-       _RL r_FWbySublim        (1:sNx,1:sNy)
-       _RL d_HEFFbySublim      (1:sNx,1:sNy)
-       _RL d_HSNWbySublim      (1:sNx,1:sNy)
- #if (defined(SEAICE_CAP_SUBLIM) || defined(SEAICE_ITD))
+ #ifdef ALLOW_SITRACER
- C     The latent heat flux which will sublimate all snow and ice
+       INTEGER iTr
- C     over one time step
+ #ifdef ALLOW_DIAGNOSTICS
-       _RL latentHeatFluxMax   (1:sNx,1:sNy)
+       CHARACTER*8   diagName
-       _RL latentHeatFluxMaxMult  (1:sNx,1:sNy,MULTDIM)
  #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)
-Line 192 
 C     actual snow thickness
+Line 155 
 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)
-Line 206 
 C     AREA_PRE :: hold sea-ice fraction
+Line 165 
 C     AREA_PRE :: hold sea-ice fraction
        _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)
+       _RL leadIceThickMin
  #endif
  C     wind speed
        _RL UG                  (1:sNx,1:sNy)
- #ifdef ALLOW_ATM_WIND
-       _RL SPEED_SQ
- #endif
- 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
+ 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
-       INTEGER it
- #ifdef SEAICE_ITD
-       INTEGER K
- #endif
-       _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)
+ #ifdef SEAICE_ITD
        _RL recip_heffActualMult(1:sNx,1:sNy,MULTDIM)
  #endif
        _RL a_QbyATMmult_cover  (1:sNx,1:sNy,MULTDIM)
-Line 251 
 C temporary variables available for the
+Line 191 
 C temporary variables available for the
        _RL r_QbyATMmult_cover  (1:sNx,1:sNy,MULTDIM)
        _RL r_FWbySublimMult    (1:sNx,1:sNy,MULTDIM)
  #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     a_QbyATM_cover :: available heat (in W/m^2) due to the interaction of
- C     ice is absent in a grid cell  (dimensionless)
+ C             the atmosphere and the ocean surface - for ice covered water
-       _RL MixedLayerTurbulenceFactor
+ C     a_QbyATM_open  :: same but for open water
+ C     r_QbyATM_cover :: residual of a_QbyATM_cover after freezing/melting processes
+ C     r_QbyATM_open  :: same but for open water
+       _RL a_QbyATM_cover      (1:sNx,1:sNy)
+       _RL a_QbyATM_open       (1:sNx,1:sNy)
+       _RL r_QbyATM_cover      (1:sNx,1:sNy)
+       _RL r_QbyATM_open       (1:sNx,1:sNy)
+ C     a_QSWbyATM_open   - short wave heat flux over ocean in W/m^2
+ 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 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)
- #ifdef ALLOW_SITRACER
+ C     The change of mean ice thickness due to out-of-bounds values following
-       INTEGER iTr
+ C     sea ice dyhnamics
-       CHARACTER*8   diagName
+       _RL d_HEFFbyNEG         (1:sNx,1:sNy)
+ 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
+ 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
+       _RL d_HEFFbyFLOODING    (1:sNx,1:sNy)
+ C     The mean ice thickness increments due to atmospheric fluxes over the open
+ 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)
+       _RL d_HSNWbyNEG         (1:sNx,1:sNy)
+       _RL d_HSNWbyATMonSNW    (1:sNx,1:sNy)
+       _RL d_HSNWbyOCNonSNW    (1:sNx,1:sNy)
+       _RL d_HSNWbyRAIN        (1:sNx,1:sNy)
+       _RL d_HFRWbyRAIN        (1:sNx,1:sNy)
+ C     a_FWbySublim :: fresh water flux implied by latent heat of
+ C                     sublimation to atmosphere, same sign convention
+ C                     as EVAP (positive upward)
+       _RL a_FWbySublim        (1:sNx,1:sNy)
+       _RL r_FWbySublim        (1:sNx,1:sNy)
+       _RL d_HEFFbySublim      (1:sNx,1:sNy)
+       _RL d_HSNWbySublim      (1:sNx,1:sNy)
+ #ifdef SEAICE_CAP_SUBLIM
+ 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)
  #endif
+ #ifdef EXF_ALLOW_SEAICE_RELAX
+ 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 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
+       CHARACTER*(max_len_mbuf) msgbuf
+ #endif
  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
-Line 286 
 C ======================================
+Line 297 
 C ======================================
        ENDIF
  C     avoid unnecessary divisions in loops
- #ifdef SEAICE_ITD
+ c#ifdef SEAICE_ITD
- CToM  SEAICE_multDim = nITD (see SEAICE_SIZE.h and seaice_readparms.F)
+ CToM this is now set by MULTDIM = nITD in SEAICE_SIZE.h
- #endif
+ 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 339 
 C conversion factors to go from precip (
+Line 352 
 C conversion factors to go from precip (
       &                       + act4*max1*max2*max3
  #endif /* ALLOW_AUTODIFF_TAMC */
  C array initializations
  C =====================
-Line 362 
 C =====================
+Line 374 
 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 389 
 C =====================
+Line 398 
 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 403 
 c
+Line 409 
 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_CAP_SUBLIM
+             latentHeatFluxMaxMult(I,J,IT) = 0.0 _d 0
+ #endif
  #ifdef SEAICE_ITD
              r_QbyATMmult_cover (I,J,IT)   = 0.0 _d 0
              r_FWbySublimMult(I,J,IT)      = 0.0 _d 0
  #endif
- #if (defined(SEAICE_CAP_SUBLIM) || defined(SEAICE_ITD))
-             latentHeatFluxMaxMult(I,J,IT) = 0.0 _d 0
- #endif
            ENDDO
           ENDDO
          ENDDO
-Line 421 
 c
+Line 427 
 c
          ENDDO
  #endif
  C =====================================================================
  C ===========PART 1: treat pathological cases (post advdiff)===========
  C =====================================================================
-Line 444 
 C ======================================
+Line 449 
 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)
            ENDDO
           ENDDO
          ENDDO
-Line 462 
 Cgf no dependency through pathological c
+Line 467 
 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 492 
 C           d_HEFFbyRLX(i,j) = 1. _d 1 *
+Line 497 
 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 area(:,:,bi,bj) = comlev1_bib
+Line 509 
 CADJ STORE area(:,:,bi,bj) = comlev1_bib
          DO J=1,sNy
           DO I=1,sNx
  #ifdef SEAICE_ITD
-           DO K=1,nITD
+           DO IT=1,nITD
-            tmpscal1=0. _d 0
             tmpscal2=0. _d 0
             tmpscal3=0. _d 0
-            tmpscal2=MAX(-HEFFITD(I,J,K,bi,bj),0. _d 0)
+            tmpscal2=MAX(-HEFFITD(I,J,IT,bi,bj),0. _d 0)
-            HEFFITD(I,J,K,bi,bj)=HEFFITD(I,J,K,bi,bj)+tmpscal2
+            HEFFITD(I,J,IT,bi,bj)=HEFFITD(I,J,IT,bi,bj)+tmpscal2
             d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2
-            tmpscal3=MAX(-HSNOWITD(I,J,K,bi,bj),0. _d 0)
+            tmpscal3=MAX(-HSNOWITD(I,J,IT,bi,bj),0. _d 0)
-            HSNOWITD(I,J,K,bi,bj)=HSNOWITD(I,J,K,bi,bj)+tmpscal3
+            HSNOWITD(I,J,IT,bi,bj)=HSNOWITD(I,J,IT,bi,bj)+tmpscal3
             d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3
-            tmpscal1=MAX(-AREAITD(I,J,K,bi,bj),0. _d 0)
+            AREAITD(I,J,IT,bi,bj)=MAX(AREAITD(I,J,IT,bi,bj),0. _d 0)
-            AREAITD(I,J,K,bi,bj)=AREAITD(I,J,K,bi,bj)+tmpscal1
-            d_AREAbyNEG(I,J)=d_AREAbyNEG(I,J)+tmpscal1
            ENDDO
  CToM      AREA, HEFF, and HSNOW will be updated at end of PART 1
  C         by calling SEAICE_ITD_SUM
  #else
            d_HEFFbyNEG(I,J)=MAX(-HEFF(I,J,bi,bj),0. _d 0)
-           d_HSNWbyNEG(I,J)=MAX(-HSNOW(I,J,bi,bj),0. _d 0)
-           AREA(I,J,bi,bj)=MAX(AREA(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
           ENDDO
          ENDDO
-Line 538 
 CADJ STORE heff(:,:,bi,bj)  = comlev1_bi
+Line 540 
 CADJ STORE heff(:,:,bi,bj)  = comlev1_bi
          DO J=1,sNy
           DO I=1,sNx
  #ifdef SEAICE_ITD
-           DO K=1,nITD
+           DO IT=1,nITD
  #endif
            tmpscal2=0. _d 0
            tmpscal3=0. _d 0
  #ifdef SEAICE_ITD
-            IF (HEFFITD(I,J,K,bi,bj).LE.siEps) THEN
+            IF (HEFFITD(I,J,IT,bi,bj).LE.siEps) THEN
-             tmpscal2=-HEFFITD(I,J,K,bi,bj)
+             tmpscal2=-HEFFITD(I,J,IT,bi,bj)
-             tmpscal3=-HSNOWITD(I,J,K,bi,bj)
+             tmpscal3=-HSNOWITD(I,J,IT,bi,bj)
-             TICES(I,J,K,bi,bj)=celsius2K
+             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,K,bi,bj) =HEFFITD(I,J,K,bi,bj) +tmpscal2
+            HEFFITD(I,J,IT,bi,bj) =HEFFITD(I,J,IT,bi,bj) +tmpscal2
-            HSNOWITD(I,J,K,bi,bj)=HSNOWITD(I,J,K,bi,bj)+tmpscal3
+            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)
-Line 580 
 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
+Line 582 
 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
          DO J=1,sNy
           DO I=1,sNx
  #ifdef SEAICE_ITD
-           DO K=1,nITD
+           DO IT=1,nITD
-            IF ((HEFFITD(i,j,k,bi,bj).EQ.0. _d 0).AND.
+            IF ((HEFFITD(I,J,IT,bi,bj).EQ.0. _d 0).AND.
-      &         (HSNOWITD(i,j,k,bi,bj).EQ.0. _d 0))
+      &         (HSNOWITD(I,J,IT,bi,bj).EQ.0. _d 0))
-      &      AREAITD(I,J,K,bi,bj)=0. _d 0
+      &      AREAITD(I,J,IT,bi,bj)=0. _d 0
            ENDDO
  #else
            IF ((HEFF(i,j,bi,bj).EQ.0. _d 0).AND.
-Line 601 
 CADJ STORE area(:,:,bi,bj)  = comlev1_bi
+Line 603 
 CADJ STORE area(:,:,bi,bj)  = comlev1_bi
          DO J=1,sNy
           DO I=1,sNx
  #ifdef SEAICE_ITD
-           DO K=1,nITD
+           DO IT=1,nITD
-            IF ((HEFFITD(i,j,k,bi,bj).GT.0).OR.
+            IF ((HEFFITD(I,J,IT,bi,bj).GT.0).OR.
-      &         (HSNOWITD(i,j,k,bi,bj).GT.0)) THEN
+      &         (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,K,bi,bj)=
+             AREAITD(I,J,IT,bi,bj)=
-      &         MAX(AREAITD(I,J,K,bi,bj),SEAICE_area_floor/float(nITD))
+      &         MAX(AREAITD(I,J,IT,bi,bj),SEAICE_area_floor/float(nITD))
             ENDIF
            ENDDO
  #else
-Line 639 
 CADJ STORE area(:,:,bi,bj)  = comlev1_bi
+Line 641 
 CADJ STORE area(:,:,bi,bj)  = comlev1_bi
            AREA(I,J,bi,bj)=MIN(AREA(I,J,bi,bj),SEAICE_area_max)
           ENDDO
          ENDDO
- #endif /* SEAICE_ITD */
+ #endif /* notSEAICE_ITD */
  #ifdef SEAICE_ITD
  CToM catch up with items 1.25 and 2.5 involving category sums AREA and HEFF
- C    first, update AREA and HEFF:
-         CALL SEAICE_ITD_SUM(bi, bj, myTime, myIter, myThid)
- C
          DO J=1,sNy
           DO I=1,sNx
  C    TICES was changed above (item 1.25), now update TICE as ice volume
  C     weighted average of TICES
+ C    also compute total of AREAITD (needed for finishing item 2.5, see below)
            tmpscal1 = 0. _d 0
            tmpscal2 = 0. _d 0
-           DO K=1,nITD
+           tmpscal3 = 0. _d 0
-            tmpscal1=tmpscal1 + TICES(I,J,K,bi,bj)*HEFFITD(I,J,K,bi,bj)
+           DO IT=1,nITD
-            tmpscal2=tmpscal2 + HEFFITD(I,J,K,bi,bj)
+            tmpscal1=tmpscal1 + TICES(I,J,IT,bi,bj)*HEFFITD(I,J,IT,bi,bj)
+            tmpscal2=tmpscal2 + HEFFITD(I,J,IT,bi,bj)
+            tmpscal3=tmpscal3 + AREAITD(I,J,IT,bi,bj)
            ENDDO
            TICE(I,J,bi,bj)=tmpscal1/tmpscal2
  C    lines of item 2.5 that were omitted:
-Line 662 
 C    in 2.5 these lines are executed bef
+Line 664 
 C    in 2.5 these lines are executed bef
  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) = AREA(I,J,bi,bj)
+           DIAGarrayA(I,J) = tmpscal3
  #endif
  #ifdef ALLOW_SITRACER
-           SItrAREA(I,J,bi,bj,1)=AREA(I,J,bi,bj)
+           SItrAREA(I,J,bi,bj,1)=tmpscal3
  #endif
           ENDDO
          ENDDO
-Line 676 
 C    and update AREA, HEFF, and HSNOW
+Line 678 
 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)
+ #ifdef SEAICE_DEBUG
+ c ToM<<< debug seaice_growth
+         WRITE(msgBuf,'(A,7F8.4)')
+      &    ' SEAICE_GROWTH: Heff increments 0, HEFFITD = ',
+      &     HEFFITD(1,1,:,bi,bj)
+         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
+      &    SQUEEZE_RIGHT , myThid)
+         WRITE(msgBuf,'(A,7F8.4)')
+      &    ' SEAICE_GROWTH: Area increments 0, AREAITD = ',
+      &     AREAITD(1,1,:,bi,bj)
+         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
+      &    SQUEEZE_RIGHT , myThid)
  #endif
+ #else
+ #ifdef SEAICE_DEBUG
+         WRITE(msgBuf,'(A,7F8.4)')
+      &    ' SEAICE_GROWTH: Heff increments 0, HEFF = ',
+      &     HEFF(1,1,bi,bj)
+         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
+      &    SQUEEZE_RIGHT , myThid)
+         WRITE(msgBuf,'(A,7F8.4)')
+      &    ' SEAICE_GROWTH: Area increments 0, AREA = ',
+      &     AREA(1,1,bi,bj)
+         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
+      &    SQUEEZE_RIGHT , myThid)
+ c ToM>>>
+ #endif
+ #endif /* SEAICE_ITD */
  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
- C        ENDIF SEAICEadjMODE.EQ.0
+ C        end SEAICEadjMODE.EQ.0 statement:
          ENDIF
  #endif
-Line 700 
 C 3) store regularized values of heff, h
+Line 729 
 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,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)
  C memorize areal and volume fraction of each ITD category
-            IF (AREA(I,J,bi,bj).GT.0) THEN
+            IF (AREA(I,J,bi,bj) .GT. ZERO) THEN
-             areaFracFactor(I,J,K)=AREAITD(I,J,K,bi,bj)/AREA(I,J,bi,bj)
+             areaFracFactor(I,J,IT)=AREAITD(I,J,IT,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
+ C           if there's 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
  C prepare SItrHEFF to be computed as cumulative sum
-         DO K=2,5
+         DO iTr=2,5
           DO J=1,sNy
            DO I=1,sNx
-            SItrHEFF(I,J,bi,bj,K)=ZERO
+            SItrHEFF(I,J,bi,bj,iTr)=ZERO
            ENDDO
           ENDDO
          ENDDO
-Line 791 
 Cgf no additional dependency of air-sea
+Line 820 
 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 821 
 CADJ STORE HEFFpreTH = comlev1_bibj, key
+Line 850 
 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/float(nITD),
-      &                              AREAITDpreTH(I,J,K))
+      &                              AREAITDpreTH(I,J,IT))
-             hsnowActualMult(I,J,K) = HSNWITDpreTH(I,J,K)/tmpscal1
+             hsnowActualMult(I,J,IT) = HSNWITDpreTH(I,J,IT)/tmpscal1
-             tmpscal2               = HEFFITDpreTH(I,J,K)/tmpscal1
+             tmpscal2               = HEFFITDpreTH(I,J,IT)/tmpscal1
-             heffActualMult(I,J,K)  = MAX(tmpscal2,SEAICE_hice_reg)
+             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 /* SEAICE_ITD */
            IF (HEFFpreTH(I,J) .GT. ZERO) THEN
-Line 864 
 cif       Do not regularize when HEFFpre
+Line 893 
 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
-Line 900 
 cif       Do not regularize when HEFFpre
+Line 929 
 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 963 
 cCADJ STORE TmixLoc = comlev1_bibj, key
+Line 993 
 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 985 
 C     Compute relative wind speed over s
+Line 1014 
 C     Compute relative wind speed over s
            ENDDO
           ENDDO
          ENDIF
- #endif /* ALLOW_ATM_WIND */
  #ifdef ALLOW_AUTODIFF_TAMC
  CADJ STORE tice(:,:,bi,bj)
-Line 1004 
 C--   Start loop over multi-categories
+Line 1032 
 C--   Start loop over multi-categories
  CToM  SEAICE_multDim = nITD (see SEAICE_SIZE.h and seaice_readparms.F)
  #endif
          DO IT=1,SEAICE_multDim
- c        homogeneous distribution between 0 and 2 x heffActual
+ C        homogeneous distribution between 0 and 2 x heffActual
  #ifndef SEAICE_ITD
-          pFac = (2.0 _d 0*real(IT)-1.0 _d 0)*recip_multDim
+          pFac = (2.0 _d 0*IT - 1.0 _d 0)*recip_multDim
+          pFacSnow = 1. _d 0
+          IF ( SEAICE_useMultDimSnow ) pFacSnow=pFac
  #endif
           DO J=1,sNy
            DO I=1,sNx
-Line 1014 
 c        homogeneous distribution betwee
+Line 1044 
 c        homogeneous distribution betwee
  CToM for SEAICE_ITD heffActualMult and latentHeatFluxMaxMult are calculated above
  C    (instead of heffActual and latentHeatFluxMax)
             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
-Line 1028 
 C    (instead of heffActual and latentHe
+Line 1059 
 C    (instead of heffActual and latentHe
  #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 1043 
 CADJ &     comlev1_bibj, key = iicekey,
+Line 1075 
 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 1079 
 C     update TICE & TICES
+Line 1109 
 C     update TICE & TICES
  #ifdef SEAICE_ITD
  C     calculate area weighted mean
  C     (although the ice's temperature relates to its energy content
- C      and hence should be averaged weighted by ice volume [heffFracFactor],
+ 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 1095 
 C     average over categories
+Line 1125 
 C     average over categories
  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 1116 
 C     (fluxes are per unit (ice surface)
+Line 1146 
 C     (fluxes are per unit (ice surface)
  # 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 1141 
 CADJ STORE a_FWbySublim    = comlev1_bib
+Line 1171 
 CADJ STORE a_FWbySublim    = comlev1_bib
  C switch heat fluxes from W/m2 to 'effective' ice meters
  #ifdef SEAICE_ITD
-         DO K=1,nITD
+         DO IT=1,nITD
           DO J=1,sNy
            DO I=1,sNx
-            a_QbyATMmult_cover(I,J,K)   = a_QbyATMmult_cover(I,J,K)
+            a_QbyATMmult_cover(I,J,IT)   = a_QbyATMmult_cover(I,J,IT)
-      &          * convertQ2HI * AREAITDpreTH(I,J,K)
+      &          * convertQ2HI * AREAITDpreTH(I,J,IT)
-            a_QSWbyATMmult_cover(I,J,K) = a_QSWbyATMmult_cover(I,J,K)
+            a_QSWbyATMmult_cover(I,J,IT) = a_QSWbyATMmult_cover(I,J,IT)
-      &          * convertQ2HI * AREAITDpreTH(I,J,K)
+      &          * convertQ2HI * AREAITDpreTH(I,J,IT)
- C and initialize r_QbyATM_cover
+ C and initialize r_QbyATMmult_cover
-            r_QbyATMmult_cover(I,J,K)=a_QbyATMmult_cover(I,J,K)
+            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,K) = ZERO
+            a_FWbySublimMult(I,J,IT) = ZERO
  #endif
-            a_FWbySublimMult(I,J,K) = SEAICE_deltaTtherm*recip_rhoIce
+            a_FWbySublimMult(I,J,IT) = SEAICE_deltaTtherm*recip_rhoIce
-      &            * a_FWbySublimMult(I,J,K)*AREAITDpreTH(I,J,K)
+      &            * a_FWbySublimMult(I,J,IT)*AREAITDpreTH(I,J,IT)
-            r_FWbySublimMult(I,J,K)=a_FWbySublimMult(I,J,K)
+            r_FWbySublimMult(I,J,IT)=a_FWbySublimMult(I,J,IT)
            ENDDO
           ENDDO
          ENDDO
-Line 1188 
 C and initialize r_QbyATM_cover/r_QbyATM
+Line 1218 
 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)
-Line 1214 
 CADJ STORE r_FWbySublim    = comlev1_bib
+Line 1244 
 CADJ STORE r_FWbySublim    = comlev1_bib
  Cgf no additional dependency through ice cover
          IF ( SEAICEadjMODE.GE.3 ) THEN
  #ifdef SEAICE_ITD
-          DO K=1,nITD
+          DO IT=1,nITD
            DO J=1,sNy
             DO I=1,sNx
-             a_QbyATMmult_cover(I,J,K)   = 0. _d 0
+             a_QbyATMmult_cover(I,J,IT)   = 0. _d 0
-             r_QbyATMmult_cover(I,J,K)   = 0. _d 0
+             r_QbyATMmult_cover(I,J,IT)   = 0. _d 0
-             a_QSWbyATMmult_cover(I,J,K) = 0. _d 0
+             a_QSWbyATMmult_cover(I,J,IT) = 0. _d 0
             ENDDO
            ENDDO
           ENDDO
-Line 1248 
 CADJ &                       key = iicek
+Line 1278 
 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 1268 
 c efficiency of turbulent fluxes : depen
+Line 1298 
 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 1283 
 c available turbulent flux
+Line 1313 
 c available turbulent flux
          CALL ZERO_ADJ_1D( sNx*sNy, r_QbyOCN, myThid)
  #endif
  C ===================================================================
  C =========PART 3: determine effective thicknesses increments========
  C ===================================================================
-Line 1296 
 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
+Line 1325 
 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)
  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
+           HSNOWITD(I,J,IT,bi,bj)  = HSNOWITD(I,J,IT,bi,bj)   - tmpscal2
       &                                                       *ICE2SNOW
-           r_FWbySublimMult(I,J,K) = r_FWbySublimMult(I,J,K) - tmpscal2
+           r_FWbySublimMult(I,J,IT)= r_FWbySublimMult(I,J,IT) - 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 1330 
 CADJ STORE r_FWbySublim    = comlev1_bib
+Line 1357 
 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)
  C         accumulate change over ITD categories
-           d_HSNWbySublim(I,J)     = d_HSNWbySublim(I,J)     - tmpscal2
+           d_HSNWbySublim(I,J)      = d_HSNWbySublim(I,J)      - tmpscal2
-           HEFFITD(I,J,K,bi,bj)    = HEFFITD(I,J,K,bi,bj)    - tmpscal2
+           HEFFITD(I,J,IT,bi,bj)    = HEFFITD(I,J,IT,bi,bj)    - tmpscal2
-           r_FWbySublimMult(I,J,K) = r_FWbySublimMult(I,J,K) - tmpscal2
+           r_FWbySublimMult(I,J,IT) = r_FWbySublimMult(I,J,IT) - tmpscal2
- C         keep total up to date, too
-           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 1351 
 C     If anything is left, it will be ev
+Line 1376 
 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,'(A,7F6.2)')
+         WRITE(msgBuf,'(A,7F8.4)')
  #ifdef SEAICE_ITD
       &    ' SEAICE_GROWTH: Heff increments 1, HEFFITD = ',
-      &     HEFFITD(20,20,:,bi,bj)
+      &     HEFFITD(1,1,:,bi,bj)
+         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
+      &    SQUEEZE_RIGHT , myThid)
+         WRITE(msgBuf,'(A,7F8.4)')
+      &    ' SEAICE_GROWTH: Area increments 1, AREAITD = ',
+      &     AREAITD(1,1,:,bi,bj)
  #else
       &    ' SEAICE_GROWTH: Heff increments 1, HEFF = ',
-      &     HEFF(20,20,bi,bj)
+      &     HEFF(1,1,bi,bj)
  #endif
          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
       &    SQUEEZE_RIGHT , myThid)
  c ToM>>>
+ #endif
  C compute ice thickness tendency due to ice-ocean interaction
  C ===========================================================
-Line 1389 
 CADJ STORE r_QbyOCN = comlev1_bibj,key=i
+Line 1419 
 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),
-            d_HEFFbyOCNonICE(I,J)= d_HEFFbyOCNonICE(I,J) + tmpscal1
+      &                               -HEFFITD(I,J,IT,bi,bj))
-            r_QbyOCN(I,J)        = r_QbyOCN(I,J)         - tmpscal1
+            d_HEFFbyOCNonICE(I,J) = d_HEFFbyOCNonICE(I,J) + tmpscal1
-            HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,bi,bj)  + tmpscal1
+            HEFFITD(I,J,IT,bi,bj) = HEFFITD(I,J,IT,bi,bj) + tmpscal1
  #ifdef ALLOW_SITRACER
             SItrHEFF(I,J,bi,bj,2) = SItrHEFF(I,J,bi,bj,2)
-      &                           + HEFFITD(I,J,K,bi,bj)
+      &                           + HEFFITD(I,J,IT,bi,bj)
  #endif
            ENDDO
           ENDDO
          ENDDO
+         DO J=1,sNy
+          DO I=1,sNx
+           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
-Line 1418 
 C          and hence weighted by fractio
+Line 1453 
 C          and hence weighted by fractio
           ENDDO
          ENDDO
  #endif /* SEAICE_ITD */
+ #ifdef SEAICE_DEBUG
  c ToM<<< debug seaice_growth
-         WRITE(msgBuf,'(A,7F6.2)')
+         WRITE(msgBuf,'(A,7F8.4)')
  #ifdef SEAICE_ITD
       &    ' SEAICE_GROWTH: Heff increments 2, HEFFITD = ',
-      &     HEFFITD(20,20,:,bi,bj)
+      &     HEFFITD(1,1,:,bi,bj)
+         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
+      &    SQUEEZE_RIGHT , myThid)
+         WRITE(msgBuf,'(A,7F8.4)')
+      &    ' SEAICE_GROWTH: Area increments 2, AREAITD = ',
+      &     AREAITD(1,1,:,bi,bj)
  #else
       &    ' SEAICE_GROWTH: Heff increments 2, HEFF = ',
-      &     HEFF(20,20,bi,bj)
+      &     HEFF(1,1,bi,bj)
  #endif
          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
       &    SQUEEZE_RIGHT , myThid)
  c ToM>>>
+ #endif
  C compute snow melt tendency due to snow-atmosphere interaction
  C ==================================================================
-Line 1440 
 CADJ STORE r_QbyATM_cover = comlev1_bibj
+Line 1482 
 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(I,J) = d_HSNWbyATMonSNW(I,J)
-           HSNOWITD(I,J,K,bi,bj)=HSNOWITD(I,J,K,bi,bj)+tmpscal2*ICE2SNOW
+      &                           + tmpscal2*ICE2SNOW
-           r_QbyATMmult_cover(I,J,K)=r_QbyATMmult_cover(I,J,K) - tmpscal2
+            HSNOWITD(I,J,IT,bi,bj)= HSNOWITD(I,J,IT,bi,bj)
- 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)
+      &                           - tmpscal2
            ENDDO
           ENDDO
          ENDDO
-Line 1477 
 Cgf no additional dependency through sno
+Line 1520 
 Cgf no additional dependency through sno
           ENDDO
          ENDDO
  #endif /* SEAICE_ITD */
+ #ifdef SEAICE_DEBUG
  c ToM<<< debug seaice_growth
-         WRITE(msgBuf,'(A,7F6.2)')
+         WRITE(msgBuf,'(A,7F8.4)')
  #ifdef SEAICE_ITD
       &    ' SEAICE_GROWTH: Heff increments 3, HEFFITD = ',
-      &     HEFFITD(20,20,:,bi,bj)
+      &     HEFFITD(1,1,:,bi,bj)
+         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
+      &    SQUEEZE_RIGHT , myThid)
+         WRITE(msgBuf,'(A,7F8.4)')
+      &    ' SEAICE_GROWTH: Area increments 3, AREAITD = ',
+      &     AREAITD(1,1,:,bi,bj)
  #else
       &    ' SEAICE_GROWTH: Heff increments 3, HEFF = ',
-      &     HEFF(20,20,bi,bj)
+      &     HEFF(1,1,bi,bj)
  #endif
          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
       &    SQUEEZE_RIGHT , myThid)
  c ToM>>>
+ #endif
  C compute ice thickness tendency due to the atmosphere
  C ====================================================
-Line 1504 
 Cgf the v1.81=>v1.82 revision would chan
+Line 1554 
 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
  #ifdef SEAICE_GROWTH_LEGACY
-            tmpscal2 =MAX(-HEFFITD(I,J,K,bi,bj),r_QbyATMmult_cover(I,J,K))
+            tmpscal2 = MAX(-HEFFITD(I,J,IT,bi,bj),
+      &                     r_QbyATMmult_cover(I,J,IT))
  #else
-            tmpscal2 =MAX(-HEFFITD(I,J,K,bi,bj),r_QbyATMmult_cover(I,J,K)
+            tmpscal2 = MAX(-HEFFITD(I,J,IT,bi,bj),
+      &                     r_QbyATMmult_cover(I,J,IT)
  c         Limit ice growth by potential melt by ocean
-      &     + AREAITDpreTH(I,J,K) * r_QbyOCN(I,J))
+      &              + AREAITDpreTH(I,J,IT) * r_QbyOCN(I,J))
  #endif /* SEAICE_GROWTH_LEGACY */
             d_HEFFbyATMonOCN_cover(I,J) = 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)
+            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
+            HEFFITD(I,J,IT,bi,bj) = HEFFITD(I,J,IT,bi,bj) + tmpscal2
  #ifdef ALLOW_SITRACER
             SItrHEFF(I,J,bi,bj,3) = SItrHEFF(I,J,bi,bj,3)
-      &                           + HEFFITD(I,J,K,bi,bj)
+      &                           + HEFFITD(I,J,IT,bi,bj)
  #endif
            ENDDO
           ENDDO
-Line 1537 
 c         Limit ice growth by potential
+Line 1589 
 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 */
-Line 1549 
 c         Limit ice growth by potential
+Line 1601 
 c         Limit ice growth by potential
  #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,'(A,7F6.2)')
+         WRITE(msgBuf,'(A,7F8.4)')
  #ifdef SEAICE_ITD
       &    ' SEAICE_GROWTH: Heff increments 4, HEFFITD = ',
-      &     HEFFITD(20,20,:,bi,bj)
+      &     HEFFITD(1,1,:,bi,bj)
+         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
+      &    SQUEEZE_RIGHT , myThid)
+         WRITE(msgBuf,'(A,7F8.4)')
+      &    ' SEAICE_GROWTH: Area increments 4, AREAITD = ',
+      &     AREAITD(1,1,:,bi,bj)
  #else
       &    ' SEAICE_GROWTH: Heff increments 4, HEFF = ',
-      &     HEFF(20,20,bi,bj)
+      &     HEFF(1,1,bi,bj)
  #endif
          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
       &    SQUEEZE_RIGHT , myThid)
  c ToM>>>
+ #endif
- 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
           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)
+            HSNOWITD(I,J,IT,bi,bj) = HSNOWITD(I,J,IT,bi,bj)
-      &     + d_HSNWbyRAIN(I,J)*areaFracFactor(I,J,K)
+      &     + d_HSNWbyRAIN(I,J)*areaFracFactor(I,J,IT)
            ENDDO
           ENDDO
          ENDDO
-Line 1611 
 C           add precip to the fresh wate
+Line 1684 
 C           add precip to the fresh wate
  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,'(A,7F6.2)')
+         WRITE(msgBuf,'(A,7F8.4)')
  #ifdef SEAICE_ITD
       &    ' SEAICE_GROWTH: Heff increments 5, HEFFITD = ',
-      &     HEFFITD(20,20,:,bi,bj)
+      &     HEFFITD(1,1,:,bi,bj)
+         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
+      &    SQUEEZE_RIGHT , myThid)
+         WRITE(msgBuf,'(A,7F8.4)')
+      &    ' SEAICE_GROWTH: Area increments 5, AREAITD = ',
+      &     AREAITD(1,1,:,bi,bj)
  #else
       &    ' SEAICE_GROWTH: Heff increments 5, HEFF = ',
-      &     HEFF(20,20,bi,bj)
+      &     HEFF(1,1,bi,bj)
  #endif
          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
       &    SQUEEZE_RIGHT , myThid)
  c ToM>>>
+ #endif
  C compute snow melt due to heat available from ocean.
  C =================================================================
-Line 1637 
 CADJ STORE r_QbyOCN = comlev1_bibj,key=i
+Line 1719 
 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),
+            tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW*areaFracFactor(I,J,IT),
-      &                  -HSNOWITD(I,J,K,bi,bj))
+      &                  -HSNOWITD(I,J,IT,bi,bj))
             tmpscal2=MIN(tmpscal1,0. _d 0)
  #ifdef SEAICE_MODIFY_GROWTH_ADJ
  Cgf no additional dependency through snow
-Line 1649 
 Cgf no additional dependency through sno
+Line 1731 
 Cgf no additional dependency through sno
  #endif
             d_HSNWbyOCNonSNW(I,J) = d_HSNWbyOCNonSNW(I,J) + tmpscal2
             r_QbyOCN(I,J)=r_QbyOCN(I,J) - tmpscal2*SNOW2ICE
-            HSNOWITD(I,J,K,bi,bj) = HSNOWITD(I,J,K,bi,bj) + tmpscal2
+            HSNOWITD(I,J,IT,bi,bj) = HSNOWITD(I,J,IT,bi,bj) + tmpscal2
            ENDDO
           ENDDO
          ENDDO
-Line 1671 
 Cgf no additional dependency through sno
+Line 1753 
 Cgf no additional dependency through sno
  #endif /* SEAICE_ITD */
  #endif /* SEAICE_EXCLUDE_FOR_EXACT_AD_TESTING */
  Cph)
+ #ifdef SEAICE_DEBUG
  c ToM<<< debug seaice_growth
-         WRITE(msgBuf,'(A,7F6.2)')
+         WRITE(msgBuf,'(A,7F8.4)')
  #ifdef SEAICE_ITD
       &    ' SEAICE_GROWTH: Heff increments 6, HEFFITD = ',
-      &     HEFFITD(20,20,:,bi,bj)
+      &     HEFFITD(1,1,:,bi,bj)
+         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
+      &    SQUEEZE_RIGHT , myThid)
+         WRITE(msgBuf,'(A,7F8.4)')
+      &    ' SEAICE_GROWTH: Area increments 6, AREAITD = ',
+      &     AREAITD(1,1,:,bi,bj)
  #else
       &    ' SEAICE_GROWTH: Heff increments 6, HEFF = ',
-      &     HEFF(20,20,bi,bj)
+      &     HEFF(1,1,bi,bj)
  #endif
          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
       &    SQUEEZE_RIGHT , myThid)
  c ToM>>>
+ #endif
  C gain of new ice over open water
  C ===============================
-Line 1696 
 CADJ STORE a_QSWbyATM_open = comlev1_bib
+Line 1785 
 CADJ STORE a_QSWbyATM_open = comlev1_bib
          DO J=1,sNy
           DO I=1,sNx
- c           Initial ice growth is triggered by open water
+ C           Initial ice growth is triggered by open water
- c           heat flux overcoming potential melt by ocean
+ 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)
-Line 1714 
 C           or 0. otherwise (no melting
+Line 1803 
 C           or 0. otherwise (no melting
  C         open water area fraction
            tmpscal0 = ONE-AREApreTH(I,J)
  C         determine thickness of new ice
- C         considering the entire open water area to refreeze
+ ctomC         considering the entire open water area to refreeze
-           tmpscal1 = tmpscal3/tmpscal0
+ ctom          tmpscal1 = tmpscal3/tmpscal0
+ C         considering a minimum lead ice thickness of 10 cm
+ C         WATCH that leadIceThickMin is smaller that Hlimit(1)!
+           leadIceThickMin = 0.1
+           tmpscal1 = MAX(leadIceThickMin,tmpscal3/tmpscal0)
+ C         adjust ice area fraction covered by new ice
+           tmpscal0 = tmpscal3/tmpscal1
  C         then add new ice volume to appropriate thickness category
-           DO K=1,nITD
+           DO IT=1,nITD
-            IF (tmpscal1.LT.Hlimit(K)) THEN
+            IF (tmpscal1.LT.Hlimit(IT)) THEN
-             HEFFITD(I,J,K,bi,bj) = HEFFITD(I,J,K,bi,bj) + tmpscal3
+             HEFFITD(I,J,IT,bi,bj) = HEFFITD(I,J,IT,bi,bj) + tmpscal3
              tmpscal3=ZERO
  C not sure if AREAITD should be changed here since AREA is incremented
  C   in PART 4 below in non-itd code
  C in this scenario all open water is covered by new ice instantaneously,
  C   i.e. no delayed lead closing is concidered such as is associated with
  C   Hibler's h_0 parameter
-             AREAITD(I,J,K,bi,bj) = AREAITD(I,J,K,bi,bj)
+             AREAITD(I,J,IT,bi,bj) = AREAITD(I,J,IT,bi,bj)
       &                           + tmpscal0
              tmpscal0=ZERO
             ENDIF
-Line 1739 
 C   Hibler's h_0 parameter
+Line 1834 
 C   Hibler's h_0 parameter
  #ifdef ALLOW_SITRACER
  #ifdef SEAICE_ITD
-         DO K=1,nITD
+         DO IT=1,nITD
           DO J=1,sNy
            DO I=1,sNx
  c needs to be here to allow use also with LEGACY branch
             SItrHEFF(I,J,bi,bj,4) = SItrHEFF(I,J,bi,bj,4)
-      &                           + HEFFITD(I,J,K,bi,bj)
+      &                           + HEFFITD(I,J,IT,bi,bj)
            ENDDO
           ENDDO
          ENDDO
  #else
          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,4)=HEFF(I,J,bi,bj)
           ENDDO
          ENDDO
  #endif
  #endif /* ALLOW_SITRACER */
+ #ifdef SEAICE_DEBUG
  c ToM<<< debug seaice_growth
-         WRITE(msgBuf,'(A,7F6.2)')
+         WRITE(msgBuf,'(A,7F8.4)')
  #ifdef SEAICE_ITD
       &    ' SEAICE_GROWTH: Heff increments 7, HEFFITD = ',
-      &     HEFFITD(20,20,:,bi,bj)
+      &     HEFFITD(1,1,:,bi,bj)
+         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
+      &    SQUEEZE_RIGHT , myThid)
+         WRITE(msgBuf,'(A,7F8.4)')
+      &    ' SEAICE_GROWTH: Area increments 7, AREAITD = ',
+      &     AREAITD(1,1,:,bi,bj)
  #else
       &    ' SEAICE_GROWTH: Heff increments 7, HEFF = ',
-      &     HEFF(20,20,bi,bj)
+      &     HEFF(1,1,bi,bj)
+         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
+      &    SQUEEZE_RIGHT , myThid)
+         WRITE(msgBuf,'(A,7F8.4)')
+      &    ' SEAICE_GROWTH: Area increments 7, AREA = ',
+      &     AREA(1,1,bi,bj)
  #endif
          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
       &    SQUEEZE_RIGHT , myThid)
  c ToM>>>
+ #endif
  C convert snow to ice if submerged.
  C =================================
-Line 1781 
 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
+Line 1888 
 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
+             tmpscal0 = (HSNOWITD(I,J,IT,bi,bj)*SEAICE_rhoSnow
-      &               +HEFFITD(I,J,K,bi,bj)*SEAICE_rhoIce)*recip_rhoConst
+      &               +  HEFFITD(I,J,IT,bi,bj) *SEAICE_rhoIce)
-             tmpscal1 = MAX( 0. _d 0, tmpscal0 - HEFFITD(I,J,K,bi,bj))
+      &                                        *recip_rhoConst
-             d_HEFFbyFLOODING(I,J) = d_HEFFbyFLOODING(I,J) + tmpscal1
+             tmpscal1 = MAX( 0. _d 0, tmpscal0 - HEFFITD(I,J,IT,bi,bj))
-             HEFFITD(I,J,K,bi,bj)  = HEFFITD(I,J,K,bi,bj)  + tmpscal1
+             d_HEFFbyFLOODING(I,J) = d_HEFFbyFLOODING(I,J)  + tmpscal1
-             HSNOWITD(I,J,K,bi,bj) = HSNOWITD(I,J,K,bi,bj) - tmpscal1
+             HEFFITD(I,J,IT,bi,bj) = HEFFITD(I,J,IT,bi,bj)  + tmpscal1
+             HSNOWITD(I,J,IT,bi,bj)= HSNOWITD(I,J,IT,bi,bj) - tmpscal1
       &                            * ICE2SNOW
             ENDDO
            ENDDO
-Line 1804 
 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
+Line 1912 
 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
             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,'(A,7F6.2)')
+         WRITE(msgBuf,'(A,7F8.4)')
  #ifdef SEAICE_ITD
       &    ' SEAICE_GROWTH: Heff increments 8, HEFFITD = ',
-      &     HEFFITD(20,20,:,bi,bj)
+      &     HEFFITD(1,1,:,bi,bj)
+         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
+      &    SQUEEZE_RIGHT , myThid)
+         WRITE(msgBuf,'(A,7F8.4)')
+      &    ' SEAICE_GROWTH: Area increments 8, AREAITD = ',
+      &     AREAITD(1,1,:,bi,bj)
  #else
       &    ' SEAICE_GROWTH: Heff increments 8, HEFF = ',
-      &     HEFF(20,20,bi,bj)
+      &     HEFF(1,1,bi,bj)
+         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
+      &    SQUEEZE_RIGHT , myThid)
+         WRITE(msgBuf,'(A,7F8.4)')
+      &    ' SEAICE_GROWTH: Area increments 8, AREA = ',
+      &     AREA(1,1,bi,bj)
  #endif
          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
       &    SQUEEZE_RIGHT , myThid)
  c ToM>>>
+ #endif
  C ===================================================================
  C ==========PART 4: determine ice cover fraction increments=========-
-Line 1851 
 C       replaces Hibler '79 scheme and l
+Line 1971 
 C       replaces Hibler '79 scheme and l
  C       because ITD accounts explicitly for lead openings and
  C       different melt rates due to varying ice thickness
  C
- C       only consider ice area loss due to total ice thickness loss
+ C       only consider ice area loss due to total ice thickness loss;
- C       ice area gain due to freezing of open water as handled above
+ C       ice area gain due to freezing of open water is handled above
  C       under "gain of new ice over open water"
  C
  C       does not account for lateral melt of ice floes
  C
  C        AREAITD is incremented in section "gain of new ice over open water" above
  C
-         DO K=1,nITD
+         DO IT=1,nITD
           DO J=1,sNy
-           DO I=1,sNx
+          DO I=1,sNx
-            IF (HEFFITD(I,J,K,bi,bj).LE.ZERO) THEN
+            IF (HEFFITD(I,J,IT,bi,bj).LE.ZERO) THEN
-             AREAITD(I,J,K,bi,bj)=ZERO
+             AREAITD(I,J,IT,bi,bj)=ZERO
             ENDIF
  #ifdef ALLOW_SITRACER
             SItrAREA(I,J,bi,bj,3) = SItrAREA(I,J,bi,bj,3)
-      &                           + AREAITD(I,J,K,bi,bj)
+      &                           + AREAITD(I,J,IT,bi,bj)
  #endif /* ALLOW_SITRACER */
            ENDDO
           ENDDO
-Line 1948 
 C apply tendency
+Line 2068 
 C apply tendency
  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
-Line 1995 
 CADJ &                       key = iicek
+Line 2115 
 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 2090 
 C set HSALT = 0 if HEFF = 0 and compute
+Line 2210 
 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 2173 
 C =================================
+Line 2292 
 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 2187 
 C compute net heat flux leaving/entering
+Line 2306 
 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 2200 
 C in principle a_QSWbyATM_cover should a
+Line 2355 
 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 2235 
 CADJ STORE d_HSNWbyATMonSNW = comlev1_bi
+Line 2396 
 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 2293 
 C ======================================
+Line 2516 
 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 2307 
 C     If r_FWbySublim>0, then it is evap
+Line 2530 
 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
  #ifdef ALLOW_SSH_GLOBMEAN_COST_CONTRIBUTION
  C--
-Line 2379 
 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
+Line 2612 
 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 2429 
 C three that actually need intermediate
+Line 2698 
 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 2468 
 C the actual Freshwater flux of sublimat
+Line 2711 
 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 2485 
 C the latent heat flux
+Line 2728 
 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 2531 
 C close bi,bj loops
+Line 2737 
 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 /* */
+ #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.4
 


changed lines


 
Added in v.1.10
 Legend:



Removed from v.1.4
 


changed lines


 
Added in v.1.10
-Removed from v.1.4
+Added in v.1.10

	ViewVC Help
Powered by ViewVC 1.1.22