_RL  hnew

        hnew  = 0.5*( (gamm-1.0)*k1local + sqrt(stuff) )

     &   (1.0 + (k1local/hnew) + (k1local*k2local/(hnew*hnew)))

        phlocal = -log10(hnew)

      _RL  hnew(nlyr)    ! new ice layer thickness (m)

         hnew(k) = hlyr

               rqh = rq * hnew(k)

                  hnew(k) = hnew(k) - etop / rq

                  hnew(k) = 0. _d 0

         qicen(k) = (hnew(k)*qicen(k)+dhi*qbot) / (hnew(k)+dhi)

         hnew(k) = hnew(k) + dhi

            if (ebot.gt.0. _d 0 .and. hnew(k).gt.0. _d 0) then

               rqh = rq * hnew(k)

                  hnew(k) = hnew(k) - ebot / rq

                  hnew(k) = 0. _d 0

         hi = hi + hnew(k)

            esurp = esurp - rhoi*qicen(k)*hnew(k)

               if (dhi.ge.hnew(k)) then

                 evap=evap-hnew(k)*rhoi/dt

                 esurp = esurp - hnew(k)*rhoi*(qicen(k)-Lfresh)

                 hnew(k)=0. _d 0

                 hq = hnew(k)*qicen(k)-dhi*Lfresh

                 hnew(k) = hnew(k) - dhi

                 qicen(k)=hq/hnew(k)

         hi = hi + hnew(k)

            esurp = esurp - rhoi*qicen(k)*hnew(k)

         rqh = rhoi*qicen(1)*hnew(1) + rhos*qsnow*dhs

         hnew(1) = hnew(1) + dhi

         qicen(1) = rqh / (rhoi*hnew(1))

            hnew(k)=hnew(k)-chi/2. _d 0

         hi = hi + hnew(k)

     &                                                 hnew, qicen

     I                            hlyr, hnew, myThid )

     I                                 hlyr, hnew, myThid )

      _RL  hnew(*)

      if (hnew(1).gt.hnew(2)) then

         f1 = (hnew(1)-hlyr)/hlyr

            qhtot = hnew(1)*qicen(1) + hnew(2)*qicen(2)

         f1 = hnew(1)/hlyr