/[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.11 by torge,
Mon Oct 22 22:18:09 2012 UTC
+revision 1.17 by torge,
Wed Apr 10 00:35:33 2013 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 99 
 C     number of surface interface layer
+Line 93 
 C     number of surface interface layer
        INTEGER kSurface
  C     IT :: ice thickness category index (MULTICATEGORIES and ITD code)
        INTEGER IT
-       _RL pFac
+ C     msgBuf      :: Informational/error message buffer
+ #ifdef ALLOW_BALANCE_FLUXES
+       CHARACTER*(MAX_LEN_MBUF) msgBuf
+ #elif (defined (SEAICE_DEBUG))
+       CHARACTER*(MAX_LEN_MBUF) msgBuf
+       CHARACTER*12 msgBufForm
+ #endif
  C     constants
+       _RL pFac
        _RL tempFrz, ICE2SNOW, SNOW2ICE
        _RL QI, QS, recip_QI
        _RL lhSublim
-Line 134 
 C     facilitate multi-category snow imp
+Line 135 
 C     facilitate multi-category snow imp
  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
-Line 165 
 C     AREA_PRE :: hold sea-ice fraction
+Line 170 
 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 leadIceThickMin
  #endif
  C     wind speed
-Line 190 
 C     temporary variables available for
+Line 194 
 C     temporary variables available for
  #ifdef SEAICE_ITD
        _RL r_QbyATMmult_cover  (1:sNx,1:sNy,MULTDIM)
        _RL r_FWbySublimMult    (1:sNx,1:sNy,MULTDIM)
+ c for lateral melt parameterization:
+       _RL latMeltFrac         (1:sNx,1:sNy,MULTDIM)
+       _RL latMeltRate         (1:sNx,1:sNy,MULTDIM)
+       _RL floeAlpha
+       _RL floeDiameter
+       _RL floeDiameterMin
+       _RL floeDiameterMax
  #endif
  C     a_QbyATM_cover :: available heat (in W/m^2) due to the interaction of
-Line 259 
 C     The change of mean ice thickness d
+Line 270 
 C     The change of mean ice thickness d
        _RL d_HEFFbyRLX         (1:sNx,1:sNy)
  #endif
+ #ifdef SEAICE_ITD
+       _RL d_HEFFbySublim_ITD         (1:sNx,1:sNy,1:nITD)
+       _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
  #ifdef ALLOW_DIAGNOSTICS
  C ICE/SNOW stocks tendencies associated with the various melt/freeze processes
        _RL d_AREAbyATM         (1:sNx,1:sNy)
-Line 281 
 C     Helper variables for diagnostics
+Line 305 
 C     Helper variables for diagnostics
        _RL HFsiGlob
        _RL FWF2HFsiTile(nSx,nSy)
        _RL FWF2HFsiGlob
-       CHARACTER*(max_len_mbuf) msgbuf
  #endif
  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
-Line 300 
 C     avoid unnecessary divisions in loo
+Line 323 
 C     avoid unnecessary divisions in loo
  c#ifdef SEAICE_ITD
  CToM this is now set by MULTDIM = nITD in SEAICE_SIZE.h
  C    (see SEAICE_SIZE.h and seaice_readparms.F)
  c     SEAICE_multDim = nITD
  c#endif
        recip_multDim        = SEAICE_multDim
        recip_multDim        = ONE / recip_multDim
-Line 335 
 C conversion factors to go from Q (W/m2)
+Line 358 
 C conversion factors to go from Q (W/m2)
  C conversion factors to go from precip (m/s) unit to HEFF (ice meters)
        convertPRECIP2HI=SEAICE_deltaTtherm*rhoConstFresh/SEAICE_rhoIce
        convertHI2PRECIP = ONE/convertPRECIP2HI
+ #ifdef SEAICE_ITD
+ c constants for lateral melt parameterization:
+ c following Steele (1992), Equ. 2
+       floeAlpha                  = 0.66 _d 0
+ c typical mean diameter used in CICE 4.1:
+ c      floeDiameter               = 300. _d 0
+ c parameters needed for variable floe diameter following Luepkes et al. (2012):
+       floeDiameterMin            = 8. _d 0
+       floeDiameterMax            = 300. _d 0
+ #endif
        DO bj=myByLo(myThid),myByHi(myThid)
         DO bi=myBxLo(myThid),myBxHi(myThid)
-Line 413 
 C =====================
+Line 446 
 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
+ c for lateral melt parameterization:
+             latMeltFrac(I,J,IT)                = 0.0 _d 0
+             latMeltRate(I,J,IT)                = 0.0 _d 0
  #endif
            ENDDO
           ENDDO
-Line 431 
 C ======================================
+Line 477 
 C ======================================
  C ===========PART 1: treat pathological cases (post advdiff)===========
  C =====================================================================
- #ifdef SEAICE_GROWTH_LEGACY
-         DO J=1,sNy
-          DO I=1,sNx
-           HEFFpreTH(I,J)=HEFFNM1(I,J,bi,bj)
-           HSNWpreTH(I,J)=HSNOW(I,J,bi,bj)
-           AREApreTH(I,J)=AREANM1(I,J,bi,bj)
-           d_HEFFbyNEG(I,J)=0. _d 0
-           d_HSNWbyNEG(I,J)=0. _d 0
- #ifdef ALLOW_DIAGNOSTICS
-           DIAGarrayA(I,J) = AREANM1(I,J,bi,bj)
-           DIAGarrayB(I,J) = AREANM1(I,J,bi,bj)
-           DIAGarrayC(I,J) = HEFFNM1(I,J,bi,bj)
-           DIAGarrayD(I,J) = HSNOW(I,J,bi,bj)
- #endif
-          ENDDO
-         ENDDO
- #ifdef SEAICE_ITD
-         DO IT=1,nITD
-          DO J=1,sNy
-           DO I=1,sNx
-            HEFFITDpreTH(I,J,IT)=HEFFITD(I,J,IT,bi,bj)
-            HSNWITDpreTH(I,J,IT)=HSNOWITD(I,J,IT,bi,bj)
-            AREAITDpreTH(I,J,IT)=AREAITD(I,J,IT,bi,bj)
-           ENDDO
-          ENDDO
-         ENDDO
- #endif
- #else /* SEAICE_GROWTH_LEGACY */
  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
  Cgf no dependency through pathological cases treatment
          IF ( SEAICEadjMODE.EQ.0 ) THEN
-Line 487 
 C           d_HEFFbyRLX(i,j) = 1. _d 1 *
+Line 502 
 C           d_HEFFbyRLX(i,j) = 1. _d 1 *
              d_HEFFbyRLX(i,j) = 1. _d 1 * siEps
             ENDIF
  #ifdef SEAICE_ITD
              AREAITD(I,J,1,bi,bj) = AREAITD(I,J,1,bi,bj)
       &                           +  d_AREAbyRLX(i,j)
              HEFFITD(I,J,1,bi,bj) = HEFFITD(I,J,1,bi,bj)
       &                           +  d_HEFFbyRLX(i,j)
-Line 506 
 CADJ STORE heff(:,:,bi,bj)  = comlev1_bi
+Line 521 
 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 531 
 C         by calling SEAICE_ITD_SUM
+Line 547 
 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
-           tmpscal2=0. _d 0
+         DO J=1,sNy
-           tmpscal3=0. _d 0
+          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 567 
 CToM        TICE will be updated at end
+Line 586 
 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 579 
 C 1.5) treat the case of area but no ice
