/[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.8 by torge,
Mon Oct 22 16:02:25 2012 UTC
+revision 1.13 by torge,
Fri Nov  2 19:15:42 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
+       CHARACTER*12 msgBufForm
  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 98 
 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 conversion factors to go from Q (W/m2) to HEFF (ice meters)
+       _RL convertQ2HI, convertHI2Q
+ C conversion factors to go from precip (m/s) unit to HEFF (ice meters)
+       _RL convertPRECIP2HI, convertHI2PRECIP
+ C     Factor by which we increase the upper ocean friction velocity (u*) when
+ C     ice is absent in a grid cell  (dimensionless)
+       _RL MixedLayerTurbulenceFactor
+ C     wind speed square
+       _RL SPEED_SQ
+ C     Regularization values squared
+       _RL area_reg_sq, hice_reg_sq
+ C     pathological cases thresholds
+       _RL heffTooHeavy
+ C     Helper variables: reciprocal of some constants
+       _RL recip_multDim
+       _RL recip_deltaTtherm
+       _RL recip_rhoIce
+ C     local value (=1/HO or 1/HO_south)
+       _RL recip_HO
+ C     local value (=1/ice thickness)
+       _RL recip_HH
+ C     facilitate multi-category snow implementation
+       _RL pFacSnow
+ C     temporary variables available for the various computations
+       _RL tmpscal0, tmpscal1, tmpscal2, tmpscal3, tmpscal4
+ #ifdef SEAICE_ITD
+       _RL tmpscal1itd(1:sNx,1:sNy), tmpscal2itd(1:sNx,1:sNy)
+       _RL tmpscal3itd(1:sNx,1:sNy)
+ #endif
+ #ifdef ALLOW_SITRACER
+       INTEGER iTr
+ #ifdef ALLOW_DIAGNOSTICS
+       CHARACTER*8   diagName
+ #endif
+ #endif /* ALLOW_SITRACER */
+ #ifdef ALLOW_AUTODIFF_TAMC
+       INTEGER ilockey
+ #endif
+ C==   local arrays ==
  C--   TmixLoc        :: ocean surface/mixed-layer temperature (in K)
        _RL TmixLoc       (1:sNx,1:sNy)
+ C     actual ice thickness (with upper and lower limit)
+       _RL heffActual          (1:sNx,1:sNy)
+ C     actual snow thickness
+       _RL hsnowActual         (1:sNx,1:sNy)
+ C     actual ice thickness (with lower limit only) Reciprocal
+       _RL recip_heffActual    (1:sNx,1:sNy)
+ C     AREA_PRE :: hold sea-ice fraction field before any seaice-thermo update
+       _RL AREApreTH           (1:sNx,1:sNy)
+       _RL HEFFpreTH           (1:sNx,1:sNy)
+       _RL HSNWpreTH           (1:sNx,1:sNy)
+ #ifdef SEAICE_ITD
+       _RL AREAITDpreTH        (1:sNx,1:sNy,1:nITD)
+       _RL HEFFITDpreTH        (1:sNx,1:sNy,1:nITD)
+       _RL HSNWITDpreTH        (1:sNx,1:sNy,1:nITD)
+       _RL areaFracFactor      (1:sNx,1:sNy,1:nITD)
+       _RL leadIceThickMin
+ #endif
+ C     wind speed
+       _RL UG                  (1:sNx,1:sNy)
+ C     temporary variables available for the various computations
+       _RL tmparr1             (1:sNx,1:sNy)
+ #ifdef SEAICE_VARIABLE_SALINITY
+       _RL saltFluxAdjust      (1:sNx,1:sNy)
+ #endif
+       _RL ticeInMult          (1:sNx,1:sNy,MULTDIM)
+       _RL ticeOutMult         (1:sNx,1:sNy,MULTDIM)
+       _RL heffActualMult      (1:sNx,1:sNy,MULTDIM)
+       _RL hsnowActualMult     (1:sNx,1:sNy,MULTDIM)
+ #ifdef SEAICE_ITD
+       _RL recip_heffActualMult(1:sNx,1:sNy,MULTDIM)
+ #endif
+       _RL a_QbyATMmult_cover  (1:sNx,1:sNy,MULTDIM)
+       _RL a_QSWbyATMmult_cover(1:sNx,1:sNy,MULTDIM)
+       _RL a_FWbySublimMult    (1:sNx,1:sNy,MULTDIM)
+ #ifdef SEAICE_ITD
+       _RL r_QbyATMmult_cover  (1:sNx,1:sNy,MULTDIM)
+       _RL r_FWbySublimMult    (1:sNx,1:sNy,MULTDIM)
+ #endif
  C     a_QbyATM_cover :: available heat (in W/m^2) due to the interaction of
  C             the atmosphere and the ocean surface - for ice covered water
  C     a_QbyATM_open  :: same but for open water
-Line 122 
 C             processes have been accoun
+Line 217 
 C             processes have been accoun
        _RL a_QbyOCN            (1:sNx,1:sNy)
        _RL r_QbyOCN            (1:sNx,1:sNy)
- C conversion factors to go from Q (W/m2) to HEFF (ice meters)
+ C     The change of mean ice thickness due to out-of-bounds values following
-       _RL convertQ2HI, convertHI2Q
+ C     sea ice dyhnamics
- C conversion factors to go from precip (m/s) unit to HEFF (ice meters)
-       _RL convertPRECIP2HI, convertHI2PRECIP
- #ifdef ALLOW_DIAGNOSTICS
- C ICE/SNOW stocks tendencies associated with the various melt/freeze processes
-       _RL d_AREAbyATM         (1:sNx,1:sNy)
-       _RL d_AREAbyOCN         (1:sNx,1:sNy)
-       _RL d_AREAbyICE         (1:sNx,1:sNy)
- #endif
- #ifdef SEAICE_ALLOW_AREA_RELAXATION
- C     ICE/SNOW stocks tendency associated with relaxation towards observation
-       _RL d_AREAbyRLX         (1:sNx,1:sNy)
- c     The change of mean ice thickness due to relaxation
-       _RL d_HEFFbyRLX         (1:sNx,1:sNy)
- #endif
- c     The change of mean ice thickness due to out-of-bounds values following
- c     sea ice dynamics
        _RL d_HEFFbyNEG         (1:sNx,1:sNy)
- c     The change of mean ice thickness due to turbulent ocean-sea ice heat fluxes
+ C     The change of mean ice thickness due to turbulent ocean-sea ice heat fluxes
        _RL d_HEFFbyOCNonICE    (1:sNx,1:sNy)
- c     The sum of mean ice thickness increments due to atmospheric fluxes over the open water
+ C     The sum of mean ice thickness increments due to atmospheric fluxes over
- c     fraction and ice-covered fractions of the grid cell
+ C     the open water fraction and ice-covered fractions of the grid cell
        _RL d_HEFFbyATMonOCN    (1:sNx,1:sNy)
- c     The change of mean ice thickness due to flooding by snow
+ C     The change of mean ice thickness due to flooding by snow
        _RL d_HEFFbyFLOODING    (1:sNx,1:sNy)
- c     The mean ice thickness increments due to atmospheric fluxes over the open water
+ C     The mean ice thickness increments due to atmospheric fluxes over the open
- c     fraction and ice-covered fractions of the grid cell, respectively
+ C     water fraction and ice-covered fractions of the grid cell, respectively
        _RL d_HEFFbyATMonOCN_open(1:sNx,1:sNy)
        _RL d_HEFFbyATMonOCN_cover(1:sNx,1:sNy)
-Line 165 
 c     fraction and ice-covered fractions
+Line 241 
 c     fraction and ice-covered fractions
        _RL d_HSNWbyRAIN        (1:sNx,1:sNy)
        _RL d_HFRWbyRAIN        (1:sNx,1:sNy)
- C
  C     a_FWbySublim :: fresh water flux implied by latent heat of
  C                     sublimation to atmosphere, same sign convention
  C                     as EVAP (positive upward)
-Line 178 
 C                     as EVAP (positive
+Line 254 
 C                     as EVAP (positive
  C     The latent heat flux which will sublimate all snow and ice
  C     over one time step
        _RL latentHeatFluxMax   (1:sNx,1:sNy)
-       _RL latentHeatFluxMaxMult  (1:sNx,1:sNy,MULTDIM)
+       _RL latentHeatFluxMaxMult(1:sNx,1:sNy,MULTDIM)
  #endif
- C     actual ice thickness (with upper and lower limit)
+ #ifdef EXF_ALLOW_SEAICE_RELAX
-       _RL heffActual          (1:sNx,1:sNy)
+ C     ICE/SNOW stocks tendency associated with relaxation towards observation
- C     actual snow thickness
+       _RL d_AREAbyRLX         (1:sNx,1:sNy)
-       _RL hsnowActual         (1:sNx,1:sNy)
+ C     The change of mean ice thickness due to relaxation
- C     actual ice thickness (with lower limit only) Reciprocal
+       _RL d_HEFFbyRLX         (1:sNx,1:sNy)
-       _RL recip_heffActual    (1:sNx,1:sNy)
- C     local value (=1/HO or 1/HO_south)
-       _RL recip_HO
- C     local value (=1/ice thickness)
-       _RL recip_HH
- C     AREA_PRE :: hold sea-ice fraction field before any seaice-thermo update
-       _RL AREApreTH           (1:sNx,1:sNy)
-       _RL HEFFpreTH           (1:sNx,1:sNy)
-       _RL HSNWpreTH           (1:sNx,1:sNy)
- #ifdef SEAICE_ITD
-       _RL AREAITDpreTH        (1:sNx,1:sNy,1:nITD)
-       _RL HEFFITDpreTH        (1:sNx,1:sNy,1:nITD)
-       _RL HSNWITDpreTH        (1:sNx,1:sNy,1:nITD)
-       _RL areaFracFactor      (1:sNx,1:sNy,1:nITD)
-       _RL 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
-       _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
-       _RL pFac
-       _RL ticeInMult          (1:sNx,1:sNy,MULTDIM)
-       _RL ticeOutMult         (1:sNx,1:sNy,MULTDIM)
-       _RL heffActualMult      (1:sNx,1:sNy,MULTDIM)
- #ifdef SEAICE_ITD
-       _RL hsnowActualMult     (1:sNx,1:sNy,MULTDIM)
-       _RL recip_heffActualMult(1:sNx,1:sNy,MULTDIM)
- #endif
-       _RL a_QbyATMmult_cover  (1:sNx,1:sNy,MULTDIM)
-       _RL a_QSWbyATMmult_cover(1:sNx,1:sNy,MULTDIM)
-       _RL a_FWbySublimMult    (1:sNx,1:sNy,MULTDIM)
  #ifdef SEAICE_ITD
-       _RL r_QbyATMmult_cover  (1:sNx,1:sNy,MULTDIM)
+       _RL d_HEFFbySublim_ITD         (1:sNx,1:sNy,1:nITD)
-       _RL r_FWbySublimMult    (1:sNx,1:sNy,MULTDIM)
+       _RL d_HSNWbySublim_ITD         (1:sNx,1:sNy,1:nITD)
+       _RL d_HEFFbyOCNonICE_ITD       (1:sNx,1:sNy,1:nITD)
+       _RL d_HSNWbyATMonSNW_ITD       (1:sNx,1:sNy,1:nITD)
+       _RL d_HEFFbyATMonOCN_ITD       (1:sNx,1:sNy,1:nITD)
+       _RL d_HEFFbyATMonOCN_cover_ITD (1:sNx,1:sNy,1:nITD)
+       _RL d_HEFFbyATMonOCN_open_ITD  (1:sNx,1:sNy,1:nITD)
+       _RL d_HSNWbyRAIN_ITD           (1:sNx,1:sNy,1:nITD)
+       _RL d_HSNWbyOCNonSNW_ITD       (1:sNx,1:sNy,1:nITD)
+       _RL d_HEFFbyFLOODING_ITD       (1:sNx,1:sNy,1:nITD)
  #endif
- C     Helper variables: reciprocal of some constants
-       _RL recip_multDim
-       _RL recip_deltaTtherm
-       _RL recip_rhoIce
- C     Factor by which we increase the upper ocean friction velocity (u*) when
- C     ice is absent in a grid cell  (dimensionless)
-       _RL MixedLayerTurbulenceFactor
- #ifdef ALLOW_SITRACER
-       INTEGER iTr
-       CHARACTER*8   diagName
- #endif
  #ifdef ALLOW_DIAGNOSTICS
- c     Helper variables for diagnostics
+ C ICE/SNOW stocks tendencies associated with the various melt/freeze processes
+       _RL d_AREAbyATM         (1:sNx,1:sNy)
+       _RL d_AREAbyOCN         (1:sNx,1:sNy)
+       _RL d_AREAbyICE         (1:sNx,1:sNy)
+ C     Helper variables for diagnostics
        _RL DIAGarrayA    (1:sNx,1:sNy)
        _RL DIAGarrayB    (1:sNx,1:sNy)
        _RL DIAGarrayC    (1:sNx,1:sNy)
        _RL DIAGarrayD    (1:sNx,1:sNy)
- #endif
+ #endif /* ALLOW_DIAGNOSTICS */
+       _RL SItflux     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
+       _RL SIatmQnt    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
+       _RL SIatmFW     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
+ #ifdef ALLOW_BALANCE_FLUXES
+       _RL FWFsiTile(nSx,nSy)
+       _RL FWFsiGlob
+       _RL HFsiTile(nSx,nSy)
+       _RL HFsiGlob
+       _RL FWF2HFsiTile(nSx,nSy)
+       _RL FWF2HFsiGlob
+       CHARACTER*(max_len_mbuf) msgbuf
+ #endif
  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
-Line 333 
 C conversion factors to go from precip (
+Line 370 
 C conversion factors to go from precip (
       &                       + act4*max1*max2*max3
  #endif /* ALLOW_AUTODIFF_TAMC */
  C array initializations
  C =====================
-Line 356 
 C =====================
+Line 392 
 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
-Line 380 
 C =====================
+Line 416 
 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 398 
 c
+Line 431 
 c
              latentHeatFluxMaxMult(I,J,IT) = 0.0 _d 0
  #endif
  #ifdef SEAICE_ITD
-             r_QbyATMmult_cover (I,J,IT)   = 0.0 _d 0
+             d_HEFFbySublim_ITD(I,J,IT)         = 0.0 _d 0
-             r_FWbySublimMult(I,J,IT)      = 0.0 _d 0
+             d_HSNWbySublim_ITD(I,J,IT)         = 0.0 _d 0
+             d_HEFFbyOCNonICE_ITD(I,J,IT)       = 0.0 _d 0
+             d_HSNWbyATMonSNW_ITD(I,J,IT)       = 0.0 _d 0
+             d_HEFFbyATMonOCN_ITD(I,J,IT)       = 0.0 _d 0
+             d_HEFFbyATMonOCN_cover_ITD(I,J,IT) = 0.0 _d 0
+             d_HEFFbyATMonOCN_open_ITD(I,J,IT)  = 0.0 _d 0
+             d_HSNWbyRAIN_ITD(I,J,IT)           = 0.0 _d 0
+             d_HSNWbyOCNonSNW_ITD(I,J,IT)       = 0.0 _d 0
+             d_HEFFbyFLOODING_ITD(I,J,IT)       = 0.0 _d 0
+             r_QbyATMmult_cover(I,J,IT) = 0.0 _d 0
+             r_FWbySublimMult(I,J,IT)   = 0.0 _d 0
  #endif
            ENDDO
           ENDDO
-Line 412 
 c
+Line 455 
 c
          ENDDO
  #endif
  C =====================================================================
  C ===========PART 1: treat pathological cases (post advdiff)===========
  C =====================================================================
-Line 453 
 Cgf no dependency through pathological c
+Line 495 
 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 483 
 C           d_HEFFbyRLX(i,j) = 1. _d 1 *
+Line 525 
 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 492 
 CADJ STORE heff(:,:,bi,bj)  = comlev1_bi
+Line 534 
 CADJ STORE heff(:,:,bi,bj)  = comlev1_bi
  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
  CADJ STORE area(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
  #endif /* ALLOW_AUTODIFF_TAMC */
+ #ifdef SEAICE_ITD
+         DO IT=1,nITD
+ #endif
          DO J=1,sNy
           DO I=1,sNx
  #ifdef SEAICE_ITD
-           DO IT=1,nITD
+           tmpscal2=0. _d 0
-            tmpscal2=0. _d 0
+           tmpscal3=0. _d 0
-            tmpscal3=0. _d 0
+           tmpscal2=MAX(-HEFFITD(I,J,IT,bi,bj),0. _d 0)
-            tmpscal2=MAX(-HEFFITD(I,J,IT,bi,bj),0. _d 0)
+           HEFFITD(I,J,IT,bi,bj)=HEFFITD(I,J,IT,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
-            d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2
+           tmpscal3=MAX(-HSNOWITD(I,J,IT,bi,bj),0. _d 0)
-            tmpscal3=MAX(-HSNOWITD(I,J,IT,bi,bj),0. _d 0)
+           HSNOWITD(I,J,IT,bi,bj)=HSNOWITD(I,J,IT,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
-            d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3
+           AREAITD(I,J,IT,bi,bj)=MAX(AREAITD(I,J,IT,bi,bj),0. _d 0)
-            AREAITD(I,J,IT,bi,bj)=MAX(AREAITD(I,J,IT,bi,bj),0. _d 0)
-           ENDDO
  CToM      AREA, HEFF, and HSNOW will be updated at end of PART 1
  C         by calling SEAICE_ITD_SUM
  #else
-Line 517 
 C         by calling SEAICE_ITD_SUM
+Line 560 
 C         by calling SEAICE_ITD_SUM
  #endif
           ENDDO
          ENDDO
+ #ifdef SEAICE_ITD
+         ENDDO
+ #endif
  C 1.25) treat the case of very thin ice:
  #ifdef ALLOW_AUTODIFF_TAMC
  CADJ STORE heff(:,:,bi,bj)  = comlev1_bibj, key = iicekey,byte=isbyte
  #endif /* ALLOW_AUTODIFF_TAMC */
-         DO J=1,sNy
-          DO I=1,sNx
  #ifdef SEAICE_ITD
-           DO IT=1,nITD
+         DO IT=1,nITD
  #endif
+         DO J=1,sNy
+          DO I=1,sNx
            tmpscal2=0. _d 0
            tmpscal3=0. _d 0
  #ifdef SEAICE_ITD
-            IF (HEFFITD(I,J,IT,bi,bj).LE.siEps) THEN
+           IF (HEFFITD(I,J,IT,bi,bj).LE.siEps) THEN
-             tmpscal2=-HEFFITD(I,J,IT,bi,bj)
+            tmpscal2=-HEFFITD(I,J,IT,bi,bj)
-             tmpscal3=-HSNOWITD(I,J,IT,bi,bj)
+            tmpscal3=-HSNOWITD(I,J,IT,bi,bj)
-             TICES(I,J,IT,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
+ CToM       TICE will be updated at end of Part 1 together with AREA and HEFF
-            ENDIF
+           ENDIF
-            HEFFITD(I,J,IT,bi,bj) =HEFFITD(I,J,IT,bi,bj) +tmpscal2
+           HEFFITD(I,J,IT,bi,bj) =HEFFITD(I,J,IT,bi,bj) +tmpscal2
-            HSNOWITD(I,J,IT,bi,bj)=HSNOWITD(I,J,IT,bi,bj)+tmpscal3
+           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 553 
 CToM        TICE will be updated at end
+Line 599 
 CToM        TICE will be updated at end
  #endif
            d_HEFFbyNEG(I,J)=d_HEFFbyNEG(I,J)+tmpscal2
            d_HSNWbyNEG(I,J)=d_HSNWbyNEG(I,J)+tmpscal3
- #ifdef SEAICE_ITD
-           ENDDO
- #endif
           ENDDO
          ENDDO
+ #ifdef SEAICE_ITD
+         ENDDO
+ #endif
  C 1.5) treat the case of area but no ice/snow:
-Line 565 
 C 1.5) treat the case of area but no ice
+Line 611 
 C 1.5) treat the case of area but no ice
  CADJ STORE heff(:,:,bi,bj)  = comlev1_bibj, key = iicekey,byte=isbyte
  CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
  #endif /* ALLOW_AUTODIFF_TAMC */
+ #ifdef SEAICE_ITD
+         DO IT=1,nITD
+ #endif
          DO J=1,sNy
           DO I=1,sNx
  #ifdef SEAICE_ITD
-           DO IT=1,nITD
+           IF ((HEFFITD(I,J,IT,bi,bj).EQ.0. _d 0).AND.
-            IF ((HEFFITD(I,J,IT,bi,bj).EQ.0. _d 0).AND.
+      &        (HSNOWITD(I,J,IT,bi,bj).EQ.0. _d 0))
-      &         (HSNOWITD(I,J,IT,bi,bj).EQ.0. _d 0))
+      &     AREAITD(I,J,IT,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.
       &        (HSNOW(i,j,bi,bj).EQ.0. _d 0)) AREA(I,J,bi,bj)=0. _d 0
  #endif
           ENDDO
          ENDDO
+ #ifdef SEAICE_ITD
+         ENDDO
+ #endif
  C 2) treat the case of very small area:
-Line 586 
 C 2) treat the case of very small area:
+Line 636 
 C 2) treat the case of very small area:
  #ifdef ALLOW_AUTODIFF_TAMC
  CADJ STORE area(:,:,bi,bj)  = comlev1_bibj, key = iicekey,byte=isbyte
  #endif /* ALLOW_AUTODIFF_TAMC */
+ #ifdef SEAICE_ITD
+         DO IT=1,nITD
+ #endif
          DO J=1,sNy
           DO I=1,sNx
  #ifdef SEAICE_ITD
-           DO IT=1,nITD
+           IF ((HEFFITD(I,J,IT,bi,bj).GT.0).OR.
-            IF ((HEFFITD(I,J,IT,bi,bj).GT.0).OR.
+      &        (HSNOWITD(I,J,IT,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
- CToM        SEAICE_area_floor*nITD cannot be allowed to exceed 1
+ C          hence use SEAICE_area_floor devided by nITD
- C           hence use SEAICE_area_floor devided by nITD
+ C          (or install a warning in e.g. seaice_readparms.F)
- C           (or install a warning in e.g. seaice_readparms.F)
+            AREAITD(I,J,IT,bi,bj)=
-             AREAITD(I,J,IT,bi,bj)=
+      &        MAX(AREAITD(I,J,IT,bi,bj),SEAICE_area_floor/float(nITD))
-      &         MAX(AREAITD(I,J,IT,bi,bj),SEAICE_area_floor/float(nITD))
+           ENDIF
-            ENDIF
-           ENDDO
  #else
            IF ((HEFF(i,j,bi,bj).GT.0).OR.(HSNOW(i,j,bi,bj).GT.0)) THEN
             AREA(I,J,bi,bj)=MAX(AREA(I,J,bi,bj),SEAICE_area_floor)
-Line 606 
 C           (or install a warning in e.g
+Line 657 
 C           (or install a warning in e.g
  #endif
           ENDDO
          ENDDO
+ #ifdef SEAICE_ITD
+         ENDDO
+ #endif
  #endif /* DISABLE_AREA_FLOOR */
  C 2.5) treat case of excessive ice cover, e.g., due to ridging:
-Line 631 
 CADJ STORE area(:,:,bi,bj)  = comlev1_bi
+Line 685 
 CADJ STORE area(:,:,bi,bj)  = comlev1_bi
  #ifdef SEAICE_ITD
  CToM catch up with items 1.25 and 2.5 involving category sums AREA and HEFF
-         DO J=1,sNy
+         DO IT=1,nITD
-          DO I=1,sNx
+          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
+            IF (IT .eq. 1) THEN
-           tmpscal2 = 0. _d 0
+             tmpscal1itd(i,j) = 0. _d 0
-           tmpscal3 = 0. _d 0
+             tmpscal2itd(i,j) = 0. _d 0
-           DO IT=1,nITD
+             tmpscal3itd(i,j) = 0. _d 0
-            tmpscal1=tmpscal1 + TICES(I,J,IT,bi,bj)*HEFFITD(I,J,IT,bi,bj)
+            ENDIF
-            tmpscal2=tmpscal2 + HEFFITD(I,J,IT,bi,bj)
+            tmpscal1itd(i,j)=tmpscal1itd(i,j) + TICES(I,J,IT,bi,bj)
-            tmpscal3=tmpscal3 + AREAITD(I,J,IT,bi,bj)
+      &                                       * HEFFITD(I,J,IT,bi,bj)
-           ENDDO
+            tmpscal2itd(i,j)=tmpscal2itd(i,j) + HEFFITD(I,J,IT,bi,bj)
-           TICE(I,J,bi,bj)=tmpscal1/tmpscal2
+            tmpscal3itd(i,j)=tmpscal3itd(i,j) + AREAITD(I,J,IT,bi,bj)
+            IF (IT .eq. nITD) THEN
+             TICE(I,J,bi,bj)=tmpscal1itd(i,j)/tmpscal2itd(i,j)
  C    lines of item 2.5 that were omitted:
  C    in 2.5 these lines are executed before "ridging" is applied to AREA
  C    hence we execute them here before SEAICE_ITD_REDIST is called
  C    although this means that AREA has not been completely regularized
  #ifdef ALLOW_DIAGNOSTICS
-           DIAGarrayA(I,J) = tmpscal3
+             DIAGarrayA(I,J) = tmpscal3itd(i,j)
  #endif
  #ifdef ALLOW_SITRACER
-           SItrAREA(I,J,bi,bj,1)=tmpscal3
+             SItrAREA(I,J,bi,bj,1)=tmpscal3itd(i,j)
  #endif
+            ENDIF
+           ENDDO
           ENDDO
          ENDDO
-Line 663 
 C    which includes ridging as in item 2
+Line 722 
 C    which includes ridging as in item 2
  C    and update AREA, HEFF, and HSNOW
          CALL SEAICE_ITD_REDIST(bi, bj, myTime, myIter, myThid)
          CALL SEAICE_ITD_SUM(bi, bj, myTime, myIter, myThid)
+ #endif /* SEAICE_ITD */
- c ToM<<< debug seaice_growth
+ #ifdef SEAICE_DEBUG
-         WRITE(msgBuf,'(A,7F8.4)')
+ #ifdef SEAICE_ITD
-      &    ' SEAICE_GROWTH: Heff increments 0, HEFFITD = ',
+         WRITE(msgBufForm,'(A,I2,A)') '(A,',nITD,'F14.10)'
-      &     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)
  #else
-         WRITE(msgBuf,'(A,7F8.4)')
+         WRITE(msgBufForm,'(A,I2,A)') '(A,  F14.10)'
+ #endif
+         WRITE(msgBuf,msgBufForm)
       &    ' SEAICE_GROWTH: Heff increments 0, HEFF = ',
+ #ifdef SEAICE_ITD
+      &     HEFFITD(1,1,:,bi,bj)
+ #else
       &     HEFF(1,1,bi,bj)
+ #endif
          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
       &    SQUEEZE_RIGHT , myThid)
-         WRITE(msgBuf,'(A,7F8.4)')
+         WRITE(msgBuf,msgBufForm)
       &    ' SEAICE_GROWTH: Area increments 0, AREA = ',
+ #ifdef SEAICE_ITD
+      &     AREAITD(1,1,:,bi,bj)
+ #else
       &     AREA(1,1,bi,bj)
+ #endif
          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
       &    SQUEEZE_RIGHT , myThid)
- c ToM>>>
+ #endif
- #endif /* SEAICE_ITD */
  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
  C        end SEAICEadjMODE.EQ.0 statement:
          ENDIF
-Line 732 
 C           if there's no ice, potential
+Line 794 
 C           if there's no ice, potential
            ENDDO
           ENDDO
          ENDDO
+ #ifdef ALLOW_SITRACER
  C prepare SItrHEFF to be computed as cumulative sum
          DO iTr=2,5
           DO J=1,sNy
-Line 747 
 C prepare SItrAREA to be computed as cum
+Line 810 
 C prepare SItrAREA to be computed as cum
           ENDDO
          ENDDO
  #endif
+ #endif /* SEAICE_ITD */
  C 4) treat sea ice salinity pathological cases
  #ifdef SEAICE_VARIABLE_SALINITY
-Line 875 
 cif       Do not regularize when HEFFpre
+Line 939 
 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 915 
 C    AND SNOW THICKNESS
+Line 979 
 C    AND SNOW THICKNESS
  #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,IT) .GT. ZERO) THEN
              latentHeatFluxMaxMult(I,J,IT) = lhSublim*recip_deltaTtherm *
-Line 971 
 cCADJ STORE TmixLoc = comlev1_bibj, key
+Line 1035 
 cCADJ STORE TmixLoc = comlev1_bibj, key
       I       TmixLoc,
       O       a_QbyATM_open, a_QSWbyATM_open,
       I       bi, bj, myTime, myIter, myThid )
- c ToM<<< debugging
-         print*,' '
-         print*,'UG = ',UG(1,1)
-         print*,'Tsurf = ',TmixLoc(1,1)
-         print*,'a_QbyATM_open = ',a_QbyATM_open(1,1)
-         print*,' '
- c ToM>>>
  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 1004 
 C     Compute relative wind speed over s
+Line 1060 
 C     Compute relative wind speed over s
            ENDDO
           ENDDO
          ENDIF
- #endif /* ALLOW_ATM_WIND */
  #ifdef ALLOW_AUTODIFF_TAMC
  CADJ STORE tice(:,:,bi,bj)
-Line 1020 
 CADJ &                       key = iicek
+Line 1075 
 CADJ &                       key = iicek
  C--   Start loop over multi-categories
  #ifdef SEAICE_ITD
- CToM  SEAICE_multDim = nITD (see SEAICE_SIZE.h and seaice_readparms.F)
+         DO IT=1,nITD
- #endif
-         DO IT=1,SEAICE_multDim
- 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
- #endif
           DO J=1,sNy
            DO I=1,sNx
- #ifndef SEAICE_ITD
  CToM for SEAICE_ITD heffActualMult and latentHeatFluxMaxMult are calculated above
  C    (instead of heffActual and latentHeatFluxMax)
+            ticeInMult(I,J,IT)=TICES(I,J,IT,bi,bj)
+            ticeOutMult(I,J,IT)=TICES(I,J,IT,bi,bj)
+            TICE(I,J,bi,bj) = ZERO
+            TICES(I,J,IT,bi,bj) = ZERO
+           ENDDO
+          ENDDO
+         ENDDO
+ #else
+         DO IT=1,SEAICE_multDim
+ C        homogeneous distribution between 0 and 2 x heffActual
+          pFac = (2.0 _d 0*IT - 1.0 _d 0)*recip_multDim
+          pFacSnow = 1. _d 0
+          IF ( SEAICE_useMultDimSnow ) pFacSnow=pFac
+          DO J=1,sNy
+           DO I=1,sNx
             heffActualMult(I,J,IT)= heffActual(I,J)*pFac
+            hsnowActualMult(I,J,IT)=hsnowActual(I,J)*pFacSnow
  #ifdef SEAICE_CAP_SUBLIM
             latentHeatFluxMaxMult(I,J,IT) = latentHeatFluxMax(I,J)*pFac
  #endif
- #endif
             ticeInMult(I,J,IT)=TICES(I,J,IT,bi,bj)
             ticeOutMult(I,J,IT)=TICES(I,J,IT,bi,bj)
             TICE(I,J,bi,bj) = ZERO
-Line 1044 
 C    (instead of heffActual and latentHe
+Line 1107 
 C    (instead of heffActual and latentHe
            ENDDO
           ENDDO
          ENDDO
+ #endif
  #ifdef ALLOW_AUTODIFF_TAMC
  CADJ STORE heffActualMult = comlev1_bibj, key = iicekey, byte = isbyte
+ CADJ STORE hsnowActualMult= comlev1_bibj, key = iicekey, byte = isbyte
  CADJ STORE ticeInMult     = comlev1_bibj, key = iicekey, byte = isbyte
  # ifdef SEAICE_CAP_SUBLIM
  CADJ STORE latentHeatFluxMaxMult
-Line 1062 
 CADJ &     comlev1_bibj, key = iicekey,
+Line 1127 
 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 1130 
 C     (fluxes are per unit (ice surface)
+Line 1193 
 C     (fluxes are per unit (ice surface)
            ENDDO
           ENDDO
          ENDDO
- c ToM<<< debugging
-         print*,' '
-         print*,'after SOLVE4TEMP: '
-         print*,'TICE  = ',TICE(1,1,bi,bj)
-         print*,'TICES = ',TICES(1,1,:,bi,bj)
-         print*,'a_QSWbyATM_cover     = ',a_QSWbyATM_cover(1,1)
-         print*,'a_QSWbyATMmult_cover = ',a_QSWbyATMmult_cover(1,1,:)
-         print*,' '
- c ToM>>>
  #ifdef SEAICE_CAP_SUBLIM
  # 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 1216 
 C and initialize r_QbyATM_cover/r_QbyATM
+Line 1270 
 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 1276 
 CADJ &                       key = iicek
+Line 1330 
 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 1296 
 c efficiency of turbulent fluxes : depen
+Line 1350 
 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)
- c ToM<<< debugging
-           if (i.eq.1 .and. j.eq.1) then
-             print *, 'salt [psu]     = ',salt(i,j,kSurface,bi,bj)
-             print *, 'theta [degC]   = ',theta(i,j,kSurface,bi,bj)
-             print *, 'tempFrz [degC] = ',tempFrz
-             print *, 'max turb flx [m]   = ',tmpscal2
-             print *, 'avail trub flx [m] = ',a_QbyOCN(i,j)
-           endif
- c ToM>>>
           ENDDO
          ENDDO
-Line 1320 
 c ToM>>>
+Line 1365 
 c ToM>>>
          CALL ZERO_ADJ_1D( sNx*sNy, r_QbyOCN, myThid)
  #endif
  C ===================================================================
  C =========PART 3: determine effective thicknesses increments========
  C ===================================================================
-Line 1342 
 C     First sublimate/deposite snow
+Line 1386 
 C     First sublimate/deposite snow
  #ifdef SEAICE_ITD
       &     MAX(MIN(r_FWbySublimMult(I,J,IT),HSNOWITD(I,J,IT,bi,bj)
       &             *SNOW2ICE),ZERO)
+           d_HSNWbySublim_ITD(I,J,IT) = - tmpscal2 * ICE2SNOW
  C         accumulate change over ITD categories
            d_HSNWbySublim(I,J)     = d_HSNWbySublim(I,J)      - tmpscal2
       &                                                       *ICE2SNOW
-           HSNOWITD(I,J,IT,bi,bj)  = HSNOWITD(I,J,IT,bi,bj)   - tmpscal2
-      &                                                       *ICE2SNOW
            r_FWbySublimMult(I,J,IT)= r_FWbySublimMult(I,J,IT) - tmpscal2
  #else
       &     MAX(MIN(r_FWbySublim(I,J),HSNOW(I,J,bi,bj)*SNOW2ICE),ZERO)
-Line 1366 
 C     If anything is left, sublimate ice
+Line 1409 
 C     If anything is left, sublimate ice
            tmpscal2 =
  #ifdef SEAICE_ITD
       &     MAX(MIN(r_FWbySublimMult(I,J,IT),HEFFITD(I,J,IT,bi,bj)),ZERO)
+           d_HEFFbySublim_ITD(I,J,IT) = - tmpscal2
  C         accumulate change over ITD categories
-           d_HSNWbySublim(I,J)      = d_HSNWbySublim(I,J)      - tmpscal2
+           d_HEFFbySublim(I,J)      = d_HEFFbySublim(I,J)      - tmpscal2
-           HEFFITD(I,J,IT,bi,bj)    = HEFFITD(I,J,IT,bi,bj)    - tmpscal2
            r_FWbySublimMult(I,J,IT) = r_FWbySublimMult(I,J,IT) - tmpscal2
  #else
       &     MAX(MIN(r_FWbySublim(I,J),HEFF(I,J,bi,bj)),ZERO)
-Line 1398 
 C     remove the fusion part for the res
+Line 1441 
 C     remove the fusion part for the res
  C       end IT loop
          ENDDO
  #endif
+ #ifdef SEAICE_DEBUG
  c ToM<<< debug seaice_growth
-         WRITE(msgBuf,'(A,7F8.4)')
+         WRITE(msgBuf,msgBufForm)
+      &    ' SEAICE_GROWTH: Hsnow increments 1, d_HSNWySublim = ',
  #ifdef SEAICE_ITD
-      &    ' SEAICE_GROWTH: Heff increments 1, HEFFITD = ',
+      &     d_HSNWbySublim_ITD(1,1,:)
-      &     HEFFITD(1,1,:,bi,bj)
+ #else
+      &     d_HSNWbySublim(1,1)
+ #endif
          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
       &    SQUEEZE_RIGHT , myThid)
-         WRITE(msgBuf,'(A,7F8.4)')
+         WRITE(msgBuf,msgBufForm)
-      &    ' SEAICE_GROWTH: Area increments 1, AREAITD = ',
+      &    ' SEAICE_GROWTH: Heff increments 1, d_HEFFbySublim = ',
-      &     AREAITD(1,1,:,bi,bj)
+ #ifdef SEAICE_ITD
+      &     d_HEFFbySublim_ITD(1,1,:)
  #else
-      &    ' SEAICE_GROWTH: Heff increments 1, HEFF = ',
+      &     d_HEFFbySublim(1,1)
-      &     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 1433 
 C          equally distributed under the
+Line 1481 
 C          equally distributed under the
  C          fractional area of each thickness category
             tmpscal1=MAX(r_QbyOCN(i,j)*areaFracFactor(I,J,IT),
       &                               -HEFFITD(I,J,IT,bi,bj))
+            d_HEFFbyOCNonICE_ITD(I,J,IT)=tmpscal1
             d_HEFFbyOCNonICE(I,J) = d_HEFFbyOCNonICE(I,J) + 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,IT,bi,bj)
- #endif
            ENDDO
           ENDDO
          ENDDO
+ #ifdef ALLOW_SITRACER
+         DO J=1,sNy
+          DO I=1,sNx
+           SItrHEFF(I,J,bi,bj,2) = HEFFpreTH(I,J)
+      &                          + d_HEFFbySublim(I,J)
+      &                          + d_HEFFbyOCNonICE(I,J)
+          ENDDO
+         ENDDO
+ #endif
          DO J=1,sNy
           DO I=1,sNx
            r_QbyOCN(I,J)=r_QbyOCN(I,J)-d_HEFFbyOCNonICE(I,J)
-Line 1459 
 C          fractional area of each thick
+Line 1512 
 C          fractional area of each thick
           ENDDO
          ENDDO
  #endif /* SEAICE_ITD */
+ #ifdef SEAICE_DEBUG
  c ToM<<< debug seaice_growth
-         WRITE(msgBuf,'(A,7F8.4)')
+         WRITE(msgBuf,msgBufForm)
+      &    ' SEAICE_GROWTH: Heff increments 2, d_HEFFbyOCNonICE = ',
  #ifdef SEAICE_ITD
-      &    ' SEAICE_GROWTH: Heff increments 2, HEFFITD = ',
+      &     d_HEFFbyOCNonICE_ITD(1,1,:)
-      &     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 = ',
+      &     d_HEFFbyOCNonICE(1,1)
-      &     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 1498 
 C     of snow. This appears to be more r
+Line 1547 
 C     of snow. This appears to be more r
  Cgf no additional dependency through snow
             IF ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
  #endif
+            d_HSNWbyATMonSNW_ITD(I,J,IT) = tmpscal2*ICE2SNOW
             d_HSNWbyATMonSNW(I,J) = d_HSNWbyATMonSNW(I,J)
       &                           + tmpscal2*ICE2SNOW
-            HSNOWITD(I,J,IT,bi,bj)= HSNOWITD(I,J,IT,bi,bj)
-      &                           + tmpscal2*ICE2SNOW
             r_QbyATMmult_cover(I,J,IT)=r_QbyATMmult_cover(I,J,IT)
       &                           - tmpscal2
            ENDDO
-Line 1524 
 Cgf no additional dependency through sno
+Line 1572 
 Cgf no additional dependency through sno
           ENDDO
          ENDDO
  #endif /* SEAICE_ITD */
+ #ifdef SEAICE_DEBUG
  c ToM<<< debug seaice_growth
-         WRITE(msgBuf,'(A,7F8.4)')
+         WRITE(msgBuf,msgBufForm)
+      &    ' SEAICE_GROWTH: Hsnow increments 3, d_HSNWbyATMonSNW = ',
  #ifdef SEAICE_ITD
-      &    ' SEAICE_GROWTH: Heff increments 3, HEFFITD = ',
+      &    d_HSNWbyATMonSNW_ITD(1,1,:)
-      &     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 = ',
+      &    d_HSNWbyATMonSNW(1,1)
-      &     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 1559 
 Cgf warming conditions, the lab_sea resu
+Line 1603 
 Cgf warming conditions, the lab_sea resu
          DO IT=1,nITD
           DO J=1,sNy
            DO I=1,sNx
+            tmpscal1 = HEFFITDpreTH(I,J,IT)
+      &              + d_HEFFbySublim_ITD(I,J,IT)
+      &              + d_HEFFbyOCNonICE_ITD(I,J,IT)
  #ifdef SEAICE_GROWTH_LEGACY
-            tmpscal2 = MAX(-HEFFITD(I,J,IT,bi,bj),
+            tmpscal2 = MAX(-tmpscal1,
       &                     r_QbyATMmult_cover(I,J,IT))
  #else
-            tmpscal2 = MAX(-HEFFITD(I,J,IT,bi,bj),
+            tmpscal2 = MAX(-tmpscal1,
       &                     r_QbyATMmult_cover(I,J,IT)
  c         Limit ice growth by potential melt by ocean
       &              + AREAITDpreTH(I,J,IT) * r_QbyOCN(I,J))
  #endif /* SEAICE_GROWTH_LEGACY */
+            d_HEFFbyATMonOCN_cover_ITD(I,J,IT) = tmpscal2
             d_HEFFbyATMonOCN_cover(I,J) = d_HEFFbyATMonOCN_cover(I,J)
       &                                 + tmpscal2
+            d_HEFFbyATMonOCN_ITD(I,J,IT) = d_HEFFbyATMonOCN_ITD(I,J,IT)
+      &                                 + tmpscal2
             d_HEFFbyATMonOCN(I,J)       = d_HEFFbyATMonOCN(I,J)
       &                                 + tmpscal2
             r_QbyATMmult_cover(I,J,IT)  = r_QbyATMmult_cover(I,J,IT)
       &                                 - 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,IT,bi,bj)
- #endif
            ENDDO
           ENDDO
          ENDDO
+ #ifdef ALLOW_SITRACER
+         DO J=1,sNy
+          DO I=1,sNx
+           SItrHEFF(I,J,bi,bj,3) = SItrHEFF(I,J,bi,bj,2)
+      &                          + d_HEFFbyATMonOCN_cover(I,J)
+          ENDDO
+         ENDDO
+ #endif
  #else /* SEAICE_ITD */
          DO J=1,sNy
           DO I=1,sNx
-Line 1591 
 c         Limit ice growth by potential
+Line 1643 
 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 1603 
 c         Limit ice growth by potential
+Line 1655 
 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,7F8.4)')
+         WRITE(msgBuf,msgBufForm)
+      &   ' SEAICE_GROWTH: Heff increments 4, d_HEFFbyATMonOCN_cover = ',
  #ifdef SEAICE_ITD
-      &    ' SEAICE_GROWTH: Heff increments 4, HEFFITD = ',
+      &     d_HEFFbyATMonOCN_cover_ITD(1,1,:)
-      &     HEFFITD(1,1,:,bi,bj)
+ #else
+      &     d_HEFFbyATMonOCN_cover(1,1)
+ #endif
          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
       &    SQUEEZE_RIGHT , myThid)
-         WRITE(msgBuf,'(A,7F8.4)')
+         WRITE(msgBuf,msgBufForm)
-      &    ' SEAICE_GROWTH: Area increments 4, AREAITD = ',
+      &    ' SEAICE_GROWTH: Heff increments 4, d_HEFFbyATMonOCN = ',
-      &     AREAITD(1,1,:,bi,bj)
+ #ifdef SEAICE_ITD
+      &     d_HEFFbyATMonOCN_ITD(1,1,:)
  #else
-      &    ' SEAICE_GROWTH: Heff increments 4, HEFF = ',
+      &     d_HEFFbyATMonOCN(1,1)
-      &     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 IT=1,nITD
           DO J=1,sNy
            DO I=1,sNx
-            HSNOWITD(I,J,IT,bi,bj) = HSNOWITD(I,J,IT,bi,bj)
+            d_HSNWbyRAIN_ITD(I,J,IT)
-      &     + d_HSNWbyRAIN(I,J)*areaFracFactor(I,J,IT)
+      &     = d_HSNWbyRAIN(I,J)*areaFracFactor(I,J,IT)
            ENDDO
           ENDDO
          ENDDO
-Line 1670 
 C           add precip to the fresh wate
+Line 1741 
 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,7F8.4)')
+         WRITE(msgBuf,msgBufForm)
+      &    ' SEAICE_GROWTH: Hsnow increments 5, d_HSNWbyRAIN = ',
  #ifdef SEAICE_ITD
-      &    ' SEAICE_GROWTH: Heff increments 5, HEFFITD = ',
+      &     d_HSNWbyRAIN_ITD(1,1,:)
-      &     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 = ',
+      &     d_HSNWbyRAIN(1,1)
-      &     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 1704 
 CADJ STORE r_QbyOCN = comlev1_bibj,key=i
+Line 1773 
 CADJ STORE r_QbyOCN = comlev1_bibj,key=i
          DO IT=1,nITD
           DO J=1,sNy
            DO I=1,sNx
+            tmpscal4 = HSNWITDpreTH(I,J,IT)
+      &              + d_HSNWbySublim_ITD(I,J,IT)
+      &              + d_HSNWbyATMonSNW_ITD(I,J,IT)
+      &              + d_HSNWbyRAIN_ITD(I,J,IT)
             tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW*areaFracFactor(I,J,IT),
-      &                  -HSNOWITD(I,J,IT,bi,bj))
+      &                  -tmpscal4)
             tmpscal2=MIN(tmpscal1,0. _d 0)
  #ifdef SEAICE_MODIFY_GROWTH_ADJ
  Cgf no additional dependency through snow
             if ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
  #endif
