_RL  hnew(nlyr)

         hnew(k) = hlyr

               rqh = rq * hnew(k)

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

                  hnew(k) = 0. _d 0

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

     &        (hnew(nlyr)+dhi)

         hnew(nlyr) = hnew(nlyr) + 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)

     &  'ThSI_CALC_TH: b-Winton: hnew, qice =', 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

        _RL  hnew

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

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

        phlocal = -log10(hnew)

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

        _RL  hnew

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

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

        phlocal = -log10(hnew)