+Line 598 
 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 600 
 C 2) treat the case of very small area:
+Line 623 
 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 620 
 C           (or install a warning in e.g
+Line 644 
 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 645 
 CADJ STORE area(:,:,bi,bj)  = comlev1_bi
+Line 672 
 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)
- C    lines of item 2.5 that were omitted:
+            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 677 
 C    which includes ridging as in item 2
+Line 709 
 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 */
  #ifdef SEAICE_DEBUG
- c ToM<<< debug seaice_growth
+ #ifdef SEAICE_ITD
-         WRITE(msgBuf,'(A,7F8.4)')
+         WRITE(msgBufForm,'(A,I2,A)') '(A,',nITD,'F14.10)'
-      &    ' 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(msgBufForm,'(A,A)') '(A,  F14.10)'
-         WRITE(msgBuf,'(A,7F8.4)')
+ #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 /* SEAICE_DEBUG */
- #endif
- #endif /* SEAICE_ITD */
  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
  C        end SEAICEadjMODE.EQ.0 statement:
          ENDIF
-Line 732 
 C 3) store regularized values of heff, h
+Line 763 
 C 3) store regularized values of heff, h
          DO IT=1,nITD
           DO J=1,sNy
            DO I=1,sNx
             HEFFITDpreTH(I,J,IT)=HEFFITD(I,J,IT,bi,bj)
             HSNWITDpreTH(I,J,IT)=HSNOWITD(I,J,IT,bi,bj)
             AREAITDpreTH(I,J,IT)=AREAITD(I,J,IT,bi,bj)
  C memorize areal and volume fraction of each ITD category
             IF (AREA(I,J,bi,bj) .GT. ZERO) THEN
              areaFracFactor(I,J,IT)=AREAITD(I,J,IT,bi,bj)/AREA(I,J,bi,bj)
             ELSE
- C           if there's no ice, potential growth starts in 1st category
+ C           if there is no ice, potential growth starts in 1st category
              IF (IT .EQ. 1) THEN
               areaFracFactor(I,J,IT)=ONE
              ELSE
               areaFracFactor(I,J,IT)=ZERO
              ENDIF
             ENDIF
            ENDDO
           ENDDO
          ENDDO
-Line 786 
 CADJ STORE heff(:,:,bi,bj)  = comlev1_bi
+Line 817 
 CADJ STORE heff(:,:,bi,bj)  = comlev1_bi
          ENDDO
  #endif /* SEAICE_VARIABLE_SALINITY */
- #endif /* SEAICE_GROWTH_LEGACY */
  #ifdef ALLOW_DIAGNOSTICS
          IF ( useDiagnostics ) THEN
           CALL DIAGNOSTICS_FILL(DIAGarrayA,'SIareaPR',0,1,3,bi,bj,myThid)
-Line 859 
 CADJ STORE HSNWpreTH = comlev1_bibj, key
+Line 888 
 CADJ STORE HSNWpreTH = comlev1_bibj, key
  #ifdef SEAICE_ITD
            IF (HEFFITDpreTH(I,J,IT) .GT. ZERO) THEN