+            d_HSNWbyOCNonSNW_ITD(I,J,IT) = tmpscal2
             d_HSNWbyOCNonSNW(I,J) = d_HSNWbyOCNonSNW(I,J) + tmpscal2
             r_QbyOCN(I,J)=r_QbyOCN(I,J) - tmpscal2*SNOW2ICE
-            HSNOWITD(I,J,IT,bi,bj) = HSNOWITD(I,J,IT,bi,bj) + tmpscal2
            ENDDO
           ENDDO
          ENDDO
-Line 1735 
 Cgf no additional dependency through sno
+Line 1808 
 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,7F8.4)')
+         WRITE(msgBuf,msgBufForm)
+      &    ' SEAICE_GROWTH: Hsnow increments 6, d_HSNWbyOCNonSNW = ',
  #ifdef SEAICE_ITD
-      &    ' SEAICE_GROWTH: Heff increments 6, HEFFITD = ',
+      &     d_HSNWbyOCNonSNW_ITD(1,1,:)
-      &     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 = ',
+      &     d_HSNWbyOCNonSNW(1,1)
-      &     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 1765 
 CADJ STORE a_QSWbyATM_open = comlev1_bib
+Line 1834 
 CADJ STORE a_QSWbyATM_open = comlev1_bib
          DO J=1,sNy
           DO I=1,sNx
