#include "SEAICE_OPTIONS.h"
#ifdef ALLOW_EXF
# include "EXF_OPTIONS.h"
#endif
#ifdef ALLOW_SALT_PLUME
# include "SALT_PLUME_OPTIONS.h"
#endif
#ifdef ALLOW_AUTODIFF
# include "AUTODIFF_OPTIONS.h"
#endif

CBOP
C     !ROUTINE: SEAICE_GROWTH_ADX
C     !INTERFACE:
      SUBROUTINE SEAICE_GROWTH_ADX
     I  ( myTime, myIter, myThid )
C     !DESCRIPTION: \bv
C     *==========================================================*
C     | SUBROUTINE seaice_growth_adx
C     | o Adjointable seaice thermodynamic code based on.
C     |   Ian Fenty's code with modifications for nonlinear free
C     |   surface budget closure by A.Bigdeli.
C     *==========================================================*
C     \ev

C     !USES:
      IMPLICIT NONE
C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "DYNVARS.h"
#include "GRID.h"
#include "FFIELDS.h"
#include "SEAICE_SIZE.h"
#include "SEAICE_PARAMS.h"
#include "SEAICE.h"
#include "SEAICE_TRACER.h"
#ifdef ALLOW_EXF
# include "EXF_PARAM.h"
# include "EXF_FIELDS.h"
#endif
#ifdef ALLOW_SALT_PLUME
# include "SALT_PLUME.h"
#endif
#ifdef ALLOW_AUTODIFF_TAMC
# include "tamc.h"
#endif

C     !INPUT/OUTPUT PARAMETERS:
C     === Routine arguments ===
C     myTime :: Simulation time
C     myIter :: Simulation timestep number
C     myThid :: Thread no. that called this routine.
      _RL myTime
      INTEGER myIter, myThid
CEOP

#if (defined ALLOW_EXF) && (defined ALLOW_ATM_TEMP)
C     !FUNCTIONS:
#ifdef ALLOW_DIAGNOSTICS
      LOGICAL  DIAGNOSTICS_IS_ON
      EXTERNAL DIAGNOSTICS_IS_ON
#endif

C     !LOCAL VARIABLES:
C     === Local variables ===
C
C unit/sign convention:

C HEFF is effective Hice thickness (m3/m2)
C HSNOW is Heffective snow thickness (m3/m2)
C HSALT is Heffective salt content (g/m2)
C AREA is the seaice cover fraction (0<=AREA<=1)

C     i,j,bi,bj :: Loop counters
      INTEGER i, j, bi, bj
C     number of surface interface layer
      INTEGER kSurface
C     IT :: ice thickness category index (MULTICATEGORIES and ITD code)
      INTEGER IT
C     msgBuf      :: Informational/error message buffer
c     CHARACTER*(MAX_LEN_MBUF) msgBuf
C     constants
      _RL tempFrz, ICE2SNOW, SNOW2ICE, surf_theta
      _RL QI, QS
      _RL RHOSW
      _RL lhSublim

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     facilitate multi-category snow implementation
      _RL pFac, pFacSnow

C     temporary variables available for the various computations
      _RL tmpscal0, tmpscal1, tmpscal2, tmpscal3

#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 hiceActual          (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_hiceActual    (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)

C     wind speed
      _RL UG                  (1:sNx,1:sNy)

C     temporary variables available for the various computations
      _RL ticeInMult          (1:sNx,1:sNy,nITD)
      _RL ticeOutMult         (1:sNx,1:sNy,nITD)
      _RL hiceActualMult      (1:sNx,1:sNy,nITD)
      _RL hsnowActualMult     (1:sNx,1:sNy,nITD)

      _RL F_io_net_mult     (1:sNx,1:sNy,nITD)
      _RL F_ia_net_mult     (1:sNx,1:sNy,nITD)
      _RL F_ia_mult     (1:sNx,1:sNy,nITD)
      _RL QSWI_mult     (1:sNx,1:sNy,nITD)
      _RL FWsublim_mult     (1:sNx,1:sNy,nITD)

#ifdef ALLOW_DIAGNOSTICS
C     Helper variables for diagnostics
c     _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 /* ALLOW_DIAGNOSTICS */
#ifdef ALLOW_SALT_PLUME
c     d(HEFF)/dt from heat fluxes in the open water fraction of the grid cell
      _RL IceGrowthRateInLeads          (1:sNx,1:sNy)

c     The fraction of salt released in leads by new ice production there
c     which is to be sent to the salt plume package
      _RL leadPlumeFraction             (1:sNx,1:sNy)
#endif

      _RL QSWO_BELOW_FIRST_LAYER (1:sNx,1:sNy)
      _RL QSWO_IN_FIRST_LAYER (1:sNx,1:sNy)
      _RL QSWO (1:sNx,1:sNy)
      _RL QSWI (1:sNx,1:sNy)

C     Sea ice growth rates (m/s)
      _RL ActualNewTotalVolumeChange     (1:sNx,1:sNy)
      _RL ActualNewTotalSnowMelt     (1:sNx,1:sNy)

      _RL NetExistingIceGrowthRate      (1:sNx,1:sNy)
      _RL IceGrowthRateUnderExistingIce (1:sNx,1:sNy)
      _RL IceGrowthRateFromSurface      (1:sNx,1:sNy)
      _RL IceGrowthRateOpenWater        (1:sNx,1:sNy)
      _RL IceGrowthRateMixedLayer       (1:sNx,1:sNy)

      _RL EnergyInNewTotalIceVolume     (1:sNx,1:sNy)
      _RL NetEnergyFluxOutOfOcean       (1:sNx,1:sNy)