- #ifdef SEAICE_GROWTH_LEGACY
-             tmpscal1          = MAX(SEAICE_area_reg/float(nITD),
-      &                              AREAITDpreTH(I,J,IT))
-             hsnowActualMult(I,J,IT) = HSNWITDpreTH(I,J,IT)/tmpscal1
-             tmpscal2               = HEFFITDpreTH(I,J,IT)/tmpscal1
-             heffActualMult(I,J,IT)  = MAX(tmpscal2,SEAICE_hice_reg)
- #else /* SEAICE_GROWTH_LEGACY */
  cif        regularize AREA with SEAICE_area_reg
             tmpscal1 = SQRT(AREAITDpreTH(I,J,IT) * AREAITDpreTH(I,J,IT)
       &                     + area_reg_sq)
  cif        heffActual calculated with the regularized AREA
             tmpscal2 = HEFFITDpreTH(I,J,IT) / tmpscal1
  cif        regularize heffActual with SEAICE_hice_reg (add lower bound)
             heffActualMult(I,J,IT) = SQRT(tmpscal2 * tmpscal2
       &                                  + hice_reg_sq)
  cif        hsnowActual calculated with the regularized AREA
             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,IT)  = AREAITDpreTH(I,J,IT) /
       &                 sqrt(HEFFITDpreTH(I,J,IT) * HEFFITDpreTH(I,J,IT)
       &                      + hice_reg_sq)
  cif       Do not regularize when HEFFpreTH = 0
            ELSE
-Line 889 
 cif       Do not regularize when HEFFpre
+Line 910 
 cif       Do not regularize when HEFFpre
            ENDIF
  #else /* SEAICE_ITD */
            IF (HEFFpreTH(I,J) .GT. ZERO) THEN
- #ifdef SEAICE_GROWTH_LEGACY
-             tmpscal1         = MAX(SEAICE_area_reg,AREApreTH(I,J))
-             hsnowActual(I,J) = HSNWpreTH(I,J)/tmpscal1
-             tmpscal2         = HEFFpreTH(I,J)/tmpscal1
-             heffActual(I,J)  = MAX(tmpscal2,SEAICE_hice_reg)
- #else /* SEAICE_GROWTH_LEGACY */
  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
-Line 903 
 Cif        regularize heffActual with SE
+Line 918 
 Cif        regularize heffActual with SE
             heffActual(I,J) = SQRT(tmpscal2 * tmpscal2 + hice_reg_sq)
  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
             recip_heffActual(I,J)  = AREApreTH(I,J) /
       &                 sqrt(HEFFpreTH(I,J)*HEFFpreTH(I,J) + hice_reg_sq)
-Line 1031 
 CADJ &                       key = iicek
+Line 1045 
 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 J=1,sNy
+           DO I=1,sNx
+ 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
- #ifndef SEAICE_ITD
           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