- c           Initial ice growth is triggered by open water
+ #ifdef SEAICE_ITD
- c           heat flux overcoming potential melt by ocean
+ C         HEFF will be updated at the end of PART 3,
+ C         hence sum of tendencies so far is needed
+           tmpscal4 = HEFFpreTH(I,J)
+      &             + d_HEFFbySublim(I,J)
+      &             + d_HEFFbyOCNonICE(I,J)
+      &             + d_HEFFbyATMonOCN(I,J)
+ #else
+ C         HEFF is updated step by step throughout seaice_growth
+           tmpscal4 = HEFF(I,J,bi,bj)
+ #endif
+ C           Initial ice growth is triggered by open water
+ C           heat flux overcoming potential melt by ocean
              tmpscal1=r_QbyATM_open(I,J)+r_QbyOCN(i,j) *
       &            (1.0 _d 0 - AREApreTH(I,J))
- c           Penetrative shortwave flux beyond first layer
+ C           Penetrative shortwave flux beyond first layer
- c           that is therefore not available to ice growth/melt
+ C           that is therefore not available to ice growth/melt
              tmpscal2=SWFracB * a_QSWbyATM_open(I,J)
  C           impose -HEFF as the maxmum melting if SEAICE_doOpenWaterMelt
  C           or 0. otherwise (no melting if not SEAICE_doOpenWaterMelt)
              tmpscal3=facOpenGrow*MAX(tmpscal1-tmpscal2,
-      &           -HEFF(I,J,bi,bj)*facOpenMelt)*HEFFM(I,J,bi,bj)
+      &           -tmpscal4*facOpenMelt)*HEFFM(I,J,bi,bj)
- c ToM<<< debugging
+ #ifdef SEAICE_ITD
-             if (I.eq.1 .and. J.eq.1) then
+ C         ice growth in open water adds to first category
-              print*,'r_QbyATM_open(I,J) = ', r_QbyATM_open(I,J)
+           d_HEFFbyATMonOCN_open_ITD(I,J,1)=tmpscal3
-              print*,'r_QbyOCN(i,j) = ', r_QbyOCN(i,j)
+           d_HEFFbyATMonOCN_ITD(I,J,1)     =d_HEFFbyATMonOCN_ITD(I,J,1)
-              print*,'1 - AREApreTH = ', (1.0 _d 0 - AREApreTH(I,J))
+      &                                    +tmpscal3
-              print*,'tmpscal1 = ', tmpscal1
+ #endif
-              print*,' '
-              print*,'SWFracB = ', SWFracB
-              print*,'a_QSWbyATM_open(I,J) = ', a_QSWbyATM_open(I,J)
-              print*,'tmpscal2 = ', tmpscal2
-              print*,' '
-              print*,'facOpenGrow = ', facOpenGrow
-              print*,'HEFF(I,J,bi,bj) = ', HEFF(I,J,bi,bj)
-              print*,'facOpenMelt = ', facOpenMelt
-              print*,'MAX = ', MAX(tmpscal1-tmpscal2,
-      &              -HEFF(I,J,bi,bj)*facOpenMelt)
-              print*,'HEFFM(I,J,bi,bj) = ', HEFFM(I,J,bi,bj)
-              print*,'tmpscal3 = ', tmpscal3
-              print*,' '
-             endif
- c ToM>>>
            d_HEFFbyATMonOCN_open(I,J)=tmpscal3
            d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal3
            r_QbyATM_open(I,J)=r_QbyATM_open(I,J)-tmpscal3