CAB   Terms for closing the budget
       _RL SItflux     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
       _RL SIaaflux    (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
       _RL SIatmFW     (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
CAB these two terms are placeholders for extra bits from advection,
C   used in budget
       _RL SIheffNeg   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
       _RL SIhsnwNeg   (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)

CAB
       _RL d_HEFFbyFlood                (1:sNx,1:sNy)

C     The energy taken out of the ocean which is not converted
C     to sea ice (Joules)
      _RL ResidualEnergyOutOfOcean      (1:sNx,1:sNy)

C     Snow accumulation rate over ice (m/s)
      _RL SnowAccRateOverIce            (1:sNx,1:sNy)

C     Total snow accumulation over ice (m)
      _RL SnowAccOverIce                (1:sNx,1:sNy)

C     The precipitation rate over the ice which goes immediately into the ocean
      _RL PrecipRateOverIceSurfaceToSea (1:sNx,1:sNy)

C     The potential snow melt rate if all snow surface heat flux convergences
C     goes to melting snow (m/s)
      _RL PotSnowMeltRateFromSurf       (1:sNx,1:sNy)

C     The potential thickness of snow which could be melted by snow surface
C     heat flux convergence (m)
      _RL PotSnowMeltFromSurf           (1:sNx,1:sNy)

C     The actual snow melt rate due to snow surface  heat flux convergence
      _RL SnowMeltRateFromSurface       (1:sNx,1:sNy)

C     The actual surface heat flux convergence used to melt snow (W/m^2)
      _RL SurfHeatFluxConvergToSnowMelt (1:sNx,1:sNy)

C     The actual thickness of snow to be melted by snow surface
C     heat flux convergence (m)
      _RL SnowMeltFromSurface           (1:sNx,1:sNy)

C     The freshwater contribution to the ocean from melting snow (m)
      _RL FreshwaterContribFromSnowMelt (1:sNx,1:sNy)

C     The freshwater contribution to (from) the ocean from melting (growing) ice (m)
      _RL FreshwaterContribFromIce      (1:sNx,1:sNy)

C     S_a : d(AREA)/dt
      _RL S_a                           (1:sNx,1:sNy)

C     S_h : d(HEFF)/dt
      _RL S_h                           (1:sNx,1:sNy)

C     S_hsnow : d(HSNOW)/dt
      _RL S_hsnow                       (1:sNx,1:sNy)

C     F_ia  - sea ice/snow surface heat flux with atmosphere (W/m^2)
C       F_ia > 0, heat loss to atmosphere
C       F_ia < 0, atmospheric heat flux convergence (ice/snow surface melt)
      _RL F_ia                          (1:sNx,1:sNy)

C     F_ia_net - the net heat flux divergence at the sea ice/snow surface
C                including sea ice conductive fluxes and atmospheric fluxes (W/m^2)
C       F_ia_net = 0, sea ice/snow surface energy balance condition met
C                     upward conductive fluxes balance surface heat loss
C       F_ia_net < 0, net heat flux convergence at ice/snow surface
C                     zero conductive fluxes and net atmospheric convergence
      _RL F_ia_net                      (1:sNx,1:sNy)

C     F_ia_net - the net heat flux divergence at the sea ice/snow surface
C                before snow is melted with any convergence (W/m^2)
C        F_ia_net < 0, some snow, if present, will melt
      _RL F_ia_net_before_snow          (1:sNx,1:sNy)

C     F_io_net - the net upward conductive heat flux through the ice+snow system
C                realized at the sea ice/snow surface
C        F_io_net > 0, heat conducting upward from ice base --> basal thickening
C        F_io_net = 0, no upward heat conduction
C                      ice/snow surface temperature > SEAICE_freeze)
      _RL F_io_net                      (1:sNx,1:sNy)

C     F_ao  - heat flux from atmosphere to ocean (W/m^2)
C        F_ao > 0
      _RL F_ao                          (1:sNx,1:sNy)

C     F_mi - heat flux from ocean to the ice (W/m^2)
      _RL F_mi                          (1:sNx,1:sNy)

      _RL FWsublim                      (1:sNx,1:sNy)

C     S_a_from_IGROW : d(AREA)/dt [from ice growth rate from open water fluxes]
      _RL S_a_IGROW                     (1:sNx,1:sNy)

C     S_a_from_IGROW : d(AREA)/dt [from ice growth rate from ocean-ice fluxes]
      _RL S_a_IGRML                     (1:sNx,1:sNy)

C     S_a_from_IGROW : d(AREA)/dt [from ice growth rate from ice-atm fluxes]
      _RL S_a_IGRNE                     (1:sNx,1:sNy)

C     Factor by which we increase the upper ocean friction velocity (u*) when

C     ice is absent in a grid cell  (dimensionless)
      _RL MLTF                          (1:sNx,1:sNy)

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C ===================================================================
C =================PART 0: constants and initializations=============
C ===================================================================

      IF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
       kSurface        = Nr
      ELSE
       kSurface        = 1
      ENDIF

C     avoid unnecessary divisions in loops
      recip_multDim        = SEAICE_multDim
      recip_multDim        = ONE / recip_multDim
C     above/below: double/single precision calculation of recip_multDim
c     recip_multDim        = 1./float(SEAICE_multDim)
      recip_deltaTtherm = ONE / SEAICE_deltaTtherm
      recip_rhoIce      = ONE / SEAICE_rhoIce

C     Cutoff for iceload
      heffTooHeavy=drF(kSurface) / 5. _d 0

C     RATIO OF SEA ICE DENSITY to SNOW DENSITY
      ICE2SNOW     = SEAICE_rhoIce/SEAICE_rhoSnow
      SNOW2ICE     = ONE / ICE2SNOW

C  Heat Flux * QI = [W] [m^3/kg] [kg/J] = [m^3/s] or [m/s] per unit area
      QI           = ONE/(SEAICE_rhoIce*SEAICE_lhFusion)
C  Heat Flux * QS = [W] [m^3/kg] [kg/J] = [m^3/s] or [m/s] per unit area
      QS           = ONE/(SEAICE_rhoSnow*SEAICE_lhFusion)

c     A reference seawater density (kg m^-3)
      RHOSW = 1026. _d 0
C     ICE LATENT HEAT CONSTANT
      lhSublim = SEAICE_lhEvap + SEAICE_lhFusion

C     regularization constants
      area_reg_sq = SEAICE_area_reg * SEAICE_area_reg
      hice_reg_sq = SEAICE_hice_reg * SEAICE_hice_reg

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
      DO bj=myByLo(myThid),myByHi(myThid)
       DO bi=myBxLo(myThid),myBxHi(myThid)

#ifdef ALLOW_AUTODIFF_TAMC
        act1 = bi - myBxLo(myThid)
        max1 = myBxHi(myThid) - myBxLo(myThid) + 1
        act2 = bj - myByLo(myThid)
        max2 = myByHi(myThid) - myByLo(myThid) + 1
        act3 = myThid - 1
        max3 = nTx*nTy
        act4 = ikey_dynamics - 1
        iicekey = (act1 + 1) + act2*max1
     &                       + act3*max1*max2
     &                       + act4*max1*max2*max3
#endif /* ALLOW_AUTODIFF_TAMC */

cph-test(
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE area = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE heff = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE hsnow = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE qnet = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE qsw = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE tices = comlev1_bibj, key = iicekey, byte = isbyte
#endif
cph-test)

C array initializations
C =====================

        DO J=1,sNy
         DO I=1,sNx

C NEW VARIABLE NAMES
          NetExistingIceGrowthRate(I,J) = 0.0 _d 0
          IceGrowthRateOpenWater(I,J)   = 0.0 _d 0
          IceGrowthRateFromSurface(I,J) = 0.0 _d 0
          IceGrowthRateMixedLayer(I,J)  = 0.0 _d 0

          EnergyInNewTotalIceVolume(I,J) = 0.0 _d 0
          NetEnergyFluxOutOfOcean(I,J)   = 0.0 _d 0
          ResidualEnergyOutOfOcean(I,J)  = 0.0 _d 0

          PrecipRateOverIceSurfaceToSea(I,J) = 0.0 _d 0

          DO IT=1,SEAICE_multDim
            ticeInMult(I,J,IT)            = 0.0 _d 0
            ticeOutMult(I,J,IT)           = 0.0 _d 0

            F_io_net_mult(I,J,IT)  = 0.0 _d 0
            F_ia_net_mult(I,J,IT)  = 0.0 _d 0
            F_ia_mult(I,J,IT)      = 0.0 _d 0

            QSWI_mult(I,J,IT)      = 0.0 _d 0
            FWsublim_mult(I,J,IT)  = 0.0 _d 0
          ENDDO

          F_io_net(I,J) = 0.0 _d 0
          F_ia_net(I,J) = 0.0 _d 0
          F_ia(I,J)     = 0.0 _d 0

          QSWI(I,J) = 0.0 _d 0
          FWsublim(I,J) = 0.0 _d 0

          QSWO_BELOW_FIRST_LAYER (I,J) = 0.0 _d 0
          QSWO_IN_FIRST_LAYER (I,J) = 0.0 _d 0
          QSWO(I,J) = 0.0 _d 0

          ActualNewTotalVolumeChange(I,J) = 0.0 _d 0
          ActualNewTotalSnowMelt(I,J)     = 0.0 _d 0

          SnowAccOverIce(I,J) = 0.0 _d 0
          SnowAccRateOverIce(I,J) = 0.0 _d 0

          PotSnowMeltRateFromSurf(I,J)    = 0.0 _d 0
          PotSnowMeltFromSurf(I,J)        = 0.0 _d 0
          SnowMeltRateFromSurface(I,J)    = 0.0 _d 0
          SurfHeatFluxConvergToSnowMelt(I,J) = 0.0 _d 0
          SnowMeltFromSurface(I,J)       = 0.0 _d 0

          FreshwaterContribFromSnowMelt(I,J) = 0.0 _d 0
          FreshwaterContribFromIce(I,J) = 0.0 _d 0

          S_a(I,J)       = 0.0 _d 0
          S_a_IGROW(I,J) = 0.0 _d 0
          S_a_IGRML(I,J) = 0.0 _d 0
          S_a_IGRNE(I,J) = 0.0 _d 0

          S_h(I,J)       = 0.0 _d 0
          S_hsnow(I,J)   = 0.0 _d 0

          MLTF(I,J)      = 0.0 _d 0

          d_HEFFbyFlood(I,J) = 0.0 _d 0

#ifdef ALLOW_SALT_PLUME
          IceGrowthRateInLeads            (I,J) = 0. _d 0
          leadPlumeFraction               (I,J) = 0. _d 0
          saltPlumeFlux             (I,J,bi,bj) = 0. _d 0
#endif

         ENDDO
        ENDDO

#if (defined (ALLOW_MEAN_SFLUX_COST_CONTRIBUTION) || defined (ALLOW_SSH_GLOBMEAN_COST_CONTRIBUTION))
        DO J=1-OLy,sNy+OLy
         DO I=1-OLx,sNx+OLx
          frWtrAtm(I,J,bi,bj)        = 0.0 _d 0
         ENDDO
        ENDDO
#endif

C =====================================================================
C  PART 1: Store ice and snow state on onset + regularize actual
C               snow and ice thickness
C =====================================================================
CAB
        DO J=1,sNy
         DO I=1,sNx
         SIheffNeg(I,J,bi,bj)=d_HEFFbyNEG(I,J,bi,bj)*SINegFac
         SIhsnwNeg(I,J,bi,bj)=d_HSNWbyNEG(I,J,bi,bj)*SINegFac
         ENDDO
        ENDDO

        DO J=1,sNy
         DO I=1,sNx

          HEFF(I,J,bi,bj) = MAX(ZERO, HEFF(I,J,bi,bj))
          AREA(I,J,bi,bj) = MAX(ZERO, AREA(I,J,bi,bj))

          IF (HEFF(I,J,bi,bj) .LE. ZERO) then
             AREA(I,J, bi,bj)  = ZERO
             HSNOW(I,J, bi,bj) = ZERO
          ELSEIF (AREA(I,J,bi,bj) .LE. ZERO) then
             HEFF(I,J,bi,bj)  = ZERO
             HSNOW(I,J,bi,bj) = ZERO
          ENDIF

          HEFFpreTH(I,J) = HEFF(I,J,bi,bj)
          HSNWpreTH(I,J) = HSNOW(I,J,bi,bj)
          AREApreTH(I,J) = AREA(I,J,bi,bj)

#ifdef ALLOW_DIAGNOSTICS
          DIAGarrayB(I,J) = AREA(I,J,bi,bj)
          DIAGarrayC(I,J) = HEFF(I,J,bi,bj)
          DIAGarrayD(I,J) = HSNOW(I,J,bi,bj)
#endif
         ENDDO
        ENDDO

#ifdef ALLOW_DIAGNOSTICS
        IF ( useDiagnostics ) THEN
         CALL DIAGNOSTICS_FILL(DIAGarrayB,'SIareaPT',0,1,3,bi,bj,myThid)
         CALL DIAGNOSTICS_FILL(DIAGarrayC,'SIheffPT',0,1,3,bi,bj,myThid)
         CALL DIAGNOSTICS_FILL(DIAGarrayD,'SIhsnoPT',0,1,3,bi,bj,myThid)
        ENDIF
#endif /* ALLOW_DIAGNOSTICS */

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE AREApreTH = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE HEFFpreTH = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE HSNWpreTH = comlev1_bibj, key = iicekey, byte = isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C     COMPUTE ACTUAL ICE/SNOW THICKNESS; USE MIN/MAX VALUES
C     TO REGULARIZE SEAICE_SOLVE4TEMP/d_AREA COMPUTATIONS

        DO J=1,sNy
         DO I=1,sNx
          IF (HEFFpreTH(I,J) .GT. ZERO) THEN
c          regularize AREA with SEAICE_area_reg
           tmpscal1 = SQRT(AREApreTH(I,J)* AREApreTH(I,J) + area_reg_sq)

c          hiceActual calculated with the regularized AREA
           tmpscal2 = HEFFpreTH(I,J) / tmpscal1

c          regularize hiceActual with SEAICE_hice_reg (add lower bound)
c           hiceActual(I,J) = SQRT(tmpscal2 * tmpscal2 + hice_reg_sq)
           hiceActual(I,J) = MAX(0.05 _d 0, tmpscal2)

c          hsnowActual calculated with the regularized AREA (no lower bound)
c          hsnowActual(I,J) = HSNWpreTH(I,J) / tmpscal1
c          actually I don't think we need to regularize this.
           hsnowActual(I,J) = HSNWpreTH(I,J) / AREApreTH(I,J)

c          regularize the inverse of hiceActual by hice_reg
           recip_hiceActual(I,J)  = AREApreTH(I,J) /
     &         sqrt(HEFFpreTH(I,J)*HEFFpreTH(I,J) + hice_reg_sq)

c         Do not regularize when HEFFpreTH = 0
          ELSE
           hiceActual (I,J)       = ZERO
           hsnowActual(I,J)       = ZERO
           recip_hiceActual(I,J)  = ZERO
          ENDIF

         ENDDO
        ENDDO

C =============================================================================
C Part 2: Precipitation as snow or rain over ice
C =============================================================================

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE SnowAccRateOverIce = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE SnowAccOverIce     = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE PrecipRateOverIceSurfaceToSea   =
CADJ &     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 */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

        DO J=1,sNy
         DO I=1,sNx
c        if we have ice and the temperature of the ice is below the freezing point
c        then the precip falls and accumulates as snow
          IF (( AREApreTH(I,J)     .GT. ZERO) .AND.
     &        ( TICES(I,J,1,bi,bj) .LT. celsius2k) ) THEN

c use either prescribed snowfall or PRECIP rate
           IF ( snowPrecipFile .NE. ' ' ) THEN
c    rate of snow accumulation in m/s over ice
c                      y [m/s] \approx 1.0 [kg/m^3] / 0.9 [m^3/kg] * x [m/s]
             SnowAccRateOverIce(I,J) = rhoConstFresh/SEAICE_rhoSnow *
     &               snowPrecip(i,j,bi,bj)

           ELSE
             SnowAccRateOverIce(I,J) = rhoConstFresh/SEAICE_rhoSnow *
     &           PRECIP(i,j,bi,bj)

           ENDIF

           PrecipRateOverIceSurfaceToSea(I,J) = ZERO

          ELSE
c            The snow/ice surface is not frozen (wet) so the precipitation
c            remains wet and runs into the ocean
             SnowAccRateOverIce(I,J) = ZERO
             PrecipRateOverIceSurfaceToSea(I,J) = PRECIP(i,j,bi,bj)
          ENDIF

c  actual change in snow thickness due to precipitation in units of
c  mean snow thickness
          SnowAccOverIce(I,J) =
     &     SnowAccRateOverIce(I,J) * SEAICE_deltaTtherm * AREApreTH(I,J)

c I,J
         ENDDO
        ENDDO

C =============================================================================
C FIND WIND SPEED
C =============================================================================

        DO j=1,sNy
         DO i=1,sNx
C ocean surface/mixed layer temperature
          TmixLoc(i,j) = theta(i,j,kSurface,bi,bj) + celsius2K
C wind speed from exf
          UG(I,J) = MAX(SEAICE_EPS, wspeed(I,J,bi,bj))
         ENDDO
        ENDDO

C =============================================================================
c           Retrieve the air-sea heat and shortwave radiative fluxes
C =============================================================================

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE qnet(:,:,bi,bj) = comlev1_bibj, key = iicekey,byte=isbyte
CADJ STORE qsw(:,:,bi,bj)  = comlev1_bibj, key = iicekey,byte=isbyte
cCADJ STORE UG = comlev1_bibj, key = iicekey,byte=isbyte
cCADJ STORE TmixLoc = comlev1_bibj, key = iicekey,byte=isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

        CALL SEAICE_BUDGET_OCEAN(
     I       UG,
     I       TmixLoc,
     O       F_ao, QSWO,
     I       bi, bj, myTime, myIter, myThid )

C =============================================================================
C      Calc air-sea fluxes in the uppermost grid cell
C =============================================================================

c--   Not all of the sw radiation is absorbed in the uppermost ocean grid cell layer.
c     Only that fraction which converges in the uppermost ocean grid cell is used to
c     melt ice.
c           SWFRACB - the fraction of incoming sw radiation absorbed in the
c                     uppermost ocean grid cell (calculated in seaice_init_vari.F)
        DO J=1,sNy
         DO I=1,sNx

c The contribution of shortwave heating is
c not included without #define SHORTWAVE_HEATING
#ifdef SHORTWAVE_HEATING
           QSWO_BELOW_FIRST_LAYER(I,J)= QSWO(I,J)*SWFRACB
           QSWO_IN_FIRST_LAYER(I,J)   = QSWO(I,J)*(ONE - SWFRACB)
#else
           QSWO_BELOW_FIRST_LAYER(I,J)= ZERO
           QSWO_IN_FIRST_LAYER(I,J)   = ZERO
#endif
           IceGrowthRateOpenWater(I,J) = QI *
     &        (F_ao(I,J) - QSWO(I,J) + QSWO_IN_FIRST_LAYER(I,J))

         ENDDO
        ENDDO

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE hsnowActual = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE hiceActual  = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE UG          = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE tices(:,:,:,bi,bj)
CADJ &     = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE salt(:,:,kSurface,bi,bj) = comlev1_bibj,
CADJ &                     key = iicekey, byte = isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C =============================================================================
c  calculate heat fluxes within ice (conduction), F_io, and across the
c  ice/atmosphere interface, F_ia
C =============================================================================

C--   Start loop over multi-categories
        DO IT=1,SEAICE_multDim

         DO J=1,sNy
          DO I=1,sNx
c record prior ice surface temperatures
           ticeInMult(I,J,IT)  = TICES(I,J,IT,bi,bj)
           ticeOutMult(I,J,IT) = TICES(I,J,IT,bi,bj)
           TICES(I,J,IT,bi,bj) = ZERO
          ENDDO
         ENDDO

c set relative thickness of ice categories
         pFac = (2.0 _d 0*IT - 1.0 _d 0)*recip_multDim
         pFacSnow = 1. _d 0

c find actual snow and ice thickness within categories categories
         IF ( SEAICE_useMultDimSnow ) pFacSnow=pFac

         DO J=1,sNy
          DO I=1,sNx
            hiceActualMult(I,J,IT)   = hiceActual(I,J) *pFac
            hsnowActualMult(I,J,IT)  = hsnowActual(I,J)*pFacSnow
          ENDDO
         ENDDO

       ENDDO /* multDim */

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE hiceActualMult = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE hsnowActualMult= comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE ticeOutMult    = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE F_io_net_mult   =
CADJ &     comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE F_ia_net_mult   =
CADJ &     comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE F_ia_mult       =
CADJ &     comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE QSWI_mult       =
CADJ &     comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE FWsublim_mult   =
CADJ &     comlev1_bibj, key = iicekey, byte = isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

c =======================================================================
c  find calculate heat fluxes within ice (conduction) and across the
c  ice/atmosphere interface for each thickness category
c =======================================================================

        DO IT=1,SEAICE_multDim
         CALL SEAICE_SOLVE4TEMP_ADX(
     I        UG, hiceActualMult(1,1,IT), hsnowActualMult(1,1,IT),
     U        ticeInMult(1,1,IT), ticeOutMult(1,1,IT),
     O        F_io_net_mult(1,1,IT),
     O        F_ia_net_mult(1,1,IT),
     O        F_ia_mult(1,1,IT),
     O        QSWI_mult(1,1,IT),
     O        FWsublim_mult(1,1,IT),
     I        bi, bj, myTime, myIter, myThid )
        ENDDO

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE hiceActualMult = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE hsnowActualMult= comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE ticeOutMult    = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE F_io_net_mult  =
CADJ &     comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE F_ia_net_mult   =
CADJ &     comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE F_ia_mult       =
CADJ &     comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE QSWI_mult       =
CADJ &     comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE FWsublim_mult   =
CADJ &     comlev1_bibj, key = iicekey, byte = isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

c =======================================================================
c  record the ice surface temperature in each category
c  and find the average of fluxes across each category

        DO IT=1,SEAICE_multDim
         DO J=1,sNy
          DO I=1,sNx

C     update TICES
            TICES(I,J,IT,bi,bj) = ticeOutMult(I,J,IT)

            F_io_net(I,J) = F_io_net(I,J) +
     &        F_io_net_mult(I,J,IT)*recip_multDim

            F_ia_net(I,J) = F_ia_net(I,J) +
     &        F_ia_net_mult(I,J,IT)*recip_multDim

            F_ia(I,J)  = F_ia(I,J)     +
     &        F_ia_mult(I,J,IT)*recip_multDim

            QSWI(I,J)  = QSWI(I,J) + QSWI_mult(I,J,IT)*recip_multDim

            FWsublim(I,J) = FWsublim(I,J) +
     &         FWsublim_mult(I,J,IT)*recip_multDim

          ENDDO
         ENDDO
        ENDDO

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE AREApreTH = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE F_io_net  = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE F_ia_net  = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE F_ia      = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE QSWI      = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE FWSublim  = comlev1_bibj, key = iicekey, byte = isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

        DO J=1,sNy
         DO I=1,sNx
c          If there is heat flux convergence at the snow surface,
c          use that energy to melt snow before melting ice.  It's
c          possible that some snow will remain after melting,
c          which will drive F_ia_net to zero, or that all of the
c          snow will be melted, leaving a nonzero F_ia_net to melt
c          some ice.
          F_ia_net_before_snow(I,J) = F_ia_net(I,J)

c         Only continue if there is snow and ice in the cell
          IF (AREApreTH(I,J) .LE. ZERO) THEN
            IceGrowthRateUnderExistingIce(I,J) = 0. _d 0
            IceGrowthRateFromSurface(I,J)      = 0. _d 0
            NetExistingIceGrowthRate(I,J)      = 0. _d 0
          ELSE
c           The growth rate (m/s) beneath existing ice is given by the upward
c           ocean-ice conductive flux, F_io_net, and QI.
            IceGrowthRateUnderExistingIce(I,J) = F_io_net(I,J)*QI

c           The rate at which snow is melted (m/s) because of surface
c           heat flux convergence.  Note, during snow melt, F_ia_net must
c           be negative (implying convergence) to make PSMRFW is positive
            PotSnowMeltRateFromSurf(I,J) = - F_ia_net(I,J)*QS

c           This is the depth of snow (m) that would be melted in one dt
            PotSnowMeltFromSurf(I,J)  =
     &        PotSnowMeltRateFromSurf(I,J) * SEAICE_deltaTtherm

c           If we can melt MORE than is actually there, then the melt
c           rate is reduced so that only that which is there
c           is melted during the time step.  In this case, not all of the
c           heat flux convergence at the surface is used to melt snow.
c           Any remaining energy will melt ice.

c           SurfHeatFluxConvergToSnowMelt is the part of the total heat
c           flux convergence which melts snow.

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CCHECK: HSNOW ACTUAL SHOULDN'T BE REGULARIZED FOR THIS
            IF (PotSnowMeltFromSurf(I,J) .GE. hsnowActual(I,J)) THEN

c           Snow melt and melt rate [m] (actual snow thickness)
              SnowMeltFromSurface(I,J)     = hsnowActual(I,J)

c             SnowMeltRateFromSurface(I,J) =
              SnowMeltRateFromSurface(I,J) =
     &           SnowMeltFromSurface(I,J) / SEAICE_deltaTtherm

              SurfHeatFluxConvergToSnowMelt(I,J) =
     &          - hsnowActual(I,J)/QS/SEAICE_deltaTtherm
           ELSE
c             In this case there will be snow remaining after melting.
c             All of the surface heat convergence will be redirected to
c             this effort.
              SnowMeltFromSurface(I,J) = PotSnowMeltFromSurf(I,J)

              SnowMeltRateFromSurface(I,J) =PotSnowMeltRateFromSurf(I,J)

              SurfHeatFluxConvergToSnowMelt(I,J) = F_ia_net(I,J)

           ENDIF

c          Reduce the heat flux convergence available to melt surface
c          ice by the amount used to melt snow
           F_ia_net(I,J) = F_ia_net(I,J) -
     &         SurfHeatFluxConvergToSnowMelt(I,J)

           IceGrowthRateFromSurface(I,J) = F_ia_net(I,J) * QI

c          The total growth rate (m/s) of the existing ice - the rate of
c          new ice accretion at the base less the rate due to surface melt
           NetExistingIceGrowthRate(I,J) =
     &            IceGrowthRateUnderExistingIce(I,J) +
     &            IceGrowthRateFromSurface(I,J)
          ENDIF

         ENDDO
        ENDDO

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE AREApreTH = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE SnowMeltRateFromSurface   = comlev1_bibj,
CADJ &                       key = iicekey, byte = isbyte
CADJ STORE IceGrowthRateUnderExistingIce  = comlev1_bibj,
CADJ &                       key = iicekey, byte = isbyte
CADJ STORE IceGrowthRateFromSurface = comlev1_bibj,
CADJ &                       key = iicekey, byte = isbyte
CADJ STORE NetExistingIceGrowthRate  = comlev1_bibj,
CADJ &                       key = iicekey, byte = isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE AREApreTH = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE theta(:,:,kSurface,bi,bj) = comlev1_bibj,
CADJ &                       key = iicekey, byte = isbyte
CADJ STORE salt(:,:,kSurface,bi,bj) = comlev1_bibj,
CADJ &                       key = iicekey, byte = isbyte
#endif
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C       Calcuate the heat fluxes from the ocean to the sea ice
        DO J=1,sNy
         DO I=1,sNx
c
c         Bound the ocean temperature to be at or above the freezing point.
          tempFrz = SEAICE_tempFrz0 +
     &              SEAICE_dTempFrz_dS * salt(I,J,kSurface,bi,bj)

          surf_theta = max(theta(I,J,kSurface,bi,bj), tempFrz)

c         inflection point
          tmpscal0 = 0.4 _d 0

c         steepness
          tmpscal1 = 7.0 _d 0

          MLTF(I,J) = ONE + (12.5 - ONE) *
     &         (ONE + EXP(  (AREApreTH(I,J) - tmpscal0)
     &                     * tmpscal1/tmpscal0))**(-ONE)

c          IF (AREApreTH(I,J) .GT. ZERO) THEN
c
c           If ice is present, MixedLayerTurbulenceFactor = 1.0, else 12.50
c            MLTF = ONE
c          ELSE
c            MLTF = 12.5 _d 0
c          ENDIF

          F_mi(I,J) = -STANTON_NUMBER * USTAR_BASE * RHOSW *
     &      HeatCapacity_Cp * (surf_theta - tempFrz) * MLTF(I,J)

          IceGrowthRateMixedLayer (I,J) = F_mi(I,J) * QI

         ENDDO
        ENDDO

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE AREApreTH   = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE F_mi        = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE IceGrowthRateMixedLayer =
CADJ &     comlev1_bibj,key = iicekey, byte = isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C       CALCULATE THICKNESS DERIVATIVES of ice (dhdt) and snow (dhsdt)
        DO J=1,sNy
         DO I=1,sNx

          S_h(I,J) =
     &       NetExistingIceGrowthRate(I,J) *  AREApreTH(I,J)
     &     + IceGrowthRateOpenWater(I,J)   * (ONE - AREApreTH(I,J))
     &     + IceGrowthRateMixedLayer(I,J)

c         Both the accumulation and melt rates are in terms
c         of actual snow thickness.  As with ice, multiplying
c         with area converts to mean snow thickness.
          S_hsnow(I,J) =  AREApreTH(I,J) *
     &      ( SnowAccRateOverIce(I,J) - SnowMeltRateFromSurface(I,J))

         ENDDO
        ENDDO

#ifdef ALLOW_SALT_PLUME
#ifdef SALT_PLUME_IN_LEADS
C       Now that we know the thickness tendency terms, we can calculate the saltPlumeFlux
        DO J=1,sNy
         DO I=1,sNx

c            It is assumed that ice production in leads can generate
c            salt plumes.  The fraction of the salt sent to the plume package
c            from ice produced in leads (defined as the open water
c            fraction) is a nonlinear function of ice concentration, AREA.
c
c            Specifically, function is a logistic curve (sigmoid) with a range and
c            domain {0,1}.  The function, f(AREA), has a single free parameter,
c            SEAICE_plumeInflectionPoint, the inflection point of the curve.
c            By construction, the function has the following properties:
c            f(1) \approx 1.0
c            f(SEAICE_plumeInflectionPoint) = 0.5
c            f(0) \approx 0.0 (when SEAICE_plumeInflectionPoint \geq 0.5)
c            f(0) > 0.0 (when SEAICE_plumeInflectionPoint < 0.5)
c
c            As AREA --> 1, the open water fraction occurs
c            in narrow leads, new ice production become spatially non-uniform,
c            and the assumptions motivating KPP no longer hold.  To treat
c            overturning in a more physically realistic way, the salt produced
c            in the leads should be sent to depth via the plume package.  To assure
c            only narrow leads generate plumes, choose a SEAICE_plumeInflectionPoint
c            of > 0.8.

c            Ensure that there is already ice present or that the total ice
c            ice tendency term is positive.  We don't want to release
c            salt if sea ice is not established in the cell.

             IF ((AREApreTH(I,J) .GT. ZERO) .OR.
     &           (S_h(I,J)           .GT. ZERO)) THEN

              leadPlumeFraction(I,J) =
     &         (ONE + EXP( ( SPinflectionPoint -
     &                       AREApreTH(I,J)
     &                     ) * 5. _d 0/
     &                     (ONE - SPinflectionPoint)
     &                   )
     &         )**(-ONE)

c            Only consider positive ice growth rate in leads for salt production
             IceGrowthRateInLeads(I,J) = max( ZERO,
     &      (ONE - AREApreTH(I,J)) * IceGrowthRateOpenWater(I,J))

              saltPlumeFlux(I,J,bi,bj) = leadPlumeFraction(I,J) *
     &            HEFFM(I,J,bi,bj)*IceGrowthRateInLeads(I,J)*
     &            ICE2WATR*rhoConstFresh*
     &            (salt(I,J,kSurface,bi,bj) - SEAICE_salt0)
             ELSE

              saltPlumeFlux(I,J,bi,bj) = ZERO

             ENDIF

         ENDDO
        ENDDO
#endif /* SALT_PLUME_IN_LEADS */
#endif /* ALLOW_SALT_PLUME */

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE AREApreTH = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE HEFFpreTH = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE S_a       = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE S_h       = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE S_hsnow   = comlev1_bibj, key = iicekey, byte = isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

c       Caculate dA/dt (S_a)
        DO J=1,sNy
         DO I=1,sNx

          S_a(I,J) =  0. _d 0

c         Caculate the ice area growth rate from the open water fluxes.
c         First, determine whether the open water growth rate is positive or
c         negative.  If positive, make sure that ice is present or that the
c         net ice thickness growth rate is positive before extending ice cover

c         this is the geometric term: Area/(2*Heff),
c         with Area/Heff regularized as Area/(Heff^2 + epsilon^2)
c         Area/(Heff^2 + ep^2) = recip_hiceActual
c         epsilon = SEAICE_hice_reg

          tmpscal0 = recip_hiceActual(I,J) / 2.0 _d 0

C    -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
C         /* IceGrowthRateOpenWater */

          S_a_IGROW(I,J) = ZERO

c         Expand ice cover if the open water growth rate is positive
          IF ( IceGrowthRateOpenWater(I,J) .GT. ZERO) THEN
            IF ((AREApreTH(I,J) .GT. ZERO)   .OR.
     &          (S_h(I,J)       .GT. ZERO))  THEN
c           Determine which hemisphere for hemisphere-dependent
c           "lead closing variable", HO
              IF ( YC(I,J,bi,bj) .LT. ZERO ) THEN
                S_a_IGROW(I,J) = (ONE - AREApreTH(I,J)) *
     &             IceGrowthRateOpenWater(I,J)/HO_south
              ELSE
                S_a_IGROW(I,J) = (ONE - AREApreTH(I,J)) *
     &             IceGrowthRateOpenWater(I,J)/HO
              ENDIF
            ENDIF

C    -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
c         Contract ice cover if the open water growth rate is negative
          ELSE
             S_a_IGROW(I,J) = tmpscal0 *
     &          IceGrowthRateOpenWater(I,J) * (ONE - AREApreTH(I,J))
          ENDIF

          S_a(I,J) = S_a(I,J) + S_a_IGROW(I,J)

C    -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
c         /* IceGrowthRateMixedLayer */
C         Contract ice if the IceGrowthRateMixedLayer is negative
          S_a_IGRML(I,J) = ZERO

          IF ( IceGrowthRateMixedLayer(I,J) .LE. ZERO)  THEN
             S_a_IGRML(I,J) = tmpscal0 * IceGrowthRateMixedLayer(I,J)
          ENDIF

          S_a(I,J) = S_a(I,J) + S_a_IGRML(I,J)

c
C    -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
C         Contract ice if the NetExistingIceGrowthRate is negative
C         /* NetExistingIceGrowthRate */
          S_a_IGRNE(I,J) = ZERO
          IF ( (NetExistingIceGrowthRate(I,J) .LE. ZERO) .AND.
     &         (HEFFpreTH(I,J)                .GT. ZERO) ) THEN

            S_a_IGRNE(I,J) =
     &       tmpscal0 * NetExistingIceGrowthRate(I,J) * AREApreTH(I,J)

          ENDIF

          S_a(I,J) = S_a(I,J) + S_a_IGRNE(I,J)

         ENDDO /* I,J */
        ENDDO /* I,J */

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE AREApreTH    = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE HEFFpreTH    = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE HSNWpreTH    = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE S_a          = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE S_h          = comlev1_bibj, key = iicekey, byte = isbyte
CADJ STORE S_hsnow      = comlev1_bibj, key = iicekey, byte = isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C       Update the area, heff, and hsnow
        DO J=1,sNy
         DO I=1,sNx
           HEFF(I,J,bi,bj)  = HEFFpreTH(I,J) +
     &           SEAICE_deltaTtherm * S_h(I,J)

           AREA(I,J,bi,bj)  = AREApreTH(I,J) +
     &           SEAICE_deltaTtherm * S_a(I,J)

           HSNOW(I,J,bi,bj) = HSNWpreTH(I,J) +
     &           SEAICE_deltaTtherm * S_hsnow(I,J)
         ENDDO
        ENDDO

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE AREApreTH        = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE HEFFpreTH        = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE HSNWpreTH        = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE AREA (:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
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 */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

        DO J=1,sNy
          DO I=1,sNx
c           Bound area, heff, and hsnow
             AREA(I,J,bi,bj)  = MIN(ONE,  AREA(I,J,bi,bj))
             AREA(I,J,bi,bj)  = MAX(ZERO, AREA(I,J,bi,bj))
             HEFF(I,J,bi,bj)  = MAX(ZERO, HEFF(I,J,bi,bj))
             HSNOW(I,J,bi,bj) = MAX(ZERO, HSNOW(I,J,bi,bj))

c            Sanity checks
             IF (HEFF(I,J,bi,bj) .LE. ZERO .OR.
     &           AREA(I,J,bi,bj) .LE. ZERO) THEN

                  AREA(I,J,bi,bj)       = 0. _d 0
                  HEFF(I,J,bi,bj)       = 0. _d 0
                  hiceActual(I,J)       = 0. _d 0
                  HSNOW(I,J,bi,bj)      = 0. _d 0
                  hsnowActual(I,J)      = 0. _d 0

             ELSE
c                 recalcuate the actual ice and snow thicknesses
                  hiceActual(I,J)  =  HEFF(I,J,bi,bj)/AREA(I,J,bi,bj)
                  hsnowActual(I,J) = HSNOW(I,J,bi,bj)/AREA(I,J,bi,bj)
             ENDIF

          ENDDO
        ENDDO

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE AREA (:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE HEFF (:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE salt(:,:,kSurface,bi,bj) = comlev1_bibj,
CADJ &                       key = iicekey, byte = isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

        DO J=1,sNy
          DO I=1,sNx

c         THE EFFECTIVE SHORTWAVE HEATING RATE
#ifdef SHORTWAVE_HEATING
            QSW(I,J,bi,bj)  =
     &         QSWI(I,J)  * (      AREApreTH(I,J)) +
     &         QSWO(I,J)  * (ONE - AREApreTH(I,J))
#else
            QSW(I,J,bi,bj) = 0. _d 0
#endif
#ifdef ALLOW_SEAICE_FLOODING
        IF(SEAICEuseFlooding) THEN

           tmpscal0 = (HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
     &              +HEFF(I,J,bi,bj)*SEAICE_rhoIce)*recip_rhoConst
          d_HEFFbyFlood(I,J) = MAX(tmpscal0 - HEFF(I,J,bi,bj), 0. _d 0)
                      HEFF(I,J,bi,bj)=HEFF(I,J,bi,bj)+d_HEFFbyFlood(I,J)
                      HSNOW(I,J,bi,bj) = HSNOW(I,J,bi,bj)-
     &                 d_HEFFbyFlood(I,J)*ICE2SNOW
c SEAICEuseFlooding
           ENDIF
c ALLOW_SEAICE_FLOODING
#endif

c           The actual ice volume change over the time step
            ActualNewTotalVolumeChange(I,J) =
     &        HEFF(I,J,bi,bj) - HEFFpreTH(I,J)

c           The net average snow thickness melt that is actually realized. e.g.
c           hsnow_orig  = 0.25 m (e.g. 1 m of ice over a cell 1/4 covered in snow)
c           hsnow_new   = 0.20 m
c           snow accum  = 0.05 m
c           melt = 0.25 + 0.05 - 0.2 = 0.1 m

c           since this is in mean snow thickness it might have been  0.4 of actual
c           snow thickness over the 1/4 of the cell which is ice covered.
            ActualNewTotalSnowMelt(I,J) =
     &                 HSNWpreTH(I,J) +
     &                 SnowAccOverIce(I,J) -
     &                 HSNOW(I,J,bi,bj)

c           The energy required to melt or form the new ice volume
            EnergyInNewTotalIceVolume(I,J) =
     &          ActualNewTotalVolumeChange(I,J)/QI
c

c     This is the net energy flux out of the ice+ocean system
c     Remember -----
c     F_ia_net : Under ice/snow surface freezing conditions,
c                vertical conductive heat flux convergence (F_c < 0) balances
c                heat flux divergence to atmosphere (F_ia > 0)
c                Otherwise, F_ia_net = F_ia (pos)
c
c     F_io_net : Under ice/snow surface freezing conditions, F_c < 0.
c                Under ice surface melting conditions, F_c = 0 (no energy flux
c                from the ice to ocean)
c
c     So if we are freezing, F_io_net = the conductive flux and there
c     is energy balance at ice surface, F_ia_net =0.  If we are melting,
c     there is a convergence of energy into the ice from above
            NetEnergyFluxOutOfOcean(I,J) = SEAICE_deltaTtherm *
     &               ( AREApreTH(I,J) *
     &                 (F_ia_net(I,J) + F_io_net(I,J) + QSWI(I,J))
     &         +     ( ONE - AREApreTH(I,J)) *  F_ao(I,J))

c     THE QUANTITY OF HEAT WHICH IS THE RESIDUAL TO THE QUANTITY OF
c     ML temperature.  If the net energy flux is exactly balanced by the
c     latent energy of fusion in the new ice created then we will not
c     change the ML temperature at all.

            ResidualEnergyOutOfOcean(I,J) =
     &             NetEnergyFluxOutOfOcean(I,J) -
     &             EnergyInNewTotalIceVolume(I,J)

C     NOW FORMULATE QNET
C     THIS QNET DETERMINES THE TEMPERATURE CHANGE
C     QNET IS A DEPTH AVERAGED HEAT FLUX FOR THE OCEAN COLUMN

            QNET(I,J,bi,bj) =
     &             ResidualEnergyOutOfOcean(I,J) / SEAICE_deltaTtherm

c    Like snow melt, if there is melting, this quantity is positive.
c    The change of freshwater content is per unit area over the entire
c    cell, not just over the ice covered bits.  This term is only used
c    to calculate freshwater fluxes for the purpose of changing the
c    salinity of the liquid cell.  In the case of non-zero ice salinity,
c    the amount of freshwater is reduced by the ratio of ice salinity
c    to water cell salinity.
            IF  (salt(I,J,kSurface,bi,bj) .GE. SEAICE_salt0 .AND.
     &           salt(I,J,kSurface,bi,bj) .GT. 0. _d 0) THEN

                 FreshwaterContribFromIce(I,J) =
     &            - ActualNewTotalVolumeChange(I,J) *
     &              SEAICE_rhoICE/rhoConstFresh
C     &             * (ONE - SEAICE_salt0/salt(I,J,kSurface,bi,bj))

            ELSE
C    If the liquid cell has a lower salinity than the specified
c    salinity of sea ice then assume the sea ice is completely fresh
                 FreshwaterContribFromIce(I,J) =
     &             -ActualNewTotalVolumeChange(I,J) *
     &              SEAICE_rhoIce/rhoConstFresh
            ENDIF
CAB
          tmpscal3 = max( 0. _d 0,
     &                    min(SEAICE_salt0,salt(I,J,kSurface,bi,bj)) )
CAB Salt in ice + ridgi - ice from snow flood
         tmpscal2 = (ActualNewTotalVolumeChange(I,J)+
     &                   SIheffNeg(I,J,bi,bj)+
     &                  0.0 )
     &                   * tmpscal3
     &                   * HEFFM(I,J,bi,bj)
     &                   / SEAICE_deltaTtherm * SEAICE_rhoIce

CAB test for getting the lines

             saltflux(I,J,bi,bj) = tmpscal2
CAB

c    The freshwater contribution from snow comes only in the form of melt
c    unlike ice, which takes freshwater upon growth and yields freshwater
c    upon melt.  This is why the the actual new average snow melt was determined.
c    In m/m^2 over the entire cell.
             FreshwaterContribFromSnowMelt(I,J) =
     &          ActualNewTotalSnowMelt(I,J)*SEAICE_rhoSnow/rhoConstFresh

c    This seems to be in m/s, original time level 2 for area
c    Only the precip and evap need to be area weighted.  The runoff
c    and freshwater contribs from ice and snow melt are already mean
c    weighted
             EmPmR(I,J,bi,bj)  = maskC(I,J,kSurface,bi,bj)*(
     &       ( EVAP(I,J,bi,bj) - PRECIP(I,J,bi,bj) )
     &          * ( ONE - AREApreTH(I,J) )
     &          - PrecipRateOverIceSurfaceToSea(I,J)*AREApreTH(I,J)
#ifdef ALLOW_RUNOFF
     &          - RUNOFF(I,J,bi,bj)
#endif
     &          - (FreshwaterContribFromIce(I,J) +
     &             FreshwaterContribFromSnowMelt(I,J)) /
     &             SEAICE_deltaTtherm ) * rhoConstFresh
     &            +( SIheffNeg(I,J,bi,bj)*SEAICE_rhoIce
     &            + SIhsnwNeg(I,J,bi,bj)*SEAICE_rhoSnow)
     &            /SEAICE_deltaTtherm * maskC(I,J,kSurface,bi,bj)

CAB Budget test , adding the energy stored in snow melt to qnet
CMoved bunch of these to line 780

       SItflux(I,J,bi,bj) = 0.0 +
     &   (AREApreTH(I,J)
     &     *(F_ia_net(I,J)
     &     + F_io_net(I,J) + QSWI(I,J))
     &  +( ONE - AREApreTH(I,J)) *  F_ao(I,J))
     &    *maskC(I,J,kSurface,bi,bj)
     &  +(SIheffNeg(I,J,bi,bj)/QI
     &  +(SIhsnwNeg(I,J,bi,bj)
     &  -ActualNewTotalSnowMelt(I,J) +SnowAccOverIce(I,J))/QS)
     &  / SEAICE_deltaTtherm * maskC(I,J,kSurface,bi,bj)

       SIaaflux(I,J,bi,bj) =  -EmPmR(I,J,bi,bj)
     &     * HeatCapacity_Cp *theta(I,J,kSurface,bi,bj)

      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
C     &     +    FWsublim(i,j)* SEAICE_rhoIce * recip_deltaTtherm
C  a_FWbySublim(
CAB

         ENDDO
        ENDDO

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE hiceActual       = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE hsnowActual      = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE AREApreTH        = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE HEFF(:,:,bi,bj)  = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE HSNOW(:,:,bi,bj) = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE AREA(:,:,bi,bj)  = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE salt(:,:,kSurface,bi,bj) = comlev1_bibj,
CADJ &                       key = iicekey, byte = isbyte
#endif /* ALLOW_AUTODIFF_TAMC */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
#ifdef ALLOW_AUTODIFF_TAMC
CADJ STORE hiceActual       = comlev1_bibj,key=iicekey,byte=isbyte
CADJ STORE hsnowActual      = comlev1_bibj,key=iicekey,byte=isbyte
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 */
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|

C Sea Ice Load on the sea surface.
C =================================
        IF ( useRealFreshWaterFlux ) THEN
         DO J=1,sNy
          DO I=1,sNx
#ifdef SEAICE_CAP_ICELOAD
           tmpscal1 = HEFF(I,J,bi,bj)*SEAICE_rhoIce
     &              + HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
           tmpscal2 = MIN(tmpscal1,heffTooHeavy*rhoConst)
#else
           tmpscal2 = HEFF(I,J,bi,bj)*SEAICE_rhoIce
     &              + HSNOW(I,J,bi,bj)*SEAICE_rhoSnow
#endif
           sIceLoad(i,j,bi,bj) = tmpscal2
          ENDDO
         ENDDO
        ENDIF

#ifdef SEAICE_DEBUG
            DO j=1,sNy
               DO i=1,sNx

               IF ( (i .EQ. SEAICE_debugPointI)   .and.
     &              (j .EQ. SEAICE_debugPointJ) ) THEN

                  print *,'ifice: myTime,myIter:',myTime,myIter

                  print  '(A,2i4,2(1x,1P3E15.4))',
     &                 'ifice i j --------------  ',i,j

                  print  '(A,2i4,2(1x,1P3E15.4))',
     &                 'ifice i j F(mi ao), rHIA  ',
     &                 i,j, F_mi(i,j), F_ao(i,j),
     &                 recip_hiceActual(i,j)

                  print  '(A,2i4,2(1x,1P3E15.4))',
     &                 'ifice i j Fi(a,ant2/1 ont)',
     &                 i,j, F_ia(i,j),
     &                 F_ia_net_before_snow(i,j),
     &                 F_ia_net(i,j),
     &                 F_io_net(i,j)

                  print  '(A,2i4,2(1x,1P3E15.4))',
     &                 'ifice i j AREA2/1 HEFF2/1 ',i,j,
     &                 AREApreTH(I,J),
     &                 AREA(i,j,bi,bj),
     &                 HEFFpreTH(I,J),
     &                 HEFF(i,j,bi,bj)

                  print  '(A,2i4,2(1x,1P3E15.4))',
     &                 'ifice i j HSNOW2/1 TMX    ',i,j,
     &                 HSNWpreTH(I,J),
     &                 HSNOW(I,J,bi,bj),
     &                 theta(I,J,kSurface,bi,bj)

                  print  '(A,2i4,2(1x,1P3E15.4))',
     &                 'ifice i j TI ATP LWD      ',i,j,
     &                 TICES(i,j,1, bi,bj) - celsius2k,
     &                 ATEMP(i,j,bi,bj) - celsius2k,
     &                 LWDOWN(i,j,bi,bj)

                  print  '(A,2i4,2(1x,1P3E15.4))',
     &                 'ifice i j S_a(tot,OW,ML,NE',i,j,
     &                 S_a(i,j),
     &                 S_a_IGROW(I,J),
     &                 S_a_IGRML(I,J),
     &                 S_a_IGRNE(I,J)

                  print  '(A,2i4,2(1x,1P3E15.4))',
     &                 'ifice i j S_a S_h S_hsnow ',i,j,
     &                 S_a(i,j),
     &                 S_h(i,j),
     &                 S_hsnow(i,j)

                  print  '(A,2i4,2(1x,1P3E15.4))',
     &                 'ifice i j IGR(ML OW ICE)  ',i,j,
     &                 IceGrowthRateMixedLayer(i,j),
     &                 IceGrowthRateOpenWater(i,j),
     &                 NetExistingIceGrowthRate(i,j)

                  print  '(A,2i4,2(1x,1P3E15.4))',
     &                 'ifice i j IVC(A ENIN)     ',i,j,
     &                 ActualNewTotalVolumeChange(i,j),
     &                 EnergyInNewTotalIceVolume(i,j)
c     &                 ExpectedIceVolumeChange(i,j),

                  print  '(A,2i4,2(1x,1P3E15.4))',
     &                 'ifice i j EF(NOS RE) QNET ',i,j,
     &                 NetEnergyFluxOutOfOcean(i,j),
     &                 ResidualEnergyOutOfOcean(i,j),
     &                 QNET(I,J,bi,bj)

                  print  '(A,2i4,3(1x,1P3E15.4))',
     &                 'ifice i j QSW QSWO QSWI   ',i,j,
     &                 QSW(i,j,bi,bj),
     &                 QSWO(i,j),
     &                 QSWI(i,j)

c                  print  '(A,2i4,3(1x,1P3E15.4))',
c     &                 'ifice i j SW(BML IML SW)  ',i,j,
c     &                 QSW_absorb_below_first_layer(i,j),
c     &                 QSW_absorb_in_first_layer(i,j),
c     &                 SWFRACB

c                  print  '(A,2i4,3(1x,1P3E15.4))',
c     &                 'ifice i j ptc(to, qsw, oa)',i,j,
c     &                 PredTempChange(i,j),
c     &                 PredTempChangeFromQSW (i,j),
c     &                 PredTempChangeFromOA_MQNET(i,j)

c                  print  '(A,2i4,3(1x,1P3E15.4))',
c     &                 'ifice i j ptc(fion,ian,ia)',i,j,
c     &                 PredTempChangeFromF_IO_NET(i,j),
c     &                 PredTempChangeFromF_IA_NET(i,j),
c     &                 PredTempChangeFromFIA(i,j)

c                  print  '(A,2i4,3(1x,1P3E15.4))',
c     &                 'ifice i j ptc(niv)        ',i,j,
c     &                 PredTempChangeFromNewIceVol(i,j)

                  print  '(A,2i4,3(1x,1P3E15.4))',
     &                 'ifice i j EmPmR EVP PRE RU',i,j,
     &                 EmPmR(I,J,bi,bj),
     &                 EVAP(I,J,bi,bj),
     &                 PRECIP(I,J,bi,bj),
     &                 RUNOFF(I,J,bi,bj)

                  print  '(A,2i4,3(1x,1P3E15.4))',
     &                 'ifice i j PRROIS,SAOI(R .)',i,j,
     &                 PrecipRateOverIceSurfaceToSea(I,J),
     &                 SnowAccRateOverIce(I,J),
     &                 SnowAccOverIce(I,J)

                  print  '(A,2i4,4(1x,1P3E15.4))',
     &                 'ifice i j SM(PM PMR . .R) ',i,j,
     &                 PotSnowMeltFromSurf(I,J),
     &                 PotSnowMeltRateFromSurf(I,J),
     &                 SnowMeltFromSurface(I,J),
     &                 SnowMeltRateFromSurface(I,J)

                  print  '(A,2i4,4(1x,1P3E15.4))',
     &                 'ifice i j TotSnwMlt       ',i,j,
     &                 ActualNewTotalSnowMelt(I,J)
c     &                 ExpectedSnowVolumeChange(I,J)

                  print  '(A,2i4,4(1x,1P3E15.4))',
     &                 'ifice i j fw(CFICE, CFSM) ',i,j,
     &                 FreshwaterContribFromIce(I,J),
     &                 FreshwaterContribFromSnowMelt(I,J)

                  print  '(A,2i4,2(1x,1P3E15.4))',
     &                 'ifice i j --------------  ',i,j

               ENDIF

               ENDDO
            ENDDO
#endif /* SEAICE_DEBUG */

C close bi,bj loops
       ENDDO
      ENDDO

        _BEGIN_MASTER( myThid )
#ifdef ALLOW_DIAGNOSTICS
      IF ( useDiagnostics ) THEN
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.
CAB
CAB        CALL DIAGNOSTICS_FILL(SItflux,

C        CALL DIAGNOSTICS_FILL(SIatmQnt,
C     &        'SIatmQnt',0,1,0,1,1,myThid)
C SIatmFW follows the same convention as empmr -- SIatmFW diag does not
CAB
C       CALL DIAGNOSTICS_FILL(SIacSubl,
C     &        'SIacSubl',0,1,3,bi,bj,myThid)

       CALL DIAGNOSTICS_FILL(SItflux,
     &        'SItflux ',0,1,0,1,1,myThid)

       CALL DIAGNOSTICS_FILL(SIaaflux,
     &        'SIaaflux',0,1,0,1,1,myThid)

        tmpscal1= - 1. _d 0
        CALL DIAGNOSTICS_SCALE_FILL(SIatmFW,
     &        tmpscal1,1,'SIatmFW ',0,1,0,1,1,myThid)
      ENDIF
#endif /* ALLOW_DIAGNOSTICS */

C        WRITE(msgBuf,'(a,a,e24.17)') 'SINegFac ',
C     &       ' = ', SINegFac
C        CALL PRINT_MESSAGE( msgBuf, standardMessageUnit,
C     &                      SQUEEZE_RIGHT, myThid )
C
        _END_MASTER( myThid )

#else  /* ALLOW_EXF and ALLOW_ATM_TEMP */
      STOP 'SEAICE_GROWTH not compiled without EXF and ALLOW_ATM_TEMP'
#endif /* ALLOW_EXF and ALLOW_ATM_TEMP */

      RETURN
      END