- #ifndef SEAICE_ITD
- 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
- #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 1058 
 C    (instead of heffActual and latentHe
+Line 1077 
 C    (instead of heffActual and latentHe
            ENDDO
           ENDDO
          ENDDO
+ #endif
  #ifdef ALLOW_AUTODIFF_TAMC
  CADJ STORE heffActualMult = comlev1_bibj, key = iicekey, byte = isbyte
-Line 1110 
 CADJ &     comlev1_bibj, key = iicekey,
+Line 1130 
 CADJ &     comlev1_bibj, key = iicekey,
  C     update TICE & TICES
  #ifdef SEAICE_ITD
  C     calculate area weighted mean
- C     (although the ice's temperature relates to its energy content
+ C     (although the ice temperature relates to its energy content
  C      and hence should be averaged weighted by ice volume,
  C      the temperature here is a result of the fluxes through the ice surface
  C      computed individually for each single category in SEAICE_SOLVE4TEMP
  C      and hence is averaged area weighted [areaFracFactor])
             TICE(I,J,bi,bj) = TICE(I,J,bi,bj)
       &          +  ticeOutMult(I,J,IT)*areaFracFactor(I,J,IT)
-Line 1173 
 CADJ STORE a_FWbySublim    = comlev1_bib
+Line 1193 
 CADJ STORE a_FWbySublim    = comlev1_bib
  C switch heat fluxes from W/m2 to 'effective' ice meters
  #ifdef SEAICE_ITD
          DO IT=1,nITD
           DO J=1,sNy
            DO I=1,sNx
             a_QbyATMmult_cover(I,J,IT)   = a_QbyATMmult_cover(I,J,IT)
-Line 1190 
 C     Negative sublimation is resublimat
+Line 1210 
 C     Negative sublimation is resublimat
             a_FWbySublimMult(I,J,IT) = SEAICE_deltaTtherm*recip_rhoIce
       &            * a_FWbySublimMult(I,J,IT)*AREAITDpreTH(I,J,IT)
             r_FWbySublimMult(I,J,IT)=a_FWbySublimMult(I,J,IT)
            ENDDO
           ENDDO
          ENDDO
          DO J=1,sNy
-Line 1246 
 CADJ STORE r_FWbySublim    = comlev1_bib
+Line 1266 
 CADJ STORE r_FWbySublim    = comlev1_bib
  Cgf no additional dependency through ice cover
          IF ( SEAICEadjMODE.GE.3 ) THEN
  #ifdef SEAICE_ITD
           DO IT=1,nITD
            DO J=1,sNy
             DO I=1,sNx
              a_QbyATMmult_cover(I,J,IT)   = 0. _d 0
-Line 1254 
 Cgf no additional dependency through ice
+Line 1274 
 Cgf no additional dependency through ice
              a_QSWbyATMmult_cover(I,J,IT) = 0. _d 0
             ENDDO
            ENDDO
           ENDDO
  #else
           DO J=1,sNy
            DO I=1,sNx
-Line 1311 
 C available turbulent flux
+Line 1331 
 C available turbulent flux
           ENDDO
          ENDDO
+ #ifdef SEAICE_ITD
+ C determine lateral melt rate at floe edges based on an
+ C average floe diameter or a floe size distribution
+ C following Steele (1992, Tab. 2)
+ C ======================================================
+         DO J=1,sNy
+          DO I=1,sNx
+           tmpscal1=(theta(I,J,kSurface,bi,bj)-tempFrz)
+           tmpscal2=sqrt(0.87 + 0.067*UG(i,j)) * UG(i,j)
+ c
+ c variable floe diameter following Luepkes et al. (2012, JGR, Equ. 26) with beta=1
+           tmpscal3=ONE/(ONE-(floeDiameterMin/floeDiameterMax))
+           floeDiameter = floeDiameterMin
+      &                 * (tmpscal3 / (tmpscal3-AREApreTH(I,J)))
+           DO IT=1,nITD
+ C following the idea of SEAICE_areaLossFormula == 1:
+            IF (r_QbyATMmult_cover(i,j,it).LT.ZERO .OR.
+      &         r_QbyATM_open(i,j) .LT.ZERO .OR.
+      &         r_QbyOCN(i,j)      .LT.ZERO) THEN
+ c lateral melt rate as suggested by Perovich, 1983 (PhD thesis)
+ c            latMeltRate(i,j,it) = 1.6 _d -6 * tmpscal1**1.36
+ c lateral melt rate as suggested by Maykut and Perovich, 1987 (JGR 92(C7)), Equ. 24
+ c            latMeltRate(i,j,it) = 13.5 _d -6 * tmpscal2 * tmpscal1**1.3
+ c further suggestion by Maykut and Perovich to avoid latMeltRate -> 0 for UG -> 0
+             latMeltRate(i,j,it) = (1.6 _d -6 + 13.5 _d -6 * tmpscal2)
+      &                          * tmpscal1**1.3
+ c factor determining fraction of area and ice volume reduction
+ c due to lateral melt
+             latMeltFrac(i,j,it) =
+      &       latMeltRate(i,j,it)*SEAICE_deltaTtherm*PI /
+      &       (floeAlpha * floeDiameter)
+             latMeltFrac(i,j,it)=max(ZERO, min(latMeltFrac(i,j,it),ONE))
+             print*,'latMelt',it,tmpscal1,latMeltRate(i,j,it),
+      &                       areaitd(i,j,it,bi,bj),latMeltFrac(i,j,it)
+            ELSE
+ c            latMeltRate(i,j,it)=0.0 _d 0
+             latMeltFrac(i,j,it)=0.0 _d 0
+             print*,'latMelt',it,' 0.0 0.0 ',
+      &                       areaitd(i,j,it,bi,bj),latMeltFrac(i,j,it)
+            ENDIF
+           ENDDO
+          ENDDO
+         ENDDO
+ #endif
  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
          CALL ZERO_ADJ_1D( sNx*sNy, r_QbyOCN, myThid)
  #endif
-Line 1336 
 C     First sublimate/deposite snow
+Line 1401 
 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 1360 
 C     If anything is left, sublimate ice
+Line 1424 
 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 1390 
 C     remove the fusion part for the res
+Line 1454 
 C     remove the fusion part for the res
          ENDDO
  #ifdef SEAICE_ITD
  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
+ #endif /* SEAICE_DEBUG */
  C compute ice thickness tendency due to ice-ocean interaction
  C ===========================================================
-Line 1424 
 CADJ STORE r_QbyOCN = comlev1_bibj,key=i
+Line 1491 
 CADJ STORE r_QbyOCN = comlev1_bibj,key=i
          DO IT=1,nITD
           DO J=1,sNy
            DO I=1,sNx
  C          ice growth/melt due to ocean heat r_QbyOCN (W/m^2) is
  C          equally distributed under the ice and hence weighted by
  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 1457 
 C          fractional area of each thick
+Line 1529 
 C          fractional area of each thick
  #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
+ #endif /* SEAICE_DEBUG */
  C compute snow melt tendency due to snow-atmosphere interaction
  C ==================================================================
-Line 1484 
 CADJ STORE r_QbyATM_cover = comlev1_bibj
+Line 1550 
 CADJ STORE r_QbyATM_cover = comlev1_bibj
  #endif /* ALLOW_AUTODIFF_TAMC */
  #ifdef SEAICE_ITD
          DO IT=1,nITD
           DO J=1,sNy
            DO I=1,sNx
  C     Convert to standard units (meters of ice) rather than to meters
-Line 1496 
 C     of snow. This appears to be more r
+Line 1562 
 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)
+            r_QbyATMmult_cover(I,J,IT)=r_QbyATMmult_cover(I,J,IT)
-      &                           + tmpscal2*ICE2SNOW
-            r_QbyATMmult_cover(I,J,IT)=r_QbyATMmult_cover(I,J,IT)
       &                           - tmpscal2
            ENDDO
           ENDDO
          ENDDO
  #else /* SEAICE_ITD */
          DO J=1,sNy
           DO I=1,sNx
-Line 1524 
 Cgf no additional dependency through sno
+Line 1589 
 Cgf no additional dependency through sno
  #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
+ #endif /* SEAICE_DEBUG */
  C compute ice thickness tendency due to the atmosphere
  C ====================================================
-Line 1556 
 Cgf the v1.81=>v1.82 revision would chan
+Line 1615 
 Cgf the v1.81=>v1.82 revision would chan
  Cgf warming conditions, the lab_sea results were not changed.
  #ifdef SEAICE_ITD
          DO IT=1,nITD
           DO J=1,sNy
            DO I=1,sNx
- #ifdef SEAICE_GROWTH_LEGACY
+            tmpscal1 = HEFFITDpreTH(I,J,IT)
-            tmpscal2 = MAX(-HEFFITD(I,J,IT,bi,bj),
+      &              + d_HEFFbySublim_ITD(I,J,IT)
-      &                     r_QbyATMmult_cover(I,J,IT))
+      &              + d_HEFFbyOCNonICE_ITD(I,J,IT)
- #else
+            tmpscal2 = MAX(-tmpscal1,
-            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,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
+           ENDDO
+          ENDDO
+         ENDDO
  #ifdef ALLOW_SITRACER
-            SItrHEFF(I,J,bi,bj,3) = SItrHEFF(I,J,bi,bj,3)
+         DO J=1,sNy
-      &                           + HEFFITD(I,J,IT,bi,bj)
+          DO I=1,sNx
+           SItrHEFF(I,J,bi,bj,3) = SItrHEFF(I,J,bi,bj,2)
+      &                          + d_HEFFbyATMonOCN_cover(I,J)
+          ENDDO
+         ENDDO
  #endif
-           ENDDO
-          ENDDO
-         ENDDO
  #else /* SEAICE_ITD */
          DO J=1,sNy
           DO I=1,sNx
- #ifdef SEAICE_GROWTH_LEGACY
-           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
       &        AREApreTH(I,J) * r_QbyOCN(I,J))
- #endif /* SEAICE_GROWTH_LEGACY */
            d_HEFFbyATMonOCN_cover(I,J)=tmpscal2
            d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal2
-Line 1608 
 C         Limit ice growth by potential
+Line 1666 
 C         Limit ice growth by potential
  #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
+ #endif /* SEAICE_DEBUG */
  C add snow precipitation to HSNOW.
  C =================================================
-Line 1668 
 C           add precip to the fresh wate
+Line 1729 
 C           add precip to the fresh wate
           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
  #else
          DO J=1,sNy
           DO I=1,sNx
-Line 1691 
 C end of IF statement snowPrecipFile:
+Line 1752 
 C end of IF statement snowPrecipFile:
  #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
+ #endif /* SEAICE_DEBUG */
  C compute snow melt due to heat available from ocean.
  C =================================================================
-Line 1724 
 CADJ STORE r_QbyOCN = comlev1_bibj,key=i
+Line 1779 
 CADJ STORE r_QbyOCN = comlev1_bibj,key=i
          DO IT=1,nITD
           DO J=1,sNy
            DO I=1,sNx
-            tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW*areaFracFactor(I,J,IT),
+            tmpscal4 = HSNWITDpreTH(I,J,IT)
-      &                  -HSNOWITD(I,J,IT,bi,bj))
+      &              + d_HSNWbySublim_ITD(I,J,IT)
+      &              + d_HSNWbyATMonSNW_ITD(I,J,IT)
+      &              + d_HSNWbyRAIN_ITD(I,J,IT)
+            tmpscal1=MAX(r_QbyOCN(i,j)*ICE2SNOW*areaFracFactor(I,J,IT),
+      &                  -tmpscal4)
             tmpscal2=MIN(tmpscal1,0. _d 0)
  #ifdef SEAICE_MODIFY_GROWTH_ADJ
  Cgf no additional dependency through snow
             if ( SEAICEadjMODE.GE.2 ) tmpscal2 = 0. _d 0
  #endif
+            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 1757 
 Cgf no additional dependency through sno
+Line 1816 
 Cgf no additional dependency through sno
  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
+ #endif /* SEAICE_DEBUG */
  C gain of new ice over open water
  C ===============================
-Line 1787 
 CADJ STORE a_QSWbyATM_open = comlev1_bib
+Line 1840 
 CADJ STORE a_QSWbyATM_open = comlev1_bib
          DO J=1,sNy
           DO I=1,sNx
+ #ifdef SEAICE_ITD
+ 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) *
-Line 1797 
 C           that is therefore not availa