- #ifdef SEAICE_ITD
- cC         open water area fraction
- c          tmpscal0 = ONE-AREApreTH(I,J)
- cC         determine thickness of new ice
- cctomC         considering the entire open water area to refreeze
- cctom          tmpscal1 = tmpscal3/tmpscal0
- cC         considering a minimum lead ice thickness of 10 cm
- cC         WATCH that leadIceThickMin is smaller that Hlimit(1)!
- c          leadIceThickMin = 0.1
- c          tmpscal1 = MAX(leadIceThickMin,tmpscal3/tmpscal0)
- cC         adjust ice area fraction covered by new ice
- c         tmpscal0 = tmpscal3/tmpscal1
- cC         then add new ice volume to appropriate thickness category
- c          DO IT=1,nITD
- c          IF (tmpscal1.LT.Hlimit(IT)) THEN
- c            HEFFITD(I,J,IT,bi,bj) = HEFFITD(I,J,IT,bi,bj) + tmpscal3
- c           tmpscal3=ZERO
- cC not sure if AREAITD should be changed here since AREA is incremented
- cC   in PART 4 below in non-itd code
- cC in this scenario all open water is covered by new ice instantaneously,
- cC   i.e. no delayed lead closing is concidered such as is associated with
- cC   Hibler's h_0 parameter
- c           AREAITD(I,J,IT,bi,bj) = AREAITD(I,J,IT,bi,bj)
- c     &                          + tmpscal0
- c           tmpscal0=ZERO
- c           ENDIF
- c          ENDDO
- ctom debugging: 1 category only
-           HEFFITD(I,J,1,bi,bj) = HEFFITD(I,J,1,bi,bj) + tmpscal3
- #else
            HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj) + tmpscal3
- #endif
           ENDDO
          ENDDO
  #ifdef ALLOW_SITRACER
- #ifdef SEAICE_ITD
-         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,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
+ #ifdef SEAICE_ITD
+           SItrHEFF(I,J,bi,bj,4)=SItrHEFF(I,J,bi,bj,3)
+      &                         +d_HEFFbyATMonOCN_open(I,J)
+ #else
            SItrHEFF(I,J,bi,bj,4)=HEFF(I,J,bi,bj)
+ #endif
           ENDDO
          ENDDO
- #endif
  #endif /* ALLOW_SITRACER */
+ #ifdef SEAICE_DEBUG
  c ToM<<< debug seaice_growth
-         WRITE(msgBuf,'(A,7F8.4)')
+         WRITE(msgBuf,msgBufForm)
+      &    ' SEAICE_GROWTH: Heff increments 7, d_HEFFbyATMonOCN_open = ',
  #ifdef SEAICE_ITD
-      &    ' SEAICE_GROWTH: Heff increments 7, HEFFITD = ',
+      &     d_HEFFbyATMonOCN_open_ITD(1,1,:)
-      &     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 = ',
+      &     d_HEFFbyATMonOCN_open(1,1)
-      &     HEFF(1,1,bi,bj)
+ #endif
          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
       &    SQUEEZE_RIGHT , myThid)