+Line 1861 
 C           that is therefore not availa
  C           impose -HEFF as the maxmum melting if SEAICE_doOpenWaterMelt
  C           or 0. otherwise (no melting if not SEAICE_doOpenWaterMelt)
              tmpscal3=facOpenGrow*MAX(tmpscal1-tmpscal2,
-      &           -HEFF(I,J,bi,bj)*facOpenMelt)*HEFFM(I,J,bi,bj)
+      &           -tmpscal4*facOpenMelt)*HEFFM(I,J,bi,bj)
+ #ifdef SEAICE_ITD
+ C         ice growth in open water adds to first category
+           d_HEFFbyATMonOCN_open_ITD(I,J,1)=tmpscal3
+           d_HEFFbyATMonOCN_ITD(I,J,1)     =d_HEFFbyATMonOCN_ITD(I,J,1)
+      &                                    +tmpscal3
+ #endif
            d_HEFFbyATMonOCN_open(I,J)=tmpscal3
            d_HEFFbyATMonOCN(I,J)=d_HEFFbyATMonOCN(I,J)+tmpscal3
            r_QbyATM_open(I,J)=r_QbyATM_open(I,J)-tmpscal3
- #ifdef SEAICE_ITD
- C         open water area fraction
-           tmpscal0 = ONE-AREApreTH(I,J)
- C         determine thickness of new ice
- ctomC         considering the entire open water area to refreeze
- 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 IT=1,nITD
-            IF (tmpscal1.LT.Hlimit(IT)) THEN
-             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,IT,bi,bj) = AREAITD(I,J,IT,bi,bj)
-      &                           + tmpscal0
-             tmpscal0=ZERO
-            ENDIF
-           ENDDO
- #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
+ #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
+ #endif /* SEAICE_DEBUG */
  C convert snow to ice if submerged.
  C =================================
- #ifndef SEAICE_GROWTH_LEGACY
  C note: in legacy, this process is done at the end
  #ifdef ALLOW_AUTODIFF_TAMC
  CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
-Line 1890 
 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
+Line 1921 
 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
  #endif /* ALLOW_AUTODIFF_TAMC */
          IF ( SEAICEuseFlooding ) THEN
  #ifdef SEAICE_ITD
           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
+            ENDDO
-             HSNOWITD(I,J,IT,bi,bj)= HSNOWITD(I,J,IT,bi,bj) - tmpscal1
+           ENDDO
-      &                            * ICE2SNOW
+          ENDDO
-            ENDDO
-           ENDDO
-          ENDDO
  #else
           DO J=1,sNy
            DO I=1,sNx
-Line 1918 
 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
+Line 1955 
 CADJ STORE hsnow(:,:,bi,bj) = comlev1_bi
           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 /* SEAICE_DEBUG */
-      &     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
+ #ifdef SEAICE_DEBUG
+ 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,
       &    SQUEEZE_RIGHT , myThid)
  c ToM>>>
- #endif
+ #endif /* SEAICE_DEBUG */
  C ===================================================================
  C ==========PART 4: determine ice cover fraction increments=========-
-Line 1969 
 CADJ STORE AREA(:,:,bi,bj) = comlev1_bib
+Line 2039 
 CADJ STORE AREA(:,:,bi,bj) = comlev1_bib
  #endif /* ALLOW_AUTODIFF_TAMC */
  #ifdef SEAICE_ITD
- C       replaces Hibler '79 scheme and lead closing parameter
+ C--     account for lateral ice growth and melt only in thinnest category
- C       because ITD accounts explicitly for lead openings and
+ C--     use HEFF, ARE, HSNOW, etc. temporarily for 1st category
- C       different melt rates due to varying ice thickness
+ C       (this way we can use same code for ITD and non-ITD case)
- C
+         DO J=1,sNy
- C       only consider ice area loss due to total ice thickness loss;
+          DO I=1,sNx
- C       ice area gain due to freezing of open water is handled above
+           HEFF(I,J,bi,bj)=HEFFITD(I,J,1,bi,bj)
- C       under "gain of new ice over open water"
+           AREA(I,J,bi,bj)=AREAITD(I,J,1,bi,bj)
- C
+           HSNOW(I,J,bi,bj)=HSNOWITD(I,J,1,bi,bj)
- C       does not account for lateral melt of ice floes
+           HEFFpreTH(I,J)=HEFFITDpreTH(I,J,1)
- C
+           AREApreTH(I,J)=AREAITDpreTH(I,J,1)
- C        AREAITD is incremented in section "gain of new ice over open water" above
+           recip_heffActual(I,J)=recip_heffActualMult(I,J,1)
- C
+          ENDDO
-         DO IT=1,nITD
+         ENDDO
-          DO J=1,sNy
+ C       all other categories only experience basal growth or melt,
-          DO I=1,sNx
+ C       i.e. change sin AREA only occur when all ice in a category is melted
-            IF (HEFFITD(I,J,IT,bi,bj).LE.ZERO) THEN
+         IF (nITD .gt. 1) THEN
-             AREAITD(I,J,IT,bi,bj)=ZERO
+          DO IT=2,nITD
-            ENDIF
+           DO J=1,sNy
- #ifdef ALLOW_SITRACER
+            DO I=1,sNx
-            SItrAREA(I,J,bi,bj,3) = SItrAREA(I,J,bi,bj,3)
+             IF (HEFFITD(I,J,IT,bi,bj).LE.ZERO) THEN
-      &                           + AREAITD(I,J,IT,bi,bj)
+              AREAITD(I,J,IT,bi,bj)=ZERO
- #endif /* ALLOW_SITRACER */
+             ELSE
-           ENDDO
+ c            melt ice laterally based on an average floe sice
-          ENDDO
+ c            following Steele (1992)
-         ENDDO
+              AREAITD(I,J,IT,bi,bj) = AREAITD(I,J,IT,bi,bj)
- #else /* SEAICE_ITD */
+      &                             * (ONE - latMeltFrac(I,J,IT))
+              AREAITD(I,J,IT,bi,bj) = max(ZERO,AREAITD(I,J,IT,bi,bj))
+             ENDIF
+            ENDDO
+           ENDDO
+          ENDDO
+         ENDIF
+ #endif
          DO J=1,sNy
           DO I=1,sNx
-Line 2003 
 C
+Line 2080 
 C
            ELSE
             recip_HO=1. _d 0 / HO
            ENDIF
- #ifdef SEAICE_GROWTH_LEGACY
+           recip_HH = recip_heffActual(I,J)
-            tmpscal0=HEFF(I,J,bi,bj) - d_HEFFbyATMonOCN(I,J)
-            recip_HH = AREApreTH(I,J) /(tmpscal0+.00001 _d 0)
- #else
-            recip_HH = recip_heffActual(I,J)
- #endif
  C gain of ice over open water : computed from
  C   (SEAICE_areaGainFormula.EQ.1) from growth by ATM
-Line 2064 
 C apply tendency
+Line 2136 
 C apply tendency
  #endif /* ALLOW_DIAGNOSTICS */
           ENDDO
          ENDDO