-         WRITE(msgBuf,'(A,7F8.4)')
+         WRITE(msgBuf,msgBufForm)
-      &    ' SEAICE_GROWTH: Area increments 7, AREA = ',
+      &    ' SEAICE_GROWTH: Heff increments 7, d_HEFFbyATMonOCN = ',
-      &     AREA(1,1,bi,bj)
+ #ifdef SEAICE_ITD
+      &     d_HEFFbyATMonOCN_ITD(1,1,:)
+ #else
+      &     d_HEFFbyATMonOCN(1,1)
  #endif
          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
       &    SQUEEZE_RIGHT , myThid)
  c ToM>>>
+ #endif
  C convert snow to ice if submerged.
  C =================================
-Line 1892 
 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
+Line 1919 
 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
           DO IT=1,nITD
            DO J=1,sNy
             DO I=1,sNx
-             tmpscal0 = (HSNOWITD(I,J,IT,bi,bj)*SEAICE_rhoSnow
+             tmpscal3 = HEFFITDpreTH(I,J,IT)
-      &               +  HEFFITD(I,J,IT,bi,bj) *SEAICE_rhoIce)
+      &               + d_HEFFbySublim_ITD(I,J,IT)
-      &                                        *recip_rhoConst
+      &               + d_HEFFbyOCNonICE_ITD(I,J,IT)
-             tmpscal1 = MAX( 0. _d 0, tmpscal0 - HEFFITD(I,J,IT,bi,bj))
+      &               + d_HEFFbyATMonOCN_ITD(I,J,IT)
+             tmpscal4 = HSNWITDpreTH(I,J,IT)
+      &               + d_HSNWbySublim_ITD(I,J,IT)
+      &               + d_HSNWbyATMonSNW_ITD(I,J,IT)
+      &               + d_HSNWbyRAIN_ITD(I,J,IT)
+             tmpscal0 = (tmpscal4*SEAICE_rhoSnow
+      &               +  tmpscal3*SEAICE_rhoIce)
+      &               * recip_rhoConst
+             tmpscal1 = MAX( 0. _d 0, tmpscal0 - tmpscal3)
+             d_HEFFbyFLOODING_ITD(I,J,IT) = tmpscal1
              d_HEFFbyFLOODING(I,J) = d_HEFFbyFLOODING(I,J)  + tmpscal1
-             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
           ENDDO
-Line 1913 
 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
+Line 1946 
 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,7F8.4)')
+         WRITE(msgBuf,msgBufForm)
+      &    ' SEAICE_GROWTH: Heff increments 8, d_HEFFbyFLOODING = ',
  #ifdef SEAICE_ITD
-      &    ' SEAICE_GROWTH: Heff increments 8, HEFFITD = ',
+      &     d_HEFFbyFLOODING_ITD(1,1,:)
-      &     HEFFITD(1,1,:,bi,bj)
+ #else
+      &     d_HEFFbyFLOODING(1,1)
+ #endif
          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
       &    SQUEEZE_RIGHT , myThid)
-         WRITE(msgBuf,'(A,7F8.4)')
+ c ToM>>>
-      &    ' SEAICE_GROWTH: Area increments 8, AREAITD = ',
+ #endif
-      &     AREAITD(1,1,:,bi,bj)
+ #ifdef SEAICE_ITD
+ C apply ice and snow thickness changes
+ C =================================================================
+          DO IT=1,nITD
+           DO J=1,sNy
+            DO I=1,sNx
+             HEFFITD(I,J,IT,bi,bj) = HEFFITD(I,J,IT,bi,bj)
+      &                            + d_HEFFbySublim_ITD(I,J,IT)
+      &                            + d_HEFFbyOCNonICE_ITD(I,J,IT)
+      &                            + d_HEFFbyATMonOCN_ITD(I,J,IT)
+      &                            + d_HEFFbyFLOODING_ITD(I,J,IT)
+             HSNOWITD(I,J,IT,bi,bj) = HSNOWITD(I,J,IT,bi,bj)
+      &                            + d_HSNWbySublim_ITD(I,J,IT)
+      &                            + d_HSNWbyATMonSNW_ITD(I,J,IT)
+      &                            + d_HSNWbyRAIN_ITD(I,J,IT)
+      &                            + d_HSNWbyOCNonSNW_ITD(I,J,IT)
+      &                            - d_HEFFbyFLOODING_ITD(I,J,IT)
+      &                            * ICE2SNOW
+            ENDDO
+           ENDDO
+          ENDDO
+ #endif
+ c ToM<<< debug seaice_growth
+         WRITE(msgBuf,msgBufForm)
+      &    ' SEAICE_GROWTH: Heff increments 9, HEFF = ',
+ #ifdef SEAICE_ITD
+      &     HEFFITD(1,1,:,bi,bj)
  #else