- #endif /* SEAICE_ITD */
+ #ifdef SEAICE_ITD
+ C       transfer 1st category values back into ITD variables
+         DO J=1,sNy
+          DO I=1,sNx
+           HEFFITD(I,J,1,bi,bj)=HEFF(I,J,bi,bj)
+           AREAITD(I,J,1,bi,bj)=AREA(I,J,bi,bj)
+           HSNOWITD(I,J,1,bi,bj)=HSNOW(I,J,bi,bj)
+          ENDDO
+         ENDDO
+ #endif
  #if (defined ALLOW_AUTODIFF_TAMC && defined SEAICE_MODIFY_GROWTH_ADJ)
  Cgf 'bulk' linearization of area=f(HEFF)
          IF ( SEAICEadjMODE.GE.1 ) THEN
  #ifdef SEAICE_ITD
           DO IT=1,nITD
            DO J=1,sNy
             DO I=1,sNx
              AREAITD(I,J,IT,bi,bj) = AREAITDpreTH(I,J,IT) + 0.1 _d 0 *
-Line 2135 
 CADJ &                       key = iicek
+Line 2216 
 CADJ &                       key = iicek
  #ifdef SEAICE_VARIABLE_SALINITY
- #ifdef SEAICE_GROWTH_LEGACY
- # ifdef ALLOW_AUTODIFF_TAMC
- CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
- # endif /* ALLOW_AUTODIFF_TAMC */
-         DO J=1,sNy
-          DO I=1,sNx
- C set HSALT = 0 if HSALT < 0 and compute salt to remove from ocean
-           IF ( HSALT(I,J,bi,bj) .LT. 0.0 ) THEN
-              saltFluxAdjust(I,J) = - HEFFM(I,J,bi,bj) *
-      &            HSALT(I,J,bi,bj) * recip_deltaTtherm
-              HSALT(I,J,bi,bj) = 0.0 _d 0
-           ENDIF
-          ENDDO
-         ENDDO
- #endif /* SEAICE_GROWTH_LEGACY */
  #ifdef ALLOW_AUTODIFF_TAMC
  CADJ STORE hsalt(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
  #endif /* ALLOW_AUTODIFF_TAMC */
-Line 2197 
 C update HSALT based on surface saltFlux
+Line 2262 
 C update HSALT based on surface saltFlux
       &         saltFlux(I,J,bi,bj) * SEAICE_deltaTtherm
            saltFlux(I,J,bi,bj) =
       &         saltFlux(I,J,bi,bj) + saltFluxAdjust(I,J)
- #ifdef SEAICE_GROWTH_LEGACY
- C set HSALT = 0 if HEFF = 0 and compute salt to dump into ocean
-           IF ( HEFF(I,J,bi,bj) .EQ. 0.0 ) THEN
-            saltFlux(I,J,bi,bj) = saltFlux(I,J,bi,bj) -
-      &          HEFFM(I,J,bi,bj) * HSALT(I,J,bi,bj) * recip_deltaTtherm
-            HSALT(I,J,bi,bj) = 0.0 _d 0
- #ifdef ALLOW_SALT_PLUME
-            saltPlumeFlux(i,j,bi,bj) = 0.0 _d 0
- #endif /* ALLOW_SALT_PLUME */
-           ENDIF
- #endif /* SEAICE_GROWTH_LEGACY */
           ENDDO
          ENDDO
  #endif /* SEAICE_VARIABLE_SALINITY */
- C =======================================================================
- C ==LEGACY PART 6 (LEGACY) treat pathological cases, then do flooding ===
- C =======================================================================
- #ifdef SEAICE_GROWTH_LEGACY
- C treat values of ice cover fraction oustide
- C the [0 1] range, and other such issues.
- C ===========================================
- Cgf note: this part cannot be heat and water conserving
- #ifdef ALLOW_AUTODIFF_TAMC
- CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
- CADJ &                       key = iicekey, byte = isbyte
- CADJ STORE heff(:,:,bi,bj) = comlev1_bibj,
- CADJ &                       key = iicekey, byte = isbyte
- #endif /* ALLOW_AUTODIFF_TAMC */
-         DO J=1,sNy
-          DO I=1,sNx
- C     NOW SET AREA(I,J,bi,bj)=0 WHERE THERE IS NO ICE
- CML   replaced "/.0001 _d 0" by "*1. _d 4", 1e-4 is probably
- CML   meant to be something like a minimum thickness
-           AREA(I,J,bi,bj)=MIN(AREA(I,J,bi,bj),HEFF(I,J,bi,bj)*1. _d 4)
-          ENDDO
-         ENDDO
- #ifdef ALLOW_AUTODIFF_TAMC
- CADJ STORE area(:,:,bi,bj) = comlev1_bibj,
- CADJ &                       key = iicekey, byte = isbyte
- #endif /* ALLOW_AUTODIFF_TAMC */
-         DO J=1,sNy
-          DO I=1,sNx
- C NOW TRUNCATE AREA
-           AREA(I,J,bi,bj)=MIN(ONE,AREA(I,J,bi,bj))
-          ENDDO
-         ENDDO
- #ifdef ALLOW_AUTODIFF_TAMC
- CADJ STORE area(:,:,bi,bj)  = comlev1_bibj,
- CADJ &                        key = iicekey, byte = isbyte
- CADJ STORE hsnow(:,:,bi,bj) = comlev1_bibj,
- CADJ &                         key = iicekey, byte = isbyte
- #endif /* ALLOW_AUTODIFF_TAMC */
-         DO J=1,sNy
-          DO I=1,sNx
-           AREA(I,J,bi,bj) = MAX(ZERO,AREA(I,J,bi,bj))
-           HSNOW(I,J,bi,bj)  = MAX(ZERO,HSNOW(I,J,bi,bj))
-           AREA(I,J,bi,bj) = AREA(I,J,bi,bj)*HEFFM(I,J,bi,bj)
-           HEFF(I,J,bi,bj) = HEFF(I,J,bi,bj)*HEFFM(I,J,bi,bj)
- #ifdef SEAICE_CAP_HEFF
- C     This is not energy conserving, but at least it conserves fresh water
-           tmpscal0         = -MAX(HEFF(I,J,bi,bj)-MAX_HEFF,0. _d 0)
-           d_HEFFbyNeg(I,J) = d_HEFFbyNeg(I,J) + tmpscal0
-           HEFF(I,J,bi,bj)  = HEFF(I,J,bi,bj)  + tmpscal0
- #endif /* SEAICE_CAP_HEFF */
-           HSNOW(I,J,bi,bj)  = HSNOW(I,J,bi,bj)*HEFFM(I,J,bi,bj)
-          ENDDO
-         ENDDO
- C convert snow to ice if submerged.
- C =================================
-         IF ( SEAICEuseFlooding ) THEN
-          DO J=1,sNy
-           DO I=1,sNx
-            tmpscal0     = (HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
-      &              +HEFF(I,J,bi,bj)*SEAICE_rhoIce)*recip_rhoConst
-            tmpscal1 = MAX( 0. _d 0, tmpscal0 - HEFF(I,J,bi,bj))
-            d_HEFFbyFLOODING(I,J)=tmpscal1
-            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 /* SEAICE_GROWTH_LEGACY */
  #ifdef ALLOW_SITRACER
          DO J=1,sNy
           DO I=1,sNx
-Line 2322 
 C compute total of "mult" fluxes for oce
+Line 2297 
 C compute total of "mult" fluxes for oce
           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
-Line 2352 
 CADJ STORE d_hsnwbyocnonsnw = comlev1_bi
+Line 2327 
 CADJ STORE d_hsnwbyocnonsnw = comlev1_bi
          DO J=1,sNy
           DO I=1,sNx
            QNET(I,J,bi,bj) = r_QbyATM_cover(I,J) + r_QbyATM_open(I,J)
- #ifndef SEAICE_GROWTH_LEGACY
- C in principle a_QSWbyATM_cover should always be included here, however
- 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_HSNWbyOCNonSNW(I,J)*SNOW2ICE
       &           + d_HEFFbyNEG(I,J)
-Line 2400 
 CADJ &                       key = iicek
+Line 2371 
 CADJ &                       key = iicek
  # endif /* ALLOW_AUTODIFF_TAMC */
  cgf 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 arbitrarily to be e.g. the ocean SST. Indeed the present seaice model
  C implies a constant ice temperature of 0degC. If melt/freeze water is exchanged
  C at a different temperature, it leads to a loss of conservation in the
  C ocean+ice system. While this is mostly a serious issue in the
  C real fresh water + non linear free surface framework, a mismatch
  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
-Line 2460 
 CML according to wave-length) fluxes but
+Line 2431 
 CML according to wave-length) fluxes but
  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)
-Line 2501 
 C linFS, rain/evap get a special treatme
+Line 2470 
 C linFS, rain/evap get a special treatme
             tmpscal2= ZERO
        ENDIF
             SItflux(I,J,bi,bj)=
       &            SIatmQnt(I,J,bi,bj)-tmpscal1-tmpscal2
            ENDDO
           ENDDO
-Line 2543 
 c and the flux leaving/entering the ocea
+Line 2512 
 c and the flux leaving/entering the ocea
       &     + a_FWbySublim(I,J) * SEAICE_rhoIce * recip_deltaTtherm
           ENDDO
          ENDDO
  #ifdef ALLOW_SSH_GLOBMEAN_COST_CONTRIBUTION
  C--
-Line 2621 
 C ======================================
+Line 2590 
 C ======================================
         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
-Line 2642 
 c to translate global mean FWF adjusteme
+Line 2611 
 c to translate global mean FWF adjusteme
         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
-Line 2754 
 CADJ STORE FWF2HFsiTile = comlev1, key=i
+Line 2723 
 CADJ STORE FWF2HFsiTile = comlev1, key=i
  # 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)
-Line 2784 
 c 3) balancing adjustments
+Line 2753 
 c 3) balancing adjustments
           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=
-Line 2807 
 c           no qnet or tflux adjustement
+Line 2776 
 c           no qnet or tflux adjustement
        ENDDO
        IF ( balancePrintMean ) THEN
         _BEGIN_MASTER( myThid )
-        WRITE(msgbuf,'(a,a,e24.17)') 'rm Global mean of ',
+        WRITE(msgBuf,'(a,a,e24.17)') 'rm Global mean of ',
       &      'SIatmFW = ', tmpscal0
         CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
       &      SQUEEZE_RIGHT , myThid)
-Line 2828 
 c           no qnet or tflux adjustement
+Line 2797 
 c           no qnet or tflux adjustement
        ENDDO
        IF ( balancePrintMean ) THEN
         _BEGIN_MASTER( myThid )
-        WRITE(msgbuf,'(a,a,e24.17)') 'rm Global mean of ',
+        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 /* */
+ #endif /* ALLOW_BALANCE_FLUXES */
  #ifdef ALLOW_DIAGNOSTICS
  c these diags need to be done outside of the bi,bj loop so that

 Legend:



Removed from v.1.11
 


changed lines


 
Added in v.1.17
 Legend:



Removed from v.1.11
 


changed lines


 
Added in v.1.17
-Removed from v.1.11
+Added in v.1.17

	ViewVC Help
Powered by ViewVC 1.1.22