-      &    ' SEAICE_GROWTH: Heff increments 8, HEFF = ',
       &     HEFF(1,1,bi,bj)
+ #endif
          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
       &    SQUEEZE_RIGHT , myThid)
-         WRITE(msgBuf,'(A,7F8.4)')
+         WRITE(msgBuf,msgBufForm)
-      &    ' SEAICE_GROWTH: Area increments 8, AREA = ',
+      &    ' SEAICE_GROWTH: Area increments 9, AREA = ',
+ #ifdef SEAICE_ITD
+      &     AREAITD(1,1,:,bi,bj)
+ #else
       &     AREA(1,1,bi,bj)
  #endif
          CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
-Line 1965 
 CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bi
+Line 2031 
 CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bi
  CADJ STORE AREA(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
  #endif /* ALLOW_AUTODIFF_TAMC */
- c#ifdef SEAICE_ITD
+ #ifdef SEAICE_ITD
- cC       replaces Hibler '79 scheme and lead closing parameter
+ C--     account for lateral ice growth and melt only in thinnest category
- cC       because ITD accounts explicitly for lead openings and
+ C--     use HEFF, ARE, HSNOW, etc. temporarily for 1st category
- cC       different melt rates due to varying ice thickness
+ C       (this way we can use same code for ITD and non-ITD case)
- cC
- cC       only consider ice area loss due to total ice thickness loss;
- cC       ice area gain due to freezing of open water is handled above
- cC       under "gain of new ice over open water"
- cC
- cC       does not account for lateral melt of ice floes
- cC
- cC        AREAITD is incremented in section "gain of new ice over open water" above
- cC
- c        DO IT=1,nITD
- c         DO J=1,sNy
- c          DO I=1,sNx
- c          IF (HEFFITD(I,J,IT,bi,bj).LE.ZERO) THEN
- c           AREAITD(I,J,IT,bi,bj)=ZERO
- c          ENDIF
- c#ifdef ALLOW_SITRACER
- c           SItrAREA(I,J,bi,bj,3) = SItrAREA(I,J,bi,bj,3)
- c     &                           + AREAITD(I,J,IT,bi,bj)
- c#endif /* ALLOW_SITRACER */
- c         ENDDO
- c        ENDDO
- c       ENDDO
- c#else /* SEAICE_ITD */
          DO J=1,sNy
           DO I=1,sNx
- ctom<<< debugging
- #ifdef SEAICE_ITD
            HEFF(I,J,bi,bj)=HEFFITD(I,J,1,bi,bj)
            AREA(I,J,bi,bj)=AREAITD(I,J,1,bi,bj)
            HSNOW(I,J,bi,bj)=HSNOWITD(I,J,1,bi,bj)
            HEFFpreTH(I,J)=HEFFITDpreTH(I,J,1)
            AREApreTH(I,J)=AREAITDpreTH(I,J,1)
            recip_heffActual(I,J)=recip_heffActualMult(I,J,1)
+          ENDDO
+         ENDDO
+ C       all other categories only experience basal growth or melt,
+ C       i.e. change sin AREA only occur when all ice in a category is melted
+         IF (nITD .gt. 1) THEN
+          DO IT=2,nITD
+           DO J=1,sNy
+            DO I=1,sNx
+             IF (HEFFITD(I,J,IT,bi,bj).LE.ZERO) THEN
+              AREAITD(I,J,IT,bi,bj)=ZERO
+             ENDIF
+            ENDDO
+           ENDDO
+          ENDDO
+         ENDIF
  #endif
- ctom>>> debugging
+         DO J=1,sNy
+          DO I=1,sNx
            IF ( YC(I,J,bi,bj) .LT. ZERO ) THEN
             recip_HO=1. _d 0 / HO_south
-Line 2070 
 C apply tendency
+Line 2126 
 C apply tendency
            d_AREAbyOCN(I,J)=
       &        HALF*recip_HH*MIN( 0. _d 0,d_HEFFbyOCNonICE(I,J) )
  #endif /* ALLOW_DIAGNOSTICS */
- ctom<<< debugging
+          ENDDO
+         ENDDO
  #ifdef SEAICE_ITD
+ C       transfer 1st category values back into ITD variables
+         DO J=1,sNy
+          DO I=1,sNx
            HEFFITD(I,J,1,bi,bj)=HEFF(I,J,bi,bj)
            AREAITD(I,J,1,bi,bj)=AREA(I,J,bi,bj)
            HSNOWITD(I,J,1,bi,bj)=HSNOW(I,J,bi,bj)
- #endif
- ctom>>> debugging
           ENDDO
          ENDDO
- c#endif /* SEAICE_ITD */
+ #endif
  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
  Cgf 'bulk' linearization of area=f(HEFF)
-Line 2111 
 C finally update total AREA, HEFF, HSNOW
+Line 2169 
 C finally update total AREA, HEFF, HSNOW
          CALL SEAICE_ITD_SUM(bi, bj, myTime,myIter,myThid)
  #endif
- c ToM<<< debugging
-         DO J=1,sNy
-          DO I=1,sNx
-           if (I.eq.1 .and. J.eq.1) then
-            print *, 'd_HEFFbyNEG(I,J)      = ', d_HEFFbyNEG(I,J)
-            print *, 'd_HEFFbyOCNonICE(I,J) = ', d_HEFFbyOCNonICE(I,J)
-            print *, 'd_HEFFbyATMonOCN(I,J) = ', d_HEFFbyATMonOCN(I,J)
-            print *, 'd_HEFFbyATMonOCN_cover(I,J) = ',
-      &               d_HEFFbyATMonOCN_cover(I,J)
-            print *, 'd_HEFFbyATMonOCN_open(I,J) = ',
-      &               d_HEFFbyATMonOCN_open(I,J)
-            print *, 'd_HEFFbyFLOODING(I,J) = ', d_HEFFbyFLOODING(I,J)
-            print *, 'd_HEFFbySublim(I,J)   = ', d_HEFFbySublim(I,J)
-           endif
-          ENDDO
-         ENDDO
- c ToM>>>
  C ===================================================================
  C =============PART 5: determine ice salinity increments=============
  C ===================================================================
-Line 2149 
 CADJ &                       key = iicek
+Line 2190 
 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 2244 
 C set HSALT = 0 if HEFF = 0 and compute
+Line 2285 
 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 2327 
 C =================================
+Line 2367 
 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 2390 
 C in principle a_QSWbyATM_cover should a
+Line 2430 
 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
- cToM<<< debugging
-         print*,'------------------'
-         print*,'OcnModFrc: QNET = ',QNET(1,1,bi,bj)
-         print*,'OcnModFrc: QSW  = ',QSW(1,1,bi,bj)
-         print*,' '
-         print*,'r_QbyATM_cover  = ', r_QbyATM_cover(1,1)
-         print*,'r_QbyATM_open   = ', r_QbyATM_open(1,1)
-         print*,'a_QSWbyATM_cover = ', a_QSWbyATM_cover(1,1)
-         print*,'d_HEFFbyOCNonICE = ', d_HEFFbyOCNonICE(1,1)
-         print*,'d_HSNWbyOCNonSNW = ', d_HSNWbyOCNonSNW(1,1)
-         print*,'d_HEFFbyNEG = ', d_HEFFbyNEG(1,1)
-         print*,'d_HSNWbyNEG = ', d_HSNWbyNEG(1,1)
-         print*,'SNOW2ICE = ',SNOW2ICE
-         print*,'maskC       = ', maskC(1,1,kSurface,bi,bj)
-         print*,'------------------'
- cToM>>>
  C switch heat fluxes from 'effective' ice meters to W/m2
  C ======================================================
-Line 2441 
 CADJ STORE d_HSNWbyATMonSNW = comlev1_bi
+Line 2471 
 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 2499 
 C ======================================
+Line 2591 
 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 2513 
 C     If r_FWbySublim>0, then it is evap
+Line 2605 
 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 2585 
 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
+Line 2687 
 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 2635 
 C three that actually need intermediate
+Line 2773 
 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 2674 
 C the actual Freshwater flux of sublimat
+Line 2786 
 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 2691 
 C the latent heat flux
+Line 2803 
 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 2737 
 C close bi,bj loops
+Line 2812 
 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.8
 


changed lines


 
Added in v.1.13
 Legend:



Removed from v.1.8
 


changed lines


 
Added in v.1.13
-Removed from v.1.8
+Added in v.1.13

	ViewVC Help
Powered by ViewVC 1.1.22